## Abstract

As fractional diffusion equations can describe the early breakthrough and the heavy-tail decay features observed in anomalous transport of contaminants in groundwater and porous soil, they have been commonly used in the related mathematical descriptions. These models usually involve long-time-range computation, which is a critical obstacle for their application; improvement of computational efficiency is of great significance. In this paper, a semi-discrete method is presented for solving a class of time-fractional diffusion equations that overcome the critical long-time-range computation problem. In the procedure, the spatial domain is discretized by the finite element method, which reduces the fractional diffusion equations to approximate fractional relaxation equations. As analytical solutions exist for the latter equations, the burden arising from long-time-range computation can effectively be minimized. To illustrate its efficiency and simplicity, four examples are presented. In addition, the method is used to solve the time-fractional advection–diffusion equation characterizing the bromide transport process in a fractured granite aquifer. The prediction closely agrees with the experimental data, and the heavy-tail decay of the anomalous transport process is well represented.

## 1. Introduction

For the contaminant transport processes in soil and groundwater, diffusion equations (such as the diffusion, advection–dispersion and advection–reaction–diffusion equations) are the traditional governing equations [1–4]. In the past several decades, however, a lot of evidence has shown that some of the critical features in contaminant transport through complex porous media cannot be described by the conventional diffusion equations [5–10]. These features include early breakthrough and heavy-tail decay of the contaminant, as well as the scale-dependent coefficients [11–13]. They have led to the increasing use of fractional diffusion equations for modelling contaminant transport in porous media.

The theoretical research on fractional diffusion equation models has achieved considerable success in physical modelling and experimental result analysis in the past decades [14–17]. Nevertheless, numerical methods for solving fractional diffusion equations are still immature for practical applications, in which the spatial problem domains are geometrically complex and large, whereas long-time-range predictions are often desired. The main reason is that the fractional derivative has a global nature, compared with the locally expressed classic derivatives. The computational cost of the time-fractional derivative term will increase dramatically with time evolution. This brings the computational challenge of approximating fractional order equations with the finite difference or the finite element methods (FEMs) [18–23]. Although the short-memory method has been proposed to tackle the long-time-range computation, it has been shown to bring accuracy deterioration in many cases [24,25]. Moreover, when computing high-Péclet-number problems, numerical schemes may produce oscillating solutions [26–30]. The problem is more annoying in fractional advection–dispersion equations.

Until now, the finite difference method has been widely used for solving anomalous diffusion equations with successes in short-time-range and small-spatial-scale problems [31,32]. Because the FEM is more suitable for modelling large and geometrically complicated spatial domains, the investigation of the FEM for fractional diffusion-type equations has attracted much attention in recent years. Fix [33] and Roop *et al.* [22] approximated the solutions of steady-state space-fractional advection–dispersion equations in two spatial dimensions using the Galerkin and least-squares FEMs. Huang *et al.* [7] proposed an unconditionally stable FEM approach to solve the one-dimensional space-fractional advection–dispersion equation, and successfully applied it to simulate atrazine transport in a saturated soil column. Deng [34] developed an FEM for the numerical resolution of the space- and time-fractional Fokker–Planck equation, with convergence order *O*(*k*^{2−α}+*h*^{μ}), where *α* and *μ* are time- and spatial-derivative orders. Zhuang *et al.* [35] investigated the Galerkin finite element approximation of symmetric space-fractional partial differential equations, and proved the stability and convergence of the proposed schemes. Zheng *et al.* [36] and Li *et al.* [37] discussed the FEM for the space-fractional advection–diffusion equation with non-homogeneous initial and boundary conditions. Yang *et al.* [38] used the transfer matrix technique to reduce the computation cost for solving the time–space fractional diffusion equation in two dimensions. Although all the aforementioned works indicate that FEMs play an important and increasing role in the applications of fractional diffusion-type equation models, efficient and simple numerical methods for fractional diffusion equations are still urgently needed.

In this paper, we introduce a semi-discrete FEM for time-fractional diffusion equations that can be expressed in the form
1.1in which *u* is the scalar unknown, *t* denotes time, *γ* is the fractional derivative order, ∇ is the gradient operator, **A** is a coefficient vector, **D** is a coefficient matrix, and *P* and *f* are scalars. Moreover, **A**, **D**, *P* and *f* are functions of the spatial variables. In the proposed semi-discrete FEM, the FEM is used for spatial discretization that reduces the time-fractional diffusion equations to approximate fractional relaxation equations (also known as temporal fractional ordinary differential equations). The FEM can conveniently be used to discretize large and geometrically complicated spatial domains. On the other hand, analytical solutions for fractional relaxation equations exist, and the burden of long-time-range computation can be significantly alleviated. The present semi-discrete method can solve one-, two- and three-dimensional problems with variable coefficients conveniently at low implementation cost.

