1 Introduction

There are many multimodeling problems in real applications of complex systems. They consist of multiple models in different regions coupled through interface conditions. The local models may be varied in type, scale, control variable, and many other physical and mathematical properties. In this paper, we focus on the coupled fluid flow and porous media flow modeled by the non-stationary Stokes-Darcy problem. There is a rich literature on the mathematical analysis, numerical methods and applications for this model, see, e.g., [2–10] and the references therein. Among them, the decoupled method might be one of the most popular approaches for solving the multimodeling problems because the decoupled method makes the existing single-model solvers applicable locally with little extra computational and software overhead. Other appealing reasons were discussed in [11, 12].

In [1], authors developed a decoupled method for the Stokes-Darcy model based on the numerical solutions from previous time level and established the corresponding error analysis. Unfortunately the estimates for \({\mathbf{u}}_{f}\) and ϕ are not optimal, namely, the order is \(\mathcal{O}(\Delta t^{\frac{1}{2}})\). These estimates may be improved to \(\mathcal{O}(\Delta t)\), as suggested by numerical experiments in [1]. It motivates us to propose some new decoupled algorithms and derive the optimal estimates of the order of \(\mathcal{O}(\Delta t)\) for both \({\mathbf{u}}_{f}\) and ϕ.

The rest of the paper is organized as follows. A coupled non-stationary Stokes-Darcy model and its weak formulation are introduced in Section 2. Numerical algorithms, including the coupled scheme and decoupled schemes, are developed in Section 3. The stabilities of these developed algorithms are provided in Section 4. Convergence is derived in Section 5 to show that these new decoupled algorithms keep the same order of approximation accuracy as the coupled method. Numerical results are reported in Section 6.

2 The non-stationary Stokes-Darcy model

Let us consider a fluid in \(\Omega_{f}\) coupled with a porous media flow in \(\Omega_{p}\), where \(\Omega_{f}, \Omega_{p}\in\mathbb{R}^{d} \) (\(d=2\mbox{ or }3\)) are bounded domains, \(\Omega_{f}\cap\Omega_{p}=\emptyset\), and \(\overline{\Omega}_{f}\cap\overline{\Omega}_{p}=\Gamma\). Denote by \(\Omega=\Omega_{f}\cup\Omega_{p}, {\mathbf{n}}_{f}\) and \({\mathbf{n}}_{p}\) the unit outward normal directions on \(\partial\Omega _{f}\) and \(\partial\Omega_{p}\), respectively. \({\boldsymbol {\tau}}_{i},i=1,\ldots,d-1\) the unit tangential vectors on the interface Γ and \({\mathbf{n}}_{f}=-{\mathbf{n}}_{p}\) on Γ.

Let \(T>0\) be a finite time. The fluid motion is governed by the Stokes equations:

$$ \left \{ \textstyle\begin{array}{@{}l@{\quad}l} \frac{\partial{\mathbf{u}}_{f}}{\partial t}-\nu\Delta{\mathbf {u}}_{f}+\nabla p_{f}={\mathbf{g}}_{f} &\mbox{in } \Omega_{f}\times (0,T],\\ \nabla\cdot{\mathbf{u}}_{f}=0 &\mbox{in } \Omega_{f}\times(0,T],\\ {\mathbf{u}}_{f}({\mathbf{x}},0)={\mathbf{u}}_{f}^{0}({\mathbf{x}}) &\mbox{in } \Omega_{f}, \end{array}\displaystyle \right . $$
(2.1)

where \({\mathbf{u}}_{f}({\mathbf{x}},t)\) represents the velocity of the fluid flow in \(\Omega_{f}\), \(p_{f}({\mathbf{x}},t)\) the kinetic pressure, \({\mathbf{g}}_{f}\) the external force, and \(\nu>0\) the kinematic viscosity.

The porous media flow motion is governed by the following equations [13, 14]:

$$ \left \{ \textstyle\begin{array}{@{}l@{\quad}l} S_{0} \frac{\partial\phi}{\partial t}+\nabla\cdot{\mathbf{q}}=g_{p} &\mbox{in } \Omega_{p}\times(0,T],\\ {\mathbf{q}}=-{\mathbf{K}}\nabla\phi &\mbox{in } \Omega_{p}\times(0,T], \ \mbox{(Darcy law)}\\ {\mathbf{u}}_{p}=\frac{{\mathbf{q}}}{n} &\mbox{in } \Omega_{p}\times(0,T],\\ \phi({\mathbf{x}},0)=\phi^{0}({\mathbf{x}}) &\mbox{in } \Omega_{p}, \end{array}\displaystyle \right . $$
(2.2)

where \(\phi({\mathbf{x}},t)\) is the piezometric head, q is the specific discharge defined as the volume of the fluid flowing per unit time through a unit cross-sectional area normal to the direction of the flow, \({\mathbf{u}}_{p}\) is the fluid velocity in \(\Omega_{p}\), \(S_{0}\) is the specific mass storativity coefficient, K is the hydraulic conductivity tensor, n is the volumetric porosity, and \(g_{p}\) is the source term. Note that \(\phi=z+\frac{p_{p}}{\rho g}\), the sum of elevation head plus pressure head, where z is the elevation from a reference level, \(p_{p}\) is the pressure in \(\Omega_{p}, \rho\) is the density of the fluid, and g is the gravity acceleration. Without loss of generality, we assume \(z=0\). Furthermore, we assume that \({\mathbf{K}}=\operatorname{diag}(K,\ldots,K)\) with \(K\in L^{\infty}(\Omega_{p}), K>0\), which implies that the porous media is homogeneous. Finally, by using Darcy’s law in (2.2), the first continuity equation in (2.2) in \(\Omega_{p}\) can be written in the parabolic form

$$ S_{0} \frac{\partial\phi}{\partial t}-\nabla\cdot({\mathbf {K}}\nabla \phi)=g_{p} \quad \mbox{in } \Omega_{p}\times(0,T]. $$
(2.3)

A key part in a mixed model is the interface coupling conditions. For the Stokes-Darcy model, the following interface conditions have been extensively studied and used in the literature [15–20]

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l@{\quad}l} {\mathbf{u}}_{f}\cdot{\mathbf{n}}_{f}+{\mathbf{u}}_{p}\cdot{\mathbf{n}}_{p}=0 &\mbox{on } \Gamma \times(0,T],\\ p_{f}-\nu{\mathbf{n}}_{f} \frac{\partial{\mathbf{u}}_{f}}{\partial{\mathbf {n}}_{f}}=\rho g\phi&\mbox{on } \Gamma\times(0,T],\\ -\nu{\boldsymbol {\tau}}_{i} \frac{\partial{\mathbf{u}}_{f}}{\partial{\mathbf{n}}_{f}} = \frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}{\mathbf{u}}_{f}\cdot{\boldsymbol {\tau}}_{i}, i=1,\ldots,d-1 &\mbox{on } \Gamma\times(0,T]. \end{array}\displaystyle \right . \end{aligned}$$
(2.4)

Here α is a positive parameter depending on the properties of the porous medium and must be experimentally determined. The first interface condition ensures the mass conservation across the interface Γ, and using the second and third equations in (2.2), it can be rewritten as

$$\begin{aligned} {\mathbf{u}}_{f}\cdot{\mathbf{n}}_{f}=\frac{{\mathbf{K}}}{n} \frac{\partial\phi}{\partial {\mathbf{n}}_{p}} \quad\mbox{on } \Gamma\times(0,T]. \end{aligned}$$

The second one is a balance of normal forces across the interface. The third one states that the slip velocity along Γ is proportional to the shear stress along Γ.

Several types of boundary conditions for this coupled model are discussed in [15]. In this paper, we consider the homogeneous Dirichlet boundary conditions for the coupled model, that is,

$$\begin{aligned} {\mathbf{u}}_{f}={\mathbf{0}} \quad\mbox{on } \partial \Omega_{f}\backslash\Gamma \quad\mbox{and} \quad \phi=0 \quad\mbox{on } \partial\Omega_{p}\backslash\Gamma. \end{aligned}$$

Denote \(W=H_{f}\times H_{p}\) and \(Q=L^{2}(\Omega_{f})\), where

$$\begin{aligned} H_{f}=\bigl\{ {\mathbf{v}}\in\bigl(H^{1}( \Omega_{f})\bigr)^{d}\mid {\mathbf{v}}=0 \quad\mbox{on } \partial \Omega _{f}\backslash\Gamma\bigr\} \end{aligned}$$

and

$$\begin{aligned} H_{p}=\bigl\{ \psi\in H^{1}(\Omega_{p})\mid \psi=0 \quad\mbox{on } \partial\Omega _{p}\backslash\Gamma\bigr\} . \end{aligned}$$

The space \(L^{2}(D)\), where \(D=\Omega_{f}\) or \(\Omega_{p}\), is equipped with the usual \(L^{2}\)-scalar product \((\cdot,\cdot)\) and \(L^{2}\)-norm \(\|\cdot\|_{L^{2}(D)}\). The spaces \(H_{f}\) and \(H_{p}\) are equipped with the following norms:

$$\begin{aligned}& \|{\mathbf{u}}_{f}\|_{H_{f}}=\|\nabla{\mathbf{u}}_{f} \|^{2}_{L^{2}(\Omega_{f})}=(\nabla {\mathbf{u}}_{f},\nabla{ \mathbf{u}}_{f})_{\Omega_{f}} \quad \forall {\mathbf{u}}_{f} \in H_{f}, \\& \|\phi\|_{H_{p}}=\|\nabla\phi\|^{2}_{L^{2}(\Omega_{p})}=(\nabla\phi, \nabla\phi )_{\Omega_{p}} \quad \forall \phi\in H_{p}. \end{aligned}$$

We equip the space W with the following norms: \(\forall u=({\mathbf {u}}_{f},\phi)\in W\)

$$\begin{aligned}& \|u\|_{0}^{2}=n({\mathbf{u}}_{f},{ \mathbf{u}}_{f})_{\Omega_{f}}+\rho gS_{0}(\phi,\phi )_{\Omega_{p}}, \\& \|u\|_{W}^{2}=n\nu(\nabla{\mathbf{u}}_{f}, \nabla{\mathbf{u}}_{f})_{\Omega_{f}}+\rho gK(\nabla\phi,\nabla \phi)_{\Omega_{p}}\thickapprox\|\nabla u\|_{0}^{2}, \end{aligned}$$

where \((\cdot,\cdot)_{D}\) refers to the scalar product \((\cdot,\cdot)\) in the corresponding domain D for \(D=\Omega_{f}\) or \(\Omega_{p}\). For simplicity, we assume \(n,\rho,g,S_{0},\nu\) and K are constants.

We also recall the Poincaré and trace inequalities that are useful in the following analysis. There exist constants \(C_{p}\) and \(C_{t}\) which only depend on Ω such that

$$\begin{aligned} \|v\|_{L^{2}(D)}\leq C_{p}\|v\|_{H^{1}(D)}, \qquad \|v\|_{L^{2}(\Gamma)}\leq C_{t}\|v\|^{1/2}_{L^{2}(D)} \|v\|^{1/2}_{H^{1}(D)}\quad \forall v\in H^{1}(D). \end{aligned}$$
(2.5)

The weak formulation for the non-stationary Stokes-Darcy problem reads as follows: Find \(u=({\mathbf{u}}_{f},\phi)\in W, p_{f}\in Q\) such that for all \(t\in(0,T]\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l@{\quad}l} ( \frac{\partial u}{\partial t},v )+a(u,v)+b(v,p_{f})=(f,v) &\forall v=({\mathbf{v}},\psi)\in W,\\ b(u,q)=0 & \forall q\in Q,\\ u({\mathbf{x}},0)=u^{0}, \end{array}\displaystyle \right . \end{aligned}$$
(2.6)

where

$$\begin{aligned}& \biggl( \frac{\partial u}{\partial t},v\biggr)=n\biggl(\frac {\partial{\mathbf{u}}_{f}}{\partial t},{\mathbf{v}} \biggr)+\rho gS_{0}\biggl(\frac{\partial\phi }{\partial t},\psi\biggr), \\& a(u,v)=a_{\Omega}(u,v)+a_{\Gamma}(u,v)=a_{\Omega_{f}}({ \mathbf{u}}_{f},{\mathbf {v}})+a_{\Omega_{p}}(\phi, \psi)+a_{\Gamma}(u,v), \end{aligned}$$

with

$$\begin{aligned}& a_{\Omega_{f}}({\mathbf{u}}_{f},{\mathbf{v}})=n \int_{\Omega_{f}}\nu\nabla{\mathbf {u}}_{f}\cdot\nabla{ \mathbf{v}} +n\sum_{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot {\mathbf{K}}{\boldsymbol {\tau}}_{i}}}({\mathbf{u}}_{f}\cdot{ \boldsymbol {\tau}}_{i}) ({\mathbf{v}}\cdot{\boldsymbol {\tau}}_{i}); \\& a_{\Omega_{p}}(\phi,\psi)=\rho g \int_{\Omega_{p}}{\mathbf{K}}\nabla\phi\cdot \nabla\psi; \qquad a_{\Gamma}(u,v)=n\rho g \int_{\Gamma}(\phi{\mathbf{v}}\cdot{\mathbf{n}}_{f}-\psi {\mathbf{u}}_{f}\cdot{\mathbf{n}}_{f}); \\& b(v,p_{f})=-n \int_{\Omega_{f}}p_{f}\operatorname{div} {\mathbf{v}}; \qquad (f,v)=n \int_{\Omega_{f}}{\mathbf{g}}_{f}\cdot{\mathbf{v}}+\rho g \int_{\Omega _{p}}g_{p}\psi. \end{aligned}$$

