Convergence Analysis of a Numerical Method for a Fractional Model of Fluid Flow in Fractured Porous Media

: The present paper is devoted to the construction and study of numerical methods for solving an initial boundary value problem for a differential equation containing several terms with fractional time derivatives in the sense of Caputo. This equation is suitable for describing the process of ﬂuid ﬂow in fractured porous media under some physical assumptions, and has an important applied signiﬁcance in petroleum engineering. Two different approaches to constructing numerical schemes depending on orders of the fractional derivatives are proposed. The semi-discrete and fully discrete numerical schemes for solving the problem are analyzed. The construction of a fully discrete scheme is based on applying the ﬁnite difference approximation to time derivatives and the ﬁnite element method in the spatial direction. The approximation of the fractional derivatives in the sense of Caputo is carried out using the L1-method. The convergence of both numerical schemes is rigorously proved. The results of numerical tests conducted for model problems are provided to conﬁrm the theoretical analysis. In addition, the proposed computational method is applied to study the ﬂow of oil in a fractured porous medium within the framework of the considered model. Based on the results of the numerical tests, it was concluded that the model reproduces the characteristic features of the ﬂuid ﬂow process in the medium under consideration.


Introduction
Comprehensive study of the fluid flow phenomenon in fractured porous media has received considerable attention due to its practical significance in variety of industrial processes, particularly in petroleum engineering. This is due to the fact that a significant percentage of world's oil reserves are concentrated in these complex formations, and therefore accurate representation of the flow behavior in the reservoir is relevant due to involvement of these deposits in production.
In recent decades, several conceptual approaches to modeling flows in fractured media have been proposed, an extensive overview of which is presented in [1,2]. The fundamental conceptual models include single-and multi-continuum models, discrete fracture matrix and discrete fracture network models. The works of many authors (e.g., [3][4][5]) are aimed at a qualitative comparison of these models. Some authors [6][7][8] are of the opinion that the fluid motion in a fractured medium cannot be adequately described within the framework of classical theory of fluid flow in porous media, assuming that the nature of the flow depends not only on the current state of the process but also on the previous history of changes in this process. This phenomenon is known as memory in petroleum engineering. Caputo [9] introduced a memory formalism using the apparatus of fractional calculus and proposed the time fractional modification of the classical Darcy's law. This idea was further developed by many authors (e.g., [6,7,[10][11][12][13][14]), and a number of fractional differential generalizations of the classical model have been developed that take into account the memory and spatial correlations effects in porous media (see a literature review in [15]). In this case, the order of the fractional derivative determines the degree of memory impact on patterns of the flow. The experimental studies [16,17] show that fractional differential equations can be effectively used to obtain a more realistic description of fluid flow in porous media.
There are several definitions of the fractional derivative and fractional integral in fractional calculus, an overview of which can be found in [18,19]. In the development of mathematical models of fluid flow in porous media under different physical assumptions, various authors have used fractional time derivatives in the sense of Caputo [7][8][9]20,21], Riemann-Liouville [6,22], Caputo-Fabrizio [23,24], Atangana-Baleanu [25] and others [26] to account for memory effects, whereas the Riemann-Liouville derivative is mainly used to describe nonlocality in the spatial direction [6,27]. In the study of dynamic processes in fractured fractal media, a fractional temporal derivative in the sense of Caputo has an advantage, since it is given by the convolution of the power law kernel and the derivative of the function. In the case of using the Riemann-Liouville time derivatives, difficulties arise in determining and physical treatment of initial conditions [28]. Despite the obvious, from a computational aspect, advantages of the recently developed Caputo-Fabrizio derivative associated with the non-singularity of the kernel, some authors doubt the correctness of its application to account for memory effects [29].
In [6], a new model of fluid flow in a fractured porous medium is derived. The proposed model is based on utilizing the fractional differential analogue of the fluid motion law describing fluid flow in natural fractured porous media in which the fractures are distributed uniformly on average over the volume. The model is also constructed under the assumption that porosity and density are functions not only of pressure but also its fractional time derivatives. However, the paper does not pay attention to the numerical implementation of the mathematical model. Furthermore, the literature review did not reveal any papers devoted to the numerical implementation or obtaining an analytical solution to this model.
In this paper, we conduct for the first time a numerical analysis of the mathematical model of fluid flow in a fractured porous medium proposed in [6] with some reasonable simplifying assumptions. The governing equation of the model contains several terms with fractional derivatives of various orders α, β ∈ (0, 2), γ ∈ (0, 1), as well as a term with the derivative of a function of some expression that contains both time fractional and integer spatial derivatives. Unlike [6], we use the Caputo fractional derivative which does not lack a physical interpretation of the initial conditions.
We consider three special cases of the fractional differential equation and propose two different implicit finite difference/finite element approaches to constructing a numerical scheme for solving the problem depending on orders of the fractional derivatives. The first approach, applied to the case α, β, γ ∈ (0, 1), is based on the use of a higher order approximation to the integer temporal derivative. The second numerical scheme, applied to the cases α ∈ (1, 2), β, γ ∈ (0, 1) and β ∈ (1, 2), α, γ ∈ (0, 1), uses a different approach and is based on utilizing of a Crank-Nicolson type scheme.
We rigorously prove the convergence of both semi-discrete and fully discrete schemes for both constructed numerical methods. The theoretical convergence order of both computational schemes is confirmed by the results of numerous computational experiments.
Further, we consider a slightly more realistic example, which is aimed at studying the flow of oil in a fractured porous medium, and present the results of computational experiments depending on orders of the fractional derivatives. In addition, the results of the numerical implementation of the classical model of fluid flow in a porous medium are presented for comparison. Based on the results of the numerical tests, it was concluded that the model reproduces the characteristic features of the process of fluid flow in a fractured porous medium.
The paper starts by providing a brief derivation of the model equation and defining the problem formulation in Section 2.1. In Section 2.2, semi-discrete formulations of the problem are determined and convergence theorems are proved. The construction of fully discrete schemes and the proof of two convergence theorems is given in Section 2.3. Section 3 describes the experimental study and presents the results of several numerical tests conducted to verify the theoretical analysis, while Section 4 provides the discussion of the obtained results.

