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.     

Transport of kinetically sorbing solutes in heterogeneous sediments with multimodal conductivity and hierarchical organization across scales

Solute transport in subsurface environments is controlled by geological heterogeneity over multiple scales. In reactive transport characterized by a low Damköhler number, it is also controlled by the rate of kinetic mass transfer. A theory for addressing the impact of sedimentary texture on the transport of kinetically sorbing solutes in heterogeneous porous formations is derived using the Lagrangian-based stochastic methodology. The resulting model represents the hierarchical organization of sedimentary textures and associated modes of log conductivity (K) for sedimentary units through a hierarchical Markov Chain. The model characterizes kinetic sorption using a spatially uniform linear reversible rate expression. Our main interest is to investigate the effect of sorption kinetics relative to the effects of K heterogeneity on the dispersion of a reactive plume. We study the contribution of each scale of stratal architecture to the dispersion of kinetically sorbing solutes in the case of a low Damköhler number. Examples are used to demonstrate the time evolution and relative contributions of the auto- and cross-transition probability terms to dispersion. Our analysis is focused on the model sensitivity to the parameters defined at each hierarchical level (scale) including the integral scales of K spatial correlation, the anisotropy ratio, the indicator correlation scales, and the contrast in mean K between facies defined at different scales. The results show that the anisotropy ratio and integral scales of K have negligible effect upon the longitudinal dispersion of sorbing solutes. Furthermore, dispersion of sorbing solutes depends mostly on indicator correlation scales, and the contrast of the mean conductivity between units at different scales.     

The role of moisture fluctuations in unsaturated transport

Within the framework of stochastic theory and the spectral perturbation techniques, three-dimensional dispersion in partially saturated soils with a finite correlation scale of log-hydraulic conductivity is analyzed. The effects of spatial variability of the moisture distribution parameter on the asymptotic spreading behavior of a unsaturated solute plume are assessed. This is accomplished by comparing two asymptotic macrodispersivities and two variance of solute concentration, obtained for a constant moisture content and spatially varied moisture, respectively.     

Modeling plume behavior for nonlinearly sorbing solutes in saturated homogeneous porous media

Transport of a sorbing solute in a two-dimensional steady and uniform flow field is modeled using a particle tracking random walk method. The solute is initially introduced from an instantaneous point source. Cases of linear and nonlinear sorption isotherms are considered. Local pore velocity and mechanical dispersion are used to describe the solute transport mechanisms at the local scale. The numerical simulation of solute particle transport yields the large scale behavior of the solute plume. Behavior of the plume is quantified in terms of the center-of-mass displacement distance, relative velocity of the center-of-mass, mass breakthrough curves, spread variance, and longitudinal skewness. The nonlinear sorption isotherm affects the plume behavior in the following way relative to the linear isotherm: (1) the plume velocity decreases exponentially with time; (2) the longitudinal variance increases nonlinearly with time; (3) the solute front is steepened and tailing is enhanced     

Stochastic study of solute transport in a nonstationary medium

A Lagrangian stochastic approach is applied to develop a method of moment for solute transport in a physically and chemically nonstationary medium. Stochastic governing equations for mean solute flux and solute covariance are analytically obtained in the first-order accuracy of log conductivity and/or chemical sorption variances and solved numerically using the finite-difference method. The developed method, the numerical method of moments (NMM), is used to predict radionuclide solute transport processes in the saturated zone below the Yucca Mountain project area. The mean, variance, and upper bound of the radionuclide mass flux through a control plane 5 km downstream of the footprint of the repository are calculated. According to their chemical sorption capacities, the various radionuclear chemicals are grouped as nonreactive, weakly sorbing, and strongly sorbing chemicals. The NMM method is used to study their transport processes and influence factors. To verify the method of moments, a Monte Carlo simulation is conducted for nonreactive chemical transport. Results indicate the results from the two methods are consistent, but the NMM method is computationally more efficient than the Monte Carlo method. This study adds to the ongoing debate in the literature on the effect of heterogeneity on solute transport prediction, especially on prediction uncertainty, by showing that the standard derivation of solute flux is larger than the mean solute flux even when the hydraulic conductivity within each geological layer is mild. This study provides a method that may become an efficient calculation tool for many environmental projects.     

