Recent advances in anomalous transport models for predicting contaminants in natural groundwater systems

1 Motivation

Despite our image of Earth as the ‘blue’ planet, global supplies of clean water are arguably the most valuable, yet fragile, natural resources. Threats to freshwater arise from contamination by agricultural, industrial, urban, and resource mismanagement activities [1]. To minimize threats, a thorough understanding of the inherently complex structures that underpin flow and transport processes in water systems is necessary. A major challenge is that hydrologic systems are natural geologic systems with an extraordinary amount of inherent heterogeneity that is impossible to fully characterize, making predictions with high confidence difficult. While all hydrologic systems face this challenge, we will restrict our discussion to one: groundwater systems. However, many of the advances that we highlight have similar impact in other hydrologic settings such as streams and rivers.

At field scales, the most common models typically use average system quantities to define parameters for the flow equation and advection–dispersion equation (ADE). Such models serve as reasonable first estimates, but ultimately fail at capturing the extreme behaviors caused by the system's structural complexity [2]. For example, the United States National Research Council (e.g. [3,4]), as part of a larger set of studies focused on groundwater contamination, highlight that court ordered pump and treat strategies can fail to adequately remediate polluted sites as much as 90% of the time. In many cases, failure can be attributed to designs with conventional models that do not adequately represent system complexity. Failing 9 of 10 times is unacceptable by any engineering standard. Significant limitations of conventional approaches highlight the critical need for better theories and models that accurately incorporate, or at the very least acknowledge, the presence and impact of heterogeneity and natural variability in hydrologic systems.

The specific nature of heterogeneities can vary considerably and is not clearly known in hydrologic systems. After all, classifying something as heterogeneous only indicates that it is not homogeneous, but tells us nothing about its nature [5]. In groundwater flows, the spatial variability of the underlying geology causes this complexity. Real aquifers exhibit structure across an astonishing range of scales, from micrometer sized pores to structures on the order of tens to hundreds of kilometers. Heterogeneity comes in a variety of forms: permeability, which controls how quickly water can flow, can vary over many orders of magnitude depending on the material making up the aquifer; fractures and other preferential flow paths can exist that transmit water much more quickly than through the pore space in consolidated media; and geochemical heterogeneity due to spatial variability in mineral surfaces allows solutes to sorb and interact with the geologic medium over a broad range of rates. The extremely broad spatial and temporal variability in a geologic system makes it impractical to characterize it by a single effective parameter. For this reason, models like the ADE will always fail — as they are implicitly derived based on the assumption of narrow distributions and local processes. Any transport model seeking to improve the predictive capability beyond conventional models must incorporate the extremely broad range of spatial variability in geologic systems. For these models to be useful, they should capture this complexity using a minimal set of effective parameters.

2 Anomalous transport models

As highlighted, ADEs inherently have great difficulty in capturing observed behaviors in real hydrologic systems. Most forms of the ADE are implicitly built on the assumption of Fickian transport; that is, diffusive and dispersive mass transfer can be modeled as proportional to the concentration gradient. While, as shown in the seminal work of GI Taylor [6], this assumption may hold at asymptotic times, such timescales may be prohibitively large for practical interest. Transport that is not well described by Fickian dispersion is typically referred to as non-Fickian or anomalous. Indeed, anomalous transport is found so commonly in natural systems well beyond hydrology, that recently, one of the pioneers of anomalous models argued it be renamed ‘ubiquitous transport’ [7]. Although our discussion is restricted to the context of hydrology, it must be noted that anomalous transport models have substantially impacted other disciplines in the natural and physical sciences also [8].

One reason Fick's law makes sense and works remarkably well in certain contexts is that for a pulse initial condition, Fickian diffusion naturally results in a Gaussian concentration profile. Diffusion/Dispersion aims to represent random jumps solutes undergo due to subscale fluctuations (molecular for diffusion, velocity variation with dispersion). By the central limit theorem, such random jumps converge to a Gaussian distribution, provided jumps are independent and identically distributed (iid) with finite mean and finite variance. When solutes can make extraordinarily large jumps or be retained for very long times, this assumption becomes questionable. Convergence to Fickian behavior, while it may ultimately occur, will be at such large scales as to be irrelevant to a practitioner.

Thankfully, a rich family of models has emerged that relaxes this assumption that particle jumps follow narrow-tailed distributions. These models take advantage of ideas such as the generalized central limit theorem to still converge to analytically tractable solutions that can capture observed anomalous transport. Among these, perhaps the most popular are fractional Advection Dispersion Equations (fADE), continuous time random walks (CTRW) and multi-rate mass transfer (MRMT). It is important to note that these are certainly not the only anomalous transport models, and perhaps not even the most theoretically sophisticated (see [9] for a comprehensive discussion), but they are parsimonious and open sharing of computational toolboxes (e.g. [10,11]) has facilitated their application. All of these models can also be expressed from Lagrangian and Eulerian perspectives, enabling users to have a clear physical understanding of what the models aim to represent. Although these models are very closely related, it is common to find that a user's particular choice is tightly coupled to their conceptual model of the system they are studying.

