Numerical Method for Fractional-Order Generalization of the Stochastic Stokes–Darcy Model

Abstract

This paper is aimed at efficient numerical implementation of the fractional-order generalization of the stochastic Stokes–Darcy model, which has important scientific, applied, and economic significance in hydrology, the oil industry, and biomedicine. The essence of this generalization of the stochastic model is the introduction of fractional time derivatives in the sense of Caputo’s definition to take into account long-term changes in the properties of media. An efficient numerical method for the implementation of the fractional-order Stokes–Darcy model is proposed, which is based on the use of a higher-order approximation formula for the fractional derivative, higher-order finite difference relations, and a finite element approximation of the problem in the spatial direction. In the paper, a rigorous theoretical analysis of the stability and convergence of the proposed numerical method is carried out, which is confirmed by numerous computational experiments. Further, the proposed method is applied to the implementation of the fractional-order stochastic Stokes–Darcy model using an ensemble technique, in which the approximation is carried out in such a way that the resulting systems of linear equations have the same coefficient matrix for all realizations. Furthermore, evaluation of the discrete fractional derivatives is carried out with the use of parallel threads. The efficiency of applying both approaches has been demonstrated in numerical tests.

1. Introduction

The model described by the joint Stokes–Darcy equations occupies an important place in computational fluid dynamics due to its wide application in problems of interaction between surface and underground flows, in predicting processes occurring in oil reservoirs with a cavernous-porous structure, in underground systems of karst aquifers, and others. In addition, this model is used in forecasting and assessing the risk of flooding of territories as a result of groundwater and surface water. In this model, the free fluid flow is described using the Stokes equation, and the flow in a porous medium is described by the Darcy equation. It is believed [1] that the joint Stokes–Darcy model is more adequate than each of the models considered separately.

This model is most fully studied under the assumption that the input parameters of the problem, such as physical quantities describing the model, media properties, exposure conditions, domain geometry, boundary conditions, and initial conditions, are known. In reality, there are many uncertainties that arise from the difficulty in determining real data due to measurement noise, internal variability in physical quantities, and the adoption of simplifying assumptions. As a result, a plausible description of the fluid flow process using classical deterministic partial differential equations seems impossible. This often gives rise to the risks of making inadequate technological and managerial decisions. In practice, knowledge about the qualitative and quantitative properties of a porous medium is often fragmentary, since the input data are formed on the basis of laboratory analysis of samples taken at several rarefied points of the territory with their subsequent interpolation.

The study of the stochastic generalization of the deterministic Stokes–Darcy model is a fairly new and relevant area in modern computational fluid dynamics. The earliest work in this direction, in our opinion, is the paper of Kumar et al. [2] published in 2018. Interest in this model has grown rapidly over the course of 3–4 years, as evidenced by numerous works published to this day [1,3,4,5,6].

A significant drawback of the stochastic Stokes–Darcy model is the neglect of an extremely important property of a porous medium—memory, which largely affects the nature of the fluid flow. The need to take into account this property lies in the fact that the porosity and permeability of the medium change in the process of fluid flow, and more noticeably lead to a deceleration in the flow rate over time. Therefore, the nature of the fluid flow is determined not only by the current state of the porous medium, but also by all its previous states.

Taking into account the importance of this direction, this paper studies a further generalization of the stochastic Stokes–Darcy model, which takes into account long-term changes in the properties of media. Accounting for memory in our work is achieved by using fractional differential calculus; namely, by replacing integer time derivatives with fractional-order derivatives in the sense of Caputo’s definition. This generalization is a natural generalization of the works of scientists in the field of fractional differential fluid flow theory [7,8,9] and is based on the hypotheses put forward in [7].

The main problem associated with the considered model is the complexity of its numerical implementation. For example, the application of the stochastic Galerkin method to this model, which is based on the expansion in generalized polynomial chaos, is difficult due to the complication of using existing solvers without their serious modification [2]. Stochastic collocation methods [10,11] do not have this disadvantage; however, the number of nodes increases with an increase in stochastic dimensions, which leads to a significant increase in computational complexity. The very popular Monte Carlo method, being easy to implement and allowing the use of well-known deterministic algorithms and having a high potential for parallelization, has a rather low convergence rate [3]. Much effort has been put into the development of fast convergent methods, such as quasi-Monte Carlo methods [12], multilevel Monte Carlo methods [2,13,14], centroidal Voronoi tilings [15], and Latin hypercube selection [16]. Despite the fact that these methods make it possible to use existing deterministic algorithms without changes, the complex form of the equations in the proposed fractional-order stochastic model leads to high computational costs. In addition, accounting for uncertainty usually leads to the need to run the solvers of these methods multiple times in order to obtain statistical moments. This drawback is partially overcome in the so-called ensemble algorithms [17], in which the problem is solved simultaneously for all realizations in the ensemble instead of solving the problem individually for each realization. In the latter work, this was achieved by approximating the nonlinear term $\mathbf{u}·\nabla \mathbf{u}$ as follows:

$${\mathbf{u}}^n·\nabla {\mathbf{u}}^{n+1}\approx {\overline{\mathbf{u}}}^n·\nabla {\mathbf{u}}^{n+1}+\left({\mathbf{u}}^n-{\overline{\mathbf{u}}}^n\right)·\nabla {\mathbf{u}}^n,$$

where ${\overline{\mathbf{u}}}^n$ is the average velocity for all realizations in the ensemble. As a result of this approximation, the sought coefficients of ${\mathbf{u}}^{n+1}$ do not depend on the current realization in the ensemble. Therefore, the systems of linear equations for all realizations have the same coefficient matrix, but different right-hand sides. This makes it possible to use efficient factorization methods for systems of linear equations, which significantly speeds up the solution of the problem.

Since the permeability fields do not have obvious separation of scales, a multiscale method [18,19] is an effective method for solving fluid flow problems because of its usage of limited global information. In particular, in [18] the problem of fluid motion in a porous medium is solved using a stochastic mixed multiscale finite element method (FEM) on a coarse grid. This approach does not require interpolation in the stochastic space and can be combined with approaches based on interpolation on a coarse grid. In addition to this, in [19], a new mixed multiscale FEM is proposed within a heterogeneous multiscale method that provides both local and global mass conservation. In [20], a stochastic FEM solution based on the perturbation method is developed. The paper [5] proposes three algorithms for quantifying uncertainty. Due to the complexity of the media under consideration, the initial data can be characterized by a weak singularity. Recent studies address this class of problems and propose Sinc-collocation methods [21] and methods based on orthogonal spline collocation [22,23] for their solution. In [3], an efficient approach for implementing the stochastic Stokes–Darcy model with a random hydraulic conductivity tensor is proposed, which is based on separation of the Stokes–Darcy system into two subproblems, as well as on assembling a common matrix of coefficients for all realizations in the ensemble. This approach significantly speeds up the solution of the resulting systems of linear equations without affecting the accuracy of calculations. In [24], a stabilized mixed method for solving the problem is proposed that does not use Lagrange multipliers. In [1,25,26], a stochastic collocation method on sparse grids is developed for the Stokes–Darcy model with random hydraulic conductivity. In [27], an ensemble domain decomposition method is proposed, the peculiarity of which is that the solution of the problem is reduced to solving a system of linear algebraic equations with a common matrix of coefficients for each deterministic numerical model.

Recent studies have focused on the construction of higher-order numerical schemes for the Stokes–Darcy equations with integer time derivatives. For example, in [6], a second-order ensemble numerical method is proposed for stochastic Stokes–Darcy equations based on the results of [3,17]. Li et al. [28] constructed a second-order fractional time-stepping method. Qin et al. [29] and Li et al. [30] proposed numerical schemes that combine a time filter and a Backward Euler scheme to increase the convergence order without increasing the amount of computation. Furthermore, [29] obtained a third-order scheme by applying a time filter to the BDF2 scheme. Chen et al. [31] proposed a third-order in time Adam–Moulton–Bashforth method. Further, recent papers have studied the nonlinear generalizations of Stokes–Darcy equations [32,33,34] by applying domain decomposition techniques, optimization methods, and mortar finite element methods. Other effective techniques should be noted, such as the optimal homotopy asymptotic method [35,36], which was successfully applied to studying nonlinear behaviors described by partial differential equations containing fractional-order time derivatives [37,38].

The present study is aimed at continuing the development of ensemble methods for solving stochastic Stokes–Darcy equations for cases when the equations contain fractional time derivatives. This work has several key differences from previous studies. First, an ensemble method for implementing a fractional-order stochastic generalization of the Stokes–Darcy model is constructed in the work. In our opinion, this is the first study of the specified generalization of the model. Secondly, in contrast to the original works [3,6] employing the ensemble technique, the constructed computational scheme has an order of convergence greater than 2 in the time variable. This is achieved by using an approximation formula of the order $O\left({\tau }^{3-\nu }\right)$ for the fractional derivative in the sense of Caputo, as well as by a special approximation of some terms. Thirdly, the paper proposes a parallel algorithm to accelerate the computation process and analyzes the acceleration rate on a test problem.

The main contribution of this study is a higher-order numerical scheme allowing for parallel execution, which makes it possible to efficiently solve the fractional-order stochastic Stokes–Darcy equations.

The structure of the article is as follows. In the next section, the formulation of a new initial boundary-value problem for the fractional-order stochastic generalization of the Stokes–Darcy equation is given. Section 2.2 proposes a method for solving this problem. Section 2.3 is devoted to the theoretical study of the solution method. Section 3 contains the results of computational experiments to confirm the results of the theoretical analysis. Finally, Section 4 discusses the results obtained.

2. Materials and Methods

2.1. Formulation of the Problem

Consider the fluid flow process in two non-intersecting bounded domains ${\mathrm{\Omega }}_f\subset {\mathbb{R}}^d$ and ${\mathrm{\Omega }}_p\subset {\mathbb{R}}^d$ $\left(d=2,3\right)$ separated by an interface I, i.e., ${\mathrm{\Omega }}_f\cap {\mathrm{\Omega }}_p=\varnothing$, ${\overline{\mathrm{\Omega }}}_f\cap {\overline{\mathrm{\Omega }}}_p=I$ and ${\mathrm{\Omega }}_f\cup {\mathrm{\Omega }}_p=\mathrm{\Omega }$ (Figure 1). It is assumed that a free surface flow occurs in ${\mathrm{\Omega }}_f$, then the flow moves into the domain ${\mathrm{\Omega }}_p$ with a porous medium. Let ${\mathbf{n}}_f$ and ${\mathbf{n}}_p$ denote outer unit normal vectors to the boundaries $\partial {\mathrm{\Omega }}_f$ and $\partial {\mathrm{\Omega }}_p$, and ${\left({\mathbf{m}}_i\right)}_{i=1,2,\dots ,d-1}$ denote an orthonormal system of tangential vectors to the interface I. Further, the fluid flow process is studied on the time interval $\overline{J}$, where $J=\left(0,T\right]$ for $T>0$.

Let $\mathbf{u}$ and p denote the fluid velocity and pressure in ${\mathrm{\Omega }}_f$, respectively; let $\varphi$ be the piezometric head in ${\mathrm{\Omega }}_p$; let $\mu$ be the fluid viscosity and $\mathbb{K}$ be the hydraulic conductivity tensor of the porous medium. Then, the fractional-order generalization of the linearized Stokes–Darcy model reads as follows:

Problem 1.

Find $\mathbf{u},p,\varphi$ in $\mathrm{\Omega }\times J$ satisfying the following conditions:

$$\begin{array}{cccc}& {\partial }_{0,t}^{\alpha }\mathbf{u}-\mu {\nabla }^2\mathbf{u}+\nabla p={\mathbf{f}}_f,\hfill & \hfill & \left(x,t\right)\in {\mathrm{\Omega }}_f\times J,\alpha \in \left(0,1\right),\hfill \end{array}$$

(1)

$$\begin{array}{cccc}& \nabla ·\mathbf{u}=0,\hfill & & \left(x,t\right)\in {\mathrm{\Omega }}_f\times J,\hfill \end{array}$$

(2)

$$\begin{array}{cccc}& \mathbf{u}=\mathbf{0},\hfill & & \left(x,t\right)\in \left(\partial {\mathrm{\Omega }}_f\setminus I\right)\times J,\hfill \end{array}$$

(3)

$$\begin{array}{cccc}& S{\partial }_{0,t}^{\beta }\varphi -\nabla ·\left(\mathbb{K}\nabla \varphi \right)=f_p,\hfill & & \left(x,t\right)\in {\mathrm{\Omega }}_p\times J,\beta \in \left(0,1\right),\hfill \end{array}$$

(4)

$$\begin{array}{cccc}& \mathbb{K}\nabla \varphi ·{\mathbf{n}}_p=0,\hfill & & \left(x,t\right)\in \left(\partial {\mathrm{\Omega }}_p\setminus I\right)\times J,\hfill \end{array}$$

(5)

$$\begin{array}{cccc}& \mathbf{u}·{\mathbf{n}}_f-\mathbb{K}\nabla \varphi ·{\mathbf{n}}_p=0,\hfill & & \left(x,t\right)\in I\times J,\hfill \end{array}$$

(6)

$$\begin{array}{cccc}& p-\mu {\mathbf{n}}_f\nabla \mathbf{u}·{\mathbf{n}}_f=g\varphi ,\hfill & & \left(x,t\right)\in I\times J,\hfill \end{array}$$

(7)

$$\begin{array}{cccc}& -{\mathbf{n}}_f·\nabla \mathbf{u}·{\mathbf{m}}_i=\frac{{\alpha }_{\mathrm{BJS}}}{\sqrt{{\mathbf{m}}_i·\mathbb{K}{\mathbf{m}}_i}}\mathbf{u}·{\mathbf{m}}_i,\hfill & & \left(x,t\right)\in I\times J,1\le i\le d-1,\hfill \end{array}$$

(8)

$$\begin{array}{cccc}& \mathbf{u}={\mathbf{u}}_0,\hfill & & \left(x,t\right)\in {\overline{\mathrm{\Omega }}}_f\times \left\{0\right\},\hfill \end{array}$$

(9)

$$\begin{array}{cccc}& \varphi ={\varphi }_0,\hfill & & \left(x,t\right)\in {\overline{\mathrm{\Omega }}}_p\times \left\{0\right\},\hfill \end{array}$$

(10)

where the fractional derivative in the sense of Caputo is defined as

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

and S is the reservoir capacity coefficient, g is the gravitational acceleration constant, and ${\alpha }_{\mathrm{BJS}}$ is a dimensionless parameter that depends on the structure of the porous medium.

A feature of this problem is that the hydraulic conductivity $\mathbb{K}$, the initial data ${\mathbf{u}}_0$, ${\varphi }_0$, and the right-hand sides ${\mathbf{f}}_f$, $f_p$ are not fully known in practice. Most often, the mean and variance of these values are available as a result of geological study. Therefore, there is a high degree of uncertainty that must be taken into account when modeling a fluid flow process in such complex media.

To account for the uncertainty effects, we utilize a commonly acceptable approach, in which a set of input data ${\left\{\left({\mathbb{K}}_{\omega },{\mathbf{u}}_{0,\omega },{\varphi }_{0,\omega },{\mathbf{f}}_{f,\omega },f_{p,\omega }\right)\right\}}_{\omega =1}^{N_{\omega }}$ is randomly generated, then Problem 1 is solved $N_{\omega }$ times with specified data. As a result, one obtains a series of problems instead of Problem 1. Thus, a set of solutions ${\left\{\left({\mathbf{u}}_{\omega },p_{\omega },{\varphi }_{\omega }\right)\right\}}_{\omega =1}^{N_{\omega }}$ is found, with the help of which the statistical moments are subsequently determined.

Therefore, the rest of this section, as well as Section 2.2, Section 2.3 and Section 2.4, are aimed at developing a numerical method for solving the problem for one of the $N_{\omega }$ realizations. Further, in Section 3, the proposed method is applied to the solution of a stochastic problem that covers $N_{\omega }$ samples of input parameters.

We introduce the following spaces:

$$\begin{array}{cc}& X_f=\left\{\mathbf{v}\in {\left(H^1\left({\mathrm{\Omega }}_f\right)\right)}^d|\mathbf{v}=\mathbf{0}\mathrm{on}\partial {\mathrm{\Omega }}_f\setminus I\right\},\hfill \end{array}$$

(11)

$$\begin{array}{cc}& Q_f=\left\{q\in L^2\left({\mathrm{\Omega }}_f\right)|{\int }_{{\mathrm{\Omega }}_f}qdx=0\right\},\hfill \end{array}$$

(12)

$$\begin{array}{cc}& Q_p=\left\{\psi \in H^1\left({\mathrm{\Omega }}_p\right)|\psi =0\mathrm{on}\partial {\mathrm{\Omega }}_p\setminus I\right\},\hfill \\ & \mathbf{X}=X_f\times Q_p,\hfill \end{array}$$

(13)

where the standard notation for Sobolev spaces $H^k\left(\mathrm{\Omega }\right)$ and Lebesgue spaces $L^p\left(\mathrm{\Omega }\right)$ is utilized. We use the notation ${\left(·,·\right)}_{\mathrm{\Omega }}$ to denote the dot product in $L^2\left(\mathrm{\Omega }\right)$, the notations ${∥\mathbf{v}∥}_{{\mathrm{\Omega }}_f}$ and ${∥\psi ∥}_{{\mathrm{\Omega }}_p}$ stand for ${∥\mathbf{v}∥}_{{\mathbf{L}}^2\left({\mathrm{\Omega }}_f\right)}$ and ${∥\psi ∥}_{L^2\left({\mathrm{\Omega }}_p\right)}$ for brevity, and the norm in $\mathbf{X}$ is defined as

$$∥\widehat{\mathbf{v}}∥={\left({∥\mathbf{v}∥}_{{\mathrm{\Omega }}_f}^2+{∥\psi ∥}_{{\mathrm{\Omega }}_p}^2\right)}^{1/2},\widehat{\mathbf{v}}\in \mathbf{X}.$$

