[2]Mingchao Cai

1]Key Laboratory of NSLSCS, Ministry of Education; Jiangsu International Joint
Laboratory of BDMCA; School of Mathematical Sciences, Nanjing Normal
University, Nanjing, 210023, Jiangsu, P.R. China

[2]Department of Mathematics, Morgan State University, Baltimore, 21251,
Maryland, USA

# An Unconditionally Stable Iterative Decoupled Algorithm for Multiple-Network
Poroelasticity

Meng Lei<EMAIL_ADDRESS><EMAIL_ADDRESS>Feng Wang
<EMAIL_ADDRESS>[  *

###### Abstract

In this work, we introduce an iterative decoupled algorithm designed for
addressing the quasi-static multiple-network poroelasticity problem. This
problem pertains to the simultaneous modeling of fluid flow and deformations
within an elastic porous medium permeated by multiple fluid networks, each
with distinct characteristics. Our approach focuses on the total-pressure-
based formulation, which treats the solid displacement, total pressure, and
network pressures as primary unknowns. This formulation transforms the
original problem into a combination of the generalized Stokes problem and the
parabolic problem, offering certain advantages such as mitigating elastic
locking effects and streamlining the discretization process. Notably, the
algorithm ensures unconditional convergence to the solution of the total-
pressure-based coupled algorithm. To validate the accuracy and efficiency of
our method, we present numerical experiments. The robustness of the algorithm
with respect to the physical parameters and the discretization parameters is
carefully investigated.

###### keywords:

multiple-network poroelasticity, iterative decoupled algorithm, total pressure

##  1 Introduction

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded polygonal domain. The quasi-
static multiple-network poroelasticity problem [1] is to find the displacement
$\boldsymbol{u}$ and the network pressures $\vec{p}$ such that

$\displaystyle-\operatorname{div}\sigma(\boldsymbol{u})+\nabla(\vec{\alpha}^{\mathsf{T}}\vec{p})=\boldsymbol{f}\quad\text{
in }\Omega\times[0,T],$ (1a)
$\displaystyle\vec{\alpha}\operatorname{div}\dot{\boldsymbol{u}}+S\dot{\vec{p}}+B\vec{p}-\operatorname{div}(K\nabla\vec{p})=\vec{g}\quad\text{
in }\Omega\times[0,T].$ (1b)

Here,
$\boldsymbol{u}=(u_{1}(\boldsymbol{x},t),u_{2}(\boldsymbol{x},t))^{\mathsf{T}}$
and
$\vec{p}=(p_{1}(\boldsymbol{x},t),p_{2}(\boldsymbol{x},t),\cdots,p_{N}(\boldsymbol{x},t))^{\mathsf{T}}$,
with a given number of networks $N$. The operators and parameters are defined
as follows. The effective stress and the strain tensor are denoted by
$\sigma(\boldsymbol{u})=2\mu\epsilon(\boldsymbol{u})+\lambda\operatorname{div}(\boldsymbol{u})\boldsymbol{I}$
and
$\epsilon(\boldsymbol{u})=\frac{1}{2}\left(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{\mathsf{T}}\right)$
, respectively. The Lamé parameters $\lambda$ and $\mu$ are expressed in terms
of the Young modulus $E$ and Poisson ratio $\nu\in[0,\frac{1}{2})$ by
$\lambda=\frac{\nu E}{(1+\nu)(1-2\nu)}$ and $\mu=\frac{E}{2(1+\nu)}$,
respectively. Column vector
$\vec{\alpha}=(\alpha_{1},\alpha_{2},\cdots,\alpha_{N})^{\mathsf{T}}$, where
$\alpha_{i}\in(0,1]$ is the Biot-Willis coefficient. Matrix
$S=\text{diag}(c_{1},c_{2},\cdots,c_{N})$, and $c_{i}\geq 0$ is the storage
coefficient. Matrix $K=\text{diag}(K_{1},K_{2},\cdots,K_{N})$, and $K_{i}$ is
the hydraulic conductivity coefficient. The coefficient matrix $B$ satisfies
$(B\vec{p})_{i}=\sum_{j=1,j\neq i}^{N}\beta_{ij}\left(p_{i}-p_{j}\right)$,
where the non-negative network transfer coefficient $\beta_{ij}$ couple the
network pressures and $\beta_{ij}=\beta_{ji},1\leq i,j\leq N,j\neq i$.
Further, $\boldsymbol{f}=(f_{1},f_{2})^{\mathsf{T}}$, where $f_{i}$ represents
the body force and $\vec{g}=(g_{1},g_{2},\cdots,g_{N})^{\mathsf{T}}$, where
$g_{i}$ represents the source in $i$th network.

The system (1) is well-posed with proper boundary and initial conditions. In
this paper, we consider the mixed boundary conditions, assuming
$\partial\Omega=\Gamma_{\boldsymbol{u},D}\cup\Gamma_{\boldsymbol{u},N}=\Gamma_{\vec{p},D}\cup\Gamma_{\vec{p},N},|\Gamma_{\boldsymbol{u},D}\cap\Gamma_{\boldsymbol{u},N}|=0,|\Gamma_{\vec{p},D}\cap\Gamma_{\vec{p},N}|=0,|\Gamma_{\boldsymbol{u},D}|>0,|\Gamma_{\vec{p},D}|>0$.
More clearly, the following boundary and initial conditions [2] are imposed.

$\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{0}\quad\text{on}\,\Gamma_{\boldsymbol{u},D},$
$\displaystyle\left(\sigma(\boldsymbol{u})-(\vec{\alpha}^{\mathsf{T}}\vec{p})\boldsymbol{I}\right)\boldsymbol{n}$
$\displaystyle=\boldsymbol{h}\quad\text{on}\,\Gamma_{\boldsymbol{u},N},$
$\displaystyle\vec{p}$
$\displaystyle=\vec{0}\quad\text{on}\,\Gamma_{\vec{p},D},$
$\displaystyle\left(K\nabla\vec{p}\right)\boldsymbol{n}$
$\displaystyle=\vec{l}\quad\text{on}\,\Gamma_{\vec{p},N},$
$\displaystyle\boldsymbol{u}(0)$
$\displaystyle=\boldsymbol{u}_{0}\quad\text{in}\,\Omega,$
$\displaystyle\vec{p}(0)$ $\displaystyle=\vec{p}_{0}\quad\text{in}\,\Omega.$

The concept of multiple-network poroelasticity has been incorporated into the
field of geomechanics [3] as a means to reflect mechanical deformation and
fluid flow within porous materials. The multiple-network poroelasticity
problem serves as an extension of Biot’s model, which has been applied to
geomechanics, geophysics, and biology. For instance, when considering the case
with $N=2$, one can obtain the Biot-Barenblatt model [4, 5], which states
consolidation processes within a fluid-saturated double-diffusion model of
fractured rock. In another context, Tully and Ventikos investigated a scenario
involving four distinct networks ($N=4$) as a macroscopic model to describe
the dynamics of fluid flows within brain tissue [6]. More recently, the
multiple-network poroelasticity model has also been employed to gain a deeper
understanding of the impact of biomechanical risk factors associated with the
early stages of Alzheimer’s disease [7].

The discretization of Biot’s equations is widely recognized as a challenging
task, due to the poroelastic locking phenomenon. Poroelastic locking is
characterized by two main features [2], (1) underestimation of the solid
deformation if the material is close to being incompressible, (2) nonphysical
pressure oscillations. Many existing studies, such as [8, 9], adopt solid
displacement and fluid pressure as the primary variables within Biot’s
consolidation model. However, it has been noted that the elastic locking
phenomenon occurs in a model based on the two-field formulation [10, 11]. To
address these challenges, researchers have employed various approaches,
including the DG method [12], stabilization techniques [13, 14], and various
three-field or four-field reformulations [15, 16, 17, 18, 19, 20, 21, 22]. In
spite of the extensive research dedicated to exploring Biot’s consolidation
model, the multiple-network poroelasticity problem has received much less
attention. Hong et al. introduced parameter-robust preconditioners [23],
proposed a fixed-stress split method [1] and presented an augmented Lagrangian
Uzawa algorithm [24] for multiple-network flux-based poroelasticity equations.
Lee et al. analyzed a mixed finite element method [2] and proposed two
partitioned numerical algorithms [25] for the multiple-network poroelasticity
problem.

By introducing an intermediate variable
$\xi=\vec{\alpha}^{\mathsf{T}}\vec{p}-\lambda\operatorname{div}\boldsymbol{u}$,
we studied a formulation which has the solid displacement, the total pressure,
and the network pressures as the primary unknowns. This formulation interprets
the original problem as a combination of the generalized Stokes problem and
the parabolic problem. The solid displacement and the total pressure are
approximated by the classical stable Stokes finite elements, while the network
pressures are discretized by the Lagrange finite elements. Importantly, this
approach eliminates the need for any assumptions regarding additional physical
parameters, such as storage coefficients $c_{i}>0,i=1,2,\cdots,N$. The article
presents an iterative decoupled algorithm for the multiple-network
poroelasticity problem (1). Unlike the fixed-stress split method for multiple-
network flux-based poroelasticity problem proposed in [1], the iterative
decoupled algorithm does not require any stabilization parameters and
unconditionally converges to the solution of the total-pressure-based coupled
algorithm. In order to assess both the precision and effectiveness of our
approach, we provide numerical experiments. Furthermore, we meticulously
examine the algorithm’s resilience concerning variations in the physical
parameters and discretization parameters.

The structure of this article is as follows. In Section 2, we present the
total-pressure-based formulation which has the solid displacement, the total
pressure, and the network pressures as the primary unknowns and the
corresponding variational problems for the multiple-network poroelasticity
equations. The iterative decoupled algorithm is given in Section 3. In Section
4, we analyze the convergence of the iterative decoupled algorithm. The error
analysis of the coupled algorithm is studied in Appendix A. Numerical
experiments are presented in Section 5, and conclusions are drawn in Section
6.

##  2 The total-pressure-based formulation and the weak form

###  2.1 Notation and preliminaries

We use $L^{2}(\Omega)$ and $L^{2}(\partial\Omega)$ to denote square-integrable
real-valued functions defined on $\Omega$ and $\partial\Omega$ respectively.
The inner product of $L^{2}(\Omega)$ and the induced norm are denoted by
$(\cdot,\cdot)$ and $|\\!|\cdot|\\!|_{L^{2}(\Omega)}$ respectively. We use
$\langle\cdot,\cdot\rangle$ to denote the $L^{2}(\partial\Omega)$ inner
product. For vector-valued or matrix-valued functions where each component
belongs to $L^{2}(\Omega)$, we use the same notation $(\cdot,\cdot)$ and
$|\\!|\cdot|\\!|_{L^{2}(\Omega)}$ to denote the inner product and norm. For
vector-valued functions where each component belongs to
$L^{2}(\partial\Omega)$, we also use the notation $\langle\cdot,\cdot\rangle$
to denote the inner product.

For a non-negative integer $m$ and $1<p<+\infty$, we introduce the following
Sobolev spaces, $W^{m,p}(\Omega)=\\{u|D^{\alpha}u\in
L^{p}(\Omega),0\leq|\alpha|\leq m,|\\!|u|\\!|_{W^{m,p}(\Omega)}<+\infty\\}$;
When $p=2$, we use $H^{m}(\Omega)$ to denote $W^{m,2}(\Omega)$. In cases
involving vector-valued functions where each component belongs to
$H^{m}(\Omega)$, we consistently employ the notation
$|\\!|\cdot|\\!|_{H^{m}(\Omega)}$ to represent the norm. Furthermore, when
$m\geq 1$, we use $H^{m}_{0,\Gamma}(\Omega)$ to denote the subspace of
$H^{m}(\Omega)$ comprising functions with a vanishing trace on
$\Gamma\subset\partial\Omega$. For ease of presentation, we assume that
$\Omega$ is a two-dimensional domain in this work.

###  2.2 The total-pressure-based formulation and the variational problem

In the following, we introduce an intermediate variable to rewrite the problem
(1) into the total-pressure-based formulation and present the corresponding
variational problem [2]. More clearly, we introduce the so-called ”total
pressure”,
$\xi=\vec{\alpha}^{\mathsf{T}}\vec{p}-\lambda\operatorname{div}\boldsymbol{u}$.
Then (1) can be rewritten as

$\displaystyle-2\mu\operatorname{div}\left(\varepsilon(\boldsymbol{u})\right)+\nabla\xi=\boldsymbol{f},$
(2a)
$\displaystyle-\operatorname{div}\boldsymbol{u}-\frac{1}{\lambda}\xi+\frac{1}{\lambda}\vec{\alpha}^{\mathsf{T}}\vec{p}=0,$
(2b)
$\displaystyle\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)\dot{\vec{p}}-\frac{1}{\lambda}\vec{\alpha}\dot{\xi}-\operatorname{div}(K\nabla\vec{p})+B\vec{p}=\vec{g}.$
(2c)

The initial and boundary conditions we provided earlier can still be applied
to the problem (2).

In order to describe the variational problem for the total-pressure-based
formulation (2), we will utilize the following functional spaces.

$\displaystyle\boldsymbol{V}$
$\displaystyle:=\\{\boldsymbol{v}\in[H^{1}(\Omega)]^{2};\boldsymbol{v}|_{\Gamma_{\boldsymbol{u},D}}=\boldsymbol{0}\\},$
$\displaystyle W$ $\displaystyle:=L^{2}(\Omega),$ $\displaystyle M$
$\displaystyle:=\\{\vec{q}\in[H^{1}(\Omega)]^{N};\vec{q}|_{\Gamma_{\vec{p},D}}=\vec{0}\\}.$

The properties of these spaces are as follows. The first Korn’s inequality
[26] holds on $\boldsymbol{V}$, that is, there exists a constant
$C_{K}=C_{K}(\Omega,\Gamma_{\boldsymbol{u},D})>0$ such that

