Characterization of mixing and spreading in a bounded stratified medium

Matheron and de Marsily [Matheron M, de Marsily G. Is the transport in porous media always diffusive? A counter-example. Water Resour Res 1980;16:901–17] studied transport in a perfectly stratified infinite medium as an idealized aquifer model. They observed superdiffusive solute spreading quantified by anomalous increase of the apparent longitudinal dispersion coefficient with the square root of time. Here, we investigate solute transport in a vertically bounded stratified random medium. Unlike for the infinite medium at asymptotically long times, disorder-induced mixing and spreading is uniquely quantified by a constant Taylor dispersion coefficient. Using a stochastic modeling approach we study the effective mixing and spreading dynamics at pre-asymptotic times in terms of effective average transport coefficients. The latter are defined on the basis of local moments, i.e., moments of the transport Green function. We investigate the impact of the position of the initial plume and the initial plume size on the (highly anomalous) pre-asymptotic effective spreading and mixing dynamics for single realizations and in average. Effectively, the system “remembers” its initial state, the effective transport coefficients show so-called memory effects, which disappear after the solute has sampled the full vertical extent of the medium. We study the impact of the intrinsic non-ergodicity of the confined medium on the validity of the stochastic modeling approach and study in this context the transition from the finite to the infinite medium.     

On boundary conditions in the lattice Boltzmann model for advection and anisotropic dispersion equation

A mirror-image method is proposed in this paper to solve the boundary conditions in the lattice Boltzmann model proposed by Zhang et al. [Adv. Water Resour. 25 (2002) 1] for the advection and anisotropic dispersion of solute transport in porous media. Three types of boundary are considered: prescribed concentration boundary, prescribed flux boundary and prescribed concentration-gradient boundary. The accuracy of the proposed method is verified against benchmark problems and finite difference method.     

Efficient flow and transport simulations in reconstructed 3D pore geometries

Upscaling pore-scale processes into macroscopic quantities such as hydrodynamic dispersion is still not a straightforward matter for porous media with complex pore space geometries. Recently it has become possible to obtain very realistic 3D geometries for the pore system of real rocks using either numerical reconstruction or micro-CT measurements. In this work, we present a finite element–finite volume simulation method for modeling single-phase fluid flow and solute transport in experimentally obtained 3D pore geometries. Algebraic multigrid techniques and parallelization allow us to solve the Stokes and advection–diffusion equations on large meshes with several millions of elements. We apply this method in a proof-of-concept study of a digitized Fontainebleau sandstone sample. We use the calculated velocity to simulate pore-scale solute transport and diffusion. From this, we are able to calculate the a priori emergent macroscopic hydrodynamic dispersion coefficient of the porous medium for a given molecular diffusion D_m of the solute species. By performing this calculation at a range of flow rates, we can correctly predict all of the observed flow regimes from diffusion dominated to convection dominated.     

Higher and lowest order mixed finite element approximation of subsurface flow problems with solutions of low regularity

In this work we study mixed finite element approximations of Richards’ equation for simulating variably saturated subsurface flow and simultaneous reactive solute transport. Whereas higher order schemes have proved their ability to approximate reliably reactive solute transport (cf., e.g. [Bause M, Knabner P. Numerical simulation of contaminant biodegradation by higher order methods and adaptive time stepping. Comput Visual Sci 7;2004:61–78]), the Raviart–Thomas mixed finite element method (RT₀) with a first order accurate flux approximation is popular for computing the underlying water flow field (cf. [Bause M, Knabner P. Computation of variably saturated subsurface flow by adaptive mixed hybrid finite element methods. Adv Water Resour 27;2004:565–581, Farthing MW, Kees CE, Miller CT. Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow. Adv Water Resour 26;2003:373–394, Starke G. Least-squares mixed finite element solution of variably saturated subsurface flow problems. SIAM J Sci Comput 21;2000:1869–1885, Younes A, Mosé R, Ackerer P, Chavent G. A new formulation of the mixed finite element method for solving elliptic and parabolic PDE with triangular elements. J Comp Phys 149;1999:148–167, Woodward CS, Dawson CN. Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media. SIAM J Numer Anal 37;2000:701–724]). This combination might be non-optimal. Higher order techniques could increase the accuracy of the flow field calculation and thereby improve the prediction of the solute transport. Here, we analyse the application of the Brezzi-Douglas-Marini element (BDM₁) with a second order accurate flux approximation to elliptic, parabolic and degenerate problems whose solutions lack the regularity that is assumed in optimal order error analyses. For the flow field calculation a superiority of the BDM₁ approach to the RT₀ one is observed, which however is less significant for the accompanying solute transport.     

