Application of continuous time random walk theory to nonequilibrium transport in soil

Journal of Contaminant Hydrology 108 (2009) 134–151

Contents lists available at ScienceDirect

Journal of Contaminant Hydrology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j c o n h yd

Application of continuous time random walk theory to nonequilibrium transport in soil Na Li, Li Ren ⁎ Department of Soil and Water Sciences, China Agricultural University, Key Laboratory of Plant–Soil Interactions, MOE, Beijing 100193, China

a r t i c l e

i n f o

Article history: Received 14 March 2009 Received in revised form 26 June 2009 Accepted 2 July 2009 Available online 12 July 2009 Keywords: Continuous time random walk (CTRW) Sorbing solute Breakthrough curve (BTC) Nonequilibrium transport Numerical experiment Prediction model

a b s t r a c t Continuous time random walk (CTRW) formulations have been demonstrated to provide a general and effective approach that quantifies the behavior of solute transport in heterogeneous media in field, laboratory, and numerical experiments. In this paper we first apply the CTRW approach to describe the sorbing solute transport in soils under chemical (or) and physical nonequilibrium conditions by curve-fitting. Results show that the theoretical solutions are in a good agreement with the experimental measurements. In case that CTRW parameters cannot be determined directly or easily, an alternative method is then proposed for estimating such parameters independently of the breakthrough curve data to be simulated. We conduct numerical experiments with artificial data sets generated by the HYDRUS-1D model for a wide range of pore water velocities (υ) and retardation factors (R) to investigate the relationship between CTRW parameters for a sorbing solute and these two quantities (υ, R) that can be directly measured in independent experiments. A series of best-fitting regression equations are then developed from the artificial data sets, which can be easily used as an estimation or prediction model to assess the transport of sorbing solutes under steady flow conditions through soil. Several literature data sets of pesticides are used to validate these relationships. The results show reasonable performance in most cases, thus indicating that our method could provide an alternative way to effectively predict sorbing solute transport in soils. While the regression relationships presented are obtained under certain flow and sorption conditions, the methodology of our study is general and may be extended to predict solute transport in soils under different flow and sorption conditions. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Quantification of chemical transport in soil/groundwater has generated considerable interest owing to concern over ground water quality. Characterizing the transport behavior of contaminants requires an understanding of physical, chemical, biological processes that influence the fate of contaminants in the subsurface (e.g., Brusseau et al., 1991b). Time-related or, nonequilibrium, sorption of organic chemicals by soils has received a significant degree of attention for quite some time since it is one of the most important processes that determine

⁎ Corresponding author. Department of Soil and Water Sciences, China Agricultural University, Beijing 100193, China. Tel.: +86 10 62731675; fax: +86 10 62733596. E-mail address: [email protected] (L. Ren). 0169-7722/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2009.07.002

the rate of movement, persistence and removal of contaminant in the environment (e.g., Brusseau and Rao, 1989b; Brusseau et al., 1991a,b; Gamerdinger et al., 1991). With respect to sorbing solute transport, asymmetric and long-tail breakthrough curves (BTCs) have been frequently observed for a variety of chemicals in a broad range of soils (e.g., Rao et al., 1979; Gamerdinger et al., 1991; Nakayama et al., 1994; Gaber et al., 1995). Such non-ideal phenomena, which have often been ascribed to multiple sources of nonequilibrium including physical nonequilibrium and sorption nonequilibrium (e.g., Brusseau et al., 1989; Brusseau and Rao, 1989a,b; Brusseau and Rao, 1990; Brusseau et al., 1991a,b), indicate a need to provide qualitative as well as quantitative understanding of nonequilibrium transport by both experiments and mathematical models. Solute transport in soils is most commonly described by the (Fickian-based) advection–dispersion equation (ADE),

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

wherein the sorption is represented by a retardation factor (Mojid and Vereecken, 2005), which can be expressed as

R

AC ðx; t Þ A2 C ðx; t Þ AC ðx; t Þ =D −υ At Ax A2 x

ð1Þ

where C(x,t) is the solution-phase solute concentration, t is time, x is distance, D is the dispersion coefficient, υ is the pore water velocity, and R is the retardation factor. For linear sorption, R = 1 + (ρbKd/θυ), where θυ is the volumetric soil water content, ρb is the soil bulk density, Kd is the distribution coefficient. We denote the average velocity of the sorbing solute by ῡ = υ/R. The ratio λ = D/υ is generally referred to as the dispersivity (Bear, 1972), which together with velocity completes the identification of the two parameters appearing in Eq. (1). The classical ADE or alternative equations have been widely used as models of chemical transport in soil– water systems for many years and have provided valuable information through their use. In real natural systems, however, transport often exhibits nonequilibrium behavior, which cannot be properly quantified by using this classical ADE (Gamerdinger et al., 1991; Cortis and Berkowitz, 2004). Improvements and extensions of the ADE have been proposed and evaluated for describing nonequilibrium transport. The two-region and the two-site transport models (van Genuchten and Wierenga, 1976; Selim et al., 1976; Cameron and Klute, 1977) have been developed for describing the sorption-related nonequilibrium (SNE) and transport-related nonequilibrium (TNE). A variety of analytical and numerical solutions for the two-site and the two-region models, which have been included into optimization packages (e.g., CFITIM, van Genuchten, 1981; CXTFIT, Parker and van Genuchten, 1984; HYDRUS-1D, Šimunek et al., 1998), are available for interpretation of chemical transport at laboratory scale or field scale. These two models have been successfully applied for quantitative studies of nonequilibrium transport (e.g., Gamerdinger et al., 1990, 1991; Gaber et al., 1995; Mao and Ren, 2004). The equivalence of the two-region and the twosite transport models has been discussed by several authors (Nkedi-Kizza et al., 1984; Selim and Amacher, 1988; van Genuchten and Wagenet, 1989). In order to further understand the chemical nonequilibrium, this two-site model was later extended to the second-order two-site (SOTS) models including SOTS-CDE and SOTS-MIM (coupled physical–chemical nonequilibrium) models (Selim and Ma, 1995), which have proven to be successful in describing atrazine retention and transport in soils (e.g., Ma and Selim, 1994a,b; Selim et al., 1999; Ma and Selim, 2005). In fact, the nature of the contaminant mobility in soil– water system is strongly influenced by the soil heterogeneity (e.g., Seuntjens et al., 2002). This physical and chemical heterogeneity existing at all scales is assumed responsible for the appearance of the early breakthrough times and long time tailing in measured BTCs (e.g., Berkowitz et al., 2001), which is often referred to as “anomalous” or “non-Fickian” transport. Recently, a general, physically based approach based on the continuous time random walk (CTRW) formalism has been introduced to account naturally for the anomalous behavior in the context of heterogeneous media by Berkowitz and Scher (1995, 1997, 1998) and Scher et al. (2002). In the CTRW framework, contaminant migration is approximated as a

135

series of transitions characterized by ψ(s,t), the probability density function of a local displacement s at transition time t. This transition probability density function incorporates all possible mechanisms governing transport such as advection, diffusion, dispersion, sorption and so on (Margolin et al., 2003; Berkowitz et al., 2006). The key physical ideal of CTRW theory is that anomalous dispersion is typically due to a wide distribution (a power law distribution) of delay times, ψ(s,t) (Berkowitz and Scher, 1995). Based on the CTRW formation, the theoretical framework enables one to determine the evolution of the solute plume, P(s,t), given a general ψ(s,t) (Berkowitz and Scher, 1998). Such a modeling approach provides an effective, upscaled transport framework based on the local-averaged behavior (Margolin et al., 2003), presenting a different understanding and insight into the solute transport behavior. To obtain more complete characterization of the transport, recently, Berkowitz et al. (2008) developed a physical model that incorporates two generic mechanisms into the CTRW framework: the complex flow field of a highly heterogeneous medium and the mass exchange between a mobile phase and a distribution of immobile states. A twotime domain form of ψ(s,t) is used to distinguish between the two different mechanisms, in which additional CTRW parameters are involved. We note, parenthetically, that a variety of fractal mobile/immobile models (MIM) (e.g., Benson et al., 2000; Metzler and Klafter, 2000; Schumer et al., 2003) and multirate mass transfer (MRMT) formulations (e.g., Pfister and Scher, 1978; Roth and Jury, 1993; Harvey and Gorelick, 1995; Haggerty and Gorelick, 1995; Haggerty et al., 2000; Carrera et al., 1998) also incorporated a spectrum of transition times, which have been considered and applied to data analyses (see review by Berkowitz et al., 2006). The relationship between CTRW and these two models has been addressed recently (Berkowitz et al., 2002; Schumer, 2002; Dentz and Berkowitz, 2003; Margolin et al., 2003; Berkowitz et al., 2006). With regard to other approaches Berkowitz et al. (2006) presented several analyses that demonstrated the ADE, fractional advection–dispersion models, and multirate mass transport and mobile–immobile models to be, in fact, specialized subsets within the CTRW framework. The CTRW theory was first applied to problems in solidstate physics such as electron movement in disorder semiconductors (e.g., Montroll and Scher, 1973; Scher and Montroll, 1975) and later introduced in the context of geological materials (Berkowitz and Scher, 1995). An initial hydrological application of the CTRW theory was to solute transport simulations in fractured and strongly heterogeneous porous media based on discrete fracture networks (Berkowitz and Scher, 1997, 1998). The plausibility of this approach to heterogeneous porous media was then demonstrated by successful application to laboratory and field experiments. For instance, it was applied to describe anomalous tracer transport in small-scale laboratory models (Berkowitz et al., 2000), to quantify a broad range of Fickian and non-Fickian transport behavior of BTCs measured in midscale laboratory flow cells (Levy and Berkowitz, 2003), and to analyze field-scale, non-Fickian transport in fractured till (Kosakowski et al., 2001). More applications to some classical experimental data demonstrated the ability of CTRW-based solutions to accurately capture the anomalous early arrival times and late time tailing even existing in the homogeneous

