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.     

Numerical simulations of non-ergodic solute transport in three-dimensional heterogeneous porous media

Numerical simulations of non-ergodic transport of a non-reactive solute plume by steady-state groundwater flow under a uniform mean velocity, , were conducted in a three-dimensional heterogeneous and statistically isotropic aquifer. The hydraulic conductivity, K(x), is modeled as a random field which is assumed to be log-normally distributed with an exponential covariance. Significant efforts are made to reduce the simulation uncertainties. Ensemble averages of the second spatial moments of the plume and the plume centroid variances were simulated with 1600 Monte Carlo (MC) runs for three variances of log K, _Y²=0.09, 0.23, and 0.46, and a square source normal to of three dimensionless lengths. It is showed that 1600 MC runs are needed to obtain stabilized results in mildly heterogeneous aquifers of _Y²0.5 and that large uncertainty may exist in the simulated results if less MC runs are used, especially for the transverse second spatial moments and the plume centroid variance in transverse directions. The simulated longitudinal second spatial moment and the plume centroid variance in longitudinal direction fit well to the first-order theoretical results while the simulated transverse moments are generally larger than the first-order values. The ergodic condition for the second spatial moments is far from reaching in all cases simulated and transport in transverse directions may reach ergodic condition much slower than that in longitudinal direction.     

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

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

Mixing cell method for solving the solute transport equation with spatially variable coefficients

The advection–dispersion equation with spatially variable coefficients does not have an exact analytical solution and is therefore solved numerically. However, solutions obtained with several of the traditional finite difference or finite element techniques typically exhibit spurious oscillation or numerical dispersion when advection is dominant. The mixing cell and semi-analytical solution methods proposed in this study avoid such oscillation or numerical dispersion when advection dominates. Both the mixing cell and semi-analytical solution methods calculate the spatial step size by equating numerical dispersion to physical dispersion. Because of the spatial variability of the coefficients the spatial step size varies in space. When the time step size Δt→0, the mixing cell method reduces to the semi-analytical solution method. The results of application to two cases show that the mixing cell and semi-analytical solution methods are better than a finite difference method used in the study. © 1998 John Wiley & Sons, Ltd.     

Microtopographical and hydrophysical controls on subsurface flow and solute transport: A continuous solute release experiment in a subarctic bog

Resource extraction and transportation activities in subarctic Canada can result in the unintentional release of contaminants into the surrounding peatlands. In the event of a release, a thorough understanding of solute transport within the saturated zone is necessary to predict plume fate and the potential impacts on peatland ecosystems. To better characterize contaminant transport in these systems, approximately 13,000 L/day of sodium chloride tracer (200 mg/L) was released into a bog in the James Bay Lowland. The tracer was pumped into a fully penetrating well (1.5 m) between July 5 and August 18, 2015. Horizontal and vertical plume development was measured via in situ specific conductance and water table depth from an adaptive monitoring network. Over the spill period, the bulk of the plume travelled a lateral distance of 100 m in the direction of the slight regional groundwater and topographical slope. The plume shape was irregular and followed the hollows, indicating preferential flow paths due to the site microtopography. Saturated transport of the tracer occurred primarily at ~25 cm below ground surface (bgs), and at a discontinuous high hydraulic conductivity layer ~125 cm bgs due to a complex and heterogeneous vertical hydraulic conductivity profile. Plume measurement was confounded by a large amount of precipitation (233 mm over the study period) that temporarily diluted the tracer in the highly conductive upper peat layer. Longitudinal solute advection can be approximated using local water table information (i.e., depth and gradient); microtopography; and meteorological conditions. Vertical distribution of solute within the peat profile is far more complex due to the heterogeneous subsurface; characterization would be aided by a detailed understanding of the site‐specific peat profile; the degree of decomposition; and the type of contaminant (e.g., reactive/nonreactive). The results of this research highlight the difficulty of tracking a contaminant spill in bogs and provide a benchmark for the characterization of the short‐term fate of a plume in these complex systems.     

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.     

