Diffusion in Relatively Homogeneous Sand Columns : A Scale-Dependent or Scale-Independent Process ?

Solute transport through homogeneous media has long been assumed to be scale-independent and can be quantified by the second-order advection-dispersion equation (ADE). This study, however, observed the opposite in the laboratory, where transport of CuSO4 through relatively homogeneous silica-sand columns exhibits sub-diffusion growing with the spatial scale. Only at a very small travel distance (approximately 10 cm) and a relatively short temporal scale can the transport be approximated by normal diffusion. This is also the only spatiotemporal scale where the fundamental concept of the “representative element volume” (which defines the scale of homogeneous cells used by the ADE-based hydrologic models) is valid. The failure of the standard ADE motivated us to apply a tempered-stable, fractional advection-dispersion equation (TS-FADE) to capture the transient anomalous dispersion with exponentially truncated power-law late-time tails in CuSO4 breakthrough curves. Results show that the tempering parameter in the TS-FADE model generally decreases with an increase of the column length (probably due to the higher probability of long retention processes), while the time index (which is a non-local parameter) remains stable for the uniformly packed columns. Transport in sand columns filled with relatively homogeneous silica sand, therefore, is scale-dependent, and the resultant transient sub-diffusion can be quantified by the TS-FADE model. Entropy 2013, 15 4377


Introduction
Fractional-derivative models have been increasingly used in the last two decades to quantify anomalous diffusion in disordered systems [1,2].The entropy theory was also combined with fractional calculus in analyzing fractional dynamics with a power-law memory [3].For example, El-Wakil and Zahran [4] applied the maximum entropy principle to reveal the structure of the probability distribution function of waiting times underlying the fractional Fokker-Planck equation.Machado [5] analyzed multi-particle systems with fractional order behavior.Ubriaco [3] applied Fisher's information theory based on his new definition of fractional entropy [6] to derive mathematical models for anomalous diffusion.Machado [7] found a power law evolution of the system energy and the entropy measures in fractional dynamical systems filled with colliding particles.Here, we focus on fractional porous media equations, which may also be related to the entropy solution theory [8].
Porous media focused on by the hydrologic community are known as complex dynamic systems containing multi-scale heterogeneity, but the application of the fractional engine is limited [9], due to probably two reasons.First, the second-order advection-dispersion equation (ADE) is believed by the hydrologic community to be valid.Most hydrologic numerical models are grid-based, where each grid is homogeneous.Transport within the grid is assumed to be scale-independent normal diffusion, which can be quantified by the ADE.The available, detailed subsurface heterogeneity can then be embedded in the model using a combination of millions of grids.Therefore, although contaminant transport through heterogeneous porous media from pore to regional scales has been well documented to be non-Fickian (as characterized by heavy tails in tracer breakthrough curves (BTC) [10,11]), the second-order ADE model with a certain resolution of velocity remains the routine modeling tool in the field of hydrology.For readers interested in this topic, we refer to the most recent comments and replies on the feasibility of the ADE-based models [12][13][14][15].Second, the fractional-derivative models do describe anomalous diffusion more efficiently than the standard ADE [9], but they also introduce additional parameters (such as the fractional order), whose linkage with medium properties may remain obscure.
Critical questions however remain for hydrologic numerical models that have been used for decades.First, does the transport through typical homogeneous porous media remain as scale-independent normal diffusion?If not, then the classical ADE-based modeling tool is questionable.This leads to the subsequent question: could the fractional-derivative model be the appropriate alternative to the ADE for transport in the deceptively simple homogeneous media?In other words, does the diffusion through homogeneous media exhibit fractional dynamics with a power-law memory in either space or time or both?Finally, for practical applications, what are the main factors affecting the fractional-derivative model parameters?These critical questions motivated this study.
The rest of the paper is organized as follows.In Section 2, we introduce the systematic laboratory transport experiments focusing on the dynamics of nonreactive tracers moving through sand columns of various lengths.Relatively uniform silica sand was used to fill the columns, forming an ideal porous medium that is much more homogeneous than the natural geological medium.In Section 3, we quantify the observed dynamics using both the classical ADE and different fractional-derivative models.The questions raised above are then discussed in Section 4. Finally, conclusions are drawn in Section 5.

