Uncertainty Associated with Semivariograms Used for Site Simulation

Received September 1992, revised March 1993, accepted March 1993.
Vol. 31, No. 5 -- GROUND WATER September-October 1993

Abstract

Introduction

At this time, the best approach is to integrate all available data from a site into a range of possible subsurface interpretations and then consider the probability of satisfactory performance of alternative remedial actions. Multiple indicator conditional simulation blends indicator kriging and stochastic simulation to statistically evaluate the range of possible subsurface geologic configurations (Journel, 1983; Journel and Isaaks, 1984; Alabert, 1987; Journel and Alabert, 1988; Johnson and Dreiss, 1989). Knowledge of the range of possible subsurface conditions aids the modeler in defining best, worst, and most likely case scenarios as well as the probability of occurrence of particular scenarios given the available data. To date, such simulations have been carried out at the level where the kriging matrix is solved, without incorporation of the uncertainty associated with definition of the semivariogram. Such an approach is based on the assumption that the specified semivariogram models are absolutely correct, but this is rarely true, particularly when one considers the limited data usually available at a typical hazardous waste site. In such a situation, use of the estimation error to evaluate the accuracy of the kriging is misleading because it appears to characterize uncertainty associated with the result but ignores the uncertainty associated with selection of the semivariogram. The result of a kriging process is based on the definition of the semivariogram. By evaluating the uncertainty in the semivariogram, the greater range of uncertainty associated with the simulated results becomes apparent. Uncertainty in the simulation process can be more completely evaluated by using methods such as Jackknifing (Shafer and Varljen, 1990; Davis, 1987), latin-hypercube sampling (McKay et al., 1979), and expert opinion (McGill, 1984) in defining the semivariogram models to be used for stochastic simulation. These methods are discussed in this paper.

Data collection is time-consuming and expensive. Data collection can be performed more efficiently by examining data as they are collected, preparing experimental semivariograms, plotting estimation errors, and using the results to select subsequent data types and locations. Some projects have used estimation errors to identify areas of' greater uncertainty which can be targeted for further data collection, thus optimizing dollars spent in site characterization (Olea, 1974; Orr and Dutton, 1983). Similarly, evaluation of experimental semivariograms as data are collected can guide the data collection program.

Because data are usually limited, the results of kriging can be misleading; depending on the parameters used to define the semivariogram, the same data can yield different results. Although kriging will produce results that honor the data, the estimated values at locations between sample sites are nonunique. The simple examples shown in Figure 1 demonstrate this point. These hypothetical two-dimensional models represent two distinctly different geologic settings that are indistinguishable by examination of only the well data, which are identical in Figure 1a and 1b. Of the 11 well borings, six are in fine-grained sediments of relatively low hydraulic conductivity (low K) and five are in coarse-grained sediments of generally high hydraulic conductivity (high K). Ideally more data should be collected, but because of cost constraints or constraints on drilling locations, this may be the only data set that can be used. Because different geologic configurations can yield distinctly different contaminant plumes (Figure 2), incorrect modeling of the site, or failure to recognize the uncertainty associated with subsurface interpretation, can result in remedial action that does not accommodate conditions at the site. An example of such a situation is illustrated by Poeter and Gaylord, 1990.

Fig. 1. Borehole data used to interpret the subsurface may not provide a unique solution. In this case, there are 11 data samples: six of fine-grained sediments with low hydraulic conductivity, and five of coarse-grained sediments with high hydraulic conductivity. Although data for each map are identical, the nature of the geology in each map is substantially different. This illustrates that there is uncertainty associated with the interpretation of the character of subsurface at locations that have not been sampled.

Fig. 2. Contaminants will migrate in different patterns within the two geologic models presented in Figure 1. It is important to evaluate the probable alternative scenarios when designing a remediation plan.

It is not sufficient to utilize a code that calculates an experimental semivariogram and selects a suitable model. For good results, the modeler must evaluate the uncertainty of the data. In many, if not most cases, there is not enough data available to clearly and absolutely define the semivariogram, but by incorporating the modeler's knowledge or expert opinion about the site, uncertainty may be reduced, possibilities limited, and reasonable results may be identifiable.

Semivariograms