It is well known [4] that \(a_{\Omega_{f}}(\cdot,\cdot), a_{\Omega _{p}}(\cdot,\cdot)\) and \(a_{\Gamma}(\cdot,\cdot)\) are continuous, and \(a(\cdot,\cdot)\) is coercive. Furthermore, \(a_{\Omega_{f}}(\cdot,\cdot )\) and \(a_{\Omega_{p}}(\cdot,\cdot)\) are symmetric,

$$\begin{aligned} a_{\Gamma}(u,v)=-a_{\Gamma}(v,u) \quad \mbox{and} \quad a_{\Gamma }(u,u)=0 \quad \forall u,v\in W. \end{aligned}$$
(2.7)

The well-posedness of the model problem (2.6) can be found in [2, 4, 5] for the stationary case. We also proceed to apply the Babuska-Brezzi theory to prove that (2.6) is well posed for the non-stationary case. After assuming the inf-sup condition for \(b(\cdot,\cdot)\), we restrict the monolithic formulation (2.6) to the null space of \(b(\cdot,\cdot)\). With the help of Riesz representation theorem, we can define an operator \(u\mapsto Au\) in a standard way by \((Au,v)=a(u,v)\) for the bilinear form \(a(\cdot ,\cdot)\). Then it follows form continuity and coercivity and the Lax-Milgram theorem that this operator is maximal monotone. As a consequence, thanks to the Hille-Yoshida theorem, we can obtain the existence of solution for the evolutionary problem.

Lemma 2.1

Assume that

$$\begin{aligned} {\mathbf{g}}_{f}\in L^{2}\bigl(0,T,L^{2}( \Omega_{f})^{d}\bigr), \qquad g_{p}\in L^{2}\bigl(0,T,L^{2}(\Omega _{p})^{d} \bigr), \quad {\mathbf{K}}\in L^{\infty}(\Omega_{p})^{d\times d}, \end{aligned}$$

and K is uniformly bounded and positive defined in \(\Omega_{p}\), i.e., there exist constants \(k_{\min},k_{\max}>0\) such that

$$\begin{aligned} k_{\min}|x|^{2}\leq{\mathbf{K}}x\cdot x\leq k_{\max}|x|^{2} \quad \textit{a.e. } x\in\Omega_{p}. \end{aligned}$$

In addition, let \({\mathbf{u}}^{0}_{f}\in L^{2}(\Omega_{f})^{d}, \phi^{0}\in L^{2}(\Omega _{p})\), then the solution \(({\mathbf{u}}_{f},p_{f},\phi)\in(L^{2}(0,T,H_{f})\cap H^{1}(0,T,L^{2}(\Omega_{f})^{d}))\times L^{2}(0,T,Q)\times L^{2}(0,T,H_{p})\) of (2.1)-(2.2) is also the solution to (2.6). Conversely the solution of (2.6) satisfies (2.1)-(2.2).

3 Numerical algorithms

Let \(W_{h}=H_{fh}\times H_{ph} \subset W\) and \(Q_{h}\subset Q\) denote the finite element subspaces. The finite element spaces \(H_{fh}\) and \(Q_{h}\) approximating velocity and pressure in the fluid flow region are assumed to satisfy the well-known discrete inf-sup condition [21]: There exists a constant \(\beta>0\) independent of h, such that there exists a \(v_{h}\in W_{h}\) for all \(q_{h}\in Q_{h}\),

$$\begin{aligned} b(v_{h},q_{h})\geq\beta\|v_{h} \|_{W}\|q_{h}\|_{Q}. \end{aligned}$$
(3.1)

Several families of finite element spaces designed for the Stokes problem are provided in [4, 21]. They all satisfy the discrete inf-sup condition (3.1) and can thus be applied for \(H_{fh}\) and \(Q_{h}\). Standard finite element approximations of \(H^{1}(\Omega_{p})\) can be applied for \(H_{ph}\) in the porous media flow region. For illustration, we assume that the finite element spaces of the first-order approximation \(\mathcal{O}(h)\) are used in the fluid flow, such as the well-known MINI elements, and the porous media flow regions, such as the linear Lagrangian elements. The corresponding inverse estimates are well known:

$$\begin{aligned} \|\nabla{\mathbf{v}}_{h}\|_{0}\leq C_{in}h^{-1} \|{\mathbf{v}}_{h}\|_{0}\quad \forall {\mathbf {v}}_{h}\in H_{fh}; \qquad \|\nabla\psi_{h} \|_{0}\leq C_{in}h^{-1}\|\psi_{h} \|_{0} \quad \forall \psi_{h}\in H_{ph}. \end{aligned}$$

We also introduce a subspace \(V_{h}\) of \(W_{h}\) defined by

$$\begin{aligned} V_{h}=\bigl\{ v_{h}\in W_{h}:b(v_{h},q_{h})=0 \ \forall q_{h}\in Q_{h}\bigr\} , \end{aligned}$$

and the projection \(R_{h}:v=({\mathbf{v}},\psi)\in W\mapsto R_{h}v=(R_{h}{\mathbf {v}},R_{h}\psi)\in V_{h}\) defined by

$$\begin{aligned} \bigl((R_{h}v,v_{h})\bigr)=\bigl((v,v_{h}) \bigr) \quad \forall v\in W,v_{h}\in V_{h}, \end{aligned}$$

where

$$\begin{aligned} \bigl((u,v)\bigr)=n\nu(\nabla{\mathbf{u}},\nabla{\mathbf{v}})_{\Omega_{f}}+ \rho gK(\nabla \phi,\psi)_{\Omega_{p}}, \quad u=({\mathbf{u}}_{f}, \phi),v=({\mathbf{v}},\psi)\in W. \end{aligned}$$

Without loss of generality, we assume a uniform mesh applied to the time interval \([0,T]\) with \(t_{m}=m\Delta t, m=0,1,\ldots,J\), where \(\Delta t=\frac{T}{J}\) is the time step.

3.1 Coupled marching schemes for the mixed model

Recall that the mixed model (2.6) is formulated as an abstract time-dependent saddle-point problem. It is natural to consider the following first-order implicit marching scheme by applying the backward divided difference for the temporal discretization and the finite element Galerkin method for the spatial discretization, which leads to the coupled backward Euler scheme:

Algorithm 3.1

Coupled backward Euler scheme (CBES)

Find \(u^{m}_{h}=({\mathbf{u}}^{m}_{fh},\phi^{m}_{h})\in W_{h}\) and \(p_{fh}^{m}\in Q_{h}\) with \(m=1,\ldots,J\), such that for all \(v_{h}=({\mathbf{v}}_{h},\psi_{h})\in W_{h}\) and \(q_{h}\in Q_{h}\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} ( \frac{u_{h}^{m}-u_{h}^{m-1}}{\Delta t},v_{h})+a(u^{m}_{h},v_{h})+b(v_{h},p^{m}_{fh})=(f^{m},v_{h}),\\ b(u^{m}_{h},q_{h})=0,\\ u_{h}^{0}=R_{h}u^{0}, \end{array}\displaystyle \right . \end{aligned}$$
(3.2)

where \(f^{m}=f(t_{m})\) and \(u^{0}=({\mathbf{u}}_{f}^{0},\phi^{0})\). Note that at each time level, CBES amounts to solving a stationary Stokes-Darcy problem and is well-posed. Theoretical analysis and numerical experiments have been provided by Mu and his co-worker. In [1], authors not only provided the upper bounds for the numerical solution \((u_{h}^{m},p_{fh}^{m})\) of (3.2), but also established the corresponding optimal error estimates.

Theorem 3.1

see [1]

For CBES (3.2), we have

$$\begin{aligned} \bigl\Vert d_{t}u_{h}^{J}\bigr\Vert _{0}^{2}+\Delta t\sum_{m=1}^{J} \bigl\Vert d_{tt}u_{h}^{m}\bigr\Vert _{0}^{2}+\Delta t\sum_{m=1}^{J} \bigl\Vert \nabla d_{t}u_{h}^{m}\bigr\Vert _{0}^{2}\leq M_{0}. \end{aligned}$$
(3.3)

Furthermore, there exist constants \(C_{*}>0\) and \(C^{*}>0\) independent of h such that if

$$\begin{aligned} C_{*}h\leq\Delta t\leq C^{*}h, \end{aligned}$$
(3.4)

then

$$\begin{aligned} \bigl\Vert \nabla d_{t}u_{h}^{J} \bigr\Vert _{0}^{2}+\Delta t\sum _{m=2}^{J}\bigl\| d_{tt}u_{h}^{m} \bigr\| _{0}^{2}\leq M_{1}. \end{aligned}$$
(3.5)

Here and below, the positive constant \(M_{i} \) (\(i=0,1,\ldots\)) is independent of Δt and h.

In order to derive error estimates, we assume the regularity \(u\in(H^{2}(\Omega_{f}))^{d}\times H^{2}(\Omega_{p})\) and \(p\in H^{1}(\Omega_{f})\), and the finite element spaces as described above of first-order approximation \(\mathcal{O}(h)\) are used for the fluid and porous media regions. For convenience, from now on, we will use \(x\lesssim y\) to denote that there exists a positive constant C, such that \(x\leq Cy\).

Theorem 3.2

see [1]

For CBES (3.2) with \(m=1,\ldots,J\), we have

$$\begin{aligned}& \bigl\Vert u(t_{m})-u^{m}_{h}\bigr\Vert _{0}\lesssim\Delta t+h^{2}, \\& \bigl\Vert \nabla\bigl(u(t_{m})-u^{m}_{h}\bigr) \bigr\Vert _{0}\lesssim\Delta t+h, \\& \bigl\Vert p_{f}(t_{m})-p^{m}_{fh} \bigr\Vert _{0}\lesssim\Delta t+h+\Delta t^{-1}h^{2}. \end{aligned}$$

From CBES (3.2), we know that the variables \({\mathbf {u}}_{fh}, p_{fh}\) and \(\phi_{h}\) are coupled together by the boundary condition. When we solve this coupled system directly, the numerical difficulties increase as the mesh size decreases. In order to solve the non-stationary Stokes-Darcy model efficiently, some decoupled algorithms will be developed in the next subsection.

3.2 Decoupled marching schemes for the mixed model

Firstly, we recall a decoupled approach based on the temporal extrapolation on the interface, which has been researched in [1].

Algorithm 3.2

Decoupled backward Euler scheme 1 (DBES1)

Find \(u^{m,h}_{3.2}=({\mathbf{u}}^{m,h}_{3.2},\phi^{m,h}_{3.2})\in W_{h}\) and \(p_{3.2}^{m,h}\in Q_{h}\) with \(m=1,\ldots,J\), such that for all \(v_{h}=({\mathbf{v}}_{h},\psi_{h})\in W_{h}\) and \(q_{h}\in Q_{h}\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} ( \frac{u_{3.2}^{m,h}-u_{3.2}^{m-1,h}}{\Delta t},v_{h})+a_{\Omega }(u^{m,h}_{3.2},v_{h})+b(v_{h},p^{m,h}_{3.2})=(f^{m},v_{h})-a_{\Gamma }(u^{m-1,h}_{3.2},v_{h}),\\ b(u^{m,h}_{3.2},q_{h})=0,\\ u_{h}^{0}=R_{h}u^{0}. \end{array}\displaystyle \right . \end{aligned}$$
(3.6)

From the coercivity of \(a_{\Omega}(\cdot,\cdot)\) and \(b(\cdot,\cdot)\) satisfies the discrete inf-sup condition, we can see that the DBES1 is well-posed. Furthermore, at each time step, the discrete model (3.6) is equivalent to two decoupled problems that correspond to a Stokes problem in \(\Omega_{f}\) and a Darcy problem in \(\Omega_{p}\), respectively, with associated boundary conditions defined by \(u_{3.2}^{m-1,h}\) from the previous time level of Γ. In order to simplify the expression, we denote \(e_{3.2}^{m}=({\mathbf {e}}_{3.2}^{m},\xi_{3.2}^{m})\) with \({\mathbf{e}}_{3.2}^{m}={\mathbf{u}}_{fh}^{m}-{\mathbf{u}}_{3.2}^{m,h}\) and \(\xi _{3.2}^{m}=\phi_{h}^{m}-\phi_{3.2}^{m,h}\). In particular, \(e_{3.2}^{0}=({\mathbf{0}},0)\).

Theorem 3.3

see [1]

Under the condition of (3.4), for DBES1 (3.6) with \(m=1,\ldots,J\), we have

$$\begin{aligned}& \bigl\| e_{3.2}^{m}\bigr\| ^{2}_{0}+\Delta t\sum _{j=1}^{m}\bigl\| \nabla e_{3.2}^{j} \bigr\| _{0}^{2}\lesssim\Delta t^{2}, \end{aligned}$$
(3.7)
$$\begin{aligned}& \bigl\| \nabla e_{3.2}^{m}\bigr\| ^{2}_{0}+2\Delta t \sum_{j=1}^{m}\bigl\| d_{t}e_{3.2}^{j} \bigr\| _{0}^{2}\lesssim\Delta t. \end{aligned}$$
(3.8)

Combining Theorems 3.2 and 3.3, we can obtain the optimal order error estimates for the decoupled numerical solution \({\mathbf{u}}_{3.2}^{m,h}\) and \(\phi_{3.2}^{m,h}\) in \(L^{2}\) norm. But the \(H^{1}\) norm for \({\mathbf{u}}_{3.2}^{m,h}\) and \(\phi_{3.2}^{m,h}\) are suboptimal, namely, the order is \(\mathcal{O}(\Delta t^{\frac{1}{2}})\). This estimate may be improved to \(\mathcal{O}(\Delta t)\), as suggested by numerical experiments in [1]. It motivates us to propose some novel decoupled algorithms and derive the optimal error estimates. Our algorithms can be described as follows:

