A Priori Estimates for the Solution of an Initial Boundary Value Problem of Fluid Flow through Fractured Porous Media

Abstract

The paper studies a model of fluid flow in a fractured porous medium in which fractures are distributed uniformly over the volume. This model includes a nonlinear equation containing several terms with fractional derivatives in the sense of Caputo of order belonging to the interval $\left(1,2\right)$. The relevance of studying this problem is determined by its practical significance in the oil industry, since most of the world’s oil reserves are in these types of reservoirs. The uniqueness of the solution to the problem in a differential form and its dependence on the initial data and the right-hand side of the equation is proved. A numerical method is proposed based on the use of the finite difference approximation for integer and fractional time derivatives and the finite element method in the spatial direction. A change of variables is introduced to reduce the order of the fractional derivatives. Furthermore, the fractional derivative is approximated by using the L1-method. The stability and convergence of the proposed numerical method are rigorously proved. The theoretical order of convergence is confirmed by the results of numerical tests for a problem of fluid flow in fractured porous media with a known exact solution.

1. Introduction

Predicting the behavior of fluid flow through naturally fractured porous media has received considerable attention due to its fundamental significance to a diverse range of industrial processes including enhanced oil recovery, carbon sequestration, and aquifer remediation. Over the last few decades, several substantially different approaches have been proposed in the effort to describe the flow and transport dynamics in these complex formations. A comprehensive review of the most important approaches has been conducted in [1,2,3,4].

The process of fluid flow in fractured porous media is characterized by anomalous kinetics, which obeys distribution laws with power-law asymptotics. This is due to the complex internal geometric structure of these media, consisting of matrix porous blocks and a system of fractures.

Fractures have a significant impact on the flow pattern in connection with its dependence on the fracture properties. Understanding the fractures propagation is known to be essential in the assessment of the fluid flow process. Numerical methods for fracture modeling include the extended finite element method [5], cracking particles method [6], cracking elements method [7], phase-field method [8], and many others. The latter has been raising much interest by virtue of its simplicity for numerical implementation. Recent studies employing this method aimed at, e.g., fluid-driven fracture modeling [9,10,11], anisotropic fracture modeling [12,13,14], and its utilization to applied problems [15,16].

Another approach to modeling flows in a fractured porous medium is to replace this medium with some homogeneous medium with memory. Fractional differential calculus is an effective tool for accounting for memory effects. In this case, equations describing the fluid flow are replaced by their fractional differential counterparts of order $\left(0,2\right)$.

Fractional differential models of fluid flow are most fully studied in the case when the order of the fractional derivative belongs to the interval $\left(0,1\right)$. For example, in [17], a modification of the classical Darcy’s law is proposed, which depends on the time fractional derivative by introducing a memory formalism to better describe the flow as well as the pressure of fluids. He [18] modified Darcy’s law to overcome the percolation flow assumption. Based on numerous computational experiments, the authors of [19] showed a significant effect of memory on the fluid flow process through a porous matrix.

Various definitions of fractional derivatives were used in the literature when constructing fractional differential models, for example, the derivative in the sense of Caputo [19,20,21,22], Caputo-Fabrizio [23,24], Riemann–Liouville [18,25,26], Atangana-Baleanu-Caputo [27], and Hilfer [28].

There are very few works in which the order of the fractional derivatives in the governing equations belongs to the interval $\left(1,2\right)$. In [25], a model of filtration in a fractured porous medium was derived using a fractional differential analog of the law of motion, as well as under the assumption that porosity and density are functions not only of pressure, but also of its fractional derivatives. The pressure equation in this model includes three terms with fractional time derivatives of the order $\alpha ,\beta ,\gamma \in \left(0,2\right)$. In [29], several numerical methods for implementing this model are proposed for three special cases, depending on orders of the fractional derivatives. The first one covers the case when $\alpha ,\beta ,\gamma \in \left(0,1\right)$. In the second case, it was assumed that $\alpha \in \left(1,2\right)$ and $\beta =\alpha -1$. In addition, when performing a theoretical analysis, a simplifying assumption was made that the coefficients for fractional derivatives were equal. In fact, this imposes serious restrictions on the use of the method in real applications.

Despite the presence of many studies in which methods for solving fractional equations of order $\alpha \in \left(1,2\right)$ are proposed, there are very few works devoted to studying the features of real physical processes for this order. For example, their applications to wave equations in problems of linear viscoelasticity [30], diffusion-wave processes [31,32], signal analysis, random walk of suspended flows, and others are known.

Construction of numerical schemes is based on the use of approximation formulas for fractional derivatives. There are discretization formulas for approximating the derivative in the sense of Caputo of order $\alpha \in \left(1,2\right)$—for example, the L2-method [33], L2C-method [34], and a method based on employing a piecewise interpolation polynomial to approximate the integrand derivative [35]. As a rule, they are more complex, and they use a multi-point pattern, which creates challenges of applying them on the first time layers. These discretization formulas are used, for example, in [34,35] to solve fractional wave equations.

Various approaches are known for the numerical solution of equations containing fractional derivatives of order $\alpha \in \left(1,2\right)$. For example, the finite sinusoidal transform method [36], the finite difference method [35,37,38,39], the local meshless method [40], the finite element method [41,42], the Galerkin spectral method [43], the collocation method using cubic B-splines [44], and others. Unfortunately, the literature review revealed very few works (for example, [43]) that provide a theoretical analysis of the stability and convergence of the proposed numerical methods for equations of the considered order.

Given the importance and high practical significance, the aim of this paper is to analyze the model of fluid flow in a fractured medium, the order of the governing equation which belongs to the interval $\left(1,2\right)$. The work proposes several results.

First, the uniqueness of the solution in a differential form and its dependence on initial data and the right-hand side of the equation is proved.

Secondly, a computational method is proposed based on the application of the finite element method in the spatial direction and the finite difference method with respect to the time variable. To reduce the computational complexity associated with the calculation of fractional derivatives, a change of variables has been introduced that reduces the order of fractional derivatives by one. In contrast to the results of [29], in this paper, we consider a more general case that does not impose restrictions on the orders of fractional derivatives and the coefficients of equations.

Thirdly, the stability and convergence of the proposed numerical method is rigorously proved. Furthermore, the theoretical order of convergence is confirmed by the results of numerical tests for a problem with a known exact solution.

The main contribution of this study consists of a rigorous theoretical analysis of the method for implementing a previously unexplored model of fluid flow in fractured porous media.

The paper is organized as follows: Section 2.1 presents the formulation of the problem under consideration. In Section 2.2, uniqueness of the solution and its continuous dependence on input data are discussed. In Section 2.3, semi-discrete and fully discrete formulations of the problem are determined. Section 2.4 discusses the stability and convergence of the numerical method proposed. Section 3 presents the results of a number of numerical tests to verify the theoretical analysis. The results obtained are discussed in Section 4.

2. Materials and Methods

2.1. Formulation of the Problem

Problem 1.

In ${\overline{Q}}_T=\overline{\mathrm{\Omega }}\times \overline{J}$, where $\mathrm{\Omega }\subset {\mathbb{R}}^d$, $d\ge 1$, $J=\left(0,T\right]$, $T>0$; consider the following initial boundary value problem:

$$\begin{array}{cccc}& {\partial }_{t}p+k_{\alpha }{\partial }_{0,t}^{\alpha +1}p+k_{\beta }{\partial }_{0,t}^{\beta +1}p-\lambda {\partial }_{0,t}^{\gamma +1}{\nabla }^{2}p=f\left(x,t,p\right),\hfill & \hfill & x\in \mathrm{\Omega },t\in J,\hfill \end{array}$$

(1)

$$\begin{array}{cccc}& p=0,\hfill & & x\in \partial \mathrm{\Omega },t\in J,\hfill \end{array}$$

(2)

$$\begin{array}{cccc}& p=p_0\left(x\right),{\partial }_{t}p=u_0\left(x\right),\hfill & & x\in \overline{\mathrm{\Omega }},t=0,\hfill \end{array}$$

(3)

where $\alpha ,\beta ,\gamma \in \left(0,1\right)$, $k_{\alpha }$, $k_{\beta }$, λ are positive constants, f is a given function, and the Caputo fractional differentiation operator ${\partial }_{0,t}^{\nu }$ is defined as [45]

$${\partial }_{0,t}^{\nu }p\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(n-\nu \right)}{\int }_0^t\frac{p^{\left(n\right)}\left(\theta \right)}{{\left(t-\theta \right)}^{1+\nu -n}}d\theta ,n-1<\nu <n,n\in \mathbb{N}.$$

(4)

Problem 1 was originated in [25] to describe one phase fluid flow in fractured porous media. However, that paper did not pay attention to the numerical implementation of the proposed model. In [29], a numerical method was constructed for three special cases of the model depending on orders of the fractional derivatives. In contrast to [29], we consider a more general case which does not impose restrictions on relations between the orders of fractional derivatives.

We assume that Problem 1 has a solution having a sufficient number of derivatives required to perform the analysis. Moreover, we assume that

$$f\left(x,t,p\right)=f_1\left(p\right)+f_2\left(x,t\right),$$

and there are positive real numbers $C_1$ and $C_2$ such that $\left|f_1\left(p\right)\right|\le C_1\left|p\right|$, $f_1^{\prime }\left(p_0\right)=0$, $\left|f_1^{\prime }\left(p\right)\right|\le C_2$. In practice, $f_1\left(p\right)$ can take into account possible nonlinear mass transfer from other continua.

The paper uses the generally accepted notation of spaces $L^q\left(\mathrm{\Omega }\right)$ and Sobolev spaces $W^{k,q}\left(\mathrm{\Omega }\right)$, and, in particular, $W^{k,2}\left(\mathrm{\Omega }\right)=H^k\left(\mathrm{\Omega }\right)$. $\left(·,·\right)$ denotes the dot product in $L^2\left(\mathrm{\Omega }\right)$.

Rewrite (1) in the following form:

