Groundwater Flow Simulation by a Stochastic Representation of Soil

A Monte Carlo technique is utilized to incorporate the uncertainty in media characteristics to the solution of a groundwater flow problem. This technique involves the repetitive solution of a significant number of equiprobable representations of the soil medium. Probability and statistics are utilized to model the soil parameters and study the significance of the output. An existing computer code was adapted and significantly modified to allow characterization of the media's hydraulic conductivity (permeability) as autocorrelated and statistically homogeneous. A first order nearest neighbor model was selected to affect the autocorrelation of this parameter within the finite difference mesh. The statistical homogeneity considers that the distribution of hydraulic conductivity values within the mesh comes from a log normal probability density function. The selection of hydraulic conductivity value at any mode of the mesh is stochastic within the framework of the autocorrelation and statistical homogeneity of the mesh aggregate. The computer code takes the stochastically generated hydraulic conductivity field and the boundary conditions and utilizing an iterative alternating direct implicit solution determines the hydraulic head values at each mode and flow rate through the medium. An array of these results are produced for each of the equiprobable representations of the soil medium. Mass transport through the region is simulated as a combination of advection and a stochastic simulation of microscopic or particle scale dispersion. A water particle is released from a preselected location along the upgradient boundary. The particle moves toward the downgradient boundary under the influence of advective forces caused by the differences in hydraulic head and the stochastic simulation of microscopic dispersion. This simulation of microscopic dispersion displaces the particle parallel to and perpendicular to the advective transport direction based on laboratory scale dispersivities. The computer code establishes arrays for the particle location at a predetermined time after start as well as the location along the downgradient boundary and total travel time upon completion of transit of the region. Uniform flow results in most of the regions considered although some alternate configurations were considered. An effective hydraulic conductivity is calculated on the basis of flow rate. After application of a shape factor this value was found to be slightly less but closest to the geometric mean of the hydraulic conductivity distribution thus confirming earlier work. An alternative effective hydraulic conductivity calculated on the basis of travel time was also determined. This value was generally less than the other effective hydraulic conductivity value but again after application of a shape factor the value was best estimated by the geometric mean. These results suggest that the mean flow rate and mean travel time may be estimated by the use of the shape factor from a flow net solution or the method of fragments and the geometric means of the hydraulic conductivity. The results of the simulations indicate that macroscopic or field scale longitudinal and lateral dispersion is significantly affected by the standard deviation of the hydraulic conductivity distribution. Region size, hydraulic gradient and time interval were found to cause lesser

