REDUCING UNCERTAINTY ASSOCIATED WITH GROUND-WATER FLOW AND TRANSPORT PREDICTIONS

Eileen P. Poeter
Department of Geology and Geological Engineering
Colorado School of Mines
1500 Illinois St.
Golden, CO 80401 USA
phone: (303)273-3829
fax: (303) 273-3859
email: epoeter@mines.edu

and

Sean A. McKenna
Department of Geology and Geological Engineering
Colorado School of Mines
Golden, Colorado 80401
Current address: Sandia National Laboratories
Geohydrology Department, P.O. Box 5800, M.S. 1324, Albuquerque, New Mexico 87185-1324.)
email: smcken@nwer.sandia.gov

Reprinted from GROUND WATER, Vol. 33, No. 6, November-December 1995

ABSTRACT

Effective evaluation of ground-water flow and transport problems requires consideration of the range of possible interpretations of the subsurface given the available, disparate types of data. Geostatistical simulation (using a modified version of ISIM3D) of hydrofacies units produces many realizations that honor the available geologic data and represent the range of subsurface interpretations of unit geometry. Hydraulic observations are utilized to accept or reject the geometric configurations of hydrofacies units and to estimate ground-water flow parameters for the units (using MODFLOWP). These realizations are employed to evaluate the uncertainty of the resulting value of the response function (ground-water flow velocity and contaminant concentration) using MT3D. The process is illustrated with a synthetic data set for which the "truth" is known, and produces a striking reduction in the distribution of predicted contaminant concentrations. The same system is evaluated three times: first with only hard data, then with both hard and soft data, and finally with only the realizations that honor the hydraulic data (i.e., those accepted after parameter estimation via inverse flow modeling). Using only hard data, the mean concentration predicted for all realizations at the point of interest is nearly two orders of magnitude lower than the true value and the standard deviation of the log of concentration is two. The addition of soft data brings the mean concentration within one order of magnitude of the true value and reduces the standard deviation of the log of concentration to one. After eliminating realizations using inverse flow modeling, the mean concentration is one-third of the true value and the standard deviation of the log of concentration less than 0.5.

INTRODUCTION

A fundamental problem with interpreting the earth's subsurface is that we only sample a small fraction (typically less than one millionth) of the material we are characterizing. The challenge is to determine the character of the material between boreholes and the continuity of high (or low) hydraulic conductivity units. For example, consider a site where there are three boreholes that intersect two hydrofacies (black and white) as shown in Figure 1. If no other information is available, all six interpretations shown (and many more that are not shown) are reasonable interpretations of the subsurface. Knowing that heterogeneity is critical to the movement of contaminants (Poeter and Gaylord, 1990), each interpretation will affect ground-water flow and contaminant transport prediction in different ways. Decisions regarding remedial actions at such a site will differ depending on the actual configuration of units in the subsurface.

Fig. 1. Alternative subsurface interpretation. a: field knowledge of the occurrence of hydrofacies in borings; b through g: a few alternative interpretations of interwell connections.

Usually we know more about a site than just the location of hydrofacies in boreholes, and this information can be used to rule out some, but not all, of the alternative interpretations. "Fusion" of these data (Olhoeft, 1992) reduces uncertainty associated with the interpretation. If each circle in Figure 2 represents the range of possible interpretations, then using all of the information together can significantly reduce the range of possible interpretations, thus reducing the uncertainty of the nature of ground-water flow and contaminant transport. In short, the sum of the parts is less than the whole! If the data circles do not overlap, the project team should strive to identify the shortcomings/errors in the data or determine what assumptions have been made that falsely constrain the interpretation.

Fig. 2. Integrated interpretation constrains predictions.

It is important to work with a range of subsurface interpretations because consideration of these alternative interpretations of the subsurface will yield a range of the response function (ground-water flow and advection-dispersion equation) values (head, flow rate, and concentration). For example, assume there is a contaminated site with two alternative remediation schemes. If there is only one deterministic picture of the subsurface, and the predicted concentration at the point and time of interest is 5 x 10^-3 ppm for remedial alternative 1 and 5 x 10^-4 ppm for remedial alternative 2, than alternative 2 appears to be the better choice. However, when the range of possible subsurface conditions is considered, it may be found that remedial alternative 2 actually has a higher probability of poor performance and the opposite selection would be made (Figure 3).