In order to construct a weak formulation of Problem 1 for the $\omega$-th implementation, we multiply Equations (1), (2) and (4) by arbitrary elements $\mathbf{v}\in X_f$, $q\in Q_f$, and $g\psi \in Q_p$, respectively, then integrate the resulting equations over the corresponding domains ${\mathrm{\Omega }}_f$ and ${\mathrm{\Omega }}_p$, taking into account the boundary conditions (3) and (5)–(8), and use the expansion of the test velocity vector $\mathbf{v}$ in terms of its normal and tangential components:

$$\mathbf{v}=\left(\mathbf{v}·{\mathbf{n}}_f\right){\mathbf{n}}_f+\sum _{i=1}^{d-1}\left(\mathbf{v}·{\mathbf{m}}_i\right){\mathbf{m}}_i.$$

(14)

Moreover, we introduce the notation $\widehat{\mathbf{u}}={\left[\mathbf{u},\varphi \right]}^{\top }$, $\widehat{\mathbf{v}}={\left[\mathbf{v},\psi \right]}^{\top }$ for brevity.

Problem 2

(Weak formulation). Find $\left({\widehat{\mathbf{u}}}_{\omega },p_{\omega }\right):J\mapsto \mathbf{X}\times Q_f$ satisfying the identity

$$\begin{array}{cc}& {\mathcal{D}}^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)+\mathcal{A}\left({\mathbb{K}}_{\omega };{\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)-\mathcal{B}\left(\widehat{\mathbf{v}},p_{\omega }\right)+\mathcal{B}\left({\widehat{\mathbf{u}}}_{\omega },q\right)\hfill \\ & +{\mathcal{A}}_I\left({\gamma }_{\omega };{\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)-\mathcal{I}\left({\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)=\left({\widehat{\mathbf{f}}}_{\omega },\widehat{\mathbf{v}}\right),\hfill \end{array}$$

(15)

for all $\left(\widehat{\mathbf{v}},q\right)\in \mathbf{X}\times Q_f$, where

$$\begin{array}{cc}& {\mathcal{D}}^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)={\left({\partial }_{0,t}^{\alpha }{\mathbf{u}}_{\omega },\mathbf{v}\right)}_{{\mathrm{\Omega }}_f}+{\left(gS{\partial }_{0,t}^{\beta }{\varphi }_{\omega },\psi \right)}_{{\mathrm{\Omega }}_p},\hfill \\ & \mathcal{A}\left({\mathbb{K}}_{\omega };{\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)={\mathcal{A}}_f\left({\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)+{\mathcal{A}}_p\left({\mathbb{K}}_{\omega };{\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right),\hfill \\ & {\mathcal{A}}_f\left({\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)={\left(\mu \nabla {\mathbf{u}}_{\omega },\nabla \mathbf{v}\right)}_{{\mathrm{\Omega }}_f},\hfill \\ & {\mathcal{A}}_p\left({\mathbb{K}}_{\omega };{\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)=g{\left({\mathbb{K}}_{\omega }\nabla {\varphi }_{\omega },\nabla \psi \right)}_{{\mathrm{\Omega }}_p},\hfill \\ & {\mathcal{A}}_I\left({\gamma }_{\omega };{\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)=\mu {\int }_I{\gamma }_{\omega }\left({\mathbf{u}}_{\omega }·{\mathbf{m}}_i\right)\left(\mathbf{v}·{\mathbf{m}}_i\right)ds,\hfill \\ & \mathcal{I}\left({\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)={\mathcal{H}}_I\left(\mathbf{v},{\varphi }_{\omega }\right)-{\mathcal{H}}_I\left({\mathbf{u}}_{\omega },\psi \right),\hfill \\ & \mathcal{B}\left(\widehat{\mathbf{v}},q\right)={\left(q,\nabla ·\mathbf{v}\right)}_{{\mathrm{\Omega }}_f},\hfill \\ & \left({\widehat{\mathbf{f}}}_{\omega },\widehat{\mathbf{v}}\right)={\left({\mathbf{f}}_{f,\omega },\mathbf{v}\right)}_{{\mathrm{\Omega }}_f}+{\left(g{f}_{p,\omega },\psi \right)}_{{\mathrm{\Omega }}_p},\hfill \end{array}$$

(16)

and the following notation is introduced:

$$\begin{array}{cc}& {\mathcal{H}}_I\left(\mathbf{v},{\varphi }_{\omega }\right)=g{\int }_I{\varphi }_{\omega }\left(\mathbf{v}·{\mathbf{n}}_f\right)ds,\hfill \\ & {\displaystyle {\gamma }_{i,\omega }=\frac{{\alpha }_{\mathrm{BJS}}}{\sqrt{{\mathbf{m}}_i·{\mathbb{K}}_{\omega }{\mathbf{m}}_i}}},{\gamma }_{\omega }=\sum _{i=1}^{d-1}{\gamma }_{i,\omega }.\hfill \end{array}$$

(17)

2.2. Numerical Method

We construct a numerical method for solving Problem 1 in two stages. First, we discretize the problem with respect to the time variable [39]. To do this, we introduce a uniform partition of the time interval ${\overline{J}}_{\tau }=\left\{t_n=n\tau ,n=0,1,\dots ,N,N\tau =T\right\}$ for the time discretization parameter $\tau >0$, and the timestamp $t=t_n$ will be referred to as the nth time layer. We also introduce the notation $f^n=f\left(·,t_n\right)$ for convenience.

In general, the construction of our computational method is based on the use of different approaches to discretize the problem on different time layers to increase the accuracy of the method. The first approach is applied to the first two time layers, $t_1$ and $t_2$, and is based on the use of a Crank–Nicolson-type approximation [40]. The second approach covers the rest of the time layers, $n\ge 3$, and is based on the use of higher-order approximation formulas. In particular, we employ the following result from [41] in order to approximate fractional derivatives starting from the third time layer.

Lemma 1.

The discrete analogue ${\Delta }_{\tau }^{\nu }{f}^n$ of the fractional derivative in the sense of Caputo ${\partial }_{0,t}^{\nu }f\left(t_n\right)$ of order $\nu \in \left(0,1\right)$ can be represented as [41]

$${\Delta }_{\tau }^{\nu }{f}^n=\frac{{\lambda }_0^{\left(\nu \right)}{\tau }^{-\nu }}{{\varsigma }_0^{\left(\nu \right)}}\left(f^n-\sum _{s=1}^n{\delta }_{n,n-s}^{\left(\nu \right)}{f}^{n-s}\right),n\ge 3,$$

(18)

where

$$\begin{array}{cc}& {\varsigma }_0^{\left(\nu \right)}=\mathrm{\Gamma }\left(3-\nu \right),{\lambda }_0^{\left(\nu \right)}=2^{1-\nu }\left(1+\frac{\nu }{2}\right),\hfill \\ & {\delta }_{3,2}^{\left(\nu \right)}=-\frac{b_1^{\left(\nu \right)}+c_2^{\left(\nu \right)}-2}{{\lambda }_0^{\left(\nu \right)}},{\delta }_{3,1}^{\left(\nu \right)}=-\frac{a_1^{\left(\nu \right)}+b_2^{\left(\nu \right)}+\frac{\nu }{2}}{{\lambda }_0^{\left(\nu \right)}},{\delta }_{3,0}^{\left(\nu \right)}=-\frac{a_2^{\left(\nu \right)}}{{\lambda }_0^{\left(\nu \right)}},\hfill \end{array}$$

and in the case of $n\ge 4$:

$$\begin{array}{cc}& {\delta }_{n,n-1}^{\left(\nu \right)}=-\frac{b_1^{\left(\nu \right)}+c_2^{\left(\nu \right)}-2}{{\lambda }_0^{\left(\nu \right)}},{\delta }_{n,n-2}^{\left(\nu \right)}=-\frac{a_1^{\left(\nu \right)}+b_2^{\left(\nu \right)}+c_3^{\left(\nu \right)}+\frac{\nu }{2}}{{\lambda }_0^{\left(\nu \right)}},\hfill \\ & {\delta }_{n,n-s}^{\left(\nu \right)}=-\frac{a_{s-1}^{\left(\nu \right)}+b_s^{\left(\nu \right)}+c_{s+1}^{\left(\nu \right)}}{{\lambda }_0^{\left(\nu \right)}},s=3,4,\dots ,n-2,\hfill \\ & {\delta }_{n,1}^{\left(\nu \right)}=-\frac{a_{n-2}^{\left(\nu \right)}+b_{n-1}^{\left(\nu \right)}}{{\lambda }_0^{\left(\nu \right)}},{\delta }_{n,0}^{\left(\nu \right)}=-\frac{a_{n-1}^{\left(\nu \right)}}{{\lambda }_0^{\left(\nu \right)}},\hfill \end{array}$$

and

$$\begin{array}{cc}& a_s^{\left(\nu \right)}=-\frac{3}{2}\left(2-\nu \right){\left(s+1\right)}^{1-\nu }+\frac{1}{2}\left(2-\nu \right)s^{1-\nu }+{\left(s+1\right)}^{2-\nu }-s^{2-\nu },\hfill \\ & b_s^{\left(\nu \right)}=2\left(2-\nu \right){\left(s+1\right)}^{1-\nu }-2{\left(s+1\right)}^{2-\nu }+2s^{2-\nu },\hfill \\ & c_s^{\left(\nu \right)}=-\frac{1}{2}\left(2-\nu \right)\left({\left(s+1\right)}^{1-\nu }+s^{1-\nu }\right)+{\left(s+1\right)}^{2-\nu }-s^{2-\nu },\hfill \end{array}$$

where the quantity $R_n^{\left(\nu \right)}={\partial }_{0,t_n}^{\nu }f-{\Delta }_{\tau }^{\nu }{f}^n$, $n\ge 2$ satisfies the estimate

$$\begin{array}{cc}\hfill & \left|R_n^{\left(\nu \right)}\right|\le \left(\frac{\sqrt{3}}{27\mathrm{\Gamma }\left(1-\nu \right)}+\frac{1}{6\mathrm{\Gamma }\left(2-\nu \right)}\right)\underset{t_0\le t\le t_n}{max}\left|{\partial }_t^{3}f\left(t\right)\right|{\tau }^{3-\nu }.\hfill \end{array}$$

The properties of the coefficients ${\delta }_{n,s}^{\left(\nu \right)}$ have been rigorously studied in [41]. The following lemma summarizes the properties that will be used while proving the stability and convergence of the proposed numerical scheme.

Lemma 2.

The coefficients ${\delta }_{n,s}^{\left(\nu \right)}$ have the following properties [41]:

(1) ${\delta }_{n,n-1}^{\left(\nu \right)}\in \left(0,\frac{4}{3}\right)$;

(2) ${\delta }_{n,n-2}^{\left(\nu \right)}\in \left(-\frac{1}{2},\frac{1}{3}\right)$;

(3) ${\delta }_{n,n-s}^{\left(\nu \right)}>0$, $s=3,4,\dots ,n$;

(4) ${\displaystyle \sum _{s=0}^{n-1}{\delta }_{n,s}^{\left(\nu \right)}=1}$.

Approximation of fractional derivatives on the first and second time layers is carried out according to another formula, which is formulated in the following lemma.

Lemma 3.

The discrete analogue ${\overline{\Delta }}_{\tau }^{\nu }{f}^n$ of the fractional derivative in the sense of Caputo ${\partial }_{0,t}^{\nu }f\left(t_n\right)$ of order $\nu \in \left(0,1\right)$ can be represented as [42]

$${\overline{\Delta }}_{\tau }^{\nu }{f}^n=\frac{{\tau }^{-\nu }}{{\overline{\varsigma }}_0^{\left(\nu \right)}}\sum _{s=1}^n{\overline{\delta }}_{n,s}^{\left(\nu \right)}\left(f^s-f^{s-1}\right),n=1,2,$$

(19)

where

$${\overline{\delta }}_{n,s}^{\left(\nu \right)}={\left(n-s+1\right)}^{1-\nu }-{\left(n-s\right)}^{1-\nu },{\overline{\varsigma }}_0^{\left(\nu \right)}=\mathrm{\Gamma }\left(2-\nu \right),$$

where the quantity ${\overline{R}}_n^{\left(\nu \right)}={\partial }_{0,t}^{\nu }f\left(t_n\right)-{\overline{\Delta }}_{\tau }^{\nu }{f}^n$ satisfies the estimate

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

By combining the results of Lemmas 1 and 3, we conclude that the following relation holds for the bilinear form ${\mathcal{D}}^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega },\widehat{\mathbf{v}}\right)$ on the nth time layer:

$${\mathcal{D}}^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)={\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)+{\left({\mathbf{r}}_{11}^n,\mathbf{v}\right)}_{{\mathrm{\Omega }}_f}+{\left(r_{12}^n,\psi \right)}_{{\mathrm{\Omega }}_p},n\ge 1,$$

(20)

where

$${\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)=\left\{\begin{array}{cc}{\left({\Delta }_{\tau }^{\alpha }{\mathbf{u}}_{\omega }^n,\mathbf{v}\right)}_{{\mathrm{\Omega }}_f}+{\left(gS{\Delta }_{\tau }^{\beta }{\varphi }_{\omega }^n,\psi \right)}_{{\mathrm{\Omega }}_p},\hfill & n\ge 3,\hfill \\ {\left({\overline{\Delta }}_{\tau }^{\alpha }{\mathbf{u}}_{\omega }^n,\mathbf{v}\right)}_{{\mathrm{\Omega }}_f}+{\left(gS{\overline{\Delta }}_{\tau }^{\beta }{\varphi }_{\omega }^n,\psi \right)}_{{\mathrm{\Omega }}_p},\hfill & n=1,2,\hfill \end{array}\right.$$

(21)

and ${\mathbf{r}}_{11}^n={\left[O\left({\tau }^{2-\alpha }\right),O\left({\tau }^{2-\alpha }\right)\right]}^{\top }$, $r_{12}^n=O\left({\tau }^{2-\beta }\right)$ when $n=1,2$; ${\mathbf{r}}_{11}^n={\left[O\left({\tau }^{3-\alpha }\right),O\left({\tau }^{3-\alpha }\right)\right]}^{\top }$, $r_{12}^n=O\left({\tau }^{3-\beta }\right)$ when $n\ge 3$.

Further, we transform the terms in (15) as shown below:

$${\mathcal{A}}_I\left({\gamma }_{\omega };{\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)={\mathcal{A}}_I\left(\sum _{i=1}^{d-1}\left({\overline{\gamma }}_i+\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right)\right);{\widehat{\mathbf{u}}}_{\omega }^n,\mathbf{v}\right)$$

$$={\mathcal{A}}_I\left(\sum _{i=1}^{d-1}{\overline{\gamma }}_i;{\widehat{\mathbf{u}}}_{\omega }^n,\mathbf{v}\right)+{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right);\sigma {\widehat{\mathbf{u}}}_{\omega }^{n-1},\mathbf{v}\right)$$

$$+\mu \sum _{i=1}^{d-1}{\int }_I\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right)\left({\mathbf{r}}_3^n·{\mathbf{m}}_i\right)\left(\mathbf{v}·{\mathbf{m}}_i\right)ds,$$

(22)

$${\mathcal{A}}_p\left({\mathbb{K}}_{\omega };{\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)={\mathcal{A}}_p\left(\overline{\mathbb{K}}+\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}}\right);{\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)$$

$$={\mathcal{A}}_p\left(\overline{\mathbb{K}};{\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)+{\mathcal{A}}_p\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}};\sigma {\widehat{\mathbf{u}}}_{\omega }^{n-1},\widehat{\mathbf{v}}\right)+g{\left(\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}}\right)\nabla r_5^n,\nabla \psi \right)}_{{\mathrm{\Omega }}_p},$$

(23)