The effect of non divergence-free velocity fields on field scale ground water solute transport

Solute plume subjected to field scale hydraulic conductivity heterogeneity shows a large dispersion/macrodispersion, which is the manifestation of existing fields scale heterogeneity on the solute plume. On the other hand, due to the scarcity of hydraulic conductivity measurements at field scale, hydraulic conductivity heterogeneity can only be defined statistically, which makes the hydraulic conductivity a random variable/function. Random hydraulic conductivity as a parameter in flow equation makes the pore flow velocity also random and the ground water solute transport equation is a stochastic differential equation now. In this study, the ensemble average of stochastic ground water solute transport equation is taken by the cumulant expansion method in order to upscale the laboratory scale transport equation to field scale by assuming pore flow velocity is a non stationary, non divergence-free and unsteady random function of space and time. Besides the stochastic explanation of macrodispersion and the velocity correction term obtained by Kavvas and Karakas (J Hydrol 179:321–351, 1996) before a new velocity correction term, which is a function of mean pore flow velocity divergence, is obtained in this study due to strict second order cumulant expansion (without omitting any term after the expansion) performed. The significance of the new velocity correction term is investigated on a one dimensional transport problem driven by a density dependent flow field.     

Do conservative solutes migrate at average pore-water velocity?

According to common understanding, the advective velocity of a conservative solute equals the average linear pore-water velocity. Yet direct monitoring indicates that the two velocities may be different in heterogeneous media. For example, at the Camp Dodge, Iowa, site the advective velocity of discrete Cl- plumes was less than one tenth of the average pore-water velocity calculated from Darcy's law using the measured hydraulic gradient, effective porosity, and hydraulic conductivity (K) from large-scale three-dimensional (3D) techniques, e.g., pumping tests. Possibly, this difference reflects the influence of different pore systems, if the K relevant to transient solute flux is influenced more by lower-K heterogeneity than a steady or quasi-steady water flux. To test this idea, tracer tests were conducted under controlled laboratory conditions. Under one-dimensional flow conditions, the advective velocity of discrete conservative solutes equaled the average pore-water velocity determined from volumetric flow rates and Darcy's law. In a larger 3D flow system, however, the same solutes migrated at approximately 65% of the average pore-water velocity. These results, coupled with direct observation of dye tracers and their velocities as they migrated through both homogeneous and heterogeneous sections of the same model, demonstrate that heterogeneity can slow the advective velocity of discrete solute plumes relative to the average pore-water velocity within heterogeneous 3D flow sytems.     

Transport of variable-density solute plumes in beach aquifers in response to oceanic forcing

A comprehensive numerical study was undertaken to investigate transport of a variable-density, conservative solute plume in an unconfined coastal aquifer subject to high and low frequency oceanic forcing. The model combined variable-density saturated flow for groundwater and solute transport, and wave hydrodynamics from a 2D Navier–Stokes solver. A sinusoidal tidal signal was specified by implementing time-varying heads at the seaward boundary. The solute plume behavior was investigated under different oceanic forcing conditions: no forcing, waves, tide, and combined waves and tide. For each forcing condition, four different injected solute densities (freshwater, brackish water, seawater, brine) were used to investigate the effects of density on the transport of the injected plume beneath and across the beach face. The plume’s low-order spatial moments were computed, viz., mass, centroid, variance and aspect ratio. The results confirmed that both tide- and wave-forcing produce an upper saline plume beneath the beach face in addition to the classical saltwater wedge. For the no-forcing and tide-only cases (during rising tides), an additional small circulation cell below the beach face was observed. Oceanic forcing affects strongly the solute plume’s flow path, residence time and discharge rate across the beach face, as well as its spreading. For the same oceanic forcing, solute plumes with different densities follow different trajectories from the source to the discharge location (beach face). The residence time and plume spreading increased with plume density. It was concluded that simulations that neglect the effect of waves or tides cannot reproduce accurately solute plume dispersion and also, in the case of coasts with small waves or tides, the solute residence time in the aquifer.     