## 2. Algorithm framework

Time-fractional diffusion equations are often used to characterize the contaminant transport processes, which exhibit typical subdiffusion features, such as heavy-tail decay of concentration and nonlinear, time-dependent mean-square displacement 〈*x*^{2}(*t*)〉∼*t*^{α} (0<*α*<1) [3]. The time-fractional diffusion problem can be expressed as
2.1Here, *u* is contaminant concentration, **A** is the generalized convective coefficient vector, **D** is the generalized diffusion coefficient matrix, *Pu* represents a reaction–absorption term, *f* is the source term and *Ω* denotes the spatial domain of the problem. Moreover, **n** is the unit outwards normal vector to the boundary, *h*_{c} is the convective coefficient, and *Γ*_{D}∪*Γ*_{N}∪*Γ*_{C}=∂*Ω* denotes the entire boundary of *Ω*. The subscripts **D**, **N** and **C** for *Γ* designate the essential (Dirichlet), natural (Neumann) and convective boundary conditions, respectively, where *Γ*_{D}, *Γ*_{N} and *Γ*_{C} are mutually exclusive. It will be assumed that **A**, **D**, *P*, *f*, *q* and are independent of time. In the problem statement, d^{γ}/*dt*^{γ} represents the Caputo fractional derivative whose definition is given as
2.2in which *Γ* is the Gamma function and *γ* is the derivative order. The Caputo fractional derivative has the following properties [39,40]:
2.3Here *E*_{γ} represents the Mittag–Leffler function with one parameter [25]:
2.4The weighted residual statements for the governing equation, natural boundary condition and convective boundary condition
2.5can be merged to form the following weak form with the help of the divergence theorem:
2.6Here the trial solution *u* equals and the weight function *ψ* vanishes on *Γ*_{D}. In the FEM, *u* and *Ψ* within each element *Ω*^{e} can be expressed, respectively, as
2.7Here and are the values of *u* at nodes away from and on *Γ*_{D}, respectively; are the value of *ψ* at nodes away from *Γ*_{D}; and are the nodal interpolation functions for and , respectively; *n* is the number of nodes in the element, and *n*_{D} is the number of nodes on *Γ*^{e}_{D}=∂*Ω*^{e}∩*Γ*_{D}; and form the row interpolation matrices **N**^{e} and , respectively; and form the vectors **U**^{e} and , respectively. By virtue of (2.7), (2.6) can be written as
2.8or
2.9in which
2.10In **F**^{e}, *Γ*^{e}_{N}=∂*Ω*^{e}∩*Γ*_{N} and *Γ*^{e}_{C}=∂*Ω*^{e}∩*Γ*_{C}. Moreover, *Ψ*, **C**, **K**, **U**, and **F** are the assembled counterparts of *Ψ*^{e}, **C**^{e}, **K**^{e}, **U**^{e}, and **F**^{e}, respectively. The arbitrary nature of *Ψ* leads to the following system equation:
2.11If
2.12the above system equation can be transformed into
2.13where
2.14

At last, the time-fractional system (2.1) leads to
2.15An exact solution of the fractional relaxation equation in (2.15) exists and can be expressed as [41,42]:
2.16where **M**=**C**^{−1}**K**, and *E*_{γ} is the Mittag–Leffler function that has been accurately evaluated by Chen [43] and Podlubny [44]. In our computations, the used value of the function will be accurate up to 10^{−12}.

It can be deduced that, if *γ*=1.0, (2.16) becomes the exponential solution of the integer-order relaxation equation. Next, we decompose the Mittag–Leffler function in (2.16) as
2.17where **B** is the modal matrix formed by the eigenvectors of −**M**. On the other hand, *Λ*_{t} is a diagonal matrix whose *i*th diagonal entries are *E*_{γ}(−*Λ*_{i}*t*^{γ}) and *Λ*_{i} is the *i*th eigenvalue of −**M**. Substituting (2.17) into (2.16), we get
2.18

Through the above manipulations, the initial–boundary value problem in (2.1) is reduced to an initial value problem through spatial finite element discretization. This reduced problem can be solved analytically in terms of the Mittag–Leffler function. This practice drastically lowers the cost associated with the long-time-range computation of the initial–boundary value problem.

For the essential boundary condition in (2.12), the equation can be transformed as 2.19which would require to satisfy a fractional relaxation equation. In particular, the typical case in which is a constant, , also satisfies (2.12).

## 3. Numerical examples

### (a) One-dimensional time-fractional diffusion equation

The following simple time-fractional diffusion problem is considered:
3.1If the diffusion coefficient *k*=*L*^{2}/*π*^{2}, the exact solution is whose value at *x*=*L*/2 is shown in figure 1 for *γ*=0.4,0.7 and 1.0.

