Diagrammatic theory of effective hydraulic conductivity

This work presents a stochastic diagrammatic theory for the calculation of the effective hydraulic conductivity of heterogeneous media. The theory is based on the mean-flux series expansion of a log-normal hydraulic conductivity medium in terms of diagrammatic representations and leads to certain general results for the effective hydraulic conductivity of three-dimensional media. A selective summation technique is used to improve low-order perturbation analysis by evaluating an infinite set of diagrammatic terms with a specific topological structure that dominates the perturbation series. For stochastically isotropic media the selective summation yeilds the anticipated exponential expression for the effective hydraulic conductivity. This expression is extended to stochastically anisotropic media. It is also shown that in the case of non homogeneous media the uniform effective hydraulic conductivity is replaced by a non-local tensor kernel, for which general diagrammatic expressions are obtained. The non-local kernel leads to the standard exponential behavior for the effective hydraulic conductivity at the homogeneous limit.     

Stochastic perturbation analysis of groundwater flow. Spatially variable soils,semi-infinite domains and large fluctuations

As is well known, a complete stochastic solution of the stochastic differential equation governing saturated groundwater flow leads to an infinite hierarchy of equations in terms of higher-order moments. Perturbation techniques are commonly used to close this hierarchy, using power-series expansions. These methods are applied by truncating the series after a finite number of terms, and products of random gradients of conductivity and head potential are neglected. Uncertainty regarding the number or terms required to yield a sufficiently accurate result is a significant drawback with the application of power series-based perturbation methods for such problems. Low-order series truncation may be incapable of representing fundamental characteristics of flow and can lead to physically unreasonable and inaccurate solutions of the stochastic flow equation. To support this argument, one-dimensional, steady-state, saturated groundwater flow is examined, for the case of a spatially distributed hydraulic conductivity field. An ordinary power-series perturbation method is used to approximate the mean head, using second-order statistics to characterize the conductivity field. Then an interactive perturbation approach is introduced, which yields improved results compared to low-order, power-series perturbation methods for situations where strong interactions exist between terms in such approximations. The interactive perturbation concept is further developed using Feynman-type diagrams and graph theory, which reduce the original stochastic flow problem to a closed set of equations for the mean and the covariance functions. Both theoretical and practical advantages of diagrammatic solutions are discussed; these include the study of bounded domains and large fluctuations.     

Renormalization group analysis of permeability upscaling

The heterogeneity of the subsurface permeability is considered as the most influential factor in determining groundwater flow and the transport of toxic contaminants. Numerical simulators cannot handle the large grids required to represent the small-scale variability of permeability, and thus explicit estimates of the large-scale behavior in terms of coarse-grained parameters are often required. Perturbation formulations of the effective permeability are based on simplifying assumptions that are valid only for certain probability distributions and weak heterogeneity. A generalized perturbation ansatz that involves higher orders has been proposed (Gelhar and Axness, 1983), but to our knowledge its validity has not been rigorously proved before in three dimensions. In this work we propose a general upscaling formulation valid for strong heterogeneity, general permeability distributions, and media with impermeable zones. We show that the effective permeability is determined by the self-energy series of the permeability fluctuations at zero frequency. Using the diagrammatic representation, we obtain a Dyson equation that involves only irreducible diagrams of the proper self-energy series. We develop a renormalization group (RG) analysis for isotropic lognormal media that proves the generalized perturbation ansatz to all orders. We show that the RG result accurately estimates laboratory permeability measurements in limestone (strong heterogeneity) and sandstone (weak heterogeneity). We also propose an explicit RG estimate for the preasymptotic effective permeability. We compare our results with an approach based on a leading order Green's function expansion (Paleologos et?al., 1996), which, however, requires intensive numerical computations. Finally, we investigate the relation between the RG expression and the algebraic means used in numerical upscaling.     

Use of the Fourier-Laplace transform and of diagrammatical methods to interpret pumping tests in heterogeneous reservoirs

Advances in computer power and in reservoir characterization allow simulation of pressure transients in complex reservoirs generated stochastically. Generally, interpretation of these transient gives useful information about the reservoir hydraulic properties: a major goal is to interpret these transients in the stochastic context. First we ensemble average the pressure over all the random permeability field realizations to derive an equation which drives the ensemble averaged pressure. We use the Fourier transform in space and the Laplace transform in time, in conjuction with a perturbation series expansion in successive powers of the permeability fluctuations to obtain an explicit solution. The Nth order term of this series involves the hydrodynamic interaction between N permeability heterogeneities and after averaging we obtain an expansion containing correlation functions of permeability fluctuations of increasing order.Next, Feynman graphs are introduced allowing a more attractive graphical interpretation of the perturbation series. Then series summation techniques are employed to reduce the graph number to be summed at each order of the fluctuation expansion. This in turn gives useful physical insights on the homogenization processes involved. In particular, it is shown that the sum of the so-called ‘one-particle irreducible graphs’ gives the kernel of a linear integro-differential equation obeyed by the ensemble average pressure. All the information about the heterogeneity structure is contained in this renormalized kernel, which is a limited range function.This equation on its own is the starting point of useful asymptotic results and approximations. In particular it is shown that interpretation of pumping tests yields the steady-state equivalent permeability after a sufficiently long time for an infinite reservoir, as expected.     

