Research Article  Open Access
Challenges in the Application of Fractional Derivative Models in Capturing Solute Transport in Porous Media: DarcyScale Fractional Dispersion and the Influence of Medium Properties
Abstract
Heterogeneous media consisting of segregated flow regions are fractionalorder systems, where the regionalscale anomalous diffusion can be described by the fractional derivative model (FDM). The standard FDM, however, first, cannot characterize the Darcyscale dispersion through repacked sand columns, and second, the link between medium properties and model parameters remains unknown. To fill these two knowledge gaps, this study applies a tempered fractional derivative model (TFDM) to capture bromide transport through laboratory repacked sand. Column transport experiments are conducted first, where glass beads and silica sand with different diameters are repacked individually. Latetime tails are observed in the breakthrough curves (BTC) of bromide even in relatively homogeneous glass beads. The TFDM can capture the observed subdiffusion, especially the latetime BTC with a transient declining rate. Results also show that both the size distribution of repacked sand and the magnitude of fluid velocity can affect subdiffusion. In particular, a wider sand size distribution or a smaller flow rate can enhance the subdiffusion, leading to a smaller time index and a higher truncation parameter in the TFDM. Therefore, the Darcyscale dispersion follows the tempered stable law, and the model parameters might be related to the soil size and flow conditions.
1. Introduction
Geological formations usually exhibit multiscale physical and/or chemical heterogeneity, which can lead to the space and/or time nonlocal dependency for solute transport (see the extensive review by Zhang et al. [1]). For example, the existence of preferential flow paths can force the fast movement of dissolved contaminants, resulting in superdiffusion. The mass flux at any position therefore depends on the flux not only at adjacent neighbors but also at a wide range of upstream zones. Such spatial nonlocal dependency can be efficiently characterized by the space fractional derivative, which is a generalization of its integerorder counterpart [2, 3]. In addition, the sorptiondesorption mechanism can cause nonMarkovian evolution of tracer mass in time, a typical subdiffusion that has been described by the time fractional derivative model (FDM) [1].
Two major knowledge gaps, however, remain for the FDM. First, while the FDM has been applied by hydrologists to simulate contaminant transport through regionalscale heterogeneous porous and fractured media for more than a decade [4], its applicability for Darcyscale dispersion remains obscure. Indeed, some studies implied that the standard FDM may not be applicable for smallscale dispersion [5], due to the discrepancy between the finite medium size and the infinite distribution of particle dynamics (i.e., jump sizes and waiting times) assumed by the standard FDM. Second, the link between medium properties and FDM parameters has not been evaluated systematically. This unknown relationship, as commented by Neuman and Tartakovsky [6], suggests a failure of the physical model itself at the Darcyscale.
This study attempts to fill the above knowledge gaps. We apply a tempered fractional derivative model (TFDM), which is a generalization of the standard FDM, to capture bromide transport through laboratory repacked sand. This way, the applicability of the fractionalorder partial differential equation on Darcyscale dispersion can be tested reliably. Hence the first knowledge gap can be filled. The combined study of laboratory experiments and stochastic analysis may also reveal the trend of major transport parameters varying with sand properties. Such trend might lead to the answer regarding the second knowledge gap.
Laboratory experiments of solute transport through sand columns were conducted extensively by the hydrology community to explore the dynamics of dissolved solutes. For example, recent experiments [7, 8] identified nonFickian diffusion for passive tracer transport through repacked laboratory columns of macroscopically homogeneous sand. Curvefitting applications further show that the nonFickian diffusion characterized by the nonsymmetric plume cannot be explained efficiently by the 2ndorder advectiondispersion equation (ADE) model based on Fickian diffusion [6, 9]. Welldesigned laboratory experiments and alternative conceptual models are needed to explore the nature of transport through sand columns that may have been missed or misinterpreted previously, and then, to build the link between medium properties (measurable in the laboratory) and model parameters (controlling diffusion).
The rest of the paper is organized as follows. In Section 2, we introduce the laboratory experiments conducted to explore the Darcyscale dispersion in repacked sand columns with different filling materials and to evaluate the influence of soil properties and fluid velocity on the transport of bromide. In Section 3, a fractional derivativebased, nonlocal transport model is developed and compared with the other nonlocal models. In Section 4, the proposed model is used to capture the observed Darcyscale transport. In Section 5, we discuss the factors that affect bromide transport. Subdiffusion dominated by either slow advection or diffusion is also discussed to further explore the physical nature of the observed anomalous transport. Conclusions are summarized in Section 6.
2. Laboratory Experiments
2.1. Experimental Setup
We conducted laboratory experiments to measure the breakthrough curve (BTC) of one conservative tracer (bromide, as NaBr) in various sand packs. A glass column with internal dimensions of 150 mm (length) × 15.9 mm (diameter) was packed with glass beads or silica sand (Figure 1(a)).
(a)
(b)
Three different types of column experiments were conducted. For the first type of experiment (denoted by Run 1 in the following), the column was filled with uniform glass beads with an average diameter of 0.4 mm (Figure 1(b)) to represent a “homogeneous” porous medium microstructure.
For the second type of experiments (denoted by Run 2), two different sizes of glass beads were well mixed and packed, forming networks with mobile and relatively immobile domains. The first group of these experiments included glass beads with average diameters of 1 and 0.2 mm, while the second group included glass beads with average diameters of 0.4 and 0.2 mm.
For the third type of experiments (denoted by Run 3), silica sand with a specific particle size distribution was used to represent “natural” soil, where the irregular shape of the sand may affect the interconnected pores and the corresponding tracer dynamics. The overall particle size distribution was obtained by combining the following size fractions obtained by sieving: 0.85~1.0 mm (representing the coarse sand), 0.35~0.425 mm (medium sand), and 0.15~0.25 mm (fine sand), respectively. For description simplicity, in the following we denote the three size fractions by 1, 0.4, and 0.2 mm size silica sand, respectively, corresponding to the glass beads with similar diameters used in the second type of experiment. The sand was then cleaned with acid before packing.
After the column was packed, the following steps were performed to obtain BTCs. A fivepoint calibration of the bromide ion selective electrode (ISE) (Orion) was performed. The potential of standard solutions was measured from the lowest to the highest bromide concentration. Deionized water (DI) was run through the column for at least 2 hours prior to tracer injection to establish the flow domain. The pulse of bromide (of volume 10 mL) was then injected into the column (oriented horizontally) at a concentration of 1 mol/L, and discrete samples were collected from the outlet using a fraction collector (Teledyne ISCO). The sampling interval at early and late times was shortened to better record the tails of the BTC, known to be critical signals of nonFickian transport. The bromide concentration was measured in all collected samples using the calibrated ISE probe with a detection limit of 10^{−5} mol/L.
The flow velocity was adjusted during the experiment using a rotary peristaltic pump and controller (Cole Palmer Masterflex) (Figure 1(a)) to evaluate the potential influence of flow rate on bromide dispersion. The average flow rate, and, hence, linear velocity were larger than those in real aquifers, to shorten the experimental periods and to provide advectiondominated transport systems.
2.2. Experiment Results
2.2.1. Run 1: Homogeneous Glass Beads
Figure 2 shows the measured BTC for this run. The BTC is symmetric for most of the times (Figure 2(a)), which can be explained by the classical Fickian diffusive model: where is the density (or the tracer concentration in this case), is the average linear velocity, and is the macroscopic dispersion coefficient.
(a)
(b)
A slight latetime tail in the BTC (Figure 2(b)) close to the detection limit implies that transport can be subdiffusive even in homogeneous media, although the subdiffusive portion can only be detected at low concentrations.
2.2.2. Run 2: Heterogeneous Glass Beads
The BTC measured in Run 2 contains a much heavier latetime BTC tail than the “homogeneous” case (Figure 3(a)) used in Run 1. The mixture of glass beads in this run has a wider size distribution of glass beads (i.e., more “heterogeneous” than the one used in Run 1). In addition, the latetime BTC tail for the larger diameter (1 + 0.2 mm) mixture is also heavier than that for the smaller diameter (0.4 + 0.2 mm) mixture (Figure 3(a)).
(a)
(b)
Note that the BTC tail declines faster for a higher fluid velocity, as shown by the glass beads with a 0.4 + 0.2 mm mixture (Figure 3(a)). The latetime BTC tail for the larger diameter (1 + 0.2 mm) mixture, however, seems less sensitive to fluid velocity (Figure 3(a)), probably due to the short duration of the experimental time. Thus, the experimental time for this specific case is apparently not long enough to capture any response of BTC to the variation of fluid velocity. Note that a dimensionless time scale is used in Figure 3. In Section 4, we will simulate the measured BTC using stochastic models, where the subtle discrepancy for BTCs with different fluid velocity might be gleaned.
2.2.3. Run 3: Heterogeneous Silica Sand
The observed BTC for Run 3 also contains an apparent latetime tail (Figure 3(b)), showing the strong subdiffusive process. The discrepancy between different mixtures, however, is not as apparent as that for glass beads.
The latetime tail in the BTC shrinks with the increase of fluid velocity (Figure 3(b)), similar to the behavior found in Run 2. The influence of fluid velocity on BTC will be viewed further using a numerical model (see Section 3.2) and a dimensional time scale (shown in Figure 4).
(a)
(b)
(c)
(d)
3. The Tempered Fractional Derivative Model
3.1. Review of Nonlocal Transport Methods
Nonlocal transport theories were developed recently to capture nonFickian diffusion, as reviewed extensively by Haggerty et al. [10], Berkowitz et al. [9], Neuman and Tartakovsky [6], and Zhang et al. [1]. The most efficient model for laboratoryscale transport was found to be the continuous time random walk (CTRW) framework by Berkowitz and Scher [11]. The CTRW framework defines empirical distributions for the transition time of solute particles after experiencing enough variations of local velocity. Levy and Berkowitz [7] found that, if the transition time has a power law tail (where ), the CTRW captures the observed nonFickian diffusion in sandboxes filled with homogeneous or heterogeneous sand, where the exponent decreases with increasing fluid velocity. Berkowitz and Scher [11] extended the CTRW model used by Levy and Berkowitz [7] by assigning a truncated power law for the transition time (see also Table 1), where is a median time for transitions between sites and is the cutoff time of the power law spectrum. From the point of view of random walks [1], the transition time also represents the time for each particle to move. Hence the standard CTRW model actually assumes that all solute particles are in motion all the time. In other words, the subdiffusion is assumed to be the result of slow advection, as also shown by Berkowitz and Scher [11].