$|\\!|\boldsymbol{u}|\\!|_{H^{1}(\Omega)}\leq
C_{K}|\\!|\epsilon(\boldsymbol{u})|\\!|_{L^{2}(\Omega)}\quad\forall\boldsymbol{u}\in\boldsymbol{V}.$

Furthermore, the following inf-sup condition [26] holds. There exists a
constant $\beta_{0}>0$ depending only on $\Omega$ and
$\Gamma_{\boldsymbol{u},D}$ such that

$\displaystyle\sup_{\boldsymbol{v}\in\boldsymbol{V}}\frac{(\operatorname{div}\boldsymbol{v},\eta)}{|\\!|\boldsymbol{v}|\\!|_{H^{1}(\Omega)}}\geq\beta_{0}|\\!|\eta|\\!|_{L^{2}(\Omega)}\quad\forall\eta\in
W.$ (3)

######  Assumption 1.

We assume that $\boldsymbol{u}_{0}\in\boldsymbol{H}^{1}(\Omega)$,
$\boldsymbol{f}\in\boldsymbol{L}^{2}(\Omega)$,
$\boldsymbol{h}\in\boldsymbol{L}^{2}(\Gamma_{\boldsymbol{u},N})$,
$\vec{p}_{0}\in[L^{2}(\Omega)]^{N}$, $\vec{g}\in[L^{2}(\Omega)]^{N}$ and
$\vec{l}\in[L^{2}({\Gamma_{\vec{p},N}})]^{N}$. We also assume that $\mu>0$,
$\lambda>0$, $K_{i}>0,i=1,\cdots,N$, $\beta_{i,j}>0,1\leq i,j\leq N,j\neq i$,
$T>0$.

For simplicity, we will assume Assumption 1 holds in the rest of our paper.
For ease of presentation, we also assume that $\boldsymbol{f}$,
$\boldsymbol{h}$, $\vec{g}$, $\vec{l}$ are independent of $t$. Given $T>0$,
$(\boldsymbol{u},\xi,\vec{p})\in\boldsymbol{V}\times W\times M$ with

$\displaystyle\boldsymbol{u}\in L^{\infty}(0,T;\boldsymbol{V}),\,\xi\in
L^{\infty}(0,T;W),\,\vec{p}\in[L^{\infty}(0,T;M)]^{N},$
$\displaystyle\dot{\vec{p}}\in[L^{2}(0,T;M^{\prime})]^{N},\,\dot{\xi}\in
L^{2}(0,T;M^{\prime}),$

is called a weak solution of problem (2a)-(2c), if there holds

$\displaystyle
2\mu(\epsilon(\boldsymbol{u}),\epsilon(\boldsymbol{v}))-(\xi,\operatorname{div}\boldsymbol{v})=(\boldsymbol{f},\boldsymbol{v})+\langle\boldsymbol{h},\boldsymbol{v}\rangle_{\Gamma_{\boldsymbol{u},N}},$
(4a)
$\displaystyle(\operatorname{div}\boldsymbol{u},\eta)+\frac{1}{\lambda}(\xi,\eta)-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}\vec{p},\eta)=0,$
(4b)
$\displaystyle\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)\dot{\vec{p}},\vec{q}\right)$
$\displaystyle-\frac{1}{\lambda}\left(\vec{\alpha}\dot{\xi},\vec{q}\right)+\left(K\nabla\vec{p},\nabla\vec{q}\right)+\left(B\vec{p},\vec{q}\right)=\left(\vec{g},\vec{q}\right)+\left\langle\vec{l},\vec{q}\right\rangle_{\Gamma_{\vec{p},N}},$
(4c)

for almost every $t\in[0,T]$.

##  3 The algorithms

In this section, based on the above total-pressure reformulation, we will
present the coupled algorithm and the iterative decoupled algorithm. We will
start with the introduction of the mixed finite element spaces for
discretizing the variational problem (4). We will also provide the
approximation properties of the corresponding projection operators.

###  3.1 Finite element spaces

Let $\mathcal{T}_{h}$ be a shape-regular and conforming triangulation of the
domain $\Omega$. It is common in the literature to set
$h:=\max_{K\in\mathcal{T}_{h}}h_{K}$. The solid displacement and total
pressure are approximated by the lowest-order Taylor-Hood elements,
i.e.,$(\boldsymbol{P}_{2},P_{1})$ Lagrange finite elements, while the network
pressures are discretized by using the $P_{1}$ Lagrange finite elements [27].
The finite element spaces on $\mathcal{T}_{h}$ are as follows.

$\displaystyle\boldsymbol{V}_{h}:=\\{\boldsymbol{v}_{h}\in\boldsymbol{V}\cap[C^{0}(\overline{\Omega})]^{2};\boldsymbol{v}_{h}|_{\Gamma_{\boldsymbol{u},D}}=\boldsymbol{0},\boldsymbol{v}_{h}|_{K}\in[P_{2}(K)]^{2},\forall
K\in T_{h}\\},$ $\displaystyle W_{h}:=\\{\eta_{h}\in W\cap
C^{0}(\overline{\Omega});\eta_{h}|_{K}\in P_{1}(K),\forall K\in T_{h}\\},$
$\displaystyle M_{h}:=\\{\vec{q}_{h}\in
M\cap[C^{0}(\overline{\Omega})]^{N};\vec{q}_{h}|_{\Gamma_{\vec{p},D}}=\vec{0},\vec{q}_{h}|_{K}\in[P_{1}(K)]^{N},\forall
K\in T_{h}\\}.$

We note that $\boldsymbol{V}_{h}\times W_{h}$ is a stable Stokes pair, i.e.,
the following discrete inf-sup condition holds. There exists a constant
$\beta_{0}^{*}>0$, independent of $h$, such that

$\displaystyle\sup_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\frac{(\operatorname{div}\boldsymbol{v}_{h},\eta_{h})}{|\\!|\boldsymbol{v}_{h}|\\!|_{H^{1}(\Omega)}}\geq\beta_{0}^{*}|\\!|\eta_{h}|\\!|_{L^{2}(\Omega)}\quad\forall\eta_{h}\in
W_{h}.$ (5)

###  3.2 The coupled algorithm

The time discretization is taken as an equidistant partition
$0=t^{0}<t^{1}<\cdots<t^{L}=T$. For simplicity, we define
$\boldsymbol{u}^{n}=\boldsymbol{u}(t^{n}),\xi^{n}=\xi(t^{n}),\vec{p}^{n}=\vec{p}(t^{n})$.
Suppose that initial values
$(\boldsymbol{u}_{h}^{0},\xi_{h}^{0},\vec{p}_{h}^{0})\in\boldsymbol{V}_{h}\times
W_{h}\times M_{h}$ are provided, we apply a backward Euler scheme for the time
discretization. The following algorithm is considered.

For all $n=1,2,\cdots,L$, given
$(\boldsymbol{u}_{h}^{n-1},\xi_{h}^{n-1},\vec{p}_{h}^{n-1})\in\boldsymbol{V}_{h}\times
W_{h}\times M_{h}$, find
$(\boldsymbol{u}_{h}^{n},\xi_{h}^{n},\vec{p}_{h}^{n})\in\boldsymbol{V}_{h}\times
W_{h}\times M_{h}$, for all
$(\boldsymbol{v}_{h},\eta_{h},\vec{q}_{h})\in\boldsymbol{V}_{h}\times
W_{h}\times M_{h}$ such that

$\displaystyle
2\mu(\epsilon(\boldsymbol{u}_{h}^{n}),\epsilon(\boldsymbol{v}_{h}))$
$\displaystyle-(\xi_{h}^{n},\operatorname{div}\boldsymbol{v}_{h})=(\boldsymbol{f},\boldsymbol{v}_{h})+\langle\boldsymbol{h},\boldsymbol{v}_{h}\rangle_{\Gamma_{\boldsymbol{u},N}},$
(6a) $\displaystyle(\operatorname{div}\boldsymbol{u}_{h}^{n},\eta_{h})$
$\displaystyle+\frac{1}{\lambda}(\xi_{h}^{n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}\vec{p}_{h}^{n},\eta_{h})=0,$
(6b)
$\displaystyle\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)\frac{\vec{p}_{h}^{n}-\vec{p}_{h}^{n-1}}{\Delta
t},\vec{q}_{h}\right)$
$\displaystyle-\frac{1}{\lambda}\left(\vec{\alpha}\frac{\xi_{h}^{n}-\xi_{h}^{n-1}}{\Delta
t},\vec{q}_{h}\right)+\left(K\nabla\vec{p}_{h}^{n},\nabla\vec{q}_{h}\right)+\left(B\vec{p}_{h}^{n},\vec{q}_{h}\right)$
(6c)
$\displaystyle=\left(\vec{g},\vec{q}_{h}\right)+\left\langle\vec{l},\vec{q}_{h}\right\rangle_{\Gamma_{\vec{p},N}}.$

###  3.3 The iterative decoupled algorithm

Now, we introduce the iterative decoupled algorithm. In each time step of the
algorithm, the previous iteration is used as the initial values, and then
solve the reaction-diffusion equation for $\vec{p}$ and the generalized Stokes
equations for $\boldsymbol{u}$ and $\xi$ alternately. Let us define a sequence
$\\{(\boldsymbol{u}_{h}^{n,k},\xi_{h}^{n,k},\vec{p}_{h}^{n,k})\\}$ with
$k\in\mathbb{N}^{*}$ being the iteration index. After initialization, i.e.,
$\boldsymbol{u}_{h}^{n,0}=\boldsymbol{u}_{h}^{n-1},\xi_{h}^{n,0}=\xi_{h}^{n-1},\vec{p}_{h}^{n,0}=\vec{p}_{h}^{n-1}$,
each iteration consists of the following two steps. For a fixed time-step
index $n$, the $k$-th iteration can be expressed as follows.
Step 1. Given $\xi_{h}^{n,k-1}\in W_{h}$, find $\vec{p}_{h}^{n,k}\in M_{h}$,
for all $\vec{q}_{h}\in M_{h}$ such that

$\displaystyle\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)\vec{p}_{h}^{n,k},\vec{q}_{h}\right)+\Delta
t\left(K\nabla\vec{p}_{h}^{n,k},\nabla\vec{q}_{h}\right)+\Delta
t\left(B\vec{p}_{h}^{n,k},\vec{q}_{h}\right)$ (7a)
$\displaystyle=\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)\vec{p}_{h}^{n-1},\vec{q}_{h}\right)+\frac{1}{\lambda}\left(\vec{\alpha}\left(\xi_{h}^{n,k-1}-\xi_{h}^{n-1}\right),\vec{q}_{h}\right)+\Delta
t\left(\vec{g},\vec{q}_{h}\right)+\Delta
t\left\langle\vec{l},\vec{q}_{h}\right\rangle_{\Gamma_{\vec{p},N}}.$ Step 2.
Given $\vec{p}_{h}^{n,k}\in M_{h}$, find
$(\boldsymbol{u}_{h}^{n,k},\xi_{h}^{n,k})\in\boldsymbol{V}_{h}\times W_{h}$,
for all $(\boldsymbol{v}_{h},\eta_{h})\in\boldsymbol{V}_{h}\times W_{h}$ such
that $\displaystyle
2\mu(\epsilon(\boldsymbol{u}_{h}^{n,k}),\epsilon(\boldsymbol{v}_{h}))-(\xi_{h}^{n,k},\operatorname{div}\boldsymbol{v}_{h})=(\boldsymbol{f},\boldsymbol{v}_{h})+\langle\boldsymbol{h},\boldsymbol{v}_{h}\rangle_{\Gamma_{\boldsymbol{u},N}},$
(7b)
$\displaystyle(\operatorname{div}\boldsymbol{u}_{h}^{n,k},\eta_{h})+\frac{1}{\lambda}(\xi_{h}^{n,k},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}\vec{p}_{h}^{n,k},\eta_{h})=0.$
(7c)

##  4 Convergence analysis

For a detailed error analysis of the coupled algorithm, we direct the readers
to Appendix A. In this section, our objective is to establish that as $k$
approaches positive infinity, the sequence
${(\boldsymbol{u}_{h}^{n,k},\xi_{h}^{n,k},\vec{p}_{h}^{n,k})}$ generated
through the iterative decoupled algorithm converges towards the solution
${(\boldsymbol{u}_{h}^{n},\xi_{h}^{n},\vec{p}_{h}^{n})}$ obtained via the
coupled algorithm. To substantiate these convergence results, we present a
crucial lemma [26].

######  Lemma 1.

For all $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$, the following inequality
holds

$\displaystyle|\\!|\operatorname{div}\boldsymbol{u}_{h}|\\!|_{L^{2}(\Omega)}\leq\sqrt{2}|\\!|\epsilon(\boldsymbol{u}_{h})|\\!|_{L^{2}(\Omega)}.$
(8)

For a complete theoretical analysis of the decoupled algorithm, we firstly
assume that $\min_{1\leq i\leq N}c_{i}>0$, and we have the following theorem.

######  Theorem 2.

Let ${(\boldsymbol{u}_{h}^{n},\xi_{h}^{n},\vec{p}_{h}^{n})}$ and
${(\boldsymbol{u}_{h}^{n,k},\xi_{h}^{n,k},\vec{p}_{h}^{n,k})}$ be the
solutions of problem (6a)-(6c) and problem (7a)-(7c), respectively. Let
$e_{\boldsymbol{u}}^{k}=\boldsymbol{u}_{h}^{n,k}-\boldsymbol{u}_{h}^{n},e_{\xi}^{k}=\xi_{h}^{n,k}-\xi_{h}^{n},e_{\vec{p}}^{k}=\vec{p}_{h}^{n,k}-\vec{p}_{h}^{n}$
denote the errors between the iterative solution in the $k$-th step and the
solution of the coupled algorithm. Then, for all $k\geq 1$, it holds that

$\displaystyle|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}\leq
C^{*}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)},$