136

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

soil columns that cannot be properly quantified by using the ADE (Berkowitz and Scher, 2001; Cortis and Berkowitz, 2004; Jiménez-Hornero et al., 2005; Berkowitz et al., 2006). Such applications to the experimental measurements show that the solutions of the CTRW turn out valid for a broad range of non-Fickian and Fickian transport behaviors of conservative chemicals in heterogeneous porous media. Nevertheless, the transport behavior of a sorbing chemical in porous media is distinctly different from that of a nonsorbing one. So far, attempts to describe sorbing solute transport in soils using the general CTRW framework (e.g., Berkowitz et al., 2002; Cortis et al., 2004; Berkowitz et al., 2006) are sparse (e.g., Deng et al., 2008; Berkowitz et al., 2008). Deng et al. (2008) applied the general CTRW framework with a single ψ(s,t) to sorbing solute transport and demonstrated that the CTRW convincingly captured and predicted the full evolution of atrazine BTCs in a repacked soil column, especially with respect to the long tailing that the two-site model failed to characterize. The two-time domain form of CTRW model proposed by Berkowitz et al. (2008) produces excellent fits to both tracer transport in a fractured shear zone and sorbing species transport through a heterogeneous porous domain. Here we focus on the CTRW framework with a single ψ(s,t) to evaluate its effectiveness for describing sorbing solute transport undergoing chemical or (and) physical nonequilibrium. This approach is justified because application of a single ψ(s,t) has been highly successful in understanding the dynamics of a host of laboratory and field observations (Berkowitz et al., 2008), and more importantly, because the purpose of the present manuscript is not to separate between the two mechanisms. Moreover, it is preferable to choose a model with a fewer number of parameters for the predictive purpose presented below. First, we apply this form of CTRW with a single ψ(s,t) to several published data sets through the commonly used curve-fitting method. Within this approach, both the effects of dispersion and sorption are incorporated into the central function ψ(s,t) (the key function of the CTRW framework). Therefore, their effects cannot be separated. Then we focus on our work to propose an estimation or prediction method to determine the CTRW parameters independently of the BTC data to be simulated. Among the work involved when using mathematical models to describe solute transport behavior, the major one is to determine the input parameters. Such parameters can commonly be obtained from measured BTCs by data-fitting. However, in circumstances where BTCs cannot be measured directly or easily, such method is not applicable. For example, in practice, due to the restrictions on the use of some reactive contaminants especially radioactive materials, experiments could not be conducted under laboratory or field conditions. Therefore, an alternative approach for obtaining the input parameters instead of curve-fitting, i.e., independently of the BTC data to be simulated is expected. With respect to sorption, a common approach is to combine the results of passive tracer tests with laboratory-derived sorption parameters to predict the behavior of sorbing solutes. Successful predictions that result from this approach have been reported by several authors (Brusseau et al., 1989; Brusseau, 1992; Brusseau and Srivastava, 1997). Motivated by some ideas mentioned above, we present an alternative way based on

numerical experiments with artificially generated BTCs to predict CTRW parameters of a sorbing solute by relating them to some measurable parameters (viz. the pore water velocities, v, which can be quantified by passive tracer tests and retardation factors, R, which can be determined by batch experiments). The HYDRUS-1D package is selected to generate the required BTC “data” since it has been widely used and proven to be a reasonable tool for simulation of transport. We note parenthetically, that investigation of the relationship between CTRW parameters and the flow rate is not new (e.g., Levy and Berkowitz, 2003; Berkowitz and Scher, 2009). These authors addressed such issue by analyzing several laboratory experimental sets of BTC data in terms of CTRW models. Here, we conduct mathematical experiments in a series of imaginary soil media of different permeability and adsorption capacities to investigate relationships not only between CTRW parameters and the pore water velocity but also between CTRW parameters and the retardation factor. Our goal here is to predict sorbing solute transports in soils by direct measurements of parameters in independent experiments or determination from existing empirical correlations instead of curve-fitting when no data can be given. The main advantage of our approach is that all input parameter values to execute the CTRW model can be estimated independently of the data to be simulated. The primary objectives of this study are (i) to evaluate CTRW approach for describing sorbing solute transport in soils under chemical or/and physical nonequilibrium conditions, (ii) to estimate the CTRW parameters by regressing them to some measured parameters, so that the sorbing solute transport behavior can be predicted when the behavior of a passive tracer in the same system has been quantified, and when the equilibrium sorption isotherms of the sorbing chemical have been measured. For the first purpose, laboratory data obtained from miscible displacement experiments reported in published papers are used to evaluate the CTRW approach. For the second purpose, a series of regression relationships based on the numerical experiments with artificially generated BTC data are established as a prediction model to estimate the CTRW parameters for a sorbing solute. Finally, to illustrate their use and assess their validity, we apply such regression relationships to several published data sets. 2. Continuous time random walk (CTRW) theory In the CTRW framework, contaminant migration in a variable velocity field is envisioned as particles executing a series of steps, or transitions, through the formation via different paths with spatially changing velocities (Cortis and Berkowitz, 2004). Generally, this kind of transitions can be characterized by a joint probability density function, ψ(s,t). Being the central function in the CTRW framework, it incorporates all possible mechanisms governing transport such as advection, diffusion, sorption, and so on, so that wide variability in single transition distances and times can lead to wide temporal tails in the breakthrough curves (BTCs) (Margolin et al., 2003). Here we present a short summary associated with the CTRW theory. Comprehensive explanations and a full derivation of the method have been presented elsewhere (Berkowitz and Scher,1998; Margolin and Berkowitz,

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

2000; Berkowitz et al., 2002; Dentz and Berkowitz, 2003; Cortis et al., 2004; Berkowitz et al., 2006). 2.1. Basic equation We present here the basic equation associated with the CTRW approach. The CTRW is a generalization of the random walk from discrete time to continuous time with the spatial state description remaining discrete Rðs; t Þ =

X Z 0

s

t 0

ψðs − s0; t − t 0ÞRðs0; t 0Þdt 0

ð2Þ

where R(s,t) is the probability per time for a particle to just reach s at time t, and ψ(s,t) is the probability density for a transition between sites (i.e. the displacement is s) with a difference of arrival times of t. Just as Cortis et al. (2004) described, here, we will focus on the decoupled transition length and time ψ(s,t): ψ(s,t) ≈ p(s)ψ(t), where p(s) denotes the transition length probability density and ψ(t) denotes the transit-time (or waiting time) distribution defined below in (3). ψðt Þu

X

ψðs; t Þ

ð3Þ

s

It is important to note that the transport Eq. (2) represents a different point of departure from the usual formulations because the advective, dispersive, and diffusive transport mechanisms are inextricably combined in (2) rather than being treated separately (Berkowitz et al., 2000). The function ψ(t) is the “heart” of the CTRW formation, dominating the principal characteristics of solute plume migration patterns (Berkowitz et al., 2001). Interesting anomalous transport may be obtained with slow-decaying ψ(t) with a diverging standard deviation, particularly a power law form (e.g., Montroll and Scher, 1973; Hatano and Hatano, 1998) ψðt Þet

−1−β

; 0 b β b 2:

ð4Þ

The relation between the power law and anomalous transport has been documented by Berkowitz et al. (2001). The background of a power law distribution of waiting times can have either a physical or chemical nature, reflecting the heterogeneity of the medium (Gheorghiu and Coppens, 2004). Many choices of the transition time distribution ψ(t) have emerged from many studies employing different methodologies and a functional representation that uses a minimum parameters is the truncated power law (TPL)(Cortis et al., 2004; Berkowitz and Scher, 2009) h    i − 1 −β −1 −1 ψðt Þ = t1 τ2 exp τ 2 C −β; τ 2

expð−t = t2 Þ ;0 b β b 2 ð1 + t = t1 Þ1 + β

ð5Þ which requires three input parameters (β, t1, and t2), and τ2 = t2/t1. The t1 parameter represents a median time for transitions between sites, and t2 is a “cutoff” time. For t1 ≪ t ≪ t2, ψ(t) ∝ t− 1− β (power law region), while for t ≫ t2 the transport behavior is Gaussian. The key factor is the interplay between β and the cutoff time t2, which has a dramatic effect on the entire shape of a migrating solute

137

plume (Berkowitz and Scher, 2009). The TPL is a complete form representing the features of the power low region and the “cutoff” of this region, allowing a transition from nonFickian behavior to Fickian behavior at longer times. The nature of non-Fickian transport derives from the interplay between the span of observation times and the width of the power low region. The ψ(t) exhibits a very wide change in behavior as a function of β, which is reported by Cortis et al. (2004). The exponent β controls the particle migration behavior and thus functionally quantifies the dispersion behavior, which characterizes the nature of the transport behavior (Fickian or nonFickian). The relative shapes of the anomalous transport curves, and the rate of advance of the peak, vary strongly as a function of β. If β ⩾ 2, the first two moments of the ψ(t) exist and the solute plume will display Gaussian behavior. For 1 b β b 2, the second moment of the ψ(t) does not exist, while for 0 b β b 1, even the first moment does not exist. 1 b β b 2 is associated to moderately dispersive systems, while 0 b β b 1 indicates a highly anomalous transport. The exponent dominates the trapping and transition time distribution and considerably influences the transport behavior. It may be expected that the value of β varies as a function of the relative degree of heterogeneity (e.g., Berkowitz et al., 2001). The issue of the dependence of β on the flow rate has been addressed by analysis of experimental BTCs with the CTRW model (Levy and Berkowitz, 2003; Berkowitz and Scher, 2009). Levy and Berkowitz (2003) performed three sets of experiments using a “homogeneous” packing, a randomly heterogeneous packing, and an exponentially correlated structure. Breakthrough experiments were reproducible for three different flow rates with each of the porous medium systems. Best fit analysis using a pure power law tail ψ(t) ~ t− 1− β indicated the same trend for each of the experiment sets, namely that the key parameter β slowly increased as the flow rate decreased. Berkowitz and Scher (2009) reanalyzed the BTC data with the full spectrum of ψ(t) and showed explicitly that changing the velocity does not result in different CTRW parameter values. The use of the full spectrum of ψ(t) is not only necessary for the transition to Fickian behavior, but also to account for the dynamics of these laboratory observations of non-Fickian transport. More details about the effects of β on the transport behavior are well discussed elsewhere (e.g., Margolin and Berkowitz, 2000; Dentz et al., 2004; Cortis et al., 2004). 2.2. Transport equation In the CTRW framework, the governing equation is the Fokker–Planck with memory equation (FPME), which is derived from the equivalence of the classic Generalized Master Equation (GME) and CTRW (Berkowitz et al., 2002; Cortis et al., 2004). The Laplace transform of the FPME for one dimension is given by " # 2 ˜ ðuÞ υ A c˜ ðs; uÞ − D A c˜ðs; uÞ uc˜ ðs; uÞ − c0 ðsÞ = −M ð6Þ ψ ψ As As2 where ˜ ðuÞ = tu M