for a particular lag distance (h), where N = number of data pairs in the search area, and z(x_i) and z(x_i + h) are all the pairs of the N samples within the lag range, h. The search area is defined using a search direction and half angle. The search direction is measured clockwise from north (or the positive z axis for a cross section) and defines parallel lines along which data of the given lag distance must fall in order to be used in the calculation of the semivariogram (Figure 3). Often data exhibit anisotropy; consequently the experimental semivariogram is calculated in a number of directions. The major axis of the anisotropy is indicated by the search direction of the semivariogram with the longest range (range is the separation distance at which the semivariogram value reaches the population variance and is discussed later). Generally, few data will lie directly along a search direction line; therefore a tolerance angle (defined as the search half angle) is used to include data that are offset from the line (Figure 3). It is useful to note that any search direction accompanied by a search half angle of 90' includes all combinations of orientations of points at each spacing, thus is appropriate when evaluating data with an isotropic distribution. Johnson and Dreiss (1989) provide useful illustrations and examples of the features associated with indicator semivariograms of lithologic distribution.

Fig. 3. Features of a semivariogram and parameters defining the search area (after Englund and Sparks, 1988).

The model semivariogram, (h), is a function representing the experimental semivariogram. The distance at which the model semivariogram meets the data set variance is defined as the range (Figure 3). The variance of the sample at a separation distance of zero is called the nugget (Figure 3). This terminology arose in the mining industry where two assays from the same gold sample would sometimes yield markedly different results due to the presence of a gold nugget in one piece of the sample while another piece includes only disseminated gold. The variance of the entire data set is referred to as the sill (Figure 3). Points falling between the 1st and 2nd lag distances are used to calculate (h) for a given h. A maximum bandwidth excludes points that lie in areas well off the search direction from being included in the semivariogram calculation. The simplest spherical model has one sill and is presented by Journel and Huijbregts (1978). A nested spherical model is fit to the experimental semivariogram, in this case with two structures, and is expressed as:

where C₁ = point variance for first nest; C₂ = point variance for second nest; a₁ = range for the first nest; a₂ = range for the second nest; C₀ = nugget; and h = sample separation.

Indicator Kriging and Stochastic Simulation

An experimental semivariogram based on the well data from Figure 1 is presented as Figure 4. For simplicity in demonstrating concepts, only two indicators were employed, one for low hydraulic conductivity materials and another for high hydraulic conductivity materials. Although both models in Figure 1 share the same experimental data, semivariograms generated using many data points selected from the two models (315 points vs 11 points) are substantially different (Figure 5). These semivariograms developed from the extensive data sets illustrate that use of only one experimental semivariogram of the raw data may lead to inaccurate simulations. The actual experimental semivariogram (Figure 4) is based on very few points, and arbitrary, simple assumptions (in this case a search direction of 0° with a 90° half-angle). When more restrictive searches were analyzed (i.e., searches with different directions and smaller half-angles), it was not possible to define anisotropy in the data. That is, there are not enough data to develop convincing statistics to indicate a distinctly longer range is obtained by orienting the semivariogram in a particular direction. A spherical model was defined for the experimental semivariogram of Figure 4, where: spherical model: 0° search direction, 90° half-angle; range = 170 feet; C₁ = 0.273; and C₀ = 0.00. This semivariogram, developed from the 11 data points, contrasts to the semivariograms developed from the extensive data sets in Figure 5. An extensive data set taken from the model presented in Figure 1 a yields a semivariogram (Figure 5a) with the following character: spherical model: 0 search direction, 90° half-angle; range = 190 feet; C₁ = 0.251; and C₀ = 0.00. This model is similar to the semivariogram model determined using the field data and simple assumptions, and though they are not identical, the simulated results would be similar.

Fig. 4. Experimental and modeled semivariograms developed from the 11 labeled data points in Figure 1. A great deal of uncertainty is associated with the modeled semivariogram because of the limited number of data.

Fig. 5. These experimental semivariograms based on 315 data points from the models in Figure 1 were determined by overlaying a regular grid (25' X 25') on each model. The distribution of high and low conductivity materials in Figure la was determined to be isotropic and is described by the model semivariogram in 5a. In Figure 1b, the distribution is anisotropic and the major and minor axes of the model semivariogram ellipsoid are shown in 5b and 5c, respectively.

The purpose of this example is to illustrate that much of the uncertainty in the kriging process is directly accountable to the definition of the modeled semivariogram. The difficulty, however, is that at an actual site, sparse data often result in unsatisfactory experimental semivariograms (Journel and Isaaks, 1984; Fogg, 1986; Alabert, 1987; Ravenne and Beucher, 1988; Johnson and Dreiss, 1989; Fogg, 1989; Rautman and Treadway, 1991). Two techniques, jackknifing and latin-hypercube sampling, can be used to address the uncertainty associated with the semivariograms. In some cases, it may also be reasonable to bias the results with expert opinion. Use of expert opinion in formulating semivariograms may lack statistical rigor, but may be necessary to limit the possibilities. Ground-water hydrologists are hired for their expertise; exercising it, as opposed to blindly following a statistical method which we know has limitations, can improve results. Of course, hydrologists must remember to keep an open mind about the nature of the subsurface at a site and not assume the presence of trends nor assume a simple pattern (such as that of Figure la as opposed to that of Figure lb) without sufficient observation.

