Block-effective macrodispersion for numerical simulations of sorbing solute transport in heterogeneous porous formations

Heterogeneity is prevalent in aquifers and has an enormous impact on contaminant transport in groundwater. Numerical simulations are an effective way to deal with heterogeneity directly by assigning different hydraulic property values to each numerical grid block. Because hydraulic properties vary on different scales, but they cannot be sampled exhaustively and the number of numerical grid blocks is limited by computational considerations, the dispersive effects of unmodeled heterogeneity need to be accounted for. Dispersion tensors can be used to model the dispersion caused by unmodeled heterogeneity. The concept of block-effective macrodispersion tensors for modeling the effects of small-scale variability on solute transport introduced by Rubin et al. [Rubin Y, Sun A, Maxwell R, Bellin A. The concept of block-effective macrodispersivity and a unified approach for grid-scale- and plume-scale-dependent transport. J Fluid Mech 1999;395:161–80] is extended in this paper for use with reactive solutes. The tensors are derived for reactive solutes with spatially variable retardation factors and for solutes experiencing spatially uniform rate-limited sorption. The longitudinal block-effective macrodispersion coefficient is largest for perfect negative correlation between the log-hydraulic conductivity and the retardation factor. Because dispersion tensors, as they are usually implemented in numerical simulations, produce symmetric spreading, the applicability of the concept depends on the portion of the plume asymmetry caused by small-scale variability. The presented results show that the concept is applicable for rate-limited sorption for block sizes of one and two integral scales.     

A meshless method to simulate solute transport in heterogeneous porous media

We derive a meshless numerical method based on smoothed particle hydrodynamics (SPH) for the simulation of conservative solute transport in heterogeneous geological formations. We demonstrate that the new proposed scheme is stable, accurate, and conserves global mass. We evaluate the performance of the proposed method versus other popular numerical methods for the simulation of one- and two-dimensional dispersion and two-dimensional advective–dispersive solute transport in heterogeneous porous media under different Pèclet numbers. The results of those benchmarks demonstrate that the proposed scheme has important advantages over other standard methods because of its natural ability to control numerical dispersion and other numerical artifacts. More importantly, while the numerical dispersion affecting traditional numerical methods creates artificial mixing and dilution, the new scheme provides numerical solutions that are “physically correct”, greatly reducing these artifacts.     

Effect of permeability variations on solute transport in highly heterogeneous porous media