While these models have been immensely successful in capturing real world observations across diverse hydrologic settings, they are not without shortcomings, among which we point out the following:

(i)
While agreement between model and measurement can be remarkably good, a common criticism is that these models are able to fit observations because of the increased number of parameters they have (for conservative transport 3–5, compared with 2). Unlike for example ADE models, it is difficult to obtain physically motivated estimates for these model parameters, and parameter estimation often becomes a fitting exercise (i.e. transport is measured and then fit, rather than actually predicted). Linking model parameters to physical characteristics of the system at hand therefore remains a central challenge. However, characterizing hydrologic systems to this level of detail has historically been difficult, due to inaccessibility of porous media.
(ii)
Most typically, the success of these models is assessed by comparing them to measured breakthrough curves (BTCs), which are concentration time series some distance downstream from an injection point. The BTC is often the only form of data that one can realistically obtain at scales of interest, but it has some implicit limitations. BTCs are an integrated measure that ultimately lumps many important processes together, making it difficult to disentangle individual controls. Experimental measurements and techniques that can help distinguish different processes are essential in providing a more physical basis for these models.
(iii)
As sophisticated as the aforementioned anomalous transport models are, they are often still built on strong assumptions. For example, while we can relax the assumption of finite mean and variance in the central limit theorem, the generalized central limit theorem still requires iid random variables. Many studies over the last decade have shown that the assumption of independence may not be adequate at the kinds of scales that practitioners are interested in, depending on the system at hand.
(iv)
All of the anomalous transport models listed above use dimension reduction to predict how solutes move downstream, that is, even though aquifers are clearly three-dimensional systems, these models are often developed in one-dimension and so any prediction that one obtains is an effective projected/mean concentration. (It is true that multi-dimensional forms of these models exist, but applications are more limited and often still typically involve some degree of reduction.) Given the integrative nature of BTCs noted above, this of course makes sense: it is risky to base predictions on a model that cannot be validated by experimental measurements. However, should one be interested in more complex nonlinear processes (e.g. mixing and chemical reactions), then a prediction of mean concentration is not sufficient and more information about subscale effects is needed.

3 Advances in anomalous transport models

3.1 Improved characterization and simulation capability

Some of the greatest advances in our understanding of flow and transport processes over the last decade can directly be attributed to techniques that enable better (visual) access to the internal structure of porous media. Innovative technologies, including micro-computed tomography (micro-CT) and nuclear magnetic resonance among others, have enabled us to obtain three-dimensional images of the internal structure of complex real geologic porous media at a scale and resolution that was previously unobtainable [12,13]. Similar parallel advances in computational resources have made accessible open source, state of the art computational fluid dynamics packages (e.g. [14]), that simulate the complete velocity field within the imaged complex porous structures. These new simulation approaches enable calculation of useful quantities such as velocity probability distributions, which could previously only be inferred indirectly by inverse modeling of BTCs and the likes. Furthermore, flow simulations can be coupled with Lagrangian and Eulerian transport codes, enabling high-resolution simulation of transport processes, including solutes undergoing mixing-driven and heterogeneous reactions.

Highlighted paper: A large number of papers using high resolution imaging have emerged from a research group at Imperial College London. Here we highlight one [15•], which we feel demonstrates the power of these techniques. The authors develop a particle-based method to simulate dissolution reactions at pore scales using voxelized three-dimensional micro-CT images. Their approach is validated against a dynamic imaging experiment where a Ketton oolite is imaged during CO2-saturated brine injection under reservoir conditions, again exploiting advances in imaging technology [16] (see Figure 1). The model results agree well with measured changes in porosity and permeability, and the spatial distribution of the dissolution front is correctly replicated. Advances on observation enabled a physically based model, capable of reproducing behavior in a highly complex setting that previously would have been entirely empirical and whose generality might be questionable.

3.2 Improved experimental techniques

Similarly, experimental breakthroughs on directly observing behaviors in the porous medium have occurred. As noted, a major limitation to date is our inability to directly observe and measure transport within the porous medium, thus typically ending up with integrated measures such as BTCs that do not enable direct inference of actual processes. A major breakthrough in this regard has been using refractive index matching (RIM) approaches between the fluid and the porous medium. This renders the solid phase transparent and allows direct visual access to the internal flow. These advances, coupled with particle tracking velocimetry (PTV) and particle imaging velocimetry (PIV) techniques, allows for direct characterization of velocity fields.

Highlighted paper: Here we highlight the work of [17•] who studied the evolution of velocity in time in porous media by experimentally tracking tracer particles moving through a transparent, 3-D synthetic sandstone. Using state of the art PTV, they measured the correlated nature of velocities along stream lines. The observed behavior is well described by a correlated CTRW model and is one of the first examples where such a model is derived from experimental rather than numerical data. The model is based on an Ornstein–Uhlenbeck [18] process to model velocity evolution and can be quite simply parameterized with the correlation length as well as the mean and standard deviation of the velocity distribution, which may be obtainable in general settings. The same methods were later used to study flow evolution in a porous medium gradually invaded by biofilm [19•] (Figure 2).

