Abstract (top)
A simple threepoint vertical grid test is developed to show that the arithmetic intergrid hydraulic conductivity mean does not have mathematical physicality for unsaturated flow in a porous medium with an exponential pressureconductivity relation in all but trivial special cases. Compared to the analytic Darcian mean for this medium, the arithmetic mean can overestimate hydraulic conductivity for dryoverwet conditions and underestimate conductivity for wetoverdry conditions. This can give rise to nonphysical separation and clumping of the mass flow, and can produce astable behavior in a model of a long, vertical fracture. A numerical experiment shows that in a vertical unsaturated flow model with space steps significantly larger than the "displacement pressure", yd, the arithmetic mean produces dispersive oscillations in the leading edge of an initialcondition wet pulse and a dry spike of excessively increased matric suction at the trailing edge. The trailing dry spike tends to increase as the logarithm of time in the model, out to approximately 2.63 years of model time, and can be delayed in onset by transient conditions. This error does not bode well for such models attempting to the predict flow of water and solutes on the scale of thousands of years. By contrast, a parallel model using the analytic Darcian mean demonstrates smooth, wellbehaved response at all the model time and space step sizes.
Introduction (top)
In the field of modeling unsaturated ground water flow, a numerous virtues are ascribed in the name of , or on the basis of using, "modern massconservative methods", as best described by Celia, et al. (1990). A developer of such methods may dismiss efforts to improve time step discretization error as appearing to "suffer from the problems associated with hbased methods" (personal communication, delivered peer review, postmarked June 11, 1993). A consultant in the nuclear waste business may claim to "routinely get mass balances of less than 0.001percent at long times (1000 years) on work I do for the Nevada Test Site and Hanford Tank Initiative" (personal communication, delivered in peer review, dated March 19, 1998).
The implication in these cases is that, if one only uses a method that gives such good mass balance over such long times, there are no other significant errors left to consider. This brings to modeling a certitude that may allow a contractor for the Yucca Mountain nuclear waste site to make sweeping statements about the effective reach of such models, without finding his claims and assumptions seriously questioned by regulatory agencies: "And our method of breaking this down was, in fact, to take a drop of water and walk it through the mountain, and a drop of water in a rain cloud until it hits a receptor  what happens to it, physically and in a process manner as it walks through the mountain" (Mr. Jack Bailey, Director of Regulatory and Licensing for the M&O operating contractor at YMP, from the transcript of the 105th Meeting of the NRC Advisory Committee on Nuclear Waste, Dec 16, 1998, Rockville MD).
This certitude may be misplaced. The fact that most of the modeling error may seem to go away when a better method is used does not guarantee that it all did. In nonlinear models it is common to find that errors are more than additive. Even in linear models, one need only look at time and spacestep discretization errors. The linear components are dominant only when the step sizes are small. Past certain step sizes, the significant errors increase as powers of the step sizes. Even in "modern massconservative" models, errors other than mass balance remain very much alive and well and waiting to increase in favorable regimes if left unaccounted.
ThreePoint Test of Physicality (top)
Consider the vertical threepoint grid in a homogeneous porous medium in the general case in Figure 1, of points (x0,y0), (x1,y1) and (x2,y2), where x is vertical position (m) positive upwards, x0 = 0, x1 = Dx and x2 = 2 Dx, y is matric suction head (m), H = xy is total head (m), Ks is the saturated conductivity (m/s) and the unsaturated hydraulic conductivity is k(y) = Ks·kr(y) (m/s). The relative hydraulic conductivity is kr(y), which is dimensionless. Let y0 and y2 be constantpressure boundary conditions, resulting in a steadystate flow, q (m/s), as in [1], where Km (dimensionless) is the mean relative hydraulic conductivity in general between x0 and x2 and rq (dimensionless) is the relative flow. This is a steadystate problem in a homogeneous fracture. We also know that the flow between x0 and x2 is equal to the flow between x1 and x2 and the flow between x0 and x1, as in [2], where kmu is the mean hydraulic conductivity between x1 and x2 and kml is the mean between x0 and x1.
From this we can derive [3], the relationship between the mean hydraulic conductivities. Take each H1 term in [2] and solve for H1 against the Km term. Then set the H1 results equal and divide out the factors of (H2H0). Notice that this does not say that Km is the harmonic mean of kr0 = kr(y0) and kr2 = kr(y2), but the harmonic mean of the means of the Dx segments. If a method of calculating Km does not satisfy [3], then it cannot correctly determine the constant, steadystate flow, q, or the relative flow, rq.
[1]
[2]
[3]