Experimental Setup
We packed glass tubes (with an inner diameter of 25 mm) using silica sand with a relatively uniform size.The corresponding median grain size packed in each tube was 0.73 mm (i.e., coarse sand), 0.35 mm (medium sand), and 0.21 mm (fine sand), respectively.The silica sand was soaked in nitric acid for 24 h, and then washed with tap water and deionized (DI) water.After drying in an oven, the sand is ready for packing.Finally, the glass columns were packed using the wet sand loading method (which can minimize the air bubbles [16]), and the resultant porosity was measured (Table 1).To test the scale effect, we built sand columns with three different lengths of 10 cm, 20 cm and 40 cm.
Table 1.Parameters used in the models.R denotes the median grain size.v and D ADE are the velocity and dispersion coefficient used in the advection-dispersion (ADE) model in Equation (1), respectively.γ is the order of the time fractional derivative in the fractional-derivative models in Equations ( 3) and (5).λ is the tempering parameter, r is the scale factor and D FADE is the dispersion coefficient in the tempered-stable, fractional advection-dispersion (TS-FADE) model in Equation (5).The transport experiment involved three main steps.First, DI water with a pH of 2.0 was run through the column (oriented vertically) for a period of ten pore volumes, and then, the background solution (i.e., tap water, in this case) was run through for another five pore volumes to build the flow domain and remove the background (concentration) effect.The vertical flow is from top to bottom.A peristaltic pump (BT100-1F, LongerPump) was used to regulate the downward flow at a specific discharge around ∼ 1 mL/min.Second, the CuSO 4 solution was added into the column continuously for 40 min at a concentration of 0.5 mmol/L, followed by three pore volumes of water for flushing.Third, discrete samples were collected from the outlet using a fraction collector (BS-100A, PuYang Scientific Instrument Research Institute, Nanjing, China).
Finally, we measured the sample concentration by: (1) adding 100 uL nitric acid; (2) diluting the solution to a volume of 5 mL using DI water; (3) passing the solution through a 0.45 µm moisture film; (4) measuring the absorbance using an atomic absorption spectrophotometer (Z-2000, Hitachi, Tokyo, Japan); and (5) converting the absorbance to the concentration.

Experimental Results
Figure 1.The CuSO 4 breakthrough curves (BTC) along a 10 cm-long sand column: the measurements (symbols) versus the best-fit solutions using the ADE model in Equation (1) (grey lines), the time-fractional advection-dispersion equation (tFADE) model in Equation (3) (the dashed line) and the TS-FADE model in Equation ( 5) (black lines).In the legend, "R" denotes the median grain size.The measured CuSO 4 BTCs are shown in Figures 1-3, for the travel distance of 10 cm, 20 cm and 40 cm, respectively.The early time tails of all BTCs (representing the early arrivals of solute) are as steep as exponential, implying that there is no fast movement along preferential flow paths.This is expected, since super-diffusive transport typically requires a heterogeneous medium with a hydraulic conductivity field exhibiting large correlation length and variance [9].The late-time tails of BTCs, however, become relatively heavier with an increase of the travel distance.In the next section, we conduct numerical analysis to reveal whether the observed transport is actually scale-dependent.

Quantifying the Observed Transport
We first apply the local transport model (i.e., the standard ADE model) to simulate the observed CuSO 4 BTCs.If the ADE fails, we will then apply the non-local transport models (i.e., the fractional-derivative models) and compare them with the ADE model.

The ADE Model
The second-order ADE model takes the form: where is the macroscopic dispersion coefficient used by the ADE model (which is not necessarily the same as the other dispersion coefficients used below) and s [M L −3 T −1 ] is the source/sink term.The ADE in Equation ( 1) can be solved analytically [17].Here, we also apply the well-known Lagrangian solver (see [18,19], among many others) to solve Equation (1).The Lagrangian solver is selected here, since it can be extended conveniently to approximate all the fractional-derivative models used below.The space and time Markov processes underlying in Equation ( 1) are: where X i [L] denotes the particle position at time T i , dt i [T ] is the operational time used by the i-th jump and ω (dimensionless) is a normally distributed random variable with zero mean and unit variance.Here, the physical/clock time increases linearly with the operational time.
The above Lagrangian scheme was validated extensively against analytical solutions (for simplicity, they are not shown here).Solutions of the ADE model in Equation ( 1) fit generally well with the measured BTCs for the shortest sand column (Figure 1), but underestimate significantly the late-time tail of the other BTCs (Figures 2 and 3).

