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.     

Investigation of flow and transport processes at the MADE site using ensemble Kalman filter

In this work the ensemble Kalman filter (EnKF) is applied to investigate the flow and transport processes at the macro-dispersion experiment (MADE) site in Columbus, MS. The EnKF is a sequential data assimilation approach that adjusts the unknown model parameter values based on the observed data with time. The classic advection–dispersion (AD) and the dual-domain mass transfer (DDMT) models are employed to analyze the tritium plume during the second MADE tracer experiment. The hydraulic conductivity (K), longitudinal dispersivity in the AD model, and mass transfer rate coefficient and mobile porosity ratio in the DDMT model, are estimated in this investigation. Because of its sequential feature, the EnKF allows for the temporal scaling of transport parameters during the tritium concentration analysis. Inverse simulation results indicate that for the AD model to reproduce the extensive spatial spreading of the tritium observed in the field, the K in the downgradient area needs to be increased significantly. The estimated K in the AD model becomes an order of magnitude higher than the in situ flowmeter measurements over a large portion of media. On the other hand, the DDMT model gives an estimation of K that is much more comparable with the flowmeter values. In addition, the simulated concentrations by the DDMT model show a better agreement with the observed values. The root mean square (RMS) between the observed and simulated tritium plumes is 0.77 for the AD model and 0.45 for the DDMT model at 328 days. Unlike the AD model, which gives inconsistent K estimates at different times, the DDMT model is able to invert the K values that consistently reproduce the observed tritium concentrations through all times.     

Diffusive mass exchange of non‐reactive substances in dual‐porosity porous systems – column experiments under saturated conditions

Diffusive mass exchange into immobile water regions within heterogeneous porous aquifers influences the fate of solutes. The percentage of immobile water is often unidentified in natural aquifers though. Hence, the mathematical prediction of solute transport in such heterogeneous aquifers remains challenging. The objective of this study was to find a simple analytical model approach that allows quantifying properties of mobile and immobile water regions and the portion of immobile water in a porous system. Therefore, the Single Fissure Dispersion Model (SFDM), which takes into account diffusive mass exchange between mobile and immobile water zones, was applied to model transport in well‐defined saturated dual‐porosity column experiments. Direct and indirect model validation was performed by running experiments at different flow velocities and using conservative tracer with different molecular diffusion coefficients. In another column setup, immobile water regions were randomly distributed to test the model applicability and to determine the portion of immobile water. In all setups, the tracer concentration curves showed differences in normalized maximum peak concentration, tailing and mass recovery according to their diffusion coefficients. These findings were more pronounced at lower flow rates (larger flow times) indicating the dependency of diffusive mass exchange into immobile water regions on tracers' molecular diffusion coefficients. The SFDM simulated all data with high model efficiency. Successful model validation supported the physical meaning of fitted model parameters. This study showed that the SFDM, developed for fissured aquifers, is applicable in porous media and can be used to determine porosity and volume of regions with immobile water. Copyright © 2015 John Wiley & Sons, Ltd.     

Preferential percolation quantified by large water content sensors with artifactual macroporous envelopes

For many scientific and practical tasks, it is important to estimate the soil–water percolation fluxes. This paper builds on measurements with large horizontal time‐domain reflectometry water content sensors in a loamy Mollisol. The sensors were installed into pre‐drilled holes and the gaps between them, and the soil was filled with a slurry of local soil with water. This gave rise to envelopes around them that contained artificial macropores. The sensors reacted to intensive rains by a rapid increase of their readings, often above the native soil's porosity, followed by an almost equally rapid decrease. The paper explores the feasibility of quantifying the rapid percolation, based on these anomalous water content peaks, and demonstrates that this is possible in principle, if the processes are simulated by a suitable model. A two‐dimensional dual porosity non‐equilibrium (mobile‐immobile) model was tried. The envelope around the sensor was modelled as an annulus with higher porosity and hydraulic conductivity, which attracts preferential flow and amplifies the percolation signal. With the model at hand, the flux hydrographs can be derived from model simulations and measured precipitation. For contrast, the Durner equilibrium dual porosity model was tried but was found little suitable. However, even the mobile‐immobile model did not perform perfectly. Simulated water contents were similar to the measured ones at some depths but not in the others, and the percolation fluxes were overestimated, compared to cumulative soil–water balance. Efforts to improve model performance were not successful. Hence, the model structure needs to be improved. Copyright © 2015 John Wiley & Sons, Ltd.     