˜ ðuÞ ψ ˜ ðuÞ 1−ψ

ð7Þ

138

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

is a memory function that accounts for the unknown heterogeneities below the level of measurement resolution, and υψ =

1X pðsÞs nt s

ð8Þ

Dψ =

1X 1 pðsÞss nt s 2

ð9Þ

"

where n is the porosity, t ̄ is the characteristic time, υψ is the transport velocity, and Dψ is the generalized dispersion coefficient. Note, both the “effective velocity” υψ and the dispersion Dψ are time-dependent, and most significantly, depend fundamentally on ψ(t). The transport velocity υψ is generally different from the pore water velocity υ in the original form of the ADE model, although they have the same dimensions. These velocities are identical when Eq. (6) reduces to the classical ADE. To obtain the normalized υψ and Dψ, the values are simply divided by L and L2 respectively, in which L is the column length, so that υψ = υψ/L, and Dψ = Dψ/L2. The dimensionless ‘dispersivity’ can then be defined as P 1 s pðsÞss αψ = P 2 s pðsÞs

ð10Þ

The flux-averaged concentration is defined as j ̃ ≡ M̃υψ (c̃ − αψ∂c ̃/∂s) (Cortis and Berkowitz, 2004). Given a unit steady flux boundary condition j ̃ = u− 1 at the inlet (s = 0) and a natural boundary condition ∂C̃/∂s = 0 at the outlet (s = 1) we can obtain the analytic solution of Eq. (6) for j ̃(u) at s = 1 j˜ðuÞ =

h i 2zexp 12 α −1 ψ + z n h  i h  io ˜ ˜ + z − α −1 u expðzÞ z + α −1 ψ + 2u = Mυψ ψ − 2u = Mυψ

ð11Þ where vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 4uα −1 1 u ψ t1 + z= ˜ αψ Mυ ψ

used here are digitized from the published literature except the data from Mao and Ren (2004). Different combinations of adjustable parameters are compared to select the one that best fits the experimental BTC. As a measure for the error between observed data and the calculated values the root mean square error (RMSE) is used

ð12Þ

Analytical solutions to the Eq. (6) for a number of common boundary conditions (Dentz et al., 2004) are explicit in Laplace space and can be easily transformed to the time domain, numerically, by means of the de Hoog algorithm (de Hoog et al., 1982). A MATLAB toolbox (version 3.0) with a collection of subroutines that solve these equations (Cortis and Berkowitz, 2005) is freely available from the website http://www.weizmann.ac.il/ESER/People/Brian/CTRW. 3. Fitting published data sets using CTRW models Based on the analyses in the previous sections, it might be expected that CTRW theory can elucidate sorbing solute transport from a different aspect. Since curve-fitting is a valid procedure for model verification (Rao et al., 1979), we now present the comparison of BTCs from the CTRW solutions to experimental data of pesticide transport in soils. All data sets

RMSE =

N  2 1X Cif −Cim N i=1

#1 = 2 ð13Þ

where Cif is the fitted concentration, and Cim is the measured data of the concentration at the corresponding time point. RMSE indicates the scatterness around 1:1 line when plotting Cif against Cim. A larger RMSE means a more scattered relationship, whereas a smaller one means more appropriate. 3.1. Data set 1 Data from Mao and Ren (2004) were taken as an example of the application of the CTRW models for describing nonequilibrium transport through saturated homogeneous soil columns when the nonequilibrium conditions were assumed to be resulted from two-site sorption. Atrazine, one of the widely used herbicide, was selected for batch sorption studies and column displacement experiments that were performed with a sandy loamy soil. The experimental procedures were discussed earlier by Mao and Ren (2004). Miscible displacement experiments were conducted in saturated homogeneous soil columns at a specific pore water velocity in duplicate. Columns, 20 cm in length and 4.984 cm (for the first column) and 4.920 cm (for the second column) in inner diameter, were used in the experiments. The columns were packed vertically in incremental steps and each increment was trapped firmly. KBr was used as a nonsorbing, conservative tracer to characterize the hydrodynamic properties of the packed soil columns. A pulse input of an aqueous solution containing atrazine and KBr was applied to the top of the soil columns at a specific pore water velocity (1.9255 cm h− 1 for the first column and 1.9642 cm h− 1 for the second column) for an hour. The flow rate of the injected solution was controlled by a peristaltic pump (P-1, Pharmacia Fine Chemicals, Sweden) to ensure that steady-state flow conditions were satisfied during the displacement experiments. Fig. 1 shows the comparison of effluent concentration distributions for atrazine and Br− through the first column (υ = 1.9255 cm h− 1). The measured BTC for tracer Br− was predominantly symmetrical, indicating there were no significant effects of physical nonequilibrium during transport. Also shown in Fig. 1 are the fitted BTCs using the CTRW solutions. Clearly, the solutions result in excellent fits to both atrazine and Br− BTCs, which capture the full evolution of plume migration. The corresponding fitted parameters of CTRW-based TPL model and the goodness of fits in terms of RMSE are listed in Table 1. To tracer Br−, as can be seen, the parameter β → 2 for the TPL form of ψ̃(u). Such result means that the evolution of the tracer Br− displays Fickian behavior. Unlike Br− BTC, in the experiment conducted with atrazine through the same column (υ = 1.9255 cm h− 1) the BTC was asymmetrical and exhibited significant tailing, as was reported previously (e.g., Gamerdinger et al., 1990; Gaber

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

139

Fig. 1. Measured breakthrough curves of Br− and atrazine in the first column (after Mao and Ren (2004)). Column length equals 20 cm. Pore water velocity υ = 1.9255 cm h− 1.

Fig. 2. Prediction (the parameters determined from the first experiment) of atrazine data of the second column (after Mao and Ren (2004)) using the TPL model. Column length equals 20 cm. Pore water velocity υ = 1.9642 cm h− 1.

et al., 1995). Given the fact that Br− BTC was not influenced by physical nonequilibrium, nonequilibrium conditions exhibited by the atrazine BTC were primarily related to sorption nonequilibrium. Table 1 shows that the fitted value of β for atrazine is 0.3011, which is very different from that of Br−. Such value of β is typical of highly dispersive (non-Fickian) transport according to Cortis et al. (2004). A comparison between the pore water velocity (v = 1.9255 cm h− 1) and the value of υ ψ = 1.0877 cm h− 1 illustrates the difference between pore water velocity and chemical velocity. Conversely, the optimized values of the CTRW parameters are used to predict the BTC at υ = 1.9642 h− 1 treatment in the second column. Reasonable prediction is obtained and shown in Fig. 2.

packed with aggregates of ≤6 mm for experiments 3-3 and 4-1. Sorption processes in both the stagnate and dynamic regions of the medium were assumed to be instantaneous (i.e., with no chemical nonequilibrium effects). Here, we analyze the measured BTCs by obtaining best fit curves using solutions based on the CTRW solutions. As before, the CTRW solutions completely capture the BTC behavior, as is seen in Fig. 3. Parameter values for the CTRW model are summarized in Table 2. For pore water velocity of 0.399, 0.4501, 1.491, and 1.535 cm h− 1, the β = 0.9784, 0.7907, 0.2778, 0.2305, respectively. Note that the value of β decreases with increasing pore water velocity, as is reported previously when modeling dispersion behavior of conservative tracer (Levy and Berkowitz, 2003). With respect to a conservative tracer, changing the flow rate in a given domain leads to fundamentally different transport behaviors; this can be attributed to the change in residence time of solute within the domain which in turn changes the degree of averaging of the system disorder (Berkowitz and Scher, 2009). With respect to a sorbing solute, more conditions may influence the residence time. Thus the results, namely, the slowly decreasing values of β with increasing pore water velocities, may be valid within some specific ranges of pore velocity.

3.2. Data set 2 The second group of data sets utilized were obtained by van Genuchten et al. (1977) from miscible displacement of the herbicide (2,4,5-T) through unsaturated, aggregated sorbing porous media when the nonequilibrium conditions were assumed to be resulted from the presence of immobile water rather than two-site sorption. 2,4,5-T was applied to uniformly packed soil columns of 30-cm long at four pore water velocities (υ = 0.450 cm h− 1 for experiment 1-4, υ = 0.399 cm h− 1 for experiment 2-4, υ = 1.491 cm h− 1 for experiment 3-3, υ = 1.535 cm h− 1 for experiment 4-1). The columns used in experiments 1-4 and 2-4 were uniformly packed with aggregates of ≤2 mm, while the columns were

3.3. Data set 3 Data sets from Gaber et al. (1995) were used as additional cases for evaluating the application of the CTRW theory to describe sorbing transport through intact soil columns when

Table 1 Transport parameters for atrazine and tracer BTCs fitted by the CTRW-based TPL model. The accuracy of the model is evaluated by the root mean square errors (RMSE). Fitted parameters (TPL) Chemical

υa

Rb

cm h− 1 Br− Atrazine a b

1.9255 1.9255

– 1.8725

