A finite element solution for the fractional advection–dispersion equation

Advances in Water Resources 31 (2008) 1578–1589

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

A finite element solution for the fractional advection–dispersion equation Quanzhong Huang a,b, Guanhua Huang a,b,*, Hongbin Zhan c a b c

Center for Agricultural Water Research in China, China Agricultural University, Beijing 100083, PR China Chinese-Israeli International Center for Research and Training in Agriculture, China Agricultural University, Beijing 100083, PR China Department of Geology and Geophysics, Texas A&M University, College Station, Texas 77845-3115, USA

a r t i c l e

i n f o

Article history: Received 11 September 2007 Received in revised form 16 June 2008 Accepted 1 July 2008 Available online 5 July 2008 Keywords: Finite element method Fractional advection–dispersion equation Caputo derivative Non-Fickian dispersion

a b s t r a c t The fractional advection–dispersion equation (FADE) known as its non-local dispersion, has been proven to be a promising tool to simulate anomalous solute transport in groundwater. We present an unconditionally stable finite element (FEM) approach to solve the one-dimensional FADE based on the Caputo definition of the fractional derivative with considering its singularity at the boundaries. The stability and accuracy of the FEM solution is verified against the analytical solution, and the sensitivity of the FEM solution to the fractional order a and the skewness parameter b is analyzed. We find that the proposed numerical approach converge to the numerical solution of the advection–dispersion equation (ADE) as the fractional order a equals 2. The problem caused by using the first- or third-kind boundary with an integral-order derivative at the inlet is remedied by using the third-kind boundary with a fractional-order derivative there. The problems for concentration estimation at boundaries caused by the singularity of the fractional derivative can be solved by using the concept of transition probability conservation. The FEM solution of this study has smaller numerical dispersion than that of the FD solution by Meerschaert and Tadjeran (J Comput Appl Math 2004). For a given a, the spatial distribution of concentration exhibits a symmetric non-Fickian behavior when b = 0. The spatial distribution of concentration shows a Fickian behavior on the left-hand side of the spatial domain and a notable non-Fickian behavior on the right-hand side of the spatial domain when b = 1, whereas when b = 1 the spatial distribution of concentration is the opposite of that of b = 1. Finally, the numerical approach is applied to simulate the atrazine transport in a saturated soil column and the results indicat that the FEM solution of the FADE could better simulate the atrazine transport process than that of the ADE, especially at the tail of the breakthrough curves. Published by Elsevier Ltd.

1. Introduction The advection–dispersion equation (ADE), which is based on the Fick’s law, is commonly used to simulate contaminant transport in porous media [1]. But increasing evidence over the past two decades shows that contaminants moving through porous media usually exhibit anomalous or non-Fickian processes, including heavy leading and/or trailing edges and scale-dependent dispersion coefficients [48,13,34,35]. Heterogeneities are regarded as the primary causes of non-Fickian migration of solute in porous media [5]. Migration of solute is sensitive to heterogeneities at all scales, and the small-scale heterogeneities can significantly affect largescale behavior. It has been proven that even in carefully packed laboratory-scale flow cells and columns, pore-scale heterogeneity exists, giving rise to variation on water velocity field. Transport

* Corresponding author. Address: Chinese-Israeli International Center for Research and Training in Agriculture, China Agricultural University, Beijing 100083, PR China. Tel.: +86 10 62737144; fax: +86 10 62737138. E-mail address: [email protected] (G. Huang). 0309-1708/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.advwatres.2008.07.002

process in such heterogeneous geological media often cannot be well described by the classical ADE [5]. Several attempts have been made to model the anomalous transport, including the mobile–immobile (MIM) model [9,44], the multiple-rate and fractal MIM [39,16], the continuous time random walk (CTRW) [5,4], the non-local stochastic transport theories [30,31], and the fractional advection–dispersion equation (FADE) [2,3]. The MIM based approach seems make a big improvement over ADE, but the disadvantage of this approach is also apparent as reviewed by many researchers [34,35,52,53,11]. The CTRW has shown significant improvement over ADE [5,4], but it still needs further development. For instance, there are still no studies using CTRW to solve the reactive transport in groundwater. The non-local stochastic transport theories are capable of modeling non-Fickian transport by introducing spatial/temporal non-local parameters. Zhang et al. [51] has shown mathematically that all variants of FADE with variable coefficients, proposed to date, are special cases of a nonlocal stochastic transport theory originally proposed by Neuman [28], obtained upon replacing Neuman’s kernels with time-localized forms having power-law spatial memory.

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589

By the same token, FADE with constant coefficients may be obtained from Neuman’s theory upon taking the underlying random velocity field to be stationary, in which case Neuman’s theory reduces to that of Cushman and Ginn [10]. But the needed statistical functions for the nonlocal stochastic theory may not always be easy to obtain [5,4]. The FADE has shown some advantages to simulate non-Fickian transport process of solute via fractional derivative in space and will be the focus of this study [34,10]. Cushman and Ginn [10] claimed that the FADE could be seen as a special case of the classical mass balance equation with a convolution-Fickian flux to model non-local dispersion. Schumer et al. [40] used a generalized central limit theorem which summed the length of particle jumps during their random walk through a heterogeneous porous medium to demonstrate the behavior of the FADE. Similar to the ADE, the FADE can be solved with either analytical or numerical methods. Analytical solutions of the FADE can only be obtained at very simply initial and boundary conditions [2,18]. Numerical solutions are especially important since the boundary and initial conditions are always complicated in practical applications. Efforts have been made for developing numerical solutions of the FADE with the finite difference (FD) method [11,22,23,25,26], the finite volume (FV) method [52,53], the random walk (RW) method [14,15,47,27], the random walk particle tracking (RWPT) method [49,50], as well as the finite element (FEM) method [38]. Existing numerical methods of the FADE still have a number of problems unsolved. For instance, there is still no study about the numerical dispersion and artificial oscillation of the FD/FV solution for FADE. The RW and/or RWPT solution to ADE suffers the lack of a fixed grid or fixed coordinates, which may lead to numerical instability and computational difficulties, particularly in non-uniform media with multiple sinks or sources and complex boundary conditions [46]. The similar problem is likely to occur in the RW and/or RWPT solution for FADE. And most numerical solutions of the FADE were derived from a fractional derivative based on the Riemann–Liouville or the Grünwald–Letnikov definitions [22,38,11,23,25,26,14,15], which may cause mass balance error as proven by Zhang et al. [53]. A straightforward way to keep mass balance is to redefine the fractional dispersive flux by the Caputo definition [53]. Other shortcomings of the Riemann–Liouville and/or the Grünwald–Letnikov formula are discussed in Section 2 of this paper. Besides, some numerical solutions of the FADE only take the left-handed fractional derivative into account [11,23,25,49,50], implying that the backward transition of solute particles that is slower than the average velocity is eliminated. Unlike the local behavior of an integer-order derivative, the fractional derivative is well known for its non-local property. The fractional derivative of a certain function at a given point not only depends on the function values at nearby points, but also depends on the characteristics of the function across the entire spatial domain [40]. Therefore, only using the left-handed fractional derivative in the FADE will loss the information from the right-hand side of the function. In this study, we present an unconditionally stable numerical scheme for the FADE based on a FEM using the Caputo definition of the fractional derivative. Both the left- and right-handed fractional derivatives, which represent the forward and backward transition of solute particles, respectively, are taken into account. The convergence of the FEM solution is analyzed by comparing the results with the numerical solution of the ADE as the fractional order a = 2. The accuracy of the FEM solution is also examined by comparing with the analytical solution of the FADE. We have paid great attention to the influence of the inlet and outlet boundaries on the simulation results, and the sensitivity of the numerical solution to the fractional order and skewness. Finally, the FEM solution is applied to simulate a laboratory experiment of atrazine movement in a saturated soil column.

1579