Multitracer test approach to characterize reactive transport in karst aquifers

A method to estimate reactive transport parameters as well as geometric conduit parameters from a multitracer test in a karst aquifer is provided. For this purpose, a calibration strategy was developed applying the two-region nonequilibrium model CXTFIT. The ambiguity of the model calibration was reduced by first calibrating the model with respect to conservative tracer breakthrough and later transferring conservative transport parameters to the reactive model calibration. The reactive transport parameters were only allowed to be within a defined sensible range to get reasonable calibration values. This calibration strategy was applied to breakthrough curves obtained from a large-scale multitracer test, which was performed in a karst aquifer of the Swabian Alb, Germany. The multitracer test was conducted by the simultaneous injection of uranine, sulforhodamine G, and tinopal CBS-X. The model succeeds to represent the tracer breakthrough curves (TBCs) of uranine and sulforhodamine G and verifies that tracer-rock interactions preferably occur in the immobile fluid region, although the fraction of this region amounts to only 3.5% of the total water. However, the model failed to account for the long tailing observed in the TBC of tinopal CBS-X. Sensitivity analyses reveal that model results for the conservative tracer transport are most sensitive to average velocity and volume fraction of the mobile fluid region, while dispersion and mass transfer coefficients are least influential. Consequently, reactive tracer calibration allows the determination of sorption sites in the mobile and immobile fluid region at small retardation coefficients.     

Relationship between dual-domain parameters and practical characterization data

Dual-domain solute transport models produce significantly improved agreement to observations compared to single-domain (advection-dispersion) models when used in an a posteriori data fitting mode. However, the use of dual-domain models in a general predictive manner has been a difficult and persistent challenge, particularly at field-scale where characterization of permeability and flow is inherently limited. Numerical experiments were conducted in this study to better understand how single-rate mass transfer parameters vary with aquifer attributes and contaminant exposure. High-resolution reference simulations considered 30 different scenarios involving variations in permeability distribution, flow field, mass transfer timescale, and contaminant exposure time. Optimal dual-domain transport parameters were empirically determined by matching to breakthrough curves from the high-resolution simulations. Numerical results show that mobile porosity increases with lower permeability contrast/variance, smaller spatial correlation length, lower connectivity of high-permeability zones, and flow transverse to strata. A nonzero non-participating porosity improves empirical fitting, and becomes larger for flow aligned with strata, smaller diffusion coefficient, and larger spatial correlation length. The non-dimensional mass transfer coefficient or Damkohler number tends to be close to 1.0 and decrease with contaminant exposure time, in agreement with prior studies. The best empirical fit is generally achieved with a combination of macrodispersion and first-order mass transfer. Quantitative prediction of ensemble-average dual-domain parameters as a function of measurable aquifer attributes proved only marginally successful.     

Simulation of less‐mobile porosity dynamics in contrasting sediment water interface porous media

Considering heterogeneity in porous media pore size and connectivity is essential to predicting reactive solute transport across interfaces. However, exchange with less‐mobile porosity is rarely considered in surface water/groundwater recharge studies. Previous research indicates that a combination of pore‐fluid sampling and geoelectrical measurements can be used to quantify less‐mobile porosity exchange dynamics using the time‐varying relation between fluid and bulk electrical conductivity. For this study, we use macro‐scale (10 s of cm) advection–dispersion solute transport models linked with electrical conduction in COMSOL Multiphysics to explore less‐mobile porosity dynamics in two different types of observed sediment water interface porous media. Modeled sediment textures contrast from strongly layered streambed deposits to poorly sorted lakebed sands and cobbles. During simulated ionic tracer perturbations, a lag between fluid and bulk electrical conductivity, and the resultant hysteresis, is observed for all simulations indicating differential loading of pore spaces with tracer. Less‐mobile exchange parameters are determined graphically from these tracer time series data without the need for inverse numerical model simulation. In both sediment types, effective less‐mobile porosity exchange parameters are variable in response to changes in flow direction and fluid flux. These observed flow‐dependent effects directly impact local less‐mobile residence times and associated contact time for biogeochemical reaction. The simulations indicate that for the sediment textures explored here, less‐mobile porosity exchange is dominated by variable rates of advection through the domain, rather than diffusion of solute, for typical low‐to‐moderate rate (approximately 3–40 cm/day) hyporheic fluid fluxes. Overall, our model‐based results show that less‐mobile porosity may be expected in a range of natural hyporheic sediments and that changes in flowpath orientation and magnitude will impact less‐mobile exchange parameters. These temporal dynamics can be assessed with the geoelectrical experimental tracer method applied at laboratory and field scales.     

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.     