Value from the tracer BTC. Value from the batch experiments.

υψ

Dψ

cm h− 1

cm2 h− 1

2.0533 1.0877

0.1832 0.2227

β

1.9335 0.3011

lg t1

lg t2

lg h

lg h

− 2.1297 0.6210

6.0326 × 10− 2 0.9342

RMSE 5.2518 × 10− 3 1.5853 × 10− 3

140

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

Fig. 3. Measured data (after van Genuchten et al. (1977)) and fitted breakthrough curves for atrazine using the TPL solutions. υ = 0.450 cm h− 1 for experiment 1-4, υ = 0.399 cm h− 1 for experiment 2-4, υ = 1.491 cm h− 1 for experiment 3-3, and υ = 1.535 cm h− 1 for experiment 4-1.

both transport-related nonequilibrium (TNE) and sorptionrelated nonequilibrium (SNE) are operative. They performed a series of miscible displacement experiments under varying pore water velocities and soil–water contents in intact soil columns (30-cm length; 15.24-cm diam.). A vacuum chamber providing constant suction of 13.3 kPa was used to establish a constant hydraulic gradient. Three water Darcy flux levels (q) of 0.04, 0.27, and 0.86 cm h− 1 were established at a constant θυ of approximately 0.39 cm3 cm− 3. Correspondingly, the pore water velocities (υ) of 0.12, 0.69, and 2.16 cm h− 1 were provided. The medium pore water velocity experiment was repeated (0.74 cm h− 1) at a higher θυ of 0.44 cm3 cm− 3. A

solution pulse of about 0.12 pore volumes containing 3H2O, analytical-grade atrazine and ring-labeled 14C atrazine in 3 mM CaCl2 matrix was applied to each column. Tritiated water (3H2O) was employed as a conservative tracer to determine the dispersion characteristics of the packed soil columns and the existence of TNE. Column-effluent data were collected and analyzed for each species to generate BTCs for 3 H2O and atrazine. In their study of tracer transport, Gaber et al. (1995) demonstrated that the tritiated water BTCs were nearly symmetrical and exhibited only minor tailing at the slow (0.12 cm h− 1) and the medium (0.69 cm h− 1) pore water

Table 2 Summary of the fitted parameters of the CTRW-based TPL model in the simulations of the experimental breakthrough curves of atrazine reported by van Genuchten et al. (1977). Fitted parameters (TPL) Exp. no.

υa

Rb

cm h− 1 1-4 2-4 3-3 4-1 a b

0.4501 0.399 1.491 1.535

2.225 2.21 2.139 2.223

Value from the tracer BTC. Value from the batch experiments.

υψ

Dψ

cm h− 1

cm2 h− 1

0.2319 0.188 0.914 0.8054

0.6582 0.8291 1.5921 2.5066

β

0.7907 0.9784 0.2778 0.2305

lg t1

lg t2

lg h

lg h

0.214 0.0936 1.6096 0.0751

0.6169 0.5449 1.9105 0.2758

RMSE 3.9986 × 10− 2 2.7311 × 10− 2 1.4201 × 10− 2 2.5222 × 10− 2

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

141

Table 3 Summary of the fitted parameters of the TPL solutions in the simulations of the experimental breakthrough curves of tracer (3H2O) reported by Gaber et al. (1995). Fitted parameters (TPL) Treatment

Fast Medium (θυ = 0.44) Medium Slow a

υa

υψ

Dψ

cm h− 1

cm h− 1

cm2 h− 1

2.16 0.74 0.69 0.12

1.9182 0.5232 0.6369 0.09821

19.196 0.5935 0.2006 0.1692

β

1.0919 1.7996 1.8290 1.6138

lg t1

lg t2

lg h

lg h

2.0356 0.5370 0.1632 0.2093

4.4235 1.3785 4.0157 1.2022

RMSE 4.1835 × 10− 3 1.6393 × 10− 3 6.1277 × 10− 3 3.6144 × 10− 3

Value from the tracer BTC.

velocities, but exhibited significant asymmetry and tailing at the highest pore water velocity (2.16 cm h− 1). The CXTFIT analysis conducted by them suggested that TNE processes were operative at the high velocity and some TNE was exhibited at the medium pore water velocity (0.74 cm h− 1) with a higher volumetric soil water content (θυ = 0.44 cm3 cm− 3) experiment. Here we analyze the tracer data in terms of the CTRW framework with the TPL model. A summary of the fitted parameters values for all of the 3H2O BTCs is given in Table 3. Values for the RMSE are also listed in this table. As can be seen in Fig. 4, the CTRW provides excellent fits to all of four BTCs, for appropriate choices of β. The values of β for the high pore water

velocity is much smaller than those for the medium and the slow pore water velocities; for pore water velocities of 2.16, 0.74, 0.69, 0.12 cm h− 1, we get values of β = 1.0919, 1.7996, 1.8290, 1.6138, respectively. Such results indicate that the transport is moderate non-Fickian at high pore water velocities, while a slightly non-Fickian transport behavior exists at slow and medium pore water velocities. This can be explained by the conclusion of Gaber et al. (1995), which stated that high pore water velocity, high volumetric water content, or both enhanced the degree of preferential flow or the fraction of immobile water in intact columns and thus contributed to the TNE observed under high and medium pore water velocities at

Fig. 4. Measured data (after the Fig. 2 of Gaber et al. (1995)) and fitted breakthrough curves for the tracer (3H2O) using the TPL solutions: (a) fast, (b) medium (θυ = 0.44 cm3 cm− 3), (c) medium, and (d) slow pore water velocity.

142

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

Fig. 5. Measured data (after the Fig. 3 of Gaber et al. (1995)) and fitted breakthrough curves for atrazine using the TPL solutions: (a) fast, (b) medium (θυ = 0.44 cm3 cm− 3), (c) medium, and (d) slow pore water velocity.

influenced primarily by SNE at the slow (0.12 cm h− 1) and the medium (0.69 cm h− 1) pore velocities, and a combination of both TNE and SNE at pore water velocities of 0.74 (θυ ≈ 0.44 cm3 cm− 3) and 2.16 cm h− 1. We now present TPL model fits to the measured BTCs of atrazine. Once again, the TPL model captures the full evolution of the BTCs completely, as can be seen in Fig. 5. The fitted parameters for the TPL model as well as the values of RMSE are listed in Table 4. As is shown, the β values of TPL model range from 0.8 to 1.5 for the atrazine BTCs, which is well below the range in

higher θυ. Moreover, similar results—the decreasing values of β with increasing pore water velocity—are found from medium to fast velocities, as discussed in section 3.2. However, the β value for the slow pore water velocity (0.12 cm h− 1) treatment does not follow this pattern. We speculate that this may be explained that the lower pore water velocity in a specific range would enhance the degree of physical nonequilibrium. As was stated by Gaber et al. (1995), the BTCs for atrazine were asymmetrical and exhibited significant tailing at all of the pore water velocities, indicating that atrazine was

Table 4 Summary of the fitted parameters of the TPL solutions in the simulations of the experimental breakthrough curves of atrazine reported by Gaber et al. (1995). Fitted parameters (TPL) Treatment

υa

Rb

cm h− 1 Fast Medium (θυ = 0.44) Medium Slow a b

2.16 0.74 0.69 0.12

Value from the tracer BTC. Value from the batch experiments.

3.74 3.49 3.82 3.85

υψ

Dψ

cm h− 1

cm2 h− 1

0.9741 0.2732 0.24 0.0316

14.591 1.0098 0.4272 0.0817

β

0.8086 1.3497 1.4656 0.8205

lg t1

lg t2

lg h

lg h

0.8720 2.4867 2.4655 1.8298

1.2017 4.9876 10.5410 2.0170

RMSE 7.5623 × 10− 4 7.6166 × 10− 4 1.2695 × 10− 3 1.3448 × 10− 3

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

which Fickian transport occurs (β N 2); the values of β = 0.8086, 1.3497, 1.4656, 0.8205 for pore water velocities of 2.16, 0.74, 0.69, 0.12 cm h− 1. For the medium and the fast pore water velocities, the higher the velocity is, the lower the value of β will be. This is similar to the analysis of tracer experiments and may be explained that an increase in the pore water velocity would result in a greater degree of SNE given a constant sorption rate (Gaber et al., 1995), which leads to a broader distribution trapping time of migrating particles and a smaller value of β. On the other hand, the β value for the slow pore water velocity (0.12 cm h− 1) treatment does not follow this pattern, which is similar to the slow pore water velocity treatment of the conservative tracer 3H2O. This again might be caused by the deeper degree of TNE. The successful fits again suggest that the CTRW model provides a good description of the sorptive transport under physical or (and) chemical nonequilibrium conditions. The application above suggests that the key parameter β of the CTRW should be a function of both physical and chemical nonequilibrium conditions. In addition, for all of the four treatments, we notice that the β for atrazine is less than that of tracer 3H2O for the same column. From the analyses of the three independent sets of experiments above, we demonstrate that the CTRW framework provides excellent fits to all of the experimental BTCs. The physical meaning of the CTRW parameters and the relationship between characteristics of media and specific parameters of the CTRW theory have been investigated previously (e.g., Berkowitz et al., 2001), which stated that the β parameter characterizes the functional nature of the dispersion and such parameter may be expected to vary as a function of the relative degree of heterogeneity in fractured and heterogeneous porous media. In our application here, all of the mechanisms including both dispersion and sorption are incorporated in the central function (ψ(t) ~ t− 1− β) of this CTRW framework, whose effects cannot be separated. Thus the single parameter β quantifies all of the mechanisms that control the transport behavior. We assume that β varies as a function of the sorption and flow properties, which will be illustrated by the numerical experiments presented below in section 4. 4. Numerical experiments Our numerical experiments aim to build the regression relationships between the CTRW parameters and two measurable quantities (i.e., pore water velocity, υ, and retardation factor, R) using artificial breakthrough curve (BTC) data. For this purpose, (i) for a range of pore water velocities and retardation factors, we generate several artificial sets of BTC data (with different pore water velocities and retardation factors) by using the HYDRUS-1D model with assumption of nonequilibrium two-site sorption, (ii) we fit the generated BTCs using the CTRW-based TPL model so as to obtain the fitted parameters of TPL model for all of the generated BTCs, (iii) we regress the fitted parameters of TPL on pore water velocities (υ) and retardation factors (R) that are used for generating the artificial BTCs. After accomplishing these tasks, it is easy to use the regression relationships obtained to estimate the CTRW parameters for a sorbing solute and thus to predict its transport behavior, when v is quantified by the