Renormalization group methods in subsurface hydrology: overview and applications in hydraulic conductivity upscaling

The renormalization group (RG) approach is a powerful theoretical framework, more suitable for upscaling strong heterogeneity than low-order perturbation expansions. Applications of RG methods in subsurface hydrology include the calculation of (1) macroscopic transport parameters such as effective and equivalent hydraulic conductivity and dispersion coefficients, and (2) anomalous exponents characterizing the dispersion of contaminants due to long-range conductivity correlations or broad (heavy-tailed) distributions of the groundwater velocity. First, we review the main ideas of RG methods and their hydrological applications. Then, we focus on the hydraulic conductivity in saturated porous media with isotropic lognormal heterogeneity, and we present an RG calculation based on the replica method. The RG analysis gives rigorous support to the exponential conjecture for the effective hydraulic conductivity [Water Resour. Res. 19 (1) (1983) 161]. Using numerical simulations in two dimensions with a bimodal conductivity distribution, we demonstrate that the exponential expression is not suitable for all types of heterogeneity. We also introduce an RG coarse-grained conductivity and investigate its applications in estimating the conductivity of blocks or flow domains with finite size. Finally, we define the fractional effective dimension, and we show that it justifies fractal exponents in the range 1−2/dα<1 (where d is the actual medium dimension) in the geostatistical power average.     

A scale‐invariant property of the water retention curve in weakly heterogeneous vadose zones

Abrupt changes of hydraulic properties in a vadose zone are modelled within a stochastic framework, which regards the saturated conductivity and parameters related to the distribution of soil pores as stationary, log‐normally distributed, random space functions. As a consequence, flow variables become random fields, and we aim at deriving an effective Richards equation. To obtain the latter, we adopt a perturbation expansion truncated at the first order (weakly heterogeneous media), which leads to the effective hydraulic conductivity and water retention curves. Overall, the effective properties are scale dependent. However, within the proposed framework, we demonstrate that the inflection point of the laboratory scale retention curve is not affected by the heterogeneity of the vadose zone. Finally, to illustrate the quantitative implications of our results, we consider a monitoring experiment at field scale, and we show how our approach leads to an effective water retention curve, which differs significantly from that which would be obtained without accounting for the above scale‐invariance property.     

Spatial derivatives and perturbation derivatives of amplitude in isotropic and anisotropic media

Explicit equations for the spatial derivatives and perturbation derivatives of amplitude in both isotropic and anisotropic media are derived. The spatial and perturbation derivatives of the logarithm of amplitude can be calculated by numerical quadratures along the rays. The spatial derivatives of amplitude may be useful in calculating the higher-order terms in the ray series, in calculating the higher-order amplitude coefficients of Gaussian beams, in estimating the accuracy of zero-order approximations of both the ray method and Gaussian beams, in estimating the accuracy of the paraxial approximation of individual Gaussian beams, or in estimating the accuracy of the asymptotic summation of paraxial Gaussian beams. The perturbation derivatives of amplitude may be useful in perturbation expansions from elastic to viscoelastic media and in estimating the accuracy of the common-ray approximations of the amplitude in the coupling ray theory.     

Seismic wave scattering inversion for fluid factor of heterogeneous media

Elastic wave inverse scattering theory plays an important role in parameters estimation of heterogeneous media. Combining inverse scattering theory, perturbation theory and stationary phase approximation, we derive the P-wave seismic scattering coefficient equation in terms of fluid factor, shear modulus and density of background homogeneous media and perturbation media. With this equation as forward solver, a pre-stack seismic Bayesian inversion method is proposed to estimate the fluid factor of heterogeneous media. In this method, Cauchy distribution is utilized to the ratios of fluid factors, shear moduli and densities of perturbation media and background homogeneous media, respectively. Gaussian distribution is utilized to the likelihood function. The introduction of constraints from initial smooth models enhances the stability of the estimation of model parameters. Model test and real data example demonstrate that the proposed method is able to estimate the fluid factor of heterogeneous media from pre-stack seismic data directly and reasonably.     

