A stabilized multiple time step method for coupled Stokes-Darcy flows and transport model

1.   Introduction

The models for the coupling of fluid flows in a porous medium domain and a free flow domain have a wide applications in environment science ^[15] and biofluid dynamics ^[18]. The fluid flow in the porous medium domain is described by Darcy equation, and the fluid flow in the free flow domain is modeled by Stokes equation. The two equations are coupled through the interface conditions, which connects the porous domain and the free flow domain.

For past years, stable and convergent numerical methods for the coupled Stokes and Darcy flows system can be found in many references, which have been deeply studied. For example, the coupled finite element methods ^[7,23,28], the domain decomposition methods ^[2,6,9,37], the conforming finite volume element method ^[25], the non-conforming finite element methods ^[32,41], the mortar finite element methods ^[1,12,21], the Lagrange multiplier methods ^[5,19,24], the mixed finite element method combining with the DG method ^[30,31], the DG method combining with mimetic finite difference method ^[26], the staggered DG method ^[42], the pseudospectral least squares method ^[22], the spectral method ^[39], the weak Galerkin methods ^[11], and many other numerical methods ^{[10,17,27,29,36]}.

The aim of this paper is to construct an efficient numerical algorithm for the Stokes-Darcy flows system coupled with the transport of a chemical. This kind of model can describe solute transport in the coupled flow region and porous media flow region, which appear in the research of human health and the environment. For example, the groundwater contamination, the pollution of groundwater by transport of contaminants through rivers. Cesmelioglu and Riviére study the existence and stability bounds of the weak solution with the fluid viscosity depending on the concentration for this model in ^[8]. For numerical methods, Vassilev and Yotov solve the flow equations through the domain decomposition method, and solve the transport equation by using local discontinuous Galerkin method ^[38]. The viscosity of the fluid in this study is assumed to be independent of the concentration. In ^[33], the authors study the numerical methods for the model with concentration-dependent viscosity, they propose a mixed weak formulation for the coupled flow problem, and use conforming piecewise linear finite element to approximate concentration.

In this paper, as an extended research of ^[33], we study the numerical methods for the Stokes-Darcy-Transport system with different time steps on coupled Stokes-Darcy flows system and the transport system. The viscosity in this study is assumed to depend on concentration. Since the Stokes-Darcy-Transport system is a multi-physics problem, the partial differential equations have different time scale reflected by the corresponding physical parameters, so it reminds us that we can use larger time step in the region with slower velocity. By using the multiple time step discrete finite element scheme, we can obtain the same optimal error estimation order as ^[33], and reduce computation, effectively improve computational efficiency. The multiple time step technique for the Stokes-Darcy was studied in ^[34,35]. In this study, we construct a stabilized mixed finite element method for the coupled Stokes-Darcy flows system by using the nonconforming piecewise linear Crouzeix-Raviart element for velocity and using piecewise constant function for pressure, and propose classical piecewise linear finite element method for the transport system. Under no assumption on the restriction about the time-step and spatial meshsize, we obtain the existence and uniqueness of the numerical solutions, and a prior error estimates. From the error analysis, we also derive the ratio of different time step sizes which should be proportional to the ratio of the physical parameters.

The rest of the article is organized as follows. In Section 2, we introduce the model problem and present the mixed weak formulation. In Section 3, we propose the multiple time step scheme with different time steps on Stokes-Darcy flows system and solute transport, and give the existence and uniqueness of the scheme by adding a stabilization term, which derive the discrete inf-sup condition. The error estimates for fluid velocity and concentration, and the relationship formula between physical parameters and the ratio of different time steps are presented in Section 4. In Section 5, we present some numerical examples to verify that the numerical results are in agreement with the theoretical analysis.

Throughout this paper we use $C$, with or without subscription, to denote a generic constant, which should have different values in different appearances.

2.   Model problem and weak formulation

Consider the model of a flow in a bounded domain $\Omega\subset R^L$($L$ = 2 or 3), consisting of a free region $\Omega_s$, where the flow is governed by the Stokes equations, and a porous medium domain $\Omega_d = \Omega\backslash \overline{\Omega}_s$, where the flow is governed by Darcy's law. The two regions are separated by a interface $\Gamma_I = \partial \Omega_s\cap\partial \Omega_d$. Let $\Gamma_l = \partial\Omega_l\backslash \Gamma_I(l = s, d)$. Each interface and boundary is assumed to be polygonal. We denote by ${\boldsymbol{n}}_{l}(l = s, d)$ the unit outward normal direction on $\partial\Omega_l(l = s, d)$, and on the interface ${\boldsymbol{n}}_s = -{\boldsymbol{n}}_d$.