Jackknifing

To circumvent this problem, a process called jackknifing is used (Shafer and Varljen, 1990; and Davis, 1987). Jackknifing is a procedure where the experimental semivariogram is calculated with one (or more) data point(s) removed from the data set. By repeating this procedure for every point in the data set, a series of n (n = number of samples) experimental semivariograms is calculated. For each lag distance there are now n values. Using these values, confidence limits can be approximately determined, for the mean at a particular lag. When these are plotted, the error bars define the possible range of the modeled semivariogram. The problem with this method is that each mean value is correlated with the other mean values calculated at that lag (the same data, except for one point, is being used); therefore, the variance calculations are not strictly correct (Davis, 1987). However, this technique is not being used to select the best semivariogram model, which it cannot do (Davis, 1987). Rather it is used to guide the modeler in optimizing further data collection or identifying a likely range of reasonable model semivariograms.

A semivariogram developed by 'jackknifing the II data points from the models in Figure 1 (using a 0° search direction with a 90° half-angle) is presented in Figure 6. By examining the error-bars, it can be seen that the modeled spherical range could vary from less than 70 feet to more than 155 feet, but is probably less than 250 feet (error-bars are set at 95% confidence). This compares favorably with the experimental semivariogram shown in Figure 5a.

Fig. 6. Jackknifing the 11 data points indicated in Figure 1 allows evaluation of uncertainty associated with the semivariogram. The vertical error-bars define the 95% confidence intervals for the mean of each lag. The variance around the mean lag is represented by the horizontal error bars. Each data point represents one instance of a jackknifed experimental semivariogram. This experimental semivariogram is based on the assumption of an isotropic material distribution.

The jackknifed experimental semivariogram (Figure 6) also suggests that 11 data points are not enough to correctly define the model semivariogram. The data are not even sufficient to determine if the drilling pattern is tight enough to be within the range of the local variance, as indicated by the fact that the upper limit of the uncertainty bars associated with the smallest sample separation Is above t total (population) variance (the sill). This suggests that further drilling (data collection) is required. Ideally, a jackknifed semivariogram will appear more like Figure 8; which is calculated from a data set consisting of many digitized points from contour lines of a smooth topographic surface. In contrast, indicator semivariograms calculated from sparse, irregularly spaced samples of discontinuous lithologic units reflect much higher uncertainty (e.g., Figures 6 and 7). Uncertain semivariograms are the norm rather than the exception as indicated by the work of Shafer and Varljen (1990) and the erratic nature of published indicator semivariograms of lithology (Journel and Isaaks, 1984; Fogg, 1986; Alabert, 1987; Ravenne and Beucher, 1988; Johnson and Dreiss, 1989; Fogg, 1989; Rautman and Treadway, 1991). The lack of variation in the experimental jackknifed semivariogram illustrated in Figure 8 allows the modeler to clearly define the model semivariogram. If the experimental jackknifed semivariogram of lithology at a site had the character of Figure 8, it could be argued that too much money was expended collecting data; the semivariogram could have been modeled adequately with fewer data. In this case, if jackknifed semivariograms had been calculated while data were being collected, the characterization program could have been terminated sooner or redirected to focus on collecting data to reduce uncertainty in poorly characterized areas of the site (as indicated by areas of high kriging estimation error), as opposed to collecting data that would further define the semivariogram, thus saving time and money.

Fig. 7. Although anisotropy cannot he identified by evaluating single semivariograms of the 11 data points, anisotropic character is hinted at when the same data are jackknifed along specified search directions. For a search direction of 45°, the range is likely to be less than 100 feet. In the 135° search direction, the range is likely to be greater than 150 feet, and possibly more than 500 feet. The anisotropy defined in Figure 5b-c cannot be determined from the 11 data points, but its possibility is indicated by the data. Symbols are described in the caption of Figure 6.

Fig. 8. When a substantial amount of data are collected, the experimental semivariogram may be clearly defined. In this jackknifed simulation, there is tittle uncertainty in the lag means, and there would be little uncertainty in defining the model semivariogram (this semivariogram was not calculated from the data presented in this paper).

