The Monte Carlo and fractional kinetics approaches to the underground anomalous subdiffusion of contaminants
1 The Monte Carlo and Fractional Kinetics approaches to the underground anomalous subdiffusion of contaminants
M. Marseguerra , A. Zoia Department of Nuclear Engineering, Polytechnic of Milan, Via Ponzio 34/3, 20133 Milan, Italy
Abstract
It is nowadays recognized that the experimental evidences of many transport phenomena must be interpreted in the framework of the so-called anomalous diffusion: this is the case e.g. of the spread of contaminant particles in porous media, whose mean squared displacement (MSD) has been experimentally reported to grow in time as t a (with ). This is in deep contrast with the linear increase which characterizes the standard Fickian diffusive processes. Anomalous transport is currently tackled within the Fractional Diffusion Equation (FDE) analytical model, which has blown new life into the concept of fractional derivatives, dormant since Leibniz’s first discovery, about 300 years ago. In this paper we focus on subdiffusion, i.e. the case a<1 , which provides a suitable explanation of the observed non-Fickian persistence of the contaminant particles near the source in various transport experiences in heterogeneous media. Unfortunately, the FDE approach requires several approximations to overcome its analytical complexities. To evaluate the relevance of these approximations, we propound the use of the Monte Carlo simulation as a suitable reference for the analytical FDE results. This comparison shows that the FDE results are less and less reliable as a becomes closer to . The approach herein proposed can be easily applied to various fields of science where anomalous diffusion phenomena may occur, from physics and chemistry to biology and finance.
1. Introduction
About eighty years ago, Richardson, investigating the diffusion of a pair of fluid elements in fully developed turbulence, observed the existence of transport phenomena characterized by a mean square displacement (MSD) proportional not to the time t (as in the usual diffusion equation established by Fick in 1855 (Fick 1855)), but instead to the third power of t (Richardson 1926). This observation opened in the late 1960 the route to the investigations on a completely new kind of diffusion, in the framework of the Fractional Kinetics: nowadays, transport phenomena characterized by a MSD proportional to t α , with „a , are called anomalous diffusion processes, even if they are so ubiquitous that it might be reasonably said that they represent the norm. If α < 1 the diffusion is slower than in the case of standard diffusion and it is called subdiffusion; if α > 1 the diffusion is faster and it is called superdiffusion. Since the Richardson’s work, a lot of experimental evidences concerning the existence and the exploitation of anomalous diffusions have been obtained in quite different fields, such as physics, chemistry, biology, economics, geology. Most of the theoretical interpretations have been established within the Fractional Kinetics approach, which generalizes the Fick’s second law, the Fokker-Planck equation and the Montroll and Weiss (Montroll, Weiss 1965) framework of the continuous time random walks (CTRW). Exhaustive references may be found e.g. in the fundamental papers by Metzler and Klafter (Metzler, Klafter 2000; Metzler, Klafter 2004). Just to give an idea of the possible applications, we explicitly mention only some important results, e.g. the subdiffusive charge transport in amorphous semiconductors (Scher, Montroll 1975), the anomalous diffusion properties of heat channels (Denisov et al. 2003), the chemical reaction processes in high- k dielectric films (de Almeida, Baumvol 2003), the motion of DNA-binding proteins along the DNA structure (Slutsky et al. 2003), the waiting time between two financial transactions (Mainardi et al. 2000), the contaminant transport in geological formations (Kosakowski 2004). Here we focus on this last subject: the underground migration of contaminants driven by the subsurface water flow, where the medium heterogeneities over many scales and the trapping of the contaminant particles in the solid environment support the occurrence of subdiffusion phenomena. The importance of this problem is well recognized Corresponding Author. Tel: +39 02 2399 6355 . Fax: +39 02 2399 6309 Email address: [email protected] (Marzio Marseguerra) Email address: [email protected] (Andrea Zoia) particularly in view of the possibility of leakage and successive migration of long-term radioactive or toxic species from underground repositories. We tackle this problem within the Monte Carlo approach, whose results are compared to the analytical ones of the fractional kinetics. In Section 2 we summarize the Fractional Kinetics approach to subdiffusion essentially along the lines suggested by the CTRW model of a generalized walker and detail the various approximations. In Section 3 we present the potentialities of the Monte Carlo method in simulating subdiffusive walkers in various environments and conditions and describe the power law pdf from which we sample the rest times. In Section 4 we report the main results of the Monte Carlo simulations of a purely subdiffusive process and compare these results to the analytical ones deriving from the Fractional Kinetics approach. In Section 5 we show one of the possible extensions of the Monte Carlo approach: the simulation of an advection-dispersion process in a heterogeneous medium composed of two zones with different physical properties. In Section 6 we finally draw some conclusion of this work. Several appendices make this paper self-consistent.
2. Brief review of Fractional Kinetics
Let us consider the continuous time random walk (CTRW) of a particle in a medium, consisting of a sequence of jump-trapping-release processes: more specifically, the transport will be described in terms of the random walk of a walker which performs a succession of instantaneous jumps of length x (drawn from a probability density function (pdf) λ (x) ), followed by time intervals t (drawn from a pdf w(t) ) during which it is trapped at the location reached in the last jump. The assumption that these two probabilities do not depend on the past fate of the walker implies that in any case the process is semi-Markovian, i.e. Markovian at the jump instants. Better specification depends on the choice of w(t) and of l (t) . Physical examples of such motion may be found e.g. in the underground transport of carriers of toxic or radioactive contaminants, subject to immobilizations in and releases from the solid environment, caused by the subsurface water flow (Kosakowski 2004), or in the transport of charge carriers in amorphous semiconductors under a constant electric field (Scher, Montroll 1975). Clearly, assuming linearity – i.e. absence of interactions between the walkers –, this phenomenological picture lends itself to be Monte Carlo (MC) simulated once the pdf’s λ (x) and w(t) and the medium parameters are assigned. However, before describing the steps of the MC simulation, we prefer to summarize the analytical approach (Metzler, Klafter 2004; Metzler, Klafter 2000; Montroll, Weiss 1965), which enlightens many features of the underlying transport process. The process may be formalized by introducing firstly a bivariate probability p(x,t)dxdt that the walker arrives at a location x at time t (actually at a location between x and x+dx at a time between t and t+dt ) and then another univariate probability P(x,t)dx that the walker is at x (more precisely, between x and x+dx ) at time t , independently of when it has arrived at x . The pdf p(x,t) may be expressed by a Chapman-Kolmogorov type integral equation in terms of the probabilities that the walker i) arrived at any location ),(' ¥-¥˛ x at any previous time tt £ ' and ii) remained immobilized at ' x for a time interval ' tt - and, finally, iii) performed an instantaneous jump from x’ to x : )()()'()'()','(''),( txxxttwtxpdxdttxp t ddl ∫∫ ¥ ¥- +--= (1) where the second term at rhs represents the initial condition of the walk, namely that the walker started at x =0 at t =0 . The other pdf P(x,t) may now be expressed as the integral of the product of the probability that the walker arrived at x at any time tt £ ' , namely p(x,t’)dxdt’ , times the probability )'( tt -f that it remained immobilized there at least for t-t’ , viz., ∫ -= t tttxpdttxP )'()',('),( f (2) where ∫ -= t twdtt )'('1)( f (3) An explicit expression for P(x,t) might be obtained by firstly performing the Laplace- )( ut fi and the Fourier- )( kx fi transforms of eq. (2), thus obtaining ),(ˆ ukP , and then by inverse transforming. We shall see that the latter step, with the aid of some approximations, leads to the diffusion or to the fractional diffusion equation, according to the choice of the pdf’s w(t) and λ (x) . In the Laplace-Fourier space the relation (2) becomes (cf. eq. (A3) of the Appendix A): )(ˆ)(1 1)(1),(ˆ kuwu uwukP l--= (4) Note that this expression is rigorous, in the sense that it contains the same information as eqs. (1-3). To proceed further, we must specify the time and jump pdf’s w and λ . In the following, we shall consider the two cases which originate the diffusive and the subdiffusive transport, respectively. Gaussian jumps and exponential waiting times: the diffusive case. It is assumed that the walker performs Gaussian jumps with mean value zero and variance Σ =2 σ . After each jump, it remains immobilized for an exponential time interval with mean value τ . The Gaussian pdf, its mean, variance and F-transform are S- S= x ex pl ; = xE ; s”S= xVar ; )(ˆ s l k ek - = (5) The exponential pdf, its mean, variance and L-transform are t t t etw - = ; t = ][ tE ; ][ t= tVar ; t uuw += (6) Moreover, let W(t) be the cumulative distribution function (cdf) of w(t) , viz. t t etW - -= . Thus, the conditional probability )(1 )(1 tWdttw -- for a rest time t+dt , given that a rest time t has already elapsed, is t dt - . Then, the probability that the rest time ends at t , so that a step of the random walk will indeed take place in ( ) dttt + , , as fi dt , is tt dtdt = -- , independently of t . Then, the assumed semi-Markovian process is actually Markovian. Since both distributions have finite first and second moments, the central limit theorem applies (Metzler, Klafter 2000; Bouchaud, Georges 1990) and we know in advance that the probability P(x,t) that the walker be at x at time t is Gaussian. Substituting the above transforms in eq. (4) yields =-+- - ukPeukPu k s t (7) To the lowest order approximation in k , i.e. for large x , we approximate sl s kek k -»= - (8) and then =+- ukPkukPu st . Inverse L-transforming and then inverse F-transforming yields the standard diffusion equation, viz., ),(),( x txPDt txP ¶¶=¶¶ (9) where tts S== D is the well known diffusion coefficient. The variance of P(x,t) is readily obtained by applying the operator ∫ ¥ ¥- dxx to the diffusion equation dxx txPxtxdtd ),()( ¶¶= 〉〈 ∫ ¥ ¥- ts Integrating twice by parts the integral at rhs yields 2 (recollect that ∫ ¥ ¥- = txdxP ), and finally Dttttx =S== 〉〈 tts (10) In words, the variance of the diffusion process increases linearly with time. The solution to the diffusion equation (9) is ( )
Dtxx eDttxP -- = p (11) It is important to note that eq. (10) is rigorous in spite of the fact that it has been obtained by utilizing the approximation (8) of the F-transform of λ (x) . Indeed, back transforming eq. (7), yields ∫ ¥ ¥- -+-=¶¶ ),'()'('),(),( txPxxdxtxPt txP lt Application of the operator ∫ ¥ ¥- dxx , after some algebra concerning the convolution integral, yields S+ 〉〈 + 〉 - 〈 = 〉〈 xxxdtd t and then eq. (10). Gaussian jumps and power law with α <1 waiting times: the subdiffusive case. In the subdiffusive motion of a walker, e.g. in a porous medium (Montroll, Weiss 1965; Margolin, Berkowitz 2004), the sequence of flights and trappings is characterized by long rests, much longer than those occurring in the diffusive case. Formally, this means that the pdf P(x,t) should have tails fatter than the Gaussian ones pertaining to the diffusive case and, correspondingly, that the pdf w(t) (from which the waiting times are sampled) should be characterized by fat tails, fatter than the exponential ones of the diffusive case. Furthermore, it is required that the distribution of the sum of the waiting times – which represents the walker’s temporal fate – should have tails not substantially variable during the whole walker motion, when new terms are continuously added. This implies that the mean of the sum of the waiting times, constituted by a sequence of several iid random variables sampled from w(t) , should have the same tails, up to a scale factor, as those of w(t) . In other words, this sum should give rise to a limiting stable distribution , similarly to the case of the distributions obeying the central-limit theorem, which roughly states that the limiting distribution of the sum of any iid random variables is Gaussian, provided that the first and second moments of the terms are finite. The upshot is that, to obtain a distribution which is stable and endowed with fatter than Gaussian tails, we should abandon the requirement of the moments finiteness. The solution to this problem is due to the Lévy and Khinchine generalization (Feller 1971) of the central-limit theorem, which takes into account the possibility that the terms of the sum obey distributions having infinite moments. This theorem specifies that a power law of the kind )/()( aa t + (cid:181) ttw , with << a , τ positive, small and t‡ t to prevent divergence, characterized by infinite first moment, satisfies both the above requirements, namely that the mean of the sum of n terms thereby sampled will have fat tails and limiting stable distribution. It is of paramount importance observing that it turns out that the subdiffusive experiments mentioned in the Introduction are actually well accounted for by power laws with suitable a<1 exponents: thus, mathematics matches physics . The Appendix B, after (Metzler, Klafter 2004; Metzler, Klafter 2000), gives some more details on the Lévy-Gnedenko generalized central-limit theorem. In view of the following comparison between Monte Carlo simulations and analytical results, the formalization of the subdiffusive transport is performed by assuming the same Gaussian distribution (5) for the jump lengths and a pdf with power law asymptotic for the waiting times, viz., )(2122)( *1221 twptpttw ==+= ttaa , ),0( t=D˛ tt ( ) ( ) )(11122)( *2112 twptpttw -=-=+= ++ aaaa attaa , ),( +¥=D˛ t tt (12) (Feller 1971, p.166) Let ,...,, XXX be a sequence of mutual independent rv obeying a common distribution R and let nn XXXS +++= ... . The distribution R is said to be stable if S n and c n X have same distribution, viz., XcS ndn = . A theorem states that only the constants a nc n = , with £< a , are possible. The constant α is called the characteristic exponent of R . where aa+= p is the probability of a waiting time less or equal to t and the )( * tw i , i=1,2 , are the pdf’s normalized on the domains t i D . Moreover, let )( * tW i be the cdf of )( * tw i , viz. ∫ == t tdvvwtW )()( t and ∫ -== t tdvvwtW t aa t *2*2 . In this case, the conditional probability )(1 )(1 ** tW dttw ii -- for a rest time t+dt , given that a rest time t has already elapsed, is tdtt -- t for tt D˛ and tdt a- for tt D˛ . Then, the probabilities that a step of the random walk will indeed take place in ( ) dttt + , , as fi dt , are tdtt -t and tdt a , respectively, both depending on t . Then, the underlying process remains semi-Markovian, i.e. Markovian only at the jump instants. The initial linear shape of w(t) has been here introduced to avoid the sudden jump from zero to ta / , which would have occurred in case of a pure power law a tta + = )/)(/()( ttw for ),( ¥˛ t t . Correspondingly, due to the smallness of t , successive jumps essentially with no rest time are possible. The mean number of such jumps is
211 1 a+=- p . The L-transform of the above pdf (12) is )()()( uwuwuw += , where )( )1(122)( t taa t u euuw u - +-+= and -+--+= ∫ ¥ --- t aat ttaaa u xu dxexueuuw )()1()1)(2( 2)( To the first order in u , and then for large t, we approximate ute ut -@ - and )2( a at a -G=@ ∫∫ ¥ --¥ -- dxexdxex xu x so that aa tt )(1)( ucucuw ++= (13) where )2)(1( 2 aa a +-= c and a a a +-G= c (14) where G(.) is the complete Gamma function. Substituting the approximate eqs. (8) and (13) in the rigorous expression (4) we get an approximate expression for the Laplace - Fourier transform of the pdf
P(x,t) that the walker be in x at time t [ ] aaa aa tts t cucuk ucucukP +- +-= -- -
112 11 )()( 1)(),(ˆ (15) Inverse L- and F- transforming, as detailed in the Appendix C, yields the required Fractional Diffusion Equation (FDE) )(1)1( 1),(),(),( xttctxPDtctxPxDtxPtc tt dtatts aaaaaa ----- -G-=¶ ¶-¶ ¶-¶ ¶ (16) Note that the presence of the integral operator a- t D accounts for the non-locality of the above equation, which reflects the semi-Markovian feature of the process. In lack of an explicit solution for (16), eq. (15) is usually approximated (Metzler, Klafter 2000; Metzler, Klafter 2004) by expanding to the lowest order a u the L-transform of the waiting time w(t) . Then, = c . As pointed out in (Margolin, Berkowitz 2004) and as we shall see in the next Section devoted to the use of the Monte Carlo simulation for validating the approximate analytical approach, the term linear in u becomes important at long times only when α is close to . By neglecting the linear term, eq. (16) has an explicit solution in the F- and L- transformed space, namely the Fox function
P(x,t) , which can be expressed in a computable form (Metzler, Klafter 2000; Metzler, Klafter 2004; Margolin, Berkowitz 2004): ( )( ) nn n tKxnntKtxP +-G -= ∑ ¥+= aaaa a (17) As in the case of the standard diffusion, the most important parameter characterizing the transport is the variance )( tx . A fractional differential equation in time for this quantity is readily obtained by applying the operator ∫ ¥ ¥- dxx to eq. (16). In the analytical derivation, the most interesting algebraic step is that by integrating twice by parts dxtxPxx ∫ ¥ ¥- ¶ ¶ ),( simply yields , so that, by eq. (C4’) of the Appendix C, [ ] ( ) a aa G= -- /22 tD t . Therefore, ( ) ( ) txDdtdcctcKtxdtd t aaaaa a ta --- +G= (18) where aaa tts S== K is the generalized diffusion coefficient. Then, by expanding to the lowest order a u the L-transform of the waiting time w(t) as done before, in eq. (18) = c and the second term at rhs disappears. Neglecting this term, the above differential equation becomes no more fractional and its solution is ( ) ( ) aaa a tKctx +G= (19) which differs from those reported in literature, e.g. (Metzler, Klafter 2000; Metzler, Klafter 2004), for the presence of the coefficient a c , which grows from the initial value for α =0 , to 2 for α =2/3 and to 6.6 for α =0.9 . The above expression (19) for the mean squared displacement of the walker at the generic time t confirms the subdiffusive character of the transport when << a - i.e. the spatial spread of the walker around its initial position grows sublinearly with time - and represents one of the most important findings of the fractional diffusion approach. To evaluate the relevance of the approximation leading to (19), we come back to equation (18), which may be analytically transformed in the Laplace space, leading to the simple algebraic expression tt s aaa
211 22 ucucux -= + (20) This expression has the important advantage of lending itself to numerical Laplace inversion. To do this, we have adopted the algorithm (Hollenbeck 1998). The curve thus obtained, which accounts also for the contribution of the linear term, will be used in Section 4 in comparison with the MC variance and the analytical approximated expression (19): we anticipate that (19) will turn out to be rather correct up to values » a , while (20) holds true over the entire range of a . Indeed, it can be observed that in eq. (20), as fi a , the divergence of the coefficient a c is exactly compensated by the one of c : by adopting (19), this compensation is lost and the variance will significantly be underestimated for large values of a .
3. The Monte Carlo approach to the continuous time random walk (CTRW)
The MC approach to the particle transport is based on the simulation of the microscopic behaviour of an ensemble of particles, evolving in phase space according to a master equation such as eq. (1) and macroscopically resulting in the behaviour of the ensemble. In the case of the anomalous transport, a distinct merit of the MC approach with respect to the analytical one of the fractional kinetics is that the simulation is based on the exact w(t) and λ (t) pdf’s governing the fate of each particle, without the necessity of resorting to the various approximations of these pdf’s as detailed in the preceding Section. Therefore, the MC simulations may be utilized to quantitatively evaluate the accuracy of the analytical results. Moreover, the MC approach lends itself to examine a number of variants of the walk not easily affordable analytically. For example, it is not necessary that the pdf’s of the waiting times and of the jump lengths be independent each from the other and therefore also correlated walks may be easily studied. Indeed, if ),( tx y is the bivariate pdf of the jump lengths and trapping times, we sample a jump length X from the marginal ∫ ¥ = ),()( dttxx yl and then a waiting time from the conditional )( ),()|()( XXtXttw lyy == . In addition, the medium in which the walker is moving may be as well anisotropic, with pdf’s containing space and time dependent parameters, e.g. as in Section 5. Finally, the transport is easily simulated in 2D and even 3D spaces. As in almost all cases, the main drawbacks of the MC approach are the requirements of long CPU times and vast memory occupation: however, as time goes by, the continuous computer improvements make these burdens less and less important. In a computer code, before starting the stochastic simulation, the extension of the phase space available for the walker, namely the available space ( ) maxmin , xx and time t max , must be fixed. Then, this rectangular domain is discretized by assigning the respective numbers of subintervals N x and N t : correspondingly a matrix P of order ( N x , N t ) is assigned, whose generic element P(i,j) at the end of the simulation will contain the number of passages of the walker in the elementary ij -th cell. To avoid distortions in the bivariate distribution P and then gross errors in the mean and variance as a function of time, it is important that the space assignment must be done so that essentially no walks will exit from the space interval before t max . The simulation of the fate of each particle is performed by sampling the successive jump lengths ∆ x and waiting times ∆ t from the pdf’s (5) and (12), respectively. Specifically, G rx S=D , where r G is a normal random variable (rv) with zero mean and unit variance, i.e. G r ~ N(0,1). Concerning ∆ t, the cdf of the w i ( ∆ t) given by eq.(12) are D+==D taa ttWptW for t £D t a tat D+-=D-+=D ttWpWptW *2*12 for t ‡D t (21) To sample a ∆ t value, a rv r U is firstly sampled from the uniform distribution, i.e. r U ~ U[0,1) and then this value is equated to W( ∆ t). Solving with respect to ∆ t yields ( ) a at +=D U rt if aa+£ U r ; ( ) a at - - +=D U rt if aa+> U r (22) In Figures 1 and we show the characteristic shapes of MC simulated walks in the x-t space, with different values of a : note that for a=.3 the motion is strongly subdiffusive and indeed the story is characterized by long rest times. The other story, with a=.99 , is instead very close to the one of a fully diffusive motion ( a=1 ) . Figures 1 and 2. Characteristic shapes of walks in phase space x-t, for different values of a . Left: a=.3 ; right: a=.99. Figure 3. Typical shapes of a Fox function (tent-like, a =.5) and a Gaussian, at the same time. At the generic step of the simulation, the walker arrives at a point (x,t) within the cell (i,j) and the new ∆ x and ∆ t are successively sampled. This implies that the walker remains in the i -th spatial interval until t’=t+ ∆ t and then it suddenly performs a spatial jump towards the new point ( ) ',' txxx D+= within the cell ( ) ji mjni ++ , with n i and m j integers. Then, a unit is added to all the elements ( ) ( ) ( )( ) jnijiji i ,1,...,1,, -++ of the matrix P and the walk continues from (x’,t’) . When a large number of stories have been processed in this way and the walker passages through the discrete cells have been registered, the normalized space-time matrix P(i,j) will represent the (discretized) distribution of the particles: in particular, the j-th column of the matrix represents the spatial distribution at time j . In Figure 3 we compare the results for
P(x,t) at a given time in the case of a Fox function (cf. eq. (17)) with a=.5 and a Gaussian function ( a=1 ): observe that the Fox function shows greater persistence of the initial condition and a smaller spread around the origin with respect to the Gaussian.
4. Results and comparisons for a subdiffusive walker
In this Section we present a comparison between the results of the analytical (and numerical) FDE approach to subdiffusion and those of the Monte Carlo simulations both for the spatial distribution ),( txP at given times t and for the time evolution of the variance for different values of the parameter a . Case 1. The parameters values are: a=0.3 , t=10 -4 , s=.5 . The number of simulated stories is and each story will end if the maximum time t max =10 is reached or if the boundaries of the spatial domain are attained during the simulation. These boundaries are allocated through an a priori estimate of the mean spread of the location x given by several )( max2 tx . Figure 4 reports the spatial distribution
P(x,t) at a given time as obtained from the MC simulation and from the analytical Fox function in its representation by series (17).
Figure 5 reports the variances as a function of time, as obtained from the analytical expression (19), from the numerically inverted curve (20) and from the MC. Since in this case a is small, we expect that the effects of the linear term neglected in (19) should be almost negligible and indeed in all cases the analytical curves are in very good agreement with the reference ones of the MC simulation. Fox Function Gaussian Figures 4 and 5. Spatial distribution P(x,t) for a given time(left) and variance evolution in time (right).
Case 2. The parameters values are: a=0.5 , t=10 -2 , s=.5 . The number of simulated stories is and each story will end if the maximum time t max =10 is reached or if the boundaries of the spatial domain are attained during the simulation. In Figures 6 and we show the spatial distribution at a given time and the variance, as above. We remark that also in this case, while a has been increased to .5 , the contribution of the linear term is negligible in the time range considered, so that the analytical results still coincide with the reference MC values. Figures 6 and 7. Spatial distribution P(x,t) for a given time(left) and variance evolution in time (right).
Case 3. The parameters values are: a=0.8 , t=1, s=.5 .
The number of simulated stories is and each story will end if the maximum time t max =10 is reached or if the boundaries of the spatial domain are attained during the simulation. Figure 8 shows the analytical and the MC spatial distributions at a given time.
Figure 9 displays an enlarged representation of one of the (symmetric) tails of the
P(x,t) shown in
Figure 8 . It appears that the analytical distribution slightly but systematically underestimates the reference MC profile. In
Figures 10 and the analytical and MC variances are plotted in log-log and linear scale as a function of time. We remark that in this case ( a close to ) the contribution of the effects of the linear term are no more negligible. However, when the variance given by expression (19) is replaced by the numerical inversion of expression (20), which contains also the linear term, the data perfectly agree with the MC ones, except at small times (cf. Figures 10 and ). We conjecture that these small discrepancies could be due to the so called “ diffusion limit ” approximations (i.e. the expansion of the exact expressions of )(ˆ k l and )( uw for small k and small u , respectively), not appearing in the reference MC approach. We note also that the discrepancy between the FDE analytical model and the reference MC (which does not suffer from approximations) is by far more evident in the variance than in the spatial distribution plots. This happens because of the large contributions of the tails of P(x,t) to the variance. Figures 8 and 9. Spatial distribution P(x,t) for a given time(left) and zoom of its right tail (right).
Figures 10 and 11. Double logarithmic (left) and linear (right) scale plot of the variance vs. time.
Figure 12. Zoom of Figure 11 left, to enlighten the (small) discrepancy between MC and numerically inverted curves
5. Generalization of the anomalous advection-dispersion model for a two-layers medium
We now generalize the stochastic model so far described to the case of the anomalous transport in heterogeneous media. This problem is of importance in different instances, since in reality the transport of particles, e.g. under the influence of an external constant force, gravitational or electric, or the underground transport of radioactive or toxic particles driven by the subsurface water flow normally occurs in heterogeneous media. In the simpler case of an homogeneous medium Analytical, without linear term MC, dots Numerical inversion, with linear term Analytical, without linear term MC (dots) and numerical inversion (solid line) Analytical, without linear term MC (dots) Numerical inversion (solid line) the problem is analytically tackled by means of the fractional Fokker-Planck equation or by means of the formally equal fractional advection equation (van Kampen 1981). Here we focus on the underground transport of particles. We assume that the diffusion of the walker occurs within a constant velocity field v, in a medium constituted by two different homogeneous layers, say L on the left and L on the right, the discontinuity occurring at d xx = . In layer i ( i =1,2) the diffusion contribution to each jump is Gaussian ),0( i S and the time intervals ∆ t between jumps obey to a power law distribution w i with parameter α i (cf. eq. (12)). In both layers the advection contribution to the walker motion is v ∆ t . Differently from the sudden-jump model adopted in the preceding Sections, we now assume that during the time between successive jumps the walker travels with constant speed, variable from jump to jump, given by the ratio between the total jump length and the waiting time. For sake of simplicity, here we assume also the same parameter a for the two zones. Moreover, to a first approximation, in the MC simulation it may be assumed that the optical path related to each Gaussian jump across the discontinuity can be subdivided proportionally to the respective standard deviations. Let us now consider a single jump starting at x in L . The contribution to the walker’s travel due to diffusion is the so-called optical path r G Σ , r G ~ N(0,1), and the jump duration is T , sampled from the power law w . If the whole jump is done within the same L , i.e. if r G Σ + vT < x d - x , the walker location at the generic t ≤ T is (note that r G may well be negative) ( ) TtvTrxtx G +S+= )( (23) Otherwise, if r G Σ + vT > x d - x , the walker crosses x d at time TvTr xxt
G d +S -= and the utilized fraction of r G is S --= vtxxTtr dG At t * the walker is at the entrance of L . The remaining portion of the path within this layer is obtained from the residual fraction of r G , viz., ( ) )( S ---=- ttvxtxTttr dG so that ( ) TttvTrxtx Gd *2 )( -+S+= Ttt ££ * If r G is so negative that <+S vTr G , then d xtx < )( , which would mathematically mean that the walker, as soon as arrived at x d driven by v , instead of proceeding further, is driven back to layer 1. Here, the advection velocity pushes it again towards layer 2 and these back and forth oscillations continue up to time T with the walker essentially at rest in x d . To take into account this possibility, the above equation is modified as follows ( ) [ ] TttvTrxtx Gd *2 ,0max)( -+S+= Ttt ££ * (24) Let us consider the other case, with the starting point x of the jump in L , i.e. d xx > . Now the situation is reversed: if dG xvTrx >+S+ , the whole jump occurs in this layer and the walker location at t ≤ T is ( ) TtvTrxtx G +S+= )( (23’) Otherwise, if x + r G Σ + vT < x d , (which means that r G is so negative that also r G Σ +vT<0 ) the walker flies backwards against the velocity field and crosses x d at time TvTr xxt
G d +S -= or TvTr xxt
G d +S -=
The utilized fraction of r G is S --= vtxxTtr dG or S +-= vtxxTtr dG The remaining of the path within L is obtained as before from the residual fraction of r G , viz., ( ) TttvTrxtx Gd *1 )( -+S+= However, it may happen that >+S vTr G so that, as in the preceding case, the walker is stuck at x d . Application of an analogous remedy finally yields ( ) TttvTrxtx Gd *1 ,0min)( -+S+= Ttt ££ * (24’) The following Figures 13 and show the spatial distribution P(x,t) at two given times as obtained through MC simulation. The walkers move across two distinct zones, separated by a discontinuity at = d x . The parameters of the two regions are the following: a =a =.5, s =5, s =15, t=.01 . The flow speed v was set to , towards the right. As expected, the most relevant consequence of the medium heterogeneity is the significant discontinuity of P(x,t) at d x , due to the ratio of the walkers' speeds during their crossings of the discontinuity. Figures 13 and 14. Space distribution at time t=1 (left) and t=3 (right) for the model presented above.
6. Conclusions
Several experimental investigations concerning transport of particles in given media have shown strong evidences for anomalous non-Fickian diffusion (Kosakowski 2004; Margolin, Berkowitz 2004; Metzler, Klafter 2000; Metzler, Klafter 2004; Scher, Montroll 1975). The interpretation of these results has given rise to various analytical models, which have finally led to the so-called Fractional Diffusion Equation (FDE) approach. Here, the motion of the particles in the surrounding medium is characterized by long sojourn times obeying distributions with infinite mean value, for which the standard central limit theorem must be replaced by the Lévy – Gnedenko generalised central limit theorem (Feller 1971; Metzler, Klafter, 2000). More specifically, the time fate of the particles is described in terms of a power law pdf with exponent a+1 : for the case here considered (i.e. ), the mean squared displacement of the particles evolves as t a , thus resulting in a displacement smaller than the linear one characteristic of Fickian diffusion. Hence the name of subdiffusion given to such phenomena. The FDE analyses are performed in the Laplace and Fourier domains and a closed-form solution in the usual x-t space can not be obtained without resorting to several approximations. Moreover, the effects of the above mentioned approximations can hardly be physically interpreted in the reciprocal domain. Of utmost importance is that a term linear in u (where u is the L-transformed time variable) is not taken into account. In this paper we propound the use of the Monte Carlo (MC) simulation, which does not suffer from any approximation (apart from those due to the finite statistics and to the x-t discretization), as a reference model in the investigation of the anomalous diffusion. The main result here attained is showing that the FDE analytical solutions are far from reality for fia . Nevertheless, the contribution of the above linear term may be taken into account by hybridizing the elegant fully analytical FDE approach with a numerical Laplace inversion. In this case, the resulting curves are quite consistent with the MC results. Another application, which follows from the MC flexibility, concerns the simulation of a more realistic anomalous advection-dispersion model in heterogeneous media, for which the analytical results are awkward to obtain. In lack of an analytical formulation, the MC approach could represent a viable way for numerical investigations of the particles behaviour at the interface of media characterized by different physical parameters. Further exploitations of the MC flexibility, not reported here, concern straightforward generalizations of the anomalous diffusion scheme, such as particles transport with pdfs correlated in time and space and transport in multidimensional media. Appendix A:
Laplace – Fourier transforms of equation (2) The Laplace transform (L-transform) with respect to time )( ut fi of eq. (1) is ∫ ∫ ∫ ¥ ¥ ¥-- +--= )()'()'()','(''),( tut xxxttwtxpdxdtdteuxp dl Inversion of the first and second integral yields )()'()','(')'('),( xxxtxpdxettdtwedtuxp t ttuut ∫ ∫ ∫ ¥ ¥ ¥ ¥---- +--= dl
The inner integral is )( uw , the L-transform of w(t) and may be taken out of the triple integral. Inversion of the remaining two integrals yields )()','(')'(')(),( xetxpdtxxdxuwuxp ut ∫ ∫ ¥ ¥- ¥ - +-= dl The second integral is p(x’,u) , the L-transform of p(x’,t) , so that we get ∫ ¥ ¥- +-= )(),'()'(')(),( xuxpxxdxuwuxp dl This expression is now Fourier transformed (F-transform) with respect to position )( kx fi , viz., )'(' +=+-= ∫∫ ¥ ¥- --¥ ¥- - kukpuwexxdxeuxpdxuwukp xxjkjkx ll Solving for ),(ˆ ukP yields )(ˆ)(1 1),(ˆ kuwukp l-= (A1) We now proceed to transform
P(x,t) as given by eqs. (2) and (3). The L-transform is
IIItwdttxpdtdtetxpdtdteuxP t t ttutut -fi-= ∫ ∫ ∫ ∫∫ ¥ --¥ - )()',(')',('),( (A2) Inversion of the integrals in the first term yields u uxpetxpdtudtetxpdtI utt ut ),()',('1)',(' ==fi ∫∫ ∫ ¥ -¥ ¥ -
Inverting the last two integrals in the second term of eq. (A2) yields ∫ ∫∫ -¥ - fi t ttut txpdttwdtdteII '' )',(')( Inverting the first and the second integrals in the above expression yields ∫∫ ∫ -¥ ¥ - fi ''''
00 '''' )',(')( ttt ut txpdtdtetwdtII
Inverting the last two integrals yields u uxpuwetxpdtetwdtu etxwdttwdtudtetxpdttwdtII utut ttutt ut ),()()',(')(1 )',(')(1)',(')( '' '''' == =fi ∫ ∫ ∫∫∫ ∫∫ ¥ ¥ -- ¥ +-¥¥ ¥ + -¥
Substitution of I and II in eq. (A2) yields ),()(1),( uxpu uwuxP -= Finally, F-transforming, we obtain the required LF-transform of
P(x,t) , viz., )(ˆ)(1 1)(1),(ˆ)(1),(ˆ kuwu uwukpu uwukP l--=-= (A3)
Appendix B: the Lévy – Gnedenko generalization of the central-limit theorem Let ),;( sm xG be the normal (Gauss) pdf with mean value µ and variance σ . We say that the random variable (rv) ξ is normal ( µ , σ ) if { } dxxGdxxx ),;(Pr smx =+££ If n values n xxx ,...,, are independently sampled from the same G , the simplest version of the central limit theorem states that the rv ( ) mmx +-= ∑ = ni in nS obeys to the same G, that is { } dxxGdxxSx n ),;(Pr sm =+££ To simplify the matter, from now on we assume =m and =s . Then the theorem states that, up to a scale factor n , the sum of n normal rv (0,1), namely ∑ = ni i x obeys to the same distribution as the one of the single terms of the sum. Formally, the distribution of the rv ξ is a stable distribution (see Note 1 , § ) with characteristic exponent α =2. Actually, the theorem is much more general but here we only mention the above simple version. When dealing with a power law pdf such as the w(t) given by eqs. (12), the situation is complicated by the lack of a finite mean value. In this case we may resort to the Gnedenko – Lévy generalization (Metzler, Klafter 2000) of the central limit theorem, which states that if the sum S n of n iid (independent-identically-distributed) random variables n xxx ,...,, , apart from an appropriate normalization, converges to some distribution w(t) in the limit ¥fi n , then w(t) is stable. Then, from the definition of stability xx a nS dni in = ∑ = ” in the limit of large n . The following Figures , obtained from a MC simulation with = N stories for the cases of n=10 and a=0.2 ( Figures 15 and , case a ) and a=0.8 ( Figures 17 and , case b ) give an idea of how much relatively small n and t values are suitable to satisfy the theorem. Figures 15 and 16. a). Monte Carlo approximation of the Levy – Gnedenko theorem for finite n and t values. Dots= n S ; solid line= x a n . Left: plot of the distributions. Right: zoom. In this case a=.2 Figures 17 and 18. b). The same as above. Left: plot of the distributions. Right: zoom. In this case a=.8
Note that the convergence is faster in this latter case, where a is closer to 1. Appendix C:
The fractional diffusion equation (FDE)
For convenience, we firstly briefly recollect some definitions and properties of the fractional differ-integration. Detailed books on this subject are e.g. (Oldham, Spanier 1974; Miller, Ross 1993). i)
Fractional integral : The Riemann-Liouville fractional integral of arbitrary order α , i.e. the integral operator a - t D , is defined as ''10 '0 )()()(1)( dttftttfD tt -- ∫ -G= aa a (C1) Note that the application of the operator to a function yields a result which has a power law memory of the previous function values. An important property of the fractional integral which follows from the fact that it is a convolution, is that its L-transform is the product of the transforms of -a t and f(t) , viz., [ ] [ ] )()()( ufuuftLtfDL t aaa --- == (C2) ii) Fractional derivative : Let [ ] += a n , where [ ] a is the integer part of α . Differentiating (C1) n £n times yields anann -- = tt DDdtd . If we change notations by performing the substitution aan fi- so that <£ a , we get the fractional derivative nanna - = tt DdtdD or, explicitly, ''10 '0 )()()( 1)( dttfttdtdtfD tt -- ∫ --G= annna an (C3) In particular, for << a , so that n = ,we get '')1(0 '''0 '0 )()()(1)()()1( 1)( dttfttdttfttdtdtfD ttt +-- ∫∫ --G=--G= aaa aa (C3’) The name of fractional derivative given to the operator a t D follows from the fact that the fractional derivative of order α of an arbitrary power of a variable – and then the fractional derivative of any function which may be expanded in a power series - is performed with the same rules of the usual derivative of integer order, viz., [ ] annaana an n - +-G +G=” ttdtdtD t )1( )1( (C4) In particular, for ν =0 , we get a result at first sight surprising: the fractional derivative of a constant is a function, namely [ ] aaaa a - -G=” tdtdD t )1( 111 (C4’) However, note that, if α is positive and integer, the fractional derivative is zero as it must be, due to the divergence of the gamma function. It is worthwhile mentioning that in the above case the fractional derivative as a function of α oscillates around the zero values taken in correspondence of the integer values of α . Let us now come back to the inversion of eq. (15). Multiplying the lhs by the denominator of the rhs yields [ ] [ ] ucucukPucukPukukPc aaaaaaaa ttts +-=-+ ------ ),(ˆ),(ˆ),(ˆ The inverse L- transform ( ) tu fi is readily obtained from eq. (C2) for the factors in square parentheses and from usual algebra for the others aaaaaaa tatts ctctkPDctkPDktkPc tt + G-=-+ ----- )(),(ˆ),(ˆ),(ˆ
Inverse F-transforming and then differentiating with respect to time yields eq. (16).
References
Bouchaud , J-P., Georges, A. (1990), Phys. Rep. , 12 de Almeida, R. M. C., Baumvol, I. J. R. (2003), Surf. Sci. Rep. , 1 Denisov, S. Klafter, J., Urback, M. (2003), Phys. Rev. Lett. , 193301 Feller W. (1971), An Introduction to Probability Theory and its Applications , Vol. 2, Wiley, New York, Fick, A. (1855), Ann. Phys (Leipzig)
Stochastic Processes in Physics and Chemistry , North Holland, Amsterdam Kosakowski, G. (2004), Journal of Contaminant Hydrology, 72, 23-46 Mainardi, F. Raberto, M. Gorenflo, R. Scalas E. (2000), Physica A
468 Margolin, G. Berkowitz , B. (2004), Physica A, , 46
Metzler, R. Klafter, J. (2000), Phys. Rep., , 1-77 Metzler, R. Klafter, J. (2004), J. Phys. A: Math. Gen. R161-208 Miller, K. S. Ross, B.
An Introduction to the Fractional Calculus and Fractional Differential Equations , (1993), Wiley, New York Montroll, E. W. Weiss, G H (1965), J. Math. Phys. , 178 Oldham, K. B. Spanier J. (1974) The Fractional Calculus , Academic Press, New York Richardson, L. F. (1926), Proc. R. Soc. , 709 Scher, H. Montroll E. W. (1975), Phys. Rev. B , 2455 Slutsky, M. Kardar, M. Mirny L.A. (2004), Phys. Rev. E ,69