The effect of aquifer heterogeneity on flow and solute transport in two-dimensional isotropic porous media was analyzed using the Monte Carlo method. The two-dimensional logarithmic permeability (ln K) was assumed to be a non-stationary random field with its increments being a truncated fractional Lévy motion (fLm). The permeability fields were generated using the modified successive random additions (SRA) algorithm code SRA3DC [http://www.iamg.org/CGEditor/index.htm]. The velocity and concentration fields were computed respectively for two-dimensional flow and transport with a pulse input using the finite difference codes of MODFLOW 2000 and MT3DMS. Two fLm control parameters, namely the width parameter (C) and the Lévy index (α), were varied systematically to examine their effect on the resulting permeability, flow velocity and concentration fields. We also computed the first- and second-spatial moments, the dilution index, as well as the breakthrough curves at different control planes with the corresponding concentration fields. In addition, the derived breakthrough curves were fitted using the continuous time random walk (CTRW) and the traditional advection-dispersion equation (ADE). Results indicated that larger C and smaller α both led to more heterogeneous permeability and velocity fields. The Lévy-stable distribution of increments in ln K resulted in a Lévy-stable distribution of increments in logarithm of the velocity (ln v). Both larger C and smaller α created sharper leading edges and wider tailing edges of solute plumes. Furthermore, a relatively larger amount of solute still remained in the domain after a relatively longer time transport for smaller α values. The dilution indices were smaller than unity and increased as C increased and α decreased. The solute plume and its second-spatial moments increased as C increased and α decreased, while the first-spatial moments of the solute plume were independent of C and α values. The longitudinal macrodispersivity was scale-dependent and increased as a power law function of time. Increasing C and decreasing α both resulted in an increase in longitudinal macrodispersivity. The transport in such highly heterogeneous media was slightly non-Gaussian with its derived breakthrough curves being slightly better fitted by the CTRW than the ADE, especially in the early arrivals and late-time tails.     

Numerical simulations of transport of non-ergodic solute plumes in heterogeneous aquifers

Transport of non-ergodic solute plumes by steady-state groundwater flow with a uniform mean velocity, μ, were simulated with Monte Carlo approach in a two-dimensional heterogeneous and statistically isotropic aquifer whose transmissivity, T, is log-normally distributed with an exponential covariance. The ensemble averages of the second spatial moments of the plume about its center of mass, <S _{i i} (t)>, and the plume centroid covariance, R _{i i} (t) (i=1,2), were simulated for the variance of Y=log T, σ _Y ²=0.1, 0.5 and 1.0 and line sources normal or parallel to μ of three dimensionless lengths, 1, 5, and 10. For σ _Y ²=0.1, all simulated <S _{i i} (t)>−S _{i i} (0) and R _{i i} (t) agree well with the first-order theoretical values, where S _{i i} (0) are the initial values of S _{i i} (t). For σ _Y ²=0.5 and 1.0 and the line sources normal to μ, the simulated longitudinal moments, <S ₁₁(t)>−S ₁₁(0) and R ₁₁(t), agree well with the first-order theoretical results but the simulated transverse moments <S ₂₂(t)>−S ₂₂(0) and R ₂₂(t) are significantly larger than the first-order values. For the same two larger values of σ _Y ² but the line sources parallel to μ, the simulated <S ₁₁(t)>−S ₁₁(0) are larger than but the simulated R ₁₁ are smaller than the first-order values, and both simulated <S ₂₂(t)>−S ₂₂(0) and R ₂₂(t) stay larger than the first-order values. For a fixed value of σ _Y ², the summations of <S _{i i} (t)>−S _{i i} (0) and R _{i i} , i.e., X _{i i} (i=1,2), remain almost the same no matter what kind of source simulated. The simulated X ₁₁ are in good agreement with the first-order theory but the simulated X ₂₂ are significantly larger than the first-order values. The simulated X ₂₂, however, are in excellent agreement with a previous modeling result and both of them are very close to the values derived using Corrsin's conjecture. It is found that the transverse moments may be significantly underestimated if less accurate hydraulic head solutions are used and that the decreasing of <S ₂₂(t)>−S ₂₂(0) with time or a negative effective dispersivity, defined as , may happen in the case of a line source parallel to μ where σ _Y ² is small.     

Numerical simulations of transport of non-ergodic solute plumes in heterogeneous aquifers

Transport of non-ergodic solute plumes by steady-state groundwater flow with a uniform mean velocity, μ, were simulated with Monte Carlo approach in a two-dimensional heterogeneous and statistically isotropic aquifer whose transmissivity, T, is log-normally distributed with an exponential covariance. The ensemble averages of the second spatial moments of the plume about its center of mass, <S _{i i} (t)>, and the plume centroid covariance, R _{i i} (t) (i=1,2), were simulated for the variance of Y=log T, σ _Y ²=0.1, 0.5 and 1.0 and line sources normal or parallel to μ of three dimensionless lengths, 1, 5, and 10. For σ _Y ²=0.1, all simulated <S _{i i} (t)>−S _{i i} (0) and R _{i i} (t) agree well with the first-order theoretical values, where S _{i i} (0) are the initial values of S _{i i} (t). For σ _Y ²=0.5 and 1.0 and the line sources normal to μ, the simulated longitudinal moments, <S ₁₁(t)>−S ₁₁(0) and R ₁₁(t), agree well with the first-order theoretical results but the simulated transverse moments <S ₂₂(t)>−S ₂₂(0) and R ₂₂(t) are significantly larger than the first-order values. For the same two larger values of σ _Y ² but the line sources parallel to μ, the simulated <S ₁₁(t)>−S ₁₁(0) are larger than but the simulated R ₁₁ are smaller than the first-order values, and both simulated <S ₂₂(t)>−S ₂₂(0) and R ₂₂(t) stay larger than the first-order values. For a fixed value of σ _Y ², the summations of <S _{i i} (t)>−S _{i i} (0) and R _{i i} , i.e., X _{i i} (i=1,2), remain almost the same no matter what kind of source simulated. The simulated X ₁₁ are in good agreement with the first-order theory but the simulated X ₂₂ are significantly larger than the first-order values. The simulated X ₂₂, however, are in excellent agreement with a previous modeling result and both of them are very close to the values derived using Corrsin's conjecture. It is found that the transverse moments may be significantly underestimated if less accurate hydraulic head solutions are used and that the decreasing of <S ₂₂(t)>−S ₂₂(0) with time or a negative effective dispersivity, defined as , may happen in the case of a line source parallel to μ where σ _Y ² is small.     

Comparison of theory and experiment for solute transport in weakly heterogeneous bimodal porous media

In this work, the influence of non-equilibrium effects on solute transport in a weakly heterogeneous medium is discussed. Three macro-scale models (upscaled via the volume averaging technique) are investigated: (i) the two-equation non-equilibrium model, (ii) the one-equation asymptotic model and (iii) the one-equation local equilibrium model. The relevance of each of these models to the experimental system conditions (duration of the pulse injection, dispersivity values…) is analyzed. The numerical results predicted by these macroscale models are compared directly with the experimental data (breakthrough curves). Our results suggest that the preasymptotic zone (for which a non-Fickian model is required) increases as the solute input pulse time decreases. Beyond this limit, the asymptotic regime is recovered. A comparison with the results issued from the stochastic theory for this regime is performed. Results predicted by both approaches (volume averaging method and stochastic analysis) are found to be consistent.     

An efficient,high-order probabilistic collocation method on sparse grids for three-dimensional flow and solute transport in randomly heterogeneous porous media

In this study, a probabilistic collocation method (PCM) on sparse grids is used to solve stochastic equations describing flow and transport in three-dimensional, saturated, randomly heterogeneous porous media. The Karhunen–Loève decomposition is used to represent log hydraulic conductivity Y=ln_Ks$Y=\ln K_s$. The hydraulic head h   and average pore-velocity v$\mathbf{v}$ are obtained by solving the continuity equation coupled with Darcy’s law with random hydraulic conductivity field. The concentration is computed by solving a stochastic advection–dispersion equation with stochastic average pore-velocity v$\mathbf{v}$ computed from Darcy’s law. The PCM approach is an extension of the generalized polynomial chaos (gPC) that couples gPC with probabilistic collocation. By using sparse grid points in sample space rather than standard grids based on full tensor products, the PCM approach becomes much more efficient when applied to random processes with a large number of random dimensions. Monte Carlo (MC) simulations have also been conducted to verify accuracy of the PCM approach and to demonstrate that the PCM approach is computationally more efficient than MC simulations. The numerical examples demonstrate that the PCM approach on sparse grids can efficiently simulate solute transport in randomly heterogeneous porous media with large variances.     

Models of sub-grid variability in numerical simulations of solute transport in heterogeneous porous formations: three-dimensional flow and effect of pore-scale dispersion

Modeling transport of contaminants in the earths subsurface relies on numerical solutions over grids with blocks larger than Darcys scale. The hydraulic conductivity is homogenized over the grid blocks and the plumes spreading is reduced as a consequence of the wiped-out variability. To compensate for this loss Rubin et al. (1999) proposed to augment mixing by block-effective dispersion coefficients, and Rubin et al. (2003) showed, by means of two dimensional simulations, how this concept can be applied in practice. In this paper, we present new solutions of the block-effective dispersion tensor for an axisymmetric exponential covariance model. In addition, we discuss the influence of pore-scale dispersion in both two- and three-dimensional applications.     

Stochastic modeling of solute transport in 3-D heterogeneous porous media with random source condition

During probabilistic analysis of flow and transport in porous media, the uncertainty due to spatial heterogeneity of governing parameters are often taken into account. The randomness in the source conditions also play a major role on the stochastic behavior in distribution of the dependent variable. The present paper is focused on studying the effect of both uncertainty in the governing system parameters as well as the input source conditions. Under such circumstances, a method is proposed which combines with stochastic finite element method (SFEM) and is illustrated for probabilistic analysis of concentration distribution in a 3-D heterogeneous porous media under the influence of random source condition. In the first step SFEM used for probabilistic solution due to spatial heterogeneity of governing parameters for a unit source pulse. Further, the results from the unit source pulse case have been used for the analysis of multiple pulse case using the numerical convolution when the source condition is a random process. The source condition is modeled as a discrete release of random amount of masses at fixed intervals of time. The mean and standard deviation of concentration is compared for the deterministic and the stochastic system scenarios as well as for different values of system parameters. The effect of uncertainty of source condition is also demonstrated in terms of mean and standard deviation of concentration at various locations in the domain.     

A master equation for reactive solute transport in porous media

The mean value of a density of a cloud of points described by a generalized Liouville equation associated with a convection dispersion equation governing adsorbing solute transport yields a joint concentration probability density. The general technique can be applied for either linear or nonlinear adsorption; here the application is restricted to linear adsorption in one-dimensional transport. The equation generated for the joint concentration probability density is in the general form of a Fokker-Planck equation, but with a suitable coordinate transformation, it is possible to represent it as a diffusion equation with variable coefficients.     

Transport upscaling using multi-rate mass transfer in three-dimensional highly heterogeneous porous media

A methodology for transport upscaling of three-dimensional highly heterogeneous formations is developed and demonstrated. The overall approach requires a prior hydraulic conductivity upscaling using an interblock-centered full-tensor Laplacian-with-skin method followed by transport upscaling. The coarse scale transport equation includes a multi-rate mass transfer term to compensate for the loss of heterogeneity inherent to all upscaling processes. The upscaling procedures for flow and transport are described in detail and then applied to a three-dimensional highly heterogeneous synthetic example. The proposed approach not only reproduces flow and transport at the coarse scale, but it also reproduces the uncertainty associated with the predictions as measured by the ensemble variability of the breakthrough curves.     

Modeling the effects of nonlinear equilibrium sorption on the transport of solute plumes in saturated heterogeneous porous media

Transport of sorbing solutes in 2D steady and heterogeneous flow fields is modeled using a particle tracking random walk technique. The solute is injected as an instantaneous pulse over a finite area. Cases of linear and Freundlich sorption isotherms are considered. Local pore velocity and mechanical dispersion are used to describe the solute transport mechanisms at the local scale. This paper addresses the impact of the degree of heterogeneity and correlation lengths of the log-hydraulic conductivity field as well as negative correlation between the log-hydraulic conductivity field and the log-sorption affinity field on the behavior of the plume of a sorbing chemical. Behavior of the plume is quantified in terms of longitudinal spatial moments: center-of-mass displacement, variance, 95% range, and skewness. The range appears to be a better measure of the spread in the plumes with Freundlich sorption because of plume asymmetry. It has been found that the range varied linearly with the travelled distance, regardless of the sorption isotherm. This linear relationship is important for extrapolation of results to predict behavior beyond simulated times and distances. It was observed that the flow domain heterogeneity slightly enhanced the spreading of nonlinearly sorbing solutes in comparison to that which occurred for the homogeneous flow domain, whereas the spreading enhancement in the case of linear sorption was much more pronounced. In the case of Freundlich sorption, this enhancement led to further deceleration of the solute plume movement as a result of increased retardation coefficients produced by smaller concentrations. It was also observed that, except for plumes with linear sorption, correlation between the hydraulic conductivity and the sorption affinity fields had minimal effect on the spatial moments of solute plumes with nonlinear sorption.     

Comparison of laboratory-scale solute transport visualization experiments with numerical simulation using cross-bedded sandstone

Using a slab of Massillon Sandstone, laboratory-scale solute tracer experiments were carried out to test numerical simulations using the Advection–Dispersion Equation (ADE). While studies of a similar nature exist, our work differs in that we combine: (1) experimentation in naturally complex geologic media, (2) X-ray absorption imaging to visualize and quantify two-dimensional solute transport, (3) high resolution transport property characterization, with (4) numerical simulation. The simulations use permeability, porosity, and solute concentration measured to sub-centimeter resolution. While bulk breakthrough curve characteristics were adequately matched, large discrepancies exist between the experimental and simulated solute concentration fields. Investigation of potential experimental errors suggests that the failure to fit solute concentration fields may lie in loss of intricate connectivity within the cross-bedded sandstone occurring at scales finer than our property characterization measurements (i.e., sub-centimeter).     

Stochastic modelling of groundwater flow and solute transport in aquifers

This paper presents an introductory overview of recently developed stochastic theories for tackling spatial variability problems in predicting groundwater flow and solute transport. Advantages and limitations of the theories are discussed. Lastly, strategies based on the stochastic approaches to predict solute transport in aquifers are recommended.     

A multidimensional streamline-based method to simulate reactive solute transport in heterogeneous porous media

We present a new streamline-based numerical method for simulating reactive solute transport in porous media. The key innovation of the method is that both longitudinal and transverse dispersion are incorporated accurately without numerical dispersion. Dispersion is approximated in a flow-oriented grid using a combination of a one-dimensional finite difference scheme and a meshless approximation. In contrast to previous hybrid alternatives to incorporate dispersion in streamline-based simulations, the proposed scheme does not require a grid and, hence, it does not introduce numerical dispersion. In addition, the proposed scheme eliminates numerical oscillations and negative concentration values even when the dispersion tensor includes the off-diagonal coefficients and the flow field is non-uniform. We demonstrate that for a set of two- and three-dimensional benchmark problems, the new proposed streamline-based formulation compares favorably to two state of the art finite volume and hybrid Eulerian–Lagrangian solvers.     

On different numerical inverse Laplace methods for solute transport problems

Numerical inversion is required when Laplace transform cannot be inverted analytically by manipulating tabled formulas of special cases. However, the numerical inverse Laplace transform is generally an ill-posed problem, and there is no universal method which works well for all problems. In this study, we selected seven commonly used numerical inverse Laplace transform methods to evaluate their performance for dealing with solute transport in the subsurface under uniform or radial flow condition. Such seven methods included the Stehfest, the de Hoog, the Honig–Hirdes, the Talbot, the Weeks, the Simon and the Zakian methods. We specifically investigated the optimal free parameters of each method, including the number of terms used in the summation and the numerical tolerance. This study revealed that some commonly recommended values of the free parameters in previous studies did not work very well, especially for the advection-dominated problems. Instead, we recommended new values of the free parameters for some methods after testing their robustness. For the radial dispersion, the de Hoog, the Talbot, and the Simon methods worked very well, regardless of the dispersion-dominated or advection-dominated situations. The Weeks method can be used to solve the dispersion-dominated problems, but not the advection-dominated problems. The Stehfest, the Honig–Hirdes, and the Zakian methods were recommended for the dispersion-dominated problems. The Zakian method was efficient, while the de Hoog method was time-consuming under radial flow condition. Under the uniform flow condition, all the methods could present somewhat similar results when the free parameters were given proper values for dispersion-dominated problems; while only the Simon method, the Weeks method, and the de Hoog method worked well for advection-dominated problems.     

Flow and transport through two-dimensional isotropic media of binary conductivity distribution. Part 2: NUMERICAL simulations and comparison with theoretical results

In the present part the results of numerical simulations of flow and transport in media made up from circular inclusions of conductivity K that are submerged in a matrix of conductivity K ₀, subjected to uniform mean velocity, are presented. This is achieved for a few values of =K/K ₀ (0.01, 0.1 and 10) and of the volume fraction n (0.05, 0.1 and 0.2). The numerical simulations (NS) are compared with the analytical approximate models presented in Part 1: the composite elements (CEA), the effective medium (EMA), the dilute system (DSA) and the first-order in the logconductivity variance (FOA). The comparison is made for the longitudinal velocity variance and for the longitudinal macrodispersivity. This is carried out for n<0.2, for which the theoretical and simulation models represent the same structure of random and independent inclusions distribution. The main result is that transport is quite accurately modeled by the EMA and CEA for low , for which _L is large, whereas in the case of =10, the EMA matches the NS for n<0.1. The first-order approximation is quite far apart from the NS for the values of examined . This material is based upon work supported by the National Science Foundation under Grant No. 0218914. Authors also wish to thank the Center of Computational Research, University at Buffalo for assistance in running numerical simulations.     

Characterizing heterogeneous soil water flow and solute transport using information measures

Heterogeneous water flow and solute transport in soils are an important phenomenon and difficult to be characterized. The objectives of this study were to investigate the heterogeneity of solute transport related to heterogeneous soil water flow using dye infiltration experiments, and to characterize heterogeneous water flow and solute transport in soils using the information theory. Field experiments of dye infiltration were performed in four plots. Various information measures were applied to characterize information content and complexity of water flow and solute transport in soils. Information contents and complexities of the maximum and apparent infiltration depths, and the mean and standard deviation of concentrations in the vertical direction of the plots were calculated. More heterogeneous processes of soil water flow and transport result in higher information/complexity values. The probability distributions of mean concentration were similar to those of the corresponding apparent infiltration depths for the plots, indicating that heterogeneity of dye concentrations was closely related to that of soil water flow. However, the range of information entropy and complexity of the water flow sequences was much narrower than that of the sequences of the concentrations. The results suggested that the transport processes were more heterogeneous than the water flow processes. Compared with the probability distributions of flow parameters, the information measures appeared to be a more versatile tool to describe flow and transport heterogeneities in soils.