Analytical Solutions of a Space-Time Fractional Derivative of Groundwater Flow Equation

Abstract

The classical Darcy law is generalized by regarding the water flow as a function of a noninteger order derivative of the piezometric head. This generalized law and the law of conservation of mass are then used to derive a new equation for groundwater flow. Two methods including Frobenius and Adomian decomposition method are used to obtain an asymptotic analytical solution to the generalized groundwater flow equation. The solution obtained via Frobenius method is valid in the vicinity of the borehole. This solution is in perfect agreement with the data observed from the pumping test performed by the institute for groundwater study on one of their boreholes settled on the test site of the University of the Free State. The test consisted of the pumping of the borehole at the constant discharge rate Q and monitoring the piezometric head for 350 minutes. Numerical solutions obtained via Adomian method are compared with the Barker generalized radial flow model for which a fractal dimension for the flow is assumed. Proposition for uncertainties in groundwater studies was given.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Academic Editor:Hossein Jafari

Institute for Groundwater Studies, Faculty of Natural and Agricultural Sciences, University of the Free State, P.O. Box 399, Bloemfontein, South Africa

Received 5 September 2013; Accepted 12 September 2013; 21 January 2014

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The real problem encounter in groundwater studies up to now is the real shape of the geological formation in which water flows in the aquifer under investigation. However, there are many fractured rock aquifers where the flow of groundwater does not fit conventional geometries [1], for example, in South Africa, the Karoo aquifers, characterized by the presence of a very few bedding parallel fractures that serve as the main conduits of water in the aquifers [2]. With a challenge to fit the solution of the groundwater flow equation with experimental data from field observation in particular, the observed drawdown see [3], for all time yields a fit that undervalues the observed drawdown at early times and overvalues it at later times. The variation of observations from theoretically predictable values is usually an indication of uncertainties in the predictable. To investigate the first possibility Botha et al. [2] developed a three-dimensional model for the Karoo aquifer on the campus of the University of the Free State. This model is based on the conventional, saturated groundwater flow equation for density-independent flow: [figure omitted; refer to PDF] where _{S 0} is the specific storativity, K _ the hydraulic conductivity tensor of the aquifer, Φ ( x , t ) the piezometric head, f ( x , t ) the strength of any sources or sinks, with x _ and t the usual spatial and time coordinates; ∇ the gradient operator, and _{∂ t} the time derivative.

This model showed that the dominant flow field in these aquifers is vertical and linear and not horizontal and radial as commonly assumed. However, more recent investigations [3] suggest that the flow is also influenced by the geometry of the bedding parallel fractures, a feature that (1) cannot account for. It is therefore possible that the equation may not be applicable to flow in these fractured aquifers.

In an attempt to circumvent this problem, Barker [4] introduced a model in which the geometry of the aquifer is regarded as a fractal. Although this model has been applied with reasonable success in the analysis of hydraulic tests from boreholes in Karoo aquifers, it introduces parameters for which no sound definition exists in the case of noninteger flow dimensions.

As a review of the derivation of (1) will show, [5] Darcy law [figure omitted; refer to PDF] is used as a keystone in the derivation of (1). This law, proposed by Darcy early in the 19th century, is relying on experimental results obtained from the flow of water through a one-dimensional sand column, the geometry of which differs completely from that of a fracture [6]. There is therefore a possibility that the Darcy law may not be valid for flow in fractured rock formations but is only a very crude idealization of the reality [6]. Nevertheless, the relative success achieved by Botha et al. [2], to describe many of the properties of Karoo aquifers, suggests also that the basic principle underlying this law may be correct: the observed drawdown is to be related to either a variation in the hydraulic conductivity of the aquifer or a change in the piezometric head. Any new form of the law should therefore be reduced to the classical form under the more common conditions. Because K _ is essentially determined by the permeability of the rocks, and not the flow pattern, the gradient term in (2) is the most likely cause for the deviation between the observed and theoretical drawdown observed in the Karoo formations [6]. In the same direction, Cloot and Botha introduced the concept of non-integer fractional derivative to investigate a radially symmetric form of (1), by replacing the classical first order derivative of the piezometric head by a complementary fractional derivative [6]. However the generalized model for groundwater flow equation was solved numerically. In this work, more general form of groundwater flow equation will be introduced; the Frobenius and Adomian decomposition will be used to give an asymptotic solution of the generalized model for groundwater flow equation. Because the concepts of fractional (or non-integer) order derivatives and complementary fractional order derivatives may not be widely known, both concepts are first briefly discussed below.