New equations for binary gas transport in porous media,Part 2: experimental validation

A previous study [Water Resour Res 39 (3) (2003) doi:10.1029/2002WR001338] questioned the validity of the traditional advection–dispersion equation for describing gas flow in porous media. In an original mathematical derivation presented in Part 1 [Adv Water Resour, this issue] we have demonstrated the theoretical existence of two novel physical phenomena which govern the macroscopic transport of gases in porous media. In this work we utilize laboratory experiments and numerical modeling in order to ascertain the importance of these novel theoretical terms. Numerical modeling results indicate that the newly derived sorptive velocity, arising from closure level coupling effects, does not contribute noticeably to the overall flux, under the conditions explored in this work. We demonstrate that the newly discovered “slip coupling” phenomenon in the mass conservation equation plays an important role in describing the physics of gas flow through porous solids for flow regimes of both environmental and industrial interest.     

Application of a mixing-ratios based formulation to model mixing-driven dissolution experiments

We address the question of how one can combine theoretical and numerical modeling approaches with limited measurements from laboratory flow cell experiments to realistically quantify salient features of complex mixing-driven multicomponent reactive transport problems in porous media. Flow cells are commonly used to examine processes affecting reactive transport through porous media, under controlled conditions. An advantage of flow cells is their suitability for relatively fast and reliable experiments, although measuring spatial distributions of a state variable within the cell is often difficult. In general, fluid is sampled only at the flow cell outlet, and concentration measurements are usually interpreted in terms of integrated reaction rates. In reactive transport problems, however, the spatial distribution of the reaction rates within the cell might be more important than the bulk integrated value. Recent advances in theoretical and numerical modeling of complex reactive transport problems [De Simoni M, Carrera J, Sanchez-Vila X, Guadagnini A. A procedure for the solution of multicomponent reactive transport problems. Water Resour Res 2005;41:W11410. doi: 10.1029/2005WR004056, De Simoni M, Sanchez-Vila X, Carrera J, Saaltink MW. A mixing ratios-based formulation for multicomponent reactive transport. Water Resour Res 2007;43:W07419. doi: 10.1029/2006WR005256] result in a methodology conducive to a simple exact expression for the space–time distribution of reaction rates in the presence of homogeneous or heterogeneous reactions in chemical equilibrium. The key points of the methodology are that a general reactive transport problem, involving a relatively high number of chemical species, can be formulated in terms of a set of decoupled partial differential equations, and the amount of reactants evolving into products depends on the rate at which solutions mix. The main objective of the current study is to show how this methodology can be used in conjunction with laboratory experiments to properly describe the key processes that occur in a complex, geochemically-active system under chemical equilibrium conditions. We model three CaCO₃ dissolution experiments reported in Singurindy et al. [Singurindy O, Berkowitz B, Lowell RP. Carbonate dissolution and precipitation in coastal environments: Laboratory analysis and theoretical consideration. Water Resour Res 2004;40:W04401. doi: 10.1029/2003WR002651, Singurindy O, Berkowitz B, Lowell RP. Correction to Carbonate dissolution and precipitation in coastal environments: laboratory analysis and theoretical consideration. Water Resour Res 2005;41:W11701. doi: 10.1029/2005WR004433], in which saltwater and freshwater were mixed in different proportions. The integrated reaction rate within the cell estimated from the experiments are modeled independently by means of (a) a state-of-the-art reactive transport code, and (b) the uncoupled methodology of [12, 13], both of which use dispersivity as a single, adjustable parameter. The good agreement between the results from both methodologies demonstrates the feasibility of using simple solutions to design and analyze laboratory experiments involving complex geochemical problems.     

