Dispersive Errors Induced by a Non-Darcian Mean

in a Model of Unsaturated Flow

Topics

• Abstract
• Introduction
• Three-Point Test of Physicality
• A Comparison of Two Means
• Numerical Experiment
• Conclusions
• References

Abstract (top)

A simple three-point 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 pressure-conductivity relation in all but trivial special cases. Compared to the analytic Darcian mean for this medium, the arithmetic mean can overestimate hydraulic conductivity for dry-over-wet conditions and underestimate conductivity for wet-over-dry conditions. This can give rise to non-physical 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 initial-condition 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, well-behaved 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 mass-conservative 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 h-based 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.001-percent 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 space-step 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 mass-conservative" models, errors other than mass balance remain very much alive and well and waiting to increase in favorable regimes if left unaccounted.

Three-Point Test of Physicality (top)

Consider the vertical three-point 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 = x-y 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 constant-pressure boundary conditions, resulting in a steady-state 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 steady-state 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 (H2-H0). 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, steady-state flow, q, or the relative flow, rq.

Figure 1: Three-Point Grid Test Points (x0,y0), (x1,y1) and (x2,y2) in a homogeneous, unsaturated porous medium where x0 = 0, x1 = Dx, x2 = 2 Dx. The total hydraulic head is H = x-y (m). The relative mean hydraulic conductivities are Km between x0 and x2, kml between x0 and x1 and kmu between x1 and x2, satisfying equations [1], [2] and [3]

Because this simple mathematical test [3] depends upon Darcy's law and the pressure-conductivity 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 steady-state system back to the properly predicted physical equilibrium.

Consider the three-point system in Figure 1, with the center point perturbed from steady-state 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 non-equilibrium volumetric water content at x1, q1 (dimensionless), relaxes with each time step, Dt, from the perturbed value, until the quantity in brackets on the right-hand-side 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 steady-state flow. Note how this quantity relates to the right-hand-side of equation [2]. If one sets the right-hand-side 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.

At steady-state equilibrium, this three-point, finite difference system is also the expression of a Laplace equation [5], derived from Darcy's law and the steady-state constant flow condition. As such, the solution of y(x) must satisfy the minimum-maximum 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 non-mathematical and grossly non-physical. This is a necessary condition, but not sufficient. As we shall see, the min-max 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 steady-state equilibrium point involves no mass-rate-of-change, the pressure-saturation relation, q(y), required for [4] is not necessary in [2], [3] or [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 pressure-conductivity 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 right-hand-side 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.

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 steady-state equilibrium point in general for this porous medium, nor for many others. Notice that the right-hand-sides 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 commonly-used 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 y1-column 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 min-max 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 min-max principle. The last column is the right-hand-side 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 Dx-segments are equal, they do not equal the flow calculated with Km in the 2·Dx-segment. Nothing using the arithmetic mean truly balances for the steady-state 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·(yd-y1)), kml = (kr0+kr1)/2, kmu = (kr1+kr2)/2, Km = (kr0+kr2)/2 and the mean of means is 2·kml·kmu/(kml+kmu), the right-hand-side of equation [3]. The y1 column tests conformity to the min-max principle, that y1 is included in the range [y0,y2]. Rows 6 and 7 show violation of the min-max principle. Rows 1,2 and 5 show the boundary of violation for the min-max principle. The last two columns on the right test the balance of equation [4] for the arithmetic mean, which fails in every row.

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 steady-state 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 (kr1-kr0)/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 right-hand-side of [4] to zero. In contrast to Table 1 and the arithmetic mean, here every row shows no violation of the min-max 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 Dx-segments 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 steady-state 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 right-hand-side of [4] to zero, kr1=kr(y1) = exp(h·(yd-y1)), 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 right-hand-side of equation [4]. The y1 column tests conformity to the min-max 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 min-max principle and Darcy's law are perfectly preserved.

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³-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 3-D log³-space into a 2-D space in which the kr1=kr2 diagonal is viewed end-on as a point.