Additional Comments About Jackknifing

Latin-Hypercube Sampling

An alternative approach to random selection of a large number of possible semivariogram models is to use latin-hypercube sampling. This reduces the number of simulations required to insure the "flavor" of all alternatives is addressed (McKay et al., 1979). For this example, one might suggest the nugget must fall within one of four equiprobable regions, and the range also must fall within one of four equiprobable regions. The actual nugget, or range within each region is randomly calculated (Figure 9). This allows 16 model semivariograms to be calculated for an isotropic model. For an anisotropic model, the direction and magnitude of the anisotropy can be restricted similarly. This, however, produces many more simulations. If the anisotropy factor between the major and minor axis is evaluated at four ratios (e.g., 1.0, 0.5, 0.25, and 0.125 or some other ratios as determined from jackknifing the data to obtain a semivariogram in the direction of the minor axis of anisotropy), the number of simulations is increased to 64. If the search directions, 0° to 180°, are divided into four directions (0°, 45°, 90°, and 135°), the number of simulations is increased to 256.

Fig. 9. Reasonable models must be selected from the shaded region in 9a to represent the "flavor" of the alternative interpretations of the data. Four model semivariograms with a nugget selected from the lower quartile of possible nugget values are shown in 9b. The ranges of the four semivariograms are selected to represent each of the quartiles of possible ranges. Sixteen models would be used to represent the distribution of semivariogram models for the isotropic case. Symbols are described in the caption of Figure 6.

Recall that the objective of this approach is not to make a single best estimate of the subsurface interpretation, but to evaluate the possible range of subsurface character based on available data. From a purely mathematical approach, this may be computationally intractable; however, incorporation of expert opinion into the process makes it possible to limit the reasonable alternatives.

Expert Opinion

The process of stochastic simulation uses probabilities to estimate a value at a grid location. Unfortunately, these probabilities are based on measured values near that location and, consequently, geologically impossible configurations can be simulated. For example, the "law of original horizontality" and the "principal of stratigraphic superposition" are readily broken. Eventually, techniques that incorporate these concepts into stochastic simulation will be developed. Until that time, such simulations must be identified, deemed unreasonable, and discarded.

Although creation of such geologic fallacies cannot be prevented with the current simulation process, the simulations can be improved by incorporation of geologic knowledge at analog sites into the process. An expert can infer more information about the site than is evident in the borehole data. For example, an expert may know that sand lenses in the area tend to be between 10 and 25 feet thick. The borehole data at the site may be too sparse to determine this range of thickness, but knowledge from analog sites in the area may render it reasonable to assign a range of 10 to 25 feet to the vertical modeled semivariogram. Although such action is not based on data from the site, knowledge of analogs adds information to (decreases uncertainty associated with) the simulation process. If the site is made of horizontally bedded alluvial deposits, there is no reason to run simulations which assume the material distribution is isotropic. In such settings, units are generally continuous for greater distances horizontally than vertically. The modeler may be able to confirm the presence of layered anisotropy by demonstrating that semivariograms with different search directions and limited half-angle and bandwidth have the potential to have different ranges. Even if the indications are sketchy due to sparsity of data, the modeler can limit the simulations to produce only reasonable interpretations given the local geology. Similarly, anisotropy may be present in lateral directions, and geologic knowledge of directional trends of lenses or channels may be used to limit the number of orientations considered for semivariograms which will, in turn, limit the number of simulations that must be undertaken.

There is little reason to evaluate solutions that are mathematically possible, but geologically improbable. Discarding geologically improbable solutions adds "bias" to the results that may have to be defended later. However, omission of the bias means that we do not use all of the information available to us. When expert opinion is used wisely, the bias is likely appropriate, and will speed the site evaluation, thus limiting exploration and analysis costs.

Results

These simulations were created using the multiple indicator conditional simulation code ISIM3D (Gómez-Hernández and Srivastava, 1990). The map area was modeled in two dimensions using a 50 X 35 grid, with 10 foot square grid cells.