Algorithm 3.3

Decoupled backward Euler scheme 2 (DBES2)

Step 1. The discrete Stokes problem in the fluid region \(\Omega_{f}\) reads as follows: Find \({\mathbf{u}}^{m,h}_{3.3}\in H_{fh}, p_{3.3}^{m,h}\in Q_{h}\) with \(m=1,\ldots,J\), such that for all \({\mathbf{v}}_{h}\in H_{fh}\) and \(q_{h}\in Q_{h}\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} n( \frac{{\mathbf{u}}^{m,h}_{3.3}-{\mathbf {u}}^{m-1,h}_{3.3}}{\Delta t},{\mathbf{v}}_{h}) +a_{\Omega_{f}}({\mathbf{u}}^{m,h}_{3.3},{\mathbf{v}}_{h})+b_{\Omega_{f}}({\mathbf {v}}_{h},p_{3.3}^{m,h})\\ \quad=(n{\mathbf{g}}^{m}_{f},{\mathbf{v}}_{h})-\int_{\Gamma}n\rho g\phi_{3.3}^{m-1,h}{\mathbf {v}}_{h}\cdot{\mathbf{n}}_{f} ,\\ b_{\Omega_{f}}({\mathbf{u}}^{m,h}_{3.3},q_{h})=0,\\ {\mathbf{u}}^{0}_{fh}=R_{h}{\mathbf{u}}^{0}_{f}. \end{array}\displaystyle \right . \end{aligned}$$

Step 2. The discrete Darcy problem on the porous media region reads as follows: Find \(\phi^{m,h}_{3.3} \in H_{ph}\) with \(m=1,\ldots,J\), such that for all \(\psi_{h}\in H_{ph}\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} \rho gS_{0}( \frac{\phi_{3.3}^{m,h}-\phi _{h}^{m-1,h}}{\Delta t},\psi_{h})+ a_{\Omega_{p}}(\phi^{m,h}_{3.3},\psi_{h})=(\rho gg^{m}_{p},\psi_{h})+\int_{\Gamma }n\rho g\psi_{h}{\mathbf{u}}^{m,h}_{3.3}\cdot{\mathbf{n}}_{f},\\ \phi^{0}_{h}=R_{h}\phi^{0}. \end{array}\displaystyle \right . \end{aligned}$$

The Steps 1 and 2 in Algorithm 3.3 can be rewritten as: Find \(u^{m,h}_{3.3}=({\mathbf{u}}^{m,h}_{3.3},\phi^{m,h}_{3.3})\in W_{h}\) and \(p_{3.3}^{m,h}\in Q_{h}\) with \(m=1,\ldots,J\), such that for all \(v_{h}=({\mathbf{v}}_{h},\psi_{h})\in W_{h}\) and \(\forall q_{h}\in Q_{h}\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} ( \frac{u_{3.3}^{m,h}-u_{3.3}^{m-1,h}}{\Delta t},v_{h})+a_{\Omega}(u^{m,h}_{3.3},v_{h})+b(v_{h},p^{m,h}_{3.3})\\ \quad =(f^{m},v_{h})-a_{\Gamma}(u^{m-1,h}_{3.3},v_{h}) +\int_{\Gamma}n\rho g\psi_{h}({\mathbf{u}}^{m,h}_{3.3}-{\mathbf {u}}^{m-1,h}_{3.3})\cdot{\mathbf{n}}_{f},\\ b(u^{m,h}_{3.3},q_{h})=0,\\ u_{h}^{0}=R_{h}u^{0}. \end{array}\displaystyle \right . \end{aligned}$$
(3.9)

Algorithm 3.4

Decoupled backward Euler scheme 3 (DBES3)

Step 1. The discrete Darcy problem on the porous media region reads as follows: Find \(\phi^{m,h}_{3.4} \in H_{ph}\) with \(m=1,\ldots,J\), such that for all \(\psi_{h}\in H_{ph}\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} \rho gS_{0}( \frac{\phi_{3.4}^{m,h}-\phi _{3.4}^{m-1,h}}{\Delta t},\psi_{h})+ a_{\Omega_{p}}(\phi^{m,h}_{3.4},\psi_{h})=(\rho gg^{m}_{p},\psi_{h})+\int_{\Gamma }n\rho g\psi_{h}{\mathbf{u}}^{m-1,h}_{3.4}\cdot{\mathbf{n}}_{f},\\ \phi^{0}_{h}=R_{h}\phi^{0}. \end{array}\displaystyle \right . \end{aligned}$$

Step 2. The discrete Stokes problem in the fluid region \(\Omega_{f}\) reads as follows: Find \({\mathbf{u}}^{m,h}_{3.4}\in H_{fh}\) and \(p_{3.4}^{m,h}\in Q_{h}\) with \(m=1,\ldots,J\), such that for all \({\mathbf{v}}_{h}\in H_{fh}\) and \(q_{h}\in Q_{h}\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} n( \frac{{\mathbf{u}}^{m,h}_{3.4}-{\mathbf {u}}^{m-1,h}_{3.4}}{\Delta t},{\mathbf{v}}_{h}) +a_{\Omega_{f}}({\mathbf{u}}^{m,h}_{3.4},{\mathbf{v}}_{h})+b_{\Omega_{f}}({\mathbf {v}}_{h},p_{3.4}^{m,h})\\ \quad=(n{\mathbf{g}}^{m}_{f},{\mathbf{v}}_{h})-\int_{\Gamma}n\rho g\phi_{3.4}^{m,h}{\mathbf {v}}_{h}\cdot{\mathbf{n}}_{f} ,\\ b_{\Omega_{f}}({\mathbf{u}}^{m,h}_{3.4},q_{h})=0,\\ {\mathbf{u}}^{0}_{fh}=R_{h}{\mathbf{u}}^{0}_{f}. \end{array}\displaystyle \right . \end{aligned}$$

The Steps 1 and 2 in Algorithm 3.4 can be rewritten as: Find \(u^{m,h}_{3.4}=({\mathbf{u}}^{m,h}_{3.4},\phi^{m,h}_{3.4})\in W_{h}\) and \(p_{3.4}^{m,h}\in Q_{h}\) with \(m=1,\ldots,J\), such that for all \(v_{h}=({\mathbf{v}}_{h},\psi_{h})\in W_{h}\) and \(q_{h}\in Q_{h}\)

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} ( \frac{u_{3.4}^{m,h}-u_{3.4}^{m-1,h}}{\Delta t},v_{h})+a_{\Omega}(u^{m,h}_{3.4},v_{h})+b(v_{h},p^{m,h}_{3.4})\\ \quad =(f^{m},v_{h})-a_{\Gamma }(u^{m-1,h}_{3.4},v_{h}) -\int_{\Gamma}n\rho g(\phi_{3.4}^{m,h}-\phi_{3.4}^{m-1,h}){\mathbf {v}}_{h}\cdot{\mathbf{n}}_{f} ,\\ b(u^{m,h}_{3.4},q_{h})=0,\\ u_{h}^{0}=R_{h}u^{0}. \end{array}\displaystyle \right . \end{aligned}$$
(3.10)

Similar to the DBES1, we can see that both Algorithms 3.3 and 3.4 are well-posed. In scheme (3.6), we use the numerical solution \(u^{m-1,h}_{h}\) from the previous time level to approximate the interface conditions. One advantage of algorithm (3.6) is that it can be used in parallelism based on the solution of previous time level. In order to improve the computational accuracy, we separate the coupled model (2.6) into two steps (one Stokes equation in \(\Omega_{f}\) and one Darcy problem in \(\Omega _{p}\)), and use the numerical solution obtained in step 1 to approximate the boundary condition of step 2 at the same time level.

4 Stability

This section is devoted to establishing the upper bounds for the solutions \(u_{3.3}^{m,h}\) and \(u_{3.4}^{m,h}\) of decoupled Algorithms 3.3 and 3.4, respectively, which will be used in the error estimates for \(e^{m}_{3.3}=({\mathbf{u}}_{fh}^{m}-{\mathbf {u}}_{3.3}^{m,h},\phi_{h}^{m}-\phi_{3.3}^{m,h})\) and \(e_{3.4}^{m}=({\mathbf{u}}_{fh}^{m}-{\mathbf{u}}_{3.4}^{m,h},\phi_{h}^{m}-\phi _{3.4}^{m,h})\) in Section 5. For convenience, let us introduce the following notations. We denote the backward divided difference operator \(d_{t}\) by

$$\begin{aligned} d_{t}u_{h}^{m}=\frac{u_{h}^{m}-u_{h}^{m-1}}{\Delta t}, \quad \mbox{for } m=1,2\ldots,J. \end{aligned}$$

When \(m=0\), we define \(d_{t}u_{h}^{0}=(d_{t}{\mathbf{u}}^{0}_{fh},d_{t}\phi_{h}^{0})\) as the solution to the following problem:

$$\begin{aligned} \bigl(d_{t}u_{h}^{0},v_{h} \bigr)+a\bigl(u_{h}^{0},v_{h}\bigr)= \bigl(f^{0},v_{h}\bigr) \quad \forall v_{h}\in V_{h}. \end{aligned}$$
(4.1)

We also denote

$$\begin{aligned} d_{tt}u_{h}^{m}=\frac{d_{t}u_{h}^{m}-d_{t}u_{h}^{m-1}}{\Delta t} \quad\mbox{for } m=1,2,\ldots,J. \end{aligned}$$

Similarly, we introduce the divide differences \(d_{t}u_{3.i}^{m,h}\) and \(d_{tt}u_{3.i}^{m,h}\) for \(m=1,\ldots,J\) and \(i=2,3,4\) for the solutions \(u_{3.i}^{m,h}\) of Algorithms 3.2, 3.3 and 3.4, respectively. For \(m=0, d_{t}u_{3.i}^{0}\) is defined in the same way as \(d_{t}u_{h}^{0}\) in (4.1).

Lemma 4.1

For DBES2 (3.9), under the condition (3.4) we have

$$\begin{aligned} \bigl\| d_{t}u_{3.3}^{J,h}\bigr\| _{0}^{2}+ \frac{1}{4}\Delta t^{2}\sum_{m=1}^{J} \bigl\| d_{tt}u_{3.3}^{m,h}\bigr\| _{0}^{2}+ \frac{3}{4}\Delta t\sum_{m=1}^{J}\bigl\| d_{t}u_{3.3}^{m,h}\bigr\| _{W}^{2} \leq M_{2}. \end{aligned}$$

Proof

We subtract the decoupled backward Euler scheme (3.9) on two adjacent time levels and notice the definition of \(d_{t}u_{3.3}^{0}\), for all \(v_{h}\in V_{h}\) and \(m=1,\ldots,J\) we have

$$\begin{aligned} &\biggl(\frac{d_{t}u_{3.3}^{m,h}-d_{t}u_{3.3}^{m-1,h}}{\Delta t},v_{h}\biggr)+a_{\Omega } \bigl(d_{t}u_{3.3}^{m,h},v_{h}\bigr) \\ &\quad=\frac{1}{\Delta t}\bigl(f^{m}-f^{m-1},v_{h} \bigr)-a_{\Gamma}\bigl(d_{t}u_{3.3}^{m-1,h},v_{h} \bigr) +\Delta t^{2} \int_{\Gamma}n\rho g\psi_{h}d_{tt}{ \mathbf{u}}_{3.3}^{m,h}\cdot{\mathbf{n}}_{f}. \end{aligned}$$
(4.2)

Taking \(v_{h}=2\Delta td_{t}u_{3.3}^{m,h}=2\Delta t(d_{t}{\mathbf {u}}_{3.3}^{m,h},d_{t}\phi_{3.3}^{m,h})\in V_{h}\) in (4.2) one finds

$$\begin{aligned} &2\Delta t\biggl(\frac{d_{t}u_{3.3}^{m,h}-d_{t}u_{3.3}^{m-1,h}}{\Delta t},d_{t}u_{3.3}^{m,h} \biggr)+2\Delta ta_{\Omega}\bigl(d_{t}u_{3.3}^{m,h},d_{t}u_{3.3}^{m,h} \bigr) \\ &\quad=\frac{1}{\Delta t}\bigl(f^{m}-f^{m-1},2\Delta td_{t}u_{3.3}^{m,h}\bigr) -2\Delta ta_{\Gamma} \bigl(d_{t}u_{3.3}^{m-1,h},d_{t}u_{3.3}^{m,h} \bigr) \\ &\qquad{}+2\Delta t^{3} \int_{\Gamma}n\rho g\cdot d_{tt}{\mathbf{u}}_{3.3}^{m,h} \cdot d_{t}\phi_{3.3}^{m,h}\cdot{ \mathbf{n}}_{f}. \end{aligned}$$

Then by using the equality \((a-b,2a)=a^{2}-b^{2}+(a-b)^{2}\) and (2.7) we have