The Time-Fractional Advection-Dispersion Equation (tFADE) Model
The time-fractional advection-dispersion equation model can capture the heavy tail of tracer BTCs [9], since it describes the heavy tailed memory in time.Using the subordination tool [20], the time-fractional advection-dispersion equation (tFADE) model can be written as [21]: where the fractional time derivative of the order of γ (0 < γ < 1) (dimensionless) (i.e., the scale index) and the scale factor, r [t γ−1 ], describe an inverse stable distribution of clock times between jump events.Here, the dispersion coefficient, D FADE , can be different from (i.e., smaller than) D ADE in Equation ( 1), since the tFADE model in Equation ( 3) accounts for the variation of transport velocity using the fractional-time derivative.The Caputo fractional derivative for time is used in this study.Similar to Equation (2), the Lagrangian solver for the tFADE model in Equation (3) contains the following two Markov processes, after using the extended Langevin approach [22]: where dL γ (dimensionless) represents a stable random variable with the maximum skewness (β = +1), unit scale (σ = 1) and zero mean shift (µ = 0).Note that the operational time, dt i , is now different from the clock time.
The above Lagrangian solver in Equation ( 4) is validated, with several examples shown in Figure 4, where the Lagrangian solutions generally match the fast Fourier transform solutions (see, also [20]) of Equation (3).
Further applications, however, show that, as opposed to the ADE model in Equation ( 1), the tFADE model in Equation (3) overestimates significantly the late-time tail of the BTCs (see Figures 1a, 2a and 3a).The tFADE model assumes a power-law distribution for the particle clock time, which leads to a late-time BTC tail heavier than the measurement.

The Tempered-Stable, Fractional Advection-Dispersion Equation (TS-FADE) Model
Truncated stable Lévy flights were first proposed by Mantegna and Stanley [24] to censor arbitrarily large jumps and capture the natural cutoff present in real physical systems.Exponentially-tempered stable processes were proposed by various researchers [25][26][27][28] as a smoother alternative, without a sharp cutoff.The tempered stable density may describe the distribution of the random clock time between jump events.
To capture the truncated power-law decline of the late-time BTC, we apply the following tempered-stable, fractional advection-dispersion equation (TS-FADE): where λ [T −1 ] is the tempering parameter.The operator, ∂ γ,λ /∂t γ,λ , denotes [29]: The Lagrangian solver developed for the tFADE model in Equation ( 3) can be used for the TS-FADE model in Equation (5), where the only change is the stable random variable, dL γ , in Equation ( 4) replaced by a tempered stable random variable, dL γ,λ .
Applications show that the TS-FADE model in Equation ( 5) can capture the observed BTCs for all cases (Figures 1-3).The best-fit model parameters are shown in Table 1.The linear-linear scale of Figures 1-3 shows that some solutions of the ADE model in Equation ( 1) are almost identical to those of the TS-FADE in Equation ( 5), while the discrepancy in the simulated BTC tails between the two models is amplified by the logarithmic scale.This is not a surprise, given that the mass in the late-time BTC tail (such as Figure 2a) represents only a small fraction (0.78%) of the total mass.
In the next section, we discuss the possible underlying mechanism and the main controlling factors for the observed transport dynamics.