$$\mathcal{I}\left({\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)=\mathcal{I}\left(\sigma {\widehat{\mathbf{u}}}_{\omega }^{n-1},\widehat{\mathbf{v}}\right)+{\mathcal{H}}_I\left(\mathbf{v},r_2^n\right)-{\mathcal{H}}_I\left({\mathbf{r}}_4^n,\psi \right),$$

(24)

where

$$\sigma {\widehat{\mathbf{u}}}_{\omega }^n={\left[\sigma {\mathbf{u}}_{\omega }^n,\sigma {\varphi }_{\omega }^n\right]}^{\top },\sigma {\mathbf{u}}_{\omega }^n=3{\mathbf{u}}_{\omega }^n-3{\mathbf{u}}_{\omega }^{n-1}+{\mathbf{u}}_{\omega }^{n-2},\sigma {\varphi }_{\omega }^n=3{\varphi }_{\omega }^n-3{\varphi }_{\omega }^{n-1}+{\varphi }_{\omega }^{n-2},$$

$${\overline{\gamma }}_i=\frac{1}{N_{\omega }}\sum _{\omega =1}^{N_{\omega }}{\gamma }_{i,\omega },\overline{\mathbb{K}}=\frac{1}{N_{\omega }}\sum _{\omega =1}^{N_{\omega }}{\mathbb{K}}_{\omega },$$

$$r_2^n=O\left({\tau }^3\right),{\mathbf{r}}_3^n={\left[O\left({\tau }^3\right),O\left({\tau }^3\right)\right]}^{\top },{\mathbf{r}}_4^n={\left[O\left({\tau }^3\right),O\left({\tau }^3\right)\right]}^{\top },r_5^n=O\left({\tau }^3\right).$$

By using the approximations (20)–(24) and introducing the notation

$${\widehat{\mathbf{u}}}_{\omega }^{n-1/2}=\frac{1}{2}\left({\widehat{\mathbf{u}}}_{\omega }^n+{\widehat{\mathbf{u}}}_{\omega }^{n-1}\right),$$

(25)

we rewrite Problem 2 in the following way.

Problem 3.

Let $\left({\widehat{\mathbf{u}}}_{\omega }^k,p_{\omega }^k\right)\in \mathbf{X}\times Q_f$, $k=0,1,\dots ,n-1$ be known; in particular, ${\widehat{\mathbf{u}}}_{\omega }^0={\left[{\mathbf{u}}_0,{\varphi }_0\right]}^{\top }$. Find $\left({\widehat{\mathbf{u}}}_{\omega }^n,p_{\omega }^n\right)\in \mathbf{X}\times Q_f$, $n\ge 1$ satisfying the following identities for any $\left(\widehat{\mathbf{v}},q\right)\in \mathbf{X}\times Q_f$:

Case I:  when $n=1,2$:

$$\begin{array}{cc}& {\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega }^{n-1/2},\widehat{\mathbf{v}}\right)+\mathcal{A}\left(\overline{\mathbb{K}};{\widehat{\mathbf{u}}}_{\omega }^{n-1/2},\widehat{\mathbf{v}}\right)-\mathcal{B}\left(\widehat{\mathbf{v}},p_{\omega }^{n-1/2}\right)+\mathcal{B}\left({\widehat{\mathbf{u}}}_{\omega }^{n-1/2},q\right)\hfill \\ & +\mathcal{I}\left({\widehat{\mathbf{u}}}_{\omega }^{n-1/2},\widehat{\mathbf{v}}\right)+{\mathcal{A}}_I\left({\gamma }_{\omega };{\widehat{\mathbf{u}}}_{\omega }^{n-1/2},\widehat{\mathbf{v}}\right)\hfill \\ & =\left({\widehat{\mathbf{f}}}_{\omega }^{n-1/2},\widehat{\mathbf{v}}\right)+{\left({\mathbf{r}}_{11}^{n-1/2},\mathbf{v}\right)}_{{\mathrm{\Omega }}_f}+{\left(r_{12}^{n-1/2},\psi \right)}_{{\mathrm{\Omega }}_p}\hfill \\ & +{\mathcal{H}}_I\left(\mathbf{v},r_2^{n-1/2}\right)+{\mathcal{H}}_I\left({\mathbf{r}}_4^{n-1/2},\psi \right)\hfill \\ & -\mu \sum _{i=1}^{d-1}{\int }_I{\gamma }_{i,\omega }\left({\mathbf{r}}_3^{n-1/2}·{\mathbf{m}}_i\right)\left(\mathbf{v}·{\mathbf{m}}_i\right)ds+g{\left({\mathbb{K}}_{\omega }\nabla r_5^{n-1/2},\nabla \psi \right)}_{{\mathrm{\Omega }}_p},\hfill \end{array}$$

(26)

Case II:  when $n\ge 3$:

$$\begin{array}{cc}& {\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)+\mathcal{A}\left(\overline{\mathbb{K}};{\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)-\mathcal{B}\left(\widehat{\mathbf{v}},p_{\omega }^n\right)+\mathcal{B}\left({\widehat{\mathbf{u}}}_{\omega }^n,q\right)\hfill \\ & +{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}{\overline{\gamma }}_i;{\widehat{\mathbf{u}}}_{\omega }^n,\widehat{\mathbf{v}}\right)+{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right);\sigma {\widehat{\mathbf{u}}}_{\omega }^{n-1},\widehat{\mathbf{v}}\right)\hfill \\ & +\mathcal{I}\left(\sigma {\widehat{\mathbf{u}}}_{\omega }^{n-1},\widehat{\mathbf{v}}\right)+{\mathcal{A}}_p\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}};\sigma {\widehat{\mathbf{u}}}_{\omega }^{n-1},\widehat{\mathbf{v}}\right)\hfill \\ & =\left({\widehat{\mathbf{f}}}_{\omega }^n,\widehat{\mathbf{v}}\right)+{\left({\mathbf{r}}_{11}^n,\mathbf{v}\right)}_{{\mathrm{\Omega }}_f}+{\left(r_{12}^n,\psi \right)}_{{\mathrm{\Omega }}_p}+{\mathcal{H}}_I\left(\mathbf{v},r_2^n\right)+{\mathcal{H}}_I\left({\mathbf{r}}_4^n,\psi \right)\hfill \\ & -\mu \sum _{i=1}^{d-1}{\int }_I\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right)\left({\mathbf{r}}_3^n·{\mathbf{m}}_i\right)\left(\mathbf{v}·{\mathbf{m}}_i\right)ds+g{\left(\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}}\right)\nabla r_5^n,\nabla \psi \right)}_{{\mathrm{\Omega }}_p}.\hfill \end{array}$$

(27)

Now we discretize the problem with respect to spatial variables. To this end, we introduce the partitions ${\mathfrak{T}}_{fh}$ and ${\mathfrak{T}}_{ph}$ in ${\mathrm{\Omega }}_f$ and ${\mathrm{\Omega }}_p$, respectively, that may consist of triangles, quadrilaterals, prisms, or parallelepipeds, and let ${\mathfrak{T}}_h={\mathfrak{T}}_{fh}\cup {\mathfrak{T}}_{ph}$. Denote by $X_{fh}\subset X_f$, $Q_{fh}\subset Q_f$, $Q_{ph}\subset Q_p$ finite element spaces with the following approximation properties:

$$\begin{array}{cccc}& \underset{{\mathbf{v}}_h\in X_{fh}}{inf}{∥\mathbf{v}-{\mathbf{v}}_h∥}_{{\mathrm{\Omega }}_f}\le C{h}^{l+1}{∥\mathbf{v}∥}_{{\mathbf{H}}^{l+1}\left({\mathrm{\Omega }}_f\right)},\hfill & \hfill & \forall \mathbf{v}\in {\left(H^{l+1}\left({\mathrm{\Omega }}_f\right)\right)}^d,\hfill \\ & \underset{{\mathbf{v}}_h\in X_{fh}}{inf}{∥\nabla \left(\mathbf{v}-{\mathbf{v}}_h\right)∥}_{{\mathrm{\Omega }}_f}\le C{h}^l{∥\mathbf{v}∥}_{{\mathbf{H}}^{l+1}\left({\mathrm{\Omega }}_f\right)},\hfill & & \forall \mathbf{v}\in {\left(H^{l+1}\left({\mathrm{\Omega }}_f\right)\right)}^d,\hfill \\ & \underset{q_h\in Q_{fh}}{inf}{∥q-q_h∥}_{{\mathrm{\Omega }}_f}\le C{h}^l{∥q∥}_{H^l\left({\mathrm{\Omega }}_f\right)},\hfill & & \forall q\in H^l\left({\mathrm{\Omega }}_f\right),\hfill \\ & \underset{{\psi }_h\in Q_{ph}}{inf}{∥\psi -{\psi }_h∥}_{{\mathrm{\Omega }}_p}\le C{h}^{m+1}{∥\psi ∥}_{H^{m+1}\left({\mathrm{\Omega }}_p\right)},\hfill & & \forall \psi \in H^{m+1}\left({\mathrm{\Omega }}_p\right),\hfill \\ & \underset{{\psi }_h\in Q_{ph}}{inf}{∥\nabla \left(\psi -{\psi }_h\right)∥}_{{\mathrm{\Omega }}_p}\le C{h}^m{∥\psi ∥}_{H^{m+1}\left({\mathrm{\Omega }}_p\right)},\hfill & & \forall \psi \in H^{m+1}\left({\mathrm{\Omega }}_p\right),\hfill \end{array}$$

(28)

and let ${\mathbf{X}}_h=X_{fh}\times Q_{ph}$.

By discarding the error terms in (26) and (27), the following discrete problem can be obtained.

Problem 4.

Let $\left({\widehat{\mathbf{u}}}_{\omega ,h}^k,p_{\omega ,h}^k\right)\in {\mathbf{X}}_h\times Q_{fh}$, $k=0,1,\dots ,n-1$ be known; in particular, let ${\widehat{\mathbf{u}}}_{\omega ,h}^0$ consist of $L^2$-projections of the initial conditions (9), (10). Find $\left({\widehat{\mathbf{u}}}_{\omega ,h}^n,p_{\omega ,h}^n\right)\in {\mathbf{X}}_h\times Q_{fh}$, $n\ge 1$ satisfying the following identities for any $\left({\widehat{\mathbf{v}}}_h,q_h\right)\in {\mathbf{X}}_h\times Q_{fh}$:

Case I:  when $n=1,2$:

$$\begin{array}{cc}& {\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{v}}}_h\right)+\mathcal{A}\left(\overline{\mathbb{K}};{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{v}}}_h\right)-\mathcal{B}\left({\widehat{\mathbf{v}}}_h,p_{\omega ,h}^{n-1/2}\right)+\mathcal{B}\left({\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},q_h\right)\hfill \\ & +\mathcal{I}\left({\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{v}}}_h\right)+{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}{\gamma }_{i,\omega };{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{v}}}_h\right)\hfill \\ & =\left({\widehat{\mathbf{f}}}_{\omega }^{n-1/2},{\widehat{\mathbf{v}}}_h\right),\hfill \end{array}$$

(29)

Case II:  when $n\ge 3$:

$$\begin{array}{cc}& {\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega ,h}^n,{\widehat{\mathbf{v}}}_h\right)+\mathcal{A}\left(\overline{\mathbb{K}};{\widehat{\mathbf{u}}}_{\omega ,h}^n,{\widehat{\mathbf{v}}}_h\right)-\mathcal{B}\left({\widehat{\mathbf{v}}}_h,p_{\omega ,h}^n\right)+\mathcal{B}\left({\widehat{\mathbf{u}}}_{\omega ,h}^n,q_h\right)\hfill \\ & +{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}{\overline{\gamma }}_i;{\widehat{\mathbf{u}}}_{\omega ,h}^n,{\widehat{\mathbf{v}}}_h\right)+{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right);\sigma {\widehat{\mathbf{u}}}_{\omega ,h}^{n-1},{\widehat{\mathbf{v}}}_h\right)\hfill \\ & +\mathcal{I}\left(\sigma {\widehat{\mathbf{u}}}_{\omega ,h}^{n-1},{\widehat{\mathbf{v}}}_h\right)+{\mathcal{A}}_p\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}};\sigma {\widehat{\mathbf{u}}}_{\omega ,h}^{n-1},{\widehat{\mathbf{v}}}_h\right)=\left({\widehat{\mathbf{f}}}_{\omega }^n,{\widehat{\mathbf{v}}}_h\right).\hfill \end{array}$$

(30)

As in the original papers [3,17], the advantage of the proposed discrete scheme (29) and (30) is that the stiffness matrix is the same for all realizations $\omega =1,2,\dots ,N_{\omega }$. When implementing the algorithm, only the right-hand sides of systems of linear algebraic equations are different. This makes it possible to use efficient methods for solving systems of equations and save computational resources. Further implementation details will be covered in Section 2.4 and Section 3.2.

2.3. Theoretical Analysis of the Numerical Method

2.3.1. Preliminaries

This section presents auxiliary results and the assumptions under which we obtain a priori estimates of the stability and convergence for the numerical scheme (29) and (30). First, let us present the following well-known inequalities, which will be used in obtaining the main results:

$${∥\mathbf{v}∥}_I\le C_{\mathrm{tr}}{∥\mathbf{v}∥}_{{\mathrm{\Omega }}_f}^{1/2}{∥\nabla \mathbf{v}∥}_{{\mathrm{\Omega }}_f}^{1/2},{∥\psi ∥}_I\le C_{\mathrm{tr}}{∥\psi ∥}_{{\mathrm{\Omega }}_p}^{1/2}{∥\nabla \psi ∥}_{{\mathrm{\Omega }}_p}^{1/2},$$

(31)

$${∥\mathbf{v}∥}_{{\mathrm{\Omega }}_f}\le C_{P,f}{∥\nabla \mathbf{v}∥}_{{\mathrm{\Omega }}_f},{∥\psi ∥}_{{\mathrm{\Omega }}_p}\le C_{P,p}{∥\nabla \psi ∥}_{{\mathrm{\Omega }}_p},$$

(32)

$$ab\le \epsilon a^2+\frac{1}{4\epsilon }{b}^2,a,b>0.$$

(33)

Assumption 1.

Problem 1 has a unique solution having a sufficient number of derivatives required to perform the analysis.

Assumption 2.

Let $k_{min}\left(x\right)$ and ${\overline{k}}_{min}\left(x\right)$ be the minimum eigenvalue of tensors and vectors. Let us assume that the following conditions hold:

$$\begin{array}{cc}& {\overline{\gamma }}_i^{min}=\underset{x\in I}{min}{\overline{\gamma }}_i\left(x\right),\hfill \end{array}$$

(34)

$$\begin{array}{cc}& {\overline{\gamma }}_i^{max}=\underset{x\in I}{max}{\overline{\gamma }}_i\left(x\right),\hfill \end{array}$$

(35)

$$\begin{array}{cc}& {\overline{\gamma }}^{max}=\underset{i=1,2,\dots ,d-1}{max}{\overline{\gamma }}_i^{max},\hfill \end{array}$$

(36)

$$\begin{array}{cc}& {\gamma }_{i,\omega }^{\prime \max }=\underset{x\in I}{max}\left|{\gamma }_{i,\omega }\left(x\right)-{\overline{\gamma }}_i\left(x\right)\right|,\hfill \end{array}$$

(37)

$$\begin{array}{cc}& {\gamma }_{\omega }^{\prime \max }=\underset{i=1,2,\dots ,d-1}{max}{\gamma }_{i,\omega }^{\prime \max }=\underset{i=1,2,\dots ,d-1}{max}\underset{x\in I}{max}\left|{\gamma }_{i,\omega }\left(x\right)-{\overline{\gamma }}_i\left(x\right)\right|,\hfill \end{array}$$

(38)

$$\begin{array}{cc}& {\kappa }_{\omega ,max}^{\prime }=\underset{x\in {\mathrm{\Omega }}_p}{max}\left|{\mathbb{K}}_{\omega }-\overline{\mathbb{K}}\right|,\hfill \end{array}$$

(39)

$$\begin{array}{cc}& \mu \left(1-ϵ-{\overline{\gamma }}_i^{max}{C}_{\mathrm{tr}}{C}_{P,f}\right)\ge {\mu }_0>0,\hfill \end{array}$$

(40)

where $ϵ>0$, and $C_{\mathrm{tr}}$ and $C_{P,f}$ are constants from (31) and (32).

The following Gronwall lemma will be used several times in obtaining the main results.

Lemma 4.

Let non-negative sequences $\left\{y_n\right\}$, $\left\{z_n\right\}$, $\left\{f_n\right\}$ satisfy the inequalities ${\displaystyle y_n\le f_n+\sum _{i=0}^{n-1}{z}_i{y}_i}$. Then, the following inequality holds:

$$y_n\le f_n+\sum _{i=0}^{n-1}{f}_i{z}_{i}exp\left(\sum _{j=i+1}^{n-1}{z}_j\right).$$

2.3.2. Stability of the Numerical Scheme

Now we turn to the proof of the stability of the proposed numerical method (29) and (30) similarly to [43]. Let us first consider the case of $n=1,2$. The symbol C with possible subscripts will denote a generic positive constant that may take different values case by case.

Lemma 5.

Let $\left({\widehat{\mathbf{u}}}_{\omega ,h}^n,p_{\omega ,h}^n\right)\in {\mathbf{X}}_h\times Q_{fh}$ be the solution of Problem 4. Then, the following inequality holds under Assumption 2:

$${\mathcal{E}}_{\alpha ,\beta }^n+{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}∥}^2\le C_{\alpha ,\beta }\left({\mathcal{E}}_{\alpha ,\beta }^0+\sum _{s=1}^n{∥{\widehat{\mathbf{f}}}_{f,\omega }^{s-1/2}∥}^2\right),n=1,2,$$

which implies the stability of the numerical method (29) and (30) with respect to the initial data and the right-hand sides of the equations on the first two time layers, where

$${\mathcal{E}}_{\alpha ,\beta }^n={\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2,$$

and $C_{\alpha ,\beta }$ is a positive constant that only depends on the orders of fractional derivatives, α and β.

Proof. 

Choose $\left({\widehat{\mathbf{v}}}_h,q_h\right)=\left({\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},p_{\omega ,h}^{n-1/2}\right)$ in the identity (29) to obtain

$${\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}\right)+\mathcal{A}\left(\overline{\mathbb{K}};{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}\right)+{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}{\gamma }_{i,\omega };{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}\right)$$

$$-\mathcal{I}\left({\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}\right)=\left({\widehat{\mathbf{f}}}_{\omega }^{n-1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}\right).$$

(41)

Let us estimate the terms of (41). Due to the expansion (19) and the notation (25), the first term on the left-hand side of (41) is estimated as

$${\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega ,h}^{1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{1/2}\right)=\frac{1}{2}\left[{\left({\overline{\Delta }}_{\tau }^{\alpha }{\mathbf{u}}_{\omega ,h}^1,{\mathbf{u}}_{\omega ,h}^{1/2}\right)}_{{\mathrm{\Omega }}_f}+{\left(gS{\overline{\Delta }}_{\tau }^{\beta }{\varphi }_{\omega ,h}^1,{\varphi }_{\omega ,h}^{1/2}\right)}_{{\mathrm{\Omega }}_p}\right]$$