2. Fractional Order Derivatives

On one hand, the concept of fractional calculus is popularly believed to have stemmed from a question raised in the year 1695 by Marquis de L Hospital (1661-1704) to Gottfried Wilhelm Leibniz (1646-1716), which sought the meaning of Leibniz's currently popular notation ^{d n} y / d ^{x n} for derivative of order n ∈ _{... 0} [: =] { 0 , 1 , 2 , ... } when n = 1 / 2 (what if n = 1 / 2 ). In his reply, dated September 30, 1695, Leibniz wrote to L'Hospital as follows: "This is an apparent paradox from which, one day, useful consequences will drawn ... ." On the other hand, the concept of fractional order derivatives for a function, say f ( x ) , is based on a generalization of the Abel integral [7]: [figure omitted; refer to PDF] where n is a nonzero positive integer [8-12] and Γ ( n ) the Gamma function [13]. This represents an integral of order n for the continuous function f ( x ) , whenever f and all its derivatives vanish at the origin, x = 0 . This result can be extended to the concept of an integral of arbitrary order c , defined as: [figure omitted; refer to PDF] where c is a positive real number and j an integer such that 0 < s ...4; 1 .

Let ρ now be the least positive integer larger than α such that α = m - ρ ; 0 < ρ ...4; 1 . Equation (4) can then be used to define the derivative of (positive) fractional order, say α , of a function f ( x ) as [figure omitted; refer to PDF] Note that these results, like Abel's integral, are only valid subject to the condition that ^{f ( k )} ( x ) |" x = 0 = 0 for k = 0,1 , 2 , ... , p [8-11].

2.1. Properties

Properties of the operator can be found in [14, 15], we mention only the following: for f ∈ _{C μ} , μ ...5; - 1 , α , β ...5; 0 and γ > - 1 : [figure omitted; refer to PDF]

3. A Generalized Mathematical Groundwater Flow Model

For the sake of clarity the generalization of the classical model for groundwater flow in the case of density-independent flow in the uniform homogeneous aquifer is considered in this paper [4]. Consider the following groundwater flow equation [figure omitted; refer to PDF] where both the specific storability, _{S 0} , and hydraulic conductivity, K , are scalar and constant quantities and n = 1,2 or 3 is the radial dimension. To be complete, the following set of initial and boundary conditions is added: [figure omitted; refer to PDF] Here Φ ( r , 0 ) = a r means that before a pumping test begins the level of water or the initial hydraulic head in the aquifer is linear function of space with a positive gradient a to be found, _{lim ... r [arrow right] ∞} ... Φ ( r , t ) = a r means that during the pumping test the level of water is not affected for a very long distance from the borehole, and Q = ( 2 ^{π n / 2} / Γ ( n / 2 ) ) ^{r b n - 1} K ^{d 3 - n} _{∂ r} Φ ( _{r b} , t ) means that the rate of pumping is proportional to the hydraulic conductivity.

