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.