where
$C^{*}=\frac{\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}{\delta+\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}$
is a positive constant less than 1,
$|\\!|\vec{\alpha}|\\!|=\left(\sum_{i=1}^{N}\alpha_{i}^{2}\right)^{\frac{1}{2}}$,
$\delta=\min_{1\leq i\leq N}c_{i}$. Moreover, it holds that

$\displaystyle\lim_{k\rightarrow+\infty}|\\!|\nabla
e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}$
$\displaystyle=\lim_{k\rightarrow+\infty}|\\!|e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}=\lim_{k\rightarrow+\infty}|\\!|e_{\boldsymbol{u}}^{k}|\\!|_{H^{1}(\Omega)}=\lim_{k\rightarrow+\infty}|\\!|e_{\boldsymbol{u}}^{k}|\\!|_{L^{2}(\Omega)}=0.$

Proof. Subtracting (6c), (6a) and (6b) from (7a), (7b) and (7c), and taking
$\boldsymbol{v}_{h}=e_{\boldsymbol{u}}^{k},\eta_{h}=e_{\xi}^{k},\vec{q}_{h}=e_{\vec{p}}^{k}$,
respectively, we have

$\displaystyle
2\mu(\epsilon(e_{\boldsymbol{u}}^{k}),\epsilon(e_{\boldsymbol{u}}^{k}))-(e_{\xi}^{k},\operatorname{div}e_{\boldsymbol{u}}^{k})=0,$
(9)
$\displaystyle(\operatorname{div}e_{\boldsymbol{u}}^{k},e_{\xi}^{k})+\frac{1}{\lambda}(e_{\xi}^{k},e_{\xi}^{k})=\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k},e_{\xi}^{k}),$
(10)
$\displaystyle\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)e_{\vec{p}}^{k},e_{\vec{p}}^{k}\right)$
$\displaystyle+\Delta t\left(K\nabla e_{\vec{p}}^{k},\nabla
e_{\vec{p}}^{k}\right)+\Delta
t\left(Be_{\vec{p}}^{k},e_{\vec{p}}^{k}\right)=\frac{1}{\lambda}\left(\vec{\alpha}e_{\xi}^{k-1},e_{\vec{p}}^{k}\right).$
(11)

Using (11), there holds

$\displaystyle\left(Se_{\vec{p}}^{k},e_{\vec{p}}^{k}\right)+\frac{1}{\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}\leq\frac{1}{\lambda}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)},$
(12)

where we have used the fact $\Delta t\left(K\nabla e_{\vec{p}}^{k},\nabla
e_{\vec{p}}^{k}\right)\geq 0,\Delta
t\left(Be_{\vec{p}}^{k},e_{\vec{p}}^{k}\right)\geq 0.$
Note that

$\displaystyle\frac{|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}}{|\\!|\vec{\alpha}|\\!|^{2}}\leq|\\!|e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}.$
(13)

Plugging (13) into (12), we have

$\displaystyle\delta\frac{|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}}{|\\!|\vec{\alpha}|\\!|^{2}}+\frac{1}{\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq\delta|\\!|e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq\left(Se_{\vec{p}}^{k},e_{\vec{p}}^{k}\right)+\frac{1}{\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq\frac{1}{\lambda}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)},$

where $\delta=\min_{1\leq i\leq N}c_{i}$.
In the above inequality, dividing both sides by
$|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}$ leads to

$\displaystyle|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}\leq\frac{\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}{\delta+\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}.$
(14)

Ignoring the nonnegative term
$\frac{1}{\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}$
in (12) and applying (14), we obtain

$\displaystyle\delta|\\!|e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}\leq\left(Se_{\vec{p}}^{k},e_{\vec{p}}^{k}\right)$
$\displaystyle\leq\frac{1}{\lambda}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}$
$\displaystyle\leq\frac{1}{\lambda}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}\frac{\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}{\delta+\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}.$

By taking square roots of both sides, we have

$\displaystyle|\\!|e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}\leq\left(\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda\delta(\lambda\delta+|\\!|\vec{\alpha}|\\!|^{2})}\right)^{\frac{1}{2}}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}.$
(15)

Similarly, ignoring the nonnegative term
$\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)e_{\vec{p}}^{k},e_{\vec{p}}^{k}\right)$
and $\Delta t\left(Be_{\vec{p}}^{k},e_{\vec{p}}^{k}\right)$ in (11), and using
(14), we derive that

$\displaystyle\Delta t\zeta|\\!|\nabla
e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}^{2}\leq\Delta t\left(K\nabla
e_{\vec{p}}^{k},\nabla e_{\vec{p}}^{k}\right)$
$\displaystyle\leq\frac{1}{\lambda}\left(\vec{\alpha}e_{\xi}^{k-1},e_{\vec{p}}^{k}\right)$
$\displaystyle\leq\frac{1}{\lambda}|\\!|e_{\xi}^{k-1}|\\!|\frac{\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}{\delta+\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)},$

where $\zeta=\min_{1\leq i\leq N}K_{i}$.
By taking the square roots of both sides, we obtain

$\displaystyle|\\!|\nabla
e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}\leq\left(\frac{\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}{(\lambda\delta+|\\!|\vec{\alpha}|\\!|^{2})\Delta
t\zeta}\right)^{\frac{1}{2}}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}.$ (16)

Summing up (9) and (10), there holds

$\displaystyle
2\mu|\\!|\epsilon(e_{\boldsymbol{u}}^{k})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{\lambda}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k},e_{\xi}^{k})\leq\frac{1}{\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}.$
(17)

Ignoring the first positive term
$2\mu|\\!|\epsilon(e_{\boldsymbol{u}}^{k})|\\!|_{L^{2}(\Omega)}^{2}$ and using
(14), we deduce that

$\displaystyle|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}\leq|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}\leq\frac{\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}{\delta+\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}.$
(18)

Applying (8) to (9) , it holds that

$\displaystyle
2\mu|\\!|\epsilon(e_{\boldsymbol{u}}^{k})|\\!|_{L^{2}(\Omega)}^{2}=(e_{\xi}^{k},\operatorname{div}e_{\boldsymbol{u}}^{k})$
$\displaystyle\leq|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}|\\!|\operatorname{div}e_{\boldsymbol{u}}^{k}|\\!|_{L^{2}(\Omega)}\leq\sqrt{2}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}|\\!|\epsilon(e_{\boldsymbol{u}}^{k})|\\!|_{L^{2}(\Omega)}.$

Using the first Korn’s inequality [26], it can be seen that

$\displaystyle|\\!|e_{\boldsymbol{u}}^{k}|\\!|_{L^{2}(\Omega)}\leq|\\!|e_{\boldsymbol{u}}^{k}|\\!|_{H^{1}(\Omega)}\leq
C_{K}|\\!|\epsilon(e_{\boldsymbol{u}}^{k})|\\!|_{L^{2}(\Omega)}\leq\frac{C_{K}}{\sqrt{2}\mu}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}.$
(19)

According to (18) and the inequality
$0<\frac{\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}{\delta+\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda}}<1$,
then

$\displaystyle\lim_{k\rightarrow+\infty}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}=0.$

Following from (15), (16) and (19), we have

$\displaystyle\lim_{k\rightarrow+\infty}|\\!|\nabla
e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}=\lim_{k\rightarrow+\infty}|\\!|e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}=\lim_{k\rightarrow+\infty}|\\!|e_{\boldsymbol{u}}^{k}|\\!|_{H^{1}(\Omega)}=\lim_{k\rightarrow+\infty}|\\!|e_{\boldsymbol{u}}^{k}|\\!|_{L^{2}(\Omega)}=0.$
(20)

$\hfill\qed$

For $\min_{1\leq i\leq N}c_{i}=0$, the main conclusions are presented in the
following remark.

######  Remark 1.

If $\min_{1\leq i\leq N}c_{i}=0$, we will show that the arguments (20) also
hold true. It means that the iterative decoupled algorithm also converges when
$\min_{1\leq i\leq N}c_{i}$ degenerates to $0$.

According to (12), we obtain

$\displaystyle|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}\leq|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}.$
(21)

Thanks to (11), (17), (21), and the Poincaré inequality, we can derive the
following:

$\displaystyle|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}$
$\displaystyle\leq|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)},$ (22)
$\displaystyle
2\mu|\\!|\epsilon(e_{\boldsymbol{u}}^{k})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{\lambda}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq\frac{1}{\lambda}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)},$
(23) $\displaystyle|\\!|e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}\leq
C_{P}|\\!|\nabla e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}$ $\displaystyle\leq
C_{P}\left(\frac{1}{\lambda\Delta
t\zeta}\right)^{\frac{1}{2}}|\\!|e_{\xi}^{k-1}|\\!|_{L^{2}(\Omega)},$ (24)

where $C_{P}$ is the constant in the Poincaré inequality.

We are going to use the method of contradiction to show that the limit of
$\\{|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}\\}$ is 0. If not, let us assume

$\lim_{k\rightarrow+\infty}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}>0.$

From (23), we know that

$\lim_{k\rightarrow+\infty}|\\!|\epsilon(e_{\boldsymbol{u}}^{k})|\\!|_{L^{2}(\Omega)}=0.$

Applying the discrete inf-sup condition, we see that

$\displaystyle\beta_{0}^{*}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}\leq\sup_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\frac{|(e_{\xi}^{k},\operatorname{div}\boldsymbol{v}_{h})|}{|\\!|\boldsymbol{v}_{h}|\\!|_{H^{1}(\Omega)}}=\sup_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\frac{2\mu|(\epsilon(e_{\boldsymbol{u}}^{k}),\epsilon(\boldsymbol{v}_{h}))|}{|\\!|\boldsymbol{v}_{h}|\\!|_{H^{1}(\Omega)}}\leq
2\mu|\\!|\epsilon(e_{\boldsymbol{u}}^{k})|\\!|_{L^{2}(\Omega)}.$ (25)

Letting $k\rightarrow\infty$ in (25), we derive that

$\lim_{k\rightarrow+\infty}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}\leq 0,$

which is a contradiction. Therefore

$\displaystyle\lim_{k\rightarrow+\infty}|\\!|e_{\xi}^{k}|\\!|_{L^{2}(\Omega)}=0.$
(26)

Combined with (19), (24) and (26), it follows that

$\displaystyle\lim_{k\rightarrow+\infty}|\\!|\nabla
e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}=\lim_{k\rightarrow+\infty}|\\!|e_{\vec{p}}^{k}|\\!|_{L^{2}(\Omega)}=\lim_{k\rightarrow+\infty}|\\!|e_{\boldsymbol{u}}^{k}|\\!|_{H^{1}(\Omega)}=\lim_{k\rightarrow+\infty}|\\!|e_{\boldsymbol{u}}^{k}|\\!|_{L^{2}(\Omega)}=0.$

##  5 Numerical experiments

Let $\Omega=[0,1]\times[0,1]$, $\Gamma_{1}=\\{(1,y);0\leq y\leq 1\\}$,
$\Gamma_{2}=\\{(x,0);0\leq x\leq 1\\}$, $\Gamma_{3}=\\{(0,y);0\leq y\leq
1\\}$, and $\Gamma_{4}=\\{(x,1);0\leq x\leq 1\\}$. The terminal time is
$T=0.01$ seconds. In our numerical experiments, we impose a pure Dirichlet
boundary condition on this domain. We initialize the mesh with a size of
$h=\frac{1}{8}$ and subsequently carry out four rounds of mesh refinement,
which involves connecting the midpoints of each triangle. We use the variable
$iters$ to represent the number of iterations utilized within the iterative
decoupled algorithm. In these tests, we intentionally employ relatively large
time step sizes to effectively showcase the efficiency of the iterative
decoupled algorithm. Specifically, we choose the values of $\Delta t$ and
$iters$ such that the overall computational expense is approximately equal to
that of the coupled algorithm. More precisely, we set the values of $\Delta t$
and $iters$ such that $T/\Delta t\times iters=50$. At the conclusion of the
simulations, we calculate the $H^{1}$ and $L^{2}$ norms at the terminal time
$T$. All the tests discussed in this paper are implemented using the open-
source software FreeFEM++ [28].

The exact solutions of (2a)-(2c) are given as follows.

$\displaystyle\boldsymbol{u}$
$\displaystyle=\left[\begin{array}[]{c}\text{sin}(2\pi y)(-1+\text{cos}(2\pi
x))+\frac{1}{\mu+\lambda}\text{sin}(\pi x)\text{sin}(\pi y)\\\ \text{sin}(2\pi
x)(1-\text{cos}(2\pi y))+\frac{1}{\mu+\lambda}\text{sin}(\pi x)\text{sin}(\pi
y)\end{array}\right]\text{sin}(t),$ $\displaystyle p_{1}$
$\displaystyle=-\text{sin}(\pi x)\text{sin}(\pi y)\text{cos}(t),$
$\displaystyle p_{2}$ $\displaystyle=-2\text{sin}(\pi x)\text{sin}(\pi
y)\text{cos}(t).$

###  5.1 Tests for the parameter $\nu$

In this subsection, we test the performance of the algorithms under different
settings of the Poisson ratio $\nu$. For other physical parameters, we set

$\displaystyle
E=1,\alpha_{1}=\alpha_{2}=1,c_{1}=c_{2}=1,K_{1}=K_{2}=1,\beta_{12}=\beta_{21}=1.$

In Table 1-Table 3, we present the numerical results of the coupled algorithm
and the iterative decoupled algorithm for the case of Poisson’s ratio $v=0.3$.
Table 2 displays the results obtained from the iterative decoupled algorithm
with 10 iterations, while Table 3 shows the results obtained with 20
iterations. By observing the numerical results of Table 1-Table 3, it can be
found that the energy-norm errors of all variables converge to the optimal
orders. Furthermore, the results in Tables 3 indicate that increasing the
number of iterations leads to enhanced accuracy in the iterative decoupled
algorithm.