To construct the spatial discretization, both linear and quadratic elements are attempted. For the linear element, the shape functions are
3.2and the element matrices are
3.3in which *h* is the nodal spacing. For the quadratic element, the shape functions are
3.4and the corresponding element matrices are
3.5

Table 1 lists the normalized errors at different time instants yielded by using 10 linear, 10 quadratic and 100 linear elements. The proposed method can achieve accurate results, no matter whether linear or quadratic elements are used. As usual, the quadratic element delivers much more accurate results than the linear element at the same nodal spacing. Another important feature of this method is that the accuracy of the numerical result at large time constants can be improved by reducing the nodal spacing *h*.

To estimate the convergence ratio of the linear element and the quadratic element, the results in table 2 evaluated at *t*=10 but different nodal spacings are prepared, and the error is
3.6It can be seen that the convergence ratio of the linear element is *O*(*h*^{2}), and the quadratic element is *O*(*h*^{4}).

In table 1, the normalized errors increase with *t*. To investigate the efficiency in tackling long-time-range diffusion problems, the normalized errors of the linear and the quadratic elements at large *t* are computed and listed in table 3. It can be seen that the normalized errors remain fairly steady with respect to *t*.

### (b) One-dimensional time-fractional convection–dispersion equation

The time-fractional advection–dispersion equation (e.g. time-fractional Fokker–Planck equation), which exhibits a heavy-tail concentration decay feature, is usually used to characterize contaminant transport in natural porous media. A simple example is
3.7Assuming *a* and *k* are constants, *a*>*k*, the exact solution of equation (3.7) can be written as
3.8which is portrayed in figure 2 for *t*<1. For the linear element, the corresponding element matrices are
3.9

The normalized errors obtained by using 10, 20 and 40 elements at *x*=*L*/2 and different *t* are shown in table 4. With only 10 elements, the errors are less than 0.1 per cent.

### (c) Two-dimensional time-fractional diffusion equation

In this example, the following two-dimensional problem is considered:
3.10in which *k* is the diffusion coefficient, and *Ω*=[0,*L*]×[0,*L*]. If *k*=1/*π*^{2} and *L*=1.0, the exact solution of the problem is .

The square problem domain is modelled by 4×4, 8×8 and 16×16 four-node square elements. The element interpolation functions are 3.11and 3.12

Following the calculation steps (2.15)–(2.18), figure 3 plots the numerical solution for *γ*=0.8 at *t*=2. The normalized errors for *γ*=0.8 at *x*=*L*/2, *y*=*L*/2 and various values of *t* are given in table 5. The errors drop with the nodal spacings. Indeed, the FEM can readily take coordinate-dependent and direction-dependent diffusion coefficients into account.

### (d) Time-fractional diffusion equation in a quarter of a circular domain

An important advantage of the FEM over the finite difference method is that the former can readily consider complex spatial domains. In this example, the following problem defined over a circular domain is considered:
3.13in which *Ω*={(*x*,*y*)|*x*>0, *y*>0, *x*^{2}+*y*^{2}<1}. If *k*=1, the exact solution of (3.13) is , in which *J*_{0} represents the zeroth-order Bessel function of the first kind. For symmetry, we need to model only a quarter of the problem domain, and a typical mesh is depicted in figure 4.

For the four-node element, the (*x*,*y*) coordinates are also interpolated with the functions given in (3.11), i.e.
3.14in which (, ) are the coordinates of the *i*th element nodes, and
3.15In the above expressions,
3.16

The matrices **C**^{e} and **K**^{e} are computed by the second-order Gauss–Legendre rule. In our computations, a quarter circle is partitioned into three, 48 and 217 elements. The corresponding numerical results at some selected spatial locations are listed in table 6 for *γ*=0.8. Table 6 indicates that the proposed method is capable of delivering accurate solution to the anomalous diffusion problem (3.13) and the accuracy can be improved by using more elements in modelling the computational domain.

## 4. Application

