Are the current models of unsaturated flow and transport adequate to predict the performance of the Yucca Mountain nuclear waste site, so as to ensure its eternal safety? Perhaps. Perhaps not. This work introduces a new way to look at some common modeling concepts that will raise valid questions.

A May 1997 report to the U.S. Dept. of Energy Yucca Mountain Site
Characterization Office, "Unsaturated zone flow model expert
elicitation project" (Coopersmith, et al., 1997) confirms
the importance of modeling. It states, "It is important to
the performance assessment to understand how much of that volume
[percolation flux expressed as the volume rate of flow per unit
area] is carried in the rock matrix and how much in the fracture
of TSw [Topopah Spring welded tuff]." It cites testimony
of the elicited experts, placing substantial reliance in existing
models to make decisions. For example, "Dr. Stephens noted
that this partitioning is quite uncertain and, as shown by the
UZ site-scale modeling, the relative components of the flux are
a function of the total percolation flux.", and, "Dr.
Campbell postulated that although the fractures probably carry
90% of the total flux, less than 1% of the total flux occurs in
fast flow, based on simulations conducted by LBNL and LANL."

This investigator is not privy to the inner workings of the models
cited in Coopersmith, et al., but nevertheless will raise valid
questions for those who are better informed of them to answer.
Warrick (1991) first pointed out that finite method models of
unsaturated flow often neglect to calculate interblock mean flows
that are entirely consistent with Darcy's law. In models where
the pressure head or matric suction, and thus the unsaturated
conductivity through the conductivity relation, are known only
at grid points, it is common to assume that flow during a time
step is constant with space between adjacent grid points. This
requires that the mean hydraulic conductivity between the grid
points be estimated and used to calculate flow by Darcy's law.

Warrick pointed out that if Darcy's law is integrated between
grid points, under the assumption of constant flow, it requires
flow that is often different by orders of magnitude from that
calculated with standard interblock conductivity means, such as
the arithmetic, geometric and upstream. Unfortunately, some soil
classes produced integration through a singularity with Warrick's
approach. Baker (1994, 1995) recast the problem as an elliptic
boundary condition solution, and demonstrated that Darcian means
for vertical unsaturated flow vary nonlinearly both with model
and medium parameters. The variation shows a complexity hitherto
unrecognized in its richness.

In 1997, Baker won a Phase I Small Business Innovation Research
contract to continue this work under contract DE-FG02-97ER82329.
This effort produced a new modeling and mathematical framework
for other investigators to use to judge the appropriateness of
interblock hydraulic conductivity means and their approximations
for unsaturated flow, particularly in the vertical. It includes
most standard interblock unsaturated conductivity means as special
cases. Topopah Spring welded tuff was used as the example medium.
This work will raise valid questions about any model of vertical
unsaturated flow that has been designed without taking its results
into account.

This will not be welcome in some quarters. Those who have staked
their reputations on the assurances of mass balance as the final
arbiter of model accuracy may not gladly admit to having missed
something. Some have already been strongly critical and dismissive,
even to the point of verbal abuse and one charge of "flim-flam".
Fortunately, the math is there for anyone to use and explore,
independent of any agenda. Unfortunately, mass balance only guarantees
that the calculated flow agrees with the assumptions, not that
the assumptions are correct. If the assumed interblock means are
off by orders of magnitude, then so are the flows, and so are
the advective components of contaminant transport. If this is
an incorrect assessment, it must be demonstrated so by the mathematics,
not by bluster.

If an interblock mean produces non-Darcian flow, it is not enough
to demonstrate that a numerical method is consistent with the
partial differential equation to guarantee a physically-based
result. Consistence only means that the numerical method approaches
the modeled partial differential equation in the limit as the
space and time steps approach zero extent. For a given distribution
of moisture and conductivity, as step size goes to zero, the adjacent
conductivities approach each other. No interblock mean can be
physically valid unless it achieves the common mean in the limit
as adjacent conductivities approach each other. In effect, all
means approach the each other in that limit. But consistence is
only a necessary condition. Along with stability, it only guarantees
convergence to a solution in linear problems, not convergence
to a correct solution, especially for nonlinear problems.

The better trick maintains consistence with the assumptions as
step sizes increase or the conductivities diverge. Warrick (1991)
demonstrated that this is often not the case with standard means.
The resulting errors have long been neglected. Not one of the
experts in the Unsaturated Zone Flow Model Expert Elicitation
Program (Coopersmith, et al., 1997) addressed them, even though
they can amount to orders of magnitude. This investigator is ignorant
of how the models used in that study may have adequately addressed
them, but the work of understanding and producing Darcian means
is difficult and recent. In this investigator's admittedly limited
cognizance, this work is the first exposition of workable, practical,
verifiable approximations of Darcian means for running models
of vertical, unsaturated flow.

There are, of course, qualifications. The scope of this work is
limited, and much remains to be done in developing the concepts.
This work only addresses the development of Darcian means for
vertical and horizontal single-phase unsaturated flow in homogeneous
media. It does not address such things as ponded heads, layered
media, multi-phase flow, hysteresis or capillary barrier effects.
It has thus far produced analytic equations for Darcian mean only
for the case of an exponential conductivity relation, and workable
approximations only for the case of a Brooks-Cory conductivity
relation.

In application to Topopah Spring welded tuff, it assumes a composite
conductivity relation, after the fashion of Pruess and Wang (1987),
where fractures are assumed to be vertical and parallel to the
flow. True Darcian means and flows (relative to saturated flow)
were developed for this case using the elliptic boundary condition
solution approach, over five decade changes in space step size
from 1 cm to l km, and a relative conductivity range of about
10-10 to 1. The relative flows for four standard means, the arithmetic,
geometric, upstream and an integral form are developed for the
same cases to demonstrate the magnitude of the non-Darcian flow
errors.