Hence Q is the discharge rate of the borehole, with radius _{r b} and d the thickness of the aquifer from which the borehole taps. In order to include explicitly the possible effect of the geometry into mathematical model the radial component of the gradient of the piezometric head, _{∂ r} Φ ( r , t ) is replaced by the Weyl-fractional derivatives of order α = m - ρ ; the fractional derivative in this paper is Caputo derivative and is defined as follows: [figure omitted; refer to PDF] This provides a generalized form of the classical equation governing the flow equation (1): [figure omitted; refer to PDF] The integrodifferential equation does contain the additional parameter α = m - ρ , 0 < ρ ...4; 1 , which can be viewed as a new physical parameter that characterizes the flow through the geological formations [6]. The same transformation generates also a more general form for the boundary condition at the borehole [6]: [figure omitted; refer to PDF] The relations (10) and (11), together with the initial condition described in (10), represent a complete set of equations for which a solution exists. The integro-differential character of the relations makes the search for analytical solution very difficult however. Nevertheless, in this paper we make use of Frobenius and Adomian decomposition methods to give an asymptotic solution.

4. Solution of the Generalized Flow Equation

4.1. Frobenius Method

In this work to perform the Frobenius method, we consider the groundwater flow governed by the following fractional Caputo-Weyl derivation partial differential of order 2 α , where α is real number, 0 < α ...4; 1 . Also, we consider the dimension of the flow to be 2. Therefore, for n = 2 (10) can then be transformed into the following equation: [figure omitted; refer to PDF] Applying the Laplace operator on both sides of the above equation, we have the following ordinary differential equation [figure omitted; refer to PDF] where s is the variable of Laplace. In this matter we choose Φ ( r , 0 ) = 0 or we choose a = 0 , meaning that the level of water is the same everywhere in the aquifer if the water is not taken out from the aquifer. If we let [figure omitted; refer to PDF] then (13) becomes [figure omitted; refer to PDF] We put p ( r ) = 1 / r and q ( r ) = - ( _{S 0} / K ) s . To meet the condition under which Frobenius method can be used, we have to prove that p ( r ) and q ( r ) are analytical around _{r b} that means we have to prove that p ( r ) and q ( r ) can be written as series. It is very obvious to see that p ( r ) and q ( r ) can be expressed as follows: [figure omitted; refer to PDF] We start here with the coefficients of p ( r ) : [figure omitted; refer to PDF] implying that 1 = ^{∑ n = 0 ∞} _{p n} ^{( r - _{r b} ) n α} r .

Putting R = r - _{r b} we have that 1 = ^{∑ n = 0 ∞} _{p n} ^{R n α} ( R + _{r b} ) and equating the coefficients of same power we obtain the following set of equations: [figure omitted; refer to PDF] Therefore, the coefficients can be given with the general following recursive formula: [figure omitted; refer to PDF] And obviously the coefficients _{q n} are given below as [figure omitted; refer to PDF] From the above expression we can see that p ( r ) and q ( r ) are analytical around _{r b} , which follows from Frobenius method that the solution of (15) can be in the form [figure omitted; refer to PDF] This solution is not convergent for a large r ; that is, the solution will diverge if we observe the drawdown at a position very far from the borehole from which the water is taken out. Therefore, we restrict our solution in the vicinity of the borehole; more precisely we investigate the solution when | r - _{r b} | < 1 . That means we pump water from the borehole and we observe the drawdown in the vicinity of the borehole.

Thus substituting (21), p ( r ) , and q ( r ) into (15) and equating the coefficients of the same power, we obtain the following recursive formula for which the coefficients _{Ψ n} ( n ...5; 0 ) , coefficients of our series: [figure omitted; refer to PDF] Here we have the following: [figure omitted; refer to PDF] However, making use of the boundary condition in (11), we have the following: [figure omitted; refer to PDF] Then the general solution of the fractional Caputo-Weyl derivation of groundwater flow of order 2 α in Laplace space is given below as: [figure omitted; refer to PDF] To observe the behavior of the solution at the borehole, which corresponds to r = _{r b} , the series solution is reduced here to the coefficient with order zero which is obtained when n = 1 / α - 1 and it's given below as [figure omitted; refer to PDF] In the following section the analytical asymptotic solution obtained in Laplace space via Frobenius method will be compared with the experimental data.