To investigate the efficiency and applicability of the proposed semi-discrete method in solving real-world problems, it is used to solve the problem of tracer solute transport in an aquifer. The experiment was conducted using a test aquifer in Nevada and the schematic is shown in figure 5 [5]. Bromide of quantity *M*=20.81 kg was used as a non-sorbing tracer solute and introduced into the injection well for a period of *T*_{0}=3.54 days at an average concentration of 3.77 kg m^{−3}. The reference point, the injection well and the extraction well are located at *r*=0, 30 m (*R*_{i}) and 60 m (*R*_{e}) along the downstream direction of the underground water flow. The radius of the extraction well is 0.127 m, the centre of the extraction well is at *r*_{c}=60.127 m, and the pumping rate is *Q*=12.4 m^{3} d^{−1}. The solute concentration in the extraction well was monitored for about 321 days, and the screened interval was *b*=35 m. A more detailed description of the experiment can be found in earlier studies [5,46–48].

Because *T*_{0} is short compared with the total time range (approx. 321 days) of measurement, the following radial initial–boundary value problem for the solute transport in the fractured granite aquifer is established in terms of the solute concentration *u* as
4.1where *υ*_{0} is the convective coefficient, *d*_{0} is the dispersion coefficient, *υ*_{0}=*ad*_{0} and *a* is the dispersivity. Moreover, *υ*_{0}/(*r*_{c}−*r*) and *d*_{0}/(*r*_{c}−*r*) have the units of [L/T^{γ}] and [L^{2}/T^{γ}], which represent the non-local aquifer properties. The initial value is normalized as *f*(*r*)=*Mδ*(*r*−(*r*_{c}−*R*_{i}))/(2*π*(*r*_{c}−*R*_{i})*bθT*_{0}), where *θ* is the hydraulic parameter. The boundary conditions imply that the solute cannot reach *r*=0 by upstream dispersion, and the solute moves by advection at *r*=*R*_{e}, which gives the wall of the extraction well [5].

In order to obtain a highly accurate numerical approximation, the quadratic element is adopted. The element matrices **C**^{e} and **K**^{e} are computed by the second-order Gauss–Legendre rule.

A comparison of the numerical predictions and the experimental data is shown in figure 6, and the heavy-tail feature characterized by the time-fractional model (4.1) with different derivative orders is shown in figure 7. It can be observed from figure 6 that the numerical result offers a good fit to most of the experimental data. Owing to the subdiffusion behaviour in the aquifer matrix and immobile water, in the experimental result, the concentration of bromide exhibits a rather slow decay at late times. Figure 7 confirms that the time-fractional radial flow model (4.1) captures the long-time behaviour with heavy tail. Figure 7 also illustrates that the heavy-tail feature becomes more remarkable with decreasing time-fractional derivative order *γ*. Hence, in this model (4.1), the time-fractional derivative order *γ* is an indicator of the non-Fickian transport caused by the complex structure of the fractured aquifer.

## 5. Discussions

By using finite elements to discretize the spatial domain, the fractional diffusion equations can be reduced to approximate fractional relaxation equations that possess analytical solutions. The semi-discrete method can not only compute time-fractional diffusion equations in long time range at low computational cost but also deliver accurate numerical predictions for complex and large spatial problem domains. The accuracy in the spatial domain can be improved by using more elements, high-order elements or elements based on advanced finite element formulations. Because the exact solution is used in the time domain, the stability and convergence conditions of the proposed method can be easily satisfied. It can be said that the proposed method is more robust than previous ones.

The main restriction on the proposed method is that the weak forms of the time diffusion equations can be transformed into the following form: 5.1In cases that the source term, physical parameters and/or boundary conditions are only weak function(s) of time, a multiple time-step method can be used.

## 6. Concluding remarks

From the formulations and examples presented, it is clear that a class of time-fractional diffusion equations can be easily computed and the heavy-tail feature can be accurately characterized by the presented semi-discrete method. Our future research work will focus on advanced finite element formulations, such as hybrid elements [49], to compute temporal–spatial fractional diffusion equations that characterize more complex contaminant transport problems.

## Acknowledgements

The first author thanks Prof. G. Pohll and Prof. M. M. Meerschaert for providing the experimental data, and Dr Q. H. Zhang for valuable discussions on finite element programming. The work described in this paper was supported by the National Natural Science Foundation of China (project no. 11202066), the National Basic Research Program of China (973 project no. 2010CB832702), the R&D Special Fund for Public Welfare Industry (Hydrodynamics, project no. 201101014) and the Opening Fund of the State Key Laboratory of Structural Analysis for Industrial Equipment (project no. GZ0902).

## Footnotes

One contribution of 14 to a Theme Issue ‘Fractional calculus and its applications’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.