143

behavior of a passive tracer in the same system and R is measured by sorption isotherms. 4.1. Determination of parameters to generate artificial BTCs Experimental evidence and theoretical results suggest that the transport of sorbing solute is influenced by the physical and chemical properties. Among these we limit our attention to two important parameters that govern the sorptive transport, i.e., the pore water velocity (υ) associated to the characteristic of the physical property and the retardation factor (R) representing the sorption. Particularly, these two parameters can be directly measured in independent experiments. Ten pore water velocity (υ) levels of 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, and 2.5 cm h− 1 are specified at a saturated water content θs of approximately 0.4 cm3 cm− 3. Under each velocity condition, gradually increasing equilibrium sorption constant (Kd) from 0.1 to 2 L mg− 1 (Kd = 0.1, 0.2, 0.3, …, 2 L mg− 1) are used to obtain increasing retardation factor R (=1 +ρbKd/θυ). The other parameters required to execute the HYDRUS-1D to generate the BTCs are assigned as follows. Soil bulk density is 1.33 g cm− 3. The depth of the soil profile (L) is 20 cm. A value of 2 cm is assumed for dispersivity λ, which is expressed by λ = 0.1L, based on the work of Lallemand-Barres and Peaudecerf (1978) (see, e.g., the reviews by Fetter, 1993). A Neumann boundary condition at the outlet is assumed for all the simulations. Here, we note that only chemical nonequilibrium is considered in the generation of BTC data sets. To generate the sorption BTCs which are primarily related to chemical nonequilibrium, it is necessary to determine some specific parameters including F (dimensionless fraction of adsorption sites with instantaneous sorption) and α (firstorder rate coefficient for two-site nonequilibrium adsorption) to execute the HYDRUS-1D model. Sensitivity analysis of the F revealed that the value of F does not influence model output significantly (e.g., Brusseau, 1992). A value of F (0.5) (reviewed by Brusseau et al., 1989) is used herein. Brusseau (1992) pointed that values for α can be determined by use of an empirical correlation equation relating the sorption rate constant (α) to the equilibrium sorption constant (Kd) instead of the individual laboratory experiments. Investigation of this empirical correlation has received a significant degree of attention (e.g., Brusseau and Rao, 1989a,b; Brusseau et al., 1989). Brusseau and Rao (1989b) compiled data from extensive literature on nonequilibrium sorption of organic solutes including a broad spectrum of nonpolar, hydrophobic organic chemicals (Type I) and polar/ionizable organic chemicals (Type II). An inverse log–log linear relationship was found between equilibrium sorption coefficient (Kd) and sorption rate constant (α) for each type. Here, since both atrazine and 2,4-D amine (selected for model validation) belong to Type II, the α–Kd relationship for the Type II log α = − 1:789 − 0:62 log Kd

ð14Þ

reported by Brusseau and Rao (1989b) is adopted to estimate the values for α. It is reported by Gamerdinger et al. (1990) that the value of α estimated with the two-site model is consistent with that predicted by this relationship. A summary of the input parameter values to generate the BTCs is given in Table 5.

144

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

Table 5 Parameter values used to generate the artificial BTC data in the numerical experiments. Parameter Definition

Value

L ρb θυ λ F

20 1.33 0.4 2 0.5

T0

Length of column, cm Soil bulk density, g cm− 3 Volumetric water content, cm3 cm− 3 Dispersivity, cm Dimensionless fraction of adsorption sites with instantaneous sorption Pulse size

1

4.2. New regression relationships First, for each pore water velocity from 0.25 to 2.5 cm h− 1, keep the chosen velocity fixed, then change the equilibrium sorption constant (Kd) from 0.1 to 2 L mg− 1 to vary the retardation factor. This gives a total of 200 BTCs artificially generated. Then, by fitting solutions of Eq. (6) to these artificial BTCs, we obtain five free parameters υψ, Dψ, β, t1 and t2 for each BTC. The comparison between the simulated BTCs and the best-fitting results reveals that each generated BTC can be nicely described by the TPL model. It is shown that all the fitted values of β are less than 2. Here we present a typical group of the generated BTCs for a range of Kd values (0.1 ≤ Kd ≤ 2 L mg− 1) at υ = 0.5 cm h− 1 in Fig. 6. Notice that the curves are asymmetrical and exhibit significant tailing at all Kd values. As is expected, the BTC becomes increasingly retarded and the C/C0 value at the breakthrough peak is reduced and shifted to longer time as the Kd increases. Also shown in Fig. 6 are the best fits to artificial BTCs using the CTRW solutions (solid line). Notice that the TPL model provides fairly good fits to the skewed BTCs. Fig. 7 shows the distributions of the fitted parameters (υψ, Dψ, β, t1, and t2) of TPL model across all (υ, R) treatments. For each pore water velocity varying from 0.25 to 2.5 cm h− 1, we plot the key parameter β against the retardation factor R (Fig. 7a). Note that the β continues to increase as R increases. We speculate that this is because the decreasing Kd causes corresponding increase in first-order rate coefficient α (recalling that Eq. (14)), which may result in increased left-handed displacement of the BTCs (i.e., increased nonequilibrium behavior) and thus result in decreased β. Hence, the direct relationship between β and R is obtained, i.e., β =e exp(fR) for υ = 0.25, 1.5, 1.75, 2, 2.25, 2.5 cm h− 1, and β =e ln R+f for υ = 0.5, 0.75, 1, 1.25 cm h− 1, where e, f are the regression coefficients. These regression equations for β and R are reported in Table 6. High coefficients of determination (r2 generally N0.9) for all the equations exhibit high correlation between β and R, indicating that the retardation factor have unnegligible influence on the key parameter β of the CTRW theory. We should also notice that the pore water velocity has significant influence on the key parameter β. The distribution of β across all pore water velocity conditions is plotted in Fig. 7b. As can be seen, the effects of pore water velocity on β are statistically significant. On average, the value of β decreases significantly when the pore water velocity increases except the υ =0.25 cm h− 1 treatments. Furthermore, to demonstrate the influence of the pore water velocity on BTC, we present another typical group of the BTCs for

a range of υ values (from 0.25 to 2.5 cm h− 1) at Kd = 0.1 L mg− 1 (R = 1.3023) in Fig. 8. Notice that increasing pore water velocities result in the increase in asymmetry and long-tailing displacement of the BTCs. Thus we may infer that increasing pore water velocities will result in decreasing value of β. Similar analyses are performed for υψ and Dψ. In Fig. 7c, for each pore water velocity varying from 0.25 to 2.5 cm h− 1, we plot the value of the parameter υψ against the retardation factor R. Notice that the υψ continues to decrease as R increases, where the relationship is log–log linear (lg υψ =a lg R+b, where a, b is the slope and intercept, respectively) at all of the pore velocities. The best approximations achieved for the slopes and intercepts of the regression lines are listed in Table 6. In Fig. 7d, we plot the value of parameter Dψ against the retardation factor R. Similar results—the continuously decreasing values of Dψ with increasing retardation factor (R)—are obtained, where the relationship is Dψ =cln R +d (where c, d are the regression coefficients listed in Table 6). Then a set of regression functions are computed to relate each regression coefficient (a, b, c, d, e, f) to the pore water velocity, which are also listed in Table 6. As is shown in Fig. 7e and f, increasing R results in a slight increased lg(t1) and lg(t2), respectively. As also can be seen, when R increases, the scatter range of both lg(t1) and lg(t2) mostly increases, especially for lg(t2). We also present the distribution of lg(t1) and lg(t2) with different pore velocities in Fig. 7g, h, respectively. The values of lg(t1) and lg(t2) decrease significantly when the pore water velocity increases from 1.25 to 1.5 cm h− 1. However, lg(t1) and lg(t2) are slightly affected for υ from 0.25 to 1.25 cm h− 1 and from 1.25 to 2.5 cm h− 1. Here, the discontinuous functions are used to represent the relationships between lg(t1) (lg(t2)) and υ. Summarizing, the CTRW parameters (β, t1, and t2) are significantly influenced by the pore water velocity, while are mildly influenced by the retardation factor. So far, the predictive model has been established. Once the pore water velocity has been quantified by passive tracer tests, and the equilibrium sorption isotherms of sorbing chemical have been measured by batch experiments, the regression

Fig. 6. Representative breakthrough curves (BTCs) generated by HYDRUS-1D. Dots: artificially generated BTC data with pore water velocity υ = 0.5 cm h− 1, Kd varying from 0.1 to 2 mg L− 1 (Kd = 0.1, 0.2, 0.3, …, 2 mg L− 1). Solid lines: best fits to artificial data using the TPL model.

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

145

Fig. 7. Relationships between TPL parameters and pore water velocity (υ) and the retardation factor (R). The data used are calculated from all artificially generated breakthrough curves. Edges of boxes represent 25 and 75% percentiles. Error bars represent 5 and 95% percentiles. The median (50% percentile) is represented by solid line in the box. The symbols (plus) refer to the maximum and minimum values.

results obtained above can be applied to predict the sorbing solute transport whose flow velocity and retardation factor are in the range of what has been used in the numerical

experiments to obtain the regression equations. Moreover, these relationships can help us to understand the effects of pore water velocity and retardation factor on the CTRW parameters.

146

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

Table 6 Regression results calculated from all artificially generated breakthrough curve data. Pore water velocity (υ)

υψ

Dψ

β

log υψ = a log R + b