Figure 1 gives a schematic representation of the geometry with $L = 2$.

In all of $\Omega$, denote the fluid velocity by ${\boldsymbol{u}}({\boldsymbol{x}}, t)$, the pressure by $p({\boldsymbol{x}}, t)$ and the concentration by $c({\boldsymbol{x}}, t)$. For any vector and scalar functions ${\boldsymbol{v}}$, $q$ defined in $\Omega$, it often plays different mathematical roles in $\Omega_s$ and $\Omega_d$, we will often need to distinguish them, especially their traces on $\Gamma_I$, thus define

Let $J = [0, T_1]$ be the time interval, and $\partial_t$ denotes the usual partial derivative $\partial/\partial t$, then in the free region $\Omega_s$, the equations of motion, continuity and mass transport can be written as

In the porous medium $\Omega_d$, the equations of motion, continuity and mass transport can be written as

where $\mu(c)$ is the concentration-dependent fluid viscosity, we will assume it by the quarter-power rule, see Eq (2.20). ${\boldsymbol{S}}$ is the deformation rate tensor and defined by ${\boldsymbol{S}}({\boldsymbol{v}}) = \frac{1}{2}(\nabla {\boldsymbol{v}}+\nabla {\boldsymbol{v}}^{T})$, ${\boldsymbol{f}}\in(L^2(\Omega))^L(L = 2\ or\ 3)$ is a term related to body forces, $d$ is the molecular diffusion coefficient, ${\boldsymbol{D}}$ is the diffusion-dispersion cofficient, $\lambda(c) = \frac{K({\boldsymbol{x}})}{\mu(c)}$, $K = \mbox{diag }\,{k_j}\in L^{\infty}(\Omega_d)^{L\times L}$ is the permeability of the medium, $q^I$ represents a source term, $q^P$ represents a sink term, $\phi\in L^{\infty}(\Omega_d)$ is the porosity of the medium, and $c^I$ is the injected concentration.

On the interface $\Gamma_I$, the conditions are imposed

Here $\gamma_j = \alpha_1/\sqrt{k_j}$, $\alpha_1$ is a parameter determined by experimental evidence. Equations (2.7), (2.10) and (2.11) represent continuity of mass flux and concentration, Eq (2.8) represents the balance of normal forces, Eq (2.9) is the Beavers-Joseph-Saffman condition. Moreover, $\tau_j$ denote a orthonormal system of tangent vectors on $\Gamma_I$.

To complete the system, we give the following boundary conditions

and initial conditions

Here

where ${\boldsymbol{I}}$ is the identity matrix.

Equations (2.1)–(2.16) consist of the coupled Stokes and Darcy flows system with an advection-diffusion equation that models transport of a chemical in which the viscosity is dependent on the concentration $c$.

To facilitate the subsequent discussion, we now make the following assumptions about some physical quantities in the system.

(1) The porosity of the medium $\phi({\boldsymbol{x}})$ and the permeability of the medium $K({\boldsymbol{x}})$ are uniformly bounded and positive defined in $\Omega_d$. There exists positive constants $\phi_{\max}$, $\phi_{\min}$, $k_{\max}$ and $k_{\min}$ such that

(2) The form of $\mu$ is assumed by the quarter-power rule

From (2.20), we know that $\mu(c)$ is bounded and monotone for concentration $c\in[0, 1]$

where $\mu_{\min} = \min\{\mu(1), \mu(0)\}$, $\mu_{\max} = \max\{\mu(1), \mu(0)\}$.

What's more, $\mu(c)$ is a Lipschitz continuous function for concentration $c\in[0, 1]$ with Lipschitz constant $\mu_L$.

Thus, from the assumptions (1) and (2), $\lambda$ is also bounded, monotone and Lipschitz continuous for concentration $c\in[0, 1]$, we can get the estimate inequality of $\lambda$ directly

(3) The source term and sink term $q^I, q^P\ge 0$ satisfy the compatibility condition

and $q^I, q^P\in L^{\infty}(J; L^2(\Omega_d))$. $c^I$ satisfies $0\le c^I\le 1\ a.e.\ in \ \Omega_d^T.$

(4) ${\boldsymbol{f}}(c)$ is a Lipschitz continuous function for concentration $c\in[0, 1]$ with Lipschitz constant $f_{L}$.

(5) The diffusion-dispersion coefficient ${\boldsymbol{D}}$ (refer to ^[16]) is taken to be