Table 1: Convergence rate of the coupled algorithm. $\nu=0.3$ and $\Delta t=2\times 10^{-4}s$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 1.230e-03 & 1.768e-02 |  | 3.652e-02 & 1.083e+00 |
16 | 3.013e-04 & 4.032e-03 | 2.03 & 2.13 | 9.105e-03 & 5.506e-01 | 2.00 & 0.976
32 | 7.536e-05 & 9.421e-04 | 2.00 & 2.10 | 2.269e-03 & 2.760e-01 | 2.00 & 0.996
64 | 1.890e-05 & 2.257e-04 | 2.00 & 2.06 | 5.670e-04 & 1.381e-01 | 2.00 & 0.999
128 | 4.766e-06 & 5.523e-05 | 1.99 & 2.03 | 1.423e-04 & 6.908e-02 | 1.99 & 1.000
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.432e-02 & 3.581e-01 |  | 2.851e-02 & 7.161e-01 |
16 | 3.681e-03 & 1.816e-01 | 1.96 & 0.98 | 7.342e-03 & 3.633e-01 | 1.96 & 0.98
32 | 9.354e-04 & 9.134e-02 | 1.98 & 0.99 | 1.868e-03 & 1.827e-01 | 1.97 & 0.99
64 | 2.403e-04 & 4.576e-02 | 1.96 & 1.00 | 4.809e-04 & 9.153e-02 | 1.96 & 1.00
128 | 6.586e-05 & 2.290e-02 | 1.87 & 1.00 | 1.327e-04 & 4.579e-02 | 1.86 & 1.00

Table 2: Convergence rate of the iterative decoupled algorithm. $\nu=0.3$, $\Delta t=2\times 10^{-3}s$, and $iters=10$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 1.229e-03 & 1.765e-02 |  | 3.667e-02 & 1.084e+00 |
16 | 3.011e-04 & 4.012e-03 | 2.03 & 2.14 | 9.146e-03 & 5.507e-01 | 2.00 & 0.977
32 | 7.577e-05 & 9.367e-04 | 1.99 & 2.10 | 2.283e-03 & 2.760e-01 | 2.00 & 0.996
64 | 1.961e-05 & 2.260e-04 | 1.95 & 2.05 | 5.738e-04 & 1.381e-01 | 1.99 & 0.999
128 | 5.921e-06 & 5.813e-05 | 1.73 & 1.96 | 1.478e-04 & 6.908e-02 | 1.96 & 1.000
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.200e-02 & 3.609e-01 |  | 2.625e-02 & 7.203e-01 |
16 | 3.036e-03 & 1.827e-01 | 1.98 & 0.98 | 6.662e-03 & 3.652e-01 | 1.98 & 0.98
32 | 7.626e-04 & 9.182e-02 | 1.99 & 0.99 | 1.678e-03 & 1.836e-01 | 1.99 & 0.99
64 | 1.920e-04 & 4.599e-02 | 1.99 & 1.00 | 4.239e-04 & 9.198e-02 | 1.98 & 1.00
128 | 4.950e-05 & 2.301e-02 | 1.96 & 1.00 | 1.100e-04 & 4.602e-02 | 1.95 & 1.00
Table 3: Convergence rate of the iterative decoupled algorithm. $\nu=0.3$, $\Delta t=4\times 10^{-3}s$, and $iters=20$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 1.184e-03 & 1.734e-02 |  | 3.656e-02 & 1.084e+00 |
16 | 2.879e-04 & 3.877e-03 | 2.04 & 2.16 | 9.114e-03 & 5.508e-01 | 2.00 & 0.977
32 | 7.230e-05 & 8.955e-04 | 1.99 & 2.11 | 2.279e-03 & 2.760e-01 | 2.00 & 0.997
64 | 1.858e-05 & 2.150e-04 | 1.96 & 2.06 | 5.782e-04 & 1.381e-01 | 1.98 & 0.999
128 | 5.148e-06 & 5.469e-05 | 1.85 & 1.97 | 1.542e-04 & 6.908e-02 | 1.91 & 1.000
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.192e-02 & 3.613e-01 |  | 2.613e-02 & 7.206e-01 |
16 | 3.004e-03 & 1.828e-01 | 1.99 & 0.98 | 6.616e-03 & 3.652e-01 | 1.98 & 0.98
32 | 7.554e-04 & 9.183e-02 | 1.99 & 0.99 | 1.668e-03 & 1.836e-01 | 1.99 & 0.99
64 | 1.919e-04 & 4.600e-02 | 1.98 & 1.00 | 4.248e-04 & 9.199e-02 | 1.97 & 1.00
128 | 5.088e-05 & 2.301e-02 | 1.91 & 1.00 | 1.136e-04 & 4.602e-02 | 1.90 & 1.00

Table 1-Table 3 is for the case that the poroelastic material is compressible.
In Table 4 through Table 6, we keep the physical parameters constant while
setting the Poisson ratio to $\nu=0.49999$. With such a Poisson ratio, the
poroelastic material exhibits nearly incompressible behavior. Both Table 4 and
Table 5 through Table 6 present results obtained from the coupled algorithm
and the iterative decoupled algorithm, respectively, each with different
iteration counts. It is noteworthy that as the Poisson ratio $\nu$ approaches
0.5, which indicates a state of near-incompressibility, both algorithms
demonstrate improved numerical error rates compared to when $\nu=0.3$. Upon
examining Table 4 through Table 6, it becomes evident that in both algorithms,
the $L^{2}$ norm errors and energy norm errors of each variable exhibit
optimal convergence rates. A comparative analysis of the numerical outcomes
presented in Table 4 and Table 6 underscores the superior error accuracy
achieved by our proposed iterative decoupled algorithm.

Table 4: Convergence rate of the coupled algorithm. $\nu=0.49999$ and $\Delta t=2\times 10^{-4}s$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 4.109e-04 & 1.668e-02 |  | 3.945e-02 & 1.090e+00 |
16 | 3.574e-05 & 3.081e-03 | 3.52 & 2.44 | 9.785e-03 & 5.607e-01 | 2.01 & 0.959
32 | 3.262e-06 & 5.696e-04 | 3.45 & 2.44 | 2.431e-03 & 2.779e-01 | 2.01 & 1.013
64 | 2.975e-07 & 1.068e-04 | 3.45 & 2.42 | 6.065e-04 & 1.386e-01 | 2.00 & 1.004
128 | 2.776e-08 & 2.076e-05 | 3.42 & 2.36 | 1.521e-04 & 6.919e-02 | 2.00 & 1.002
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.469e-02 & 3.599e-01 |  | 2.911e-02 & 7.197e-01 |
16 | 3.743e-03 & 1.825e-01 | 1.97 & 0.98 | 7.414e-03 & 3.651e-01 | 1.97 & 0.98
32 | 9.426e-04 & 9.180e-02 | 1.99 & 0.99 | 1.867e-03 & 1.836e-01 | 1.99 & 0.99
64 | 2.365e-04 & 4.599e-02 | 1.99 & 1.00 | 4.684e-04 & 9.198e-02 | 1.99 & 1.00
128 | 5.945e-05 & 2.301e-02 | 1.99 & 1.00 | 1.177e-04 & 4.602e-02 | 1.99 & 1.00
Table 5: Convergence rate of the iterative decoupled algorithm. $\nu=0.49999$, $\Delta t=2\times 10^{-3}s$, and $iters=10$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 4.109e-04 & 1.668e-02 |  | 3.942e-02 & 1.090e+00 |
16 | 3.574e-05 & 3.081e-03 | 3.52 & 2.44 | 9.783e-03 & 5.607e-01 | 2.01 & 0.959
32 | 3.262e-06 & 5.696e-04 | 3.45 & 2.44 | 2.437e-03 & 2.779e-01 | 2.01 & 1.013
64 | 2.975e-07 & 1.068e-04 | 3.45 & 2.42 | 6.144e-04 & 1.386e-01 | 1.99 & 1.004
128 | 2.776e-08 & 2.076e-05 | 3.42 & 2.36 | 1.607e-04 & 6.919e-02 | 1.94 & 1.002
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.469e-02 & 3.599e-01 |  | 2.909e-02 & 7.198e-01 |
16 | 3.743e-03 & 1.825e-01 | 1.97 & 0.98 | 7.413e-03 & 3.651e-01 | 1.97 & 0.98
32 | 9.449e-04 & 9.180e-02 | 1.99 & 0.99 | 1.872e-03 & 1.836e-01 | 1.99 & 0.99
64 | 2.396e-04 & 4.599e-02 | 1.98 & 1.00 | 4.746e-04 & 9.198e-02 | 1.98 & 1.00
128 | 6.281e-05 & 2.301e-02 | 1.93 & 1.00 | 1.244e-04 & 4.602e-02 | 1.93 & 1.00

Table 6: Convergence rate of the iterative decoupled algorithm. $\nu=0.49999$, $\Delta t=4\times 10^{-3}s$, and $iters=20$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 4.096e-04 & 1.657e-02 |  | 3.927e-02 & 1.090e+00 |
16 | 3.552e-05 & 3.046e-03 | 3.53 & 2.44 | 9.738e-03 & 5.606e-01 | 2.01 & 0.959
32 | 3.224e-06 & 5.575e-04 | 3.46 & 2.45 | 2.429e-03 & 2.779e-01 | 2.00 & 1.012
64 | 2.911e-07 & 1.027e-04 | 3.47 & 2.44 | 6.168e-04 & 1.386e-01 | 1.98 & 1.004
128 | 2.669e-08 & 1.944e-05 | 3.45 & 2.40 | 1.657e-04 & 6.918e-02 | 1.90 & 1.002
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.462e-02 & 3.600e-01 |  | 2.903e-02 & 7.199e-01 |
16 | 3.723e-03 & 1.826e-01 | 1.97 & 0.98 | 7.388e-03 & 3.651e-01 | 1.97 & 0.98
32 | 9.410e-04 & 9.180e-02 | 1.98 & 0.99 | 1.867e-03 & 1.836e-01 | 1.98 & 0.99
64 | 2.403e-04 & 4.599e-02 | 1.97 & 1.00 | 4.768e-04 & 9.198e-02 | 1.97 & 1.00
128 | 6.476e-05 & 2.301e-02 | 1.89 & 1.00 | 1.285e-04 & 4.602e-02 | 1.89 & 1.00

###  5.2 Tests for the parameter $K_{i}$

In this subsection, we conduct experiments to assess the performance of the
iterative decoupled algorithm by varying the values of the parameter $K_{i}$.
It’s worth noting that we previously examined the case where $K_{1}=K_{2}=1$
in our earlier tests, but now we set $K_{1}=K_{2}=10^{-6}$.

For the other physical parameters, we maintain the following values:

$\displaystyle
E=1,\nu=0.3,\alpha_{1}=\alpha_{2}=1,c_{1}=c_{2}=1,\beta_{12}=\beta_{21}=1.$

The results presented in Table 8 through Table 9 illustrate that increasing
the number of iterations results in improved accuracy when using the iterative
decoupled algorithm. Upon analyzing the numerical results in Table 3 and Table
9, it becomes apparent that the energy-norm errors for all variables converge
to their optimal orders. This suggests that the accuracy of the iterative
decoupled algorithm remains robust and is minimally affected by variations in
hydraulic conductivity.

Table 7: Convergence rate of the coupled algorithm. $K_{1}=K_{2}=10^{-6}$ and $\Delta t=2\times 10^{-4}s$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 1.021e-03 & 1.674e-02 |  | 3.674e-02 & 1.092e+00 |
16 | 2.283e-04 & 3.321e-03 | 2.16 & 2.33 | 9.210e-03 & 5.561e-01 | 2.00 & 0.973
32 | 5.444e-05 & 6.833e-04 | 2.07 & 2.28 | 2.326e-03 & 2.772e-01 | 1.99 & 1.004
64 | 1.334e-05 & 1.477e-04 | 2.03 & 2.21 | 5.852e-04 & 1.384e-01 | 1.99 & 1.002
128 | 3.325e-06 & 3.369e-05 | 2.00 & 2.13 | 1.475e-04 & 6.914e-02 | 1.99 & 1.001
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.210e-02 & 3.662e-01 |  | 2.584e-02 & 7.288e-01 |
16 | 3.037e-03 & 1.863e-01 | 1.99 & 0.98 | 6.453e-03 & 3.690e-01 | 2.00 & 0.98
32 | 7.711e-04 & 9.268e-02 | 1.98 & 1.01 | 1.622e-03 & 1.844e-01 | 1.99 & 1.00
64 | 1.945e-04 & 4.620e-02 | 1.99 & 1.00 | 4.069e-04 & 9.214e-02 | 1.99 & 1.00
128 | 4.911e-05 & 2.306e-02 | 1.99 & 1.00 | 1.024e-04 & 4.605e-02 | 1.99 & 1.00
Table 8: Convergence rate of the iterative decoupled algorithm. $K_{1}=K_{2}=10^{-6}$, $\Delta t=2\times 10^{-3}s$, and $iters=10$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 1.025e-03 & 1.676e-02 |  | 3.688e-02 & 1.092e+00 |
16 | 2.301e-04 & 3.326e-03 | 2.16 & 2.33 | 9.248e-03 & 5.561e-01 | 2.00 & 0.973
32 | 5.618e-05 & 6.874e-04 | 2.03 & 2.27 | 2.337e-03 & 2.772e-01 | 1.98 & 1.004
64 | 1.588e-05 & 1.543e-04 | 1.82 & 2.16 | 5.926e-04 & 1.384e-01 | 1.98 & 1.002
128 | 7.569e-06 & 4.856e-05 | 1.07 & 1.67 | 1.557e-04 & 6.915e-02 | 1.93 & 1.001
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.217e-02 & 3.663e-01 |  | 2.592e-02 & 7.288e-01 |
16 | 3.054e-03 & 1.863e-01 | 1.99 & 0.98 | 6.474e-03 & 3.690e-01 | 2.00 & 0.98
32 | 7.753e-04 & 9.269e-02 | 1.98 & 1.01 | 1.630e-03 & 1.844e-01 | 1.99 & 1.00
64 | 1.971e-04 & 4.621e-02 | 1.98 & 1.00 | 4.130e-04 & 9.215e-02 | 1.98 & 1.00
128 | 5.370e-05 & 2.308e-02 | 1.88 & 1.00 | 1.093e-04 & 4.606e-02 | 1.92 & 1.00
Table 9: Convergence rate of the iterative decoupled algorithm. $K_{1}=K_{2}=10^{-6}$, $\Delta t=4\times 10^{-3}s$, and $iters=20$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 1.021e-03 & 1.664e-02 |  | 3.675e-02 & 1.091e+00 |
16 | 2.289e-04 & 3.291e-03 | 2.16 & 2.34 | 9.226e-03 & 5.560e-01 | 1.99 & 0.973
32 | 5.500e-05 & 6.757e-04 | 2.06 & 2.28 | 2.341e-03 & 2.772e-01 | 1.98 & 1.004
64 | 1.389e-05 & 1.474e-04 | 1.99 & 2.20 | 6.010e-04 & 1.384e-01 | 1.96 & 1.002
128 | 3.871e-06 & 3.579e-05 | 1.84 & 2.04 | 1.636e-04 & 6.914e-02 | 1.88 & 1.001
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.209e-02 & 3.663e-01 |  | 2.588e-02 & 7.288e-01 |
16 | 3.038e-03 & 1.863e-01 | 1.99 & 0.98 | 6.470e-03 & 3.690e-01 | 2.00 & 0.98
32 | 7.750e-04 & 9.269e-02 | 1.97 & 1.01 | 1.634e-03 & 1.844e-01 | 1.99 & 1.00
64 | 1.994e-04 & 4.621e-02 | 1.96 & 1.00 | 4.184e-04 & 9.215e-02 | 1.97 & 1.00
128 | 5.432e-05 & 2.307e-02 | 1.88 & 1.00 | 1.139e-04 & 4.606e-02 | 1.88 & 1.00