Stochastic analysis of bounded unsaturated flow in heterogeneous aquifers: Spectral/perturbation approach

This paper describes a stochastic analysis of steady state flow in a bounded, partially saturated heterogeneous porous medium subject to distributed infiltration. The presence of boundary conditions leads to non-uniformity in the mean unsaturated flow, which in turn causes non-stationarity in the statistics of velocity fields. Motivated by this, our aim is to investigate the impact of boundary conditions on the behavior of field-scale unsaturated flow. Within the framework of spectral theory based on Fourier–Stieltjes representations for the perturbed quantities, the general expressions for the pressure head variance, variance of log unsaturated hydraulic conductivity and variance of the specific discharge are presented in the wave number domain. Closed-form expressions are developed for the simplified case of statistical isotropy of the log hydraulic conductivity field with a constant soil pore-size distribution parameter. These expressions allow us to investigate the impact of the boundary conditions, namely the vertical infiltration from the soil surface and a prescribed pressure head at a certain depth below the soil surface. It is found that the boundary conditions are critical in predicting uncertainty in bounded unsaturated flow. Our analytical expression for the pressure head variance in a one-dimensional, heterogeneous flow domain, developed using a nonstationary spectral representation approach [Li S-G, McLaughlin D. A nonstationary spectral method for solving stochastic groundwater problems: unconditional analysis. Water Resour Res 1991;27(7):1589–605; Li S-G, McLaughlin D. Using the nonstationary spectral method to analyze flow through heterogeneous trending media. Water Resour Res 1995; 31(3):541–51], is precisely equivalent to the published result of Lu et al. [Lu Z, Zhang D. Analytical solutions to steady state unsaturated flow in layered, randomly heterogeneous soils via Kirchhoff transformation. Adv Water Resour 2004;27:775–84].     

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.     

Modeling non-equilibrium mass transport in biologically reactive porous media

We develop a one-equation non-equilibrium model to describe the Darcy-scale transport of a solute undergoing biodegradation in porous media. Most of the mathematical models that describe the macroscale transport in such systems have been developed intuitively on the basis of simple conceptual schemes. There are two problems with such a heuristic analysis. First, it is unclear how much information these models are able to capture; that is, it is not clear what the model's domain of validity is. Second, there is no obvious connection between the macroscale effective parameters and the microscopic processes and parameters. As an alternative, a number of upscaling techniques have been developed to derive the appropriate macroscale equations that are used to describe mass transport and reactions in multiphase media. These approaches have been adapted to the problem of biodegradation in porous media with biofilms, but most of the work has focused on systems that are restricted to small concentration gradients at the microscale. This assumption, referred to as the local mass equilibrium approximation, generally has constraints that are overly restrictive. In this article, we devise a model that does not require the assumption of local mass equilibrium to be valid. In this approach, one instead requires only that, at sufficiently long times, anomalous behaviors of the third and higher spatial moments can be neglected; this, in turn, implies that the macroscopic model is well represented by a convection–dispersion–reaction type equation. This strategy is very much in the spirit of the developments for Taylor dispersion presented by Aris (1956). On the basis of our numerical results, we carefully describe the domain of validity of the model and show that the time-asymptotic constraint may be adhered to even for systems that are not at local mass equilibrium.     

Semi-analytical solutions for solute transport in fractured porous media using a strip source of finite width

