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.     

Transect study on solute transport in a macroporous soil

Solute transport experiments using a non-reactive tracer were conducted on short, undisturbed, saturated columns of a sandy loam soil. All columns, 20 cm in diameter and 20 cm long, were collected along a transect of 35 m. Most of the soil columns had pre-existing macropores. The columns were leached at a steady flow-rate under ponding conditions. The resulting breakthrough curves (BTCs) showed a large heterogeneity. Several of the BTCs displayed early breakthrough and long tailing. All the data were interpreted in terms of dimensional time moments, the classical convection-dispersion equation (CDE) and the mobile-immobile transport model (MIM). Experimental time moments were found to vary significantly among the different BTCs. Analysis of the time moments also revealed that the variance of the field-scale BTC was several times larger than the average of the local-scale variance. The pore water velocity v and dispersion coefficient D were obtained by fitting the CDE to the local-scale BTCs, resulting in an average dispersivity of 7·4 cm. Frequency distributions for the CDE parameters v and D were equally well described by a normal or log-normal probability density function (pdf). When a log-normal pdf for D is considered, the variance of the log_e transformed D values (σ_{ln D}²) was found to be 2·1. For the MIM model, two additional parameters were fitted: the fraction of mobile water, θ_m/θ, and the first-order mass transfer coefficient, α. The MIM was more successful in describing the data than the CDE transport model. For the MIM model, the average dispersivity was about 2 cm. The MIM parameters v, D and θ_m/θ were best described by a log-normal pdf rather than a normal pdf. Only the parameter α was better described by a normal pdf. Mobile water fractions, θ_m/θ ranged from 0·01 to 0·98, with a mean of 0·43 (based on a log-normal pdf). When the CDE and MIM were applied to the data, the fitted pore water velocities, v, compared favourably with the effective pore water velocities, v_eff, obtained from moment analysis.     

Modeling multiaquifer wells with MODFLOW

Multiaquifer wells, i.e., wells that are open across more than one aquifer, can have a profound effect on the hydraulics of a ground water system. These wells change the physical system by establishing direct hydraulic links between otherwise isolated strata. Several methods are available to simulate multiaquifer wells in the context of comprehensive ground water flow simulators. In this paper, we review four methods to represent multiaquifer wells with the widely used code MODFLOW. These methods include a specialized code developed, but never formally released, by the U.S. Geological Survey (USGS), the Multi-Aquifer Well (MAW1) Package. An expanded implementation of the techniques in the MAW1 Package has been incorporated in the Multi-Node Well Package released recently by the USGS (Halford and Hanson 2002). We examine the performance of the methods in the context of a benchmarking study against the analytical solutions of Papadopulos (1966) and Sokol (1963). Our results demonstrate that results obtained with the MAW1 Package closely match exact solutions for pumping and nonpumping conditions, using both coarse and refined grids.     

Effect of snow structure on water flow and solute transport

Water flow through a melting snow pack modifies its structure and stability and affects the release of water and nutrients into soils and surface waters. Field and laboratory observations indicate a large spatial variability on various scales of the liquid water content and flow, a dominant system feature currently not included in numerical models. We investigated experimentally water and dye tracer movement through microstructurally different snow pack horizons and the persistence of preferential flow paths. Naturally rounded snow of varying grain size was artificially packed to obtain well known conditions by sieving it into rectangular bins. Surface melt was induced with infrared lamps. The flow paths were visualized with tracers and liquid water content was monitored with time domain reflectometry probes. Vertical cuts through the snow pack were imaged. The dye tracer patterns allowed the two flow regimes ‘matrix flow’ and ‘preferential flow’ to be distinguished. Matrix flow is apparently dominated by film and capillary flow in the unsaturated snow matrix. The capillary barrier effect at a boundary between a fine over a coarse textured layer on matrix flow in snow was confirmed. In contrast, preferential flow appears as well‐defined flow fingers that advance from 0·1 to 1 cm s^?1. During a melt phase, the advancing flow fingers enlarge and are only partially time invariant. It remains to be shown whether the continuum concept, including the Darcy–Buckingham law is apt to describe the extremely non‐linear nature of water flow and the travel time of solutes in snow under conditions of melt water percolation. Probably, snow packs that include faceted crystals and large variations in bulk density, feature more pronounced capillary barriers and preferential flow triggering, but also stronger impeding of fingers by lateral dispersion. Further, triggering and persistence of preferential flow is complicated by the usually transient infiltration rate. Copyright © 2004 John Wiley & Sons, Ltd.     

Impact of spatial rainfall resolution on point-source solute transport modelling

