Identifying the Fractional Orders in Anomalous Diffusion Models from Real Data

An attempt is made to identify the orders of the fractional derivatives in a simple anomalous diffusion model, starting from real data. We consider experimental data taken at the Columbus Air Force Base in Mississippi. Using as a model a one-dimensional fractional diffusion equation in both space and time, we fit the data by choosing several values of the fractional orders and computing the infinite-norm “errors”, representing the discrepancy between the numerical solution to the model equation and the experimental data. Data were also filtered before being used, to see possible improvements. The minimal discrepancy is attained correspondingly to a fractional order in time around 0.6 and a fractional order in space near 2. These results may describe well the memory properties of the porous medium that can be observed.


Introduction
In this paper, we provide some justification of the widespread use of simple linear parabolic fractional partial differential equations (fPDEs) to model a number of physical phenomena, such as, for instance, diffusion of fluids in porous media [1] and in particular fluid flow of tracers through porous media.We do that by identifying the fractional orders of the aforementioned fPDE from experimental data.
Several laboratory experiments have been performed concerning one-dimensional flows, to determine the flow rate through uniformly packed columns filled in by certain porous media, where a constant pressure was applied to the lower boundary, see [2,3].These measurements were done using diffusivity meters.The experimental setup was designed aiming at determining the memory properties of the observed dynamical behavior.The latter could be considered as an anomalous diffusion, and was well reproduced by a suitable fPDE modeling of the constitutive laws involving pressure gradient-flux and pressure-density variations.In such experiments, one may use water as fluid and a variety of porous media.Clearly, one of the key-factors was the size (or better, the size distribution) of the specific medium's grains.In the laboratory experiments of [2,3], special attention was paid to detect memory effects, experimentally identifying the fractional orders in time in the aforementioned constitutive equations.
On the theoretical side, there are works in the literature on several kinds of inverse problems for fractional differential equations [4][5][6][7], but very little attention seems to have been devoted to determining the fractional orders.We can mention two papers by Bondarenko and Ivaschenko, concerning only some special problems for time fractional diffusion equations [8,9], and the recent contribution in [10] where the focus is on identifying the order in simple fractional ordinary differential equations, exploiting asymptotics and special (Mittag-Leffler) functions.See also [11,12].
In this paper, in contrast to what was done in [2,3], we compare certain real, i.e., experimental data taken at the Columbus Air Force Base in Mississippi, USA [13,14], with those provided by the numerical solution of a certain fPDE, by varying both fractional orders in time and in space, assessing the mismatch for certain fixed values of the orders of fractional derivatives as well as of the initial conditions.The real data we used here refer to the MADE-2 experiment, where some tritium concentration was injected, representing contaminant transport in a highly heterogeneous aquifer [13,14].In such papers, the 1D data were actually obtained by suitably processing typical 2D or 3D concentration data sets.This concentration depends, in general, on both time and place (measured in days and meters, respectively).We considered three fixed times of observation, when they were obtained.Our approach can be considered in the spirit of the modern attitude of "learning from data", in such a way to infer the governing equations of the physical model, as put forth at the present time by machine learning techniques, within the modern Data Science [15].However, the idea here is grounded on basic physical principles.
Here is the plan of the paper.In Section 2, we state the problem with some emphasis on the numerical method we adopted, in Section 3 we describe the discretization method used to solve numerically the model equation, and in Section 4 we illustrate our method through some numerical examples.Conclusions and future developments are given in a final section, at the end of the paper.

