Shlomo P Neuman

Interests

Teaching

My teaching philosophy is rooted in the belief that hydrology requires the application of rigorous thinking to complex but often poorly defined scientific and engineering problems. I therefore construct my courses around a rigorous theoretical framework and exploration of its applicability to real-world problems. As a result, my courses are an amalgam of theoretical and applied. In the area of theory, I place equal emphases on qualitative understanding of physical processes and their mathematical treatment. In the applied area, I emphasize methodology and illustrate it, where possible, by case studies and synthetic examples. My courses are reputed to be challenging on both the conceptual and mathematical levels. Though I believe that every hydrology student should be able to rise to this challenge few dare to do so. This, and the fact that my courses are designated as electives, explains why enrollment in them is generally low.

Research

Neumans research centers around theoretical, computational and field analyses of fluid flow and solute transport in porous and fractured geologic media. It reflects a philosophy according to which hydrology requires the application of rigorous thinking to complex but often poorly defined scientific and engineering problems. Neuman therefore conducts his research within a rigorous theoretical framework which he nevertheless renders applicable to real-world problems. His research results are thus an amalgam of science and engineering, of theoretical and applied. They have led to major advances in our conceptual understanding of flow and transport processes under a broad range of field conditions, within a wide range of hydrogeologic environments, that are of both scientific significance and practical relevance; in our ability to describe and generalize such processes mathematically; in our ability to simulate and analyze them computationally; and in our ability to measure and observe them in the field.

Courses

No activities entered.

Scholarly Contributions

Journals/Publications

• Guadagnini, A., Neuman, S. P., Schaap, M. G., & Riva, M. (2014). Anisotropic statistical scaling of soil and sediment texture in a stratified deep vadose zone near Maricopa, Arizona. Geoderma, 214-215, 217-227.
  More info

  Abstract: We analyze the statistical scaling properties of vertical and horizontal increments in soil and sediment texture data measured to a depth of 15m over an area of 3600m2 in a vadose zone near Maricopa, Arizona. The data include sand, silt and clay fractions, their principal components and logit transforms. The statistical scaling properties we uncover are difficult to detect with standard geostatistical methods. They include (a) pronounced intermittency (rough, irregular spatial variability) and antipersistence (tendency for large and small values to alternate rapidly); (b) symmetric, non-Gaussian frequency distributions characterized by heavy tails that often decay with separation distance or lag; (c) nonlinear power-law scaling of sample structure functions (statistical moments of absolute increments) in a midrange of lags, with breakdown in such scaling at small and large lags; (d) extended power-law scaling (linear relations between log structure functions of successive orders) at all lags; (e) nonlinear scaling of power-law exponent with order of sample structure function; and (f) pronounced anisotropy in these behaviors. Similar statistical scaling behaviors are known to be exhibited by a wide variety of earth, environmental and other variables (including ecological, biological, physical, astrophysical and financial). The literature has traditionally interpreted this to imply that the variables are multifractal, which however explains neither the observed breakdown in power-law scaling at small and large lags nor extended power-law scaling. We offer an alternative interpretation that is simpler and consistent with all the above phenomena. Our interpretation views the data as samples from stationary, anisotropic sub-Gaussian random fields subordinated to truncated fractional Brownian motion (tfBm) or truncated fractional Gaussian noise (tfGn). Such sub-Gaussian fields are mixtures of Gaussian fields with random variances. Truncation of monofractal fBm (which is non-stationary) and fGn entails filtering out components below the measurement or resolution scale of the data and above the scale of their sampling domain. Our novel interpretation of the data allows us to obtain maximum likelihood estimates of all parameters characterizing the underlying truncated sub-Gaussian fields. These parameters in turn make it possible to downscale or upscale all statistical moments to situations entailing smaller or larger measurement or resolution and sampling scales, respectively. They also allow one to perform conditional or unconditional Monte Carlo simulations of random field realizations corresponding to these scales. © 2013 Elsevier B.V.
• Neuman, S. P. (2014). Comment on "The New Potential for Understanding Groundwater Contaminant Transport" by Paul W. Hadley and Charles Newell. Ground Water, 52(5), 653-656.
• Neuman, S. P., Xue, L., Zhang, D., & Guadagnini, a. (2014). Multimodel Bayesian Analysis of Groundwater Data Worth. Water Resources Research, 50(11), 8481-8496.
• Panzeri, M., Riva, M., Guadagnini, A., & Neuman, S. P. (2014). Comparison of Ensemble Kalman Filter groundwater-data assimilation methods based on stochastic moment equations and Monte Carlo simulation. Advances in Water Resources, 66, 8-18.
  More info

  Abstract: Traditional Ensemble Kalman Filter (EnKF) data assimilation requires computationally intensive Monte Carlo (MC) sampling, which suffers from filter inbreeding unless the number of simulations is large. Recently we proposed an alternative EnKF groundwater-data assimilation method that obviates the need for sampling and is free of inbreeding issues. In our new approach, theoretical ensemble moments are approximated directly by solving a system of corresponding stochastic groundwater flow equations. Like MC-based EnKF, our moment equations (ME) approach allows Bayesian updating of system states and parameters in real-time as new data become available. Here we compare the performances and accuracies of the two approaches on two-dimensional transient groundwater flow toward a well pumping water in a synthetic, randomly heterogeneous confined aquifer subject to prescribed head and flux boundary conditions. © 2014 Elsevier Ltd.
• Guadagnini, A., Neuman, S. P., Schaap, M. G., & Riva, M. (2013). Anisotropic statistical scaling of vadose zone hydraulic property estimates near Maricopa, Arizona. Water Resources Research, 49(12), 8463-8479.
  More info

  Abstract: Fluid flow and mass transport in the vadose zone are strongly influenced by spatial variability in soil hydraulic properties. It has become common to characterize spatial variability geostatistically and to solve corresponding flow and/or transport problems stochastically. The typical (though not only) approach is to treat log saturated hydraulic conductivity, Y=log10 Ks, and perhaps other medium properties as statistically homogeneous, isotropic, or anisotropic multivariate-Gaussian random fields with unique variance and autocorrelation scale. A growing body of literature suggests that Y as well as many other variables may possess heavy-tailed, non-Gaussian distributions, and/or scale-dependent statistical parameters representing a multiscale, hierarchical structure. Elsewhere, we have demonstrated that these important phenomena are difficult to detect, and are not fully tractable, with standard geostatistical methods. They are however detectable and tractable with a novel geostatistical method of analysis which treats such variables as samples from sub-Gaussian random fields subordinated to truncated fractional Brownian motion (tfBm) or truncated fractional Gaussian noise (tfGn). Each such field is a mixture of Gaussian components having random variances. The purpose of this paper is to explore the extent to which hydraulic parameters of unsaturated soils exhibit heavy-tailed distributions and statistical scaling of the above type. As a test bed we have selected an experimental site near Maricopa, Arizona, at which soil texture data have been collected in boreholes and trenches to a depth of 15 m over an area of 3600 m2. Since hydraulic parameters are difficult to measure to such depths we estimate them using a neural network pedotransfer model. Input data include percent sand, silt, and clay. Outputs include estimates of saturated and residual volumetric water content, saturated hydraulic conductivity, and (due to their wide applicability) parameters of the van Genuchten-Mualem constitutive model of unsaturated soil property variations with capillary pressure head. We find indeed that vertical and horizontal increments of our hydraulic parameter estimates exhibit heavy-tailed frequency distributions and statistical scaling which conform closely to our novel geostatistical framework. Its generality and applicability to the Maricopa site has far reaching implications vis-à-vis the geostatistical characterization of vadose zone properties and the stochastic modeling of flow and transport in unsaturated media at other sites. © 2013. American Geophysical Union. All Rights Reserved.
• Guadagnini, A., Neuman, S. P., Schaap, M. G., & Riva, M. (2013). Statistical and scaling analyses of neural network soil property inputs/outputs at an arizona field site. SIMULTECH 2013 - Proceedings of the 3rd International Conference on Simulation and Modeling Methodologies, Technologies and Applications, 489-494.
  More info

  Abstract: Analyses of flow and transport in the shallow subsurface require information about spatial and statistical distributions of soil hydraulic properties (water content and permeability, their dependence on capillary pressure) as functions of scale and direction. Measuring these properties is relatively difficult, time consuming and costly. It is generally much easier, faster and less expensive to collect and describe the makeup of soil samples in terms of textural composition (e.g. per cent sand, silt, clay and organic matter), bulk density and other such pedological attributes. Over the last two decades soil scientists have developed a set of tools, known collectively as pedotransfer functions (PTFs), to help translate information about the spatial distribution of pedological indicators into corresponding information about soil hydraulic properties. One of the most successful PTFs is the nonlinear Rosetta neural network model developed by one of us. Among remaining open questions are the extents to which spatial and statistical distributions of Rosetta hydraulic property outputs, and their scaling behavior, reflect those of Rosetta pedological inputs. We address the last question by applying Rosetta, coupled with a novel statistical scaling analysis recently proposed by three of us, to soil sample data from an experimental site in southern Arizona, USA.
• Panzeri, M., Riva, M., Guadagnini, A., & Neuman, S. P. (2013). Data assimilation and parameter estimation via ensemble Kalman filter coupled with stochastic moment equations of transient groundwater flow. Water Resources Research, 49(3), 1334-1344.
  More info

  Abstract: The ensemble Kalman filter (EnKF) is a powerful tool for assimilating data in earth system models. The approach allows real time Bayesian updating of system states and parameters as new data become available. This paper focuses on EnKF data assimilation in models of groundwater flow through complex geologic media. It has become common to treat the hydraulic conductivity of such media as correlated random fields conditioned on measured conductivity (medium property) and/or hydraulic head (system state) values. This renders the conductivity nonstationary and the corresponding conditional flow equations stochastic. Solving these equations and coupling them with EnKF generally entails computationally intensive Monte Carlo (MC) simulation. We propose to circumvent the need for MC through a direct solution of approximate nonlocal (integrodifferential) equations that govern the space-time evolution of conditional ensemble means (statistical expectations) and covariances of hydraulic heads and fluxes. We illustrate and explore our approach on synthetic two-dimensional examples in which a well pumps water from a randomly heterogeneous aquifer subject to prescribed head and flux boundary conditions. Embedding the solution in EnKF provides sequential updates of conductivity and head estimates throughout the space-time domain of interest. We demonstrate the computational feasibility and accuracy of our methodology, showing that hydraulic conductivity estimates are more sensitive to early than to later head values and improve with increasing assimilation frequency at early time. ©2013. American Geophysical Union. All Rights Reserved.
• Riva, M., Neuman, S. P., & Guadagnini, A. (2013). Sub-Gaussian model of processes with heavy-tailed distributions applied to air permeabilities of fractured tuff. Stochastic Environmental Research and Risk Assessment, 27(1), 195-207.
  More info

  Abstract: Earth and environmental variables are commonly taken to have multivariate Gaussian or heavy-tailed distributions in space and/or time. This is based on the observation that univariate frequency distributions of corresponding samples appear to be Gaussian or heavy-tailed. Of particular interest to us is the well-documented but heretofore little noticed and unexplained phenomenon that whereas the frequency distribution of log permeability data often seems to be Gaussian, that of corresponding increments tends to exhibit heavy tails. The tails decay as powers of -α where 1 lt; α lt; 2 is either constant or grows monotonically toward an asymptote with increasing separation distance or lag. We illustrate the latter phenomenon on 1-m scale log air permeabilities from pneumatic tests in 6 vertical and inclined boreholes completed in unsaturated fractured tuff near Superior, Arizona. We then show theoretically and demonstrate numerically, on synthetically generated signals, that whereas the case of constant α is consistent with a collection of samples from truncated sub-Gaussian fractional Lévy noise, a random field (or process) subordinated to truncated fractional Gaussian noise, the case of variable α is consistent with a collection of samples from truncated sub-Gaussian fractional Lévy motion (tfLm), a random field subordinated to truncated fractional Brownian motion. Whereas the first type of signal is relatively regular and characterized by Lévy index α, the second is highly irregular (punctuated by spurious spikes) and characterized by the asymptote of α values associated with its increments. We describe a procedure to estimate the parameters of univariate distributions characterizing such signals and apply it to our log air permeability data. The latter are found to be consistent with a collection of samples from tfLm with α slightly smaller than 2, which is easily confused with a Gaussian field (characterized by constant α = 2). The irregular (spiky) nature of this signal is typical of observed fractured rock properties. We propose that distributions of earth and environmental variable be inferred jointly from measured values and their increments in a way that insures consistency between these two sets of data. © 2012 Springer-Verlag.
• Riva, M., Neuman, S. P., Guadagnini, A., & Siena, M. (2013). Anisotropic scaling of Berea sandstone log air permeability statistics. Vadose Zone Journal, 12(3).
  More info

  Abstract: Minipermeameter Berea sandstone data and their increments are found to be consistent with sub-Gaussian random fields subordinated to truncated fractional Gaussian noise. Their statistics and scaling behavior are shown to depend on direction. We analyze directional dependence and scaling of air log permeability statistics characterizing minipermeameter measurements by Tidwell and Wilson on the faces of a 30 by 30 by 30 cm3 size unsaturated block of Berea sandstone. We find distinct differences between the statistics and scaling behaviors of data measured on faces parallel and normal to bedding, and of incremental measurements in three orthogonal directions along these faces. Whereas the distribution of data and their increments parallel to bedding are heavy tailed, those of increments normal to bedding are Gaussian. Order q sample structure functions of increments parallel to bedding scale as powers ξ(q) of directional lag, sd, over limited ranges of sd. Using moment and extended self-similarity methods of analysis we find ξ(q) to be generally nonlinear in q, tending to be concave in q on faces normal to bedding and convex on faces parallel to bedding. Whereas the literature attributes nonlinear scaling of ξ(q) with q to multifractals or fractional Laplace motions, we find the data to be consistent with sub-Gaussian random fields subordinated to truncated (monofractal, self-affine) fractional Gaussian noise. The increments exhibit negative statistical dependence (antipersistence) that varies with direction parallel to bedding and is more pronounced on faces parallel than normal to bedding. © Soil Science Society of America 5585 Guilford Rd., Madison, WI 53711 USA. All rights reserved.
• Dan, L. u., Ming, Y. e., Neuman, S. P., & Xue, L. (2012). Multimodel Bayesian analysis of data-worth applied to unsaturated fractured tuffs. Advances in Water Resources, 35, 69-82.
  More info

  Abstract: To manage water resource and environmental systems effectively requires suitable data. The worth of collecting such data depends on their potential benefit and cost, including the expected cost (risk) of failing to take an appropriate decision. Evaluating this risk calls for a probabilistic approach to data-worth assessment. Recently we [39] developed a multimodel approach to optimum value-of-information or data-worth analysis based on model averaging within a maximum likelihood Bayesian framework. Adopting a two-dimensional synthetic example, we implemented our approach using Monte Carlo (MC) simulations with and without lead order approximations, finding that the former approach was almost equally accurate but computationally more efficient. Here we apply our methodology to pneumatic permeability data from vertical and inclined boreholes drilled into unsaturated fractured tuff near Superior, Arizona. In an attempt to improve computational efficiency, we introduce three new approximations that require less computational effort and compare results with those obtained by the original Monte Carlo method. The first approximation disregards uncertainty in model parameter estimates, the second does so for estimates of potential new data, and the third disregards both uncertainties. We find that only the first approximation yields reliable quantitative assessments of reductions in predictive uncertainty brought about by the collection of new data. We conclude that, whereas parameter uncertainty may sometimes be disregarded for purposes of analyzing data worth, the same does not generally apply to uncertainty in estimates of potential new data. © 2011 Elsevier Ltd.
• Guadagnini, A., Neuman, S. P., & Riva, M. (2012). Numerical investigation of apparent multifractality of samples from processes subordinated to truncated fBm. Hydrological Processes, 26(19), 2894-2908.
  More info

  Abstract: We investigate numerically apparent multi-fractal behavior of samples from synthetically generated processes subordinated to truncated fractional Brownian motion (tfBm) on finite domains. We are motivated by the recognition that many earth and environmental (including hydrological) variables appear to be self-affine (monofractal) or multifractal with Gaussian or heavy-tailed distributions. The literature considers self-affine and multifractal types of scaling to be fundamentally different, the first arising from additive and the second from multiplicative random fields or processes. It has been demonstrated theoretically by one of us that square or absolute increments of samples from Gaussian/Lévy processes subordinated to tfBm exhibit apparent/spurious multifractality at intermediate ranges of separation lags, with breakdown in power-law scaling at small and large lags as is commonly exhibited by real data. A preliminary numerical demonstration of apparent multifractality by the same author was limited to Gaussian fields having nearest neighbor autocorrelations and led to rather noisy results. Here, we adopt a new generation scheme that allows us to investigate apparent multifractal behaviors of samples taken from a broad range of processes including Gaussian with and without symmetric Lévy and log-normal (as well as potentially other) subordinators. Our results shed new light on the nature of apparent multifractality, which has wide implications vis-a-vis the scaling of many hydrological as well as other earth and environmental variables. © 2011 John Wiley & Sons, Ltd.
• Guadagnini, A., Riva, M., & Neuman, S. P. (2012). Extended power-law scaling of heavy-tailed random air-permeability fields in fractured and sedimentary rocks. Hydrology and Earth System Sciences, 16(9), 3249-3260.
  More info

  Abstract: We analyze the scaling behaviors of two field-scale log permeability data sets showing heavy-tailed frequency distributions in three and two spatial dimensions, respectively. One set consists of 1-m scale pneumatic packer test data from six vertical and inclined boreholes spanning a decameters scale block of unsaturated fractured tuffs near Superior, Arizona, the other of pneumatic minipermeameter data measured at a spacing of 15 cm along three horizontal transects on a 21 m long and 6 m high outcrop of the Upper Cretaceous Straight Cliffs Formation, including lower-shoreface bioturbated and cross-bedded sandstone near Escalante, Utah. Order q sample structure functions of each data set scale as a power xi(q)of separation scale or lag, s, over limited ranges of s. A procedure known as extended self-similarity (ESS) extends this range to all lags and yields a nonlinear (concave) functional relationship between xi(q) and. Whereas the literature tends to associate extended and nonlinear power-law scaling with multifractals or fractional Laplace motions, we have shown elsewhere that (a) ESS of data having a normal frequency distribution is theoretically consistent with (Gaussian) truncated (additive, self-affine, monofractal) fractional Brownian motion (tfBm), the latter being unique in predicting a breakdown in power-law scaling at small and large lags, and (b) nonlinear power-law scaling of data having either normal or heavy-tailed frequency distributions is consistent with samples from sub-Gaussian random fields or processes subordinated to tfBm or truncated fractional Gaussian noise (tfGn), stemming from lack of ergodicity which causes sample moments to scale differently than do their ensemble counterparts. Here we (i) demonstrate that the above two data sets are consistent with sub-Gaussian random fields subordinated to tfBm or tfGn and (ii) provide maximum likelihood estimates of parameters characterizing the corresponding Lávy stable subordinators and tfBm or tfGn functions. © 2012 Author(s).
• Mishra, P. K., Vesselinov, V. V., & Neuman, S. P. (2012). Radial flow to a partially penetrating well with storage in an anisotropic confined aquifer. Journal of Hydrology, 448-449, 255-259.
  More info

  Abstract: Drawdowns generated by extracting water from large diameter (e.g. water supply) well are affected by wellbore storage. We present an analytical solution in Laplace transformed space for drawdown in a uniform anisotropic aquifer caused by withdrawing water at a constant rate from partially penetrating well with storage. The solution is back transformed into the time domain numerically. When the pumping well is fully penetrating our solution reduces to that of . Papadopulos and Cooper (1967); Hantush (1964) when the pumping well has no wellbore storage; . Theis (1935) when both conditions are fulfilled and . Yang (2006) when the pumping well is partially penetrating, has finite radius but lacks storage. Newly developed solution is then used to explore graphically the effects of partial penetration, wellbore storage and anisotropy on time evolutions of drawdown in the pumping well and in observation wells. We concluded after validating the developed analytical solution using synthetic pumping test. © 2012 Elsevier B.V..
• Neuman, S. P., & Mishra, P. K. (2012). Comments on "a revisit of drawdown behavior during pumping in unconfined aquifers" by D. Mao, L. Wan, T.-C. J. Yeh, C.-H. Lee, K.-C. Hsu, J.-C. Wen, and W. Lu. Water Resources Research, 48(2).
• Neuman, S. P., Xue, L., Ming, Y. e., & Dan, L. u. (2012). Bayesian analysis of data-worth considering model and parameter uncertainties. Advances in Water Resources, 36, 75-85.
  More info

  Abstract: The rational management of water resource systems requires an understanding of their response to existing and planned schemes of exploitation, pollution prevention and/or remediation. Such understanding requires the collection of data to help characterize the system and monitor its response to existing and future stresses. It also requires incorporating such data in models of system makeup, water flow and contaminant transport. As the collection of subsurface characterization and monitoring data is costly, it is imperative that the design of corresponding data collection schemes be cost-effective, i.e., that the expected benefit of new information exceed its cost. A major benefit of new data is its potential to help improve one's understanding of the system, in large part through a reduction in model predictive uncertainty and corresponding risk of failure. Traditionally, value-of-information or data-worth analyses have relied on a single conceptual-mathematical model of site hydrology with prescribed parameters. Yet there is a growing recognition that ignoring model and parameter uncertainties render model predictions prone to statistical bias and underestimation of uncertainty. This has led to a recent emphasis on conducting hydrologic analyses and rendering corresponding predictions by means of multiple models. We describe a corresponding approach to data-worth analyses within a Bayesian model averaging (BMA) framework. We focus on a maximum likelihood version (MLBMA) of BMA which (a) is compatible with both deterministic and stochastic models, (b) admits but does not require prior information about the parameters, (c) is consistent with modern statistical methods of hydrologic model calibration, (d) allows approximating lead predictive moments of any model by linearization, and (e) updates model posterior probabilities as well as parameter estimates on the basis of potential new data both before and after such data become actually available. We describe both the BMA and MLBMA versions theoretically and implement MLBMA computationally on a synthetic example with and without linearization. © 2011 Elsevier Ltd.
• Siena, M., Guadagnini, A., Riva, M., & Neuman, S. P. (2012). Extended power-law scaling of air permeabilities measured on a block of tuff. Hydrology and Earth System Sciences, 16(1), 29-42.
  More info

  Abstract: We use three methods to identify power-law scaling of multi-scale log air permeability data collected by Tidwell and Wilson on the faces of a laboratory-scale block of Topopah Spring tuff: method of moments (M), Extended Self-Similarity (ESS) and a generalized version thereof (GESS). All three methods focus on q-th-order sample structure functions of absolute increments. Most such functions exhibit power-law scaling at best over a limited midrange of experimental separation scales, or lags, which are sometimes difficult to identify unambiguously by means of M. ESS and G-ESS extend this range in a way that renders power-law scaling easier to characterize. Our analysis confirms the superiority of ESS and G-ESS over M in identifying the scaling exponents, ξ(q), of corresponding structure functions of orders q, suggesting further that ESS is more reliable than G-ESS. The exponents vary in a nonlinear fashion with q as is typical of real or apparent multifractals. Our estimates of the Hurst scaling coefficient increase with support scale, implying a reduction in roughness (anti-persistence) of the log permeability field with measurement volume. The finding by Tidwell and Wilson that log permeabilities associated with all tip sizes can be characterized by stationary variogram models, coupled with our findings that log permeability increments associated with the smallest tip size are approximately Gaussian and those associated with all tip sizes scale show nonlinear variations in ξ(q) with q, are consistent with a view of these data as a sample from a truncated version (tfBm) of self-affine fractional Brownian motion (fBm). Since in theory the scaling exponents, ξ(q), of tfBm vary linearly with q we conclude that nonlinear scaling in our case is not an indication of multifractality but an artifact of sampling from tfBm. This allows us to explain theoretically how power-law scaling of our data, as well as of non-Gaussian heavy-tailed signals subordinated to tfBm, are extended by ESS. It further allows us to identify the functional form and estimate all parameters of the corresponding tfBm based on sample structure functions of first and second orders. © Author(s) 2012.
• Dan, L. u., Ming, Y. e., & Neuman, S. P. (2011). Dependence of Bayesian Model Selection Criteria and Fisher Information Matrix on Sample Size. Mathematical Geosciences, 43(8), 971-993.
  More info

  Abstract: Geostatistical analyses require an estimation of the covariance structure of a random field and its parameters jointly from noisy data. Whereas in some cases (as in that of a Matérn variogram) a range of structural models can be captured with one or a few parameters, in many other cases it is necessary to consider a discrete set of structural model alternatives, such as drifts and variograms. Ranking these alternatives and identifying the best among them has traditionally been done with the aid of information theoretic or Bayesian model selection criteria. There is an ongoing debate in the literature about the relative merits of these various criteria. We contribute to this discussion by using synthetic data to compare the abilities of two common Bayesian criteria, BIC and KIC, to discriminate between alternative models of drift as a function of sample size when drift and variogram parameters are unknown. Adopting the results of Markov Chain Monte Carlo simulations as reference we confirm that KIC reduces asymptotically to BIC and provides consistently more reliable indications of model quality than does BIC for samples of all sizes. Practical considerations often cause analysts to replace the observed Fisher information matrix entering into KIC with its expected value. Our results show that this causes the performance of KIC to deteriorate with diminishing sample size. These results are equally valid for one and multiple realizations of uncertain data entering into our analysis. Bayesian theory indicates that, in the case of statistically independent and identically distributed data, posterior model probabilities become asymptotically insensitive to prior probabilities as sample size increases. We do not find this to be the case when working with samples taken from an autocorrelated random field. © 2011 International Association for Mathematical Geosciences.