###  5.3 Tests for the parameter $c_{i}$

In this subsection, we will evaluate the influence of the specific storage
coefficient $c_{i}$ on the iterative decoupled algorithm. We now let
$c_{1}=c_{2}=0$. For other physical parameters, we set

$\displaystyle
E=1,\nu=0.3,\alpha_{1}=\alpha_{2}=1,K_{1}=K_{2}=1,\beta_{12}=\beta_{21}=1.$

Comparing the results in Table 2 with those in Table 11, it is evident that in
cases where $c_{1}=c_{2}=0$, the error rates for all variables exhibit a
slight degradation when utilizing the iterative decoupled algorithm. However,
a closer examination of Table 12 reveals that as the number of iterations
increases to 20, the error rates for all variables in the iterative decoupled
algorithm tend to decrease, and the energy-norm errors attain their optimal
convergence rates. These experimental findings corroborate the assertion made
in Remark 1.

Table 10: Convergence rate of the coupled algorithm. $c_{1}=c_{2}=0$ and $\Delta t=2\times 10^{-4}s$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 8.027e-04 & 1.657e-02 |  | 3.108e-02 & 1.084e+00 |
16 | 1.860e-04 & 3.729e-03 | 2.11 & 2.15 | 7.778e-03 & 5.506e-01 | 2.00 & 0.977
32 | 4.606e-05 & 8.620e-04 | 2.01 & 2.11 | 1.941e-03 & 2.760e-01 | 2.00 & 0.996
64 | 1.152e-05 & 2.047e-04 | 2.00 & 2.07 | 4.845e-04 & 1.381e-01 | 2.00 & 0.999
128 | 2.883e-06 & 4.967e-05 | 2.00 & 2.04 | 1.210e-04 & 6.908e-02 | 2.00 & 1.000
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 9.448e-03 & 3.622e-01 |  | 2.296e-02 & 7.210e-01 |
16 | 2.417e-03 & 1.829e-01 | 1.97 & 0.99 | 5.890e-03 & 3.652e-01 | 1.96 & 0.98
32 | 6.097e-04 & 9.184e-02 | 1.99 & 0.99 | 1.488e-03 & 1.836e-01 | 1.98 & 0.99
64 | 1.529e-04 & 4.600e-02 | 2.00 & 1.00 | 3.734e-04 & 9.198e-02 | 1.99 & 1.00
128 | 3.827e-05 & 2.301e-02 | 2.00 & 1.00 | 9.345e-05 & 4.602e-02 | 2.00 & 1.00
Table 11: Convergence rate of the iterative decoupled algorithm. $c_{1}=c_{2}=0$, $\Delta t=2\times 10^{-3}s$, and $iters=10$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 8.829e-04 & 1.673e-02 |  | 3.264e-02 & 1.084e+00 |
16 | 2.080e-04 & 3.760e-03 | 2.09 & 2.15 | 8.127e-03 & 5.507e-01 | 2.01 & 0.977
32 | 5.719e-05 & 8.770e-04 | 1.86 & 2.10 | 2.011e-03 & 2.760e-01 | 2.01 & 0.996
64 | 2.944e-05 & 2.442e-04 | 0.96 & 1.84 | 4.972e-04 & 1.381e-01 | 2.02 & 0.999
128 | 2.702e-05 & 1.421e-04 | 0.12 & 0.78 | 1.529e-04 & 6.909e-02 | 1.70 & 1.000
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 1.011e-02 & 3.618e-01 |  | 2.383e-02 & 7.208e-01 |
16 | 2.554e-03 & 1.828e-01 | 1.99 & 0.98 | 6.081e-03 & 3.652e-01 | 1.97 & 0.98
32 | 6.388e-04 & 9.184e-02 | 2.00 & 0.99 | 1.527e-03 & 1.836e-01 | 1.99 & 0.99
64 | 1.722e-04 & 4.600e-02 | 1.89 & 1.00 | 3.834e-04 & 9.199e-02 | 1.99 & 1.00
128 | 8.820e-05 & 2.302e-02 | 0.97 & 1.00 | 1.191e-04 & 4.603e-02 | 1.69 & 1.00
Table 12: Convergence rate of the iterative decoupled algorithm. $c_{1}=c_{2}=0$, $\Delta t=4\times 10^{-3}s$, and $iters=20$. $1/h$ |  $L^{2}\&H^{1}$ errors of $\boldsymbol{u}$ | Orders |  $L^{2}\&H^{1}$ errors of $\xi$ | Orders
---|---|---|---|---
8 | 7.340e-04 & 1.620e-02 |  | 3.063e-02 & 1.084e+00 |
16 | 1.643e-04 & 3.555e-03 | 2.16 & 2.19 | 7.616e-03 & 5.508e-01 | 2.01 & 0.977
32 | 4.024e-05 & 8.060e-04 | 2.03 & 2.14 | 1.890e-03 & 2.760e-01 | 2.01 & 0.997
64 | 9.963e-06 & 1.890e-04 | 2.01 & 2.09 | 4.675e-04 & 1.381e-01 | 2.02 & 0.999
128 | 2.411e-06 & 4.533e-05 | 2.05 & 2.06 | 1.130e-04 & 6.908e-02 | 2.05 & 1.000
$1/h$ |  $L^{2}\&H^{1}$ errors of $p_{1}$ | Orders |  $L^{2}\&H^{1}$ errors of $p_{2}$ | Orders
8 | 9.229e-03 & 3.627e-01 |  | 2.263e-02 & 7.213e-01 |
16 | 2.323e-03 & 1.830e-01 | 1.99 & 0.99 | 5.769e-03 & 3.653e-01 | 1.97 & 0.98
32 | 5.801e-04 & 9.186e-02 | 2.00 & 0.99 | 1.451e-03 & 1.836e-01 | 1.99 & 0.99
64 | 1.437e-04 & 4.600e-02 | 2.01 & 1.00 | 3.619e-04 & 9.199e-02 | 2.00 & 1.00
128 | 3.457e-05 & 2.301e-02 | 2.06 & 1.00 | 8.866e-05 & 4.602e-02 | 2.03 & 1.00

##  6 Conclusion

In this paper, we extend the iterative decoupled algorithm of Biot’s
consolidation model proposed in [26] to the multiple-network poroelasticity
problem. Our analysis affirms that the sequence
${(\boldsymbol{u}_{h}^{n,k},\xi_{h}^{n,k},\vec{p}_{h}^{n,k})}$ produced by the
iterative decoupled algorithm converges to the solution
${(\boldsymbol{u}_{h}^{n},\xi_{h}^{n},\vec{p}_{h}^{n})}$ acquired through the
coupled algorithm. Furthermore, the iterative decoupled algorithm demonstrates
unconditional stability and convergence characteristics. We validate the
soundness of our analysis through numerical experiments.

##  Appendix A Error analysis of the coupled algorithm

This appendix provides the optimal order error estimate for the coupled
algorithm. Specifically, we demonstrate that the temporal error is of order
$O(\Delta t)$; The energy-norm errors for $\boldsymbol{u}$ and $\vec{p}$ are
of order $O(h^{2})$ and order $O(h)$, respectively; The $L^{2}$ norm errors
for $\xi$ and $\vec{p}$ are of order $O(h^{2})$. On the continuous level, we
refer readers to [2] for the energy estimate of the variational problem (4).
In the following lemma, we provide the energy estimate for the coupled
algorithm.

######  Lemma A.1.

Let $\\{(\boldsymbol{u}_{h}^{n},\xi_{h}^{n},\vec{p}_{h}^{n})\\}$ be the
solutions of problem (6a)-(6c), then the following identity holds

$\displaystyle J_{h}^{l}+S_{h}^{l}=J_{h}^{0}\quad
for\,l\in\\{1,2,\cdots,L\\}.$ (A.1)

where

$\displaystyle
J_{h}^{l}:=\mu|\\!|\epsilon(\boldsymbol{u}_{h}^{l})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}\vec{p}_{h}^{l}-\xi_{h}^{l}|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|\hat{S}\vec{p}_{h}^{l}|\\!|_{L^{2}(\Omega)}^{2}-(\boldsymbol{f},\boldsymbol{u}_{h}^{l})-\langle\boldsymbol{h},\boldsymbol{u}_{h}^{l}\rangle_{\Gamma_{\boldsymbol{u},N}}.$
$\displaystyle S_{h}^{l}:=$ $\displaystyle\Delta t\sum_{n=1}^{l}\bigg{[}\Delta
t\left(\mu|\\!|d_{t}\epsilon(\boldsymbol{u}_{h}^{n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|d_{t}(\vec{\alpha}^{\mathsf{T}}\vec{p}_{h}^{n}-\xi_{h}^{n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|d_{t}(\hat{S}\vec{p}_{h}^{n})|\\!|_{L^{2}(\Omega)}^{2}\right)$
$\displaystyle+|\\!|\hat{K}\nabla\vec{p}_{h}^{n}|\\!|_{L^{2}(\Omega)}^{2}+\sum_{1\leq
i<j\leq
N}\beta_{i,j}|\\!|p_{i,h}^{n}-p_{j,h}^{n}|\\!|_{L^{2}(\Omega)}^{2}-(\vec{g},\vec{p}_{h}^{n})-\left\langle\vec{l},\vec{p}_{h}^{n}\right\rangle_{\Gamma_{\vec{p},N}}\bigg{]}.$

Here, we denote $d_{t}\eta^{n}:=\frac{\eta^{n}-\eta^{n-1}}{\Delta t}$, where
$\eta$ can be a vector or a scalar,
$\hat{S}=\text{diag}(\sqrt{c_{1}},\sqrt{c_{2}},\cdots,\sqrt{c_{N}}),\hat{K}=\text{diag}(\sqrt{K_{1}},\sqrt{K_{2}},\cdots,\sqrt{K_{N}})$.
Moreover, there holds

$\displaystyle|\\!|\xi_{h}^{l}|\\!|_{L^{2}(\Omega)}\leq
C(|\\!|\epsilon(\boldsymbol{u}_{h}^{l})|\\!|_{L^{2}(\Omega)}+|\\!|\boldsymbol{f}|\\!|_{L^{2}(\Omega)}+|\\!|\boldsymbol{h}|\\!|_{L^{2}(\Gamma_{\boldsymbol{u},N})}),$
(A.2)

where $C$ is a positive constant.

Proof. Setting $\boldsymbol{v}_{h}=d_{t}\boldsymbol{u}_{h}^{n}$ in (6a),
$\eta_{h}=\xi_{h}^{n}$ in (6b), $\vec{q}_{h}=\vec{p}_{h}^{n}$ in (6c), we have

$\displaystyle
2\mu(\epsilon(\boldsymbol{u}_{h}^{n}),d_{t}\epsilon(\boldsymbol{u}_{h}^{n}))-(\xi_{h}^{n},\operatorname{div}d_{t}\boldsymbol{u}_{h}^{n})$
$\displaystyle=d_{t}(\boldsymbol{f},\boldsymbol{u}_{h}^{n})+d_{t}\langle\boldsymbol{h},\boldsymbol{u}_{h}^{n}\rangle_{\Gamma_{\boldsymbol{u},N}},$
(A.3)
$\displaystyle(\operatorname{div}d_{t}\boldsymbol{u}_{h}^{n},\xi_{h}^{n})+\frac{1}{\lambda}(d_{t}\xi_{h}^{n},\xi_{h}^{n})$
$\displaystyle=\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}d_{t}\vec{p}_{h}^{n},\xi_{h}^{n}),$
(A.4)
$\displaystyle\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)d_{t}\vec{p}_{h}^{n},\vec{p}_{h}^{n}\right)-\frac{1}{\lambda}\left(\vec{\alpha}d_{t}\xi_{h}^{n},\vec{p}_{h}^{n}\right)$
$\displaystyle+\left(K\nabla\vec{p}_{h}^{n},\nabla\vec{p}_{h}^{n}\right)+\left(B\vec{p}_{h}^{n},\vec{p}_{h}^{n}\right)$
(A.5)
$\displaystyle=(\vec{g},\vec{p}_{h}^{n})+\left\langle\vec{l},\vec{p}_{h}^{n}\right\rangle_{\Gamma_{\vec{p},N}}.$