4.2. Numerical Results

In order to examine the validation of this solution, the above asymptotic solution is compared with the experimental data from the pumping test performed by the Institute for Groundwater Studies on one of their borehole settled on the campus test site of the University of the Free State. The test consisted of the pumping of the borehole at the constant discharge rate Q and monitoring the piezometric head for 350 minutes. The first step in the section is to discretize the range of Laplace transform since the exact Laplace transform cannot be obtained in practice. This is done as follows: [figure omitted; refer to PDF] For _{t i} ...4; t ...4; _{t i + 1} we approximate Φ ( r , t ) = ( Φ ( r , _{t i} ) + Φ ( _{t i + 1} ) ) / 2 where Φ ( r , _{t i} ) for 0 ...4; i ...4; n ; then the results obtained in the fields and Laplace transform become [figure omitted; refer to PDF] Using this numerical scheme, the physical data was transformed into Laplace space. A comparison between these values and asymptotic computed data can only be provided in the in real space not in Laplace space. Since it is not worth concluding the validity of this solution in Laplace space, the inverse Laplace transform is applied in (26). To test the validity of this solution in real space and applying the inverse Laplace transform on (26), (29) is obtained [figure omitted; refer to PDF]

The above solution is compared graphically to the experimental data from the pumping test performed by the institute for groundwater study on one of their borehole settled on the test site of the University of the Free State. The small difference observed in the above graph (Figure 1) is due to uncertainties in measurement and this will be discussed in Section 6 of this work. The aquifer parameters used in this models are recorded in Table 1, the observed data from field observation will be attached to this paper. Although this solution is in agreement with the experimental data, there will be the need to investigate the case where the observation can be done for a long distance. In the next section another approach will be introduced to solve the space-time-fractional derivative of groundwater flow equation, and this method is the Adomian decomposition method.

Table 1: Aquifer parameters used in the model.

Parameters	Values	Units
Q	83.1 0 - 3	(m3 s-1 )
r b	0.025	(m)
K	13.1 0 - 4	(m s-1 )
d	1	(m)
S 0	13.1 0 - 4	(m-1 )
Φ 0	0	(m)
2 α = 1.46

Figure 1: Comparison of experimental data from real world observation with solution of fractional groundwater flow equation.

[figure omitted; refer to PDF]

5. Adomian Decomposition Method

5.1. Example 1

This section is concerned with the groundwater equation with time- and space-fractional derivatives of the form [figure omitted; refer to PDF] Subject to the initial and boundary conditions described in (8), the level of water is assumed to be the same throughout the aquifer before the pumping so that the gradient a described in (8) is zero. Furthermore, it is assumed that a fractional change in drawdown is constant for t = 0 meaning ^{∂ r α} Φ ( r , 0 ) = constant .

The method used here is based on applying the operator ^{∫ 0 t} d τ on both sides of (30) to obtain [figure omitted; refer to PDF] The Adomian decomposition method [16, 17] assumes a series solution for (31) to be [figure omitted; refer to PDF] where the components _{Φ n} ( r , t ) are determined recursively. Substituting (32) into both sides of (30) gives [figure omitted; refer to PDF] Following the decomposition method, the recursive relations are introduced as [figure omitted; refer to PDF] It is worth noting that if the component _{Φ 0} ( r , t ) is defined, then the remaining components n ...5; 1 can be completely determined such that each term is determined by using the previous terms, and the series solutions are thus entirely determined. Finally, the solution Φ ( r , t ) is approximated by the truncated series [figure omitted; refer to PDF] However, the inclusion of boundary conditions in fractional differential equations introduces additional difficulties. The Adomian decomposition method can handle these difficulties by using the space-fractional operator ^{∂ r α} and the initial conditions only. The method provides the solution in the form of a rapidly convergent series that may lead to the exact solution in the case of integer derivatives ( α = 1 , μ = 2 ) and to an efficient numerical solution with high accuracy for 0 ...4; α < 1 . The convergence of the decomposition series has been investigated in [18, 19].