Statement of the Problem
Consider the one-dimensional anomalous diffusion equation, fractional in both, space and time, where u = u(x, t) represents some mass concentration (of tritium), 0 D γ t u is the Caputo fractional derivative in time of u, ∂ α u ∂|x| α is the Riemann-Liouville fractional derivative in space of u, a > 0, 0 for some x, x, T > 0, along with the boundary conditions u(x, t) = u 1 , u(x, t) = u 2 , and the initial condition u(x, 0) = u 0 (x), for some constants u 1 , u 2 , and some u 0 (x), while f (x, t) is a certain external source (if any).
We solved such a problem numerically using finite differences of the kind of Grünwald-Letnikov (GL) for our approximations [16].We started considering a mesh in the space-time domain [x, x] × [0, T], denoting with u n j the numerical approximation of the solution at the mesh points, (x j , t n ), i.e., u n j := u(x j , t n ).We then replaced the continuous fractional derivatives operators, in both, time and space, in Equation (1), with their GL approximation.We thus obtain a difference equation whose solution (u n j ) will provide estimated values of the exact solution, u(x, t), at the mesh points.For a given fractional differential operator, there are several possible forms of difference operators, and hence there are several different finite difference methods capable to solve our problem.
We should recall, at this point, that many authors have already considered similar problems, but most often adopting a fractional order only either in time or in space, see, e.g., [17][18][19][20].When fluid flow through porous media is studied, usually, only fractional derivatives in time were considered.
Here we take a uniform spatial step-size, ∆x := x i+1 − x i , and consider the space fractional differential operator in the sense of Riemann-Liouville, that we will approximate by means of "GL differences".Fractional time derivatives in the sense of Caputo will be discretized again by GL differences, with a constant time step-size, ∆t := t n+1 − t n .
We thus obtain a finite difference algorithm, ending up with a linear system like where and A and S represent the discretized time and space operators (two matrices), respectively, obtained through the GL discretization, see [16].
A code was then written to compute U n , by minimizing, at each time step, the ∞ norms (with respect to space) of the discrepancies between the experimental and the numerically computed values obtained by solving (numerically) the aforementioned anomalous diffusion equation.We will denote with v n j the experimental data, measured in the field, and with the vector of the measured data.
Relatively recently, in [2], experimental results have been obtained to estimate the physical parameters governing the dynamical behavior of the flow of tracers through a certain given porous medium, in particular to observe the aforementioned memory effects described by parameters such as, for instance, the order of the fractional derivative with respect to time.Memory effects may also be due to the size of the tank containing the porous medium in the laboratory, the particle size distribution, the stability of the initial particle distribution.
The function to minimize is then the error, i.e., the discrepancy between the numerically computed value, U n , and the value of u(x, t) obtained by using the real experimental data, V n , that is

Discretization of Partial Derivatives in Time and Space
Let us consider the nodes (j∆x, n∆t), for a given fixed j, j ∈ {1, 2, . . ., M}, n = 0, 1, 2, . . ., N, corresponding to all time levels at the j-th space discretization node.Then, the approximation scheme can be represented through a matrix that we will call B γ N [16], so that we can write where u γ j,n represents the numerical approximation of the time derivative of order γ at the (j, n) node of the mesh.
In order to obtain a simultaneous approximation of the time derivative of order γ and space derivatives of order α of u(x, t), at each grid point shown in Figure 1, we rearrange all function values, u n j , at the discretization nodes, in a column vector, Figure 1.Discretization nodes for partial derivatives, see [16].
We can arrange all the numerical values of the solution in a matrix, which contains the numerical results of every discretized value of time and space.Such a matrix represents the iteration matrix which yields the discretization vector U γ t of the fractional partial derivative of order γ with respect to time, and can be written as the Kronecker product of the matrix B γ N , corresponding to the fractional ordinary derivative of order γ, and the identity matrix, I M , see [16].Recall that here N is the number of time steps and M is the number of space discretization nodes.The corresponding stencil is depicted in Figure 1.
Similarly, the matrix governing the transformation of the vector U n into the vector U α x of the fractional partial space derivative of order α can be obtained as the Kronecker product of the unit matrix I N (where N is the number of the time steps) in the mesh.The matrix R α N , which corresponds to a derivative of order α [21,22] is given by All this is illustrated in Figure 1, where the nodes in black and gray (corresponding to all discretization nodes, from the leftmost to the rightmost nodes) are used to approximate the fractional-order derivative at the node in gray.With these approximations for the fractional partial derivatives made with respect to both variables, we can discretize at once the fractional diffusion equation in (1) (see Figure 2), that is without any forcing term, obtaining Note that the system's matrix is banded upper triangular.1), see [16].Here, the Riemann-Liouville derivative is also denoted with D α x .
In closing this section, we observe the following.There are currently in the literature still unclear issues concerning imposing local boundary values to nonlocal differential equations such as fPDEs, see, e.g., [23].When the domain is bounded, in many cases it seems reasonable to extend the Dirichlet data to all space, outside the domain.Here we will adopt this strategy.

Numerical Examples
In this section, we present some numerical examples, where we show that the identification of both fractional orders of derivative, in space and time, can be accomplished, from the knowledge of the available real data.We attempt to model a certain real phenomenon either taking into account a possible source term or not.The real data were also filtered prior to using them.