Mass transport through the region is simulated as a combination of advection and a stochastic simulation of microscopic or particle scale dispersion.
A water particle is released from a preselected location along the upgradient boundary. The particle moves toward the downgradient boundary under the influence of advective forces caused by the differences in hydraulic head and the stochastic simulation of microscopic dispersion. This simulation of microscopic dispersion displaces the particle parallel to and perpendicular to the advective transport direction based on laboratory scale dispersivities. The computer code establishes arrays for the particle location at a predetermined time after start as well as the location along the downgradient boundary and total travel time upon completion of transit of the region.
Uniform flow results in most of the regions considered although some alternate configurations were considered. An effective hydraulic conductivity is calculated on the basis of flow rate. After application of a shape factor this value was found to be slightly less but closest to the geometric mean of the hydraulic conductivity distribution thus confirming earlier work. An alternative effective hydraulic conductivity calculated on the basis of travel time was also determined. This value was generally less than the other effective hydraulic conductivity value but again after application of a shape factor the value was best estimated by the geometric mean. These results suggest that the mean flow rate and mean travel time may be estimated by the use of the shape factor from a flow net solution or the method of fragments and the geometric means of the hydraulic conductivity.    Table 1 Section II Table 1   Table 2 Appendix A A. I Appendix B Table 1 LIST                                      Perhaps the technique most widely used today is that of flow nets (10,11). Its greatest advantage is its simplicity and ease in use. However, this technique is limited in its ability to vary soil hydraulic conductivity and groundwater surface or free surface.
A similarity with all of these techniques is the lack of flexibility in handling boundary conditions and/or varying soil property characteristics. When these techniques are utilized, the consequences of the required averaging and other approximations are often unknown.
With the recent interest in this area from the multitude of scientific disciplines and with the advancements in computer technology in the last decade, major changes in the way these flow problems are analyzed have A review of these methods may be found elsewhere (24,35).
The FEM and FDM methods have allowed not only more complex boundary conditions and input parameter characterization but also allowed introduction of new types of considerations. Of these, stochastic generation of a correlated hydraulic conductivity distribution to characterize the flow regions, stochastic simulation of microscopic dispersion and the utilization of probability concepts or theory to analyze input and output are major considerations of this work.
The flow region considered in this study is two dimensional. It should be pointed out that the term hydraulic conductivity is utilized rather that the more common geotechnical engineering term of permeability to avoid confusion with workers in other disciplines. The term permeability is reserved for the medium's flow characteristic without the effects of the pore fluid's characteristics.
This is comparable with the work in the geologic and petroleum industries and with similar terminology in such areas as heat, electrical and chemical flow.
The most significant of the soil parameters governing groundwater flow is hydraulic conductivity. The range of values for hydraulic conductivity of normal soil deposits is as much as 13 orders of magnitude.
This variation is large for a soil property. By comparison, porosity may only vary between zero and one, unit weights by less than one and dispersivity by two orders of magnitude .
The variation within one single soil type can also be large. A variation of three orders of magnitude within a soil designation, such as sand or silt, is possible . The engineer has, by necessity, to categorize similar soil encountered in soil explorations into various soil strata. The variation in hydraulic conductivity in one soil stratum may be as much as one order of magnitude. This suggests that knowing the order of magnitude of hydraulic conductivity in a soil stratum may be enough and to expect to know the value at one particular point in a soil mass to within one order of magnitude is somewhat presumptive. Indeed we should expect considerable variation within even a "homogeneous" soil.
It is this potential for wide variation and therefore the implied uncertainty which has led engineers and geologists to turn to probability and statistical techniques to characterize hydraulic conductivity. The value of hydraulic conductivity in situ is considered to be a random variable in this methodology. This implies that the exact value at any point cannot be predicted with certainty a priori. There is uncertainty in any estimate prior to experimentation.
A random variable may be described by a probability density function (PDF) which simply relates the value of a variable with its corresponding probability. A Monte Carlo solution is accomplished by generating a significant number of equally probable representations of the flow medium from the PDF, repetitively solving each for the unknown desired results, and analyzing the results to determine the "most probable" results as well as a measure of the magnitude of variation of that result.
There are several well studied PDF's which have been found to adequately characterize soil properties (6,20,24). The Gaussian or normal distribution is one of the most commonly studies PDF's. This distribution is completely defined by the mean and standard deviation. This distribution has been found to model well those processes which are additive of many processes or random variables. This is the premise of the Central Limit Theorem. An example of a soil property which has been found to follow this PDF is porosity.
Another PDF which has been found to model some soil random variables is the log normal distribution . This distribution represents well those processes which are the result of multiplying other random variables .
Hydraulic conductivity has been found to be adequately characterized by a log normal PDF (13,29,32). A useful identity is that the PDF of the log of a random variable with a lognormal PDF is normal .
The physical processes which result in soil deposits may be quite random and conditions predominating at one time and place will be modified elsewhere.
It should be expected, however, to f i nd similarities at adjacent locations when compared to locations separated by some distance. This similarity is often quantified by calculation of the autocorrelation of the random variable.
Techniques to accomplish this spatial autocorrelation include a nearest neighbor model, trend analysis, spectral density and variations of these. The nearest neighbor models develop equations for the dependency of the variable at adjacent locations and solve these equations for model development. This technique is adopted here so that further discussion will follow .
The trend analysis technique was developed at Northwestern University by Krumbein and Whitten (23,42). The technique assumes that the random variable at any point may be characterized as the summation of a trend component (dependent on location) and a random component. Koch and Link (22) present an excellent introduction to this technique. Tabba and Yong (36) have developed a working model and in a companion paper (37) demonstrate its application to a real problem. This technique interpolates and extrapolates information on the random variables to all points in the flow media.
The spectral density technique recognizes the Fourier transform pair consisting of the spectral density function of the random variable S(f), and the autocorrelation function, R(x). S(f) describes the frequency information for the process and R(x) the autocorrelation. This approach may describe the variations in soil properties more precisely than the other. However, its inherent complexity and the need for large amounts of data has discouraged its use for real problems in groundwater flow.
The effects of dispersion during advective transport is becoming more critical as interest in tracking movement of pollutants and predicting concentrations in groundwater systems increases.
Others have studied the problem of dispersion in groundwater flow stochastically (17,35,41 (27). An alternate technique by contrast uses spectral analysis assuming a vertical variability of hydraulic conductivity (17).

MODEL INPUT DEVELOPMENT
The methodology presented in this paper adopts the general premise that a soil stratum may be considered statistically homogeneous with the hydraulic conductivity having a lognormal PDF. A result of this premise is that the variation in the soil property may be represented by a curve similar to those presented in Figure 2. The upper curves illustrate PDF's of hydraulic conductivity for two different magnitudes of standard deviation but with the same mean value. These curves have lognormal PDF's. The lower curves illustrate the change in form of these PDF's when the log of hydraulic conductivity for these same two distributions are plotted. These curves have normal PDF's. It should be noted that a convenient aspect of the lognormal distribution is that negative values of the parameter are avoided.
The method selected to achieve the auto-correlation for the hydraulic conductivity values at adjacent nodes in the finite difference mesh is the stochastic nearest neighbor-model originally adapted by Smith (31) for use in groundwater flow applications.  (5,7,9,25,42). A computer code was written to generate statistically independent realizations of hydraulic conductivity for each of the required Monte Carlo simulations.
Smith (31) presents a detailed development of the stochastic nearest neighbor model and much of this information is also available elsewhere (32,33,34). A swnmary is presented here.
The finite difference method utilized to solve the groundwater flow problem requires as input hydraulic conductivities at regularly spaced grid points. These points or nodes are located at the center of a block of constant hydraulic conductivity, the so-called block centered approach. Figure 3    Smith (31) assumed that the ~ term is divided by the number of adjacent nodes. Figure 3 also presents the situation along the upper boundary.
Here this approach results in the following equation: Since a log normal distribution of hydraulic conductivity was adopted, the {s}matrix generated for each simulation will have a mean value of O, a standard deviation of 1 and a normal distribution. A transformation will later be made so that the hydraulic conductivity values will be log normally distributed. Smith (31)  To alleviate this limitation, the flow region may be considered a compo-site of two or more "blocks" that contain statistically homogeneous and internally correlated hydraulic conductivity meshes. For example, the shaved 5 x 5 region may be systematically added to 3 other statistically equivalent regions resulting in a 10 x 10 region. The computer time to add these regions together is much less than that required for inverting the larger matrix.
The program utilized to solve the flow problem was developed by Reiter (29).
This code was adapted from other available work for the two dimensional cross-sectional model (38).

MASS TRANSPORT SIMULATION
A computer code was developed to simulate transport across the hydraulic field resulting from a steady state solution. Transport results from hydraulic advection, a stochastic simulation of microscopic dispersion and dispersion resulting from variations in hydraulic conductivity. The major objective of this effort was to analyze the effects of hydraulic conductivity variations on dispersion.
Others (17,35) have indicated that the variation in porosity is much less than the variation in hydraulic conductivity and probably does not significantly influence dispersion. The value of the porosity has therefore been assumed constant throughout the region.
The finite difference method assumes that the hydraulic head, hydraulic conductivity, porosity and therefore the resulting pore water velocity are constant within each individual block in the flow region. These velocities are calculated for each conductivity simulation and for each block within the region by the method suggested by Sauty (30).
Microscopic dispersion is incorporated into the model in a manner first described by Ahlstrom, et al (1) and then modified in part by Smith and Schwartz (35). In this method, dispersion is treated as a stochastic phenomenon. The water particle is assumed to undergo a random displacement in the direction of flow and another random displacement perpendicular to flow. These displacements are calculated from the following two equations: In these equations, DL and DT are the displacements in the longitudinal and transverse directions, respectively, relative to the direction of advective transport. DLC and DTC are the longitudinal and transverse dispersion coefficients. RANOL and RANOT are random numbers normally distributed with a mean of zero and a standard deviation of one then divided by six to keep the range of values in the -0.5 to +0.5 area.
Kelly (21) developed a direct relationship, which has been adopted for this study, between grain size and microscopic dispersivity. This relationship may be expressed as: o = 0.08 n 50 (12) Where o = dispersivity, in meters, and = mean particle size in millimeters The longitudinal dispersion coefficient, DLC, is taken equal to the dis- Since the advective transport direction will not always coincide with the X and Y axes the displacements in the X and Y directions are obtained by a transformation.
The computer code simulates water movement from the upgradient side of the region to the downgradient side. A particle is released from the upgradient boundary and moves across the region as dictated by the advective displacements and stochastic dispersive displacements. When an impervious boundary is encountered, the particle is reflected back into the region. Eventually the water particle encounters the downgradient boundary. The location, travel time and other data are recorded and available for statistical analysis.
Monte Carlo simulations result in a great deal of output. It was necessary to develop a large data collection and analysis system to deal with the data. Data from this study was stored in disk data files as it was produced and saved for use and analysis later. The Statistical Analysis System (SAS) has been used in the analysis of the data.
The data output to disk files was varied but consists of information on the hydraulic conductivity means and correlation for each simulation,  (15) where Q is the flow quantity, h 1 is the total head loss within the region, L is the total length over which that head loss occurs, H is the height of the flow region and S is a shape factor computed in accordance with a technique suggested by Warren and Price (40) which will be presented later.
The value of S is a shape factor which may be directly compared with the shape factor, ND/NF' in a flow net solution or the form factor in the K . Figure 4 presents a typical layout of a flow region with approxq imately horizontal uniform flow from left to right. The flow is not precisely uniform due to the dispersion present in the flow simulation.
The velocity of flow, v, may be calculated by the following formula: The mean velocity of all flow particles traversing the flow region, v may be calculated by dividing the length of the particle travel, LT' by the mean travel time, t, as below: The effective hydraulic conductivity on the basis of travel time, Kt, has been calculated by setting v equal to v in the above two equa- tions, introducing the height of the region by inserting in both the numerator and denominator and solving for Kt resulting in the following: (18) where LT is the total length of particle travel and for the uniform flow case will equal L, n is the porosity, t is the mean travel time and L, H / and h 1 are as defined above. St is a shape factor dependent on boundary / conditions and is computed in a manner similar to S and will be discussed later.
St is a modified shape factor which may be compared to the other shape factors, S, and the shape factor in the flow net solution or by the method of fragments.
Twenty-five variations of flow region size and characteristics were studied to determine the relationship between these calculated means. A typical plot of these results is presented in Figure 5 and the results are summarized in Table 1.    The results of these computations support the observation that the geo- Warren and Price (40) have indicated that the geometric mean of the hydraulic conductivity is a good approximation of the effective hydraulic conductivity even for non-uniform flow provided a shape factor is utilized to adjust the value for the shape of the flow region. They This effective value of hydraulic conductivity should then be divided into the singular known value and the result is the shape factor. The shape factor is then multiplied by the effective hydraulic conductivity before comparing with the geometric mean . Examples of these computations are presented elsewhere (2) .

Region
The results of this study tend to support this technique. For the cases considered, the adjusted values of effective hydraulic conductivity generally fall between the geometric and harmonic means and is closest to the geometric mean as presented in Table 1.
The importance of this finding is that the mean flow rate may be estimated by rearranging equation 15 to solve for the flow rate, Q, as below: The geometric mean should be used for the value of K and the shape factor, S, may be determined by a flow net solution or the method of fragments.  Initially consider the mean and standard deviation of the hydraulic conductivity formulations.    Fried (15) gives an expression for the longitudinal dispersion coefficient, D 1 , assuming that the development of particle concentration at some distant point is normally distributed. He illustrates that this coefficient is directly related to the variance of the travel times.
The equation developed for a long travel distance relative to the dispersion is: A record was maintained of the starting and ending y coordinate of each particle. A normalized ending y coordinate was then computed along with a cumulative mean and cumulative mean standard deviation as each new record of water transport was added.
The same was computed for the y coordinate at the estimated time to transit half the region. A typical plot of these results is presented in Figure 11. These results also indicate that the 400 water particles are sufficient to reach a consistent mean value for these variables.
Normality Check on Travel Times The ensemble of travel times was utilized to plot a breakthrough curve. A check for the goodness of fit of the assumption of a normal distribution for the ensemble of travel time records was made.

A SAS
program was written to test the null hypothesis that the data could be a random sample from a normal distribution. A Kolomogorov D statistic was computed and compared with critical values. A Kolomogorov statistic was utilized rather than the more common chi square check since it is more easily adapted to continuous uncategorized data.   At the ends of the flow region the head is held constant and therefore the standard deviation is zero at these nodes. At nodes further away, there is less influence from the boundary and therefore the standard deviation and variance are greater.
The number of Monte Carlo simulations required to reach a fairly consistent value of the hydraulic head was also investigated. Figure 14 presents a typical plot of the results of one such investigation at an interior node. This plot indicates that the initial simulation had a value very close to the resulting mean value and after three simulations        The relationship between this quantity and the standard deviation of the hydraulic conductivity is of more significance. Figure 17 presents this relationship which indicates a direct relationship. This is significant since the uncertainty in the hydraulic conductivity value directly affects our uncertainty in the value of the hydraulic head value within the region.

Macroscopic Dispersivity Sensitivity
The macroscopic or field scale longitudinal dispersivity has been shown to be a function of the standard deviation of the travel time. The       The procedure to develop the site specific soil parameters is essentially unchanged. Field programs are utilized to measure insitu soil parameters and to secure samples for laboratory testing. The field and laboratory tests will be analyzed so that a probability density function of the hydraulic conductivity may be formulated. Other parameters such as porosity, and laboratory dispersivity characteristics must also be developed.   IN l40 (26) where ta 12 is the value of the cumulative distribution function of the I student t distribution at a /2 significance, s is the value of the standard deviation and N is the sample size.
The mean value would then be expected to lie between the limits of the sample mean (4.57) plus or minus the error (0.07) i.e., 4.50 and 4.64.
The effective hydraulic conductivity, K , was calculated by dividing the q flow by the mean hydraulic gradient and area. The flow rate may then be estimated with confidence limits of 95% by multiplying by these constants.
In so doing the flow limits of 38 to 39.5 ft 3 /day were estimated with 95% confidence limits.
The same methodology may be utilized to develop confidence limits on the travel times. The mean and standard deviation of the travel time is 1576.2 and 146.9 days. The magnitude of error of the mean with 95% confidence may be estimated as: The mean value of travel time is therefore expected to lie between 1562 and 1591 days.
The standard deviation of the travel time may also be estimated with confidence by the expressions given below (26): (28) where Za/ 2 is the value of the cumulative density function at a signif- The water particles are expected to follow a normal distribution so that a hypothetical breakthrough curve may be developed for the extremes of these parameters. The result of such an exercise is presented in Figure   23. This allows the engineer to estimate not only the mean arrival time with confidence but also the distribution of travel times.   The use of a computer model which characterizes flow in a porous medium as a stochastic process is demonstrated on two field problems.
The computer model allows the hydraulic conductivity (permeability) to be varied within the blocks of the finite difference mesh according to a predetermined autocorrelation structure and a log normal probability density function.
Mass transport is simulated as a combination of advection and microscopic or particle scale dispersion. The model allows development of confidence in the flow rate and travel times within the region.

INTRODUCTION
Geotechnical engineers are asked to make predictions of groundwater flow in a variety of boundary conditions and soil media. Sometimes there is considerable data available, but more often the problem is so large, the boundary conditions and soil medium so complex, that the resources available can not be expected to identify in detail the conditions insitu.
Nevertheless, predictions are required.
The geotechnical engineer has  generally be included by performing a transformation.

Anisotropy can
Perhaps the technique most widely used today is the flow net. This technique has been available for nearly half a century and been utilized extensively (3,4). Its greatest advantage is its simplicity and ease of use. The technique suffers from its limited capacity to vary the soil's hydraulic conductivity and freewater surface, however.
Numerical methods allow more realistic problem formulations. In most cases, however, the problem must be simplified into regions where the soil media have constant hydraulic conductivity. The companion paper (2) discusses the problems with this type of approach as well as other more realistic modifications.
The computer model discussed in the companion paper (2) requires input information normally required in any of the foregoing techniques as well as some additional information.
Initially the soil medium must be divided into similar strata. Boundary conditions must also be evaluated and locations for constant head and no flow or impervious boundaries identified.
The hydraulic conductivity of the soil medium is presumed to be characterized stochastically by a log normal probability density function (PDF).
An estimate of the mean and standard deviation of the hydraulic conductivity must be determined. The method of estimating these values is to compile all the data available and assign weighting factors for the various sources and test procedures used. Perhaps more weight will be given to insitu tests than lab tests or correlations from other measured parameters. An estimate of the integral scale or distance over which the value of hydraulic conductivity is positively correlated must also be determined. This may be done by plotting the discrete values of hydraulic conductivity from samples taken in close proximity vertically and laterally, computing correlation coefficients and the resulting integral scale. (10) Other data required are median grain size, the ratio of lateral to longitudinal microscopic or particle scale dispersivity and soil porosity.
The median grain size may be determined from grain size distribution curves. Soil porosity may be estimated from standard penetration tests and soil descriptions. Figure 1 presents a correlation between standard penetration test results and porosity for different soils types. This figure is adapted from Navfac DM 7 (8). On the right ordinate scale the results of information were superimposed correlating density descriptions with standard penetration results from Peck (9).
The methodology to be described herein and presented in detail in the companion paper (2) can best be illustrated by example. The following two examples will demonstrate its use in solving real flow problems.   be set up. This would increase seepage into Lake Mishnock perhaps increasing water levels on a lake already ringed with homes only a few feet above water levels. Other homes between the reservoir and the lake might be similarly affected.    (7) presents a similar hydrogeologic section oriented perpendicular to that in Figure 4 and presents a section across the dike and along the location of potential leakage from the reservoir.
The major problem confronting the design engineer in this area is to quantify potential leakage from the reservoir if an impervious dike were constructed along the proposed alignment. Seepage losses beneath an impervious dike must be estimated with confidence. In addition, design alternatives to reduce leakage to acceptable levels need to be developed  An idealized soil profile assuming a singular stratum of soil was sufficient to characterize the soil profile was developed and is presented in Figure 6.
The soil properties were estimated from the results of a series of laboratory tests. The mean and standard deviation of hydraulic conductivity and mean particle size were estimated from a series of grain size analyses on this soil. The autoregressive parameters were assumed equivalent to the results of field studies based on similar stratified soil by Smith (11). The porosity was estimated from soil descriptions and standard penetration test results.
A second soil profile which characterizes the soil conditions by using two soil strata is presented in Figure 7. The deeper soil stratum has a higher hydraulic conductivity and thus may act as a conduit for flow beneath the dike. Its thickness is limited to one row of the finite difference mesh located along the bottom of the section. . The flow region and boundary conditions in this example do not result in uniform flow so a shape factor must be applied to the effective hydraulic conductivities. The technique to determine this shape factor was included in the companion paper (2).
The results of 40 Monte Carlo simulations for these two flow situations are summarized in Table 1. In addition to the two soil profiles already discussed is a third profile, profile lA, which has the same layout but the soil is assumed to have a constant value for the hydraulic conduc- LECEMD : A -1 oas . I -2 088. ETC .

STAATINO X COORDINATE
LEOEND : A -1 aas. I -:;z 088. ETC .  Confidence limits can therefore not be predicted with accuracy. Even a quick observation of Figure 11, however, would certainly lead one to believe that the mean value of 56. 5 cubic ft/ day is a fairly good estimate of the flow rate.  simulations indicate that approximately 45% of the water exiting at the river will pass through the landfill and be contaminated.
The mean x and y coordinates after 191 days are 218.8 feet and 54.5 feet, respectively.      1 A . .

3.
The presence of a high hydraulic conductivity stratum at depth may significantly increase the flow through the region even when the flow is generally parallel to the stratum.

4.
Estimates of confidence limits may be developed from the results of these analyses when normality is confirmed and assuming that the soil properties adequately characterize the in situ soil.

5.
Even when the distribution of travel times and flow rates cannot be accepted as conforming to a normal distribution, confidence may be developed in the predictions of these results by reviewing the histograms and cumulative mean results.

6.
More work is needed to develop autoregressive parameters for sediments.
It is not likely that these parameters will vary over a large range but more actual field studies are required to develop these parameters.  Perhaps the most widely used technique today is that of flow nets, see , Cedergren (1977). Its greatest advantage is its simplicity and ease in use. This technique is also limited in it's ability to vary soil hydraulic conductivity and groundwater free surface.
Scale models in sand tanks have also been utilized to study the flow of To understand the basic concepts, a simplified discussion in two dimensions follows. The second derivative is approximated as the difference between adjacent first derivatives. Thus for the case considered in Figure A Hydraulic conductivities in natural deposits of soil have been found to vary by thirteen orders of magnitude. Freeze and Cherry (1979) indicate that the hydraulic conductivity of some gravels may be as high as 100 centimeters per second and that of some clays as low as 10 -ll centimeters per second. The hydraulic conductivity of some rocks may even be lower than this.
It is not common for an engineering parameter to vary by this much. It has been found that the same soil type may vary by as much as six orders of magnitude and the same stratum by more than one order of magnitude. This information suggests that simply knowing the order of magnitude of the hydraulic conductivity will be useful and that it is unlikely that we can determine the value at any point to within much closer than an order of magnitude. We also can expect to see a relatively large variation, perhaps an order of magnitude or more, within a "homogeneous" deposit.  ing field pumping tests and outflow tests. The results of these many attempts suggest that no method is very accurate and there are advantages and disadvantages for each of the methods adopted.
Estimates utilizing grain size distr i bution have been tried by many investigators.
Perhaps the most commonly used equation is the Hazen formula given below. (7) In this equation d 10 is the particle grain size which 90 percent of the soil is larger and 10 percent is smaller .
The constant, C, has been estimated to be from 1. 0 to 1 . 5 if d 10 is in millimeters and K is in centimeters per second. Although this equation was developed originally for fairly uniform sands, it has been found to be adequate for most sands and gravels . The major advantage of this equation is its relative simplicity .   Leonards (1962) ity of the grain sizes in the soil. A soil which is well sorted contains mostly one particle size and the hydraulic conductivity will tend to be higher. A soil which is not well sorted will be well graded, the smaller particles filling the voids between the larger particles and thus lower hydraulic conductivity. An example of these is one developed by Masch and Denny (1966) given in the form of a semi-log plot in  In this formula, p is the density and u is the viscosity of the fluid, g is the acceleration of gravity and d is the representative grain size. m Other less significant factors such as particle shape and packing have also been suggested, Todd (1980), for correlation with hydraulic conductivity. The effects of these other factors are relatively insignificant when compared with those already discussed, and therefore are not generally accounted for. In light of the fact that the value of hydraulic conductivity is highly variable, many investigators have chosen to consider this parameter as a random variable. A random variable is a variable whose precise magnitude cannot be determined with certainty, a priori . This is certainly the case for hydraulic conductivity . Probability aspects will be discussed later in this appendix.
Another input parameter has already been mentioned . This is porosity which is defined as the volume of voids divided by the total soil volume, expressed in percent. Porosity is related to hydraulic conductivity but its maximum range of values is only from zero to 100%.
Additionally, typical values for soil are only within the range of 15 to 60 percent. Porosity is utilized to calculate the velocity of groundwater flow since the Darcy velocity is not the actual water particle velocity but really a homogeneized term.
The other remaining input parameters are the microscopic or particle scale longitudinal and transverse dispersivities. The dispersivity is a property of the soil which relates to the spreading of the water particles about a mean direction of movement. This term attempts to incorporate the spreading due to the tortuous path followed as the water particle traverses the flow region. Dispersivity has been found to be an approximate function of the grain size of the porous media, Kelly (1982). Kelly pointed out the difficulty in scaling up sand column dispersivities from the laboratory to a field problem. A discussion on dispersion will follow in a later section.

VI. Stochastic Soil Modeling
Recent work in the study of groundwater flow indicates that engineers and geologists now recognize that it is oftentimes dangerous to characterize soil properties simply with some average value or some other presumably "conservative" estimate. Even the addition of multi-strata flow regions and sensitivity analyses on input parameters are not always satisfactory.
Recently investigators have turned to the use of probability and statistics to facilitate their understanding of these complex problems.
Simple statistics have, of course, long been utilized to characterize soil properties, i.e., the average or mean value. However, it has been only recently that the randomness of the data, to some extent quantified by its standard deviation, has been incorporated into the problem analysis e.g. Harr (1977). Other more powerful concepts are available and are now also being used.
There are many good introductory texts dealing with probability and statistics for engineers and scientists, e.g . , Benjamin & Cornell (1970), Miller and Freund (1977) and Walpole & Myers (1978) . It is not the purpose of this appendix to introduce these concepts since they are readily available elsewhere, however, some discussion will follow on their usage in the study of groundwater flow.
It is very unlikely that a soil property such as hydraulic conductivity or porosity will be uniform in situ.
Variations are inevitable. However, when analyzing groundwater flow, the investigator must, by necessity, recognize similarities and divide the region into like strata to facilitate analysis.
Soil media involves numerous characteristics which may be considered as random rather than deterministic variables. This implies that the exact value of these characteristics or properties at any point is not known and cannot be predicted with certainty, a priori. There is uncertainty in any estimate prior to experimentation.
The principal difference between stochastic and deterministic soil modeling is that the former recognizes variability and attempts to quantify its effects.  which has been found to adequately characterize the characteristics of some soil properties. e.g., Lumb (1974) and Harr (1977). This PDF is for the Gaussian, or normal distribution, which is also the most commonly studied distribution. This distribution is completely defined by two parameters, the mean and the standard deviation.
A unique aspect of the normal distribution proven by probability theory is the Central Limit Theorem. The hypothesis of this theorem is that the normal distribution will approximate a random variable which is itself the sum of many other random variables. Soil properties which are treated as random variables need not be categorized into a specific probability distri~ution with its known PDF. Of particular concern with the strategy of matching data with a known PDF are problems in the tails of the distribution, i.e. , the extreme values of the distribution.
There are several accepted methods to verify that an assumed PDF properly models the data collected. The probability and statistics texts cited earlier treat these methods in detail. These methods may be divided into two main types, graphical and computational. Graphical techniques may be used to roughly determine if the assumed PDF is adequate.
Comparison of histograms generated from collected data and from the assumed PDF such as those in Figure A.8 may be used as a qualitative measure of acceptance .  In spite of these problems fitting data into these known probability distributions, numerous soil properties have been categorized by various investigators. Table A. I compiles the results of several of these attempts.
Although this review should not be considered exhaustive, it does indicate a consensus for some of these soil properties.
Of particular interest to this work is that hydraulic conductivity has been found to have a log-normal PDF and porosity a normal PDF.
With the extreme variability possible for input parameters and output, a measure of confidence in the predictions is desirable. Research work in this area has resulted in several methods for estimating confidence limits on the mean values and estimates of the error in these estimates. Miller and Freund (1977) show that for a normal distribution, the maximum error, E, in the prediction of the mean value of the population is: s/ vN (10) In this inequality Z6. / 2 is a tabulated value for the normal distribution, such that the area under this distribution to its right is equal Also & is the significance level, s is the standard deviation of the sample and N is the number of observations in the sample.
Similarly, they also develop estimates of confidence limits for the population mean, X as follows: . . .
. . .  In this inequality, x is the sample mean, + z~/2 (11) cr is the population standard x deviation and za/2 is the value of (x -X) at which the CDF of the normal distribution equals a/2.
Here B and C are undetermined constants depending on the data. If we plot the data with the dependent variable as the ordinate and the coordinate x as the abcissa, then C will equal the intercept and B the slope of the line. The "best fit" for this line is characteristically determined by the method of least squares.
Thee term may now be evaluated from the residuals, i.e., the variation of the data points from the fitted line. It is generally assumed to have a normal distribution with a mean of zero and a variance equal to the mean value of the squared distance for the "best fit" line.
This method may get far more intricate as the parameter to be modeled requires it. Two and three dimensional models are possible. The coordinate terms may include not only the first power of these coordinates, but also higher powers and cross products in order that the trends be adequately mapped. Tabba and Yong (198la) have developed a working model for a quasi homogenous soil utilizing the trend analysis technique. They combine the bias in the parameter measurement method with the random variation component so that the mean of the random variation will no longer always be zero. Tabba and Yong (198lb)  A thorough presentation of this complex subject may be found in random vibration texts such as Newland (1975) and Bendat and Piersol (1971). J S (f) cos 2 1T xf (x) df (20) c;,oo S (f) = J R (x) cos 2 1T f x dx (21) -co In these expressions, f is the frequency defined as the inverse of the distance, x, the distance between data points or the lag distance.
A random process which is large in the space domain will be small in the frequency domain. This is illustrated in Figure In this equation K is the value of the constant hydraulic conductivity c and the K is the effective hydraulic conductivity for this homogeneous q soil.

VII. MASS TRANSPORT SIMULATION
Mass transport in groundwater through a stochastically generated flow medium has been studied by various investigators, e.g., Smith and Schwartz (1980) and Gelhar, et. al. (1979). The classical development for mass transport in groundwater in a porous medium is given in the literature by various researchers. Freeze and Cherry (1979) indicate that mass transport may involve as many as six mechanisms: advection, dispersion, adsorption, chemical reaction, biological transformation and radioactive decay. The first two are present in virtually all cases while the remaining four may or may not be significant factors in the groundwater transport.
Advection is the movement of the water caused by hydraulic gradients.
Dispersion is the movement of the water resulting from the spreading of the water particles due to diffusion and the natural hetergeneity at the pore scale of the flow medium.   Todd (1980) which may be used to demonstrate this concept. Initially, a wave front is released from the upper portal and movement and distortion is traced. Ogata (1970) presented the solution to this differential equation in terms of the complementary error function, erfc, For all practical purposes the second term in the brackets is negligible and may be neglected.
The Central Limit theorem has been utilized to justify the assumption that the breakthrough curve illustrated in Figure A.13 should be normally distributed, e.g., Bear (1972), Fried (1975). This had been suggested for the displacement relationship provided the water particle is allowed to undergo a significant number of uncorrelated displacements. Fried (1975) presents a method to calculate the macroscopic or field scale dispersion coefficients given the results of a tracer transport experiment and utilizing the assumption of a normal distribution.
Smith and Schwartz (1980) utilized a nearest neighbor stochastic process model of hydraulic conductivity to characterize a uniform flow field .
Water particles were released on one side and recovered along the downgradient boundary. The particles were allowed to move by hydraulically driven advection and a stochastic simulation of microscopic dispersion.
Their results indicated that normality for the breakthrough curve and a constant macroscopic dispersivity could not be confirmed . Gelhar, et. al. (1979) studied the mass transport problem and found that for large distances the macroscopic dispersivity did appear to reach a constant value and appeared to be a function of the properties of the medium. It has been an objective of this work to develop a procedure to more realistically characterize the medium of groundwater flow for a field problem and to arrange the output so that definitive conclusions may be drawn utilizing statistical models.
It is hoped that this will result in an improvement in the confidence that an engineer may have in the predictions of flow quantity and travel times.

B.2 General Technique
In this work, a two dimensional characterization has been chosen because it is the more common and simpler representation. A three dimensional model could have been adopted, but it would have been very difficult to implement to solve any real problems because of the tremendous amount of computation required.
Although the computer has made the numerical methods possible, in practice, three-dimensional solutions are still not commonly done.
A rectangular flow region was chosen for study so that a general problem might be solved. Figure  characterization, region shape and boundary conditions were also studied and will be discussed later.
The problem considered was one of a steady state flow situation. This condition was chosen because many real problems may be considered to be approximated by steady state conditions. The autocorrelation and dispersion techniques may be adapted to transient solutions as well.
The general procedure adopted to study this groundwater flow problem is detailed below: 1.
Select a region size and shape, develop boundary conditions and soil media parameters which approximate the real situation.

2.
Generate a stochastic representation of the flow media.

Utilizing a finite difference program with the Iterative
Alternating Direct Implicit Procedure (IADI), solve for the heads and flow quantity for the steady state condition.

5.
Repeat the above process until a sufficient number of Monte Carlo simulations of the region and water transport has been completed.

6.
Tabularize the data and perform statistical checks of significance on output so that definitive conclusions may be drawn.

B.3 Problem Definition
The initial step requires a definition of the region size and the soil properties. This step is required in any modeling of groundwater flow whether analytical, numerical or otherwise.
An idealized representation of the region in terms of finite blocks of homogeneous soil media is required. The hydraulic conductivity, porosity and head are assumed to be constant everywhere within each fin i te block.
As is often the case, a good deal of experience is required to make the selection of the number of blocks to represent the region. The · {E} matrix is a column matrix of random numbers with a preselected mean and standard deviation. The basis for selecting the mean and standard deviation will be discussed later. ,J An alternative treatment of these boundary nodes was considered. This alternative approach treats these areas in the same way as finite difference grid manipulations at the boundaries resulting in the following equation .
This approach results in greater weight being assigned to the K value horizontally along a soil layer for the end boundaries and less along the upper and lower boundaries. There would be no difference in the equations written at the corners of the flow region. Comparison of results computed by each treatment will be discussed later. A log normal distribution of hydraulic conductivity is adopted here .
Appendix A discusses the background for this assumption. As a consequence of this assumption, the {E }matrix generated for each simula-   Figure B.4 -Weighting Matrix for 7 By 7 Region tion will have a mean value of 0, a standard deviation of 1 and will follow a normal distribution . A log transformation will later be made so that the hydraulic conductivity values will be log normally distributed. Smith (1978) demonstrated that if the { E} matrix is multiplied by a constant term, n, the desired variance in the hydraulic conductivity values results.
He showed that this constant term, n , is a function of the desired variance of the conductivity distribution and the square root of a quantity G defined as: (8) where p (1) and p (1) are the lag 1 correlation coefficients of K in the x y x and y directions respectively. P (2) and P (2) are the lag 2 x y correlation coefficients of K in the x and y directions. p (1, 1) is the lag 1 in the x direction and lag 1 in the y direction correlation coefficients of K.
Because of the assumption of a log normal distribution of hydraulic conductivity, a transformation of the mean and standard deviation is re- Following the multiplication by the filter matrix and addition of the mean of the ln K, a transformation back to K is made. This transformation erodes some of the correlation although this has not been found to be significant (Smith (1978)).
This method does not produce a distribution of hydraulic conductivity which fits exactly the desired correlation structure. This is due to the randomness inserted by the random number generator, the exponential transformation and boundary effects. However each generated distribution is correlated and approximates the desired correlation structure.
When the ensemble of all Monte Carlo simulations is studied, a closer approximation to the desired correlation is obtained.
The computer code developed to perform these operations referred to as Inverse is included in Appendix D.
A step-by-step description of the procedure is given below. It should be noted that this computer code utilizes selected subroutines from the International Mathematics and Statistical Libraries (IMSL).
These subroutines perform necessary tasks for which new subroutines could have been developed. It is considered unnecessary to do that, however, when efficient, documented routines are available.
An example of the use of such a routine may be found in Step 4. The EPSILN matrix is generated by invoking the GGNQF command. This step generates a random number from a normal distribution with a mean of zero and a standard deviation of one. This procedure is done within a DO loop as many times as necessary to assign a value at each node.
A second example of the use of the IMSL library may be found in Step 6.
The subroutine LGINF is called to invert the (I-W) matrix and thereby determine the FILTR Matrix.
The IMSL library has several routines which may be utilized to invert matrices.
LGINF was selected because it is a generalized solution which could easily handle the problem.
Other more efficient and/or accurate inversion routines could have been utilized to reduce required computer time or increase the accuracy of the solution. The reduction of required computer time to invert the matrix was not a major objective of this study although some efficiency might be desirable for subsequent work. The results cited later will demonstrate that the accuracy of the solution is sufficient.
Considering the difficulty at the boundary of maintaining the autocorrelation structure, it was decided to "shave off" the outer boundary nodes all around the block region. To alleviate this limitation, the flow region can be considered a composite of two or more "blocks" that contain statistically homogeneous and internally correlated hydraulic conductivity meshes. The shaved 5 x 5 region may be systematically added to 3 other statistically equivalent regions to develop in a 10 x 20 region. We will see that the computer time to add these regions together is much less than that required for inverting a larger matrix. Indeed there are real limits to the size of the matrix that may realistically be inverted.

B.5 Numerical Solution of Finite Difference Technique
The finite difference solution utilized to solve the flow problem was that computer code developed by Reiter (1981). This code was adapted from other available work, e.g.  for the two dimensional cross-sectional model. The methodology utilized to determine the hydraulic head distribution has been described earlier.
The results of each solution of a stoch-astically generated realization of hydraulic conductivity is a steady state flow field. These results provide all the necessary input data to simulate advective transport.
As cited in Appendix A, the variation in porosity is much less than the variation in hydraulic conductivity. Porosity has been assumed constant throughout the region similar to the assumption in Gelhar, et. al. In these boundary blocks, the gradient is averaged over only one block length.
An example is the upper boundary block where dh/dy = 0 as in Figure B.5, which is simulated by setting H 0 = H 2 . The velocity in the X direction is calculated as before but in the Y direction, the equation reduces to the following: (11) Similar equations are written for other boundary blocks.
The advective transport SX and SY, is then calculated for a given time Microscopic dispersion is included in the model in a manner first described by Ahlstrom, et al (1977) and then modified in part by Smith and Schwartz (1980). In this method, dispersion is treated as a stochastic phenomenon. The water particle is assumed to undergo a random displacement in the direction of flow and another random displacement perpendicular to flow. These displacements are calculated from the following two equations from Ahlstrom, et al (1977). Since the advective transport direction will not always coincide with the X and Y axes the displacements in the X and Y directions are obtained by a transformation, ie, THETA = ARCTAN (SY/SX) (19) The transformation is computed such that DSX and DSY are the dispersive displacements in the X and Y directions, respectively. Figure The computer code simulates water movement from the upgradient side of the block or region to the downgradient side. A check is made to determine if the particle is displaced vertically out of the block as well.
If not, then the next time interval is begun.
However, if vertical movement has occurred into the next block, a check is made for block skipping as before. If vertical movement results in contact with a boundary, the particle is reflected back into the flow region. No additional adjustment in displacements is made in the event of horizontal and vertical block changes, a new time interval is simply begun.
In the event that the particle was displaced vertically out of the block System (SAS) to assist in the analysis of the data from this study.
A listing of the computer code referred to as STKBLKNP, which was adapted from that developed by Reiter (1981) is included in Appendix D.
The code may be subdivided into 11 steps as listed below with the line numbers referring to the listing.  ln T = 0.01075 B + 0.23 (22) where T is the time to complete the inverse and B is the number of nodes in the mesh . The practical limit on region size based on this time requirement will depend on the computer system available; however, the largest region considered in this study was a 12 x 22 region.
where Q is the flow quantity, ~ is the total head loss within the region, L is the total length over which that head loss occurs, H is the The effective hydraulic conductivity on the basis of travel time, Kt has been calculated by setting v equal to v in the above two equations, introducing the height of the region, H, by inserting in both the numerator and denominator and solving for Kt resulting in the following: (26) = where LT is the total length of particle travel, n is the porosity, t is I the mean travel time and L, H and ~ are as defined above. St is a I shape factor computed in the same general manner as S above except that mean travel time rather than flow rate is used. St is a modified shape factor which may be compared to the shape factors from a flow net, the method of fragments or to S calculated above.
Over thirty variations of flow region size and characteristics were studied to determine the relationship between these calculated means.      Values listed for K and K have been adjusted by use of the shape factor S . The results of this study support the technique suggested by Warren & Price (1961) for calculation of a shape factor. For the cases considered, the computed values of effective hydraulic conductivity generally fall between the geometric and harmonic means and is closest to the geometric mean as presented in Figure B12.  Irrespective of this finding, there was still concern over the autocorrelation of hydraulic conductivity in the boundary blocks. This ..   Initially consider the mean and standard deviation of the hydraulic conductivity formulations.    The macroscopic or field scale lateral dispersion coefficient, DT' has been shown by Fried (1975) to be a function of the scatter about the mean ending location after traversing a region .
following equation: Fried presented the (31) where a is the standard deviation of the normalized location pery pendicular to the direction of flow and U and x are as defined earlier .
A record was maintained of the starting and ending y coordinate of each water particle.
A normalized ending y coordinate was then computed along with a cumulative mean and cumulative mean of the standard deviation as each new record of water transport was added. The same was computed for the y coordinate at the estimated time to transport half the region. A typical plot of these results is presented in Figure   B.22. These results also indicate that 400 particles are sufficient to reach a consistent mean value for these variables.

B. 8.5 Normality Check on Travel Times
The ensemble of travel times was utilized to plot a breakthrough curve.     increases. This is probably due to the larger numbers in the calculation however, since when divided by the mean value there is no general trend.  The procedure to develop the site specific soil parameters is essentially unchanged. Field programs are utilized to measure insitu soil parameters and to secure samples for laboratory testing The field and laboratory tests will be analyzed so that a probability density function of the hydraulic conductivity may be formulated. Other parameters such as porosity, and laboratory dispersivi ty characteristics must also be developed.
These  (32) where t&/ 2 is the value of the cumulative distribution function of the I student t distribution at the a/2 significance level, s is the value of the standard deviation and N is the sample size.
The mean value would then be expected to lie between the limits of the sample mean (4.57) plus or minus the error (0.07) i.e., 4.50 and 4.64.
The effective hydraulic conductivity, K , was calculated by dividing the q flow by the mean hydraulic gradient and area. The flow rate may then be estimated with confidence limits of 95% by multiplying by these constants.
In so doing the flow limits of 38 .1 to 39. 9 3 ft /day were estimated with 95% confidence limits.
The water particles are expected to follow a normal distribution so that a hypothetical breakthrough curve may be developed for the extremes of these parameters. The result of such an exercise is presented in Figure   B.34. This allows the engineer to estimate not only the mean arrival time with confidence but also the distribution of travel times.

B.8.9 Macroscopic Dispersivity
The results of the various flow situations was analyzed for characterizations of macroscopic dispersion coefficients and dispersivity.