Richard L. Naff

Nonreactive and reactive solute transport in three‐dimensional heterogeneous porous media: Mean displacement, plume spreading, and uncertainty

The field-scale transport of reactive and nonreactive solutes by groundwater in a statistically anisotropic aquifer is examined by means of high-resolution, three-dimensional numerical solutions of the steady state flow and transient advection-dispersion equations. The presence of physical and chemical heterogeneities in the aquifer media is modeled with the use of a geostatistical description of the hydraulic conductivity and chemical distribution coefficient. The geostatistical parameters describing the spatial variations of the log-transformed hydraulic conductivity fields, ln [K(x)], are chosen to resemble those of the Borden aquifer. For sorbing solutes the spatial variation in the log-transformed distribution coefficient fields, ln [Kd(x)], is generated so as to have a geometric mean and variance similar to that estimated by Durant (1986) for the organic chemical tetrachlorethane with Borden sand. It is further assumed that the ln [K(x)] and ln [Kd(x)] fields are inversely correlated to each other and that they possess the same spatial correlation structure. Five realizations of media which incorporate these characteristics are generated by means of the Fourier spectral technique of Robin et al. (1993). It is shown that joint K(x) and Kd(x) variability can impart a large-scale “pseudokinetic” behavior, in that the ensemble mean bulk retardation factor can increase with time and plume displacement distance, even though sorption is modeled as being linear and instantaneous at the scale of any single heterogeneity and the flow field is assumed to be steady state. From realization to realization, however, the apparent bulk retardation factor can either increase dramatically at early time in a manner similar to that observed by Roberts et al. (1986) during the Borden tracer test or it can decrease. At large time the ensemble mean velocity of the centroid of the reactive plume is close to a value given by the mean fluid velocity divided by the arithmetic mean of the locally variable retardation factors. Local-scale transport nonidealities, such as intraparticle diffusion and relatively rapid kinetic sorption, are shown to have minimal influence on the plume centroid velocity and macrodispersivity, relative to the effect of aquifer heterogeneity. The results of the numerical simulations further demonstrate that the first-order stochastic analyses of Gelhar and Axness (1983), Dagan (1988), and Naff (1990) tend to overestimate the actual field-scale longitudinal spreading of a nonreactive solute. This result is believed to occur because these analyses include the artificial effect of plume centroid dispersion about the ensemble mean position, in addition to the actual spreading of each plume realization about its center of mass. When the effect of plume centroid dispersion is added to the numerical simulation results, a reasonable agreement with the solution of Naff (1990), which takes into consideration local-scale dispersion, is obtained for the longitudinal macrodispersivity. In the transverse direction, the ensemble mean spreading agrees reasonably well with the solution of Dagan (1988). It is also shown that the longitudinal macrodispersivity of a reactive solute can be enhanced relative to that of a nonreactive one. Its actual value, however, is overpredicted by the first-order solution of Garabedian (1987) for the reasons given above. Also explored, in a preliminary fashion, are some issues related to prediction uncertainty for both reactive and nonreactive solutes migrating through heterogeneous aquifers.

Radial flow in heterogeneous porous media: An analysis of specific discharge

A perturbation solution to three-dimensional radial flow in heterogeneous porous media is reported on; the analytical solution is in terms of the first and second moments of the gradient in the radial direction and the specific discharge in this direction. The solution requires evaluation, numerically in the case of this paper, of some rather difficult integrals which compose its basic form. By proper adjustment of the statistical anisotropy of the medium, a solution for the two-dimensional flow problem can be emulated; the two-dimensional solution for the mean and variance in specific discharge resembles asymptotes reported earlier by another investigator. In general, the effective hydraulic conductivity, outside a few length scales of the well itself, is found to be nearly constant and dependent upon the value selected for the statistical anisotropy factor; the variance in specific discharge behaves nearly as the inverse square of the distance from the well in all cases.

HIGH-RESOLUTION MONTE CARLO SIMULATION OF FLOW AND CONSERVATIVE TRANSPORT IN HETEROGENEOUS POROUS MEDIA 1. METHODOLOGY AND FLOW RESULTS