where $\widetilde{\tau}$ is a positive constant in $(0, 1)$, and represents the tortuosity of the porous medium, $d_l$ and $d_t$ are the longitudinal and transverse dispersion coefficients, respectively. For ${\boldsymbol{u}} = (u_1, \cdots, u_L)$, $|{\boldsymbol{u}}| = \sqrt{u^2_1+\cdots+u_L^2}$ and the matrices ${\boldsymbol{E}}, {\boldsymbol{E}}^{\bot}$ are given by

Usually $d_l$ is considerably larger than $d_t$, hence we assume $d_l > d_t$.

From the analysis of ^[16], the definitions (2.17) and (2.23), we know that for ${\boldsymbol{u}}, {\boldsymbol{v}}\in C(\overline{\Omega_d})^L$, there holds

Before giving the suitable weak formulations of the problems (2.1)$-$(2.16), we introduce some useful notations.

The Sobolev space $W^{m, n}$ and $L^q(J; W^{m, n}(\Omega))$ are defined in the usual way with the usual norm $\|\cdot\|_{W^{m, n}}$ and $\|\cdot\|_{L^q(J; W^{m, n}(\Omega))}$, where $0\le m < \infty$, $0\le n \le\infty$, $0\le q \le\infty$. When $n = 2$, we simply substitute $H^m(\Omega)$ for $W^{m, 2}(\Omega)$ with $\|\cdot\|_{m, \Omega} = \|\cdot\|_{W^{m, 2}(\Omega)}$, $|\cdot|_{m, \Omega} = |\cdot|_{W^{m, 2}(\Omega)}$. In particular, when $m = 0$, we have $L^2(\Omega) = H^{0}(\Omega)$, with $\|\cdot\|$ for $\|\cdot\|_{0, \Omega}$. Let $(\cdot, \cdot)_{\Omega}$ denote $L^2(\Omega)$, $L^2(\Omega)^L$ or $L^2(\Omega)^{L\times L}$ inner product or duality pairing. Also, $\|\cdot\|_{l}, |\cdot|_{l}, l = s, d$, will be the same with $\Omega$ replaced by $\Omega_l, l = s, d$. The $L^2(\Gamma_I)$ inner product or duality pairing is denoted by $\langle\cdot, \cdot\rangle_{\Gamma_I}$.

The special vector-function space is defined by

To present a variational form of the coupled problem, the spaces for the velocity, pressure and concentration are defined as follows

equipped with the norm

$Q = \{q\in L^2(\Omega):\int_{\Omega}qdx = 0\}$, with the norm $\|\cdot\|_Q = \|\cdot\|$, and $W = H^1(\Omega)$ with the norm $\|\cdot\|_W = \|\cdot\|_1$. Note that the vector valued functions in ${\boldsymbol{V}}$ have (weakly) continuous normal components on $\Gamma_I$ ^[4].

We also consider the convection term of the concentration equation as in ^[33]. For $z\in W$ it is clear that

So we have

Using (2.26), we now propose the following weak formulation of the coupled problems (2.1)–(2.16): find ${\boldsymbol{u}}(t)\in {\boldsymbol{V}}$, $p(t)\in Q$, $c(t)\in W$ such that for a.e. $t\in J$

where

Interface conditions (2.7)–(2.9) and (2.11) are posed weakly in the above variational form, while (2.10) is treated as an essential condition. Due to (2.18), $(\cdot, \cdot)_{\overline{\phi}}$ is an equivalent scalar product on $L^2(\Omega)$ and $\|c\|_{\overline{\phi}} = (c, c)^{1/2}_{\overline{\phi}}$ defines an equivalent norm on $L^2({\Omega})$.

3.   Finite element discretization and multiple time step method

In this section, we consider the finite element discretization of the coupled problem and propose a multiple time step method. We will use the nonconforming Crouzeix-Raviart piecewise linear finite element approximation for velocity, piecewise constant approximation for pressure, and continuous piecewise linear polynomial approximation for concentration, the penalizing term is the same as ^[32]. It should be noted that, the time step for velocity and pressure is different from concentration, which will increase the computing efficiency. The existence and uniqueness of a finite element solution of the discrete problem will also be established in this section. It should be mentioned that the existence and uniqueness is unconditionally satisfied, which benefit from using Crouzeix-Raviart finite element method to solve coupling problems.