A numerical method of moments for solute transport in physically and chemically nonstationary formations: linear equilibrium sorption with random K<Subscript>d</Subscript>

A Lagrangian perturbation method is applied to develop a method of moments for solute flux through a three-dimensional nonstationary flow field. The flow nonstationarity stems from medium nonstationarity and internal and external boundaries of the study domain. The solute flux is described as a space-time process where time refers to the solute flux breakthrough through a control plane (CP) at some distance downstream of the solute source and space refers to the transverse displacement distribution at the CP. The analytically derived moment equations for solute transport in a nonstationarity flow field are too complicated to solve analytically, a numerical finite difference method is implemented to obtain the solutions. This approach combines the stochastic model with the flexibility of the numerical method to boundary and initial conditions. The developed method is applied to study the effects of heterogeneity and nonstationarity of the hydraulic conductivity and chemical sorption coefficient on solute transport. The study results indicate all these factors will significantly influence the mean and variance of solute flux.     

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.     

Dilution of non-reactive tracers in variably saturated sandy structures

Based on a dye tracer experiment in a sand tank we addressed the problem of local dispersion of conservative tracers in the unsaturated zone. The sand bedding was designed to have a defined spatial heterogeneity with a strong anisotropy. We estimated the parameters that characterize the local dispersion and dilution from concentration maps of a high spatial and temporal resolution obtained by image analysis. The plume spreading and mixing behavior was quantified on the basis of the coefficient of variation of the concentration and of the dilution index. The heterogeneous structure modified the flow pattern depending on water saturation. The shape of the tracer plumes revealed the structural signature of the sand bedding at low saturation only. In this case pronounced preferential flow was observed. At higher flow rates the structure remained hidden by a spatially almost homogeneous behavior of the plumes. In this context, we mainly discuss the mechanism of re-distributing a finite mass of inert solutes over a large volume, due to macro- and micro-heterogeneities of the structure.     

A Critical Analysis of Transverse Dispersivity Field Data

Transverse dispersion, or tracer spreading orthogonal to the mean flow direction, which is relevant e.g, for quantifying bio-degradation of contaminant plumes or mixing of reactive solutes, has been studied in the literature less than the longitudinal one. Inferring transverse dispersion coefficients from field experiments is a difficult and error-prone task, requiring a spatial resolution of solute plumes which is not easily achievable in applications. In absence of field data, it is a questionable common practice to set transverse dispersivities as a fraction of the longitudinal one, with the ratio 1/10 being the most prevalent. We collected estimates of field-scale transverse dispersivities from existing publications and explored possible scale relationships as guidance criteria for applications. Our investigation showed that a large number of estimates available in the literature are of low reliability and should be discarded from further analysis. The remaining reliable estimates are formation-specific, span three orders of magnitude and do not show any clear scale-dependence on the plume traveled distance. The ratios with the longitudinal dispersivity are also site specific and vary widely. The reliability of transverse dispersivities depends significantly on the type of field experiment and method of data analysis. In applications where transverse dispersion plays a significant role, inference of transverse dispersivities should be part of site characterization with the transverse dispersivity estimated as an independent parameter rather than related heuristically to longitudinal dispersivity.     

On the effects of preferential or barrier flow features on solute plumes in permeable porous media