A novel platform to evaluate the dampening of water and solute transport in an experimental channel under unsteady flow conditions

This paper presents a novel platform to study the dampening of water and solute transport in an experimental channel under unsteady flow conditions, where literature data are scarce. We address the question about what could be the smallest size of experimental platform that is useful for research, project studies, and teaching activities and that allows to do rational experiments characterized by small space occupation, short experimental duration, high measurement precision, high quality and reproducible experimental curves, low water and energy consumption, and the possibility to test a large variety of hydrograph scenarios. Whereas large scale hydraulic laboratories have focused their studies on sediment transport, our platform deals with solute transport. The objectives of our study are (a) building a platform that allows to do rational experiments, (b) enriching the lack of experimental data concerning water and solute transport under unsteady state conditions, and (c) studying the dampening of water and solute transport. We studied solute transport in a channel with lateral gain and lateral loss under different experimental configurations, and we show how the same lateral loss flow event can lead to different lateral loss mass repartitions under different configurations. In order to characterize water and solute dampening between the input and the output of the channel, we calculate dampening ratios based on peak coordinates of time flow curves and time mass curves and that express the decrease of peak amplitude and the increase of peak occurrence time between the input and output curves. Finally, we use a solute transport model coupling the diffusive wave equation for water transfer and the advection–diffusion equation for solute transport in order to simulate the experimental data. The simulations are quite good with a Nash–Sutcliffe efficiency NSE > 0.98 for water transfer and 0.84 < NSE < 0.97 for solute transport. This platform could serve hydrological modellers because it offers a variety of measured parameters (flow, water height, and solute concentration), at a fine time step under unsteady flow conditions.     

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

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

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.     

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.     

Calculating the macrodispersion coefficient of the ensemble averaged solute transport equation in the discrete domain

Recognizing the spatial heterogeneity of hydraulic parameters, many researchers have studied the solute transport by both groundwater and channel flow in a stochastic framework. One of the methodologies used to up‐scale the stochastic solute transport equation, from a point‐location scale to a grid scale, is the cumulant expansion method combined with the calculus for the time‐ordered exponential and the calculus for the Lie operator. When the point‐location scale transport equation is scaled up to the grid scale, using the cumulant expansion method, a new dispersion coefficient emerges in the dispersive term of the solute transport equation in addition to the molecular dispersion coefficient. This velocity driven dispersion is called ‘macrodispersion’. The macrodispersion coefficient is the integral function of the time‐ordered covariance of the random velocity field. The integral is calculated over a Lagrangian trajectory of the flow. The Lagrangian trajectory depends on the following: (i) the spatial origin of the particle; (ii) the time when the macrodispersion is calculated; and (iii) the mean velocity field along the trajectory itself. The Lagrangian trajectory is a recursive function of time because the location of the particle along the trajectory at a particular time depends on the location of the particle at the previous time. This recursive functional form of the Lagrangian trajectory makes the calculation of the macrodispersion coefficient difficult. Especially for the unsteady, spatially non‐stationary, non‐uniform flow field, the macrodispersion coefficient is a highly complex expression and, so far, calculated using numerical methods in the discrete domains. Here, an analytical method was introduced to calculate the macrodispersion coefficient in the discrete domain for the unsteady and steady, spatially non‐stationary flow cases accurately and efficiently. This study can fill the gap between the theory of the ensemble averaged solute transport model and its numerical implementations. Copyright © 2011 John Wiley & Sons, Ltd.     

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.     

Variable-density groundwater flow and solute transport in porous media containing nonuniform discrete fractures

