A simple and efficient random walk solution of multi-rate mobile/immobile mass transport equations

We extend the particle-tracking method to simulate general multi-rate mass transfer (MRMT) equations. Previous methods for single-rate equations used two-state Markov chains and found that the time a particle spends in the mobile state between waiting time epochs is random and exponentially distributed. Using Bochner’s subordination technique for Markov processes, we find that the random mobile times are still exponential for the stochastic process that corresponds to the MRMT equations. The random times in the immobile phase have a distribution that is directly related to the memory function of the MRMT equation. This connection allows us to interpret the MRMT memory function as the rate at which particles of a certain age, measured by residence time in the immobile zone, exit to become mobile once again. Because the exact distributions of mobile and immobile times are known from the MRMT equations, they can be simulated very simply and efficiently using random walks.     

Experiment and Simulation of Non-Reactive Solute Transport in Porous Media

To more accurately predict the migration behavior of pollutants in porous media, we conduct laboratory scale experiments and model simulation. Aniline (AN) is used in one-dimensional soil column experiments designed under various media and hydrodynamic conditions. The advection-dispersion equation (ADE) and the continuous-time random walk (CTRW) were used to simulate the breakthrough curves (BTCs) of the solute transport. The results show that the media and hydrodynamic conditions are two important factors affecting solute transport and are related to the degree of non-Fickian transport. The simulation results show that CTRW can more effectively describe the non-Fickian phenomenon in the solute transport process than ADE. The sensitive parameter in the CTRW simulation process is , which can reflect the degree of non-Fickian diffusion in the solute transport. Understanding the relationship of with velocity and media particle size is conducive to improving the reactive solute transport model. The results of this study provide a theoretical basis for better prediction of pollutant transport in groundwater.     

Local Equilibrium and Retardation Revisited

In modeling solute transport with mobile‐immobile mass transfer (MIMT), it is common to use an advection‐dispersion equation (ADE) with a retardation factor, or retarded ADE. This is commonly referred to as making the local equilibrium assumption (LEA). Assuming local equilibrium, Eulerian textbook treatments derive the retarded ADE, ostensibly exactly. However, other authors have presented rigorous mathematical derivations of the dispersive effect of MIMT, applicable even in the case of arbitrarily fast mass transfer. We resolve the apparent contradiction between these seemingly exact derivations by adopting a Lagrangian point of view. We show that local equilibrium constrains the expected time immobile, whereas the retarded ADE actually embeds a stronger, nonphysical, constraint: that all particles spend the same amount of every time increment immobile. Eulerian derivations of the retarded ADE thus silently commit the gambler's fallacy, leading them to ignore dispersion due to mass transfer that is correctly modeled by other approaches. We then present a particle tracking simulation illustrating how poor an approximation the retarded ADE may be, even when mobile and immobile plumes are continually near local equilibrium. We note that classic “LEA” (actually, retarded ADE validity) criteria test for insignificance of MIMT‐driven dispersion relative to hydrodynamic dispersion, rather than for local equilibrium.     

Perspective on theories of non-Fickian transport in heterogeneous media

Subsurface fluid flow and solute transport take place in a multiscale heterogeneous environment. Neither these phenomena nor their host environment can be observed or described with certainty at all scales and locations of relevance. The resulting ambiguity has led to alternative conceptualizations of flow and transport and multiple ways of addressing their scale and space–time dependencies. We focus our attention on four approaches that give rise to nonlocal representations of advective and dispersive transport of nonreactive tracers in randomly heterogeneous porous or fractured continua. We compare these approaches theoretically on the basis of their underlying premises and the mathematical forms of the corresponding nonlocal advective–dispersive terms. One of the four approaches describes transport at some reference support scale by a classical (Fickian) advection–dispersion equation (ADE) in which velocity is a spatially (and possibly temporally) correlated random field. The randomness of the velocity, which is given by Darcy’s law, stems from random fluctuations in hydraulic conductivity (and advective porosity though this is often disregarded). Averaging the stochastic ADE over an ensemble of velocity fields results in a space–time-nonlocal representation of mean advective–dispersive flux, an approach we designate as stnADE. A closely related space–time-nonlocal representation of ensemble mean transport is obtained upon averaging the motion of solute particles through a random velocity field within a Lagrangian framework, an approach we designate stnL. The concept of continuous time random walk (CTRW) yields a representation of advective–dispersive flux that is nonlocal in time but local in space. Closely related to the latter are forms of ADE entailing fractional derivatives (fADE) which leads to representations of advective–dispersive flux that are nonlocal in space but local in time; nonlocality in time arises in the context of multirate mass transfer models, which we exclude from consideration in this paper. We describe briefly each of these four nonlocal approaches and offer a perspective on their differences, commonalities, and relative merits as analytical and predictive tools.     

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.     

