Peter K. Kitanidis

Quasi-Linear Geostatistical Theory for Inversing

A quasi-linear theory is presented for the geostatistical solution to the inverse problem. The archetypal problem is to estimate the log transmissivity function from observations of head and log transmissivity at selected locations. The unknown is parameterized as a realization of a random field, and the estimation problem is solved in two phases: structural analysis, where the random field is characterized, followed by estimation of the log transmissivity conditional on all observations. The proposed method generalizes the linear approach of Kitanidis and Vomvoris (1983). The generalized method is superior to the linear method in cases of large contrast in formation properties but informative measurements, i.e., there are enough observations that the variance of estimation error of the log transmissivity is small. The methodology deals rigorously with unknown drift coefficients and yields estimates of covariance parameters that are unbiased and grid independent. The applicability of the methodology is demonstrated through an example that includes structural analysis, determination of best estimates, and conditional simulations.

The concept of the Dilution Index

In many applications, it is important to make the distinction between spreading and dilution of a plume in groundwater. Spreading is associated with the stretching and deformation of a contaminant plume, whereas dilution is associated with the increase in volume of the fluid occupied by the solute. The dilution and spreading of a Gaussian plume in a homogeneous porous medium with constant velocity are related in a simple fashion and are both characterized by the same parameters, the dispersion coefficients. However, the geological formations of interest in field applications are heterogeneous, and the plumes are irregular in shape. The dispersion coefficients that are deduced from tracer tests usually measure an overall rate at which a tracer plume spreads about its centroid and depend critically on the heterogeneity of the formation. These macroscopic dispersion coefficients are not reliable measures of the rate at which the maximum concentration is reduced because in heterogeneous formations the rates of dilution and spreading can be quite different. The main objective of this work is to introduce a new macroscopic measure of dilution, the dilution index E. Examples serve to demonstrate the usefulness of the measure. A general expression for the rate of dilution of a tracer plume is derived. The exact rate of increase of the dilution index under the idealized conditions of constant dispersion coefficients and a Gaussian plume is computed, and a lower bound is found to the same quantity for non-Gaussian plumes. For the general heterogeneous case the analysis demonstrates that the instantaneous rate of increase of ln E is proportional to the small-scale dispersion coefficients, everything else being the same. The rate of increase of ln E depends also on the degree of irregularity in the shape of the plume. Thus, in the long term, geologic heterogeneity should increase the rate of dilution because spatial variability in the flow velocity tends to deform plumes and make them less regular.

A comparison of seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow

This paper describes the first major attempt to compare seven different inverse approaches for identifying aquifer transmissivity. The ultimate objective was to determine which of several geostatistical inverse techniques is better suited for making probabilistic forecasts of the potential transport of solutes in an aquifer where spatial variability and uncertainty in hydrogeologic properties are significant. Seven geostatistical methods (fast Fourier transform (FF), fractal simulation (FS), linearized cokriging (LC), linearized semianalytical )LS), maximum likelihood (ML), pilot point (PP), and sequential self-calibration (SS)) were compared on four synthetic data sets. Each data set had specific features meeting (or not) classical assumptions about stationarity, amenability to a geostatistical description, etc. The comparison of the outcome of the methods is based on the prediction of travel times and travel paths taken by conservative solutes migrating in the aquifer for a distance of 5 km. Four of the methods, LS, ML, PP, and SS, were identified as being approximately equivalent for the specific problems considered. The magnitude of the variance of the transmissivity fields, which went as high as 10 times the generally accepted range for linearized approaches, was not a problem for the linearized methods when applied to stationary fields; that is, their inverse solutions and travel time predictions were as accurate as those of the nonlinear methods. Nonstationarity of the “true” transmissivity field, or the presence of “anomalies” such as high-permeability fracture zones was, however, more of a problem for the linearized methods. The importance of the proper selection of the semivariogram of the log10 (T) field (or the ability of the method to optimize this variogram iteratively) was found to have a significant impact on the accuracy and precision of the travel time predictions. Use of additional transient information from pumping tests did not result in major changes in the outcome. While the methods differ in their underlying theory, and the codes developed to implement the theories were limited to varying degrees, the most important factor for achieving a successful solution was the time and experience devoted by the user of the method.

Unbiased minimum-variance linear state estimation

A method is developed for linear estimation in the presence of unknown or highly non-Gaussian system inputs. The state update is determined so that it is unaffected by the unknown inputs. The filter may not be globally optimum in the mean square error sense. However, it performs well when the unknown inputs take extreme or unexpected values. In many geophysical and environmental applications, it is performance during these periods which counts the most. The application of the filter is illustrated in the real-time estimation of mean areal precipitation.

Prediction by the method of moments of transport in a heterogeneous formation

Abstract Natural geologic formations are highly heterogeneous. It is impossible and often unnecessary to describe in deterministic terms the spatial variability of their properties. However, the hydrogeologic parameters may be represented in probabilistic terms. Prediction of solute transport may then be defined as the derivation of the probabilistic properties of concentration. This work deals with the first two integral (or spatial) moments of solute concentration in a heterogeneous formation of infinite extent. The first moment is the vector of the mean position of the centroid of the plume and, in a generalized sense, represents advection. The second moment is the matrix of dyadics of the mean squared displacement about the average position of the centroid of the plume and, again in a generalized sense, represents dispersion. Assuming that the mixing at the laboratory scale is Fickian, with random but time-invariant velocities and lab-scale dispersion matrices, the differential equations satisfied by the first two moments are derived. An analytical first-order (or small perturbation) solution is obtained for stationary velocity, and compared with a numerical solution.