$$=\frac{1}{4}\left[\frac{{\tau }^{-\alpha }{\overline{\delta }}_{1,1}^{\alpha }}{{\overline{\varsigma }}_0^{\left(\alpha \right)}}\left({∥{\mathbf{u}}_{\omega ,h}^1∥}_{{\mathrm{\Omega }}_f}^2-{∥{\mathbf{u}}_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2\right)+gS\frac{{\tau }^{-\beta }{\overline{\delta }}_{1,1}^{\beta }}{{\overline{\varsigma }}_0^{\left(\beta \right)}}\left({∥{\varphi }_{\omega ,h}^1∥}_{{\mathrm{\Omega }}_p}^2-{∥{\varphi }_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_p}^2\right)\right]$$

(42)

and

$${\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega ,h}^{3/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{3/2}\right)=\frac{1}{2}\left[{\left({\overline{\Delta }}_{\tau }^{\alpha }{\mathbf{u}}_{\omega ,h}^2,{\mathbf{u}}_{\omega ,h}^{3/2}\right)}_{{\mathrm{\Omega }}_f}+{\left({\overline{\Delta }}_{\tau }^{\alpha }{\mathbf{u}}_{\omega ,h}^1,{\mathbf{u}}_{\omega ,h}^{3/2}\right)}_{{\mathrm{\Omega }}_f}\right]$$

$$+\frac{1}{2}\left[{\left(gS{\overline{\Delta }}_{\tau }^{\beta }{\varphi }_{\omega ,h}^2,{\varphi }_{\omega ,h}^{3/2}\right)}_{{\mathrm{\Omega }}_p}+{\left(gS{\overline{\Delta }}_{\tau }^{\beta }{\varphi }_{\omega ,h}^1,{\varphi }_{\omega ,h}^{3/2}\right)}_{{\mathrm{\Omega }}_p}\right]$$

$$=\frac{{\tau }^{-\alpha }}{4{\overline{\varsigma }}_0^{\left(\alpha \right)}}\left[{\overline{\delta }}_{2,2}^{\alpha }{\left({\mathbf{u}}_{\omega ,h}^2-{\mathbf{u}}_{\omega ,h}^1,{\mathbf{u}}_{\omega ,h}^2+{\mathbf{u}}_{\omega ,h}^1\right)}_{{\mathrm{\Omega }}_f}+\left({\overline{\delta }}_{2,1}^{\alpha }+{\overline{\delta }}_{1,1}^{\alpha }\right){\left({\mathbf{u}}_{\omega ,h}^1-{\mathbf{u}}_{\omega ,h}^0,{\mathbf{u}}_{\omega ,h}^{3/2}\right)}_{{\mathrm{\Omega }}_f}\right]$$

$$+\frac{gS{\tau }^{-\beta }}{4{\overline{\varsigma }}_0^{\left(\beta \right)}}\left[{\overline{\delta }}_{2,2}^{\beta }{\left({\varphi }_{\omega ,h}^2-{\varphi }_{\omega ,h}^1,{\varphi }_{\omega ,h}^2+{\varphi }_{\omega ,h}^1\right)}_{{\mathrm{\Omega }}_p}+\left({\overline{\delta }}_{2,1}^{\beta }+{\overline{\delta }}_{1,1}^{\beta }\right){\left({\varphi }_{\omega ,h}^1-{\varphi }_{\omega ,h}^0,{\varphi }_{\omega ,h}^{3/2}\right)}_{{\mathrm{\Omega }}_p}\right]$$

$$\ge \frac{{\tau }^{-\alpha }}{4{\overline{\varsigma }}_0^{\left(\alpha \right)}}\left[\frac{1}{2}{\overline{\delta }}_{2,2}^{\alpha }{∥{\mathbf{u}}_{\omega ,h}^2∥}_{{\mathrm{\Omega }}_f}^2-C_{\alpha }\left({∥{\mathbf{u}}_{\omega ,h}^1∥}_{{\mathrm{\Omega }}_f}^2+{∥{\mathbf{u}}_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2\right)\right]$$

$$+\frac{gS{\tau }^{-\beta }}{4{\overline{\varsigma }}_0^{\left(\beta \right)}}\left[\frac{1}{2}{\overline{\delta }}_{2,2}^{\beta }{∥{\varphi }_{\omega ,h}^2∥}_{{\mathrm{\Omega }}_p}^2-C_{\beta }\left({∥{\varphi }_{\omega ,h}^1∥}_{{\mathrm{\Omega }}_f}^2-{∥{\varphi }_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2\right)\right],$$

(43)

where $C_{\nu }$ is a positive constant depending only on $\nu$. For the second term on the left-hand side of (41) we have:

$$\mathcal{A}\left(\overline{\mathbb{K}};{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}\right)\ge \mu {∥\nabla {\mathbf{u}}_{\omega ,h}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2+g{k}_{min}{∥\nabla {\varphi }_{\omega ,h}^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2.$$

(44)

Further, by virtue of Assumption 2, relation (14), and inequalities (31) and (32), we arrive at the inequality

$$\begin{array}{cc}\hfill {\mathcal{A}}_I\left(\sum _{i=1}^{d-1}{\gamma }_{i,\omega };{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}\right)& \le \mu {\overline{\gamma }}_i^{max}\sum _{i=1}^{d-1}{∥{\mathbf{u}}_{\omega ,h}^{n-1/2}·{\mathbf{m}}_i∥}_I^2\le \mu {\overline{\gamma }}_i^{max}{∥{\mathbf{u}}_{\omega ,h}^{n-1/2}∥}_I^2\hfill \\ \hfill & \le \mu {\overline{\gamma }}_i^{max}{C}_{\mathrm{tr}}{C}_{P,f}{∥\nabla {\mathbf{u}}_{\omega ,h}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2.\hfill \end{array}$$

(45)

Finally, using (32) and (33) yields

$$\left({\widehat{\mathbf{f}}}_{\omega }^{n-1/2},{\widehat{\mathbf{u}}}_{\omega ,h}^{n-1/2}\right)\le C_{P,f}{∥{\mathbf{f}}_{f,\omega }^{n-1/2}∥}_{{\mathrm{\Omega }}_f}{∥\nabla {\mathbf{u}}_{\omega ,h}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}+C_{P,p}g{∥f_{p,\omega }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}{∥\nabla {\varphi }_{\omega ,h}^{n-1/2}∥}_{{\mathrm{\Omega }}_p}$$

$$\le ϵ\mu {∥\nabla {\mathbf{u}}_{\omega ,h}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2+ϵg{k}_{min}{∥\nabla {\varphi }_{\omega ,h}^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2+C_1\left({∥{\mathbf{f}}_{f,\omega }^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2+{∥f_{p,\omega }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2\right),$$

(46)

where $ϵ$ is a sufficiently small number and $C_1$ is a constant depending on the domains ${\mathrm{\Omega }}_f$ and ${\mathrm{\Omega }}_p$.

Making use of the inequalities (42)–(46) and Assumption (40), we obtain from the identity (41) that

$${\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^1∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^1∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^{1/2}∥}^2$$

$$+C_{\alpha ,\beta }\left({\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_p}^2+{∥{\mathbf{f}}_{f,\omega }^{1/2}∥}_{{\mathrm{\Omega }}_f}^2+{∥f_{p,\omega }^{1/2}∥}_{{\mathrm{\Omega }}_p}^2\right)$$

(47)

and

$${\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^2∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^2∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^{3/2}∥}^2\le {\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^1∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^1∥}_{{\mathrm{\Omega }}_f}^2$$

$$+{\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2+C_{\alpha ,\beta }\left({∥{\mathbf{f}}_{f,\omega }^{3/2}∥}_{{\mathrm{\Omega }}_f}^2+{∥f_{p,\omega }^{3/2}∥}_{{\mathrm{\Omega }}_p}^2\right).$$

(48)

By estimating the first two terms on the right-hand side of (48) with the use of inequality (47), we obtain

$${\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^2∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^2∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^{3/2}∥}^2$$

$$\le C_{\alpha ,\beta }\left({\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2+\sum _{s=1}^2\left({∥{\mathbf{f}}_{f,\omega }^{s-1/2}∥}_{{\mathrm{\Omega }}_f}^2+{∥f_{p,\omega }^{s-1/2}∥}_{{\mathrm{\Omega }}_p}^2\right)\right).$$

(49)

Thus, combining (47) and (49) yields the assertion of the lemma.    □

Before we obtain the main result, let us prove the following auxiliary lemma that will be used while estimating the terms with fractional derivatives in the case of $n\ge 3$.

Lemma 6.

The following inequality holds for the bilinear form ${\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}^n,{\widehat{\mathbf{u}}}^n\right)$, $n\ge 3$:

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}^n,{\widehat{\mathbf{u}}}^n\right)& \ge {\tau }^{-\alpha }\left[{∥{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2-{\chi }_{\alpha }\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\alpha \right)}\right|{∥{\mathbf{u}}^{n-s}∥}_{{\mathrm{\Omega }}_f}^2\right]\hfill \\ \hfill & +{\tau }^{-\beta }\left[{∥{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2-{\chi }_{\beta }\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\beta \right)}\right|{∥{\varphi }^{n-s}∥}_{{\mathrm{\Omega }}_p}^2\right],\hfill \end{array}$$

(50)

where ${\chi }_{\alpha }=\frac{{\left({\lambda }_0^{\left(\alpha \right)}\right)}^2}{2{\varsigma }_0^{\left(\alpha \right)}\left({\lambda }_0^{\left(\alpha \right)}-{\varsigma }_0^{\left(\alpha \right)}\right)}$, ${\chi }_{\beta }=\frac{{\left(gS{\lambda }_0^{\left(\beta \right)}\right)}^2}{2{\varsigma }_0^{\left(\beta \right)}\left(gS{\lambda }_0^{\left(\beta \right)}-{\varsigma }_0^{\left(\beta \right)}\right)}$.

Proof. 

Choose $\widehat{\mathbf{v}}={\widehat{\mathbf{u}}}^n$ in (21) and use the expansion (18) to obtain

$$\begin{array}{cc}& {\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}^n,{\widehat{\mathbf{u}}}^n\right)=\frac{{\lambda }_0^{\left(\alpha \right)}{\tau }^{-\alpha }}{{\varsigma }_0^{\left(\alpha \right)}}{\left({\mathbf{u}}^n-\sum _{s=1}^n{\delta }_{n,n-s}^{\left(\alpha \right)}{\mathbf{u}}^{n-s},{\mathbf{u}}^n\right)}_{{\mathrm{\Omega }}_f}\hfill \\ & +\frac{gS{\lambda }_0^{\left(\beta \right)}{\tau }^{-\beta }}{{\varsigma }_0^{\left(\beta \right)}}{\left({\varphi }^n-\sum _{s=1}^n{\delta }_{n,n-s}^{\left(\beta \right)}{\varphi }^{n-s},{\varphi }^n\right)}_{{\mathrm{\Omega }}_p}.\hfill \end{array}$$

Further, applying the Cauchy inequality and the inequality (33) yields

$${\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}^n,{\widehat{\mathbf{u}}}^n\right)\ge \frac{{\lambda }_0^{\left(\alpha \right)}{\tau }^{-\alpha }}{{\varsigma }_0^{\left(\alpha \right)}}\left[{∥{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2-\frac{1}{4{\epsilon }_{\alpha }}\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\alpha \right)}\right|{∥{\mathbf{u}}^{n-s}∥}_{{\mathrm{\Omega }}_f}^2-{\epsilon }_{\alpha }{∥{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\alpha \right)}\right|\right]$$

$$+\frac{gS{\lambda }_0^{\left(\beta \right)}{\tau }^{-\beta }}{{\varsigma }_0^{\left(\beta \right)}}\left[{∥{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2-\frac{1}{4{\epsilon }_{\beta }}\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\beta \right)}\right|{∥{\varphi }^{n-s}∥}_{{\mathrm{\Omega }}_p}^2-{\epsilon }_{\beta }{∥{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\beta \right)}\right|\right].$$

(51)

Since it follows from Lemma 2 that ${\delta }_{n,s}^{\left(\nu \right)}>0$ for $s\ne n-2$, and $\left|{\delta }_{n,n-2}^{\left(\nu \right)}\right|<\frac{1}{2}$, it is easy to see that ${\displaystyle \sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\nu \right)}\right|<2}$. In this case, it suffices to choose ${\epsilon }_{\alpha }=\frac{1}{2}\left(1-\frac{{\varsigma }_0^{\left(\alpha \right)}}{{\lambda }_0^{\left(\alpha \right)}}\right)$, ${\epsilon }_{\beta }=\frac{1}{2}\left(1-\frac{{\varsigma }_0^{\left(\beta \right)}}{gS{\lambda }_0^{\left(\beta \right)}}\right)$ in (51) to obtain an assertion of the lemma.    □

Now we are ready to prove the stability of the proposed numerical method.

Theorem 1.

Let $\left({\widehat{\mathbf{u}}}_{\omega ,h}^n,p_{\omega ,h}^n\right)\in {\mathbf{X}}_{fh}\times Q_{ph}$ be the solution of Problem 4. Then, the following inequality holds under Assumption 2:

$${\mathcal{E}}_{\alpha ,\beta }^n+{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^n∥}^2\le C_{\alpha ,\beta }\left({\mathcal{E}}_{\alpha ,\beta }^0+\sum _{s=0}^n{∥{\widehat{\mathbf{f}}}_{f,\omega }^s∥}^2\right),n\ge 1,$$

which implies the stability of the numerical method with respect to the initial data and the right-hand sides of Equations (1) and (4), where

$${\mathcal{E}}_{\alpha ,\beta }^n={\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2,$$

and $C_{\alpha ,\beta }$ is a positive constant that only depends on the orders of fractional derivatives, α and β.

Proof. 

Choose $\left({\widehat{\mathbf{v}}}_h,q_h\right)=\left({\widehat{\mathbf{u}}}_{\omega ,h}^n,p_{\omega ,h}^n\right)$ in the identity (30) to obtain

$$\begin{array}{cc}\hfill & \sum _{k=1}^7{\mathcal{T}}_k\equiv {\mathcal{D}}_{\tau }^{\alpha ,\beta }\left({\widehat{\mathbf{u}}}_{\omega ,h}^n,{\widehat{\mathbf{u}}}_{\omega ,h}^n\right)+\mathcal{A}\left(\overline{\mathbb{K}};{\widehat{\mathbf{u}}}_{\omega ,h}^n,{\widehat{\mathbf{u}}}_{\omega ,h}^n\right)+{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}{\overline{\gamma }}_i;{\widehat{\mathbf{u}}}_{\omega ,h}^n,{\widehat{\mathbf{u}}}_{\omega ,h}^n\right)\hfill \\ \hfill & +{\mathcal{A}}_I\left(\sum _{i=1}^{d-1}\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right);\sigma {\widehat{\mathbf{u}}}_{\omega ,h}^{n-1},{\widehat{\mathbf{u}}}_{\omega ,h}^n\right)-\mathcal{I}\left(\sigma {\widehat{\mathbf{u}}}_{\omega ,h}^{n-1},{\widehat{\mathbf{u}}}_{\omega ,h}^n\right)\hfill \\ \hfill & +{\mathcal{A}}_p\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}};\sigma {\widehat{\mathbf{u}}}_{\omega ,h}^{n-1},{\widehat{\mathbf{u}}}_{\omega ,h}^n\right)-\left({\widehat{\mathbf{f}}}_{\omega ,h}^n,{\widehat{\mathbf{u}}}_{\omega ,h}^n\right)=0.\hfill \end{array}$$

(52)

Let us estimate the terms ${\mathcal{T}}_1$, ${\mathcal{T}}_2$, …, ${\mathcal{T}}_7$ similarly to Lemma 5 with the use of Lemma 6, Assumption 2, and inequalities (31)–(33):

$${\mathcal{T}}_1\ge {\tau }^{-\alpha }\left[{∥{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2-{\chi }_{\alpha }\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\alpha \right)}\right|{∥{\mathbf{u}}_{\omega ,h}^{n-s}∥}_{{\mathrm{\Omega }}_f}^2\right]+{\tau }^{-\beta }\left[{∥{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2-{\chi }_{\beta }\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\beta \right)}\right|{∥{\varphi }_{\omega ,h}^{n-s}∥}_{{\mathrm{\Omega }}_p}^2\right],$$

$${\mathcal{T}}_2\ge \mu {∥\nabla {\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+g{k}_{min}{∥\nabla {\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2,$$

$${\mathcal{T}}_3\le \mu {\overline{\gamma }}^{max}{C}_{\mathrm{tr}}{C}_{P,f}{∥\nabla {\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2,$$

$$\begin{array}{cc}\hfill {\mathcal{T}}_4+{\mathcal{T}}_5& \le \mu {\gamma }_{\omega }^{\prime \max }\left|\sum _{i=1}^{d-1}{\int }_I\left(\sigma {\mathbf{u}}_{\omega ,h}^{n-1}·{\mathbf{m}}_i\right)\left({\mathbf{u}}_{\omega ,h}^n·{\mathbf{m}}_i\right)ds\right|\hfill \\ \hfill & +g{∥\sigma {\varphi }_h^{n-1}∥}_I{∥{\mathbf{u}}_{\omega ,h}^n·{\mathbf{n}}_f∥}_I+g{∥{\varphi }_h^n∥}_I{∥\sigma {\mathbf{u}}_{\omega ,h}^{n-1}·{\mathbf{n}}_f∥}_I\hfill \\ \hfill & \le \mu {\gamma }_{\omega }^{\prime \max }{∥\sigma {\mathbf{u}}_{\omega ,h}^{n-1}∥}_I{∥{\mathbf{u}}_{\omega ,h}^n∥}_I+g{∥\sigma {\varphi }_h^{n-1}∥}_I{∥{\mathbf{u}}_{\omega ,h}^n∥}_I+g{∥{\varphi }_h^n∥}_I{∥\sigma {\mathbf{u}}_{\omega ,h}^{n-1}∥}_I\hfill \\ \hfill & \le \frac{\mu ϵ}{2C_{\mathrm{tr}}{C}_{P,f}}{∥{\mathbf{u}}_{\omega ,h}^n∥}_I^2+\frac{g{k}_{min}ϵ}{3C_{\mathrm{tr}}{C}_{P,p}}{∥{\varphi }_h^n∥}_I+C{∥\sigma {\mathbf{u}}_{\omega ,h}^{n-1}∥}_I^2+C\left({∥\sigma {\mathbf{u}}_{\omega ,h}^{n-1}∥}_I^2+{∥\sigma {\varphi }_{\omega ,h}^{n-1}∥}_I^2\right)\hfill \\ \hfill & \le \frac{1}{2}\mu ϵ{∥\nabla {\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+\frac{1}{3}g{k}_{min}ϵ{∥\nabla {\varphi }_h^n∥}_{{\mathrm{\Omega }}_p}^2+C\sum _{s=n-3}^{n-1}{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^s∥}^2,\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{T}}_6& \le g{\int }_{{\mathrm{\Omega }}_p}{\kappa }_{\omega }^{\prime }\left(x\right){\left|\nabla \left(\sigma {\varphi }_{\omega ,h}^{n-1}\right)\right|}_2{\left|\nabla {\varphi }_{\omega ,h}^n\right|}_{2}dx\hfill \\ \hfill & \le C{∥\nabla \left(\sigma {\varphi }_{\omega ,h}^{n-1}\right)∥}_{{\mathrm{\Omega }}_p}{∥\nabla {\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}\hfill \\ \hfill & \le C\sum _{s=n-3}^{n-1}{∥\nabla {\varphi }_{\omega ,h}^s∥}_{{\mathrm{\Omega }}_p}^2+\frac{1}{3}ϵg{k}_{min}{∥\nabla {\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2,\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{T}}_7& \le {∥{\mathbf{f}}_{f,\omega }^n∥}_{{\mathrm{\Omega }}_f}{∥{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}+g{∥f_{p,\omega }^n∥}_{{\mathrm{\Omega }}_p}{∥{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}\hfill \\ \hfill & \le \frac{1}{2}\mu ϵ{∥\nabla {\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+\frac{1}{3}ϵg{k}_{min}{∥\nabla {\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2+C{∥{\widehat{\mathbf{f}}}_{\omega }^n∥}^2,\hfill \end{array}$$

where ${|·|}_2$ denotes the two-norm of a vector.

Then, taking into account the obtained inequalities, we obtain from (52) that

$${\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2+\mu \left(1-ϵ-{\overline{\gamma }}^{max}{C}_{\mathrm{tr}}{C}_{P,f}\right){∥\nabla {\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2$$

$$+g{k}_{min}\left(1-ϵ\right){∥\nabla {\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2\le {\tau }^{-\alpha }{\chi }_{\alpha }\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\alpha \right)}\right|{∥{\mathbf{u}}_{\omega ,h}^{n-s}∥}_{{\mathrm{\Omega }}_f}^2$$

$$+{\tau }^{-\beta }{\chi }_{\beta }\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\beta \right)}\right|{∥{\varphi }_{\omega ,h}^{n-s}∥}_{{\mathrm{\Omega }}_p}^2+C\sum _{s=n-3}^{n-1}{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^s∥}^2+C{∥{\widehat{\mathbf{f}}}_{\omega }^n∥}^2.$$

Given the assumption (40), it follows that

$${\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^n∥}^2$$

$$\le C_{\alpha ,\beta }^{\prime }\sum _{s=0}^{n-1}\left|{\delta }_{n,s}^{\left(\alpha \right)}\right|\left({\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^s∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^s∥}_{{\mathrm{\Omega }}_p}^2\right)$$

$$+C_{\alpha ,\beta }^{\prime }\sum _{s=n-3}^{n-1}\left|{\delta }_{n,s}^{\left(\alpha \right)}\right|{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^s∥}^2+C{∥{\widehat{\mathbf{f}}}_{\omega }^n∥}^2,$$

(53)

where $C_{\alpha ,\beta }^{\prime }$ is a positive constant depending on constants ${\displaystyle c_{\alpha ,\beta }^{\prime }=max\left\{\underset{1\le s\le n}{max}\frac{\left|{\delta }_{n,s}^{\left(\alpha \right)}\right|}{\left|{\delta }_{n,s}^{\left(\beta \right)}\right|},1\right\}}$ and ${\displaystyle c_{\alpha }^{″}=max\left\{\underset{1\le s\le 3}{max}\frac{1}{\left|{\delta }_{n,n-s}^{\left(\alpha \right)}\right|},1\right\}}$. Further, by estimating the terms ${∥{\mathbf{u}}_{\omega ,h}^s∥}_{{\mathrm{\Omega }}_f}^2$ and ${∥{\varphi }_{\omega ,h}^s∥}_{{\mathrm{\Omega }}_p}^2$ for $s=1,2$ with the use of Lemma 5, we conclude that

$${\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^n∥}^2$$

$$\le C_{\alpha ,\beta }^{\prime }\sum _{s=3}^{n-1}\left|{\delta }_{n,s}^{\left(\alpha \right)}\right|\left({\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^s∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^s∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\widehat{\mathbf{u}}}_{\omega ,h}^s∥}^2\right)$$

$$+C_{\alpha }{\tau }^{-\alpha }{∥{\mathbf{u}}_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_f}^2+C_{\beta }{\tau }^{-\beta }{∥{\varphi }_{\omega ,h}^0∥}_{{\mathrm{\Omega }}_p}^2+C\left({∥{\widehat{\mathbf{f}}}_{\omega }^n∥}^2+\sum _{s=1}^2{∥{\widehat{\mathbf{f}}}_{\omega }^{n-1/2}∥}^2\right).$$

Therefore, by virtue of Lemma 4, we obtain the assertion of the theorem.    □

2.3.3. Convergence of the Numerical Scheme

Consider the projection operators ${\mathrm{\Phi }}_h:Q_p\to Q_{ph}$, ${\mathrm{\Pi }}_h:X_f\to X_{fh}$, ${\mathrm{\Xi }}_h:Q_f\to Q_{fh}$ satisfying the identities

$$\begin{array}{cccc}& \left(\nabla \left(\psi -{\mathrm{\Phi }}_h\psi \right),\nabla {\psi }_h\right)=0,\hfill & \hfill & \forall {\psi }_h\in Q_{ph},\hfill \end{array}$$

(54)

$$\begin{array}{cccc}& \left(\nabla ·\left(\mathbf{v}-{\mathrm{\Pi }}_h\mathbf{v}\right),q_h\right)=0,\hfill & & \forall q_h\in Q_{fh},\hfill \end{array}$$

(55)

$$\begin{array}{cccc}& \left(\nabla ·{\mathbf{v}}_h,q_h-{\mathrm{\Xi }}_h{q}_h\right)=0,\hfill & & \forall {\mathbf{v}}_h\in X_{fh}.\hfill \end{array}$$

(56)

Let us introduce the notation

$$\begin{array}{cc}& {\mathbf{u}}_{\omega }^n-{\mathbf{u}}_{\omega ,h}^n=\left({\mathbf{u}}_{\omega }^n-{\mathrm{\Pi }}_h{\mathbf{u}}_{\omega }^n\right)+\left({\mathrm{\Pi }}_h{\mathbf{u}}_{\omega }^n-{\mathbf{u}}_{\omega ,h}^n\right)={\xi }_{\mathbf{u}}^n+{\theta }_{\mathbf{u}}^n,\hfill \\ & p_{\omega }^n-p_{\omega ,h}^n=\left(p_{\omega }^n-{\mathrm{\Xi }}_h{p}_{\omega }^n\right)+\left({\mathrm{\Xi }}_h{p}_{\omega }^n-p_{\omega ,h}^n\right)={\xi }_p^n+{\theta }_p^n,\hfill \\ & {\varphi }_{\omega }^n-{\varphi }_{\omega ,h}^n=\left({\varphi }_{\omega }^n-{\mathrm{\Phi }}_h{\varphi }_{\omega }^n\right)+\left({\mathrm{\Phi }}_h{\varphi }_{\omega }^n-{\varphi }_{\omega ,h}^n\right)={\xi }_{\varphi }^n+{\theta }_{\varphi }^n,\hfill \end{array}$$

(57)

where the subscript $\omega$ is temporarily omitted.

To increase the convergence order in the time variable on the first and second time layers, a smaller time step in the intervals $\left(0,t_1\right)$ and $\left(t_1,t_2\right)$ is chosen. Specifically, assuming that L is the smallest integer such that ${\displaystyle L\ge {\tau }^{\frac{2}{\nu -4}}}$, we introduce the substep ${\tau }_1$ such that

$${\tau }_1=\frac{\tau }{L}={\tau }^{\frac{6-\nu }{4-\nu }},$$

(58)

where $\nu =max\left\{\alpha ,\beta \right\}$.

Let us first consider the case of $n=1,2$.

Lemma 7.

Let $\left({\mathbf{u}}_{\omega ,h}^n,p_{\omega ,h}^n,{\varphi }_{\omega ,h}^n\right)\in X_{fh}\times Q_{fh}\times Q_{ph}$ be the solution of Problem 4, and$\left({\mathbf{u}}_{\omega }^n,p_{\omega }^n,{\varphi }_{\omega }^n\right)\in X_f\times Q_f\times Q_p$ be the solution of Problem 3. Then, the following estimate holds for $n=1,2$ under Assumptions 1 and 2:

$${∥{\mathbf{u}}_{\omega }^n-{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+{∥{\varphi }_{\omega }^n-{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2+{\tau }^{\frac{6{\nu }_1-{\nu }_1^2}{4-{\nu }_1}}{∥\nabla \left({\mathbf{u}}_{\omega }^{n-1/2}-{\mathbf{u}}_{\omega ,h}^{n-1/2}\right)∥}_{{\mathrm{\Omega }}_f}^2$$

$$+{\tau }^{\frac{6{\nu }_1-{\nu }_1^2}{4-{\nu }_1}}{∥\nabla \left({\varphi }_{\omega }^{n-1/2}-{\varphi }_{\omega ,h}^{n-1/2}\right)∥}_{{\mathrm{\Omega }}_p}^2$$

$$\le C_{\alpha ,\beta }\left({\tau }^{\frac{6{\nu }_1-{\nu }_1^2}{4-{\nu }_1}}\left(h^{2l+1}+h^{2m+1}\right)+h^{2l+2}+h^{2m+2}+{\tau }^{6-{\nu }_1}\right),$$

where ${\nu }_1=max\left\{\alpha ,\beta \right\}$ and l and m are defined in (28).

Proof. 

Consider the difference between the identities (26) and (29) with the substep ${\tau }_1$. Then, use the notation (57) and choose $\left({\mathbf{v}}_h,q_h,{\psi }_h\right)=\left({\theta }_{\mathbf{u}}^{n-1/2},{\theta }_p^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)$ to obtain

$$\begin{array}{cc}& {\left({\overline{\Delta }}_{{\tau }_1}^{\alpha }{\theta }_{\mathbf{u}}^{n-1/2},{\theta }_{\mathbf{u}}^{n-1/2}\right)}_{{\mathrm{\Omega }}_f}+g{\left(S{\overline{\Delta }}_{{\tau }_1}^{\beta }{\theta }_{\varphi }^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)}_{{\mathrm{\Omega }}_p}+{\left(\mu \nabla {\theta }_{\mathbf{u}}^{n-1/2},\nabla {\theta }_{\mathbf{u}}^{n-1/2}\right)}_{{\mathrm{\Omega }}_f}\hfill \\ & +g{\left({\mathbb{K}}_{\omega }\nabla {\theta }_{\varphi }^{n-1/2},\nabla {\theta }_{\varphi }^{n-1/2}\right)}_{{\mathrm{\Omega }}_p}+\mu \sum _{i=1}^{d-1}{\int }_I{\gamma }_{i,\omega }{\left({\theta }_{\mathbf{u}}^{n-1/2}·{\mathbf{m}}_i\right)}^{2}ds\hfill \\ & =-{\left({\overline{\Delta }}_{{\tau }_1}^{\alpha }{\xi }_{\mathbf{u}}^{n-1/2},{\theta }_{\mathbf{u}}^{n-1/2}\right)}_{{\mathrm{\Omega }}_f}-g{\left(S{\overline{\Delta }}_{{\tau }_1}^{\beta }{\xi }_{\varphi }^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)}_{{\mathrm{\Omega }}_p}\hfill \\ & +{\mathcal{H}}_I\left({\xi }_{\mathbf{u}}^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)-{\mathcal{H}}_I\left({\theta }_{\mathbf{u}}^{n-1/2},{\xi }_{\varphi }^{n-1/2}\right)\hfill \\ & +{\mathcal{H}}_I\left({\theta }_{\mathbf{u}}^{n-1/2},r_2^{n-1/2}\right)+{\mathcal{H}}_I\left({\mathbf{r}}_4^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)\hfill \\ & -\mu \sum _{i=1}^{d-1}{\int }_I{\gamma }_{i,\omega }\left({\xi }_{\mathbf{u}}^{n-1/2}·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^{n-1/2}·{\mathbf{m}}_i\right)ds\hfill \\ & +{\left({\mathbf{r}}_{11}^{n-1/2},{\theta }_{\mathbf{u}}^{n-1/2}\right)}_{{\mathrm{\Omega }}_f}+{\left(r_{12}^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)}_{{\mathrm{\Omega }}_p}+g{\left({\mathbb{K}}_{\omega }\nabla r_5^{n-1/2},\nabla {\theta }_{\varphi }^{n-1/2}\right)}_{{\mathrm{\Omega }}_p}\hfill \\ & +\mu \sum _{i=1}^{d-1}{\int }_I{\gamma }_{i,\omega }\left({\mathbf{r}}_3^{n-1/2}·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^{n-1/2}·{\mathbf{m}}_i\right)ds.\hfill \end{array}$$

(59)

Estimate the terms in (59) using the technique applied in Lemma 5. For the first and second terms on the left-hand side of (59), we have:

$${\left({\overline{\Delta }}_{{\tau }_1}^{\alpha }{\theta }_{\mathbf{u}}^{1/2},{\theta }_{\mathbf{u}}^{1/2}\right)}_{{\mathrm{\Omega }}_f}+g{\left(S{\overline{\Delta }}_{{\tau }_1}^{\beta }{\theta }_{\varphi }^{1/2},{\theta }_{\varphi }^{1/2}\right)}_{{\mathrm{\Omega }}_p}=\frac{1}{4}\left(\frac{{\tau }_1^{-\alpha }{\overline{\delta }}_{1,1}^{\left(\alpha \right)}}{{\overline{\varsigma }}_0^{\left(\alpha \right)}}{∥{\theta }_{\mathbf{u}}^1∥}_{{\mathrm{\Omega }}_p}^2+\frac{{\tau }_1^{-\beta }{\overline{\delta }}_{1,1}^{\left(\beta \right)}}{{\overline{\varsigma }}_0^{\left(\beta \right)}}gS{∥{\theta }_{\varphi }^1∥}_{{\mathrm{\Omega }}_p}^2\right),$$

$${\left({\overline{\Delta }}_{{\tau }_1}^{\alpha }{\theta }_{\mathbf{u}}^{3/2},{\theta }_{\mathbf{u}}^{3/2}\right)}_{{\mathrm{\Omega }}_f}+g{\left(S{\overline{\Delta }}_{{\tau }_1}^{\beta }{\theta }_{\varphi }^{3/2},{\theta }_{\varphi }^{3/2}\right)}_{{\mathrm{\Omega }}_p}$$

$$\ge \frac{{\tau }_1^{-\alpha }}{8{\overline{\varsigma }}_0^{\left(\alpha \right)}}\left({\overline{\delta }}_{2,2}^{\left(\alpha \right)}{∥{\theta }_{\mathbf{u}}^2∥}_{{\mathrm{\Omega }}_f}^2-C_{\alpha }{∥{\theta }_{\mathbf{u}}^1∥}_{{\mathrm{\Omega }}_f}^2\right)+\frac{gS{\tau }_1^{-\beta }}{8{\overline{\varsigma }}_0^{\left(\beta \right)}}\left({\overline{\delta }}_{2,2}^{\left(\beta \right)}{∥{\theta }_{\varphi }^2∥}_{{\mathrm{\Omega }}_p}^2-C_{\beta }{∥{\theta }_{\varphi }^1∥}_{{\mathrm{\Omega }}_p}^2\right).$$

For the third, fourth, and fifth terms, we have:

$${\left(\mu \nabla {\theta }_{\mathbf{u}}^{n-1/2},\nabla {\theta }_{\mathbf{u}}^{n-1/2}\right)}_{{\mathrm{\Omega }}_f}+g{\left({\mathbb{K}}_{\omega }\nabla {\theta }_{\varphi }^{n-1/2},\nabla {\theta }_{\varphi }^{n-1/2}\right)}_{{\mathrm{\Omega }}_p}$$

$$\ge \mu {∥\nabla {\theta }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2+g{k}_{min}{∥\nabla {\theta }_{\varphi }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2,$$

$$-\mu \sum _{i=1}^{d-1}{\int }_I{\gamma }_{i,\omega }{\left({\theta }_{\mathbf{u}}^{n-1/2}·{\mathbf{m}}_i\right)}^{2}ds\le \mu {\overline{\gamma }}_i^{max}{C}_{\mathrm{tr}}{C}_{P,f}{∥\nabla {\theta }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2.$$

Considering that ${\xi }_{\mathbf{u}}^0=\mathbf{0}$ and ${\xi }_{\varphi }^0=0$, the first term on the right-hand side of (59) is estimated as follows:

$$-{\left({\overline{\Delta }}_{{\tau }_1}^{\alpha }{\xi }_{\mathbf{u}}^{1/2},{\theta }_{\mathbf{u}}^{1/2}\right)}_{{\mathrm{\Omega }}_f}-g{\left(S{\overline{\Delta }}_{{\tau }_1}^{\beta }{\xi }_{\varphi }^{1/2},{\theta }_{\varphi }^{1/2}\right)}_{{\mathrm{\Omega }}_p}$$

$$\le \frac{ϵ{\tau }_1^{-\alpha }{\overline{\delta }}_{1,1}^{\alpha }}{4{\overline{\varsigma }}_0^{\left(\alpha \right)}}{∥{\theta }_{\mathbf{u}}^1∥}_{{\mathrm{\Omega }}_f}^2+\frac{ϵ{\tau }_1^{-\beta }gS{\overline{\delta }}_{1,1}^{\beta }}{4{\overline{\varsigma }}_0^{\left(\beta \right)}}{∥{\theta }_{\varphi }^1∥}_{{\mathrm{\Omega }}_p}^2+C_{\alpha }{\tau }_1^{-\alpha }{∥{\xi }_{\mathbf{u}}^1∥}_{{\mathrm{\Omega }}_f}^2+C_{\beta }{\tau }_1^{-\beta }{∥{\xi }_{\varphi }^1∥}_{{\mathrm{\Omega }}_p}^2,$$

$$-{\left({\Delta }_{\tau }^{\alpha }{\xi }_{\mathbf{u}}^{3/2},{\theta }_{\mathbf{u}}^{3/2}\right)}_{{\mathrm{\Omega }}_f}-g{\left(S{\Delta }_{\tau }^{\beta }{\xi }_{\varphi }^{3/2},{\theta }_{\varphi }^{3/2}\right)}_{{\mathrm{\Omega }}_p}$$

$$\le \frac{ϵ{\overline{\delta }}_{2,2}^{\alpha }{\tau }_1^{-\alpha }}{8{\overline{\varsigma }}_0^{\left(\alpha \right)}}{∥{\theta }_{\mathbf{u}}^2∥}_{{\mathrm{\Omega }}_f}^2+\frac{ϵ{\overline{\delta }}_{2,2}^{\alpha }{\tau }_1^{-\alpha }}{8{\overline{\varsigma }}_0^{\left(\alpha \right)}}{∥{\theta }_{\mathbf{u}}^1∥}_{{\mathrm{\Omega }}_f}^2+C_{\alpha }{\tau }_1^{-\alpha }\sum _{s=1}^2{∥{\xi }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2$$

$$+\frac{gSϵ{\overline{\delta }}_{2,2}^{\beta }{\tau }_1^{-\beta }}{8{\overline{\varsigma }}_0^{\left(\beta \right)}}{∥{\theta }_{\varphi }^2∥}_{{\mathrm{\Omega }}_p}^2+\frac{gSϵ{\overline{\delta }}_{2,2}^{\beta }{\tau }_1^{-\beta }}{8{\overline{\varsigma }}_0^{\left(\beta \right)}}{∥{\theta }_{\varphi }^1∥}_{{\mathrm{\Omega }}_p}^2+C_{\beta }{\tau }_1^{-\beta }\sum _{s=1}^2{∥{\xi }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2.$$

The remaining terms are estimated as follows:

$${\mathcal{H}}_I\left({\xi }_{\mathbf{u}}^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)-{\mathcal{H}}_I\left({\theta }_{\mathbf{u}}^{n-1/2},{\xi }_{\varphi }^{n-1/2}\right)+{\mathcal{H}}_I\left({\theta }_{\mathbf{u}}^{n-1/2},r_2^{n-1/2}\right)+{\mathcal{H}}_I\left({\mathbf{r}}_4^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)$$

$$\le 2ϵ\mu {∥\nabla {\theta }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2+2ϵg{k}_{min}{∥\nabla {\theta }_{\varphi }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2$$

$$+C\left({∥{\xi }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}{∥\nabla {\xi }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}+{∥{\xi }_{\varphi }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}{∥\nabla {\xi }_{\varphi }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}\right)$$

$$+C\left({∥r_2^{n-1/2}∥}_{{\mathrm{\Omega }}_p}{∥\nabla r_2^{n-1/2}∥}_{{\mathrm{\Omega }}_p}+{∥{\mathbf{r}}_4^{n-1/2}∥}_{{\mathrm{\Omega }}_f}{∥\nabla {\mathbf{r}}_4^{n-1/2}∥}_{{\mathrm{\Omega }}_f}\right),$$

$$\mu \sum _{i=1}^{d-1}{\int }_I{\gamma }_{i,\omega }\left({\xi }_{\mathbf{u}}^{n-1/2}·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^{n-1/2}·{\mathbf{m}}_i\right)ds\le C{∥{\xi }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}{∥\nabla {\xi }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}+ϵ\mu {∥\nabla {\theta }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2,$$

$${\left({\mathbf{r}}_{11}^{n-1/2},{\theta }_{\mathbf{u}}^{n-1/2}\right)}_{{\mathrm{\Omega }}_f}\le C{∥{\mathbf{r}}_{11}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2+ϵ\mu {∥\nabla {\theta }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2,$$

$${\left(r_{12}^{n-1/2},{\theta }_{\varphi }^{n-1/2}\right)}_{{\mathrm{\Omega }}_p}\le C{∥r_{12}^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2+ϵg{k}_{min}{∥\nabla {\theta }_{\varphi }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2,$$

$$g{\left({\mathbb{K}}_{\omega }\nabla r_5^{n-1/2},\nabla {\theta }_{\varphi }^{n-1/2}\right)}_{{\mathrm{\Omega }}_p}\le C{∥\nabla r_5^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2+ϵg{k}_{min}{∥\nabla {\theta }_{\varphi }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2,$$

$$\mu \sum _{i=1}^{d-1}{\int }_I{\gamma }_{i,\omega }\left({\mathbf{r}}_3^{n-1/2}·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^{n-1/2}·{\mathbf{m}}_i\right)ds\le C{∥\nabla {\mathbf{r}}_3^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2+ϵ\mu {∥\nabla {\theta }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2.$$

By taking into account the obtained inequalities and approximation properties (28), we obtain after elementary transformations that

$${\tau }_1^{-\alpha }{∥{\theta }_{\mathbf{u}}^1∥}_{{\mathrm{\Omega }}_f}^2+{\tau }_1^{-\beta }{∥{\theta }_{\varphi }^1∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\theta }_{\mathbf{u}}^{1/2}∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^{1/2}∥}_{{\mathrm{\Omega }}_p}^2$$

$$\le C_{\alpha ,\beta }\left({\tau }_1^{-\alpha }{h}^{2l+2}+{\tau }_1^{-\beta }{h}^{2m+2}+h^{2l+1}+h^{2m+1}+{\tau }_1^{4-2\alpha }+{\tau }_1^{4-2\beta }\right)$$

(60)

and

$${\tau }_1^{-\alpha }{∥{\theta }_{\mathbf{u}}^2∥}_{{\mathrm{\Omega }}_f}^2+{\tau }_1^{-\beta }{∥{\theta }_{\varphi }^2∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\theta }_{\mathbf{u}}^{3/2}∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^{3/2}∥}_{{\mathrm{\Omega }}_p}^2$$

$$\le C_{\alpha ,\beta }\left({\tau }_1^{-\alpha }{∥{\theta }_{\mathbf{u}}^1∥}_{{\mathrm{\Omega }}_f}^2+{\tau }_1^{-\beta }{∥{\theta }_{\varphi }^1∥}_{{\mathrm{\Omega }}_p}^2+h^{2m+1}+h^{2l+1}+{\tau }_1^{-\alpha }{h}^{2l+2}\right)$$

$$+C_{\alpha ,\beta }\left({\tau }_1^{-\beta }{h}^{2m+2}+{\tau }_1^{4-2\alpha }+{\tau }_1^{4-2\beta }\right).$$

(61)

Combining the estimates (60) and (61) yields

$${∥{\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+{∥{\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}+{\tau }_1^{{\nu }_1}\left({∥\nabla {\theta }_{\mathbf{u}}^{n-1/2}∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^{n-1/2}∥}_{{\mathrm{\Omega }}_p}^2\right)$$

$$\le C_{\alpha ,\beta }\left({\tau }_1^{{\nu }_1}\left(h^{2l+1}+h^{2m+1}\right)+h^{2l+2}+h^{2m+2}+{\tau }_1^{4-{\nu }_1}\right),$$

where ${\nu }_1=max\left\{\alpha ,\beta \right\}$. Finally, taking into account the relation (58), we obtain the assertion of the lemma.    □

We now prove the main result of this section.

Theorem 2.

Suppose that $0<\alpha \le \beta <1$ and denote $\nu =\beta -\alpha$. Let $\left({\mathbf{u}}_{\omega ,h}^n,p_{\omega ,h}^n,{\varphi }_{\omega ,h}^n\right)\in X_{fh}\times Q_{fh}\times Q_{ph}$ be the solution of Problem 4, and $\left({\mathbf{u}}_{\omega }^n,p_{\omega }^n,{\varphi }_{\omega }^n\right)\in X_f\times Q_f\times Q_p$ be the solution of Problem 3, $n\ge 1$. Then, the following estimate is valid under Assumptions 1 and 2:

$${\tau }^{\nu }{∥{\mathbf{u}}_{\omega }^n-{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+{∥{\varphi }_{\omega }^n-{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2+{\tau }^{\beta }{∥\nabla \left({\mathbf{u}}_{\omega }^{n-1/2}-{\mathbf{u}}_{\omega ,h}^{n-1/2}\right)∥}_{{\mathrm{\Omega }}_f}^2$$

$$+{\tau }^{\beta }{∥\nabla \left({\varphi }_{\omega }^{n-1/2}-{\varphi }_{\omega ,h}^{n-1/2}\right)∥}_{{\mathrm{\Omega }}_p}^2\le C_{\alpha ,\beta }\left({\tau }^{\frac{2\beta \left(5-2\beta \right)}{4-\beta }}\left(h^{2m+1}+h^{2l+1}\right)+{\tau }^{\nu -1}{h}^{2l+2}\right)$$

$$+C_{\alpha ,\beta }\left({\tau }^{-1}{h}^{2m+2}+{\tau }^{\beta }{h}^{2l}+{\tau }^{\beta }{h}^{2m}+{\tau }^{6-2\alpha +\beta }+{\tau }^{6-\beta }\right),$$

(62)

which implies the convergence of the approximate solution to the solution of the differential problem, where l and m are defined in (28). In the case of $0<\beta \le \alpha <1$ and ${\nu }_1=\alpha -\beta$:

$${∥{\mathbf{u}}_{\omega }^n-{\mathbf{u}}_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{{\nu }_1}{∥{\varphi }_{\omega }^n-{\varphi }_{\omega ,h}^n∥}_{{\mathrm{\Omega }}_p}^2+{\tau }^{\alpha }{∥\nabla \left({\mathbf{u}}_{\omega }^{n-1/2}-{\mathbf{u}}_{\omega ,h}^{n-1/2}\right)∥}_{{\mathrm{\Omega }}_f}^2$$

$$+{\tau }^{\alpha }{∥\nabla \left({\varphi }_{\omega }^{n-1/2}-{\varphi }_{\omega ,h}^{n-1/2}\right)∥}_{{\mathrm{\Omega }}_p}^2\le C_{\alpha ,\beta }\left({\tau }^{\frac{2\alpha \left(5-2\alpha \right)}{4-\alpha }}\left(h^{2m+1}+h^{2l+1}\right)+{\tau }^{{\nu }_1-1}{h}^{2l+2}\right)$$

$$+C_{\alpha ,\beta }\left({\tau }^{-1}{h}^{2m+2}+{\tau }^{\alpha }{h}^{2l}+{\tau }^{\alpha }{h}^{2m}+{\tau }^{6+\alpha -2\beta }+{\tau }^{6-\alpha }\right).$$

(63)

Proof. 

Subtract (30) from (27), use the notation (57), then take $\left({\mathbf{v}}_h,q_h,{\psi }_h\right)=\left({\theta }_{\mathbf{u}}^n,{\theta }_p^n,{\theta }_{\varphi }^n\right)$ to obtain

$$\begin{array}{cc}& {\left({\Delta }_{\tau }^{\alpha }{\theta }_{\mathbf{u}}^n,{\theta }_{\mathbf{u}}^n\right)}_{{\mathrm{\Omega }}_f}+g{\left(S{\Delta }_{\tau }^{\beta }{\theta }_{\varphi }^n,{\theta }_{\varphi }^n\right)}_{{\mathrm{\Omega }}_p}+{\left(\mu \nabla {\theta }_{\mathbf{u}}^n,\nabla {\theta }_{\mathbf{u}}^n\right)}_{{\mathrm{\Omega }}_f}+g{\left(\overline{\mathbb{K}}\nabla {\theta }_{\varphi }^n,\nabla {\theta }_{\varphi }^n\right)}_{{\mathrm{\Omega }}_p}\hfill \\ & +\mu \sum _{i=1}^{d-1}{\int }_I{\overline{\gamma }}_i\left({\theta }_{\mathbf{u}}^n·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^n·{\mathbf{m}}_i\right)ds\hfill \\ & +{\mathcal{H}}_I\left({\theta }_{\mathbf{u}}^n,3{\theta }_{\varphi }^{n-1}-3{\theta }_{\varphi }^{n-2}+{\theta }_{\varphi }^{n-3}\right)-{\mathcal{H}}_I\left(3{\theta }_{\mathbf{u}}^{n-1}-3{\theta }_{\mathbf{u}}^{n-2}+{\theta }_{\mathbf{u}}^{n-3},{\theta }_{\varphi }^n\right)\hfill \\ & +\mu \sum _{i=1}^{d-1}{\int }_I\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right)\left(\left(3{\theta }_{\mathbf{u}}^{n-1}-3{\theta }_{\mathbf{u}}^{n-2}+{\theta }_{\mathbf{u}}^{n-3}\right)·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^n·{\mathbf{m}}_i\right)ds\hfill \\ & +g{\left(\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}}\right)\nabla \left(3{\theta }_{\varphi }^{n-1}-3{\theta }_{\varphi }^{n-2}+{\theta }_{\varphi }^{n-3}\right),\nabla {\theta }_{\varphi }^n\right)}_{{\mathrm{\Omega }}_p}\hfill \\ & =-{\left({\Delta }_{\tau }^{\alpha }{\xi }_{\mathbf{u}}^n,{\theta }_{\mathbf{u}}^n\right)}_{{\mathrm{\Omega }}_f}-g{\left(S{\Delta }_{\tau }^{\beta }{\xi }_{\varphi }^n,{\theta }_{\varphi }^n\right)}_{{\mathrm{\Omega }}_p}\hfill \\ & +{\mathcal{H}}_I\left(3{\xi }_{\mathbf{u}}^{n-1}-3{\xi }_{\mathbf{u}}^{n-2}+{\xi }_{\mathbf{u}}^{n-3},{\theta }_{\varphi }^n\right)-{\mathcal{H}}_I\left({\theta }_{\mathbf{u}}^n,3{\xi }_{\varphi }^{n-1}-3{\xi }_{\varphi }^{n-2}+{\xi }_{\varphi }^{n-3}\right)\hfill \\ & +\mu \sum _{i=1}^{d-1}{\int }_I\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right)\left(\left(3{\xi }_{\mathbf{u}}^{n-1}-3{\xi }_{\mathbf{u}}^{n-2}+{\xi }_{\mathbf{u}}^{n-3}\right)·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^n·{\mathbf{m}}_i\right)ds\hfill \\ & -g{\left(\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}}\right)\nabla \left(3{\xi }_{\varphi }^{n-1}-3{\xi }_{\varphi }^{n-2}+{\xi }_{\varphi }^{n-3}\right),\nabla {\theta }_{\varphi }^n\right)}_{{\mathrm{\Omega }}_p}\hfill \\ & -\mu \sum _{i=1}^{d-1}{\int }_I{\overline{\gamma }}_i\left({\xi }_{\mathbf{u}}^n·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^n·{\mathbf{m}}_i\right)ds\hfill \\ & +{\left(r_{12}^n,{\theta }_{\varphi }^n\right)}_{{\mathrm{\Omega }}_p}+{\mathcal{H}}_I\left({\mathbf{r}}_4^n,{\theta }_{\varphi }^n\right)+{\mathcal{H}}_I\left({\theta }_{\mathbf{u}}^n,r_2^n\right)+g{\left(\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}}\right)\nabla r_5^n,\nabla {\theta }_{\varphi }^n\right)}_{{\mathrm{\Omega }}_p}\hfill \\ & +{\left({\mathbf{r}}_{11}^n,{\theta }_{\mathbf{u}}^n\right)}_{{\mathrm{\Omega }}_f}+\mu \sum _{i=1}^{d-1}{\int }_I\left({\gamma }_{i,\omega }-{\overline{\gamma }}_i\right)\left({\mathbf{r}}_3^n·{\mathbf{m}}_i\right)\left({\theta }_{\mathbf{u}}^n·{\mathbf{m}}_i\right)ds.\hfill \end{array}$$

(64)

Let us estimate the terms in (64). By using the technique applied in Lemma 7, one obtains

$${\tau }^{-\alpha }\left(1-ϵ\right){∥{\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }\left(1-ϵ\right){∥{\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2$$

$$+\mu \left(1-{\overline{\eta }}^{max}{C}_{\mathrm{tr}}{C}_{P,f}-7ϵ\right){∥\nabla {\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+g{k}_{min}\left(1-7ϵ\right){∥\nabla {\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2$$

$$\le {\tau }^{-\alpha }{\chi }_{\alpha }\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\alpha \right)}\right|{∥{\theta }_{\mathbf{u}}^{n-s}∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{\chi }_{\beta }\sum _{s=1}^n\left|{\delta }_{n,n-s}^{\left(\beta \right)}\right|{∥{\theta }_{\varphi }^{n-s}∥}_{{\mathrm{\Omega }}_p}^2$$

$$+C\sum _{s=n-3}^{n-1}\left({∥\nabla {\theta }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2\right)+C_{\alpha }{\tau }^{-\alpha }\sum _{s=0}^n{∥{\xi }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2+C_{\beta }{\tau }^{-\beta }\sum _{s=0}^n{∥{\xi }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2$$

$$+C\sum _{s=n-3}^{n-1}{∥\nabla {\xi }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2+C\sum _{s=n-3}^{n-1}{∥\nabla {\xi }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2+C{∥{\xi }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}{∥\nabla {\xi }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}$$

$$+C\left({∥{\mathbf{r}}_{11}^n∥}_{{\mathrm{\Omega }}_f}^2+{∥r_{12}^n∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla r_2^n∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\mathbf{r}}_3^n∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\mathbf{r}}_4^n∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla r_5^n∥}_{{\mathrm{\Omega }}_p}^2\right).$$

Applying Assumption (40) and using the technique similarly to obtaining inequality (53), we arrive at the inequality

$${\tau }^{-\alpha }{∥{\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2$$

$$\le C_{\alpha ,\beta }\sum _{s=0}^2\left|{\delta }_{n,s}^{\left(\alpha \right)}\right|\left({\tau }^{-\alpha }{∥{\theta }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\theta }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2\right)$$

$$+C_{\alpha ,\beta }^{\prime }\sum _{s=3}^{n-1}\left|{\delta }_{n,s}^{\left(\alpha \right)}\right|\left({\tau }^{-\alpha }{∥{\theta }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\theta }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\theta }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2\right)$$

$$+C_{\alpha ,\beta }\left({\tau }^{-\alpha -1}{h}^{2l+2}+{\tau }^{-\beta -1}{h}^{2m+2}+h^{2l}+h^{2m}+{\tau }^{6-2\alpha }+{\tau }^{6-2\beta }\right).$$

By utilizing the result of Lemma 7 and taking into account the fact that $\alpha \le \beta$, we obtain

$${\tau }^{-\alpha }{∥{\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2$$

$$\le C_{\alpha ,\beta }\sum _{s=3}^{n-1}\left|{\delta }_{n,s}^{\left(\alpha \right)}\right|\left({\tau }^{-\alpha }{∥{\theta }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\theta }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\theta }_{\mathbf{u}}^s∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^s∥}_{{\mathrm{\Omega }}_p}^2\right)$$

$$+C_{\alpha ,\beta }\left({\tau }^{\frac{6\beta -{\beta }^2}{4-\beta }}\left(h^{2m+1}+h^{2l+1}\right)+{\tau }^{-\alpha -1}{h}^{2l+2}+{\tau }^{-\beta -1}{h}^{2m+2}+h^{2l}+h^{2m}+{\tau }^{6-2\alpha }+{\tau }^{6-2\beta }\right).$$

Finally, applying Lemma 4, we obtain

$${\tau }^{-\alpha }{∥{\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+{\tau }^{-\beta }{∥{\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2+{∥\nabla {\theta }_{\mathbf{u}}^n∥}_{{\mathrm{\Omega }}_f}^2+{∥\nabla {\theta }_{\varphi }^n∥}_{{\mathrm{\Omega }}_p}^2$$

$$\le C_{\alpha ,\beta }\left({\tau }^{\frac{6\beta -{\beta }^2}{4-\beta }}\left(h^{2m+1}+h^{2l+1}\right)+{\tau }^{-\alpha -1}{h}^{2l+2}+{\tau }^{-\beta -1}{h}^{2m+2}+h^{2l}+h^{2m}+{\tau }^{6-2\alpha }+{\tau }^{6-2\beta }\right),$$

which implies the estimate (62). Inequality (63) is proved in a similar way.    □

2.4. Implementation of the Discrete Scheme

Assume that the permeability tensor ${\mathbb{K}}_{\omega }$, fractional derivative orders $\alpha ,\beta$, initial data ${\mathbf{u}}_{0,\omega },{\varphi }_{0,\omega }$, functions $S,{\mathbf{f}}_{f,\omega },f_{p,\omega }$, constants $\mu ,{\alpha }_{\mathrm{BJS}}$, and discretization parameter $\tau$ are known. Let us present the following Algorithm 1 for solving the $\omega$th realization of the discrete scheme (29) and (30).

Algorithm 1. Implementation of the $\omega$th realization.

1. For the first $\left(n=1\right)$ and second $\left(n=2\right)$ time layers:

1.1. Introduce subtime layers $\left\{t_{n,k}=n\tau +k{\tau }_1,k=0,1,\dots ,L,n=0,1\right\}$ within the time intervals $\left[0,t_1\right]$ and $\left[t_1,t_2\right]$ according to (58).

1.2. Calculate the coefficients ${\overline{\delta }}_{n,s}^{\left(\nu \right)}$, $\nu \in \left\{\alpha ,\beta \right\}$, $s=1,2,\dots ,n$ for the current n.

1.3. Find velocity ${\mathbf{u}}_{\omega ,h}^n\in X_f$ and pressure $p_{\omega ,h}^n\in Q_f$ in ${\mathrm{\Omega }}_f$ according to (29) with ${\widehat{\mathbf{v}}}_h={\left[{\mathbf{v}}_h,0\right]}^{\top }$, or equivalently, $$\begin{array}{c}\frac{{\tau }^{-\alpha }{\overline{\delta }}_{n,n}^{\alpha }}{{\overline{\varsigma }}_0^{\left(\alpha \right)}}{\left({\mathbf{u}}_{\omega ,h}^n,{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}+{\left(\mu \nabla {\mathbf{u}}_{\omega ,h}^n,\nabla {\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}+c_I\left({\mathbf{v}}_h,{\varphi }_{\omega ,h}^n\right)\hfill \\ -\mu \sum _{i=1}^{d-1}{\int }_I{\eta }_{i,\omega }\left({\mathbf{u}}_{\omega ,h}^n·{\mathbf{m}}_i\right)\left({\mathbf{v}}_h·{\mathbf{m}}_i\right)ds-{\left(p_{\omega ,h}^n,\nabla ·{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}\hfill \\ -\frac{{\tau }^{-\alpha }{\overline{\delta }}_{n,n}^{\alpha }}{{\overline{\varsigma }}_0^{\left(\alpha \right)}}{\left({\mathbf{u}}_{\omega ,h}^{n-1},{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}+\frac{{\tau }^{-\alpha }}{{\overline{\varsigma }}_0^{\left(\alpha \right)}}\sum _{s=1}^{n-1}\left({\overline{\delta }}_{n,s}^{\alpha }+{\overline{\delta }}_{n,s+1}^{\alpha }\right){\left({\mathbf{u}}_{\omega ,h}^s-{\mathbf{u}}_{\omega ,h}^{s-1},{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}\hfill \\ +{\left(\mu \nabla {\mathbf{u}}_{\omega ,h}^{n-1},\nabla {\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}+c_I\left({\mathbf{v}}_h,{\varphi }_{\omega ,h}^{n-1}\right)-{\left(p_{\omega ,h}^{n-1},\nabla ·{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}\hfill \\ -\mu \sum _{i=1}^{d-1}{\int }_I{\eta }_{i,\omega }\left({\mathbf{u}}_{\omega ,h}^{n-1}·{\mathbf{m}}_i\right)\left({\mathbf{v}}_h·{\mathbf{m}}_i\right)ds={\left(2{\mathbf{f}}_{f,\omega }^{n-1/2},{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f},\hfill \\ {\left(q_h,\nabla ·{\mathbf{u}}_{\omega ,h}^n\right)}_{{\mathrm{\Omega }}_f}+{\left(q_h,\nabla ·{\mathbf{u}}_{\omega ,h}^{n-1}\right)}_{{\mathrm{\Omega }}_f}=0\hfill \end{array}$$ (65) for all $\left({\mathbf{v}}_h,q_h\right)\in X_f\times Q_f$.

1.4. Find piezometric head ${\varphi }_{\omega ,h}^n\in Q_p$ in ${\mathrm{\Omega }}_p$ according to (29) with ${\widehat{\mathbf{v}}}_h={\left[\mathbf{0},{\psi }_h\right]}^{\top }$ and $q_h=0$, or equivalently, $$\begin{array}{}\\ \frac{g{\tau }^{-\beta }{\overline{\delta }}_{n,n}^{\left(\beta \right)}}{{\overline{\varsigma }}_0^{\left(\beta \right)}}{\left(S{\varphi }_{\omega ,h}^n,{\psi }_h\right)}_{{\mathrm{\Omega }}_p}+g{\left({\mathbb{K}}_{\omega }\nabla {\varphi }_{\omega ,h}^n,\nabla {\psi }_h\right)}_{{\mathrm{\Omega }}_p}-c_I\left({\mathbf{u}}_{\omega ,h}^n,{\psi }_h\right)\hfill \\ +g{\left({\mathbb{K}}_{\omega }\nabla {\varphi }_{\omega ,h}^{n-1},\nabla {\psi }_h\right)}_{{\mathrm{\Omega }}_p}-\frac{g{\tau }^{-\beta }{\overline{\delta }}_{n,n}^{\left(\beta \right)}}{{\overline{\varsigma }}_0^{\left(\beta \right)}}{\left(S{\varphi }_{\omega ,h}^{n-1},{\psi }_h\right)}_{{\mathrm{\Omega }}_p}-c_I\left({\mathbf{u}}_{\omega ,h}^{n-1},{\psi }_h\right)\hfill \\ +\frac{g{\tau }^{-\beta }}{{\overline{\varsigma }}_0^{\left(\beta \right)}}\sum _{s=1}^{n-1}\left({\overline{\delta }}_{n,s}^{\left(\beta \right)}+{\overline{\delta }}_{n,s+1}^{\left(\beta \right)}\right){\left({\varphi }_{\omega ,h}^s-{\varphi }_{\omega ,h}^{s-1},S{\psi }_h\right)}_{{\mathrm{\Omega }}_p}=2g{\left(f_{p,\omega }^{n-1/2},{\psi }_h\right)}_{{\mathrm{\Omega }}_p}\hfill \end{array}$$ (66) for all ${\psi }_h\in Q_p$.

2. For the other time layers $\left(n=3,4,\dots ,N\right)$:

2.1. Calculate the coefficients ${\delta }_{n,s}^{\left(\nu \right)}$, $s=0,1,\dots ,n$, $\nu \in \left\{\alpha ,\beta \right\}$, the sums $$\sum _{s=1}^n{\delta }_{n,n-s}^{\left(\alpha \right)}{\left({\mathbf{u}}_{\omega ,h}^{n-s},{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}\mathrm{and}\sum _{s=1}^n{\delta }_{n,n-s}^{\left(\beta \right)}{\left(S{\varphi }_{\omega ,h}^{n-s},{\psi }_h\right)}_{{\mathrm{\Omega }}_p}$$ (67) entering (68) and (69) below, as well as finite element functions $\sigma {\mathbf{u}}_{\omega ,h}^{n-1}$, $\sigma {\varphi }_{\omega ,h}^{n-1}$.

2.2. Find ${\mathbf{u}}_{\omega ,h}^n\in X_f$ and $p_{\omega ,h}^n\in Q_f$ in ${\mathrm{\Omega }}_f$ according to (30) with ${\widehat{\mathbf{v}}}_h={\left[{\mathbf{v}}_h,0\right]}^{\top }$, or equivalently, $$\begin{array}{c}\frac{{\lambda }_0^{\left(\alpha \right)}{\tau }^{-\alpha }}{{\varsigma }_0^{\left(\alpha \right)}}{\left({\mathbf{u}}_{\omega ,h}^n,{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}+{\left(\mu \nabla {\mathbf{u}}_{\omega ,h}^n,\nabla {\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}-\mu \sum _{i=1}^{d-1}{\int }_{\mathrm{\Gamma }}{\overline{\eta }}_i\left({\mathbf{u}}_{\omega ,h}^n·{\mathbf{m}}_i\right)\left({\mathbf{v}}_h·{\mathbf{m}}_i\right)ds\hfill \\ +c_I\left({\mathbf{v}}_h,\sigma {\varphi }_{\omega ,h}^{n-1}\right)-\frac{{\lambda }_0^{\left(\alpha \right)}{\tau }^{-\alpha }}{{\varsigma }_0^{\left(\alpha \right)}}\sum _{s=1}^n{\delta }_{n,n-s}^{\left(\alpha \right)}{\left({\mathbf{u}}_{\omega ,h}^{n-s},{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}-{\left(p_{\omega ,h}^n,\nabla ·{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f}\hfill \\ -\mu \sum _{i=1}^{d-1}{\int }_{\mathrm{\Gamma }}\left({\eta }_{i,\omega }-{\overline{\eta }}_i\right)\left(\sigma {\mathbf{u}}_{\omega ,h}^{n-1}·{\mathbf{m}}_i\right)\left({\mathbf{v}}_h·{\mathbf{m}}_i\right)ds={\left({\mathbf{f}}_{f,\omega }^n,{\mathbf{v}}_h\right)}_{{\mathrm{\Omega }}_f},\hfill \\ {\left(q_h,\nabla ·{\mathbf{u}}_{\omega ,h}^n\right)}_{{\mathrm{\Omega }}_f}=0\hfill \end{array}$$ (68) for all $\left({\mathbf{v}}_h,q_h\right)\in X_f\times Q_f$.

2.3. Find ${\varphi }_{\omega ,h}^n\in Q_p$ in ${\mathrm{\Omega }}_p$ according to (30) with ${\widehat{\mathbf{v}}}_h={\left[\mathbf{0},{\psi }_h\right]}^{\top }$ and $q_h=0$, or equivalently, $$\begin{array}{c}g\frac{{\lambda }_0^{\left(\beta \right)}{\tau }^{-\beta }}{{\varsigma }_0^{\left(\beta \right)}}{\left(S{\varphi }_{\omega ,h}^n,{\psi }_h\right)}_{{\mathrm{\Omega }}_p}+g{\left(\overline{\mathbb{K}}\nabla {\varphi }_{\omega ,h}^n,\nabla {\psi }_h\right)}_{{\mathrm{\Omega }}_p}\hfill \\ -g\frac{{\lambda }_0^{\left(\beta \right)}{\tau }^{-\beta }}{{\varsigma }_0^{\left(\beta \right)}}\sum _{s=1}^n{\delta }_{n,n-s}^{\left(\beta \right)}{\left(S{\varphi }_{\omega ,h}^{n-s},{\psi }_h\right)}_{{\mathrm{\Omega }}_p}-c_I\left(\sigma {\mathbf{u}}_{\omega ,h}^{n-1},{\psi }_h\right)\hfill \\ +g{\left(\left({\mathbb{K}}_{\omega }-\overline{\mathbb{K}}\right)\nabla \left(\sigma {\varphi }_{\omega ,h}^{n-1}\right),\nabla {\psi }_h\right)}_{{\mathrm{\Omega }}_p}=g{\left(f_{p,\omega }^n,{\psi }_h\right)}_{{\mathrm{\Omega }}_p}\hfill \end{array}$$ (69) for all ${\psi }_h\in Q_p$.

Further implementation of the finite element method is well known (see, for example, [44]).

3. Results

3.1. Verification of the Convergence Order

In this section, we present the results of numerical experiments to illustrate their agreement with theoretical studies. The aim of the first numerical tests is to determine the relation between the empirical convergence order and the orders of fractional derivatives, as well as to compare the empirical convergence order with the theoretical convergence order obtained in Theorem 2. Consider the domains ${\mathrm{\Omega }}_f=\left\{\left(x_1,x_2\right)\in {\mathbb{R}}^2:0<x_1<\pi ,-1<x_2<0\right\}$ and ${\mathrm{\Omega }}_p=\left\{\left(x_1,x_2\right)\in {\mathbb{R}}^2:0<x_1<\pi ,0<x_2<1\right\}$. Consider Problem 1 in $\overline{\mathrm{\Omega }}={\overline{\mathrm{\Omega }}}_f\cup {\overline{\mathrm{\Omega }}}_p$ with the input data: ${\alpha }_{\mathrm{BJS}}={10}^{-4}$, $\mu ={10}^{-6}$, $g=S=0.01$, $\mathbb{K}=\mathrm{diag}\left(k_{11},k_{22}\right)$, where the right-hand sides of the equations and the initial conditions are determined from the following exact solution:

$$\begin{array}{cc}\hfill & \mathbf{u}={\left[u_1,u_2\right]}^{\top },\hfill \\ \hfill & u_1=k_{22}{t}^6\left(\left(x_1^2-\pi x_1-2\right)sin{x}_1+\left(2x_1-\pi \right)cos{x}_1+\pi \right)\hfill \\ \hfill & \times \left(\left(-{\pi }^2{\left(x_2^2-1\right)}^2+12x_2^2-4\right)sin\left(\pi x_2\right)+8\pi x_2\left(x_2^2-1\right)cos\left(\pi x_2\right)\right),\hfill \\ \hfill & u_2=-t^6{k}_{22}{x}_1\left(x_1-\pi \right)\left(x_2^2-1\right)cos{x}_1\left(4x_{2}sin\left(\pi x_2\right)+\pi \left(x_2^2-1\right)cos\left(\pi x_2\right)\right),\hfill \\ \hfill & p=t^6{x}_1\left(x_1-\pi \right)cos{x}_1\left(\left(\left(g+\mu k_{22}{\pi }^2\right){\left(x_2^2-1\right)}^2-4\mu k_{22}\left(3x_2^2-1\right)\right)sin\left(\pi x_2\right)\right.\hfill \\ \hfill & \left.-8\pi \mu k_{22}{x}_2\left(x_2^2-1\right)cos\left(\pi x_2\right)\right),\hfill \\ \hfill & \varphi =t^6{x}_1\left(x_1-\pi \right)cos{x}_1{\left(x_2^2-1\right)}^{2}sin\left(\pi x_2\right).\hfill \end{array}$$

Let us describe the procedure for finding the empirical convergence order of the numerical scheme. We consider a set of time discretization parameters ${\tau }_i$, $i=1,2,\dots$ shown in the first column of

Table 1. $L^2$ errors and convergence orders obtained for random tensors ${\mathbb{K}}_{\omega }$ when $\alpha =0.98$ and $\beta =0.6$.

$\mathit{\tau }\approx \mathit{h}$	${\mathit{k}}_{11}={\mathit{k}}_{22}=0.05$		${\mathit{k}}_{11}={\mathit{k}}_{22}=0.01$		${\mathit{k}}_{11}={\mathit{k}}_{22}=0.005$
	${∥{\mathbf{u}}_{\mathbf{1}}-{\mathbf{u}}_{\mathbf{1},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{f}}}$	Order	${∥{\mathbf{u}}_{\mathbf{2}}-{\mathbf{u}}_{\mathbf{2},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{f}}}$	Order	${∥{\mathbf{u}}_{\mathbf{3}}-{\mathbf{u}}_{\mathbf{3},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{f}}}$	Order
$1/5$	$9.4221\times {10}^{-2}$	–	$2.0260\times {10}^{-2}$	–	$1.1025\times {10}^{-2}$	–
$1/10$	$1.8790\times {10}^{-2}$	2.33 $\left(\simeq 2.30\right)$	$4.0021\times {10}^{-3}$	2.34 $\left(\simeq 2.30\right)$	$2.1937\times {10}^{-3}$	2.33 $\left(\simeq 2.30\right)$
$1/20$	$3.7822\times {10}^{-3}$	2.31 $\left(\simeq 2.30\right)$	$7.9993\times {10}^{-4}$	2.32 $\left(\simeq 2.30\right)$	$4.3934\times {10}^{-4}$	2.32 $\left(\simeq 2.30\right)$
$1/40$	$7.6509\times {10}^{-4}$	2.31 $\left(\simeq 2.30\right)$	$1.6205\times {10}^{-4}$	2.30 $\left(\simeq 2.30\right)$	$8.9009\times {10}^{-5}$	2.30 $\left(\simeq 2.30\right)$
$\mathit{\tau }\approx \mathit{h}$	${∥{\mathit{\varphi }}_{\mathbf{1}}-{\mathit{\varphi }}_{\mathbf{1},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{p}}}$	Order	${∥{\mathit{\varphi }}_{\mathbf{2}}-{\mathit{\varphi }}_{\mathbf{2},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{p}}}$	Order	${∥{\mathit{\varphi }}_{\mathbf{3}}-{\mathit{\varphi }}_{\mathbf{3},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{p}}}$	Order
$1/5$	$4.9494\times {10}^{-1}$	–	$3.9874\times {10}^{-1}$	–	$3.4272\times {10}^{-1}$	–
$1/10$	$8.6492\times {10}^{-2}$	2.52 $\left(\simeq 2.49\right)$	$7.0000\times {10}^{-2}$	2.51 $\left(\simeq 2.49\right)$	$5.9839\times {10}^{-2}$	2.52 $\left(\simeq 2.49\right)$
$1/20$	$1.5203\times {10}^{-2}$	2.51 $\left(\simeq 2.49\right)$	$1.2358\times {10}^{-2}$	2.50 $\left(\simeq 2.49\right)$	$1.0553\times {10}^{-2}$	2.50 $\left(\simeq 2.49\right)$
$1/40$	$2.6973\times {10}^{-3}$	2.49 $\left(\simeq 2.49\right)$	$2.1954\times {10}^{-3}$	2.49 $\left(\simeq 2.49\right)$	$1.8773\times {10}^{-3}$	2.49 $\left(\simeq 2.49\right)$

and generate a finite element mesh with a diameter $h_i$ such that $h_i\approx {\tau }_i$. Further, we solve Problem 1 several times with the specified discretization parameters and evaluate the corresponding errors $E_i$. The empirical convergence orders are then calculated using the formula

$${\mathrm{\Theta }}_i=\frac{log\left(E_i/E_{i-1}\right)}{log\left(h_i/h_{i-1}\right)},i=2,3,\dots$$

We employ Taylor–Hood elements in all numerical tests, which are known to satisfy the LBB condition.

The stochastic nature of the problem requires that it be solved multiple times with different input data. However, we confine ourselves to two simple cases in this section. The first case assumes that only the components of the hydraulic conductivity tensor $\mathbb{K}$ are random parameters. To be precise, consider a set of tensors ${\left\{{\mathbb{K}}_{\omega }\right\}}_{\omega =1}^{N_{\omega }}$ with constant diagonal components $k_{11}=k_{22}=k$, then solve Problem 1 $N_{\omega }$ times for each ${\mathbb{K}}_{\omega }$. The results of calculating the errors for $N_{\omega }=3$, as well as the empirical order of convergence, are shown in Table 1.

It is assumed throughout the section that the orders of the fractional derivatives satisfy the relation $\beta \le \alpha$. This is due to the fact that the influence of memory effects is usually insignificant in the domain of free fluid flow; therefore, $\alpha$ should be taken close to 1. A completely different situation is usually observed in a porous medium, since its properties change with time, which can significantly affect the nature of the flow in this medium [45]. Thus, we consider the case $\alpha =0.98$ and $\beta =0.6$ for definiteness in this numerical test. For convenience, the theoretical orders of convergence predicted in Theorem 2 are given in parentheses for this choice of the parameters.

The second case assumes that only orders of fractional derivatives are random parameters. This assumption is due to the fact that the properties of media can be uncertain. More precisely, the degree of influence of memory on the process of flow cannot be determined unambiguously. Therefore, the study considers the simplest case, when the parameter $\alpha =0.98$ is fixed, and the parameter $\beta$ takes on different values ${\beta }_{\omega }$, $\omega =1,2,\dots ,N_{\omega }$. The corresponding calculation results are shown in

Table 2. $L^2$ errors and convergence orders obtained for a fixed $\alpha =0.98$, random parameters ${\beta }_{\omega }$, and $k_{11}=k_{22}=0.01$.

$\mathit{\tau }\approx \mathit{h}$	${\mathit{\beta }}_1=0.55$		${\mathit{\beta }}_2=0.60$		${\mathit{\beta }}_3=0.65$
	${∥{\mathbf{u}}_{\mathbf{1}}-{\mathbf{u}}_{\mathbf{1},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{f}}}$	Order	${∥{\mathbf{u}}_{\mathbf{2}}-{\mathbf{u}}_{\mathbf{2},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{f}}}$	Order	${∥{\mathbf{u}}_{\mathbf{3}}-{\mathbf{u}}_{\mathbf{3},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{f}}}$	Order
$1/5$	$1.8414\times {10}^{-2}$	–	$2.0260\times {10}^{-2}$	–	$2.2289\times {10}^{-2}$	–
$1/10$	$3.6877\times {10}^{-3}$	2.32 $\left(\simeq 2.28\right)$	$4.0021\times {10}^{-3}$	2.34 $\left(\simeq 2.30\right)$	$4.3719\times {10}^{-3}$	2.35 $\left(\simeq 2.33\right)$
$1/20$	$7.4884\times {10}^{-4}$	2.30 $\left(\simeq 2.28\right)$	$7.9993\times {10}^{-4}$	2.32 $\left(\simeq 2.30\right)$	$8.6350\times {10}^{-4}$	2.34 $\left(\simeq 2.33\right)$
$1/40$	$1.5442\times {10}^{-4}$	2.28 $\left(\simeq 2.28\right)$	$1.6205\times {10}^{-4}$	2.30 $\left(\simeq 2.30\right)$	$1.7174\times {10}^{-4}$	2.33 $\left(\simeq 2.33\right)$
$\mathit{\tau }\approx \mathit{h}$	${∥{\mathit{\varphi }}_{\mathbf{1}}-{\mathit{\varphi }}_{\mathbf{1},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{p}}}$	Order	${∥{\mathit{\varphi }}_{\mathbf{2}}-{\mathit{\varphi }}_{\mathbf{2},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{p}}}$	Order	${∥{\mathit{\varphi }}_{\mathbf{3}}-{\mathit{\varphi }}_{\mathbf{3},\mathit{h}}∥}_{{\mathrm{\Omega }}_{\mathit{p}}}$	Order
$1/5$	$4.0568\times {10}^{-1}$	–	$3.9874\times {10}^{-1}$	–	$3.9081\times {10}^{-1}$	–
$1/10$	$7.1219\times {10}^{-2}$	2.51 $\left(\simeq 2.49\right)$	$7.0000\times {10}^{-2}$	2.51 $\left(\simeq 2.49\right)$	$6.8609\times {10}^{-2}$	2.51 $\left(\simeq 2.49\right)$
$1/20$	$1.2590\times {10}^{-2}$	2.50 $\left(\simeq 2.49\right)$	$1.2358\times {10}^{-2}$	2.50 $\left(\simeq 2.49\right)$	$1.2128\times {10}^{-2}$	2.50 $\left(\simeq 2.49\right)$
$1/40$	$2.2411\times {10}^{-3}$	2.49 $\left(\simeq 2.49\right)$	$2.1954\times {10}^{-3}$	2.49 $\left(\simeq 2.49\right)$	$2.1589\times {10}^{-3}$	2.49 $\left(\simeq 2.49\right)$

.

It follows from the results of both computational experiments that the empirical convergence order is in complete agreement with the theoretical order of convergence predicted in Theorem 2.

3.2. Application of the Numerical Method to the Implementation of the Fractional-Order Stochastic Stokes–Darcy Model

In the case of the stochastic Stokes–Darcy model with integer time derivatives, it is known [3] that the assembly of the stiffness matrix is the main bottleneck of the computational process due to the need to solve the problem multiple times with random input parameters. Therefore, the well-known ensemble method significantly reduces the computational complexity of the problem, since the same stiffness matrix is used for all implementations of the ensemble.

As applied to the fractional-order generalization of the stochastic Stokes–Darcy model, this statement is not entirely valid. This is due to the fact that the evaluation of the discrete fractional derivatives (18) and (19) has a greater computational complexity than the assembly of the stiffness matrix. This complexity increases significantly with an increasing time layer number. Let us demonstrate this statement with a simple example.

Consider a set of time discretization parameters ${\tau }_i$ equal to $1/10$, $1/20$, $1/40$, $1/80$, and $1/160$, generate finite element meshes such that ${\tau }_i\approx h_i$, then solve the problem studied in Section 3.1 with random hydraulic conductivity tensors ${\mathbb{K}}_{\omega }$, $\omega =1,2,\dots ,N_{\omega }$. Note that the time required for evaluation of the discrete fractional derivatives by (18) and (19), as well as the construction of finite element functions for the terms $\sigma {\mathbf{u}}_h^n$ and $\sigma {\varphi }_h^n$, took nearly 10.9 milliseconds for $h\approx \tau =1/10$ (

Table 3. Evaluation timings of discrete fractional derivatives in the case of $N_{\omega }=10$.

$\mathit{\tau }\approx \mathit{h}$	Degrees of Freedom in $\left({\mathit{u}}_{\mathit{h}}^{\mathit{n}};{\mathit{\varphi }}_{\mathit{h}}^{\mathit{n}}\right)$	Evaluation of ${\mathit{\delta }}_{\mathit{n},\mathit{s}}^{\left(\mathit{\nu }\right)}$ and ${\overline{\mathit{\delta }}}_{\mathit{n},\mathit{s}}^{\left(\mathit{\nu }\right)}$	Construction of FE Functions $\mathit{\sigma }{\mathbf{u}}_{\mathit{h}}^{\mathit{n}}$ and $\mathit{\sigma }{\mathit{\varphi }}_{\mathit{h}}^{\mathit{n}}$	Construction of FE Functions with Discrete Fractional Derivatives by (18) and (19)	Total
$1/10$	$\left(2524;1262\right)$	48.1 $\mathsf{\mu }$s	3.78 ms	6.98 ms	10.9 ms
$1/20$	$\left(10,004;5002\right)$	218 $\mathsf{\mu }$s	19.5 ms	618 ms	638 ms
$1/40$	$\left(40,164;20,082\right)$	599 $\mathsf{\mu }$s	84.4 ms	7.23 s	7.32 s
$1/80$	$\left(160,964;80,482\right)$	2 ms	2.71 s	111 s	114 s
$1/160$	$\left(643,204;321,602\right)$	7.55 ms	37.1 s	1902 s	1939 s

), which is less than 0.2% of the total CPU time. However, the calculation of the specified data on a finer grid, $h\approx \tau =1/160$, required approximately 76% of the total CPU time.

Therefore, we extend Algorithm 1 such that the calculation of the discrete fractional derivatives is performed simultaneously for all implementations of the ensemble before the implementation of each time layer along with the ensemble algorithm proposed in [3,17]. The extended algorithm is presented in Algorithm 2.

Algorithm 2. Implementation of the stochastic problem.

1. Set the deterministic input data of the problem $\left(\mu ,S,{\alpha }_{\mathrm{BJS}}\right)$, the number of random samples $N_{\omega }$, as well as the mean and variance of random parameters used to generate random physical data.

2. Generate random parameters such as ${\mathbb{K}}_{\omega }$, ${\mathbf{u}}_{0,\omega },{\varphi }_{0,\omega }$, ${\mathbf{f}}_{f,\omega }$, and $f_{p,\omega }$ for $\omega =1,2,\dots ,N_{\omega }$.

3. For each sample number $\omega$, find ${\mathbf{u}}_{\omega ,h}^n$, $p_{\omega ,h}^n$, and ${\varphi }_{\omega ,h}^n$ at the first two time layers according to Steps 1.1–1.4 of Algorithm 1.

4. For the other time layers $\left(n=3,4,\dots ,N\right)$:

4.1. Calculate in parallel the coefficients ${\delta }_{n,s}^{\left(\nu \right)}$, $s=0,1,\dots ,n$, $\nu \in \left\{\alpha ,\beta \right\}$, the sums (67), as well as finite element functions $\sigma {\mathbf{u}}_{\omega ,h}^{n-1}$, $\sigma {\varphi }_{\omega ,h}^{n-1}$ for all $\omega =1,2,\dots ,N_{\omega }$.

4.2. Compute the entries of the stiffness matrix and right-hand side for $\omega =1$ according to (68). Evaluate the right-hand side entries for the rest sample numbers $\omega =2,3,\dots ,N_{\omega }$. Then solve the resulting $N_{\omega }$ systems of equations with identical matrices.

4.3. Compute the entries of the stiffness matrix and right-hand side for $\omega =1$ according to (69). Evaluate the right-hand side entries for the other sample numbers $\omega =2,3,\dots ,N_{\omega }$. Then solve the resulting $N_{\omega }$ systems of equations with identical matrices.

As a result of Steps 4.2 and 4.3, the solution on the nth time layer is obtained for all realizations $\omega =1,2,\dots ,N_{\omega }$.

5. Based on the obtained $N_{\omega }$ solutions of the stochastic problem, evaluate statistical moments.

Note that Steps 4.2 and 4.3 in Algorithm 2 are also executed in parallel, since the corresponding integral identities do not depend on each other.

Now let us analyze the results of some numerical tests conducted using the presented algorithm. The numerical tests were carried out on a workstation with a 24-core Intel(R) Xeon(R) E5-2650 processor with a clock speed of 2.20 GHz and 64 GB of RAM.

We considered $N_{\omega }=1,2,5,10,50,100,500,$ and 1000 samples of random permeability tensors. The discretization step was chosen to be $\tau =0.1$ and triangulation was generated as in Section 3.1. Firstly, the problem was solved individually for all samples without using the ensemble approach. The results of measuring the total time required to solve the problems are given in the column “Serial Algorithm” of

Table 4. Discrete fractional derivative evaluation timing (in seconds) in the case of $\tau =1/10$.

${\mathit{N}}_{\mathit{\omega }}$	Serial Algorithm	Parallel Ensemble Algorithm	Acceleration Rate
1	2.68	–	–
2	6.81	3.71	1.83
5	13.45	5.15	2.61
10	24.52	6.51	3.76
50	120.15	28.04	4.28
100	276.73	53.45	5.18
500	1470.68	223.14	6.59
1000	2823.43	419.69	6.73

. It can be seen that the implementation of the fractional-order stochastic Stokes–Darcy model took around 47 min for 1000 random samples.

Then, we combined the ensemble approach and the technique of parallel evaluation of the discrete fractional derivatives. It is evident from the results presented in the column “Parallel Ensemble Algorithm” of Table 4 that the calculation time reduced to nearly 7 min for the same number of samples. Thus, an almost seven-fold acceleration was achieved.

4. Discussion

Now let us make a few comments on the results obtained and outline future research directions. As far as we know, a fractional-order stochastic generalization of the Stokes–Darcy model is considered for the first time in our paper. Therefore, in the current work, the main emphasis is on the construction of an efficient finite element method for the implementation of this model. Note that the constructed method has a higher convergence order in time than in original papers [3,6] due to the use of a special approximation of fractional derivatives and some of the terms in a weak formulation. This was demonstrated on a test problem with a known exact solution. The use of the approximation formula for the fractional derivative presented in Lemma 1 did not allow achieving a higher convergence order, unlike [29,31], where the equations contained integer time derivatives. This limitation can be removed by utilizing higher-order approximations of fractional derivatives (see, for example [46]). The authors believe that this study deserves a separate paper.

Since a significant part of the paper is devoted to a theoretical study of the constructed method, little attention has been paid to the features of the fractional-order stochastic Stokes–Darcy model. In particular, in subsequent works, the authors intend to study the influence of the orders of fractional derivatives on the fluid flow pattern. In addition, it will be interesting to study the influence of the spatial fractional derivative orders on the fluid flow in fractal porous media. In a previous study [47], the authors constructed a numerical method for a model space-fractional partial differential equation and proposed a parallel algorithm for its implementation. In future papers, the authors intend to conduct a comprehensive analysis of the time–space fractional generalization of the Stokes–Darcy model and address its parallel implementation.

When implementing a parallel ensemble approach on a test example, an almost seven-fold acceleration of the calculation was obtained. This result is close to the result obtained in [3], where the ensemble approach saved a massive 88% of CPU time in a stochastic problem with integer time derivatives. In contrast to the above work, in our problem most of the CPU time is spent on calculating the discrete fractional derivative with an increase in the number of time layers and the number of degrees of freedom. In subsequent papers, parallelization issues using other parallel architectures will be considered in more detail. This will make it possible to apply the proposed numerical method to model more realistic fluid flows in the case of uncertainties based on the fractional-order stochastic model.