Because this simple mathematical test [3] depends upon Darcy's law and the pressureconductivity relation of the porous medium, kr(y), it is also a test of physicality. As we shall see, if a method of estimating the intergrid hydraulic conductivity mean, Km, cannot satisfy this test, it cannot correctly bring a perturbed steadystate system back to the properly predicted physical equilibrium.
Consider the threepoint system in Figure 1, with the center point perturbed from steadystate equilibrium, keeping the boundary conditions at x0 and x2 constant. Perturbed, it becomes a transient flow problem. In the finite difference expression of Richards' equation [4], the nonequilibrium volumetric water content at x1, q1 (dimensionless), relaxes with each time step, Dt, from the perturbed value, until the quantity in brackets on the righthandside become zero as t goes to infinity. Depending on the form of kmu and kml, this equilibrium value of q1 may or may not agree with a real physical case of steadystate flow. Note how this quantity relates to the righthandside of equation [2]. If one sets the righthandside to zero and solves for y1, this predicts the eventual equilibrium point of the transient flow problem. One can also do the same thing with equation [3], which must provide the same estimate of y1 for the equilibrium point.
[4]
At steadystate equilibrium, this threepoint, finite difference system is also the expression of a Laplace equation [5], derived from Darcy's law and the steadystate constant flow condition. As such, the solution of y(x) must satisfy the minimummaximum principle for elliptic boundary value problems (Duchateau and Zachmann, 1986). This means that every value y1(x1) must be between the boundary conditions y0 and y2, y1 Î [y0,y2]. If it violates this principle, the method of estimating the mean that produces it, Km, kml and kmu, is both nonmathematical and grossly nonphysical. This is a necessary condition, but not sufficient. As we shall see, the minmax principle can be satisfied when [3] is not, or when satisfying [3] and [4] separately do not produce the same equilibrium value of y1. Note that since the steadystate equilibrium point involves no massrateofchange, the pressuresaturation relation, q(y), required for [4] is not necessary in [2], [3] or [5].
[5]
A Comparison of Two Means (top)
Suppose for the moment that Km is the arithmetic mean, Ka = (kr0+kr2)/2, and that kml = (kr0+kr1)/2 and kmu = (kr1+kr2)/2, and that the homogenous porous medium in Figure 1 is characterized by the pressureconductivity relation in [6]. This crudely approximates a fracture in Topopah Spring welded volcanic tuff (Schenker, et al., 1995), but assuming the tuff is impervious to flow. If one solves [3] for kr1, one gets [7]. But if one uses the equilibrium point of the righthandside of [4] and solves [8] for y1, one gets [9]. The only way that the equilibrium values of y1 calculated from [7] and [9] can be the same is for the term involving Dx in [9] to be zero. This can only be true if Dx = 0 or kr0 = kr2.
[6] , where h = 6.4 (1/m), and yd = 0.08 (m), both are fitting parameters for a medium crudely approximating a fracture in impervious rock.
[7]
[8]
[9]
These two cases are trivial, and not even the case for Dx = 0 is physically correct for this medium, as we shall see. Thus, the arithmetic mean cannot produce the correct steadystate equilibrium point in general for this porous medium, nor for many others. Notice that the righthandsides of equations [2] and [4] have a dependence on Dx through the total heads, H1 and H2. This implies that the mean that satisfies [3] must also have a dependence on Dx. None of the commonlyused simple means, such as the arithmetic, geometric or upstream, have this dependence.
Table 1 shows solutions to equation [4] as time goes to infinity using the arithmetic mean, kml = (kr0+kr1)/2, kmu = (kr1+kr2)/2, and Km = (kr0+kr2)/2, with y1 = 2 m and y2 = 3 m. The kr1 column is calculated by substituting y1 in to [6]. In rows 3, 4, 6 and 7 the y1column is calculated by iterating [8]. In the remaining rows, y1 is set to y2 and Dx is iterated in [4] to a solution to obtain the boundary value of y1 for violation of the minmax principle. It is apparent that above a certain value of Dx for this problem y1 will settle to a value outside the range of [y0,y2] in violation of the minmax principle. The last column is the righthandside of [3], the harmonic mean of kml and kmu, which shows a violation of consistency with Km and Darcy's law in every row. This means that while the flows calculated with kml and kmu in the Dxsegments are equal, they do not equal the flow calculated with Km in the 2·Dxsegment. Nothing using the arithmetic mean truly balances for the steadystate solution in this medium, except for the trivial case of y0 = y1 = y2 or kr0 = kr1 = kr2.
The variables h, Dx and y1 are as described above, with y0 = 2 m and y2 = 3 m. The variable, y1, is determined by the solution to equation [8], kr1=kr(y1) = exp(h·(ydy1)), kml = (kr0+kr1)/2, kmu = (kr1+kr2)/2, Km = (kr0+kr2)/2 and the mean of means is 2·kml·kmu/(kml+kmu), the righthandside of equation [3]. The y1 column tests conformity to the minmax principle, that y1 is included in the range [y0,y2]. Rows 6 and 7 show violation of the minmax principle. Rows 1,2 and 5 show the boundary of violation for the minmax principle. The last two columns on the right test the balance of equation [4] for the arithmetic mean, which fails in every row.