Short-Duration of Normal Diffusion in Relatively Homogeneous Porous Media
The above laboratory and numerical tests reveal that normal diffusion may only exist for a short travel distance (i.e., ∼10 cm) in the relatively homogeneous sand columns.The short-duration of normal diffusion may be due to two reasons.First, there may be small-scale variations in the packing of the sands, which leads to micro-structure (such as aggregates) in the macroscopic homogeneous medium.The silica sand used in our experiments is not perfectly uniform, but has a relatively narrow size distribution, which also helps to build internal structures.Solute particles diffusing into the sand matrix or a dead-end water zone may be delayed and, therefore, build the late-time tail of the BTC.If the transport is a non-dissipative process, the microscopic scale heterogeneity may control the macroscopic dynamics [30].With an increase of the spatial scale, more (and perhaps larger) aggregates may be formed, and the solute transport is delayed further, resulting in scale-dependent sub-diffusion.It is well-known that the solute particle jumps can be regarded as instantaneous [1], and hence, the clock time [expressed, for example, by Equation (4b)] between jump events represents the random waiting time for random-walking particles.At a small scale in both space and time, most particles have not experienced large retention periods yet, and the transport exhibits initial behavior similar to normal diffusion.Normal diffusion, therefore, may only be an approximation of real-world anomalous diffusion at a small spatiotemporal scale where the anomaly has not apparently developed yet.
Another explanation is the fractal geometry of silica sand (see, for example, [31,32], among others), which tends to generate anomalous diffusion.This explanation, however, is difficult to validate directly.Some researchers also argued that the fractal properties of soil might not be so obvious [30].The qualitative link between fractal properties of sand (such as texture and surface area) and fractional dynamics remains to be shown [33].In addition, a recent study [34] found that uniform glass beads (which may contain aggregates or relatively immobile flow zones when they are packed in glass tubes) without any multi-fractality can also lead to sub-diffusion.
It is noteworthy that the limited duration of normal diffusion may be even shorter in natural geological formations with intrinsic multi-scale heterogeneity.Mixing and structured sands in the field can enhance (i.e., super-diffusion) or decelerate (i.e., sub-diffusion) the motion of solute particles, which can appear much earlier than those in the laboratory sand columns.
The short-duration of normal diffusion challenges the ADE-based hydrologic modeling, with one example discussed blow.

Challenge on the Definition of the Representative Element Volume
The scale-limited normal diffusion constrains seriously the size of the representative elementary volume (REV) [35].The conventional local ADE is believed to be applicable at the scale of REV, so that the heterogeneity of a large-scale medium must be adequately represented at the REV scale.For a regional scale (i.e., kilometers) model, the size of each homogeneous cell is usually larger than one meter, which is at least one order of magnitude larger than the valid scale of normal diffusion revealed by this study.A finer resolution with a small REV, however, can lead to a prohibitive computational burden.
Several recent studies also identified a very small REV or even could not find the scale for REV.For example, Yoon and McKenna [36] found that the REV may exist at the length of 0.25 cm, while local-heterogeneity features below the REV should still be quantified in numerical modeling.Klise et al. [37] conducted an unprecedented laboratory experiment by taking a thin slab of Massillon sandstone and exhaustively sampling the permeability (k) via air permeability sampling.The 30.5 × 30.5 × 2.1 cm slab was measured for k every 0.33 cm, yielding 17,328 measurements.Each sample support volume was on the order of 0.45 cm 3 .The finely discretized ADE, however, could not capture the observed early or late time tails of the tracer BTC [37].Major et al. [38] further found that sub-grid dispersion (below the support volume) is non-Fickian, and the non-local transport model is needed to capture the observed transport.

Fractional Dynamics for Tracer Transport in Relatively Homogeneous Sand Columns
The underlying dynamics for tracer transport in relatively homogeneous sand columns is transient sub-diffusion, due to probably the physical and chemical properties in the transport process.The finite retention capacity of sand matrix, probably due to the limited thickness of matrix and the non-negligible diffusive displacement of solute [39], acts as an upper bound for tracer particle waiting times.This is similar to what we observed in solute retention in alluvial aquifers [39].Fractional dynamics in porous media therefore may depend on both the physical properties of the media and the chemical properties of the tracer.
The waiting time for solute particles between jump events exhibits multi-fractal scaling, which evolves in space.Numerical fitting in Section 3.3 shows that the particle waiting time distribution at a given spatial scale is a power-law function transferring gradually to become exponential.According to [40], multi-fractality can arise from a linear, additive process, whose increments have power-law tails with a variable truncation.Note that, here, the multi-fractal waiting times also grow in space, since the tempering parameter, λ (which defines the rate of exponential tempering of the power-law tail), changes with the travel distance.

