Shlomo P. Neuman

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.

Prediction of steady state flow in nonuniform geologic media by conditional moments: exact nonlocal formalism, effective conductivities, and weak approximation

We consider the effect of measuring randomly varying local hydraulic conductivities K(x) on ones ability to predict steady state flow within a bounded domain, driven by random source and boundary functions. More precisely, we consider the prediction of local hydraulic head h(x) and Darcy flux q(x) by means of their unbiased ensemble moments 〈h(x)〉κ and 〈q(x)〉κ conditioned on measurements of K(x). These predictors satisfy a deterministic flow equation in which 〈q(x)〉κ = −κ(x)∇〈h(x)〉κ + rκ(x), where κ(x) is a relatively smooth unbiased estimate of K(x) and rκ(x) is a “residual flux.” We derive a compact integral expression for rκ(x) which is rigorously valid for a broad class of K(x) fields, including fractals. It demonstrates that 〈q(x)〉κ is nonlocal and non-Darcian so that an effective hydraulic conductivity does not generally exist. We show analytically that under uniform mean flow the effective conductivity may be a scalar, a symmetric or a nonsymmetric tensor, or a set of directional scalars which do not form a tensor. We demonstrate numerically that in two-dimensional mean radial flow it may increase from the harmonic mean of K(x) near interior and boundary sources to the geometric mean far from such sources. For cases where rκ(x) can neither be expressed nor approximated by a local expression, we propose a weak (integral) approximation (closure) which appears to work well in media with pronounced heterogeneity and improves as the quantity and quality of K(x) measurements increase. The nonlocal deterministic flow equation can be solved numerically by standard methods; our theory shows clearly how the scale of grid discretization should relate to the scale, quantity, and quality of available data. After providing explicit approximations for the second moments of head and flux prediction errors, we conclude by discussing practical methods to compute κ(x) from noisy measurements of K(x) and to calculate required second moments of the associated estimation errors when K(x) is lognormal.

Eulerian‐Lagrangian Theory of transport in space‐time nonstationary velocity fields: Exact nonlocal formalism by conditional moments and weak approximation

A unified Eulerian-Lagrangian theory is presented for the transport of a conservative solute in a random velocity field that satisfies locally ∇ · v(x, t) = f(x, t), where f(x, t) is a random function including sources and/or the time derivative of head. Solute concentration satisfies locally the Eulerian equation ∂c(x, t)/∂t + ∇ · J(x, t) = g(x, t), where J(x, t) is advective solute flux and g(x, t) is a random source independent of f(x, t). We consider the prediction of c(x, t) and J(x, t) by means of their unbiased ensemble moments 〈c(x, t)〉ν and 〈J(x, t)〉ν conditioned (as implied by the subscript) on local hydraulic measurements through the use of the latter in obtaining a relatively smooth unbiased estimate ν(x, t) of v(x, t). These predictors satisfy ∂〈c(x, t)〉v/∂t + ∇ · 〈J(x, t)〉ν = 〈g(x, t)〉ν, where 〈J(x, t)〉ν = ν(x, t)〈c(x, t)〉ν + Qν(x, t) and Qν(x, t) is a dispersive flux. We show that Qν, is given exactly by three space-time convolution integrals of conditional Lagrangian kernels αν with ∇·Qν, βν with ∇〈c〉ν, and γν with 〈c〉ν for a broad class of v(x, t) fields, including fractals. This implies that Qν(x, t) is generally nonlocal and non-Fickian, rendering 〈c(x, t)〉ν non-Gaussian. The direct contribution of random variations in f to Qν depends on 〈c〉ν rather than on ∇〈c〉ν,. We elucidate the nature of the above kernels; discuss conditions under which the convolution of βν and ∇〈c〉 becomes pseudo-Fickian, with a Lagrangian dispersion tensor similar to that derived in 1921 by Taylor; recall a 1952 result by Batchelor which yields an exact expression for 〈c〉ν at early time; use the latter to conclude that linearizations which predict that 〈c〉ν bifurcates at early time when the probability density function of v is unimodal cannot be correct; propose instead a weak approximation which leads to a nonlinear integro-differential equation for 〈c〉ν due to an instantaneous point source and which improves with the quantity and quality of hydraulic data; demonstrate that the weak approximation is analogous to the “direct interaction” closure of turbulence theory; offer non-Fickian and pseudo-Fickian weak approximations for the second conditional moment of the concentration prediction error; demonstrate that it yields the so-called “two-particle covariance” as a special case; conclude that the (conditional) variance of c does not become infinite merely as a consequence of disregarding local dispersion; and discuss how to estimate explicitly the cumulative release of a contaminant across a “compliance surface” together with the associated estimation error.

Generalized scaling of permeabilities: Validation and effect of support scale