$$\begin{aligned} &\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| _{0}^{2}-\bigl\| d_{t}u_{3.3}^{m-1,h} \bigr\| _{0}^{2}+\Delta t^{2}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| _{0}^{2}+2\Delta t\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| _{W}^{2} \\ &\quad\leq2n \int_{\Omega_{f}}\bigl({\mathbf{g}}_{f}^{m}-{ \mathbf{g}}_{f}^{m-1}\bigr)\cdot d_{t}{\mathbf {u}}_{3.3}^{m,h} +2\rho g \int_{\Omega_{p}}\bigl(g_{p}^{m}-g_{p}^{m-1} \bigr)\cdot d_{t}\phi_{3.3}^{m,h} \\ &\qquad{}-2\Delta ta_{\Gamma}\bigl(d_{t}u_{3.3}^{m-1,h},d_{t}u_{3.3}^{m,h} \bigr) +2\Delta t^{3} \int_{\Gamma}n\rho g\cdot d_{tt}{\mathbf{u}}_{3.3}^{m,h} \cdot d_{t}\phi_{3.3}^{m,h}\cdot{ \mathbf{n}}_{f} \\ &\quad\leq Cn\nu^{-1} \int_{t_{m-1}}^{t_{m}}\|{\mathbf{g}}_{ft} \|_{0}^{2}\,dt +C\rho gK^{-1} \int_{t_{m-1}}^{t_{m}}\bigl\| g_{pt}^{m} \bigr\| _{0}^{2}\,dt +n\nu\Delta t\bigl\| \nabla d_{t}{ \mathbf{u}}_{3.3}^{m,h}\bigr\| _{0}^{2} \\ &\qquad{}+\rho gK\Delta t\bigl\| \nabla d_{t}\phi_{3.3}^{m,h} \bigr\| _{0}^{2} +2\Delta ta_{\Gamma}\bigl(d_{t}u_{3.3}^{m,h}-d_{t}u_{3.3}^{m-1,h},d_{t}u_{3.3}^{m,h} \bigr) \\ &\qquad{}+2\Delta t^{3} \int_{\Gamma}n\rho g\cdot d_{tt}{\mathbf{u}}_{3.3}^{m,h} \cdot d_{t}\phi_{3.3}^{m,h}\cdot{ \mathbf{n}}_{f}. \end{aligned}$$
(4.3)

For the last two terms in right-hand side (4.3), thanks to the trace theorem, the inverse and Cauchy inequalities, and (3.4), we can treat them as follows:

$$\begin{aligned} &2\Delta ta_{\Gamma}\bigl(d_{t}u_{3.3}^{m,h}-d_{t}u_{3.3}^{m-1,h},d_{t}u_{3.3}^{m,h} \bigr) \\ &\quad=2\Delta t^{2}a_{\Gamma}\bigl(d_{tt}u_{3.3}^{m,h},d_{t}u_{3.3}^{m,h} \bigr) \\ &\quad\leq2\Delta t^{2}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| _{L^{2}(\Gamma)}\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| _{L^{2}(\Gamma)} \\ &\quad\leq2C_{t}^{2}\Delta t^{\frac{7}{4}} \bigl\| d_{tt}u_{3.3}^{m,h}\bigr\| ^{1/2}_{0}\bigl\| d_{tt}u_{3.3}^{m,h}\bigr\| ^{1/2}_{W} \bigl\| d_{t}u_{3.3}^{m,h}\bigr\| ^{1/2}_{0} \cdot\Delta t^{\frac{1}{4}}\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| ^{1/2}_{W} \\ &\quad\leq\frac{1}{4}\Delta t\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| ^{2}_{W}+\frac {3}{4}\bigl(2C_{t}^{2} \bigr)^{\frac{4}{3}}\Delta t^{\frac{7}{3}}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| ^{2/3}_{0}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| ^{2/3}_{W}\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| ^{2/3}_{0} \\ &\quad\leq\frac{1}{4}\Delta t\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| ^{2}_{W}+\frac {3}{4}\bigl(2C_{t}^{2} \bigr)^{\frac{4}{3}}\Delta t^{2}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| ^{2/3}_{0}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| ^{2/3}_{W}\cdot\Delta t^{\frac{1}{3}}\bigl\| d_{t}u_{3.3}^{m,h}\bigr\| ^{2/3}_{0} \\ &\quad\leq\frac{1}{4}\Delta t\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| ^{2}_{W}+\frac {2}{3C_{in}}\biggl(\frac{3}{4} \biggr)^{\frac{3}{2}}\Delta t^{3}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| _{0}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| _{W} \\ &\qquad{}+\frac{2^{4}C_{t}^{8}C_{in}}{3}\Delta t\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| _{0}^{2} \\ &\quad\leq\frac{1}{4}\Delta t\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| ^{2}_{W}+\frac{1}{2}\Delta t^{2} \bigl\| d_{tt}u_{3.3}^{m,h}\bigr\| ^{2}_{0} + \frac{2^{4}C_{t}^{8}C_{in}}{3}\Delta t\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| _{0}^{2}, \end{aligned}$$
(4.4)
$$\begin{aligned} &2\Delta t^{3} \int_{\Gamma}n\rho g\cdot d_{tt}{\mathbf{u}}_{3.3}^{m,h} \cdot d_{t}\phi_{3.3}^{m,h}\cdot{ \mathbf{n}}_{f} \\ &\quad\leq2\Delta t^{3}n\rho g\bigl\| d_{tt}{\mathbf{u}}_{3.3}^{m,h} \bigr\| _{L^{2}(\Gamma)}\bigl\| d_{t}\phi_{3.3}^{m,h} \bigr\| _{L^{2}(\Gamma)} \\ &\quad\leq2C_{t}^{2}\Delta t^{3}n\rho g \bigl\| d_{tt}{\mathbf{u}}_{3.3}^{m,h}\bigr\| ^{1/2}_{L^{2}(\Omega_{f})}\bigl\| d_{tt}{\mathbf{u}}_{3.3}^{m,h} \bigr\| ^{1/2}_{H^{1}(\Omega_{f})} \bigl\| d_{t}\phi_{3.3}^{m,h} \bigr\| ^{1/2}_{L^{2}(\Omega_{p})}\bigl\| d_{t}\phi_{3.3}^{m,h} \bigr\| ^{1/2}_{H^{1}(\Omega_{p})} \\ &\quad\leq2C_{t}^{2}C_{in}\Delta t^{2}n \rho g\bigl\| d_{tt}{\mathbf{u}}_{3.3}^{m,h}\bigr\| _{L^{2}(\Omega_{f})} \bigl\| d_{t}\phi_{3.3}^{m,h} \bigr\| _{L^{2}(\Omega_{p})} \\ &\quad\leq\frac{1}{4}n\Delta t^{2}\bigl\| d_{tt}{ \mathbf{u}}_{3.3}^{m,h}\bigr\| ^{2}_{L^{2}(\Omega _{f})}+ \frac{n\rho g\Delta tC_{t}^{4}C_{in}^{2}}{S_{0}}\rho gS_{0}\Delta t \bigl\| d_{t} \phi_{3.3}^{m,h}\bigr\| ^{2}_{L^{2}(\Omega_{p})}. \end{aligned}$$
(4.5)

Combining (4.4)-(4.5) with (4.3) and summing m form 1 to J we have

$$\begin{aligned} &\bigl\| d_{t}u_{3.3}^{J,h} \bigr\| _{0}^{2}+\frac{1}{4}\Delta t^{2}\sum _{m=1}^{J}\bigl\| d_{tt}u_{3.3}^{m,h} \bigr\| _{0}^{2}+\frac{3}{4}\Delta t\sum _{m=1}^{J}\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| _{W}^{2} \\ &\quad\leq Cn\nu^{-1} \int_{0}^{T}\|{\mathbf{g}}_{ft} \|_{0}^{2}\,dt +C\rho gK^{-1} \int_{0}^{T}\bigl\| g_{pt}^{m} \bigr\| _{0}^{2}\,dt+\bigl\| d_{t}u_{3.3}^{0} \bigr\| _{0}^{2} \\ &\qquad{}+\biggl(\frac{2^{4}C_{t}^{8}C_{in}}{3}+\frac{n\rho g\Delta tC_{t}^{4}C_{in}^{2}}{S_{0}}\biggr)\Delta t\sum _{m=1}^{J}\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| _{0}^{2}. \end{aligned}$$

Applying the Gronwall lemma, we obtain the desired result. □

Lemma 4.2

Under the condition of (3.4), for the decoupled Algorithm  3.4 with \(m=1,2,\ldots,J\) we have

$$\begin{aligned} \bigl\Vert d_{t}u_{3.4}^{J,h}\bigr\Vert _{0}^{2}+\frac{1}{4}\Delta t^{2}\sum _{m=1}^{J}\bigl\Vert d_{tt}u_{3.4}^{m,h} \bigr\Vert _{0}^{2}+\frac{3}{4}\Delta t\sum _{m=1}^{J}\bigl\Vert d_{t}u_{3.4}^{m,h} \bigr\Vert _{W}^{2} \leq M_{3}. \end{aligned}$$

Proof

For all \(v_{h}\in V_{h}\), the following error equation can be obtained by using (3.10) on two adjacent time levels

$$\begin{aligned} &\biggl(\frac{d_{t}u_{3.4}^{m,h}-d_{t}u_{3.4}^{m-1,h}}{\Delta t},v_{h}\biggr)+a_{\Omega } \bigl(d_{t}u_{3.4}^{m,h},v_{h}\bigr) \\ &\quad=\frac{1}{\Delta t}\bigl(f^{m}-f^{m-1},v_{h} \bigr)-a_{\Gamma }\bigl(d_{t}u_{3.4}^{m-1,h},v_{h} \bigr) -\Delta t^{2} \int_{\Gamma}n\rho gd_{tt}\phi_{3.4}^{m,h}{ \mathbf{v}}_{h}\cdot{\mathbf{n}}_{f}. \end{aligned}$$
(4.6)

Taking \(v_{h}=2\Delta td_{t}u_{3.4}^{m,h}=2\Delta t(d_{t}{\mathbf {u}}_{3.4}^{m,h},d_{t}\phi_{3.4}^{m,h})\in V_{h}\) in (4.6) we obtain

$$\begin{aligned} &\biggl(\frac{d_{t}u_{3.4}^{m,h}-d_{t}u_{3.4}^{m-1,h}}{\Delta t},2\Delta td_{t}u_{3.4}^{m,h} \biggr)+2\Delta ta_{\Omega}\bigl(d_{t}u_{3.4}^{m,h},d_{t}u_{3.4}^{m,h} \bigr) \\ &\quad=\frac{1}{\Delta t}\bigl(f^{m}-f^{m-1},2\Delta td_{t}u_{3.4}^{m,h}\bigr) -2\Delta ta_{\Gamma} \bigl(d_{t}u_{3.4}^{m-1,h},d_{t}u_{3.4}^{m,h} \bigr) \\ &\qquad{}-2\Delta t^{3} \int_{\Gamma}n\rho g d_{tt}\phi_{3.4}^{m,h}d_{t}{ \mathbf {u}}_{3.4}^{m,h}\cdot{\mathbf{n}}_{f}. \end{aligned}$$

As a consequence we have

$$\begin{aligned} &\bigl\| d_{t}u_{3.4}^{m,h} \bigr\| _{0}^{2}-\bigl\| d_{t}u_{3.4}^{m-1,h} \bigr\| _{0}^{2}+\Delta t^{2}\bigl\| d_{tt}u_{3.4}^{m,h} \bigr\| _{0}^{2}+2\Delta t\bigl\| d_{t}u_{3.4}^{m,h} \bigr\| _{W}^{2} \\ &\quad\leq2n \int_{\Omega_{f}}\bigl({\mathbf{g}}_{f}^{m}-{ \mathbf{g}}_{f}^{m-1}\bigr)\cdot d_{t}{\mathbf {u}}_{3.4}^{m,h} +2\rho g \int_{\Omega_{p}}\bigl(g_{p}^{m}-g_{p}^{m-1} \bigr)\cdot d_{t}\phi_{3.4}^{m,h} \\ &\qquad{}-2\Delta ta_{\Gamma}\bigl(d_{t}u_{3.4}^{m-1,h},d_{t}u_{3.4}^{m,h} \bigr) -2\Delta t^{3} \int_{\Gamma}n\rho g\cdot d_{t}{\mathbf{u}}_{3.4}^{m,h} \cdot d_{tt}\phi_{3.4}^{m,h}\cdot{ \mathbf{n}}_{f} \\ &\quad\leq Cn\nu^{-1} \int_{t_{m-1}}^{t_{m}}\|{\mathbf{g}}_{ft} \|_{0}^{2}\,dt +C\rho gK^{-1} \int_{t_{m-1}}^{t_{m}}\bigl\| g_{pt}^{m} \bigr\| _{0}^{2}\,dt +n\nu\Delta t\bigl\| \nabla d_{t}{ \mathbf{u}}_{3.4}^{m,h}\bigr\| _{0}^{2} \\ &\qquad{}+\rho gK\Delta t\bigl\| \nabla d_{t}\phi_{3.4}^{m,h} \bigr\| _{0}^{2} +2\Delta ta_{\Gamma}\bigl(d_{t}u_{3.4}^{m,h}-d_{t}u_{3.4}^{m-1,h},d_{t}u_{3.4}^{m,h} \bigr) \\ &\qquad{}-2\Delta t^{3} \int_{\Gamma}n\rho g\cdot d_{t}{\mathbf{u}}_{3.4}^{m,h} \cdot d_{tt}\phi_{3.4}^{m,h}\cdot{ \mathbf{n}}_{f}. \end{aligned}$$
(4.7)

For the last term in (4.7), using trace and inverse inequalities and (3.4) yields