Let $\mathscr{T}_h$ be a family of triangulations of $\Omega$ with nondegenerate elements, that means $\Omega$ is completely triangulated into triangles($L = 2$) or tetrahedron ($L = 3$). For any $T\in\mathscr{T}_h$, we denote the diameter of $T$ as $h_T$, $h = \max\limits_{T\in\mathscr{T}_h}h_T$ and the diameter of the sphere inscribed in $T$ as $\rho_T$. $\mathscr{T}_h$ is regular in the sense of Ciarlet ^[13], that is, there exists a constant $\sigma$ independent of $h$ and $T$ such that

Remark 3.1. $\Omega$ is completely triangulated into triangles ($L = 2$) or tetrahedron ($L = 3$), so each interface and boundary are assumed to be polygonal. This approximation makes it impossible for the numerical discrete scheme to fully match Figure 1. Figure 1 is intended to demonstrate the general model that satisfies the Eqs (2.1)$-$(2.16). Polygonal simplification was carried out during discretization to facilitate approximate calculation and solution. Further research is needed on the approximation of general smooth curves.

Assume each element $T\in\mathscr{T}_h$ is in either $\Omega_s$ or $\Omega_d$, denote $\mathscr{T}_h^s$ and $\mathscr{T}_h^d$ as the corresponding induced triangulations of $\Omega_s$ and $\Omega_d$. For any $T\in\mathscr{T}_h$, $\mathcal{E}(T)$ is denoted as the set of its edges ($L = 2$) or face ($L = 3$), and set $\mathcal{E}_h = \mathop{\cup}\limits_{T\in\mathscr{T}_h}\mathcal{E}(T)$. We split $\mathcal{E}_h$ into the following form

where $\Omega_s^{+} = \Omega_s\cup\Gamma_s$, $\mathcal{E}_h(\mathscr{S}) = \{E\in\mathcal{E}_h :E\subset \mathscr{S}\}$.

With every edge $E\in\mathcal{E}_h$ we associate a unit vector ${\boldsymbol{n}}_E$ such that ${\boldsymbol{n}}_E$ is orthogonal to $E$. For any $E\in\mathcal{E}_h$ and any piecewise continuous function $\varphi$, we denote by $[\varphi]_E$ the jump of $\varphi$ across $E$ in the direction ${\boldsymbol{n}}_E$:

Define the nonconforming Crouzeix-Raviart piecewise linear finite element space

and piecewise constant function space

where $P_m(T)$ is the space of the restrictions to $T$ of all polynomials of degree less than or equal to $m$. And define the classical Galerkin finite element space

it is clear that $W_h\subset W$.

Denote the operator $div_h\in \mathscr{L}({\boldsymbol{V}}_h, Q_h)\cap\mathscr{L}({\boldsymbol{V}}, Q)$ by $div_h{\boldsymbol{v}}|_T = \nabla\cdot{\boldsymbol{v}}|_T, \ \forall T\in \mathscr{T}_h.$ Then define some bilinear forms

and penalty term

with

where $(\cdot, \cdot)_T$ is $L^2(T)$, $(L^2(T))^L$ or $(L^2(T))^{L\times L}$ inner product or duality pairing, $h_E$ is the length ($L = 2$) or the diameter ($L = 3$) of E. Note that each edge ($L = 2$) or face ($L = 3$) of $\mathcal{E}_h$ only correlates with one jump term of $j({\boldsymbol{u}}, {\boldsymbol{v}})$.

Now we give the numerical scheme with different time steps. Assume that to each time level $\tau_m$ for concentration, there exists a time level $t_{n_m}$ for velocity and pressure. For simplicity, define uniform time levels as

where $\Delta t = \frac{T_1}{N}, \ \Delta \tau = \frac{T_1}{M}$, and $N = r\cdot M$, which means $n_m = r\cdot m$. Here $r$ is the ratio of the different time steps, we will give the relationship between the value of the ratio and the physical parameters of the equations in next section.

The relationship of the time steps is given in Figure 2.

Let

and

Set the initial approximations as

where the operators ${\boldsymbol{r}}_h$ and $\Xi_h$ will be defined in (3.17) and (4.3), the numerical algorithm of the fully discrete finite element scheme with different time steps is proposed as follows:

Here $\overline{C}_h^m = \min\{1, \max\{0, C_h^m\}\}\in[0, 1]$. Since $c\in[0, 1]$, so it is clear that

Each algebraic system of numerical schemes (3.8)–(3.10) in the numerical algorithm is linear. And the time step for velocity and pressure is different from concentration, which depending on the ratio $r$, so we can choose appropriate $r$ to improve computing efficiency.

Define the norm of ${\boldsymbol{V}}_h$:

and the space $Q_h$ is equipped with the norm $\|\cdot\|$.