Summing up (A.3)-(A.5), we get

$\displaystyle
2\mu(\epsilon(\boldsymbol{u}_{h}^{n}),d_{t}\epsilon(\boldsymbol{u}_{h}^{n}))-d_{t}(\boldsymbol{f},\boldsymbol{u}_{h}^{n})-d_{t}\langle\boldsymbol{h},\boldsymbol{u}_{h}^{n}\rangle_{\Gamma_{\boldsymbol{u},N}}+\frac{1}{\lambda}(d_{t}\xi_{h}^{n},\xi_{h}^{n})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}d_{t}\vec{p}_{h}^{n},\xi_{h}^{n})$
(A.6)
$\displaystyle+\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)d_{t}\vec{p}_{h}^{n},\vec{p}_{h}^{n}\right)-\frac{1}{\lambda}\left(\vec{\alpha}d_{t}\xi_{h}^{n},\vec{p}_{h}^{n}\right)+\left(K\nabla\vec{p}_{h}^{n},\nabla\vec{p}_{h}^{n}\right)+\left(B\vec{p}_{h}^{n},\vec{p}_{h}^{n}\right)-(\vec{g},\vec{p}_{h}^{n})$
$\displaystyle-\left\langle\vec{l},\vec{p}_{h}^{n}\right\rangle_{\Gamma_{\vec{p},N}}=0.$

Note that, the following identities hold.

$\displaystyle
2(\eta_{h}^{n},d_{t}\eta_{h}^{n})=d_{t}|\\!|\eta_{h}^{n}|\\!|_{L^{2}(\Omega)}^{2}$
$\displaystyle+\Delta t|\\!|d_{t}\eta_{h}^{n}|\\!|_{L^{2}(\Omega)}^{2},$
$\displaystyle\frac{1}{\lambda}(d_{t}\xi_{h}^{n},\xi_{h}^{n})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}d_{t}\vec{p}_{h}^{n},\xi_{h}^{n})$
$\displaystyle+\frac{1}{\lambda}\left(\vec{\alpha}\vec{\alpha}^{\mathsf{T}}d_{t}\vec{p}_{h}^{n},\vec{p}_{h}^{n}\right)-\frac{1}{\lambda}\left(\vec{\alpha}d_{t}\xi_{h}^{n},\vec{p}_{h}^{n}\right)$
$\displaystyle=\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}\vec{p}_{h}^{n}-\xi_{h}^{n},d_{t}(\vec{\alpha}^{\mathsf{T}}\vec{p}_{h}^{n}-\xi_{h}^{n})),$
$\displaystyle\left(B\vec{p}_{h}^{n},\vec{p}_{h}^{n}\right)=\sum_{1\leq
i<j\leq N}\beta_{i,j}|\\!|$ $\displaystyle
p_{i,h}^{n}-p_{j,h}^{n}|\\!|_{L^{2}(\Omega)}^{2}.$

According to the above identities, (A.6) can be rewritten as

$\displaystyle
d_{t}\bigg{(}\mu|\\!|\epsilon(\boldsymbol{u}_{h}^{n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}\vec{p}_{h}^{n}-\xi_{h}^{n}|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|\hat{S}\vec{p}_{h}^{n}|\\!|_{L^{2}(\Omega)}^{2}-(\boldsymbol{f},\boldsymbol{u}_{h}^{n})$
(A.7)
$\displaystyle-\langle\boldsymbol{h},\boldsymbol{u}_{h}^{n}\rangle_{\Gamma_{\boldsymbol{u},N}}\bigg{)}+\Delta
t\left(\mu|\\!|d_{t}\epsilon(\boldsymbol{u}_{h}^{n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|d_{t}(\vec{\alpha}^{\mathsf{T}}\vec{p}_{h}^{n}-\xi_{h}^{n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|d_{t}(\hat{S}\vec{p}_{h}^{n})|\\!|_{L^{2}(\Omega)}^{2}\right)$
$\displaystyle+|\\!|\hat{K}\nabla\vec{p}_{h}^{n}|\\!|_{L^{2}(\Omega)}^{2}+\sum_{1\leq
i<j\leq
N}\beta_{i,j}|\\!|p_{i,h}^{n}-p_{j,h}^{n}|\\!|_{L^{2}(\Omega)}^{2}-(\vec{g},\vec{p}_{h}^{n})-\left\langle\vec{l},\vec{p}_{h}^{n}\right\rangle_{\Gamma_{\vec{p},N}}=0.$

Applying the summation operator $\Delta t\sum_{n=1}^{l}$ to both sides of the
(A.7), we get (A.1). In addition, using (5) and (6a), we has (A.2).
$\hfill\qed$

As a preliminary step for a priori estimate of the discrete formulation (6),
we introduce a set of auxiliary projection operators. In particular, we define
projection operators

$\Pi_{h}^{\boldsymbol{V}}:\boldsymbol{V}\rightarrow\boldsymbol{V}_{h},~{}~{}~{}~{}\Pi_{h}^{W}:W\rightarrow
W_{h},$

as follows. For any $(\boldsymbol{u},\xi)\in\boldsymbol{V}\times W$, we define
its projection
$(\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u},\Pi_{h}^{W}\xi)\in\boldsymbol{V}_{h}\times
W_{h}$ as the unique discrete solution to the Stokes problem:

$\displaystyle
2\mu(\epsilon(\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u}),\epsilon(\boldsymbol{v}_{h}))-(\Pi_{h}^{W}\xi,\operatorname{div}\boldsymbol{v}_{h})$
$\displaystyle=2\mu(\epsilon(\boldsymbol{u}),\epsilon(\boldsymbol{v}_{h}))-(\xi,\operatorname{div}\boldsymbol{v}_{h})\quad\forall\boldsymbol{v}_{h}\in\boldsymbol{V}_{h},$
(A.8a)
$\displaystyle-(\operatorname{div}\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u},\eta_{h})$
$\displaystyle=-(\operatorname{div}\boldsymbol{u},\eta_{h})\quad\forall\eta_{h}\in
W_{h}.$ (A.8b)

We also define the projection operator

$\Pi_{h}^{M}:M\rightarrow M_{h},$

as follows. For any $\vec{p}\in M$, we define its projection
$\Pi_{h}^{M}\vec{p}\in M_{h}$ as the unique discrete solution to the weighted
elliptic problem

$\displaystyle(K\nabla\Pi_{h}^{M}\vec{p},\nabla\vec{q}_{h})=(K\nabla\vec{p},\nabla\vec{q}_{h})\quad\forall\vec{q}_{h}\in
M_{h}.$ (A.9)

Here, we list the properties of the operators
$(\Pi_{h}^{\boldsymbol{V}},\Pi_{h}^{W},\Pi_{h}^{M})$. For all
$(\boldsymbol{u},\xi,\vec{p})\in[H^{3}(\Omega)]^{2}\times
H^{2}(\Omega)\times[H^{2}(\Omega)]^{N}$, there holds [2]

$\displaystyle|\\!|\boldsymbol{u}-\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u}|\\!|_{H^{1}(\Omega)}+|\\!|\xi-\Pi_{h}^{W}\xi|\\!|_{L^{2}(\Omega)}$
$\displaystyle\leq
Ch^{2}\left(|\\!|\boldsymbol{u}|\\!|_{H^{3}(\Omega)}+|\\!|\xi|\\!|_{H^{2}(\Omega)}\right),$
(A.10) $\displaystyle|\\!|\vec{p}-\Pi_{h}^{M}\vec{p}|\\!|_{H^{1}(\Omega)}$
$\displaystyle\leq Ch|\\!|\vec{p}|\\!|_{H^{2}(\Omega)},$ (A.11)
$\displaystyle|\\!|\vec{p}-\Pi_{h}^{M}\vec{p}|\\!|_{L^{2}(\Omega)}$
$\displaystyle\leq Ch^{2}|\\!|\vec{p}|\\!|_{H^{2}(\Omega)}.$ (A.12)

For convenience, we introduce the standard decomposition of the errors into
projection and discretization errors in the subsequent analysis,

$\displaystyle
e_{\boldsymbol{u}}^{n}=\boldsymbol{u}^{n}-\boldsymbol{u}_{h}^{n}=(\boldsymbol{u}^{n}-\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u}^{n})+(\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u}^{n}-\boldsymbol{u}_{h}^{n}):=e_{\boldsymbol{u}}^{I,n}+e_{\boldsymbol{u}}^{h,n},$
$\displaystyle
e_{\xi}^{n}=\xi^{n}-\xi_{h}^{n}=(\xi^{n}-\Pi_{h}^{W}\xi^{n})+(\Pi_{h}^{W}\xi^{n}-\xi_{h}^{n}):=e_{\xi}^{I,n}+e_{\xi}^{h,n},$
$\displaystyle
e_{\vec{p}}^{n}=\vec{p}^{n}-\vec{p}_{h}^{n}=(\vec{p}^{n}-\Pi_{h}^{M}\vec{p}^{n})+(\Pi_{h}^{M}\vec{p}^{n}-\vec{p}_{h}^{n}):=e_{\vec{p}}^{I,n}+e_{\vec{p}}^{h,n}.$

######  Lemma A.2.

Let $\\{(\boldsymbol{u}_{h}^{n},\xi_{h}^{n},\vec{p}_{h}^{n})\\}$ be the
solutions of problem (6a)-(6c), then the following identity holds

$\displaystyle E_{h}^{l}+\Delta t\sum_{n=1}^{l}\left(|\\!|\hat{K}\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}+\left(Be_{\vec{p}}^{h,n},e_{\vec{p}}^{h,n}\right)\right)$
(A.13) $\displaystyle+(\Delta
t)^{2}\sum_{n=1}^{l}\left(\mu|\\!|d_{t}\epsilon(e_{\boldsymbol{u}}^{h,n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|d_{t}(\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{h,n}-e_{\xi}^{h,n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|d_{t}\hat{S}e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}\right)$
$\displaystyle=E_{h}^{0}+\Delta
t\sum_{n=1}^{l}\Bigg{[}\left(\operatorname{div}(d_{t}\boldsymbol{u}^{n}-\dot{\boldsymbol{u}}^{n}),e_{\xi}^{h,n}\right)+\frac{1}{\lambda}\left(d_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n},e_{\xi}^{h,n}\right)$
$\displaystyle-\frac{1}{\lambda}\left(\vec{\alpha}^{\mathsf{T}}(d_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}),e_{\xi}^{h,n}\right)+\left(B\left(\Pi_{h}^{M}\vec{p}^{n}-\vec{p}^{n}\right),e_{\vec{p}}^{h,n}\right)$
$\displaystyle-\frac{1}{\lambda}\left(\vec{\alpha}\left(d_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n}\right),e_{\vec{p}}^{h,n}\right)+\bigg{(}\bigg{(}S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\bigg{)}\left(d_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}\right),e_{\vec{p}}^{h,n}\bigg{)}\Bigg{]},$

where

$E_{h}^{l}:=\mu|\\!|\epsilon(e_{\boldsymbol{u}}^{h,l})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{h,l}-e_{\xi}^{h,l}|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|\hat{S}e_{\vec{p}}^{h,l}|\\!|_{L^{2}(\Omega)}^{2}.$

Proof. Letting $t=t_{n},\boldsymbol{v}=\boldsymbol{v}_{h}$ in (4a),
subtracting (6a) from equation (4a) and combining with (A.8), there holds

$\displaystyle
2\mu(\epsilon(e_{\boldsymbol{u}}^{h,n}),\epsilon(\boldsymbol{v}_{h}))-(e_{\xi}^{h,n},\operatorname{div}\boldsymbol{v}_{h})=0,$
(A.14)

Taking the derivative of (4b) with respect to time $t$, and letting
$t=t_{n},\eta=\eta_{h}$, we have

$\displaystyle(\text{div}\dot{\boldsymbol{u}}^{n},\eta_{h})+\frac{1}{\lambda}(\dot{\xi}^{n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}\dot{\vec{p}}^{n},\eta_{h})=0.$
(A.15)

Using (6b), it holds that

$\displaystyle(\operatorname{div}d_{t}\boldsymbol{u}^{n},\eta_{h})+\frac{1}{\lambda}(d_{t}\xi^{n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}d_{t}\vec{p}^{n},\eta_{h})=0.$
(A.16)

By subtracting (A.15) from (A.16), we obtain

$\displaystyle(\operatorname{div}(d_{t}\boldsymbol{u}^{n}-\dot{\boldsymbol{u}}^{n}),\eta_{h})+\frac{1}{\lambda}(d_{t}\xi^{n}-\dot{\xi}^{n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}(d_{t}\vec{p}^{n}-\dot{\vec{p}}^{n}),\eta_{h})=0,$
(A.17)

The combination of (4b), (6b) and (A.8) implies that

$\displaystyle(\operatorname{div}e_{\boldsymbol{u}}^{h,n},\eta_{h})+\frac{1}{\lambda}(e_{\xi}^{I,n}+e_{\xi}^{h,n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}(e_{\vec{p}}^{I,n}+e_{\vec{p}}^{h,n}),\eta_{h})=0,$
(A.18)
$\displaystyle(\operatorname{div}e_{\boldsymbol{u}}^{h,n-1},\eta_{h})+\frac{1}{\lambda}(e_{\xi}^{I,n-1}+e_{\xi}^{h,n-1},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}(e_{\vec{p}}^{I,n-1}+e_{\vec{p}}^{h,n-1}),\eta_{h})=0.$
(A.19)

Using (A.19) and (A.18), we get