3.3 Newer models that relax previous assumptions

About a decade ago, Le Borgne et al. [20,21] introduced a model that now goes by the name the spatial Markov model (SMM). It is closely related to CTRWs and fADEs, but differs in that it imposes correlation between successive jumps, relaxing the assumption of independent increments. The authors found that velocity correlations in space were sufficiently short range that successive velocities could be represented by a Markov process. They originally studied it in heterogeneous porous media at geologic scales, but since then it has been broadly shown to work well in fractured media (e.g. [22]), pore scale settings (e.g. [23]) and beyond. It was also shown that particles display highly intermittent behavior, alternating between quiescent periods of low velocities and small accelerations and energetic periods of large velocities and large accelerations. Such intermittent behavior necessitates a correlated model with the specific characteristics of the SMM [24]. Another development we follow with interest has been the PhEDEX model (Phase Exposure-Dependent EXchange) [25], which builds on the ideas of MRMT. At the root of this model is the idea that solute concentration evolves with respect to two separate times: time-mobile as well as exposure time (i.e. when a particle is exposed to another process — such as immobilization). One of the main points of interest of the PhEDEX is that it casts the problem in a way that clearly separates transport and delay mechanisms. In many other anomalous transport models these mechanisms are lumped together via complex convolutions with memory functions. This paves the way for physical parameterization of the model processes, although to our knowledge, evidence to demonstrate these attributes is not yet experimentally available.

Highlighted paper: Here we highlight the work of [26], who applied the SMM to tracer experiments in fractured media at a field site in France (see Figure 3). What stands out about this paper is that it is, to our knowledge, the first application of the SMM to field data that did not rely on high resolution numerical simulations to measure travel time distributions and parameterize correlation effects. The authors introduce a relatively parsimonious model (conceptually similar to the one used by [17•] noted above) with which they successfully capture field experiment data of characteristic anomalous transport behavior in a real complex setting. This work paves the way for applications in practical settings, which to date have not been achievable. While other approaches with real experimental data are possible (e.g. [27,28]) the parsimony of [26] is elegant and appealing.

3.4 Upscaling in three-dimensions

As noted, dimension reduction is common for many anomalous transport models. However, hydrologic systems are three-dimensional and much important information can be lost when not considering the full three-dimensional concentration field. Increasingly, authors are considering spreading in multiple dimensions (e.g. [29]), but do not always consider the full coupling between longitudinal and transverse dimensions. A proper description of this coupling is necessary to capture the pronounced spatial coherence of advective processes, which strongly determines mixing behaviors. Part of the challenge lies in the complexity of flow and transport, as well as in available techniques to quantify large-scale behavior in three spatial dimensions. [30] demonstrated this by studying particle trajectories in a Doddington sandstone sample, showing transport processes are strongly correlated in all three directions such that full parameterization would require correlation descriptions using a nine-dimensional set of transition matrices. While they showed this approach to be effective, it is impractical for real porous media. As an alternative, novel trajectory-based methods are emerging to upscale simulated particle positions in a physically consistent manner (e.g. [31,32]).

Highlighted paper: The recent paper we highlight here is [33•]. Using high resolution trajectories obtained from simulations in a Doddington sandstone, they proposed a novel training trajectory method, where trajectories are cut into small fragments that are then stitched together into much longer trajectories that ensure continuity of velocity magnitude and direction. Using this method, the authors fully reconstructed BTCs and dilution profiles obtained from much more costly direct numerical simulations. The method, inspired by training image methods in geostatistics, is in our view radically different from previous approaches and has great potential for application in other settings and at other scales that are difficult to upscale in three-dimensions.

3.5 Connecting spreading and mixing behaviors

A wide range of nonlinear processes in hydrologic systems depend directly on local concentrations of interacting solutes and concentration gradients, including kinetic/equilibrium reactions and biological activity [34,35]. Most anomalous transport models, due to dimension reduction, cannot explicitly account for concentration variations, but instead provide a measure of plume spreading and mean concentrations. While mixing and spreading are fundamentally different, recent studies show that they are tightly linked (e.g. [36]), illustrating the potential to utilize existing transport theory to predict mixing and reactions. One promising advance is lamella theory, which conceptualizes a mixing interface as a line that distorts into lamellar structures as it stretches and folds in a heterogeneous velocity field [37] (see Figure 4). The lamellae eventually coalesce by diffusion, resulting in a progression through multiple mixing regimes over time. By recognizing these different processes and regimes, the lamellar framework predicts the global evolution of mixing based on limited information relating to structural heterogeneity and plume spreading characteristics. The framework has also been used to predict how chemical reactions progress in a complex porous medium [38,39]