2. Theory 2.1. Basic equation The FADE is characterized by its non-local transport behavior over large distances via a convolutional fractional-derivative. Zhang et al. [49,50] proposed three forms of the FADE, which are the ADE with a fractional flux (FF–ADE) and the ADE with a fractional divergence (FD–ADE) and the ADE with a fully fractional divergence (FFD–ADE). They found that the FF–ADE and the FD– ADE result in very similar plume behavior if the diffusion coefficient (D) varied linearly with distance, providing that its gradient was relatively small. The FF–ADE arises when the non-locality occurs only in the form of concentration gradients. This form is the same as the fractional nonlinear Fokker–Planck equation proposed by Tsallis and Lenzi [43]. When the mean groundwater velocity (v) and the dispersion coefficient (D) are constant, the FF-ADE reduces to the one used by Benson et al. [2]. One-dimensional fractional flux advection–dispersion equation with a constant D and v can be written as

oC oC o ¼ v þ D ot ox ox

    a1 ! 1 þ b oa1 C 1b o C ; þ 2 oxa1 2 oðxÞa1

ð1Þ

where C is the concentration of contaminant, t is the duration of the contaminant evolving in space, x is the spatial coordinate, a is the order of the fractional derivative, b (1 6 b 6 1) is the skewness of the transport process, which controls the bias of the dispersion. The parameter b reflects the relative weight of forward versus backward transition probability. If b is smaller than zero, the dispersion is skewed backward describing a slowly evolving contaminant plume followed by a heavy tail. If b is greater than zero, the dispersion is skewed forward describing a fast evolving contaminant plume followed by a light tail. If the solute forward transition has the same probability as the backward transition, b equals to zero. This is the symmetric case for the solute transport. Three formulas are commonly used to describe the fractional derivative. These formulas are the Riemann–Liouville formula, the Grünwald–Letnikov formula, and the Caputo formula [29,37,6,21]. The Riemann–Liouville derivative within the interval [0, L] has the form as

8 a 1 > < o oxf ðxÞ a ¼ CðnaÞ a > 1 : o f ðxÞa ¼ Cðn aÞ oðxÞ

i  fÞn1a f ðfÞdf ; 0 6 f < x; hR i L on ðf  xÞn1a f ðfÞdf ; x < f 6 L; x oxn

on oxn

hR

x ðx 0

ð2Þ

where n is the smallest integer larger than the fractional order a , a oa f ðxÞ i.e.,n  1 < a < n. o oxf ðxÞ a and oðxÞa are the left and right-handed fractional derivatives of order a, respectively. C(z) is the gamma function. L is the size of the domain. If one assumes that the function f(x) is smooth enough, then the Riemann–Liouville derivatives of Eq. (2) are equivalent to the Grünwald–Letnikov derivatives, which are

8 n1 Rx P ðx0Þka f ðkÞ ð0Þ > n1a ðnÞ oa f ðxÞ 1 > f ðfÞdf; > < oxa ¼ Cðkþ1aÞ þ CðnaÞ 0 ðx  fÞ

0 6 f < x;

k¼0

n1 > RL P > k ðLxÞka f ðkÞ ðLÞ n1a ðnÞ oa f ðxÞ ð1Þn > f ðfÞdf; x < f 6 L: : oðxÞa ¼ ð1Þ Cðkþ1aÞ þ CðnaÞ x ðf  xÞ k¼0

ð3Þ The equivalency between Eqs. (2) and (3) can be verified by the method of integration by parts [37,21]. Several numerical schemes for solving the FADE have been developed on the basis of the definitions of the Riemann–Liouville and/or the Grünwald–Letnikov [22,23,25]. However, Caputo [7], Caputo and Mainardi [8], Kilbas et al. [20], Heymans and Podlubny [17], Li and Deng [21], and Zhang et al. [53] have pointed out that

1580

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589

the definition of the fractional derivative as shown in Eq. (2) has several shortcomings. First, it may cause hyper-singular improper integral. Second, the derivative of a constant is not zero. Third, the fractional derivative involved in the initial condition is often ill-defined. Furthermore, Zhang et al. [53] have proven that the Riemann–Liouville definition may cause a mass balance problem, and suggested to use the Caputo definition to approximate the fractional derivative. Podlubny [37], Gorenflo et al. [15], and Butzer and Westphal [6] have pointed out that the Caputo fractional derivative represents a sort of regularization in the time origin for the Riemann–Liouville fractional derivative and satisfies the requirement of being zero when applied to a constant. Besides, the Caputo method does not use the fractional-order derivative in the initial condition, thus is convenient in physical and engineering applications where the initial conditions are usually expressed in terms of the integer-order derivatives. Therefore, in the following analysis, we use the Caputo formula to describe the fractional derivatives, which are oa f ðxÞ oxa oa f ðxÞ oðxÞa

1 ¼ Cðn aÞ 1 ¼ Cðn aÞ

Rx 0

n

ðx  fÞna1 o oxf ðfÞ n df;

RL x

ðf 

n xÞna1 o oxf ðfÞ n df;

0 6 f < x; x < f 6 L:

ð4Þ

1 ¼ Cð2 aÞ

oa1 f ðxÞ oðxÞa1

1 ¼ Cð2 aÞ

Rx 0

ðx  fÞð1aÞ ofoxðfÞ df;

RL x

ðf  xÞð1aÞ ofoxðfÞ df:

Applying the method of integration by parts, one can easily obtain the following result

!  Z xjþ1 Z x Dð1 þ bÞ o ðx  fÞð1aÞ C 0 ðfÞdf /j dx 2Cð2  aÞ ox xj1 0  Z xjþ1 Z x Dð1 þ bÞ b 0 ðfÞdf /0 ðxÞdx: ðx  fÞ1a C  j 2Cð2  aÞ xj1 0

F3 ¼ 