Factors Affecting Sub-Diffusion in Relatively Homogeneous Sand Columns
In our experiments, the spatial scale affects the tempering parameter, λ, but the time fractional order, γ (which is also the BTC or waiting-time tail index), remains stable in these uniformly packed sand columns.For a given tracer, the index, γ, is a non-local parameter characterizing the overall retention capacity of the medium.In general, a stationary medium with strong immobile regions may be characterized by a small, constant γ [9].The tempering parameter, λ, on the other hand, records the extreme retention event, which, therefore, may change with the local variation of immobile zone properties.
In addition, the sand size also significantly affects the tailing behavior of transport, which can also be captured by adjusting the tempering parameter in the TS-FADE model in Equation (5).For coarse and medium sand, λ decreases with an increase of the travel distance (see Table 1).For fine sand, however, λ fluctuates with the spatial scale, where the relative amplitude of fluctuation is smaller than that for the coarser sand.The discrepancy might imply that the immobile zones formed by fine sand have relatively less variability (such as properties of segregates) than that for coarser sands.This hypothesis needs further experimental and numerical tests in a future study.
It is noteworthy that, for a short travel distance in relatively homogeneous media, the parameters in the TS-FADE model in Equation ( 5) may not be unique when fitting the measured BTC.This is because the late-time tailing of transport cannot fully develop in a limited time; see, for example, Figure 1c.Even if the observation time is long enough to capture the full range behavior of late-time tailing, the measured BTC may still remain incomplete, due to the concentration detection limit of tracers.Caution is therefore needed when quantifying a small-scale transport using the TS-FADE model in Equation (5).

Possible Influence of the Small Diameter of the Sand Columns on Anomalous Diffusion
The repacked sand columns in this study have a relatively small diameter (25 mm), which may affect solute transport in three ways.First, a narrower column may force sand to be packed tightly, resulting in more dead zones for flow that can enhance the trapping of solutes.Second, solute particles may be trapped between the narrow glass tube and the filled sands, and solutes may also sorb on to the glass tube.Third, a narrower sand column provides less spatial inter-connectivity that is necessary for super-diffusion.The first two impacts tend to enhance sub-diffusion, while the last one may constrain the generation of super-diffusion.
Further analysis, however, shows that the above possible impacts might be minimized in this study.For example, for a short column (Figure 1), the classical ADE model captures the observed BTCs, showing that the narrow sand column with a short length does not necessarily cause sub-diffusion.In addition, although none of the observed BTCs exhibit heavy early-time tails (the BTC with a power-law early-time tail is one of the typical characteristics of super-diffusion), this does not mean that the sand column with a small diameter must constrain super-diffusion.Previous studies, such as Herrick et al. [41] and Kohlbecker et al. [42], showed that the heavy-tailed and long-range-dependent hydraulic conductivity (K) field is needed to generate heavy-tailed solute displacement.In this study, the repacked silica-sand columns do not contain such a heterogeneous K field.Indeed, to the best of our knowledge, no previous studies showed that the column experiments packed with relatively homogeneous sand can build heavy-tailed super-diffusive transport.Therefore, the small diameter of the sand columns may have limited impact on the anomalous transport observed in this study.We will check this hypothesis in a future study, by using glass tubes with various inner diameters.Additional factors, such as variable-density flow and undisturbed soils may also be considered to explore the possible influence of the lateral dimension of sand columns on tracer transport.
It is also noteworthy that glass tubes with a centimeter-scale inner diameter have also been used to study various aspects of transport, especially in the last two years.For example, Ngueleu et al. [43] used a glass column with an inner diameter of 24 mm (and a length of 150 mm) to minimize the sorption of lindane onto the equipment.Sagee et al. [44] conducted column experiments for silver nanoparticle transport in closed, cylindrical columns with an inner diameter of 31 mm (and a length of 100 mm).Gouet-Kaplan et al. [45] packed sands in an acrylic vertical column with an internal diameter of 30 mm (to a height of 150 mm) to study the mixture of water.Zhang et al. [34] measured bromide transport through a horizontal glass tube with internal dimensions of 15.9 mm (diameter) (and 150 mm in length).Historic and well-known column experiments also used columns with similar sizes.For example, Gramling et al. [46] monitored bimolecular reactions through a thin glass chamber with the smallest thickness being 18 mm.Raje and Kapoor [47] measured reactive kinetics across a glass column with a diameter of 45 mm (and a length of 180 mm).We also emphasize that the focus of the above column experiments differs from this study.

