Normal and Anomalous Transport across an Interface: Monte Carlo and Analytical Approach
1 Normal and Anomalous Transport across an Interface: Monte Carlo and Analytical Approach
M. Marseguerra , A. Zoia Department of Nuclear Engineering, Polytechnic of Milan, Via Ponzio 34/3, 20133 Milan, Italy
Abstract
We present a Monte Carlo scheme to simulate particles going across an interface separating two layers of a medium characterized by different physical properties, together with an analytical formulation of the same problem, for both normal diffusive and subdiffusive regimes. We relate the Monte Carlo simulation parameters to the coefficients and boundary conditions appearing in the companion analytical equations. Under suitable physical hypotheses on the constraints to be imposed on such parameters, we show that the Monte Carlo simulation results are in good agreement with the corresponding analytical solutions. In particular, we remark that – while in the normal diffusion case the conservation of particle local velocities across the interface leads to a smoothly varying concentration profile – in the subdiffusive case the same condition leads to a neat jump in resident concentration.
1. Introduction
A general approach to the analysis of transport phenomena is based on Continuous-Time Random Walk (CTRW) (Metzler, Klafter, 2000a; Metzler, Klafter, 2004; Scalas et al, 2004; Margolin, Berkowitz, 2004; Weiss, 1994; Bouchaud, Georges, 1990; Klafter at al., 1996), in which the travel of a particle (a walker) in a medium is modelled as a series of jumps of random lengths, separated by random waiting times. The theory of CTRW with algebraically decaying probability distribution functions (pdf’s) has been originally introduced in Physics in a series of seminal papers by Weiss, Scher, Montroll and co-workers (Weiss, 1994; Montroll, Weiss, 1965; Scher, Montroll, 1975) in the late 1960s to explain evidences of anomalous diffusion occurring in the drift-diffusion processes in amorphous semiconductors. The diffusion is called anomalous if the mean squared displacement (MSD) is not linearly proportional to time t as in the standard Fickian case, but to powers of t larger (superdiffusion) or smaller (subdiffusion) than unity. More recently, anomalous diffusion has turned out to be quite ubiquitous in almost every field of science (see e.g. (Metzler, Klafter, 2000a; Metzler, Klafter, 2004) for a detailed review) and the CTRW model has been applied with success to interpret the experimental results and to make predictions on the evolution of the examined systems (Margolin, Berkowitz, 2004; Weiss, 1994; Bouchaud, Georges, 1990; Klafter at al., 1996; Scher, Montroll, 1975). Such applications concern, among others, e.g. the behaviour of chaotic Hamiltonian systems – with application to the transport of charged particles in turbulent plasma – (Shlesinger et al., 1993; Zaslavsky, 2002), the evolution of financial markets (Mainardi et al., 2000), the dynamics of ad-atoms on the surface of a solid (Vega et al., 2002) or the underground transport of contaminant (toxic and/or radioactive) particles in presence of rock fractures and porosity (Berkowitz, Scher, 2001; Scher et al. 1991; Berkowitz, Scher, 1998; Berkowitz et al., 2001; Cortis, Berkowitz, 2004; Berkowitz, Scher, 1997; Levy, Berkowitz, 2003; Margolin, Berkowitz, 2002). Much attention has been paid to this last topic, as many laboratory-scale experiments as well as direct field measurements have evidenced subdiffusive and superdiffusive behaviours of migrating particles. It can be shown that exact transport formulation of CTRW, under the so called “diffusion limit”, i.e. considering the asymptotic behaviour of the walkers far from the origin both in space and in time, can be successfully approximated by the so-called Fractional Diffusion Equation (FDE). FDE approach allows to obtain closed-form analytical solutions (Scalas et al., 2004; Metzler, Klafter 2000a; Metzler, Klafter, 2004) and to deal with boundary conditions problems (Metzler, Klafter 2000b; Metzler, Klafter 2000a; Metzler, Klafter, 2004), by making use of non-integer order (pseudo-) derivatives and of their definition in the Laplace and Fourier transformed spaces. Moreover, the CTRW model with exponentially decaying pdf’s naturally takes into account Fickian Corresponding Author. Tel: +39 02 2399 6355 . Fax: +39 02 2399 6309 Email address: [email protected] (Marzio Marseguerra) Email address: [email protected] (Andrea Zoia) diffusion and under the previous asymptotic expansion reduces to the normal diffusion equation, which may be considered as a particular case of the more general FDE. The present paper focuses on the Monte Carlo approach to the diffusive as well as subdiffusive transport of particles in a medium constituted by two regions characterized by different physical properties. A large number of “walkers” are followed along their simulated travels through the medium: their ensemble behaviour yields the flux f (x,t) , the net current J(x,t) and the resident concentration
P(x,t) . The core of the present analysis is the investigation of the effects of the interface, which introduces a macroscopic heterogeneity in the traversed medium. The discontinuity has the effect of modifying the distributions of the jump lengths and of the flight times across the discontinuity itself. Moreover, while in the standard CTRW model the walker is assumed to be trapped for the whole waiting time at the starting point of each jump and then suddenly transferred to the new sojourn location at infinite speed, here we will assume a uniform linear motion between successive interactions. In a previous paper we have presented a Monte Carlo simulation scheme by which this general problem may be tackled, both in the normal diffusion and subdiffusion case (Marseguerra, Zoia, 2006). In particular, we showed that the resident concentration reveals neat jumps at the interface when the simulation parameters are arbitrarily chosen. In parallel with the Monte Carlo approach, in this paper we present an analytical scheme to model diffusion in two-layered media and compare the results both for the normal diffusive and subdiffusive regimes. In particular, we show how the Monte Carlo simulation parameters are related to the coefficients and boundary conditions of the companion equations. The paper is structured as follows: in Section 2 we introduce the statement of the problem, i.e. simulating particles crossing an interface. In Section 3 we summarize the Monte Carlo simulation scheme we will adopt to this aim. Then, in Section 4 we briefly recall how flux, current and resident concentration may be collected during Monte Carlo simulation. Section 5 is devoted to the development of the analytical formulation for the problem of normal diffusion and subdiffusion across an interface in a two-layered 1D medium. Results are presented and discussed in Section 6, where Monte Carlo simulation findings are compared to the curves obtained from the analytical formulation for different choices of the parameters, motivated by physical considerations. Conclusions are finally drawn in Section 7.
2. Statement of the problem: walking across a discontinuity
Both normal and anomalous transport can be successfully described in the framework of the CTRW approach, in which a walker diffuses by performing a sequence of semi-Markovian jumps in space and time. For sake of simplicity, the transport is here thought to occur on a one-dimensional support. We assume that the transport takes place in a medium constituted by two different layers, say L on the left and L on the right, the discontinuity occurring at x=x d . Each layer is characterized by different physical properties, i.e. by different parameters of the jump lengths and waiting times distributions. In principle, the pdf’s for waiting times and space lengths of the walker could depend on each other: however, throughout the paper we will assume that the two pdf’s are decoupled. The semi-Markovianity assumption amounts to saying that, along the walker’s trajectory, each jump starts without memory of the way in which the jumping point has been reached. The process is regenerated at each jump start along the trajectory and it is then ruled by space and time probability density functions (pdf’s) which – besides depending on the features of the surrounding medium – in general depend also on the time and space elapsed from the jump origin. As a special case, normal diffusion is recovered if the process is Markovian. In each layer i (i=1,2) the walker undergoes a succession of jumps, each performed in a time interval ∆ t drawn from a pdf w i ( ∆ t) : the jump lengths are given by a diffusive stochastic contribution ∆ x independently drawn from a pdf λ i ( ∆ x) . When the particles trajectories start in a layer and end in the other one, thus passing through the discontinuity, the situation is similar to that of the time paradox described by Feller (Feller, 1971), who put in evidence that, even in the trivial case of a succession of exponential time intervals with mean value τ , the distribution of the intervals which went across a fixed time point was “different”, i.e. no more exponential and with a mean value τ . Similarly, in our transport problem, it turns out that the space and time distributions of the jumps across the heterogeneity are different and may even cause an important discontinuity at x d in the shape of the resident concentration P(x,t) . Moreover, the present analysis stands out from the usual one in the CTRW context (Metzler, Klafter, 2000a; Metzler, Klafter, 2004; Bouchaud, Georges, 1990; Scalas et al., 2004; Klafter at al., 1996; Weiss, 1994), in which it is assumed that the walker, after a spatial jump, rests in the reached location for the time interval ∆ t and then performs another instantaneous jump (at infinite speed) to the successive location. Instead, we will assume that at each jump the walker moves linearly at constant local velocity – given by the ratio of the jump length to the flight time – between successive interactions: under this assumption, we will have a spectrum of velocities as a result of the stochastic variability of both the jump lengths and the flight times.
3. The Monte Carlo scheme for walker’s jumps across a discontinuity
This section summarizes and extends the Monte Carlo scheme proposed in (Marseguerra, Zoia, 2006), which describes how the walkers cross an interface between different media, in presence of an advective field V and of a bias m . The presence of a velocity V , which leads to the so-called Galilei invariant advection scheme (Metzler, Klafter, 2000a), increases the length of each jump of duration t by a quantity Vt : for the jumps across the interface, this space is partitioned in proportion to the times spent in the first and the second layer, respectively. On the other hand, the presence of a bias, which leads to the so-called Galilei variant advection scheme (Metzler, Klafter, 2000a), increases each jump length by a fixed contribution m : for the jumps crossing the interface, m is partitioned in proportion to the spaces traversed by the walker in the two layers. A. Normal diffusion The diffusion contribution to each jump length in layer L i is assumed to be Gaussian with pdf λ ( ∆ x|0, σ i ) and cdf Λ ( ∆ x|0, σ i ) , respectively (to simplify the notations we shall write x and t instead of ∆ x and ∆ t ), viz. i xii ex s spsl - = and dzzx x ii ),0|(),0|( ∫ ¥- =L sls . (1) As for the distribution of the flight times, we adopt exponential distributions with means τ i . The pdf’s and cdf’s are w(t| τ i ) and W(t| τ i ) , respectively, viz. i tii etw t tt - = and i ti etW t t - -= (2) Let us now consider a walker’s jump starting at a generic point Lx ˛ at time t , as determined by sampling two random numbers x R and t R , both uniform in [0,1). If the jump is totally performed in L , the duration and the length of the jump are )1ln()|( tt RRWt --== - tt and m++= ** VtxX diff (3) respectively, where ),0|( s xdiff Rx - L= (4) is the diffusion contribution to the jump length. Suppose now that, along its trajectory, the walker arrives at a point x at a known distance *01 XxxX d <-= from the discontinuity and that the jump starting in x crosses the discontinuity after an (unknown) time t VtxX +=- rrr , so that ( ) VtxVtx VtxVtx iiiii +++ +=+-= mrrm . Then: Recollect that ),0|(),|( smsm bxbx -L=L and ),0|(),|( smsm RaaR -- L+=L >" ba . ( ) iiiiii VtxAVtxX +=++= m , (16) where VtxVtxA ++++= m . Note that the two portions of each jump across the interface are travelled by each walker with local velocities += +== VtxAVtxAtXu diff *11111 and += +== VtxAVtxAtXu diff *211222222 ttss (17) Thus, each walker does not change its local velocity going across the discontinuity provided that tsts = . The ratio σ i / τ i has the dimensions of a velocity and it will be here called parametric velocity . Thus, eq. (17) states that in the jumps across the discontinuity the local velocity does not change provided the parametric velocity does not change either. We shall see below that all the present Monte Carlo results support the conjecture that the persistence of the velocity of the walkers across the discontinuity represents a key feature for characterizing the macroscopic quantities – resident concentration, flux, current – near the discontinuity. In the jumps across the discontinuity, the presence of the bias burdens the computations since µ must be divided in proportions to X and X , this latter being unknown. The problem may be tackled by firstly computing µ by solving the 2 nd degree algebraic equation which results from the substitution in eq. (15) of t taken from eq. (5) and X = (1- µ )X / µ . Once µ has been determined – or in absence of bias – one computes in succession t from eq. (5), x from eq. (6), t from eq. (9) and x from eq. (14’). Then, the discontinuity is crossed at (t +t , x +X ) and the starting point in L for the next jump is (t +t +t , x +X +X ) . We now consider a jump starting at x in L . As before, two uniform random numbers x R and t R are sampled and the jump duration and the diffusive contribution to the length are )|( t t RWt - = and ),0|( s xdiff Rx - L= , respectively. The total jump length would be m++ * Vtx diff if the new location were still in L . The analysis is similar to the one performed for Lx ˛ , with reversed σ , τ and σ , τ . B. Anomalous diffusion The space travelled in each jump is ruled by the same Gaussian distribution as in the previous case of normal diffusion. Therefore, the walker jumps differ only for the time behaviour. In particular, expressions (6) are still valid, provided that t * is now sampled from the new jump time pdf. As for the flight times, we adopt a power law distribution with pdf and cdf w(t| τ i , α i ) and W(t| τ i , α i ) , respectively, viz. tptw iiii tat = and ),|( = iiii tptW tat for i t t££ (18’) ii tptw iiiii aa taat + -= )1(),|( and ii tpdzzwpptW iiit iiiii at tatat --=-+= ∫ )1(1),|()1(),|( for ¥<£ t i t , (18’’) where )2/(),|( iiiiii Wp aaatt +== is the probability of a waiting time smaller or equal to τ i . If the jump is totally performed in L , its duration is ),|( at t RWt - = (19) and x * and x diff are given by (3) and (4). Let us consider a jump across the discontinuity. The time t spent in L and the diffusion contribution x are given by (5) and (6). The time t spent in L is given by the difference between ),|( at tp RWt - = , the travel time if it were totally effectuated in L , and P t , solution of ),|(),|( atat tWtW P = , the time corresponding to the amount of t R not utilized in L . Then, ( ) ( ) ),|,|(,| atatat tWWRWt t -- -= (20) and the total time duration is T=t +t . Note that, in absence of discontinuity, i.e. if ttt == and aaa == , we have ( ) ),|,|( ttWW = - atat and the expression for T becomes the standard one for a jump totally evolving in an homogeneous medium, i.e. ),|( at t RWT - = . Moreover, if sss == , we have ( )( ) **1 /||/ ttxttx ddiffddiff =LL - ss and the expression for X becomes the standard one for a jump totally evolving in an homogeneous medium, i.e. )|( s x RX - L= . Because of the above choice of the time distribution, we must further distinguish the following cases: i) ),min( aa£ t R : then the time cdf’s for the two layers are both parabolic and t Rt at += , ( ) tttRt t -++=++-+= aatttaatat (21’) The walker’s velocity is +++== VtxAtXu diff *212112222 aattss . ii) ),max( aa‡ t R : then the time cdf’s for the two layers are both power law and ( ) /111* a t t Rt -= , ( ) ( ) ( ) -= --= - aaaaaaaaa t tttt tttRt t (21’’) The walker’s velocity is +--== Vtt tttxAtXu diff /1/* 1**2/112222 )()( aaaaaa ttss . iii) ),max(),min( aaaa ££ t R : then one time cdf is parabolic and the other power law. The expressions for the times spent in the two layers must be taken from case i) and the other from case ii). From the above expressions it appears that, in case of a jump across the discontinuity, the two velocities are equal not only provided that tsts = as in the normal diffusion case, but also that aa = . The computational steps are similar to those described for the normal diffusion. However, in this case the presence of the bias adds more complications, since in order to compute µ we are lead to a 2 nd order equation only for the case i), while in the other case we arrive at a transcendental equation in )/(11 aa m + , to be solved numerically. As a final remark, we would like to point out a possible source of errors in a Monte Carlo game, resulting from the presence of V or µ . Let us consider e.g. a jump starting in L , composed by a backward x diff (originated by an < x R in (4)) and by a large Vt or µ which might overcome x diff and lead the walker inside L . When the walker arrives at the discontinuity, it must adopt the value of the backward diffusion contribution which results from the L rules: this new value, still negative, may be large enough to exceed the forward component due to V or µ and drive the walker backwards in L . In such cases the walkers can only remain stuck at the discontinuity for all the remaining jump times. A spurious spike will then appear at the discontinuity and the remedy consists in eliminating the whole trajectory. We now consider a jump starting at x in L . As before, two uniform random numbers x R and t R are sampled and the jump duration and the diffusive contribution to the length are ),|( at t RWt - = and )|( s xdiff Rx - L= , respectively. The total jump length would be m ++ * Vtx diff if the new location were still in L . The analysis is similar to that for Lx ˛ , with reversed a , σ , τ and a , , σ , τ . 4. Monte Carlo estimates of flux, current and concentration Resident concentration is physically defined as the mass of tracer particles per unit volume of flowing substance contained in an elementary volume of the medium at a given time. Flux is defined as the total mass of tracer particles passing through a given cross section during an elementary time interval (regardless of their direction), while current is the net mass of tracer particles passing through a given cross section during an elementary time interval in a given direction (Kreft, Zuber, 1978; Cortis, Berkowitz, 2004). With reference to contaminant tracer particles transport, it should be remarked that current is the most commonly measured variable, e.g. in experiments determining breakthrough curves at the end of a test column where tracer particles are free to diffuse (moreover, at the end of the column flux and current coincide, as there is no back-flow from outside). In a Monte Carlo simulation, the phase space { } xt r , is discretized in cells where weights are accumulated. When a walker moves across a sequence of cells, the contribution of its trajectory to the mean flux in a given cell at a given (discrete) time is obtained by accumulating the length of the travel spanned within the cell at that time. Current is collected in a similar way, by simply giving to each spanned length contribution a positive sign if the walker is crossing the cell in the positive x-axis direction and a negative sign otherwise. Analogously, the contribution of the trajectory to the mean resident concentration in a given cell at a given (discrete) position is obtained by accumulating the time length spent within the cell at that position (Briesmeister, 2000). In our case of one-dimensional transport, the phase space is discretized with a rectangular grid composed of identical cells of area dtdx . Thus, the contribution to the flux in each cell is constant and equal to dx and similarly the contribution to the concentration is constant and equal to dt . The contribution to the current is either +dx or –dx , depending on the particle direction. 5. Analytical formulation of the problem A general approach to model normal diffusion in a two-layered medium characterized by an interface located at x d >0 may be derived in the following way: we assume mass conservation in the form =¶ ¶+¶ ¶ txJxtxPt , (23) where P(x,t) is the resident concentration and J(x,t) is the net current exiting the reference volume, positive when forwardly directed. Here we will restrict our analysis to the case of dispersion alone, in absence of an advective field V or a bias m . We will first consider the normal diffusive case, for which the following general equation may be adopted: ),()(2 )(),( txPxvxxxtxPt ¶ ¶¶ ¶=¶ ¶ s (24) where the parameters s and t have the same meaning as in Section 3 and )( )()( xxxv ts= is the parametric velocity at x . A feature of expression (24) is that it generalizes Fick’s first law by splitting the classical diffusion coefficient D(x) in two factors, namely x s outside the gradient and )( xv inside. This decomposition has been adopted as a particular case of the general diffusion operator presented in (Lejay, Martinez, 2006; Etoré, 2006): as a further motivation, we observe that this choice, together with the boundary conditions (28) and (29) below, gives rise to analytical solutions which are always in good agreement with the Monte Carlo simulation results for any choice of the parameters. A comparison between the above equation and (23) shows that the corresponding form of the current J(x,t) is: ),()(2 )(),( txPxvxxtxJ ¶ ¶-= s (25) As known, this analytical formulation is strictly valid far both in time and in space from the source and corresponds to the asymptotic limit of the exact underlying CTRW transport model. The process is assumed to start at the spatial origin at t=0 , namely, )(),0(lim xtP t d= fi , and, if the 1D domain has an infinite extension, the natural boundary conditions at infinity for this Cauchy problem are =–¥ tP . For sake of simplicity, we further assume that s (x) and v(x) , the parameters which characterize the medium, are piecewise constant in each layer. If the interface occurs at x d , we thus have )(,)(,)( vxvxx === ttss for d xx < and )(,)(,)( vxvxx === ttss for d xx > , respectively, so that the problem may be decomposed as follows: ),(2),( txPxxvtxPt ¶ ¶¶ ¶=¶ ¶ s for d xx < (26) ),(2),( txPxxvtxPt ¶ ¶¶ ¶=¶ ¶ s for d xx > , (27) with the following conditions at the interface x d : +- ¶ ¶-=¶ ¶- dd xx txPxvtxPxv ),(2),(2 ss (28) and +- = dd xx txPvtxPv ),(),( . (29) Condition (28) stems from mass conservation, while boundary condition (29) is chosen so that the argument of the inner spatial derivative appearing in (24) be continuous at the interface. Moreover, this latter condition clarifies the relationship existing between the local and parametric velocities and concentration: when vv = , the local velocities of the walkers will be preserved at the interface (see (17)) and the resident concentration will be continuous. Indeed, as a particular case, if we require the Fick’s first law to hold true on the whole domain, in the form ),(2 )(),()(),( txPxvxtxPxxDtxJ ¶ ¶-=¶ ¶-= s , (30) then it must be vv = at the interface and the resident concentrations will be equal at x d . Thus, if the parametric velocity v(x) is conserved crossing the boundary, one recovers the celebrated solution originally proposed by Carslaw and Jaeger and reconsidered by Uffink (Carslaw, Jaeger, 1959; Uffink, 1985; Uffink, 1990; LaBolle et al. 1996; LaBolle et al. 1998; LaBolle et al. 2000), formulated in terms of the method of images. Defining iii vD s= , i=1,2 , and supposing that DD > , we obtain: dd xxxx txPtxPtxP >< += ),(),(),( , (31) where tDxxtDxxx dd etDRetDtxP 44 1),( ---< += pp for d xx < (32) tD xxxx dd etDRtxP -+-< -= b p for d xx > . (33) In order to satisfy the boundary condition +- = dd xx txPtxP ),(),( , (34) the parameters R , which plays the role of a reflection coefficient, and b must have the following forms (Carslaw, Jaeger, 1959; Uffink, 1985; Uffink, 1990; LaBolle et al. 1996; LaBolle et al. 1998; LaBolle et al. 2000): 21 21 DD DDR +-= and DD =b . (35) This situation, i.e. continuity of resident concentration, is the most commonly found to occur in physical experiments of diffusion through a neat interface between two different media (Carslaw, Jaeger, 1959; Uffink, 1985; Uffink, 1990; LaBolle et al. 1996; LaBolle et al. 1998; LaBolle et al. 2000; Hoteit et al., 2000; Parlange et al., 1984). However, other evidences point out that boundary condition (29) with vv „ may take place, too (Ovaskainen, Cornell, 2003; Van Genuchten et al., 1984; Van Genuchten, Parker, 1984; Schwartz et al., 1999; Novakowski, 1992a; Novakowski, 1992b; Leij et al., 1991; Barrat, Chiaruttini, 2003): as a consequence, the resident concentration will present a neat jump at the interface. The general solution to this problem is here obtained by suitably generalizing Carslaw and Jaeger’s solution, requiring that at the interface +- = dd xx txPvtxPv ),(),( . In order to satisfy this boundary condition, it must be 21 21 tt tt +-” R in (32) and in (33). We remark that condition (28) does not depend on R if we assume P(x,t) in the form (31). The case of anomalous diffusion in a two-layered infinite medium is here derived as a generalization of the normal diffusion. Assuming as usual mass conservation (23), we propose to write the current as follows: ),()(2 )(),( )(1,0 txPxvxxtxJ xt aa s - ¶¶ ¶-= , (36) where a- ¶ t is the Riemann-Liouville pseudo-differential operator and )( )( )()( x x xxv aa t s= . Expression (36) could be regarded as a generalized fractional Fick’s law (Metzler, Klafter, 2000a). As we shall see below, this assumption leads to an analytical formulation which is in good agreement with the Monte Carlo results. Analogously as in the previous case, we assume for sake of simplicity that the parameters characterizing the physical properties are piecewise constant in the two layers, so that )(,)(,)(,)( aa aattss vxvxxx ==== for d xx < and )(,)(,)(,)( aa aattss vxvxxx ==== for d xx > . Substituting (36) in (23) yields ),()(2 )(),( )(1,0 txPxvxxxtxPt xt aa s - ¶¶ ¶¶ ¶=¶ ¶ , (37) which, when the medium is homogeneous, takes the form of the well-known Fractional Diffusion Equation (FDE) with generalized diffusion coefficient a ts = a D (Kilbas et al., 2006; Podlubny, 1999; Oldham, Spanier, 1974; Miller, Ross, 1993; Mainardi et al., 2005; Gorenflo et al., 2004; Gorenflo et al., 2002; Metzler, Klafter, 2000a; Metzler, Klafter, 2004). Moreover, we remark that the inner x ¶ ¶ in (37) operates on a quantity which has the dimensions of the product of a speed times a concentration, similarly as in the case of Fickian diffusion. Exactly as in the case of normal diffusion, the FDE equation is strictly valid only far in time and in space from the source and therefore it suitably describes the asymptotic behaviour of the walkers with respect to the exact CTRW transport formulation (Scalas et al., 2004). Again, the problem may be decomposed on two separate domains, where ),(2),( txPxxvtxPt t ¶ ¶¶ ¶¶=¶ ¶ - aa s for d xx < (38) ),(2),( txPxxvtxPt t ¶ ¶¶ ¶¶=¶ ¶ - aa s for d xx > , (39) with the following conditions at the interface x d : By definition, { } ∫ -- -¶ ¶G”¶ tt dttt txfttxf ')'( )',()(1),( aa a for any sufficiently well-behaved function f (See e.g. Kilbas et al., 2006; Podlubny, 1999). A similar approach may be adopted to describe also superdiffusive transport and in this case it would lead to a non-local fractional Fick’s law for the walkers (Paradisi et al., 2001). +- ¶ ¶¶-=¶ ¶¶- -- dd xtxt txPxvtxPxv ),(2),(2 aaaa ss (40) and +- -- ¶=¶ dd xtxt txPvtxPv ),(),( aaaa . (41) Boundary conditions problems for Fractional Diffusion Equations have been carefully and thoroughly examined by several authors (see e.g. Metzler, Klafter, 2000b; Metzler, Klafter, 2000a; Krepysheva et al. 2006a; Krepysheva et al., 2006b) both for absorbing and reflective boundaries. Yet, to the authors’ best knowledge, the analytical approach to subdiffusion through a two-layered medium and its connection with the underlying Monte Carlo microscopic dynamics are new and unexplored subjects. Boundary conditions (40) and (41) generalize expressions (28) and (29). The fundamental solution of the Cauchy problem for (37) when the medium is homogeneous and we require =–¥ tP is given by the Fox’s H function, so that P(x,t)=H(x,t) (Metzler, Klafter, 2000a; Metzler, Klafter, 2004). Thus, we propose to generalize solution (31) for the two-layered medium by applying the method of images to the Fox’s function. Then, defining ,, iii vD aa s= , i=1,2 , and requiring aa DD > , the Laplace transform of P(x,t) reads ( ) ( ) ( ) ( )( ) ddd xxdxxdxx uxuxHuRuxxHuRuxHuxP ><< -+-+-+= ,1)()(1,2)(,),( aaa b . (42) Here we have denoted the Laplace transform of P by its argument u and the index of the H functions indicates whether the parameters of the right or left layer are to be chosen. In order to satisfy boundary condition (41), the coefficients )( u a b and )( uR a (which generalize b and R appearing in (32) and (33), respectively) are: )( aa aaa b uD uDu = and 21 21 )()( )()()( 21 21 aa aaa tt tt uu uuuR +-= . (43) Note that the inverse transform of (42) is given by the convolution of the H i functions with )( uR a . As in the normal diffusion case, condition (40) does not depend on )( uR a , provided that P(x,u) has the functional form (42). As a particular case, when aaa == , a b and a R do not depend on u and may be identically rewritten as aaa b DD = and aa aaa tt tt 21 21 +-= R , (44) which have a similar structure as the corresponding b and R of the normal diffusion case. Moreover, when fia the expressions for normal diffusion are recovered. As a general remark, the analytical formulation of the transport across a discontinuity, both for anomalous and normal diffusion, is expected to be valid only in the so-called “diffusion limit” approximation, i.e. t>> t and s>> x for each layer separately. This is a consequence of the fact that both the normal diffusion equation and the FDE are asymptotic approximations of the exact CTRW transport model, as mentioned before (Metzler, Klafter, 2000a; Metzler, Klafter, 2004; Scalas et al., 2004; Margolin, Berkowitz, 2004; Marseguerra, Zoia, 2006). 6. Results By taking inspiration from the usual Monte Carlo and analytical approaches to transport phenomena, throughout this Section we assume that the Monte Carlo simulation parameters adopted in the scheme of Section 3 coincide with the coefficients denoted by the same symbols appearing in the analytical formulation of Section 5, for both normal and anomalous diffusion. Then, we proceed to compare the Monte Carlo simulation results for resident concentration and net current with the corresponding analytical curves, under different physical assumptions, leading to different choices of the parameters. As said before, the boundary conditions adopted in (24) and in (37) are those which can best take into account the interface conditions implicitly assumed in the Monte Carlo formulation. In all the following Figures, the analytical curves are plotted as solid lines, whereas Monte Carlo results as dots. The larger statistical fluctuations of the net current J(x,t) with respect to those of the concentration P(x,t) are due to the fact that J is obtained as a difference between forward and backward contributions. We begin with the normal diffusion case. Concerning the Monte Carlo simulation, we assume at first that walkers maintain their local velocities while crossing the boundary. This assumption is motivated by the fact that in Nature diffusing particles are often found to conserve their speed while crossing an interface separating two media with different physical properties, e.g. in the case of neutrons (Weinberg, Wigner, 1958). By making use of expression (17), this condition imposes the equality of the parametric velocities, namely tsts = . Within the analytical approach, this assumption on the microscopic behaviour of the single walkers leads to a continuity of macroscopic resident concentration at the interface, namely +- = dd xx txPtxP ),(),( (see eq. (29)). Figures 1a and 1b report resident concentration P(x,t) and current J(x,t) as obtained from Monte Carlo simulation and from analytical formulation as given by (31) and (25). The parameters were set as follows: x d =5, t =0.1, t =0.01, s =0.707, s =0.0707, so that D =2.5 and D =0.25. It appears that the two approaches are in good agreement. Any other choice of the parameters implying different parametric velocities (and different local velocities, as well) on the two sides of the interface will give rise to neat bumps in resident concentration profiles: this observation is substantiated by some physical evidences pointing out that macroscopic jumps in resident concentration are possible, even if not so common, e.g. in contaminant transport in test columns (Van Genuchten et al., 1984; Van Genuchten, Parker, 1984; Schwartz et al., 1999; Novakowski, 1992a; Novakowski, 1992b; Leij et al., 1991) or even in temperature profiles (Barrat, Chiaruttini, 2003). An example of resident concentration and current as obtained from Monte Carlo simulation with parameters x d =5 t =0.1, t =0.03, s =0.707, s = s /7 (so that D =2.5, D =0.17 and vv „ ) is shown in Figures 2a and 2b together with the corresponding analytical curves as given by (31) and (25). A good agreement between analytical and Monte Carlo results is evident also in this case. −20 −15 −10 −5 0 5 1000.511.522.533.54 x 10 −3 Figure 1a. Normal diffusion: resident concentration profile P(x,t) at t=6 when local parametric velocities are equal. −15 −10 −5 0 5 10−0.05−0.04−0.03−0.02−0.0100.010.020.030.040.05 Figure 1b. Normal diffusion: net current J(x,t) at t=1.5 when parametric velocities are equal. −20 −15 −10 −5 0 5 1001234567 x 10 −3 Figure 2a. Normal diffusion: resident concentration profile P(x,t) at t=6 when parametric velocities are different. −15 −10 −5 0 5 10−0.04−0.03−0.02−0.0100.010.020.030.04 Figure 2b. Normal diffusion: net current J(x,t) at t=2 when parametric velocities are different. As for the anomalous diffusion case, the same physical assumption of particles preserving their local velocities across the boundary imposes two separate constraints on the Monte Carlo parameters, namely the equality of parametric velocities tsts = and also aa = (cf. (22)). Differently from the normal diffusion case, these requirements do not automatically imply a smoothly varying resident concentration profile: Figure 3a, where x d =5, t =10 -4 , t =10 -5 , s =0.7, s = s /10, a = a =0.5 (so that D a, =24.5 and D a, =0.775), reveals a neat bump in resident concentration at the interface. Figure 3b shows the net current for the same parameters: the analytical formulations (42) and (36) compare well with the Monte Carlo results, but for slight differences at the interface. With reference to the analytical formulation (42), we remark that imposing resident concentration continuity at the interface requires that the generalized parametric velocities )( )( )()( x x xxv aa t s= should be equal on both sides and aa = . If the same choice is adopted within Monte Carlo simulation, resident concentration profile is shown to smoothly vary across the boundary, while local velocities are not preserved. Figures 4a and 4b compare the Monte Carlo results with the corresponding analytical curves (42) and (36) when the parameters are x d =5, t =10 -3 , t =1.768*10 -4 , s =0.4, s = s /4, a = a =0.8 (so that v a, =v a, =100.5, D a, =20.1 and D a, =5). −30 −25 −20 −15 −10 −5 0 5 10 1500.0020.0040.0060.0080.010.0120.014 Figure 3a. Anomalous diffusion: resident concentration profile P(x,t) at t=14 with equal parametric velocities. −35 −30 −25 −20 −15 −10 −5 0 5 10 15−8−6−4−2024681012 x 10 −3 Figure 3b. Anomalous diffusion: net current J(x,t) at t=1.2 with equal parametric velocities. −15 −10 −5 0 5 100123456 x 10 −3 Figure 4a. Anomalous diffusion: resident concentration profile P(x,t) at t=0.8 with equal generalized velocities. −15 −10 −5 0 5 10−0.04−0.03−0.02−0.0100.010.020.030.04 Figure 4b. Anomalous diffusion: net current J(x,t) at t=0.38 with equal generalized velocities. We remark that the previous comparisons between the Monte Carlo simulations and the analytical formulation require that the asymptotic regime, i.e. the diffusion limit, has been attained. This observation is particularly pertinent in correspondence of a bump at the interface: we conjecture that the interface represents a sort of “virtual source” for the walkers and that the attainment of the asymptotic regime must be faster in the Fickian than in the anomalous case. This conclusion is supported by the results shown in Figures 5 and 6, where slight deviations between Monte Carlo and analytical curves at the interface are present only in the anomalous diffusion case (Figure 6) and are completely negligible for normal diffusion (Figure 5). This fact is thought to be related to the “memory effect” which characterizes non-Markovian process such as subdiffusion. −3 Figure 5. Normal diffusion: zoom of the bump appearing in Figure 2a at two different times, t=3 (o) and t=6 (*). Figure 6. Anomalous diffusion: zoom of the bump appearing in Figure 3a at two different times, t=2 (o) and t=14(*). Finally, note that in the analysis of field or laboratory data the choice of the kind of underlying kinetics (normal or anomalous diffusion regimes) and consequently of the appropriate parameters is usually performed by fitting procedures. However, with reference e.g. to transport of contaminant tracer particles, in most cases resident concentration P(x,t) is not the directly measurable variable, so that experimental evidences are mainly based on breakthrough curves, i.e. on the net current J(x,t) as a function of time collected at the end of an experimental setup (Berkowitz, Scher, 2001; Scher et al. 1991; Berkowitz, Scher, 1998; Berkowitz et al., 2001; Cortis, Berkowitz, 2004; Berkowitz, Scher, 1997; Levy, Berkowitz, 2003; Margolin, Berkowitz, 2002). 7. Conclusions In this paper we have presented a Monte Carlo approach to the transport of particles across an interface separating two media characterized by different physical properties. In particular, this scheme details how to switch from the jump lengths and waiting times distributions of one side to those of the other side. This scheme can take into account both normal and anomalous diffusion by suitably varying the underlying waiting times distribution. Guided by the Monte Carlo results, we have proposed an analytical formulation of the same problem: the most crucial point has been to establish a correspondence between the parameters of the Monte Carlo simulation and the coefficients appearing in the analytical solutions and also between the Monte Carlo rules and the boundary conditions at the interface. For the case of normal diffusion, most physical evidences suggest that particles preserve their local velocities while crossing the interface: we adopted this condition in the Monte Carlo scheme and we derived the corresponding analytical constraints. Under this assumption, resident concentration is found to smoothly vary across the boundary. On the contrary, when the Monte Carlo parameters are such that local velocities vary across the interface and the coefficients of the companion analytical solution are chosen accordingly, resident concentration is found to present a neat jump: this situation, although less common than the previous one, has been reported to occur in some physical contexts. As for the anomalous diffusion case, we found that imposing that particles preserve their local velocities while crossing the interface does not automatically lead to a smoothly varying resident concentration profile: instead, a neat bump is evident. In order to obtain a smoothly varying resident concentration profile, the analytical formulation suggests to impose the equality of the generalized velocities across the boundary. In both Fickian and anomalous diffusion cases, the resident concentration and current profiles given by the analytical and Monte Carlo formulations turn out to be in good agreement for any choice of the parameters. Acknowledgements This work was partially supported by the Italian Ministry of Education, University and Research (MIUR), through a project titled Transport of toxic and/or radioactive contaminants through natural and artificial porous media: models and experiments . References Barrat, J.-L., Chiaruttini, F., Kapitza resistance at the liquid solid interface. Molecular Physics, 101, (2003) 1605 Oldham, K. B. Spanier J., The Fractional Calculus, Academic Press, New York (1974) Berkowitz, B. Kosakowski, G., Margolin, G., Scher, H., Application of continuous time random walk theory to tracer test measurements in fractured and heterogeneous media, Ground Water 39 (2001) 593-604. Berkowitz, B., Scher, H., Anomalous transport in random fracture networks, Phys. Rev. Lett. 79 (20) 1997 Berkowitz, B., Scher, H., The role of probabilistic approaches to transport theory in heterogeneous media, Transp. Porous Media, 42 (2001) 241-2631 Berkowitz, B., Scher, H., Theory of anomalous chemical transport in random fracture networks, Phys. Rev. E 57 (1998) 5858-5869 Bouchaud, J.P., Georges, A., Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications, Phys. Rep., 195 (1990) 12 Briesmeister, J.F., Ed., MCNP – A general Monte Carlo N-particle transport code, Version 4C. LA-13709-M (2000), p. 2-79 Carslaw, H.S., Jaeger, J.C., Conduction of heat in solids, Clarendon, Oxford, 1959 Cortis, A. Berkowitz, B., Anomalous transport in “classical” soil and sand columns, Soil Sci. Soc. Am. J. 68 (2004) 1539-1548 Etoré, P., On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients. Elec. J. of Prob. 11 (2006) 249-275 Feller, W., An Introduction to Probability Theory and its Applications, Vol. 2, Wiley, New York (1971) Gorenflo, R., Mainardi, F., Moretti, D., Pagnini, G., Paradisi, P., Discrete random walks model for space-time fractional diffusion. Chemical Physics 284 (2004) 521-541 Gorenflo, R., Mainardi, F., Moretti, D., Pagnini, G., Paradisi, P., Fractional diffusion: probability distributions and random walk models. Physica A 305 (2002) 106-112 Hoteit, H., Mose, R., Younes, A., Lehmann, F., Ackerer, Ph. Three-dimensional modelling of mass transfer in porous media using the mixed hybrid finite elements and the random walk methods. Mathematical Geology 34 (4) 2002 Kilbas, A. A., Srivastava, H. M., Trujillo, J. J., Theory and Applications of Fractional Differential Equations. Elsevier, Amsterdam (North-Holland Mathematics Studies No 204) (2006) Klafter, J., Shlesinger, M.G., Zumofen, G., Beyond Brownian Motion, Phys. Today 49 (2) (1996) 33 Kreft, A., Zuber, A., On the physical meaning of the dispersion equation and its solutions for different initial boundary conditions. Chem. Eng. Sci. 33 (1978) 1471-1480 Krepysheva, N., Di Pietro, L., Neel, M.C., Fractional diffusion and reflective boundary conditions, Physica A 368 (2006a) 355-361 Krepysheva, N., Di Pietro, L., Neel, M.C., Space fractional advection-diffusion and reflective boundary conditions, Phys. Rev. E 73 (2006b) 021104//1-9 LaBolle, E.M., Fogg, G.E., Tompson, A.F.B., Random-Walk simulation of transport in heterogeneous porous media: Local mass conservation problem and implementation methods. Water Resources Research 32 (3) 583-593 (1996) LaBolle, E.M., Quastel, J., Fogg, G.E., Diffusion theory for transport in porous media: Transition-probability densities of diffusion processes corresponding to advection-dispersion equations. Water Resources Research 34 (7) 1685-1693 (1998) LaBolle, E.M., Quastel, J., Fogg, G.E., Gravner, J., Diffusion processes in composite porous media and their numerical integration by random walks: generalized stochastic differential equations with discontinuous coefficients. Water Resources Research 36 (3) 651-662 (2000) Leij, F.J., Dane, J. H., van Genuchten, M. Th., Mathematical Analysis of One-Dimensional Solute Transport in a Layered Soil Profile. Soil Sci. Soc. Am. J. 55 (1991) Lejay, A., Martinez, M., A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients. Ann. Appl. Probab. 16 (1) (2006) 107–139 Levy, M., Berkowitz, B., Measurement and analysis of non-Fickian dispersion in heterogeneous porous media, J. Contaminant Hydrology, 64 (2003) 203-226 Mainardi, F. Raberto, M. Gorenflo, R. Scalas E., Fractional calculus and continuous- time finance II: the waiting-time distribution. Physica A 287 (2000) 468 Mainardi, F., Pagnini, G., Saxena, R.K., The Fox H functions in fractional diffusion. J. of Comp. and Appl. Mathematics, 178 (2005) 321-331 Margolin, G. Berkowitz , B., Continuous Time Random Walk revisited: first passage time and spatial distribution. Physica A 334 (2004) 46 Margolin, M., Berkowitz, B., Spatial behaviour of anomalous transport, Phys. Rev. E 65 (3) 031101 (2002) Marseguerra, M., Zoia, A., Monte Carlo simulation of anomalous transport in presence of a discontinuity and of an advection field. Submitted to Physica A Marseguerra, M., Zoia, A., The Monte Carlo and Fractional Kinetics approaches to the underground anomalous subdiffusion of contaminants. Annals of Nuclear Energy 33 (3) (2006) 223-235 Metzler, R. Klafter, J., The random walk’s guide to anomalous diffusion: a Fractional Dynamics approach. Phys. Rep., 339, (2000a) 1-77 Metzler, R., Klafter, J., Boundary value problems for fractional diffusion equations, Physica A, 278 (2000b), 107-125 Metzler, R. Klafter, J., The restaurant at the end of the random walk: recent developments in the description of anomalous transport by Fractional Dynamics. J. Phys. A: Math. Gen. 37 (2004) R161-208 Miller, K. S. Ross, B. An Introduction to the Fractional Calculus and Fractional Differential Equations, Wiley, New York (1993) Montroll E. W. and Weiss G. H., Random walks on lattices. II. J. Math. Phys. 6, (1965) 167-181 Novakowski, K.S., An evaluation of boundary conditions for one-dimensional solute transport 1. Mathematical development. Water Resources Research 28 (1992) 2399-2410 Novakowski, K.S., An evaluation of boundary conditions for one-dimensional solute transport 2. Column experiments. Water Resources Research 28 (1992) 2411-2423 Ovaskainen, O., Cornell, S.J., Biased movement at a boundary and conditional occupancy times for diffusion processes. J. Appl. Prob. 40 (2003) 557-580 Paradisi, P., Cesari, R., Mainardi, F., Tampieri, F., The fractional Fick's law for non-local transport processes, Physica A 293 (2001) 130-142 Parlange, J.-Y., Barry, D.A., Starr, J.L., Comments on “Boundary conditions for displacement experiments through short laboratory soil columns”. Soil Sci. Soc. Am. J. 48 (1984) Podlubny , I., Fractional Differential Equations. Academic Press, San Diego, 1999 Scher, H. Montroll E. W., Anomalous transit-time dispersion in amorphous solids. Phys. Rev. B 12, (1975) 2455 Scher, H., Margolin, G. Metzler, R., Klafter, J., Berkowitz, B. The dynamical foundation of fractal stream chemistry. Phys. Today, Jan. 1991 (1991) 26-34 Scalas, E., Gorenflo, R., Mainardi, F., Uncoupled continuous-time random walks: Solution and limiting behaviour of the master equation, Physical Review E 69 (2004) 011107/1-8 Schwartz, R.C., McInnes, K.J., Juo, A.S.R., Wilding, L.P., Reddel, D.L., Boundary effects on solute transport in finite soil columns. Water Resources Research 35 (3) 671-681 (1999) Shlesinger, M.F., Zaslavsky, G.M, Klafter, J., Strange Kinetics, Nature 363 (1993) 31 Uffink, G.J.M., A random-walk method for the simulation of macrodispersion in a stratified aquifer. Proceedings of IAHS Symposia, 65 26-34, 18 th IUGG Assembly, Hamburg, 1985 Uffink, G.J.M., Analysis of dispersion by the random walk method. PhD dissertation, Delft University, The Netherlands, 1990 Van Genuchten, M.Th., Parker, J.C., Boundary conditions for displacement experiments through short laboratory soil columns. Soil Sci. Soc. Am. J. 48 (1984) 703-70819