A Discrete Method Based on the CE-SE Formulation for the Fractional Advection-Dispersion Equation

We obtain a numerical algorithm by using the space-time conservation element and solution element (CE-SE) method for the fractional advection-dispersion equation.The fractional derivative is defined by the Riemann-Liouville formula.We prove that the CESE approximation is conditionally stable under mild requirements. A numerical simulation is performed for the one-dimensional case by considering a benchmark with a discontinuous initial condition in order to compare the results with the analytical solution.


Introduction
Recently, the analytical and numerical framework of fractional differential equations has attracted much attention because of its extended use in the modeling of nonlocal transport phenomena which appear in many fields like magnetic plasma, turbulence, and flows in porous media [1][2][3][4][5].A nonlocal transport can lead an anomalous diffusion because of large tracer displacements (super-diffusion) or trapping structure vortices (sub-diffusion).The standard local diffusion description is linked to Gaussian processes and diffusion differential models are obtained.However, a nonlocal transport flow presents a non-Gaussian particle motion and in this case fractional diffusion models are derived [6][7][8].One of the most considered proposals for the modeling of the mass flow in media with anomalous diffusion is the fractional advection-dispersion equation (FADE); see, for example, [7][8][9][10].
Different definitions for the fractional derivative can be used to describe the FADE, and the most common ones are the Caputo, Riemann-Liouville, and Grünwald-Letnikov derivatives [11][12][13].Here, we consider the one-dimensional FADE given by where  and  represent the space and the time, respectively,  = (, ) is the solute concentration,  is the fractional dispersion,  is the skewness of the diffusive transport taking values in the interval [−1, 1], and  is the fractional order.In our case we consider 1 <  < 2 and by definition of the ceil function we get  = 2.We can rewrite the fractional equation (1) in a more compact form as follows: by defining the fractional operator ∇   as The analytical solution of (3) is only known under particular initial and boundary conditions, so numerical solutions are calculated and studied for most of the problems.In the last years, a lot of works have been focused on the development of numerical algorithms for fractional differential equations.The first efforts are based on extensions of standard methods used for the nonfractional case like the finite difference method (FDM) [14][15][16][17], the finite volume method (FVM) [18,19], and the finite element method (FEM) [20,21].More recently other methods have been developed, for example, the spectral method [22], the collocation method [23], or a combination of both [24].Generally, the FDM is constructed using the shifted Grünwald formula for the discretization of the fractional derivative, the FVM under a conservative formulation is based on the discretization of the integral using the Riemann-Liouville definition or the Caputo definition of the fractional derivative, and the FEM and collocation methods are variational techniques based on the discretization of the integral on some interpolating polynomial that approximates the solution.Therefore, all these techniques generate algorithms with dense coefficient matrices for which it is complicated to analyze their stability and a high computational cost to solve them is required.
As an alternative to avoid the discretization of the integral associated with the fractional derivative, in this paper we propose a numerical algorithm based on the space-time conservation element and solution element (CE-SE) method developed by Chang and To [25].For this method, the solution is approximated by a first-order Taylor expansion which must fulfill the integral form of the equation in the conservation elements and the differential form in the solution elements.Therefore, the space-time domain is discretized two times by: a rhomboid mesh and a rectangular mesh.The theoretical framework of the CE-SE method for the nonfractional case is well established and it has been applied to a large number of real problems with successful results [26][27][28][29][30][31][32][33].An important property of the fractional differential operators is the nonlocality which is achieved with the discretization of the integral associated with the fractional derivative.In our case we seek to keep the nonlocality by providing information of the numerical solution at nearby points in both meshes and using nonlocal approximations of the standard derivatives.It is important to remark that we consider this work as a first step to open a venue for a different way to approximate the fractional differential equations like the FADE equation.
The paper is organized as follows.In Section 2, we give the discretization of the space-time domain for the CE-SE method, its extension to the fractional advection-dispersion equation (4), and the von Neumann stability analysis of the proposed algorithm.In Section 3, we present numerical experiments of the fractional CE-SE method for the FADE with an initial condition type shock and it is compared with the analytical solution of the problem.Finally, this work is completed with a section of conclusions and an appendix with the calculation of the fractional Riemann-Liouville of a linear polynomial of two variables.By definition, SE(, ) is the interior of the space-time region bounded by the opened rhombus in Figure 2. In each SE, the flow variables are assumed to be continuous with possible discontinuities at the interface between the two rhombuses.The second mesh is given by nonoverlapping rectangular regions, CE(, ), which fill the computational domain, see Figure 2. In each CE, a local space-time flux balance is enforced so a global flux conservation relation is satisfied over the entire space-time domain.Note that each CE consists of three different SEs and if we compute the solution in a integer time step then the spatial index is running in the half-integer points and vice versa, see Figure 1.[25,34] for the one-dimension FADE.First, the solution  = (, ) of ( 3) is approximated in each SE(,  + 1/2) by a first-order Taylor expansion of two variables which is denoted by û(, ) and  is given by