Then giving the definition of ${\boldsymbol{Z}}_h$, which is the subspace of ${\boldsymbol{V}}_h$:

we have the following lemma.

Lemma 3.1. if ${\boldsymbol{v}}_h\in{\boldsymbol{Z}}_h$, then $div_h{\boldsymbol{v}}_h = 0$.

Proof. Since $div_h{\boldsymbol{v}}_h\in Q_h$ for ${\boldsymbol{v}}_h\in{\boldsymbol{V}}_h$, if ${\boldsymbol{v}}_h\in Z_h$, take $q_h = div_h{\boldsymbol{v}}_h$, we have $(div_h{\boldsymbol{v}}_h, div_h{\boldsymbol{v}}_h)_{\Omega} = 0$, the lemma follows.□

Lemma 3.2. The bilinear form $A_h(\cdot; \cdot, \cdot)$ is coercive on $Z_h$: there is a positive constant $\alpha > 0$ such that

where $\alpha\propto \frac{\mu_{\min}}{k_{\max}}$ depends on $\sigma$, and is independent of $h$, $\Delta t$ and $\Delta \tau$.

Proof. Follow Lemma 3.1, inequality (2.21) and the Korn's inequality for piecewise $H^1$ vector fields ^[3], we have

where $C$ is a positive constant which depends only on $\sigma$ and $\mu_{\min}$. Since $\gamma_j = \alpha_1/\sqrt{k_j}$ and follow inequality (2.22), we have $\alpha\propto \frac{\mu_{\min}}{k_{\max}}$.□

Lemma 3.3. The bilinear forms $A_h(\cdot; \cdot, \cdot)$ and $b_h(\cdot, \cdot)$ are Bounded in ${\boldsymbol{V}}_h+{\boldsymbol{V}}$: there is a positive constant $A$ such that

where $A\propto \frac{\mu_{\max}}{k_{\min}}+\mu_{\max}$ is independent of $h$, $\Delta t$ and $\Delta \tau$.

Proof. The lemma follows from Hölder inequality.□

Define the space

we give the definition of the Crouzeix-Raviart interpolation operator ${\boldsymbol{r}}_h:{\boldsymbol{X}}\to{\boldsymbol{V}}_h, \ \forall t\in J$

From the interpolation theory for ${\boldsymbol{r}}_h$^[14], we have the standard Bramble-Hilbert estimation

where $R$ is independent of $h$.

From ^[20] and the definition of ${\boldsymbol{r}}_h{\boldsymbol{v}}$, we know that the discrete inf-sup condition is satisfied. The proof can be found in ^[33].

Lemma 3.4. There exists a constant $\beta > 0$ depending on $\sigma$, $\mu_{\min}$ such that

Lemma 3.5. The nonlinear form $d(\cdot; \cdot, \cdot)$ is positive definite, that is

Proof. From the definition (2.36), we have

proof of positive definiteness of $d(; ,)$.□

From Lemmas 3.2–3.4 and the abstract theory of mixed problem ^[4,21], we derive the existence and uniqueness of velocity ${\boldsymbol{u}}$ and pressure $p$ in both free region $\Omega_s$ and porous medium domain $\Omega_d$. From Lemma 3.5 and the fact that the numerical scheme (3.10) is linear, we can also get the existence and uniqueness of concentration $c$. Thus we drive the existence and uniqueness of the the fully discrete finite element scheme with different time steps.

Theorem 3.2. There exists unique solution $({\boldsymbol{u}}^n_h, p^n_h)\in{\boldsymbol{V}}_h\times Q_h, \ n = 1, 2, \dots, N$, and $C_h^m\in W_h, \ m = 1, 2, \dots, M$ satisfying the numerical schemes (3.8)$-$(3.10) in ${\boldsymbol{Algorithm}}^*$.

4.   Error estimates

In this section, we derive the error estimates for ${\boldsymbol{Algorithm}}^*$ and the relationship between the value of $r$ and the physical parameters of the equations. We first define two projection operators and give corresponding estimation.

Define the projection operator $\Pi_h: Q\to Q_h, \ \forall t\in J$ by

then

with $P$ depending only on $\sigma$ and $L$, so that

Define the elliptic projection operator $\Xi_h: W\to W_h, \ \forall t\in J$ by (one can also find similar definition in ^[16,40])

From the theory of finite element methods for elliptic problem, when $c$ is sufficiently smooth, there hold

where the constants $C$ and $C_c$ are independent of $h$.

From the trace theorem, we have

and from ^[13], the following approximation property satisfied