$\displaystyle(\operatorname{div}(d_{t}e_{\boldsymbol{u}}^{h,n}),\eta_{h})+\frac{1}{\lambda}(d_{t}e_{\xi}^{I,n}+d_{t}e_{\xi}^{h,n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}(d_{t}e_{\vec{p}}^{I,n}+d_{t}e_{\vec{p}}^{h,n}),\eta_{h})=0.$
(A.20)

Combining (A.17) with (A.20), we obtain

$\displaystyle(\operatorname{div}(d_{t}e_{\boldsymbol{u}}^{h,n}),\eta_{h})+\frac{1}{\lambda}(d_{t}e_{\xi}^{h,n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}d_{t}e_{\vec{p}}^{h,n},\eta_{h})=(\operatorname{div}(d_{t}\boldsymbol{u}^{n}-\dot{\boldsymbol{u}}^{n}),\eta_{h})$
(A.21)
$\displaystyle+\frac{1}{\lambda}(d_{t}\xi^{n}-\dot{\xi}^{n},\eta_{h})-\frac{1}{\lambda}(d_{t}e_{\xi}^{I,n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}(d_{t}\vec{p}^{n}-\dot{\vec{p}}^{n}),\eta_{h})+\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}d_{t}e_{\vec{p}}^{I,n},\eta_{h}).$

The above equation (A.21) can also be written as

$\displaystyle(\operatorname{div}(d_{t}e_{\boldsymbol{u}}^{h,n}),\eta_{h})+\frac{1}{\lambda}(d_{t}e_{\xi}^{h,n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}d_{t}e_{\vec{p}}^{h,n},\eta_{h})=(\operatorname{div}(d_{t}\boldsymbol{u}^{n}-\dot{\boldsymbol{u}}^{n}),\eta_{h})$
(A.22)
$\displaystyle+\frac{1}{\lambda}(d_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n},\eta_{h})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}(d_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}),\eta_{h}).$

In (4c), we let $t=t_{n}$ and $\vec{q}=\vec{q}_{h}$, then subtract equation
(6c) from it. Combining with (A.8) and (A.9), we have

$\displaystyle\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)d_{t}e_{\vec{p}}^{h,n},\vec{q}_{h}\right)-\frac{1}{\lambda}\left(\vec{\alpha}d_{t}e_{\xi}^{h,n},\vec{q}_{h}\right)+\left(K\nabla
e_{\vec{p}}^{h,n},\nabla\vec{q}_{h}\right)$ (A.23)
$\displaystyle+\left(Be_{\vec{p}}^{h,n},\vec{q}_{h}\right)=\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)\left(d_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}\right),\vec{q}_{h}\right)-\frac{1}{\lambda}\left(\vec{\alpha}\left(d_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n}\right),\vec{q}_{h}\right)$
$\displaystyle+\left(B\left(\Pi_{h}^{M}\vec{p}^{n}-\vec{p}^{n}\right),\vec{q}_{h}\right).$

In (A.14), (A.22) and (A.23), letting
$\boldsymbol{v}_{h}=d_{t}e_{\boldsymbol{u}}^{h,n}$, $\eta_{h}=e_{\xi}^{h,n}$,
$\vec{q}_{h}=e_{\vec{p}}^{h,n}$, then adding them up, one can get

$\displaystyle
2\mu(\epsilon(e_{\boldsymbol{u}}^{h,n}),\epsilon(d_{t}e_{\boldsymbol{u}}^{h,n}))+\frac{1}{\lambda}(d_{t}e_{\xi}^{h,n},e_{\xi}^{h,n})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}d_{t}e_{\vec{p}}^{h,n},e_{\xi}^{h,n})+\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)d_{t}e_{\vec{p}}^{h,n},e_{\vec{p}}^{h,n}\right)$
$\displaystyle-\frac{1}{\lambda}\left(\vec{\alpha}d_{t}e_{\xi}^{h,n},e_{\vec{p}}^{h,n}\right)+\left(K\nabla
e_{\vec{p}}^{h,n},\nabla
e_{\vec{p}}^{h,n}\right)+\left(Be_{\vec{p}}^{h,n},e_{\vec{p}}^{h,n}\right)=(\operatorname{div}(d_{t}\boldsymbol{u}^{n}-\dot{\boldsymbol{u}}^{n}),e_{\xi}^{h,n})$
$\displaystyle+\frac{1}{\lambda}(d_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n},e_{\xi}^{h,n})-\frac{1}{\lambda}(\vec{\alpha}^{\mathsf{T}}(d_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}),e_{\xi}^{h,n})-\frac{1}{\lambda}\left(\vec{\alpha}\left(d_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n}\right),e_{\vec{p}}^{h,n}\right)$
$\displaystyle+\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)\left(d_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}\right),e_{\vec{p}}^{h,n}\right)+\left(B\left(\Pi_{h}^{M}\vec{p}^{n}-\vec{p}^{n}\right),e_{\vec{p}}^{h,n}\right).$

The above formula can also be written as

$\displaystyle
d_{t}\left(\mu|\\!|\epsilon(e_{\boldsymbol{u}}^{h,n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{h,n}-e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|\hat{S}e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}\right)$
(A.24) $\displaystyle+|\\!|\hat{K}\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}+\left(Be_{\vec{p}}^{h,n},e_{\vec{p}}^{h,n}\right)$
$\displaystyle+\Delta
t\left(\mu|\\!|d_{t}\epsilon(e_{\boldsymbol{u}}^{h,n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|d_{t}(\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{h,n}-e_{\xi}^{h,n})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|d_{t}\hat{S}e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}\right)$
$\displaystyle=\left(\operatorname{div}(d_{t}\boldsymbol{u}^{n}-\dot{\boldsymbol{u}}^{n}),e_{\xi}^{h,n}\right)+\frac{1}{\lambda}\left(d_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n},e_{\xi}^{h,n}\right)-\frac{1}{\lambda}\left(\vec{\alpha}^{\mathsf{T}}(d_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}),e_{\xi}^{h,n}\right)$
$\displaystyle+\left(\left(S+\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}\right)\left(d_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}\right),e_{\vec{p}}^{h,n}\right)-\frac{1}{\lambda}\left(\vec{\alpha}\left(d_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n}\right),e_{\vec{p}}^{h,n}\right)$
$\displaystyle+\left(B\left(\Pi_{h}^{M}\vec{p}^{n}-\vec{p}^{n}\right),e_{\vec{p}}^{h,n}\right).$

Applying the summation operator $\Delta t\sum_{n=1}^{l}$ to both sides, we
obtain (A.13). $\hfill\qed$

In the rest part of the analysis, we will use the following discrete Gronwall
inequality [29].

######  Lemma A.3.

Let $\tau,\,B$ and $a_{j},\,b_{j},\,c_{j},\,\gamma_{j},\,\forall j\geq 1$ be
non-negative numbers such that

$\displaystyle a_{n}+\sum_{j=1}^{n}b_{j}\leq
B+\tau\sum_{j=1}^{n}c_{j}+\tau\sum_{j=1}^{n}\gamma_{j}a_{j}.$ (A.25)

If $\tau\gamma_{j}\leq 1,j=1,\cdots,n$, there holds

$\displaystyle a_{n}+\sum_{j=1}^{n}b_{j}\leq
C\left(B+\tau\sum_{j=1}^{n}c_{j}\right),$ (A.26)

where $C=e^{\tau\sum_{j=1}^{n}\frac{\gamma_{j}}{1-\tau\gamma_{j}}}$.

######  Theorem A.4.

Let $\\{(\boldsymbol{u}_{h}^{n},\xi_{h}^{n},\vec{p}_{h}^{n})\\}$ be the
solution of problem (6a)-(6c), then the following error estimate holds

$\displaystyle E_{h}^{l}+\frac{1}{2}\Delta t\sum_{n=1}^{l}|\\!|\hat{K}\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}\leq C_{1}(\Delta
t)^{2}+C_{2}h^{4}+C_{3}\Delta th^{4},$ (A.27)

where

$\displaystyle
C_{1}=C_{1}\left(|\\!|\ddot{\boldsymbol{u}}|\\!|_{L^{2}(0,t_{l};H^{1}(\Omega))},|\\!|\ddot{\xi}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))},|\\!|\ddot{\vec{p}}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))}\right),$
$\displaystyle
C_{2}=C_{2}\left(|\\!|\dot{\boldsymbol{u}}|\\!|_{L^{2}(0,t_{l};H^{3}(\Omega))},|\\!|\dot{\xi}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))},|\\!|\dot{\vec{p}}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))}\right),$
$\displaystyle
C_{3}=C_{3}\left(|\\!|\vec{p}|\\!|_{L^{\infty}(0,t_{l};H^{2}(\Omega))}\right).$

Proof. According to the Cauchy-Schwarz inequality, the Young inequality with a
$\epsilon_{1}>0$, and Taylor’s formula, we obtain

$\displaystyle\Delta
t\sum_{n=1}^{l}\left(\operatorname{div}(\text{d}_{t}\boldsymbol{u}^{n}-\dot{\boldsymbol{u}}^{n}),e_{\xi}^{h,n}\right)$
(A.28)
$\displaystyle=\sum_{n=1}^{l}\left(\operatorname{div}(\boldsymbol{u}^{n}-\boldsymbol{u}^{n-1}-\Delta
t\dot{\boldsymbol{u}}^{n}),e_{\xi}^{h,n}\right)$
$\displaystyle\leq\sum_{n=1}^{l}|\\!|\boldsymbol{u}^{n}-\boldsymbol{u}^{n-1}-\Delta
t\dot{\boldsymbol{u}}^{n}|\\!|_{H^{1}(\Omega)}|\\!|e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}$
$\displaystyle\leq\left(\sum_{n=1}^{l}\frac{(\Delta
t)^{3}}{3}\int_{t_{n-1}}^{t_{n}}|\\!|\ddot{\boldsymbol{u}}|\\!|_{H^{1}(\Omega)}^{2}\text{d}s\right)^{\frac{1}{2}}\left(\sum_{n=1}^{l}|\\!|e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}\right)^{\frac{1}{2}}$
$\displaystyle\leq\left(\frac{(\Delta
t)^{2}}{3}\int_{t_{0}}^{t_{l}}|\\!|\ddot{\boldsymbol{u}}|\\!|_{H^{1}(\Omega)}^{2}\text{d}s\right)^{\frac{1}{2}}\left(\Delta
t\sum_{n=1}^{l}|\\!|e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}\right)^{\frac{1}{2}}$
$\displaystyle\leq\frac{(\Delta
t)^{2}}{6\epsilon_{1}}|\\!|\ddot{\boldsymbol{u}}|\\!|_{L^{2}(0,t_{l};H^{1}(\Omega))}^{2}+\frac{\epsilon_{1}\Delta
t}{2}\sum_{n=1}^{l}|\\!|e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}.$

Similarly, the following term can be bounded by

$\displaystyle\Delta
t\sum_{n=1}^{l}\frac{1}{\lambda}(\text{d}_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n},e_{\xi}^{h,n})$
(A.29)
$\displaystyle\leq\frac{1}{\lambda^{2}\epsilon_{2}}\left(Ch^{4}|\\!|\dot{\xi}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))}^{2}+Ch^{4}|\\!|\dot{\boldsymbol{u}}|\\!|_{L^{2}(0,t_{l};H^{3}(\Omega))}^{2}+\frac{(\Delta
t)^{2}}{3}|\\!|\ddot{\xi}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))}^{2}\right)$
$\displaystyle+\frac{\epsilon_{2}\Delta
t}{2}\sum_{n=1}^{l}|\\!|e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}.$

In consideration of (A.12), there holds

$\displaystyle\Delta
t\sum_{n=1}^{l}\frac{1}{\lambda}\left(\vec{\alpha}^{\mathsf{T}}(\text{d}_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}),e_{\xi}^{h,n}\right)$
(A.30)
$\displaystyle\leq\frac{1}{\lambda^{2}\epsilon_{3}}|\\!|\vec{\alpha}|\\!|^{2}\left(Ch^{4}|\\!|\dot{\vec{p}}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))}^{2}+\frac{(\Delta
t)^{2}}{3}|\\!|\ddot{\vec{p}}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))}^{2}\right)+\frac{\epsilon_{3}\Delta
t}{2}\sum_{n=1}^{l}|\\!|e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}.$

By use of estimate (A.12) and the Poincaré inequality, we get

$\displaystyle\Delta
t\sum_{n=1}^{l}\left(S\left(\text{d}_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}\right),e_{\vec{p}}^{h,n}\right)$
(A.31)
$\displaystyle\leq\frac{|\\!|S|\\!|^{2}}{\epsilon_{4}}\left(Ch^{4}|\\!|\dot{\vec{p}}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))}^{2}+\frac{(\Delta
t)^{2}}{3}|\\!|\ddot{\vec{p}}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))}^{2}\right)+\frac{\Delta
t\epsilon_{4}}{2}\sum_{n=1}^{l}|\\!|\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}.$

Likewise, based on the Cauchy-Schwarz inequality, we derive