Following the recursive formula equation (35) and using the fact that ^{∂ r α} Φ ( r , 0 ) = constant = a , the equations below are obtained: [figure omitted; refer to PDF] The component _{Φ 4} ( r , t ) was also determined and will be used, but for brevity it is not listed. In this matter five components of the decomposition series (30) were obtained for which Φ ( r , t ) was evaluated to have the following expansion: [figure omitted; refer to PDF] Applying the boundary condition yields [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

5.2. Example 2

Consider the groundwater flow equation with time- and space-fractional derivatives of the form [figure omitted; refer to PDF] subject to the initial condition described in (11). Furthermore we suppose that the gradient a = 0 , meaning that the water level is everywhere in the aquifer at ^{∂ r μ} Φ ( r , 0 ) = constant and the fractional change in drawdown is a constant and boundary condition: [figure omitted; refer to PDF] Following the discussion earlier, we have the below recursive formula: [figure omitted; refer to PDF] It follows from the recursive formula that [figure omitted; refer to PDF] The component _{Φ 4} ( r , t ) was also determined and will be used, but for brevity it is not listed. In this matter five components of the decomposition series (41) were obtained of which Φ ( r , t ) was evaluated to have the following expansion: [figure omitted; refer to PDF] Applying the boundary condition yields [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

5.3. Example 3

Consider time fractional derivative of the groundwater flow equation with time-fractional derivatives of the form [figure omitted; refer to PDF] Subject to the initial and boundary conditions described in (8) it is assumed that a ...0; 0 . Following the discussion presented earlier, we obtained the below recursive formula: [figure omitted; refer to PDF] The component _{Φ 4} ( r , t ) was also determined and will be used, but for brevity not listed. In this matter five components of the decomposition series (31) were obtained of which Φ ( r , t ) was evaluated to have the following expansion [figure omitted; refer to PDF] Applying the boundary conditions yield to [figure omitted; refer to PDF] where [figure omitted; refer to PDF] The normalized solutions with Φ ( r , 0 ) of (50), (51) and (52) are illustrated graphically in Figures 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 13 for different values of various parameters; these solutions are compared to the solution proposed by Barker [4]. These graphs show the behaviour of the drawdown during the pumping test, first as a function of space and time from the borehole to the point of observation, secondly as function of time for a fixed distance from the borehole, and finally as a function of space for a fixed time. Based on the assumption that groundwater in an aquifer flows through an equipotential surface that are projections of l -dimensional spheres onto two-dimensional space, Barker [4] obtained the following analytical solution for an infinite aquifer with a line source: [figure omitted; refer to PDF] where Γ is the incomplete gamma function and λ the dimension of the flow which equals the special dimension, n being an integer, and the other quantities all have the same meaning as before. Although this model has been applied with reasonable success in the analysis of hydraulic tests from boreholes in Karoo aquifers, it introduces parameters for which no sound definition exists in the case of non-integer flow dimensions. The main raison of comparison of these results obtained via the Adomian decomposition methods with those of Barker's fractal radial flow model is to establish a possible relationship between the fractional order of the derivative and the parameter fractal introduced earlier by Barker.

Figure 2: Solution of (30) (piezometric head).

[figure omitted; refer to PDF]

Figure 3: Solution of (41) (piezometric head).

[figure omitted; refer to PDF]

Figure 4: Barker's solution.

[figure omitted; refer to PDF]

Figure 5: Solution of (48) (piezometric head).

[figure omitted; refer to PDF]

Figure 6: Solution of (30) (drawdown).

[figure omitted; refer to PDF]

Figure 7: Solution of (41) (drawdown).

[figure omitted; refer to PDF]

Figure 8: Solution of (48) (drawdown).

[figure omitted; refer to PDF]

Figure 9: Barker's solution.

[figure omitted; refer to PDF]

Figure 10: Solution of (30) as function of space (piezometric head).

[figure omitted; refer to PDF]

Figure 11: Solution of (41) as function of space (piezometric head).

[figure omitted; refer to PDF]

Figure 12: Barker's solution.

[figure omitted; refer to PDF]

Figure 13: Solution of (48) as function of space (piezometric head).

[figure omitted; refer to PDF]

Table 2 shows the theoretical values of the discharge rate and aquifer parameters used in the numerical simulations.

Table 2: Theoretical values.

Parameters	Values	Units	Parameters	Values	Units
Q	10.8	(m3 s-1 )	α	0.005	None
r b	0.15	(m)	β	1.05	None
K	0.65	(m s-1 )	μ	1.005	None
d	40	(m)	λ	0.012	None
S 0	0.75	(m-1 )	α 1	0.5	None
Φ 0	0	(m)	α 2	0.005	None
			μ 1	0.025	None

6. Discussion and Propositions

Although the analytical solution obtained via Frobenius method fit the experimental data or have described successfully the events taking place in the vicinity of the borehole, on one hand, and the analytical solutions obtained via Adomian decomposition method was successfully compared to the solution proposed by Barker, on the other hand, the problem of choosing an appropriate geometry for the geological system in which the flow occurs still remains a challenge in groundwater studies. We personally believe that we do not propose solution to a problem because it going to be useful for the generation in which we are, but we propose solution because we hope that they will one day be useful for mankind; therefore the following proposition may not be useful for this generation because we may not have the adequate technology to perform the steps involved but it will be in the future. We think that describing the groundwater flow with one equation for the whole aquifer is unrealistic, because from one point of the aquifer to another properties including geology and geometry change, therefore the flow. Assumption such as homogeneity, isotropy, uniform thickness over the area under investigation and so on render the study of groundwater uncertain. In order to study the geology or the geometry of an aquifer, we have to divide the aquifer in small portions, from north to south and from east to west. The geology and the geometry of each portion can then be studied. If the results of the study reveal that the portion under investigation is for instance a hyperboloid, since there is no exact solution to groundwater flow model by hyperboloid flow and there are solutions to circular flow, then a suitable transformation can be done including transformation of hyperboloid coordinates to Cartesian coordinates, then to Cartesian coordinates to cylindrical coordinate such that this new coordinate can now be put in the equation describing groundwater flow, and the solution to the new equation can then be investigated. Henceforth knowing the real geology and geometry of this portion, the real paths flow will be known. Then having a good knowledge of each small portion including its geometry and geology, the real geometry and geology of the aquifer can be not exactly but more accurately known.

For groundwater remediation it will be possible to know where the maximum, minimum, and average chemical concentrations are found in the aquifer. We have no proof for this but we believe that this proposition will be useful in reducing uncertainties in groundwater study. It is believed that the field test gives the characteristic of an aquifer, but we believe that the field test gives both uncertainties and characteristics of the aquifer; therefore quantify uncertainties in this measurement lead us to the real picture of aquifer characteristics, henceforth we propose that the studies in groundwater should focus on both uncertainties and fields observations, because what is known is bounded by what is not known; knowing what is not known give a real picture to what was known, and it follows that the knowledge of uncertainties in groundwater study will give a clear picture of what we already know in groundwater.

7. Conclusion

The classical Darcy Law has been generalized using the concept of complementary fractional order derivatives of Weyl fractional derivative. This leads to the formulation of a new generalized form of the groundwater flow equation [6]. The applications of Adomian decomposition and Frobenius methods were extended to obtain explicit and numerical solutions of the space-time fractional groundwater flow. The two methods were very clearly efficient and powerful techniques in finding the solutions of the proposed equations. The solution obtained via Frobenius takes into account the events taking place in the vicinity of the borehole during the pumping test whereas the solution obtained via Adomian decomposition methods takes into account the events that take place far from the point where water is pumped out, that is, in the borehole. The solution obtained via Frobenius method was in perfect agreement with the observed data obtained from the pumping test performed by the Institute for Groundwater Studies on one of their borehole settled on the campus test site of the University of the Free State. The Adomian decomposition method requires less computational work than existing approaches while supplying quantitatively reliable results. The obtained results demonstrate the reliability of the algorithms and their wider applicability to fractional evolution equations. A comparison of these results obtained via the Adomian decomposition methods with those of Barker's fractal radial flow model suggests that there exists a relation between the fractional order of the derivative and the non-integral dimension of the flow.

Conflict of Interests

The authors declare that there is no conflict of interests for this paper.

References

[1] J. H. Black, J. A. Barber, D. J. Noy, "Crosshole investigations: the method, theory and analysis of crosshole sinusoidal pressure tests in fissured rocks," Stripa Projects Internal Reports , no. 86-03, SKB, Stockholm, Sweden

[2] J. F. Botha, I. Verwey van Voort, J. J. P. Viviers, W. P. Collinston, J. C. Loock, "Karoo aquifers. Their geol­ogy, geometry and physical behaviour," WRC Report , no. 487/1/98, Water Research Commission, Pretoria, 1998.

[3] G. J. van Tonder, J. F. Botha, W.-H. Chiang, H. Kunstmann, Y. Xu, "Estimation of the sustainable yields of boreholes in fractured rock formations," Journal of Hydrology , vol. 241, no. 1-2, pp. 70-90, 2001.

[4] J. A. Barker, "A generalized radial flow model for hydraulic tests in fractured rock," Water Resources Research , vol. 24, no. 10, pp. 1796-1804, 1988.

[5] J. Bear Dynamics of Fluids in Porous Media , of American Else­vier Environmental Science Series, Elsevier, New York, NY, USA, 1972.

[6] A. Cloot, J. F. Botha, "A generalised groundwater flow equation using the concept of non-integer order derivatives," Water SA , vol. 32, no. 1, pp. 1-7, 2006.

[7] R. Courant, F. John Introduction to Calculus and Analysis , vol. 2, John Wiley & Sons, New York, NY, USA, 1974.

[8] Y. Cherruault, G. Adomian, "Decomposition methods: a new proof of convergence," Mathematical and Computer Modelling , vol. 18, no. 12, pp. 103-106, 1993.

[9] I. Podlubny Fractional Differential Equations , Academic Press, San Diego, Calif, USA, 1999.

[10] I. Podlubny, "Geometric and physical interpretation of fractional integration and fractional differentiation," Fractional Calculus & Applied Analysis , vol. 5, no. 4, pp. 367-386, 2002.

[11] K. Adolfsson, "Nonlinear fractional order viscoelasticity at large strains," Nonlinear Dynamics , vol. 38, no. 1-4, pp. 233-246, 2004.

[12] O. P. Agrawal, "Application of fractional derivatives in thermal analysis of disk brakes," Nonlinear Dynamics , vol. 38, pp. 191-206, 2004.

[13] G. Afken Mathematical Methods for Physicists , Academic Press, London, UK, 1985.

[14] K. S. Miller, B. Ross An Introduction to the Fractional Calculus and Fractional Differential Equations , John Wiley & Sons, New York, NY, USA, 1993.

[15] K. B. Oldham, J. Spanier The Fractional Calculus , Academic Press, New York, NY, USA, 1974.

[16] G. Adomian, "A review of the decomposition method in applied mathematics," Journal of Mathematical Analysis and Applications , vol. 135, no. 2, pp. 501-544, 1988.

[17] A. Atangana, "New class of boundary value problems," Information Science Letters , vol. 1, no. 2, pp. 67-76, 2012.

[18] G. Adomian Solving Frontier Problems of Physics: The Decomposition Method , Kluwer Academic, Dodrecht, The Netherlands, 1994.

[19] Y. Cherruault, "Convergence of Adomian's method," Kybernetes , vol. 18, no. 2, pp. 31-38, 1989.