Molecular diffusion, however, may also cause the subdiffusive effect, as suggested by the physical process of multiple rate mass transfer [10]. After solutes whose transport is controlled by advection are flushed out, diffusion out of the relatively immobile domains causes later arrivals and the apparent latetime tail of a breakthrough curve. Repacked soils in the laboratory can contain immobile or stagnant regions, where the effect of diffusion on subdiffusion should not be neglected. In this study, we check this mechanism and compare it with the slow advectionrelated subdiffusion.
3.2. The TFDM for DiffusionDominated Subdiffusion
The tempered stable model proposed by Meerschaert et al. [5] is a concise version of the multirate mass transfer model with a finite number of rate coefficients. It contains the least number of parameters and can be computationally efficient, if solved appropriately. Hence we choose it as the appropriate model for the diffusiondominated subdiffusion (i.e., subdiffusion due to the effect of slow diffusion of solute particles).
In our representation, we propose the following tempered fractional derivative model, or TFDM, by generalizing the current time fractional derivative models [1] and the tempered stable model [5]where and denote the resident concentration in the total and mobile phase, respectively, denotes the capacity coefficient, is the truncation parameter, ( in this study) is the scale index in time characterizing the powerlaw slope of waiting times, () is the space index characterizing the displacement of the plume front, and accounts for the initial condition. When , the TFDM (2a) and (2b) reduces to the tempered stable model proposed by Meerschaert et al. [5]. In addition, the RiemannLiouville fractional derivative is used for both the space and time fractional derivatives. This type of fractional derivative is selected because the corresponding Langevin method is known [12] and used for numerical approximations in this study.
Model (2a) and (2b) can be derived using either the fractal mobile/immobile (FMI) approach [4] or the subordination approach [5]. In the FMI approach, the generalized transport equations for the total and mobile concentrations are where the symbol denotes convolution, represents a generalized memory function, denotes the initial mass (which can be normalized to 1), and the operator describes the flux due to advection and dispersion. When the memory function and the nonlocal dispersive flux are used, the above model reduces to (2a) and (2b). In the 2nd approach (subordination), the total concentration can be expressed as , where denotes the operational time and and are densities of random walking particles in the mobile and immobile phases, respectively [12]. The governing equation for the total concentration therefore is model (2a), if the first density is governed by the motion process and the second density is governed by the waiting time process [12]. Similar arguments (see [5]) lead to the governing equation for the mobile concentration, which is (2b) in this case.
The two fractional derivative terms in model (2a) and (2b) have specific physical meanings that may help to explain the experimental data. First, the time fractional term is used to distinguish solute particle status (i.e., mobile versus immobile). In particular, the time drift term () is assigned to mobile particles, with the waiting time (represented by the two remaining terms on the left hand side of (2a) and (2b)) for immobile particles. The physical time therefore increases linearly when particles are in motion, and then it has a positive dispersive component (because ) for each immobile particle. The evolution of physical time due to drift and dispersion in time is analogous to the advective and dispersive displacement for solute particles in space—the advective term accounts for the mean solute displacement, while the dispersive term adds random noise (either positive or negative) caused by the deviation of local velocities. The distinction for solute particle status is required for field applications, where the mobile concentration or mass can differ significantly from that in the immobile or total phase [13], especially at early or late times. Second, the space fractional derivative term on the right hand side of (2a) and (2b) describes the possible fast motion through preferential flow paths. In practical applications, both the early and late arrivals can be critical [1]. Heavy leading edges of tracer plumes have been observed for regionalscale transport (see, e.g., Adams and Gelhar [14]). These observations lead to a specific question: will any smallscale preferential flow path in saturated repacked sand generate the leading edge of a plume? The BTCs observed above provide the first hand material to answer this question. It is also noteworthy that the CTRW framework used by Levy and Berkowitz [7] and Berkowitz and Scher [11] does not have the above two properties. Further comparison of the TFDM (2a) and (2b) and the standard CTRW framework can be seen in Table 1.
The model in (2a) and (2b) can be approximated by a spatiotemporal Lagrangian solver proposed by Zhang et al. [12]. The approximation for the classical fractional derivative models is combined with the exponential rejection method proposed by Baeumer and Meerschaert [15] (that can generate the tempered stable random variables) to form a fully Lagrangian solver for (2a) and (2b). The resultant particle tracking scheme is similar to the one proposed by Zhang and Papelis [16], where the time and space fractional terms were separated (in different models).
4. Applications
Model (2a) and (2b) now can be used to capture the BTCs described in Section 2. Note that the resident concentration (i.e., solution of (2a) and (2b)) needs to be transformed to its flux counterpart (i.e., the BTC). We first fit Run 1 (with homogeneous glass beads) by adjusting the dispersion coefficient , capacity coefficient , truncation parameter , and the two scale indexes and . Note that the five parameters have different impacts on the BTC. For example, affects the peak concentration, dominates the early tail, controls the slope of the latetime tail, affects the mass partition for particles at different phases, and controls the transition time from the power law tail to the exponential one. This helps us find quickly the bestfit value for each parameter. The average linear velocity (3.56 cm/min) was measured in the laboratory. Results (Figure 2) show that TFDM (2a) and (2b) can fit the observed BTC, with the bestfit parameters cm^{2}/min, min^{−0.01}, min^{−1}, , and . The standard FDM, however, slightly overestimates the latetime BTC tail (Figure 2).
For Run 2 using 1 + 0.2 mm glass beads, we first fit the BTC using the TFDM (2a) and (2b) with the slowest flow velocity ( cm/min, see the black line in Figure 4(a)). The bestfit parameters, including cm^{2}/min, min^{−0.1}, min^{−1}, , and , were then used to predict the BTC for the other two cases with larger velocities. Predictions of model (2a) and (2b) generally match the measured BTCs, while the classical secondorder ADE (shown by the dashed line in Figure 4(a)) underestimates the BTC latetime tail.
For Run 2 using 0.4 + 0.2 mm glass beads, the bestfit parameters for the BTC with the slowest flow velocity (3.89 cm/min) are cm^{2}/min, min^{−0.02}, min^{−1}, , and . Model prediction, however, overestimates the latetime tail for BTC with velocities and (shown by the solid lines in Figure 4(b)). The relatively lighter BTC tail due to a larger velocity can be captured by a slightly larger truncation parameter and/or a smaller capacity coefficient (see the dotted lines in Figure 4(b)). For example, the bestfit parameters for the BTC with a larger velocity (5.09 cm/min) are min^{−0.02}, min^{−1}, and .
The same conclusion is found for Run 3 with silica sand (Figures 4(c) and 4(d)), where a higher velocity corresponds to a larger and/or a smaller . Note, however, that model parameters are not sensitive to the size distribution for this run, which is consistent to the measurements described in Section 2.
5. Discussion
5.1. Factors Affecting the Subdiffusion of Bromide and the Model Parameters
Laboratory experiments show that the subdiffusion increases by increasing the range of the sand size distribution, especially when the column is filled with glass beads. A wider sand size distribution tends to enhance subdiffusion, because broader distributions of particle diameter more readily form immobile regions. Bromide transport through silica sand is not as sensitive to the size distribution as the glass beads, which might be due to either the strong influence of the irregular shape of silica sand on the mass exchange between mobile and relatively immobile regions and/or the relatively large flow rate required for the laboratory experiments that counterbalances the size effect. Future studies are needed to explore further the influence of sand shape and low flow rate on subdiffusion.
Water flow rates across the sand column also affect subdiffusion. As the flow rate in nonaggregated material increases, the percentage of the total domain dominated by diffusive transport decreases. Hence the increase of fluid velocity likely decreases the contribution of diffusion to the arrival times of solute particles. In particular, if the observation period is short, the latetime subdiffusive behavior affected by the flow rate may not be detected. Hence the total experimental period should be as long as possible, to identify the full behavior of subdiffusion at late times.
The TFDM parameters can efficiently capture the subtle variation of subdiffusion. For example, the temporal scale index increases (representing the decrease of subdiffusion) with the decrease of the size range of mixed glass beads, given the relatively declining contribution of immobile regions to subdiffusion. Meanwhile, the capacity coefficient decreases, also illustrating the decline of subdiffusion. Though our results are compelling, they are not sufficient to build a quantitative relationship between the TFDM parameters and medium heterogeneity. To establish a purely predictive physical model, substantial effort involving laboratory experiments, analytical analysis, and numerical evaluations is still needed.
In addition, the standard FDM tends to overestimate the latetime BTC tail (see, e.g., the dashed line in Figure 2), since it assumes an infinite waiting time distribution. In a typical sand column at the Darcyscale, the maximum trapping period (also known as the residence time) of solute particles may be finite. In other words, the waiting time distribution may have an upper limit. Such limit can be captured efficiently by the TFDM using the truncation parameter.
Finally, the bestfit space scale index in model (2a) and (2b) is limited to 2.0 for all the observed BTCs in this study, implying that the dynamics for solute particles in mobile time are limited to Brownian motion (with a drift). As shown in Figure 4, the early tail of the BTC is as steep as an exponential function. The lack of an apparent leading edge confirms the difference between the repacked sand and realworld soils; clearly, simulating the realworld fast motion paths using repacked sand is difficult, if not impossible. In contrast, the realworld subdiffusive behavior can be captured by laboratory experiments, most likely due to the insensitivity of subdiffusion to the exact location of immobile regions [1].
5.2. The Slow Advection Dominated Subdiffusive Model and Its Limitations in Capturing RealWorld Transport
To further understand the subdiffusive process, we simulate again the measured BTCs by assuming that the observed subdiffusion is driven by slow advection. A time fractional derivative model can be built to describe an advectiondominated subdiffusion. Assuming a CTRW with independent jump sizes and waiting times (or the elapsed time during two subsequent jumps), the corresponding scaling limit is which is the reduced form of the TFDM (2a) and (2b) without the time drift term and capacity coefficient. Model (4) can also be derived directly by assuming that the drift in time is zero in the TFDM model (2a) and (2b).
Applications show that model (4) captures most BTCs for Run 2 with 1 + 0.2 mm glass beads (Figure 5(a)). The bestfit parameters are cm/min, cm^{2}/min, min^{−1}, , and , for the BTC with the smallest fluid velocity (4.66 cm/min). Theoretically, the space scale index (capturing superdiffusion) and the time scale index (representing the degree of subdiffusion) are independent. The truncation parameter however might be related to , since both of them describe the waiting time distribution of tracer particles and they all depend on properties of relatively immobile domains. A future study with extensive laboratory experiments is needed to reveal the quantitative relationship among model parameters. It is also noteworthy that both the time scale index and the truncation parameter are larger than those for model (2a) and (2b). The relatively large and in model (4) have to be used to capture the relatively steep power law latetime tail of the BTC (with a slope ~−6.5 in a loglog plot). In model (2a) and (2b), the steep latetime BTC is explained by the relatively small time ratio that particles spend in the immobile and mobile phases (so that particles exit the immobile phase quickly and form the steep latetime tail of BTC), which can be captured conveniently by a small value of capacity coefficient . Therefore, model (2a) and (2b) can describe a wide range of BTCs with various latetime tails, while model (4) has limited capability due to the lack of the controlling parameter . In addition, the bestfit velocity in model (4) is smaller than the measured average linear velocity, an artificial effect due to more jumps than that in model (2a) and (2b) (because model (4) misses the actual mobile time; see also Zhang et al. [1]).
(a)
(b)
Further applications show that model (4) cannot fit any other BTCs. For example, model (4) misses the BTC tail for Run 2 with 0.4 + 0.2 mm glass beads (Figure 5(b)). Fit 1 and fit 2 shown in Figure 5(b) represent two fitting results using (4). The parameters used are min^{−1} and for fit 2 and min^{−1} and for fit 1. The scale index approaches the maximum limit 1, and the truncation parameter cannot improve the fit at all.
It is possible to extend model (4) to capture BTCs similar to the one shown in Figure 5(b), by using the time fractional derivative model with two time scales proposed by Meerschaert et al. [17] (so that can be larger than 1). This extension, however, still has two serious limitations. First, the CTRW model differs significantly for different ranges of index and it represents different physical processes [17]. Second, model (4) cannot capture the mobile mass decline, no matter the range of . Therefore, the TFDM model (2a) and (2b) is superior to its simplified version (such as model (4)) in capturing realworld subdiffusion.
5.3. Applicability and Limitation of the TFDM (2a) and (2b)
The TFDM model (2a) and (2b) may also be used to capture open channel flow and transport, such as the transport of dye in rivers [1, 4, 10]. Although advection is the dominant factor for transport in surface systems, the observed latetime tail of the BTC for a dye is caused by the molecular diffusion during the mass exchange between open channel and hyporheic zone [4] or the many relatively immobile domains in natural rivers [1]. The mechanism for the diffusionrelated subdiffusion is therefore similar to that discussed above for porous media. We will check the applicability of the model (2a) and (2b) in surface dynamic processes in a future study.
The TFDM (2a) and (2b) may be applied to dynamic processes observed in the other fractionalorder systems in multiple disciplines, such as sedimentation engineering and chemical engineering. For example, quantifying anomalous dynamics of suspended and bedload sediment transport in natural rivers remains a significant challenge in river morphology studies, due to the stochastic nature of sediment transport in a complex system with multiscale intrinsic heterogeneity. The TFDM (2a) and (2b) may capture the random process of sediment transport, especially the intermittent mobile and immobile dynamics. In addition, anomalous kinetics is well documented for chemical reactions, where the nonFickian motion of reactant molecules (the main reason why the diffusionlimited anomalous kinetics deviate significantly from the thermodynamic law) may be efficiently simulated by the TFDM (2a) and (2b) with a particlebased scheme.
One of the major limitations of the TFDM (2a) and (2b) is the representative scale of model parameters. Nonlocal transport models are upscaling tools that replace detailed medium heterogeneity information with memory kernels in space and/or time. How to define the representative scale and how to delineate the effective range for the model index in nonstationary systems remain to be shown. This study shows that the representative scale for the tempered fractional diffusion is no less than the laboratory scale. Will the variableorder or the distributedorder fractional derivatives (where the order of the fractional derivative is no longer a constant) capture the evolution of heterogeneity, and should the nonlocal transport models be conditioning on local system properties measured at each representative scale? These remain open questions. In addition, the TFDM (2a) and (2b) differs significantly from the standard fractional derivative models because the former is scale dependent. Can the TFDM (2a) and (2b) capture the scaling behavior for transport observed in practical engineering processes? We will focus on these questions in a future study.
6. Conclusions
The fractional engine is a promising tool to capture anomalous dispersion in heterogeneous media, but major challenges do exist, including the Darcyscale fractional dispersion and the influence of medium heterogeneity. In this study, laboratory experiments were combined with stochastic model analysis to explore the applicability of the fractional engine in capturing Darcyscale dispersion in sand columns filled with various materials and to explore the potential link between medium properties and model parameters. The following four main conclusions are drawn.(1)The tempered fractional derivative model can capture subdiffusion at the Darcyscale. The physical model distinguishes solute status, contains the least number of parameters, and can be extended conveniently to capture advanced transport processes. Most importantly, the TFDM can characterize the transient decline rate of the latetime BTC, probably due to the finite distribution of particle waiting times, while the standard FDM tends to overestimate the latetime tail of BTC.(2)Both the sand particle size distribution and the fluid velocity can affect the Darcyscale subdiffusion. All of the measured BTCs of bromide contain an apparent latetime tail, which is heavier for a wider particle size distribution of sand or a smaller fluid velocity. These two properties can enhance the relative contribution of diffusion to the latetime arrivals of solute particles. Hence both medium properties and flow conditions can affect subdiffusion, which is consistent with the conclusion of Berkowitz and Scher [11]. However, to build a quantitative relationship between the two properties and model parameters, additional laboratory and numerical experiments are needed.(3)Diffusioncontrolled subdiffusion is possible. In heterogeneous or even homogeneous sand columns, subdiffusive transport due to molecular diffusion occurs even for a large Peclet number. The relatively immobile regions formed during soil repacking cause diffusioncontrolled subdiffusion, following the physical process of multirate mass transfer. The diffusioncontrolled subdiffusion can be apparent in undisturbed soils, where the immobile zones are almost inevitable and the corresponding mass transfer rate varies significantly in space.(4)There are serious limitations in applying the slowadvectiondominated subdiffusive model, such as model (4), to capture realworld subdiffusion due to mass exchange. In the slow advectiondominated subdiffusive model, the transient anomalous mass decline cannot be captured. The bestfit velocity differs from the measurement, and the time scale index has a wider range and represents different physical processes. These limitations cause high uncertainty in predicting nonFickian transport.
Acknowledgments
This work was supported by the Desert Research Institute (DRI) and the National Science Foundation (NSF) under Grant no. DMS1025417. Part of the experiments was conducted by Wallace Atterberry supported by the first author. This paper does not necessary reflect the views of NSF or DRI. The authors thank three anonymous reviewers for helpful comments that improved the paper.
References
 Y. Zhang, D. A. Benson, and D. M. Reeves, “Time and space nonlocalities underlying fractionalderivative models: distinction and literature review of field applications,” Advances in Water Resources, vol. 32, no. 4, pp. 561–581, 2009. View at: Publisher Site  Google Scholar
 R. Metzler and J. Klafter, “The random walk's guide to anomalous diffusion: a fractional dynamics approach,” Physics Reports, vol. 339, no. 1, pp. 1–77, 2000. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Metzler and J. Klafter, “The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics,” Journal of Physics A, vol. 37, no. 31, pp. R161–R208, 2004. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Schumer, D. A. Benson, M. M. Meerschaert, and B. Baeumer, “Fractal mobile/immobile solute transport,” Water Resources Research, vol. 39, no. 10, Article ID W1296, 2003. View at: Publisher Site  Google Scholar
 M. M. Meerschaert, Y. Zhang, and B. Baeumer, “Tempered anomalous diffusion in heterogeneous systems,” Geophysical Research Letters, vol. 35, no. 17, Article ID L17403, 2008. View at: Publisher Site  Google Scholar
 S. P. Neuman and D. M. Tartakovsky, “Perspective on theories of nonFickian transport in heterogeneous media,” Advances in Water Resources, vol. 32, no. 5, pp. 670–680, 2009. View at: Publisher Site  Google Scholar
 M. Levy and B. Berkowitz, “Measurement and analysis of nonFickian dispersion in heterogeneous porous media,” Journal of Contaminant Hydrology, vol. 64, no. 34, pp. 203–226, 2003. View at: Publisher Site  Google Scholar
 M. Bromly and C. Hinz, “NonFickian transport in homogeneous unsaturated repacked sand,” Water Resources Research, vol. 40, no. 7, Article ID W07402, 2004. View at: Publisher Site  Google Scholar
 B. Berkowitz, A. Cortis, M. Dentz, and H. Scher, “Modeling Nonfickian transport in geological formations as a continuous time random walk,” Reviews of Geophysics, vol. 44, no. 2, Article ID RG2003, 2006. View at: Publisher Site  Google Scholar
 R. Haggerty, S. A. McKenna, and L. C. Meigs, “On the latetime behavior of tracer test breakthrough curves,” Water Resources Research, vol. 36, no. 12, pp. 3467–3479, 2000. View at: Publisher Site  Google Scholar
 B. Berkowitz and H. Scher, “Exploring the nature of nonFickian transport in laboratory experiments,” Advances in Water Resources, vol. 32, no. 5, pp. 750–755, 2009. View at: Publisher Site  Google Scholar
 Y. Zhang, M. M. Meerschaert, and B. Baeumer, “Particle tracking for timefractional diffusion,” Physical Review E, vol. 78, no. 3, Article ID 036705. View at: Publisher Site  Google Scholar
 C. Harvey and S. M. Gorelick, “Ratelimited mass transfer or macrodispersion: which dominates plume evolution at the Macrodispersion Experiment (MADE) site?” Water Resources Research, vol. 36, no. 3, pp. 637–650, 2000. View at: Publisher Site  Google Scholar
 E. E. Adams and L. W. Gelhar, “Field study of dispersion in a heterogeneous aquifer, 2. Spatial moments analysis,” Water Resources Research, vol. 28, no. 12, pp. 3293–3307, 1992. View at: Publisher Site  Google Scholar
 B. Baeumer and M. M. Meerschaert, “Tempered stable Lévy motion and transient superdiffusion,” Journal of Computational and Applied Mathematics, vol. 233, no. 10, pp. 2438–2448, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 Y. Zhang and C. Papelis, “Particletracking simulation of fractional diffusionreaction processes,” Physical Review E, vol. 84, no. 6, Article ID 066704, 2011. View at: Publisher Site  Google Scholar
 M. M. Meerschaert, Y. Zhang, and B. Baeumer, “Particle tracking for fractional diffusion with two time scales,” Computers & Mathematics with Applications, vol. 59, no. 3, pp. 1078–1086, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
Copyright
Copyright © 2013 Yong Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.