Transient and steady-state analytical solutions are derived to investigate solute transport in a fractured porous medium consisting of evenly spaced, parallel discrete fractures. The solutions incorporate a finite width strip source, longitudinal and transverse dispersion in the fractures, source decay, aqueous phase decay, one-dimensional diffusion into the matrix, sorption to fracture walls, and sorption within the matrix. The solutions are derived using Laplace and Fourier transforms, and inverted by interchanging the order of integration and utilizing a numerical Laplace inversion algorithm. The solutions are verified for simplified cases by comparison to solutions derived by Batu [Batu V. A generalized two-dimensional analytical solution for hydrodynamic dispersion in bounded media with the first-type condition at the source. Wat Resour Res 1989;25(6):1125] and Sudicky and Frind [Sudicky EA, Frind EO. Contaminant transport in fractured porous media: analytical solutions for a system of parallel fractures. Wat Resour Res 1982;18(6):1634]. The application of the solutions to a fractured sandstone demonstrates that narrower source widths and larger values of transverse dispersivity both lead to lower downstream concentrations in the fractures and shorter steady-state plumes. The incorporation of aqueous phase decay and source concentration decay both lead to lower concentrations and shorter plumes, with even moderate amounts of decay significantly shortening the persistence of contamination.     

Conditional simulations of water–oil flow in heterogeneous porous media

This study is an extension of the stochastic analysis of transient two-phase flow in randomly heterogeneous porous media (Chen et al. in Water Resour Res 42:W03425, 2006), by incorporating direct measurements of the random soil properties. The log-transformed intrinsic permeability, soil pore size distribution parameter, and van Genuchten fitting parameter are treated as stochastic variables that are normally distributed with a separable exponential covariance model. These three random variables conditioned on given measurements are decomposed via Karhunen–Loève decomposition. Combined with the conditional eigenvalues and eigenfunctions of random variables, we conduct a series of numerical simulations using stochastic transient water–oil flow model (Chen et al. in Water Resour Res 42:W03425, 2006) based on the KLME approach to investigate how the number and location of measurement points, different random soil properties, as well as the correlation length of the random soil properties, affect the stochastic behavior of water and oil flow in heterogeneous porous media.     

Analytical solutions to the transient,unsaturated transport of water and contaminants through horizontal porous media

We present a range of analytical solutions to the combined transient water and solute transport for horizontal flow. We adopt the concept of a scale and time dependent dispersivity used for contaminant transport in aquifers and apply it to transient, unsaturated horizontal flow to develop similarity solutions for both constant solute concentration and solute flux boundary conditions. Through the use of a specific form of the water profile as used by Brutsaert [Water Resour Res 1968:4;785], the solute profiles can be reduced to a simple quadrature. We also derive a solution for the instantaneous injection of water and solute into a horizontal media for an arbitrary dispersivity. It is found that the solute concentration remains constant in both space and time as the water redistributes, suggesting that the solute does not disperse relative to the water.     

Chaotic advection at the pore scale: Mechanisms,upscaling and implications for macroscopic transport

The macroscopic spreading and mixing of solute plumes in saturated porous media is ultimately controlled by processes operating at the pore scale. Whilst the conventional picture of pore-scale mechanical dispersion and molecular diffusion leading to persistent hydrodynamic dispersion is well accepted, this paradigm is inherently two-dimensional (2D) in nature and neglects important three-dimensional (3D) phenomena. We discuss how the kinematics of steady 3D flow at the pore scale generate chaotic advection—involving exponential stretching and folding of fluid elements—the mechanisms by which it arises and implications of microscopic chaos for macroscopic dispersion and mixing. Prohibited in steady 2D flow due to topological constraints, these phenomena are ubiquitous due to the topological complexity inherent to all 3D porous media. Consequently 3D porous media flows generate profoundly different fluid deformation and mixing processes to those of 2D flow. The interplay of chaotic advection and broad transit time distributions can be incorporated into a continuous-time random walk (CTRW) framework to predict macroscopic solute mixing and spreading. We show how these results may be generalised to real porous architectures via a CTRW model of fluid deformation, leading to stochastic models of macroscopic dispersion and mixing which both honour the pore-scale kinematics and are directly conditioned on the pore-scale architecture.     

Transport in chemically and mechanically heterogeneous porous media. I: Theoretical development of region-averaged equations for slightly compressible single-phase flow