$$\begin{aligned} &2\Delta t^{3} \int_{\Gamma}n\rho g\cdot d_{t}{\mathbf{u}}_{3.4}^{m,h} \cdot d_{tt}\phi_{3.4}^{m,h}\cdot{ \mathbf{n}}_{f} \\ &\quad\leq2\Delta t^{3}n\rho g\bigl\| d_{t}{\mathbf{u}}_{3.4}^{m,h} \bigr\| _{L^{2}(\Gamma)}\bigl\| d_{tt}\phi_{3.4}^{m,h} \bigr\| _{L^{2}(\Gamma)} \\ &\quad\leq2C_{t}^{2}\Delta t^{3}n\rho g \bigl\| d_{t}{\mathbf{u}}_{3.4}^{m,h}\bigr\| ^{1/2}_{L^{2}(\Omega_{f})}\bigl\| d_{t}{\mathbf{u}}_{3.4}^{m,h} \bigr\| ^{1/2}_{H^{1}(\Omega_{f})} \bigl\| d_{tt}\phi_{3.4}^{m,h} \bigr\| ^{1/2}_{L^{2}(\Omega_{p})}\bigl\| d_{tt}\phi _{3.4}^{m,h} \bigr\| ^{1/2}_{H^{1}(\Omega_{p})} \\ &\quad\leq2C^{2}_{t}C_{in}\Delta t^{2}n \rho g\bigl\| d_{t}{\mathbf{u}}_{3.4}^{m,h}\bigr\| _{L^{2}(\Omega_{f})} \bigl\| d_{tt}\phi_{3.4}^{m,h} \bigr\| _{L^{2}(\Omega_{p})} \\ &\quad\leq\frac{1}{4}\rho gS_{0}\Delta t^{2} \bigl\| d_{tt}\phi_{3.4}^{m,h}\bigr\| ^{2}_{L^{2}(\Omega_{p})}+ \frac{4n\rho gC_{t}^{4}C_{in}^{2}\Delta t}{S_{0}}n\Delta t \bigl\| d_{t}{\mathbf{u}}_{3.4}^{m,h} \bigr\| ^{2}_{L^{2}(\Omega_{f})}. \end{aligned}$$
(4.8)

Combining (4.4), (4.8) with (4.7), summing m from 1 to J, and applying the Gronwall lemma, we complete the proof. □

5 Convergence analysis

In this section we present the error estimates for the decoupled Algorithms 3.3 and 3.4. As we have mentioned before, the suboptimal \(H^{1}\)-norm error estimates for \({\mathbf{u}}_{3.2}^{m,h}\) and \(\phi_{3.2}^{m,h}\) have been derived by Mu and Zhu in [1]. The estimate (3.8) can be improved to \(\mathcal{O}(\Delta t)\), as suggested by numerical experiments. Firstly, we improve the estimate (3.8) to the optimal order.

Lemma 5.1

Let \(({\mathbf{u}}_{fh}^{m},\phi_{h}^{m})\) and \(({\mathbf {u}}_{3.2}^{m,h},\phi_{3.2}^{m,h})\) be the solutions of the discrete models (3.2) and (3.6), respectively. Denote \(e_{3.2}^{m}=({\mathbf{e}}_{3.2}^{m},\xi_{3.2}^{m})=({\mathbf{u}}_{fh}^{m}-{\mathbf {u}}_{3.2}^{m,h},\phi_{h}^{m}-\phi_{3.2}^{m,h})\) we have

$$\begin{aligned} \bigl\Vert \nabla e_{3.2}^{m}\bigr\Vert ^{2}_{0}+\Delta t\sum_{j=1}^{m} \bigl\| d_{t}e_{3.2}^{j}\bigr\| _{0}^{2} \lesssim\Delta t^{2}. \end{aligned}$$

Proof

For all \(v_{h}\in V_{h}\) and \(e_{3.2}^{m}\) with \(m=1,\ldots,J\), the following error equation can be obtained by comparing the discrete models (3.2): and (3.6)

$$\begin{aligned} \biggl(\frac{e_{3.2}^{m}-e_{3.2}^{m-1}}{\Delta t},v_{h}\biggr)+a_{\Omega} \bigl(e_{3.2}^{m},v_{h}\bigr) +a_{\Gamma} \bigl(u_{h}^{m}-u^{m-1}_{3.2},v_{h} \bigr)=0. \end{aligned}$$
(5.1)

Taking \(v_{h}=2(e_{3.2}^{m}-e_{3.2}^{m-1})\in V_{h}\) in (5.1), we have

$$\begin{aligned} &\biggl(\frac{e_{3.2}^{m}-e_{3.2}^{m-1}}{\Delta t},2\bigl(e_{3.2}^{m}-e_{3.2}^{m-1} \bigr)\biggr)+a_{\Omega} \bigl(e_{3.2}^{m},2 \bigl(e_{3.2}^{m}-e_{3.2}^{m-1}\bigr)\bigr) \\ &\quad=-a_{\Gamma} \bigl(u_{h}^{m}-u^{m-1,h}_{3.2},2 \bigl(e_{3.2}^{m}-e_{3.2}^{m-1}\bigr)\bigr). \end{aligned}$$
(5.2)

For the right-hand side term of (5.2), we can estimate it as follows:

$$\begin{aligned} &a_{\Gamma}\bigl(u_{h}^{m}-u^{m-1,h}_{3.2},2 \bigl(e_{3.2}^{m}-e_{3.2}^{m-1}\bigr)\bigr) \\ &\quad=2\Delta ta_{\Gamma} \bigl(u_{h}^{m}-u_{h}^{m-1}+u_{h}^{m-1}-u^{m-1,h}_{3.2},d_{t}e_{3.2}^{m} \bigr) \\ &\quad=2\Delta t^{2}a_{\Gamma}\bigl(d_{t}u_{h}^{m},d_{t}e_{3.2}^{m} \bigr) +2\Delta ta_{\Gamma}\bigl(e_{3.2}^{m-1},d_{t}e_{3.2}^{m} \bigr) \\ &\quad\leq2\Delta t^{2}\bigl\| d_{t}u_{h}^{m} \bigr\| _{L^{2}(\Gamma)}\bigl\| d_{t}e_{3.2}^{m} \bigr\| _{L^{2}(\Gamma )}+2\Delta t\bigl\| e_{3.2}^{m-1}\bigr\| _{L^{2}(\Gamma)} \bigl\| d_{t}e_{3.2}^{m}\bigr\| _{L^{2}(\Gamma)} \\ &\quad\leq2C_{t}^{2}\Delta t^{2} \bigl\| d_{t}u_{h}^{m}\bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla d_{t}u_{h}^{m}\bigr\| ^{\frac{1}{2}}_{0} \bigl\| d_{t}e_{3.2}^{m}\bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{\frac{1}{2}}_{0} \\ &\qquad{}+2C_{t}^{2}\Delta t\bigl\| e_{3.2}^{m-1} \bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla e_{3.2}^{m-1} \bigr\| ^{\frac{1}{2}}_{0} \bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla d_{t}e_{3.2}^{m} \bigr\| ^{\frac{1}{2}}_{0}. \end{aligned}$$
(5.3)

For the first term of the right-hand side in (5.3), using the Young inequality, we obtain

$$\begin{aligned} &\Delta t^{2}\bigl\| d_{t}u_{h}^{m} \bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla d_{t}u_{h}^{m} \bigr\| ^{\frac{1}{2}}_{0} \bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla d_{t}e_{3.2}^{m} \bigr\| ^{\frac{1}{2}}_{0} \\ &\quad=\Delta t^{\frac{3}{2}}\bigl\| d_{t}u_{h}^{m} \bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla d_{t}u_{h}^{m} \bigr\| ^{\frac{1}{2}}_{0} \bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{\frac{1}{2}}_{0}\cdot\Delta t^{\frac{1}{2}}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{\frac{1}{2}}_{0} \\ &\quad\leq\frac{3}{4}\Delta t^{2}\bigl\| d_{t}u_{h}^{m} \bigr\| ^{\frac{2}{3}}_{0}\bigl\| \nabla d_{t}u_{h}^{m} \bigr\| ^{\frac{2}{3}}_{0} \bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{\frac{2}{3}}_{0}+\frac{1}{4}\Delta t^{2}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{2}_{0} \\ &\quad=\frac{3}{4}\Delta t^{\frac{5}{3}}\bigl\| d_{t}u_{h}^{m} \bigr\| ^{\frac{2}{3}}_{0}\bigl\| \nabla d_{t}u_{h}^{m} \bigr\| ^{\frac{2}{3}}_{0}\cdot \Delta t^{\frac{1}{3}}\bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{\frac{2}{3}}_{0}+\frac {1}{4}\Delta t^{2}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{2}_{0} \\ &\quad\leq\frac{2}{3}\biggl(\frac{3}{4}\biggr)^{\frac{3}{2}}\Delta t^{\frac {5}{2}}\bigl\| d_{t}u_{h}^{m} \bigr\| _{0}\bigl\| \nabla d_{t}u_{h}^{m} \bigr\| _{0}+ \frac{1}{3}\Delta t\bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{2}_{0}+\frac{1}{4}\Delta t^{2}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{2}_{0} \\ &\quad=\frac{2}{3}\biggl(\frac{3}{4}\biggr)^{\frac{3}{2}} \Delta t^{2}\bigl\| d_{t}u_{h}^{m} \bigr\| _{0}\cdot\Delta t^{\frac{1}{2}}\bigl\| \nabla d_{t}u_{h}^{m} \bigr\| _{0}+ \frac{1}{3}\Delta t\bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{2}_{0}+\frac{1}{4}\Delta t^{2}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{2}_{0}. \end{aligned}$$

For the second term of the right-hand side in (5.3), by the Young inequality we get

$$\begin{aligned} &\Delta t\bigl\| e_{3.2}^{m-1}\bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla e_{3.2}^{m-1}\bigr\| ^{\frac{1}{2}}_{0} \bigl\| d_{t}e_{3.2}^{m}\bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{\frac{1}{2}}_{0} \\ &\quad=\Delta t^{\frac{1}{2}}\bigl\| e_{3.2}^{m-1} \bigr\| ^{\frac{1}{2}}_{0}\bigl\| \nabla e_{3.2}^{m-1} \bigr\| ^{\frac{1}{2}}_{0} \bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{\frac{1}{2}}_{0}\cdot\Delta t^{\frac{1}{2}}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{\frac{1}{2}}_{0} \\ &\quad\leq\frac{3}{4}\Delta t^{\frac{2}{3}}\bigl\| e_{3.2}^{m-1} \bigr\| ^{\frac {2}{3}}_{0}\bigl\| \nabla e_{3.2}^{m-1} \bigr\| ^{\frac{2}{3}}_{0} \bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{\frac{2}{3}}_{0}+\frac{1}{4}\Delta t^{2}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{2}_{0} \\ &\quad=\frac{3}{4}\Delta t^{\frac{1}{3}}\bigl\| e_{3.2}^{m-1} \bigr\| ^{\frac{2}{3}}_{0}\bigl\| \nabla e_{3.2}^{m-1} \bigr\| ^{\frac{2}{3}}_{0} \cdot\Delta t^{\frac{1}{3}}\bigl\| d_{t}e_{3.2}^{m} \bigr\| ^{\frac{2}{3}}_{0}+\frac {1}{4}\Delta t^{2}\bigl\| \nabla d_{t}e_{3.2}^{m}\bigr\| ^{2}_{0} \\ &\quad\leq\frac{2}{3}\biggl(\frac{3}{4}\biggr)^{\frac{3}{2}} \Delta t^{\frac {1}{2}}\bigl\| e_{3.2}^{m-1}\bigr\| _{0}\bigl\| \nabla e_{3.2}^{m-1}\bigr\| _{0} +\frac{1}{3}\Delta t \bigl\| d_{t}e_{3.2}^{m}\bigr\| ^{2}_{0}+ \frac{1}{4}\Delta t^{2}\bigl\| \nabla d_{t}e_{3.2}^{m} \bigr\| ^{2}_{0}. \end{aligned}$$

On the other hand, for the left-hand side of (5.2), using the trick of [1] we have

$$\begin{aligned} &\biggl( \frac{e_{3.2}^{m}-e_{3.2}^{m-1}}{\Delta t},2\bigl(e_{3.2}^{m}-e_{3.2}^{m-1} \bigr)\biggr)+a_{\Omega} \bigl(e_{3.2}^{m},2 \bigl(e_{3.2}^{m}-e_{3.2}^{m-1}\bigr)\bigr) \\ &\quad=2\Delta t\bigl(d_{t}e_{3.2}^{m},d_{t}e_{3.2}^{m} \bigr)+a_{\Omega }\bigl(e_{3.2}^{m},e_{3.2}^{m} \bigr)-a_{\Omega}\bigl(e_{3.2}^{m-1},e_{3.2}^{m-1} \bigr)+\Delta t^{2}a_{\Omega}\bigl(d_{t}e_{3.2}^{m},d_{t}e_{3.2}^{m} \bigr) \\ &\quad=2\Delta t\bigl\| d_{t}e_{3.2}^{m} \bigr\| _{0}^{2}+\bigl\| \nabla e_{3.2}^{m} \bigr\| _{0}^{2} \\ &\qquad{}+n\sum_{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.2}^{m} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2}-\bigl\| \nabla e_{3.2}^{m-1}\bigr\| _{0}^{2} +\Delta t^{2}\bigl\| \nabla d_{t}e_{3.2}^{m} \bigr\| _{0}^{2} \\ &\qquad{}-n\sum_{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.2}^{m-1} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2} +n\Delta t^{2} \sum_{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl(d_{t}{ \mathbf{e}}_{3.2}^{m-1}\cdot{\boldsymbol {\tau}}_{i} \bigr)^{2}. \end{aligned}$$
(5.4)

Combining (5.3), (5.4) with (5.2), we obtain

$$\begin{aligned} &\Delta t\bigl\| d_{t}e_{3.2}^{m} \bigr\| _{0}^{2}+\bigl\| \nabla e_{3.2}^{m} \bigr\| _{0}^{2}-\bigl\| \nabla e_{3.2}^{m-1} \bigr\| _{0}^{2} \\ &\qquad{}+n\sum_{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.2}^{m} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2} -n\sum _{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.2}^{m-1} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2} \\ &\quad\leq\frac{2}{3}\biggl(\frac{3}{4}\biggr)^{\frac{3}{2}} \bigl[\Delta t^{2}\bigl\| d_{t}u_{h}^{m} \bigr\| _{0}\cdot\Delta t^{\frac{1}{2}}\bigl\| \nabla d_{t}u_{h}^{m} \bigr\| _{0}+ \Delta t^{\frac{1}{2}}\bigl\| e_{3.2}^{m-1} \bigr\| _{0}\bigl\| \nabla e_{3.2}^{m-1}\bigr\| _{0} \bigr]. \end{aligned}$$
(5.5)