The Fractional CE-SE Method. Now, we construct an explicit algorithm based on the CE-SE method
and  +1/2  =   (  ,  +1/2 ) for all  and  integers, we need to know these coefficients to obtain the approximation (5).
Figure 2: The graphic to the left shows a conservation element CE(, ) and the graphic to the right shows the associated solution element SE(, ).Assuming that they are known in the time step , we are going to calculate the coefficients  +1/2  and  +1/2  imposing that the numerical approximation (5) verifies the integral form of (3) in each CE(,  + 1/2).We split the rectangle CE(,  + 1/2) in two subrectangles and solve the integral form in each subrectangle; see Figure 3. Consider for  = 1, 2. Applying the Green theorem to the integral of the left hand side of ( 6) we have that where   : [0, 4] → R 2 is the parameterization of CE  given by: Substituting the values of   () and    () for each  and replacing û and its partial derivatives by the coefficients 's, 's, and 's in (7), we obtain Now, we rewrite the integral of the right hand side of ( 6) as a function of the coefficients 's and 's.For that, it is necessary to calculate the Riemann-Liouville fractional derivative ∇  of the function û defined in (5).This calculation is carry out in the Appendix.Let [, ] be a spatial interval, and for any  ∈ N 2 we write   −  = Δ and  −   = ( − )Δ ( the right-end index of N 2 ) and substituting them in (A.5), we get where Since the approximation û defined in (5) must fulfill the FADE (1) in each SE, then the coefficient  can be obtained by for  =  ± 1/2.Denoting by  = (Δ/Δ) the Courant number and by  = (Δ/(Δ)  ) the anomalous diffusion number, we can rewrite the previous formula for the coefficient    as Using (15) for  =  ± 1/2 and substituting ( 9) and ( 11) into (6) for the case of  = 1 and for  = 2 substituting ( 10) and ( 12) into ( 6), we have the following discrete system: The system ( 16) can be written in matrix form as follows: where , , and  are the coefficient matrices that are given by The CE-SE method is a family of schemes [35], and the developed algorithm ( 16) is called the a-scheme which is the core algorithm of the CE-SE family.Next, we are going to prove the convergence of this fractional a-scheme verifying its consistency with the linear FADE (1) and obtaining its stability conditions.For the stability analysis, we use the von Neumann technique [36] which is a local analysis.Thus, the boundary effects are not considered which is a reasonable assumption since in most cases the numerical instability appears in the interior nodes of the discretization of the domain.

Stability and Convergence
Analysis.We start the von Neumann stability assuming that the solution and its partial derivative with respect to  can be expressed as a sum of eigenmodes at each node of the mesh: where  is the spatial wave number and  = () is a complex number.Now, we substitute (19) in the discrete system (17) and obtain and multiplying both sides by  −  −  yields with Taking the square of Euclidean norm in both sides of ( 21), we have where  is the amplification factor of each mode in the expansions (19).By calculating we have where  ±1/2 = 1− ±1/2 ,   = 1+  ,  ±1/2 = − + 4 ±1/2 , and   =  − 4  .The stability of the method is achieved if the magnitude of the amplification factor is less than 1.By definition, the values of the positive coefficients   and  ±1/2 are close; thus the following condition is fulfilled By definition of the coefficient   , we get ( 2 +1/2 +  2 −1/2 + 2) ≈ 2( 2  + 1); thus || < 1 if we take Δ small enough.Now, we analyze the local truncation error to obtain the consistency of the fractional CE-SE method (17) with the FADE (1).Given that we assume the flow variables ,   , and   to be continuous in each SE and considering that the numerical solution û is a first-order Taylor approximation (5), we conclude that the fractional a-scheme ( 17) is second-order accurate in space and time.Given the above results, we have the following theorem using the Lax Equivalence theorem [36] which states that For a uniformly solvable linear finite difference scheme which is consistent with a well-posed linear evolution problem, the stability is a necessary and sufficient condition for its convergence.17) is convergent for a Δ small enough.