In this, the first of two papers concerned with the use of numerical simulation to examine flow and transport parameters in heterogeneous porous media via Monte Carlo methods, various aspects of the modelling effort are examined. In particular, the need to save on core memory causes one to use only specific realizations that have certain initial characteristics; in effect, these transport simulations are conditioned by these characteristics. Also, the need to independently estimate length scales for the generated fields is discussed. The statistical uniformity of the flow field is investigated by plotting the variance of the seepage velocity for vector components in the x, y, and z directions. Finally, specific features of the velocity field itself are illuminated in this first paper. In particular, these data give one the opportunity to investigate the effective hydraulic conductivity in a flow field which is approximately statistically uniform; comparisons are made with first- and second-order perturbation analyses. The mean cloud velocity is examined to ascertain whether it is identical to the mean seepage velocity of the model. Finally, the variance in the cloud centroid velocity is examined for the effect of source size and differing strengths of local transverse dispersion.

High‐resolution Monte Carlo simulation of flow and conservative transport in heterogeneous porous media: 2. Transport results

In this, the second of two papers concerned with the use of numerical simulation to examine flow and transport parameters in heterogeneous porous media via Monte Carlo methods, results from the transport aspect of these simulations are reported on. Transport simulations contained herein assume a finite pulse input of conservative tracer, and the numerical technique endeavors to realistically simulate tracer spreading as the cloud moves through a heterogeneous medium. Medium heterogeneity is limited to the hydraulic conductivity field, and generation of this field assumes that the hydraulic-conductivity process is second-order stationary. Methods of estimating cloud moments, and the interpretation of these moments, are discussed. Techniques for estimation of large-time macrodispersivities from cloud second-moment data, and for the approximation of the standard errors associated with these macrodispersivities, are also presented. These moment and macrodispersivity estimation techniques were applied to tracer clouds resulting from transport scenarios generated by specific Monte Carlo simulations. Where feasible, moments and macrodispersivities resulting from the Monte Carlo simulations are compared with first- and second-order perturbation analyses. Some limited results concerning the possible ergodic nature of these simulations, and the presence of non-Gaussian behavior of the mean cloud, are reported on as well.

Shape Functions for Velocity Interpolation in General Hexahedral Cells

Numerical methods for grids with irregular cells require discrete shape functions to approximate the distribution of quantities across cells. For control-volume mixed finite-element (CVMFE) methods, vector shape functions approximate velocities and vector test functions enforce a discrete form of Darcys law. In this paper, a new vector shape function is developed for use with irregular, hexahedral cells (trilinear images of cubes). It interpolates velocities and fluxes quadratically, because as shown here, the usual Piola-transformed shape functions, which interpolate linearly, cannot match uniform flow on general hexahedral cells. Truncation-error estimates for the shape function are demonstrated. CVMFE simulations of uniform and non-uniform flow with irregular meshes show first- and second-order convergence of fluxes in the L2 norm in the presence and absence of singularities, respectively.

Test functions for three-dimensionalcontrol-volume mixed finite element methods on irregular grids

Numerical methods based on unstructured grids, with irregular cells, usually require discrete shape functions to approximate the distribution of quantities across cells. For control-volume mixed finite element methods, vector shape functions are used to approximate the distribution of velocities across cells and vector test functions are used to minimize the error associated with the numerical approximation scheme. For a logically cubic mesh, the lowest-order shape functions are chosen in a natural way to conserve intercell fluxes that vary linearly in logical space. Vector test functions, while somewhat restricted by the mapping into the logical reference cube, appear to accept a wider class of possibilities. Ideally, an error minimization procedure to select the test function from an acceptable class of candidates would be the best procedure. Lacking such a procedure, we first investigate the effect of possible test functions on the pressure distribution over the control volume; we look for test functions which allow for the elimination of intermediate pressures on cell faces. From these results, we select three forms for the test function for use in a control-volume mixed method code and subject them to an error analysis for different forms of grid irregularity; errors are reported in terms of the discrete L2 norm of the velocity error. Of these three forms, one appears to produce optimal results for most forms of grid irregularity.

An Eulerian scheme for the second-order approximation of subsurface transport moments