Transport processes in heterogeneous porous media are often treated in terms of one-equation models. Such treatment assumes that the velocity, pressure, temperature, and concentration can be represented in terms of a single large-scale averaged quantity in regions having significantly different mechanical, thermal, and chemical properties. In this paper we explore the process of single-phase flow in a two-region model of heterogeneous porous media. The region-averaged equations are developed for the case of a slightly compressible flow which is an accurate representation for a certain class of liquid-phase flows. The analysis leads to a pair of transport equations for the region averaged pressures that are coupled through a classic exchange term, in addition to being coupled by a diffusive cross effect. The domain of validity of the theory has been identified in terms of a series of length and timescale constraints.In Part II the theory is tested, in the absence of adjustable parameters, by comparison with numerical experiments for transient, slightly compressible flow in both stratified and nodular models of heterogeneous porous media. Good agreement between theory and experiment is obtained for nodular and stratified systems, and effective transport coefficients for a wide range of conditions are presented on the basis of solutions of the three closure problems that appear in the theory. Part III of this paper deals with the principle of large-scale mechanical equilibrium and the region-averaged form of Darcy's law. This form is necessary for the development and solution of the region-averaged solute transport equations that are presented in Part IV. Finally, in Part V we present results for the dispersion tensors and the exchange coefficient associated with the two-region model of solute transport with adsorption.     

Exploring the nature of non-Fickian transport in laboratory experiments

Observations of non-Fickian transport in sandbox experiments [Levy M, Berkowitz B. Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. J Contam Hydrol 2003;64:203–26] were analyzed previously using a power law tail ψ(t) ∼ t^−1−β with 0 < β < 2 for the spectrum of transition times comprising a tracer plume migration. For each sandbox medium a choice of β resulted in an excellent fit to the breakthrough curve (BTC) data, and the value of β decreased slowly with increasing flow velocity. Here, the data are reanalyzed with the full spectrum of ψ(t) gleaned from analytical calculations [Cortis A, Chen Y, Scher H, Berkowitz B. Quantitative characterization of pore-scale disorder effects on transport in “homogeneous” granular media. Phys Rev E 2004;10(70):041108. doi: 10.1103/PhysRevE.70.041108], numerical simulations [Bijeljic B, Blunt MJ. Pore-scale modeling and continuous time random walk analysis of dispersion in porous media. Water Resour Res 2006;42:W01202. doi: 10.1029/2005WR004578] and permeability fields [Di Donato G, Obi E-O, Blunt MJ. Anomalous transport in heterogeneous media demonstrated by streamline-based simulation. Geophys Res Lett 2003;30:1608–12s. doi: 10.1029/2003GL017196]. We represent the main features of the full spectrum of transition times with a truncated power law (TPL), ψ(t) ∼ (t₁ + t)^−1−βexp(−t/t₂), where t₁ and t₂ are the limits of the power law spectrum. An excellent fit to the entire BTC data set, including the changes in flow velocity, for each sandbox medium is obtained with a single set of values of t₁, β, t₂. The influence of the cutoff time t₂ is apparent even in the regime t < t₂. Significantly, we demonstrate that the previous apparent velocity dependence of β is a result of choosing a pure power law tail for ψ(t). The key is the change in the log–log slope of the TPL form of ψ(t) with a shifting observational time window caused by the change in the mean velocity. Hence, the use of the full spectrum of ψ(t) is not only necessary for the transition to Fickian behavior, but also to account for the dynamics of these laboratory observations of non-Fickian transport.     

An ELLAM approximation for advective–dispersive transport with nonlinear sorption

We consider an Eulerian–Lagrangian localized adjoint method (ELLAM) applied to nonlinear model equations governing solute transport and sorption in porous media. Solute transport in the aqueous phase is modeled by standard advection and hydrodynamic dispersion processes, while sorption is modeled with a nonlinear local-equilibrium model. We present our implementation of finite volume ELLAM (FV-ELLAM) and finite element (FE-ELLAM) discretizations to the reactive transport model and evaluate their performance for several test problems containing self-sharpening fronts.     