Formulation of the Problem
In the classical theory of fluid flow in porous media, the governing equations of the single-phase isothermal flow of a viscous fluid through an isotropic porous medium include the continuity equation and Darcy's law: where φ and k are the porosity and absolute permeability of the medium; p, ρ, andμ are the pressure, density, and viscosity of the fluid, respectively; U is the velocity vector; f is the density of mass sources. For simplicity, the influence of gravitational forces is neglected. One approach to describe fluid flow through a fractured porous medium is to replace the above medium with some model homogeneous porous medium with power-law memory. Due to the fact that porosity depends on the pressure of the fluid and on the stress-strain state of a medium exhibiting viscoelastic properties, it is concluded in [6] that porosity is a function not only of pressure but also of its fractional derivative or fractional integral: where ∂ᾱ 0,t refers to the Caputo fractional differentiation operator which for positive real numbers ν such that n − 1 < ν < n, n = 1, 2, 3, ... is defined as [18] The negative order Caputo fractional differentiation operator ∂ −ν 0,t is equivalent to the Caputo fractional integration operator, I ν 0,t , which can be defined as [54,55] where g(t) is a function defined in [54]. In addition, we use the following fractional differential generalization of the classical state equation [56]: To account for the effect of fractures on the fluid flow process, we use the following generalized motion law which describes fluid flow in natural fractured porous media in which the fractures are distributed uniformly on average over the volume [6]: where F (z) is a given function. Substituting (3)-(6) into (1), we obtain the following nonlinear fractional differential equation describing the flow of a viscoelastic fluid in a fractured porous medium: wheref 0 = generalized fractional differential counterparts [6]. A complete analysis of the resulting Equation (7) is quite difficult due to the presence of a large number of parameters. Let us introduce some simplifying assumptions that are not overly restrictive in practical applications. First, we assume that the flow is onedimensional. Following [6], we neglect the last term in the left-hand side of (7) due to its smallness compared to the fourth term. Furthermore, the porosity function φ is assumed to be linear with respect to its arguments: where b 1 , b 2 are some positive constants. Let us further assume that the fluid density explicitly depends only on the fractional derivative of the pressure, and there exists a positive constant quantity b 3 such that φc fβ = b 3 . Taking into account the above assumptions, (7) reduces to the following form: when α, β, γ ∈ (0, 1); whenᾱ, γ ∈ (0, 1); whenβ, γ ∈ (0, 1).
The following additional assumptions will be used throughout the paper: Assumption 1. Problem 1 has a unique solution having sufficient number of derivatives required to perform the analysis.

Assumption 2.
In Case I, F is of the form where χ is a given function, µ > 0 is a constant. In Cases II and III, F is of the form Assumption 3. In Cases II and III, we assume thatc φα =c f β .
We adopt the standard notations of spaces L q (Ω) and Sobolev spaces W k,q (Ω) with the special case q = 2 denoted by H k (Ω) [57]. In addition, (·, ·) denotes the scalar product in L 2 (Ω) and · k denotes the norm in H k (Ω). Now we define the variational formulation of Problem 1. Since the sought solution belongs to different functional spaces depending on orders of the fractional derivatives,ᾱ, β andγ, let us consider the variational formulations of Problem 1 separately in accordance with Cases I, II, and III.