Dψ = c log R + d

β = f(R) a

q/θυ

a

b

r2

c

d

r2

e

f

r2

0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5

− 1.0347 − 0.9332 − 0.9017 − 0.8794 − 0.8617 − 0.8730 − 0.8616 − 0.8505 − 0.8386 − 0.8324

− 0.5807 − 0.2761 − 0.1024 0.0230 0.1153 0.1323 0.2592 0.3432 0.3621 0.4378

0.9994 0.9986 0.9990 0.9982 0.9980 0.9980 0.9971 0.9966 0.9966 0.9968

− 0.0969 − 0.2628 − 0.3919 − 0.4711 − 0.9874 − 0.9355 − 0.9405 − 1.0117 − 1.2224 − 1.3966

0.3646 0.6358 0.9864 1.2656 2.2521 2.2755 2.4120 2.5991 3.1644 3.6152

0.8272 0.9506 0.9909 0.9767 0.9515 0.9930 0.9947 0.9907 0.9925 0.9961

0.1689 0.9020 0.8935 0.8545 0.6785 0.0486 0.0883 0.1968 0.0096 0.0089

0.2287 0.0624 − 0.1836 − 0.2843 − 0.1859 0.2862 0.1559 0.1431 0.2981 0.1980

0.9663 0.9260 0.9563 0.9210 0.9083 0.9447 0.9259 0.9430 0.9749 0.9812

a = 0.0787 ln(υ) − 0.8965(r 2 = 0.9263); b = 0.4331 ln(υ) + 0.0176(r 2 = 0.9945). c = − 0.5544υ − 0.0094(r 2 = 0.9312); d = 1.4177υ + 0.0077(r 2 = 0.9671).  − 0:67υ2 + 0:8887υ + 0:6198; 0:5 V υ V 1:25 e= 9:0315υ4 − 71:624υ 3 + 210:06υ2 − 269:95υ + 128:35; 1:5 V υ V 2:5 :  1:3776υ2 − 2:749υ + 1:0952; 0:5 V υ V 1:25 f = −5:0475υ4 + 38:393υ3 − 107:56υ 2 + 131:45υ − 58:904; 1:5 V υ V 2:5 : 8 < e lnðRÞ + f ; 0:5 V υ V 1:25 a β = f ðRÞ = : e expð fRÞ; else

4.3. Model evaluation using literature data To illustrate the application of these regression relationships, we consider several published experimental data sets. All of the selected experiments were conducted in packed columns during steady-state water flow involving different sorbing herbicides. Breakthrough behavior curves (BTCs) measured in these experiments are predicted in terms of TPL model based on the input parameters (υψ, Dψ, β, t1, and t2) calculated using the regression equations obtained above. Recall that the regression results can be applied to predict sorbing solute transport whose pore water velocity is in the range of 0.5 and 1.25 cm h− 1 or 1.5 and 2.5 cm h− 1. The first data set used to verify the regression relationships is from Mao and Ren (2004) discussed in section 3.1. The pore water velocity is calculated by the measured q and θs, where υ = q/θs. The R value is obtained from batch experiments, where R = 1 + ρKd/θυ. Independent predicted BTC (Fig. 9) is generated with the TPL model based on the input parameters calculated from the regression equations listed in Table 6. The use of such estimated parameter set results in a BTC shift to the left and lower peak concentration compared to the measured BTC although the general shape of the measured data is relatively matched. In fact, the prediction deviation in some degree is caused by overprediction of the average particle velocity υψ. However, from the view point of practical application, one may be satisfied with such prediction based on the model parameters from independent source. The second group of data sets is that of van Genuchten et al. (1977) discussed in section 3.2. Based on the range of pore water velocities, only the data set from experiment 4-1 (υ = 1.535 cm h− 1, aggregate diameter ≤ 6 mm) can be directly predicted by our regression results. The R value (2.223) was calculated with the “linearized” Kd value ob-

tained from a batch isotherm experiment. With the values of υ and R, the TPL model parameters are estimated by the regression equations listed in Table 6. As can be seen in Fig. 10a, the results indicate a good prediction for the measured data. The parameters for the experiment 4-1 are used to predict the BTC obtained from a different column experiment performed under similar conditions (experiment 3-3, υ = 1.491 cm h− 1, aggregated diameter ≤ 6 mm). Thus all parameter values for experiment 3-3 are obtained from independent source. The resultant independent prediction is presented in Fig. 10b, which provides a reasonable match to the measured data. Then, some parameters (β, t1, and t2) for the experiment 4-1 are also used to predict the BTC of experiment 1-4 (υ = 0.45 cm h− 1, aggregate diameter ≤ 2 mm), a system with a velocity and an aggregate size one third that of experiment 4-1. Note that all parameters except for the υψ and Dψ are obtained from experiment 4-1. The values of υ ψ and D ψ are obtained from the regression equations based on the velocity and retardation factor of experiment 1-4 itself. The predicted and the measured BTCs for experiment 1-4 are presented in Fig. 10c. These values of β, t1, and t2 for the experiment 4-1 are also used to predict the BTC of experiment 2-4 (υ = 0.399 cm h− 1, aggregate diameter ≤ 2 mm) (Fig. 10d). Values for υψ and Dψ are determined in the same manner of independent prediction for experiment 1-4. That the same values for β, t1, and t2 can be used to predict experiments with different velocities and aggregate sizes suggests that the transport process are relatively insensitive to changes in velocity and aggregate size, at least over the range of values associated with these four experiments, as reported by Brusseau et al. (1989). The third group of data sets is from Gaber et al. (1995) discussed in section 3.3. The first data set selected for use is from the highest pore water velocity (υ = 2.16 cm h− 1) treatment. As can be seen in Fig. 11a, the use of parameter set

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

Fig. 8. Representative breakthrough curves generated by HYDRUS-1D. R = 1.3023 (Kd = 0.1 mg L− 1), υ varies from 0.25 to 2.5 cm h− 1 (υ = 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2, 2.25, 2.5 cm h− 1).

calculated from the regression equations listed in Table 6 results in a BTC shift to the right and higher peak concentrations. This may be due to the failure to use the empirical correlation equation reported in Eq. (14) that relates the sorption rate coefficient to the equilibrium sorption coefficient. Such failure to predict the atrazine BTCs of Gaber et al. (1995) reflects the shortcoming of our proposed method. Despite that, our predictive model can provide better predictions for medium pore water velocity treatments, as can be seen in Fig. 11b, c. The regression relationships are further tested for its capability to predict transport of pesticides with high input concentrations. Rao and Davidson (1979) measured BTCs from miscible displacement of two pesticides through three soils for two input concentrations of each pesticide. The two pesticides used were: 2,4-D (2,4-dichlorophenoxyacetic acid) and atrazine. BTCs were measured for 2,4-D amine displacement at two input concentrations (C0 = 50 and 5000 μg mL− 1) using Webster (column length L =15.90 cm) and Cecil (L = 15.79 cm) soils, and two concentrations (C0 =5 and 50 μg mL− 1) of atrazine using Eustis soil (L = 15.79 and 14.90 cm). These BTC data were later used to evaluate two conceptual models for describing the nonequilibrium adsorption–desorption of pesticides in soils under steady-state water flow conditions (Rao et al., 1979). They observed that pesticide mobility was greater for the high input concentration as indicated by the left-hand shifts of the BTC. As can be seen in Fig. 12, the TPL model provides reasonable predictions of these digitized data sets although it overpredicts the tailing of the BTCs. Additionally, in contrast to the high concentration treatment (Fig. 12a), the proposed approach provides better predictions for the two lower concentration treatments (Fig. 12b, c). From the analyses of these independent experimental data sets, two points can be given. First, the predictive method presented provides an alternative approach to predict a sorbing solute transport when its nonequilibrium conditions are similar to that used in the numerical experiments to obtain the regression relationships, i.e., the α–Kd empirical relationship (Eq. (14)) and the value for F are relatively valid for the transport processes to be described.

147

Next, although each published data set has different column lengths that are different from the condition in the numerical experiments (L = 20 cm), the reasonable prediction can be obtained in most cases. This suggests that, within some range of length scale (in our study, the column length varies in the range of 10–30 cm), the TPL model may be easily used to predict breakthrough curves for different values of L, once the relevant parameters have been obtained for one breakthrough curve. In other words, dimensional parameters of the CTRW model determined through optimization are reasonably independent of column length, at least within the range employed above (Berkowitz et al., 2001). With regard to the numerical experiments, our ultimate goal is to estimate CTRW parameters for a sorbing solute independently of the BTC data to be simulated. For this purpose, we relate CTRW parameters for a sorbing solute to some measurable quantities (viz. υ and R) using artificially generated BTCs. The HYDRUS-1D package is selected to generate the required BTC “data” since it has been widely used and proven to be a reasonable tool for simulation of transport. We assume that the generated BTCs can approximate real experimental results. With respect to generation of BTCs data using HYDRUS-1D model, we note two points. First, the artificial BTCs are based on a two-site model and their decay need not display power law behavior. However, these generated BTCs are apparently asymmetric and long tailed, behind which the process can be captured by the CTRW model. In addition, we fit the generated BTCs with the CTRW model using a full spectrum of ψ(t), which is a compact form expressing the entire shape of a migrating solute plume: the power law region and the “cutoff” of this region allowing a transition from non-Fickian behavior to Fickian behavior at longer times. Therefore, it is expected that the CTRW model will produce good fits to these numerical BTCs. Simulation results in section 4.2 above are consistent with this expectation. Second, it is interesting and necessary to investigate the relationship between CTRW parameters and some measurable parameters from sufficient BTCs which are subject to

Fig. 9. Measured (after Mao and Ren (2004)) and predicted breakthrough curves in the second column with L = 20 cm, υ = 1.9642 cm h− 1, and R = 1.8253. TPL parameter values calculated from the regression equations are υψ = 1.2038 cm h− 1, Dψ = 2.0895 cm2 h− 1, β = 0.2487, t1 = 3.05 h, and t2 = 4.17 h.