Electrokinetic Migration of Nitrate Through Heterogeneous Granular Porous Media

This study investigates and quantifies the influence of physical heterogeneity in granular porous media, represented by materials with different hydraulic conductivity, on the migration of nitrate, used as an amendment to enhance bioremediation, under an electric field. Laboratory experiments were conducted in a bench‐scale test cell under a low applied direct current using glass bead and clay mixes and synthetic groundwater to represent ideal conditions. The experiments included bromide tracer tests in homogeneous settings to deduce controls on electrokinetic transport of inorganic solutes in the different materials, and comparison of nitrate migration under homogeneous and heterogeneous scenarios. The results indicate that physical heterogeneity of subsurface materials, represented by a contrast between a higher‐hydraulic conductivity and lower‐hydraulic conductivity material normal to the direction of the applied electric field exerts the following controls on nitrate migration: (1) a spatial change in nitrate migration rate due to changes in effective ionic mobility and subsequent accumulation of nitrate at the interface between these materials; and (2) a spatial change in the voltage gradient distribution across the hydraulic conductivity contrast, due to the inverse relationship with effective ionic mobility. These factors will contribute to higher mass transport of nitrate through low hydraulic conductivity zones in heterogeneous porous media, relative to homogeneous host materials. Overall electrokinetic migration of amendments such as nitrate can be increased in heterogeneous granular porous media to enhance the in situ bioremediation of organic contaminants present in low hydraulic conductivity zones.     

A geostatistical approach to an inverse problem: Identification of geometry and estimate of equivalent conductivities for highly heterogeneous porous media with the differential system method

A methodology for identifying the geometry of different materials in highly heterogeneous porous media in discrete inverse problems (DIP) is described. It applies a geostatistical approach within the differential system method (DSM). DSM calculates conductivity values along an integration path beginning at a point with known conductivity. In aquifers with zero source terms, DSM completely describes the conductivity field through a spatially distributed parameter depending on hydraulic head gradients and integration path. A factor analysis of the structural components of this parameter (i.e. coregionalisation analysis) was carried out to identify the geometry of different materials, corresponding to distinct statistically homogeneous areas. The equivalent conductivity values for homogeneous areas were estimated.This approach was applied for a synthetic aquifer. The identification of geometry was accurate and the estimates of equivalent parameters were good, compared with reference values. The accuracy of the results depended on errors in hydraulic gradients, compared with conductivity gradients.     

Binary upscaling—the role of connectivity and a new formula

Although recognized as important, measures of connectivity (i.e. the existence of high-conductivity paths that increase flow and allow for early solute arrival) have not yet been incorporated into methods for upscaling hydraulic conductivities of porous media. We present and evaluate a binary upscaling formula that utilizes connectivity information. The upscaled hydraulic conductivity (K) of binary media is determined as a function of the proportions and conductivities of the two materials, the geometry of the inclusions, and the mean distance between them. The use of a phase interchange theorem renders the formula equally applicable to two-dimensional media with inclusions of low K and high K as compared with the matrix. The new upscaling formula is tested on two-dimensional binary random fields spanning a broad range of spatial correlation structures and conductivity contrasts. The computed effective conductivities are compared to what is obtained using self-consistent effective medium theory, the coated ellipsoids approximation, and to a streamline approach. It is shown that, although simple, the proposed formula performs better than available methods for binary upscaling. The use of connectivity information leads to significantly improved behavior close to the percolation threshold. The proposed upscaling formula depends exclusively on parameters that are obtainable from field investigations.     

Nonlinear ray perturbation theory with its applications to ray tracing and inversion in anisotropic media

A comprehensive approach, based on the general nonlinear ray perturbation theory (Druzhinin, 1991), is proposed for both a fast and accurate uniform asymptotic solution of forward and inverse kinematic problems in anisotropic media. It has been developed to modify the standard ray linearization procedures when they become inconsistent, by providing a predictable truncation error of ray perturbation series. The theoretical background consists in a set of recurrent expressions for the perturbations of all orders for calculating approximately the body wave phase and group velocities, polarization, travel times, ray trajectories, paraxial rays and also the slowness vectors or reflected/transmitted waves in terms of elastic tensor perturbations. We assume that any elastic medium can be used as an unperturbed medium. A total 2-D numerical testing of these expressions has been established within the transverse isotropy to verify the accuracy and convergence of perturbation series when the elastic constants are perturbed. Seismological applications to determine crack-induced anisotropy parameters on VSP travel times for the different wave types in homogeneous and horizontally layered, transversally isotropic and orthorhombic structures are also presented. A number of numerical tests shows that this method is in general stable with respect to the choice of the reference model and the errors in the input data. A proof of uniqueness is provided by an interactive analysis of the sensitivity functions, which are also used for choosing optimum source/receiver locations. Finally, software has been developed for a desktop computer and applied to interpreting specific real VSP observations as well as explaining the results of physical modelling for a 3-D crack model with the estimation of crack parameters.     