The log-log 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. Wet-over-dry conditions are left of the center point and dry-over-wet are on the right. Note that as Dx goes to zero, so does u in [10], to the limit of Kv = (kr1-kr2)/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.

Figure 2a: Ka/Kg and Kv/Kg versus kr1/kr2 The circles are Ka/Kg. The lines are Kv/Kg for Dx (or Dx) = 0.04, 0.08, 0.16, 0.32, 0.64 and 1.28 m	Figure 2b: Ka/Kg and Kv/Kg versus kr1/kr2 The circles are Ka/Kg. The lines are Kv/Kg for Dx (or Dx) = 1.28, 1.829, 2.56, 3.657, 5.12, 8.533 and 12.8 m

These plots thus indicate some consequences for an unsaturated flow model. If the water profile has a wet-over-dry 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 dry-over-wet, then the flow below the slope will be too fast. This creates a situation where the vertical mass flow can be artificially separated into non-physical 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 pressure-conductivity relation [6] and the pressure-saturation 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 no-flow (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.

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*(H1-H0)*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⁵) 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⁶), 1.024(10⁶), 2.048(10⁶) and 4.096(10⁶) 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 initial-condition 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 non-physical. 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.

Figure 3a: Arithmetic Mean Model Initial pulse of y = -1 m between 0.35xl and 0.45xl, relaxing and translating with time, for np = 40 at t = 0, 10, 0.512(106), 1.024(106), 2.048(106) and 4.096(106) s, in the example fracture	Figure 3b: Darcian Mean Model for np = 40 running in parallel with the same time steps in the same fracture and the same output times.

As the models progress in time, the arithmetic mean model develops a persistent non-physical spike in matric suction at the top edge of the pulse (left on graph). This is apparently a direct result of violating the min-max principle, as demonstrated in the three-point grid test in the analysis above. The arithmetic mean model also develops severe oscillations in the peak of the pulse, producing many non-physical peaks that are consistent with the concept of mass clumping due to a differential error in hydraulic conductivity between wet-over-dry and dry-over-wet conditions. By contrast, the Darcian mean model is very smooth and well-behaved.

Figure 4a: Arithmetic Mean Model Convergence to the fine-grid solution, for np = 40, 100, 200 and 400, or Dx = 12.8, 5.12, 2.56 and 1.28 m.	Figure 4b: Darcian Mean Model Convergence to the fine-grid solution, for np = 40, 100, 200 and 400, or Dx = 12.8, 5.12, 2.56 and 1.28 m. Parallel time step run with Arithmetic mean model.

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⁶)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 fine-grid solution, but the Darcian mean model shows superior error and stability characteristics.

Non-Darcian flow errors are not apparent in this example for vertical space step sizes below where the arithmetic mean actually violates the min-max 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 non-Darcian flow errors are tolerable in many cases, as the wide-spread 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 van-Genuchten-style 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 non-physical 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.

Figure 5: Trailing Edge Matric Suction Spike for arithmetic mean model for np = 40, 60, 100, 140, 200, 280 and 400, Dx = 12.8 to 1.28m; The grid point with the largest spike value at t = 8.39(108) s is tracked from t = 10 to 8.39(108) s.

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 no-flow 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 non-physical 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 non-Darcian 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 steady-state 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 steady-state equilibrium, then it cannot return a perturbed system to a physically correct steady-state equilibrium. Further, it may even violate the minimum-maximum principle for elliptic boundary value problems.

This papers uses a simple three-point grid test to show the non-physicality 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 minimum-maximum principle as well. It may be that in many cases the non-Darcian 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 min-max principle in steady-state flow, it produces non-physical oscillations similar to those of numerical dispersion. It also shows the ability to generate non-physical dry and wet spikes due to its underestimation of flow for wet-over-dry conditions and overestimation of flow for dry-over-wet 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 well-behaved at the space step sizes studied. It produces no apparent non-physical 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 non-Darcian means. The apparently nonlinear, astable nature of non-Darcian 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 non-Darcian 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)