Fig. 3. Evaluation of alternative remedial actions using a range of interpretations (indicated by the frequency distribution) as opposed to one deterministic interpretation (indicated by the arrow) of the geology may result in selection of a different remedial action.

A current approach to this problem is geostatistical simulation using hard data, soft data, and spatial statistics to describe the character of the subsurface (McKenna and Poeter, 1994). Hard data are data with negligible uncertainty, such as direct measurements of hydraulic conductivity or observations of lithology. Soft data are information with non-negligible uncertainty, such as indirect measurements gathered in a geophysical survey and expert opinion regarding geologic fabric or structure. Hard and soft data are used to develop indicator semivariograms for the hydrofacies (represented by integer indicators) and to simulate numerous realizations [using multiple indicator conditional stochastic simulations, MICSS (Gomez-Hernandez and Srivastava, 1990)] of the hydrofacies configuration which honor the data and the statistical character of the subsurface. The hydrofacies within each realization are populated with hydraulic parameters and a forward flow and transport model is run to yield the distribution of predicted contaminant concentration over the model domain for a period of time. This process is extended through the use of inverse flow modeling which can be used to: (1) improve geologic interpretation, somewhat narrowing the distribution of predicted concentrations; and (2) eliminate realizations that do not fit the hydrologic data, significantly reducing uncertainty.

HYPOTHETICAL PROBLEM

To facilitate demonstration of the process, a scenario requiring a decision concerning establishment of a water-supply well near a stream (Figure 4) is constructed. Alternatively, but at much greater cost, water could be supplied by pipeline. Use of water from the nearby stream is not an option due to poor water quality, which improves by infiltration through the aquifer and dilution with ground water. Knowledge of the probability distribution of concentration at the proposed well location resulting from contaminants entering the groundwater from a nearby landfill will aid in the decision.

Fig. 4. Synthetic model used to illustrate uncertainty reduction procedure.

This study explores the uncertainty associated with predicted contaminant concentration calculated using three approaches and compares those uncertainties with the known concentration. Uncertainty of the geometric configuration and hydraulic conductivity of five hydrofacies are considered. The approaches include: (1) using hard geologic data from the domestic wells in MICSS to produce geologic realizations which are populated with hydraulic conductivities based on hydrofacies type and used in forward flow and transport modeling to predict the probability distribution of contaminant concentration at the proposed well location; (2) repetition of the first approach with addition of soft information regarding continuity of alluvial fan hydrofacies gained through inverse ground-water flow modeling; and (3) repetition of the second approach incorporating hydraulic data through application of inverse flow modeling. The MICSS code used is ISIM3D (Gomez-Hernandez and Srivastava, 1990) as modified by McKenna (1994).

SYNTHETIC SETTING

A synthetic data set is used to demonstrate the process because, with a synthetic data set, the true conditions are known. Reduction of uncertainty is only useful if the distribution of predicted concentration narrows on the truth. The synthetic data set represents a simplified system of fluvial hydrofacies, with a north to south hydraulic head gradient, no-flow conditions at the eastern, western and underlying boundaries, recharge on the surface, and discharge to a stream. The model domain is 6000 m long, 2000 m wide, 60 to 80 m thick (depending on the elevation of the water table in the uppermost block), and is represented by a grid of 30 x 20 x 6 cells. A landfill and stream are located as shown in Figure 4(a).