$\displaystyle\Delta
t\sum_{n=1}^{l}\left(\frac{1}{\lambda}\vec{\alpha}\vec{\alpha}^{\mathsf{T}}(\text{d}_{t}\Pi_{h}^{M}\vec{p}^{n}-\dot{\vec{p}}^{n}),e_{\vec{p}}^{h,n}\right)$
(A.32)
$\displaystyle\leq\frac{|\\!|\vec{\alpha}\vec{\alpha}^{\mathsf{T}}|\\!|^{2}}{\lambda^{2}\epsilon_{5}}\left(Ch^{4}|\\!|\dot{\vec{p}}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))}^{2}+\frac{(\Delta
t)^{2}}{3}|\\!|\ddot{\vec{p}}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))}^{2}\right)+\frac{\Delta
t\epsilon_{5}}{2}\sum_{n=1}^{l}|\\!|\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}.$ $\displaystyle\Delta
t\sum_{n=1}^{l}\frac{1}{\lambda}(\vec{\alpha}(\text{d}_{t}\Pi_{h}^{W}\xi^{n}-\dot{\xi}^{n}),e_{\vec{p}}^{h,n})$
(A.33)
$\displaystyle\leq\frac{|\\!|\vec{\alpha}|\\!|^{2}}{\lambda^{2}\epsilon_{6}}\bigg{(}Ch^{4}|\\!|\dot{\xi}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))}^{2}+Ch^{4}|\\!|\dot{\boldsymbol{u}}|\\!|_{L^{2}(0,t_{l};H^{3}(\Omega))}^{2}+\frac{(\Delta
t)^{2}}{3}|\\!|\ddot{\xi}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))}^{2}\bigg{)}$
$\displaystyle+\frac{\Delta t\epsilon_{6}}{2}\sum_{n=1}^{l}|\\!|\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}.$ $\displaystyle\Delta
t\sum_{n=1}^{l}\left(Be_{\vec{p}}^{I,n},e_{\vec{p}}^{h,n}\right)\leq\frac{Cl|\\!|B|\\!|^{2}h^{4}}{2\epsilon_{7}}\Delta
t|\\!|\vec{p}|\\!|_{L^{\infty}(0,t_{l};H^{2}(\Omega))}^{2}+\frac{\Delta
t\epsilon_{7}}{2}\sum_{n=1}^{l}|\\!|\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}.$ (A.34)

By using the discrete inf-sup condition, we see that

$\displaystyle|\\!|e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}\leq\frac{1}{\beta_{0}^{*}}\sup_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\frac{(e_{\xi}^{h,n},\operatorname{div}\boldsymbol{v}_{h})}{|\\!|\boldsymbol{v}_{h}|\\!|_{H^{1}(\Omega)}}$
$\displaystyle=\frac{1}{\beta_{0}^{*}}\sup_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\frac{2\mu(\epsilon(e_{\boldsymbol{u}}^{h,n}),\epsilon(\boldsymbol{v}_{h}))}{|\\!|\boldsymbol{v}_{h}|\\!|_{H^{1}(\Omega)}}\leq\frac{2\mu}{\beta_{0}^{*}}|\\!|\epsilon(e_{\boldsymbol{u}}^{h,n})|\\!|_{L^{2}(\Omega)}.$

So we can choose $\epsilon_{i}$, $i\in\\{1,2,3\\}$, small enough such that

$\displaystyle\frac{\epsilon_{i}\Delta
t}{2}\sum_{n=1}^{l}|\\!|e_{\xi}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}\leq\frac{\Delta
t\mu}{3}\sum_{n=1}^{l}|\\!|\epsilon(e_{\boldsymbol{u}}^{h,n})|\\!|_{L^{2}(\Omega)}^{2}.$
(A.35)

Similarly, select $\epsilon_{i}$, $i\in\\{4,5,6,7\\}$, small enough to
guarantee

$\displaystyle\frac{\Delta t\epsilon_{i}}{2}\sum_{n=1}^{l}|\\!|\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}\leq\frac{\Delta
t}{8}\sum_{n=1}^{l}|\\!|\hat{K}\nabla
e_{\vec{p}}^{h,n}|\\!|_{L^{2}(\Omega)}^{2}.$ (A.36)

Combining (A.13), (A.28)-(A.36) and the discrete Gronwall inequality, we see
(A.27). $\hfill\qed$

######  Theorem A.5.

Let $\\{(\boldsymbol{u}_{h}^{n},\xi_{h}^{n},\vec{p}_{h}^{n})\\}$ be the
solution of problem (6a)-(6c), then the following error estimate holds

$\displaystyle\mu|\\!|\epsilon(e_{\boldsymbol{u}}^{l})|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{l}|\\!|_{L^{2}(\Omega)}^{2}+\frac{1}{2}|\\!|\hat{S}e_{\vec{p}}^{l}|\\!|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq C_{1}(\Delta t)^{2}+C_{2}h^{4}+C_{3}\Delta th^{4},$ (A.37)
$\displaystyle\Delta t\sum_{n=1}^{l}|\\!|\hat{K}\nabla
e_{\vec{p}}^{n}|\\!|_{L^{2}(\Omega)}^{2}$ $\displaystyle\leq C_{1}(\Delta
t)^{2}+C_{2}h^{4}+C_{3}\Delta th^{2},$ (A.38)
$\displaystyle|\\!|e_{\boldsymbol{u}}^{l}|\\!|_{H^{1}(\Omega)}^{2}$
$\displaystyle\leq C_{1}(\Delta t)^{2}+C_{2}h^{4}+C_{3}\Delta th^{4},$ (A.39)
$\displaystyle|\\!|e_{\xi}^{l}|\\!|_{L^{2}(\Omega)}^{2}$ $\displaystyle\leq
C_{1}(\Delta t)^{2}+C_{2}h^{4}+C_{3}\Delta th^{4},$ (A.40)

where

$\displaystyle
C_{1}=C_{1}\big{(}|\\!|\ddot{\boldsymbol{u}}|\\!|_{L^{2}(0,t_{l};H^{1}(\Omega))},|\\!|\ddot{\xi}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))},|\\!|\ddot{\vec{p}}|\\!|_{L^{2}(0,t_{l};L^{2}(\Omega))}\big{)},$
$\displaystyle
C_{2}=C_{2}\big{(}|\\!|\boldsymbol{u}|\\!|_{L^{\infty}(0,t_{l};H^{3}(\Omega))},|\\!|\xi|\\!|_{L^{\infty}(0,t_{l};H^{2}(\Omega))},|\\!|\vec{p}|\\!|_{L^{\infty}(0,t_{l};H^{2}(\Omega))},|\\!|\dot{\boldsymbol{u}}|\\!|_{L^{2}(0,t_{l};H^{3}(\Omega))},$
$\displaystyle|\\!|\dot{\xi}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))},|\\!|\dot{\vec{p}}|\\!|_{L^{2}(0,t_{l};H^{2}(\Omega))}\big{)},$
$\displaystyle
C_{3}=C_{3}\big{(}|\\!|\vec{p}|\\!|_{L^{\infty}(0,t_{l};H^{2}(\Omega))}\big{)}.$

Proof. By using the Cauchy-Schwarz inequality and discrete inf-sup condition,
we have

$\displaystyle|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{h,l}|\\!|_{L^{2}(\Omega)}$
$\displaystyle\leq|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{h,l}-e_{\xi}^{h,l}|\\!|_{L^{2}(\Omega)}+|\\!|e_{\xi}^{h,l}|\\!|_{L^{2}(\Omega)}$
(A.41)
$\displaystyle\leq|\\!|\vec{\alpha}^{\mathsf{T}}e_{\vec{p}}^{h,l}-e_{\xi}^{h,l}|\\!|_{L^{2}(\Omega)}+\frac{2\mu}{\beta_{0}}|\\!|\epsilon(e_{\boldsymbol{u}}^{h,l})|\\!|_{L^{2}(\Omega)}.$

Combining (A.27), (A.41) and (A.10)-(A.12), we obtain (A.37), (A.38) and
(A.40). According to the first Korn’s inequality and (A.37), we derive (A.39).
$\hfill\qed$

## References

  *   * Hong et al. [2020] Hong, Q., Kraus, J., Lymbery, M., Wheeler, M.F.: Parameter-robust convergence analysis of fixed-stress split iterative method for multiple-permeability poroelasticity systems. Multiscale modeling & simulation 18(2), 916–941 (2020)
  * Lee et al. [2019] Lee, J.J., Piersanti, E., Mardal, K.-A., Rognes, M.E.: A mixed finite element method for nearly incompressible multiple-network poroelasticity. SIAM journal on scientific computing 41(2), 722–747 (2019)
  * Bai et al. [1993] Bai, M., Elsworth, D., Roegiers, J.-C.: Multiporosity/multipermeability approach to the simulation of naturally fractured reservoirs. Water Resources Research 29(6), 1621–1633 (1993)
  * Nordbotten et al. [2010] Nordbotten, J.M., Rahman, T., Repin, S.I., Valdman, J.: A posteriori error estimates for approximate solutions of the Barenblatt-Biot poroelastic model. Computational methods in applied mathematics 10(3), 302–314 (2010)
  * Showalter and Momken [2002] Showalter, R., Momken, B.: Single-phase flow in composite poroelastic media. Mathematical methods in the applied sciences 25(2), 115–139 (2002)
  * Tully and Ventikos [2011] Tully, B., Ventikos, Y.: Cerebral water transport using multiple-network poroelastic theory: application to normal pressure hydrocephalus. Journal of Fluid Mechanics 667, 188–215 (2011)
  * Guo et al. [2018] Guo, L., Vardakis, J.C., Lassila, T., Mitolo, M., Ravikumar, N., Chou, D., Lange, M., Sarrami-Foroushani, A., Tully, B.J., Taylor, Z.A., et al.: Subject-specific multi-poroelastic model for exploring the risk factors associated with the early stages of Alzheimer’s disease. Interface focus 8(1), 20170019 (2018)
  * Both et al. [2017] Both, J.W., Borregales, M., Nordbotten, J.M., Kumar, K., Radu, F.A.: Robust fixed stress splitting for Biot’s equations in heterogeneous media. Applied Mathematics Letters 68, 101–108 (2017)
  * Mikelić et al. [2014] Mikelić, A., Wang, B., Wheeler, M.F.: Numerical convergence study of iterative coupling for coupled flow and geomechanics. Computational Geosciences 18, 325–341 (2014)
  * Cai and Zhang [2017] Cai, M., Zhang, G.: Comparisons of some iterative algorithms for Biot equations. International journal of evolution equations 10(3-4), 267 (2017)
  * Lee et al. [2017] Lee, J.J., Mardal, K.-A., Winther, R.: Parameter-robust discretization and preconditioning of Biot’s consolidation model. SIAM Journal on Scientific Computing 39(1), 1–24 (2017)
  * Phillips and Wheeler [2009] Phillips, P.J., Wheeler, M.F.: Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach. Computational Geosciences 13, 5–12 (2009)
  * Rodrigo et al. [2016] Rodrigo, C., Gaspar, F., Hu, X., Zikatanov, L.: Stability and monotonicity for some discretizations of the Biot’s consolidation model. Computer Methods in Applied Mechanics and Engineering 298, 183–204 (2016)
  * Rodrigo et al. [2018] Rodrigo, C., Hu, X., Ohm, P., Adler, J.H., Gaspar, F.J., Zikatanov, L.: New stabilized discretizations for poroelasticity and the Stokes’ equations. Computer Methods in Applied Mechanics and Engineering 341, 467–484 (2018)
  * Hong and Kraus [2017] Hong, Q., Kraus, J.: Parameter-robust stability of classical three-field formulation of Biot’s consolidation model. arXiv preprint arXiv:1706.00724 (2017)
  * Hu et al. [2017] Hu, X., Rodrigo, C., Gaspar, F.J., Zikatanov, L.T.: A nonconforming finite element method for the Biot’s consolidation model in poroelasticity. Journal of Computational and Applied Mathematics 310, 143–154 (2017)
  * Ju et al. [2020] Ju, G., Cai, M., Li, J., Tian, J.: Parameter-robust multiphysics algorithms for Biot model with application in brain edema simulation. Mathematics and computers in simulation 177, 385–403 (2020)
  * Lee [2016] Lee, J.J.: Robust error analysis of coupled mixed methods for Biot’s consolidation model. Journal of Scientific Computing 69(2), 610–632 (2016)
  * Oyarzúa and Ruiz-Baier [2016] Oyarzúa, R., Ruiz-Baier, R.: Locking-free finite element methods for poroelasticity. SIAM Journal on Numerical Analysis 54(5), 2951–2973 (2016)
  * Qi et al. [2021] Qi, W., Seshaiyer, P., Wang, J.: Finite element method with the total stress variable for Biot’s consolidation model. Numerical Methods for Partial Differential Equations 37(3), 2409–2428 (2021)
  * Wang et al. [2022] Wang, F., Cai, M., Wang, G., Zeng, Y.: A mixed virtual element method for Biot’s consolidation model. Computers $\&$ Mathematics with Applications 126, 31–42 (2022)
  * Yi and Bean [2017] Yi, S.-Y., Bean, M.L.: Iteratively coupled solution strategies for a four-field mixed finite element method for poroelasticity. International Journal for Numerical and Analytical Methods in Geomechanics 41(2), 159–179 (2017)
  * Hong et al. [2019] Hong, Q., Kraus, J., Lymbery, M., Philo, F.: Conservative discretizations and parameter-robust preconditioners for Biot and multiple-network flux-based poroelasticity models. Numerical Linear Algebra with Applications 26(4), 2242 (2019)
  * Hong et al. [2020] Hong, Q., Kraus, J., Lymbery, M., Philo, F.: Parameter-robust uzawa-type iterative methods for double saddle point problems arising in biot’s consolidation and multiple-network poroelasticity models. Mathematical Models and Methods in Applied Sciences 30(13), 2523–2555 (2020)
  * Lee [2019] Lee, J.J.: Unconditionally stable second order convergent partitioned methods for multiple-network poroelasticity. arXiv preprint arXiv:1901.06078 (2019)
  * Gu et al. [2023] Gu, H., Cai, M., Li, J.: An iterative decoupled algorithm with unconditional stability for biot model. Mathematics of Computation 92(341), 1087–1108 (2023)
  * Taylor and Hood [1973] Taylor, C., Hood, P.: A numerical solution of the Navier-Stokes equations using the finite element technique. Computers $\&$ Fluids 1(1), 73–100 (1973)
  * Hecht [2012] Hecht, F.: New development in FreeFem++. Journal of numerical mathematics 20(3-4), 251–266 (2012)
  * Ahmed et al. [2018] Ahmed, N., Linke, A., Merdon, C.: On really locking-free mixed finite element methods for the transient incompressible Stokes equations. SIAM Journal on Numerical Analysis 56(1), 185–209 (2018)