Numerical simulations of unsaturated solute transport from point sources were carried out using HYDRUS-1D. Three different soil types cropped with spring wheat were considered at three different locations in Sweden: Malmö, Norrköping and Petisträsk. Two types of rainfall data were used, point-scale raingauge measurements and a gauge-adjusted weather radar product at four spatial resolutions, 2 × 2, 6 × 6, 10 × 10 and 14 × 14 km². The results showed that differences in the mean solute transport depths were small and not significant, with the exception of Petisträsk. Maximum transport depths were in most cases significantly larger using raingauge data compared to radar data. The results showed that using areal-averaged rainfall input will give solute transport estimations close to those using point-measured data. This shows the great potential in using radar-measured rainfall data in small-scale hydrological applications.     

Density effects on solute release from streambeds

Contaminants that entered the streambed during previous surface water pollution events can be released to the stream, causing secondary pollution of the stream and impacting its eco-environmental condition. By means of laboratory experiments and numerical simulations, we investigated density effects on the release of solute from periodic bedforms. The results show that solute release from the upper streambed is driven by bedform-induced convection, and that density effects generally inhibit the solute release from the lower streambed. Density gradients modify the pore water flow patterns and form circulating flows in the area of lower streambed. The formation of circulating flows is affected by density gradients associated with the solute concentration and horizontal pressure gradients induced by stream slope. The circulating flows near the bottom of the streambed enhance mixing of the hyporheic zone and the ambient flow zone.     

Effect of anisotropy on solute transport in degraded fen peat soils

Peat soils are heterogeneous, anisotropic porous media. Compared to mineral soils, there is still limited understanding of physical and solute transport properties of fen peat soils. In this study, we aimed to explore the effect of soil anisotropy on solute transport in degraded fen peat. Undisturbed soil cores, taken in vertical and horizontal direction, were collected from one drained and one restored fen peatland both in a comparable state of soil degradation. Saturated hydraulic conductivity (K _s) and chemical properties of peat were determined for all soil cores. Miscible displacement experiments were conducted under saturated steady state conditions using potassium bromide as a conservative tracer. The results showed that (1) the K _s in vertical direction (K _sv) was significantly higher than that in horizontal direction (K_sh), indicating that K _s of degraded fen peat behaves anisotropically; (2) pronounced preferential flow occurred in vertical direction with a higher immobile water fraction and a higher pore water velocity; (3) the 5% arrival time (a proxy for the strength of preferential flow) was affected by soil anisotropy as well as study site. A strong correlation was found between 5% arrival time and dispersivity, K _s and mobile water fraction; (4) phosphate release was observed from drained peat only. The impact of soil heterogeneity on phosphate leaching was more pronounced than soil anisotropy. The soil core with the strongest preferential flow released the highest amount of phosphate. We conclude that soil anisotropy is crucial in peatland hydrology but additional research is required to fully understand anisotropy effects on solute transport.     

Imaging hyporheic zone solute transport using electrical resistivity

Traditional characterization of hyporheic processes relies upon modelling observed in‐stream and subsurface breakthrough curves to estimate hyporheic zone size and infer exchange rates. Solute data integrate upstream behaviour and lack spatial coverage, limiting our ability to accurately quantify spatially heterogeneous exchange dynamics. Here, we demonstrate the application of near‐surface electrical resistivity imaging (ERI) methods, coupled with experiments using an electrically conductive stream tracer (dissolved NaCl), to provide in situ imaging of spatial and temporal dynamics of hyporheic exchange. Tracer‐labelled water in the stream enters the hyporheic zone, reducing electrical resistivity in the subsurface (to which subsurface ERI is sensitive). Comparison of background measurements with those recording tracer presence provides distributed characterization of hyporheic area (in this application, ∼0·5 m²). Results demonstrate the first application of ERI for two‐dimensional imaging of stream‐aquifer exchange and hyporheic extent. Future application of this technique will greatly enhance our ability to quantify processes controlling solute transport and fate in hyporheic zones, and provide data necessary to inform more complete numerical models. Copyright © 2010 John Wiley & Sons, Ltd.     

Multicomponent solute transport with sorption and soluble complexation

A numerical algorithm for predicting the migration of multiple subsurface pollutants has been developed accounting for dispersion, convection, soluble complexation and solid phase accumulations (sorption). The basis of the model is a finite element solution of the mass transport equation. The essence of the algorithm is the treatment of the sorption terms as implicit functions of the total soluble concentrations facilitating a general and modularized treatment of solution phase and solid phase chemistry. Examples are presented illustrating the effects of soluble complexation and competitive sorption on the transport of multicomponent solutions.     

Discretizing the fracture-matrix interface to simulate solute transport