Construction of the Semi-Discrete Schemes
To define a semi-discrete formulation of Problem 1 with respect to time, it is necessary to approximate both integer and fractional derivatives in time. To this end, we introduce a uniform partition of the time interval [0, T] by points t n = nτ, τ > 0, n = 0, 1, . . . , N such that Nτ = T. Let p n denote a semi-discrete approximation of the function p at t = t n . Additionally, we sometimes use the notation p(·, t) = p(t) for brevity.
We use the following approximation formula for the fractional derivative in the sense of Caputo. where Moreover, the following estimate holds for r ν n = ∂ ν 0,t p(t n ) − ∆ ν 0,t p n : Some obvious properties of the coefficients δ ν n,s presented in Lemma 1 are as below.

Lemma 2.
The following properties hold for the coefficients δ ν n,s , s = 1, 2, . . . , n: In this paper, we propose two different approaches to constructing numerical schemes depending on the orders of the fractional derivatives. The first approach is applied to Case I and is based on the use of the fractional derivative approximation (15) and a higher order approximation of the integer temporal derivative which is defined as Let us define a semi-discrete formulation of Problem 1 corresponding to Case I as follows.

Problem 5.
Let p i n−1 i=0 , p i ∈ H 1 0 (Ω) be given where p 0 = p 0 (x), x ∈ Ω and n ≥ 1. Find p n ∈ H 1 0 (Ω) satisfying the following identities for all v ∈ H 1 0 (Ω) and α, β, γ ∈ (0, 1): The second proposed numerical scheme is applied to Cases II and III and is based on the use of the fractional derivative approximation (15) and the Crank-Nicolson type scheme. To this end, introduce the notations Then the semi-discrete formulations corresponding to Cases II and III is stated as follows.

Error Estimates of the Semi-Discrete Schemes
First, we prove a few auxiliary lemmas. Lemmas 3 and 4 below are used to obtain estimates of terms containing fractional and integer time derivatives. Lemma 3. The following inequality holds for ∆ ν 0,t p n , 0 < ν < 1, n ≥ 1: Proof. It follows from Lemma 1 that Using the Cauchy-Schwarz inequality and Lemma 2 taking into account the elementary inequality ab ≤ 1 2 a 2 + b 2 , we obtain The lemma is proved.

Lemma 4.
Let the sequence p i n i=0 , p i ∈ L 2 (Ω) be given. Then Proof. This lemma can be verified directly.
In addition, we formulate the discrete analogue of Gronwall's lemma which will be used several times throughout the paper.

Lemma 5 ([58]
). If {a n } and {b n } are two positive sequences, {c n } is a monotone positive sequence, and they satisfy the inequalities a 0 + b 0 ≤ c 0 , a n + b n ≤ c n + λ n−1 ∑ i=0 a i , λ > 0, then the following estimate holds: a n + b n ≤ c n exp(nλ), n = 0, 1, 2, ... Now we prove the main results of the section. The symbol C below denotes a positive constant that does not depend on τ and may take different values.
Taking into account inequalities (24)-(26), we obtain from (22) and (23) that where Sum (29) with respect to n from 2 to n, make use of Young inequality to get c φα δ α n,n +c f β δ β n,n π n 2 0 + µδ γ n,n π n where σ = min{α, β}. Using Lemma 1 as well as the facts that and δ ν n,n −1 < Γ(2 − ν)T ν , we arrive at By applying Lemma 5, we conclude that Estimate the right-hand side of (28) using the Cauchy-Schwarz and Young inequalities and apply Lemma 1 to obtain Making use of (31), we get Combining (33) and (34), we arrive at the inequality Finally, by extracting the square root of both sides of the last inequality and applying the inequality |a i | and noticing that 1 √ n ≥ τ T , we arrive at the statement of the theorem.
Before we prove the convergence theorem for Case II, we state the following lemma.

Remark 1.
The convergence result of the semi-discrete scheme corresponding to Case III is similar to Theorem 2.

Construction of the Fully Discrete Schemes
Let K h denote a uniform partition on Ω with a mesh size h. For l ∈ N, let P l (Ω) denote the space of polynomials of degree no greater than l on Ω. Define the finite element space as follows: Define the projection operator Q h : The operator Q h has the following well known property: Based on the semi-discrete scheme (17), (18) constructed in Section 2.2, we propose the following fully discrete scheme corresponding to Case I.
Similarly, we propose the following two fully discrete schemes that correspond to Cases II and III and are based on the semi-discrete schemes (20) and (21).