148

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

Fig. 10. Predicted breakthrough curves for 2,4,5-T data from van Genuchten et al. (1977). Column length L = 30 cm for all experiments. (a) υ = 1.535 cm h− 1, R = 2.223, υψ = 0.8014 cm h− 1, Dψ = 1.4965 cm2 h− 1; (b) υ = 1.491 cm h− 1, R = 2.139, υψ = 0.8014 cm h− 1, Dψ = 1.4965 cm2 h− 1; (c) υ = 0.4501 cm h− 1, R = 2.225, υψ = 0.2181 cm h− 1, Dψ = 0.4387 cm2 h− 1; (d) υ = 0.399 cm h− 1, R = 2.21, υψ = 0.1932 cm h− 1, Dψ = 0.3905 cm2 h− 1. In all four prediction curves, TPL parameter values are β = 0.03, t1 = 1.56 h, and t2 = 2.02 h, which are calculated based on the pore water velocity and retardation factor of experiment 4-1.

chemical or (and) physical nonequilibrium. In our work, the HYDRUS-1D package considering nonequilibrium sorption is selected to generate the required BTC “data”. The only nonideal behavior in these numerically generated BTCs is due to the nonequilibrium sorption. In fact, we should have generated BTC using a multi-process, nonequilibrium model which considers simultaneous physical and chemical nonequilibrium, i.e., adsorption sites are subdivided into equilibrium and nonequilibrium sites, and the water phase is subdivided into mobile and immobile regions. The reasons for this are as follows: (1) the two-region and two-site models are mathematically equivalent in nondimensional form for the linear sorption isotherm case (van Genuchten, 1981; van Genuchten and Wagenet, 1989). Hence, under the stated conditions in our numerical experiments, the models for both chemical and physical mechanisms are equivalent. This means that a single model (e.g., two-site model) may be used to represent any one of the two nonequilibrium mechanisms (see also Brusseau and Rao, 1989b); (2) all of the experimental BTCs used to validate our regression relationships are subject to chemical nonequilibrium, and furthermore, chemical nonequilibrium conditions for some of

these BTCs are key; and (3) last but not least, we generate BTCs exhibiting early breakthrough times with long time tailing, regardless of whether physical or chemical heterogeneity is responsible for such non-ideal behaviors. For the three reasons above, we build the regression equations using artificial BTCs undergoing chemical nonequilibrium conditions, and then apply these equations to several groups of BTCs which may exhibit both chemical or/and physical nonequilibrium conditions. 5. Conclusion In this study, firstly, we evaluate the continuous time random walk (CTRW) theory for its capability to characterize the breakthrough behavior of sorbing chemicals undergoing chemical (or) and physical nonequilibrium. Data-fitting results show that the CTRW approach convincingly captures the full evolution of pesticide BTCs in saturated or unsaturated columns. Then, we investigate the CTRW parameter space for sorbing solute transport behaviors under a wide range of pore water velocity (υ) and retardation factor (R) conditions and relate the parameters of CTRW model to these

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

149

to assess the transport of sorbing solute through soil. In particular, numerical experimental studies on the effects of pore water velocity and retardation factor on BTCs help to

Fig. 11. Measured (after Gaber et al. (1995)) and predicted breakthrough curves for soil columns with L=30 cm. (a) υ=2.16 cm h−1, R=3.74, υψ =0.7452 cm h−1, Dψ =1.4779 cm2 h−1, β=0.2541, t1 =2.03 h, t2 =2.77 h; (b) υ=0.74 cm h− 1 (θ=0.44), R=3.49, υψ =0.2442 cm h−1, Dψ =0.5323 cm2 h−1, β=0.9534, t1 =86.30 h, t2 =208.55 h; (c) υ=0.69 cm h− 1, R=3.82, υψ =0.208 cm h−1, Dψ =0.4606 cm2 h− 1, β=1.0793, t1 =79.36 h, t2 =241.10 h.

two measurable parameters by regression analysis. The obtained regression relationships can be easily applied for estimating the CTRW parameters for sorbing solutes and thus for predicting their transport behaviors in soils. Predictions of several experimental data sets show reasonable results, thus indicating that our method could be an alternative approach

Fig. 12. Measured (after Rao and Davidson (1979) and Rao et al. (1979)) and predicted breakthrough curves. (a) L = 15.90 cm, υ = 0.544 cm h− 1, R = 1.5559, υ ψ = 0.3738 cm h − 1, D ψ = 0.6415 cm 2 h − 1, β = 0.4075, t1 = 59.72 h, t2 = 305.84 h; (b) L = 15.79 cm, υ = 0.666 cm h− 1, R = 3.0551, υ ψ = 0.2462 cm h− 1, D ψ = 0.529 cm2 h− 1, β = 0.8967, t 1 = 76.07 h, t2 = 256.68 h; (c) L = 14.90 cm, υ = 0.745 cm h− 1, R = 2.3212, υψ = 0.3579 cm h− 1, Dψ = 0.7202 cm2 h− 1, β = 0.5808, t1 = 87.00 h, t2 = 205.45 h.

150

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151

shed additional light on the relationship between the soil properties and specific parameters of the CTRW theory. A total of 200 artificial BTCs and 21 experimental BTCs have been analyzed in this paper. The analyses of actual measurements and artificial data suggest that both the flow and the sorption conditions may largely influence the key parameter β when modeling transport of sorbing solutes in soil. Higher velocity which appears to indicate higher degree of nonequilibrium may result in lower β value. Thus we expect that β would act as a function of the degree of chemical or physical nonequilibrium. A small value for β may indicate a relatively higher degree of nonequilibrium, and vice verse. More comprehensive interpretation of the physical meaning of the CTRW parameters and their relationship with the nonequilibrium parameters are needed in the future. Acknowledgments This study was partly supported by the National Key Basic Research Special Funds (NKBRSF, 2009CB118607) and by the National Natural Science Foundation of China (grant no. 50779064). We acknowledge the Agricultural Research Foundation of China–Israel (no. SIARF2001-05). We are grateful to Dr. Brian Berkowitz for many useful discussions and his generous help in answering many questions about the CTRW theory. We also thank the anonymous reviewers for their valuable comments and helpful suggestions, which led to an improvement of this paper.

References Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York. Benson, D.A., Wheatcraft, S.W., Meerschaert, M.M., 2000. The fractional-order governing equation of Lévy motion. Water Resour. Res. 36 (6), 1413–1423. Berkowitz, B., Scher, H., 1995. On characterization of anomalous dispersion in porous and fractured media. Water Resour. Res. 31 (6), 1461–1466. Berkowitz, B., Scher, H., 1997. Anomalous transport in random fracture networks. Phys. Rev. Lett. 79 (20), 4038–4041. Berkowitz, B., Scher, H., 1998. Theory of anomalous chemical transport in random fracture networks. Phys. Rev. E 57 (5), 5858–5869. Berkowitz, B., Scher, H., 2001. The role of probabilistic approaches to transport theory in heterogeneous porous media. Transp. Porous Media 42, 241–263. Berkowitz, B., Scher, H., 2009. Exploring the nature of non-Fickian transport in laboratory experiments. Adv. Water Resour. 32, 750–755. Berkowitz, B., Scher, H., Silliman, S.E., 2000. Anomalous transport in laboratoryscale, heterogeneous porous media. Water Resour. Res. 36 (1), 149–158. Berkowitz, B., Kosakowski, G., Margolin, G., Scher, H., 2001. Application of continuous time random walk theory to tracer test measurements in fractured and heterogeneous porous media. Ground Water 39 (4), 593–604. Berkowitz, B., Klafter, J., Metzler, R., Scher, H., 2002. Physical pictures of transport in heterogeneous media: advection–dispersion, random-walk, and fractional derivative formulations. Water Resour. Res. 38 (10), 1191. doi:10.1029/2001WR001030. Berkowitz, B., Cortis, A., Dentz, M., Scher, H., 2006. Modeling non-Fickian transport in geological formations as a continuous time random walk. Rev. Geophys. 44, RG2003. doi:10.1029/2005RG000178. Berkowitz, B., Emmanuel, S., Scher, H., 2008. Non-Fickian transport and multiplerate mass transfer in porous media. Water Resour. Res. 44, W03402. Brusseau, M.L., 1992. Transport of rate-limited sorbing solutes in heterogeneous porous media: application of a one-dimensional multi-factor nonideality model to field data. Water Resour. Res. 28 (9), 2485–2497. Brusseau, M.L., Rao, P.S.C., 1989a. Sorption nonideality during organic contaminant transport in porous media. CRC Crit. Rev. Environ. Control 19 (1), 33–99. Brusseau, M.L., Rao, P.S.C., 1989b. The influence of sorbate–organic matter interactions on sorption nonequilibrium. Chemosphere 18 (9/10), 1691–1706.

