Note: This text may not appear properly unless the user has the TrueType font "SYMBOL" installed on his computer.
What is, as one massbalance partisan so thoughtfully put it, "this absurd parameter called the Darcian mean"? The statement is deceptively simple. It is the interblock conductivity mean for unsaturated flow that is consistent with Darcy’s law [1] and the modeling assumption of constant flow between grid centers [2] [Figure 1].

This is a common assumption in finite method modeling of conservative flow in continuum mechanics. In effect, we hold flow constant over the space between two grid points and the time between to time points, and then add the mass changes back in at the end. Mass balance only means that the mass changes agree with the flow, not that the flow is correct.
If the conductivities in adjacent cells are significantly different, and known only at the grid points in the cell centers, that still leaves the question of what average, or mean, conductivity to use in Darcy’s law to calculate flow between the cells. Practitioners of finite element methods sometimes make the mistake in presuming that finite difference techniques can only assume that the properties in the cells are constant throughout the cells. This assumption leads naturally to the harmonic mean, at least for horizontal flow, which tends to predict much lower flows than justified by physical reality.
This need not be the case. Warrick (1991) solved for Darcian means by integrating [1] to get [3], which has to be solved by iterating rq over the possible range until the equation balances. Then the substitution back into [4], the difference form of [1], produces the vertical Darcian mean, Kv. Warrick found that compared to these flows, the flows produced by using standard means, such as the arithmetic, geometric and upstream, can be off by orders of magnitude.
[1] , 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).
[2]
[3] , where
Dx = x_{2}x_{1}, the vertical separation between grid cells centered at x_{2} > x_{1} (m), with the vertical axis positive upwards, and Y_{1} and Y_{2} are the matric suctions at x_{1} and x_{2}, respectively.
[4] , where
Kv = the Darcian relative conductivity mean for vertical flow (dimensionless).
Unfortunately, Warrick’s (1991) approach presents a problem when the integrand of [3] contains a singularity. It is more useful to integrate [1] to get [5], then to solve [4] and [5] for Kv [6]. The components of Kv, are the mean Darcian conductivity for gravity or advective flow [7], Kx, and the mean Darcian conductivity for horizontal and diffusive flow [8], Kh. Warrick’s attention to tilted flow is neglected here. When flow is purely horizontal, then the gravity component disappears and Kv becomes Kh.
When finite differences are used to solve Richards’ equation for unsaturated vertical flow [9], the righthandside can be expressed using the Kx and Kh components [10]. Using true Darcian means reduces the error associated with standard approximations, like the arithmetic mean, and produces solutions more consistent with Richards’ equation. In many cases, Kh is analytically integrable for the chosen form of kr(Y), but there are considerable difficulties in knowing kr(x) in general.
[5]
[6]
[7]
[8]
[9]
[10] ,
where Dx = (x_{i+1}x_{i1})/2
Consider that all approximations of Kx and Kh in [10] can be mathematically reduced to the equivalent form of Kv in [4] and [6]. If the proper combination is not chosen, the effective Kv can be negative, producing nonphysical backwards flow. It is not enough to claim, for example, that only finite element methods can solve this problem by varying physical properties between grid points. The variation has to be consistent with both Darcy’s law and the modeling assumptions to produce Darcian flow.
One way to determine kr(x) between two vertical grid cells, centered at x1 and x2, is to combine [1] and [2] to form the elliptic boundary condition problem [11]. All one need know is the boundary conditions, and the relationship between kr and Y . It can be solved using finite differences and Newton or Picard iteration, until the extreme spread of rq across all the solution fine grid elements falls below a given tolerance. The distribution kr(x) can then be numerically integrated to produce Kx. If kr(Y ) is analytically integrable, then Kh develops from know quantities. Kv is of course determined after solving [11] for the mean rq by back substitution into [4].
[11] ,
where kr(x1) = kr1, kr(x2) = kr2, Y (x1) = Y1, and Y (x2) = Y2, given the relationship kr(Y ).
As noted before, this problem is deceptively simple. This investigator has only been able to find one explicit analytic solution for [11], for the case of an exponential relative conductivity relation, kr(y ) = exp(a (yyd)). All of the rest seem to be, at best, implicit. Using Newton’s iteration method to solve for [11] contains many pitfalls. Either bad initial conditions or bad grid point distributions can cause failure, as well as several other convergence problems, none of which will be discussed in detail here. The problems get worse for the kind of composite function presented in the next section to represent kr(y) for Topopah Spring welded tuff. They can include very sharp fronts, sensitivity to computer arithmetic processing unit bit resolution, small differences in large numbers between the advective and diffusive components, and sensitivity to space steps being, of all things, too small.
One of the toughest problems is getting initial conditions for a Newton’s method solver. This investigator has found it useful to recast the problem as an ordinary differential equation shooting method problem, which comes very close to Warrick’s (1991) original method. If one sets the continuum and discrete forms of Darcy’s law, [1] and [4], equal, one can write the distribution of y(x) entirely in terms of x and y [12].
[12]
Kv^ is the shooting parameter, the estimate of Kv to satisfy the boundary conditions of [11]. It seems to work best to shoot upstream, against the flow, iterating Kv^ until both boundary conditions are satisfied. Using a method like adaptive RungeKutta integration, this provides a rough solution to [11], which can be improved with Newton’s method. With this information, other investigators can verify the results developed here.
The reader should recall and realize a few things. First, this is a steadystate problem, solved between two grid points in a model of Richards’ equation, for the duration of one time step. Because of the computational effort involved, this approach is not practical to use in a running model. We do it here to learn about the nature of Darcian means, and to develop adequate and computationallyefficient approximations for them.
Second, this is three or more approximations removed from physical reality. Here we use a numerical method to approximate the solution to a partial differential equation that describes an ideal soils that approximates a real soil. We do this because we have no better way.
Third, this is a small but critical part of the overall problem. It will help us to understand how to treat matrix and fracture flow differently in a model of unsaturated flow, but it will not resolve all difficulties. However, because of the potential for orderofmagnitude errors that it represents, it cannot lightly dismissed until it has been investigated.
REFERENCE:
Warrick, A.W., 1991, Numerical approximation of Darcian flow through unsaturated soil, Water Resour. Res., 27(6):12151222.