• Guadagnini, A., & Neuman, S. P. (2011). Extended power-law scaling of self-affine signals exhibiting apparent multifractality. Geophysical Research Letters, 38(13).
  More info

  Abstract: Many earth and environmental variables appear to scale as multiplicative (multifractal) processes with spatial or temporal increments possessing Gaussian or heavy-tailed distributions. The behavior, characterized by power-law scaling, is typically limited to intermediate ranges of separation scales (lags) considered, in the case of fully developed turbulence, to be dominated by inertia. It has been established empirically that, in numerous cases (e.g. turbulence, diffusion-limited aggregates, natural images, kinetic surface roughening, fluvial turbulence, sand wave dynamics, Martian topography, river morphometry, gravel-bed mobility, barometric pressure, low-energy cosmic rays, cosmic microwave background radiation, metal-insulator transition, irregularities in human heartbeat time series, turbulence in edge magnetized plasma of fusion devices and turbulent boundary layers of the Earth's magnetosphere), this range of lags can be enlarged significantly, at both ends of the spectrum, via a procedure known as Extended Self-Similarity (ESS). We demonstrate numerically that a similar procedure extends the power-law scaling range over which additive (self-affine) signals exhibit apparent multifractality. We conclude that signals appearing to exhibit either standard or extended (such as those listed) multifractal scaling may potentially represent self-affine processes. Copyright 2011 by the American Geophysical Union.
• Mishra, P. K., & Neuman, S. P. (2011). Saturated-unsaturated flow to a well with storage in a compressible unconfined aquifer. Water Resources Research, 47(5).
  More info

  Abstract: Mishra and Neuman (2010) developed an analytical solution for flow to a partially penetrating well of zero radius in a compressible unconfined aquifer that allows inferring its saturated and unsaturated hydraulic properties from responses recorded in the saturated and/or unsaturated zones. Their solution accounts for horizontal as well as vertical flows in each zone. It represents unsaturated zone constitutive properties in a manner that is at once mathematically tractable and sufficiently flexible to provide much improved fits to standard constitutive models. In this paper we extend the solution of Mishra and Neuman [2010] to the case of a finite diameter pumping well with storage; investigate the effects of storage in the pumping well and delayed piezometer response on drawdowns in the saturated and unsaturated zones as functions of position and time; validate our solution against numerical simulations of drawdown in a synthetic aquifer having unsaturated properties described by the van Genuchten [1980]-Mualem [1976] model; use our solution to analyze 11 transducer-measured drawdown records from a seven-day pumping test conducted by University of Waterloo researchers at the Canadian Forces Base Borden in Ontario, Canada; validate our parameter estimates against manually-measured drawdown records in 14 other piezometers at Borden; and compare (a) our estimates of aquifer parameters with those obtained on the basis of all these records by Moench [2008], (b) on the basis of 11 transducer-measured drawdown records by Endres et al. [2007], (c) our estimates of van Genuchten-Mualem parameters with those obtained on the basis of laboratory drainage data from the site by Akindunni and Gillham [1992], and (d) our corresponding prediction of how effective saturation varies with elevation above the initial water table under static conditions with a profile based on water contents measured in a neutron access tube at a radial distance of about 5 m from the center of the pumping well. We also use our solution to analyze 11 transducer-measured drawdown records from a 7 day pumping test conducted by University of Waterloo researchers at the Canadian Forces Base Borden in Ontario, Canada. We validate our parameter estimates against manually measured drawdown records in 14 other piezometers at Borden. We compare our estimates of aquifer parameters with those obtained on the basis of all these records by Moench (2008) and on the basis of 11 transducer-measured drawdown records by Endres et al. (2007), and we compare our estimates of van Genuchten-Mualem parameters with those obtained on the basis of laboratory drainage data from the site by Akindunni and Gillham (1992); finally, we compare our corresponding prediction of how effective saturation varies with elevation above the initial water table under static conditions with a profile based on water contents measured in a neutron access tube at a radial distance of about 5 m from the center of the pumping well. © 2011 by the American Geophysical Union.
• Neuman, S. P. (2011). Apparent multifractality and scale-dependent distribution of data sampled from self-affine processes. Hydrological Processes, 25(11), 1837-1840.
  More info

  Abstract: It has been previously demonstrated theoretically and numerically by the author that square or absolute increments of data sampled from fractional Brownian/Lévy motion (fBm/fLm), or of incremental data sampled from fractional Gaussian/Lévy noise (fGn/fLn), exhibit apparent/spurious multifractality. Here, we generalize these previous development in a way that (a) rigorously subordinates (truncated) fLn to fGn or, in a statistically equivalent manner, (truncated) fLm to fBm; (b) extends the analysis to a wider class of subordinated self-affine processes; (c) provides a simple way to generate such processes and (d) explains why the distribution of corresponding increments tends to evolve from heavy tailed at small lags (separation distances or scales) to Gaussian at larger lags. © 2011 John Wiley & Sons, Ltd.
• Riva, M., Panzeri, M., Guadagnini, A., & Neuman, S. P. (2011). Role of model selection criteria in geostatistical inverse estimation of statistical data- and model-parameters. Water Resources Research, 47(7).
  More info

  Abstract: We analyze theoretically the ability of model quality (sometimes termed information or discrimination) criteria such as the negative log likelihood NLL, Bayesian criteria BIC and KIC and information theoretic criteria AIC, AICc, and HIC to estimate (1) the parameter vector θ of the variogram of hydraulic log conductivity (Y = ln K), and (2) statistical parameters σ hE2 and σ YE2 proportional to head and log conductivity measurement error variances, respectively, in the context of geostatistical groundwater flow inversion. Our analysis extends the work of Hernandez et al. (2003, 2006) and Riva et al. (2009), who developed nonlinear stochastic inverse algorithms that allow conditioning estimates of steady state and transient hydraulic heads, fluxes and their associated uncertainty on information about conductivity and head data collected in a randomly heterogeneous confined aquifer. Their algorithms are based on recursive numerical approximations of exact nonlocal conditional equations describing the mean and (co)variance of groundwater flow. Log conductivity is parameterized geostatistically based on measured values at discrete locations and unknown values at discrete "pilot points." Optionally, the maximum likelihood function on which the inverse estimation of Y at pilot points is based may include a regularization term reflecting prior information about Y. The relative weight λ = σ hE2/σ YE2 assigned to this term and its components σ hE2 and σ YE2, as well as θ are evaluated separately from other model parameters to avoid bias and instability. This evaluation is done on the basis of criteria such as NLL, KIC, BIC, HIC, AIC, and AICc. We demonstrate theoretically that, whereas all these six criteria make it possible to estimate σ hE2, KIC alone allows one to estimate validly θ and σ YE2 (and thus λ. We illustrate this discriminatory power of KIC numerically by using a differential evolution genetic search algorithm to minimize it in the context of a two-dimensional steady state groundwater flow problem. We find that whereas σ hE2, σ YE2, and the integral scale of Y can be estimated on the basis of a zero-order mean flow equation, the sill of the Y-variogram is estimated more accurately by a second-order approximation of flow. This notwithstanding, KIC prefers the simpler zero-order moment over the more complex second-order version. © 2011 by the American Geophysical Union.
• Wang, X., Neuman, S. P., Strack, O. D., Verruijt, A., Jamali, M., Seymour, B., Bear, J., Cheng, A. H., Chen, C., Kuang, X., & Jiao, J. J. (2011). Methods to Derive the Differential Equation of the Free Surface Boundary. Ground Water, 49(2), 133-143.
  More info

  PMID: 21244417;
• Ming, Y. e., Dan, L. u., Neuman, S. P., & Meyer, P. D. (2010). Comment on "inverse groundwater modeling for hydraulic conductivity estimation using Bayesian model averaging and variance window" by Frank T.-C. Tsai and Xiaobao Li. Water Resources Research, 46(2).
• Ming, Y. e., Meyer, P. D., Lin, Y., & Neuman, S. P. (2010). Quantification of model uncertainty in environmental modeling. Stochastic Environmental Research and Risk Assessment, 24(6), 807-808.
• Morales-Casique, E., Neuman, S. P., & Vesselinov, V. V. (2010). Maximum likelihood Bayesian averaging of airflow models in unsaturated fractured tuff using Occam and variance windows. Stochastic Environmental Research and Risk Assessment, 24(6), 863-880.
  More info

  Abstract: We use log permeability and porosity data obtained from single-hole pneumatic packer tests in six boreholes drilled into unsaturated fractured tuff near Superior, Arizona, to postulate, calibrate and compare five alternative variogram models (exponential, exponential with linear drift, power, truncated power based on exponential modes, and truncated power based on Gaussian modes) of these parameters based on four model selection criteria (AIC, AICc, BIC and KIC). Relying primarily on KIC and cross-validation we select the first three of these variogram models and use them to parameterize log air permeability and porosity across the site via kriging in terms of their values at selected pilot points and at some single-hole measurement locations. For each of the three variogram models we estimate log air permeabilities and porosities at the pilot points by calibrating a finite volume pressure simulator against two cross-hole pressure data sets from sixteen boreholes at the site. The traditional Occam's window approach in conjunction with AIC, AICc, BIC and KIC assigns a posterior probability of nearly 1 to the power model. A recently proposed variance window approach does the same when applied in conjunction with AIC, AICc, BIC but spreads the posterior probability more evenly among the three models when used in conjunction with KIC. We compare the abilities of individual models and MLBMA, based on both Occam and variance windows, to predict space-time pressure variations observed during two cross-hole tests other than those employed for calibration. Individual models with the largest posterior probabilities turned out to be the worst or second worst predictors of pressure in both validation cases. Some individual models predicted pressures more accurately than did MLBMA. MLBMA was far superior to any of the individual models in one validation test and second to last in the other validation test in terms of predictive coverage and log scores. © 2010 The Author(s).
• Neuman, S. P. (2010). Apparent/spurious multifractality of data sampled from fractional Brownian/Lévy motions. Hydrological Processes, 24(15), 2056-2067.
  More info

  Abstract: Many earth and environmental variables appear to be self-affine (monofractal) or multifractal with spatial (or temporal) increments having exceedance probability tails that decay as powers of -α where 1 < α ≤ 2. The literature considers selfaffine and multifractal modes of scaling to be fundamentally different, the first arising from additive and the second from multiplicative random fields or processes. We demonstrate theoretically that data having finite support, sampled across a finite domain from one or several realizations of an additive Gaussian field constituting fractional Brownian motion (fBm) characterized by α = 2, give rise to positive square (or absolute) increments which behave as if the field was multifractal when in fact it is monofractal. Sampling such data from additive fractional Lévy motions (fLm) with 1 < α < 2 causes them to exhibit spurious multifractality. Deviations from apparent multifractal behaviour at small and large lags are due to nonzero data support and finite domain size, unrelated to noise or undersampling (the causes cited for such deviations in the literature). Our analysis is based on a formal decomposition of anisotropic fLm (fBm when α = 2) into a continuous hierarchy of statistically independent and homogeneous random fields, or modes, which captures the above behaviour in terms of only E + 3 parameters where E is Euclidean dimension. Although the decomposition is consistent with a hydrologic rationale proposed by Neuman (2003), its mathematical validity is independent of such a rationale. Our results suggest that it may be worth checking how closely would variables considered in the literature to be multifractal (e.g. experimental and simulated turbulent velocities, some simulated porous flow velocities, landscape elevations, rain intensities, river network area and width functions, river flow series, soil water storage and physical properties) fit the simpler monofractal model considered in this paper (such an effort would require paying close attention to the support and sampling window scales of the data). Parsimony would suggest associating variables found to fit both models equally well with the latter. © 2010 John Wiley & Sons, Ltd.
• Kuhlman, K. L., & Neuman, S. P. (2009). Laplace-transform analytic-element method for transient porous-media flow. Journal of Engineering Mathematics, 64(2), 113-130.
  More info

  Abstract: A unified theory of the Laplace-transform analytic-element method (LT-AEM) for solving transient porous-media flow problems is presented. LT-AEM applies the analytic-element method (AEM) to the modified Helmholtz equation, the Laplace-transformed diffusion equation. LT-AEM uses superposition and boundary collocation with Laplace-space convolution to compute flexible semi-analytic solutions from a small collection of fundamental elements. The elements discussed are derived using eigenfunction expansions of element shapes in their natural coordinates. A new formulation for a constant-strength line source is presented in terms of elliptical coordinates and complex-parameter Mathieu functions. Examples are given illustrating how leaky and damped-wave hydrologic problems can be solved with little modification using existing LT-AEM techniques. © Springer Science+Business Media B.V. 2008.
• Morales-Casique, E., & Neuman, S. P. (2009). Laplace-transform finite element solution of nonlocal and localized stochastic moment equations of transport. Communications in Computational Physics, 6(1), 131-161.
  More info

  Abstract: Morales-Casique et al. (Adv. Water Res., 29 (2006), pp. 1238-1255) developed exact first and second nonlocal moment equations for advective-dispersive transport in finite, randomly heterogeneous geologic media. The velocity and concentration in these equations are generally nonstationary due to trends in heterogeneity, conditioning on site data and the influence of forcing terms. Morales-Casique et al. (Adv. Water Res., 29 (2006), pp. 1399-1418) solved the Laplace transformed versions of these equations recursively to second order in the standard deviation σy of (natural) log hydraulic conductivity, and iteratively to higher-order, by finite elements followed by numerical inversion of the Laplace transform. They did the same for a space-localized version of the mean transport equation. Here we recount briefly their theory and algorithms; compare the numerical performance of the Laplace-transform finite element scheme with that of a high-accuracy ULTIMATE-QUICKEST algorithm coupled with an alternating split operator approach; and review some computational results due to Morales-Casique et al. (Adv. Water Res., 29 (2006), pp. 1399-1418) to shed light on the accuracy and computational efficiency of their recursive and iterative solutions in comparison to conditional Monte Carlo simulations in two spatial dimensions. ©2009 Global-Science Press.
• Neuman, S. P. (2009). Comment on 'Siegel D. 2008. Reductionist hydrogeology: Ten fundamental principles. Hydrological processes 22: 4967-4970'. Hydrological Processes, 23(11), 1676-1677.
• Neuman, S. P. (2009). Hortonian scaling and effect of cutoffs on statistics of self-affine river networks. Water Resources Research, 45(12).
  More info

  Abstract: The planar structure of river networks exhibits fractal properties. In particular, available data indicate a tendency for drainage areas A to be distributed according to a power law; their boundaries and main channels to form self-affine curves; the characteristic lengths of drainage areas to be independently self-affine in the one-dimensional space Z of total channel length, rendering the network elongated (and anisotropic); the inverse network density A/Z to be either constant or self-affine in Z, depending on whether or not the network fills the two-dimensional space of A; and (as found in a recent study) areas of given size, vegetative cover, and mean steady state soil moisture, weighted by their flow distance from the basin outlet, to be self-affine in the one-dimensional space of this distance. Various theoretical and semiempirical relationships have been proposed among exponents defining some of these and other scale dependencies. Expressions have been proposed for ensemble moments of A and some length measures associated with finite size river basins. We present a new model that views any self-affine basin property, Y(X), as belonging to an infinite hierarchy of mutually uncorrelated, statistically homogeneous random functions defined on elementary subbasins, each of which is characterized by a unique integral (spatial correlation) scale λ. We cite a mathematically rigorous hydrologic rationale for our model and use it in conjunction with published scaling relations to obtain the probability density function of λ; to derive analytical expressions for ensemble analogs of Horton's scaling laws; to deduce from them that streams of any Horton-Strahler order ω are associated with integral scales λω ≤ λ ≤ λω+1, where the ratio λω+1/λω is a constant independent of ω; to develop analytical expressions relating statistical moments of Y(X) to arbitrary lower and upper cutoff scales that may (but need not) be taken to represent data support and maximum watershed size; to describe ways of estimating the corresponding parameters; and to provide a theoretical basis for the heretofore unexplained observation that self-affine amplitude fluctuations of basin boundaries and main channels (as measured by their variance), having a common Hurst exponent, are larger in the former than in the latter. Copyright 2009 by the American Geophysical Union.
• Neuman, S. P., & Tartakovsky, D. M. (2009). Perspective on theories of non-Fickian transport in heterogeneous media. Advances in Water Resources, 32(5), 670-680.
  More info

  Abstract: Subsurface fluid flow and solute transport take place in a multiscale heterogeneous environment. Neither these phenomena nor their host environment can be observed or described with certainty at all scales and locations of relevance. The resulting ambiguity has led to alternative conceptualizations of flow and transport and multiple ways of addressing their scale and space-time dependencies. We focus our attention on four approaches that give rise to nonlocal representations of advective and dispersive transport of nonreactive tracers in randomly heterogeneous porous or fractured continua. We compare these approaches theoretically on the basis of their underlying premises and the mathematical forms of the corresponding nonlocal advective-dispersive terms. One of the four approaches describes transport at some reference support scale by a classical (Fickian) advection-dispersion equation (ADE) in which velocity is a spatially (and possibly temporally) correlated random field. The randomness of the velocity, which is given by Darcy's law, stems from random fluctuations in hydraulic conductivity (and advective porosity though this is often disregarded). Averaging the stochastic ADE over an ensemble of velocity fields results in a space-time-nonlocal representation of mean advective-dispersive flux, an approach we designate as stnADE. A closely related space-time-nonlocal representation of ensemble mean transport is obtained upon averaging the motion of solute particles through a random velocity field within a Lagrangian framework, an approach we designate stnL. The concept of continuous time random walk (CTRW) yields a representation of advective-dispersive flux that is nonlocal in time but local in space. Closely related to the latter are forms of ADE entailing fractional derivatives (fADE) which leads to representations of advective-dispersive flux that are nonlocal in space but local in time; nonlocality in time arises in the context of multirate mass transfer models, which we exclude from consideration in this paper. We describe briefly each of these four nonlocal approaches and offer a perspective on their differences, commonalities, and relative merits as analytical and predictive tools. © 2008 Elsevier Ltd. All rights reserved.
• Riva, M., Guadagnini, A., Neuman, S. P., Janetti, E. B., & Malama, B. (2009). Inverse analysis of stochastic moment equations for transient flow in randomly heterogeneous media. Advances in Water Resources, 32(10), 1495-1507.
  More info

  Abstract: We present a nonlinear stochastic inverse algorithm that allows conditioning estimates of transient hydraulic heads, fluxes and their associated uncertainty on information about hydraulic conductivity (K) and hydraulic head (h) data collected in a randomly heterogeneous confined aquifer. Our algorithm is based on Laplace-transformed recursive finite-element approximations of exact nonlocal first and second conditional stochastic moment equations of transient flow. It makes it possible to estimate jointly spatial variations in natural log-conductivity (Y = ln K), the parameters of its underlying variogram, and the variance-covariance of these estimates. Log-conductivity is parameterized geostatistically based on measured values at discrete locations and unknown values at discrete "pilot points". Whereas prior values of Y at pilot point are obtained by generalized kriging, posterior estimates at pilot points are obtained through a maximum likelihood fit of computed and measured transient heads. These posterior estimates are then projected onto the computational grid by kriging. Optionally, the maximum likelihood function may include a regularization term reflecting prior information about Y. The relative weight assigned to this term is evaluated separately from other model parameters to avoid bias and instability. We illustrate and explore our algorithm by means of a synthetic example involving a pumping well. We find that whereas Y and h can be reproduced quite well with parameters estimated on the basis of zero-order mean flow equations, all model quality criteria identify the second-order results as being superior to zero-order results. Identifying the weight of the regularization term and variogram parameters can be done with much lesser ambiguity based on second- than on zero-order results. A second-order model is required to compute predictive error variances of hydraulic head (and flux) a posteriori. Conditioning the inversion jointly on conductivity and hydraulic head data results in lesser predictive uncertainty than conditioning on conductivity or head data alone. © 2009 Elsevier Ltd. All rights reserved.
• Meyer, P. D., Ming, Y. e., Neuman, S. P., & Rockhold, M. L. (2008). Application of maximum likelihood Bayesian model averaging to groundwater flow and transport at the Hanford Site 300 area. IAHS-AISH Publication, 64-69.
  More info

  Abstract: A methodology to systematically and quantitatively assess model predictive uncertainty was applied to saturated zone uranium transport at the 300 Area of the US Department of Energy Hanford Site in Washington State, USA. The methodology extends Maximum Likelihood Bayesian Model Averaging (MLBMA) to account jointly for uncertainties due to the conceptual-mathematical basis of models, model parameters, and the scenarios to which the models are applied. Conceptual uncertainty was represented by postulating four alternative models of hydrogeology and uranium adsorption. Parameter uncertainties were represented by estimation covariances resulting from the joint calibration of each model to observed heads and uranium concentration. Posterior model probability was dominated by one model. Results demonstrated the role of model complexity and fidelity to observed system behaviour in determining model probabilities, as well as the impact of prior information. Two scenarios representing alternative future behaviour of the Columbia River adjacent to the site were considered. Predictive simulations carried out with the calibrated models illustrated the computation of model- and scenario-averaged predictions and how results can be displayed to clearly indicate the individual contributions to predictive uncertainty of the model, parameter, and scenario uncertainties. The application demonstrated the practicability of applying a comprehensive uncertainty assessment to large-scale, detailed groundwater flow and transport modelling.
• Ming, Y. e., Meyer, P. D., & Neuman, S. P. (2008). On model selection criteria in multimodel analysis. Water Resources Research, 44(3).
  More info

  Abstract: Hydrologic systems are open and complex, rendering them prone to multiple conceptualizations and mathematical descriptions. There has been a growing tendency to postulate several alternative hydrologic models for a site and use model selection criteria to (1) rank these models, (2) eliminate some of them, and/or (3) weigh and average predictions and statistics generated by multiple models. This has led to some debate among hydrogeologists about the merits and demerits of common model selection (also known as model discrimination or information) criteria such as AIC, AICc, BIC, and KIC and some lack of clarity about the proper interpretation and mathematical representation of each criterion. We examine the model selection literature to find that (1) all published rigorous derivations of AIC and AICc require that the (true) model having generated the observational data be in the set of candidate models; (2) though BIC and KIC were originally derived by assuming that such a model is in the set, BIC has been rederived by Cavanaugh and Neath (1999) without the need for such an assumption; and (3) KIC reduces to BIC as the number of observations becomes large relative to the number of adjustable model parameters, implying that it likewise does not require the existence of a true model in the set of alternatives. We explain why KIC is the only criterion accounting validly for the likelihood of prior parameter estimates, elucidate the unique role that the Fisher information matrix plays in KIC, and demonstrate through an example that it imbues KIC with desirable model selection properties not shared by AIC, AICc, or BIC. Our example appears to provide the first comprehensive test of how AIC, AICc, BIC, and KIC weigh and rank alternative models in light of the models' predictive performance under cross validation with real hydrologic data. Copyright 2008 by the American Geophysical Union.
• Morales-Casique, E., Neuman, S. P., & Vesselinov, V. V. (2008). Maximum likelihood Bayesian averaging of air flow models in unsaturated fractured tuff. IAHS-AISH Publication, 70-75.
  More info

  Abstract: MLBMA is a maximum likelihood (ML) version of Bayesian model averaging (BMA) that renders it compatible with ML methods of model calibration and thus applicable to cases where prior information about the parameter may be unavailable. We explore the role of prior information in MLBMA by applying it to air flow during a cross-hole pneumatic injection test in unsaturated fractured tuff with and without reliance on packer-test data from six boreholes. We parameterize log air permeability and porosity geostatistically using pilot points and estimate them by calibrating a finite volume pressure simulator (FEHM) against cross-hole pressure data by means of a parallelized version of PEST considering several alternative variogram models. We assess the predictive capabilities of each model based on various model selection criteria and discuss future plans to generate corresponding predictions via MLBMA, cross-validate them against pressure data from the same cross-hole test, and validate them against data from another such test. Copyright © 2008 IAHS Press.
• Neuman, S. P. (2008). Multiscale relationships between fracture length, aperture, density and permeability. Geophysical Research Letters, 35(22).
  More info

  Abstract: Fractured rocks exhibit a hierarchical structure which renders their attributes scale-dependent. In particular available data indicate a tendency for fracture length scales to be distributed according to a power law, average fracture aperture to be given by a power of the fracture length scale, and fracture density as well as log permeability to behave as random fractals. To date, no consistent theoretical relationship has been developed between fracture type (as categorized, for example, by length scale and/or aperture) and corresponding fractal attributes (such as density and log permeability). We develop and explore multiscale relationships between these fracture categories and attributes on the basis of a theory recently proposed by the author, which allows linking them in a formal way. Copyright 2008 by the American Geophysical Union.
• Neuman, S. P. (2008). Shlomo p. Neuman: a brief autobiography.. Ground water, 46(1), 164-169.
  More info

  PMID: 18181876;
• Neuman, S. P., Riva, M., & Guadagnini, A. (2008). On the geostatistical characterization of hierarchical media. IAHS-AISH Publication, 3-8.
  More info

  Abstract: The subsurface exhibits hierarchical structures with multiscale hydraulic and transport properties. The structure and properties are often well characterized by stationary variograms. We propose that the latter is an artefact of sampling over finite windows, which disappears when one uses truncated power variograms (TPVs). We cite evidence that parameters of traditional stationary variograms vary with support and window scales; demonstrate the ability of TPVs to capture these scale variations; show that stationary variograms are often difficult to distinguish from TPVs; note that TPVs are unique in their ability to represent multiscale random fields having either Gaussian or heavy-tailed symmetric Levy stable probability distributions; detail the way in which TPVs allow conditioning on multiscale measurements via co-kriging; and illustrate these capabilities on multiscale hydraulic data from Tubingen, Germany. Copyright © 2008 IAHS Press.
• Neuman, S. P., Riva, M., & Guadagnini, A. (2008). On the geostatistical characterization of hierarchical media. Water Resources Research, 44(2).
  More info

  Abstract: The subsurface consists of porous and fractured materials exhibiting a hierarchical geologic structure, which gives rise to systematic and random spatial and directional variations in hydraulic and transport properties on a multiplicity of scales. Traditional geostatistical moment analysis allows one to infer the spatial covariance structure of such hierarchical, multiscale geologic materials on the basis of numerous measurements on a given support scale across a domain or "window" of a given length scale. The resultant sample variogram often appears to fit a stationary variogram model with constant variance (sill) and integral (spatial correlation) scale. In fact, some authors, who recognize that hierarchical sedimentary architecture and associated log hydraulic conductivity fields tend to be nonstationary, nevertheless associate them with stationary "exponential-like" transition probabilities and variograms, respectively, the latter being a consequence of the former. We propose that (1) the apparent ability of stationary spatial statistics to characterize the covariance structure of nonstationary hierarchical media is an artifact stemming from the finite size of the windows within which geologic and hydrologic variables are ubiquitously sampled, and (2) the artifact is eliminated upon characterizing the covariance structure of such media with the aid of truncated power variograms, which represent stationary random fields obtained upon sampling a nonstationary fractal over finite windows. To support our opinion, we note that truncated power variograms arise formally when a hierarchical medium is sampled jointly across all geologic categories and scales within a window; cite direct evidence that geostatistical parameters (variance and integral scale) inferred on the basis of traditional variograms vary systematically with support and window scales; demonstrate the ability of truncated power models to capture these variations in terms of a few scaling parameters; show that exponential and truncated power variograms are often difficult to distinguish from each other, which helps explain why hierarchical data may appear to fit the former; note that truncated power models are unique in their ability to represent multiscale random fields having either Gaussian or heavy-tailed symmetric Levy stable probability distributions; detail the way in which these models allow conditioning the spatial statistics of such fields on multiscale measurements via cokriging; and illustrate these capabilities on multiscale hydraulic data from an unconfined aquifer near Tübingen, Germany. Copyright 2008 by the American Geophysical Union.
• Riva, M., Guadagnini, A., & Neuman, S. P. (2008). Characterizing the spatial variability of transmissivity using stochastic type-curve and numerical inverse analyses of data from a sequence of pumping tests. IAHS-AISH Publication, 39-44.
  More info

  Abstract: We discuss two recent methods of characterizing the spatial variability of a random (natural) log transmissivity field on the basis of observed space-time variations in hydraulic head: a graphical stochastic type-curve method and a geostatistical method of inverting (ensemble) mean flow equations. While both methods allow estimating the unconditional variance and integral (correlation) scale of log transmissivities, geostatistical inversion is computationally more intensive, but also provides tomographic images of how log transmissivity estimates and their variance vary in space. We apply the two approaches to synthetic scenarios and to measured late time (quasi-steady state) drawdowns from a sequence of transient pumping tests in an unconfined aquifer near Tübingen, Germany. Copyright © 2008 IAHS Press.
• Tartakovsky, A. M., & Neuman, S. P. (2008). Effects of Peclet number on pore-scale mixing and channeling of a tracer and on directional advective porosity. Geophysical Research Letters, 35(21).
  More info

  Abstract: Field tracer tests indicate that advective porosity, the quantity relating advective velocity to Darcy flux, may exhibit directional dependence. Theory (S. P. Neuman, 2005) suggests that directional variations in advective porosity could arise from incomplete mixing of an inert tracer between directional flow channels within a sampling volume of soil or rock that may be hydraulically isotropic or anisotropic, rendering advective porosity a function of Peclet number. Whereas Neuman was able to quantify the process for idealized intersecting fracture sets, his theory is too abstract to allow doing so for more complex systems of fractures or pores. We use Smoothed Particle Hydrodynamics (e.g., A. M. Tartakovsky and P. Meakin, 2005, 2006) to simulate advective-diffusive transport of a passive tracer in the interior of an idealized pore space between impermeable circular or elliptical solid grains of random size, arranged randomly within a square domain. Our numerical experiments reveal a transition from complete mixing of tracer at small Peclet numbers to the gradual evolution of pore-scale channels of tracer which allow less and less mixing as Peclet numbers grow. Advective porosity, computed as the ratio between total water flux across the domain boundaries and migration rate of the tracer center of mass, is found to diminish asymptotically from a value close to the total porosity at small Peclet numbers to a constant dominated entirely by advection at large Peclet numbers. When total porosity is kept constant, channeling is less pronounced and mixing more prominent when flow takes place parallel to the longer axes of elliptical grains as compared to when the grains are circular. On the other hand, channeling is more pronounced and mixing less so when flow takes place perpendicular to the longer elliptical grain axes. Correspondingly, the rate at which advective porosity diminishes with increasing Peclet number is slower and the asymptote is larger when flow takes place parallel to the longer axes of elliptical grains than when the grains are circular. However, the rate is faster and the asymptote smaller when flow takes place perpendicular to the longer elliptical grain axes. Copyright 2008 by the American Geophysical Union.
• Neuman, S. P., Blattstein, A., Riva, M., Tartakovsky, D. M., Guadagnini, A., & Ptak, T. (2007). Type curve interpretation of late-time pumping test data in randomly heterogeneous aquifers. Water Resources Research, 43(10).
  More info

  Abstract: The properties of heterogeneous media vary spatially in a manner that can seldom be described with certainty. It may, however, be possible to describe the spatial variability of these properties in terms of geostatistical parameters such as mean, integral (spatial correlation) scale, and variance. Neuman et al. (2004) proposed a graphical method to estimate the geostatistical parameters of (natural) log transmissivity on the basis of quasi-steady state head data when a randomly heterogeneous confined aquifer is pumped at a constant rate from a fully penetrating well. They conjectured that a quasi-steady state, during which heads vary in space-time while gradients vary only in space, develops in a statistically homogeneous and horizontally isotropic aquifer as it does in a uniform aquifer. We confirm their conjecture numerically for Gaussian log transmissivities, show that time-drawdown data from randomly heterogeneous aquifers are difficult to interpret graphically, and demonstrate that quasi-steady state distance-drawdown data are amenable to such interpretation by the type curve method of Neuman et al. The method yields acceptable estimates of statistical log transmissivity parameters for fields having either an exponential or a Gaussian spatial correlation function. These estimates are more robust than those obtained using the graphical time-drawdown method of Copty and Findikakis (2003, 2004a). We apply the method of Neuman et al. (2004) simultaneously to data from a sequence of pumping tests conducted in four wells in an aquifer near Tübingen, Germany, and compare our transmissivity estimate with estimates obtained from 312 flowmeter measurements of hydraulic conductivity in these and eight additional wells at the site. We find that (1) four wells are enough to provide reasonable estimates of lead log transmissivity statistics for the Tübingen site using this method, and (2) the time-drawdown method of Cooper and Jacob (1946) underestimates the geometric mean transmissivity at the site by 30-40%. Copyright 2007 by the American Geophysical Union.
• Tartakovsky, G. D., & Neuman, S. P. (2007). Three-dimensional saturated-unsaturated flow with axial symmetry to a partially penetrating well in a compressible unconfined aquifer. Water Resources Research, 43(1).
  More info

  Abstract: A new analytical solution is presented for the delayed response process characterizing flow to a partially penetrating well in an unconfined aquifer. The new solution generalizes that of Neuman (1972, 1974) by accounting for unsaturated flow above the water table. Three-dimensional, axially symmetric flow in the unsaturated zone is described by a linearized version of Richards' equation in which hydraulic conductivity and water content vary exponentially with incremental capillary pressure head relative to its air entry value (defining the interface between the saturated and unsaturated zones). Unsaturated soil properties are characterized by an exponent κ having the dimension of inverse length or, equivalently, a dimensionless exponent κD = κb, where b is initial saturated thickness. Our treatment of the unsaturated zone is similar to that of Kroszynski and Dagan (1975), who, however, have ignored internal (artesian) aquifer storage. According to Kroszynski and Dagan, aquifers that are not excessively shallow have values of κD (their parameter a) much greater than 10. We find that in such typical cases, unsaturated flow has little impact on early and late dimensionless time drawdown a short distance below the water table. Unsaturated flow causes drawdown to increase slightly at intermediate dimensionless time values that represent transition from an early artesian-dominated to a late water-table-dominated flow regime. Delayed drainage from the unsaturated zone becomes less and less important as κD increases; as κD → ∞, this effect dies out, and drawdown is controlled entirely by delayed decline in the water table as in the model of Neuman. The unsaturated zone has a major impact on drawdown at intermediate time and a significant impact at early and late times, in the atypical case of κD ≤ 1, becoming the dominant factor as κD approaches zero (the soil water retention capacity becomes very large and/or saturated thickness becomes insignificant). Our new solution was used to analyze field data from a pumping test conducted by Moench et al. (2001) in a glacial outwash deposit at Cape Cod, Massachusetts. The solution was fitted individually and simultaneously to time-drawdown data from 20 piezometers and observation wells and simultaneously to data from three piezometers in each of two clusters at various depths and distances from the pumping well, with very good results. Our parameter estimates of hydraulic conductivities from the simultaneous fit are similar to those obtained previously by Moench (2004), but our estimates of specific yield and storage are smaller and larger, respectively, while our estimate of κ is not comparable with his estimates of three empirical parameters. Copyright 2007 by the American Geophysical Union.
• Yanhua, L. i., & Neuman, S. P. (2007). Flow to a well in a five-layer system with application to the Oxnard Basin. Ground Water, 45(6), 672-682.
  More info

  PMID: 17973745;Abstract: Nearly 40 years ago, Neuman (1968) developed an analytical solution for drawdown in a system of three aquifers separated by two aquitards when one of the aquifers is pumped at a constant rate. Whereas the simpler case of two aquifers separated by one aquitard has been presented by Neuman and Witherspoon (1969a), the full five-layer solution has not been previously evaluated. We do so here using numerical inversion of its Laplace-transformed version and present selected results graphically in dimensionless form. The solution demonstrates that the effect of pumping propagates across all five layers, adding emphasis to a question previously raised by Neuman and Witherspoon about the validity of leaky aquifer theories that disregard drawdowns in unpumped aquifers. A large-scale, long-term pumping test spanning three aquifers separated by two aquitards near Oxnard, California, has been conducted and analyzed by Neuman and Witherspoon (1972). They evaluated the vertical hydraulic diffusivities of the aquitards using the Neuman-Witherspoon ratio method and their specific storage values on the basis of laboratory consolidation tests. We reinterpret the Oxnard pumping test by coupling the five-layer analytical solution of Neuman (1968) with the parameter estimation code PEST (Doherty 2002) and validate our results against drawdowns from a subsequent pumping test at the site. Our parameter estimates compare favorably with those of Neuman and Witherspoon (1972). © 2007 National Ground Water Association.
• Hernandez, A. F., Neuman, S. P., Guadagnini, A., & Carrera, J. (2006). Inverse stochastic moment analysis of steady state flow in randomly heterogeneous media. Water Resources Research, 42(5).
  More info

  Abstract: Nonlocal stochastic moment equations have been used successfully to analyze steady state and transient flow in randomly heterogeneous media conditional on measured values of medium properties. We present a nonlinear geostatistical inverse algorithm for steady state flow that makes it possible to further condition such analyses on measured values of state variables. Our approach accounts for all scales of spatial variability resolvable by the computational grid. It is based on recursive finite element approximations of exact nonlocal first and second conditional moment equations. Hydraulic conductivity is parameterized geostatistically on the basis of measured values at discrete locations and unknown values at discrete "pilot points." Prior estimates of pilot point values are obtained (optionally) by generalized kriging. Posterior estimates at pilot points and (optionally) at measurement points are obtained by calibrating mean flow against measured values of head. The parameters are projected onto a computational grid via kriging. Maximum likelihood calibration allows estimating not only hydraulic but also (optionally) unknown variogram parameters with or without prior information about the former. The approach yields covariance matrices for parameter estimation as well as head and flux prediction errors, the latter being obtained a posteriori. We implement our inverse approach on highly and mildly heterogeneous media under superimposed mean uniform and convergent flows in a bounded two-dimensional domain. Our examples illustrate how conductivity and head data act separately and jointly to reduce parameter estimation errors and to model predictive uncertainty. We also evaluate the functional form of the log conductivity variogram and its parameters using likelihood-based model discrimination criteria. Copyright 2006 by the American Geophysical Union.
• Morales-Casique, E., Neuman, S. P., & Guadagnini, A. (2006). Non-local and localized analyses of non-reactive solute transport in bounded randomly heterogeneous porous media: Theoretical framework. Advances in Water Resources, 29(8), 1238-1255.
  More info

  Abstract: Solute transport in randomly heterogeneous media is described by stochastic transport equations that are typically solved via Monte Carlo simulation. A promising alternative is to solve a corresponding system of statistical moment equations directly. We present exact (though not closed) implicit conditional first and second moment equations for advective-dispersive transport in finite domains. The velocity and concentration are generally non-stationary due to possible trends in heterogeneity, conditioning on data, temporal variations in velocity, fluid and/or solute sources, initial and boundary conditions. Our equations are integro-differential and include non-local parameters depending on more than one point in space-time. To allow solving these equations, we close them by perturbation and develop recursive moment equations in Laplace space for the special case of steady state flow, to second order in σY where σY2 is a measure of (natural) log hydraulic conductivity variance. We also propose a higher-order iterative closure. Our recursive equations and iterative closure are formally valid for mildly heterogeneous media, or well-conditioned strongly heterogeneous media in which the random component of heterogeneity is relatively small. The non-local moment equations suggest (and a companion paper [Morales Casique E, Neuman SP, Guadagnini A. Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: computational analysis. Adv Water Resour, submitted for publication] demonstrates numerically) that, in general, transport cannot be validly described by means of Fick's law with a (constant or variable) macrodispersion coefficient. We show nevertheless that, under a limited set of conditions, the mean transport equation can be localized to yield a familiar-looking advection-dispersion equation with a conditional macrodispersion tensor that varies generally in space-time. In a companion paper [Morales Casique E, Neuman SP, Guadagnini A. Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: computational analysis. Adv Water Resour, submitted for publication] we present a high-accuracy computational algorithm for our iterative non-local and recursive localized moment equations, assessing their accuracy and computational efficiency in comparison to unconditional and conditional Monte Carlo simulations. © 2005 Elsevier Ltd. All rights reserved.
• Morales-Casique, E., Neuman, S. P., & Guadagnini, A. (2006). Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: Computational analysis. Advances in Water Resources, 29(9), 1399-1418.
  More info

  Abstract: In a companion paper we presented exact (though implicit and not closed) nonlocal conditional first and second moment equations for nonreactive advective-dispersive transport under both steady state and transient flow regimes in bounded, randomly heterogeneous porous domains. To allow solving our nonlocal equations we developed recursive moment equations in Laplace space for the special case of steady state flow to second order in σY (a measure of the standard deviation of natural log hydraulic conductivity Y, which is generally nonhomogeneous) and proposed a higher-order iterative closure scheme. We also showed that, under a limited set of conditions, the mean transport equation can be localized to yield a familiar-looking advection-dispersion equation with a conditional macrodispersion tensor that generally varies in space-time. The purpose of this paper is to explore the behavior and assess the accuracy and computational efficiency of our moment solutions in comparison to conditional and unconditional Monte Carlo simulations. To do so, we present a high-accuracy computational algorithm for our iterative nonlocal and recursive localized moment equations and corresponding computational results in two spatial dimensions conditional on measurements of Y. Our algorithm solves the moment equations by finite elements in Laplace-transformed space and inverts the solution numerically back into the time domain. Conditional results obtained with our iterative algorithm compare well with Monte Carlo simulations for σY2 = 0.3 and Peclet number Pe = 100 defined in terms of the integral scale of Y, and for Pe = 10 in the unconditional case. As σY2, Pe and time increase the quality of our iterative moment solution deteriorates. We show that this is due to our disregarding velocity moments of order higher than two and propose that including such moments should render our iterative solution workable over a wider range of these parameters. Second-order recursive nonlocal and space-localized results are considerably less accurate than those obtained with our iterative nonlocal algorithm. Even though our moment solution does not require computing (space-time localized) macrodispersion coefficients, we nevertheless do so to examine the influence of boundaries and conditioning on their behavior. Our results support an earlier observation by the second author [Morales Casique E, Neuman SP, Guadagnini A. Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: theoretical framework. Adv Water Resour, in press., Neuman SP. Universal scaling of hydraulic conductivities and dispersivities in geologic media. Water Resour Res 1990;26(8):1749-58], based on world-wide tracer test results, that the rate at which apparent longitudinal dispersivity increases with scale diminishes with conditioning. In preliminary runs conducted on a relatively small grid without optimizing our algorithms and without parallelization, the moment solutions required considerably less computer time than did the Monte Carlo simulations. © 2005 Elsevier Ltd. All rights reserved.
• Neuman, S. P. (2006). Blueprint for perturbative solution of flow and transport in strongly heterogeneous composite media using fractal and variational multiscale decomposition. Water Resources Research, 42(SUPPL.).
  More info

  Abstract: This paper offers a blueprint for extending perturbative ensemble moment (and potentially other, e.g., Gaussian or perturbative Karhunen-Loéve) solutions of stochastic flow and transport problems in composite, randomly heterogeneous porous media to cases in which medium heterogeneity is arbitrarily large. The proposed approach is based on a novel idea of combined fractal and variational multiscale decomposition. It considers log hydraulic conductivity to be a composite random field with piecewise statistically homogeneous increments characterized by a truncated power variogram, i.e., a truncated fractal. This consideration is supported either exactly or approximately by a growing amount of observational data. It allows decomposing the field into two or more mutually uncorrelated components having relatively small variances associated with a hierarchy of spatial correlation scales. Combining this fractal decomposition with a variational multiscale decomposition formalism, motivated by the multiscale nature of turbulence, allows decomposing flow and transport problems into two or more coupled sets of variational problems each of which is amenable to perturbative solution regardless of how heterogeneous the original field was. The problems are defined on a hierarchy of computational grids having discretization intervals proportional to the correlation scales. The proposed approach does not require introducing effective parameters on any scale; instead, it allows resolving flow and transport ensemble moments fully on multiple scales of spatial resolution. The approach nevertheless allows subgrid closure (another idea inspired by the literature on turbulence) through subgrid refinement or the use of effective subgrid parameters. It additionally offers the critically important capability of conditioning flow and transport analyses on multiscale measurements. Details are presented in the context of steady state single-phase flow in a composite medium. Copyright 2006 by the American Geophysical Union.
• Neuman, S. P. (2006). Comment on "Spatial correlation of permeability in cross-stratified sediment with hierarchical architecture" by Robert W. Ritzi, Zhenxue Dai, David F. Dominic, and Yoram N. Rubin. Water Resources Research, 42(5).
• Neuman, S. P. (2006). Longitudinal dispersivity data and implications for scaling behavior: Comment. Ground Water, 44(2), 139-141.
  More info

  PMID: 16556193;
• Vrugt, J. A., & Neuman, S. P. (2006). Introduction to the special section in Vadose Zone Journal: Parameter identification and uncertainty assessment in the unsaturated zone. Vadose Zone Journal, 5(3), 915-916.
• Cintoli, S., Neuman, S. P., & Federico, V. D. (2005). Generating and scaling fractional Brownian motion on finite domains. Geophysical Research Letters, 32(8), 1-4.
  More info

  Abstract: Power variograms of statistically isotropic or anisotropic fractal fields (common in earth science) are weighted integrals of variograms representing statistically homogeneous fields (modes) having mutually uncorrelated increments. Large- and small-scale cutoffs were previously assumed proportional to length scales of the sampling window and data support. We verify this assumption numerically for two-dimensional isotropic fractional Brownian motion (fBm). It was previously concluded semi-empirically that, for Hurst coefficient H = 0.25, the constant of proportionality is μ = 1/3. We confirm this but find μ to vary with mode type and H. We find that due to lack of ergodicity, sample fBm variograms generated on finite windows exhibit directional dependence and differ sharply between realizations. Many realizations are required to obtain an average sample variogram resembling the theoretical power model, especially for persistent fields. We propose generating fBm on finite windows using truncated power variograms and provide guidance for doing so effectively. Copyright 2005 by the American Geophysical Union.
• Ming, Y. e., Neuman, S. P., Meyer, P. D., & Pohlmann, K. (2005). Sensitivity analysis and assessment of prior model probabilities in MLBMA with application to unsaturated fractured tuff. Water Resources Research, 41(12), 1-14.
  More info

  Abstract: [1] Previous application of maximum likelihood Bayesian model averaging (MLBMA, Neuman (2002, 2003)) to alternative variogram models of log air permeability data in fractured tuff has demonstrated its effectiveness in quantifying conceptual model uncertainty and enhancing predictive capability (Ye et al., 2004). A question remained how best to ascribe prior probabilities to competing models. In this paper we examine the extent to which lead statistics of posterior log permeability predictions are sensitive to prior probabilities of seven corresponding variogram models. We then explore the feasibility of quantifying prior model probabilities by (1) maximizing Shannon's entropy H (Shannon, 1948) subject to constraints reflecting a single analyst's (or a group of analysts') prior perception about how plausible each alternative model (or a group of models) is relative to others, and (2) selecting a posteriori the most likely among such maxima corresponding to alternative prior perceptions of various analysts or groups of analysts. Another way to select among alternative prior model probability sets, which, however, is not guaranteed to yield optimum predictive performance (though it did so in our example) and would therefore not be our preferred option, is a minimum-maximum approach according to which one selects a priori the set corresponding to the smallest value of maximum entropy. Whereas maximizing H subject to the prior perception of a single analyst (or group) maximizes the potential for further information gain through conditioning, selecting the smallest among such maxima gives preference to the most informed prior perception among those of several analysts (or groups). We use the same variogram models and log permeability data as Ye et al. (2004) to demonstrate that our proposed approach yields the least amount of posterior entropy (residual uncertainty after conditioning) and enhances predictive model performance as compared to (1) the noninformative neutral case in which all prior model probabilities are set equal to each other and (2) an informed case that nevertheless violates the principle of parsimony. Copyright 2005 by the American Geophysical Union.
• Neuman, S. P. (2005). Comment on the "Nothing older than three years".. Ground water, 43(3), 294-.
  More info

  PMID: 15882319;
• Neuman, S. P. (2005). Erratum: "On the tensorial nature of advective porosity" (Advances Water Resources (2005) vol. 28 (2) (149-159) 10.1016/j.advwatres.2004.09.008). Advances in Water Resources, 28(8), 884-.
• Neuman, S. P. (2005). On the tensorial nature of advective porosity. Advances in Water Resources, 28(2), 149-159.
  More info

  Abstract: Field tracer tests indicate that advective porosity, the quantity relating advective velocity to Darcy flux, may exhibit directional dependence. Hydraulic anisotropy explains some but not all of the reported directional results. The present paper shows mathematically that directional variations in advective porosity may arise simply from incomplete mixing of an inert tracer between directional flow channels within a sampling (or support) volume ψ of soil or rock that may be hydraulically isotropic or anisotropic. In the traditional fully homogenized case, our theory yields trivially a scalar advective porosity equal to the interconnected porosity φ, thus explaining neither the observed directional effects nor the widely reported experimental finding that advective porosity is generally smaller than φ. We consider incomplete mixing under conditions in which the characteristic time tD of longitudinal diffusion along channels across ω is much shorter than the characteristic time tH required for homogenization through transverse diffusion between channels. This may happen where flow takes place preferentially through relatively conductive channels and/or fractures of variable orientation separated by material that forms a partial barrier to diffusive transport. Our solution is valid for arbitrary channel Peclet numbers on a correspondingly wide range of time scales tD ≤ t ≪ tH. It shows that the tracer center of mass is advected at a macroscopic velocity which is generally not collinear with the macroscopic Darcy flux and exceeds it in magnitude. These two vectors are related through a second-rank symmetric advective dispersivity tensor Φ. If the permeability k of ω is a symmetric positive-definite tensor, so is Φ. However, the principal directions and values of these two tensors are generally not the same; whereas those of k are a fixed property of the medium and the length-scale of ω, those of Φ depend additionally on the direction and magnitude of the applied hydraulic gradient. When the latter is large, diffusion has negligible effect on Φ and one may consider tracer mass to be distributed between channels in proportion to the magnitude of their Darcy flux. This is made intuitive through a simple example of an idealized fracture network. Our analytical formalism reveals the properties of Φ but is too schematic to allow predicting the latter accurately on the basis of realistic details about the void structure of ω and tracer mass distribution within it. Yet knowing the tensorial properties of Φ is sufficient to allow determining it indirectly on the basis of ω-scale hydraulic and tracer data, including concentrations that represent homogenized samples extracted from (or sensed externally across) an ω-scale plume. © 2004 Elsevier Ltd. All rights reserved.
• Neuman, S. P. (2005). Trends, prospects and challenges in quantifying flow and transport through fractured rocks. Hydrogeology Journal, 13(1), 124-147.
  More info

  Abstract: Among the current problems that hydrogeologists face, perhaps there is none as challenging as the characterization of fractured rock (Faybishenko and Benson 2000). This paper discusses issues associated with the quantification of flow and transport through fractured rocks on scales not exceeding those typically associated with single- and multi-well pressure (or flow) and tracer tests. As much of the corresponding literature has focused on fractured crystalline rocks and hard sedimentary rocks such as sandstones, limestones (karst is excluded) and chalk, so by default does this paper. Direct quantification of flow and transport in such rocks is commonly done on the basis of fracture geometric data coupled with pressure (or flow) and tracer tests, which therefore form the main focus. Geological, geophysical and geochemical (including isotope) data are critical for the qualitative conceptualization of flow and transport in fractured rocks, and are being gradually incorporated in quantitative flow and transport models, in ways that this paper unfortunately cannot describe but in passing. The hydrogeology of fractured aquifers and other earth science aspects of fractured rock hydrology merit separate treatments. All evidence suggests that rarely can one model flow and transport in a fractured rock consistently by treating it as a uniform or mildly nonuniform isotropic continuum. Instead, one must generally account for the highly erratic heterogeneity, directional dependence, dual or multicomponent nature and multiscale behavior of fractured rocks. One way is to depict the rock as a network of discrete fractures (with permeable or impermeable matrix blocks) and another as a nonuniform (single, dual or multiple) continuum. A third way is to combine these into a hybrid model of a nonuniform continuum containing a relatively small number of discrete dominant features. In either case the description can be deterministic or stochastic. The paper contains a brief assessment of these trends in light of recent experimental and theoretical findings, ending with a short list of prospects and challenges for the future. © Springer-Verlag 2005.
• Amir, O., & Neuman, S. P. (2004). Gaussian Closure of Transient Unsaturated Flow in Random Soils. Transport in Porous Media, 54(1), 55-77.
  More info

  Abstract: The Gaussian closure approximation, previously used by the authors to solve steady state stochastic unsaturated flow problems in randomly heterogeneous soils, is extended here to transient flow. The method avoids linearizing the governing flow equations or the soil constitutive relations. It places no theoretical limit on the variance of constitutive parameters and applies to a broad class of soils with flow properties that scale according to a linearly separable model. Closure is obtained by treating the dimensionless pressure head ψ as a multivariate Gaussian function. It yields a system of coupled nonlinear differential equations for the first and second moments of ψ. We apply the Gaussian closure technique to the problem of transient infiltration into a randomly stratified soil. In each layer, hydraulic conductivity and water content vary exponentially with ψ. Elsewhere we show that application of the technique to other constitutive relations is straightforward. Our solution for the mean and variance of ψ in a one-dimensional layer with random conductivity compares well with Monte Carlo results over a wide range of parameters, provided that the spatial variability of the constitutive exponent is small. The solution provides considerable insight into the behavior of the transient unsaturated stochastic flow problem.
• Furman, A., & Neuman, S. P. (2004). Laplace transform analytic element method for transient flow problems. Developments in Water Science, 55(PART 1), 371-379.
  More info

  Abstract: Recently we [12] proposed a new Laplace transform analytic element method (LT-AEM) for the solution of transient groundwater flow problems. Laplace transformation of the transient flow equations leads to a time-independent modified Helmholtz equation and associated conditions at discontinuities. The latter are solved by AEM and the results transformed numerically back into the time domain. Though neither continuity of potential nor continuity of flux are satisfied automatically at internal discontinuities, both are satisfied approximately via least squares at an overdetermined system of control points in a manner similar to that of Janković [17] and Barnes and Janković [5]. LT-AEM preserves all advantages of the AEM in Laplace space, most importantly its mathematical elegance and grid-free nature. Solution in Laplace space and numerical back transformation into the time domain are done independently for any given time and are thus amenable to parallel computation on multiple processors. This renders the method particularly well suited for cases where a high-accuracy solution is required at a relatively small number of discrete space-time locations. LT-AEM requires a new family of analytic elements associated transformation of known analytic solutions in the time domain, or developed directly in the Laplace domain. We use both methods to develop a number of analytic elements for LT-AEM. We then illustrate the method on transient flow in a two-dimensional confined aquifer containing various inhomogeneities and time-dependent sources. © 2004 Elsevier B.V.
• Ming, Y. e., Neuman, S. P., & Meyer, P. D. (2004). Maximum likelihood Bayesian averaging of spatial variability models in unsaturated fractured tuff. Water Resources Research, 40(5), W051131-W0511317.
  More info

  Abstract: Hydrologic analyses typically rely on a single conceptual-mathematical model. Yet hydrologic environments are open and complex, rendering them prone to multiple interpretations and mathematical descriptions. Adopting only one of these may lead to statistical bias and underestimation of uncertainty. Bayesian model averaging (BMA) [Hoeting et al., 1999] provides an optimal way to combine the predictions of several competing models and to assess their joint predictive uncertainty. However, it tends to be computationally demanding and relies heavily on prior information about model parameters. Neuman [2002, 2003] proposed a maximum likelihood version (MLBMA) of BMA to render it computationally feasible and to allow dealing with cases where reliable prior information is lacking. We apply MLBMA to seven alternative variogram models of log air permeability data from single-hole pneumatic injection tests in six boreholes at the Apache Leap Research Site (ALRS) in central Arizona. Unbiased ML estimates of variogram and drift parameters are obtained using adjoint state maximum likelihood cross validation [Samper and Neuman, 1989a] in conjunction with universal kriging and generalized least squares. Standard information criteria provide an ambiguous ranking of the models, which does not justify selecting one of them and discarding all others as is commonly done in practice. Instead, we eliminate some of the models based on their negligibly small posterior probabilities and use the rest to project the measured log permeabilities by kriging onto a rock volume containing the six boreholes. We then average these four projections and associated kriging variances, using the posterior probability of each model as weight. Finally, we cross validate the results by eliminating from consideration all data from one borehole at a time, repeating the above process and comparing the predictive capability of MLBMA with that of each individual model. We find that MLBMA is superior to any individual geostatistical model of log permeability among those we consider at the ALRS.
• Ming, Y. e., Neuman, S. P., Guadagnini, A., & Tartakovsky, D. M. (2004). Nonlocal and localized analyses of conditional mean transient flow in bounded, randomly heterogeneous porous media. Water Resources Research, 40(5), W051041-W0510419.
  More info

  Abstract: We consider the numerical prediction of transient flow in bounded, randomly heterogeneous porous media driven by random sources, initial heads, and boundary conditions without resorting to Monte Carlo simulation. After applying the Laplace transform to the governing stochastic flow equations, we derive exact nonlocal (integro-differential) equations for the mean and variance-covariance of transformed head and flux, conditioned on measured values of log conductivity Y = 1n K. Approximating these conditional moment equations recursively to second order in the standard deviation σY of Y, we solve them by finite elements for superimposed mean-uniform and convergent flows in a two-dimensional domain. An alternative conditional mean solution is obtained through localization of the exact moment expressions. The nonlocal and localized solutions are obtained using a highly efficient parallel algorithm and inverted numerically back into the time domain. A comparison with Monte Carlo simulations demonstrates that the moment solutions are remarkably accurate for strongly heterogeneous media with σY2 as large as 4. The nonlocal solution is only slightly more accurate than the much simpler localized solution, but the latter does not yield information about predictive uncertainty. The accuracy of each solution improves markedly with conditioning. A preliminary comparison of computational efficiency suggests that both the nonlocal and localized solutions for mean head and its variance require significantly less computer time than is required for Monte Carlo statistics to stabilize when the same direct matrix solver is used for all three (we do not presently know how using iterative solvers would have affected this conclusion). This is true whether the Laplace inversion and Monte Carlo simulations are conducted sequentially or in parallel on multiple processors and regardless of problem size. The underlying exact and recursive moment equations, as well as the proposed computational algorithm, are valid in both two and three dimensions; only the numerical implementation of our algorithm is two-dimensional.
• Neuman, S. P., & Bredehoft, J. D. (2004). "From Models to Performance Assessment: The Conceptualization Problem," by John D. Bredehoft, September-October 2003 issue, v. 41, no. 5: 571-577. Ground Water, 42(5), 785-787.
  More info

  PMID: 15457802;
• Neuman, S. P., Guadagnini, A., & Riva, M. (2004). Type-curve estimation of statistical heterogeneity. Water Resources Research, 40(4), W042011-W042017.
  More info

  Abstract: The analysis of pumping tests has traditionally relied on analytical solutions of groundwater flow equations in relatively simple domains, consisting of one or at most a few units having uniform hydraulic properties. Recently, attention has been shifting toward methods and solutions that would allow one to characterize subsurface heterogeneities in greater detail. On one hand, geostatistical inverse methods are being used to assess the spatial variability of parameters, such as permeability and porosity, on the basis of multiple cross-hole pressure interference tests. On the other hand, analytical solutions are being developed to describe the mean and variance (first and second statistical moments) of flow to a well in a randomly heterogeneous medium. We explore numerically the feasibility of using a simple graphical approach (without numerical inversion) to estimate the geometric mean, integral scale, and variance of local log transmissivity on the basis of quasi steady state head data when a randomly heterogeneous confined aquifer is pumped at a constant rate. By local log transmissivity we mean a function varying randomly over horizontal distances that are small in comparison with a characteristic spacing between pumping and observation wells during a test. Experimental evidence and hydrogeologic scaling theory suggest that such a function would tend to exhibit an integral scale well below the maximum well spacing. This is in contrast to equivalent transmissivities derived from pumping tests by treating the aquifer as being locally uniform (on the scale of each test), which tend to exhibit regional-scale spatial correlations. We show that whereas the mean and integral scale of local log transmissivity can be estimated reasonably well based on theoretical ensemble mean variations of head and drawdown with radial distance from a pumping well, estimating the log transmissivity variance is more difficult. We obtain reasonable estimates of the latter based on theoretical variation of the standard deviation of circumferentially averaged drawdown about its mean.
• Yao, T., Wierenga, P. J., Graham, A. R., & Neuman, S. P. (2004). Neutron probe calibration in a vertically stratified vadose zone. Vadose Zone Journal, 3(4), 1400-1406.
  More info

  Abstract: A procedure is described for neutron probe calibration in deep multilayer vadose zones. Calibration equations relating neutron count ratios (CR) to soil water content were developed for the upper 2.5 m of layered soil profiles using soil texture, water content, and neutron probe data. These equations were extended down to a depth of 10 m by relying on neutron probe data to delineate soil texture zones at depth. Data from two constant-flux field infiltration experiments were used to verify the calibration procedure. The water balance for each of nine soil profiles within the infiltrated area was computed using up to four separate calibration equations for each soil profile. The use of two or more texture-based calibration equations greatly improved the agreement between water applied at the surface and water measured within each of the nine layered profiles. Best results were obtained using four such equations, three corresponding to distinct soil horizons and one to interfaces between them (soil texture boundaries). The calibration procedure proposed here is applicable to deep soil profiles when soil texture and bulk density data are only available for the upper soil profile. © Soil Science Society of America.
• Ando, K., Kostner, A., & Neuman, S. P. (2003). Stochastic continuum modeling of flow and transport in a crystalline rock mass: Fanay-Augères, France, revisited. Hydrogeology Journal, 11(5), 521-535.
  More info

  Abstract: A stochastic discrete-fracture model was used by Cacas et al,a,b to interpret flow measurements and transport experiments in a fractured crystalline rock mass at Fanay-Augères. They considered continuum models to be incapable of properly interpreting small-scale measurements or tracer tests in fractured systems, which, in their view, require three-dimensional modeling of numerous discrete channels; in their opinion, continuum modeling applies only to average flow on a relatively large scale. Cacas et al. considered their discrete fracture model to have been validated by its demonstrated ability to reproduce selected experimental results. In this paper, flow and transport at Fanay-Augères are modeled by viewing the fractured rock as a stochastic continuum in a manner originally proposed by Neumanc,d. The stochastic continuum approach obviates the need for detailed information about fracture geometry or assumptions about how individual fractures control flow and transport. All it requires is the delineation of a few dominant features, which can be embedded into the stochastic continuum model as heterogeneous porous slabs. Though a fault zone has been identified at the Fanay-Augères experimental site, it has been modeled neither by Cacas et al. nor in this paper. In fact, in this paper, a larger selection of experimental results than those considered by Cacas et al. are reproduced merely by modeling the rock as a statistically homogeneous continuum in two dimensions. These results demonstrate that a continuum approach may be well suited for the analysis of flow and transport in fractured rock. This does not constitute a validation of the continuum approach, just as the results of Cacas et al. fall short of validating the discrete fracture approach. Instead, the two sets of results illustrate jointly the well-established principle that an open system, especially one as complex as fractured hydrogeologic environments tend to be, cannot be described uniquely on the basis of sparse data and need not be described in great detail to capture its salient behavior by a model.
• Furman, A., & Neuman, S. P. (2003). Laplace-transform analytic element solution of transient flow in porous media. Advances in Water Resources, 26(12), 1229-1237.
  More info

  Abstract: A Laplace-transform analytic element method (LT-AEM) is described for the solution of transient flow problems in porous media. Following Laplace transformation of the original flow problem, the analytic element method (AEM) is used to solve the resultant time-independent modified Helmholtz equation, and the solution is inverted numerically back into the time domain. The solution is entirely general, retaining the mathematical elegance and computational efficiency of the AEM while being amenable to parallel computation. It is especially well suited for problems in which a solution is required at a limited number of points in space-time, and for problems involving materials with sharply contrasting hydraulic properties. We illustrate the LT-AEM on transient flow through a uniform confined aquifer with a circular inclusion of contrasting hydraulic conductivity and specific storage. Our results compare well with published analytical solutions in the special case of radial flow. © 2003 Elsevier Ltd. All rights reserved.
• Guadagnini, A., Riva, M., & Neuman, S. P. (2003). Three-dimensional steady state flow to a well in a randomly heterogeneous bounded aquifer. Water Resources Research, 39(3), SBH11-SBH116.
  More info

  Abstract: We consider flow in a confined aquifer of uniform thickness due to a well of zero radius that fully penetrates the aquifer and discharges at a constant rate. If the lateral extent of the aquifer is infinite, a steady state flow regime never develops. It is, however, well known that if the aquifer is additionally uniform, a quasi-steady state region extends from the well out to a cylindrical surface whose radius expands as the square root of time. On the expanding surface, head is uniform and time invariant. Inside this surface, head at any time is described by a steady state solution. A rigorous analysis of the analogous situation in a randomly heterogeneous aquifer would require the solution of a three-dimensional transient stochastic flow problem in an aquifer of infinite lateral extent. Here we take a different approach by developing a three-dimensional steady solution for mean flow to a well in a randomly heterogeneous aquifer with a cylindrical prescribed head boundary. In analogy to the uniform case we expect our solution to approximate a quasi-steady state region whose radius is initially small in comparison to the horizontal correlation scale of log conductivity but grows with time to become eventually much larger. We treat log conductivity as a statistically homogeneous random field characterized by a Gaussian spatial covariance function that may have different horizontal and vertical correlation scales. Our solution consists of analytical expressions for the ensemble mean and variance of head in the aquifer to second order in the standard deviation of log conductivity. It is based on recursive approximations of exact nonlocal moment equations that are free of distributional assumptions and so apply to both Gaussian and non-Gaussian log conductivity fields. The analytical solution is supported by numerical Monte Carlo simulations. It clarifies the manner in which relationships between the horizontal and vertical scales of the quasi-steady state region and those of statistical anisotropy impact the statistical moments of drawdown and the equivalent and apparent hydraulic conductivities of the aquifer. Both conductivities are shown to exhibit a scale effect by growing with distance from the well within a radius of one to two horizontal integral scales from it.
• Hernandez, A. F., Neuman, S. P., Guadagnini, A., & Carrera, J. (2003). Conditioning mean study state flow on hydraulic head and conductivity through geostatistical inversion. Stochastic Environmental Research and Risk Assessment, 17(5), 329-338.
  More info

  Abstract: Nonlocal moment equations allow one to render deterministically optimum predictions of flow in randomly heterogeneous media and to assess predictive uncertainty conditional on measured values of medium properties. We present a geostatistical inverse algorithm for steady-state flow that makes it possible to further condition such predictions and assessments on measured values of hydraulic head (and/or flux). Our algorithm is based on recursive finite-element approximations of exact first and second conditional moment equations. Hydraulic conductivity is parameterized via universal kriging based on unknown values at pilot points and (optionally) measured values at other discrete locations. Optimum unbiased inverse estimates of natural log hydraulic conductivity, head and flux are obtained by minimizing a residual criterion using the Levenberg-Marquardt algorithm. We illustrate the method for superimposed mean uniform and convergent flows in a bounded two-dimensional domain. Our examples illustrate how conductivity and head data act separately or jointly to reduce parameter estimation errors and model predictive uncertainty.
• Illman, W. A., & Neuman, S. P. (2003). Steady-state analysis of cross-hole pneumatic injection tests in unsaturated fractured tuff. Journal of Hydrology, 281(1-2), 36-54.
  More info

  Abstract: Numerous single-hole and cross-hole pneumatic injection tests have been conducted in unsaturated fractured tuff at the Apache Leap Research Site (ALRS) near Superior, Arizona. Steady-state analyses of single-hole tests conducted by Guzman et al. [Summary of Air Permeability Data From Single-Hole Injection Tests in Unsaturated Fractured Tuffs at the Apache Leap Research Site: Results of Steady-State Test Interpretation (1996)] have yielded values of air permeability at various locations throughout the tested rock volume on nominal scales from 0.5 to 3.0 m. Analyses of transient data from single-hole and larger-scale cross-hole tests were performed by Illman and Neuman [Ground Water 38 (2000) 899; Water Resour. Res. 37 (2001) 583] using type-curves and by Vesselinov and Neuman [Ground Water 39 (2001) 685] and Vesselinov et al. [Water Resour. Res. 37 (2001a) 3001; Water Resour. Res. 37 (2001b) 3041] using numerical inversion. These have yielded bulk fracture permeabilities and porosities on scales ranging from a few meters to several tens of meters. We complement the latter results by steady-state analysis of a relatively large set of cross-hole tests. Steady-state analysis (a) allows interpreting cross-hole tests that were not interpreted using type-curves due to weak or noisy signals, (b) is relatively easy to perform, and (c) yields results comparable with those obtained from transient analyses. We analyze the results statistically and discuss their implications vis-a-vis the pneumatic properties of unsaturated fractured tuff at the ALRS. Our results strengthen the evidence for a previously surmised permeability scale effect at the site. © 2003 Elsevier B.V. All rights reserved.
• Neuman, S. P. (2003). Maximum likelihood Bayesian averaging of uncertain model predictions. Stochastic Environmental Research and Risk Assessment, 17(5), 291-305.
  More info

  Abstract: Hydrologic analyses typically rely on a single conceptual-mathematical model. Yet hydrologic environments are open and complex, rendering them prone to multiple interpretations and mathematical descriptions. Adopting only one of these may lead to statistical bias and underestimation of uncertainty. A comprehensive strategy for constructing alternative conceptual-mathematical models of subsurface flow and transport, selecting the best among them, and using them jointly to render optimum predictions under uncertainty has recently been developed by Neuman and Wierenga (2003). This paper describes a key formal element of this much broader and less formal strategy that concerns rendering optimum hydrologic predictions by means of several competing deterministic or stochastic models and assessing their joint predictive uncertainty. The paper proposes a Maximum Likelihood Bayesian Model Averaging (MLBMA) method to accomplish this goal. MLBMA incorporates both site characterization and site monitoring data so as to base the outcome on an optimum combination of prior information (scientific knowledge plus data) and model predictions. A preliminary example based on real data is included in the paper.
• Neuman, S. P. (2003). Relationship between juxtaposed, overlapping, and fractal representations of multimodal spatial variability. Water Resources Research, 39(8), SBH41-SBH47.
  More info

  Abstract: It has been shown by Neuman [1990], Di Federico and Neuman [1997], and Di Federico et al. [1999] that the power variogram of a statistically isotropic or anisotropic fractal field can be constructed as a weighted integral from zero to infinity of exponential or Gaussian variograms of overlapping, homogeneous random fields (modes) having mutually uncorrelated increments and variance proportional to a power 2H of the integral (spatial correlation) scale where H is the Hurst coefficient. Low- and high-frequency cutoffs are related to length scales of the sampling window (domain) arid data support (sample volume), respectively. Intermediate cutoffs account for lacunarity due to gaps in the multiscale hierarchy, created by a hiatus of modes associated with discrete ranges of scales. Alternative mathematical representations of multimodal spiatial variability were formulated by various authors in which space is filled by a discrete number of juxtaposed (mutually exclusive) materials or categories, each having its own architecture and attributes. The spatial distribution of categories is characterized by indicator random variables and their attributes by random fields. This paper focuses on expressions developed by Lu and Zhang [2002] in which the indicator variables and their attributes are mutually uncorrelated while each is autocorrelated and cross correlated within and between categories. Upon rewriting their expressions for statistically homogeneous and anisotropic media, it is demonstrated mathematically that in the limit as the categories stretch to occupy each point in space (overlap) in a way that preserves their local architecture, their attributes become mutually uncorrelated. Categories are said to overlap if they are found in fixed proportions within a representative, sampling volume centered about any mathematical point throughout space; the idea is analogous to the well-known dual continuum concept in which two categories, most commonly fractures and porous blocks, are considered to overlap. In reality, the categories do not overlap but are considered to preserve their relative architectural arrangement throughout space. The variogram of an attribute, sampled jointly over overlapping statistically homogeneous and anisotropic categories, is the sum of the variograms of statistically homogeneous and anisotropic components of this attribute sampled over individual categories, weighted by their volumetric proportions. It is further demonstrated that when the overlapping categories represent an infinite or truncated hierarchy of modes whose attributes have exponential or Gaussian variograms, such that the product of category probability density function (pdf) and attribute variance is proportional to a power 2H + 1 of the attribute integral scale, they represent a fractal field characterized by a power or truncated power variogram as in the work of Di Federico and Neuman. Any category pdf that satisfies this condition, and any attribute pdf having a variance that satisfies the same condition, are admissible; neither of them need to be Gaussian (to form fractional Brownian motion) or Levy-like (to approximate fractional Levy motion) for the attribute to form a random fractal. Intermediate cutoffs render the latter model discrete, as is that of Lu and Zhang. Hence the two multimodal models are mutually consistent and complementary.
• Neuman, S. P., & Federico, V. D. (2003). Multifaceted nature of hydrogeologic scaling and its interpretation. Reviews of Geophysics, 41(3), 4-1 - 4-31.
  More info

  Abstract: Geology is ubiquitously heterogeneous, exhibiting both discrete and continuous spatial variations on a multiplicity of scales. It is therefore natural to expect that hydrogeologic and other geophysical variables would do likewise. We present evidence that hydrogeologic variables exhibit isotropic and directional dependencies on scales of measurement (data support), observation (extent of phenomena such as a dispersing plume), sampling window (domain of investigation), spatial correlation (structural coherence), and spatial resolution (descriptive detail). We then show that it is possible to interpret these multiple scale dependencies within a unified theoretical framework. This and similar theoretical frameworks may be applicable to a wider range of geophysical scale issues.
• Tartakovsky, A. M., Neuman, S. P., & Lenhard, R. J. (2003). Immiscible front evolution in randomly heterogeneous porous media. Physics of Fluids, 15(11), 3331-3341.
  More info

  Abstract: The evolution of a sharp interface between two immiscible fluids in a randomly heterogeneous porous medium is investigated analytically using a stochastic moment approach. The displacing fluid is taken to be at constant saturation and to have a much larger viscosity than does the displaced fluid, which is therefore effectively static. Capillary pressure at the interface is related to porosity and permeability via the Leverett J-function. Whereas porosity is spatially uniform, permeability forms a spatially correlated random field. Displacement is governed by stochastic integro-differential equations defined over a three-dimensional domain bounded by a random interface. The equations are expanded and averaged in probability space to yield leading order recursive equations governing the ensemble mean and variance of interface position, rate of propagation and pressure gradient within the displacing fluid. Solutions are obtained for one-dimensional head- and flux-driven displacements in statistically homogeneous media and found to compare well with numerical Monte Carlo simulations. The manner in which medium heterogeneity affects the mean pressure gradient is indicative of how it impacts the stability of the mean interface. Capillary pressure at the interface is found to have a potentially important effect on its mean dynamics and stability. © 2003 American Institute of Physics.
• Wang, W., Neuman, S. P., Yao, T., & Wierenga, P. J. (2003). Simulation of large-scale field infiltration experiments using a hierarchy of models based on public, generic, and site data. Vadose Zone Journal, 2(3), 297-312.
  More info

  Abstract: Comparative simulations of a large-scale field infiltration experiment at the Maricopa Agricultural Center (MAC) near Phoenix, AZ, were conducted using a hierarchy of models based on public, generic, and site data joined with pedotransfer functions and an inverse procedure. Our purpose was to investigate the ability of simple models and relatively inexpensive data to reproduce and predict reliably the time evolution of water content profiles in nine 10-m-deep neutron monitoring boreholes at the site. By relying solely on public sources of information one might conclude that soil at the MAC site is uniform to a depth of 16 m, with a water table at about 22 m. Upon collecting soil samples at the site, we learned the soils are layered and laterally discontinuous, with a perched water table at about 13 m. To identify the least level of complexity required to simulate infiltration at the MAC, we compared models that consider one- and two-dimensional flow in a uniform soil, a soil consisting of uniform layers, and a stratified soil with laterally distinct zones. There is a paucity of hydraulic characterization data for the site. To investigate the feasibility of obtaining hydraulic parameter estimates for the models on the basis of soil type, we ascribed uniform properties to each layer or zone using mean values of three generic databases. To improve these estimates, we ascribed variable soil hydraulic properties to individual soil samples using regression and neural network pedotransfer functions based on soil type and bulk density; we then used them to obtain Bayesian updates of mean hydraulic properties in each layer or zone. We used the various models and mean parameter estimates to simulate water contents during one of several infiltration experiments at the MAC. None of the results compared well with measured water contents, although one set of generic mean parameter values yielded much better results than the other two. Estimating parameters using pedotransfer functions and Bayesian updating does not lead to improved simulations. Only when the parameters are estimated by means of an inverse procedure does one notice a significant improvement in model fit. We comparedand ranked the various models and parameter estimates using likelihood-based model discrimination criteria and confirmed our choice of best model by successfully simulating flow during an earlier infiltration experiment. © Soil Science Society of America.
• Hernandez, A. F., Neuman, S. P., Guadagnini, A., & Carrera, J. (2002). Conditioning steady state mean stochastic flow equations on head and hydraulic conductivity measurement. Acta Universitatis Carolinae, Geologica, 46(2-3), 158-162.
  More info

  Abstract: Nonlocal moment equations allow one to render deterministically optimum predictions of flow in randomly heterogeneous media and to assess predictive uncertainty conditional on measured values of medium properties. We present a geostatistical inverse algorithm for steady state flow that makes it possible to further condition such predictions and assessments on measured values of hydraulic head (and/or flux). Our algorithm is based on recursive finite-element approximations of exact first and second conditional moment equations. Hydraulic conductivity is parameterized via universal kriging based on unknown values at pilot points and (optionally) measured values at other discrete locations. We illustrate the method for superimposed mean uniform and convergent flows in a bounded two-dimensional domain. Our examples illustrate how conductivity and head data act separately or jointly to reduce parameter estimation errors and model predictive uncertainty.
• Hernandez, A. F., Neuman, S. P., Guadagnini, A., & Carrera, J. (2002). Conditioning steady-state mean stochastic flow equations on head and hydraulic conductivity measurements. IAHS-AISH Publication, 122-128.
  More info

  Abstract: Non-local moment equations allow one to obtain deterministically optimum predictions of flow in randomly heterogeneous media and to assess predictive uncertainty conditional on measured values of medium properties. We present a geostatistical inverse algorithm for steady-state flow that makes it possible to further condition such predictions and assessments on measured values of hydraulic head (and/or flux). Our algorithm is based on recursive finite-element approximations of exact first and second conditional moment equations. Hydraulic conductivity is parameterized via universal kriging based on unknown values at pilot points and (optionally) measured values at other discrete locations. We illustrate the method for superimposed mean uniform and convergent flows in a bounded two-dimensional domain. Our examples illustrate how conductivity and head data act separately or jointly to reduce parameter estimation errors and model predictive uncertainty.
• Hyun, Y., Neuman, S. P., Vesselinov, V. V., Illman, W. A., Tartakovsky, D. M., & Federico, V. D. (2002). Theoretical interpretation of a pronounced permeability scale effect in unsaturated fractured tuff. Water Resources Research, 38(6), 281-288.
  More info

  Abstract: Numerous single-hole and cross-hole pneumatic injection tests have been conducted in unsaturated fractured tuff at the Apache Leap Research Site (ALRS) near Superior, Arizona. Single-hole tests have yielded values of air permeability at various locations throughout the tested rock volume on a nominal scale of ∼1 m. Cross-hole tests have yielded equivalent air permeabilities (and air-filled porosities) for a rock volume characterized by a length scale of several tens of meters. Cross-hole tests have also provided high-resolution tomographic estimates of how air permeability (and air-filled porosity), defined over grid blocks having a length scale of 1 m, vary throughout a similar rock volume. The results have revealed a highly pronounced scale effect in permeability (and porosity) at the ALRS. We examine the extent to which the permeability scale effect is amenable to interpretation by a recent stochastic scaling theory, which treats the rock as a truncated random fractal.
• Neuman, S. P. (2002). Accounting for conceptual model uncertainty via maximum likelihood Bayesian model averaging. IAHS-AISH Publication, 303-313.
  More info

  Abstract: Analyses of groundwater flow and transport typically rely on a single conceptual model of site hydrogeology. Yet hydrogeological environments are open and complex, rendering them prone to multiple interpretations. Adopting only one of these may lead to statistical bias and underestimation of uncertainty. A comprehensive strategy for constructing alternative conceptual-mathematical models, selecting the best among them, and using them jointly to render optimum predictions under uncertainty is being developed by the author. This paper proposes a maximum likelihood Bayesian model averaging approach, MLBMA, to rendering optimum predictions by means of several competing models and assessing their joint predictive uncertainty.
• Neuman, S. P. (2002). Accounting for conceptual model uncertainty via maximum likelihood bayesian model averaging. Acta Universitatis Carolinae, Geologica, 46(2-3), 529-534.
  More info

  Abstract: Analyses of groundwater flow and transport typically rely on a single conceptual model of site hydrogeology. Yet hydrogeologic environments are open and complex, rendering them prone to multiple interpretations. Adopting only one of these may lead to statistical bias and underestimation of uncertainty. A comprehensive strategy for constructing alternative conceptual-mathematical models, selecting the best among them, and using them jointly to render optimum predictions under uncertainty is being developed by the author. This paper proposes a Maximum Likelihood Bayesian Model Averaging approach, MLBMA, to rendering optimum predictions by means of several competing models and assessing their joint predictive uncertainty.
• Zhiming, L. u., Neuman, S. P., Guadagnini, A., & Tartakovsky, D. M. (2002). Conditional moment analysis of steady state unsaturated flow in bounded, randomly heterogeneous soils. Water Resources Research, 38(4), 91-916.
  More info

  Abstract: We consider steady state unsaturated flow in bounded, randomly heterogeneous soils under the influence of random boundary and source terms. Our aim is to predict pressure heads and fluxes without resorting to Monte Carlo simulation, upscaling, or linearization of the constitutive relationship between unsaturated hydraulic conductivity and pressure head. We represent this relationship through Gardner's exponential model while treating its exponent α as a random constant and saturated hydraulic conductivity Ks as a spatially correlated random field. We linearize the steady state unsaturated flow equations by means of the Kirchhoff transformation and integrate them in probability space to obtain exact integro-differential equations for the conditional mean and variance-covariance of transformed pressure head and flux. After approximating these equations recursively to second order in the standard deviation σY of Y = ln Ks, we solve them by finite elements for superimposed mean uniform and divergent flows in the vertical plane, with and without conditioning on measured Y values. Comparison with Monte Carlo solutions demonstrates that whereas our nonlocal solution is nominally restricted to mildly nonuniform media with σY2 ≪ 1, it yields remarkably accurate results for strongly nonuniform media with σY2 at least as large as 2. This accords well with a previous theoretical analysis, which shows that the solution may remain asymptotic for values of σY2 as large as 2.
• Amir, O., & Neuman, S. P. (2001). Gaussian closure of one-dimensional unsaturated flow in randomly heterogeneous soils. Transport in Porous Media, 44(2), 355-383.
  More info

  Abstract: We propose a new method for the solution of stochastic unsaturated flow problems in randomly heterogeneous soils which avoids linearizing the governing flow equations or the soil constitutive relations, and places no theoretical limit on the variance of constitutive parameters. The proposed method applies to a broad class of soils with flow properties that scale according to a linearly separable model provided the dimensionless pressure head ψ has a near-Gaussian distribution. Upon treating ψ as a multivariate Gaussian function, we obtain a closed system of coupled nonlinear differential equations for the first and second moments of pressure head. We apply this Gaussian closure to steady-state unsaturated flow through a randomly stratified soil with hydraulic conductivity that varies exponentially with j where ψ = (1/α)ψ is dimensional pressure head and α is a random field with given statistical properties. In one-dimensional media, we obtain good agreement between Gaussian closure and Monte Carlo results for the mean and variance of ψ over a wide range of parameters provided that the spatial variability of α is small. We then provide an outline of how the technique can be extended to two- and three-dimensional flow domains. Our solution provides considerable insight into the analytical behavior of the stochastic flow problem.
• Guadagnini, A., & Neuman, S. P. (2001). Recursive Conditional Moment Equations for Advective Transport in Randomly Heterogeneous Velocity Fields. Transport in Porous Media, 42(1), 37-67.
  More info

  Abstract: Flow and transport parameters such as hydraulic conductivity, seepage velocity, and dispersivity have been traditionally viewed as well-defined local quantities that can be assigned unique values at each point in space-time. Yet in practice these parameters can be deduced from measurements only at selected locations where their values depend on the scale (support volume) and mode (instruments and procedure) of measurement. Quite often, the support of the measurements is uncertain and the data are corrupted by experimental and interpretive errors. Estimating the parameters at points where measurements arc not available entails an additional random error. These errors and uncertainties render the parameters random and the corresponding How and transport equations stochastic. The stochastic flow and transport equations can be solved numerically by conditional Monte Carlo simulation. However, this procedure is computationally demanding and lacks well-established convergence criteria. An alternative to such simulation is provided by conditional moment equations, which yield corresponding predictions of flow and transport deterministically. These equations are typically integro-differential and include nonlocal parameters that depend on more than one point in space-time. The traditional concept of a REV (representative elementary volume) is neither necessary nor relevant for their validity or application. The parameters are nonunique in that they depend not only on local medium properties but also on the information one has about these properties (scale, location, quantity, and quality of data). Darcy's law and Fick's analogy are generally not obeyed by the flow and transport predictors except in special cases or as localized approximations. Such approximations yield familiar-looking differential equations which, however, acquire a non-traditional meaning in that their parameters (hydraulic conductivity, seepage velocity, dispersivity) and state variables (hydraulic head, concentration) are information-dependent and therefore, inherently nonunique. Nonlocal equations contain information about predictive uncertainty, localized equations do not. We have shown previously (Guadagnini and Neuman, 1997, 1998, 1999a, b) how to solve conditional moment equations of steady-state flow numerically on the basis of recursive approximations similar to those developed for transient flow by Tartakovsky and Neuman (1998, 1999). Our solution yields conditional moments of velocity, which are required for the numerical computation of conditional moments associated with transport. In this paper, we lay the theoretical groundwork for such computations by developing exact integro-differential expressions for second conditional moments, and recursive approximations for all conditional moments, of advective transport in a manner that complements earlier work along these lines by Neuman (1993).
• Illman, W. A., & Neuman, S. P. (2001). Type curve interpretation of a cross-hole pneumatic injection test in unsaturated fractured tuff. Water Resources Research, 37(3), 583-603.
  More info

  Abstract: Cross-hole pneumatic injection tests have been conducted in 16 vertical and inclined boreholes in unsaturated fractured tuffs at the Apache Leap Research Site (ALRS) near Superior, Arizona. Their purpose was to characterize the bulk pneumatic properties and connectivity of fractures at the site on scales ranging from meters to several tens of meters. We describe the design, conduct, and type curve interpretation of one of these tests. Our cross-hole type curves are modified after Hsieh and Neuman [1985] to consider single-phase airflow and extended to consider the effects of storage and skin in monitoring intervals. Cross-hole type curves of pressure derivatives and recovery are included for improved pneumatic characterization of the site. We analyze recorded pressures in each monitoring interval separately from those in other intervals while treating the fractured rock as a uniform, isotropic porous continuum. Each record yields an equivalent directional air permeability and air-filled porosity for fractures that connect the corresponding monitoring and injection intervals. Both parameters are found to vary considerably from one monitoring interval to another, reflecting the nonuniform nature of pneumatic rock properties at the ALRS. The geometric mean of these equivalent permeabilities is found to be larger by a factor of 50 than that obtained from single-hole pneumatic injection tests on a nominal scale of 1-3 m (single-hole tests yield only limited information about porosities, which therefore cannot be meaningfully compared with cross-hole results). Our results find support in numerical inverse modeling of cross-hole tests at the site by Vesselinov [2000; see also Illman et al., 1998; Vesselinov et al., 2000]. Vesselinov [2000] has further demonstrated [see also Chen et al., 2000] that a similar scale effect is exhibited by fracture porosity at the ALRS and that both scale effects disappear when cross-hole tests at the site are interpreted by means of a numerical inverse model, which resolves heterogeneities down to a scale of 1 m. When considered jointly with these and other studies of the site, our analysis implies that the pneumatic pressure behavior of unsaturated fractured tuffs at the ALRS can be described quite accurately by means of linearized single-phase airflow equations; this behavior can be interpreted by treating the rock as a continuum on scales ranging from meters to tens of meters: the continuum is representative primarily of interconnected fractures; as these fractures are filled primarily with air, their pneumatic permeabilities and porosities are close to the bulk intrinsic properties of fractures at the site: these intrinsic properties vary randomly with location and direction across the ALRS, and they depend strongly on the scale at which they are determined.
• Riva, M., Guadagnini, A., Neuman, S. P., & Franzetti, S. (2001). Radial flow in a bounded randomly heterogeneous aquifer. Transport in Porous Media, 45(1), 139-193.
  More info

  Abstract: Flow to wells in nonuniform geologic formations is of central interest to hydrogeologists and petroleum engineers. There are, however, very few mathematical analyses of such flow. We present analytical expressions for leading statistical moments of vertically averaged hydraulic head and flux under steady-state flow to a well that pumps water from a bounded, randomly heterogeneous aquifer. Like in the widely used Thiem equation, we prescribe a constant pumping rate deterministically at the well and a constant head at a circular outer boundary of radius L. We model the natural logarithm Y = In T of aquifer transmissivity T as a statistically homogeneous random field with a Gaussian spatial correlation function. Our solution is based on exact nonlocal moment equations for multidimensional steady state flow in bounded, randomly heterogeneous porous media. Perturbation of these nonlocal equations leads to a system of local recursive moment equations that we solve analytically to second order in the standard deviation of Y. In contrast to most stochastic analyses of flow, which require that log transmissivity be multivariate Gaussian, our solution is free of any distributional assumptions. It yields expected values of head and flux, and the variance-covariance of these quantities, as functions of distance from the well. It also yields an apparent transmissivity, Ta, defined as the negative ratio between expected flux and head gradient at any radial distance. The solution is supported by numerical Monte Carlo simulations, which demonstrate that it is applicable to strongly heterogeneous aquifers, characterized by large values of log transmissivity variance. The two-dimensional nature of our solution renders it useful for relatively thin aquifers in which vertical heterogeneity tends to be of minor concern relative to that in the horizontal plane. It also applies to thicker aquifers when information about their vertical heterogeneity is lacking, as is commonly the case when measurements of head and flow rate are done in wells that penetrate much of the aquifer thickness. Potential uses include the analysis of pumping tests and tracer test conducted in such wells, the statistical delineation of their respective capture zones, and the analysis of contaminant transport toward fully penetrating wells.
• Vesselinov, V. V., & Neuman, S. P. (2001). Numerical inverse interpretation of single-hole pneumatic tests in unsaturated fractured tuff. Ground Water, 39(5), 685-695.
  More info

  PMID: 11554246;Abstract: A numerical inverse method was used to interpret simultaneously multirate injection and recovery data from single-hole pneumatic tests in unsaturated fractured tuff at the Apache Leap Research Site near Superior, Arizona. Our model represents faithfully the three-dimensional geometry of boreholes at the site, and accounts directly for their storage and conductance properties by treating them as high-permeability and high-porosity cylinders of finite length and radius. It solves the airflow equations in their original nonlinear form and yields information about air permeability, air-filled porosity and dimensionless borehole storage coefficient. Some of this is difficult to accomplish with analytical type-curves. Air permeability values obtained by our inverse method agree well with those obtained by steady-state and type-curve analyses.
• Vesselinov, V. V., Neuman, S. P., & Illman, W. A. (2001). Three-dimensional numerical inversion of pneumatic cross-hole tests in unsaturated fractured tuff 1. Methodology and borehole effects. Water Resources Research, 37(12), 3001-3017.
  More info

  Abstract: We describe 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. The model combines a finite volume flow simulator, FEHM, an automatic mesh generator, X3D, a parallelized version of an automatic parameter estimator, PEST, and a geostatistical package, GSTAT. The tests are simulated by considering single-phase airflow through a porous continuum, which represents primarily interconnected fractures at the site. The simulator solves the airflow equations in their original nonlinear form and accounts directly for the ability of all packed-off borehole intervals to store and conduct air through the system. Computations are performed in parallel on a supercomputer using 32 processors. We analyze pneumatic cross-hole test data, previously conducted by our group at ALRS, 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 the rock volume being tested, having length scales of the order of meters to tens of meters. 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, a concept originally proposed by one of us in 1987. The first paper of this two-part series describes the field data, the model, and the effect of boreholes on pressure propagation through the rock. The second paper implements our approach on selected cross-hole test data from ALRS.
• Vesselinov, V. V., Neuman, S. P., & Illman, W. A. (2001). Three-dimensional numerical inversion of pneumatic cross-hole tests in unsaturated fractured tuff 2. Equivalent parameters, high-resolution stochastic imaging and scale effects. Water Resources Research, 37(12), 3019-3041.
  More info

  Abstract: 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.
• Amir, O., & Neuman, S. P. (2000). Gaussian analysis of one-dimensional unsaturated flow in randomly heterogeneous soils. Special Paper of the Geological Society of America, 348, 105-117.
  More info

  Abstract: We propose a new method for the solution of stochastic unsaturated flow problems in randomly heterogeneous soils that applies to both steady-state and transient regimes, avoids linearizing the soil constitutive relations, and places no theoretical limit on the variance of constitutive parameters. The proposed method applies to a broad class of soils with flow properties that scale according to the linearly separable model of L. Vogel and coworkers, provided that the dimensionless pressure head ψ has a near-Gaussian distribution. Upon treating ψ as a multivariate Gaussian function, we obtain a closed system of coupled nonlinear differential equations for the first and second moments of pressure head. We apply this Gaussian closure to one-dimensional steady-state unsaturated flow through a randomly stratified soil with hydraulic conductivity that varies exponentially with αψ where ψ=1/α ψ is dimensional pressure head and α is a random field with given statistical properties. Upon disregarding correlation between αand ψ, we obtain good agreement between Gaussian closure and Monte Carlo results for the mean of ψover a wide range of parameters; the agreement between variances is good far from the soil bottom where pressure head is prescribed, but not as good near the bottom boundary. We expect that the latter problem will be remedied once we account for correlation between α and ψ. Our solution provides considerable insight into the analytical behavior of the stochastic flow problem.
• Chen, G., Neuman, S. P., & Hills, R. G. (2000). Wetting front instability and fingering in random permeability fields. Special Paper of the Geological Society of America, 348, 53-59.
  More info

  Abstract: Numerical simulations are used to investigate wetting front instability and finger evolution in the vertical plane during surface-water infiltration into a randomly heterogeneous soil. As is common in the literature, the natural log permeability of the soils is taken to be a random, multivariate Gaussian function of space. The mean (expectation) of this function is constant and its fluctuations about the mean are statistically homogeneous and isotropic with constant variance and autocorrelation scale. The wetting front is allowed to diffuse in accord with Richards' equation of flow in a partially saturated soil. The infiltrating water is labeled with a passive tracer, the 50% concentration contour of which is taken to indicate the position of the wetting front. The variogram of the contour is introduced to characterize the length and width of the fingers. Results show that whereas mean finger length grows at a nonlinear rate that diminishes monotonically with time, mean finger width remains constant in time. Both mean finger growth rate and mean finger width increase with the variances and correlation scale of permeability.
• Illman, W. A., & Neuman, S. P. (2000). Type-curve interpretation of multirate single-hole pneumatic injection tests in unsaturated fractured rock. Ground Water, 38(6), 899-911.
  More info

  Abstract: More than 270 multirate single-hole pneumatic injection tests were conducted by Guzman et al. (1996) within six shallow vertical and inclined boreholes in unsaturated fractured tuff at the Apache Leap Research Site (ALRS) near Superior, Arizona. The authors used steady-state formulae to obtain air permeability values for the rock based solely on late pressure data from each test. We developed pressure and pressure-derivative type-curves for the interpretation of transient data from these tests, which account for storage in the test interval due to air compressibility and skin effect due to damage to the surrounding rock. We applied our type-curves to pressure data from more than 40 of the tests to obtain information about air permeability, skin factor, phenomenology, and dimensionality of the flow regime on a nominal scale of 1 m in the immediate vicinity of each test interval. Our air permeabilities agreed well with those determined previously by means of steady- state formulae but correlated poorly with fracture density data. Nonlinear effects due to air compressibility were pronounced, but skin effect was generally small or nonexistent. Nonlinear two-phase flow and inertia were evident but did not have a significant impact on type-curve interpretation of the test data. Most of these data fit a continuum model of spherical airflow around the test interval, but a few fit a continuum model of radial flow or suggest the dominance of a single fracture. We take this to imply that fractures at the site are pneumatically well connected and can, for the most part, be treated as a porous continuum when analyzing airflow on scales of 1 m or more. This is supported by the finding that 1 m scale air permeability data from single-hole tests at the ALRS are amenable to geostatistical analysis (Chen et al. 2000), which views them as a sample from a random field defined over a continuum. That the same holds true on larger scale is indicated by our ability to interpret cross-hole pneumatic injection tests at the site by means of type-curves and inverse methods (Illman et al. 1998), which treat the entire rock mass as a porous continuum.
• Lu, Z., Neuman, S. P., Guadagnini, A., & Tartakovsky, D. M. (2000). Direct solution of unsaturated flow in randomly heterogeneous soils. Computational methods in water resources - Volume 2 - Computational methods,surface water systems and hydrology, 785-792.
  More info

  Abstract: We consider steady state unsaturated flow in bounded randomly heterogeneous soils under the influence of random forcing terms. Our aim is to predict pressure heads and fluxes without resorting to Monte Carlo simulation, upscaling or linearization of the constitutive relationship between unsaturated hydraulic conductivity and pressure head. We represent this relationship through Gardner's exponential model, treating its exponent α as a random constant and saturated hydraulic conductivity, K(s), as a spatially correlated random field. This allows us to linearize the steady state unsaturated flow equations by means of the Kirchhoff transformation, integrate them in probability space, and obtain exact integro-differential equations for the conditional mean and variance-covariance of transformed pressure head and flux. We solve the latter for flow in the vertical plane, with a point source, by finite elements to second-order of approximation. Our solution compares favorably with conditional Monte Carlo simulations, even for soils that are strongly heterogeneous.
• Neuman, S. P., & Guadagnini, A. (2000). A new look at traditional deterministic flow models and their calibration in the context of randomly heterogeneous media. IAHS-AISH Publication, 213-221.
  More info

  Abstract: We examine the tradition of modelling subsurface flow deterministically from a stochastic viewpoint. In particular, we show that traditional deterministic flow equations do apply to randomly heterogeneous media, albeit in an approximate manner, provided they are interpreted in a non-traditional manner. Our paper explains why parameter estimates obtained by traditional inverse methods tend to vary as one modifies the database. It also makes clear that the traditional Monte Carlo method of assessing uncertainty in the output of a calibrated deterministic model generally overestimates the predictive capabilities of the model. The only valid way to assess predictive uncertainty is by means of a stochastic model.
• Riva, M., Guadagnini, A., Neuman, S. P., & Tartakovsky, D. M. (2000). Three-dimensional steady state flow to a well in a randomly heterogeneous aquifer. IAHS-AISH Publication, 131-136.
  More info

  Abstract: We consider three-dimensional steady state flow towards a well that fully penetrates a randomly heterogeneous aquifer confined between horizontal no-flow boundaries, and bounded laterally by a cylindrical, deterministically prescribed constant head boundary. The well is represented by a line sink that produces water at a deterministically prescribed constant rate Q for unit aquifer thickness. The log hydraulic conductivity, Y = InK, of the aquifer is multivariate Gaussian, statistically homogeneous with a Gaussian spatial autocorrelation function. We develop an analytical solution for the variance of hydraulic head as a function of dimensionless vertical and horizontal locations within the aquifer, variance σ(Y)2 of Y and dimensionless ratios between the principal spatial correlation scales. Our analysis is based on the non-local theory first proposed for steady state and transient flows in bounded, randomly heterogeneous media by Neuman and Orr (1993), Neuman et al. (1996), Guadagnini and Neuman (1999a,b) and Tartakovsky and Neuman (1998, 1999). In particular, we develop and solve analytically recursive closure approximations of the governing non-local moment equations to second order in σ(Y) by means of an appropriate Green's function. We evaluate our analytical solutions by means of Gaussian quadratures for the special case of the isotropic Y field.
• Vesselinov, V. V., & Neuman, S. P. (2000). Numerical inverse interpretation of multistep transient single-hole pneumatic tests in unsaturated fractured tuffs at the Apache Leap Research Site. Special Paper of the Geological Society of America, 348, 175-190.
  More info

  Abstract: More than 270 single-hole multiple-step air-injection (pneumatic) tests have been conducted by A.G. Guzman and coworkers in six shallow vertical and slanted boreholes in unsaturated fractured tuffs at the Apache Leap Research Site near Superior, Arizona. Guzman and coworkers used steady-state formulae for single-phase air flow in a uniform, isotropic porous continuum to interpret late data from each step of an injection test. W.A. Illman and coworkers used transient type curves for single-phase air flow in a similar continuum to analyze all data from the first step of several injection tests; however, this analysis did not allow reliable identification of air-filled porosity and the dimensionless borehole storage coefficient. This chapter describes numerical inversion of multistep and recovery data from some of these same tests, based on similar assumptions, by means of a three-dimensional finite-volume code (FEHM) coupled with a parameter-estimation code (PEST). Our numerical inverse model accounts directly for the geometry, flow properties, and storage capabilities of open borehole intervals by treating them as high-permeability and high-porosity cylinders of finite length and radius. It also allows interpreting multiple injection-step and recovery data simultaneously, and yields information about air permeability, air-filled porosity, and the dimensionless borehole storage coefficient. Some of this is difficult to accomplish with the analytical type-curve method. Air permeability values obtained by our inverse method agree well with those obtained by steady-state and type-curve analyses.
• Vesselinov, V. V., Neuman, S. P., & Illman, W. A. (2000). Three-dimensional inverse modelling of pneumatic tests in unsaturated fractured rocks. IAHS-AISH Publication, 243-249.
  More info

  Abstract: A three-dimensional inverse model was used to analyse air pressure data from single- and cross-hole air injection tests in unsaturated fractured tuffs at the Apache Leap Research Site in Arizona. The model simulates single-phase airflow in a uniform, isotropic continuum. It treats open borehole intervals as high-permeability and high-porosity cylinders of finite length and radius. Single-hole test data yield local-scale permeabilities, porosities and borehole storage coefficients. Cross-hole test data yield larger-scale permeabilities, porosities and their spatial distributions.
• Vesselinov, V. V., Neuman, S. P., Illman, W. A., & Zyvoloski, G. A. (2000). Three-dimensional inversion of pneumatic tests in fractured rocks. Computational methods in water resources - Volume 1 - Computational methods for subsurface flow and transport, 307-314.
  More info

  Abstract: A three-dimensional geostatistically-based numerical inverse model was developed for the interpretation of cross-hole pneumatic tests in unsaturated fractured tuffs at the Apache Leap Research Site (ALRS) in Arizona, USA. The model combines a finite-volume flow simulator, FEHM, an automatic mesh generator, X3D, a parallelized version of an automatic parameter estimator, PEST, and a geostatistical code, GSTAT. The tests are simulated by considering singel-phase airflow through an equivalent stochastic porous continuum. The simulator accounts directly for the ability of all packed-off borehole intervals to store and conduct air through the system. Computations are performed in parallel on a supercomputer using 32 processors. We analyze data from several pneumatic cross-hole tests simultaneously to assess the spatial distribution of air permeability and air-filled porosity throughout the test volume. The analysis amounts to three-dimensional pneumatic 'tomography' or stochastic imaging of the rock, a concept originally proposed in connection with hydraulic cross-hole tests in fractured crystalline rocks by Neuman (1987).
• Federico, V. D., Neuman, S. P., & Tartakovsky, D. M. (1999). Anisotropy, lacunarity, and upscaled conductivity and its autocovariance in multiscale random fields with truncated power variograms. Water Resources Research, 35(10), 2891-2908.
  More info

  Abstract: It has been shown by Di Federico and Neuman [1997, 1998a, b] that observed multiscale behaviors of subsurface fluid flow and transport variables can be explained within the context of a unified stochastic framework, which views hydraulic conductivity as a random fractal characterized by a power variogram. Such a random field is statistically nonhomogeneous but possesses homogeneous spatial increments. Di Federico and Neuman [1997] have demonstrated that the power variogram and associated spectra of a statistically isotropic fractal field can be constructed as a weighted integral from zero to infinity (an infinite hierarchy) of exponential or Gaussian variograms and spectra of mutually uncorrelated fields (modes) that are homogeneous and isotropic. We show in this paper that the same holds true when the field and its constituent modes are statistically anisotropic, provided the ratios between principal integral (spatial correlation) scales are the same for all modes. We then analyze the effect of filtering out (truncating) modes of low, high, and intermediate spatial frequency from this infinite hierarchy in the real and spectral domains. A low-frequency cutoff renders the truncated hierarchy homogeneous. The integral scales of the lowest- and highest-frequency cutoff modes are related to length scales of the sampling window (domain) and data support (sample volume), respectively. Taking the former to be proportional to the latter renders expressions for the integral scale and variance of the truncated field dependent on window and support scale (in a manner previously shown to be consistent with observations in the isotropic case). It also allows (in principle) bridging across scales at a specific locale, as well as among locales, by adopting either site-specific or generalized variogram parameters. The introduction of intermediate cutoffs allows us to account, in a straightforward manner, for lacunarity due to gaps in the multiscale hierarchy created by the absence of modes associated with discrete ranges of scales (for example, where textural and structural features are associated with distinct ranges of scale, such as fractures having discrete ranges of trace length and density, which dissect the rock into matrix blocks having corresponding ranges of sizes). We explore mathematically and graphically the effects that anisotropy and lacunarity have on the integral scale, variance, covariance, and spectra of a truncated fractal field. We then develop an expression for the equivalent hydraulic conductivity of a box-shaped porous block, embedded within a multiscale log hydraulic conductivity field, under mean-uniform flow. The block is larger than the support scale of the field but is smaller than a surrounding sampling window. Consequently, its equivalent hydraulic conductivity is a random variable whose variance and spatial autocorrelation function, conditioned on a known mean value of support-scale conductivity across the window, are given explicitly by our multiscale theory.
• Guadagnini, A., & Neuman, S. P. (1999). Nonlocal and localized analyses of conditional mean steady state flow in bounded, randomly nonuniform domains 1. Theory and computational approach. Water Resources Research, 35(10), 2999-3018.
  More info

  Abstract: We consider the effect of measuring randomly varying hydraulic conductivities K(x) on one's 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, (c) and (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 (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 in K(x); and show how to solve these equations to first order in σ2(y) 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 (c) ≃ - K(c)(x) Δ (h(x)>(c) in which K(c) (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.
• Guadagnini, A., & Neuman, S. P. (1999). Nonlocal and localized analyses of conditional mean steady state flow in bounded, randomly nonuniform domains 2. Computational examples. Water Resources Research, 35(10), 3019-3039.
  More info

  Abstract: In paper 1 [Guadagnini and Neuman, this issue] of this two-part series we presented exact nonlocal equations for first and second conditional ensemble moments of hydraulic head and flux under steady state groundwater flow in bounded, randomly nonuniform domains in the presence of uncertain source and boundary terms. We derived recursive approximations for these equations based on expansion in powers of a small parameter σ(y), which represents the standard estimation error of log hydraulic conductivity; developed a two-dimensional finite element scheme for their solution, under the action of deterministic forcing terms, to first order in σ(y)/2; showed how to localize the conditional mean flow equations so as to yield familiar-looking, Darcian differential equations in which, however, all quantities are conditional on data and therefore nonunique; and proposed that traditional deterministic flow models be reinterpreted in light of these localized equations. We stressed that all quantities which enter into our nonlocal and localized conditional moment equations are defined on a consistent support scale thereby obviating the need for upscaling (from support to grid scale) or the a priori definition of a computational grid; that such a grid can therefore be defined a posteriori based on the degree of smoothness one expects the conditional moments to exhibit; and that conditioning points are often spaced widely enough to render these moments relatively smooth and hence resolvable on a coarser grid than is necessary for Monte Carlo simulations. This notwithstanding, we use in this paper a high-resolution grid to explore superimposed mean-uniform and convergent flows in unconditional and conditional log hydraulic conductivity domains and to compare the relative accuracies of corresponding nonlocal, localized, and Monte Carlo finite element solutions. Whereas our nonlocal solution is nominally restricted to mildly nonuniform media with σ(y)/2
• Tartakovsky, D. M., & Neuman, S. P. (1999). Erratum: Transient flow in bounded randomly heterogeneous domains, 1, Exact conditional moment equations and recursive approximations (Water Resources Research (1998) 34:1 (1-12)). Water Resources Research, 35(6), 1947-.
• Tartakovsky, D. M., & Neuman, S. P. (1999). Extension of 'Transient flow in bounded randomly heterogeneous domains, 1, Exact conditional moment equations and recursive approximations'. Water Resources Research, 35(6), 1921-1925.
  More info

  Abstract: In a recent paper we developed exact nonlocal conditional moment equations for transient flow in bounded domains driven by random forcing terms (sources, initial head, and boundary conditions). Whereas our conditional mean equations took into account the randomness of forcing terms, our conditional second moment equations did not. We do so in this brief addendum.
• Tartakovsky, D. M., Neuman, S. P., & Zhiming, L. u. (1999). Conditional stochastic averaging of steady state unsaturated flow by means of Kirchhoff transformation. Water Resources Research, 35(3), 731-745.
  More info

  Abstract: We consider the effect of measuring randomly varying soil hydraulic properties on one's ability to predict steady state unsaturated flow subject to random sources and/or initial and boundary conditions. Our aim is to allow optimum unbiased prediction of system states (pressure head, water content and fluxes deterministically, without upscaling and without linearizing the constitutive characteristics of the soil. It has been shown by Neuman et al. [1999] that such prediction is possible by means of first ensemble moments of system states and fluxes, conditioned on measured values of soil properties, when the latter scale in a linearly separable fashion as proposed by Vogel et al. [1991]; the uncertainty associated with such predictions can be quantified by means of the corresponding conditional second moments. The derivation of moment equations for soils whose properties do not scale in the above manner requires linearizing the corresponding constitutive relations, which may lead to major inaccuracies when these relations are highly nonlinear, as is often the case in nature. When the scaling parameter of pressure head is a random variable independent of location, the steady state unsaturated flow equation can be linearized by means of the Kirchhoff transformation for gravity-free flow. Linearization is also possible in the presence of gravity when hydraulic conductivity varies exponentially with pressure head. For the latter case we develop exact conditional first- and second-moment equations which are nonlocal and therefore non-Darcian. We solve these equations analytically by perturbation for unconditional vertical infiltration and compare our solution with the results of numerical Monte Carlo simulations. Our analytical solution demonstrates in a rigorous manner that the concept of effective hydraulic conductivity does not apply to ensemble-averaged unsaturated flow except when gravity is the sole driving force.
• Federico, V. D., & Neuman, S. P. (1998). Flow in multiscale log conductivity fields with truncated power variograms. Water Resources Research, 34(5), 975-987.
  More info

  Abstract: In a previous paper we offered an interpretation for the observation that the log hydraulic conductivity of geologic media often appears to be statistically homogeneous but with variance and integral scale which grow with the size of the region (window) being sampled. We did so by demonstrating that such behavior is typical of any random field with a truncated power (semi)variogram and that this field can be viewed as a truncated hierarchy of mutually uncorrelated homogeneous fields with either exponential or Gaussian spatial autocovariance structures. The low- and high-frequency cutoff scales λ(l) and λ(u) are related to the length scales of the sampling window and data support, respectively. We showed how this allows the use of truncated power variograms to bridge information about a multiscale random field across windows of different sizes, either at a given locale or between different locales. In this paper we investigate mean uniform steady state groundwater flow in unbounded domains where the log hydraulic conductivity forms a truncated multiscale hierarchy of Gaussian fields, each associated with an exponential autocovariance. We start by deriving an expression for effective hydraulic conductivity, as a function of the Hurst coefficient H and the cutoff scales in one-, two-, and three-dimensional domains which is qualitatively consistent with experimental data. We then develop leading-order analytical expressions for two- and three-dimensional autocovariance and cross-covariance functions of hydraulic head, velocity, and log hydraulic conductivity versus H, λ(l) and λ(u); examine their behavior; and compare them with those corresponding to an exponential log hydraulic conductivity autocovariance. Our results suggest that it should be possible to bridge information about hydraulic heads and groundwater velocities across windows of disparate scales. In particular, when λ(l) >> λ(u), the variance of head is infinite in two dimensions and grows in proportion to λ(l)/(2+2H) in three dimensions, while the variance and longitudinal integral scale of velocity grow in proportion to λ(l)/(2H) and λ(l), respectively, in both cases.
• Federico, V. D., & Neuman, S. P. (1998). Transport in multiscale log conductivity fields with truncated power variograms. Water Resources Research, 34(5), 963-973.
  More info

  Abstract: Di Federico and Neuman [this issue] investigated mean uniform steady state groundwater flow in an unbounded domain with log hydraulic conductivity that forms a truncated multiscale hierarchy of statistically homogeneous and isotropic Gaussian fields, each associated with an exponential autocovariance. Here we present leading-order expressions for the displacement covariance, and dispersion coefficient, of an ensemble of solute particles advected through such a flow field in two or three dimensions. Both quantities are functions of the mean travel distance s, the Hurst coefficient H, and the low- and high-frequency cutoff integral scales λ(l) and λ(u). The latter two are related to the length scales of the sampling window (region under investigation) and sample volume (data support), respectively. If one considers transport to be affected by a finite domain much larger than the mean travel distance, so that s
• Guadagnini, A., & Neuman, S. P. (1998). Deterministic solution of stochastic groundwater flow equations by nonlocal finite elements. International Conference on Computational Methods in Water Resources, CMWR, 1, 347-354.
  More info

  Abstract: We consider the effect of measuring randomly varying hydraulic conductivities K(x) on the prediction of groundwater flow in a bounded porous domain under uncertainty. Hydraulic head is governed by a stochastic Poisson equation subject to random source and boundary terms. We present a system of exact nonlocal deterministic equations for optimum unbiased predictors of these quantities and for measures of corresponding prediction errors. We then develop recursive approximations for these equations and solve them to leading order in the variance of ln K(x) by nonlocal Galerkin finite elements. Our results compare well with Monte Carlo simulations of mean uniform and convergent flows in media with large variance and arbitrary correlation.
• Hsu, K. -., & Neuman, S. P. (1998). Erratum: Second-order expressions for velocity moments in two- and three-dimensional statistically anisotropic media (Water Resources Research (1997) 33:4 (625-637)). Water Resources Research, 34(1), 155-.
• Liu, H. H., Molz, F. J., Neuman, S. P., & Federico, V. D. (1998). Comment on 'Scaling of random fields by means of truncated power variograms and associated spectra' by Vittorio Di Federico and Shlomo P. Neuman. Water Resources Research, 34(11), 3207-3210.
• Tartakovsky, D. M., & Neuman, S. P. (1998). Transient effective hydraulic conductivities under slowly and rapidly varying mean gradients in bounded three-dimensional random media. Water Resources Research, 34(1), 21-32.
  More info

  Abstract: We have shown elsewhere [Tartakovsky and Neuman, this issue (a)] that in randomly heterogeneous media, the ensemble mean transient flux is generally nonlocal in space-time and therefore non-Darcian. We have also shown [Tartakovsky and Neuman, this issue (b)] that there are special situations in which this flux can be localized so as to render it Darcian in real or transformed domains. Each such situation gives rise to an effective hydraulic conductivity which relates mean gradient to mean flux at any point in real or transformed space-time. In this paper we develop first-order analytical expressions for effective hydraulic conductivity under three-dimensional transient flow through a box-shaped domain due to a mean hydraulic gradient that varies slowly in space and time. When the mean gradient varies rapidly in time, the Laplace transform of the mean flux is local but its real-time equivalent includes a temporal convolution integral; we develop analytical expressions for the real-time kernel of this convolution integral. The box is embedded within a statistically homogeneous natural log hydraulic conductivity field that is Gaussian and exhibits an anisotropic exponential spatial correlation structure. By the effective hydraulic conductivity of a finite box in such a field we imply the ensemble mean (expected value) of all random equivalent conductivities that one could associate with the box under these conditions. We explore the influence of domain size, time, and statistical anisotropy on effective conductivity and include a simple new formula for its variation with statistical anisotropy ratio in an infinite domain under steady state.
• Tartakovsky, D. M., & Neuman, S. P. (1998). Transient flow in bounded randomly heterogeneous domains 1. Exact conditional moment equations and recursive approximations. Water Resources Research, 34(1), 1-12.
  More info

  Abstract: We consider the effect of measuring randomly varying local hydraulic conductivities K(x) on one's ability to predict transient flow within bounded domains, driven by random sources, initial head; and boundary conditions. Our aim is to allow optimum unbiased prediction of local hydraulic heads h(x, t) and Darcy fluxes q(x, t) by means of their ensemble moments, (c) and (c), conditioned on measurements of K(x). We show that these predictors satisfy a compact deterministic flow equation which contains a space-time integrodifferential 'residual flux' term. This term renders (c) nonlocal and non-Darcian so that the concept of effective hydraulic conductivity looses meaning in all but a few special cases. Instead, the residual flux contains kernels that constitute nonlocal parameters in space-time that are additionally conditional on hydraulic conductivity data and thus nonunique. The kernels include symmetric and nonsymmetric second-rank tensors as well as vectors. We also develop nonlocal equations for second conditional moments of head and flux which constitute measures of predictive uncertainty. The nonlocal expressions cannot be evaluated directly without either a closure approximation or high-resolution conditional Monte Carlo simulation. To render our theory workable, we develop recursive closure approximations for the moment equations through expansion in powers of a small parameter which represents the standard estimation error of natural log K(x). These approximations are valid to arbitrary order for either mildly heterogeneous or well-conditioned strongly heterogeneous media. They allow, in principle, evaluating the conditional moments numerically on relatively coarse grids, without upscaling, by standard methods such as finite elements.
• Tartakovsky, D. M., & Neuman, S. P. (1998). Transient flow in bounded randomly heterogeneous domains 2. Localization of conditional mean equations and temporal nonlocality effects. Water Resources Research, 34(1), 13-20.
  More info

  Abstract: In randomly heterogeneous porous media one cannot predict flow behavior with certainty. One can, however, render optimum unbiased predictions of such behavior by means of conditional ensemble mean hydraulic heads and fluxes. We have shown in paper 1 [Tartakovsky and Neuman, this issue] that under transient flow, these optimum predictors are governed by nonlocal equations. In particular, the conditional mean flux is generally nonlocal in space-time and therefore non-Darcian. As such, it cannot be associated with an effective hydraulic conductivity except in special cases. Here we explore analytically situations under which localization is possible so that Darcy's law applies in real, Laplace, and/or infinite Fourier transformed spaces, approximately or exactly, with or without conditioning. We show that the corresponding conditional effective hydraulic conductivity tensor is generally nonsymmetric. An alternative to Darcy's law in each case, valid under mean no-flow conditions along Neumann boundaries, is a quasi-Darcian form that includes only a symmetric tensor which, however, does not constitute a bona fide effective hydraulic conductivity. Both lack of symmetry and differences between Darcian and quasi-Darcian forms disappear to first (but not necessarily higher) order of approximation in the (conditional) variance of natural log hydraulic conductivity. We adopt such an approximation to investigate analytically the effect of temporal nonlocality on one- and three-dimensional mean flows in infinite, statistically homogeneous media. Our results show that temporal nonlocality may manifest itself under either monotonic or oscillatory time variations in the mean hydraulic gradient. The effect of temporal nonlocality increases with the variance of log hydraulic conductivity and is more pronounced in one dimension than in three.
• Zimmerman, D. A., Marsily, G. D., Gotway, C. A., Marietta, M. G., Axness, C. L., Beauheim, R. L., Bras, R. L., Carrera, J., Dagan, G., Davies, P. B., Gallegos, D. P., Galli, A., Gómez-Hernández, J., Grindrod, P., Gutjahr, A. L., Kitanidis, P. K., Lavenue, A. M., McLaughlin, D., Neuman, S. P., , RamaRao, B. S., et al. (1998). A comparison of seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow. Water Resources Research, 34(6), 1373-1413.
  More info

  Abstract: 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.
• Federico, V. D., & Neuman, S. P. (1997). Scaling of random fields by means of truncated power variograms and associated spectra. Water Resources Research, 33(5), 1075-1085.
  More info

  Abstract: An interpretation is offered for the observation that the log hydraulic conductivity of geologic media often appears to be statistically homogeneous but with variance and integral scale which grow with domain size. We first demonstrate that the power (semi)variogram and associated spectra of random fields, having homogeneous isotropic increments, can be constructed as weighted integrals from zero to infinity (an infinite hierarchy) of exponential or Gaussian variograms and spectra of mutually oncorrelated homogeneous isotropic fields (modes). We then analyze the effect of filtering out (truncating) high- and low-frequency modes from this infinite hierarchy in the real and spectral domains. A low-frequency cutoff renders the truncated hierarchy homogeneous with an autocovariance function that varies monotonically with separation distance in a manner not too dissimilar than that of its constituent modes. The integral scales of the lowest- and highest-frequency modes (cutoffs) are related, respectively, to the length scales of the sampling window (domain) and data support (sample volume). Taking each relationship to be one of proportionality renders our expressions for the integral scale and variance of a truncated field dependent on window and support scales in a manner consistent with observations. The traditional approach of truncating power spectral densities yields autocovariance functions that oscillate about zero with finite (in one and two dimensions) or vanishing (in one dimension) integral scales. Our hierarchical theory allows bridging across scales at a specific locale, by calibrating a truncated variogram model to data observed on a given support in one domain and predicting the autocovariance structure of the corresponding multiscale field in domains that are either smaller or larger. One may also venture (we suspect with less predictive power) to bridge across both domain scales and locales by adopting generalized variogram parameters derived on the basis of juxtaposed hydraulic and tracer data from many sites.
• Hsu, K. -., Zhang, D., & Neuman, S. P. (1997). Erratum: Higher-order effects on flow and transport in randomly heterogeneous porous media (Water Resources Research (1996) 32(3) (571-582)). Water Resources Research, 33(4), 911-.
• Hsu, K., & Neuman, S. P. (1997). Second-order expressions for velocity moments in two- and three-dimensional statistically anisotropic media. Water Resources Research, 33(4), 625-637.
  More info

  Abstract: Second-order expressions are derived for the mean and covariance of steady state seepage velocity under mean uniform flows in infinite two- and three-dimensional domains. The order of approximation is defined in terms of the variance σ2 of a statistically homogeneous and anisotropic natural log hydraulic conductivity field Y with a Gaussian spatial autocorrelation function. Results show that second-order mean velocity either exceeds or is close to its first-order counterpart, depending on anisotropy. Head fluctuations of order larger than σ affect second-order velocity moments to the same extent as do head fluctuations of order σ in virtually all cases, hence neglecting the former renders the results nonasymptotic. Velocity variances are generally larger when approximated consistently to second than to first order. The ratio between second- and first-order variance approximations is larger in three than in two dimensions, larger for transverse than for longitudinal velocity, and increases with σ2. Anisotropy has a significant effect on second-order velocity variance. Second-order effects have the greatest influence on longitudinal velocity variance at extreme anisotropy ratios and on transverse velocity variance in isotropic domains.
• Chen, G., & Neuman, S. P. (1996). Wetting front instability in randomly stratified soils. Physics of Fluids, 8(2), 353-369.
  More info

  Abstract: A probabilistic criterion is derived for the onset of wetting front instability during surface water infiltration into a randomly stratified soil. It is based on the common assumption that the natural log hydraulic conductivity of the soil is a random, multivariate Gaussian function of space. Whereas the mean (expectation) of this function may exhibit a drift, its fluctuations about the mean are statistically homogeneous with constant variance and autocorrelation scale. The wetting front is taken to form a sharp boundary. Closed-form expressions for the probability of instability, and for the mean critical wave number, are obtained either directly or via a first-order reliability method. Monte Carlo simulations are used to verify these analytical solutions as well as to determine the mean maximum rate of incipient finger growth and corresponding mean wave number. The effects of applied pressure gradient, capillary pressure head at the wetting front, and statistical parameters of the hydraulic conductivity field on instability and incipient finger growth are investigated for a wide range of these variables. © 1996 American Institute of Physics.
• Hsu, K., Zhang, D., & Neuman, S. P. (1996). Higher-order effects on flow and transport in randomly heterogeneous porous media. Water Resources Research, 32(3), 571-582.
  More info

  Abstract: A higher-order theory is presented for steady state, mean uniform saturated flow and nonreactive solute transport in a random, statistically homogeneous natural log hydraulic conductivity field Y. General integral expressions are derived for the spatial covariance of fluid velocity to second order in the variance σ2 of Y in two and three dimensions. Integrals involving first-order (in σ) fluctuations in hydraulic head are evaluated analytically for a statistically isotropic two-dimensional Y field with an exponential autocovariance. Integrals involving higher-order head fluctuations are evaluated numerically for this same field. Complete second- order results are presented graphically for σ2 = 1 and σ2 = 2. They show that terms involving higher-order head fluctuations are as important as those involving lower-order ones. The velocity variance is larger when approximated to second than to first order in σ2. Discrepancies between second- and first-order approximations of the velocity autocovariance diminish rapidly with separation distance and are very small beyond two integral scales. Transport requires approximation at two levels: the flow level at which velocity statistics are related to those of Y, and the advection level at which macrodispersivities are related to velocity fluctuations. Our results show that a second-order flow correction affects transport to a greater extent than does a second-order correction to advection. Asymptotically. the second-order transverse macrodispersivity tends to zero as does its first- order counterpart. An approximation of advection alone based on Corrsin's conjecture, coupled with either a first- or a second-order flow approximation, leads to a transverse macrodispersivity which is significantly larger than that obtained by standard perturbation and tends to a nonzero asymptote. Published Monte Carlo results yield macrodispersivities that lie significantly below those predicted by first- and second-order theories. Considering that Monte Carlo simulations often suffer from sampling and computational errors, that standard perturbation approximations are theoretically valid only for σ2 < 1, and that Corrsin's conjecture represents the leading term in a renormalization group perturbation which contains contributions from an infinite number of high-order terms, we find it difficult to tell which of these approximations is closest to representing transport in strongly heterogeneous media with σ2 ≤ 1.
• Neuman, S. P., & Chen, G. (1996). On instability during immiscible displacement in fractures and porous media. Water Resources Research, 32(6), 1891-1894.
  More info

  Abstract: Existing linear theories for the onset of instability during immiscible fluid displacement in fractures and porous media are generalized to account for arbitrary fracture, mcan flow, and fluid interface orientations. The generalized analysis shows that contrary to common belief, the interface is not perpendicular (and hence incipient fingers do not extend in parallel) to the mean direction of displacement except in special cases.
• Neuman, S. P., Orr, S., Tartakovsky, D., Wallstrom, T. C., & Winter, C. L. (1996). Erratum: Prediction of steady state flow in nonuniform geologic media by conditional moments: Exact nonlocal formalism, effective conductivities, and weak approximation (Water Resources Research (1993) 29:2 (341- 364)). Water Resources Research, 32(5), 1479-.
• Paleologos, E. K., Neuman, S. P., & Tartakovsky, D. (1996). Effective hydraulic conductivity of bounded, strongly heterogeneous porous media. Water Resources Research, 32(5), 1333-1341.
  More info

  Abstract: We develop analytical expressions for the effective hydraulic conductivity K(e) of a three-dimensional, heterogeneous porous medium in the presence of randomly prescribed head and flux boundaries. The log hydraulic conductivity Y forms a Gaussian, statistically homogeneous and anisotropic random field with an exponential autocovariance. By effective hydraulic conductivity of a finite volume in such a field, we mean the ensemble mean (expected value) of all random equivalent conductivities that one could associate with a similar volume under uniform mean flow. We start by deriving a first-order approximation of an exact expression developed in 1993 by Neuman and Orr. We then generalize this to strongly heterogeneous media by invoking the Landau-Lifshitz conjecture. Upon evaluating our expressions, we find that K(e) decreases rapidly from the arithmetic mean K(A) toward an asymptotic value as distance between the prescribed head boundaries increases from zero to about eight integral scales of Y. The more heterogeneous is the medium, the larger is K(e) relative to its asymptote at any given separation distance. Our theory compares well with published results of spatially power-averaged expressions and with a first- order expression developed intuitively by Kitanidis in 1990.
• Zhang, D., & Neuman, S. P. (1996). Effect of local dispersion on solute transport in randomly heterogeneous media. Water Resources Research, 32(9), 2715-2723.
  More info

  Abstract: Recently, an exact Eulerian-Lagrangian theory of advective transport in spacetime random velocity fields was developed by one of us. We present a formal extension of this theory that accounts for anisotropic local dispersion. The resultant (conditional) mean transport equation is generally nonlocal in space-time. To assess the impact of local dispersion on the prediction of transport under uncertainty, we adopt a first-order pseudoFickian approximation for this equation. We then solve it numerically by Galerkin finite elements for two-dimensional transport from an instantaneous square source in a uniform (unconditional) mean flow field subject to isotropic local dispersion. We use a higher-order approximation to compute explicitly the standard deviation and coefficient of variation of the predicted concentrations. Our theory shows (in an exact manner), and our numerical results demonstrate (under the above closure approximations), that the effect of local dispersion on first and second concentration moments varies monotonically with the magnitude of the local dispersion coefficient. When this coefficient is small relative to macrodispersion, its effect on the prediction of nonreactive transport under uncertainty can, for all practical purposes, be disregarded. This is contrary to some recent assertions in the literature that local dispersion must always be taken into account, no matter how small.
• Zhang, D., & Neuman, S. P. (1996). Head and velocity covariances under quasi-steady state flow and their effects on advective transport. Water Resources Research, 32(1), 77-83.
  More info

  Abstract: Most stochastic models of solute transport assume flow to be at steady state. However, even locally uniform gradients tend to show seasonal fluctuations in magnitude and direction. Such transients affect the prediction of flow and plume migration and spread. The question is how and to what extent. We address this question by developing low-order approximations for autocovariances and cross covariances of velocity, head, and log hydraulic conductivity under quasi-steady state flow, then using a first-order Lagrangian approach to examine their effect on advective transport. Our results show that whereas periodic temporal fluctuations in the magnitude of the mean velocity may either enhance or reduce dispersion, similar fluctuations in its direction always cause longitudinal dispersion to decrease and transverse dispersion to increase.
• Neuman, S. P. (1995). On advective transport in fractal permeability and velocity fields. Water Resources Research, 31(6), 1455-1460.
  More info

  Abstract: Advective transport is considered in a steady state random velocity field with homogeneous increments. Such a field is self-affine with a power law dyadic semivariogram γ(s) proportional to d2ω, where d is distance and ω is a Hurst coefficient. It is characterized by a fractal dimension D = E + 1 - ω, where E is topological dimension. As the mean and variance of such a field are undefined, they are conditioned on measurement at some point x0. A tracer is introduced at another point y0 and elementary theoretical considerations are invoked to demonstrate that its conditional mean dispersion is local at all times. Its conditional mean concentration and variance are given explicitly by well-established expressions which, have not been previously recognized as being valid in fractal fields. -from Author
• Zhang, D., & Neuman, S. P. (1995). A stochastically-derived, deterministic model of solute transport conditioned on hydraulic data. Models for assessing and monitoring groundwater quality. Proc. symposium, Boulder, 1995, 227, 185-191.
  More info

  Abstract: In this paper we describe a stochastically-derived, deterministic model of solute transport in randomly heterogeneous geologic media. The model allows for local dispersion, linear equilibrium sorption, and radioactive decay or biodegradation. It utilizes an analytical-numerical approach to deterministically predict the spacetime evolution of concentrations conditioned on measurements of hydraulic conductivities (transmissivities) and /or hydraulic heads. We solve the conditional transport problem analytically at early time, and use a pseudo-Fickian approximation at later time. We illustrate the model by considering instantaneous point sources in a two-dimensional domain with a mildly fluctuating, statistically homogeneous, lognormal transmissivity field. Though we take the unconditional mean velocity to be uniform, conditioning on log transmissivity and hydraulic head data renders the velocity field statistically nonhomogeneous with reduced variances and correlation scales, renders the predicted plume irregular and non-Gaussian, and generally reduces both predictive dispersion and uncertainty.
• Zhang, D., & Neuman, S. P. (1995). Eulerian-Lagrangian analysis of transport conditioned on hydraulic data 1. Analytical-numerical approach. Water Resources Research, 31(1), 39-51.
  More info

  Abstract: The special case of steady state flow in a mildly fluctuating, statistically homogeneous, lognormal hydraulic conductivity field is considered. The unconditional mean velocity is taken to be uniform but conditioning on measurements of log hydraulic conductivity (or transmissivity) and/or hydraulic head is allowed. This renders the velocity field nonstationary. The conditional transport problem is solved analytically and expressed later in pseudo-Fickian form. The deterministic pseudo-Fickian equations involved a conditional, space-time dependent dispersion tensor which is evaluated numerically along mean "partial' trajectories. These equations lend themselves to accurate solution by standard Galerkin finite elements onto a relatively coarse grid. The final step is an explicit numerical computation of lower bounds on conditional concentration prediction variance-covariance (and coefficeint of variation), travel time distribution, cumulative mass release across a "compliance surface', the associated error, and plume spatial moments. -from Authors
• Zhang, D., & Neuman, S. P. (1995). Eulerian-Lagrangian analysis of transport conditioned on hydraulic data 2. Effects of log transmissitivity and hydraulic head measurements. Water Resources Research, 31(1), 53-63.
  More info

  Abstract: Velocity is conditioned on log transmissivity and/or hydraulic head data via cokriging. An early time analytical solution is combined with a pseudo-Fickian Galerkin finite element scheme for later time to obtain conditional predictions of concentration and lower bounds on its variance and coefficient of variation. The pseudo-Fickian scheme involves a conditional dispersion tensor which depends on information and varies in space-time. Hence the predicted plumes travel along curved trajectories and attain irregular, non-Gaussian shapes. A measure of uncertainty is computed for the original source location of a solute "particle' whose position at some late time is known from sampling. Spatial maps of this "particle origin covariance' provide vivid images of preferential flow paths and exclusion zones identified by the available data. These concepts are illustrated and results on instantaneous point and area sources are given. -from Authors
• Zhang, D., & Neuman, S. P. (1995). Eulerian-Lagrangian analysis of transport conditioned on hydraulic data 3. Spatial moments, travel time distribution, mass flow rate, and cumulative release across a compliance surface. Water Resources Research, 31(1), 65-75.
  More info

  Abstract: In paper two the effects of log transmissivity and hydraulic head data on conditional predictions of concentration due to instantaneous point and nonpoint sources are examined. In this paper the same approach is used directly to estimate mass flow rate across a "compliance surface', cumulative mass release, and the probability distribution of travel times across this surface and the associated uncertainty. Contrary to some other methods in the literature, the approach requires neither a special theory for travel times nor a prior assumption about their probability distribution. It is shown how one can compute explicitly the second spatial moment of the conditional mean plume about its center of mass, the conditional mean second spatial moment of the actual plume about its center of mass, and the conditional covariance of the plume center of mass. These quantities and the effect of conditioning on some of them are illustrated by considering instantaneous point, line, and area sources in a two-dimensional, statistically homogeneous, isotropic, mildly varying log transmissivity field under uniform prior (unconditional) mean flow. -from Authors
• Zhang, D., & Neuman, S. P. (1995). Eulerian-Langrangian analysis of transport conditioned on hydraulic data 4. Uncertain initial plume state and non-Gaussian velocities. Water Resources Research, 31(1), 77-88.
  More info

  Abstract: The analysis to date assumed that the initial state of the plume (the source term or initial concentration) is known with certainty and that the groundwater velocity field is Gaussian. In reality, the initial state of the plume is almost never known with certainty, especially when this state is inferred by sampling a plume at some arbitrary time t0 following introduction of the solute into the subsurface (as is the case at many contaminated sites). Likewise, Monte Carlo simulations have shown that in Gaussian log hydraulic conducitivity or log transmissivity fields, of the kind often indicated by in situ hydraulic test data, the longitudinal velocity becomes rapidly lognormal as the variance of the heterogeneities increases. In this paper it is shown how to handle both of these complications by means of an analytical-numerical method of computation. By revisiting previous examples, the effects that uncertainty in plume initial state and non-Gaussian velocities may have on predicted plume concentations and uncertainty due to instantaneous point and nonpoint sources with and without conditioning on measurements of log transmissivity or hydraulic head in two dimensions are explored. -from Authors
• Neuman, S. P. (1994). Generalized scaling of permeabilities: validation and effect of support scale. Geophysical Research Letters, 21(5), 349-352.
  More info

  Abstract: 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 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. The present paper cites evidence for such behavior over scales from 10 cm to 45 km based directly on permeability and transmissivity data. The paper then investigates theoretically the implications of such power-law behavior on the equivalent permeability of a block of rock having a characteristic length. -from Author
• Orr, S., & Neuman, S. P. (1994). Operator and integro-differential representations of conditional and unconditional stochastic subsurface flow. Stochastic Hydrology and Hydraulics, 8(2), 157-172.
  More info

  Abstract: Operator representations of stochastic subsurface flow equations allow writing their solutions implicitly or explicitly in terms of integro-differential expressions. Most of these representations involve Neumann series that must be truncated or otherwise approximated to become operational. It is often claimed that truncated Neumann series allow solving groundwater flow problems in the presence of arbitrarily large heterogeneities. Such claims have so far not been backed by convincing computational examples, and we present an analysis which suggests that they may not be justified on theoretical grounds. We describe an alternative operator representation due to Neuman and Orr (1993) which avoids the use of Neumann series yet accomplishes a similar purpose. It leads to a compact integro-differential form which provides considerable new insight into the nature of the solution. When written in terms of conditional moments, our new representation contains local and nonlocal effective parameters that depend on scale and information. As such, these parameters are not unique material properties but may change as more is learned about the flow system. © 1994 Springer-Verlag.
• Neuman, S. P. (1993). Eulerian-Lagrangian theory of transport in space-time nonstationary velocity fields: exact nonlocal formalism by conditional moments and weak approximation. Water Resources Research, 29(3), 633-645.
  More info

  Abstract: 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. -from Author
• Neuman, S. P., & Orr, S. (1993). Prediction of steady state flow in nonuniform geologic media by conditional moments: exact nonlocal formalism, effective conductivities, and weak approximation. Water Resources Research, 29(2), 341-364.
  More info

  Abstract: Considers the effect of measuring randomly varying local hydraulic conductivities K(x) on one's ability to predict steady state flow within a bounded domain, driven by random source and boundary functions. The prediction of local hydraulic head h(x) and Darcy flux q(x) are considered 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)]K 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. -from Authors
• Galeati, G., Gambolati, G., & Neuman, S. P. (1992). Coupled and partially coupled Eulerian-Lagrangian model of freshwater-seawater mixing. Water Resources Research, 28(1), 149-165.
  More info

  Abstract: The problem of density-dependent transport of salt in unconfined coastal aquifers is solved numerically by means of an implicit Eulerian-Lagrangian finite element formulation. The method has been used to study the effect of dewatering on seawater intrusion within a vertical cross section through an aquifer in southern Italy, related to the construction of a thermoelectric power plant. To investigate the extent to which the dependence of fluid density on salt concentration affects the numerical solution, the flow and advection-dispersion equations were solved in coupled (iterative), partially coupled (noniterative) and completely decoupled modes. Partial coupling was found to yield results very close to those obtained by full coupling but at great savings in computer time; the less rigorous decoupled approach led to results substantially different from those obtained through coupling and partial coupling. Effects of aquifer heterogeneity and the construction of a cutoff wall on seawater intrusion are discussed. -from Authors
• Neuman, S. P., Orr, S., Levin, O., & Paleologos, E. (1992). Theory and high-resolution finite element analysis of 2-D and 3-D effective permeabilities in strongly heterogeneous porous media. Finite Elements in Water Resources, Proceedings of the International Conference, 2, 117-136.
  More info

  Abstract: Starting from the premise that Darcy's law applies with a random hydraulic conductivity field K(x) defined on a support ω, we consider steady state flow in a larger domain Ω driven by statistically independent random source and boundary functions. Writing K(x) as the sum of a slowly varying unbiased estimate κ(x) and a zero mean estimation error K′(x), the conditional ensemble moments 〈h(x)〉κ≡〈h(x)|κ(x)〉 and 〈q(x)〉κ≡〈q(x)|κ(x)〉 become conditionally unbiased predictors of the hydraulic head h(x) and Darcy flux q(x), respectively; 〈q(x)〉κ satisfies a standard continuity equation driven by ensemble mean source and boundary functions; and 〈q(x)〉κ = -κ(x)▽〈h(x)〉κ+rκ(x) where rκ(x) is a residual flux. We present an exact formal expression for 〈r(x)〉κ which demonstrates that 〈q(x)〉κ is nonlocal (depends on head gradients at points other than x) with a well-defined kernel, and non-Darcian (there is no effective or equivalent conductivity valid for arbitrary directions of conditional mean flow), except in special cases. In some of these cases the effective conductivity is κ(x), in some it is a symmetric or a nonsymmetric second-rank tensor, and in yet other situations it exists only as a set of directional scalars but not as a tensor. We describe a weak approximation for 〈r(x)〉κ which improves with the quality of the estimate κ(x), and an ancillary approximation for the effective hydraulic conductivity tensor in an infinite domain. We then use high-resolution finite element Monte Carlo simulation to verify aspects of these approximations under uniform 3-D and radial 2-D flows in strongly heterogenous, statistically homogeneous and isotropic porous media.
• Vickers, B. C., Neuman, S. P., Sully, M. J., & Evans, D. D. (1992). Reconstruction and geostatistical analysis of multiscale fracture apertures in a large block of welded tuff. Geophysical Research Letters, 19(10), 1029-1032.
  More info

  Abstract: Reconstruction of fracture apertures in a block of welded tuff has revealed two scales of autocorrelation, one on the order of millimeters, the other on the order of tens of centimeters. The spatial and frequency distributions of apertures changed drastically as a result of minor nondestructive shear. -Authors
• Dagan, G., & Neuman, S. P. (1991). Nonasymptotic behavior of a common Eulerian approximation for transport in random velocity fields. Water Resources Research, 27(12), 3249-3256.
  More info

  Abstract: A particular Eulerian approximation which has become popular in recent years derives from the premise that, in mildly fluctuating velocity fields, terms involving the products of fluctuating quantities can be disregarded in comparison to terms which involve only one such quantity. This leads to a non-Fickian low-order approximation for the ensemble mean concentration (C), where the dispersive flux is not proportional to ∇(C), but is given instead by a convolution integral. The spatial moments of (C) based on this approximation are shown to be in conflict with those obtained from a first-order Lagrangian analysis. Explores this apparent contradiction and concludes that the reason lies in the nonasymptotic nature of the Eulerian approximation whereby terms neglected are of the same order as terms retained. -from Authors
• Wittmeyer, G., & Neuman, S. P. (1991). Monte Carlo experiments with robust estimation of aquifer model parameters. Advances in Water Resources, 14(5), 252-272.
  More info

  Abstract: Automatic calibration procedures for fitting groundwater models to measured data have almost exclusively used L2 error norms for estimation criteria. Quadratic estimation criteria result from applying the method of maximum likelihood under the assumption of jointly normal errors in the measured head data which, for the case of uncorrelated errors with constant variance, results in the ordinary least squares criterion. Quadratic estimation criteria are extremely sensitive to errors which violate the Gaussian assumption and pose an additional source of instability for the inverse problem. We consider the use of statistically robust estimation procedures to reduce these destabilizing effects. In particular, we test the robust M-estimator of Huber with a series of Monte Carlo experiments, compute the sample mean squared error, squared bias and variance of the estimated log-hydraulic conductivities for realizations of head measurements contaminated with mixed normal noise and compare these measures to their counterparts for the OLS estimator. We show that Huber's M-estimator achieves the same estimation variance as the OLS estimator with reduced need for prior information to regularize the problem. © 1991.
• Gambolati, G., Galeati, G., & Neuman, S. P. (1990). Eulerian-Lagrangian finite element model for coupled groundwater transport. Computational Methods in Subsurface Hydrology, 341-347.
  More info

  Abstract: The paper proposes a model for the coupled and partially uncoupled analysis of groundwater flow and salt dispersion and convection based on a Eulerian-Lagrangian approach for the contamination problem combined with a traditional finite element method for the flow problem. Unlike the schemes developed by Pinder and Cooper and Sanford and Konikow, the present approach implements a 'single step reverse tracking' technique similar to the one introduced by Neuman and Baptista et al.. The overall model is implemented in a FORTRAN code called CUNEO.
• Guzman, A. G., Wilson, L. G., Neuman, S. P., & Osborn, M. D. (1990). Simulating effect of channel changes on stream infiltration. Journal of Hydraulic Engineering, 115(12), 1631-1645.
  More info

  Abstract: We use a saturated-unsaturated flow model (UNSAT2) to study the effect of soil cementing and channel widening on infiltration and recharge in two idealized profiles along the Rillito River of southern Arizona. The channel sections associated with the profiles are 202-ft (61.6-m) and 326-ft (99.4m) wide. A synthetic 10-yr runoff in the river provides head values for generating infiltration before and after bank protection in these two sections. Results from the two base simulations are combined to produce four additional cases. Soil cementing the banks of the river without widening the channel reduces infiltration and recharge by 15% for the 326-ft wide (99.4-m) channel and by 25% for the 202-ft wide (61.6-m) channel. Increasing channel width together with bank protection increases infiltration 25% for our two sections. Increasing channel width without bank protection, however, increases infiltration by 46% for our conditions. Channel narrowing with bank protection markedly increases infiltration rates, about 50% for the two sections.
• Neuman, S. P. (1990). Adjoint Petrov-Galerkin method with optimum weight and interpolation functions defined on multi-dimensional nested grids. Computational Methods in Surface Hydrology, 347-356.
  More info

  Abstract: An adjoint Petrov-Galerkin method is proposed which solves multidimensional advection-dispersion, diffusion, and Poisson equations by means of weight (test) functions that form numerical solutions of adjoint state equations on a sequence of nested grids. The cells of these grids may have straight or curved boundaries. The weights are numerical relatives of 'optimal test functions' derived analytically by Celia et al. but possess important additional attributes. Solutions corresponding to a given member of the sequence are ideally suited for parallel computation. All solutions with the possible exception of those corresponding to the finest local grids are free of advective terms and contain diffusive terms only along intercell boundaries, a property which enhances accuracy and makes the method especially well suited for the treatment of nonuniform velocity fields in heterogeneous media. While the dependent variable is evaluated only at the nodes of the coarset grid, its spatial variation is described with a resolution corresponding to the finest grids. This is achieved in part by associating with each member of the sequence a novel family of 'optimal interpolation (trial) functions' which depend solely on the corresponding weight functions and reflect the general asymmetry of the latter. Both the weight and interpolation functions can vary with time to maintain optimality when the velocity and parameter fields are not constant.
• Neuman, S. P. (1990). Universal scaling of hydraulic conductivities and dispersivities in geologic media. Water Resources Research, 26(8), 1749-1758.
  More info

  Abstract: An interpretation is offered for the observation that dispersivities increase with scale. Apparent longitudinal dispersivity data from a variety of hydrogeologic settings are assumed to represent a continuous hierarchy of log hydraulic conductivity fields with mutually uncorrelated increments, each field having its own exponential autocovariance, associated integral scale, and variance that increases as a power of scale. Such a hierarchy is shown theoretically to form a self-similar random field with homogenous increments. Porous and fractured media appear to follow the same idealized scaling rule for both flow and transport, raising a question about the validity of many distinctions commonly drawn between such media. Finally, the data suggest that conditioning transport models through calibration against hydraulic measurements has the effect of filtering out large-scale modes from the hierarchy.
• Neuman, S. P., & Gardner, D. A. (1989). Determination of aquitard/aquiclude hydraulic properties from arbitrary water-level fluctuations by deconvolution. Ground Water, 27(1), 66-76.
  More info

  Abstract: A deconvolution method is proposed for the estimation of aquitard/aquiclude vertical hydraulic diffusivities on the basis of water-level measurements in piezometers completed both in the confining layer(s) and in neighboring aquifers. The water levels may fluctuate in an arbitrary manner which is not necessarily caused by a prescribed pumping regime (as required in conventional pumping tests). A simple formula from the literature can be used to correct the measured water levels for time lag in the piezometer response. If the specific storage of the confining material is known from laboratory consolidation tests, the hydraulic diffusivity estimates can be translated into vertical hydraulic conductivities. The method is illustrated with the aid of a case study from southern California and is shown to yield robust parameter estimates.
• Samper, F. J., & Neuman, S. P. (1989). Estimation of spatial covariance structures by adjoint state maximum likelihood cross validation, 1. Theory. Water Resources Research, 25(3), 351-362.
  More info

  Abstract: A series of three papers describes a cross-validation method to estimate the spatial covariance structure of intrinsic or nonintrinsic random functions from point or spatially averaged data that may be corrupted by noise. Any number of relevant parameters, including nugget effect, can be estimated. The theory, described in this paper, is based on a maximum likelihood approach which treats the cross-validation errors as Gaussian. Various a posteriori statistical tests are used to verify this hypothesis and to show that in many cases, correlation between these errors is weak. The log likelihood criterion is optimized through a combination of conjugate gradient algorithms. An adjoint state theory is used to efficiently calculate the gradient of the estimation criterion, optimize the step size downgradient, and compute a lower bound for the covariance matrix of the estimation errors. Issues related to the identifiability, stability, and uniqueness of the resulting adjoint state maximum likelihood cross-validation method are discussed. -from Authors
• Samper, F. J., & Neuman, S. P. (1989). Estimation of spatial covariance structures by adjoint state maximum likelihood cross validation, 2. Synthetic experiments. Water Resources Research, 25(3), 363-371.
  More info

  Abstract: Uses synthetic data to investigate the properties of the adjoint state maximum likelihood cross-validation (ASMLCV) method presented in paper 1. More than 40 synthetic experiments are performed to compare various conjugate gradient algorithms; investigate the manner in which computer time varies with ASMLCV parameters; study the effect of sample size and choice of kriging points on ASMLCV estimates; evaluate the ability of various model structure identification criteria to help select the most appropriate semivariogram model among given alternatives; study the conditions required for parameter identifiability, uniqueness, and stability; quantify the statistics of cross-validation errors; test hypotheses concerning the distribution and auto-correlation of these errors; and illustrate the computation of approximate quality indicators for ASMLCV parameter estimates. -from Authors
• Samper, F. J., & Neuman, S. P. (1989). Estimation of spatial covariance structures by adjoint state maximum likelihood cross validation, 3. Application to hydrochemical and isotopic data. Water Resources Research, 25(3), 373-384.
  More info

  Abstract: Presents applications of the adjoint state maximum likelihood cross-validation (ASMLCV) method to aquifer data. The Madrid basin in Spain serves as the source of information about 11 hydrochemical variables (pH, electrical conductivity, silica content, and the concentration of major ions) and two isotopes (oxygen 18 and carbon 14). With the exception of oxygen 18 and silica, the variables appear to be free of a horizontal drift. No discernible directional effects are seen. All variables exhibit a large nugget effect which is indicative of background noise. Concludes that more detailed and careful sampling in three dimensions is required if groundwater quality information is to become less prone to such noise and thereby more useful in the context of quantitative hydrogeological analyses. -from Authors
• Adar, E. M., & Neuman, S. P. (1988). Estimation of spatial recharge distribution using environmental isotopes and hydrochemical data, II. Application to Aravaipa Valley in Southern Arizona, U.S.A.. Journal of Hydrology, 97(3-4), 279-302.
  More info

  Abstract: In Part I* * Adar et al. (this volume). a mathematical model was developed to estimate the spatial distribution of annual recharge rates into an aquifer using environmental isotopes and hydrochemical data. In Part II the model is applied to real data from the Aravaipa Valley in southern Arizona. These data suffer from a paucity of hydrologic information but contain numerous isotopic and hydrochemical analyses. The model is able to extract from the available data both qualitative and quantitative information about recharge into, and flow through the shallow unconfined aquifer in the valley. The model shows that this aquifer is replenished primarily by lateral recharge from the pediments and the tributary Stowe Gulch Basin, as well as by upward leakage from the underlying confined aquifer. Recharge from the Aravaipa Creek is relatively minor. Lateral recharge from the eastern pediments into the unconfined aquifer exceeds that from the western pediments by a factor of two. The western pediments are a major source of recharge for the confined aquifer. These results are tentative because independent means to verify them are currently lacking. We are encouraged, however, that the mass balance errors of water and dissolved constituents do not exceed 4%, and that the flow rate obtained from Darcy's law at one location within the unconfined aquifer in which transmissivity is known, agrees with model prediction. While chemical reactions are ignored in the model, some of the isotopic and chemical species show sufficient spatial variability due to mixing to justify the approach. This is supported by analyses of equilibrium water-mineral interactions conducted by WATEQF. © 1988.
• Adar, E. M., Neuman, S. P., & Woolhiser, D. A. (1988). Estimation of spatial recharge distribution using environmental isotopes and hydrochemical data, I. Mathematical model and application to synthetic data. Journal of Hydrology, 97(3-4), 251-277.
  More info

  Abstract: A mathematical model is proposed to estimate the spatial distribution of annual recharge rates into an aquifer using environmental isotopes and hydrochemical data. The aquifer is divided into cells within which the isotopes and dissolved constituents are assumed to undergo complete mixing. For each mixing cell mass-balance equations expressing the conservation of water, isotopes and dissolved chemicals are written. These equations are solved simultaneously for unknown rates of recharge into the various cells by quadratic programming. The degree to which individual dissolved constituents may be considered conservative is tested a-priori by means of a chemical equilibrium model such as wateqf. Constituents which do not pass this test are either disregarded or suitably assigned a small weight in the quadratic program. In Part I, the model is applied to synthetic data corrupted by random noise and its sensitivity to input errors is examined. Part II* * Adar and Neuman (this volume). describes an application of the model to real data from the Aravaipa Valley in southern Arizona. © 1988.
• Cady, R., & Neuman, S. P. (1988). Three-Dimensional Adaptive Eulerian-Lagrangian Finite Element Method for Advection-Dispersion. Developments in Water Science, 36(C), 183-193.
  More info

  Abstract: The adaptive Eulerian-Lagrangian finite element method of Neuman' is improved and extended to three-dimensional domains. Accurate results are obtained for the entire range of Peclet numbers from very small to infinity and Courant numbers at least up to 15. Difficulties resolved include oscillations, numerical smearing, the clipping of peaks. and grid orientation effects. © 1988.
• Neuman, S. P., & Depner, J. S. (1988). Use of variable-scale pressure test data to estimate the log hydraulic conductivity covariance and dispersivity of fractured granites near Oracle, Arizona. Journal of Hydrology, 102(1-4), 475-501.
  More info

  Abstract: Application of modern stochastic theories to nonuniform geologic media requires, among other things, knowing the principal directions and correlation lengths of the log hydraulic conductivity covariance. To infer this information by standard geostatistical methods directly from point measurements of hydraulic conductivity in three-dimensional space is seldom possible due to a ubiquitous scarcity of data. To show how this difficulty can be overcome we use an example of fractured granitic rocks near Oracle, Arizona. At this site, hydraulic conductivities were measured on a small scale by packer tests and on a larger scale by cross-hole tests. The test results suggest that the measured hydraulic conductivities can be viewed as quantities defined over a continuum which are scalar on a local scale but exhibit large scale anisotropy due to spatial nonuniformity and the effect of fracture orientations. Based on this observation, the two sets of hydraulic conductivity data are related theoretically and the resulting expression is used to estimate the anisotropic covariance function of log hydraulic conductivities at the Oracle site. We then employ these parameters to calculate a Fickian dispersivity tensor for the Oracle granite by relying on the recent theory of Neuman et al. (1987a). Our results highlight the need to extend available stochastic theories to the domain of strongly nonuniform media the properties of which fluctuate with large amplitudes. © 1988.
• Farmer, I. W., Daemen, J. J., Desai, C. S., Glass, C. E., & Neuman, S. P. (1987). Rock mechanics. Proceedings of the 28th US symposium, Tucson, June-July 1987.. Rock mechanics. Proceedings of the 28th US symposium, Tucson, June-July 1987..
  More info

  Abstract: The symposium covered a wide range of applications with no specific theme but with more emphasis on fluid flow through rock masses and coupled flow; rock fracture and thermomechanical behaviour; rock characterization and remote techniques; field and laboratory testing; constitutive models; hydrofracture; well bore and tunnel stability; rock fragmentation and excavation; mine stability and rock mechanics in geologic repository design. Case histories and scale model studies are included. The ca.140 papers are arranged under 17 subject headings.-from Editors
• Neuman, S. P. (1987). ON METHODS OF DETERMINING SPECIFIC YIELD.. Ground Water, 25(6), 679-684.
  More info

  Abstract: The discussion relies on a pumping test conducted recently by Nwankwor and coworkers from the University of Waterloo in an unconfined aquifer at Borden, Ontario, Canada. A water-balance method promulgated by these workers and used earlier by others is invalid because it overlooks a major component of the water budget. This method suggests erroneously that specific yield increase with time during a pumping test. It is shown that these specific yield values are consistent with water-balance considerations when all the components of the water budget are properly taken into account. Whereas the larger specific yields usually obtained from laborator drainage experiments may be well suited for the evaluation of ground-water reserves that can potentially be recovered from an unconfined aquifer over long periods of time, they are not directly relevant to the problem of relating ground-water level variations to pumpage which is characterized by a shorter time scale. Additional results are discussed.
• Neuman, S. P., Samper, J., & Hernandez, M. (1987). Reply to comments by G. Gambolati and G. Galeati. Mathematical Geology, 19(3), 259-266.
• Neuman, S. P., Winter, C., & Newman, C. M. (1987). STOCHASTIC THEORY OF FIELD-SCALE FICKIAN DISPERSION IN ANISOTROPIC POROUS MEDIA.. Water Resources Research, 23(3), 453-466.
  More info

  Abstract: A three-dimensional theory is described for field-scale Fickian dispersion in anisotropic porous media due to the spatial variability of hydraulic conductivities. The study relies partly on earlier work by the authors the attributes of which are briefly reviewed. It leads to results which differ in important ways from earlier theoretical conclusions about dispersion in anisotropic media. We express the dispersion tensor D as the sum of a local component d and a field-scale component DELTA . The local component is assumed to be independent of velocity (which is most appropriate if it represents molecular diffusion) and its principal terms are taken to act parallel and normal to the mean velocity vector mu . The field-scale component is written as alpha mu , where alpha is a dispersivity tensor and mu equals vertical mu vertical . We show that at large Peclet numbers P, the dispersivity tensor reduces to a single principal component parallel to the mean velocity, regardless of how mu is oriented. This result, valid for arbitrary covariance functions of log-hydraulic conductivity, differs from that of L. F. Gelhar and C. L. Axness (1983), according to whom the asymptotic dispersivity tensor may possess more than one nonzero eigen value.
• Carrera, J., & Neuman, S. P. (1986). Estimation of aquifer parameters under transient and steady state conditions: 2. Uniqueness, stability, and solution algorithms.. Water Resources Research, 22(2), 211-227.
  More info

  Abstract: Under what conditions is the aquifer inverse problem well-posed? After defining the terms uniqueness, identifiability, and stability, theoretical considerations and synthetic examples are used to demonstrate that ill-posedness can be mitigated by including prior information about the parameters in the estimation criterion to be minimized. At the same time, the inclusion of such information is shown to be insufficient to guarantee uniqueness and stability in all cases. Several test problems in the recent literature, which have resulted in pessimistic conclusions about the solvability of the aquifer inverse problem, are shown to be ill-posed.-from Authors
• Carrera, J., & Neuman, S. P. (1986). Estimation of aquifer parameters under transient and steady state conditions: 3. Application to synthetic and field data.. Water Resources Research, 22(2), 228-242.
  More info

  Abstract: Illustrates and explores various features of the methodology proposed in papers 1 and 2 by applying it to a synthetic test case and to a set of field data from the southwestern United States.- from Authors
• Carrera, J., & Neuman, S. P. (1986). Estimation of aquifer parameters under transient and steady state conditions: I. Maximum likelihood method incorporating prior information.. Water Resources Research, 22(2), 199-210.
  More info

  Abstract: The parameters include values and directions of principal hydraulic conductivities (or transmissivities) in anisotropic media, specific storage (or storativity), interior and boundary recharge or leakage rates, coefficients of head-dependent interior and boundary sources, and boundary heads. Concludes by suggesting that certain model identification criteria developed in the time series literature, all of which are based on the maximum likelihood concept, might be useful for selecting the best groundwater model (or the best method of parameterizing a particular model) among a number of given alternatives.-from Authors
• Hsieh, P. A., & Neuman, S. P. (1985). Field determination of the three-dimensional hydraulic conductivity tensor of anisotropic media 1. Theory.. Water Resources Research, 21(11), 1655-1665.
  More info

  Abstract: The first paper presents theoretical expressions describing transient and steady state head response in monitoring intervals or arbitrary lengths and orientations, to constant-rate injection into (or withdrawal from) intervals having similar or different lengths and orientations. The conditions under which these intervals can be treated mathematically as points are investigated by an asymptotic analysis. The effect of planar no-flow and constant-head boundaries on the response is the theory of images. The second paper describes the field methodology and shows how the proposed approach works in the case of fractured granitic rocks.-from Authors
• Hsieh, P. A., Neuman, S. P., Stiles, G. K., & Simpson, E. S. (1985). Field determination of the three-dimensional hydraulic conductivity tensor of anisotropic media 2. Methodology and application to fractured rocks.. Water Resources Research, 21(11), 1667-1676.
  More info

  Abstract: The analytical solutions developed in the first paper can be used to interpret the results of cross-hole tests conducted in anisotropic porous or fractured media. Test results from a granitic rock near Oracle in southern Arizona are presented to illustrate how the method works for fractured rocks. At the site, the Oracle granite is shown to respond as a near-uniform, anisotropic medium, the hydraulic conductivity of which is strongly controlled by the orientations of major fracture sets. The cross-hole test results are shown to be consistent with the results of more than 100 single- hole packer tests conducted at the site. -from Authors
• Neuman, S. P., & Carrera, J. (1985). Maximum-likelihood adjoint-state finite-element estimation of groundwater parameters under steady- and nonsteady-state conditions. Applied Mathematics and Computation, 17(4), 405-432.
  More info

  Abstract: A method is presented for estimating the hydraulic parameters of groundwater flow models under steady- and nonsteady-state conditions. The estimation problem is posed in the framework of maximum-likelihood theory by means of a log-likelihood criterion that includes prior estimates of the parameters. To allow for an incomplete knowledge of the covariances of the prior head and parameter errors, these covariances are expressed in terms of a few unknown statistical parameters that may be estimated jointly with the hydraulic parameters. Computational efficiency is achieved by evaluating the gradient of the estimation criterion with an adjoint-state finite-element scheme and using a combination of conjugate-gradient algorithms, coupled with Newton's method for determining the step size to be taken at each iteration. Model structure identification criteria developed in the time-series literature (all of which utilize the maximum-likelihood concept) are shown to be useful for selecting the best way to parametrize a groundwater flow region when a number of alternative schemes of parametrization are given. The paper also demonstrates the potential utility of the proposed estimation method for the optimum design of space-time measurement networks. A case study dealing with three-dimensional flow in a multiaquifer system is briefly discussed. © 1985.
• Carrera, J., & Neuman, S. P. (1984). ADJOINT STATE FINITE ELEMENT ESTIMATION OF AQUIFER PARAMETERS UNDER STEADY-STATE AND TRANSIENT CONDITIONS.. Finite Elements in Water Resources, Proceedings of the International Conference, 199-211.
• Neuman, S. P. (1984). ADAPTIVE EULERIAN-LAGRANGIAN FINITE ELEMENT METHOD FOR ADVECTION-DISPERSION.. International Journal for Numerical Methods in Engineering, 20(2), 321-337.
  More info

  Abstract: A new adaptive finite element method is proposed for the advection-dispersion equation using an Eulerian-Lagrangian formulation. The method is based on a decomposition of the concentration field in two parts, one advective and one dispersive, in a rigorous manner that does not leave room for ambiguity. The advective component of steep concentration fronts is tracked forward with the aid of moving particles clustered around each front. Away from such fronts the advection problem is handled by an efficient modified method of characteristics called single-step reverse particle tracking. When a front dissipates with time, its forward tracking stops automatically and the corresponding cloud of particles is eliminated. The dispersion problem is solved by an unconventional Lagrangian finite element formulation on a fixed grid which involves only symmetric and diagonal matrices.
• Neuman, S. P., & Jacobson, E. A. (1984). Analysis of nonintrinsic spatial variability by residual kriging with application to regional groundwater levels. Journal of the International Association for Mathematical Geology, 16(5), 499-521.
  More info

  Abstract: 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. © 1984 Plenum Publishing Corporation.
• Neuman, S. P., Walter, G. R., Bentley, H. W., Ward, J. J., & Gonzalez, D. D. (1984). DETERMINATION OF HORIZONTAL AQUIFER ANISOTROPY WITH THREE WELLS.. Ground Water, 22(1), 66-72.
  More info

  Abstract: Existing methods for the determination of horizontal aquifer anisotropy by means of pumping tests require a minimum of four wells, one for water withdrawal and three for drawdown observations. This paper shows how the same methods can be used to determine anisotropy with as few as three wells, if at least two of them can be pumped in sequence. A field example is included. A method of analyzing data from more wells than the above minimum, by least squares, is also described.
• Winter, C. L., Newman, C. M., & Neuman, S. P. (1984). PERTURBATION EXPANSION FOR DIFFUSION IN A RANDOM VELOCITY FIELD.. SIAM Journal on Applied Mathematics, 44(2), 411-424.
  More info

  Abstract: A formal perturbation expansion is presented for the large-scale effective drift velocity and diffusion matrix of a medium with stationary random velocity, V(x), and constant nonrandom diffusion matrix, a, on a small scale. If mu denotes the mean of V(x), and V(x) is expressed as mu plus epsilon U(x), then the expansion is in powers of epsilon with mu , U and a fixed. Explicit expressions up to second order in epsilon are obtained which generalize the standard formulae to the case a does not equal 0 0.
• Neuman, S. P., & Adar, E. (1983). Estimation of recharge in a small southern Arizona basin by means of hydrological, hydrochemical, and environmental isotope data.. Array.
  More info

  Abstract: A mathematical model is proposed to estimate the ratios between various components of recharge into an alluvial aquifer system. The model is similar in principal to one developed earlier by Woolhiser, et al (1982) the river channels. It is based on long-term average mass balance expressions applied to water and various dissolved chemical constituents, coupled with a quadratic programming optimization scheme. The suitability of this model to the semiarid Aravaipa basin in southern Arizona is explored. Potential components of recharge into the basin are identified qualitatively on the basis of hydro-geological, hydrological, hydrochemical, and hydroisotope data. A method to obtain quantitative estimates of these components is outlined. Actual application of the quantitative model to the Aravaipa basin will be the subject of a future report.-from STAR, 22(20), 1984
• Clifton, P. M., & Neuman, S. P. (1982). EFFECTS OF KRIGING AND INVERSE MODELING ON CONDITIONAL SIMULATION OF THE AVRA VALLEY AQUIFER IN SOUTHERN ARIZONA.. Water Resources Research, 18(4), 1215-1234.
  More info

  Abstract: The Avra Valley aquifer in southern Arizona is modeled stochastically at three levels of uncertainty. The highest level of uncertainty occurs when log transmissivity estimates are based on measured values of this parameter but without regard to the geographic location of each measurement point. The resulting steady state hydraulic heads in the aquifer, computed by unconditional simulation with the aid of a multivariate normal random number generator coupled with a finite element model, have a relatively large variance. This variance can be reduced by conditioning the log transmissivity estimates on the spatial arrangement of the data by means of kriging.
• NEUMAN, S. P., PRELLER, C., & NARASIMHAN, T. N. (1982). ADAPTIVE EXPLICIT-IMPLICIT QUASI THREE-DIMENSIONAL FINITE ELEMENT MODEL OF FLOW AND SUBSIDENCE IN MULTIAQUIFER SYSTEMS.. Water Resource Res, V 18(N 5), 1551-1561.
  More info

  Abstract: A QUASI THREE-DIMENSIONAL FINITE ELEMENT MODEL IS PRESENTEDFOR THE ANALYSIS OF GROUNDWATER FLOW AND LAND SUBSIDENCE DUE TO PUMPAGE IN MULTIAQUIFER SYSTEMS. IN THE MODEL, AQUIFERSARE SIMULATED WITH THE AID OF TWO-DIMENSIONAL HORIZONTAL FINITE ELEMENT GRIDS. EACH AQUIFER IS CONNECTED TO ITS NEIGHBORS ABOVE AND BELOW BY ONE-DIMENSIONAL VERTICAL FINITE ELEMENT STRINGS WHICH ALLOW LEAKAGE TO TAKE PLACE ACROSS AQUITARDSAND AQUICLUDES. LAND SUBSIDENCE IS MODELED BY VARYING THE VOID RATIO OF EACH VERTICAL ELEMENT ACCORDING TO A NONLINEAR VERSION OF TERZAGHI'S ONE-DIMENSIONAL CONSOLIDATION THEORY.
• Neuman, S. P. (1982). Statistical characterization of aquifer heterogeneities: An overview. Special Paper of the Geological Society of America, 189, 81-102.
  More info

  Abstract: In recent years, ground-water hydrologists have become increasingly aware of the need to describe the spatial variability of aquifer characteristics in statistical terms. This has led to the development of new theoretical models whose parameters, representing the aquifer characteristics, are treated as stochastic variables rather than deterministic functions of space. In this paper, the state of the art in stochastic modeling is reviewed, and the experience and new knowledge gained with these models are summarized. Many of the stochastic models developed to date allow the parameters to fluctuate with equal amplitude at every point in space, including points at which the material properties have actually been measured. A more recent trend has been to try to reduce the variance of the computed hydraulic head values by conditioning the model on measured values of aquifer transmissivities. It is argued that a further reduction in this variance could be effected by conditioning the model not only upon measurements of the aquifer characteristics, but also upon historical data relating to the prevailing flow regime. This additional conditioning can be achieved by estimating the model parameters with the aid of inverse methods that are compatible with the stochastic interpretation of spatial variability. Two such inverse methods are described in this paper. It is suggested that in the future, the output from inverse models should be used as input into stochastic models of groundwater flow. In addition to the need for interphasing inverse models with stochastic models, additional research is required to improve the reliability, versatility, and computational efficiency of these models.
• Neuman, S. P., & Sorek, S. (1982). EULERIAN-LAGRANGIAN METHODS FOR ADVECTION-DISPERSION.. Finite Elements in Water Resources, Proceedings of the International Conference, 14. 41-14. 68.
• Neuman, S. P., Resnick, S. D., & Peebles, R. W. (1982). Developing a new deconvolution technique to model rainfall-runoff in arid environments.. Array.
  More info

  Abstract: Noisy rainfall runoff records of short duration are encountered in arid environments where rainfall often occurs in short isolated bursts and the data are measured with a considerable margin of error. Known noise was superimposed on synthetic rainfall runoff data to examine the ability of the Neuman-de Marsily deconvolution method to estimate the correct impulse response of the system when the data include only a single storm event. It is demonstrated that the deconvolution model leads to reliable estimates and can be used in the presence of realistic noise levels.-from STAR, 21(17), 1983
• Simpson, E. S., Neuman, S. P., Thompson, G. M., Hsieh, P., & Jones, J. (1982). HYDRAULIC AND TRACER EXPERIMENTS IN FRACTURED GRANITE FOR HIGH LEVEL WASTE SITE EVALUATION.. Proceedings of the Symposium on Waste Management, 371-379.
• Stephens, D. B., & Neuman, S. P. (1982). Free surface and saturated-unsaturated analyses of borehole infiltration tests above the water table. Advances in Water Resources, 5(2), 111-116.
  More info

  Abstract: Constant head borehole infiltration tests are widely used for the in situ evaluation of saturated hydraulic conductivities of unsaturated soils above the water table. The formulae employed in analysing the results of such tests disregard the fact that some of the infiltrating water may flow under unsaturated conditions. Instead, these formulae are based on various approximations of the classical free surface theory which treats the flow region as if it were fully saturated and enclosed within a distinct envelope, the so-called 'free surface'. A finite element model capable of solving free surface problems is used to examine the mathematical accuracy of the borehole infiltration formulae. The results show that in the hypothetical case where unsaturated flow does not exist, the approximate formulae are reasonably accurate within·a practical range of borehole conditions. To see what happens under conditions closer to those actually encountered in the field, the effect of unsaturated flow on borehole infiltration is investigated by means of two different numerical models: a mixed explicit-implicit finite element model, and a mixed explicit-implicit integrated finite difference model. Both of these models give nearly identical results; however, the integrated finite difference model is considerably faster than the finite element model. The relatively low computational efficiency of the finite element scheme is attributed to the large number of operations required in order to re-evaluate the conductivity (stiffness) matrix at each iteration in this highly non-linear saturated-unsaturated flow problem. The saturated-unsaturated analysis demonstrates that the classical free surface approach provides a distorted picture of the flow pattern in the soil. Contrary to what one would expect on the basis of this theory, only a finite region of the soil in the immediate vicinity of the borehole is saturated, whereas a significant percentage of the flow takes place under unsaturated conditions. As a consequence of disregarding unsaturated flow, the available formulae may underestimate the saturated hydraulic conductivity of fine grained soils by a factor of two, three, or more. Our saturated-unsaturated analysis leads to an improved design of borehole infiltration tests and a more accurate method for interpreting the results of such tests. The analysis also shows how one can predict the steady state rate of infiltration from data collected during the early transient period of the test. © 1982.
• Stephens, D. B., & Neuman, S. P. (1982). VADOSE ZONE PERMEABILITY TESTS: STEADY STATE RESULTS.. ASCE J Hydraul Div, 108(HY5), 640-659.
  More info

  Abstract: Free surface and saturated-unsaturated flow models are applied to the problem of constant head borehole infiltration tests for determining saturated hydraulic conductivity above a water table. Based on saturated-unsaturated simulations in four different homogeneous, isotropic soils, an empirical solution is developed for deep water table conditions to estimate saturated hydraulic conductivity with improved accuracy.
• Stephens, D. B., & Neuman, S. P. (1982). VADOSE ZONE PERMEABILITY TESTS: SUMMARY.. ASCE J Hydraul Div, 108(HY5), 623-639.
  More info

  Abstract: Solutions to determine saturated hydraulic conductivity above a water table by means of constant head borehole infiltration tests are summarized and critiqued. Analytical solutions are available for open and partly cased boreholes and deep and shallow water table conditions. The formulas are based on various simplifying assumptions petaining to soil properties and boundary conditions.
• Stephens, D. B., & Neuman, S. P. (1982). VADOSE ZONE PERMEABILITY TESTS: UNSTEADY FLOW.. ASCE J Hydraul Div, 108(HY5), 660-677.
  More info

  Abstract: Saturated-unsaturated flow models are used to simulate transient conditions in four soils during constant head borehole infiltration tests in the vadose zone. Simple graphical procedures are developed to overcome field problems associated with large water volume requirements or long time required to reach steady state by using early time data to estimate the steady infiltration rate.
• Winter, C. L., & Neuman, S. P. (1982). EVALUATION OF COEFFICIENTS IN THE MEGASCOPIC CONVECTION-DISPERSION EQUATION TO BE USED IN FAR FIELD MODELING OF HIGH LEVEL WASTE SITES.. Proceedings of the Symposium on Waste Management, 333-340.
• Neuman, S. P. (1981). A Eulerian-Lagrangian numerical scheme for the dispersion-convection equation using conjugate space-time grids. Journal of Computational Physics, 41(2), 270-294.
  More info

  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. © 1981.
• Neuman, S. P. (1980). Adjoint-state finite element equations for parameter estimation.. Preprints of papers. Volume 1. 3rd International Conference on finite elements in water resources, May 1980, Univ. of Mississsppi, Oxford, (Science and Education Administration, US Dept of Agriculture),, 2.66-2.75.
  More info

  Abstract: The equations have been used by Neuman (1980) in conjunction with a statistical parameter estimation procedure to determine the transmissivities of a hypothetical aquifer, and are currently being used to determine the transmissivities of the Cortaro Aquifer in Southern Arizona. -L.F.Musk
• Neuman, S. P. (1980). STATISTICAL APPROACH TO THE INVERSE PROBLEM OF AQUIFER HYDROLOGY - 3. IMPROVED SOLUTION METHOD AND ADDED PERSPECTIVE.. Water Resources Research, 16(2), 331-346.
  More info

  Abstract: A new statistically based approach to the problem of estimating spatially varying aquifer transmissivities on the basis of steady water level and flux data. The new method is based on a variational theory developed by G. Chavent which is extended to the case of generalized nonlinear least squares. The method is implemented numerically by a finite element scheme. The inverse problem is posed in terms of log transmissivities instead of transmissivities and is solved by a Fletcher-Reeves conjugate gradient algorithm in conjunction with Newton's method for determining the step size to be taken at each iteration. Two theoretical examples are included to demonstrate the ability of the new method to deal with artificial noise of a relatively large amplitude, derived from a given stochastic model.
• Neuman, S. P., Fogg, G. E., & Jacobson, E. A. (1980). STATISTICAL APPROACH TO THE INVERSE PROBLEM OF AQUIFER HYDROLOGY - 2. CASE STUDY.. Water Resources Research, 16(1), 33-58.
  More info

  Abstract: The paper describes a case study in which the new method is applied to actual field data from the Cortaro Basin in southern Arizona. The estimated transmissivities are shown to compare favorably with those obtained earlier for the same basin by an ad hoc trial-and-error procedure. Both sets of transmissivity values have been used successfully in conjunction with a mixed explicit-implicit finite element model to reproduce 25 years of water level variations. The finite element model has revealed some important features of the local hydrogeological regime.
• Stephens, D. B., & Neuman, S. P. (1980). Free surface and saturated-unsaturated analyses of borehole infiltration tests above the water table.. Preprints of papers. Volume 1. 3rd International Conference on finite elements in water resources, May 1980, Univ. of Mississippi, Oxford, (Science and Education Administration, US Dept of Agriculture), 2.229-2.238.
  More info

  Abstract: In this paper we will examine the accuracy of the Glover and Nasberg-Terletskata formulae under the assumption that the free surface concept is a valid representation of the real flow system. This will be don by comparing these formulae with results obtained from a finite model of free surface flow which does not suffer from the limitations of the analytical approach. -from Authors
• Fogg, G. E., Simpson, E. S., & Neuman, S. P. (1979). Aquifer modeling by numerical methods applied to an Arizona ground-water basin ( Tucson Basin).. Arizona University, Tucson, Department of Hydrology and Water Resources, Technical Report, 32.
  More info

  Abstract: FLUMP, a recently developed mixed explicit-implicit finite-element program, was calibrated against a data base obtained from a portion of the Tucson Basin aquifer, Arizona, and represents its first application to a real-world problem. -from US Govt Reports Announcements, 25, 1979
• Neuman, S. P. (1979). Perspective on 'delayed yield'.. WATER RESOUR. RES., 15(4 , Aug. 1979), 899-908.
  More info

  Abstract: The phenomenon of 'delayed yield' has received wide attention in the literature during recent years. Considerable work has been devoted to the development of physically based mathematical models capable of successfully reproducing observed aquifer behaviour. As a result of this effort, there appears to be at present a consensus of opinion among the main contributors to the theory of delayed yield as to an acceptable physical and mathematical interpretation of this phenomemon. In spite of this apparent consensus of opinion, there is still a considerable amount of controversy in the literature about some fundamental aspects of unconfined well hydraulics. The outstanding issues concern the physical assumptions underlying various delayed yield models, changes in the flow regime of an unconfined aquifer as a function of radial distance from the pumping well, the relationship between Boulton's ALPHA and other physical characteristics of the aquifer, the validity of conditions required in order to represent unconfined flow by integrodifferential equations, and the behaviour of unconfined aquifers at large values of time. The present paper offers a perspective from the author's viewpoint on these and other related issues. (A)
• Neuman, S. P., & Yakowitz, S. (1979). STATISTICAL APPROACH TO THE INVERSE PROBLEM OF AQUIFER HYDROLOGY - 1. THEORY.. Water Resources Research, 15(4), 845-860.
  More info

  Abstract: A new statistically based approach to the problem of estimating spatially varying aquifer transmissivities on the basis of steady state water level data is presented. The method involves solving a family of generalized nonlinear regression problems and then selecting one particular solution from this family by means of a comparative analysis of residuals. A linearized error analysis of the solution is included. This analysis allows one to estimate the convariance of the transmissivity estimates as well as the square error of the estimates of hydraulic head. An efficient explicit numerical scheme for solving the inverse problem in an approximate manner when errors in water level data are sufficiently small is also described.
• Neuman, S. P., & Yakowitz, S. J. (1979). Identification of aquifer parameters using numerical techniques. Completion report.. Array.
  More info

  Abstract: A new statistically based approach to the problem of estimating spatially varying aquifer transmissivities on the basis of steady state water level and flux data has been developed. This analysis allows one to estimate the covariance of the transmissivity estimates as well as the square error of the estimates of hydraulic head. -from US Govt Reports Announcements, 11, 1980
• Narasimhan, T. N., Neuman, S. P., & Witherspoon, P. A. (1978). FINITE ELEMENT METHOD FOR SUBSURFACE HYDROLOGY USING A MIXED EXPLICIT-IMPLICIT SCHEME.. Water Resources Research, 14(5), 863-877.
  More info

  Abstract: The mixed explicit-implicit Galerkin finite element method developed previously by the authors is shown to be suited for a wide class of problems arising in subsurface hydrology. These problems include confined saturated flow, unconfined flow under free surface conditions subject to the Dupuit assumption, flow in aquifers which are partly confined and partly unconfined, axisymmetric flow to a well with storage, and flow in saturated-unsaturated soils. Five examples are presented to demonstrate the versatility and power of this new approach. A purely physical derivation of the finite element equations which does not rely on the Galerkin formalism is also included.
• Dasberg, S., & Neuman, S. P. (1977). Peat hydrology in the Hula Basin, Israel: I. Properties of peat. Journal of Hydrology, 32(3-4), 219-239.
  More info

  Abstract: The results of laboratory and field studies concerned with various mechanical and hydraulic properties of peat in the Hula Basin are presented. Particular emphasis has been placed on the measurement of bulk density, porosity and void ratio, hydraulic conductivity, water retention under unsaturated conditions, compressibility, release of water due to consolidation and the content of organic matter. It is shown that the properties of peat in the basin change drastically when it becomes partially desaturated. As a result of this change, one can distinguish between a permanently saturated layer below the zone of water-table fluctuation having relatively uniform characteristics and an overlying partially desaturated layer the properties of which vary gradually with depth. The saturated layer is currently undergoing a process of consolidation and its properties may therefore be expected to change in the future. One of the conclusions reached in this work is that compressibility is very important and must not be neglected in dealing with fluid flow through peats. This is true for confined as well as unconfined (i.e., water-table) conditions. © 1977.
• Neuman, S. P. (1977). Theoretical derivation of Darcy's law. Acta Mechanica, 25(3-4), 153-170.
  More info

  Abstract: Darcy's 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 Darcy's 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 Darcy's law. © 1977 Springer-Verlag.
• Neuman, S. P., & Dasberg, S. (1977). Peat hydrology in the Hula Basin, Israel: II. Subsurface flow regime. Journal of Hydrology, 32(3-4), 241-256.
  More info

  Abstract: The subsurface flow regime in the drained portion of the Hula Basin has been investigated on a local and regional scale. The results are based on data obtained largely from a three-dimensional network of piezometers and observation wells covering the entire drained area, and another three-dimensional network covering a representative experimental plot. The data include measurements of piezometric levels under various conditions, analyses of stable and radioactive environmental isotopes in water samples from various depths and chemical analyses of water. It is shown that the position of the water table is controlled by the local topography. The groundwater at depths exceeding 2 m below the surface is subject to a mild regional hydraulic gradient from north to south, and a much more pronounced vertical gradient upward. The upwelling groundwater is probably derived from recharge over the surrounding mountains. When the groundwater reaches the water table it mixes with waters originating at the surface, discharges into the channels and ditches which cross the drained area, and is carried downstream into the Jordan River. One of the important conclusions of this study is that water originating at the surface has practically no chance of reaching a depth exceeding 2-3 m. © 1977.
• Neuman, S. P., & Krzysztofowicz, R. (1977). An iterative algorithm for interactive multiobjective programming. Advances in Water Resources, 1(1), 1-14.
  More info

  Abstract: A new iterative algorithm for interactive multiobjective programming is proposed. The algorithm is based on the Lagrange multiplier technique of generating noninferior solutions, and it is shown to converge under certain conditions. It reduces a complex multiobjective problem into a sequence of two-objective problems which the decision maker can handle more easily. The number of two-objective problems with which the decision maker is confronted, as well as the total number of noninferior solutions that must be generated, increase more or less linearly with the number of objectives. Computational efficiency is further enhanced by avoiding the need for regression. The decision maker interacts with the model directly in the functional space, and he is not required to translate his judgment of relative worth into numbers. Due to the iterative nature of the algorithm, the decision maker can articulate his preferences in a progressive manner. Furthermore, he may modify his attitude at any stage of the computation, based on partial results, without adversely affecting the quality of the solution. An example problem previously solved by other methods, including the surrogate worth trade-off approach, is used to illustrate the new algorithm. © 1977.
• Neuman, S. P., Narasimhan, T. N., & Edwards, A. L. (1977). MIXED EXPLICIT-IMPLICIT ITERATIVE FINITE ELEMENT SCHEME FOR DIFFUSION-TYPE PROBLEMS: I. THEORY.. International Journal for Numerical Methods in Engineering, 11(2), 309-323, 325.
  More info

  Abstract: A Galerkin finite element formulation of diffusion processes based on a diagonal capacity matrix is analysed from the standpoint of local stability and convergence. The theoretical analysis assumes that the conductance matrix is locally diagonally dominant, and it is shown that one can always construct a finite element network of linear triangles satisfying this condition. Time derivatives are replaced by finite differences, leading to a mixed explicit-implicit system of algebraic equations which can be efficiently solved by a point iterative technique. In this work the accelerated point iterative method is adopted and is shown to converge when the conductance matrix is locally diagonally dominant. Several examples are included in Part II of this paper to demonstrate the efficiency of the new approach.
• Neuman, S. P., Narasimhan, T. N., & Witherspoon, P. A. (1977). APPLICATION OF MIXED EXPLICIT-IMPLICIT FINITE ELEMENT METHOD TO NONLINEAR DIFFUSION-TYPE PROBLEMS.. Array, 1. 153-1. 186.
  More info

  Abstract: The point iterative nature of the solution process makes the mixed approach ideally suited for the treatment of quasilinear problems in which the conductances, capacitances, source terms, and boundary conditions vary with the dependent variable or with time. To demonstrate the accuracy and flexibility of the new approach, the paper considers three different types of problems with varying degrees of nonlinearity that are commonly encountered in subsurface hydrology; saturated unconfined flow governed by the Boussinesq equation; unsaturated flow governed by the diffusion equation; and saturated-unsaturated flow governed by the Richards equation.
• Neuman, S. P., & Marsily, G. d. (1976). IDENTIFICATION OF LINEAR SYSTEMS RESPONSE BY PARAMETRIC PROGRAMMING.. Water Resources Research, 12(2), 253-262.
  More info

  Abstract: The problem is formulated as a multicriterion decision process under uncertainty in a manner analogous to that previously described by one of the authors in connection with the inverse problem of groundwater hydrology. Parametric programming serves as a means of generating a continuous set of alternative solutions to the identification problem together with a bicriterion function representing these alternatives. The shape of this bicriterion curve is used as a guide in selecting a particular solution. The method is illustrated by two examples, and the results are compared with those obtained by another approach developed previously by one of the authors.
• Feddes, R. A., Kowalik, P., Neuman, S. P., & Bresler, E. (1975). FINITE DIFFERENCE AND FINITE ELEMENT SIMULATION OF FIELD WATER UPTAKE BY PLANTS.. Hydrol Sci BULL Sci Hydrol, 21(1), 81-98.
  More info

  Abstract: The problem of nonsteady flow of water in a soil-plant system can be described by adding a sink term to the continuity equation for soil water flow. The sink term is defined in two different ways. Firstly, it is considered to be dependent on the hydraulic conductivity of the soil, on the difference in pressure head between the soil and the root-soil interface and some root effectiveness function. Secondly, the sink is taken to be a prescribed function of the soil water content. The partial differential equation applying to the first problem is solved by both a finite difference and a finite element technique, that applying to the second problem by a finite difference approach. The purpose of this paper is to verify the numerical methods against field measurements, to compare the results obtained by the three numerical methods and to show how the finite element method can be applied to complex but realistic two-dimensional flow situations. Two examples are given.
• Feddes, R. A., Neuman, S. P., & Bresler, E. (1975). FINITE ELEMENT ANALYSIS OF TWO-DIMENSIONAL FLOW IN SOILS CONSIDERING WATER UPTAKE BY ROOTS - 2. FIELD APPLICATIONS.. Proc Soil Sci Soc Am, 39(2), 231-237.
  More info

  Abstract: The finite element model developed in Part I is verified against field data and the results are compared with those obtained from a one-dimensional finite difference model. Two field examples are given: vertical flow in a clay soil with red cabbage, and two-dimensional flow in a complex layered and anisotropic situation. The results indicate that the finite element model is well suited for the treatment of complex and realistic field problems.
• Neuman, S. P. (1975). ANALYSIS OF PUMPING TEST DATA FROM ANISOTROPIC UNCONFINED AQUIFERS CONSIDERING DELAYED GRAVITY RESPONSE.. Water Resources Research, 11(2), 329-342.
  More info

  Abstract: Previously a new analytical model was proposed by the author for the delayed response process characterizing flow to a well in an unconfined aquifer. The new approach differs from that of Boulton (1954, 1963, 1970, 1973) and Boulton and Pontin (1971) in that it is based only on well-defined physical parameters of the aquifer system. As such, it can be used to develop methods for determining the hydraulic properties of anisotropic unconfined aquifers from field drawdown data. Two methods of analysis are described, one based on the matching of field data with theoretical type curves and the other based on the semilogarithmic relationship between drawdown and time. Owing to the reversible nature of the delayed response process as represented by the analytical model, data from recovery tests can be used to determine aquifer transmissivity.
• Neuman, S. P. (1975). ROLE OF SUBJECTIVE VALUE JUDGEMENT IN PARAMETER IDENTIFICATION.. Array, 59-84.
  More info

  Abstract: A theoretical analysis of the inverse problem associated with a continuous model of groundwater flow indicates that, in general, the solution to this problem is nonunique and is extremely sensitive to minor errors in the model and the data. A similar difficulty is encountered when the flow equations are discretized in space by a numerical scheme such as the finite difference or finite element methods. In order to obtain a meaningful solution it is necessary to introduce additional criteria reflecting the modeler's subjective view point as to the most plausible shape of the solution surface. The absolute need to rely on such prior assumptions in dealing with the inverse problem means that information obtained from direct field and laboratory tests on the system to be simulated is a vital importance to the parameter identification process. A post-optimal sensitivity analysis is suggested as a means of ascertaining the reliability of each individual parameter obtained in this manner. Refs.
• Neuman, S. P., Feddes, R. A., & Bresler, E. (1975). FINITE ELEMENT ANALYSIS OF TWO-DIMENSIONAL FLOW IN SOILS CONSIDERING WATER UPTAKE BY ROOTS - 1. THEORY.. Proc Soil Sci Soc Am, 39(2), 224-230.
  More info

  Abstract: The problem of two-dimensional nonsteady flow of water in unsaturated and partly saturated porous media is solved by a Galerkin-type finite element approach. Special emphasis is placed on the treatment of atmospheric boundaries, water uptake by plants, irregularly shaped boundaries, soil heterogeneity, and soil anisotropy. The apparent advantages of this method over more conventional numerical techniques are discussed.
• Neuman, S. P., Feddes, R. A., & Bresler, E. (1975). GALERKIN METHOD OF SIMULATING WATER UPTAKE BY PLANTS.. Array, 325-348.
  More info

  Abstract: This particular finite element method is used to simulate the nonsteady flow of water in a multidimensional soil-plant environment. The approach can handle nonuniform soils having irregularly shaped boundaries and arbitrary degrees of local anisotropy. The soil can be unsaturated, saturated, or partly saturated. Flow can take place in a two-dimensional region or in a three-dimensional region having a radial symmetry. Particular emphasis is placed upon the simulation of atmospheric boundaries and water uptake by plant roots. The nonlinear atmospheric boundary conditions are handled by a unique iterative procedure relying on the ease with which fluxes normal to the boundary are assigned in the finite element approach. An example is included describing two-dimensional flow in a complex field situation. Refs.
• Feddes, R. A., Bresler, E., & Neuman, S. P. (1974). FIELD TEST OF A MODIFIED NUMERICAL MODEL FOR WATER UPTAKE BY ROOT SYSTEMS.. Water Resources Research, 10(6), 1199-1206.
  More info

  Abstract: Data obtained from careful water balance studies on water uptake by the roots of red cabbage are compared with results obtained from a modified numerical model of M. N. Nimah and R. J. Hanks. In the modified model the air dry moisture content at the soil surface may vary with time depending on meteorological conditions. The maximum possible rate of evapotranspiration is calculated by considering both meteorological conditions and crop properties. Data are quoted to suggest that the coefficient of the root sink may sometimes vary exponentially with depth. A period of 7 weeks was simulated, and the calculated weekly moisture profiles did not agree completely with those measured in the field. On the other hand, the calculated cumulative rates of evaporation and transpiration were in excellent agreement with the field data. When the original model was used without the suggested modifications, the agreement of these rates with the field data was not as good, an indication that some of these modifications actually improve the predictive capabilities of the model.
• NEUMAN, S. P. (1974). GALERKIN APPROACH TO UNSATURATED FLOW IN SOILS.. Array, 1974-.
  More info

  Abstract: AN ITERATIVE GALERKIN-TYPE FINITE ELEMENT METHOD IS EMPLOYED TO SOLVE THE QUASILINEAR EQUATIONS OF NON-STEADY WATER FLOW IN UNSATURATED AND PARTLY-SATURATED POROUS MEDIA.THE PARAMETERS IN THESE EQUATIONS AS WELL AS THE BOUNDARY CONDITIONS ALONG SEEPAGE FACES ARE HIGHLY NONLINEAR FUNCTIONS OF THE DEPENDENT VARIABLE.OWING TO THIS STRONG NONLINEARITY, PAST ATTEMPT TO SOLVE SUCH PROBLEMS WITH FINITE DIFFERENCE METHODS WERE NEVER COMPLETELY SUCCESSFUL.THE FINITE ELEMENT METHOD APPEARS TO BE COMPUTATIONALLY MORE EFFICIENT THAN THE CONVENTIONAL FINITE DIFFERENCE APPROACH IN OVERCOMING DIFFICULTIES ARISING FROM THE NONLINEAR NATURE OF THE PROBLEM.EXPERIENCE WITH THE FINITE ELEMENT ALGORITHM SHOWS THAT CONVERGENCE IS OFTEN NEAR QUADRATIC.THIS ALGORITHM IS CAPABLE OF HANDLING FLOW IN NONUNIFORM SOILS HAVING COMPLEX BOUNDARIES, ARBITRARY ANISOTROPY, AND AN UNLIMITED NUMBER OF SEEPAGE FACES.FLOW TO A FULLY OR PARTIALLY PENETRATING WELL OF A FINITE DIAMETER TAKING INTO ACCOUNT WELL STORAGE AS WELL AS PUMP CHARACTERISTICS, CAN ALSO BE TREATED.THE RESULTS INDICATE THAT THE CLASSICAL CONCEPT OF A 'FREE SURFACE' IS UNSUITABLE FOR MOST PROBLEMS INVOLVING TRANSIENT SEEPAGE THROUGH POROUS MEDIA.(A).
• Neuman, S. P. (1974). EFFECT OF PARTIAL PENETRATION ON FLOW IN UNCONFINED AQUIFERS CONSIDERING DELAYED GRAVITY RESPONSE.. Water Resources Research, 10(2), 303-312.
  More info

  Abstract: A theory is advances to account for the effect of a well partially penetrating a homogeneous anisotropic unconfined aquifer. Field and laboratory evidence is quoted to suggest that the elastic storage properties of unconfined aquifers may often be much more pronounced that those of deep-seated confined aquifers composed of similar materials. The need for further research into the mechanical properties of unconfined aquifers and the role of unsaturated flow in the response process is emphasized.
• Neuman, S. P. (1974). GALERKIN APPROACH TO UNSATURATED FLOW IN SOILS.. Array, 517-522.
  More info

  Abstract: An iterative Galerkin-type finite element method is employed to solve the quasilinear equations of nonsteady water flow in unsaturated and partly-saturated porous media. The parameters in these equations as well as the boundary conditions along seepage faces are highly nonlinear functions of the dependent variable. The finite element method appears to be computationally more efficient than the conventional finite difference approach in overcoming difficulties arising from the nonlinear nature of the problem. Experience with the finite element algorithm shows that convergence is often near quadratic. This algorithm is capable of handling flow in nonuniform soils having complex boundaries, arbitrary anisotropy, and an unlimited number of seepage faces. Flow to a fully or partially penetrating well of a finite diameter, taking into account well storage as well as pump characteristics, can also be treated.
• NEUMAN, S. P. (1973). SATURATED - UNSATURATED SEEPAGE BY FINITE ELEMENTS.. J. HYDRAUL. DIV., PROC., ASCE., 99(HY12), DECEMBER, 1973.
  More info

  Abstract: A GALERKIN-TYPE FINITE ELEMENT METHOD IS EMPLOYED TO SOLVE THE QUASILINEAR PARTIAL DIFFERENTIAL EQUATIONS OF TRANSIENT SEEPAGE IN SATURATED-UNSATURATED POROUS MEDIA.THE RESULTING COMPUTER PROGRAM IS CAPABLE OF HANDLING NONUNIFORM FLOW REGIONS HAVING COMPLEX BOUNDARIES AND ARBITRARY DEGREES OF LOCAL ANISOTROPY.FLOW CAN TAKE PLACE IN A VERTICAL PLANE, IN A HORIZONTAL PLANE, OR IN A THREE-DIMENSIONAL SYSTEM WITH RADIAL SYMMETRY.AN ARBITRARY NUMBER OF SEEPAGE FACES CAN BE CONSIDERED SIMULTANEOUSLY, AND THE POSITIONS OF THE EXIT POINTS ON THESE BOUNDARIES ARE ADJUSTED AUTOMATICALLY DURING EACH TIME STEP.TWO EXAMPLES, ONE OF SEEPAGE THROUGH AN EARTH DAM WITH A SLOPING CORE AND HORIZONTAL DRAINAGE BLANKET, AND THE OTHER OF SEEPAGE THROUGH A LAYERED MEDIUM CUT BY A COMPLEX TOPOGRAPHY, ARE INCLUDED.THESE EXAMPLES INDICATE THAT THE CLASSICAL CONCEPT OF A FREE SURFACE IS NOT ALWAYS APPLICABLE WHEN DEALING WITH TRANSIENT SEEPAGE THROUGH SOILS.(A)
• Neuman, S. P. (1973). SATURATED-UNSATURATED SEEPAGE BY FINITE ELEMENTS.. ASCE J Hydraul Div, 99(HY12), 2233-.
  More info

  Abstract: A Galerkin-type finite element method is employed to solve the quasilinear partial differential equations of transient seepage in saturated-unsaturated porous media. The resulting computer program is capable of handling nonuniform flow regions having complex boundaries and arbitrary degrees of local anisotropy. Flow can take place in a vertical plane, in a horizontal plane, or in a three-dimensional system with radial symmetry.
• Neuman, S. P., & Witherspoon, P. A. (1970). Comments on 'a correspondence principle for the theory of leaky aquifers' by Ismael Herrera and Germán E. Figueroa V.. Water Resources Research, 6(3), 1009-1010.
• Neuman, S. P., & Witherspoon, P. A. (1970). Finite element method analyzing steady seepage with a free surface. Water Resources Research, 6(3), 889-897.
  More info

  Abstract: A new iterative approach to steady seepage of ground water with a free surface has been developed using the finite element method. This approach eliminates a number of difficulties that were inherent in the iterative procedures previously u sed to solve this problem and rapid convergenceis now assured in all cases. The method is applicable to heterogeneous porous media with complex geometric boundaries and arbitrary degrees of anisotropy. It can handle problems where the free surface is discontinuous and where portions of the free surface are vertical or near vertical. In addition, infiltration or evapotranspiration at the free surface can be handled with ease. Several examples are included to demonstrate the power of this new approach and to show how it can apply to a wider variety of free surface problems than has been possible before.
• Neuman, S. P., & Witherspoon, P. A. (1970). Variational principles confined and unconfined flow of ground water. Water Resources Research, 6(5), 1376-1382.
  More info

  Abstract: Finite element techniques are being employed to an increasing degree in solving complicated groundwater flow problems, including problems with a free surface. These methods usually rely on the development of appropriate variational principles. This paper presents several such principles for both confined and unconfined flow.
• NEUMAN, S. P., & WITHERSPOON, P. A. (1969). APPLICABILITY OF CURRENT THEORIES OF FLOW IN LEAKY AQUIFERS. Water Resources Research, 5(4), 817-829.
  More info

  Abstract: Applicability of current theories of flow in leaky aquifers has been investigated using new analytic and numeric solutions. In general it has been customary to use one or both of the following assumptions; storage in the aquitards is negligible, and drawdowns in the unpumped aquifers remain zero. At small values of time the first assumption amounts to neglecting the effects of leakage completely. At large values of time, neglecting drawdowns in the unpumped aquifers can also lead to serious errors unless the transmissibilities of these aquifers are significantly greater than that of the pumped aquifer. In addition, present methods of analyzing leaky aquifersusually rely only on drawdown data from the pumped aquifer. It is concluded that current methods of analyzing field data from leaky aquifers need to be improved.
• NEUMAN, S. P., & WITHERSPOON, P. A. (1969). THEORY OF FLOW IN A CONFINED TWO AQUIFER SYSTEM. Water Resources Research, 5(4), 803-816.
  More info

  Abstract: Complete analytic solution has been developed for the problem of flow to a well in a confined infinite radial system composed of two aquifers that are separated by an aquitard. The well, which is represented by a line sink, completely penetrates only one of the aquifers and discharges at constant rate. The present solution differs from previous work in that the effects of storage in the aquitard and drawdown in the unpumped aquifer are both considered. Asymptotic solutions for small values of time are also included. The analytic solutions have been evaluated numerically for selected cases. The validity of these solutions has been verified independently using the finite element method of analysis.