Theorem 1. The fractional CE-SE scheme (
It is well known that the a-scheme of the CE-SE family is a nondissipative algorithm.We are interested in flow problems with discontinuities, so we need to introduce some numerical dissipation in the approximation.In the next section, we modify the a-scheme using another way to approximate the derivative   , and the modified algorithm is known as a-scheme.

An Approximation of the Derivative for Discontinuities.
In order to introduce information about the behavior of the solution in the calculation of the partial derivative of û with respect to , we propose another approximation for the coefficient  +1/2  .By adding both equations of the system (16), we can rewrite the fractional CE-SE algorithm as follows: As a first approximation, we consider the coefficient  +1/2

𝑗
given in [25,28]: where the ratio , with  a positive constant, can be considered as a sloper limiter coefficient which measures the smoothness of the function by using the forward difference and the backward difference: with û(  ,  +1/2 ) =  +1/2  .As the points ( ±1/2 ,  +1/2 ) are not solution nodes, we can obtain û( ±1/2 ,  +1/2 ) as follows: The estimation (27) for the partial spatial derivative has been used in nonfractional differential problems with nonsmooth solutions obtaining successful results [25,26,28].In the boundary nodes, the coefficient  +1/2 0 is calculated from the forward difference  + and the coefficient  +1/2  is obtained with the backward difference  − .Remark 2. Notice that we have used a weighted finite difference approximation to approximate the coefficient  in the step  + 1/2, but other approximations can be used.We think that this can be a way to introduce nonlocal information in the calculation of the numerical solution.

Numerical Results
For the fractional case it is important to remark that there are few benchmark problems with known analytical solution.So in this section, we consider (3) with a jump condition of the shock type: with  0 a positive constant function and analyze the performance of the fractional CE-SE method ( 26)-( 27) by comparing with the exact solution, given by Sousa [37].So, taking  = 0 and the boundary conditions (−∞, ) =  0 and (∞, ) = 0, the exact solution is where  = ( − )/(| cos(()/2)|) 1/ and the cumulative probability function is defined by with For all  ̸ = 0 and  ≤ 0, the cumulative probability function is obtained using the identity   (−) = 1 −   ().Next, we show the numerical results of the fractional a--scheme for the FADE (1) with initial condition (30) using the same test problem given in [37].
Example 3. We consider the solution of ( 1) with (30) in the interval [ −4, 4] for  0 = 1,  = 0.5, and  = 0.2 and the boudary conditions are (−4, ) = 1 and (4, ) = 0.In Figures 4 and 5, we show the profile of the analytical solution (31) and the CE-SE approximation ( 26)-( 27) at  = 1 for several values of the parameter .Moreover, we compare the global error of the fractional a- scheme versus some finite difference methods that appear in [37] for a fixed time   using the  ∞ norm given by where   is the difference between the exact and the numerical solutions; see Table 1.Notice that for the three schemes as  is closer to one the numerical error is larger than in the cases as  is closer to two, since for smaller  the solution loses differentiability.This fact is visible at the extremes of the shock in Figure 4, where we display the CE-SE behavior for  = 1.2 and 1.4.Although the three methods present similar numerical errors, the CE-SE method has the advantage of holding a low computational cost compared to the other two methods.Finally in Figure 5, we show the accuracy of the fractional CE-SE method as the step size is diminished.

Conclusions
In this work we have constructed a discrete method under the CE-SE formulation for the fractional advection-diffusion equation described by the Riemann-Liouville derivative.The new fractional method avoids the fact that the calculation at the new point depends on the function values at points of further regions at that point.We have also proved the convergence of the fractional CE-SE method and show its performance for a benchmark problem with a shock initial condition.Future work will focus on the analysis of different proposals to approximate the partial spatial derivative of the solution to include the nonlocal information of the flow.We also consider extending the fractional CE-SE method to two dimensions in order to cover a wider range of flow problems important.(A.5)

Figure 3 :
Figure 3: In both graphics (a) and (b) different subdivisions of the conservation element, CE(,  + 1/2), are shown.In (a) the two regions CEs that the CE(,  + 1/2) is divided into are shown whereas in (b) the four SEs regions are shown.

Figure 4 :Figure 5 :
Figure 4: Numerical solution (red line) of methods (26)-(27) taking Δ = 0.2 and the exact solution (blue line) of the FADE.The graphic to the left shows the results for  = 1.2 and the graphic to the right shows the results for  = 1.4,for both cases Δ = 0.2.