Example 1
Consider the one-dimensional fractional diffusion equation characterized by fractional derivatives in both, space and time (measured in meters and days, respectively), satisfied by the mass concentration of tritium, u = u(x, t), measured in picoCuries per milliliter, with a = 10 −5 , 0 < γ ≤ 1, 1 < α ≤ 2, T = 328, under the homogeneous boundary conditions u(0, t) = u(1, t) = 0 (that will be extended to the entire real line), and the initial condition u 0 (x) := x(1 − x).Some values for a can be found in the literature, often much smaller than the value we used here [24,25], but the value we used seemed to be more appropriate in the present case as a modeling choice.Both, the initial profile and the homogeneous boundary values are inferred by inspecting the available real data.The values of u, obtained solving numerically the problem in ( 15) should be compared with the experimental data, see Figure 3.Such data come from tests made at the Columbus Air Force Base in Mississippi, see [13,14].We fitted them against the results of the numerical integration of Equation ( 15) (i.e., the numerically approximated values u n j of u(x, t)), for several values of the parameters involved in the model.All data have been normalized to the spatial range [0, 1], by dividing by their largest value.In Tables 1-3 below, the discrepancies between the numerical solution and the experimental data are shown, for several chosen values of γ and α.The notation 1.1(0.1)2.0 means taking all values between 1.1 and 2.0 at steps of size 0.1, that is considering all values 1.1, 1.2, 1.3, . . ., 1.9.In all tables, M = 14 is the number of points of the space-grid, and N = 328 is that of time steps.In Figure 4, the error versus α and γ is plotted in log scale.

Example 2
Consider the one-dimensional fractional diffusion equation fractional in both, space and time, satisfied by the tritium concentration, u = u(x, t), with a = 10 −5 , 0 < γ ≤ 1, 1 < α ≤ 2, T = 328, under the homogeneous boundary conditions u(0, t) = u(1, t) = 0 (that will be extended to the entire real line), and the initial condition u(x, 0) := u 0 (x) ≡ 0.Here δ are the delta functions.Indeed, introducing such a forcing term should more realistically account for describing a source located initially in the middle of the space domain; see the experimental data in Figure 3. Better results are in fact obtained doing so.
In the numerical simulations, the delta functions on the right-hand side of ( 15) have been implemented using suitable gaussian functions.We chose the gaussian functions exp(− , and exp(− t 2 2•0.025 2 ).We were interested in using peaked functions but with a prescribed maximum value, i.e., the value 1, in view of the data normalization.
In Tables 5-7 below, we show the discrepancies between the numerical solution and the experimental data, for several values that we chose for γ and α.
Table 5. γ = 0.1(0.1)0.4 and α = 1.1(0.1)2.0.(γ, α) We can see that, again, the best results are achieved for (γ, α) = (0.6, 2.0).Using as before N = 328 and M = 14, we obtained about which is appreciably smaller than that of Example 1.Note however that having a smaller value for the minimal discrepancy is only meaningful in view of a reliable identification of the optimal parameter γ and α.In fact, different data coming from the same experiments will provide, in general, different discrepancy.What matters is that again the same values for γ and α have been inferred.An error surface is shown, as a function of α and γ in Figure 7.
In Figure 8, the error is plotted as a function of α, for a fixed value of γ (γ = 0.6).

Conclusions and Future Developments
The aim of this work is that of identifying in a practical way the time and space orders of the fractional derivatives in a certain anomalous diffusion model, on the basis of the knowledge of real data.We used as a model a one-dimensional fractional diffusion equation in both, space and time, and fit the data by testing several possible, reasonable values for the fractional orders.We observed that, in the specific case we considered, the minimal discrepancy was attained correspondingly to a time fractional order around 0.6-0.7,but very close to 0.6, and a space fractional order around 1.9-2.0,but very close to 2.0.This means that the porous medium under investigation behaves approximately as being governed by a Laplacian in space, but the delays and memory effects are significantly described by a time fractional order very different from 1 (and appreciably below 1).Some modifications of the basic model we adopted in this paper might be used to handle some other set of data, pertaining to other experiments or sites.Improving or refining our results can also be considered, making the model more sophisticated and realistic, for instance, by introducing some terms on the right-hand side. of the equation, to represent external forces applied to the medium, or advection terms to account for the flow of the contaminant in heterogenous porous media.Furthermore, a viscoelastic behavior could be described coupling the model equation with other specific equations [28].In this way, the elasticity and viscosity coefficients might enter the model, and hence the relevant physical dynamics could be better described.
Clearly, further refining the spatial grid would improve the accuracy of the numerical solution, but it is impractical and costly to obtain real data at many locations and at many times, for comparison).Changing the initial condition, or introducing a reaction term in the model Equation (15) could also be considered to improve the model.The case of higher space dimensions would definitely be of great interest.

Figure 4 .
Figure 4. Error surface in log scale, i.e., the log of the error U n − V n ∞ , is shown vs. γ and α, for a = 10 −5 .

Figure 7 .
Figure 7. Error surface in log scale, for a = 10 −5 , i.e., that is, the log of the error U n − V n ∞ , vs. γ and α.