Problem 10.
Let p i h n−1 i=0 , p i h ∈ V h be given such that p 0 h = Q h p 0 and ∆ t p 1/2 h = Q h u 0 . Find p n h ∈ V h , n = 2, 3, . . . , N satisfying the following identity for any v h ∈ V h : whereβ, γ ∈ (0, 1).

Error Estimates of the Fully Discrete Schemes
, p i h ∈ V h be the solution of Problem 8 and p be the solution of Problem 1. Then p n h ∈ V h satisfies the following inequality for α, β, γ ∈ (0, 1) and for all τ such that τ < 1/6 under Assumptions 1 and 2: Proof. Let p(t n ) − p n h = (p(t n ) − Q h p n ) + Q h p n − p n h = ψ n + ξ n . By substracting (44) from (18) and taking v h = ξ n , we have By using Lemma 3 under Assumption 2, we obtain where the notation (30) with is used. Next, using Lemmas 1 and 2, Cauchy-Schwarz and Young inequalities gives where t n−1 ≤ ζ s ≤ t n . In a similar way, one gets To estimate the last term in the left-hand side of (47), apply Cauchy-Schwarz and Young inequalities, as well as the inequality b a f dx Then taking into account (48)-(51), it follows from (47) that Sum (52) with respect to n from 2 to n to obtain By assuming that τ satisfies the condition τ < 1/6, and by applying Lemma 5, we have By considering the substraction of (17) and (43), taking v h = ξ 1 and using the same technique of estimating the terms in the resulting identity similarly to obtaining the inequalities (48)- (51), we arrive at the inequality By combining (53) and (54), we obtain Using the triangle inequality, we have Then taking into account Theorem 1, (42) and (55), we arrive at the assertion of the theorem. Lemma 7. Let {p n } N n=0 be the solution of Problem 6. Then Proof. By using Lemma 1, it is easy to show that where z(t n ) = ∂ t p(t n ) +c φα ∂ᾱ 0,t (∂ t p(t n ) + p(t n )). Then using the triangle inequality and property (42) yields (57).
Now we discuss the convergence of the fully discrete scheme corresponding to Case II.
be the solution of Problem 9 and p be the solution of Problem 1. Then p n h ∈ V h satisfies the following inequality forᾱ, γ ∈ (0, 1) and for all τ such that τ < 2 under Assumptions 1-3: Proof. Denote ξ n = Q h p n − p n h . Using (20), (45) and (42), we have where R n−1/2 h is defined in (58). By taking v h = ∆ t ξ n−1/2 + ξ n−1/2 and applying the technique used to obtain inequalities (37)-(40), we arrive at where Θᾱ n = Hence, taking into account (60)-(63), it follows from (59) that where Θ n =c φα Θᾱ n + µΘ γ n . By summing (64) with respect to n from 1 to n, then using Cauchy-Schwarz and Young inequalities and applying the same technique as in obtaining the inequality (32) and noticing that Assuming that τ satisfies the condition τ < 2, we conclude that Finally, by applying Lemma 5 and taking into account (42), (56), (65) and Theorem 2, we arrive at the statement of the theorem.

Remark 2.
The convergence result of the fully discrete scheme corresponding to Case III is similar to Theorem 4.

Verification of the Convergence Order
To verify the theoretical convergence estimates obtained in Theorems 3 and 4 for the fully discrete schemes, a number of computational experiments were carried out on the example of two model problems with known exact solutions. Numerical algorithms were implemented in FreeFEM++.
The purpose of the computational experiments carried out for Example 1 below was to determine the dependence of the empirical convergence order on orders of the fractional derivatives, and compare the empirical convergence order with the theoretical convergence order obtained in Theorem 3.