According to the definition of \(e_{3.2}^{0}\) we know that \(e_{3.2}^{0}=0\). From (3.3) as \(J=1\) one finds

$$\begin{aligned} \bigl\| d_{t}u_{h}^{1}\bigr\| _{0}^{2}+ \Delta t\bigl\| \nabla d_{t}u_{h}^{1} \bigr\| _{0}^{2}\leq M_{0}. \end{aligned}$$

From (5.5) as \(m=1\)

$$\begin{aligned} &\Delta t\bigl\| d_{t}e_{3.2}^{1}\bigr\| _{0}^{2}+ \bigl\| \nabla e_{3.2}^{1}\bigr\| _{0}^{2} +n\sum _{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.2}^{1} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2} \\ &\quad\leq \frac{\sqrt{3}}{4} \Delta t^{2}\bigl\| d_{t}u_{h}^{1} \bigr\| _{0}\cdot\Delta t^{\frac{1}{2}}\bigl\| \nabla d_{t}u_{h}^{1} \bigr\| _{0}\lesssim\Delta t^{2}. \end{aligned}$$

Then we see that Lemma 5.1 holds for all \(m\leq J-1\), that is to say,

$$\begin{aligned} \bigl\| \nabla e_{3.2}^{J-1}\bigr\| ^{2}_{0}+\Delta t \sum_{m=1}^{J-1}\bigl\| d_{t}e_{3.2}^{m} \bigr\| _{0}^{2}\lesssim\Delta t^{2}. \end{aligned}$$

Now, we take \(m=J\) in (5.5) and use Theorem 3.1 and the trace theorem to obtain

$$\begin{aligned} &\Delta t\bigl\| d_{t}e_{3.2}^{J}\bigr\| _{0}^{2}+ \bigl\| \nabla e_{3.2}^{J}\bigr\| _{0}^{2} +n\sum _{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.2}^{J} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2} \\ &\quad\leq\frac{\sqrt{3}}{4} \bigl[\Delta t^{2}\bigl\| d_{t}u_{h}^{J} \bigr\| _{0}\cdot\Delta t^{\frac{1}{2}}\bigl\| \nabla d_{t}u_{h}^{J} \bigr\| _{0}+ \Delta t^{\frac{1}{2}}\bigl\| e_{3.2}^{J-1} \bigr\| _{0}\bigl\| \nabla e_{3.2}^{J-1}\bigr\| _{0} \bigr] \\ &\qquad{}+\bigl\| \nabla e_{3.2}^{J-1}\bigr\| _{0}^{2} +n\sum_{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.2}^{m-1} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2} \lesssim\Delta t^{2}. \end{aligned}$$

Then we finish the proof of Lemma 5.1. □

For the decoupled Algorithm 3.3, we have the following error estimates.

Lemma 5.2

Let \(({\mathbf{u}}_{fh},\phi_{h})\) and \(({\mathbf {u}}_{3.3}^{m,h},\phi_{3.3}^{m,h})\) be defined by the discrete models (3.2) and (3.9), denote \(e_{3.3}^{m}=({\mathbf{e}}_{3.3}^{m},\xi_{3.3}^{m})=({\mathbf{u}}_{fh}^{m}-{\mathbf {u}}_{3.3}^{m,h},\phi_{h}^{m}-\phi_{3.3}^{m,h})\), under the condition of (3.4) we have

$$\begin{aligned} \bigl\Vert e_{3.3}^{J}\bigr\Vert _{0}^{2}+ \frac{3}{4}\sum_{m=1}^{J}\Delta t\bigl\Vert e_{3.3}^{m}\bigr\Vert _{W}^{2} \lesssim\Delta t^{2}, \qquad \bigl\Vert \nabla e_{3.3}^{J} \bigr\Vert ^{2}_{0}+\Delta t\sum _{m=1}^{J}\bigl\Vert d_{t}e_{3.3}^{m} \bigr\Vert _{0}^{2}\lesssim\Delta t^{2}. \end{aligned}$$

Proof

By comparing the discrete models (3.2) and (3.9), we have the following error equation for all \(v_{h}\in V_{h}\) and \(e_{3.3}^{m}\) with \(m=1,\ldots,J\):

$$\begin{aligned} &\biggl(\frac{e_{3.3}^{m}-e_{3.3}^{m-1}}{\Delta t},v_{h}\biggr)+a_{\Omega} \bigl(e_{3.3}^{m},v_{h}\bigr) +a_{\Gamma} \bigl(u_{h}^{m}-u^{m-1,h}_{3.3},v_{h} \bigr) \\ &\quad =- \int_{\Gamma}n\rho g\psi_{h}\bigl({\mathbf{u}}^{m,h}_{3.3}-{ \mathbf {u}}^{m-1,h}_{3.3}\bigr)\cdot{\mathbf{n}}_{f}. \end{aligned}$$
(5.6)

Taking \(v_{h}=2\Delta te_{3.3}^{m}\) in (5.6), we have

$$\begin{aligned} &\biggl(\frac{e_{3.3}^{m}-e_{3.3}^{m-1}}{\Delta t},2\Delta te_{3.3}^{m} \biggr)+a_{\Omega}\bigl(e_{3.3}^{m},2\Delta te_{3.3}^{m}\bigr) \\ &\quad=-a_{\Gamma}\bigl(u_{h}^{m}-u^{m-1,h}_{3.3},2 \Delta te_{3.3}^{m}\bigr) -2\Delta t \int_{\Gamma}n\rho g\xi_{3.3}^{m}\bigl({ \mathbf{u}}^{m,h}_{3.3}-{\mathbf {u}}^{m-1,h}_{3.3} \bigr)\cdot{\mathbf{n}}_{f}. \end{aligned}$$
(5.7)

Note that the left-hand side of (5.7) can be rewritten as

$$\begin{aligned} &\bigl\| e_{3.3}^{m}\bigr\| _{0}^{2}- \bigl\| e_{3.3}^{m-1}\bigr\| _{0}^{2}+\Delta t^{2}\bigl\| d_{t}e_{3.3}^{m} \bigr\| _{0}^{2} +2\Delta t\bigl\| e_{3.3}^{m} \bigr\| _{W}^{2} \\ &\quad{}+2n\Delta t\sum_{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.3}^{m} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2}. \end{aligned}$$
(5.8)

For the terms of the right-hand side in (5.7), following the trick used in [1], we have

$$\begin{aligned} &a_{\Gamma}\bigl(u_{h}^{m}-u^{m-1,h}_{3.3},2 \Delta te_{3.3}^{m}\bigr) \\ &\quad\leq\frac{1}{2}\Delta t\bigl(\bigl\| e_{3.3}^{m} \bigr\| _{W}^{2}-\bigl\| e_{3.3}^{m-1} \bigr\| _{W}^{2}\bigr) +\Delta t^{2}\bigl\| d_{t}e_{3.3}^{m} \bigr\| _{0}^{2}+\widehat{C}\Delta t^{3} \bigl\| d_{t}u_{h}^{m}\bigr\| _{W}^{2}, \end{aligned}$$
(5.9)
$$\begin{aligned} &2\Delta t \int_{\Gamma}n\rho g\xi_{3.3}^{m}\bigl({ \mathbf{u}}^{m,h}_{3.3}-{\mathbf {u}}^{m-1,h}_{3.3} \bigr)\cdot{\mathbf{n}}_{f} \\ &\quad\leq2\Delta t^{2}n\rho g\bigl\| \xi_{3.3}^{m} \bigr\| _{L^{2}(\Gamma)}\bigl\| d_{t}{\mathbf {u}}^{m,h}_{3.3} \bigr\| _{L^{2}(\Gamma)} \\ &\quad\leq2C_{t}^{2}\Delta t^{\frac{7}{4}}n\rho g\bigl\| \xi_{3.3}^{m}\bigr\| ^{1/2}_{L^{2}(\Omega_{p})} \bigl\| d_{t}{\mathbf{u}}^{m,h}_{3.3}\bigr\| ^{1/2}_{L^{2}(\Omega_{f})} \bigl\| \nabla d_{t}{\mathbf {u}}^{m,h}_{3.3} \bigr\| ^{1/2}_{L^{2}(\Omega_{f})} \cdot\Delta t^{\frac{1}{4}} \bigl\| \nabla \xi_{3.3}^{m}\bigr\| ^{1/2}_{L^{2}(\Omega_{p})} \\ &\quad\leq\frac{1}{4}\Delta t\bigl\| \nabla\xi_{3.3}^{m} \bigr\| ^{2}_{L^{2}(\Omega_{p})}+\frac {3}{4}\bigl(2C_{t}^{2}n \rho g\bigr)^{\frac{4}{3}}\Delta t^{2} \bigl\| d_{t}{ \mathbf{u}}^{m,h}_{3.3}\bigr\| ^{2/3}_{L^{2}(\Omega_{f})}\bigl\| \nabla d_{t}{\mathbf {u}}^{m,h}_{3.3}\bigr\| ^{2/3}_{L^{2}(\Omega_{f})} \\ &\qquad{} \cdot\Delta t^{\frac{1}{3}}\bigl\| \xi_{3.3}^{m} \bigr\| ^{2/3}_{L^{2}(\Omega_{p})} \\ &\quad\leq\frac{1}{4}\Delta t\bigl\| \nabla\xi_{3.3}^{m} \bigr\| ^{2}_{L^{2}(\Omega_{p})}+\frac {1}{3}\Delta t\bigl\| \xi_{3.3}^{m} \bigr\| ^{2}_{L^{2}(\Omega_{p})} \\ &\qquad{}+\frac{\sqrt{3}}{4}\bigl(2C_{t}^{2}n\rho g \bigr)^{2}\Delta t^{\frac{5}{2}} \bigl\| d_{t}{ \mathbf{u}}^{m,h}_{3.3}\bigr\| _{L^{2}(\Omega_{f})}\cdot\Delta t^{\frac{1}{2}}\bigl\| \nabla d_{t}{\mathbf{u}}^{m,h}_{3.3} \bigr\| _{L^{2}(\Omega_{f})}. \end{aligned}$$
(5.10)

Combining (5.8)-(5.10) with (5.7) and summing it over m from 1 to J we obtain

$$\begin{aligned} &\bigl\| e_{3.3}^{J}\bigr\| _{0}^{2} + \frac{3}{4}\sum_{m=1}^{J}\Delta t \bigl\| e_{3.3}^{m}\bigr\| _{W}^{2}+\frac{1}{2} \Delta t\bigl\| e_{3.3}^{m}\bigr\| _{W}^{2} \\ &\quad\leq\Delta t^{2}\Biggl(\frac{\sqrt{3}}{4}(2n\rho g)^{2}\bigl\| d_{t}{\mathbf {u}}^{m,h}_{3.3} \bigr\| _{L^{2}(\Omega_{f})} \cdot\Delta t\sum_{m=1}^{J} \bigl\| \nabla d_{t}{\mathbf{u}}^{m,h}_{3.3} \bigr\| _{L^{2}(\Omega_{f})} \\ &\qquad{}+\widehat{C}\Delta t\sum_{m=1}^{J} \bigl\| d_{t}u_{h}^{m}\bigr\| _{W}^{2} \Biggr) +\frac{1}{3}\Delta t\sum_{m=1}^{J} \bigl\| e_{3.3}^{m}\bigr\| ^{2}_{0}. \end{aligned}$$
(5.11)

Thanks to Theorem 3.1, Lemma 4.1, and the Gronwall lemma, one finds

$$\begin{aligned} \bigl\Vert e_{3.3}^{J}\bigr\Vert _{0}^{2} +\frac{3}{4}\sum _{m=1}^{J}\Delta t\bigl\Vert e_{3.3}^{m} \bigr\Vert _{W}^{2}+\frac{1}{2}\Delta t\bigl\Vert e_{3.3}^{m}\bigr\Vert _{W}^{2} \lesssim \Delta t^{2}. \end{aligned}$$
(5.12)

Taking \(v_{h}=2(e_{3.3}^{m}-e_{3.3}^{m-1})\in V_{h}\) in (5.6), we get

$$\begin{aligned} &\biggl(\frac{e_{3.3}^{m}-e_{3.3}^{m-1}}{\Delta t},2\bigl(e_{3.3}^{m}-e_{3.3}^{m-1} \bigr)\biggr)+2a_{\Omega} \bigl(e_{3.3}^{m},e_{3.3}^{m}-e_{3.3}^{m-1} \bigr) \\ &\quad=-a_{\Gamma} \bigl(u_{h}^{m}-u^{m-1,h}_{3.3},2 \bigl(e_{3.3}^{m}-e_{3.3}^{m-1}\bigr)\bigr) \\ &\qquad{}-2 \int_{\Gamma}n\rho g\bigl(\xi_{3.3}^{m}- \xi_{3.3}^{m-1}\bigr)\cdot\bigl({\mathbf {u}}^{m,h}_{3.3}-{ \mathbf{u}}^{m-1,h}_{3.3}\bigr)\cdot{\mathbf{n}}_{f}. \end{aligned}$$
(5.13)