Despite that discrete flow features (DFFs, e.g. fractures and faults) are common features in the subsurface, few studies have explored the influence of DFFs on solute plumes in otherwise permeable rocks (e.g. sandstone, limestone), compared to low-permeability rock settings (e.g. granite and basalt). DFFs can provide preferential flow pathways (i.e. ‘preferential flow features’; PFFs), or can act to impede flow (i.e. ‘barrier flow features’; BFFs). This research uses a simple analytical expression and numerical modelling to explore how a single DFF influences the steady-state distributions of solute plumes in permeable aquifers. The analysis quantifies the displacement and widening (or narrowing) of a steady-state solute plume as it crosses a DFF in idealised, 1 × 1 m moderately permeable rock aquifers. Previous research is extended by accounting for DFFs as 2D flow features, and including BFF situations. A range of matrix-DFF permeability ratios (0.01 to 100) and DFF apertures (0.25 mm to 2 cm), typical of sedimentary aquifers containing medium-to-large fractures, are considered. The results indicate that for the conceptual models considered here, PFFs typically have a more significant influence on plume distributions than BFFs, and the impact of DFFs on solute plumes generally increases with increasing aperture. For example, displacement of peak solute concentration caused by DFFs exceeds 20 cm in some PFF cases, compared to a maximum of 0.64 cm in BFF cases. PFFs widen plumes up to 9.7 times, compared to a maximum plume widening of 2.0 times in BFF cases. Plumes crossing a PFF are less symmetrical, and peak solute concentrations beneath PFFs are up to two orders of magnitude lower than plumes in BFF cases. This study extends current knowledge of the attenuating influence of DFFs in otherwise permeable rocks on solute plume characteristics, through evaluation of 2D flow effects in DFFs for a variety of DFF apertures, and by considering BFF situations.     

Solute transport in heterogeneous media: The impact of anisotropy and non-ergodicity in risk assessment

Conceptual model selection is a key issue in risk assessment studies. We analyze the effect of a number of conceptual aspects related to solute transport in two-dimensional heterogeneous media. The main issues addressed are non-ergodicity, anisotropy in the correlation structure of the transmissivity field, and dispersion at the local scale. In particular, we study the development of a solute plume when mean flow is oriented at an angle with respect to the principal directions of anisotropy. The study is carried out in a Lagrangian framework using Monte Carlo analysis. Of special interest is the evolution of individual plumes. A number of aspects are analyzed, namely the location of the center of mass for each plume and the different ways to compute the angles that the main axes of the plume develop with respect to the direction of the mean flow. Stochastic theories based upon ergodicity conclude that the plume gets oriented in the mean flow direction. In our non-ergodic simulations, the mean of the offset angles, for each individual plume in each particular realization, is offset from the mean flow direction towards the direction of maximum anisotropy. If, instead, the analysis is performed on the ensemble plume (superposition of all different simulations), it is then found oriented closer to the direction of the mean flow than the average offset angle for the different plumes considered separately. This last result adds an extra word of caution to the use of ensemble averaged values in solute transport studies. Serious implications for risk assessment follow from the conceptual model adopted. First, in any single realization there will a large uncertainty in locating the plume at any given time; second, real dilution would be less than what would be expected if the macrodispersion values obtained for ergodic conditions were applied; third, the volume that is affected by a non-zero concentration is smaller than that predicted from macrodispersion concepts; fourth, the orientation of the plume does not correspond to that of the mean flow; and fifth, accounting for local dispersion helps reducing uncertainty.     

Effective dispersion in conditioned transmissivity fields