Brusseau, M.L., Rao, P.S.C., 1990. Modeling solute transport in structured soils: a review. Geoderma 46, 169–192. Brusseau, M.L., Srivastava, R., 1997. Non-ideal transport of reactive solutes in heterogeneous porous media: 2. Quantitative analysis of the Borden natural-gradient field experiment. J. Contam. Hydrol. 28, 115–155. Brusseau, M.L., Jessup, R.E., Rao, P.S.C., 1989. Modeling the transport of solutes influenced by multiprocess nonequilibrium. Water Resour. Res. 25 (9), 1971–1988. Brusseau, M.L., Jessup, R.E., Rao, P.S.C., 1991a. Nonequilibrium sorption of organic chemicals: elucidation of rate-limiting processes. Environ. Sci. Technol. 25 (1), 134–142. Brusseau, M.L., Larsen, T., Christensen, T.H., 1991b. Rate-limited sorption and nonequilibrium transport of organic chemicals in low organic carbon aquifer materials. Water Resour. Res. 27 (6), 1137–1145. Cameron, D.R., Klute, A., 1977. Convective–dispersive solute transport with a combined equilibrium and kinetic adsorption model. Water Resour. Res. 13 (1), 183–188. Carrera, J., Sánchez-Vila, X., Benet, I., Medina, A., Galarza, G., Guimerà, J., 1998. On matrix diffusion: formulations, solution methods, and qualitative effects. Hydrogeol. J. 6, 178–190. Cortis, A., Berkowitz, B., 2004. Anomalous transport in “classical” soil and sand columns. Soil Sci. Soc. Am. J. 68, 1539–1548. Cortis, A., Berkowitz, B., 2005. Computing “anomalous” contaminant transport in porous media: the CTRW MATLAB toolbox. Ground Water 43 (6), 947–950. Cortis, A., Gallo, C., Scher, H., Berkowitz, B., 2004. Numerical simulation of nonFickian transport in geological formations with multiple-scale heterogeneities. Water Resour. Res. 40, W04209. doi:10.1029/2003WR002750. de Hoog, F.R., Knight, J.H., Stokes, A.N.,1982. An improved method for numerical inversion of Laplace transforms. S.I.A.M., J. Sci. Stat. Comput. 3, 357–366. Deng, J.C., Jiang, X., Zhang, X.X., Hu, W.P., Crawford, J.W., 2008. Continuous time random walk model better describes the tailing of atrazine transport in soil. Chemosphere 71, 2150–2157. Dentz, M., Berkowitz, B., 2003. Transport behavior of a passive solute in continuous time random walks and multirate mass transfer. Water Resour. Res. 39 (5), 1111. doi:10.1029/2001WR001163. Dentz, M., Cortis, A., Scher, H., Berkowitz, B., 2004. Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport. Adv. Water Resour. 27 (2), 155–173. Fetter, C.W., 1993. Contaminant Hydrogeology. Prentice-Hall, Inc., A Simon & Schuster Company, Upper Saddle River, New Jersey. 71. Gaber, H.M., Inskeep, W.P., Comfort, S.D., Wraith, J.M., 1995. Nonequilibrium transport of atrazine through large intact soil cores. Soil Sci. Soc. Am. J. 59, 60–67. Gamerdinger, A.P., Wagenet, R.J., van Genuchten, M.Th., 1990. Application of two-site/two-region models for studying simultaneous nonequilibrium transport and degradation of pesticides. Soil Sci. Soc. Am. J. 54, 957–963. Gamerdinger, A.P., Lemley, A.T., Wagenet, R.J., 1991. Nonequilibrium sorption and degradation of three 2-Chloro-s-Triazine herbicides in soil–water systems. J. Environ. Qual. 20, 815–822. Gheorghiu, S., Coppens, M.-O., 2004. Heterogeneity explains features of “anomalous” thermodynamics and statistics. PNAS 101 (45), 15852–15856. Haggerty, R., Gorelick, S.M., 1995. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Resour. Res. 31 (10), 2383–2400. Haggerty, R., McKenna, S.A., Meigs, L.C., 2000. On the late time behavior of tracer test breakthrough curves. Water Resour. Res. 36 (12), 3467–3479. Harvey, C.F., Gorelick, S.M., 1995. Temporal moment-generating equations: modeling transport and mass transfer in heterogeneous aquifers. Water Resour. Res. 31 (8), 1895–1911. Hatano, Y., Hatano, N., 1998. Dispersive transport of ions in column experiments: an explanation of long-tailed profiles. Water Resour. Res. 34 (5), 1027–1033. Jiménez-Hornero, F.J., Giráldez, J.V., Laguna, A., Pachepsky, Y., 2005. Continuous time random walks for analyzing the transport of a passive tracer in a single fissure. Water Resour. Res. 41, W04009. doi:10.1029/ 2004WR003852. Kosakowski, G., Berkowitz, B., Scher, H., 2001. Analysis of field observations of tracer transport in a fractured till. J. Contam. Hydrol. 47, 29–51. Lallemand-Barres, P., Peaudecerf, P.,1978. Recherche des relations entre la valeur de la dispersivite macroscopique d'un milieu aquifere, ses autres caracteristiques et les conditions de mesure, etude bibliographique. Bulletin, Bureau de Recherches Geologiques et Minieres Sec. 3/4, 277–287. Levy, M., Berkowitz, B., 2003. Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. J. Contam. Hydrol. 64, 203–226. Ma, L.W., Selim, H.M., 1994a. Predicting atrazine adsorption–desorption in soils: a modified second-order kinetic model. Water Resour. Res. 30 (2), 447–456. Ma, L.W., Selim, H.M., 1994b. Predicting the transport of atrazine in soils: second-order and multireaction approaches. Water Resour. Res. 30 (12), 3489–3498.

N. Li, L. Ren / Journal of Contaminant Hydrology 108 (2009) 134–151 Ma, L.W., Selim, H.M., 2005. Predicting pesticide transport in Mulch-Amended soils: a two-compartment model. Soil Sci. Soc. Am. J. 69, 318–327. Mao, M., Ren, L., 2004. Simulating nonequilibrium transport of atrazine through saturated soil. Ground Water 42 (4), 500–508. Margolin, G., Berkowitz, B., 2000. Application of continuous time random walks to transport in porous media. J. Phys. Chem., B 104 (16), 3942–3947. Margolin, G., Dentz, M., Berkowitz, B., 2003. Continuous time random walk and multirate mass transfer modeling of sorption. Chem. Phys. 295, 71–80. Metzler, R., Klafter, J., 2000. The random walk's guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 339 (1), 1–77. Mojid, M.A., Vereecken, H., 2005. On the physical meaning of retardation factor and velocity of a nonlinearly sorbing solute. J. Hydrology 302, 127–136. Montroll, E.W., Scher, H., 1973. Random walks on lattices, IV, Continuous-time walks and influence of adsorbing boundaries. J. Stat. Phys. 9, 101–135. Nakayama, S., Vandergraaf, T.T., Kumata, M., 1994. Experimental study on nuclides migration under deep geological condition (in Japanese). Radioact. Waste Res. 1, 67–76. Nkedi-Kizza, P., Biggar, J.W., Selim, H.M., van Genuchten, M.Th., Wierenga, P.J., Davidson, J.M., Nielson, D.R., 1984. On the equivalence of two conceptual models for describing ion exchange during transport through an aggregated oxisol. Water Resour. Res. 20 (8), 1123–1130. Parker, J.C., van Genuchten, M.Th, 1984. Determining Transport Parameters from Laboratory and Field Tracer Experiments. InVirginia Agricultural Experiment Station, Blacksburg, VA. Bulletin 84-3. Pfister, G., Scher, H., 1978. Non-Gaussian transient transport in disordered solids. Adv. Phys. 27, 747–798. Rao, P.S.C., Davidson, J.M., 1979. Adsorption and movement of selected pesticides at high concentrations in soils. Water Res. 13, 375–380. Rao, P.S.C., Davidson, J.M., Jessup, R.E., Selim, H.M., 1979. Evaluation of conceptual models for describing nonequilibrium adsorption–desorption of pesticides during steady-flow in soils. Soil Sci. Soc. Am. J. 43, 22–28. Roth, K., Jury, W.A., 1993. Linear transport models for adsorbing solutes. Water Resour. Res. 29 (4), 1195–1203. Scher, H., Montroll, E.W., 1975. Anomalous transit-time dispersion in amorphous solids. Phys. Rev. B 12 (6), 2455–2477.

151

Scher, H., Margolin, G., Berkowitz, B., 2002. Towards a unified framework for anomalous transport in heterogeneous media. Chem. Phys. 284, 349–359. Schumer, R., 2002. Fractional Derivatives, Continuous Time Random Walks, and Anomalous Solute Transport. Ph.D. Thesis, Univ. of Nev., Reno. Schumer, R., Benson, D.A., Meerschaert, M.M., 2003. Fractal mobile/immobile solute transport. Water Resour. Res. 39 (10), 1296. doi:10.1029/ 2003WR002141. Selim, H.M., Ma, L.W., 1995. Transport of reactive solutes in soils: a modified two-region approach. Soil Sci. Soc. Am. J. 59, 75–82. Selim, H.M., Amacher, M.C., 1988. A second-order kinetic approach for modeling solute retention and transport in soils. Water Resour. Res. 24 (12), 2061–2075. Selim, H.M., Davidson, J.M., Mansell, R.S., 1976. Evaluation of a two-site adsorption–desorption model for describing solute transport in soils. Proceedings of the Summer Computer Simulation Conference, Washington, DC, pp. 444–448. Selim, H.M., Ma, L.W., Zhu, H.X., 1999. Predicting solute transport in soils: second-order two-site models. Soil Sci. Soc. Am. J. 63, 768–777. Seuntjens, P., Mallants, D., Šimunek, J., Patyn, J., Jacques, D., 2002. Sensitivity analysis of physical and chemical properties affecting field-scale cadmium transport in a heterogeneous soil profile. J. Hydrol. 264, 185–200. Šimunek, J., Šejna, M., van Genuchten, M.T.h., 1998. The HYDRUS-1D software package for simulating water flow, heat, and solute transport in onedimensional variably saturated media. Version 2.0. IGWMC-TPS-70, International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado. van Genuchten, M.T.h., 1981. Non-equilibrium transport parameters from miscible displacement experiments. Res. Rep. 119, U.S. Salinity Lab., USDA-ARS, Riverside, CA. van Genuchten, M.T.h., Wierenga, P.J., 1976. Mass transfer studies in sorbing porous media: I. Analytic solutions. Soil Sci. Soc. Am. J. 40, 473–480. van Genuchten, M.T.h., Wagenet, R.J., 1989. Two-site/two-region models for pesticide transport and degradation: theoretical development and analytical solutions. Soil Sci. Soc. Am. J. 53, 1303–1310. van Genuchten, M.Th., Wierenga, P.J., O'Connor, G.A., 1977. Mass transfer studies in sorbing porous media: III. Experimental evaluation with 2,4,5-T. Soil Sci. Soc. Am. J. 41, 278–285.