Investigation of Small‐Scale Preferential Flow with a Forced‐Gradient Tracer Test

A new tracer experiment (referred to as MADE‐5) was conducted at the well‐known Macrodispersion Experiment (MADE) site to investigate the influence of small‐scale mass‐transfer and dispersion processes on well‐to‐well transport. The test was performed under dipole forced‐gradient flow conditions and concentrations were monitored in an extraction well and in two multilevel sampler (MLS) wells located at 6, 1.5, and 3.75 m from the source, respectively. The shape of the breakthrough curve (BTC) measured at the extraction well is strongly asymmetric showing a rapidly arriving peak and an extensive late‐time tail. The BTCs measured at seven different depths in the two MLSs are radically different from one another in terms of shape, arrival times, and magnitude of the concentration peaks. All of these characteristics indicate the presence of a complex network of preferential flow pathways controlling solute transport at the test site. Field‐experimental data were also used to evaluate two transport models: a stochastic advection‐dispersion model (ADM) based on conditional multivariate Gaussian realizations of the hydraulic conductivity field and a dual‐domain single‐rate (DDSR) mass‐transfer model based on a deterministic reconstruction of the aquifer heterogeneity. Unlike the stochastic ADM realizations, the DDSR accurately predicted the magnitude of the concentration peak and its arrival time (within a 1.5% error). For the multilevel BTCs between the injection and extraction wells, neither model reproduced the observed values, indicating that a high‐resolution characterization of the aquifer heterogeneity at the subdecimeter scale would be needed to fully capture 3D transport details.     

A general approach to advective–dispersive transport with multirate mass transfer

A numerical solution that is significantly more general than other semi-analytical solutions is presented for governing equations describing advective–dispersive transport with multirate mass transfer between mobile and immobile domains. The new solution approach is general in the sense that it does not impose any restrictive assumption on the spatial or temporal variability of advective and dispersive processes in the mobile domain. A single integro-differential equation (IDE) is developed for the concentration in the mobile domain by separating the concentration in the immobile domain from the set of two partial differential equations. The solution to the IDE requires the evaluation of a temporal integral of the concentration in the mobile domain, which is a function of the Laplace transform of the distribution of the mass transfer rate coefficient. The Laplace transform is not limited to flow fields with known constant velocities. The solutions for one- and two-dimensional examples obtained using the new approach agree with those obtained by existing semi-analytical and numerical approaches.     

Lessons Learned from 25 Years of Research at the MADE Site

Field studies at well‐instrumented research sites have provided extensive data sets and important insights essential for development and testing of transport theories and mathematical models. This paper provides an overview of over 25 years of research and lessons learned at one of such field research sites on the Columbus Air Force Base in Mississippi, commonly known as the Macrodispersion Experiment (MADE) site. Since the mid‐1980s, field data from the MADE site have been used extensively by researchers around the world to explore complex contaminant transport phenomena in highly heterogeneous porous media. Results from field investigations and modeling analyses suggested that connected networks of small‐scale preferential flow paths and relative flow barriers exert dominant control on solute transport processes. The classical advection‐dispersion model was shown to inadequately represent plume‐scale transport, while the dual‐domain mass transfer model was found to reproduce the primary observed plume characteristics. The MADE site has served as a valuable natural observatory for contaminant transport studies where new observations have led to better understanding and improved models have sprung out analysis of new data.     

Changes in water table level influence solute transport in uniform porous media