Stochastic analysis of spatial variability in unconfined groundwater flow

In this paper, spatial variability in steady one-dimensional unconfined groundwater flow in heterogeneous formations is investigated. An approach to deriving the variance of the hydraulic head is developed using the nonlinear filter theory. The nonlinear governing equation describing the one-dimensional unconfined groundwater flow is decomposed into three linear partial differential equations using the perturbation method. The linear and quadratic frequency response functions are obtained from the first- and second-order perturbation equations using the spectral method. Furthermore, under the assumption of the exponential covariance function of log hydraulic conductivity, the analytical solutions of both the spectrum and the variance of the hydraulic head produced from the linear system are derived. The results show that the variance derived herein is less than that of Gelhar (1977). The reason is that the log transmissivity is linearized in Gelhars work. In addition, the analytical solutions of both the spectrum and the variance of the hydraulic head produced from the quadratic system are derived as well. It is found that the correlation scale and the trend in mean of log hydraulic conductivity are important to the dimensionless variance ratio.     

Truncated multiGaussian fields and effective conductance of binary media

Truncated Gaussian fields provide a flexible model for defining binary media with dispersed (as opposed to layered) inclusions. General properties of excursion sets on these truncated fields are coupled with a distance-based upscaling algorithm and approximations of point process theory to develop an estimation approach for effective conductivity in two-dimensions. Estimation of effective conductivity is derived directly from knowledge of the kernel size used to create the multiGaussian field, defined as the full-width at half maximum (FWHM), the truncation threshold and conductance values of the two modes. Therefore, instantiation of the multiGaussian field is not necessary for estimation of the effective conductance. The critical component of the effective medium approximation developed here is the mean distance between high conductivity inclusions. This mean distance is characterized as a function of the FWHM, the truncation threshold and the ratio of the two modal conductivities. Sensitivity of the resulting effective conductivity to this mean distance is examined for two levels of contrast in the modal conductances and different FWHM sizes. Results demonstrate that the FWHM is a robust measure of mean travel distance in the background medium. The resulting effective conductivities are accurate when compared to numerical results and results obtained from effective media theory, distance-based upscaling and numerical simulation.     

Spatial weighting functions: transient hydraulic tests and heterogeneous media

To improve understanding of property measurements in heterogeneous media, an energy-based weighting function concept is developed. In (assumed) homogeneous media, the instrument spatial weighting function (ISWF) depends only on the energy dissipation distribution set up by the measurement procedure and it reduces to simply inverse sample volume (uniform weighting) for 1-D parallel flow case (ideal permeameter). For 1-D transient flow in homogeneous media, such as with slug tests, the ISWF varies with position and time, with 95% of the total weighting contained within 115 well radii, even late in the test. In the heterogeneous case, the determination of the ISWF is connected to the problem of determining an equivalent hydraulic conductivity (K), where the criterion for equivalence is based on equal energy dissipation rate rather than equal volume discharge. The discharge-based equivalent K (K(E)) and the energy-based equivalent K in heterogeneous media (K(eh)) are not equal in general, with K(eh) typically above the nodal arithmetic mean K. The possibly more fundamental problem is that as one makes K measurements in heterogeneous media at different locations or on different cores of heterogeneous materials, the ISWF will be heterogeneity dependent, implying that the averaging process resulting in the equivalent K value also varies with position. If the testing procedure is transient, then the averaging process varies with time. This suggests a fundamental ambiguity in the interpretation of hydraulic conductivity measurements in heterogeneous media that may impact how we approach modeling and prediction in a practical sense (Molz 2003). Further research is suggested.     

Perturbation‐based moveout approximations in anisotropic media

The moveout approximations play an important role in seismic data processing. The standard hyperbolic moveout approximation is based on an elliptical background model with two velocities: vertical and normal moveout. We propose a new set of moveout approximations based on a perturbation series in terms of anellipticity parameters using the alternative elliptical background model defined by vertical and horizontal velocities. We start with a transversely isotropic medium with a vertical symmetry axis. Then, we extend this approach to a homogeneous orthorhombic medium. To define the perturbation coefficients for a new background, we solve the eikonal equation with horizontal velocities in transversely isotropic medium with a vertical symmetry axis and orthorhombic media. To stabilise the perturbation series and improve the accuracy, the Shanks transform is applied for all the cases. We select different parameterisations for both velocities and anellipticity parameters for an orthorhombic model. From the comparison in traveltime error, the new moveout approximations result in better accuracy comparing with the standard perturbation‐based methods and other approximations.     