Analytical solutions of one-dimensional macrodispersion in stratified porous media by the quadrupole method: convergence to an equivalent homogeneous porous medium

In this article, the quadrupole method is implemented in order to simulate the effects of heterogeneities on one dimensional advective and diffusive transport of a passive solute in porous media. Theoretical studies of dispersion in heterogeneous stratified media can bring insight into transport artefacts linked to scale effects and apparent dispersion coefficients. The quadrupole method is an efficient method for the calculation of transient response of linear systems. It is based here on the Laplace transform technique. The analytical solutions that can be derived by this method assists understanding of upscaled parameters relevant to heterogeneous porous media.First, the method is developed for an infinite homogeneous porous medium. Then, it is adapted to a stratified medium where the fluid flow is perpendicular to the interfaces. The first heterogeneous medium studied is composed of two semi-infinite layers perpendicular to the flow direction each having different transport properties. The concentration response of the medium to a Dirac injection is evaluated. The case studied emphasises the importance in the choice of the boundary conditions.In the case of a periodic heterogeneous porous medium, the concentration response of the medium is evaluated for different numbers of unit-cells. When the number of unit cells is great enough, depending on the transport properties of each layer in the unit cell, an equivalent homogeneous behaviour is reached. An exact determination of the transport properties (equivalent dispersion coefficient) of the equivalent homogeneous porous medium is given.     

Experimental and numerical study of the validity of Hele–Shaw cell as analogue model for variable-density flow in homogeneous porous media

Several laboratory experiments were conducted to identify the validity domain under which a Hele–Shaw cell may serve as a suitable analogue for variable-density flow in homogeneous porous media. These experiments are concerned with the injection into a Hele–Shaw cell of a salt solute at different concentrations and flow rates. The experimental data analysis highlighted two types of mixing zone shape: with and without ‘fingers’. A semi-empirical criterion based on the ratio between gravitational and injected velocities was used to forecast the change from one shape to another. The experimental data were then analysed using numerical solutions of the classical Hele–Shaw equations by taking into account an anisotropic dispersion tensor whose components depend on fluid density gradients. The good agreement between experimental and numerical results clearly shows that the validity of the concentration-dependent dispersion tensor strongly depends on the local Péclet number variation. For Péclet numbers lower 50, the Hele–Shaw cell can be considered as an analogous model of a homogeneous and isotropic 2D porous medium. It can be successfully used to study, at the laboratory scale, the gravitational instability effects induced by flow and transport phenomena into a porous medium.     

Numerical modelling of mass transport in hydrogeologic environments: performance comparison of the Laplace Transform Galerkin and Arnoldi modal reduction schemes

The Laplace Transform Galerkin (LTG) method and the Arnoldi modal reduction method (AMRM) have been implemented in finite element schemes designed to solve mass transport problems in porous media by Sudicky [Sudicky, E.A., Water Resour. Res., 25(8) (1989) 1833-46] and Woodbury et al [Woodbury, A.D., Dunbar, W.S., & Nour-Omid, B., Water Resour. Res., 26(10) (1990) 2579-90]. In this work, a comparative analysis of the two methods is performed with attention focused on efficiency and accuracy. The analysis is performed over one- and two-dimensional domains composed of homogeneous and heterogeneous material properties. The results obtained using homogeneous material properties indicate that for a given mesh design the LTG method maintains a higher degree of accuracy than does the AMRM. However, in terms of efficiency, the Arnoldi attains a pre-defined level of accuracy faster than does the LTG method. It is also shown that for problems involving homogeneous material properties the solution obtained using the LTG method on a coarse mesh is comparable in terms of solution time and accuracy to that obtained using the AMRM on a fine mesh. Comparisons similar to those performed using homogeneous material properties are also performed for the case where the hydraulic conductivity field is heterogeneous. For this case, the level of accuracy achieved by the AMRM and the LTG method are similar. However, as with the analysis involving homogeneous material properties, the AMRM is found to be more efficient than the LTG method. It is also shown that for heterogeneous material properties, use of the LTG method under high grid Peclet conditions can be potentially problematic.