We analyze the impact of conditioning to measurements of hydraulic transmissivity on the transport of a conservative solute. The effects of conditioning on solute transport are widely discussed in the literature, but most of the published works focuses on the reduction of the uncertainty in the prediction of the plume dispersion. In this study both ensemble and effective plume moments are considered for an instantaneous release of a solute through a linear source normal to the mean flow direction, by taking into account different sizes of the source. The analysis, involving a steady and spatially inhomogeneous velocity field, is developed by using the stochastic finite element method. Results show that conditioning reduces the ensemble moment in comparison with the unconditioned case, whereas the effective dispersion may increase because of the contribution of the spatial moments related to the lack of stationarity in the flow field. As the number of conditioning points increases, this effect increases and it is significant in both the longitudinal and transverse directions. Furthermore, we conclude that the moment derived from data collected in the field can be assessed by the conditioned second-order spatial moment only with a dense grid of measured data, and it is manifest for larger initial lengths of the plume. Nevertheless, it seems very likely that the actual dispersion of the plume may be underestimated in practical applications.     

Eulerian spatial moments for solute transport in three-dimensional heterogeneous,dual-permeability media

A Eulerian analytical method is developed for nonreactive solute transport in heterogeneous, dual-permeability media where the hydraulic conductivities in fracture and matrix domains are both assumed to be stochastic processes. The analytical solution for the mean concentration is given explicitly in Fourier and Laplace transforms. Instead of using the fast fourier transform method to numerically invert the solution to real space (Hu et al., 2002), we apply the general relationship between spatial moments and concentration (Naff, 1990; Hu et al., 1997) to obtain the analytical solutions for the spatial moments up to the second for a pulse input of the solute. Owing to its accuracy and efficiency, the analytical method can be used to check the semi-analytical and Monte Carlo numerical methods before they are applied to more complicated studies. The analytical method can be also used during screening studies to identify the most significant transport parameters for further analysis. In this study, the analytical results have been compared with those obtained from the semi-analytical method (Hu et al., 2002) and the comparison shows that the semi-analytical method is robust. It is clearly shown from the analytical solution that the three factors, local dispersion, conductivity variation in each domain and velocity convection flow difference in the two domains, play different roles on the solute plume spreading in longitudinal and transverse directions. The calculation results also indicate that when the log-conductivity variance in matrix is 10 times less than its counterpart in fractures, it will hardly influence the solute transport, whether the conductivity field is matrix is treated as a homogeneous or random field.     

Assessment of uncertainty associated with the estimation of well catchments by moment equations

Non-local stochastic moment equations are used successfully to analyze groundwater flow in randomly heterogeneous media. Here we present a moment equations-based approach to quantify the uncertainty associated with the estimation of well catchments. Our approach is based on the development of a complete second order formalism which allows obtaining the first statistical moments of the trajectories of conservative solute particles advected in a generally non-uniform groundwater flow. Approximate equations of moments of particles’ trajectories are then derived on the basis of a second order expansion in terms of the standard deviation of the aquifer log hydraulic conductivity. Analytical expressions are then obtained for the predictors of locations of mean stagnation points, together with their associated uncertainties. We implement our approach on heterogeneous media in bounded two-dimensional domains, with and without including the effect of conditioning on hydraulic conductivity information. The impact of domain size, boundary conditions, heterogeneity and non-stationarity of hydraulic conductivity on the prediction of a well catchment is explored. The results are compared against Monte Carlo simulations and semi-analytical solutions available in the literature. The methodology is applicable to both infinite and bounded domains and is free of distributional assumptions (and so applies to both Gaussian and non-Gaussian log hydraulic conductivity fields) and formally includes the effect of conditioning on available information.     

Stochastic analysis of unsaturated transport in soils with fractal log-conductivity distribution

Within the framework of stochastic theory and the spectral perturbation techniques, three-dimensional dispersion in partially saturated soils with fractal log hydraulic conductivity distribution is analyzed. Our analysis is focused on the impact of fractal dimension of log hydraulic conductivity distribution, local dispersivity, and unsaturated flow parameters, such as the soil poresize distribution parameter and the moisture distribution parameter, on the spreading behavior of solute plume and the concentration variance. Approximate analytical solutions to the stochastic partial differential equations are derived for the variance of asymptotic solute concentration and asymptotic macrodispersivities.