1 








2 








3 








4 








5 








6 








7 








The analytic solution (Baker, et al., 1999a) for [5] and [6] together is [10], expressed as kmu. The form is exactly the same for Km and kml, with appropriate substitutions for Dx, kr1 and kr0. This form satisfies both [2] and [3] by definition. Because it satisfies Darcy's law and the steadystate constant flow condition [5] over any range of x for this type of porous medium, this kind of mean is called a Darcian mean (Baker, 1994, 1995). The equilibrium value of y1 is shown in equation [11], and is considerably different from the form in [9]. In the limit, as Dx goes to zero, kml goes to (kr1kr0)/ln(kr1/kr0). The corresponding limit for y1 is apparent from [11].
[10] , u = h·Dx
[11]
Table 2 shows results similar to Table 1, except that the Darcian mean for the conductivity relation [6] is used for Km, kml and kmu. Let the Darcian mean in [10] be Kv(kr1, kr2, Dx). Then Km = Kv(kr0, kr2, 2·Dx), kml = Kv(kr0, kr1, Dx) and kmu = Kv(kr1, kr2, Dx). In every row of Table 2, Dx is the same as in Table 1 and y1 is iterated to a solution of setting the righthandside of [4] to zero. In contrast to Table 1 and the arithmetic mean, here every row shows no violation of the minmax principle, or of consistency between Km and the harmonic mean of kml and kmu in the last column. In every case, the flows in both Dxsegments and their combination, calculated with kml, kmu and Km are equal, by the definition of Darcian means. There is no logical reason to suppose that if a method of estimating intergrid conductivity mean fails these tests for steadystate flow that it will do better in transient flow.
The variables h, Dx and y1 are as described above, with y0 = 2 m and y2 = 3 m. The variable, y1, is determined by the solution to equation [11], or setting the righthandside of [4] to zero, kr1=kr(y1) = exp(h·(ydy1)), kml = Kv(kr0,kr1,Dx), kmu = Kv(kr1,kr2,Dx), Km = Kv(kr0,kr2,2·Dx) and the mean of means is 2·kml·kmu/(kml+kmu), the righthandside of equation [4]. The y1 column tests conformity to the minmax principle, that y1 is included in the range [y0,y2]. The last two columns on the right test the balance of equation [3]. In each case, the minmax principle and Darcy's law are perfectly preserved.









1 








2 








3 








4 








5 








6 








7 