Changes in the water table level result in variable water saturation and variable hydrological fluxes at the interface between the unsaturated and saturated zone. This may influence the transport and fate of contaminants in the subsurface. The objective of this study was to examine the impact of a decreasing and an increasing water table on solute transport. We conducted tracer experiments at downward flow conditions in laboratory columns filled with two different uniform porous media under static and transient flow conditions either increasing or decreasing the water table. Tracer breakthrough curves were simulated using a mobile–immobile transport model. The resulting transport parameters were compared to identify dominant transport processes. Changes in the water table level affected dispersivities and mobile water fractions depending on the direction of water table movement and the grain size of the porous media. In fine glass beads, the water flow velocity was similar to the decline rate of the water table, and the mobile water fraction was decreased compared with steady‐state saturated conditions. However, immobile water was negligible. In coarse glass beads, water flow was faster because of fingered flow in the unsaturated part, and the mobile water fraction was smaller than in the fine material. Here, a rising water table led to an even smaller mobile water fraction and increased solute spreading because of diffusive interaction with immobile water. We conclude that changes of the water table need to be considered to correctly simulate transport in the subsurface at the transition of the unsaturated–saturated zone. Copyright © 2014 John Wiley & Sons, Ltd.     

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.     

Effects of uncertainty of lithofacies, conductivity and porosity distributions on stochastic interpretations of a field scale tracer test

We investigate the importance of selecting two different methodologies for the determination of hydraulic conductivity from available grain-size distributions on the stochastic modeling of the depth-averaged breakthrough curve observed during a forced-gradient tracer test experiment. The latter was performed in the Lauswiesen alluvial aquifer, located near the city of Tübingen, Germany, by injecting NaBr into a well at a distance of about 50 m from a pumping well. We also examine the joint effect of the choice of the transport model adopted to describe solute transport at the site and the way the spatial distribution of porosity is assessed. In the absence of direct measurements of porosity, we consider: (a) the model used by Riva et al. (J Contam Hydrol 88:92–118, 2006; J Contam Hydrol 101:1–13, 2008), which relates the natural logarithms of effective porosity and conductivity through an empirical, experimentally-based, linear relationship derived for a nearby experimental site; and (b) a model based on a commonly used relationship linking the total porosity to the coefficient of uniformity of grain size distributions. Transport is described in terms of a purely advective process and/or by including mass exchange processes between mobile and immobile regions. Modeling of flow and transport is performed within a Monte Carlo framework, upon conceptualizing the aquifer as a random composite medium. Our results indicate that the model adopted to describe the correlation between conductivity and porosity and the way grain-sieve information are incorporated to depict the heterogeneous distribution of hydraulic conductivity can have relevant effects in the interpretation of the data at the site. All the conceptual models employed to describe the structural heterogeneity of the system and transport features can reasonably reproduce the global characteristics of the experimental depth-averaged breakthrough curve. Specific details, such as the peak concentration and the time of first arrival, can be better reproduced by a double porosity transport model when a correlation between conductivity and porosity based on grain size information at the site is considered. The best prediction of the late-time behavior of the measured breakthrough curves, in terms of the observed heavy tailing, is offered by directly linking porosity distribution to the spatial variability of particle size information.     

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.     

Exploring an ‘ideal hill': how lithology and transport mechanisms affect the possibility of a steady state during weathering and erosion

We present a model of chemical reaction within hills to explore how evolving porosity (and by inference, permeability) affects flow pathways and weathering. The model consists of hydrologic and reactive-transport equations that describe alteration of ferrous minerals and feldspar. These reactions were chosen because previous work emphasized that oxygen- and acid-driven weathering affects porosity differently in felsic and mafic rocks. A parameter controlling the order of the fronts is presented. In the absence of erosion, the two reaction fronts move at different velocities and the relative depths depend on geochemical conditions and starting composition. In turn, the fronts and associated changes in porosity drastically affect both the vertical and lateral velocities of water flow. For these cases, estimates of weathering advance rates based on simple models that posit unidirectional constant-velocity advection do not apply. In the model hills, weathering advance rates diminish with time as the Darcy velocities decrease with depth. The system can thus attain a dynamical steady state at any erosion rate where the regolith thickness is constant in time and velocities of both fronts become equal to one another and to the erosion rate. The slower the advection velocities in a system, the faster it attains a steady state. For example, a massive rock with relatively fast-dissolving minerals such as diabase – where solute transport across the reaction front mainly occurs by diffusion – can reach a steady state more quickly than granitoid rocks in which advection contributes to solute transport. The attainment of a steady state is controlled by coupling between weathering and hydrologic processes that force water to flow horizontally above reaction fronts where permeability changes rapidly with depth and acts as a partial barrier to fluid flow. Published 2020. This article is a U.S. Government work and is in the public domain in the USA.     

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.