$$\left\{\begin{array}{c}{\partial }_{t}p=u,\hfill \\ k_{\alpha }{\partial }_{0,t}^{\alpha }u+k_{\beta }{\partial }_{0,t}^{\beta }u-\lambda {\partial }_{0,t}^{\gamma }{\nabla }^{2}u+u=f_1\left(p\right)+f_2\left(x,t\right).\hfill \end{array}\right.$$

Then, the corresponding variational formulation of Problem 1 reads as below:

Problem 2.

Find $\left\{u,p\right\}:J\mapsto H_0^1\left(\mathrm{\Omega }\right)\times L^2\left(\mathrm{\Omega }\right)$ such that for all $w\in L^2\left(\mathrm{\Omega }\right)$ and $v\in H_0^1\left(\mathrm{\Omega }\right)$:

$$\left({\partial }_{t}p,w\right)=\left(u,w\right),$$

(5)

$$k_{\alpha }\left({\partial }_{0,t}^{\alpha }u,v\right)+k_{\beta }\left({\partial }_{0,t}^{\beta }u,v\right)+\lambda \left({\partial }_{0,t}^{\gamma }\nabla u,\nabla v\right)+\left(u,v\right)=\left(f_1\left(p\right),v\right)+\left(f_2,v\right),$$

(6)

where $\alpha ,\beta ,\gamma \in \left(0,1\right)$.

2.2. Uniqueness of the Solution and Its Continuous Dependence on Input Data

Let us present a few well-known lemmas:

Lemma 1

([46]). The following inequality holds for any function $g\left(t\right)$ absolutely continuous on $\left[0,T\right]$:

$$\left({\partial }_{0,t}^{\nu }g,g\right)\ge \frac{1}{2}{\partial }_{0,t}^{\nu }{∥g∥}_{L^2\left(\mathrm{\Omega }\right)}^2,0<\nu <1.$$

Lemma 2

([47]). Let $\alpha \left(t\right)$, $\gamma \left(t\right)$, and $\lambda \left(t\right)$ be three non-negative functions satisfying the inequality

$$\alpha \left(t\right)+\beta \left(t\right)\le c+{\int }_0^t\lambda \left(s\right)ds+{\int }_0^t\alpha \left(s\right)\gamma \left(s\right)ds$$

for all $t\in \overline{J}$, where $\beta \left(t\right)$ is a non-negative function on $\overline{J}$. Then,

$$\alpha \left(t\right)+\beta \left(t\right)\le \left(c+{\int }_0^t\lambda \left(s\right)ds\right)exp\left({\int }_0^t\gamma \left(s\right)ds\right).$$

First, let us prove the following result.

Theorem 1.

The following inequality holds for the solution of Problem 2:

$${∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\left({∥p_0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{\int }_0^t{∥f_2∥}_{L^2\left(\mathrm{\Omega }\right)}^{2}ds\right),C>0,$$

(7)

which implies the uniqueness and continuous dependence of the solution on input data.

Proof. 

Let us choose $\left(w,v\right)=\left(p,u\right)$ in (5) and (6) to obtain

$$\left({\partial }_{t}p,p\right)=\left(u,p\right),$$

(8)

$$k_{\alpha }\left({\partial }_{0,t}^{\alpha }u,u\right)+k_{\beta }\left({\partial }_{0,t}^{\beta }u,u\right)+\lambda \left({\partial }_{0,t}^{\gamma }\nabla u,\nabla u\right)+\left(u,u\right)=\left(f_1\left(p\right),u\right)+\left(f_2,u\right).$$

(9)

By estimating the terms in (8) and (9) by applying Cauchy inequality and Lemma 1, it is easy to obtain that

$$\frac{d}{dt}{∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2+k_{\alpha }{\partial }_{0,t}^{\alpha }{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2+k_{\beta }{\partial }_{0,t}^{\beta }{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+\lambda {\partial }_{0,t}^{\gamma }{∥\nabla u∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2+{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\left({∥f_2∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right).$$

(10)

Denote $\nu =max\left\{\alpha ,\beta \right\}$, $\mu =min\left\{\alpha ,\beta \right\}$ and consider the two cases. For the first case, $\gamma <\nu$, apply the fractional integration operator ${\partial }_{0,t}^{-\nu }$,

$${\partial }_{0,t}^{-\nu }g=\frac{1}{\mathrm{\Gamma }\left(\nu \right)}{\int }_0^t{\left(t-\tau \right)}^{\nu -1}g\left(x,\tau \right)d\tau ,t\ge 0,$$

to both sides of (10) to have

$${\partial }_{0,t}^{1-\nu }{∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2+k_{\nu }{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2+k_{\mu }{\partial }_{0,t}^{\mu -\nu }{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\lambda {\partial }_{0,t}^{\gamma -\nu }{∥\nabla u∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2$$

$$+{\partial }_{0,t}^{-\nu }{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C{\partial }_{0,t}^{-\nu }\left({∥f_2∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right).$$

(11)

By the definition of the fractional integral, the third, fourth, and fifth terms on the left-hand side of (11) are non-negative. Therefore, it follows from (11) that

$$\frac{d}{dt}{∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\left({∥f_2∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right).$$

(12)

By integrating (12) from 0 to t and applying Lemma 2, we obtain (7).

Similar computations for the second case, $\gamma >\nu$, arrive at the inequality

$${\partial }_{0,t}^{1-\gamma }{∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2+k_{\alpha }{\partial }_{0,t}^{\alpha -\gamma }{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2+k_{\beta }{\partial }_{0,t}^{\beta -\gamma }{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+\lambda {∥\nabla u∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2+{\partial }_{0,t}^{-\gamma }{∥u∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C{\partial }_{0,t}^{-\gamma }\left({∥f_2∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)$$

which also yields (7). □

2.3. Formulation of the Discrete Problem

To construct a numerical method, we first introduce a uniform partition of the time interval $\overline{J}$ by points $t_n=n\tau$, $\tau >0$, $n=0,1,\dots ,N$, such that $N\tau =T$. Let $p^n$ denote a semi-discrete approximation of the function p with respect to time at $t=t_n$.

To define a semi-discrete formulation of Problem 1, we use the following approximation formula for the fractional derivative in the sense of Caputo.

Lemma 3.

The discrete analog ${\mathrm{\Delta }}_{0,t}^{\nu }{p}^n$ of the fractional derivative in the sense of Caputo ${\partial }_{0,t}^{\nu }p\left(t_n\right)$ of order $0<\nu <1$ can be represented as [48]

$${\mathrm{\Delta }}_{0,t}^{\nu }{p}^n=\sum _{s=1}^n{\delta }_{n,s}^{\nu }\left(p^s-p^{s-1}\right),$$

where

$${\delta }_{n,s}^{\nu }=\frac{{\tau }^{-\nu }}{\mathrm{\Gamma }\left(2-\nu \right)}\left[{\left(n-s+1\right)}^{1-\nu }-{\left(n-s\right)}^{1-\nu }\right].$$

Moreover, the following estimate holds for $r_n^{\nu }={\partial }_{0,t}^{\nu }p\left(t_n\right)-{\mathrm{\Delta }}_{0,t}^{\nu }{p}^n$:

$$\left|r_n^{\nu }\right|\le \frac{5-\nu }{8\left(1-\nu \right)}\underset{0\le t\le t_n}{max}\left|{\partial }_t^{2}p\left(t\right)\right|{\tau }^{2-\nu }.$$

Lemma 4.

The following properties hold for the coefficients ${\delta }_{n,s}^{\nu }$ presented in Lemma 3:

(a)

${\delta }_{n,s}^{\nu }>0$, $s=1,2,\dots ,n$;

(b)

${\delta }_{n,s}^{\nu }<{\delta }_{n,s+1}^{\nu }$, $s=1,2,\dots ,n-1$;

(c)

${\delta }_{n,s}^{\nu }={\delta }_{n-1,s-1}^{\nu }$;

(d)

${\displaystyle \sum _{s=1}^n{\delta }_{n,s}^{\nu }=\frac{n{t}_n^{-\nu }}{\mathrm{\Gamma }\left(2-\nu \right)}}$.

Approximate the derivative of the function p at $t=t_n$ as follows:

$${\partial }_{t}p\left(t_n\right)=\left\{\begin{array}{cc}{\displaystyle \frac{3p^n-4p^{n-1}+p^{n-2}}{2\tau }+\frac{{\tau }^2}{3}{\partial }_t^3\left({\zeta }_n\right),}\hfill & n\ge 2,{\zeta }_n\in \left(t_{n-2},t_n\right),\hfill \\ {\displaystyle \frac{p^1-p^0}{\tau }+\frac{\tau }{2}{\partial }_t^{2}p\left({\zeta }_1\right),}\hfill & n=1,{\zeta }_1\in \left(t_0,t_1\right).\hfill \end{array}\right.$$

(13)

Let us define a semi-discrete formulation of Problem 1 as below:

Problem 3.

Let ${\left\{\left(u^i,p^i\right)\right\}}_{i=0}^{n-1}$, $\left(u^i,p^i\right)\in H_0^1\left(\mathrm{\Omega }\right)\times L^2\left(\mathrm{\Omega }\right)$ be given, in particular, $p^0=p_0\left(x\right)$ and $u^0=u_0\left(x\right)$. Find $\left(u^n,p^n\right)\in H_0^1\left(\mathrm{\Omega }\right)\times L^2\left(\mathrm{\Omega }\right)$ such that for all $w\in L^2\left(\mathrm{\Omega }\right)$ and $v\in H_0^1\left(\mathrm{\Omega }\right)$: when $n=1$:

$$\left(\frac{p^1-p^0}{\tau },w\right)=\left(u^1,w\right),$$

(14)

$$k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{u}^1,v\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{u}^1,v\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla u^1,\nabla v\right)$$

$$+\left(u^1,v\right)=\left(f_1\left(p^0\right),v\right)+\left(f_2^1,v\right),$$

(15)

when $n\ge 2$:

$${\displaystyle \left(\frac{3p^n-4p^{n-1}+p^{n-2}}{2\tau },w\right)=\left(u^n,w\right),}$$

$$k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{u}^n,v\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{u}^n,v\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla u^n,\nabla v\right)+\left(u^n,v\right)$$

(16)

$$=\left(2f_1\left(p^{n-1}\right)-f_1\left(p^{n-2}\right),v\right)+\left(f_2^n,v\right),$$

(17)

where $\alpha ,\beta ,\gamma \in \left(0,1\right)$.

Now, we introduce a quasi-uniform triangulation ${\mathcal{K}}_h$ in $\overline{\mathrm{\Omega }}$ with the mesh parameter h. For a non-negative integer l, let $P_l\left(\mathrm{\Omega }\right)$ denote the space of polynomials of degree at most l on $\overline{\mathrm{\Omega }}$. Define the finite element spaces

$$V_h=\left\{v_h\in H_0^1\left(\mathrm{\Omega }\right)\cap C^0\left(\overline{\mathrm{\Omega }}\right)|{\left.v_h\right|}_e\in P_1\left(e\right),\forall e\in {\mathcal{K}}_h\right\},$$

$$W_h=\left\{w_h\in L^2\left(\mathrm{\Omega }\right)|{\left.w_h\right|}_e\in P_0\left(e\right),\forall e\in {\mathcal{K}}_h\right\}.$$

Define the projection operator $Q_h:H_0^1\left(\mathrm{\Omega }\right)\to V_h$ satisfying

$$\left(\nabla \left(Q_{h}p-p\right),\nabla p_h\right)=0\forall p\in H_0^1\left(\mathrm{\Omega }\right),p_h\in V_h$$

and $L^2$-projection operator ${\mathrm{\Pi }}_h:L^2\left(\mathrm{\Omega }\right)\to W_h$ such that

$$\left({\mathrm{\Pi }}_{h}p-p,p_h\right)=0\forall p\in L^2\left(\mathrm{\Omega }\right),p_h\in W_h.$$

The projection operators have the following properties:

$${∥p-Q_{h}p∥}_{L^2\left(\mathrm{\Omega }\right)}+h{∥p-Q_{h}p∥}_{H^1\left(\mathrm{\Omega }\right)}\le C{h}^2{∥p∥}_{H^2\left(\mathrm{\Omega }\right)}\forall p\in H_0^1\left(\mathrm{\Omega }\right)\cap H^2\left(\mathrm{\Omega }\right),$$

$${∥{\mathrm{\Pi }}_{h}p∥}_{L^2\left(\mathrm{\Omega }\right)}\le C{∥p∥}_{L^2\left(\mathrm{\Omega }\right)},\forall p\in L^2\left(\mathrm{\Omega }\right),$$

$${∥p-{\mathrm{\Pi }}_{h}p∥}_{L^2\left(\mathrm{\Omega }\right)}\le Ch{∥p∥}_{H^1\left(\mathrm{\Omega }\right)},\forall p\in H^1\left(\mathrm{\Omega }\right).$$

Let us define the following finite element procedure:

Problem 4.

Let ${\left\{\left(u_h^i,p_h^i\right)\right\}}_{i=0}^{n-1}$, $\left(u_h^i,p_h^i\right)\in V_h\times W_h$ be given, in particular, $p_h^0$ and $u_h^0$ be the $L^2$-projections of $p_0$ and $u_0$. Find $\left(u_h^n,p_h^n\right)\in V_h\times W_h$ such that, for all $w_h\in W_h$ and $v_h\in V_h$: when $n=1$:

$${\displaystyle \left(\frac{p_h^1-p_h^0}{\tau },w_h\right)=\left(u_h^1,w_h\right),}$$

(18)

$$k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{u}_h^1,v_h\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{u}_h^1,v_h\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla u_h^1,\nabla v_h\right)$$

$$+\left(u_h^1,v_h\right)=\left(f_1\left(p_h^0\right),v_h\right)+\left(f_2^1,v_h\right),$$

(19)

when $n\ge 2$:

$${\displaystyle \left(\frac{3p_h^n-4p_h^{n-1}+p_h^{n-2}}{2\tau },w_h\right)=\left(u_h^n,w_h\right),}$$

(20)

$$k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{u}_h^n,v_h\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{u}_h^n,v_h\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla u_h^n,\nabla v_h\right)+\left(u_h^n,v_h\right)$$

$$=\left(2f_1\left(p_h^{n-1}\right)-f_1\left(p_h^{n-2}\right),v_h\right)+\left(f_2^n,v_h\right),$$

(21)

where $\alpha ,\beta ,\gamma \in \left(0,1\right)$.

2.4. Stability and Convergence of the Numerical Method

First, let us present a few auxiliary lemmas.

Lemma 5

([29]). The following inequality holds for ${\mathrm{\Delta }}_{0,t}^{\nu }{p}^n$, $n\ge 1$:

$$\left({\mathrm{\Delta }}_{0,t}^{\nu }{p}^n,p^n\right)\ge {\mathrm{\Theta }}_n^{\nu }-{\mathrm{\Theta }}_{n-1}^{\nu }-\frac{1}{2}{\delta }_{n,1}^{\nu }{∥p^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2,$$

where ${\displaystyle {\mathrm{\Theta }}_i^{\nu }=\frac{1}{2}\left({\delta }_{i,1}^{\nu }{∥p^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{\delta }_{i,2}^{\nu }{∥p^2∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\dots +{\delta }_{i,i}^{\nu }{∥p^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)}$, $i\ge 1$ and ${\mathrm{\Theta }}_0^{\nu }=0$.

Lemma 6

([29]). Let the sequence ${\left\{p^i\right\}}_{i=0}^n$, $p^i\in L^2\left(\mathrm{\Omega }\right)$ be given. Then,

$$\left(\frac{p^n-p^{n-1}}{\tau },p^n\right)=\frac{1}{2\tau }\left({∥p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2-{∥p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p^n-p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right),n\ge 1,$$

$$\left(\frac{3p^n-4p^{n-1}+p^{n-2}}{2\tau },p^n\right)=\frac{1}{4\tau }\left({∥p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2p^n-p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2-{∥p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right.$$

$$\left.-{∥2p^{n-1}-p^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p^n-2p^{n-1}+p^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right),n\ge 2.$$

In addition, we formulate a discrete analogue of Gronwall’s lemma, which will be used several times throughout the paper.

Lemma 7

([49]). If $\left\{a_n\right\}$ and $\left\{b_n\right\}$ are two positive sequences, $\left\{c_n\right\}$ is a monotone positive sequence, and they satisfy the inequalities

$$a_0+b_0\le c_0,a_n+b_n\le c_n+\lambda \sum _{i=0}^{n-1}{a}_i,\lambda >0;$$

then, the following estimate holds:

$$a_n+b_n\le c_{n}exp\left(n\lambda \right),n=0,1,2,\dots$$

Now, we prove the stability of the constructed numerical method.

Theorem 2.

The discrete scheme (18)–(21) is unconditionally stable with respect to the initial data and right-hand side of the equation for all $\tau <1/4$, and the following estimate holds:

$${∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+2\tau {∥u_h^n∥}_{H^1\left(\mathrm{\Omega }\right)}^2\le C\left({∥f_2∥}_{L^{\infty }\left(J;L^2\left(\mathrm{\Omega }\right)\right)}^2+{∥p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥u_h^0∥}_{H^1\left(\mathrm{\Omega }\right)}^2\right).$$

Proof. 

For $n=1$, let us choose $w_h=p_h^1$ in (18) and $v_h=u_h^1$ in (20).

$$\left(\frac{p_h^1-p_h^0}{\tau },p_h^1\right)=\left(u_h^1,p_h^1\right),$$

$$k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{u}_h^1,u_h^1\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{u}_h^1,u_h^1\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla u_h^1,\nabla u_h^1\right)+\left(u_h^1,u_h^1\right)=\left(f_1\left(p_h^0\right),u_h^1\right)+\left(f_2^1,u_h^1\right).$$

Making use of Lemmas 5 and 6, we obtain

$$\frac{1}{2\tau }\left({∥p_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2-{∥p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p_h^1-p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)\le {∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}{∥p_h^1∥}_{L^2\left(\mathrm{\Omega }\right)},$$

(22)

$$k_{\alpha }\left({\mathrm{\Theta }}_1^{\alpha }-{\mathrm{\Theta }}_0^{\alpha }-\frac{1}{2}{\delta }_{1,1}^{\alpha }{∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)+k_{\beta }\left({\mathrm{\Theta }}_1^{\beta }-{\mathrm{\Theta }}_0^{\beta }-\frac{1}{2}{\delta }_{1,1}^{\beta }{∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)$$

$$+\lambda \left({\mathrm{\Theta }}_1^{\gamma }-{\mathrm{\Theta }}_0^{\gamma }-\frac{1}{2}{\delta }_{1,1}^{\gamma }{∥\nabla u_h^0∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right)+{∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$\le \mu {∥p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}+{∥f_2^1∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)},$$

(23)

where

$${\mathrm{\Theta }}_n=k_{\alpha }{\mathrm{\Theta }}_n^{\alpha }+k_{\beta }{\mathrm{\Theta }}_n^{\beta }+\lambda {\mathrm{\Theta }}_n^{\gamma },$$

(24)

$${\mathrm{\Theta }}_n^{\nu }=\frac{1}{2}\sum _{s=1}^n{\delta }_{n,s}^{\nu }{∥u_h^s∥}_{L^2\left(\mathrm{\Omega }\right)}^2,\nu \in \left\{\alpha ,\beta \right\},{\mathrm{\Theta }}_n^{\gamma }=\frac{1}{2}\sum _{s=1}^n{\delta }_{n,s}^{\gamma }{∥\nabla u_h^s∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2.$$

(25)

Combining (22) and (23) and multiplying the resulting inequality by $2\tau$, we obtain:

$${∥p_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p_h^1-p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+2\tau {\mathrm{\Theta }}_1+2\tau {∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$\le 2\tau \mu {∥p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}+2\tau {∥f_2^1∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}+{∥p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+2\tau {∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}{∥p_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+2\tau {\mathrm{\Theta }}_0+\tau \left(k_{\alpha }{\delta }_{1,1}^{\alpha }+k_{\beta }{\delta }_{1,1}^{\beta }\right){∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau \lambda {\delta }_{1,1}^{\gamma }{∥\nabla u_h^0∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2.$$

(26)

By applying Cauchy and Young inequalities to the terms on the right-hand side of (26), we obtain that, for all $\tau$ satisfying $\tau <1/2$:

$${∥p_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p_h^1-p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+2\tau {\mathrm{\Theta }}_1+\tau {∥u_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$\le C\tau {∥f_2^1∥}_{L^2\left(\mathrm{\Omega }\right)}+{∥p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C{∥u_h^0∥}_{H^1\left(\mathrm{\Omega }\right)}^2,$$

(27)

where $C=Tmax\left\{k_{\alpha }{\delta }_{1,1}^{\alpha }+k_{\beta }{\delta }_{1,1}^{\beta },\lambda {\delta }_{1,1}^{\gamma }\right\}$.

For $n\ge 2$, let us choose $w_h=p_h^n$ in (19) and $v_h=u_h^n$ in (20):

$$\left(\frac{3p_h^n-4p_h^{n-1}+p_h^{n-2}}{2\tau },p_h^n\right)=\left(u_h^n,p_h^n\right),$$

$$k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{u}_h^n,u_h^n\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{u}_h^n,u_h^n\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla u_h^n,\nabla u_h^n\right)+\left(u_h^n,u_h^n\right)$$

$$=\left(2f_1\left(p_h^{n-1}\right)-f_1\left(p_h^{n-2}\right),u_h^n\right)+\left(f_2^n,u_h^n\right).$$

Making use of Lemmas 5 and 6, we obtain

$$\frac{1}{4\tau }\left({∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2p_h^n-p_h^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2-{∥p_h^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right.$$

$$\left.-{∥2p_h^{n-1}-p_h^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p_h^n-2p_h^{n-1}+p_h^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)\le {∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}{∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)},$$

(28)

$$k_{\alpha }\left({\mathrm{\Theta }}_n^{\alpha }-{\mathrm{\Theta }}_{n-1}^{\alpha }-\frac{1}{2}{\delta }_{n,1}^{\alpha }{∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)+k_{\beta }\left({\mathrm{\Theta }}_n^{\beta }-{\mathrm{\Theta }}_{n-1}^{\beta }-\frac{1}{2}{\delta }_{n,1}^{\beta }{∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)$$

$$+\lambda \left({\mathrm{\Theta }}_n^{\gamma }-{\mathrm{\Theta }}_{n-1}^{\gamma }-\frac{1}{2}{\delta }_{n,1}^{\gamma }{∥\nabla u_h^0∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right)+{∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$\le {∥2f_1\left(p_h^{n-1}\right)-f_1\left(p_h^{n-2}\right)∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}+{∥f_2∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}.$$

(29)

Combining (28) and (29) and multiplying the resulting inequality by $4\tau$, we obtain:

$${∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2p_h^n-p_h^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p_h^n-2p_h^{n-1}+p_h^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {\mathrm{\Theta }}_n+4\tau {∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$\le 4\tau {∥2f_1\left(p_h^{n-1}\right)-f_1\left(p_h^{n-2}\right)∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}+4\tau {∥f_2∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+{∥p_h^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2p_h^{n-1}-p_h^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}{∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+4\tau {\mathrm{\Theta }}_{n-1}+2\tau \left(k_{\alpha }{\delta }_{n,1}^{\alpha }+k_{\beta }{\delta }_{n,1}^{\beta }\right){∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+2\tau \lambda {\delta }_{n,1}^{\gamma }{∥\nabla u_h^0∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2.$$

Sum the last inequality with respect to n from 2 to n to obtain

$${∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {\mathrm{\Theta }}_n+4\tau {∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le 4\tau {∥2f_1\left(p_h^{n-1}\right)-f_1\left(p_h^{n-2}\right)∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+4\tau \sum _{i=2}^{n-1}{∥2f_1\left(p_h^{i-1}\right)-f_1\left(p_h^{i-2}\right)∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^i∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+4\tau {∥f_2^n∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}+4\tau \sum _{i=2}^{n-1}{∥f_2^i∥}_{L^2\left(\mathrm{\Omega }\right)}{∥u_h^i∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+{∥p_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2p_h^1-p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}{∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}+4\tau \sum _{i=2}^{n-1}{∥u_h^i∥}_{L^2\left(\mathrm{\Omega }\right)}{∥p_h^i∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+4\tau {\mathrm{\Theta }}_1+2\tau \sum _{i=2}^n\left(k_{\alpha }{\delta }_{i,1}^{\alpha }+k_{\beta }{\delta }_{i,1}^{\beta }\right){∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+2\tau \lambda \sum _{i=2}^n{\delta }_{i,1}^{\gamma }{∥\nabla u_h^0∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2.$$

By applying Cauchy and Young inequalities to the first term on the right-hand side, we obtain that for all $\tau$ such that $\tau <1/4$:

$${∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {\mathrm{\Theta }}_n+\tau {∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\tau \sum _{i=0}^{n-1}{∥p_h^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C\tau \sum _{i=2}^n{∥f_2^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+C\tau \sum _{i=2}^{n-1}\left({∥u_h^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p_h^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right)+{∥p_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2p_h^1-p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+4\tau {\mathrm{\Theta }}_1+2\tau \sum _{i=2}^n\left(k_{\alpha }{\delta }_{i,1}^{\alpha }+k_{\beta }{\delta }_{i,1}^{\beta }\right){∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+2\tau \lambda \sum _{i=2}^n{\delta }_{i,1}^{\gamma }{∥\nabla u_h^0∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2.$$

Making use of the discrete Gronwall’s lemma, we obtain

$${∥p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {\mathrm{\Theta }}_n+\tau {∥u_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\left(\tau \sum _{i=2}^n{∥f_2^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥p_h^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2p_h^1-p_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2\right.$$

$$\left.+2\tau {\mathrm{\Theta }}_1+2\tau \sum _{i=2}^n\left(k_{\alpha }{\delta }_{i,1}^{\alpha }+k_{\beta }{\delta }_{i,1}^{\beta }\right){∥u_h^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+2\tau \lambda \sum _{i=2}^n{\delta }_{i,1}^{\gamma }{∥\nabla u_h^0∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right).$$

(30)

Combining (27) and (30), we arrive at the assertion of the theorem. □

Now, we prove the main results of the section:

Theorem 3.

Let ${\left\{\left(p^i,u^i\right)\right\}}_{i=0}^N$, $\left(p^i,u^i\right)\in H_0^1\left(\mathrm{\Omega }\right)\times L^2\left(\mathrm{\Omega }\right)$ be the solution of Problem 3. Then, the following inequality holds for $\alpha ,\beta ,\gamma \in \left(0,1\right)$, for all $\left(p^n,u^n\right)\in H_0^1\left(\mathrm{\Omega }\right)\times L^2\left(\mathrm{\Omega }\right)$ and $\tau <1/4$:

$${∥p\left(t_n\right)-p^n∥}_{L^2\left(\mathrm{\Omega }\right)}\le C\left({\tau }^{2-\alpha }+{\tau }^{2-\beta }+{\tau }^{2-\gamma }\right),$$

(31)

$${∥u\left(t_n\right)-u^n∥}_{H^1\left(\mathrm{\Omega }\right)}\le C\tau .$$

(32)

Proof. 

Denote ${\sigma }^n=u\left(t_n\right)-u^n$, ${\pi }^n=p\left(t_n\right)-p^n$. Substract (14)–(15) and (16)–(17) from (5)–(6) at $t=t_n$, take $w={\pi }^n$, $v={\sigma }^n$ and use (6) and (13) to obtain: when $n=1$:

$$\left(\frac{{\pi }^1-{\pi }^0}{\tau },{\pi }^1\right)+\frac{\tau }{2}\left({\partial }_t^{2}p\left({\zeta }_1\right),{\pi }^1\right)=\left({\sigma }^1,{\pi }^1\right),$$

$$k_{\alpha }\left(r_n^{\alpha },{\sigma }^1\right)+k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{\sigma }^1,{\sigma }^1\right)+k_{\beta }\left(r_n^{\beta },{\sigma }^1\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{\sigma }^1,{\sigma }^1\right)$$

$$+\lambda \left(r_n^{\gamma },\nabla {\sigma }^1\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla {\sigma }^1,\nabla {\sigma }^1\right)+\left({\sigma }^1,{\sigma }^1\right)=\left(f_1\left(p\left(t_0\right)\right)-f_1\left(p^0\right),{\sigma }^1\right),$$

when $n\ge 2$:

$$\left(\frac{3{\pi }^n-4{\pi }^{n-1}+{\pi }^{n-2}}{2\tau },{\pi }^n\right)+\frac{{\tau }^2}{3}\left({\partial }_t^{3}p\left({\zeta }_n\right),{\pi }^n\right)=\left({\sigma }^n,{\pi }^n\right),$$

$$k_{\alpha }\left(r_n^{\alpha },{\sigma }^n\right)+k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{\sigma }^n,{\sigma }^n\right)+k_{\beta }\left(r_n^{\beta },{\sigma }^n\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{\sigma }^n,{\sigma }^n\right)$$

$$+\lambda \left(r_n^{\gamma },\nabla {\sigma }^n\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla {\sigma }^n,\nabla {\sigma }^n\right)+\left({\sigma }^n,{\sigma }^n\right)$$

$$=2\left(f_1\left(p\left(t_{n-1}\right)\right)-f_1\left(p^{n-1}\right),{\sigma }^n\right)-\left(f_1\left(p\left(t_{n-2}\right)\right)-f_1\left(p^{n-2}\right),{\sigma }^n\right),$$

where $t_{i-1}\le {\zeta }_i\le t_i$.

By using Lemma 5, we have:

when $n=1$:

$${\displaystyle {∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le 2\tau {∥{\sigma }^1∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}+{\tau }^2{∥{\partial }_t^{2}p\left({\zeta }_1\right)∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)},}$$

(33)

$${\mathrm{\Theta }}_1+{∥{\sigma }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le {\mathrm{\Theta }}_0+k_{\alpha }{∥r_1^{\alpha }∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\sigma }^1∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+k_{\beta }{∥r_1^{\beta }∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\sigma }^1∥}_{L^2\left(\mathrm{\Omega }\right)}+\lambda {∥r_1^{\gamma }∥}_{L^2\left(\mathrm{\Omega }\right)}{∥\nabla {\sigma }^1∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}+C{∥{\pi }^0∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\sigma }^1∥}_{L^2\left(\mathrm{\Omega }\right)},$$

when $n\ge 2$:

$${\displaystyle {∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2{\pi }^n-{\pi }^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2-{∥{\pi }^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2-{∥2{\pi }^{n-1}-{\pi }^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2}$$

$$+{∥{\pi }^n-2{\pi }^{n-1}+{\pi }^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\frac{4{\tau }^3}{3}{∥{\partial }_t^{3}p\left({\zeta }_n\right)∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$\le 4\tau {∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)},$$

(34)

$${\mathrm{\Theta }}_n+{∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le {\mathrm{\Theta }}_{n-1}+k_{\alpha }{∥r_n^{\alpha }∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}$$

$$+k_{\beta }{∥r_n^{\beta }∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}+\lambda {∥r_n^{\gamma }∥}_{L^2\left(\mathrm{\Omega }\right)}{∥\nabla {\sigma }^n∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}$$

$$+C{∥{\pi }^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}+C{∥{\pi }^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)},$$

(35)

where

$${\mathrm{\Theta }}_n=k_{\alpha }{\mathrm{\Theta }}_n^{\alpha }+k_{\beta }{\mathrm{\Theta }}_n^{\beta }+\lambda {\mathrm{\Theta }}_n^{\gamma },$$

$$\begin{array}{ccc}{\displaystyle {\mathrm{\Theta }}_n^{\alpha }=\frac{1}{2}\sum _{s=1}^n{\delta }_{n,s}^{\alpha }{∥{\sigma }^s∥}_{L^2\left(\mathrm{\Omega }\right)}^2},\hfill & & r_n^{\alpha }={\partial }_{0,t}^{\alpha }u\left(t_n\right)-{\mathrm{\Delta }}_{0,t}^{\alpha }u\left(t_n\right),\hfill \\ {\displaystyle {\mathrm{\Theta }}_n^{\beta }=\frac{1}{2}\sum _{s=1}^n{\delta }_{n,s}^{\beta }{∥{\sigma }^s∥}_{L^2\left(\mathrm{\Omega }\right)}^2},\hfill & & r_n^{\beta }={\partial }_{0,t}^{\beta }u\left(t_n\right)-{\mathrm{\Delta }}_{0,t}^{\beta }u\left(t_n\right),\hfill \\ {\displaystyle {\mathrm{\Theta }}_n^{\gamma }=\frac{1}{2}\sum _{s=1}^n{\delta }_{n,s}^{\gamma }{∥\nabla {\sigma }^s∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2},\hfill & & r_n^{\gamma }={\partial }_{0,t}^{\gamma }\nabla u\left(t_n\right)-{\mathrm{\Delta }}_{0,t}^{\gamma }\nabla u\left(t_n\right).\hfill \end{array}$$

By adding (34) and (35) multiplied by $\tau$, we obtain:

$${∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2{\pi }^n-{\pi }^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥{\pi }^n-2{\pi }^{n-1}+{\pi }^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n+\frac{\tau }{2}{∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$\le {∥{\pi }^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2{\pi }^{n-1}-{\pi }^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_{n-1}+2\tau {∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+C\tau \left({∥r_n^{\alpha }∥}_{L^2\left(\mathrm{\Omega }\right)}+{∥r_n^{\beta }∥}_{L^2\left(\mathrm{\Omega }\right)}\right){∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}+\lambda \tau {∥r_n^{\gamma }∥}_{L^2\left(\mathrm{\Omega }\right)}{∥\nabla {\sigma }^n∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}$$

$$+\frac{4{\tau }^3}{3}{∥{\partial }_t^{3}p\left({\zeta }_n\right)∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}+C\tau {∥{\pi }^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C\tau {∥{\pi }^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2.$$

(36)

Sum (36) with respect to n from 2 to n to obtain

$${∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n+\frac{\tau }{2}\sum _{i=2}^n{∥{\sigma }^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le 5{∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_1$$

$$+C\tau {∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C\tau \sum _{i=0}^{n-1}{∥{\pi }^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\frac{\tau }{4}\left(\left(k_{\alpha }+k_{\beta }\right){∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\lambda {∥\nabla {\sigma }^n∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right)$$

$$+\frac{{\tau }^2}{2}\sum _{i=2}^{n-1}\left(\left(k_{\alpha }+k_{\beta }\right){∥{\sigma }^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\lambda {∥\nabla {\sigma }^i∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right)+C\left({\tau }^{4-2\alpha }+{\tau }^{4-2\beta }+{\tau }^{4-2\gamma }\right),$$

and notice that ${\displaystyle \frac{1}{2}\left(k_{\alpha }+k_{\beta }\right){∥{\sigma }^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\frac{\lambda }{2}{∥\nabla {\sigma }^i∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\le {\mathrm{\Theta }}_i}$ and choose ${\tau }_1$ such that the condition $1-C{\tau }_1>0$ is satisfied, we conclude that for all $\tau <{\tau }_1$:

$${∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n+\tau {∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le 10{∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+2\tau {\mathrm{\Theta }}_1+6\tau \sum _{i=2}^{n-1}{∥{\pi }^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+6{\tau }^2\sum _{i=2}^{n-1}{\mathrm{\Theta }}_i+C\left({\tau }^{4-2\alpha }+{\tau }^{4-2\beta }+{\tau }^{4-2\gamma }\right).$$

Making use of Lemma 7, we have

$${∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n+\tau {∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\left({∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_1+{\tau }^{4-2\alpha }+{\tau }^{4-2\beta }+{\tau }^{4-2\gamma }\right).$$

(37)

Summing (33) and (35) multiplied by $\tau$ and applying Young inequality, we arrive at

$${∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_1\le C\left({\tau }^4+{\tau }^{5-2\alpha }+{\tau }^{5-2\beta }+{\tau }^{5-2\gamma }\right).$$

(38)

Combining (37) and (38), we obtain

$${∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n\le C\left({\tau }^{4-2\alpha }+{\tau }^{4-2\beta }+{\tau }^{4-2\gamma }\right)$$

which implies the inequality (31).

To obtain the inequality (32), we derive from (36) the following inequality:

$${∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n+\frac{\tau }{2}\sum _{i=2}^n{∥{\sigma }^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le 5{∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_1+3\tau {∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+3\tau \sum _{i=0}^{n-1}{∥{\pi }^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\frac{\tau }{4}\left(\left(k_{\alpha }+k_{\beta }\right){∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\lambda {∥\nabla {\sigma }^n∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right)$$

$$+\frac{1}{2}{\tau }^{2-max\left\{\alpha ,\beta ,\gamma \right\}}\sum _{i=2}^{n-1}\left(\left(k_{\alpha }+k_{\beta }\right){∥{\sigma }^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\lambda {∥\nabla {\sigma }^i∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right)+C\left({\tau }^{4-\alpha }+{\tau }^{4-\beta }+{\tau }^{4-\gamma }\right).$$

Utilizing the technique as in obtaining inequality (37) and applying the Gronwall lemma, we arrive at

$${∥{\pi }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n+\tau {∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\left({∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_1+{\tau }^{4-\alpha }+{\tau }^{4-\beta }+{\tau }^{4-\gamma }\right).$$

(39)

Furthermore, considering the case $n=1$ in the same manner, we obtain the following estimate instead of (38):

$${∥{\pi }^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_1\le C{\tau }^4.$$

(40)

By combining (39) and (40), and using the definition of ${\mathrm{\Theta }}_n$, we arrive at the inequality

$${∥{\pi }^n∥}_0^2+\sum _{s=1}^n\left(\frac{{\tau }^{1-\alpha }}{2}{∥{\sigma }^s∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\frac{{\tau }^{1-\beta }}{2}{∥{\sigma }^s∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\frac{{\tau }^{1-\gamma }}{2}{∥\nabla {\sigma }^s∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right)$$

$$\le C\left({\tau }^{4-\alpha }+{\tau }^{4-\beta }+{\tau }^{4-\gamma }\right).$$

(41)

By taking the square root of both sides of (41), utilizing the elementary inequality ${\displaystyle {\left(\sum _{i=1}^m{a}_i^2\right)}^{1/2}\ge \frac{1}{\sqrt{m}}\sum _{i=1}^m\left|a_i\right|}$ and using the fact that ${\displaystyle \frac{1}{\sqrt{n}}\ge \sqrt{\frac{\tau }{T}}}$, we obtain

$${\tau }^{1-\frac{1}{2}max\left\{\alpha ,\beta ,\gamma \right\}}\left({∥{\sigma }^n∥}_{L^2\left(\mathrm{\Omega }\right)}+{∥\nabla {\sigma }^n∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}\right)\le C{\tau }^{2-\frac{1}{2}max\left\{\alpha ,\beta ,\gamma \right\}},$$

which yields the inequality (32). □

Theorem 4.

Let ${\left\{\left(p_h^i,u_h^i\right)\right\}}_{i=0}^N$, $p_h^i\in W_h$, $u_h^i\in V_h$ be the solution of Problem 4. Then, the following inequality holds for $\alpha ,\beta ,\gamma \in \left(0,1\right)$, for all $\left(p_h^n,u_h^n\right)\in W_h\times V_h$ and sufficiently small τ:

$${∥p\left(t_n\right)-p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}\le C\left(h^2+{\tau }^{1/2}{h}^2+{\tau }^{2-\alpha }+{\tau }^{2-\beta }+{\tau }^{2-\gamma }\right),$$

(42)

$${∥u\left(t_n\right)-u_h^n∥}_{H^1\left(\mathrm{\Omega }\right)}\le C\left(\tau +{\tau }^{-1}{h}^2\right).$$

(43)

Proof. 

Denote

$$p^n-p_h^n=\left(p^n-{\mathrm{\Pi }}_h{p}^n\right)+\left({\mathrm{\Pi }}_h{p}^n-p_h^n\right)={\psi }_p^n+{\xi }_p^n,$$

$$u^n-u_h^n=\left(u^n-Q_h{u}^n\right)+\left(Q_h{u}^n-u_h^n\right)={\psi }_u^n+{\xi }_u^n.$$

Consider the difference of Equations (14)– (17) and (18)–(21):

when $n=1$:

$$\left(\frac{{\xi }_p^1-{\xi }_p^0}{\tau },w_h\right)=\left({\psi }_u^1+{\xi }_u^1,w_h\right),$$

(44)

$$k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }\left({\psi }_u^1+{\xi }_u^1\right),v_h\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }\left({\psi }_u^1+{\xi }_u^1\right),v_h\right)$$

$$+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla {\xi }_u^1,\nabla v_h\right)+\left({\psi }_u^1+{\xi }_u^1,v_h\right)=\left(f_1\left(p^0\right)-f_1\left(p_h^0\right),{\sigma }^1\right),$$

when $n\ge 2$:

$$\left(\frac{3{\xi }_p^n-4{\xi }_p^{n-1}+{\xi }_p^{n-2}}{2\tau },w_h\right)=\left({\psi }_u^n+{\xi }_u^n,w_h\right),$$

(45)

$$k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }\left({\psi }_u^n+{\xi }_u^n\right),v_h\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }\left({\psi }_u^n+{\xi }_u^n\right),v_h\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla {\xi }_u^n,\nabla v_h\right)$$

$$+\left({\psi }_u^n+{\xi }_u^n,v_h\right)=2\left(f_1\left(p^{n-1}\right)-f_1\left(p_h^{n-1}\right),v_h\right)-\left(f_1\left(p^{n-2}\right)-f_1\left(p_h^{n-2}\right),v_h\right).$$

(46)

Choose $\left(w_h,v_h\right)=\left({\xi }_p^n,{\xi }_u^n\right)$ in (45), (46) and add the resulting equations multiplied by $4\tau$:

$${∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2{\xi }_p^n-{\xi }_p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2-{∥{\xi }_p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2-{∥2{\xi }_p^{n-1}-{\xi }_p^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+{∥{\xi }_p^n-2{\xi }_p^{n-1}+{\xi }_p^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {∥{\xi }_u^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+4\tau \left[k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{\xi }_u^n,{\xi }_u^n\right)+k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{\xi }_u^n,{\xi }_u^n\right)+\lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla {\xi }_u^n,\nabla {\xi }_u^n\right)\right]$$

$$+4\tau k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{\psi }_u^n,{\xi }_u^n\right)+4\tau k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{\psi }_u^n,{\xi }_u^n\right)+4\tau \left[\left({\psi }_u^n,{\xi }_u^n\right)-\left({\psi }_u^n,{\xi }_p^n\right)-\left({\xi }_u^n,{\xi }_p^n\right)\right]$$

$$\le C\tau {∥{\psi }_p^{n-1}+{\xi }_p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\xi }_u^n∥}_{L^2\left(\mathrm{\Omega }\right)}+C\tau {∥{\psi }_p^{n-2}+{\xi }_p^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}{∥{\xi }_u^n∥}_{L^2\left(\mathrm{\Omega }\right)}.$$

(47)

Applying Lemma 6, we obtain:

$$4\tau k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{\xi }_u^n,{\xi }_u^n\right)+4\tau k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{\xi }_u^n,{\xi }_u^n\right)+4\tau \lambda \left({\mathrm{\Delta }}_{0,t}^{\gamma }\nabla {\xi }_u^n,\nabla {\xi }_u^n\right)$$

$$\ge 4\tau \left({\mathrm{\Theta }}_n-{\mathrm{\Theta }}_{n-1}\right)-2\tau \left[k_{\alpha }{\delta }_{n,1}^{\alpha }{∥{\xi }_u^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+k_{\beta }{\delta }_{n,1}^{\beta }{∥{\xi }_u^0∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\lambda {\delta }_{n,1}^{\gamma }{∥\nabla {\xi }_u^0∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2\right],$$

(48)

where ${\mathrm{\Theta }}_n=k_{\alpha }{\mathrm{\Theta }}_n^{\alpha }+k_{\beta }{\mathrm{\Theta }}_n^{\beta }+\lambda {\mathrm{\Theta }}_n^{\gamma }$ with

$${\displaystyle {\mathrm{\Theta }}_n^{\overline{\nu }}=\frac{1}{2}\sum _{s=1}^n{\delta }_{n,s}^{\overline{\nu }}{∥{\xi }_u^s∥}_{L^2\left(\mathrm{\Omega }\right)}^2},\overline{\nu }\in \left\{\alpha ,\beta \right\};{\mathrm{\Theta }}_n^{\gamma }=\frac{1}{2}\sum _{s=1}^n{\delta }_{n,s}^{\gamma }{∥\nabla {\xi }_u^s∥}_{{\mathbf{L}}^2\left(\mathrm{\Omega }\right)}^2.$$

Furthermore, it follows from (47) that

$$4\tau k_{\alpha }\left({\mathrm{\Delta }}_{0,t}^{\alpha }{\psi }_u^n,{\xi }_u^n\right)+4\tau k_{\beta }\left({\mathrm{\Delta }}_{0,t}^{\beta }{\psi }_u^n,{\xi }_u^n\right)$$

$$\le 4\tau \left({\left(\frac{k_{\alpha }{T}^{1-\alpha }}{\mathrm{\Gamma }\left(2-\alpha \right)}\right)}^2+{\left(\frac{k_{\beta }{T}^{1-\beta }}{\mathrm{\Gamma }\left(2-\beta \right)}\right)}^2\right)\underset{1\le s\le n}{max}{∥{\partial }_t{\psi }_u\left({\zeta }_s\right)∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {∥{\xi }_u^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2.$$

To evaluate the remaining terms in the left-hand side of (47), we use Cauchy and Young inequalities:

$$4\tau \left({\psi }_u^n,{\xi }_u^n\right)-4\tau \left({\psi }_u^n,{\xi }_p^n\right)-4\tau \left({\xi }_u^n,{\xi }_p^n\right)$$

$$\le C\tau h^2{∥{\xi }_u^n∥}_{L^2\left(\mathrm{\Omega }\right)}+C\tau h^2{∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}+\tau {∥{\xi }_u^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C\tau {∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2.$$

(49)

Then, taking into account the inequalities (48)–(49), we obtain from (47) that

$${∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2{\xi }_p^n-{\xi }_p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {\mathrm{\Theta }}_n\le {∥{\xi }_p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+{∥2{\xi }_p^{n-1}-{\xi }_p^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+4\tau {\mathrm{\Theta }}_{n-1}+C\tau h^2{∥{\xi }_u^n∥}_{L^2\left(\mathrm{\Omega }\right)}+C\tau h^2{∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}+C\tau {∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2$$

$$+C\tau {∥{\xi }_p^{n-1}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C\tau {∥{\xi }_p^{n-2}∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C\tau h^4.$$

(50)

Sum the inequality (50) with respect to n from 2 to n, apply Cauchy and Young inequalities and notice that ${∥{\xi }_u^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2\le 2{\mathrm{\Theta }}_i$ to obtain

$$\left(1-C\tau \right){∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\left(4-C\tau \right)\tau {\mathrm{\Theta }}_n\le 5{∥{\xi }_p^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+4\tau {\mathrm{\Theta }}_1$$

$$+C{\tau }^2\sum _{i=2}^{n-1}{∥{\xi }_u^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C\tau \sum _{i=2}^{n-1}{∥{\xi }_p^i∥}_{L^2\left(\mathrm{\Omega }\right)}^2+C{h}^4+C\tau h^4.$$

Hence, for sufficiently small $\tau$, applying Lemma 7, we obtain

$${∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n\le C\left({∥{\xi }_p^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_1+h^4+\tau h^4\right).$$

(51)

Applying the same technique of estimation of terms to (44), (46) as in obtaining the inequality (51), we obtain

$${∥{\xi }_p^1∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_1\le C\left(h^4+\tau h^4\right).$$

(52)

Combining (51) and (52), we obtain

$${∥{\xi }_p^n∥}_{L^2\left(\mathrm{\Omega }\right)}^2+\tau {\mathrm{\Theta }}_n\le C\left(h^4+\tau h^4\right).$$

(53)

By taking into account the results of Theorem 3, we arrive at the inequality (42).

Finally, combining (53) with the estimate (32) in Theorem 3, we obtain the inequality

$${∥u\left(t_n\right)-u_h^n∥}_{H^1\left(\mathrm{\Omega }\right)}\le C\left({\tau }^{-1}{h}^2+{\tau }^{-1/2}{h}^2+\tau \right),$$

which implies the inequality (43). The theorem is proved. □

3. Results

To verify the theoretical convergence estimates obtained in Theorem 4 for the fully discrete scheme, a number of computational experiments was carried out on the example of a model problem for the equation describing fluid flow in fractured porous media. The implementation of the algorithms is carried out using FreeFEM++.

The purpose of the computational experiments carried out for the example below is to determine the dependence of actual convergence order on the orders of fractional derivatives, and compare it with the theoretical convergence order obtained in Theorem 4 in case of $\alpha ,\beta ,\gamma \in \left(0,1\right)$. For convenience, let us denote $\overline{\alpha }=\alpha +1$, $\overline{\beta }=\beta +1$, $\overline{\gamma }=\gamma +1$.

Example 1.

In $Q_T=\overline{\mathrm{\Omega }}\times \overline{J}$, where $\mathrm{\Omega }={\left(0,1\right)}^2$, $J=\left(0,1\right]$, consider the model of fluid flow in fractured porous media consisting of the fractional differential equation

$${\partial }_{t}p+{\partial }_{0,t}^{\overline{\alpha }}p+{\partial }_{0,t}^{\overline{\beta }}p-{\partial }_{0,t}^{\overline{\gamma }}{\nabla }^{2}p=p^2-\frac{12t^3{t}^{-\overline{\gamma }}\left(x_1^2-x_1+x_2^2-x_2\right)}{\mathrm{\Gamma }\left(2-\overline{\gamma }\right)\left({\overline{\gamma }}^2-5\overline{\gamma }+6\right)}$$

$$+\frac{6t^{3-\overline{\alpha }}{x}_1{x}_2\left(1-x_1\right)\left(1-x_2\right)}{\mathrm{\Gamma }\left(2-\overline{\alpha }\right)\left(2-\overline{\alpha }\right)\left(3-\overline{\alpha }\right)}+\frac{6t^{3-\overline{\beta }}{x}_1{x}_2\left(1-x_1\right)\left(1-x_2\right)}{\mathrm{\Gamma }\left(2-\overline{\beta }\right)\left(2-\overline{\beta }\right)\left(3-\overline{\beta }\right)}$$

$$-t^6{x}_1^2{x}_2^2{\left(1-x_1\right)}^2{\left(1-x_2\right)}^2+3t^2{x}_1{x}_2\left(1-x_1\right)\left(1-x_2\right),\left(x,t\right)\in \mathrm{\Omega }\times J,$$

subject to the following initial and boundary conditions:

$$p\left(x,0\right)=0,x\in \overline{\mathrm{\Omega }},$$

$$p\left(x,t\right)=0,x\in \partial \mathrm{\Omega },t\in J,$$

where $\overline{\alpha },\overline{\beta },\overline{\gamma }\in \left(1,2\right)$, $x=\left(x_1,x_2\right)$.

The exact solution of the problem is $p\left(x,t\right)=t^3{x}_1{x}_2\left(1-x_1\right)\left(1-x_2\right)$.

By assuming $h=C\tau$ for some positive real number C, Theorem 4 implies the temporal convergence orders

$$\begin{array}{cc}\hfill E_N^1& \equiv \underset{n\in \left\{0,1,\dots ,N\right\}}{max}{∥p\left(t_n\right)-p_h^n∥}_{L^2\left(\mathrm{\Omega }\right)}=O\left({\tau }^{3-\nu }\right),\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill E_N^2& \equiv \underset{n\in \left\{0,1,\dots ,N\right\}}{max}{∥u\left(t_n\right)-u_h^n∥}_{H^1\left(\mathrm{\Omega }\right)}=O\left(\tau \right)\hfill \end{array}$$

(55)

where $\nu =max\left\{\overline{\alpha },\overline{\beta },\overline{\gamma }\right\}$. Similar considerations are valid for the spatial convergence orders in the case when $\tau =C_{1}h$, $C_1>0$. Therefore, we considered a set of time steps in the range from ${\tau }_1=1/7$ to ${\tau }_7=1/448$ showed in the first column of

Table 1. $L^2$-errors, $E_{N_i}^1$, $i\ge 1$ calculated by (54) and convergence orders, $R_i^1$, $i\ge 2$, calculated by (56) for different values of ${\tau }_i\approx h_i$, $N_i=1/{\tau }_i$. The numbers in parentheses indicate the theoretical convergence order.

$\mathit{\tau }\approx \mathit{h}$	$\overline{\mathit{\alpha }}=\overline{\mathit{\beta }}=1.1$		$\overline{\mathit{\alpha }}=\overline{\mathit{\beta }}=1.5$		$\overline{\mathit{\alpha }}=\overline{\mathit{\beta }}=1.9$
	${\mathit{L}}^2$-Error	Order	${\mathit{L}}^2$-Error	Order	${\mathit{L}}^2$-Error	Order
$\overline{\gamma }=1.1$
$1/7$	$7.1923\times {10}^{-4}$	-	$8.4733\times {10}^{-4}$	-	$1.4621\times {10}^{-3}$	-
$1/14$	$1.8335\times {10}^{-4}$	1.97 $\left(\simeq 1.90\right)$	$2.4162\times {10}^{-4}$	1.81 $\left(\simeq 1.50\right)$	$5.6768\times {10}^{-4}$	1.36 $\left(\simeq 1.10\right)$
$1/28$	$4.7978\times {10}^{-5}$	1.93 $\left(\simeq 1.90\right)$	$6.9894\times {10}^{-5}$	1.79 $\left(\simeq 1.50\right)$	$2.2887\times {10}^{-4}$	1.31 $\left(\simeq 1.10\right)$
$1/56$	$1.2550\times {10}^{-5}$	1.93 $\left(\simeq 1.90\right)$	$2.0544\times {10}^{-5}$	1.77 $\left(\simeq 1.50\right)$	$9.6903\times {10}^{-5}$	1.24 $\left(\simeq 1.10\right)$
$1/112$	$3.2735\times {10}^{-6}$	1.94 $\left(\simeq 1.90\right)$	$6.1861\times {10}^{-6}$	1.73 $\left(\simeq 1.50\right)$	$4.3006\times {10}^{-5}$	1.17 $\left(\simeq 1.10\right)$
$1/224$	$8.5643\times {10}^{-7}$	1.93 $\left(\simeq 1.90\right)$	$1.9090\times {10}^{-6}$	1.70 $\left(\simeq 1.50\right)$	$1.9399\times {10}^{-5}$	1.15 $\left(\simeq 1.10\right)$
$1/448$	$2.2434\times {10}^{-7}$	1.93 $\left(\simeq 1.90\right)$	$6.0342\times {10}^{-7}$	1.66 $\left(\simeq 1.50\right)$	$8.8776\times {10}^{-6}$	1.13 $\left(\simeq 1.10\right)$
$\overline{\gamma }=1.5$
$1/7$	$1.8809\times {10}^{-3}$	-	$2.0215\times {10}^{-3}$	-	$2.5844\times {10}^{-3}$	-
$1/14$	$6.5152\times {10}^{-4}$	1.53 $\left(\simeq 1.50\right)$	$6.9763\times {10}^{-4}$	1.53 $\left(\simeq 1.50\right)$	$9.6132\times {10}^{-4}$	1.43 $\left(\simeq 1.10\right)$
$1/28$	$2.2333\times {10}^{-4}$	1.54 $\left(\simeq 1.50\right)$	$2.3892\times {10}^{-4}$	1.55 $\left(\simeq 1.50\right)$	$3.6367\times {10}^{-4}$	1.40 $\left(\simeq 1.10\right)$
$1/56$	$7.7119\times {10}^{-5}$	1.53 $\left(\simeq 1.50\right)$	$8.2526\times {10}^{-5}$	1.53 $\left(\simeq 1.50\right)$	$1.4199\times {10}^{-4}$	1.36 $\left(\simeq 1.10\right)$
$1/112$	$2.7128\times {10}^{-5}$	1.51 $\left(\simeq 1.50\right)$	$2.9060\times {10}^{-5}$	1.51 $\left(\simeq 1.50\right)$	$5.7823\times {10}^{-5}$	1.30 $\left(\simeq 1.10\right)$
$1/224$	$9.5040\times {10}^{-6}$	1.51 $\left(\simeq 1.50\right)$	$1.0187\times {10}^{-5}$	1.51 $\left(\simeq 1.50\right)$	$2.3874\times {10}^{-5}$	1.28 $\left(\simeq 1.10\right)$
$1/448$	$3.3411\times {10}^{-6}$	1.51 $\left(\simeq 1.50\right)$	$3.5833\times {10}^{-6}$	1.51 $\left(\simeq 1.50\right)$	$1.0074\times {10}^{-5}$	1.24 $\left(\simeq 1.10\right)$
$\overline{\gamma }=1.9$
$1/7$	$6.3612\times {10}^{-3}$	-	$6.3824\times {10}^{-3}$	-	$6.7604\times {10}^{-3}$	-
$1/14$	$2.7133\times {10}^{-3}$	1.23 $\left(\simeq 1.10\right)$	$2.7023\times {10}^{-3}$	1.24 $\left(\simeq 1.10\right)$	$2.8581\times {10}^{-3}$	1.24 $\left(\simeq 1.10\right)$
$1/28$	$1.1950\times {10}^{-3}$	1.18 $\left(\simeq 1.10\right)$	$1.1844\times {10}^{-3}$	1.19 $\left(\simeq 1.10\right)$	$1.2537\times {10}^{-3}$	1.19 $\left(\simeq 1.10\right)$
$1/56$	$5.4010\times {10}^{-4}$	1.15 $\left(\simeq 1.10\right)$	$5.3351\times {10}^{-4}$	1.15 $\left(\simeq 1.10\right)$	$5.6562\times {10}^{-4}$	1.15 $\left(\simeq 1.10\right)$
$1/112$	$2.5186\times {10}^{-4}$	1.10 $\left(\simeq 1.10\right)$	$2.4818\times {10}^{-4}$	1.10 $\left(\simeq 1.10\right)$	$2.6363\times {10}^{-4}$	1.10 $\left(\simeq 1.10\right)$
$1/224$	$1.1655\times {10}^{-4}$	1.11 $\left(\simeq 1.10\right)$	$1.1463\times {10}^{-4}$	1.11 $\left(\simeq 1.10\right)$	$1.2194\times {10}^{-4}$	1.11 $\left(\simeq 1.10\right)$
$1/448$	$5.4151\times {10}^{-5}$	1.11 $\left(\simeq 1.10\right)$	$5.3182\times {10}^{-5}$	1.11 $\left(\simeq 1.10\right)$	$5.6646\times {10}^{-5}$	1.11 $\left(\simeq 1.10\right)$

and

Table 2. $H^1$-errors, $E_{N_i}^2$, $i\ge 1$ calculated by (55) and convergence orders, $R_i^2$, $i\ge 2$, calculated by (56) for different values of ${\tau }_i\approx h_i$, $N_i=1/{\tau }_i$. The numbers in parentheses indicate the theoretical convergence order.

$\mathit{\tau }\approx \mathit{h}$	$\overline{\mathit{\alpha }}=\overline{\mathit{\beta }}=1.1$		$\overline{\mathit{\alpha }}=\overline{\mathit{\beta }}=1.5$		$\overline{\mathit{\alpha }}=\overline{\mathit{\beta }}=1.9$
	${\mathit{H}}^1$-Error	Order	${\mathit{H}}^1$-Error	Order	${\mathit{H}}^1$-Error	Order
$\overline{\gamma }=1.1$
$1/7$	$9.2985\times {10}^{-2}$	-	$9.3026\times {10}^{-2}$	-	$9.3255\times {10}^{-2}$	-
$1/14$	$4.1026\times {10}^{-2}$	1.18 $\left(\simeq 1.00\right)$	$4.1034\times {10}^{-2}$	1.18 $\left(\simeq 1.00\right)$	$4.1122\times {10}^{-2}$	1.18 $\left(\simeq 1.00\right)$
$1/28$	$1.9187\times {10}^{-2}$	1.10 $\left(\simeq 1.00\right)$	$1.9189\times {10}^{-2}$	1.10 $\left(\simeq 1.00\right)$	$1.9223\times {10}^{-2}$	1.10 $\left(\simeq 1.00\right)$
$1/56$	$9.2678\times {10}^{-3}$	1.05 $\left(\simeq 1.00\right)$	$9.2681\times {10}^{-3}$	1.05 $\left(\simeq 1.00\right)$	$9.2818\times {10}^{-3}$	1.05 $\left(\simeq 1.00\right)$
$1/112$	$4.6342\times {10}^{-3}$	1.00 $\left(\simeq 1.00\right)$	$4.6342\times {10}^{-3}$	1.00 $\left(\simeq 1.00\right)$	$4.6398\times {10}^{-3}$	1.00 $\left(\simeq 1.00\right)$
$1/224$	$2.2968\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$	$2.2968\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$	$2.2992\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$
$1/448$	$1.1434\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$	$1.1434\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$	$1.1444\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$
$\overline{\gamma }=1.5$
$1/7$	$9.3673\times {10}^{-2}$	-	$9.3833\times {10}^{-2}$	-	$9.4468\times {10}^{-2}$	-
$1/14$	$4.1217\times {10}^{-2}$	1.18 $\left(\simeq 1.00\right)$	$4.1258\times {10}^{-2}$	1.19 $\left(\simeq 1.00\right)$	$4.1506\times {10}^{-2}$	1.19 $\left(\simeq 1.00\right)$
$1/28$	$1.9238\times {10}^{-2}$	1.10 $\left(\simeq 1.00\right)$	$1.9248\times {10}^{-2}$	1.10 $\left(\simeq 1.00\right)$	$1.9342\times {10}^{-2}$	1.10 $\left(\simeq 1.00\right)$
$1/56$	$9.2809\times {10}^{-3}$	1.05 $\left(\simeq 1.00\right)$	$9.2836\times {10}^{-3}$	1.05 $\left(\simeq 1.00\right)$	$9.3187\times {10}^{-3}$	1.05 $\left(\simeq 1.00\right)$
$1/112$	$4.6375\times {10}^{-3}$	1.00 $\left(\simeq 1.00\right)$	$4.6382\times {10}^{-3}$	1.00 $\left(\simeq 1.00\right)$	$4.6512\times {10}^{-3}$	1.00 $\left(\simeq 1.00\right)$
$1/224$	$2.2977\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$	$2.2978\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$	$2.3028\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$
$1/448$	$1.1436\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$	$1.1436\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$	$1.1455\times {10}^{-3}$	1.01 $\left(\simeq 1.00\right)$
$\overline{\gamma }=1.9$
$1/7$	$1.0462\times {10}^{-1}$	-	$1.0487\times {10}^{-1}$	-	$1.0654\times {10}^{-1}$	-
$1/14$	$4.6312\times {10}^{-2}$	1.18 $\left(\simeq 1.00\right)$	$4.6339\times {10}^{-2}$	1.18 $\left(\simeq 1.00\right)$	$4.7081\times {10}^{-2}$	1.18 $\left(\simeq 1.00\right)$
$1/28$	$2.1550\times {10}^{-2}$	1.10 $\left(\simeq 1.00\right)$	$2.1536\times {10}^{-2}$	1.11 $\left(\simeq 1.00\right)$	$2.1869\times {10}^{-2}$	1.11 $\left(\simeq 1.00\right)$
$1/56$	$1.0314\times {10}^{-2}$	1.06 $\left(\simeq 1.00\right)$	$1.0301\times {10}^{-2}$	1.06 $\left(\simeq 1.00\right)$	$1.0451\times {10}^{-2}$	1.07 $\left(\simeq 1.00\right)$
$1/112$	$5.0944\times {10}^{-3}$	1.02 $\left(\simeq 1.00\right)$	$5.0861\times {10}^{-3}$	1.02 $\left(\simeq 1.00\right)$	$5.1542\times {10}^{-3}$	1.02 $\left(\simeq 1.00\right)$
$1/224$	$2.4989\times {10}^{-3}$	1.03 $\left(\simeq 1.00\right)$	$2.4945\times {10}^{-3}$	1.03 $\left(\simeq 1.00\right)$	$2.5251\times {10}^{-3}$	1.03 $\left(\simeq 1.00\right)$
$1/448$	$1.2319\times {10}^{-3}$	1.02 $\left(\simeq 1.00\right)$	$1.2297\times {10}^{-3}$	1.02 $\left(\simeq 1.00\right)$	$1.2434\times {10}^{-3}$	1.02 $\left(\simeq 1.00\right)$

. Then, for every ${\tau }_i$, the triangulation of $\overline{\mathrm{\Omega }}$ with the diameter $h_i$ such that $\left|h_i-{\tau }_i\right|<5\times {10}^{-4}$ was generated. We calculated the corresponding errors $E_{N_i}^1$ and $E_{N_i}^2$ for each ${\tau }_i$, $N_i=1/{\tau }_i$, $i\ge 1$ and the convergence orders by the formula

$$R_i^m=\frac{log{E}_{N_i}^m-log{E}_{N_{i-1}}^m}{log{\tau }_i-log{\tau }_{i-1}},m=1,2,i\ge 2.$$

(56)

In the experiments, several combinations of the fractional derivative orders, $\overline{\alpha }$, $\overline{\beta }$ and $\overline{\gamma }$, from the set $\left\{1.1,1.5,1.9\right\}$ were considered, and the values of $\overline{\alpha }$ and $\overline{\beta }$ were assumed to be equal for simplicity. The calculation results are shown in Table 1 and Table 2, respectively. In the ”Order“ column, the theoretical convergence orders are indicated in parentheses. Corresponding errors plot are depicted in Figure 1 and Figure 2.

It follows from Table 1 that the convergence order significantly depended on the orders of fractional derivatives. To be more precise, it is clearly seen that the empirical convergence order, $E_N^1$, was around 1.90 in case of $\overline{\alpha }=\overline{\beta }=\overline{\gamma }=1.1$, and it significantly reduced with increasing order of the fractional derivatives. The convergence order was close to 1.10 when any of $\overline{\alpha },\overline{\beta }$ or $\overline{\gamma }$ was equal to 1.90. This behavior confirms the theoretical order $O\left({\tau }^{3-\nu }\right)$, $\nu =max\left\{\overline{\alpha },\overline{\beta },\overline{\gamma }\right\}$ predicted in Theorem 4.

Similarly, the results presented in Table 2 demonstrate that empirical convergence order for $E_N^2$ approached 1.00 as the discretization parameter decreased. Thus, it can be concluded that numerical results conducted for Example 1 fully confirmed the theoretically predicted convergence order obtained in Theorem 4.

4. Discussion

Having carried out a theoretical analysis of the method for solving the problem under consideration as well as analyzing the results of computational experiments, the following conclusions can be drawn:

(1) Based on Theorem 1, it is concluded that the solution of the problem under consideration is unique and continuously depends on the input data.

(2) A computational method constructed to solve the problem under consideration, in the case of $\overline{\alpha },\overline{\beta },\overline{\gamma }\in \left(1,2\right)$, allows one to obtain an approximate solution with an order of accuracy $O\left(h^2+h^2{\tau }^{1/2}+{\tau }^{3-\overline{\alpha }}+{\tau }^{3-\overline{\beta }}+{\tau }^{3-\overline{\gamma }}\right)$. The results of computational experiments carried out for different orders of fractional derivatives and mesh configurations fully confirm the results of theoretical analysis. Since the proposed method makes it possible to obtain an approximate solution for an arbitrary choice of orders of fractional derivatives, $\overline{\alpha },\overline{\beta },\overline{\gamma }\in \left(1,2\right)$, this significantly improves the result obtained in [29].

(3) In subsequent studies, the authors intend to consider more realistic examples of the problems of fluid flow in a fractured porous medium on the base of the constructed method. Specifically, a more extensive comparison of the simulation results with the results of similar works will be carried out.

(4) The outlines obtained in this paper can be used when considering other classes of problems for fractional differential equations of the order belonging to the interval $\left(1,2\right)$. These studies include, for example, multidimensional models of fluid flow in fractured porous media with a fractal fracture geometry and models of single-phase fluid flow with spatial fractional derivatives.