Variations in fluid density can greatly affect fluid flow and solute transport in the subsurface. Heterogeneities such as fractures play a major role for the migration of variable-density fluids. Earlier modeling studies of density effects in fractured media were restricted to orthogonal fracture networks, consisting of only vertical and horizontal fractures. The present study addresses the phenomenon of 3D variable-density flow and transport in fractured porous media, where fractures of an arbitrary incline can occur. A general formulation of the body force vector is derived, which accounts for variable-density flow and transport in fractures of any orientation. Simulation results are presented that show the verification of the new model formulation, for the porous matrix and for inclined fractures. Simulations of variable-density flow and solute transport are then conducted for a single fracture, embedded in a porous matrix. The simulations show that density-driven flow in the fracture causes convective flow within the porous matrix and that the high-permeability fracture acts as a barrier for convection. Other simulations were run to investigate the influence of fracture incline on plume migration. Finally, tabular data of the tracer breakthrough curve in the inclined fracture is given to facilitate the verification of other codes.     

Modeling contaminant transport in homogeneous porous media with fractional advection-dispersion equation

The newly developed Fractional Advection-Dispersion Equation (FADE), which is FADE was extended and used in this paper for modelling adsorbing contaminant transport by adding an adsorbing term. A parameter estimation method and its corresponding FORTRAN based program named FADEMain were developed on the basis of Nonlinear Least Square Algorithm and the analytical solution for one-dimensional FADE under the conditions of step input and steady state flow. Data sets of adsorbing contaminants Cd and NH4+-N transport in short homogeneous soil columns and conservative solute NaCI transport in a long homogeneous soil column, respectively were used to estimate the transport parameters both by FADEMain and the advection-dispersion equation (ADE) based program CXTFIT2.1. Results indicated that the concentration simulated by FADE agreed well with the measured data. Compared to the ADE model, FADE can provide better simulation for the concentration in the initial lower concentration part and the late higher concentration part of the breakthrough curves for both adsorbing contaminants. The dispersion coefficients for ADE were from 0.13 to 7.06 cm2/min, while the dispersion coefficients for FADE ranged from 0.119 to 3.05 cm1.856/min for NaCI transport in the long homogeneous soil column. We found that the dispersion coefficient of FADE increased with the transport distance, and the relationship between them can be quantified with an exponential function. Less scale-dependent was also found for the dispersion coefficient of FADE with respect to ADE.     

Flow and transport through two-dimensional isotropic media of binary conductivity distribution. Part 1: NUMERICAL methodology and semi-analytical solutions

Flow and transport take place in a heterogeneous medium made up from inclusions of conductivity K submerged in a matrix of conductivity K ₀. We consider two-dimensional isotropic media, with circular inclusions of uniform radii, that are placed at random and without overlap in the matrix. The system is completely characterized by the conductivity contrast =K/K ₀ and by the volume fraction n. The flow is uniform in the mean, of velocity U=const. The derivation of the velocity field is achieved by a numerical method of high accuracy, based on analytical elements. Approximate analytical solutions are derived by a few methods: composite elements, effective medium, dilute systems and first-order approximation in logconductivity variance. The latter was employed by Rubin (1995), while the dilute system approximation was used by Eames and Bush (1999) and Dagan and Lessoff (2001). Transport is solved in a Lagrangean framework, with trajectories determined numerically from the velocity field, by particle tracking. Results for the velocity variance and for the longitudinal macrodispersivity, for a few values of and n, are presented in Part 2.     

Finite volume model for reactive transport in fractured porous media with distance- and time-dependent dispersion

Abstract

There are very few studies of fractured porous media that use distance- and time-dependent dispersion models, and, to the best of our knowledge, none which compare these with constant dispersion models. Therefore, in this study, the behaviour of temporal and spatial concentration profiles with distance- and time-dependent dispersion models is investigated. A hybrid finite volume method is used to solve the governing equations for these dispersion models. The developed numerical model is used to study the effects of matrix diffusion coefficient, groundwater velocity and matrix and fracture retardation factor on concentration profiles in the application of constant, distance-dependent and time-dependent dispersion models. In addition, an attempt is made to evaluate the applicability of these dispersion models by using the models to simulate experimental data. It was found that a better fit to the observed data is obtained in the case of distance- and time-dependent dispersion models as compared to the constant dispersion model. Thus, these numerical experiments indicate that distance- and time-dependent dispersion models have better simulation potential than the constant dispersion model.     

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

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

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

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