A Practical Modeling Framework for Non‐Fickian Transport and Multi‐Species Sequential First‐Order Reaction

Many studies indicate that small‐scale heterogeneity and/or mobile–immobile mass exchange produce transient non‐Fickian plume behavior that is not well captured by the use of the standard, deterministic advection‐dispersion equation (ADE). An extended ADE modeling framework is presented here that is based on continuous time random walk theory. It can be used to characterize non‐Fickian transport coupled with simultaneous sequential first‐order reactions (e.g., biodegradation or radioactive decay) for multiple degrading contaminants such as chlorinated solvents, royal demolition explosive, pesticides, and radionuclides. To demonstrate this modeling framework, new transient analytical solutions are derived and are inverted in Laplace space. Closed‐form, steady‐state, multi‐species analytical solutions are also derived for non‐Fickian transport in highly heterogeneous aquifers with linear sorption–desorption and matrix diffusion for use in spreadsheets. The solutions are general enough to allow different degradation rates for the mobile and immobile zones. The transient solutions for multi‐species transport are applied to examine the effects of source remediation on the natural attenuation of downgradient plumes of both parent and degradation products in highly heterogeneous aquifers. Results for representative settings show that the use of the standard, deterministic ADE can over‐estimate cleanup rates and under‐predict the cleanup timeframe in comparison to the extended ADE analytical model. The modeling framework and calculations introduced here are also applied for a 30 year groundwater cleanup program at a site in Palm Bay, Florida. The simulated plume concentrations using the extended ADE exhibited agreement with observed long concentration tails of trichloroethene, cis 1,2 DCE, and VC that remained above cleanup goals.     

Interpretation and nonuniqueness of CTRW transition distributions: Insights from an alternative solute transport formulation

The continuous time random walk (CTRW) has both an elegant mathematical theory and a successful record at modeling solute transport in the subsurface. However, there are some interpretation ambiguities relating to the relationship between the discrete CTRW transition distributions and the underlying continuous movement of solute that have not been addressed in existing literature. These include the exact definition of “transition”, and the extent to which transition probability distributions are unique/quantifiable from data. Here, we present some theoretical results which address these uncertainties in systems with an advective bias. Simultaneously, we present an alternative, reduced parameter CTRW formulation for general advective transport in heterogeneous porous media, which models early- and late-time transport by use of random transition times between sparse, imaginary planes normal to flow. We show that even in the context of this reduced-parameter formulation there is nonuniqueness in the definitions of both transition lengths and waiting time distributions, and that neither may be uniquely determined from experimental data. For practical use of this formulation, we suggest Pareto transition time distributions, leading to a two-degree-of-freedom modeling approach. We then demonstrate the power of this approach in fitting two sets of existing experimental data. While the primary focus is the presentation of new results, the discussion is designed to be pedagogical and to provide a good entry point into practical modeling of solute transport with the CTRW.     

Lattice Boltzmann Models for Flow and Transport in Saturated Karst

Flow and transport simulation in karst aquifers remains a significant challenge for the ground water modeling community. Darcy's law–based models cannot simulate the inertial flows characteristic of many karst aquifers. Eddies in these flows can strongly affect solute transport. The simple two-region conduit/matrix paradigm is inadequate for many purposes because it considers only a capacitance rather than a physical domain. Relatively new lattice Boltzmann methods (LBMs) are capable of solving inertial flows and associated solute transport in geometrically complex domains involving karst conduits and heterogeneous matrix rock. LBMs for flow and transport in heterogeneous porous media, which are needed to make the models applicable to large-scale problems, are still under development. Here we explore aspects of these future LBMs, present simple examples illustrating some of the processes that can be simulated, and compare the results with available analytical solutions. Simulations are contrived to mimic simple capacitance-based two-region models involving conduit (mobile) and matrix (immobile) regions and are compared against the analytical solution. There is a high correlation between LBM simulations and the analytical solution for two different mobile region fractions. In more realistic conduit/matrix simulation, the breakthrough curve showed classic features and the two-region model fit slightly better than the advection-dispersion equation (ADE). An LBM-based anisotropic dispersion solver is applied to simulate breakthrough curves from a heterogeneous porous medium, which fit the ADE solution. Finally, breakthrough from a karst-like system consisting of a conduit with inertial regime flow in a heterogeneous aquifer is compared with the advection-dispersion and two-region analytical solutions.     