Simulations resulting from use of the modeled semivariogram using the extensive data set (Figure 5a) are presented in Figure 10. These simulations differ significantly from the model because only the 11 data points were used to condition the simulation (if the 315 data points used to develop the semivariogram were also used as conditioning data, the simulations would be much more similar to the model). Simulations presented in Figure 11 are based on a model semivariogram (a₁ = 115', C₁ = 0.25, C₀ = 0.0) sampled from the jackknifed experimental semivariogram shown in Figure 6. Although neither simulation (Figure 11a or 11b) is identical to the model in Figure la, they are reasonable approximations considering the limited data. The simulation in Figure 11b is particularly close to the model of Figure 1a. The appearance of the resulting simulation is rather insensitive to the choice of range (compare Figure 11 with Figure 10 which was generated using a range of 190'). Both the experimental semivariograms (Figure 5a and Figure 6) were developed based on an assumption of isotropic material distribution. The simulations in Figure 10a and Figure 11a are nearly identical because the same random path (same random number seed) was used to generate all of the(a)simulations in Figures 10-13. A different path was used to generate the (b) simulations. These simulations bear little resemblance to the model in Figure 1b which is a viable interpretation of the data from the 11 field measurements. This inability to represent the full range of possible interpretations is not unexpected.

If expert opinion indicated that the site would be expected to exhibit the locally observed NW-SE trend of high and low hydraulic conductivity deposits, then the simulations presented in Figures 10 and 11 could be assumed to be less probable. They would be superseded by the probability f occurrence of anisotropic representations of the site. If such expert opinion were not available, the two alternative configurations would have to be considered equally likely to occur. The two simulations in Figure 12 were generated in the latin-hypercube sampling process, using one of the semivariograms that would fall in the shaded area in Figure 9a with a range between 120 and 180 feet (third quartile estimate of range), a nugget between 0.0 and 0.061 (first quartile estimate of the nugget), an anisotropy factor of (minor to major axis) 0.125, a major axis orientation of 135°, and using different random paths through the grid. Although they are not identical to the model in Figure 1b, they mimic its nature. When using the range, sill and nugget terms identified by the semivariogram developed from the extensive data set (Figures 5b and 5c), the simulation results (Figure 13) are not significantly different from the simulation results (Figure 12) obtained using the jackknifed semivariogram (Figure 7), indicating that an extensive field sampling would not improve the character of these simulations but might improve the certainty of occurrence of units with a 135° orientation. That is, more data will improve the certainty of the semivariogram having a given orientation whereas the jackknife approach only indicates the possibility of units having that orientation. Of course, a larger data set improves conditioning of the simulations.

Fig. 10. These two simulations were generated assuming isotropy and using the model semivariogram developed from the extensive data set and illustrated in Figure 5a. The solutions are a reasonable approximation of the map in Figure la.