A synthetic "true stratigraphy" is constructed consisting of six layers of geologic materials designed to mimic a generic fluvial geology. Five hydrofacies are used in the synthetic data set: (1) medium-grained alluvial fan material, (2) fine-grained overbank deposits, (3) fine to medium-grained splay/chute deposits, (4) medium-grained channel deposits, and (5) coarse-grained channel deposits [Figure 4(a)]. All of the hydrofacies have a constant value of effective porosity and lognormally distributed (but not spatially correlated) hydraulic conductivities with a standard deviation of 1/6 in natural log space (thus K varies by a factor of three within any one facies and there is no significant overlap of K's between facies). These simplified hydrofacies properties are not representative of a field situation; however, a complex hypothetical distribution is not necessary to illustrate the process presented herein. The percentage that each hydrofacies represents in the synthetic system and their hydraulic properties are presented in Table 1.

Facies Volume
Percent Mean K^*
(m/d) Effective
porosity
Alluvial fan 6 50.00 0.25
Overbank 61 0.05 0.08
Splay/chute 11 5.00 0.15
Medium channel 15 50.00 0.25
Coarse channel 11 500.00 0.30
Table 1. Facies Statistics for Synthetic Model

**Table 1. Facies Statistics for Synthetic Model**
Facies	Volume Percent	Mean K^* (m/d)	Effective porosity
Alluvial fan	6	50.00	0.25
Overbank	61	0.05	0.08
Splay/chute	11	5.00	0.15
Medium channel	15	50.00	0.25
Coarse channel	11	500.00	0.30

^* All K distributions have _lnK = 1/6

The alluvial fan hydrofacies are limited to the western edge of the alluvial valley, and are discontinuous. The main channel hydrofacies are oriented north-south and are further east in sequentially younger layers. Coarse-grained channel hydrofacies are discontinuous and occur on the outside bends of meanders. The medium-grained channel hydrofacies are more continuous. Fine to medium-grained channel hydrofacies are located along the main channels to represent splay deposits. Some medium-grained hydrofacies are located within the floodplain area to represent chute cutoff oxbow lakes. Fine to medium-grained hydrofacies of similar distribution represent neck cut-off oxbow lakes. This definition and spatial distribution of hydrofacies is a gross simplification of fluvial hydrofacies as presented by Allen (1965), Bridge and Leeder (1979), and Walker and Cant (1984).

A "true" steady state flow system is obtained with heads fixed at 100 m on the northern and 90 m on the southern boundary. An average recharge of 2 x 10^-4m/d is applied to the top of the domain with three times that rate applied at the margin of the valley (where overland flow from surrounding uplands infiltrates the system), and a zero rate near the stream which is a discharge area. The 31 domestic wells [Figure 4(a)] penetrate to various depths and pump at low rates (their discharges are normally distributed with a mean of 1 m³/d and a standard deviation of 1/6 m³/d). A head dependent flux boundary represents the stream [Figure (4a)]. The USGS code MODFLOWP (Hill, 1992) is utilized for the forward and inverse ground-water flow modeling. The forward simulation of the synthetic system using MODFLOW yields heads at 37 wells (31 domestic and 6 monitoring wells near the landfill) and ground-water discharge to the stream between gaging stations. Multilayer heads from MODFLOW were composited by averaging head of each layer weighted by transmissivity of the layer.

A "true" event of a breached landfill with recharge through the landfill at a dimensionless concentration C_o = 1.0 is simulated for a conservative contaminant with constant dispersivities of 5 m, 1.5 m, and 0.25 m in the longitudinal, transverse horizontal and vertical directions, respectively for all hydrofacies. These dispersivities reflect advective variation at small scale within each hydrofacies, while larger scale dispersion is attained through the use of heterogeneity in the flow model. The simplification of similar small-scale dispersivity for each facies does not interfere with the illustration of the analysis process presented herein. The MT3D code is used with the MOC (Method of Characteristics) option and the fourth-order Runge-Kutta algorithm for particle tracking to simulate transport (Zheng, 1991). Time steps on the order of 10 days yield mass balances on the order of 1%. Forward transport modeling yields the "true" dimensionless concentration of 1.2 x 10^-2 (i.e. C/C_o = 0.0012) at the proposed well location 50 years after the landfill began leaking. For the purpose of discussion we assume the regulatory level for this contaminant is C/C_o = 1.0 x 10^-3.

INVERSE MODELING

Performance of inverse modeling on the true stratigraphy is evaluated before exploring the possibility of using inverse modeling to eliminate some of the stratigraphic realizations. Because we have a synthetic truth to work with, "perfect" input can be provided (i.e. exact zonations, heads, flows, and initial values for the parameters) to the inverse model, in which case the model quickly converges without substantial change of the parameter values. However, stochastic simulation of parameter values within the hydrofacies is not considered here because it is not necessary in order to illustrate the uncertainty reduction process, thus the inverse model estimates only one value of K for each hydrofacies (the mean of the lognormal distribution of K is used as the starting value for the inverse run). In addition, because use data are not typically available for domestic wells, each is assigned a discharge of 1 m³/d. Consequently, the inverse model undertakes a number of iterations attempting to match the head distribution that resulted from a system of heterogeneous K within hydrofacies to heads produced by a system in which each hydrofacies is represented as internally homogeneous. If desired, MODFLOWP can accommodate stochastic variation of K within zones.

Variance of field observations must be defined for an inverse modeling procedure. The specified variance must reflect not only measurement error but also error related to the fact that the grid block represents average head over an area while the head measured in the well represents the head over the open interval at one location. Interpolation of head between grid blocks is used to determine the head that will be compared with the measurement; however, such interpolation does not consider the flow equations, thus a variance greater than that associated with measurement error must be specified for observations. In short, the selection of accuracy and confidence interval is somewhat arbitrary in both field and synthetic studies. How reasonably the user assesses this uncertainty is reflected by the code output. For this study, true head observations are specified with 90% confidence that the measured head is within 1.2 m of the head it represents in the model. The one flow observation (groundwater seepage into the stream between gage stations) is defined with 90% confidence that the measurement is within 400 m³/day of the true flows which are 10,500 and 15,300 m³/d at the up and downstream stations, respectively.

PARAMETER ESTIMATION

Initially, hydraulic conductivities of all five hydrofacies were independently estimated by the inverse model along with vertical anisotropy, recharge, and streambed conductance. The number of estimated parameters was eventually reduced to five for reasons given below. The data do not provide a basis for estimating K of all hydrofacies independently (perhaps due to the minimal number of wells in the alluvial fan facies), therefore the alluvial fan and medium-grained channel deposits (hydrofacies 1 and 4) are represented by one estimate because of their hydraulic similarity. The K of overbank deposits (hydrofacies 2) and the recharge rate are correlated with a coefficient of 1.0. This is a typical situation encountered in many ground-water models but often goes unrecognized because calibration is undertaken by trial-and-error procedures; thus, correlation coefficients are not calculated. In this model, the correlation was addressed by fixing the recharge rate at the "true" value and estimating K of hydrofacies 2. The parameter estimation is insensitive to the streambed conductivity, hence it is omitted from the estimation process. This situation is also common in ground-water modeling but typically goes unnoticed because sensitivities are not calculated for all estimated parameters. Given these adjustments, the parameter estimation is performed on 5 parameters: the vertical anisotropy of K and the K of the five hydrofacies with hydrofacies 1 and 4 grouped together. Prior information on the estimated parameters is defined as the mean of the true distribution with a standard deviation of 0.5 of the ln of K to constrain the estimates of K to reasonable values and a standard deviation of 0.2 to constrain estimates of the vertical anisotropy factor. Standard deviations on prior information can occasionally be obtained from statistical analysis of many field tests (e.g., based on numerous specific capacity tests conducted within an area); however, often the value is based on qualitative judgement and is employed as a tool to prevent the code from estimating parameters substantially different from those measured in aquifer tests or from what is deemed reasonable for a given lithology. Such decisions are typically made, but not quantified, when trial-and-error parameter estimation is undertaken.

GEOSTATISTICAL SIMULATION AND CONCENTRATION DISTRIBUTIONS

Semivariograms are developed with hydrofacies data from the 37 wells and stochastic simulations are conducted with the semivariograms and hard data using a modified version (McKenna, 1994) of ISIM3D (Gomez-Hernandez and Srivastava, 1990) in the public domain software UNCERT (Wingle et al., 1994), resulting in 100 realizations. Flow is simulated in these configurations using "field estimates" of K equal to the "true" mean ln of K for each hydrofacies. Using MT3D, the range of concentrations at the proposed well location after 50 years of contaminant leakage is broad and the "truth" (indicated by the arrow) falls within the distribution (Figure 5A). The vertical line at the center of the gray bars indicates the mean of the concentration distribution for all of the realizations and the gray bars represent one, two, and three standard deviations. The distribution straddles the assumed regulatory level (log 0.001 = -3). There is no clear case of the range of possibilities falling entirely above or below the regulatory level and the decision of whether to use a well or a pipeline is not obvious. One can only conclude that the uncertainty must be reduced. Either additional analyses or data collection must be undertaken to narrow the distribution. Given that analysis is less expensive than data collection, it is desirable to use the available data to reduce the uncertainty as much as possible. Once this is achieved, the results can be used to design further data collection if the distribution is still not narrow enough.

Fig. 5. Histograms of predicted concentrations at the proposed well location 50 yrs after breach of the landfill. A based on hard data and field estimated hydraulic parameters; B: based on hard and soft geologic data and field estimated hydraulic parameters; and C: after eliminating configurations from B and estimating parameters using inverse modeling.

When inverse flow modeling (MODFLOWP) is used to estimate the hydraulic parameters of the system in each of these initial 100 realizations, the results do not produce the proper relative order of hydraulic conductivity for the hydrofacies. The relative order of hydraulic conductivity is determined from hydraulic testing of the hydrofacies in the field and observations of the material character: coarse-grained channel deposits have the highest hydraulic conductivity, followed by medium-grained channel deposits and alluvial fan deposits, fine-grained splay/chute deposits, and overbank deposits. Parameter estimations from this initial inverse modeling consistently indicate that the hydraulic conductivity of the alluvial fan hydrofacies is lower than that of the overbank hydrofacies.

It is concluded that the simulated alluvial fan hydrofacies are too well-connected in these initial 100 realizations. As a consequence, the inverse modeling procedure is estimating a low K for the alluvial fan hydrofacies in order to fit the "measured" heads which exhibit higher gradients produced by the "true" discontinuous zones of high K. Inspection of the realizations reveals the alluvial fan hydrofacies to be either completely continuous along the left boundary of the domain or connected to the high hydraulic conductivity hydrofacies of the continuous channel hydrofacies near the centerline of the domain in all of these first 100 realizations [Figure (4b)].

INCORPORATION OF EXPERT OPINION AS SOFT DATA

Based on this observation, soft data are added to the simulation process in the form of a decreasing probability of occurrence of alluvial fan material with distance from the surface manifestation of each alluvial fan. For each 10 feet of depth directly below the alluvial fan, probability of occurrence of alluvial fan hydrofacies is specified as 95%, 80%, 50%, 10%, and 0%. On the periphery of the fan a 95% probability of occurrence 100 feet to the east and 200 feet to the north or south is specified (distances are limited by grid spacing in the model) with a zero probability at greater distances. A 10% probability of occurrence is specified at locations underlying these peripheral zones. The same semivariograms are used with both the hard and soft data to simulate 400 realizations. Again, using MT3D, the concentration at the proposed well location after 50 years of contaminant leakage is calculated for each realization. Uncertainty is slightly reduced (Figure 5B); that is, precision is improved (note the narrower gray bars). The distribution still includes the "truth", so accuracy is maintained. However, the distribution of predicted concentrations is still broad, straddles the regulatory level, and does not provide a definitive basis for a decision.

ELIMINATION OF REALIZATIONS USING INVERSE FLOW MODELING

Many of the 400 realizations are eliminated using inverse flow modeling based on four criteria. The realization is eliminated if the flow model does not converge to a fairly loose tolerance (0.1 m) on heads at any point during the parameter estimation procedure (30% are eliminated by this criterion). The realization is eliminated if the flow model does not converge to a fairly loose tolerance (1% change in parameter value) on changes in parameters between iterations within the specified maximum number of 20 iterations (20 iterations is large for estimating 5 parameters) (20% eliminated). Examination indicates that realizations eliminated for lack of convergence were progressing in a direction of unreasonable parameter values. The realization is also eliminated if the relative order of hydraulic conductivities is not as expected (7% eliminated). Examination suggests that these realizations have hydrofacies connections that are not similar to the true connectivity of hydrofacies. Finally, the realization is eliminated if the converged parameter estimation is a relatively poor fit to the field measurements (heads and stream flow) based on the value of the sum of squares of weighted residuals (40% eliminated). Sum of squares values range from approximately 180 to about 1200 for this problem. Realizations yielding values greater than 600 (half the maximum value) are eliminated; this is an arbitrary cutoff. Only 2.5% of the realizations remain after these eliminations. Contaminant transport is simulated in these realizations with MT3D using the hydraulic parameters estimated with MODFLOWP. The uncertainty is substantially reduced (Figure 5C) and the distribution mean is near the truth, indicating that inverse flow modeling is a promising approach to incorporating hydrologic data into the simulation process to reduce uncertainty. Although the elimination criteria are subjective, we emphasize that it is preferable to eliminate many acceptable realizations in order to avoid retaining any unacceptable realizations; thus the elimination criteria should be stringent.

SUMMARY

The problem presented herein illustrates a situation in which the use of hard data and multiple indicator conditional stochastic simulation leads to a large uncertainty. Inverse flow modeling identifies too much continuity in hydrofacies based on inappropriate estimates of relative hydraulic conductivity values for the hydrofacies. This observation is used to incorporate expert opinion on the probability of the extent of the alluvial fan hydrofacies as soft data in the simulations, resulting in some reduction of the uncertainty in contaminant concentrations. Finally, inverse modeling eliminates geologic realizations that do not yield a close fit to the hydrologic data, substantially reducing uncertainty.

Given that hydrogeologists must make decisions based on incomplete information, it is necessary that these decisions incorporate a recognition of uncertainty. Stochastic modeling of geologic uncertainty provides a means of quantitatively addressing uncertainty in the response function (advection and dispersion). As shown in Figure 2, this paper presents incorporation of disparate types of information for reduction of uncertainty. This "data fusion" forces the resulting predictions of flow and transport to be consistent with all available data. Furthermore, this work demonstrates two techniques for incorporating disparate types of data into site assessment: use of soft data in stochastic simulation of geologic units and inverse flow modeling. In this case soft data are in the form of expert opinion on geologic conditions (which is identified through inconsistent results of inverse flow modeling). The same approach can be used to incorporate geophysical data (McKenna, 1994; McKenna and Poeter, 1994) and other types of expert knowledge about the site. Inverse modeling allows the incorporation of hydrologic data into uncertainty assessment through automated flow model calibration.

The first step in the demonstration of data fusion in this paper illustrates the danger of ignoring some of the available data. Using only hard data, to generate stochastic geostatistical simulation yields realizations of the subsurface which are implausible. These implausible realizations are most likely the result of biased sample data and violations of the underlying geostatistical assumptions of stationarity and ergodicity. Without the incorporation of hydrologic data into the modeling process, these implausibilities might not have been identified and the resulting response function would have remained imprecise. Analyses and data collection must focus on reducing uncertainty (increasing precision) while capturing the truth in the distribution of predicted concentration (maintaining accuracy).

ACKNOWLEDGMENTS

This work was partially supported by the U.S. Bureau of Reclamation under agreement number 1-FC-81-17500.

REFERENCES

Allen J.R.L., 1965, A review of the origin and characteristics of recent alluvial sediments, Sedimentology, v. 5, pp. 89-191.

Bridge, J. S. and M. R. Leeder. 1979. A simulation model of alluvial stratigraphy. Sedimentology. v. 26, pp. 617-644.

Gomez-Hernandez J. J. and R. M. Srivastava. 1990. ISIM-3D; An ANSI-C three dimensional multiple indicator conditional simulation program. Computers in Geoscience. v. 16, no. 4, pp. 395-414.

Hill, M. C., 1992. A computer program (MODFLOWP) for estimating parameters of a transient, three-dimensional, ground-water flow model using nonlinear regression. USGS OFR 91-484.

McKenna, S. A., 1994. Utilization of soft data for uncertainty reduction in groundwater flow and transport modeling. Colorado School of Mines PhD dissertation T-4291.

McKenna S. A. and E. P. Poeter, 1994, Simulating geological uncertainty with imprecise data for groundwater flow and advective transport modeling. In: Stochastic Modeling and Geostatistics: Case Histories and Practical Examples. J. Yarus and R. Chambers, eds. American Assoc. Petroleum Geologists Special Publication, Computer Applications in Geology, no. 3.

Olhoeft, G.R., 1992. Data fusion. Chapter 3 in Workshop on Non-invasive Geophysical Site Characterization. Compiled by E. Van Eeckhout and C. Calef. Los Alamos Report LA-12311-C UC-903. pp. 4-6.

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.

Poeter, E.P., and S.A. McKenna. 1994. Geostatistical simulation and inverse flow modeling to reduce uncertainty associated with flow and transport predictions. Proceedings of the International Ground Water Modeling Center's 2nd conference, Fort Collins, CO.

Walker, R.G. and Douglas J. Cant. 1984. Sandy fluvial systems. In: Hydrofacies Models, 2nd ed., R.G. Walker, ed. Geoscience Canada. Reprint Series 1, pp. 71-90.

Wingle, W.L., S.A. McKenna and E.P. Poeter, 1995, UNCERT: A Geostatistical Uncertainty Analysis Package Applied to Groundwater Flow and Contaminant Transport Modeling, Colorado School of Mines.

Zheng, C. 1991, MT3D, a Modular Three-Dimensional Transport Model User's Manual. S. S. Papadopulos & Assoc., Rockville, MD