The permeabilities and dispersivities of geologic media are known to vary with the scale of observation. Particularly well documented is the consistent increase in apparent longitudinal dispersivity with the mean travel distance of a tracer. This has been previously interpreted by the author to imply that the permeabilities of many geologic media scale, on the average, according to the power-law semivariogram γ(s)=c ν s where c is a constant and s is distance. Tracer test data support this conclusion indirectly at least over scales from 10 cm to 3,500 m. The present paper cites evidence for such behavior over scales from 10 cm to 45 km based directly on permeability and transmissivity data

Subsurface flow and transport : a stochastic approach

Part I. Introduction: 1. Stochastic modeling of flow and transport: the broad perspective Gedeod Dagan Part II. Subsurface Characterization and Parameter Estimation: 2. Characterization of geological heterogeneity Mary P. Anderson 3. Application of geostatistics in subsurface hydrology Javier Samper 4. Formulations and computational issues of the inversion of random fields Jesus Carrera, Agustin Medina, Carl Axness and Tony Zimmerman Part III. Flow Modeling and Aquifer Management: 5. Groundwater in heterogeneous formations Peter Kitanidis 6. Aspects of numerical methods in multiphase flows Richard E. Ewing 7. Incorporating uncertainty into aquifer management models Steve Gorelick Part IV. Transport in Heterogeneous Aquifers: 8. Transport of inert solutes by groundwater: recent developments and current issues Yoram Rubin 9. Transport of reactive solutes Vladimir Cvetkovic 10. Nonlocal reactive transport with physical and chemical heterogeneity: linear nonequilibrium sorption with random rate coefficients Bill X. Hu, Fei-Wen Deng and John Cushman 11. Perspectives on field scale application of stochastic subsurface hydrology Lynn W. Gelhar Part V. Fractured Rocks and Unsaturated Soils: 12. Component characterization: an approach to fracture hydrology Jane C. S. Long, Christine Doughty, Akhil Datta-Gupta, Kevin Hestir and Don Vasco 13. Stochastic analysis of solute transport in partially saturated heterogeneous soils David Russo 14. Field scale modeling of immiscible organic chemical spills Jack Parker Part VI. A View to the Future: 15. Stochastic approach to subsurface flow and transport: a view to the future Shlomo Neuman.

Nonlocal and localized analyses of conditional mean steady state flow in bounded, randomly nonuniform domains: 1. Theory and computational approach

We consider the effect of measuring randomly varying hydraulic conductivitiesK(x) on ones ability to predict numerically, without resorting to either Monte Carlo simulation or upscaling, steady state flow in bounded domains driven by random source and boundary terms. Our aim is to allow optimum unbiased prediction of hydraulic heads h(x) and fluxes q(x) by means of their ensemble moments, 〈h(x)〉c and 〈q(x)〉c, respectively, conditioned on measurements of K(x). These predictors have been shown by Neuman and Orr [1993a] to satisfy exactly an integrodifferential conditional mean flow equation in which 〈q(x)〉c is nonlocal and non-Darcian. Here we develop complementary integrodifferential equations for second conditional moments of head and flux which serve as measures of predictive uncertainty; obtain recursive closure approximations for both the first and second conditional moment equations through expansion in powers of a small parameter σY which represents the standard estimation error of ln K(x); and show how to solve these equations to first order in σY2 by finite elements on a rectangular grid in two dimensions. In the special case where one treats K(x) as if it was locally homogeneous and mean flow as if it was locally uniform, one obtains a localized Darcian approximation 〈q(x)〉c ≈ −Kc(x)∇〈h(x)〉c in which Kc(x) is a space-dependent conditional hydraulic conductivity tensor. This leads to the traditional deterministic, Darcian steady state flow equation which, however, acquires a nontraditional meaning in that its parameters and state variables are data dependent and therefore inherently nonunique. It further explains why parameter estimates obtained by traditional inverse methods tend to vary as one modifies the database. Localized equations yield no information about predictive uncertainty. Our stochastic derivation of these otherwise standard deterministic flow equations makes clear that uncertainty measures associated with estimates of head and flux, obtained by traditional inverse methods, are generally smaller (often considerably so) than measures of corresponding predictive uncertainty, which can be assessed only by means of stochastic models such as ours. We present a detailed comparison between finite element solutions of nonlocal and localized moment equations and Monte Carlo simulations under superimposed mean-uniform and convergent flow regimes in two dimensions. Paper 1 presents the theory and computational approach, and paper 2 [Guadagnini and Neuman, this issue] describes unconditional and conditional computational results.

A Eulerian-Lagrangian numerical scheme for the dispersion-convection equation using conjugate space-time grids