Multi-Rate Mass Transfer (MRMT) models for general diffusive porosity structures

We determine the relevance of Multi-Rate Mass Tansfer (MRMT) models (Haggerty and Gorelick, 1995) to general diffusive porosity structures. To this end, we introduce Structured INteracting Continua (SINC) models as the combination of a finite number of diffusion-dominated interconnected immobile zones exchanging with an advection-dominated mobile domain. It directly extends Multiple INteracting Continua framework (Pruess and Narasimhan, 1985) by introducing a structure in the immobile domain, coming for example from the dead-ends of fracture clusters or poorly-connected dissolution patterns. We demonstrate that, whatever their structure, SINC models can be made equivalent in terms of concentration in the mobile zone to a unique MRMT model. We develop effective shape-free numerical methods to identify its few dominant rates, that comply with any distribution of rates and porosities. We show that differences in terms of macrodispersion are not larger than 50% for approximate MRMT models with only one rate (double porosity models), and drop down to less than 0.1% for five rates MRMT models. Low-dimensional MRMT models accurately approach transport in structured diffusive porosities at intermediate and long times and only miss early responses.     

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.     

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.     

Predicting the tails of breakthrough curves in regional-scale alluvial systems

The late tail of the breakthrough curve (BTC) of a conservative tracer in a regional-scale alluvial system is explored using Monte Carlo simulations. The ensemble numerical BTC, for an instantaneous point source injected into the mobile domain, has a heavy late tail transforming from power law to exponential due to a maximum thickness of clayey material. Haggerty et al.'s (2000) multiple-rate mass transfer (MRMT) method is used to predict the numerical late-time BTCs for solutes in the mobile phase. We use a simple analysis of the thicknesses of fine-grained units noted in boring logs to construct the memory function that describes the slow decline of concentrations at very late time. The good fit between the predictions and the numerical results indicates that the late-time BTC can be approximated by a summation of a small number of exponential functions, and its shape depends primarily on the thicknesses and the associated volume fractions of immobile water in "blocks" of fine-grained material. The prediction of the late-time BTC using the MRMT method relies on an estimate of the average advective residence time, t(ad). The predictions are not sensitive to estimation errors in t(ad), which can be approximated by L/v , where v is the arithmetic mean ground water velocity and L is the transport distance. This is the first example of deriving an analytical MRMT model from measured hydrofacies properties to predict the late-time BTC. The parsimonious model directly and quantitatively relates the observable subsurface heterogeneity to nonlocal transport parameters.     

Applicability of the Dual‐Domain Model to Nonaggregated Porous Media

More theoretical analysis is needed to investigate why a dual‐domain model often works better than the classical advection‐dispersion (AD) model in reproducing observed breakthrough curves for relatively homogeneous porous media, which do not contain distinct dual domains. Pore‐scale numerical experiments presented here reveal that hydrodynamics create preferential flow paths that occupy a small part of the domain but where most of the flow takes place. This creates a flow‐dependent configuration, where the total domain consists of a mobile and an immobile domain. Mass transfer limitations may result in nonequilibrium, or significant differences in concentration, between the apparent mobile and immobile zones. When the advection timescale is smaller than the diffusion timescale, the dual‐domain mass transfer (DDMT) model better captures the tailing in the breakthrough curve. Moreover, the model parameters (mobile porosity, mean solute velocity, dispersivity, and mass transfer coefficient) demonstrate nonlinear dependency on mean fluid velocity. The studied case also shows that when the Peclet number, Pe, is large enough, the mobile porosity approaches a constant, and the mass transfer coefficient can be approximated as proportional to mean fluid velocity. Based on detailed analysis at the pore scale, this paper provides a physical explanation why these model parameters vary in certain ways with Pe. In addition, to improve prediction in practical applications, we recommend conducting experiments for parameterization of the DDMT model at a velocity close to that of the relevant field sites, or over a range of velocities that may allow a better parameterization.     

Influence of small-scale heterogeneities on contaminant transport in fractured crystalline rock