Determination of the effective hydraulic conductivity for heterogeneous porous media using a numerical spectral approach: 1. Method

Using the numerical spectral method described in the companion paper (Dykaar and Kitanidis, this issue) the effective conductivity of a three-dimensional, Isotropic, stationary, lognormally distributed hydraulic conductivity is computed. Six cases were investigated for variances in log conductivity ranging between 1 and 6. The results show that averaging volumes of at least 30 integral scales are required before the effective conductivity reaches its asymptotic value. The results support Matherons (1967) hypothesis that Ke = KG exp (σY2/6), where Ke is the effective hydraulic conductivity, and KG and σY2 are the geometric mean and variance, respectively, of the hydraulic conductivity field. We also find that for a Gaussian covariance function, heterogeneities smaller than about 1.3 integral scales do not significantly contribute to the effective conductivity. In two dimensions, averaging volumes of more than 80 integral scales are required before the effective conductivity reaches the analytic infinite-domain result of the geometric mean, while the effective conductivity is insensitive to heterogeneities less than 2 integral scales in size. The method is also applied to data from a shale and sandstone formation, and the results are compared with those from other methods.

Characterization of mixing and dilution in heterogeneous aquifers by means of local temporal moments

Breakthrough curves of a conservative tracer in a heterogeneous two-dimensional aquifer are analyzed by means of their temporal moments. The average velocity and the longitudinal macrodispersion coefficient of the equivalent one-dimensional aquifer obtained through cross-sectional averaging of concentration can be defined from the first and second central moments of a breakthrough curve integrated over the outflow boundary of the domain. On the basis of an integrated breakthrough curve, one cannot distinguish between actual solute dilution, which involves concentration reduction, and variability of arrival times among parts of the plume at different cross-sectional positions. Analyzing the temporal moments of breakthrough curves at a “point” within the domain gives additional information about the dilution of the tracer. From these local first and second central moments an apparent seepage velocity υa and an apparent dispersivity of mixing αa can be derived. For short travel distances, αa equals the local-scale longitudinal dispersivity. It increases with the travel distance but much more slowly than the macrodispersivity. At the large-distance limit, αa may eventually reach the level of macrodispersivity. Lenses of high conductivity where groundwater flow converges are identified as regions of preferential enhanced mixing. The spatial distribution of these regions causes a high degree of variability of αa within a domain, indicating a high degree of uncertainty in the quantification of dilution at early stages. In an accompanying paper [Cirpka and Kitanidis, this issue] the results of conservative tracer transport are utilized for the study of mixing-controlled reactive transport.

Maximum likelihood parameter estimation of hydrologic spatial processes by the Gauss-Newton method

Abstract The Gauss-Newton method for maximum likelihood (ML) parameter estimation of spatially correlated hydrologic fields is described with emphasis on computational aspects associated with its implementation. The essence of the Gauss-Newton method is that, for purposes of applying Newtons optimization method, the Hessian of the function to be minimized is replaced by its average value eliminating the need to calculate second derivatives. In ML estimation, the parameters are determined through minimization of the negative log-likelihood function and the average of the Hessian is known as the information matrix. The information matrix provides a lower bound to the covariance matrix of parameter estimation error. Methodologies to account for near-singular information matrix and to enforce lower or upper bounds in the parameters are described. A general method to eliminate the bias in covariance parameter estimates caused by unknown drift parameters is also presented. Some simple examples are presented to illustrate the potential advantages of Gauss-Newton. The conditions under which this method should perform well are also discussed.

A geostatistical approach to contaminant source identification

A geostatistical approach to contaminant source estimation is presented. The problem is to estimate the release history of a conservative solute given point concentration measurements at some time after the release. A Bayesian framework is followed to derive the best estimate and to quantify the estimation error. The relation between this approach and common regularization and interpolation schemes is discussed. The performance of the method is demonstrated for transport in a simple one-dimensional homogeneous medium, although the approach is directly applicable to transport in two- or three-dimensional domains. The methodology produces a best estimate of the release history and a confidence interval. Conditional realizations of the release history are generated that are useful in visualization and risk assessment. The performance of the method with sparse data and large measurement error is examined. Emphasis is placed on formulating the estimation method in a computationally efficient manner. The method does not require the inversion of matrices whose size depends on the grid size used to resolve the solute release history. The issue of model validation is addressed.

Concentration fluctuations and dilution in aquifers

The concentration of solute undergoing advection and local dispersion in a random hydraulic conductivity field is analyzed to quantify its variability and dilution. Detailed numerical evaluations of the concentration variance sc are compared to an approximate analytical description, which is based on a characteristic variance residence time (VRT), over which local dispersion destroys concentration fluctuations, and effective dispersion coefficients that quantify solute spreading rates. Key features of the analytical description for a finite size impulse input of solute are (1) initially, the concentration fields become more irregular with time, i.e., coefficient of variation, CV 5 sc/^c&, increases with time (^c& being the mean concentration); (2) owing to the action of local dispersion, at large times (t . VRT), s c is a linear combination of ^c& 2 and (›^c&/› xi) 2 , and the CV decreases with time (at the center, CV > (N) 1/2 VRT/t, N being the macroscopic dimensionality of the plume); (3) at early time, dilution and spreading can be severely disconnected; however, at large time the volume occupied by solute approaches that apparent from its spatial second moments; and (4) in contrast to the advection-local dispersion case, under advection alone, the CV grows unboundedly with time (at the center, CV } t N/4 ), and spatial second moment is increasingly disconnected from dilution, as time progresses. The predicted large time evolution of dilution and concentration fluctuation measures is observed in the numerical simulations.