Two-dimensional nonlinear time fractional reaction–diffusion equation in application to sub-diffusion process of the multicomponent fluid in porous media

In the present article, an efficient operational matrix based on the famous Laguerre polynomials is applied for the numerical solution of two-dimensional non-linear time fractional order reaction–diffusion equation. An operational matrix is constructed for fractional order differentiation and this operational matrix converts our proposed model into a system of non-linear algebraic equations through collocation which can be solved by using the Newton Iteration method. Assuming the surface layers are thermodynamically variant under some specified conditions, many insights and properties are deduced e.g., nonlocal diffusion equations and mass conservation of the binary species which are relevant to many engineering and physical problems. The salient features of present manuscript are finding the convergence analysis of the proposed scheme and also the validation and the exhibitions of effectiveness of the method using the order of convergence through the error analysis between the numerical solutions applying the proposed method and the analytical results for two existing problems. The prominent feature of the present article is the graphical presentations of the effect of reaction term on the behavior of solute profile of the considered model for different particular cases.


Introduction
It is not justified to categorize the fractional calculus theory as a young science. The origin of fractional calculus is as old as classical calculus itself. During the past few decades it has become the focus of interest of many disciplines of science and technology which provides an excellent and efficient tool for modeling and describing various scientific and complex engineering phenomena such as aerodynamics, polymer rheology, electrodynamics of complex medium, fluiddynamic traffic model [23,33,34,44,56] etc. In the last few decades there have been a lot of research on the applications of fractional calculus theory to various scientific fields ranging from the physics of diffusion and advection phenomena to control system of finance and economics. Fractional differentiation has a lot of advantages on the simulation of dynamical systems and physical phenomena as compared to integer order differentiation, due to its non-Markovian and non-local behaviors. In many cases it is not possible to model the known equations in fractional order forms. For this some basic physical postulates are to be satisfied before giving its shape into fractional order system. Therefore every equation can not be generalised simply by replacing the integer order derivative by fractional order derivative.
From literature survey, it is seen that the fundamental ideas of the algorithms correlated to the ideas proposed by authors of [3,57] and Bhrawy [6] have been used to develop the efficient and accurate algorithms for the purpose of solving partial differential equations [32]. Additionally, for finding the numerical solutions of linear fractional differential equations (FDEs), Bhrawy et al. [7] have developed an operational matrix with the popular Laguerre polynomials for the fractional order integration and have revised the generalized Laguerre polynomials on the semi-infinite intervals. In [1], the authors have introduced an operational matrix for the fractional order derivative for solving the linear and non-linear FDEs with given initial conditions. A significant portion of fractional calculus theory is the two-dimensional fractional differential equations and integral equations which have further been studied in the articles [43,[45][46][47][48]. In [12], the authors have solved the two-dimensional fractional percolation equation in a non-homogeneous porous med-ium. In the present article, an accurate and efficient numerical technique has been proposed to solve two-dimensional fractional order reaction-diffusion equation arising in porous media.
Several kinds of flow problems like percolation flow problems, seepage flow problems have been examined and discussed in many aspects of research fields including groundwater hydraulics, seepage hydraulics, fluid dynamics and groundwater dynamics in the porous media [13,26]. The two-dimensional flow of the solute in the porous medium is governed under the hypothesis of continuity and Darcy's law is given by J. Bear and A. Verruijt [4] in the equation where A 0 denotes the specific storativity, which is approximated by the effective porosity or specific yield; P denotes the pressure head, the intrinsic permeability of the medium is given by K ¼ kqg l , which is the hydraulic conductivity having the components in the x, y directions as k x and k y ; q denotes the water density ; l be the viscosity; and g the gravitational acceleration; X is the percolation domain.
Since there is a large digression from the ordinary Gaussian diffusion, therefore the underlying superposition is not appropriate in tracing the movement of the solute concentration in a porous medium which is nonhomogeneous. Additionally, a supplementary realistic model can be obtained by the consideration of the application of fractional order density gradient to recapture the non-homogeneity of the porous medium. The Modified Darcy law was proposed by the authors of [22] with Riemann-Liouville fractional order derivative. Under the assumption of the continuity of seepage flow and the Darcy law, the percolation equation (1) is deduced into the the following forms as these assumptions are not valid for real seepage flow. The differential system with ordinary derivative does not quantify the flow of fluid with spatial path dependency and/or time memory. This limitation in the application of ordinary differential system has motivated the authors of [22] to model the fractional Darcy law.
In the special case, the usual Darcy law can be obtained by putting a 1 ¼ a 2 ¼ 1: Einstein's theory of Brownian motion reveals that the mean square displacement (MSD) of a particle moving randomly is proportional to time, which also can be justified for the case of a simple integer order linear diffusion equation. But as the research on fractional calculus progresses, it is found that the MSD for the anomalous case i.e., for time fractional diffusion equation, increases slowly with time. For the MSD is hX 2 ðtÞi $ t a ; where 0\a\1 is the anomalous diffusion exponent. The equation represents an evolution equation which generates the fractional Brownian motion, a generalization of Brownian motion. Thus it is seen that for the diffusion model, if the integer order time derivative is replaced by the fractional order time derivative, it changes the fundamental concept of evolution of foundation of physics. The physical meaning of the fractional order time derivative related to the statistics is the waiting time in accordance with the Montroll-Weiss theory. Hilfer and Anton [25] have showed that Montroll-Weiss continuous time random walk (CTRW) with a Mittag-Leffler waiting time density is equivalent to a fractional order master equation. Later, Hilfer [24] explained that this underlying CTRW of the model is connected to the time fractional diffusion equation in the asymptotic sense of long time and large distance. Thus random walk approach is needed to simulate diffusive phenomena of a fractional order equation. Gorenflo et al. [18] stated that the time fractional order diffusion equation generates a class of symmetric densities whose moments of order 2m are proportional to the ma power of time. Thus classes of non-Markovian stochastic processes can be obtained, which exhibit slow anomalous diffusions. By using fractional order Fokker-Plank equation approach, Metzler et al. [35] have shown that anomalous diffusion is based upon Boltzman statistics. Many researchers have used fractional equations during description of Levy flights or diverging diffusion. The tool is very powerful in modeling multi scale problems, characterized by wide time or length scale. The fractional order differential operator has the characteristic of non-local property, which states that the future state not only depends upon the present state but also upon all of the history of its previous states. Due to the greater flexibilities, the fractional order models have gained popularity to investigate the dynamical system. The anomalous diffusion process seeks the attention of many researchers [5,8,14]. Anomalous diffusion can be easily seen within complex systems like diffusion process in porous media. The fluid particles undergoes in sub-diffusion process if a\1. A subclass of anomalous diffusive system is the fractional sub-diffusion equation (FSDE), which can be illustrated from the standard parabolic PDE by replacing the first-order time derivative with a fractional derivative of order a; 0\a\1 (Moodie and Tait [37]).
In the real world problems, the FSDE has been widely applied in several fields of research [11,20,30]. So finding the solutions of these equations have become increasingly interesting and popular. The authors of [17,36,39] have provided a typical explanation of the anomalous diffusion process based on the continuous time random walk. The process of anomalous diffusion for real complex physical systems can be well explained by fractional order diffusion models.
The variations of the solute concentration with respect to the column length and time depend upon the several species and factors in fluid. It is very hard to find the overshoots of probability density function for the several cases like conservative and non-conservative systems. In addition to this, interactions of the solutes with medium or other species play a significant role in the flow of the fluid through porous media. In order to analyze the diffusion of the solute in porous media in two-dimension, an attempt has taken in the present article to develop a two-dimensional subdiffusion physical model (4) with the prescribed initial and boundary conditions (5)(6)(7)(8)(9), where the fractional derivative is taken in Caputo sense. The fractional time derivative corresponds to long-time heavy tail decay. where a is an arbitrary fractional order, k denotes the reaction coefficient and c 0 be a constant. For the approximation of solute concentration we have used Kronecker product which is discussed later. When the ground water is mostly adulterate, then the resuscitate is considered to be very difficult and more expensive. A very careful approach and attention is very much necessary for describing the boundary conditions, problem domain and model parameters for using the numerical approach of groundwater model of the field problems. Hydrology is an interdisciplinary part of science and engineering, in which the topic of solute transport through the groundwater is included.
The generalized solute transport model in the porous media is the well known reaction-advectiondispersion equation (RADE) [15,27]. This equation has the combined effects of advection, reaction and dispersion process due to which the concentration of the solute is dispersed while transported down along the stream and sometimes it reacts with the medium through which it moves. From the mass balance principles, the advection-dispersion equation can be easily derived. The modeling of solute transport is very helpful for the prediction of solute concentration in rivers, streams, aquifers and lakes, which can be represented mathematically as where R denotes the reaction term of the species u, and the symbol r is defined in R n in terms of partial derivative operators as r ¼ ð o ox 1 ; . . .; o ox n Þ. In the right-hand side of the equation (10), the first term describes the dispersion phenomena, the second term represents the advection process and the last term is for the reaction kinetics. If there is no reaction between the solute and the medium through which the solute moves, and also there is no kind of radioactive decay then this type of system is called as conservative system, otherwise it is non-conservative. In the case of non-conservative system, the last term of the equation (10) is encountered. If there is only diffusion process which is responsible for the movement of the solute concentration, then the above equation is known as diffusion equation.
In this article a drive has been taken to develop the operational matrix with Laguerre polynomials to use in collocation method to find the numerical solution of two-dimensional fractional order non-linear reactiondiffusion equation. We have introduced an efficient way to approximate the non-linear multiple-order fractional boundary value problem with the help of Laguerre spectral collocation method to obtain the numerical solution with the use of Laguerre operational matrix.
The article is organized in following ways. Section 2 contains some important properties and definitions of fractional calculus. A brief overview of the conservation principle is given in Sect. 3. In Sect. 4, some outlines of Laguerre polynomials are given. In Sect. 5, an overview of Kronecker product and its properties are given. In Sects. 6 and 7, the Laguerre operational matrix for fractional order derivative in Caputo sense is given to solve our proposed model. The detail of convergence analysis of the proposed numerical scheme is given in Sect. 8. In Sect. 9, few examples to illustrate the accuracy and applicability of the proposed method are considered and comparisons between the numerical results obtained by our proposed method with the existing analytical solutions are shown through the graphs and tables. Numerical results and discussion of the considered two-dimensional nonlinear model are given in Sect. 10, which is followed by the outcomes of the overall research work given in the section Conclusion.

Preliminaries
In this section of the manuscript, the basic definitions and properties of variable order derivatives are given. The Riemann-Liouville's derivative has some drawbacks in the modeling the mathematical modeling for some real and physical phenomena. Fractional order derivative operator D g of the given order g [ 0 of a function f(t) in the Caputo sense is given by [29,41] if r À 1\g\r: In the above expression r is the integer number. Few properties of Caputo fractional derivatives are: where C is an arbitrary constant.
if e 2 N [ f0g; e ! dgeore 6 2 N; e [ bgc; where dge be the ceiling function and bgc be the floor function.

The general conservation principle
In this section of the article, a fluid continuum in flowing state is considered. In general a multicomponent fluid has been considered in such a way that it follows the continuum approach and each component in itself is a continuum filling of the entire space, as in Fig. 1. Let H a be the initial amount of the the extensive fluid property (e.g., mass) is contained in a volume V within a portion of space at time t and enclosed by a surface S. Consider the position of the system at time t þ 4t, in the time interval 4t, the system moves and gets deformed. At time t, the system is of surface area S and volume V, whereas the system at time t þ 4t is of volume V 0 and surface S 0 . These are composed of the regions: The region V 2 is common to V and V 0 . In Lagrange point of view, the temporal rate of change ðD=DtÞH a of H a for the moving system is given by The general conservation principle of a property for the species 0 a 0 is given by where I a is the temporal rate per unit volume and V H a is the velocity field.

Mass conservation of species
Mass conservation of the fluid across a porous material involves the fundamental principle that total mass flux entered minus total mass flux out equals to the increase in mass amount stored by the medium. This means that the total mass of fluid is conserved. For a system of multicomponent fluid of a species 0 a 0 , substituting h a ¼ q a and V H a ¼ V a in above equation, we get where I a is the rate at which mass of the species 0 a 0 produced per unit volume during the chemical reaction in the system. From the Eulerian point of view, above equation is the equation of continuity (or equation of mass conservation) [21].
If we assume the case of stationary fluid and I a ¼ 0 then from above, we deduce the Fick's second law for diffusion for binary fluid as where D ab is the diffusion coefficient.

Laguerre polynomials and its some properties
The involvement of Laguerre polynomial in several aspects of engineering and applied mathematics seeks the attention of many researchers and establishes it as a strong tool for finding the numerical solution of Fig. 1 Fluid continuum in the flowing process integer as well as fractional order partial differential equations (PDEs). The l th degree Laguerre polynomial is defined as [40] L l ðtÞ ¼ 1 l! e t o l t ðt l e Àt Þ; l ¼ 0; 1; :::; The Laguerre polynomial of degree l on the interval I is simplified from above as l!ðÀ1Þ k ðÀk þ lÞ!ðk!Þ 2 x k ; k ¼ 0; 1; :::; The following recurrence relations are satisfied by Laguerre polynomial The set of all Laguerre polynomials form an orthogonal system in L 2 r ðIÞ i.e., Z where rðtÞ ¼ e Àt is weight function and d lm is the Kronecker delta function.

Kronecker Product and its some properties
In this section some basic and fundamental properties of the well known Kronecker product are defined [9,19]. Furthermore, we try to give some overview of the applications and properties of the Kronecker Product, which is an operation between two arbitrary size matrices represented by the symbol .
Definition The Kronecker product of the matrix A of order p Â q and B of order m Â n is again pm Â qn order block matrix, which is denoted by A B and is defined by Moreover, If A and B represent some linear transformations on vector spaces i.e., W 1 ! V 1 and W 2 ! V 2 respectively, then A B represents the tensor product of these two maps, W 1 W 2 ! V 1 V 2 .

Basic properties
Since Kronecker product is a special case of tensor product, therefore it can satisfy the property like associativity and bi-linearity. Also, it has a lot of interesting and effective properties, few are stated below: where I, J and K are arbitrary matrices and k is a scalar.

Laguerre operational matrix for fractional order differentiation
The collection of Laguerre polynomials form a basis set for L 2 r ðIÞ: So any function fðtÞ 2 L 2 r ðIÞ can be expanded as fðtÞ ¼ X 1 g¼0 b g L g ðtÞ; m ¼ 0; 1; 2; :::; ð28Þ where the unknown constants b g can be calculated as In the terms of first ðr þ 1Þ Laguerre polynomials the above approximation is reduced to where B T ¼ ½b 0 ; b 1 ; . . .; b r is the matrix of unknown constants and rðtÞ ¼ ½L 0 ðtÞ; L 1 ðtÞ; . . .; L r ðtÞ T is the Laguerre vector. Again a function of two variables fðx; tÞ 2 L 2 ½0; 1 can be expressed in terms of Laguerre polynomials as fðx; tÞ ' f r ðx; tÞ ¼ X r where the unknown constant matrix B ¼ ½b g;h is known as Laguerre coefficient.
In similar manner, we can approximate an arbitrary two-dimensional function fðx; y; tÞ 2 L 2 ½0; 1 in terms of Laguerre polynomial as fðx; y; tÞ ' f r ðx; y; tÞ ¼ X r where the unknown constant matrix B ¼ ½b g;h;i is the Laguerre coefficient. These coefficients can be determine by using the initial and boundary conditions. In order to derive the novel operational matrix we have provided the following lemma.
Lemma 1 Consider the fractional derivative operator D l of fractional order l in Caputo sense. Let rðtÞ be the Laguerre vector defined as in equation (30), then we have following results: C 0 D l t L r ðtÞ ¼ 0; r ¼ 0; 1; . . .; dle À 1; 0\l 1: Proof This result can be easily proved by using the basic properties of fractional Caputo derivative given in the equation (13). h We are going to derive the operational matrix of orthogonal Laguerre polynomials [40]. On solving the equations (28)- (29) and equation (20), we get the following expression as As set of Laguerre polynomials is a basis so we can write the term t gÀl as where the constant coefficient b h is given by ðÀ1Þ a h!Cðg À l þ a þ 1Þ ðh À aÞ!ða!Þ 2 : ð36Þ Using the equations (35)-(36) into the equation (34), we have C 0 D l t L l ðxÞ ¼ X r h¼0 P l ðl; hÞL h ðtÞ; l ¼ dle; :::; r; In the above equation the term P l ðl; hÞ is given by the following equation ðl À gÞ!Cðg À l þ 1Þg!ðh À aÞ!ða!Þ 2 : The vector form of the equation (37) is written as where H l is an operational matrix with fractional derivatives of r þ 1 Â r þ 1 order and is given as In above equations I be the ðr þ 1Þ Â ðr þ 1Þ identity matrix. Now from prescribed initial and boundary conditions (5)-(9) with the aid of equation (44) Now, we collocate Eq. (4) with the help of equations (48)(49)(50)(51)(52) at points x b = b r for b ¼ 0; 1; 2; . . .; r, y r = b r for b ¼ 0; 1; 2; . . .; r and t r = b r for b ¼ 0; 1; 2; . . .; r.
Using these collocation points and fractional operational matrix, a nonlinear system of algebraic equations are found. On simplifying this system we can calculate the Laguerre coefficients matrix B which can further be used in Eq. (44) to find the approximate numerical solution.

Error bound of the approximation
In this section the upper bound of the error of proposed numerical approximation is estimated. ; if À 1\h 0; x ! 0; l ¼ 0; 1; 2; ::: where x; y; t ! 0; f ; g; m ¼ 0; 1; 2; . . .; where v mbr ¼ P r m¼0 S a ðb; mÞ: Proof In view of equation (43) We can write a order partial derivative of u(x, y, t) and u r ðx; y; tÞ w.r.to t as From the above equations, we can write Using (37), the above equation reduces to Applying equation (54), we get Using the formulae S1 and S2 given in [31] and subtracting the truncated series from the infinite series, bounding each term in the difference, and summing the bounds complete the proof of the theorem. Laguerre fractional operational matrix has been applied on some existing fractional order problems to find their approximate numerical solutions. The obtained numerical results have been compared with the existing analytical solutions through finding the absolute error for approximation to illustrate the capability and validity of the proposed numerical scheme. All the numerical calculations are carried out by using the software MATHEMATICA 11.3. The absolute error is defined by [16] as E abs ðx; y; tÞ ¼j u exact À u num ðx; y; tÞ j : Maximum error for the problem in The order of convergence for two successive approximation N 1 and N 2 is also calculated to validate the higher order accuracy of the numerical scheme which is given as where Max E abs ðNÞ is the maximum absolute error during the approximation of degree N.
Example 1 Let us consider a two-dimensional nonlinear non-homogeneous fractional order advectiondiffusion equation as with f ðx; y; tÞ ¼ Àx ffi ffi t p À y ffi ffi t p À x 2 y 2 t 3 À 0:682734xy t 0:3 and boundary conditions uðx; y; 0Þ ¼ 0; uð0; y; tÞ ¼ 0; For the above example, the exact solution is uðx; y; tÞ ¼ xy ffi ffi t p for the fractional order parameters a ¼ 0:8. The absolute errors between the results obtained by our proposed numerical scheme and the exact solutions are exhibited through Tables 1 and 2 for y ¼ 0:5; 0 x 1 and x ¼ 0:5; 0 y 1 respectively for the order of approximation N ¼ 5; 7; 9 at t ¼ 1. The tables clearly confirm that the order of convergence for our proposed method increases as the degree of Laguerre polynomials in x and y increases. Again absolute error decreases with the increase of the order of approximation of the polynomial N, which clearly shows the effectiveness of our proposed numerical method. Again the data of CPU time shows that it takes minimum time to obtain accurate result showing that the method is computationally effective to solve the two-dimensional fractional order PDEs.
The absolute error is graphycally shown through Fig. 2 for various values of y and t for a fixed value of x ¼ 0:5. Figure 2 shows that the numerical solution is approximately accurate to the exact solution even for small value N ¼ 5, the better convergence can be achieved increasing N. Thus comparing the numerical solution with the exact solution, we can say that the proposed method is very much efficient.
The absolute error between the numerical solution and the exact solution is depicted through Fig. 3 at t ¼ 1 for various values of x and y. Figure 3 shows that the numerical solution is approximately accurate to the exact solution even for N ¼ 5. The better convergence can be achieved with increase in N.
Example 3 Consider following two-dimensional non-linear non-homogeneous fractional order advection-reaction-diffusion equation as where f ðx; y; tÞ ¼ 1:91116t 1:1 x 2 y 2 þ t 4 x 4 y 4 þ t 2 ðÀ2xy 2 À 2y 2 þ x 2 ðÀ2 À y 2 ÞÞ À 2x 2 yt 2 with the following initial and boundary conditions The exact solution for this concerned example is uðx; y; tÞ ¼ x 2 y 2 t 2 for fractional order time parameter a ¼ 0:9. The absolute errors, order of convergence and the CPU time required for the PDE with our proposed method are given in Tables 5 and 6 for y ¼ 0:5; 0 x 1 and x ¼ 0:5; 0 y 1 respectively for the order of approximation N ¼ 5; 7; 9 at t ¼ 1. It is clear from the tables that our proposed scheme is computationally effective even for less CPU time.
The absolute error between the numerical solution and the exact solution is depicted through Fig. 4 at t ¼ 1 for various values of x and y. Figure 4 shows that the numerical solution is approximately accurate to the exact solution even for N ¼ 5. The better convergence can be achieved with increase in N.

Results and discussion for proposed model
After being a justification of the accuracy and effectiveness of the method, the authors have been motivated to apply their proposed numerical scheme to solve the concerned two-dimensional nonlinear time fractional order reaction-diffusion model (4) under the prescribed initial and boundary conditions (5)- (9). The variations of the solute  The effects of reaction term on the solution profile for various values of the fractional order time derivative are found for k ¼ À1 and k ¼ 0; which are shown through Figs. 5 and 6 respectively. It is seen that the overshoots of sub-diffusion are decreased as the system advances towards fractional order from standard order. It is also seen from the figures that dampings are found due the presence of sink term (k ¼ À1). Figure 7 reveals that it consumes less time to stabilize the probability density function u(x, y, t) at t ¼ 0:5 for the case a ¼ 0:7 due to the effect of sink term k ¼ À1 as compared to the case of k ¼ 0 and also the presence of source term k ¼ 1, which clearly shows that physical importance of the presence of sink term in the system to enhance the stability region.   Four objectives have been accomplished through the present scientific contribution. First one is finding the solute concentration u(x, y, t) for two-dimensional non-linear fractional order reaction-diffusion equation employing the powerful and efficient technique Laguerre operational matrix. Second one is the pictorial presentations of the nature of overshoots in the form of sub-diffusion due to presence of reaction terms for different particular cases. The third one is the graphical presentations of the damping nature of the solute concentration when the system advances towards fractional order from standard order in the presence of sink term. The most important part of the present study is the exhibition of the order of convergence during validation of the proposed method while applying on the existing problems having analytical solutions. The FDE model has already been applied to describe heat transfer in ceramic polycrystalline composite materials made of two phases [49,51,38], in metal matrix composites [42,50,10,54,55] or in ceramic porous materials [52,53]. The authors are optimist that the numerical solution by the proposed algorithm can be extended to a wide range i.e., for solving 3D nonlinear time-space fractional order diffusion equations which will be the topic of our future research.