Fig. 11. These two simulations were generated assuming isotropy and using a latin-hypercube sample from the jackknifed model semivariogram (C₀ = 0.0, C₁ = 0.25, a₁ = 115') developed from the H data points and illustrated in Figure 6. The solutions are a reasonable approximation of the map in Figure la, and are very similar to those generated in Figure 10. Much of the reason that the simulations in Figures 10 and 11 are similar is that the same random path through the grid was used to simulate 10a and 11a, and another path was used to simulate 10b and 11b.

Fig. 12. These two stochastic simulations were generated assuming anisotropy using the jackknifed model semivariogram based on the 11 data points and illustrated in Figure 6. The latin-hypercube technique was applied, and these are two simulations of a potential 256, as described in the text. Even though the geologic models presented in Figure 1 are different, use of jackknifing and latin hypercube sampling can produce both configurations from limited data. These solutions are a reasonable approximation of the map in Figure 1b. Unfortunately, the method will not indicate whether these simulations or the simulations in Figures 10 and 11 are the most likely because the data are not sufficient to draw such a conclusion.

Fig. 13. These two simulations were generated assuming anisotropy using the extensive model semivariogram based on the extensive data set and illustrated in Figures 5b-c. The solutions are a reasonable approximation of the map in Figure 1b, and are very similar to those generated in Figure 12, indicating that extensive data are more important to determining the character of the semivariogram than they are to conditioning the simulation.

Conclusions

By jackknifing the data to determine a reasonable range of model semivariograms, and using latin-hypercube sampling and incorporating expert opinion to limit the required simulations, the uncertainty associated with the subsurface interpretation can be more completely assessed utilizing a reasonable amount of simulations. Unfortunately, the uncertainty may be so great that little can be concluded about the site. However, this is important information because it indicates that more data must be collected before conclusions are made about the site. Given one data sample, one can begin to make interpretations of the site, but quantifying the uncertainty associated with those interpretations is important.

The method presented herein is useful when a significant amount of uncertainty is associated with the experimental semivariogram. If the uncertainty is small, the process only adds unnecessary work. For small data sets, where there is significant uncertainty, this process may be the only way to correctly assess the potential variability of the subsurface, and evaluate potential flow paths for contaminants.

Although the application considered herein Pertains to indicator conditional simulation, evaluation of the uncertainty associated with a semivariogram is important whenever a semivariogram is used. Many commonly used contouring codes are based on kriging and have an algorithm that automatically defines the semivariogram. This semivariogram may assume the same semivariogram model for all data sets, and may use rather arbitrary procedures for determining the semivariogram parameters (i.e. the nugget, sill, and range). The lack of uncertainty analysis in defining the semivariograms used in such contouring codes not only leaves the user without alternative interpretations, but also may provide one of the less likely interpretations.

Acknowledgments

References

Davis, B. M. 1987. Uses and abuses of cross-validation in geostatistics. Mathematical Geology. v. 19, no. 3, pp. 241-248.

Englund, E. and A. Sparks. 1988. GEO-EAS. U.S. Environmental Protection Agency, Environmental Monitoring Systems Laboratory. EPA600/4-88/033.

Fogg, G. E. 1986. Stochastic analysis of aquifer interconnectedness, with a test case in the Wilcox Group, East Texas. Ph.D. dissertation, Univ. of Texas at Austin. University Microfilms International, Ann Arbor, MI. p. 216.

Fogg, G. E. 1989. Stochastic analysis of aquifer interconnectedness: Wilcox Group, East Texas. Bureau of Economic Geology, Austin, TX. Report of Investigations No. 189. Gómez-Hernández, J. J. and R. M. Srivastava. 1990. ISIM3D: An ANSI-C three dimensional multiple indicator conditional simulation program. Computers in Geoscience. v. 16, no. 4, pp. 395@.

Johnson, N. M. and S. J. Dreiss. 1989. Hydrostatic interpretation using indicator geostatistics. Water Resources Research. v. 25, no. 12, pp. 2501-2510.

Journel, A. G. and Ch. J. Huijbregts. 1978. Mining Geostatistics. Academic Press, London. p. 207.

Journel, A. G. 1983. Nonparametric estimation of spatial distributions. Mathematical Geology. v. 15, no. 3, pp. 445-468.

Journel, A. G. and E. H. Isaaks. 1984. Conditional indicator simulation: Application to a Saskatchewan uranium deposit. Mathematical Geology. v. 16, no. 7, pp. 685-718.

Journel, A. G. and F. G. Alabert. 1988. Focusing on spatial connectivity of extreme-valued attributes: Stochastic indicator models of reservoir heterogeneities. Society of Petroleum Engineers 63rd Annual Technical Conference SPE 18324, Richardson, TX. pp. 621-632.

McKay, M. D., R. J. Beckman, and W. J. Conover. 1979. A comparison of three methods for searching of input variables in the analysis of output from a computer code. Technometrics. v. 21, no. 2, pp. 239-245.

McGill, R. E. 1984. An Introduction to Risk Analysis. 2nd edition. Penn Well Publishing Co., Tulsa, OK. 274 pp.

Olea, R. A. 1974. Optimal contour mapping using universal kriging. Journal of Geophysical Research. v. 79, no. 5, pp. 695-702.

Orr, E. D. and A. R. Dutton. 1983. An application of geostatistics to determine regional ground-water flow in the San Andreas Formation, Texas and New Mexico. Ground Water. v. 21, no. 5, pp. 619-624.

Poeter, E. P. and D. R. Gaylord. 1990. Influence of aquifer heterogeneity on contaminant transport at the Hanford site. Ground Water. v. 28, no. 6, pp. 900-909.

Rautman, C. A. and A. H. Treadway. 199 1. Geologic uncertainty in a regulatory environment: An example from the potential Yucca Mountain nuclear waste repository site. Environmental Geology and Water Sciences. v. 18, no. 3, pp. 171-184.

Ravenne, C. and H. Beucher. 1988. Recent developments in description of sedimentary bodies in a fluvio deltaic reservoir and their 3D conditional simulations. Society of Petroleum Engineers 183 1 0, Richardson, TX. pp. 463-476.

Shafer, J. M. and M. D. Varljen. 1990. Approximation of confidence limits on sample semivariograms from single realizations of spatially correlated random fields. Water Resources Research. v. 26, no. 8, pp. 1787-1802.

Wingle, W. L. and E. P. Poeter. 1992. Evaluation of uncertainty associated with contaminant migration in ground water A technically feasible approach. Proceedings of the Fifth NGWA Conference on Solving Ground Water Problems with Models. pp. 687-695.