We present a sequence of purely advective transport models that demonstrate the influence of small-scale geometric inhomogeneities on contaminant transport in fractured crystalline rock. Special weight is placed on the role of statistically generated variable fracture apertures. The fracture network geometry and the aperture distribution are based on information from an in situ radionuclide retardation experiment performed at Grimsel test site (Swiss Alps). The obtained breakthrough curves are fitted with the advection dispersion equation and continuous-time random walks (CTRW). CTRW is found to provide superior fits to the late-arrival tailing and is also found to show a good correlation with the velocity distributions obtained from the hydraulic models. The impact of small-scale heterogeneities, both in fracture geometry and aperture, on transport is shown to be considerable.     

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     

Computational issues in the determination of solute discharge moments and implications for comparison to analytical solutions

Solute discharge moments (mean and variance) are computed using numerical modeling of flow and advective transport in two-dimensional heterogeneous aquifers and are compared to theoretical results. The solute discharge quantifies the temporal evolution of the total contaminant mass crossing a certain compliance boundary. In addition to analyzing the solute discharge moments within a classical absolute dispersion framework, we also analyze relative dispersion formulation, whereby plume meandering (deviation from mean flow path caused by velocity variations at scales larger than plume size) is removed. This study addresses some important issues related to the computation of solute discharge moments from random walk particle tracking experiments, and highlights some of the important differences between absolute and relative dispersion frameworks. Relative dispersion formulation produces maximum uncertainty that coincides with the peak mean discharge. Absolute dispersion, however, results in earlier arrival of the uncertainty peak as compared to the first moment peak. Simulations show that the standard deviation of solute discharge in a relative dispersion framework requires increasingly large temporal sampling windows to smooth out some of the large fluctuations in breakthrough curves associated with advective transport. Using smoothing techniques in particle tracking to distribute the particle mass over a volume rather than at a point significantly reduces the noise in the numerical simulations and removes the need to use large temporal windows. Same effect can be obtained by adding a local dispersion process to the particle tracking experiments used to model advective transport. The effect of the temporal sampling window bears some relevance and important consequences for evaluating risk-related parameters. The expected value of peak solute discharge and its standard deviation are very sensitive to this sampling window and so will be the risk distribution relying on such numerical models.     

Exact analytical solutions for two-dimensional advection-dispersion equation in cylindrical coordinates subject to third-type inlet boundary condition

Exact analytical solutions for two-dimensional advection-dispersion equation (ADE) in cylindrical coordinates subject to the third-type inlet boundary condition are presented in this study. The finite Hankel transform technique in combination with the Laplace transform method is adopted to solve the two-dimensional ADE in cylindrical coordinates. Solutions are derived for both continuous input and instantaneous slug input. The developed analytical solutions are compared with the solutions for first-type inlet boundary condition to illustrate the influence of the inlet condition on the two-dimensional solute transport in a porous medium system with a radial geometry. Results show significant discrepancies between the breakthrough curves obtained from analytical solutions for the first-type and third-type inlet boundary conditions for large longitudinal dispersion coefficients. The developed solutions conserve the solute mass and are efficient tools for simultaneous determination of the longitudinal and transverse dispersion coefficients from a laboratory-scale radial column experiment or an in situ infiltration test with a tracer.     

Gas tracer transport through a heterogeneous fracture zone under two phase flow conditions: Model development and parameter sensitivity

Large amounts of gas can result from anaerobic corrosion of metals and from chemical and biological degradation of organic substances in underground repositories for radioactive waste. Gas generation may lead to the formation of a gas phase bubble and to the migration of radioactive gaseous species. Transport occurs in, at least, in two forms: (1) gas bubble, migration is controlled by advection, dispersion and diffusion in the gas phase, and (2) within water pockets, the dissolved species migrate mainly by diffusion. We consider a two-dimensional system representing an isolated heterogeneous fractured zone. A dipole gas flow field is generated and gas tracers are injected. The delay in the breakthrough curves is studied. A simple method is used to solve the gas species transport equations in multiphase conditions. This method is based on a formal analogy between the equations of gas transport in a two phase system and the equations of solute tracer transport in water saturated systems. We perform a sensitivity analysis to quantify the relevance of the various transport mechanisms. We find that gas tracer migration is very sensitive to gas tracer solubility, which affects gas tracer transport of both mobile and immobile zones, and shows high sensitivity to diffusion in the gas phase, to heterogeneity and to gas pressure, but the largest sensitivity was observed with respect to injection borehole properties, i.e. borehole volume and water filled fraction.