where $n = 2, \infty$, $\varphi\in U^{2, n} = \left\{\varphi\in W:\varphi_s\in W^{2, n}(\Omega_s), \varphi_d\in W^{2, n}(\Omega_d)\right\}$.

For estimating the approximation error, we assume that the exact solutions satisfy the following smoothness conditions:

Remark 4.1. The assumptions (4.9)$-$(4.11) of smoothness for exact solutions may not be sufficient for some special practical problems, but they are essential for error estimation analysis. Studying error estimation analysis approach that satisfy low smoothness of exact solutions will be our future research work.

The smoothness assumption of ${\boldsymbol{u}}$ implies $j({\boldsymbol{u}}, {\boldsymbol{v}}_h) = 0$, thus

At a certain time point $n, n = 1, 2, \dots, N$, let Eqs (2.1) and (2.4) do inner product with ${\boldsymbol{v}}_h\in V_h$ respectively, and add up the results. By using Green's formula on each $T\in\mathscr{T}_h$, the interface conditions (2.8), (2.9) and Eq (4.12), we obtain

and let Eqs (2.2) and (2.5) do inner product with $q_h\in Q_h$ respectively, we obtain

Split the error as

then subtract (4.13) and (4.14) from (3.8) and (3.9), respectively, we can obtain the following equations by inserting ${\boldsymbol{r}}_h{\boldsymbol{u}}$ and $\Pi_h p$

where (4.16) is obtained from the definition (3.16).

Lemma 4.1. Suppose that the analytical solution satisfies properties (4.9)$-$(4.11), then there exists a positive constant $\hat{C}_1$ independent of $h$, $\Delta t$ and $\Delta \tau$, and exists a positive constant $\hat{C}_2$ independent of $h$, $\Delta t$, $\Delta \tau$, and independent of physical parameters $\mu_{\max}$, $\mu_{\min}$, $k_{\min}$ and $k_{\max}$, such that for $0\le l\le M-1$,

Proof. Taking ${\boldsymbol{v}}_h = 2\Delta t e_{{\boldsymbol{u}}, h}^{n+1}$ in (4.15), and combining with (4.16), we have

As inner product has the following property

using (3.13), (4.19) and summing over (4.18) with $n = n_m, n_m+1, \dots, n_{m+1}-1$, we obtain

Now we estimate the terms $S_1, S_2, \dots, S_7$ individually by using Cauchy-Schwarz inequality.

Integrating by parts about $t$, we obtain

where $\partial_t^2 {\boldsymbol{u}} = \partial^2 {\boldsymbol{u}}/\partial t^2$, and combine with (3.17), we obtain

From (3.11), (4.7), (4.12) and the assumption that $\mu$ is Lipschitz continuous, we obtain

where $C$ is independent of $h$, $\Delta t$ and independent of physical parameters $\mu_{\max}, \ \mu_{\min}, \ k_{\min}, \ k_{\max}$.

We estimate the last three items of (4.23) in turn, since for $n = n_m, \dots, n_{m+1}-1$,

and $n_{m+1}-n_m = r$, we have

Thus

Take (4.25)$-$(4.27) back to (4.23), the estimate of $T_3$ is

From (3.14) and (3.17), we obtain

Similarly, from (3.15) and (4.2)

In the same way as estimating $T_3$, follow the the assumption that ${\boldsymbol{f}}(c)$ is Lipschitz continuous, we obtain

Using the proof of Lemma 5.3 in ^[37], we have

where $C$ is a constant depending only on $\sigma$.

Taking the estimates of $S_1, S_2, \dots, S_7$ back to (4.20) and summing over with $m = 0, 1, \dots, l$, meanwhile, since $n_0 = 0$, we obtain

Using the property ${\boldsymbol{u}}_{h, s}^{0} = {\boldsymbol{r}}_h {\boldsymbol{u}}_{s, 0}$ and applying (3.17), we complete the proof.□

Now we consider the error estimate for concentration. At a certain time point $n_{m+1}, m = 0, 1, \dots, M$, subtracting (3.10) from (2.29) and using (4.3) and (2.37), we have

Lemma 4.2. Suppose that the analytical solution satisfies properties (4.9)$-$(4.11), if $\Delta t$ and $\Delta \tau$ are sufficiently small, then there exists a positive constant $\hat{C}_3$ independent of $h$, $\Delta t$ and $\Delta \tau$, and exists a positive constant $\hat{C}_4$ independent of $h$, $\Delta t$ and $\Delta \tau$, and independent of physical parameters $\mu_{\max}$, $\mu_{\min}$, $k_{\min}$, $k_{\max}$, $D_{\min}$ and $\phi_{\max}$, such that for $0\le l\le M-1$

