An iterative cokriging-like technique for ground-water flow modeling

An iterative cokriging-like technique for ground-water flow modeling

Yeh, T-C Jim


Geologic formations are inherently heterogeneous. Importance of the heterogeneity in the transport of contaminants in porous media has been the focus of research during the last two decades. Yeh (1992) provided an overview of several stochastic approaches developed in the last few years for modeling solute transport in heterogeneous aquifers and classified them into two categories: homogeneous (effective parameters) and heterogeneous approaches.

The effective parameter approach assumes that the heterogeneous geologic formation can be homogenized to obtain effective parameters with which one can predict the ensemble behavior of the flow and transport processes. Examples of such studies include those by Gelhar and Axness (1983), and Dagan (1988) for saturated porous media, and Yeh et al. (1985a, b, and c) and Mantoglou and Gelhar (1987a, b, and c), among others for unsaturated media. Although these studies have contributed to a better understanding of flow and transport in heterogeneous aquifers, the major drawback of this approach is that it predicts only the ensemble behavior of the aquifer, which can be quite different from that of the particular realization we encounter in reality.

The heterogeneous approach is designed to consider the nature of spatial variability of hydrologic properties of the aquifer with a limited amount of data. Methods in this approach generally consist of geostatistics, Monte Carlo simulation, and conditional simulation.

Geostatistics is a mathematical interpolation and extrapolation tool which uses the spatial statistics of the data set (e.g., aquifer hydraulic property) to estimate the property at locations where samples are not available. It consists of estimation of variograms (spatial statistics) and kriging. Kriging is the best linear unbiased predictor which provides the estimates of the property at unsampled locations but retains the sample values at locations where the values are known. That is, kriging is a special type of conditional expectation. In the case of simulation of flow in large-scale aquifers where only a limited amount of hydraulic conductivity data are collected, geostatistics is generally found to be useful (Gutjahr, 1981; and Clifton and Neuman, 1982).

Based on the same principle, cokriging takes advantage of the spatial cross correlation of two hydrologic variables, in addition to the spatial correlation of the variable itself, to predict the value of one or both variables at unsampled locations. A typical example in ground-water hydrology is the use of cokriging, based on some head and transmissivity values at some sample locations, to estimate the unknown transmissivity and/or hydraulic head values at other locations (e.g., Ahmed and Marsily, 1993; Kitanidis and Vomvoris, 1983; Hoeksema and Kitanidis, 1984; Gutjahr and Wilson, 1989; and Hoeksema and Kitanidis, 1989). The resultant hydraulic head and transmissivity fields are then used to simulate the transport of solutes. Since the conductivity and head values at the locations where samples are not available are estimated by smoothed values analogous to a conditional expectation, the effect of variability around the expected value on solute transport is neglected. Therefore, the uncertainty associated with the predicted concentration field cannot be addressed. However, by eliminating the determination of the uncertainty, which requires intensive computational efforts, the cokriging approach is computationally economic and is often considered as a more practical and more realistic approach than the effective parameter approach.

Although hydraulic head and transmissivity fields derived from cokriging have been found reasonable, there is no guarantee that the hydraulic head and transmissivity estimates satisfy the principle of conservation of mass because of the linear assumption used in the estimation procedure. Therefore, these conductivity and head estimates may produce erroneous velocity fields and, in turn, an erroneous concentration distribution.

The objective of the paper is to present our newly developed technique (iterative cokriging-like method) which combines the use of classical cokriging, a numerical model for flow, and an iterative scheme to improve our estimates of transmissivity, hydraulic head, and velocity fields. Several examples are used to demonstrate the robustness of our new iterative approach.


Classical Cokriging

Detailed discussion of the theory of cokriging can be found in many literatures (e.g., Marsily, 1986). Here we will present a brief description. Consider that Z sub 1 (x) and Z sub 2 (x) are two correlated second-order stationary random fields (i.e., transmissivity and hydraulic head, respectively) with known covariances and means. The estimation of Z sub 1 (and Z sub 2 , if necessary) using cokriging is again based on the best, linear, and unbiased estimator technique. That is,

[Equation (1) omitted]

where Z sub 1 *(x sub 0 ) is the estimate of Z sub 1 at location x sub 0 , Z sub 1 (x sub j ) are the measured values of Z sub 1 at locations x sub 1 , x sub 2 ,… x sub n . Similarly, Z sub 2 (x sub n+k ) are the measured values of Z sub 2 at locations x sub n+1 , x sub n+2 ,… x sub n+m . Note that Z sub 1 and Z sub 2 need not be measured at the same location. lambda sub 1j and lambda sub 2k are cokriging weights corresponding to Z sub 1 and Z sub 2 , respectively.