Eq. (11) contains integrations at the interval [0, x], which can be written as [0, xj1] [ [xj1, x] for x 2 [xj1, xj] or [0, xj] [ [xj, x] for Rx b 0 ðfÞdf can be x 2 [xj, xj+1]. Thus the integration term 0 ðx  fÞ1a C R Pj1 R xi x b 0 ðfÞdf þ b 0 ðfÞdf ðxj  fÞ1a C approximated by i¼1 xi1 ðxi  fÞ1a C xj1 Rx Pj R xi 1a b 0 1a b 0 for x 2 [xj1, xj] or C ðfÞdf þ xj ðxj  fÞ C ðfÞdf i¼1 xi1 ðxi  fÞ b 0 ðfÞ in each sub-domain for x 2 [xj, xj+1]. The differential term C 0 b ½xi1 ; xi  can be evaluated by C ðfÞ ¼ C i1 u0i1 þ C i u0i . The detailed derivation of the third term on the left-hand side of Eq. (7) is given in Appendix. The final form is

 C i ðj  i  2Þ3a þ 4ðj  i  1Þ3a  6ðj  iÞ3a i¼1    þ4ðj  i þ 1Þ3a  ðj  i þ 2Þ3a þ C j1 6 þ 25a  33a   o þC j 4  23a  C jþ1 ð12Þ þ

ð5Þ

Substituting Eq. (5) into Eq. (1), one has

 Z x oC oC o 1þb ðx  fÞð1aÞ C 0 ðfÞdf þv D ot ox ox 2Cð2  aÞ 0  Z L 1b þ ðf  xÞð1aÞ C 0 ðfÞdf ¼ 0: 2Cð2  aÞ x

ð6Þ

L

ð7Þ

ð8Þ

Using a uniform spatial discretization with a grid space of D, one has:

F1 ¼

Z

xjþ1

xj1

F2 ¼

Z

xjþ1

xj1

b oC oC j D; / dx  ot j ot b oC C jþ1  C j1 v/ dx  v : ox j 2

  þ4ði  j þ 1Þ3a  ði  j þ 2Þ3a þ C N ðN  j þ 1Þ3a o 3ðN  jÞ3a þ 3ðN  j  1Þ3a  ðN  j  2Þ3a 

ð9Þ ð10Þ

ð13Þ

Dð1bÞ D1a . 2Cð4aÞ

where GM ¼ Above Galerkin method is used only for approximating the spatial derivatives. For the time derivatives, we use the backward Euler scheme. Setting the time step as s, and substituting Eqs. (9), (10), (12), (13) into Eq. (7), one can obtain

!

 sv s  þ  þ WA C kþ1 j1 2D D D i¼1     s sv s þ 1 þ WB C jkþ1 þ þ WE C kþ1 jþ1 2D D D ! N2 X s ¼ C kj j ¼ 1; 2; . . . ; N  1; þ GM C kþ1 W i þ W N C kþ1 i N D i¼jþ2

s

b ðx; tÞ ¼ PN C j ðtÞ/ ðxÞ  Cðx; tÞ, N is the number of sub-dowhere C j j¼0 main within [0, L], j is the node index and /j is the weighting function and the linear basis function associated with node j. The first two terms in Eq. (7) are the same as the FEM scheme of the ADE. The linear basis function is [19,42]

8 xx   < x xj1 x 2 xj1 ; xj ; j j1 /j ¼ xjþ1 x   : x 2 xj ; xjþ1 : xjþ1 xj

n   F 4 ¼ GM C j1 þ C j ð4  23a Þ þ C jþ1 6 þ 25a  33a N1  X C i ði  j  2Þ3a þ 4ði  j  1Þ3a  6ði  jÞ3a þ jþ2

The FEM has been widely used in the simulation of groundwater flow and solute transport [19,42]. Among them, a commonly used FEM is the Galerkin method, which assumes that the dependent variable, i.e. the concentration function C(x, t) can be approximated b ðx; tÞ, with which the basis function serves as the by a finite series C weighting function [19,42]. According to the Galerkin method, the FEM scheme for the FADE can be expressed as

" Z x  b b oC oC Dð1 þ bÞ o ðx  fÞð1aÞ C 0 ðfÞdf þv  ot ox 2Cð2  aÞ ox 0 0 Z L  Dð1  bÞ o  ðf  xÞð1aÞ C 0 ðfÞdf /j dx ¼ 0; 2Cð2  aÞ ox x

j2 X

1a in which GMþ ¼ 2Dð1þbÞ . As show in Appendix, the fourth term Cð4aÞ D on the left-hand side of Eq. (7) can be expressed as

2.2. Finite element solution of the FADE

Z

ð11Þ

n   F 3 ¼ GMþ þC 0 ðj  2Þ3a þ 3ðj  1Þ3a  3ðjÞ3a þ ðj þ 1Þ3a

For the case of 1 < a < 2, one has oa1 f ðxÞ oxa1

The third term on the left-hand side of Eq. (7) contains integraR  x ðx  fÞð1aÞ C 0 ðfÞdf /j dx in the sub-domain of [xj1, xj+1]. 0

tion of oxo

þ

GM

W þ0 C kþ1 0

þ

j2 X

C kþ1 W þi i

ð14Þ where,

W þi ¼ ðj  i  2Þ3a þ 4ðj  i  1Þ3a  6ðj  iÞ3a þ 4ðj  i þ 1Þ3a  ðj  i þ 2Þ3a ;  W i ¼ ði  j  2Þ3a þ 4ði  j  1Þ3a  6ði  jÞ3a þ 4ði  j þ 1Þ3a  ði  j þ 2Þ3a ;   WA ¼ GMþ 6 þ 4ð2Þ3a  ð3Þ3a  GM ;     WB ¼ GMþ 4  ð2Þ3a þ GM 4  ð2Þ3a ;   WE ¼ GMþ þ GM 6 þ 4ð2Þ3a  ð3Þ3a ; W þ0 ¼ ðj  2Þ3a þ 3ðj  1Þ3a  3ðjÞ3a þ ðj þ 1Þ3a for j P 2; W N ¼ ðN  j  2Þ3a þ 3ðN  j  1Þ3a  3ðN  jÞ3a þ ðN  j þ 1Þ3a :

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589

Eq. (14) is a linear system with N  1 equations and N + 1 variþ ables. W þ 0 and W i act as the weighting coefficients for i < j  1,  and W whereas W  N i act as the weighting coefficients for i > j + 1, respectively. The superscript k is the step number at time t. This linear system can be solved for given initial and boundary conditions. The three weighting coefficients of the series denoted by j  1, j, and j + 1 in Eq. (14) are dominating. Similar results were given by Deng et al. [11] and Zhang et al. [52] who called the dominating three terms of the series a Gaussian core. The Gaussian core exhibits a Gaussian behavior, while the remaining terms in the equations designate a non-Gaussian tail. It is easy to see from Fig. 1 that the weighting coefficients (Wi) decrease rapidly when the nodes departure away from the node j. The spatial domain was divided into 100 finite elements with 101 nodes in space, and a linear system of the FEM scheme has 101 variables. The 50th weighting coefficient W50 is very close to 3, while W49 and W51 are very close

3

1581

to 1. These three weighting coefficients are significantly greater than the others, with which their absolute values are much smaller than unity. Similar result was also obtained by Schumer et al. [40] and Diethelm et al. [12]. Fig. 2 shows the matrix elements (Ai,j) composing the linear system, where the variable j stands for the jth equation, and the variable i stands for the ith coefficient of the jth equation. The number of nodes is 100 in this figure. It can be easily found that the tri-diagonal elements of the matrix are the pivotal elements. This means that the linear system expressed by Eq. (14) is expected to be numerically stable and has a unique  þ  solution. It is very interesting to note that W þ i , W i , W 0 and W N in Eq. (14) converge to zero for a = 2, WA, WB and WE in Eq. (14) equal to  2DD, DD and  2DD, respectively, and the numerical scheme of the FADE then becomes identical to that of the ADE. Eq. (14) is a simpler case of FADE, a more complicated form may be obtained by using non-constant v and D in Eq. (1) as proposed by Huang et al. [18] and Zhang et al. [51]. It has been demonstrated by Zhang et al. [51] that a FADE with variable parameters could be more suitable to the measured data at MADE site when using the FD method. Similar results could be expected if we use variable parameters in the FEM scheme for FADE. 2.3. Boundary conditions

Wi

2

We can solve the linear algebraic equations given by Eq. (14) with appropriate inlet and outlet boundary conditions, which could be the prescribed concentration (the first-kind), the prescribed flux (the second-kind), or the concentration-dependent flux (the third-kind) conditions [52,42]. Among these conditions, the prescribed concentration is relative simple and widely used as the inlet boundary condition with both the ADE [42,45] and the FADE [2,22,18,23,25]. This condition can be written as

1

0

C in ¼ C 0 ; -1 0

20

40

60

80

100

i Fig. 1. Distribution of the weighting coefficients (Wi, i = 0 ,. . . , N) of the finite element scheme based on the Caputo definition of the fractional derivative with a = 1.6, N = 100.

ð15Þ

where C0 and Cin are the concentrations at the inlet and the first node from upper stream, respectively. On the other hand, contaminant transport at the outlet is considered as a process dominated by advection, thus the zero concentration gradient or the so-called Danckwerts’ condition is usually considered there

oC=oxjx¼L ¼ 0:

ð16Þ

Fig. 2. Visualization of the weighting matrix elements Ai,j (i, j = 0 ,. . . , N) discussed in Section 3.2 for the linear system Eq. (14) and boundary conditions Eqs. (18) and (19).

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589

The use of the Danckwerts’ condition at the outlet has been justified by many investigators including van Genuchten and Parker [45]. However, improper treatment of the inlet boundary may cause significant mass balance errors on the simulation results for both the ADE [45] and the FADE [53]. It has been evident that using the third-kind condition to treat the inlet boundary will satisfy the mass conservation and may be more close to the real situations [45]. Therefore, it is used in the following numerical analysis

vC in  D

oC ¼ vC 0 ; ox x¼0

ð17Þ

Due to the fractional flux used in the FADE, the non-local behavior of transport may occur at the inlet and outlet boundaries as well. Therefore, using a third-kind boundary condition with a fractional-order derivative at the inlet may be desirable and consistent with the FADE. Zhang et al. [52,53] used such a condition as the inlet boundary and a zero gradient concentration with a fractionalorder derivative as the outlet boundary. To compare the effect of different kinds of boundary conditions on the FEM solution of the FADE, we follow the lines of Zhang et al. [52,53] to deal with the boundary conditions. The third-kind inlet boundary with a fractional-order derivative can be expressed as

    a1 ! 1 þ b oa1 C 1b o C vC in  D D 2 oxa1 2 oðxÞa1

¼ vC 0 ;

ð18Þ

1.2

1

0.8

C / C0

1582

0.6

IB1 IB2 IB3

0.4

0.2

0 0

100

200

300

400

x (cm) Fig. 3. Comparison of the spatial distribution of concentrations with the first-kind condition (IB1), the third-kind condition with an integral-order derivative (IB2), and the third-kind condition with a fractional-order derivative (IB3), respectively as the inlet boundary, and a zero concentration gradient with an integral derivative as the outlet boundary condition for a step input source.

x¼0

and the fractional-order zero gradient outlet boundary condition can be written as

    a1 ! 1 þ b oa1 C 1b o C D D 2 oxa1 2 oðxÞa1

¼ 0:

ð19Þ

x¼L

We evaluate these boundary conditions at x0 + D and at xN  D, respectively. So the second term on the left-hand side of Eq. (18)

 C 1 C 0 2a D according to the Caputo can be simplified to D 1þb 2DCð3aÞ 2 fractional derivative. After applying the Caputo definition to the third term of Eq. (18), the fractional-order, third-kind inlet boundary condition can be rearranged as

ð1  GLÞC 0 þ ðGL  GUÞC 1 N1   X 2a þ GU C i ði  2Þ2a þ 2ði  1Þ2a  i i¼2

  þ GUC N ðN  1Þ2a  ðN  2Þ2a ¼ C in ;

ð20Þ

1a where GU ¼ 2CDð1bÞ D1a and GL ¼ 2Dð1þbÞ . ð3aÞv Cð3aÞv D With the same procedure, the second term on the left-hand side  C N C N1 D2a by using the Capof Eq. (19) can be simplified as D 1b 2DCð3aÞ 2 uto definition. The outlet boundary condition denoted by Eq. (19) consequently becomes

   GL C 0 ðN  1Þ2a  ðN  2Þ2a 

N2 X

  C i ðN  iÞ2a þ 2ðN  i  1Þ2a  ðN  i  2Þ2a

!

i¼1

 ð1  GL þ GUÞC N1 þ ð1 þ GUÞC N ¼ 0:

ð21Þ

3. Results and discussion

third-kind boundary with an integral-order derivative (IB2) as given by Eq. (17), and (3) the third-kind boundary with a fractional-order derivative (IB3) as given by Eq. (18). All the outlet boundary conditions used to obtain Fig. 3 are zero concentration gradient with an integral-order derivative as defined in Eq. (16). A step input was used with v = 0.599 cm/min, D = 0.3 cm1.6/min, a = 1.6, b = 0. The length of domain is 400 cm and the number of spatial nodes is 400. The steady-state is reached at approximately 1400 min. The same parameters of v, D, a , b, domain length, and number of nodes will be used in Figs. 4–11 if not otherwise specified. The dimensionless steady-state concentration should be unity in the domain, regardless of what kind of boundary conditions used. However, as shown in Fig. 3, except for the boundary locations, the simulated normalized concentrations with both the IB1 and IB2 boundaries cannot reach unity, whereas the simulated concentrations with the IB3 boundary can approach unity. This implies that the IB3 boundary is the only correct one whereas the IB1 and IB2 boundaries are incorrect because of their obvious violation of mass balance. Actually, the problems with the use of IB1 and IB2 have long been recognized in solute transport with ADE, as seen from the work of van Genuchten and Parker [45] and Parker and van Genuchten [36]. It has also been observed experimentally by Novakowski [32,33] and Schwartz et al. [41]. Therefore, it is not surprise to see the similar problems occur in the present FADE study. Fig. 4 shows the BTCs obtained with these three different boundary conditions and it indicates that the numerical schemes with the IB1 and IB2 boundaries yield asymptotic values smaller than unity, while the numerical scheme with the IB3 boundary predicts an asymptotic value that is very close to unity in Fig. 4. In conclusion, the IB3 condition should be used for the inlet boundary.

3.1. Effect of boundary conditions 3.1.1. Effect of the inlet boundary condition Fig. 3 shows the effect of different inlet boundary conditions on solute spatial distribution under steady-state condition. Three kinds of inlet boundary conditions used to obtain this figure are: (1) the first-kind boundary (IB1) as given by Eq. (15), (2) the

3.1.2. Effect of the outlet boundary condition We also conducted simulations to investigate the effect of outlet boundary condition on the spatial distribution of concentration. The inlet boundary condition used in this simulation is the IB3 boundary. Two different types of outlet boundary conditions, i.e. the zero gradient concentration with an integral-order derivative

1583

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589 1.2

1.2

1

1

IB1 IB2 IB3

0.8

C / C0

C / C0

0.8

0.6

0.6

x =100 0.4

x=200

x=300

x=400

0.4

OB1 OB2 0.2

0.2

0

0 0

200

400

600

800

0

1000

200

400

Time (h)

1.2

1.2

1

1

0.8

0.8

OB1 OB2

0.6

0.4

0.2

0.2

0 200

1000

300

400

x (cm) Fig. 5. Comparison of the spatial distribution of concentrations with a zero gradient concentration with an integral derivative (OB1), and a zero gradient concentration with a fractional derivative (OB2) as the outlet boundary conditions, respectively, and the third-kind condition with a fractional derivative (IB3) as the inlet boundary for a step input.

(OB1), and the zero gradient concentration with a fractional-order derivative (OB2), are considered in the simulation. As shown in Fig. 5, except for the boundary locations, the spatial distribution of concentration is nearly the same for the two outlet boundary conditions at steady-state. This implies that the outlet boundary condition only have slight effect on the simulation result at large times. Fig. 6 shows the simulated BTCs at the distances of 100 cm, 200 cm, 300 cm, and 400 cm with the IB3 inlet boundary and the OB1 and OB2 outlet boundaries. The normalized concentration at the outlet point x = 400 cm is systematically lower than unity, whereas the normalized concentrations at the distances of 100 cm, 200 cm and 300 cm seem reasonable. This means that the FEM solution with the IB3 inlet boundary and either the OB1

Original Improved

0.6

0.4

100

800

Fig. 6. Comparison of the breakthrough curves at 100 cm, 200 cm, 300 cm and 400 cm from the left boundary with the zero gradient concentration with an integral derivative (OB1) and the zero gradient concentration with a fractional derivative (OB2) as the outlet boundary conditions, respectively, and the third-kind condition with a fractional derivative (IB3) as the inlet boundary for a step input source.

C / C0

C / C0

Fig. 4. Comparison of the breakthrough curves at 400 cm (at the right boundary) with the first-kind condition (IB1), the third-kind condition with an integral-order derivative (IB2), and the third-kind condition with a fractional-order derivative (IB3) conditions, respectively as the inlet boundary, and a zero concentration gradient with an integral derivative as the outlet boundary condition for a step input source.

0

600

Time (h)

0 0

100

200

300

400

x (cm) Fig. 7. Comparison of the difference of spatial distribution of concentration at steady-state between the original FEM and the improved FEM by the theory of transition probability conservation.

or OB2 outlet boundary is acceptable for simulating solute transport at locations other than the outlet location. As shown in Figs. 3 and 5, for the steady state situations, at the inlet boundary locations the concentration is larger than unity, while at the outlet boundary locations the concentration is small than unity (see the breakthrough curve at x = 400 in Fig. 4). The reason for the over and under estimation problems associated with BTCs at both the inlet and outlet is caused by the singularity of the fractional derivatives at the boundary. Actually, there is an Abel kernel (1/(x  f)a1) in the Caputo definition of fractional derivative (as shown in Eq. (4)), the integral expressions of the fractional derivative are weak singular when x is close to zero or L [12]. The weak singularity becomes worse at two boundary intervals [0, x0] and [xN1, xN] when D ? 0. This then leads to that the concentrations are over estimated at the inlet boundary and under estimated at the outlet boundary by using the FADE model.

1584

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589 1.2

a

0.6

1

0.8

0.6

x=100

x=200

x=300

Alpha=2.0 Alpha=1.8 Alpha=1.5 Alpha=1.2

C / C0

C / C0

0.4

x=400

0.4 0.2

N=200 N=400 N=800

0.2

0 0

200

400

600

800

0

1000

-200

Time (h)

0

200

x-vt (cm)

Fig. 8. Comparison of the breakthrough curves at 100 cm, 200 cm, 300 cm and 400 cm from the left boundary for different numbers of spatial nodes (N) 200, 400, 800, with the zero gradient concentration with a fractional derivative (OB2) as the outlet boundary conditions, respectively, and the third-kind condition with a fractional derivative (IB3) as the inlet boundary.

b

1

Alpha=2.0 Alpha=1.8 Alpha=1.5 Alpha=1.2

0.5

1.2

C / C0

0.2

0.1

1

0.05

C / C0

0.8

0.6

FEM N=800 FEM N=200 FD N=800 FD N=200 Analytical

0.4

0.01 0.005 -200

0

200

x-vt (cm)

0.2

Fig. 10. Impact of a on the breakthrough curves at 300 cm from the left boundary with v = 0.599 cm/min, D = 0.3 cma/min, b = 0, (a) linear axes, and (b) semi-log axes.

0

200

400

600

800

1000

Time (h) Fig. 9. Comparison of the analytical solution with the finite element solution at 300 cm from the left boundary of a 400 cm long domain with v = 0.599 cm/min, D = 0.3 cm1.5/min, b = 0, a = 1.5, and 200 and 800 spatial nodes (N).

3.2. Alternative method for evaluation of concentration at the boundaries Gorenflo et al. [14,15] used random walk method to solve the fractional diffusion equation numerically. They first used a finite difference approach to discretize the fractional diffusion equation in space and time. Then the weighting coefficients of the FD scheme are considered as the transition probability of particles undergoing random walk. Unlike the integral diffusion, the fractional diffusion has the transition probability of jump from one lattice point to others not only occurring in its nearby lattice points, but also occurring in all lattice points of the spatial domain [14,15,40]. It has been proven by Schumer et al. [40] and Gorenflo et al. [14,15] that for one sided fractional derivative (e.g. b = 1 or b = 1), the summation of the transition probability of particles from one lattice point to others should be equal to unity, i.e. the

summation of the weighting coefficients of the FD or FEM should be equal to unity for one-sided fractional derivative. The FADE considered in this paper is a two-sided fractional derivative with weight (1 + b)/2 for the right-handed derivative and (1  b)/2 for the left-handed derivative, so the summation of the weighting coefficients of the FD or FEM should be also equal to unity. The linear system defined in Eq. (14) together with the boundary conditions can be expressed by a linear equation group of ACk+1 = Ck, in which A is the weighting matrix with elements of Ai,j, i,j = 1, 2, . . . , N, and N is the dimension of matrix, C is the concentration vector, and k represents the k-th time interval. Thus for j – 0 and N, one can obtain N X

Ai;j ¼ 1:

ð22Þ

i¼0

The weak singularity may lead to noticeable deviation of the left hand side of Eq. (22) from unity. To move around the problem, it is recommended not to evaluate the integration near the singular points directly. Instead, one can force the left hand side of Eq.  (22) to be unity to calculate the updated W þ 0 and W N , denoted as 0 0 W 0 and W N , respectively. Let b = 1, one can obtain

1585

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589

a

b

0.4

Beta=1.0 Beta=0.0 Beta=-1.0

1

Alpha=1.6 Beta=1 Alpha=2.0

0.5

0.2

C / C0

C / C0

0.1 0.05

0.2

0.01

NonFickian

Fickian

0.005

0.001 0

200

400

600

800

-200

1000

c

1

Alpha=1.6 Beta=0.0 Alpha=2.0

0.5

d

1

Alpha=1.6 Beta=-1.0 Alpha=2.0

0.5

0.2

0.1

0.1

C / C0

C / C0

0.2

0.05

0.05

0.01

0.01

0.005

0.005

NonFickian

NonFickian 0.001

NonFickian

Fickian

0.001 -200

200

0

x-vt (cm)

Time (h)

0

-200

200

0

200

x-vt (cm)

x-vt (cm)

Fig. 11. Impact of b on the breakthrough curves at 300 cm from the left boundary with v = 0.599 cm/min, D = 0.3 cm1.6/min, a = 1.6, and (a) b = 1, 0, 1, (b) b = 1 with comparing the result of a = 2, (c) b = 0 with comparing the result of a = 2, and (d) b = 1, with comparing the result of a = 2, respectively.

W 00 ¼ 1 

j2 X

  W þi  6 þ 4ð2Þ3a  ð3Þ3a þ ð4  ð2Þ3a Þ:

ð23Þ

i¼1

AC = kC. Assuming Ci is the maximum absolute values of the elePN ment in vector C, then one has j¼0 Ai;j C j ¼ kC i , the eigenvalue k can be expressed as

Let b = 1, one can obtain

W 0N ¼ 1 

N1 X





W þi  6 þ 4ð2Þ3a  ð3Þ3a þ ð4  ð2Þ3a Þ:

k ¼ Ai;i þ ð24Þ

i¼jþ2

As D ? 0, it can be easily verified that distinct difference between W 00 and W þ 0 occurs when j is less than 3, and notable difference between W 0N and W  N also occurs when j is greater than N  3. Fig. 7 compares the difference of spatial distribution of concentration at steady-state between the original FEM and the improved FEM by the theory of transition probability conservation with IB3 as the inlet boundary condition. It is distinct to see from Fig. 7 that the improper estimations at the boundary points in the original FEM scheme are remedied by the improved FEM scheme. 3.3. Stability analysis We used the method as in Meerschaert and Tadjeran [25] to evaluate the stability of the proposed FEM scheme. Let k be an eigenvalue of the matrix A, and C is a nonzero vector, then we have

N Cj X Ai;j ; C i j¼0;j–i

ð25Þ

where Ai,i is the ith diagonal element of the matrix A. From Eq. (22), one can obtain

1þ

s D

WB þ

N X

Ai;j ¼ 1

for j – 0 and N. Eq. (26) means DD1a ð4  ð2Þ3a Þ > 0, 2Cð4aÞ s WB þ C j PN j¼0;j–i Ai;j > 0. D Ci

k¼1þ

ð26Þ

i¼0;i–j

s D

WB þ

s WB þ PN

D

then

PN

j¼0;j–i Ai;j

¼ 0.

j¼0;j–i Ai;j <0.

Because Since

Cj Ci

WB ¼

6 1,

so

Thus one has

N Cj X Ai;j P 1 C i j¼0;j–i

ð27Þ

for i – 0 and N . From Section 2.3, it also can be found that k P 1 for i = 0 or N.

1586

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589

As shown in Eq. (27), we conclude that every eigenvalue of the matrix A satisfies jkj P 1. So A is invertible, and every eigenvalue (g) of the inverse matrix A1 satisfies jgj 6 1. Therefore, the spectral radius of the inverse matrix q(A1) 6 1, an error e0 in the initial concentration vector C0 will result in an error e1 in the concentration vector C1, and is given by e1 = A1e0, so the error is bounded ke1k 6 ke0k. Thus the proposed FEM scheme for the FADE is unconditionally stable consequently. Numerical experiments have also been carried out to investigate the stability of the FEM solution for the FADE. Fig. 8 shows the simulated BTCs at the distances of 100 cm, 200 cm, 300 cm, and 400 cm, respectively, by using the improved FEM scheme. The number of nodes for temporal discretization is 900. The boundary conditions as shown in Eqs. (18) and (19) are adopted in the numerical simulations. Three different nodes numbers 200, 400, 800 are used for discretizing the space domain in order to see the sensitivity of results to the number of nodes. From Fig. 8, it can be found that the BTCs are almost the same for different nodes number. This again implies that the proposed FEM scheme for the FADE is stable considering the nodes number. Furthermore, Fig. 8 shows a correctly prediction of BTCs at the outlet point for the improved FEM scheme comparing to Fig. 6. The accuracy of the FEM scheme is also evaluated by comparing the numerical results with the analytical solution [2] of the FADE, and their difference is characterized with the numerical error [23].

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N N uX X dC ¼ t ðC n ðti Þ  Cðti ÞÞ2 = Cðti Þ2 ; i¼0

ð28Þ

i¼0

where dC is the numerical error, Cn(ti) is the numerical solution of Eq. (16) at time ti, and C(ti) is the analytical solution of the symmetrical FADE at ti. Three a values of 1.3, 1.6 and 1.9 and five node numbers of 100, 200, 400, 800 and 1600 are used to evaluate the accuracy of the FEM solution. All the other parameters remain the same as those used for obtaining Fig. 8. Table 1 lists the results of the accuracy analysis. The numerical error decreases with the increasing number of the spatial nodes. The numerical errors for the case with a equaling 1.3 are larger than those of other two cases with a values equaling 1.6 and 1.9. This means that the convergence of the FEM solution is sensitive to the a value, and the FEM solution converges more rapidly as a value becomes greater. Similar to those numerical schemes of Zhang et al. [52], Deng et al. [11], and Lynch et al. [23], the FEM solution is close to a second-order approximation when a approaches 2, and is close to a first-order approximation when a approaches 1, respectively. It is well known that the numerical schemes of the ADE are subject to numerical dispersion [19]. For investigating the numerical dispersion of the FADE, the comparison between the analytical solution [2] and the numerical solutions of the FADE at the distance of 300 cm from the left boundary of 400 cm long domain

Table 1 Errors between the analytical solution and the numerical solution on different a values for the nodes number changing from N = 100 to N = 1600 Nodes number

dC

a = 1.3

a = 1.6

a = 1.9

100 200 400 800 1600

0.150470 0.128019 0.120889 0.117267 0.115142

0.062985 0.048299 0.042324 0.039406 0.037838

0.046403 0.033317 0.028047 0.025869 0.024899

The concentration is calculated at the 300 cm from the left boundary of a 400 cm long domain. v = 0.599, D = 0.3 cma/min, b = 0, and the number of nodes for temporal discretization is 900.

with a = 1.5 is shown in Fig. 9 as an example. The values of v, D, and b used in Fig. 9 are the same as those used in Fig. 3. In this figure, the FD solution of the FADE [25,26] and the FEM solution of this study for two different numbers of spatial nodes (N = 200 and 800) are included. As can be seen from this figure, the numerical dispersion exists in both the FD and FEM solutions but will be reduced when the number of nodes increases. The FEM solution of this study has smaller numerical dispersion than that of the FD solution [25,26] at a given number of nodes. However, the reduction of the numerical dispersion for the FEM solution is not as significant as that for the FD solution when the number of nodes increases. To completely remove the numerical dispersion, one has to use a better numerical solver such as the method of characteristics (MOC) and its several improved forms or the total variation diminishing (TVD) which has been excellently explained by Zheng and Bennett [54] and many other investigators. 3.4. Effect of the fractional order a The parameter a exhibits the extent of non-Fickian behavior of the BTCs. Fig. 10 shows the simulated results for an instantaneous source input with v = 0.599 cm/min, D = 0.3 cma/min, b = 0, and a values of 1.3, 1.5, 1.8 and 2.0, respectively. It can be found that the BTC follows a Gaussian distribution for the case of a = 2. The BTCs exhibit non-Fickian behavior with heavy leading edges and tails and relatively large dispersions for a smaller than 2, and the extent of non-Fickian behavior increases with the decreasing a value. Similar results were also found by Benson et al. [2] and Schumer et al. [40]. 3.5. Effect of the skewness parameter b The parameter b reflects the relative weight of forward versus backward of non-Gaussian transition probability. Fig. 11a shows the BTCs with an instantaneous source input for b = 1, 0, 1, respectively. The parameters used in the numerical experiment are v = 0.599 cm/min, D = 0.3 cm1.6/min, and a = 1.6. It can be found that the model predicts the BTCs with a both heavy leading edge and tail for b = 0, a heavy leading edge and a light tail for b = 1, and a light leading edge and a heavy tail for b = 1, respectively. Comparing to the Fickian transport, the BTC for b = 1 shows a more notable non-Fickian behavior on the right-hand side, and an asymptotic Fickian behavior on the left-hand side, as shown in Fig. 11b. The BTC for b = 0 shows a remarkable non-Fickian evolving shape on both sides of the spatial domain, as shown in Fig. 11c. The BTC for b = 1 shows a more notable non-Fickian behavior on the left-hand side, and an asymptotic Fickian behavior on the righthand side, as shown in Fig. 11d. 3.6. Model application A series of experiments have been conducted by Mao and Ren [24] to investigate the adsorption property and transport process of atrazine in homogenous sandy loamy soil column. The soil column is 20 cm in length and 5.0 cm in diameter. The pore velocity of the experiment is 1.925 cm/h when the flow reaches steady state. During a 1 h pulse input, the atrazine with a concentration of 220.5 mg/L was applied to the inlet of the soil column under steady-state flow condition. The BTCs of atrzine show a remarkable retardation followed by a heavy tail, as shown in Fig. 12. More than 98% of the atrazine was recovered during the experiment. The isotherm batch experiment [24] showed that the retardation factor (R) of the atrazine in sandy loamy soil is 1.825. Both numerical solutions of the ADE and the FADE were used to fit the transport process of the atrazine. We used the following root

1587

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589 25

Measured ADE FADE

20

mg / L

15

10

5

0

10

20

30

40

(2) The simulation results can be improved by using the thirdkind boundary with a fractional-order derivative as the inlet boundary condition. The simulation result is insensitive to the choice of either a fractional-order or an integral-order derivative at the outlet boundary, provided that a zero concentration gradient (the Danckwerts’ condition) is applied there. (3) For a given a, the spatial distribution of concentration exhibits a symmetric non-Fickian behavior when b = 0. The spatial distribution of concentration shows an asymptotic Fickian behavior on the left-hand side of the spatial domain and a notable non-Fickian behavior on the right-hand side of the spatial domain when b = 1, whereas when b = 1 the spatial distribution of concentration is the opposite of that of b = 1. (4) Comparing to the numerical solution of the ADE, the FEM solution of the FADE can better simulate the atrazine transport process, especially at the tail of the breakthrough curve.

Time (h) Fig. 12. Comparison of the breakthrough curves estimated by the FADE, the ADE, and the measured data of atrazine at the end of the sand column.

mean square error (RMSE), mean relative error (MRE), and model efficiency (EF) to evaluate the model fitting quality

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X RMSE ¼ t ðPi  Oi Þ2 ; N i¼1 MRE ¼

N 1X ABSðPi  Oi Þ  100%; N i¼1 Oi N P

EF ¼ 1 

ð29Þ

ð30Þ

ðPi  Oi Þ2

i¼1 N P

;

ð31Þ

ðOi  OÞ2

i¼1

where N is the total number of observations, Oi and Pi are, respectively the observed and predicted values of the ith observation, and O is the mean of the observed values (i = 1 to N). The fitted results for the ADE are v = 1.900 cm/min, D = 1.304 cm2/min, and R = 1.825, with RMSE = 1.29, MRE = 69.44 and EF = 0.969, respectively. Whereas the fitted results for the FADE are v = 1.920 cm/min, D = 1.304 cm1.6/min, a = 1.6, b = 1, and R = 1.825, with RMSE = 0.75, MRE = 42.6 and EF = 0.990, respectively. The values of RMSE and MRE for the FADE fitting are lower than those for the ADE, and the EF value for the FADE is higher than that for the ADE. This indicates that the FADE can give a better fitting to the BTC of the atrazine than the ADE. Especially, the FADE simulates the tail better (shown in Fig. 12). 4. Conclusions We developed a FEM scheme to solve the FADE based on the Caputo definition of the fractional derivative coupling with considering its singularity at boundaries. The stability and accuracy of the FEM solution was examined, and the sensitivity of the FEM solution to the fractional order a and the skewness parameter b was then analyzed. Finally, the numerical approach was applied to simulate the atrazine transport in a saturated soil column. Several conclusions can be drawn from this study: (1) The proposed FEM solution of the FADE is unconditionally stable. The problems for concentration estimation at boundaries caused by the singularity of the fractional derivative can be solved by using the concept of transition probability conservation.

Acknowledgements This research was partly supported by the National Natural Science Foundation of China (Grant numbers: 50639040, 50779067, 50479011), the Program for New Century Excellent Talents in University (Grant number: NCET-05-0125), the Program for Changjiang Scholars and Innovative Research Team in University (IRT0657) and the Program of Beijing Key Subject of Hydrology and Water Resources. The authors acknowledge Prof. Li Ren and Dr. Meng Mao for providing the laboratory experimental data. The constructive comments from two anonymous reviewers, the Editor, and Dr. Xiaoxian Zhang are also gratefully acknowledged, which help us improve the quality of the manuscript. Appendix Numerical approximation for the third and forth terms on the left-hand side of Eq. (7) is listed as following. Applying the method of integration by parts and making use of the property of the linear basis function, one can easily obtain.

F3 ¼

Dð1 þ bÞ 2Cð2  aÞ

Z

xjþ1

Z

xj1

0

x

 b 0 ðfÞdf /0 ðxÞdx: ðx  fÞ1a C j

ðA1Þ

Eq. (A1) can further be expressed as

(Z j1 Z xi xj X Dð1 þ bÞ b 0 ðfÞdf F3 ¼ ðx  fÞ1a C 2Cð2  aÞ xj1 xi1 i¼1 ! Z x 1a b 0 þ ðx  fÞ C ðfÞdf /0j ðxÞdx xj1

þ

Z

xjþ1 xj

j Z X i¼1

)

 /0j ðxÞdx :

xi

xi1

ðx  fÞ

1a b 0

C ðfÞdf þ

Z

!

x

ðx  fÞ

1a b 0

C ðfÞdf

xj

ðA2Þ

b 0 ðfÞ ¼ C i1 u0 þ C i u0 for f 2 [xi1,xi] into Eq. (A2), Substituting C i1 i and with a uniform spatial discretizationD, one has

(Z Z j1 xj X Dð1 þ bÞ C i  C i1 xxi1 1a F3 ¼ t dt 2Cð2  aÞD xj1 i¼1 D xxi  Z xxj1 C j  C j1 t 1a dt dx þ D 0 ! ) Z Z Z xjþ1 X j C i  C i1 xxi1 1a C jþ1  C j xxj 1a  t dt þ t dt dx : D D xj xxi 0 i¼1 ðA3Þ

1588

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589

After integration, one obtains

F3 ¼

Dð1 þ bÞ 2Cð4  aÞD

( j1  X ðC i  C i1 Þ ðxj  xi1 Þ3a  ðxj1  xi1 Þ3a 2

ðxj  xi Þ

þ ðxj1  xi Þ

3a



þ ðC j  C j1 ÞD

)  ðxjþ1  xi Þ3a þ ðxj  xi Þ3a  ðC jþ1  C j ÞD3a :

ðA4Þ

Eq. (A4) can be expressed in a simplified form for easy programming

2Cð4  aÞD

n  C 0 ðxjþ1  3DÞ3a þ 3ðxjþ1  2DÞ3a 2

3ðxjþ1  DÞ3a þ ðxjþ1 Þ3a 

j2 X



 C i ðxjþ1  xiþ3 Þ3a þ 4ðxjþ1  xiþ2 Þ3a  6ðxjþ1  xiþ1 Þ3a

i¼1

 þ4ðxjþ1  xi Þ3a  ðxjþ1  xi1 Þ3a     þC j1 6ðDÞ3a þ 4ð2DÞ3a  ð3DÞ3a þ C j 4D3a  ð2DÞ3a o ðA5Þ þC jþ1 ðD3a Þ : Numerical approximation for the forth term on the left-hand side of Eq. (7) can be expressed as

F4 ¼

Z

Dð1  bÞ 2Cð2  aÞ

xjþ1

Z

L

x

xj1

 b 0 ðfÞdf /0 ðxÞdx: ðf  xÞ1a C j

ðA6Þ

b 0 ðfÞ ¼ C i1 u0 þ C i u0 for f 2 [xi1,xi], Eq. (6) becomes With C i1 i

F4 ¼

Dð1bÞ 2Cð2 aÞD

þ

N1 X C iþ1 C i i¼j



(Z

Z

xjþ1

xj

 Z C j C j1 xj x 1a t dt D xj1 0 ! Z xj

xiþ1 x

D

t 1a dt dx

xi x

C jþ1 C j D

Z

xjþ1 x

t1a dt þ

0

Z N1 X C iþ1 C i i¼jþ1

D

xiþ1 x

!

)

t 1a dt dx :

xi x

ðA7Þ After integration of Eq. (A7), one has

F4 ¼

(

Dð1  bÞ 2Cð4  aÞD

3a

ðxi  xj Þ

ðC j  C j1 ÞD3a 

2

j1  X ðC i  C i1 Þ ðxiþ1  xj Þ3a i¼1

3a

 ðxiþ1  xj1 Þ

ðC jþ1  C j ÞD3a þ

j X

þ ðxi  xj1 Þ3a

ðxiþ1  xj Þ



 ðC i  C i1 Þ ðxiþ1  xjþ1 Þ3a

i¼1 3a

3a

 ðxi  xjþ1 Þ

þ ðxi  xj Þ

3a

n   C j1 ðD3a Þ þ C j 4D3a  ð2DÞ3a

2Cð4  aÞD   þC jþ1 6ðDÞ3a þ 4ð2DÞ3a  ð3DÞ3a j2 X

 C i ðxiþ1  xjþ3 Þ3a þ 4ðxiþ1  xjþ2 Þ3a  6ðxiþ1  xjþ1 Þ3a

i¼1

 ðC i  C i1 Þ ðxjþ1  xi1 Þ3a  ðxj  xi1 Þ3a

j X

Dð1 þ bÞ

2

þ

3a

i¼1

F3 ¼

Dð1  bÞ

i¼1

3a



F4 ¼



) :

ðA8Þ

Eq. (A8) can also be expressed in a simplified form for easy programming

 þ4ðxiþ1  xj Þ3a  ðxiþ1  xj1 Þ3a  þC N ðxN  xj  2DÞ3a þ 3ðxN  xj  DÞ3a  3ðxN  xj Þ3a )  3a : ðA9Þ þðxN  xj þ DÞ

References [1] Bear J. Hydraulics of groundwater. New York, USA: McGraw-Hill; 1979. [2] Benson DA, Wheatcraft SW, Meerschaert MM. The fractional-order governing equation of Lévy motion. Water Resour Res 2000;36(6):1413–23. [3] Benson DA, Wheatcraft SW, Meerschaert MM. Application of a fractional advective–dispersion equation. Water Resour Res 2000;36(6):1403–12. [4] Berkowitz B. Characterizing flow and transport in fractured geological media: a review. Adv Water Resour 2002;25(8–12):861–84. [5] Berkowitz B, Cortis A, Dentz M, Scher H. Modeling non-Fickian transport in geological formations as a continuous time random walk. Rev Geophys 2006;44(2):RG2003. doi:10.1029/2005RG000178. [6] Butzer PL, Westphal U. An introduction to fractional calculus. Singapore: World Scientific; 2000. [7] Caputo M. Linear models of dissipation whose Q is almost frequency independent-II. Geophys J Roy Astr Soc 1967;13(5):529–39. [8] Caputo M, Mainardi F. A new dissipation model based on memory mechanism. Pure Appl Geophys 1971;91(8):134–47. [9] Coats KH, Smith BD. Dead-end pore volume and dispersion in porous media. Soc Petrol Eng J 1964;4:73–84. [10] Cushman JH, Ginn TR. Fractional advection–dispersion equation: A classical mass balance with convolution-Fickian flux. Water Resour Res 2000;36(12): 3763–6. [11] Deng Z, Singh VP, Bengtsson L. Numerical solution of fractional advection– dispersion equation. J Hydraul Eng 2004;130(5):422–31. [12] Diethelm K, Ford NJ, Freed AD, Luchko Y. Algorithms for the fractional calculus: a selection of numerical methods. Comput Methods Appl Mech Eng 2005;194: 743–73. [13] Ellsworth TR, Shouse PJ, Skaggs TH, Jobes JA, Fargerlund J. Solute transport in unsaturated soil: experimental design, parameter estimation, and model discrimination. Soil Sci Soc Am J 1996;60(2):397–407. [14] Gorenflo R, de Fabritiis G, Mainardi F. Discrete random walk models for symmetric Lévy–Feller diffusion process. Physica A 1999;269(1):79–89. [15] Gorenflo R, Mainardi F, Moretti D, Pagnini G, Paradisi P. Discrete random walk models for space-time fractional diffusion. Chem Phys 2002;284:521–41. [16] Haggerty R, Harvey CF, Freiherr von Schwerin C, Meigs LC. What controls the apparent timescale of solute mass transfer in aquifers and soils? A comparison of experimental results. Water Resour Res 2004;40:W01510. doi:10.1029/ 2002WR001716. [17] Heymans N, Podlubny I. Physical interpretation of initial conditions for fractional differential equations with Riemann–Liouville fractional derivatives. Rheol Acta 2006;45(5):765–71. [18] Huang G, Huang Q, Zhan H. Evidence of one-dimensional scale-dependent fractional advection–dispersion. J Contam Hydrol 2006;85:53–71. [19] Huyakorn PS, Pinder GF. Computational methods in subsurface flow. New York, USA: Academic Press; 1983. [20] Kilbas A, Srivastava HM, Trujillo JJ. Theory and applications of fractional differential equations. North-Holland mathematics studies 204. Boston, USA: Elsevier; 2006. [21] Li C, Deng W. Remarks on fractional derivatives. Appl Math Comput 2007;187(2):777–87. [22] Liu F, Anh V, Turner I. Numerical solution of the space fractional Fokker–Planck equation. J Comput Appl Math 2004;166(1):209–19. [23] Lynch VE, Carreras BA, del-Castillo-Negrete D, Ferreira-Mejias KM, Hicks HR. Numerical methods for the solution of partial differential equations of fractional order. J Comput Phys 2003;192(2):406–21. [24] Mao M, Ren L. Simulating nonequilibrium transport of atrazine through saturated soil. Ground Water 2004;42(4):500–8. [25] Meerschaert MM, Tadjeran C. Finite difference approximations for fractional advection–dispersion flow equations. J Comput Appl Math 2004;172(1): 65–77. [26] Meerschaert MM, Tadjeran C. Finite difference approximations for two-sided space-fractional partial differential equations. Appl Numer Math 2006;56(1): 80–90. [27] Metzler R, Klafter J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys Rep 2000;339(1):1–77.

Q. Huang et al. / Advances in Water Resources 31 (2008) 1578–1589 [28] Neuman SP. Eulerian–Lagrangian theory of transport in spacetime nonstationary velocity fields: exact nonlocal formalism by conditional moments and weak approximation. Water Resour Res 1993;29(3):633–45. [29] Miller KS, Ross B. An introduction to the fractional calculus and fractional differential equations. New York, USA: Wiley; 1993. [30] Morales Casique E, Neuman SP, Guadagnini A. Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: theoretical framework. Adv Water Resour 2006;29:1238–55. [31] Morales Casique E, Neuman SP, Guadagnini A. Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: computational analysis. Adv Water Resour 2006;29:1399–418. [32] Novakowski K. An evaluation of boundary conditions for one-dimensional solute transport 1. Math Develop Water Resour Res 1992;28(9):2399–410. [33] Novakowski K. An evaluation of boundary conditions for one-dimensional solute transport 2. Column Exp Water Resour Res 1992;28(9):2411–23. [34] Pachepsky Y, Benson DA, Rawls W. Simulating scale-dependent solute transport in soils with the fractional advective–dispersive equation. Soil Sci Soc Am J 2000;64:1234–43. [35] Pang L, Hunt B. Solutions and verification of a scale-dependent dispersion model. J Contam Hydrol 2001;53:21–39. [36] Parker J, Van Genuchten M. Flux-averaged and volume-averaged concentrations in continuum approaches to solute transport. Water Resour Res 1984;20(7):866–72. [37] Podlubny I. Fractional differential equations. New York, USA: Academic Press; 1999. [38] Roop JP. Computational aspects of FEM approximation of fractional advection– dispersion equations on bounded domains in R2. J Comput Appl Math 2006;193(1):243–68. [39] Schumer R, Benson DA, Meerschaert MM, Baeumer B. Fractal mobile/immobile solute transport. Water Resour Res 2003;39(10):1296. doi:10.1029/ 2003WR002141. [40] Schumer R, Benson DA, Meerschaert MM, Wheatcraft SW. Eulerian derivation of the fractional advection–dispersion equation. J Contam Hydrol 2001;48:69–88. [41] Schwartz R, McInnes K, Juo A, Wilding L, Reddell D. Boundary effects on solute transport in finite soil columns. Water Resour Res 1999;35(3):671–81. [42] Šimu˚nek J, Huang K, van Genuchten MTh. The HYDRUS code for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-

[43] [44] [45]

[46]

[47]

[48] [49]

[50]

[51]

[52]

[53]

[54]

1589

saturated media. Version 6.0, Research Rep. 144. US Salinity Laboratory, Riverside, CA, USA; 1998. Tsallis C, Lenzi EK. Anomalous diffusion: non-linear fractional Fokker–Plancker equation. Chem Phys 2002;284:341–7. van Genuchten MT, Wierenga PJ. Mass transfer studies in sorbing porous media: 1. Anal Solutions Soil Sci Soc Am J 1976;40(4):473–80. van Genuchten MTh, Parker JC. Boundary-conditions for displacement experiments through short laboratory columns. Soil Sci Soc Am J 1984;48: 703–8. Yeh GT. A Lagrangian–Eulerian Method with Zoomable hidden fine-mesh approach to solving advection–dispersion equations. Water Resour Res 1990;26(6):1133–44. Zhang H, Liu F, Anh V. Numerical approximation of Lévy–Feller diffusion equation and its probability interpretation. J Comput Appl Math 2007;206(2): 1098–115. Zhang R, Huang K, Xiang J. Solute movement through homogeneous and heterogeneous soil columns. Adv Water Resour 1994;17(5):317–24. Zhang Y, Benson DA, Meerschaert MM, Scheffler HP. On using random walks to solve the space-fractional advection–dispersion equations. J Stat Phys 2006;123(1). doi:10.1007/s10955-006-9042-x. Zhang Y, Benson DA, Meerschaert MM, LaBolle EM, Scheffler HP. Random walk approximation of fractional-order multiscaling anomalous diffusion. Phys Rev E 2006;74:026706. Zhang Y, Benson DA, Meerschaert MM, LaBolle EM. Space-fractional advection–dispersion equations with variable parameters: diverse formulas, numerical solutions, and application to the macrodispersion experiment site data. Water Resour Res 2007;43:W05439. doi:10.1029/2006WR004912. Zhang X, Crawford JW, Deeks LK, Stutter MI, Bengough AG, Young IM. A mass balance based numerical method for the fractional advection– dispersion equation: theory and application. Water Resour Res 2005;41(7):W07029. Zhang X, Lv M, Crawford JW, Young IM. The impact of boundary on the fractional advection–dispersion equation for solute transport in soil: defining the fractional dispersive flux with the Caputo derivatives. Adv Water Resour 2007;30(5):1205–17. Zheng C, Bennett GD. Applied contaminant transport modeling. 2nd ed. New York, USA: John Wiley & Sons; 2002.