Influence of heterogeneity on unsaturated hydraulic properties: (1) local heterogeneity and scale effect

Spatial heterogeneity is ubiquitous in nature, which may significantly affect the soil hydraulic property curves. The models of a closed‐form functional relationship of soil hydraulic property curves (e.g. VG model or exponential model) are valid at point or local scale based on a point‐scale hydrological process, but how do scale effects of heterogeneity have an influence on the parameters of these models when the models are used in a larger scale process? This paper uses a two‐dimensional variably saturated flow and solute transport finite element model (VSAFT2) to simulate variations of pressure and moisture content in the soil flume under a constant head boundary condition. By changing different numerical simulation block sizes, a quantitative evaluation of parameter variations in the VG model, resulting from the scale effects, is presented. Results show that the parameters of soil hydraulic properties are independent of scale in homogeneous media. Parameters of α and n in homogeneous media, which are estimated by using the unsaturated hydraulic conductivity curve (UHC) or the soil water retention curve (WRC), are identical. Variations of local heterogeneities strongly affect the soil hydraulic properties, and the scale affects the results of the parameter estimations when numerical experiments are conducted. Furthermore, the discrepancy of each curve becomes considerable when moisture content becomes closer to a dry situation. Parameters estimated by UHC are totally different from the ones estimated by WRC. Copyright © 2012 John Wiley & Sons, Ltd.     

Interpretation of spring recession curves

Recession curves contain information on storage properties and different types of media such as porous, fractured, cracked lithologies and karst. Recession curve analysis provides a function that quantitatively describes the temporal discharge decay and expresses the drained volume between specific time limits (Hall 1968). This analysis also allows estimating the hydrological significance of the discharge function parameters and the hydrological properties of the aquifer. In this study, we analyze data from perennial springs in the Judean Mountains and from others in the Galilee Mountains, northern Israel. All the springs drain perched carbonate aquifers. Eight of the studied springs discharge from a karst dolomite sequence, whereas one flows out from a fractured, slumped block of chalk. We show that all the recession curves can be well fitted by a function that consists of two exponential terms with exponential coefficients alpha1 and alpha2. These coefficients are approximately constant for each spring, reflecting the hydraulic conductivity of different media through which the ground water flows to the spring. The highest coefficient represents the fast flow, probably through cracks, or quickflow, whereas the lower one reflects the slow flow through the porous medium, or baseflow. The comparison of recession curves from different springs and different years leads to the conclusion that the main factors that affect the recession curve exponential coefficients are the aquifer lithology and the geometry of the water conduits therein. In normal years of rainy winter and dry summer, alpha1 is constant in time. However, when the dry period is longer than usual because of a dry winter, alpha1 slightly decreases with time.     

Generalization of TOPMODEL for a power law transmissivity profile

A generalization of the TOPMODEL equations for a power law vertical profile of hydraulic conductivity is introduced. The exponential profile of TOPMODEL is obtained as a limit case of the new general form. © 1997 John Wiley & Sons, Ltd.     

FracKfinder: A MATLAB Toolbox for Computing Three-Dimensional Hydraulic Conductivity Tensors for Fractured Porous Media

Fractures in porous media have been documented extensively. However, they are often omitted from groundwater flow and mass transport models due to a lack of data on fracture hydraulic properties and the computational burden of simulating fractures explicitly in large model domains. We present a MATLAB toolbox, FracKfinder, that automates HydroGeoSphere (HGS), a variably saturated, control volume finite-element model, to simulate an ensemble of discrete fracture network (DFN) flow experiments on a single cubic model mesh containing a stochastically generated fracture network. Because DFN simulations in HGS can simulate flow in both a porous media and a fracture domain, this toolbox computes tensors for both the matrix and fractures of a porous medium. Each model in the ensemble represents a different orientation of the hydraulic gradient, thus minimizing the likelihood that a single hydraulic gradient orientation will dominate the tensor computation. Linear regression on matrices containing the computed three-dimensional hydraulic conductivity (K) values from each rotation of the hydraulic gradient is used to compute the K tensors. This approach shows that the hydraulic behavior of fracture networks can be simulated where fracture hydraulic data are limited. Simulation of a bromide tracer experiment using K tensors computed with FracKfinder in HGS demonstrates good agreement with a previous large-column, laboratory study. The toolbox provides a potential pathway to upscale groundwater flow and mass transport processes in fractured media to larger scales.