Now consider just the upper segment of Figure 1, from x1 to x2. Let the Darcian mean in [10] be designated by Kv, and let the arithmetic mean be designated by Ka. Note that if both means are divided by the geometric mean, Kg = (kr1·kr2)^{1/2}, they both become functions of the ratio, kr1/kr2. It happens that in the log^{3}space of (log(kr1), log(kr2) and log(Kv or Ka)), that log(Kg) is parallel to the kr1=kr2 diagonal and log(kr1/kr2) is perpendicular to it. Dividing both means by Kg effectively collapses the 3D log^{3}space into a 2D space in which the kr1=kr2 diagonal is viewed endon as a point.
The loglog plots in Figures 2a,b illustrate this. Kv/Kg plots as lines for Dx = 0.04, 0.08, 0.16, 0.32, 0.64, and 1.28 m in 2a and Dx = 1.28, 1.829, 2.56, 3.657, 5.12, 8.533 and 12.8 m in 2b. Ka/Kg plots as circles. Wetoverdry conditions are left of the center point and dryoverwet are on the right. Note that as Dx goes to zero, so does u in [10], to the limit of Kv = (kr1kr2)/ln(kr1/kr2), a symmetric curve in this space. But as Dx goes to +infinity, Kv goes to kr2. For the special and very trivial case of Dx goes to zero and kr1 goes to kr2, Kv goes to Ka, as the crossing point in the center of the plots suggests.




These plots thus indicate some consequences for an unsaturated flow model. If the water profile has a wetoverdry slope, then the flow above that slope may be retarded below normal by use of the arithmetic mean for large Dx. But if the slope is dryoverwet, then the flow below the slope will be too fast. This creates a situation where the vertical mass flow can be artificially separated into nonphysical clumps purely by the numerical characteristics of the intergrid hydraulic conductivity mean.
Numerical Experiment (top)
Now consider a numerical experiment for the fracture described by pressureconductivity relation [6] and the pressuresaturation relation [12] and their parameters given in this paper. A homogeneous fracture of depth, xl = 512 m, is modeled with a finite difference form of the Richards' unsaturated flow equation. The modified Picard method (Celia, et al., 1990) is used with an adaptive time step (Baker, et al., 1998). In this case, just to make x equal to depth, x will be positive downwards. The upper boundary is set at a matric suction (2.95823 m) such that kr = 10^{8}. The lower boundary has a noflow (impermeable) condition and is set up such that the depth, xl, is constant no matter how many grid points, 0 to np, are used in the model.
[12] , where q = volumetric water content (m^{3}/m^{3}, effectively dimensionless), qr = residual water content, qs = saturated water content, and b and yd are fitting parameters
The initial condition of all the grid points in the model will be a matric suction such that kr = 10^{8}, except for the points from 0.35*xl to 0.45*xl, which shall be set to a positive pressure head of 1 m (y = 1 m). The number of points, np, in the model will be an even multiple of 10, such as, 40, 60, 100, 140, 200, 280, 400, giving space steps of Dx = 12.8, 8.533, 5.12, 3.657, 2.56 1.829 and 1.28 m for xl = 512 m. The same errors occur in smaller reaches, but this larger scale has been chosen to make the dispersive nature of some oscillations more apparent.
Because the mass inflow and outflow at the boundaries of this experiment will be orders of magnitude smaller than the mass flow in the interior, the relative global mass balance, as defined by Celia, et al., will not be used. Instead, each time step will be calculated to converge to an error, rmb, in the equivalent depth of water of 10^{10} m. The error, rmb = sumth  flx, where sumth = is the trapezoidal integral of the mass change in the model during one time step, and flx = Km_{1/2}*(H1H0)*Dt/Dx is the mass flow into the upper boundary of the model in one time step. Km_{1/2} is the intergrid hydraulic conductivity mean between the x0 = 0 upper boundary and the first grid point at x1 = Dx.
In this case, two models will be run in parallel for comparison. One uses the arithmetic mean, Ka, for the intergrid hydraulic conductivity mean. The other uses the Darcian mean in [10], Kv. The model using the Darcian mean will be set to adjust the time steps so that it converges to rmb < 10^{10} m in 40 iterations or less. The model using the arithmetic mean uses exactly the same time steps, but is allowed to converge in 80 iterations or less. If the either model does not converge in the allotted number of iterations, then the time step is reset and reduced, and both models are rerun for that time step. If the Darcian mean model converges in less than 40 iteration, the time step is increased slightly. The maximum time step is limited to 5(10^{5}) s. In this way, any difference between the two models involving time step as well as space step discretization error is removed. Both models are equally affected.
Figure 3a shows the results for the arithmetic mean model for np = 40, Dx = 12.8 m, at t = 0, 10, 0.512(10^{6}), 1.024(10^{6}), 2.048(10^{6}) and 4.096(10^{6}) s. Figure 3b shows the results in the same run for the Darcian mean model. Note the different vertical scales, which are necessary due to the radically different responses. This is the model of an initialcondition pressure pulse that should be both relaxing and drifting downwards in the fracture. In the first ten seconds, the arithmetic mean model produces excess negative matric suction (positive pressure) heads that are nonphysical. The Darcian mean model, by contrast, relaxes completely to the just under saturation near a matric suction of 0.08 m, with no apparent change in pulse shape.