For the first term of the right-hand side in (5.13), we can treat it as (5.3). For the second term, we have

$$\begin{aligned} &2 \int_{\Gamma}n\rho g\bigl(\xi_{3.3}^{m}- \xi_{3.3}^{m-1}\bigr)\cdot\bigl({\mathbf {u}}^{m,h}_{3.3}-{ \mathbf{u}}^{m-1,h}_{3.3}\bigr)\cdot{\mathbf{n}}_{f} \\ &\quad=2\Delta t^{2} \int_{\Gamma}n\rho gd_{t}\xi_{3.3}^{m} \cdot d_{t}{\mathbf {u}}^{m,h}_{3.3}\cdot{ \mathbf{n}}_{f} \\ &\quad\leq2C_{t}^{2}n\rho g\Delta t^{2} \bigl\| d_{t}\xi_{3.3}^{m}\bigr\| _{L^{2}(\Omega_{p})}^{1/2} \bigl\| d_{t}\xi_{3.3}^{m}\bigr\| ^{1/2}_{H^{1}(\Omega_{p})} \bigl\| d_{t}{\mathbf{u}}^{m,h}_{3.3}\bigr\| _{L^{2}(\Omega_{f})}^{1/2} \bigl\| d_{t}{\mathbf {u}}^{m,h}_{3.3}\bigr\| ^{1/2}_{H^{1}(\Omega_{f})} \\ &\quad=2C_{t}^{2}n\rho g\Delta t^{\frac{3}{2}} \bigl\| d_{t}\xi_{3.3}^{m}\bigr\| _{L^{2}(\Omega_{p})}^{1/2} \bigl\| d_{t}{\mathbf{u}}^{m,h}_{3.3}\bigr\| _{L^{2}(\Omega_{f})}^{1/2} \bigl\| d_{t}{\mathbf {u}}^{m,h}_{3.3}\bigr\| ^{1/2}_{H^{1}(\Omega_{f})} \cdot\Delta t^{\frac{1}{2}}\bigl\| d_{t}\xi_{3.3}^{m} \bigr\| ^{1/2}_{H^{1}(\Omega_{p})} \\ &\quad\leq\frac{1}{4}\Delta t^{2}\bigl\| \nabla d_{t} \xi_{3.3}^{m}\bigr\| _{L^{2}(\Omega _{p})}^{2} \\ &\qquad{}+\frac{3}{4} \bigl(2C_{t}^{2}n\rho g\bigr)^{\frac{4}{3}}\Delta t^{2} \bigl\| d_{t}\xi_{3.3}^{m} \bigr\| _{L^{2}(\Omega_{p})}^{2/3} \bigl\| d_{t}{\mathbf{u}}^{m,h}_{3.3} \bigr\| _{L^{2}(\Omega_{f})}^{2/3}\bigl\| d_{t}{\mathbf {u}}^{m,h}_{3.3} \bigr\| ^{2/3}_{H^{1}(\Omega_{f})} \\ &\quad\leq\frac{1}{4}\Delta t^{2}\bigl\| \nabla d_{t} \xi_{3.3}^{m}\bigr\| _{L^{2}(\Omega _{p})}^{2}+\frac{1}{3} \Delta t \bigl\| d_{t}\xi_{3.3}^{m}\bigr\| _{L^{2}(\Omega_{p})}^{2} \\ &\qquad{}+\frac{\sqrt{3}}{4}\bigl(2C_{t}^{2}n\rho g \bigr)^{2}\Delta t^{2} \bigl\| d_{t}{ \mathbf{u}}^{m,h}_{3.3}\bigr\| _{L^{2}(\Omega_{f})}\cdot\Delta t^{\frac{1}{2}}\bigl\| d_{t}{\mathbf{u}}^{m,h}_{3.3} \bigr\| _{H^{1}(\Omega_{f})}. \end{aligned}$$
(5.14)

Combining (5.3), (5.14) with (5.13) we obtain

$$\begin{aligned} &\frac{2}{3}\Delta t\bigl\| d_{t}e_{3.3}^{m} \bigr\| _{0}^{2}+\bigl\| \nabla e_{3.3}^{m} \bigr\| _{0}^{2}-\bigl\| \nabla e_{3.3}^{m-1} \bigr\| _{0}^{2} \\ &\qquad{}+n\sum_{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.3}^{m} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2} -n\sum _{i=1}^{d-1} \int_{\Gamma}\frac{\alpha}{\sqrt{{\boldsymbol {\tau}}_{i}\cdot{\mathbf{K}}{\boldsymbol {\tau }}_{i}}}\bigl({\mathbf{e}}_{3.3}^{m-1} \cdot{\boldsymbol {\tau}}_{i}\bigr)^{2} \\ &\quad\leq\frac{\sqrt{3}}{4} \bigl[\bigl(1+\bigl(2C_{t}^{2}n \rho g\bigr)^{2}\bigr)\Delta t^{2}\bigl\| d_{t}u_{3.3}^{m,h} \bigr\| _{0}\cdot\Delta t^{\frac{1}{2}}\bigl\| \nabla d_{t}u_{3.3}^{m,h} \bigr\| _{0}+ \Delta t^{\frac{1}{2}}\bigl\| e_{3.3}^{m-1} \bigr\| _{0}\bigl\| \nabla e_{3.3}^{m-1}\bigr\| _{0} \bigr]. \end{aligned}$$

By the induction method and Lemma 4.1, similar to Lemma 5.1, we complete the proof. □

For the decoupled Algorithm 3.4, we have the following error estimates.

Lemma 5.3

Let \(({\mathbf{u}}_{fh},\phi_{h})\) and \(({\mathbf {u}}_{3.4}^{m,h},\phi_{3.4}^{m,h})\) be defined by the discrete models (3.2) and (3.10). Denote \(e_{3.4}^{m}=({\mathbf{e}}_{3.4}^{m},\xi_{3.4}^{m})=({\mathbf{u}}_{fh}^{m}-{\mathbf {u}}_{3.4}^{m,h},\phi_{h}^{m}-\phi_{3.4}^{m,h})\), under the condition of (3.4) we have

$$\begin{aligned} \bigl\| e_{3.4}^{J}\bigr\| _{0}^{2}+ \frac{3}{4}\sum_{m=1}^{J}\Delta t \bigl\| e_{3.4}^{m}\bigr\| _{W}^{2}\lesssim\Delta t^{2}, \qquad \bigl\| \nabla e_{3.4}^{J} \bigr\| ^{2}_{0}+\Delta t\sum_{m=1}^{J} \bigl\| d_{t}e_{3.4}^{m}\bigr\| _{0}^{2} \lesssim\Delta t^{2}. \end{aligned}$$

Proof

By comparing the discrete models (3.2) and (3.10), we have the following error equation for all \(v_{h}\in V_{h}\) and \(e_{3.4}^{m}\) with \(m=1,\ldots,J\):

$$\begin{aligned} &\biggl(\frac{e_{3.4}^{m}-e_{3.4}^{m-1}}{\Delta t},v_{h}\biggr)+a_{\Omega } \bigl(e_{3.4}^{m},v_{h}\bigr)+a_{\Gamma} \bigl(u_{h}^{m}-u^{m-1,h}_{3.4},v_{h} \bigr) \\ &\quad=- \int_{\Gamma}n\rho g\bigl(\phi_{3.4}^{m,h}- \phi_{3.4}^{m-1,h}\bigr)\cdot{\mathbf {v}}_{h}\cdot{ \mathbf{n}}_{f}. \end{aligned}$$
(5.15)

Choosing \(v_{h}=2\Delta te_{3.4}^{m}\) in (5.15), we have

$$\begin{aligned} &\biggl(\frac{e_{3.4}^{m}-e_{3.4}^{m-1}}{\Delta t},2\Delta te_{3.4}^{m} \biggr)+a_{\Omega}\bigl(e_{3.4}^{m},2\Delta te_{3.4}^{m}\bigr) \\ &\quad=-a_{\Gamma}\bigl(u_{h}^{m}-u^{m-1,h}_{3.4},2 \Delta te_{3.4}^{m}\bigr) -2\Delta t \int_{\Gamma}n\rho g\bigl(\phi_{3.4}^{m,h}- \phi_{3.4}^{m-1,h}\bigr)\cdot {\mathbf{e}}_{3.4}^{m} \cdot{\mathbf{n}}_{f}. \end{aligned}$$
(5.16)

For the last term of the right-hand side in (5.16) we have

$$\begin{aligned} &2\Delta t \int_{\Gamma}n\rho g\bigl(\phi_{3.4}^{m,h}-\phi _{3.4}^{m-1,h}\bigr)\cdot{\mathbf{e}}_{3.4}^{m} \cdot{\mathbf{n}}_{f} \\ &\quad\leq2\Delta t^{2}n\rho g\bigl\| {\mathbf{e}}_{3.4}^{m} \bigr\| _{L^{2}(\Gamma)}\bigl\| d_{t}\phi ^{m,h}_{3.4} \bigr\| _{L^{2}(\Gamma)} \\ &\quad\leq2C_{t}^{2}\Delta t^{\frac{7}{4}}n\rho g \bigl\| {\mathbf{e}}_{3.4}^{m}\bigr\| ^{1/2}_{L^{2}(\Omega_{p})} \bigl\| d_{t}\phi^{m,h}_{3.4}\bigr\| ^{1/2}_{L^{2}(\Omega_{f})} \bigl\| \nabla d_{t}\phi ^{m,h}_{3.4}\bigr\| ^{1/2}_{L^{2}(\Omega_{f})} \cdot\Delta t^{\frac{1}{4}} \bigl\| \nabla{\mathbf{e}}_{3.4}^{m} \bigr\| ^{1/2}_{L^{2}(\Omega_{p})} \\ &\quad\leq\frac{1}{4}\Delta t\bigl\| \nabla{\mathbf{e}}_{3.4}^{m} \bigr\| ^{2}_{L^{2}(\Omega _{p})} \\ &\qquad{}+\frac{3}{4}\bigl(2C_{t}^{2}n \rho g\bigr)^{\frac{4}{3}}\Delta t^{2} \bigl\| d_{t} \phi^{m,h}_{3.4}\bigr\| ^{2/3}_{L^{2}(\Omega_{f})}\bigl\| \nabla d_{t}\phi ^{m,h}_{3.4}\bigr\| ^{2/3}_{L^{2}(\Omega_{f})} \cdot\Delta t^{\frac{1}{3}}\bigl\| {\mathbf{e}}_{3.4}^{m} \bigr\| ^{2/3}_{L^{2}(\Omega_{p})} \\ &\quad\leq\frac{1}{4}\Delta t\bigl\| \nabla{\mathbf{e}}_{3.4}^{m} \bigr\| ^{2}_{L^{2}(\Omega _{p})}+\frac{1}{3}\Delta t\bigl\| {\mathbf{e}}_{3.4}^{m} \bigr\| ^{2}_{L^{2}(\Omega_{p})} \\ &\qquad{}+\frac{\sqrt{3}}{4}\bigl(2C_{t}^{2}n\rho g \bigr)^{2}\Delta t^{2} \bigl\| d_{t}\phi^{m,h}_{3.4} \bigr\| _{L^{2}(\Omega_{f})}\cdot\Delta t\bigl\| \nabla d_{t}\phi ^{m,h}_{3.4} \bigr\| _{L^{2}(\Omega_{f})}. \end{aligned}$$
(5.17)

Combining (5.8)-(5.9), (5.17) with (5.16), one finds

$$\begin{aligned} &\bigl\| e_{3.4}^{m}\bigr\| _{0}^{2}- \bigl\| e_{3.4}^{m-1}\bigr\| _{0}^{2} + \frac{3}{4}\Delta t\bigl\| e_{3.4}^{m}\bigr\| _{W}^{2}+ \frac{1}{2}\Delta t\bigl(\bigl\| e_{3.4}^{m}\bigr\| _{W}^{2}-\bigl\| e_{3.4}^{m} \bigr\| _{W}^{2}\bigr) \\ &\quad\leq\Delta t^{2}\biggl(\frac{\sqrt{3}}{4}\bigl(2C_{t}^{2}n \rho g\bigr)^{2} \bigl\| d_{t}\phi^{m,h}_{3.4} \bigr\| _{L^{2}(\Omega_{f})}\cdot\Delta t\bigl\| \nabla d_{t}\phi ^{m,h}_{3.4} \bigr\| _{L^{2}(\Omega_{f})} +\widehat{C}\Delta t\bigl\| d_{t}u_{h}^{m} \bigr\| _{W}^{2}\biggr) \\ &\qquad{}+\frac{1}{3}\Delta t \bigl\| e_{3.4}^{m}\bigr\| ^{2}_{0}. \end{aligned}$$
(5.18)

Summing (5.18) over m from 1 to J and using Theorem 3.1, Lemma 4.2, and the Gronwall lemma, we have

$$\begin{aligned} \bigl\| e_{3.4}^{J}\bigr\| _{0}^{2} + \frac{3}{4}\sum_{m=1}^{J}\Delta t \bigl\| e_{3.4}^{m}\bigr\| _{W}^{2}+\frac{1}{2} \Delta t\bigl\| e_{3.4}^{m}\bigr\| _{W}^{2} \lesssim \Delta t^{2}. \end{aligned}$$
(5.19)

Next, we take \(v_{h}=2(e_{3.4}^{m}-e_{3.4}^{m-1})\) in (5.15) and obtain