This article examines the required spatial discretization perpendicular to the fracture-matrix interface (FMI) for numerical simulation of solute transport in discretely fractured porous media. The discrete-fracture, finite-element model HydroGeoSphere ( Therrien et al. 2005 ) and a discrete-fracture implementation of MT3DMS ( Zheng 1990 ) were used to model solute transport in a single fracture, and the results were compared to the analytical solution of Tang et al. (1981) . To match analytical results on the relatively short timescales simulated in this study, very fine grid spacing perpendicular to the FMI of the scale of the fracture aperture is necessary if advection and/or dispersion in the fracture is high compared to diffusion in the matrix. The requirement of such extremely fine spatial discretization has not been previously reported in the literature. In cases of high matrix diffusion, matching the analytical results is achieved with larger grid spacing at the FMI. Cases where matrix diffusion is lower can employ a larger grid multiplier moving away from the FMI. The very fine spatial discretization identified in this study for cases of low matrix diffusion may limit the applicability of numerical discrete-fracture models in such cases.     

An evaluation of conditioning data for solute transport prediction

The large and diverse body of subsurface characterization data generated at a field research site near Oyster, Virginia, provides a unique opportunity to test the impact of conditioning data of various types on predictions of flow and transport. Bromide breakthrough curves (BTCs) were measured during a forced-gradient local-scale injection experiment conducted in 1999. Observed BTCs are available at 140 sampling points in a three-dimensional array within the transport domain. A detailed three-dimensional numerical model is used to simulate breakthrough curves at the same locations as the observed BTCs under varying assumptions regarding the character of hydraulic conductivity spatial distributions, and variable amounts and types of conditioning data. We present comparative results of six cases ranging from simple (deterministic homogeneous models) to complex (stochastic indicator simulation conditioned to cross-borehole geophysical observations). Quantitative measures of model goodness-of-fit are presented. The results show that conditioning to a large number of small-scale measurements does not significantly improve model predictions, and may lead to biased or overly confident predictions. However, conditioning to geophysical interpretations with larger spatial support significantly improves the accuracy and precision of model predictions. In all cases, the effects of model error appear to be significant in relation to parameter uncertainty.     

Temporal variability of solute transport under vadose zone conditions

The heterogeneity of the solute flux field in the horizontal plane at the field scale has been documented in several field studies. On the other hand, little information is available on the persistence of certain solute transport scenarios over consecutive infiltration cycles. This study was initiated to analyse the recurrence of solute leaching behaviour as estimated in two soil column tests emphasizing the preferential flow phenomenon. Twenty-four small-sized soil samples were subjected to two consecutive unsaturated steady-state flow leaching experiments with bromide as tracer. Observed breakthrough curves (BTCs) were analysed by the method of moments and by the advection–dispersion equation (ADE) to classify solute behaviour. Frequency distributions of the parameters indicating the solute velocity were heavily skewed or bimodal, reflecting the broad variability of the leaching scenarios, including some with pronounced preferential solute breakthrough. Exclusion of the preferential flow columns from our calculations revealed an average amount of 37% of immobile water. The large-scale BTCs derived from assembling the individual concentration courses of each run showed similar features, such as an early bromide breakthrough. However, two distinct apices, viz. one preferential and one matrix, were observed only in the first run, whereas the concentration decrease between the peaks was missing from the second run. A change in soil structure with continuous leaching was presumed to modify the interplay of the various flow domains, thereby altering the spreading of the BTCs. Correlation analysis between parameters of both tests suggests that preferential transport conditions are likely to occur at the same locations in the field over several infiltration cycles, whereas the ‘classical’ or expected matrix flow is time variant and therefore seems to be hardly predictable. © 1998 John Wiley & Sons, Ltd.     

Positive solution of two-dimensional solute transport in heterogeneous aquifers

The transport of contaminants in aquifers is usually represented by a convection-dispersion equation. There are several well-known problems of oscillation and artificial dispersion that affect the numerical solution of this equation. For example, several studies have shown that standard treatment of the cross-dispersion terms always leads to a negative concentration. It is also well known that the numerical solution of the convective term is affected by spurious oscillations or substantial numerical dispersion. These difficulties are especially significant for solute transport in nonuniform flow in heterogeneous aquifers. For the case of coupled reactive-transport models, even small negative concentration values can become amplified through nonlinear reaction source/sink terms and thus result in physically erroneous and unstable results. This paper includes a brief discussion about how nonpositive concentrations arise from numerical solution of the convection and cross-dispersion terms. We demonstrate the effectiveness of directional splitting with one-dimensional flux limiters for the convection term. Also, a new numerical scheme for the dispersion term that preserves positivity is presented. The results of the proposed convection scheme and the solution given by the new method to compute dispersion are compared with standard numerical methods as used in MT3DMS.     

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.     

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.     