Proof. Taking $z_h = 2\Delta \tau e_{c, h}^{n_{m+1}}$ in (4.34), we have

using (2.24) and (4.19), we obtain

Using Cauchy-Schwarz inequality, now we estimate $H_1, H_2, \dots, H_7$ term by term.

Similarly as the estimates for $S_1$ and $S_2$, we have

Using (4.6) and (2.25), we have

Here

From the discussion of (4.24), the second term of (4.41) can be estimated as follows

Taking (4.41) and (4.42) back to (4.40), we get the estimate of $H_3$

Similarly, $H_4$ and $H_5$ can also be estimated as follows

The last two terms $H_6$ and $H_7$ can be estimated directly

Collecting each estimate above and summing over with $m = 0, 1, \dots, l$, due to the fact that $e_{c, h}^{n_0} = 0$ and $(\cdot, \cdot)_{\overline{\phi}}$ is an equivalent scalar product on $L^2(\Omega)$, we have

If $\Delta t$ and $\Delta \tau$ are sufficiently small, applying the discrete Gronwalll inequality and using (4.17) in Lemma 4.1, we can obtain

Here positive constants $\tilde{C}_1$ and $\tilde{C}_2$ are independent of $h$, $\Delta t$ and $\Delta \tau$. Combining with the estimate for $c-\Xi_h c$ in (4.4) and the setting $C_h^{m} = c_h^{n_m}$, we get (4.35). □

Combining the above two lemmas, we give the error estimate conclusion solely in terms of $\Delta t$ under appropriate ratio $r$, which determine time steps for each equation.

Theorem 4.2. Suppose that the analytical solution satisfies properties (4.9)$-$(4.11). Let

where the positive constant $\hat{C}_5$ is independent of physical parameters $\mu_{\max}$, $\mu_{\min}$, $k_{\min}$, $k_{\max}$, $D_{\min}$ and $\phi_{\max}$. If $\Delta t$ is sufficiently small, then there exists a positive constant $\hat{C}_6$ independent of $h$ and $\Delta t$, and exists a positive constant $\hat{C}_7$ independent of $h$ and $\Delta t$, and independent of physical parameters $\mu_{\max}$, $\mu_{\min}$, $k_{\min}$, $k_{\max}$, $D_{\min}$ and $\phi_{\max}$, such that for $0\le l\le M-1$

Proof. From Lemmas 4.1 and 4.2, $\alpha\propto \frac{\mu_{\min}}{k_{\max}}$ in Lemma 3.2 and $\Delta \tau = r\cdot\Delta t$, we have

where $\hat{C}_8 = \max\{\hat{C}_2, \hat{C}_4\}$. We derive the estimate (4.51) by taking $\hat{C}_7 = \hat{C}_8(1+(\hat{C}_5)^2)$. □

5.   Numerical experiments

In this section we give some numerical experiments to demonstrate the error estimates results obtained in the previous section. For simplicity, take the unit square domain $\Omega = [0, 1]\times[0, 1]$, and choose $\Omega_s = [0, 1/2]\times[0, 1]$ and $\Omega_d = [1/2, 1]\times[0, 1]$ with the interface $\Gamma_I = \{1/2\}\times(0, 1)$. The time interval is $J = [0, 2]$. Unless specified otherwise, for all the numerical experiments the values of the parameters are assigned as $\phi = 0.3, \ \gamma = d = q^I = 1$, $K = \frac{1}{2}I$, and

which can be obtained that $D_{min} = 1$. $\mathscr{T}_h$ is taken as a uniform grid. The upper and lower bounds of $\mu$ in (2.20) are $\mu(1) = 1.2$ and $\mu(0) = 0.1$. Since it is difficult to construct the exact solutions that satisfy the entire coupled Stokes and Darcy flows with mass transfer (2.1)$-$(2.16), especially the interface conditions. To solve this, we generalize the interface conditions (2.8), (2.9) and (2.11) to include nonhomogeneous terms.

where $\eta_1$, $\eta_2$ and $\eta_3$ are given functions on $\Gamma_I$ according to the analytical solutions. The variational form for this modified system will only include two additional terms $-\langle\eta_1, {\boldsymbol{v}}_s\cdot{\boldsymbol{n}}_s\rangle_{\Gamma_I}+\langle\eta_2, {\boldsymbol{v}}_s\cdot{\boldsymbol{\tau}}_I\rangle_{\Gamma_I}$ to the right-hand side of (2.27), and one additional term $-\langle\eta_3, z\rangle_{\Gamma_I}$ to the right-hand side of (2.29). All the right-hand terms and boundary conditions are selected according to the analytical solution.

