_{1}

^{*}

The space-time fractional advection dispersion equations are linear partial pseudo-differential equations with spatial fractional derivatives in time and in space and are used to model transport at the earth surface. The time fractional order is denoted by *β*∈* * and

The development which has happened on the last twentyfive years on the fractional calculus opened many new applications on many fields such as physics, hydrodynamics, chemistry, financial mathematics, and some other fields. Actually a growing number of articles and books which are interesting on this field and its applications have appeared in these last 25 years (see for example: [1-5] and see also my thesis [

The behaviour of particles in transport under the earth surface is an important problem. For examples, the transport of solute and contaminant particles in surface and subsurface water flows, the behaviour of soil particles and associated soil particles, and the transport of sediment particles and sediment-borne substances in turbulent flow. There are many other examples in this field. The classical advection dispersion equation, ade, has been used to formulate such problems. The generalized fractional advection-dispersion equation, fade, has recently gotten an increasing interest from many scientists because it has many applications specially on studying the transport of passive tracers carried by fluid flow in a porous medium, see Benson, Meerschaert et al. [8-12]. In their work they gave applications and experimental results for the space-fade.

There is no unique solution for the space-time fractional diffusion processes but there are some attempts using different forms of the hyper geometric functions, as for example: in [

I am interested in this paper to find the approximate solutions of the space-time fractional advection equation, space-time fade, by adopting the backward GrünwaldLetnikov Scheme joined with the common finite difference methods. The space-time fade is considered as a diffusion process under the action of a constant force in a fractal medium with a memory. I study and numerically investigate the effect of the time fractional on the path of the particle motion as well as the effect of the spacefractional order for the three cases as: 0 < α < 1, 1 < α < 2, and. I compare between all these cases numerically. My numerical results are consistent with the results of [

The generalized fade reads

Here a, and b are positive constants representing the dispersion coefficient, and the average fluid velocity and it acts as the drift term to the right respectively. My aim is to give the approximation solutions of the space-time fad equations for all values of and. I study also the convergence of the approximation solutions to the solutions of the corresponding analytical solutions of the space-time fad equations in the Fourier-Laplace domain.

The used time-fractional derivative operator

is called Caputo fractional operator, see [

This equation is important for solving the fractional differential equations because it show the dependence on the initial conditions. Here is called the Riesz space-fractional differentiation operator. I adopt here the notation introduced by [

and must not be confused with a power of the first order differential operator (see [

For the proof of convergence in distribution I need to use the Fourier transform of which reads

while

This means, in the Zaslavski’ s notations,

From (2.2) - (2.4), one easily sees that in the case

Since, we can setwhich proves that the Riesz derivative is a symmetric fractional generalization of the second derivative. For more information about the Fourier transform and the pseudo-differential operators as semi groups of linear operators, see e.g. [26,27]. In my paper, I discuss the approximate solution of the Equation (2.1) for all values of and, to do so, I descretize and by the grid with,

. Here, and are the steps in space and in time, respectively, and is the number of steps at the x direction. Treating as a density of an extensive quantity (like mass, charge, solute concentration, or probability), the approximation of the collected quantity presents in a spatial cell at the instant by a clump,

For, I introduce the column vector

where. To proceed on the proof of convergence in distribution, one needs to use the method of generating functions, see [

for the sequence of clumps

. Using the initial conditions, I introduce the function

and applying the Fourier-transform, we obtain

Now, introduce the bivariate (two-fold) generating function

as a function of, where, for the sequence

Introduce the function and apply the Laplace-transform, one gets

From Equations (2.7) and (2.10), we deduce that if we replace by and by, in Equation (2.8), we get the Fourier-Laplace transform of the bivariate sequence which is obtained by collecting all the sequences in (2.9). This means

Our aim now is to prove that is related asymptotically to the Fourier-Laplace transform of which represents the analytical solution of Equation (2.1) for any values of and, and for a fixed and, as. So far, I will prove

for each case.

I describe in this section the classical partial differential ade and its proof of convergence in distribution. It is well known that the classical ade is a partial differential equation describing the solute transport in aquifers and it reads

here is the solute concentration. The conditions imposed on the solution are

With the initial condition. The classical ade is also interpreted as a deterministic equation with the probability function which describes the particle spreading away from the plume center of mass. The stochastic process described by Equation (3.1) is a Brownian motion with a constant drift [

If, then one has the diffusion of a free particle, that is, a particle in which no forces other than those due to the molecules of the surrounding medium are actingwhich reads. It has the solution. Consequently by using the Galilei transformation of the independent variables to, then the solution of (3.1), asis. In the Fourier domain

and hence fourth in the FourierLaplace domain, see [

Now descretizing (3.1) by the central symmetric difference in space and forward in time, one gets

Introduce the scaling relation

and for the positivity of all the coefficients of, one must put. Now let, one can write as

The discrete solution at Equation (3.5) describes also a random walk with sojourn probability of a particle at the point at the instant and it may jump either to the points, , or at the time instant, see [

The transition probabilities, and in Equation (3.6) satisfy the essential condition

Now one can use these transition probabilities to constitute a tridiagonal, P matrix, in which

. Therefore, Equation (3.5) is written in the following matrix form

Introduce the row vector, defined as

In order to find the explicit discrete solution of Equation (3.1), I have to take the transpose of each sides of the matrix Equation (3.7) and rewrite it as

and for the numerical calculations, it is convenient to write the stochastic matrix in the form

here I is the unit matrix and is a matrix whose rows are summed to zero. In Section 7, I give the evolution of for different values of. Now I am going to prove that the discrete solution at Equation (3.5) converges to the Fourier-Laplace transform of Equation (3.1). Rewrite Equation (3.5) as

Then multiplying both sides by and summing over all, to get

Multiplying both sides by and summing over all, one gets

The choice of the initial condition of the column vector satisfying that, guarantees that

Now, replace by and by, in Equation (3.12), to get

Now after using Taylor expansion and taking the limits as and, one gets

Compare this equation with Equation (3.2), then Equation (2.12), is satisfied for the classical ade.

In this section, I replace the first-order time derivative in Equation (3.1) by the Caputo fractional derivative,

, with. Then this generalized time-fractional advection-dispersion equation, fade, reads

For more information about the Caputo fractional derivative and its relations to the Riemann-Liouville, see [

To descretize, I utilize the backward Grünwald-Letnikov scheme which has been successfully utilizing at [19-24] for modelling and simulating the timefractional diffusion processes and the time-fractional Fokker-Planck equations.

Join this discretization with the common symmetric finite difference for and, then one has

Now introduce the scaling parameter

and solve for, one gets

where for ease of writing, I use and which has been originally introduced in [

with, and all,. Finally, and satisfy the relation

where, see [

to be positive, it is required that. Rewrite Equation (4.6) in the following form

This equation can be interpreted as a random walk with a memory, see [

where is matrix and is not a stochastic matrix as its rows are summed to, where. A useful numerical method to ease the computation is to write the matrix as

where is the same matrix defined in the last section, i.e. it does not depend on the value of then it does not depend on the time. In Section, I compare the evolution of for different values of and.

Now, I want to prove that the discrete solution (4.6), in the Fourier-Laplace domain as and, converges to the Fourier-Laplace transform of Equation (4.1). To do so, I have to adopt the initial condition which satisfies that, then. Then multiply each sides of Equation (4.4) by and sum over all, to get

Now proceed further, multiply each sides by and sum over all, to get

For manipulating the R.H.S., one needs to use the rule of multiplication of two sequences and, see Feller [

After using this rule, and put, then apply Taylor expansion, and take the limit as, one gets

Now, substitute, then again apply Taylor expansion, and take the limit as, you get

then multiply each sides by, one gets

Then I have proved the required aim for the timefractional ade.

In this section I consider the space-time fade, Equation (2.1). It is known that the space-fractional ade arises when velocity variations are heavy tailed and describe particle motion that accounts for variation in the flow field over entire system. The time fractional ade arises as a result of power law particle residence time distributions and describe particle motion with memory in time, see

[

where the symmetric Riesz Potential operator is defined as

where

The Fourier-Laplace transformation of Equation (2.1) reads

When descretizing the Riesz fractional operator

one must use a suitable finite difference scheme and exclude the case. To do so, I use the approximation of the inverse operators by the Grünwald-Letnikov scheme, see Oldham & Spanier [

where denotes the approximating Grünwald-Letnikov scheme which reads, see [4,18,19]

a)

b)

The shift in the index in Equation (5.5) is required to obtain a scheme with all coefficients are non-negative in the final formula for which gives schemes for simulating particle paths which results after replacing the second order space-derivative in Equation (4.1) by the Feller operator [

one must distinguish the descretization of with respect to the value of, as follows:

while

Now we adjoin the descretization of, with the descretization of, and with the finite sequence

. In what follows, I give the descretization of the space-time fractional ade for each case.

In order to ensure that all the coefficients ofI have to descritize by using the symmetric central scheme, so the discretization of Equation (2.1) reads

Adopting the scaling relation

and using Equation (5.1), one can solve Equation (5.9) for to get

To ensure that the coefficients of all, it requires that

Let us write Equation (5.11) in the form of a random walk, in which the walker is sitting at at and jumps to at, where

see [

The first two terms at the LHS of this equation represent the memory part and the other terms represent the diffusion under the drift term. I like to write as it represents the transition to the next step at the same point, represents jumping one step to the left, represents jumping one step to the right, and similarly jumping steps to the left or right. By using the identity (4.7), and the identity

then it is easily to prove that the summation of all the transition probabilities is one.

Now for the numerical calculations, and since there is a symmetric random walk, I adjust and ignore all the transitions outside this interval. So it is convenient to write the last equation in the form of a matrix form

Here is an elegant fifth diagonal matrix with its elements are computed as

where and. In the special case as, one recovers the well-known three point jumps of the classical ade because as, and hence the matrix be the same as matrix defined in Equation (3.7). In Section 7, I will plot the path of the particle representing by (5.13). Now, I will prove that the discrete solution at (5.9) converges in the Fourier-Laplace domain to the Fourier-Laplace transform of Equation (2.1) as. To do so, multiplying both sides of Equation (5.9) by where and sum over all, where it is easy to prove

then one gets

The second step is to multiply both sides of this equation by and sum over all,

then put and and use the previous results. The only new part on the proof is as, one can easily prove by using the fundamentals of complex analysis that

So far

Finally multiply each side by and substitute the value of, and compare with (5.2), one gets

So, I proved the desired aim.

I have to descritize by using the backward difference scheme, in order to ensure that all the coefficients of are positive for. Then the descretization of the space-time fade in this case reads

After using the scaling relation (5.10), I am going to separate the coefficients related to from the last summation and rearrange Equation (5.16) as

To ensure that all the coefficients of are positive, the scaling relation must satisfy

It is known as, the value ofand that ensure that. I can use the same symbols of the last subsection to write this equation in the form of the random walk as

Again the summation of all the transition probabilities of this equation over all is one. As in the previous, Equation (5.17) can be written in the same matrix form (5.13) where the diagonal elements of the matrix are defined as (see (5.19) below)

In the numerical calculations, I plot the path of the particle for different values of. I choose the values of according to condition of each case. The situation as is the same as the previous cases. Now, I am going to prove the convergence of the discrete solution (5.18) in the Fourier-Laplace domain to the solution of the space-time fade as. To do so, one has to shift the indices of and rewrite Equation (5.18) as

Again multiply each side by and sum over all j and use all the identities of the last section, to get

Replace, where, it is clear that

So one gets

Now, multiply both sides by, , and sum over all, to get

Finally, replace by, solve for, and multiply both sides by, you will get the desired aim, Equation (5.2). Then the discrete solution converges to the solution of its corresponding space-time fade in the Fourier-Laplace domain.

The Cauchy fractional ade needs a special treatment. I am going to prove first the convergence of the discrete solution. Therefore, the Fourier-Laplace transformation of the space-time fade (2.1), as, reads

This case is related to the Cauchy distribution and one cannot use the Grünwald-Letnikov discretization of

at Equations (5.7) and (5.8) because the denominator is zero and for. Instead of GrünwaldLetnikov discretization, one must use the descretization used in [

[

, , in Equations (5.7) and (5.8) by for, and for. Therefore by using the scaling parameter

one gets

As the previous cases, can be written in the form of a random walk with a memory as

To have all the coefficients of are positive, one should restrict the values as. Following

[

The numerical result of this model is discussed in the next section. Noting that the coefficient and play here a significant role. To prove the convergence, multiply each side of Equation (6.3) by and sum over all, then use the identity

to get

Replace, z by, and use the identity

Finally replace, by and take the limit as, you get, Equation (6.1).

In this section, I give the numerical approximate solutions for Equation (2.1). I give the evolution of with different values of such, different values of the space fractional order and different values of the time fractional order. I fix the values of while and, with the initial condition as it must satisfy

. Since, then the iteration index while is calculated from the scaling parameter of the specified model and its values are varying according to the restriction put on. Since I used the explicit discrete scheme, therefore, sometimes I need a huge number of steps for calculating specially for. I calculate most of the numerical results for and but as and 1 < α < 2 I used to ensure that all the elements must be, because is probability transition matrices. To ease comparing the numerical results, I wrote the values of, and their corresponding with the values of at the figure. The classical case, i.e. as and, is plotted at Figures 1 and 2 and one can observe how rapidly the paths diffuse as the time increases. The time-fractional ade is simulated at Figures 3 and 4, i.e. for and. Figures 5, 6 are corresponding to the space-time fractional ade as 0 < α < 1 and. While Figures 7 and 8 are denoted for 0 < α < 1 and. Figures 9 and 10 are corresponding to 1 < α < 2 and at Figures 11 and 12 for 1 < α < 2 and. Finally, I plot the space-time fractional ade as and is plotted at Figures 13 and 14. The numerical results of this paper as and are consistent with the results at [

discussed by many authors. The pates as are different than the paths corresponding to. The pathes for need a huge number of time steps to calculate them as I use the explicit difference methods.