Using PHREEQC to simulate solute transport in fractured bedrock

The geochemical computer model PHREEQC can simulate solute transport in fractured bedrock aquifers that can be conceptualized as dual-porosity flow systems subject to one-dimensional advective-dispersive transport in the bedrock fractures and diffusive transport in the bedrock matrix. This article demonstrates how the physical characteristics of such flow systems can be parameterized for use in PHREEQC, it provides a method for minimizing numerical dispersion in PHREEQC simulations, and it compares PHREEQC simulations with results of an analytical solution. The simulations assumed a dual-porosity conceptual model involving advective-reactive-dispersive transport in the mobile zone (bedrock fracture) and diffusive-reactive transport in the immobile zone (bedrock matrix). The results from the PHREEQC dual-porosity transport model that uses a finite-difference approach showed excellent agreement compared with an analytical solution.     

Approximate solution for solute transport during spherical-flow push-pull tests

An approximate analytical solution to the advection-dispersion equation was derived to describe solute transport during spherical-flow conditions in single-well push-pull tests. The spherical-flow case may be applicable to aquifer tests conducted in packed intervals or partially penetrating wells. Using results of two-dimensional numerical simulations, we briefly illustrate the applicability of the derived spherical-flow solution and provide a comparison with its cylindrical-flow counterpart. Good agreement between simulated extraction-phase breakthrough curves and the spherical-flow solution was found when the length of the injection/extraction region was small compared to both aquifer thickness and maximum solute frontal position at the end of the injection phase. On the other hand, discrepancies between simulated breakthrough curves and the spherical-flow solution increased with increasing anisotropy in hydraulic conductivities. Several inherent limitations embedded in its derivation such as assumptions of isotropy and homogeneity warrant the cautious use of the spherical-flow solution.     

Laplace transform solutions for solute transport in fissured aquifers

The main processes affecting the migration of a solute in a fissured aquifer will be advection and dispersion in the fissures, diffusion into the porous matrix; and adsorption. This paper considers solute transport in an idealized fissured aquifer consisting of slabs of saturated rock-matrix separated by equally spaced, planar fissures. The solution of the transport equations is developed as far as Laplace transforms of the solute concentrations in the fissure and matrix water. Numerical inversion of the transforms is used to investigate characteristic behaviour of the model for a number of special cases.     

Continuous time random walks for non-local radial solute transport

This study formulates and analyzes continuous time random walk (CTRW) models in radial flow geometries for the quantification of non-local solute transport induced by heterogeneous flow distributions and by mobile–immobile mass transfer processes. To this end we derive a general CTRW framework in radial coordinates starting from the random walk equations for radial particle positions and times. The particle density, or solute concentration is governed by a non-local radial advection–dispersion equation (ADE). Unlike in CTRWs for uniform flow scenarios, particle transition times here depend on the radial particle position, which renders the CTRW non-stationary. As a consequence, the memory kernel characterizing the non-local ADE, is radially dependent. Based on this general formulation, we derive radial CTRW implementations that (i) emulate non-local radial transport due to heterogeneous advection, (ii) model multirate mass transfer (MRMT) between mobile and immobile continua, and (iii) quantify both heterogeneous advection in a mobile region and mass transfer between mobile and immobile regions. The expected solute breakthrough behavior is studied using numerical random walk particle tracking simulations. This behavior is analyzed by explicit analytical expressions for the asymptotic solute breakthrough curves. We observe clear power-law tails of the solute breakthrough for broad (power-law) distributions of particle transit times (heterogeneous advection) and particle trapping times (MRMT model). The combined model displays two distinct time regimes. An intermediate regime, in which the solute breakthrough is dominated by the particle transit times in the mobile zones, and a late time regime that is governed by the distribution of particle trapping times in immobile zones. These radial CTRW formulations allow for the identification of heterogeneous advection and mobile-immobile processes as drivers of anomalous transport, under conditions relevant for field tracer tests.     

A stochastic interpretation of the tailing effect in solute transport

 The advection-dispersion equation (ADE) is inadequate for describing tails in solute breakthrough curves. Re-examination of solute breakthrough curves from one-dimensional experiments in porous media and channel flow literature shows a consistent discrepancy compared with solutions to the ADE. The leading tail of breakthrough curves is sharper, and the trailing tail is longer and smoother, than best fitting, least-squares ADE solutions. A random particle simulation exercise shows that the ADE may firstly be erroneous because of the assumption of time steps over which random solute movements are considered independent. Definition of such time steps hinges upon the slowest random movements, such as those predominantly by molecular diffusion. A second potential source of error is the highly skewed nature of the inverse distribution of underlying, micro-scale velocities, which causes slow convergence to normality under the central limit theorem.