Abstract
Heterogeneous media consisting of segregated flow regions are fractional-order systems, where the regional-scale anomalous diffusion can be described by the fractional derivative model (FDM). The standard FDM, however, first, cannot characterize the Darcy-scale 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. Late-time 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 late-time 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 Darcy-scale 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 integer-order counterpart [2, 3]. In addition, the sorption-desorption mechanism can cause non-Markovian 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 regional-scale heterogeneous porous and fractured media for more than a decade [4], its applicability for Darcy-scale dispersion remains obscure. Indeed, some studies implied that the standard FDM may not be applicable for small-scale 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 Darcy-scale.
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 fractional-order partial differential equation on Darcy-scale 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 non-Fickian diffusion for passive tracer transport through repacked laboratory columns of macroscopically homogeneous sand. Curve-fitting applications further show that the non-Fickian diffusion characterized by the nonsymmetric plume cannot be explained efficiently by the 2nd-order advection-dispersion equation (ADE) model based on Fickian diffusion [6, 9]. Well-designed 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 Darcy-scale 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 derivative-based, nonlocal transport model is developed and compared with the other nonlocal models. In Section 4, the proposed model is used to capture the observed Darcy-scale 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 five-point 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 non-Fickian 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 advection-dominated 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 late-time 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 late-time 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 late-time 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 late-time 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 late-time 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 late-time 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 non-Fickian 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 laboratory-scale 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 non-Fickian 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 cut-off 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 late-time 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 advection-related subdiffusion.
3.2. The TFDM for Diffusion-Dominated 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 diffusion-dominated 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 power-law 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 Riemann-Liouville 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 regional-scale transport (see, e.g., Adams and Gelhar [14]). These observations lead to a specific question: will any small-scale 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 late-time 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 best-fit 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 best-fit parameters cm2/min, min−0.01, min−1, , and . The standard FDM, however, slightly overestimates the late-time 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 best-fit parameters, including cm2/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 second-order ADE (shown by the dashed line in Figure 4(a)) underestimates the BTC late-time tail.
For Run 2 using 0.4 + 0.2 mm glass beads, the best-fit parameters for the BTC with the slowest flow velocity (3.89 cm/min) are cm2/min, min−0.02, min−1, , and . Model prediction, however, overestimates the late-time 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 best-fit 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 late-time 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 late-time 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 Darcy-scale, 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 best-fit 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 real-world soils; clearly, simulating the real-world fast motion paths using repacked sand is difficult, if not impossible. In contrast, the real-world 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 Real-World 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 advection-dominated 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 best-fit parameters are cm/min, cm2/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 late-time tail of the BTC (with a slope ~−6.5 in a log-log plot). In model (2a) and (2b), the steep late-time 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 late-time 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 late-time tails, while model (4) has limited capability due to the lack of the controlling parameter . In addition, the best-fit 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 real-world 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 late-time 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 diffusion-related 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 fractional-order 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 non-Fickian motion of reactant molecules (the main reason why the diffusion-limited anomalous kinetics deviate significantly from the thermodynamic law) may be efficiently simulated by the TFDM (2a) and (2b) with a particle-based 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 variable-order or the distributed-order 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 Darcy-scale 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 Darcy-scale 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 Darcy-scale. 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 late-time BTC, probably due to the finite distribution of particle waiting times, while the standard FDM tends to overestimate the late-time tail of BTC.(2)Both the sand particle size distribution and the fluid velocity can affect the Darcy-scale subdiffusion. All of the measured BTCs of bromide contain an apparent late-time 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 late-time 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)Diffusion-controlled 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 diffusion-controlled subdiffusion, following the physical process of multirate mass transfer. The diffusion-controlled 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 slow-advection-dominated subdiffusive model, such as model (4), to capture real-world subdiffusion due to mass exchange. In the slow advection-dominated subdiffusive model, the transient anomalous mass decline cannot be captured. The best-fit 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 non-Fickian transport.
Acknowledgments
This work was supported by the Desert Research Institute (DRI) and the National Science Foundation (NSF) under Grant no. DMS-1025417. 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.