Abstract A new numerical scheme is proposed for the dispersion-convection equation which combines the utility of a fixed grid in Eulerian coordinates with the computational power of the Lagrangian method. Convection is formally decoupled from dispersion in a manner which does not leave room for ambiguity. The resulting convection problem is solved by the method of characteristics on a grid fixed in space. Dispersion is handled by finite elements on a separate grid which may, but need not, coincide wit the former at selected points in spacetime. Information is projected from one grid to another by local interpolation. The conjugate grid method is implemented by linear finite elements in conjunction with piecewise linear interpolation functions and applied to five problems ranging from predominant dispersion to pure convection. The results are entirely free of oscillations. Numerical dispersion exists but can be brought under control either by reducing the spatial increment, or by increasing the time step size, of the grid used to solve the convection problem. Contrary to many other methods, best results are often obtained when the Courant number exceeds 1.

Three‐dimensional numerical inversion of pneumatic cross‐hole tests in unsaturated fractured tuff: 2. Equivalent parameters, high‐resolution stochastic imaging and scale effects

In paper 1 of this two-part series we described a three-dimensional numerical inverse model for the interpretation of cross-hole pneumatic tests in unsaturated fractured tuffs at the Apache Leap Research Site (ALRS) near Superior, Arizona. Our model is designed to analyze these data in two ways: (1) by considering pressure records from individual borehole monitoring intervals one at a time, while treating the rock as being spatially uniform, and (2) by considering pressure records from multiple tests and borehole monitoring intervals simultaneously, while treating the rock as being randomly heterogeneous. The first approach yields a series of equivalent air permeabilities and air- filled porosities for rock volumes having length scales ranging from meters to tens of meters, represented nominally by radius vectors extending from injection to monitoring intervals. The second approach yields a high-resolution geostatistical estimate of how air permeability and air-filled porosity, defined on grid blocks having a length scale of 1 m, vary spatially throughout the tested rock volume. It amounts to three-dimensional pneumatic “tomography” or stochastic imaging of the rock. Paper 1 described the field data, the model, and the effect of boreholes on pressure propagation through the rock. This second paper implements our inverse model on pressure data from five cross-hole tests at ALRS. We compare our cross-hole test interpretations by means of the two approaches with earlier interpretations by means of type curves and with geostatistical interpretations of single-hole test data. The comparisons show internal consistency between all pneumatic test interpretations and reveal a very pronounced scale effect in permeability and porosity at ALRS.

Theoretical derivation of Darcy's law

SummaryDarcys law for anisotropic porous media is derived from the Navier-Stokes equation by using a formal averaging procedure. Particular emphasis is placed upon the proof that the permeability tensor is symmetric. In addition, it is shown that there is a one-to-one relationship between the local and macroscopic velocity fields. This leads to the interesting phenomenological observation that the local velocity vector at any given point must always lie either on a fixed line or in a fixed plane. All of this holds true for an incompressible homogeneous Newtonian fluid moving slowly through a rigid porous medium with uniform porosity under isothermal and steady state conditions. The question whether Darcys law is applicable under nonsteady or compressible flow conditions, or when the medium has nonuniform porosity, is also discussed. Finally, it is shown that the Hagen-Poiseuille equation, as well as the expression describing Couette flow between parallel plates, can be derived from the equations presented in this work and may thus be viewed as special cases of Darcys law.ZusammenfassungDas Darcysche Gesetz für anisotrope poröse Werkstoffe wird von den Navier-Stokes-Gleichungen durch formale Mittelwertbildung abgeleitet. Insbesondere betont wird der Beweis der Symmetrie des Permeabilitätstensors. Weiter wird gezeigt, daß eine eineindeutige Beziehung zwischen lokalen und makroskopischen Geschwindigkeitsfeldern existiert. Dies führt zur interessanten phänomenologischen Beobachtung, daß in jedem Punkt der lokale Geschwindigkeitsvektor entweder auf einer festen Geraden oder in einer festen Ebene liegt. All dies gilt für inkompressible homogene Newtonsche Flüssigkeiten, die sich langsam, stationär unter isothermen Bedingungen durch einen starren porösen Körper gleichförmiger Porösität bewegen. Die Frage, ob das Darcysche Gesetz für instationäre Strömungen oder kompressible Fälle oder für ungleichförmige Porösität gilt, wird ebenfalls diskutiert. Abschließend wird gezeigt, daß die Hagen-Poiseuille-Gleichung und der Ausdruck für die Couette-Strömung zwischen parallelen Platten von der in dieser Arbeit angegebenen Gleichung abgeleitet und daher als Spezialfälle des Darcyschen Gesetzes betrachtet werden können.

Analysis of nonintrinsic spatial variability by residual kriging with application to regional groundwater levels

A method for obtaining pointwise or spatially averaged estimates of a nonintrinsic function is introduced based on residual kriging. The method relies on a stepwise iterative regression process for simultaneously estimating the global drift and residual semivariogram. Estimates of the function are then obtained by solving a modified set of simple kriging equations written for the residuals. The modification consists of replacing the true variogram in the kriging equations by the variogram of the residual estimates as obtained from the iterative regression process. The method is illustrated by considering groundwater levels in an Arizona aquifer. The results are compared with those obtained for the aquifer by the generalized covariance package BLUEPACK-3D.