As the models progress in time, the arithmetic mean model develops a persistent nonphysical spike in matric suction at the top edge of the pulse (left on graph). This is apparently a direct result of violating the minmax principle, as demonstrated in the threepoint grid test in the analysis above. The arithmetic mean model also develops severe oscillations in the peak of the pulse, producing many nonphysical peaks that are consistent with the concept of mass clumping due to a differential error in hydraulic conductivity between wetoverdry and dryoverwet conditions. By contrast, the Darcian mean model is very smooth and wellbehaved.




Figures 4a and b show the convergence of the arithmetic mean and Darcian mean models for np = 40, 100, 200 and 400, or Dx = 12.8, 5.12, 2.56 and 1.28 m, at t = 4.096(10^{6})s. As the step size decreases, the trailing edge suction spike and the leading edge oscillations in the arithmetic mean model decrease. Both models converge to the same finegrid solution, but the Darcian mean model shows superior error and stability characteristics.
NonDarcian flow errors are not apparent in this example for vertical space step sizes below where the arithmetic mean actually violates the minmax principle. But in another example (Baker, 1999b), of infiltration into a fracture to less then 1 cm, with space steps from 1.5 mm to 21 mm, and an adaptive grid set to maintain a constant ratio between adjacent grid conductivities, using the arithmetic mean produced errors in the wetting front position of up to 18.75%, compared to 0.36% for an approximate Darcian mean. It may be that nonDarcian flow errors are tolerable in many cases, as the widespread use and popularity of the arithmetic mean suggest, but this cannot be certified unless they are actually accounted.
The oscillations in the leading edge of the pulse in the arithmetic mean model are reminiscent of numerical dispersion in hyperbolic systems. But classical numerical dispersion is created by the differing speeds of propagation of different frequency components of the pulse. Here the differing speeds of propagation are generated directly by errors in the intergrid conductivity mean, and depend as well on the slope of the pulse. This kind of oscillation has been seen previously in fracture flow infiltration using a vanGenuchtenstyle conductivity relation in Baker, et al. (1999a).
Figure 5 shows how the matric suction spike evolves as a function of time and model vertical space step size, Dx. The trend, out to 83,886,100 s in model time, is for the nonphysical spike to increase logarithmically in time, once it starts to develop. The plot for np = 40, Dx = 12.8 m, is atypical, possibly because of increasing space step discretization error. Note that the plots for 1.829 and 2.56 m start to decrease before rising above the initial conditions behind the pulse. The reasons for delayed onset and the apparent logarithmic increase of the dry spike are not fully understood at this time.