Even in this limited case, there are some interesting nonlinear
surprises. Both standard and Darcian means can predict flow between
two grid points that can jump up in magnitude as one of the conductivities
decreases. For this kind of composite conductivity relation Darcian
means predict local regimes where intergrid flow stays significantly
constant with changes in adjacent conductivities, where standard
means do not. And these effects can change dramatically with model
space step size.

These results, by themselves, have significant implications for
any model of unsaturated flow in Topopah Spring tuff. Several
investigators in the AGU Geophysical Monograph 42, "Flow
and Transport Through Unsaturated Fractured Rock" (1987),
Reda, Eaton and Bixler, and Peters, et al., studied relatively
small core samples of Yucca Mountain rock, presumably with the
intent of transferring the results to larger scales. The work
presented here demonstrates that transferring the results of models
with sub-millimeter space steps to field-scale steps of tens to
hundreds of meters, as noted in Coopersmith, et al., (1997), may
be neither mathematically nor physically valid, unless the appropriate
changes in Darcian means are taken into account.

Further, some of the experts in Coopersmith, et al. (1997), make
some recommendation of inverse modeling to estimate field-scale
parameters, while at least one notes with dismay that few actual
measurements of unsaturated hydraulic parameters are available
for rocks involved. The common pitfall in inverse modeling and
parameter estimation is that, in the presence of large noise and
uncertainty, there are many possible model parameters that will
fit a particular set of time-series data. But if the accepted
fit does not truly correspond to the actual physical system, the
model will break down when the boundary conditions and forcing
functions change. This investigator has had the experience of
using limited, noisy data to fit a model of infiltration to a
dry year, only to have it fail in the following wet year. This
can only be exacerbated if the model is flawed with non-Darcian
flows.

The reader should understand that the composite conductivity relation
chosen here to represent Topopah Spring welded tuff is only meant
to be a vehicle for understanding and a connection with previous
work, not an absolute truth. The work from which its parameters
were derived, Schenker, et al. (1995), notes that the expected
value of fracture spacing in this tuff is 0.740m, with a coefficient
of variation of 3.147, a minimum of 3*10-6m and a maximum of 19.878m.
This calls into question the validity in continuum mechanics of
any use of this relation in a representative elementary volume
with any dimension smaller than several meters to tens of meters.

This investigator believes that one should use linked models of
the matrix and fracture flow with the elliptic boundary condition
solution approach to properly calculate Darcian means for this
fractured medium. Then the issue of the proper representative
elementary volume will be better addressed, as will issues such
as layered media and capillary barriers. The reader is cautioned
to understand that such work will likely produce more surprises.
This investigator decided to develop the means make adequate approximations
to Darcian means on the grid scale of the elliptic boundary condition
approach before tackling those problems. This will allow the calculation
of Darcian means to bootstrap from previous approximations that
will help the convergence and stability of the solutions.

Those approaches are presented here for the cases of an exponential
conductivity relation and a Brooks-Corey relation. The van Genuchten-style
relation has been left for another time. Then the analytic solution
for the exponential relation and the approximation for the Brooks-Corey
relation are applied piecewise to the cases developed for the
composite relation to see if they improve the flow error over
standard means.

More than one reviewer, including one of the Phase I SBIR proposal,
has noted the difficulty of this material. This is a new approach
in an area long thought to be fully mined of ideas. Although a
large number of graphs have been included to help with understanding,
the reader who is not fully prepared to set aside preconceptions
and delve deeply into the math may find this presentation unrewarding.

Although the concept of Darcian means will help to show how fracture
and matrix flow must be handled differently in the same model,
particularly as model scales change, it is only a small part of
the puzzle of unsaturated flow in fractured media. As small, perhaps,
as an o-ring on a space shuttle solid rocket booster. Only time
and effort and math will tell.

REFERENCES:

- 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, 99p.

Baker, D.L., 1995, Darcian weighted interblock conductivity means for vertical unsaturated flow, GROUND WATER, May-Jun 1995, 33(3):385-390.

Coopersmith, K.J., R.C. Perman, R.R. Young, M. Pendleton and J.L. Younker, 1997, Unsaturated zone flow model expert elicitation project, document WBS 1.2.5.7 / SCPB: NA / QA: L of the Civilian Radioactive Waste Management System and Operating Contractor, Prepared for the U.S. Dept. of Energy, Yucca Mountain Site Characterization Office, P.O. Box 30307, North Las Vegas NV 89193-8608 under Contract DE-AC01-91RW00134, 171p.

*Flow and Transport Through Unsaturated Fractured Rock*, 1987, D.D. Evans and T.J. Nicholson, Des, Geophysical Monograph 42, American Geophysical Union, Washington DC, 187p.

Pruess, K. and J.S.Y. Wang, 1987, Numerical modeling of isothermal and nonisothermal flow in unsaturated fractured rock - a review, Flow and Transport Through Unsaturated Fractured Rock, D.D. Evans and T.J. Nicholson, Des, Geophysical Monograph 42: 11-21, American Geophysical Union, Washington DC.

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 total-system performance assessments, Sandia Report SAND94-0244 * UC-814, Prepared by Sandia National Laboratories, Albuquerque NM 87185 and Livermore CA 84550 for the United States Department of Energy under Contract DE-AC04-94AL85000, 87p + appendices.

Warrick, A.W., 1991, Numerical approximation of Darcian flow through unsaturated soil, Water Resour. Res., 27(6):1215-1222.