Notation wise, We define the following symbols to represent the computational errors.

With these notations, we compute the convergence rate of the approximate solutions by

with meshsizes $h_1$ and $h_2$. $E^{h_i}, i = 1, 2$ are any of the errors described in (5.1).

We will use two different numerical examples to verify the convergence rates of the multiple time step discrete finite element schemes (3.8)$-$(3.10). In each example we also compare the errors, convergence rates and CPU calculation time for concentration (unit:second) of the single time step discrete finite element scheme (which means $r = 1$, $M = N$, $\Delta t = \Delta \tau$) and the multiple time step discrete finite element scheme.

We list and plot all the errors in (5.1) and calculate the corresponding convergence rates for $h$, which is set to $2^{-m}, m = 2, 3, 4, 5$. The number of time steps is $N = 100$, and $r$ in the multiple time step discrete finite element scheme is set to be $r = 20$. To make this setting valid, we take $\hat{C}_5 \approx 19.73$, so that $M = 5$.

Example 5.1. In this example, the analytical velocity, pressure and concentration in Stokes domain and Darcy domain are listed as follows

The numerical results, convergence rates and CPU calculation time of Example 5.1 are listed in Tables 1 and 2. The convergence rates are plotted with respect of nodes on each direction in Figure 3. The numerical velocity quiver, numerical pressure distribution and numerical concentration distribution of the multiple time step discrete finite element scheme with $h = 1/32$ are plotted in Figures 4 and 5.

Example 5.2. In this example, the analytical velocity, pressure and concentration in Stokes domain and Darcy domain are listed as follows

The numerical results, convergence rates and CPU calculation time of Example 5.2 are listed in Tables 3 and 4. The convergence rates are plotted with respect of nodes on each direction in Figure 6. The numerical velocity quiver, numerical pressure distribution and numerical concentration distribution of the multiple time step discrete finite element scheme with $h = 1/32$ are plotted in Figures 7 and 8.

From Tables 1–4, Figures 3 and 4, we find that the numerical results are consistent with the theoretical analysis. The convergence rates for the pressure and velocity on $\Omega$ in $l^{\infty}(L^2)$ norm and the velocity on $\Omega_s$ in $l^{\infty}(H^1)$ norm are first-order, and the convergence rate for the velocity on $\Omega_s$ in $l^{\infty}(L^2)$ norm and the concentration on $\Omega$ in $l^{\infty}(L^2)$ norm are at least first-order. All the convergence rates are optimal. About the superconvergence for the velocity on $\Omega_s$ in $l^{\infty}(L^2)$ norm and the concentration on $\Omega$ in $l^{\infty}(L^2)$ norm, we have no approach to give the corresponding analysis temporarily, the theoretical investigation of the phenomena will be our future work.

From the comparison of Tables 1 and 2, Tables 3 and 4, we find that, the errors are similar in every discretization parameter, and the convergence rates are the same, but the multiple time step discrete finite element scheme costs less CPU time. And from Figures 4 and 5 and Figures 7 and 8, we can see that, it can clearly reflect the modification of velocity, pressure and concentration by using multiple time step discrete finite element scheme with fewer time steps for concentration. The comparison verify that the multiple time step method is useful to increase computational efficiency.

6.   Conclusions

This article presents a stabilized finite element algorithm with different time steps on different physical variables for the coupled Stokes-Darcy flows system with the solution transport (2.1)–(2.16). We use nonconforming piecewise linear Crouzeix-Raviart element and piecewise constant to approximate velocity and pressure in the coupled Stokes-Darcy flows system, and use conforming piecewise linear finite element to approximate concentration in the transport system. The time derivatives are discretized with different step sizes in these two systems. The existence and uniqueness of the approximate solution are proved unconditionally satisfied. The order of convergence and error estimates are given, which also provide a guidance on the ratio of time step sizes with respect to the ratio of the physical parameters. We present two numerical examples to verify that the numerical results are in agreement with the theoretical analysis. The results obtained through our numerical examples indicate that we can obtain optimal error estimation order by using the multiple time step discrete finite element scheme provided in this article and reduce computation, effectively improve computational efficiency.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The author would like to thank the referees for the helpful suggestions.

The work is supported by the National Natural Science Foundation of China No.12001308 and the Natural Science Foundation of Shandong Province No.ZR2021QF040.

Conflict of interest

The authors declare no conflicts of interest in this paper.