Although it does not show well, note that even for np = 400, Dx = 1.28 m, the grid point at the trailing edge of the initial pulse rises from the initial condition of 2.95823 m to 2.99623 m at the end of the run. There is no physical reason for it to do so; the gravity flow into the fracture is the same as in the fracture to the top edge of the relaxing pulse. If the pulse were diffusing upwards, the trend would be in the opposite direction. If the pulse had reached the noflow lower boundary and the fracture were filling with water, due to the upper boundary inflow of 4.74(10^{10}) m/s, the trend would be in the opposite direction. The model end time is about 2.63 years, and the nonphysical spike for the Dx = 1.28 m case is just beginning to show. This does not bode well for simulation models that use the arithmetic mean, or any other significantly nonDarcian mean, to predict flow over scales of thousands of years.
Conclusions (top)
Previous work (Baker, et al., 1999a,b) has defined Darcian intergrid hydraulic conductivity means as those necessary to the mathematically and physically correct solution of steadystate unsaturated flow problems, accounting both for the unsaturated flow characteristics of the porous medium and the model vertical space step size (for models with gravity). By definition, a Darcian mean provides the correct intergrid mean at all vertical scales. This paper demonstrates the principle that if a method of estimating an intergrid unsaturated conductivity mean does not correctly solve for the steadystate equilibrium, then it cannot return a perturbed system to a physically correct steadystate equilibrium. Further, it may even violate the minimummaximum principle for elliptic boundary value problems.
This papers uses a simple threepoint grid test to show the nonphysicality of the arithmetic mean in regard to a trial homogeneous porous medium for which an analytic form of the Darcian mean is known to exist. Unlike the Darcian mean, the arithmetic mean, like many simple means, has no provision to account for a model's vertical space step size. Except for trivial cases, such as pure gravity flow, the arithmetic mean is shown to violate basic mathematical rules of flow. At vertical space step sizes significantly larger than the "displacement pressure", yd, it is shown to violate the minimummaximum principle as well. It may be that in many cases the nonDarcian flow errors due to simple means, like the arithmetic mean, may be tolerable, as demonstrated by their wide spread use and popularity. But unless those errors are actually accounted, their tolerability cannot be legitimately certified.
A relaxing wet pulse in a long, unsaturated, vertical fracture is modeled in parallel models using the Darcian mean and the arithmetic mean with exactly the same initial and boundary conditions and time steps. At the larger vertical space step sizes, where the arithmetic mean violates the minmax principle in steadystate flow, it produces nonphysical oscillations similar to those of numerical dispersion. It also shows the ability to generate nonphysical dry and wet spikes due to its underestimation of flow for wetoverdry conditions and overestimation of flow for dryoverwet conditions. Although the reasons are not yet well understood, the onset of the dry spike above the wet pulse can be delayed for smaller space step sizes, and can apparently grow as the logarithm of time once it starts. By contrast, the model using the Darcian mean is smooth and wellbehaved at the space step sizes studied. It produces no apparent nonphysical behavior beyond that reasonably expected from space and time step discretization errors.
The does not bode well for models attempting to predict flow on the scales of thousands of years while using nonDarcian means. The apparently nonlinear, astable nature of nonDarcian flow error and the delay of onset due to changing transient conditions may make this kind of error very difficult to separate from other errors and effects. There is no longer any question that nonDarcian flow errors exist. The only question now is how much and how bad. The simple tests provided in this paper begin to give an indication of this, but more study is needed to fully understand the problem.
References (top)
Baker, D.L., 1994. Improved algorithms for finite difference modeling of Richards' equation. Ph.D. Dissertation in Soil Physics, Agronomy Dept., Colorado State University, Ft. Collins, CO 80523, Summer 1994; Also available through University Microfilms, Inc.
Baker, D.L. 1995. Darcian weighted interblock conductivity means for vertical unsaturated flow. GROUND WATER, v. 33, no. 3, pp. 385390.
Baker, D.L., M.E. Arnold and H.D. Scott, 1998, Some comparisons of new time steps for massconservative infiltration models, not yet published, draft available at: http://www.aquarien.com.
Baker, D.L., M.E. Arnold and H.D. Scott. 1999a. Some analytical and approximate Darcian means, GROUND WATER, v.37, no.4, pp. 532538. (click here for preprint draft)
Baker, D.L., 1999b, Some methods to examine the validity of Darcian means, removed from site Sept 2004.
Celia, M.A., E.T. Bouloutas and R.L. Zarba. 1990. A general massconservative numerical solution for the unsaturated flow equation. Water Resources Research 26 (7), 14831496.
DuChateau, P. and D.W. Zachmann. 1986. Schaum's Outline of Theory and Problems of Partial Differential Equations. McGrawHill, New York, 241p, ISBN 0070178976.
Schenker, A.R., D.C. Guerin, T.H. Robey, C.A. Rautman and R.W. Bernard. 1995. Stochastic hydrogeological units and hydrogeologic properties development for totalsystem performance assessments, Sandia Report SAND940244 * UC814, Prepared by Sandia National Laboratories, Albuquerque NM 87185 and Livermore CA 84550 for the United States Department of Energy under Contract DEAC0494AL85000, 87p + appendices.