Reason to Select the Fractional Models with Temporal Derivatives
The introduction of the time-fractional derivative in Equations ( 3) and ( 5) is motivated by the observed late time tailing of BTCs.Other time-nonlocal transport models can also capture the delayed transport, including the well-known multi-rate mass transport (MRMT) model [48] and the continuous time random walk (CTRW) framework [10], which have been used widely by the hydrologic community.The time-fractional derivative models in Equations ( 3) and ( 5), which can describe the time nonlocal transport processes, such as sub-diffusion, are specific and simplified MRMT models with power-law distributed mass exchange rates.Model in Equation (5) may also be functionally equivalent to the CTRW framework with a truncated power-law memory function in Equation [10], and model in Equation ( 5) requires fewer parameters (e.g., only the tempering parameter, λ) to capture the nuance of an exponentially truncated power-law tail in the BTC.Discrepancy between fractional models and the other time-nonlocal models is not the main focus of this study.We leave this discussion for future work.
The transport process observed in this study is not super-diffusive, but sub-diffusive, since the late time tailing of BTC suggests a retardation process (a typical behavior of sub-diffusion), and there is no sign of fast displacement for CuSO 4 .This process might not be related obviously to the small diameter of the sand columns, since the sand column with a much larger diameter (i.e., 300 mm) can also lead to sub-diffusion (see, for example, [49]), and the above discussion implies the minimal impact of column diameter on transport.
It is noteworthy that we did not use the fractional models with spatial derivatives, since they describe processes different from our observations.Particularly, the space fractional advection-dispersion equation with maximally positive skewness captures super-diffusion with heavy power-law early-time tails in BTCs [9], which is not observed in our column experiments.The space fractional advection-dispersion equation with maximally negative skewness does capture sub-diffusion, but the solute particles must travel backward and reach the upstream boundary [9].This power-law backward movement is not apparent, if not impossible, in our laboratory tests.

Conclusions
This study evaluates the dynamics of nonreactive tracer transport in relatively homogeneous media.The fundamental assumption in typical hydrologic models is that normal diffusion in homogeneous cells is scale-independent.To check this assumption, we conducted laboratory transport experiments and explored whether the dynamics of CuSO 4 transport through silica sand columns varies with the travel distance.The measured BTCs were then interpreted using both the standard ADE model and the fractional-derivative models.The combined study of laboratory tests and stochastic analysis leads to the following three major conclusions.
(1) Normal diffusion and the representative element volume are only valid at small scales.
(2) The TS-FADE model can capture the scale-dependent sub-diffusive transport through relatively homogeneous media.
(3) The tempering parameter in the TS-FADE model generally decreases with an increase of the column length (due probably to the higher probability of long retention precesses), while the time index (which is a non-local parameter) remains stable.

Figure 2 .
Figure 2. The CuSO 4 BTC along a 20 cm-long sand column: the measurements (symbols) versus the best-fit solutions using the ADE model in Equation (1) (grey lines), the tFADE model in Equation (3) [the dashed lines in (a) and (d)] and the TS-FADE model in Equation (5) (black lines).

Figure 3 .
Figure 3.The CuSO 4 BTC along a 40 cm-long sand column: the measurements (symbols) versus the best-fit solutions using the ADE model in Equation (1) (grey lines), the tFADE model in Equation (3) (the dashed lines in (a) and (d)) and the TS-FADE model in Equation (5) (black lines).

Figure 4 .
Figure 4.The Lagrangian solutions (symbols, denoted as "RW" (representing Random Walk) in the legend) of the tFADE model in Equation (3) versus the the fast Fourier transform (FFT) solutions (lines).(b) is the log-log plot of (a).An instantaneous point source is injected at x = 0.The velocity v = 1, dispersion coefficient D = 0 and the control plane is located at x = 10.In the legend, C R denotes the resident concentration and C F denotes the flux concentration[23].