The minimal variance criterion requires that

E[(Z sub 1 *(x sub 0 ) – Z sub 1 (x sub 0 ) sup 2 ] = minimal (2)

Similar to the development of kriging equation, one can obtain the cokriging system of equations for second-order stationary processes (Z sub 1 and Z sub 2 ) as follows:

[Equations (3)&(4) omitted]

for l = 1,…, m. The autocovariances of Z sub 1 and Z sub 1 and cross covariances between Z sub 1 and Z sub 2 are denoted by C sub 11 and C sub 22 , and C sub 12 and C sub 21 , respectively. By solving the system of equations (3) and (4), the cokriging weights can be obtained. Then, equation (1) can be utilized to obtain the estimates of either or both variables at any location.

Prior to the solution of the system of equations, the autocovariances and cross covariances as functions of the separation distance have to be known. This information can be obtained from either field data or some theoretical analyses. In this study, the autocovariances and cross covariances are derived from a perturbation-spectral analysis of steady ground-water flow in unbounded, two-dimensional, and second-order stationary random transmissivity field (Mizell et al., 1982).

If we assume that lnT (natural log of transmissivity, T) and hydraulic head, h, are second-order stationary stochastic processes and assume the mean flow is in x sub 1 direction in a two-dimensional aquifer, the spectrum of h under steady flow in unbounded porous media can be derived to be (Mizell et al., 1982):

[Equation (5) omitted]

where S sub h is the head spectrum, and S sub t is the lnT spectrum. Similarly, the cross spectrum of lnT and h can be derived as

[Equation (6) omitted]

where S sub fh is the cross spectrum of lnT and h, J sub 1 is the mean gradient in x sub 1 direction, and k sub 1 and k sub 2 are the wave numbers corresponding to x sub 1 and x sub 2 . Using (5) and (6), the spectrum of h and the cross spectrum between lnT and h can be obtained if the spectrum, S sub f , is given. It should be pointed out that equations (5) and (6) are first-order approximations which will deviate from the true solution when the variance of lnT or h becomes large.

In this paper a modified Whittle spectrum function (Mizell et al., 1982) was used since the spectrum corresponding to an exponential covariance is not stochastically differentiable. The modified Whittle spectrum function takes the form

[Equation (7) omitted]

The corresponding covariance, the inverse Fourier transform of (7), can be expressed as

[Equation (8) omitted]

where alpha = pi/4 lambda; C sub f (zeta) is the covariance function for lnT; zeta denotes the magnitude of the separation distance; [character omitted] represents the variance of transmissivity; the correlation scale is denoted by lambda and K sub 1 and K sub 0 are the modified Bessel function of second kind, first order and zeroth order, respectively. Similarly, the covariance of h and the cross covariance between lnT and h can be obtained analytically by taking the inverse Fourier transform of (7) and (8), respectively. The program for the classical cokriging used in this study is a modified version of a Co-simulator of conductivity and head using a Fast Fourier Transform (COFFT) developed by Gutjahr et al. (1992).

Iterative Cokriging-Like Method

The head and transmissivity fields obtained using classical cokriging are, in general, satisfactory for cases where the perturbations of h and lnT are small. However, these head and transmissivity fields may produce an erroneous velocity field due to the nature of the first-order approximation and the linear predictor embedded in classical cokriging. Similar problems were reported in the conditional simulation of flow and solute transport in heterogeneous porous media (Harter et al., 1992). In addition, they do not consider any boundary effects.

To overcome these difficulties, an iterative cokriging-like method is developed, which combines classical cokriging, a numerical flow model, and an iterative scheme. A schematic flow chart, illustrating the algorithm of the iterative approach is given in Figure 1.[Figure 1 omitted] The first step of this iterative cokriging-like approach is to generate transmissivity and head fields using the classical cokriging technique with some observed transmissivity data [i.e., T*(x sub i ) at sample locations, x sub i , where i = 1, n] and observed head values [h*(x sub j+n ) at sample locations, x sub j+n , where j = 1, 2,…, m]. Then, as the second step, the cokriged T field is input to a two-dimensional numerical model (MMOC2, Yeh et al., 1993) for steady-state flow with the observed head values at the sample location as internal constant-head boundary conditions. Because of the nonlinear nature of the flow system, the cokriged head field is treated as a guess solution to the numerical problem, which allows a rapid convergence of the numerical solution if an iterative matrix solver is used. Application of this approach to unsaturated flow problems shows extremely robust results in terms of CPU time saving (Harter and Yeh, 1993). Note that the head field obtained from this numerical solution satisfies specified boundary conditions and the principle of conservation of mass. Therefore, the velocity field is well-behaved.

Although the head field, h sub new (x), obtained from the numerical model preserves the measured head values at sample locations, the influence of the measured head values on the estimation of the adjacent transmissivity values is not considered. To incorporate this influence or update the transmissivity field, at the third step of our iterative approach, classical cokriging is utilized again to generate new transmissivity fields, T sub new (x) based on the measured transmissivity data [i.e., T*(x sub i )] and all the head values, h sub new (x), in the flow domain as given conditions. This newly cokriged transmissivity field, T sub new (x), is then input to the numerical model to solve for a new head field. Such an iterative process continues until the maximum change in the head field is less than a specified tolerance, epsilon. Through the iteration, the iterative approach incorporates the measured head and transmissivity values and produces transmissivity and hydraulic head fields that are consistent with the governing flow equation. However, one should recognize the fact that this approach may not be consistent with the statistical conditional expectation, and we will call our approach an iterative cokriging-like approach.

Results and Discussion

To demonstrate the performance of our iterative cokriging-like approach, we generated several hypothetical two-dimensional heterogeneous aquifers with different degrees of heterogeneity (measured by [character omitted] the variance of lnT) and then used them as the real-world analogues. Here the spectral technique (Gutjahr et al., 1992) was used to generate the random lnT field. Subsequently, head and flow fields in these hypothetical aquifers were obtained by solving the flow equation using MMOC2 with specified head boundary conditions. The resultant head fields were assumed to represent the “real-world” analogues of the hydraulic head fields corresponding to the true transmissivity fields. The hypothetical aquifer is illustrated in Figure 2, which has 16 X 16 square elements with length equal to 0.2. No flow boundaries were imposed on the upper and lower boundaries and constant heads were specified on the left and right boundaries, creating a mean hydraulic head gradient equal to 0.05 in the x sub 1 direction.

Once the hypothetical aquifer was created and the corresponding head field simulated, a random sampling scheme was then employed to determine the sample locations where transmissivity and head data in the hypothetical aquifer were collected. Using these transmissivity and head data sets, classical cokriging and our iterative cokriging-like approaches were utilized to estimate both the transmissivity and hydraulic head fields for the entire aquifer. Subsequently, the velocity fields were calculated by the application of Darcy’s law. Finally, the estimated transmissivity, hydraulic head, and velocity fields for different degrees of variability of lnT were compared to those of the “real-world” analogues.

In this paper, four cases are examined. The first three cases correspond to three levels of degrees of heterogeneities (i.e., [character omitted] equals 0.25, 2.0, and 4.0 for Cases 1, 2, and 3, respectively). Ten observations of lnT and 20 observations of h were used in these three cases. A fourth case (Case 4) was designed to examine the effect of the number of head measurements. Specifications for Case 4 are identical to those of Case 2 with the exception that the number of head observations was increased to 40. For all these four cases, the correlation scale of the lnT field, lambda (lambda = /4 alpha), was assumed to be five element lengths and to be isotropic. The value of epsilon (the tolerance level) was set to 0.00001.

Figures 3a, b, and c show the transmissivity field corresponding to the real-world analogue, the one estimated by classical cokriging approach, and one by the iterative cokriging-like approach, respectively in Case 1.[Figure 3 omitted] The head fields corresponding to these transmissivity fields are plotted in Figures 3d, e, and f. The resultant velocity fields using these transmissivity and head fields are illustrated in Figures 3g, h, and i. Note that the velocity vector is represented by arrows and the solid lines represent the stream lines associated with the velocity field. As illustrated in these figures, the resultant fields using either classical or the iterative cokriging-like technique are very similar, although slight differences exist.

Superiority of the iterative cokriging-like technique becomes apparent when the variance of lnT increases. Figures 4a, b, and c show the transmissivity field corresponding to the Case 2 real-world analogue, the one estimated by classical cokriging, and one by the iterative cokriging-like method, respectively.[Figure 4 omitted] The corresponding head and resultant velocity fields are illustrated in Figures 4d, e, and f and Figures 4g, h, and i, respectively. As shown in Figures 4g, h, and i, the estimated transmissivity and hydraulic head fields, derived from the iterative cokriging-like approach, produce a velocity field that is much closer to the true one than that based on classical cokriging approach.

The robustness of the iterative cokriging-like approach is even more evident in Figure 5 for Case 3 where the variance of lnT is increased to 4.0.[Figure 5 omitted] Figure 5h shows that classical cokriging approach produces an unacceptable velocity field and thus unacceptable streamlines. This unacceptable velocity field can be attributed to the head covariance and cross covariance functions between lnT and h derived from the perturbation-spectral approach which relies on the assumption of small perturbations. In addition, the limitation of the linear predictor approach used by classical cokriging may augment such a problem. However, as illustrated in Figure 5i, the resultant velocity field calculated using the transmissivity and hydraulic head fields estimated from our iterative cokriging-like approach is satisfactory in terms of mass conservation and accuracy as compared to the true velocity field (Figure 5g). This is also true for the head field. Figures 5d, e, and f show that the head field estimated by the iterative cokriging-like method is much closer to the true one, compared to that by the classic approach, even though the same amount of measured head and transmissivity data were used. Comparisons of the transmissivity fields estimated by classical cokriging (Figure 5b) and that by the iterative cokriging-like method (Figure 5c) also reveal some improvements of transmissivity estimates. These improvements are essentially a result of the iterative updating of the transmissivity field using the improved head field from the numerical model.

The result for Case 4 is depicted in Figures 6a through 6i.[Figure 6 omitted] In this case, the increase of the number of head observations (from 20 to 40) improves the estimation of the head and velocity based on the iterative cokriging-like method. The improvement of the lnT fields is not as significant as that of the head and velocity fields (for comparison, see Figure 4).

Some important aspects of the effect of heterogeneities on ground-water flow are unveiled by the results of this study. As expected, cokriging produces smooth transmissivity fields, omitting small-scale or local heterogeneities, as compared to the true transmissivity fields. However, the estimated velocity field based on the smoothed lnT field, regardless of the magnitude of the variance of lnT, closely resembles the true velocity field, implying the effect of local heterogeneities is negligible in terms of general flow pattern. This phenomenon can be attributed to the fact that local scale heterogeneity is responsible for local mixing effect (or dispersion) and thus the distribution of travel times. On the other hand, the cokriged lnT representing the mean values controls the general flaw pattern. This finding seems in good agreement with our findings in field tracer experiments in Georgetown, South Carolina (Yeh et al., 1992), which state that the bulk behavior of the observed three-dimensional plume appears to be governed by some “important” heterogeneity (i.e., stratification and some extremely low permeable inclusions). If this is the case, the strategy for tackling field problems should focus on the delineation of large and “important” heterogeneities. Macrodispersion concept (Gelhar and Axness, 1983) may thus be suitable for representing the effect of small-scale heterogeneity in the prediction of solute transport in field-scale aquifers.

Finally, the convergence of the iterative scheme is generally rapid but the rate of convergence depends an the tolerance and the variance of lnT. For all the cases examined, the number of iterations required to reach a tolerance of 0.00001 in head range from 10 to 50 iterations.


A new and physically correct iterative technique is developed, which overcomes the limitations of classical cokriging approach. These limitations stem from the fact that classical cokriging is a linear estimator and the theoretical covariance function of head and cross covariance function between lnT and h are limited to small perturbations. Such limitations restrict the application of classical cokriging to many real field problems where variability of lnT is expected to be large. Our new approach, combining classical cokriging technique, a numerical model, and an iterative scheme relaxes the limitations and makes the geostatistical method more attractive.

Although the head measurement does not improve the estimate of the transmissivity field substantially, a large number of hydraulic head measurements will improve the delineation of flow directions and paths. The general flow pattern appears controlled by the general transmissivity pattern and the small-scale and local variability in transmissivity has little effect on the flow direction.

Finally, we believe the paper presents a practical tool that can be easily adopted by many practitioners to tackle realistic field ground-water flow problems, although more theoretical work is still needed to generalize the approach for different flow scenarios.


This research is funded by the Subsurface Science Program, Environmental Sciences Division, U.S. Department of Energy, under Grant No. DE-FG02-91ER61199.


Ahmed, S. and G. de Marsily. 1993. Cokriged estimation of aquifer transmissivity as an indirect solution of the inverse problem: A practical approach. Water Resour. Res. v. 29, no. 2, pp. 512-530.

Clifton, P. M. and S. P. Neuman. 1982. Effects of kriging and inverse modeling on conditional simulation of the Avra valley aquifer in southern Arizona. Water Resour. Res. v. 18, no. 4, pp. 1215-1234.

Dagan, G. 1988. Time-dependent macrodispersion for solute transport in anisotropic heterogeneous aquifers. Water Resour. Res. v. 24, no. 9, pp. 1491-1500.

Gelhar, L. W. and C. L. Axness. 1983. Three-dimensional stochastic analysis of macrodispersion in aquifers. Water Resour. Res. v. 19, no. 1, pp. 161-181.

Gutjahr, A. L. 1981. Kriging in stochastic hydrology: Assessing the worth of data. AGU Chapman Conference on Spatial Variability in Hydrologic Modeling, Fort Collins, CO.

Gutjahr, A. L. and J. L. Wilson. 1989. Co-kriging for stochastic models. Transport in Porous Media. v. 4, no. 6, pp. 585-598.

Gutjahr, A., Q. Bai, and S. Hatch. 1992. Conditional simulation applied to contaminant flow modeling. Tech. Report. Dept. of Math., New Mexico Inst. of Min. and Tech. Feb.

Harter, T., T.-C.J. Yeh, and A. L. Gutjahr. 1992. Conditional simulation of unsaturated flow-fields for assessing contaminant transport in heterogeneous soils. Eos Trans. v. 73, no. 43. p. 172, AGU Fall Meeting.

Harter, T. and T.-C.J. Yeh. 1993. An efficient method for simulating steady unsaturated flow in random porous media: Using an analytical perturbation solution as initial guess to a numerical model. Water Resour. Res. v. 29, no. 12, pp. 4139-4149.

Hoeksema, R. J. and P. K. Kitanidis. 1984. An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling. Water Resour. Res. v. 20, no. 7, pp. 1003-1020.

Hoeksema, R. J. and P. K. Kitanidis. 1989. Prediction of transmissivities. heads, and seepage velocities using mathematical modeling and geostatistics. Adv. Water Resour. v. 12, pp. 90-101.

Kitanidis, P. K. and E. G. Vomvoris. 1983. A geostatistical approach to the inverse problem in groundwater modeling (steady-state) and one-dimensional simulation. Water Resour. Res. v. 19, pp. 677-690.

Mantoglou, A. and L. W. Gelhar. 1987a. Stochastic modeling of large scale transient unsaturated flow systems. Water Resour. Res. v. 23, no. 1, pp. 37-46.

Mantoglou, A. and L. W. Gelhar. 1987b. Capillary tension head variance, mean soil moisture content, and effective specific soil moisture capacity of transient unsaturated flow in stratified soils. Water Resour. Res. v. 23, no. 1, pp. 47-56.

Mantoglou, A. and L. W. Gelhar. 1987c. Effective hydraulic conductivity of transient unsaturated flow in stratified soils. Water Resour. Res. v. 23. no. 1, pp. 57-68.

Marsily, G. de. 1986. Quantitative Hydrogeology. Academic Press, San Diego, CA. 440 pp.

Mizell, S. A., A. L. Gutjahr, and L. W. Gelhar. 1982. Stochastic analysis of spatial variability in two-dimensional steady groundwater flow assuming stationary and nonstationary heads. Water Resour. Res. v. 18, no. 4, pp. 1053-1067.

Yeh, T.-C.J., L. W. Gelhar, and A. L. Gutjahr. 1985a. Stochastic analysis of unsaturated flow in heterogeneous soils: Part 1, Statistically isotropic media. Water Resour. Res. v. 21, no. 4, pp. 447-456.

Yeh, T.-C.J., L. W. Gelhar, and A. L. Gutjahr. 1985b. Stochastic analysis of unsaturated flow in heterogeneous soils: Part 2, Statistically isotropic media. Water Resour. Res. v. 21, no. 4, pp. 457-464.

Yeh, T.-C.J., L. W. Gelhar, and A. L. Gutjahr. 1985c. Stochastic analysis of unsaturated flow in heterogeneous soils: Part 3, Observations and applications. Water Resour. Res. v. 21, no. 4, pp. 465-471.

Yeh, T.-C.J., J. Mas-Pla, T. M. Williams, and J. F. McCarthy. 1992. Three-dimensional simulation of a chloride plume in a coastal sandy aquifer, Georgetown site. SC. EOS Trans. v. 73, no. 43.

Yeh, T.-C.J. 1992. Stochastic modeling of groundwater flow and solute transport in aquifers. Hydrological Processes. v. 6, pp. 369-395.

Yeh, T.-C.J., R. Srivastava, A. Guzman, and T. Harter. 1993. A numerical model for water flow and chemical transport in variably saturated porous media. Ground Water. v. 31, no. 4, pp. 634-644.

Copyright Ground Water Publishing Company Jan 1995

Provided by ProQuest Information and Learning Company. All rights Reserved.