$$\begin{aligned} &\biggl(\frac{e_{3.4}^{m}-e_{3.4}^{m-1}}{\Delta t},2\bigl(e_{3.4}^{m}-e_{3.4}^{m-1} \bigr)\biggr)+a_{\Omega} \bigl(e_{3.4}^{m},2 \bigl(e_{3.4}^{m}-e_{3.4}^{m-1}\bigr)\bigr) \\ &\quad=-a_{\Gamma} \bigl(u_{h}^{m}-u^{m-1,h}_{3.4},2 \bigl(e_{3.4}^{m}-e_{3.4}^{m-1}\bigr)\bigr) \\ &\qquad{}-2 \int_{\Gamma}n\rho g\bigl(\phi_{3.4}^{m,h}- \phi_{3.4}^{m-1,h}\bigr)\cdot\bigl({\mathbf {e}}_{3.4}^{m}-{ \mathbf{e}}_{3.4}^{m-1}\bigr)\cdot{\mathbf{n}}_{f}. \end{aligned}$$
(5.20)

For the first term in the right-hand side of (5.20), we can treat it as in (5.3). For the second term, we have

$$\begin{aligned} &2 \int_{\Gamma}n\rho g\bigl(\phi_{3.4}^{m,h}- \phi_{3.4}^{m-1,h}\bigr)\cdot\bigl({\mathbf {e}}_{3.4}^{m}-{ \mathbf{e}}_{3.4}^{m-1}\bigr)\cdot{\mathbf{n}}_{f} \\ &\quad=2 \Delta t^{2} \int_{\Gamma}n\rho gd_{t}\phi_{3.4}^{m,h} \cdot d_{t}{\mathbf {e}}_{3.4}^{m}\cdot{ \mathbf{n}}_{f} \\ &\quad\leq2C_{t}^{2}\Delta t^{2}n\rho g \bigl\| d_{t}\phi_{3.4}^{m,h}\bigr\| ^{1/2}_{L^{2}(\Omega _{p})} \bigl\| d_{t}\phi_{3.4}^{m,h}\bigr\| ^{1/2}_{H^{1}(\Omega_{p})} \bigl\| d_{t}{\mathbf{e}}_{3.4}^{m}\bigr\| ^{1/2}_{L^{2}(\Omega_{f})} \bigl\| d_{t}{\mathbf{e}}_{3.4}^{m}\bigr\| ^{1/2}_{H^{1}(\Omega_{f})} \\ &\quad\leq\frac{1}{4}\Delta t^{2}\bigl\| \nabla d_{t}{ \mathbf{e}}_{3.4}^{m}\bigr\| _{L^{2}(\Omega _{p})}^{2}+ \frac{1}{3}\Delta t \bigl\| d_{t}{\mathbf{e}}_{3.4}^{m} \bigr\| _{L^{2}(\Omega_{p})}^{2} \\ &\qquad{}+\frac{\sqrt{3}}{4}\bigl(2C_{t}^{2}n\rho g \bigr)^{2}\Delta t^{2} \bigl\| d_{t}\phi_{3.4}^{m,h} \bigr\| _{L^{2}(\Omega_{f})}\cdot\Delta t^{\frac{1}{2}}\bigl\| d_{t} \phi_{3.4}^{m,h}\bigr\| _{H^{1}(\Omega_{f})}. \end{aligned}$$
(5.21)

Combining (5.3), (5.19), (5.20) with (5.21), using the induction method we complete the proof. □

Finally, combining Theorem 3.2 with Lemmas 5.1-5.3, we have the following conclusion for the decoupled Algorithms 3.2-3.4.

Corollary

Under the condition (3.4), for the decoupled algorithms DEBS1, DEBS2, and DEBS3 with \(m=1,\ldots,J\), we have

$$\begin{aligned}& \bigl\| u(t_{m})-u_{3.i}^{m,h}\bigr\| _{0}\lesssim \Delta t+h^{2},\quad i=2,3,4, \\& \bigl\| \nabla\bigl(u(t_{m})-u_{3.i}^{m,h}\bigr) \bigr\| _{0}\lesssim\Delta t+h,\quad i=2,3,4. \end{aligned}$$

6 Numerical experiments

In order to gain insights on the established theoretical results in the previous section, we present some numerical tests in this section. Our main interest is to verify the performances of the decoupled Algorithms 3.3 and 3.4. In our experiments, let the domain Ω be composed of \(\Omega _{f}=[0,1]\times(1,2]\) and \(\Omega_{p}=[0,1]\times(0,1]\) with the interface \(\Gamma=[0,1]\times\{1\}\). The model parameters \(\rho, g, n\), and α are simply set to 1. The boundary conditions and right-hand side functions in the model are selected such that the exact solution is given by

$$\begin{aligned} \left \{ \textstyle\begin{array}{@{}l} u=(x^{2}(y-1)^{2}+y)\cos t,\\ v=(- \frac{2}{3}x(y-1)^{3}+2-\pi\sin(\pi x))\cos t,\\ p_{f}=(2-\pi\sin(\pi x))\cdot\sin(\frac{1}{2}\pi y)\cos t,\\ \phi=(2-\pi\sin(\pi x))\cdot(1-y-\cos(\pi y))\cos t, \end{array}\displaystyle \right . \end{aligned}$$

where the components of \({\mathbf{u}}_{f}\) are denoted by \((u,v)\) for convenience.

In order to show the prominent features of the decoupled Algorithms 3.3 and 3.4, we compare the numerical results of algorithms (3.9) and (3.10) with the coupled method (3.2) and DBES1 (3.6). The finite element spaces are constructed using the well-known MINI elements for the Stokes problem and linear Lagrangian elements for the Darcy flow. For the coupled scheme, the GMRES routine is used to solve the coupled system. For the decoupled schemes, a Gauss Lower and Upper triangular matrix factorization is implemented to solve the positive definite matrix subsystems. In the following tables, we will use \(\|\cdot\|_{0}\) to denote the \(L^{2}\)-norm and \(\|\cdot\|_{1}\) to denote the \(H^{1}\)-norm.

First, we compare the errors, convergence rates and CPU times for both the coupled scheme and the decoupled algorithms. In Tables 1-4, we consider these schemes at time \(T=1.0\), with varying mesh h but fixed time step Δt. These numerical algorithms achieve similar precision, although the coupled scheme is slightly more accurate than the decoupled schemes. However, the coupled scheme takes much more CPU time than the decoupled algorithms. On the other hand, we consider both the coupled and the decoupled algorithms with varying time step Δt but fixed mesh \(h=\frac{1}{32}\). In Tables 5-8, four schemes almost get the same accuracy, but the decoupled schemes need much less CPU time than the coupled scheme. Stated succinctly, the decoupled schemes are comparable with the coupled scheme but cheaper and more efficient.

Table 1 The convergence performance and CPU time of the coupled Algorithm 3.1 at time \(\pmb{T=1.0}\) , with varying mesh h but fixed time step \(\pmb{\Delta t=0.01}\)
Full size table
Table 2 The convergence performance and CPU time of the decoupled Algorithm  3.2 at time \(\pmb{T=1.0}\) , with varying mesh h but fixed time step \(\pmb{\Delta t=0.01}\)
Full size table
Table 3 The convergence performance and CPU time of the decoupled Algorithm  3.3 at time \(\pmb{T=1.0}\) , with varying mesh h but fixed time step \(\pmb{\Delta t=0.01}\)
Full size table
Table 4 The convergence performance and CPU time of the decoupled Algorithm  3.4 at time \(\pmb{T=1.0}\) , with varying mesh h but fixed time step \(\pmb{\Delta t=0.01}\)
Full size table
Table 5 The convergence performance and CPU time of the coupled Algorithm  3.1 at time \(\pmb{T=1.0}\) , with varying time step Δ t but fixed mesh \(\pmb{h=\frac{1}{32}}\)
Full size table
Table 6 The convergence performance and CPU time of the decoupled Algorithm  3.2 at time \(\pmb{T=1.0}\) , with varying time step Δ t but fixed mesh \(\pmb{h=\frac{1}{32}}\)
Full size table
Table 7 The convergence performance and CPU time of the decoupled Algorithm  3.3 at time \(\pmb{T=1.0}\) , with varying time step Δ t but fixed mesh \(\pmb{h=\frac{1}{32}}\)
Full size table
Table 8 The convergence performance and CPU time of the decoupled Algorithm  3.4 at time \(\pmb{T=1.0}\) , with varying time step Δ t but fixed mesh \(\pmb{h=\frac{1}{32}}\)
Full size table

Next, we focus on the decoupled schemes and examine the orders of convergence with respect to the mesh size h or the time step Δt. Following [1], we introduce the following approach to examine the orders of convergence with respect to the time step Δt or the mesh size h due to the approximation errors \(\mathcal{O}(\Delta t^{\gamma})+\mathcal{O}(h^{\mu})\). For example, assuming

$$\begin{aligned} v_{h}^{\Delta t}\thickapprox v(x,t_{m})+C_{1}(x,t_{m}) \Delta t^{\gamma}+C_{2}(x,t_{m})h^{\mu}, \end{aligned}$$

thus we have

$$\begin{aligned}& \rho_{v,h,j}=\frac{\|v_{h}^{\Delta t}(x,t_{m})-v_{\frac{h}{2}}^{\Delta t}(x,t_{m})\|_{j}}{\|v_{\frac{h}{2}}^{\Delta t}(x,t_{m})-v_{\frac {h}{4}}^{\Delta t}(x,t_{m})\|_{j}}\thickapprox\frac{4^{\mu}-2^{\mu}}{2^{\mu}-1} , \\& \rho_{v,\Delta t,j}=\frac{\|v_{h}^{\Delta t}(x,t_{m})-v_{h}^{\frac{\Delta t}{2}}(x,t_{m})\|_{j}}{\|v_{h}^{\frac{\Delta t}{2}}(x,t_{m})-v_{h}^{\frac{\Delta t}{4}}(x,t_{m})\|_{j}}\thickapprox\frac{4^{\gamma}-2^{\gamma}}{2^{\gamma}-1}. \end{aligned}$$

Here, v can take \({\mathbf{u}}_{f},p_{f},\phi\) and j can be 0 or 1. While \(\rho_{v,h,j},\rho_{v,\Delta t,j}\) approach 4.0 or 2.0, the convergence order will be 2.0 or 1.0, respectively.

In Tables 9-11, we study the convergence order with a fixed time step \(\Delta t=0.01\) and varying spacing \(h=\frac{1}{2},\frac{1}{4},\frac {1}{8},\frac{1}{16}\). Observe that \(\rho_{{\mathbf{u}}_{f},h,0},\rho_{\phi ,h,0}\) close to 4.0 and \(\rho_{{\mathbf{u}}_{f},h,1},\rho_{p_{f},h,0},\rho_{\phi,h,1}\) approach to 2.0, which suggest that the error estimates \(\mathcal{O}(h^{2})\) for the \(L^{2}\)-norm of \({\mathbf{u}}^{m,h}_{3.i}\) and \(\phi^{m,h}_{3.i}\) (\(i=2,3,4\)), and \(\mathcal{O}(h)\) for the \(H^{1}\)-norm of \({\mathbf{u}}^{m,h}_{3.i}\) and \(\phi ^{m,h}_{3.i}\) (\(i=2,3,4\)) in space for three decoupled algorithms. However, in Tables 12-14, we study the convergence order with a fixed spacing \(h=\frac{1}{32}\) and varying time step \(\Delta t=0.1,0.05,0.025,0.0125\). The numerical experiments strongly suggest that the orders of convergence in time are \(\mathcal{O}(\Delta t)\), which implies that the error estimates for both the \(L^{2}\)-norm and \(H^{1}\)-norm of \({\mathbf{u}}^{m,h}_{3.i}\) and \(\phi^{m,h}_{3.i} \) (\(i=2,3,4\)) are optimal. Our numerical results confirm the established theoretical analysis very well. Furthermore, we observe that Algorithm 3.3 may be the best one among four algorithms to treat the non-stationary Stokes-Darcy model due to the fact that this algorithm not only keeps good accuracy but also takes less computational time.

Table 9 Convergence orders of \(\pmb{\mathcal {O}(h^{\mu})}\) of the decoupled Algorithm  3.2 at time \(\pmb{T=1.0}\) , with varying mesh h but fixed time step \(\pmb{\Delta t=0.01}\)
Full size table
Table 10 Convergence orders of \(\pmb{\mathcal {O}(h^{\mu})}\) of the decoupled Algorithm  3.3 at time \(\pmb{T=1.0}\) , with varying mesh h but fixed time step \(\pmb{\Delta t=0.01}\)
Full size table
Table 11 Convergence orders of \(\pmb{\mathcal {O}(h^{\mu})}\) of the decoupled Algorithm  3.4 at time \(\pmb{T=1.0}\) , with varying mesh h but fixed time step \(\pmb{\Delta t=0.01}\)
Full size table
Table 12 Convergence orders of \(\pmb{\mathcal {O}(\Delta t^{\gamma})}\) of the decoupled Algorithm  3.2 at time \(\pmb{T=1.0}\) , with varying time step Δ t but fixed mesh \(\pmb{h=\frac{1}{32}}\)
Full size table
Table 13 Convergence orders of \(\pmb{\mathcal {O}(\Delta t^{\gamma})}\) of the decoupled Algorithm  3.3 at time \(\pmb{T=1.0}\) , with varying time step Δ t but fixed mesh \(\pmb{h=\frac{1}{32}}\)
Full size table
Table 14 Convergence orders of \(\pmb{\mathcal {O}(\Delta t^{\gamma})}\) of the decoupled Algorithm  3.4 at time \(\pmb{T=1.0}\) , with varying time step Δ t but fixed mesh \(\pmb{h=\frac{1}{32}}\)
Full size table