Example 1. Consider the fractional differential equation
, subject to the following initial and boundary conditions: where α, β, γ ∈ (0, 1). The exact solution of the problem is p( First, we studied the convergence behavior of the numerical solution with respect to the time step. Successively halving the time step in the range from τ = 1/10 to τ = 1/160 for a fixed number of elements M = 20,000, we calculated the corresponding errors E τ in the L 2 -norm and the convergence order by the formula The choice of a relatively large M was motivated by the intention to reduce the influence of the error with respect to the spatial variable in this numerical test. Several combinations of the fractional derivative orders from the set {0.1, 0.5, 0.9} were considered, and the values α and β were assumed to be equal for simplicity.
The calculation results are presented in Table 1. In the "Order" column, the theoretical convergence orders are indicated in parentheses.
Error plots are shown in Figure 1. It follows from Table 1 that the convergence order with respect to the time step τ significantly depended on the orders of the fractional derivatives α, β, γ. Specifically, it is easy to see that the empirical convergence order in L 2 -norm was around 1.90 in case of α = β = γ = 0.1 and it decreased significantly with increasing order of the fractional derivatives. The convergence order was close to 1.10 when any of α, β or γ was equal to 0.9. This behavior agreed well with the theoretical order O τ 2−α + τ 2−β + τ 2−γ .
In order to validate the convergence behavior with respect to the mesh size h, we fixed the time step τ = 1/20,000 and considered different values of the mesh size in the range from h = 1/5 to h = 1/40. The corresponding calculation results are presented in Table 2 and Figure 2.   Table 2 shows that the spatial convergence order was close enough to 2.00 for all considered values of the fractional derivative orders, and therefore this result agreed well with Theorem 3. One can notice that the spatial convergence order decreased slightly as any of α, β or γ approached 1. This can be explained by the influence of the error with respect to the time step due to an insufficiently small value of τ, τ = 1/20,000, and the fact that the temporal convergence order was around 1 in this case.
Thus, it can be concluded that the numerical results conducted for Example 1 fully confirmed the theoretically predicted convergence order obtained in Theorem 3.
Similarly, the purpose of the computational experiments carried out for Example 2 below was to compare the empirical convergence order with the theoretical convergence order obtained in Theorem 4. The numerical tests were mainly aimed at validating the convergence order with respect to the time step τ. We considered different values of τ in the range from 1/10 to 1/160 with a fixed mesh size h = 1/40,000. Example 2. Consider the fractional differential equation x ∈ Ω = (0, 1), t ∈ (0, 1] subject to the following initial and boundary conditions: whereᾱ, γ ∈ (0, 1). The exact solution of the problem is p(x, t) = t 2 sin(πx).
In the experiments, several combinations of the fractional derivative ordersᾱ and γ from the set {0.1, 0.5, 0.9} were considered. The results of the computational experiments are presented in Table 3 and error plots are shown in Figure 3.   Let us recall that the motion law (6) is suitable for describing the process of fluid flow in fractured porous media with the uniformly distributed fractures over the volume [6]. Figure 4 shows a deceleration of the process of fluid flow, and it can be clearly seen that an increase in the orders of the fractional derivative led to a greater deceleration. This behavior confirmed the fact that the porous medium had a retarding effect on the flow process, and the orders of the fractional derivatives determined the degree of memory impact on patterns of the flow.

Discussion
Thus, having carried out a theoretical analysis of the constructed numerical methods and having analyzed the results of the computational experiments, the following conclusions can be outlined.
(1) Both computational methods constructed for solving Problem 1 in the special casesᾱ,β ∈ (−1, 0), γ ∈ (0, 1) andᾱ, γ ∈ (0, 1),β =ᾱ − 1 (orβ, γ ∈ (0, 1),ᾱ =β − 1) allow obtaining an approximate solution with the second order of accuracy in the spatial direction and an order of 2 − ν with respect to the temporal variable, where ν is the largest of the fractional derivative orders. The results of computational experiments conducted for various orders of the fractional derivatives and grid configurations fully confirm the results of the theoretical analysis.
(2) Based on the results of the numerical tests presented in Section 3.2, we conclude that the model under study reproduces the characteristic features of the process of fluid flow in a fractured porous medium. In particular, the test carried out within the framework of the studied model showed that the fluid flow process in a fractured porous medium proceeds more slowly than in a continuous medium. This observation is also consistent with the conclusions of the paper [8].
In subsequent papers, the authors intend to consider the initial boundary value problem for the Equation (7) with a more general form of the function F . Specifically, more realistic examples of fluid flow problems in a fractured porous medium will be considered, and a more extensive comparison of the simulation results with the results of similar works will be carried out.
(3) Since this is the first work devoted to the numerical implementation of the model under consideration, the authors restricted themselves to constructing numerical schemes of nearly second order of convergence. The results obtained can be improved by using higher-order approximations of the Caputo fractional derivative (see, for example, [38,51]). A separate paper will be devoted to this study. (4) The conclusions obtained in this study can be used in solving other classes of fractional differential equations. The results of the presented work will also form the basis for further research in this direction. For example, a natural continuation of our work is the study of multi-dimensional models of fluid flow in a fractured porous medium with fractal geometry of fractures and models of single-phase fluid flow with spatial fractional derivatives.