The moments of a conservative tracer cloud migrating in a mean uniform flow field are estimated using an operator approximation scheme; results are presented for the second, third, and fourth central moments in the mean flow direction. It is assumed that the spatially variable flow field, and therefore the tracer migration problem itself, is amenable to a probabilistic description; the effects of local dispersion on cloud migration are neglected in this study. Variation in the flow field is assumed to be the result of spatial variation in the hydraulic conductivity; spatial variation in porosity is assumed negligible. The operator approximation scheme, as implemented in this study, is second-order correct, which requires a second-order correct approximation of the velocity field correlation structure. Because estimation of the velocity correlation structure is decidedly the most difficult aspect of second-order analysis, an ad hoc extension of the imperfectly stratified approximation developed earlier is implemented for this purpose. The first-order approximation resulting from the operator expansion scheme is equivalent to small perturbation Eulerian results presented earlier (Naff, 1990, 1992). The infinite-order approximation resulting from this scheme is equivalent to the exponential operator results obtained by Van Kampen (1976).

A note on conservative transport in anisotropic, heterogeneous porous media in the presence of small‐amplitude transients

The late-time macrodispersion coefficients are obtained for the case of flow in the presence of a small-scale deterministic transient in a three-dimensional anisotropic, heterogeneous medium. The transient is assumed to affect only the velocity component transverse to the mean flow direction and to take the form of a periodic function. For the case of a highly stratified medium, these late-time macrodispersion coefficients behave largely as the standard coefficients used in the transport equation. Only in the event that the medium is isotropic is it probable that significant deviations from the standard coefficients would occur.

Efficient solver for mixed finite element method and control-volume mixed finite element method in 3-D on Hexahedral grids

An efficient solver for the mixed finite element (MFE) method on 3-D distorted hexahedral cells has been implemented. The method relies on being able to define a divergence-free subspace. The divergence-free subspace is used to eliminate the pressure variable from the set of equations, resulting in a symmetric positive definite system. The pre-conditioned conjugate gradient (PCG) method is used to solve the reduced system. A two-level additive Schwarz method is implemented as the preconditioner. We show that the same efficient solver can be used for the control-volume mixed finite element (CVMFE) method. The CVMFE method uses a discontinous vector test space which leads to more accurate velocities in the presence of discontinous anisotropic hydraulic conductivities on distorted hexahedral grids.

A Comparison of Preconditioning Techniques for Parallelized PCG Solvers for the CCFD Problem

Parallel algorithms for solving sparse symmetric matrix systems that might result from the cell-centered finite difference (CCFD) scheme are compared. Parallelization is based in partitioning the mass matrix such that each partition is controlled by a separate process. These processes may then be distributed among a networked cluster of processors using the standard Message-Passing Interface (MPI). MPI software allows for multiple, simultaneous processes to coordinate and exchange information for some central purpose. In this study, partitioning of the mass matrix is based in decomposing the CCFD domain into non-overlapping subdomains; each subdomain corresponds to a partition of the mass matrix. The subdomain partitions are numbered alternately, using a red/black numbering scheme. Partitions are linked by the CCFD coefficients corresponding to cell edges that coincide with subdomain boundaries internal to the domain. A major portion of this work examines the best way to handle connectivity information between partitions. Algorithms considered in this study are based in the preconditioned conjugate gradient scheme (PCG) and differ only in the preconditioning used. Parallelization of the PCG solver entails running essentially identical conjugate-gradient loops on separate processes for every subdomain partition. These loops must exchange information globally to calculate inner products and locally for sharing connectivity information between partitions. The classic incomplete Cholesky preconditioner with zero fill (IC(0)) is used as the standard for comparison. Another preconditioning scheme considered is based principally in an approximate block Gaussian (BG) iterative solution to the problem. In both these preconditioners, arrays containing connectivity information are passed between processes corresponding to adjacent partitions. Because of this need to pass connectivity information, the incomplete Cholesky preconditioner is limited, in practice, to a zero fill application. This limitation can be partially alleviated by using a BG iteration as a preconditioner, which requires approximate solves of individual matrix partitions. These approximate solves can be carried out using any number of preconditioners, including IC(0); however, BG iteration involves an additional level of approximation, giving the result that BG iteration with approximate IC(0) solves is less efficient than using simple IC(0) as the preconditioner in the parallel PCG algorithm. Preconditioners based in the Jacobi scheme are also considered as they require no connectivity information. These preconditioners represent a trade off between lower communication costs at the expense of increased work to obtain convergence. We are presently exploring variants of these methods for use as preconditioners in the parallelized conjugate gradient scheme; we expect to report on the results of this work.