# Coupling nonconforming and enriched Galerkin methods for robust
discretization and fast solvers of poroelasticity problems

Jeonghun J. Lee1, Jacob Moore2 1,2 Department of Mathematics, Baylor
University, Waco, Texas, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Date: April, 2023)

###### Abstract.

In this paper we propose a new finite element discretization for the two-field
formulation of poroelasticity which uses the elastic displacement and the pore
pressure as primary variables. The main goal is to develop a numerical method
with small problem sizes which still achieve key features such as parameter-
robustness, local mass conservation, and robust preconditionor construction.
For this we combine a nonconforming finite element and the interior over-
stabilized enriched Galerkin methods with a suitable stabilization term.
Robust a priori error estimates and parameter-robust preconditioner
construction are proved, and numerical results illustrate our theoretical
findings.

###### Key words and phrases:

Biot’s consolidation model, error analysis, preconditioning

###### 2000 Mathematics Subject Classification:

Primary: 65N12, 65N15

##  1\. Introduction

The poroelasticity theory is a theory on fluid flows in deformable porous
media. The interactive behavior of fluid flows and elastic deformation results
in a system of partial differential equations coupling elasticity and porous
media flow equations. Since the first poroelasticity models were derived by
Biot [3, 4], mathematical properties and numerical methods for the models have
been widely explored. We refer to [38, 35] for well-posedness of linear quasi-
static poroelasticity problems. For numerical methods there are a number of
previous results with various formulations and finite element spaces [33, 34,
45, 26, 28, 29, 28, 40, 41]. The main focus of more recent studies is
developing numerical methods which provide good features of numerical
solutions. The features of interest are robustness of numerical solutions for
nearly incompressible solid matrix and saturating fluids, local mass
conservation, and fast solvers which show robust performance for realistic
physical parameter ranges. Here we list some recent results [42, 16, 17, 32,
12, 5, 19, 27, 9, 13, 11, 31, 10, 14, 6, 22] but we admit that this list of
references is by no means complete.

The goal of this paper is to develop a numerical method for poroelasticity
with small problem sizes such that numerical solutions satisfy the three
features: 1) the error bounds are robust for a nearly incompressible elastic
matrix, for arbitrarily small constrained specific storage coefficients, for
time interval, 2) local mass conservation of numerical solutions holds, 3)
preconditioners robust for mesh refinements as well as the physical parameter
values are available. For this goal we use a formulation which has the
displacement of poroelastic matrix $\boldsymbol{u}$ and the fluid pressure $p$
as unknowns.

We remark that there are a couple of previous works on robust numerical
methods using the two-field formulation. In [10], primal hybridized
discontinuous Galerkin methods are used for the displacement and pressure, and
it is proved that the methods give robust numerical solutions for nearly
incompressible materials. In [22], the enriched Galerkin method for linear
elasticty in [43] and the enriched Galerkin method for the Poisson equations
in [21] are coupled. To the best of our knowledge, construction of parameter-
robust preconditioning for these methods are not available. In contrast, we
show that the numerical method in this paper provides parameter-robust
preconditioner construction in addition to robust error estimates and local
mass conservation.

For finite element discretization of the two variables we use the
Mardal–Tai–Winther (MTW) element, a nonconforming $H^{1}$ finite element
developed in [23] for $\boldsymbol{u}$ and the interior over-stabilized
enriched Galerkin method for the pore pressure $p$ in [18]. For this
combination of discretization methods, we prove the a priori error estimates
for semidiscrete solutions without Grönwall inequalities and show a
preconditioner construction guidance with the operator preconditioning
approach in [25].

The paper is organized as follows. In Section 2 we introduce symbols,
definitions, and the partial differential equation formulation of Biot’s
consolidation model. In Section 3, finite element methods for spatial
discretization and auxiliary results for the stability of discretization are
presented. The a priori error analysis of semidiscrete solutions is proved in
Section 4. An abstract form of parameter-robust preconditioners is proposed in
Section 5, and numerical results for convergence of errors and for efficiency
of preconditioners are presented in Section 6. Conclusive remarks are given in
Section 7.

##  2\. Preliminaries

###  2.1. Notations

Let $\Omega$ be a bounded domain in ${\mathbb{R}}^{n}$, $n=2,3$ with a
polygonal or polyhedral boundary. For a nonnegative real number $m$,
$H^{m}(\Omega)$ and $H^{m}(\Omega;{\mathbb{R}}^{n})$ denote the standard
${\mathbb{R}}$ and ${\mathbb{R}}^{n}$-valued Sobolev spaces based on the
$L^{2}$ norm. For a Banach space $\mathcal{X}$ and $(a,b)\subset{\mathbb{R}}$,
$C^{0}(a,b;\mathcal{X})$ is the set of functions
$f:(a,b)\rightarrow\mathcal{X}$ which are continuous in $t\in(a,b)$. For an
integer $m\geq 1$ we define

$\displaystyle C^{m}(a,b;\mathcal{X}):=\\{f\,|\,\partial^{i}f/\partial
t^{i}\in C^{0}(a,b\,;\mathcal{X}),\,0\leq i\leq m\\},$

where $\partial f/\partial t$ is the Fréchet derivative of $f$ (see e.g.,
[44]). Recall the Bochner spaces by the norm

$\displaystyle\|f\|_{L^{p}(a,b;\mathcal{X})}=\begin{cases}\left(\int_{a}^{b}\|f\|_{\mathcal{X}}^{p}ds\right)^{1/p},\quad
1\leq p<\infty,\\\ \operatorname{esssup}_{t\in(a,b)}\|f\|_{\mathcal{X}},\quad
p=\infty.\end{cases}$

If a time interval $J\subset{\mathbb{R}}$ is clear in context, then we use
$L^{p}\mathcal{X}$ to denote $L^{p}(J;\mathcal{X})$ for simplicity. We define
the Bochner-Sobolev spaces $W^{k,p}(J;\mathcal{X})$ for nonnegative integer
$k$ and $1\leq p\leq\infty$ as the closure of $C^{k}(J;\mathcal{X})$ with the
norm $\|f\|_{W^{k,p}\mathcal{X}}=\sum_{i=0}^{k}\|\partial^{i}f/\partial
t^{i}\|_{L^{p}\mathcal{X}}$. We adopt the convention
$\|f,g\|_{\mathcal{X}}=\|f\|_{\mathcal{X}}+\|g\|_{\mathcal{X}}$ for
simplicity, and $\dot{u}$, $\ddot{u}$ is used to denote $\partial u/\partial
t$, $\partial^{2}u/\partial t^{2}$.

In this paper we use $\mathcal{T}_{h}$ to denote a shape-regular triangulation
of $\Omega$ for which $h$ is the maximum diameter of triangles/tetrahedra and
$\mathcal{E}_{h}$ is the corresponding set of edges/faces. For a set $G$ with
a positive $n$-dimensional Lebesgue measure and for functions $f,g\in
L^{2}(\Omega)$ we define $\left(f,{g}\right)_{G}$ by

$\displaystyle\left({f},{g}\right)_{G}:=\int_{G}{f}\cdot{g}\,dx.$

For a finite-dimensional vector space ${\mathbb{X}}$ on ${\mathbb{R}}$ and
functions $f,g\in L^{2}(\Omega;{\mathbb{X}})$, $\left({f},{g}\right)_{G}$ is
similarly defined using the natural inner product on ${\mathbb{X}}$. For
$e\in\mathcal{E}_{h}$ and functions ${f},{g}\in L^{2}(e)$ we define

$\displaystyle\langle{f},{g}\rangle_{e}=\int_{e}{f}\cdot{g}\,ds,\qquad\langle{f},{g}\rangle=\sum_{e\in\mathcal{E}_{h}}\langle{f},{g}\rangle_{e}.$

$\left<f,g\right>_{e}$ and $\left<f,g\right>$ are defined similarly for
$f,g\in L^{2}(e;{\mathbb{X}})$.

For an integer $k\geq 0$ and for each $K\in\mathcal{T}_{h}$,
$\mathcal{P}_{k}(K)$ is the space of polynomials of degree $\leq k$ on $K$,
and $\mathcal{P}_{k}(\mathcal{T}_{h})$ denotes the space

$\displaystyle\mathcal{P}_{k}(\mathcal{T}_{h})=\left\\{q\in
L^{2}(\Omega)\;:\;q|_{K}\in\mathcal{P}_{k}(T),\;K\in\mathcal{T}_{h}\right\\}.$

For a vector space ${\mathbb{X}}$, we use $\mathcal{P}_{k}(G;{\mathbb{X}})$
and $\mathcal{P}_{k}(\mathcal{T}_{h};{\mathbb{X}})$ to denote the space of
${\mathbb{X}}$-valued polynomials with same conditions.

For a nonnegative function (or a positive semidefinite tensor) $w$ on
$\Omega$, $\left\|q\right\|_{0,w}$ denotes
$\left\|q\right\|_{0,w}=\left(wq,q\right)_{\Omega}^{1/2}$.

###  2.2. The Biot’s consolidation model

Throughout this paper we limit our interest on quasistatic consolidation
problems. In our formulation of the model, $\boldsymbol{u}$ is the
displacement of porous media, $p$ is the pore pressure, $\boldsymbol{f}$ is
the body force, $g$ is the mass change rate of fluid. The governing equations
of the quasistatic consolidation model with an isotropic elastic porous medium
are

(1a)
$\displaystyle-\operatorname{div}\left(2\mu\epsilon(\boldsymbol{u})+(\lambda\operatorname{div}\boldsymbol{u}-\alpha
p)\mathbb{I}\right)$ $\displaystyle=\boldsymbol{f},$ (1b) $\displaystyle
s_{0}\dot{p}+\alpha\operatorname{div}\dot{\boldsymbol{u}}-\operatorname{div}(\underline{\boldsymbol{\kappa}}\nabla
p)$ $\displaystyle=g,$

where $\mu$ and $\lambda$ are the Lamé coefficients, $s_{0}\geq 0$ is the
constrained specific storage coefficient, $\underline{\boldsymbol{\kappa}}$ is
the hydraulic conductivity tensor, $\alpha>0$ is the Biot–Willis coefficient
which is close to 1, and $\mathbb{I}$ is the identity matrix.

We assume that $\mu$ is uniformly bounded above and below with positive
constants. We assume $\lambda$ has a uniformly positive lower bound but
$\lambda$ may not have a uniform upper bound and $\lambda=+\infty$ corresponds
to the incompressibility of the solid matrix. We assume that there are
constants $c_{0},c_{1}$ such that

$\displaystyle 0\leq c_{0}\leq s_{0}(x)\leq c_{1},\qquad x\in\Omega.$

We remark that $s_{0}$ is related to $\alpha$, the porosity $\phi$, and the
bulk moduli of the solid and fluid. Under the assumption that $\phi$ is
uniform with $0<\phi<\alpha$, if the solid is not incompressible, then
$s_{0}\geq C/\lambda$ holds with a constant $C$ of scale 1. However, $s_{0}$
may vanish on a subdomain if $\lambda=+\infty$ on the subdomain and the fluid
is incompressible. The hydraulic conductivity tensor
$\underline{\boldsymbol{\kappa}}=\underline{\boldsymbol{\kappa}}(x)$ is
positive definite with uniform lower and upper bounds
$\kappa_{0},\kappa_{1}>0$, i.e.,

$\displaystyle\kappa_{0}|\xi|^{2}\leq\xi^{T}\underline{\boldsymbol{\kappa}}(x)\xi\leq\kappa_{1}|\xi|^{2},\qquad\forall\;0\not=\xi\in{\mathbb{R}}^{n},\quad\text{a.e.}\;x\in\Omega.$

On details of deriving these equations from physical modelling, we refer to
standard references on poroelasticity theory, e.g., [39].

For well-posedness of the problem, the equations (1) need appropriate boundary
and initial conditions. We assume that there are partitions of
$\partial\Omega$ which are

$\displaystyle\partial\Omega=\Gamma_{p}\cup\Gamma_{f},\qquad\partial\Omega=\Gamma_{d}\cup\Gamma_{t},\qquad|\Gamma_{d}|,|\Gamma_{p}|>0$

where $|\Gamma|$ is the $(n-1)$-dimensional Lebesgue measure of $\Gamma$. We
also assume that boundary conditions are given as

(2) $\displaystyle p(t)=0\text{ on
}\Gamma_{p},\quad-\underline{\boldsymbol{\kappa}}\nabla
p(t)\cdot\boldsymbol{n}=0\text{ on }\Gamma_{f},\quad\boldsymbol{u}(t)=0\text{
on }\Gamma_{d},\quad\underline{\boldsymbol{\sigma}}(t)\boldsymbol{n}=0\text{
on }\Gamma_{t},$

for all $t\in(0,T]$ where $\boldsymbol{n}$ is the outward unit normal vector
field on $\partial\Omega$ and
$\underline{\boldsymbol{\sigma}}:=2\mu\epsilon(\boldsymbol{u})+(\lambda\operatorname{div}\boldsymbol{u}-\alpha
p)\mathbb{I}$, the Cauchy stress tensor. Here we only consider the homogeneous
boundary condition for simplicity but our method can be easily extended to
problems with nonhomogeneous boundary conditions. We also assume that given
initial data $p(0),\boldsymbol{u}(0)$ and $\boldsymbol{f}(0)$ satisfy the
compatibility condition (1a). Well-posedness of this system under these
assumptions can be found in [36].

##  3\. Finite element discretization

In this section we present a variational formulation and finite element
discretization of (1), and prove auxiliary results for the stability and error
analyses.

###  3.1. Variational formulation of continuous problem

Let us define function spaces

$\displaystyle\boldsymbol{V}=\left\\{\boldsymbol{v}\in{H}^{1}(\Omega;{\mathbb{R}}^{n})\;:\;\boldsymbol{v}|_{\Gamma_{d}}=0\right\\},\quad
Q=\\{q\in H^{1}(\Omega)\;:\;q|_{\Gamma_{p}}=0\\}.$

Then, by the integration by parts, we can show that a solution of (1),
$(\boldsymbol{u}(t),p(t))$, satisfies the following variational formulation:

(3a)
$\displaystyle\left(2\mu\epsilon(\boldsymbol{u}(t)),\epsilon(\boldsymbol{v})\right)_{\Omega}+\left(\lambda\operatorname{div}\boldsymbol{u}(t),\operatorname{div}\boldsymbol{v}\right)_{\Omega}-\left(\alpha
p(t),\operatorname{div}\boldsymbol{v}\right)_{\Omega}$
$\displaystyle=\left(\boldsymbol{f}(t),\boldsymbol{v}\right)_{\Omega},$ (3b)
$\displaystyle-\left(\alpha\operatorname{div}\dot{\boldsymbol{u}}(t),q\right)_{\Omega}-\left(s_{0}\dot{p}(t),q\right)_{\Omega}-\left(\underline{\boldsymbol{\kappa}}\nabla
p(t),\nabla q\right)_{\Omega}$ $\displaystyle=-\left(g(t),q\right)_{\Omega}$

for all $(\boldsymbol{v},q)\in\boldsymbol{V}\times Q$ and for all $t\geq 0$.

For well-posedness we need initial conditions
$(\boldsymbol{u}(0),p(0))\in\boldsymbol{V}\times Q$. By the compatibility
assumption,

(4)
$\displaystyle\left(2\mu\epsilon(\boldsymbol{u}(0)),\epsilon(\boldsymbol{v})\right)_{\Omega}+\left(\lambda\operatorname{div}\boldsymbol{u}(0),\operatorname{div}\boldsymbol{v}\right)_{\Omega}-\left(\alpha
p(0),\operatorname{div}\boldsymbol{v}\right)_{\Omega}$
$\displaystyle=\left(\boldsymbol{f}(0),\boldsymbol{v}\right)_{\Omega}$

for all $\boldsymbol{v}\in\boldsymbol{V}$.

In summary the variational problem of (1) is the following: For given initial
condition $(\boldsymbol{u}(0),p(0))\in\boldsymbol{V}\times Q$ satisfying (4)
find $(\boldsymbol{u},p)\in C^{1}(0,T;\boldsymbol{V})\times C^{1}(0,T;Q)$ such
that (3) holds for all $t\geq 0$.

###  3.2. Finite elements for spatial discretization

In this subsection we introduce finite element methods for spatial
discretizations of $\boldsymbol{u}$ and $p$.

For discretization of $\boldsymbol{u}$ we use the vector-valued nonconforming
$H^{1}$ elements by Mardal, Tai, Winther ([23, 37]). Here we only summarize
key features of the finite elements, and we refer to the original articles for
detailed definitions of the finite elements.

Let $\boldsymbol{V}_{h}$ be the subspace of the Mardal-Tai-Winther element
with an additional condition that

(5) $\displaystyle\text{all DOFs associated to the edges/faces in
}\Gamma_{d}\text{ vanish}.$

By definition one can see that $\boldsymbol{V}_{h}\subset
H(\operatorname{div},\Omega)$ and
$\operatorname{div}\boldsymbol{V}_{h}=\mathcal{P}_{0}(\mathcal{T}_{h})$ but
$\boldsymbol{V}_{h}\not\subset H^{1}(\Omega;{\mathbb{R}}^{n})$. If we define
$\|\boldsymbol{v}\|_{1,h}$ by
$\|\boldsymbol{v}\|_{1,h}^{2}:=\|\nabla\boldsymbol{v}\|_{0}^{2}+\|\boldsymbol{v}\|_{0}^{2}$
with $\nabla$ and $\epsilon$, by the element-wise gradient and symmetric
gradient operators, then a discrete Korn inequality

$\displaystyle\|\boldsymbol{v}\|_{1,h}\geq C\|\epsilon(\boldsymbol{v})\|_{0}$

holds for $C>0$ independent of $h$ (see [24]). In addition, an inf-sup
condition

(6)
$\displaystyle\inf_{q\in\mathcal{P}_{0}(\mathcal{T}_{h})}\sup_{\boldsymbol{v}\in\boldsymbol{V}_{h}}\frac{\left(\operatorname{div}\boldsymbol{v},q\right)_{\Omega}}{\|\boldsymbol{v}\|_{1,h}\|q\|_{0}}\geq
C>0$

holds with $C$ independent of $h$.

For discretization of $p$ we use the interior over-stabilized enriched
Galerkin methods developed in [18]. We first define $Q_{h}^{c}$, $Q_{h}^{0}$,
$Q_{h}$ by

$\displaystyle Q_{h}^{c}:=\mathcal{P}_{1}(\mathcal{T}_{h})\cap
H^{1}(\Omega),\quad Q_{h}^{0}:=\mathcal{P}_{0}(\mathcal{T}_{h}),\quad
Q_{h}:=Q_{h}^{c}+Q_{h}^{0}.$

Let use define
$a_{\boldsymbol{u}}:(\boldsymbol{V}+\boldsymbol{V}_{h})\times(\boldsymbol{V}+\boldsymbol{V}_{h})\rightarrow\mathbb{R}$
and $\|\cdot\|_{\boldsymbol{V}_{h}}$ by

(7) $\displaystyle
a_{\boldsymbol{u}}\left(\boldsymbol{v},\tilde{\boldsymbol{v}}\right)$
$\displaystyle:=\left(2\mu\epsilon(\boldsymbol{v}),\epsilon(\tilde{\boldsymbol{v}})\right)+\left(\lambda\operatorname{div}\boldsymbol{v},\operatorname{div}\tilde{\boldsymbol{v}}\right),\qquad\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}=\left(a_{\boldsymbol{u}}\left(\boldsymbol{v},\boldsymbol{v}\right)\right)^{1/2}.$

To define the bilinear form for $p$ with the interior over-stabilized enriched
Galerkin method, we define the jump and average operators for functions which
have well-defined traces on edges/faces. More specifically, for
$e\in\mathcal{E}_{h}^{\partial}$, we define

$\displaystyle\left[\\!\left[{q}\right]\\!\right]|_{e}=q|_{e}\boldsymbol{n},\quad\left\\{\\!\\!\left\\{{q}\right\\}\\!\\!\right\\}|_{e}=q|_{e},\qquad\left[\\!\left[{\boldsymbol{q}}\right]\\!\right]|_{e}=\boldsymbol{q}|_{e}\cdot\boldsymbol{n},\quad\left\\{\\!\\!\left\\{{\boldsymbol{q}}\right\\}\\!\\!\right\\}|_{e}=\boldsymbol{q}|_{e}$

where $q$ and $\boldsymbol{q}$ are scalar and ${\mathbb{R}}^{n}$-vector-valued
functions such that their traces on edge/face $e$ are well-defined as
integrable functions, and $\boldsymbol{n}$ is the outward unit normal vector
field on $e$. For $e\in\mathcal{E}_{h}^{0}$, $T_{+}$ and $T_{-}$ are the two
elements such that $e=\partial T_{+}\cap\partial T_{-}$. If $q$ is a scalar
function on $T_{+}\cup T_{-}$, then we use $q^{\pm}$ and $\nabla q^{\pm}$ to
denote the restrictions of $q$ and $\nabla q$ on $T_{+}$ and $T_{-}$. We use
$\boldsymbol{n}^{\pm}$ to denote the unit outward normal vector fields of
$T_{\pm}$ and we again omit the restriction $|_{e}$ if it is clear in context.
Then the jumps and the averages of $q$ and $\boldsymbol{q}$ on
$e\in\mathcal{E}_{h}^{i}$ are defined by

$\displaystyle\left[\\!\left[{q}\right]\\!\right]|_{e}$
$\displaystyle=q^{+}\boldsymbol{n}^{+}+q^{-}\boldsymbol{n}^{-},$
$\displaystyle\left\\{\\!\\!\left\\{{q}\right\\}\\!\\!\right\\}|_{e}$
$\displaystyle=\frac{1}{2}(q^{+}+q^{-}),$
$\displaystyle\left[\\!\left[{\boldsymbol{q}}\right]\\!\right]|_{e}$
$\displaystyle=\boldsymbol{q}^{+}\cdot\boldsymbol{n}^{+}+\boldsymbol{q}^{-}\cdot\boldsymbol{n}^{-},$
$\displaystyle\left\\{\\!\\!\left\\{{\boldsymbol{q}}\right\\}\\!\\!\right\\}|_{e}$
$\displaystyle=\frac{1}{2}(\boldsymbol{q}^{+}+\boldsymbol{q}^{-}).$

If $e\in\mathcal{E}_{h}$ is clear in context, we will use
$\left[\\!\left[{q}\right]\\!\right]$ instead of
$\left[\\!\left[{q}\right]\\!\right]|_{e}$ for the jump of $q$ on
$e\in\mathcal{E}_{h}$. The same simplification will apply to other quantities.

For $q\in Q_{h}$ we use $q^{c}$ and $q^{0}$ to denote the components in
$Q_{h}^{c}$ and $Q_{h}^{0}$ which are unique up to a constant addition. For
$q,\tilde{q}\in Q_{h}$ we now define

(8)
$a_{p}\left(q,\tilde{q}\right):=\left(\underline{\boldsymbol{\kappa}}\nabla
q^{c},\nabla\tilde{q^{c}}\right)_{\Omega}-\langle\underline{\boldsymbol{\kappa}}\left\\{\\!\\!\left\\{{\nabla
q^{c}}\right\\}\\!\\!\right\\},\left[\\!\left[{\tilde{q}^{0}}\right]\\!\right]\rangle-\langle\left[\\!\left[{q^{c}}\right]\\!\right],\left\\{\\!\\!\left\\{{\underline{\boldsymbol{\kappa}}\nabla\tilde{q}}\right\\}\\!\\!\right\\}\rangle\\\
+\gamma\underline{\boldsymbol{\kappa}}_{avg}\sum_{e\in\mathcal{E}_{h}^{0}}\langle
h_{e}^{-1-\beta}\left[\\!\left[{q^{0}}\right]\\!\right],\left[\\!\left[{\tilde{q}^{0}}\right]\\!\right]\rangle_{e}+\gamma\underline{\boldsymbol{\kappa}}_{avg}\sum_{e\in\mathcal{E}_{h}^{\partial}}\langle
h_{e}^{-1}(q^{c}+q^{0}),(\tilde{q}^{c}+\tilde{q}^{0})\rangle_{e}$

with $\beta\geq 1$ and a sufficiently large $\gamma>0$. This bilinear form is
well-defined independent of the choices of $q^{c},\tilde{q}^{c}$, $q^{0}$,
$\tilde{q}^{0}$ up to constant addition because all terms except the last term
in this bilinear form vanishes if $q^{c},\tilde{q}^{c}$, $q^{0}$,
$\tilde{q}^{0}$ are constant functions on $\Omega$, and the last term is
independent of the decompositions of
$q^{c},q^{0},\tilde{q}^{c},\tilde{q}^{0}$. Note that there exists a unique
decomposition $q=q^{c}+q^{0}\in Q_{h}^{c}+Q_{h}^{0}$ for $q\in Q_{h}$ such
that $q^{0}\in Q_{h}^{0}$ is mean-value zero on $\Omega$. We will use this
decomposition frequently in the rest of this paper for analysis. If $\gamma>0$
is sufficiently large depending on the shape regularity of the triangulation
$\mathcal{T}_{h}$, then $a_{p}(\cdot,\cdot)$ is a coercive bilinear form on
$Q_{h}$. Note that (8) is the enriched Galerkin (EG) method for second order
elliptic equations in [21] if $\beta=0$. In this paper we assume $\beta\geq 1$
which is the interior over-stabilized enriched Galerkin methods (IOS-EG) in
[18]. Convergence rates of numerical solutions are same in EG and IOS-EG
methods. However, as is discussed in [18], IOS-EG methods provide more robust
fast solver algorithms. In this paper we will also show this by numerical
experiments later.

We also define a stabilization bilinear form $S:Q_{h}\times
Q_{h}\rightarrow{\mathbb{R}}$ by

(9) $\displaystyle
S(q,\tilde{q})=\gamma\left(\lambda^{-1}\sum_{e\in\mathcal{E}_{h}^{0}}\langle
h_{e}^{-1}\left[\\!\left[{q}\right]\\!\right],\left[\\!\left[{\tilde{q}}\right]\\!\right]\rangle_{e}+\sum_{e\in\mathcal{E}_{h}^{0}}\langle
C_{1}h_{e}^{-1}\left[\\!\left[{q}\right]\\!\right],\left[\\!\left[{\tilde{q}}\right]\\!\right]\rangle_{e}\right)$

with a sufficiently large $\gamma>0$. This stabilization term will be used to
discretize the time derivative term of $p$ in order to guarantee the stability
of discretization independent of the lower bound of $s_{0}\geq 0$ and robust
preconditioning. Note that the edge/face integrals in $S(\cdot,\cdot)$ are
defined only on the interior edges/faces $\mathcal{E}_{h}^{0}$, so only mean-
value zero components of $q^{0}$ and $\tilde{q}^{0}$ contribute to
$S(q,\tilde{q})$. We will use this fact in our stability analysis later.

Let us define the norm $\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}$ for
$\boldsymbol{v}\in\boldsymbol{V}_{h}$ by

$\displaystyle\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}^{2}:=\int_{\Omega}\left(2\mu\epsilon(\boldsymbol{v}):\epsilon(\boldsymbol{v})+\lambda(\operatorname{div}\boldsymbol{v})^{2}\right)\,dx$

where $\epsilon(\boldsymbol{v}):\epsilon(\boldsymbol{v})$ is the Frogenius
inner product on matrices. For $q\in Q_{h}$ we define $\|q\|_{0,\lambda^{-1}}$
and $\|q\|_{h}$ by

(10) $\displaystyle\|q\|_{0,\lambda^{-1}}^{2}$
$\displaystyle:=\int_{\Omega}\lambda^{-1}q^{2}\,dx,$ (11)
$\displaystyle\|q\|_{h}^{2}$
$\displaystyle:=\int_{\Omega}\underline{\boldsymbol{\kappa}}\nabla
q^{c}\cdot\nabla
q^{c}\,dx+\sum_{e\in\mathcal{E}_{h}^{0}}\gamma\underline{\boldsymbol{\kappa}}_{avg}h_{e}^{-1-\beta}\int_{e}\left[\\!\left[{q}\right]\\!\right]^{2}\,ds$
$\displaystyle\quad+\sum_{e\in\mathcal{E}_{h}^{\partial}}\gamma\underline{\boldsymbol{\kappa}}_{avg}h_{e}^{-1}\int_{e}\left[\\!\left[{q^{c}}\right]\\!\right]^{2}\,ds+\sum_{e\in\mathcal{E}_{h}^{\partial}}\gamma\underline{\boldsymbol{\kappa}}_{avg}h_{e}^{-1}\int_{e}\left[\\!\left[{q^{0}}\right]\\!\right]^{2}\,ds.$

In the rest of this subsection we introduce a few results for stability of the
discretization. These results will be used to prove well-posedness as well as
the a priori error analysis of the discretized system.

We begin with defining an auxiliary finite element space $R_{h}$ and a linear
operator $I_{h}:Q_{h}^{c}\rightarrow R_{h}$. For a vertex $v_{i}$ of
$\mathcal{T}_{h}$ let $\chi_{i}$ be the characteristic function on the support
of the shape function of an element in $Q_{h}^{c}$ associated to the vertex
degree of freedom at $v_{i}$. Let $R_{h}$ be the subspace of $Q_{h}^{0}$
defined by

$\displaystyle R_{h}=\\{q^{0}\in
Q_{h}^{0}\,:\,q^{0}=\sum_{i=1}^{N}c_{i}\chi_{i},c_{i}\in{\mathbb{R}}\\}$

where $N$ is the number of vertices of $\mathcal{T}_{h}$. We define
$I_{h}:Q_{h}^{c}\rightarrow R_{h}$ by

$\displaystyle I_{h}q^{c}=\sum_{i=1}^{N}q^{c}(v_{i})\chi_{i}.$

The following result can be obtained with the same argument in [15, Lemma 1].

######  Lemma 3.1.

Given $\boldsymbol{v}\in\boldsymbol{V}_{h}$ and $q^{c}\in Q_{h}^{c}$

(12)
$\displaystyle(n+1)\left(\operatorname{div}\boldsymbol{v},q^{c}\right)_{\Omega}=\left(\operatorname{div}\boldsymbol{v},I_{h}q^{c}\right)_{\Omega}.$

###### Proof.

Since $\operatorname{div}\boldsymbol{V}_{h}=Q_{h}^{0}$ we can see that (12)
holds for any nodal basis $q^{c}\in Q_{h}^{c}$, so the conclusion follows. ∎

The following norm equivalence result is also proved in [15]. The proof
follows from a standard scaling argument, so we state the result here without
proof.

######  Lemma 3.2.

There exist $C_{1},C_{2}>0$ independent of $h$ such that

$\displaystyle C_{1}\|q^{c}\|_{0}\leq\|I_{h}q^{c}\|_{0}\leq
C_{2}\|q^{c}\|_{0},\quad q^{c}\in Q_{h}^{c}.$

By combining Lemma 3.1 and Lemma 3.2, we can prove an inf-sup condition with
$\lambda$-weighted norms.

######  Lemma 3.3.

There exist $C>0$ independent of $h$ such that

$\displaystyle\inf_{q^{c}\in
Q_{h}^{c}}\sup_{\boldsymbol{v}\in\boldsymbol{V}_{h}}\frac{\left(\operatorname{div}\boldsymbol{v},q^{c}\right)_{\Omega}}{\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}\|q^{c}\|_{0,\lambda^{-1}}}\geq
C.$

###### Proof.

The inf-sup condition

$\displaystyle\inf_{q^{c}\in
Q_{h}^{c}}\sup_{\boldsymbol{v}\in\boldsymbol{V}_{h}}\frac{\left(\operatorname{div}\boldsymbol{v},q^{c}\right)_{\Omega}}{\|\boldsymbol{v}\|_{1,h}\|q^{c}\|_{0}}\geq
C$

follows immediately by Lemma 3.1 and Lemma 3.2. As an equivalent condition of
this inf-sup condition, for any given $q^{c}\in Q_{h}^{c}$ there exists
$\boldsymbol{w}\in\boldsymbol{V}_{h}$ such that

$\displaystyle\left(\operatorname{div}\boldsymbol{w},\tilde{q}^{c}\right)_{\Omega}=\left(q^{c},\tilde{q}^{c}\right)_{\Omega}\quad\forall\tilde{q}^{c}\in
Q_{h}^{c},\qquad\|\boldsymbol{w}\|_{1,h}\leq C_{1}\|q^{c}\|_{0}$

with $C_{1}>0$ independent of $h$. We define
$\boldsymbol{v}=\lambda^{-1}\boldsymbol{w}$. Then,

(13)
$\displaystyle\left(\operatorname{div}\boldsymbol{v},q^{c}\right)_{\Omega}=\int_{\Omega}\lambda^{-1}(q^{c})^{2}\,dx=\|q^{c}\|_{0,\lambda^{-1}}^{2}.$

In addition, since $\mu\leq C\lambda$ with a uniform constant $C>0$, we can
obtain that

$\displaystyle\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}^{2}$
$\displaystyle=\int_{\Omega}\left(2\mu\epsilon(\boldsymbol{v}):\epsilon(\boldsymbol{v})+\lambda(\operatorname{div}\boldsymbol{v})^{2}\right)\,dx$
$\displaystyle=\lambda^{-2}\int_{\Omega}\left(2\mu\epsilon(\boldsymbol{w}):\epsilon(\boldsymbol{w})+\lambda(\operatorname{div}\boldsymbol{w})^{2}\right)\,dx$
$\displaystyle\leq\left(2C_{\mu}+1\right)\lambda^{-1}\|\boldsymbol{w}\|_{1,h}^{2}$
$\displaystyle\leq\left(2C_{\mu}+1\right)C_{1}\|q^{c}\|_{0,\lambda^{-1}}^{2}.$

Therefore, the assertion holds by (13) and the above estimate. ∎

By a discrete Poincare inequality in [7] we can obtain the following lemma.

######  Lemma 3.4.

There exists $C_{P}>0$ independent of $h$ such that

(14) $\displaystyle\|q^{0}\|_{0,s_{0}}+\|q^{0}\|_{0,\lambda^{-1}}\leq
C_{P}S(q^{0},q^{0})^{\frac{1}{2}}.$

holds for all $q^{0}\in Q_{h}^{0}$ with $\int_{\Omega}q^{0}\,dx=0$.

We introduce another interpolation operator $I_{h}^{c}:Q_{h}^{0}\rightarrow
Q_{h}^{c}$ for analysis.

######  Lemma 3.5.

There exists $I_{h}^{c}:Q_{h}^{0}\to Q_{h}^{c}$, a linear operator, such that

(15) $\displaystyle\|q^{0}-I_{h}^{c}q^{0}\|_{0,\lambda^{-1}}\leq
CS\left(q,q\right)^{\frac{1}{2}}$

with $C>0$ independent of $h$.

###### Proof.

We define $I_{h}^{c}$ by the reconstruction operator in [8] and consider
$I_{h}^{c}q^{0}$ for $q^{0}\in Q_{h}^{0}$. For a simplex $K$ in
$\mathcal{T}_{h}$ let $M_{K}$ be the union of all (closed) simplices in
$\mathcal{T}_{h}$ adjacent to $K$ and $\text{int}(M_{K})$ is the interior of
$M_{K}$. It is known that

$\displaystyle\int_{T}(q^{0}-I_{h}^{c}q^{0})^{2}\,dx\leq
Ch_{K}^{2}\sum_{e\subset\text{int}(M_{K})}h_{e}^{-1}\int_{e}\left[\\!\left[{q^{0}}\right]\\!\right]|_{e}^{2}\,ds$

with $C$ independent of the mesh sizes [8, Theorem 6]. By the shape regularity
assumption $h_{K}\sim h_{e}$ in the above inequality with constants
independent of the mesh sizes, the conclusion follows because
$\\{M_{K}\\}_{K\in\mathcal{T}_{h}}$ covers $\Omega$ with only uniformly
bounded number of finite overlapping. ∎

######  Lemma 3.6.

There exists $C_{1},C_{2}>0$ independent of $h$ satisfying the following: For
any given $q\in Q_{h}$ one can find
$\boldsymbol{v}\not=0\in\boldsymbol{V}_{h}$ satisfying

(16)
$\displaystyle\frac{\left(\operatorname{div}\boldsymbol{v},q\right)_{\Omega}}{\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}}\geq
C_{1}\|q\|_{0,\lambda^{-1}}-C_{2}S\left(q,q\right)^{\frac{1}{2}}.$

###### Proof.

Let $q=q^{c}+q^{0}\in Q_{h}^{c}+Q_{h}^{0}$ be given and consider the
decomposition $q=q^{c}+I_{h}^{c}q^{0}-I_{h}^{c}q^{0}+q^{0}$. By the inf-sup
condition of $\boldsymbol{V}_{h}$ and $Q_{h}^{c}$ in Lemma 3.7 and a scaling
with $\lambda$, there exists $\boldsymbol{w}\in\boldsymbol{V}_{h}$ such that

(17)
$\displaystyle\left(\operatorname{div}\boldsymbol{w},q^{c}+I_{h}^{c}q^{0}\right)_{\Omega}$
$\displaystyle=\|q^{c}+I_{h}^{c}q^{0}\|_{0,\lambda^{-1}}^{2},$ (18)
$\displaystyle\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}$ $\displaystyle\leq
C_{3}\|q^{c}+I_{h}^{c}q^{0}\|_{0,\lambda^{-1}}$

with a constant $C_{3}$ independent of $h$. By the triangle inequality

(19) $\displaystyle\|q^{c}+I_{h}^{c}q^{0}\|_{0,\lambda^{-1}}$
$\displaystyle=\|q+(I_{h}^{c}q^{0}-q^{0})\|_{0,\lambda^{-1}}$ (20)
$\displaystyle\leq\|q\|_{0,\lambda^{-1}}+\|I_{h}^{c}q^{0}-q^{c}\|_{0,\lambda^{-1}}.$

By the definition of $\boldsymbol{w}$ and the Cauchy-Schwarz inequality,

$\displaystyle\left(\operatorname{div}\boldsymbol{w},q\right)_{\Omega}$
$\displaystyle=\left(\operatorname{div}\boldsymbol{w},q^{c}+I_{h}^{c}q^{0}-I_{h}^{c}q^{0}+q^{0}\right)_{\Omega}$
$\displaystyle=\|q^{c}+I_{h}^{c}q^{0}\|_{\lambda^{-1}}^{2}-\left(\operatorname{div}\boldsymbol{w},I_{h}^{c}q^{0}-q^{0}\right)_{\Omega}$
$\displaystyle\geq\|q^{c}+I_{h}^{c}q^{0}\|_{\lambda^{-1}}^{2}-\|\operatorname{div}\boldsymbol{w}\|_{0,\lambda}\|I_{h}^{c}q^{0}-q^{0}\|_{0,\lambda^{-1}}$

Since
$\|\operatorname{div}\boldsymbol{w}\|_{0,\lambda}\leq\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}\leq
C_{3}\|q^{c}+I_{h}^{c}q^{0}\|_{0,\lambda^{-1}}$, dividing by
$\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}$ gives

$\displaystyle\frac{\left(\operatorname{div}\boldsymbol{w},q\right)_{\Omega}}{\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}}$
$\displaystyle\geq\frac{1}{C_{3}}\|q^{c}+I_{h}^{c}q^{0}\|_{0,\lambda^{-1}}-\|I_{h}^{c}q^{0}-q^{0}\|_{0,\lambda^{-1}}.$

By (19) and Lemma 3.5, we get

$\displaystyle\frac{\left(\operatorname{div}\boldsymbol{w},q\right)_{\Omega}}{\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}}$
$\displaystyle\geq\frac{1}{C_{3}}\|q\|_{0,\lambda^{-1}}-\frac{C_{3}+1}{C_{3}}\|I_{h}^{c}q^{0}-q^{0}\|_{0,\lambda^{-1}}$
$\displaystyle\geq\frac{1}{C_{3}}\|q\|_{0,\lambda^{-1}}-\frac{C_{4}(C_{3}+1)}{C_{3}}S\left(q,q\right)^{\frac{1}{2}}.$

The conclusion holds with $C_{1}=1/C_{3}$ and $C_{2}=C_{4}(C_{3}+1)/C_{3}$. ∎

Let $B:(\boldsymbol{V}_{h}\times Q_{h})\times(\boldsymbol{V}_{h}\times
Q_{h})\rightarrow{\mathbb{R}}$ be defined by

$\displaystyle B((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle=a_{\boldsymbol{u}}(\boldsymbol{v},\tilde{\boldsymbol{v}})-\left(\alpha
q,\operatorname{div}\tilde{\boldsymbol{v}}\right)_{\Omega}-\left(\alpha\tilde{q},\operatorname{div}\boldsymbol{v}\right)_{\Omega}-\left(s_{0}q,\tilde{q}\right)_{\Omega}-S\left(q,\tilde{q}\right).$

For simplicity we define $\|(\boldsymbol{v},q)\|_{\boldsymbol{V}_{h}\times
Q_{h}}$ by

$\displaystyle\|(\boldsymbol{v},q)\|_{\boldsymbol{V}_{h}\times
Q_{h}}:=\left(\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}^{2}+\|q\|_{0,\lambda^{-1}}^{2}+\|q\|_{0,s_{0}}^{2}\right)^{\frac{1}{2}}.$

Then, we prove the following stability result.

######  Lemma 3.7.

There exists $C>0$ independent of $h$ such that

$\displaystyle\inf_{(\tilde{\boldsymbol{v}},\tilde{q})\in\boldsymbol{V}_{h}\times
Q_{h}}\sup_{(\boldsymbol{v},q)\in\boldsymbol{V}_{h}\times
Q_{h}}\frac{B((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))}{\|(\boldsymbol{v},q)\|_{\boldsymbol{V}_{h}\times
Q_{h}}\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\boldsymbol{V}_{h}\times
Q_{h}}}\geq C>0.$

###### Proof.

For given $(\tilde{\boldsymbol{v}},\tilde{q})\in\boldsymbol{V}_{h}\times
Q_{h}$ it suffices to show that there exists
$(\boldsymbol{v},q)\in\boldsymbol{V}_{h}\times Q_{h}$ such that

$\displaystyle\|(\boldsymbol{v},q)\|_{\boldsymbol{V}_{h}\times Q_{h}}$
$\displaystyle\leq
C_{3}\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\boldsymbol{V}_{h}\times Q_{h}},$
$\displaystyle B((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle\geq
C_{4}\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\boldsymbol{V}_{h}\times
Q_{h}}^{2}$

for some constants $C_{3},C_{4}>0$ independent of $h$.

Suppose that
$(0,0)\not=(\tilde{\boldsymbol{v}},\tilde{q})\in\boldsymbol{V}_{h}\times
Q_{h}$ is given. Let $\boldsymbol{w}$ be an element in $\boldsymbol{V}_{h}$
such that

(21)
$\displaystyle\frac{\left(\operatorname{div}\boldsymbol{w},\tilde{q}\right)_{\Omega}}{\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}}\geq
C_{1}\|\tilde{q}\|_{0,\lambda^{-1}}-C_{2}S(\tilde{q},\tilde{q})^{\frac{1}{2}}$

by Lemma 3.6. By rescaling we can assume that
$\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}=\|q\|_{0,\lambda^{-1}}$ without loss
of generality, so we can obtain

(22)
$\displaystyle{\left(\operatorname{div}\boldsymbol{w},\tilde{q}\right)_{\Omega}}\geq
C_{1}^{\prime}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}-C_{2}^{\prime}S(\tilde{q},\tilde{q})$

from (21) with Young’s inequality.

Take $(\boldsymbol{v},q)$ as
$\boldsymbol{v}=\tilde{\boldsymbol{v}}+\delta_{1}\boldsymbol{w}$ and
$q=-\tilde{q}$ with $\delta_{1}>0$ which will be determined later. Then,

(23) $B((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))\\\
=a_{\boldsymbol{u}}(\tilde{\boldsymbol{v}},\tilde{\boldsymbol{v}})+\delta_{1}a_{\boldsymbol{u}}(\boldsymbol{w},\tilde{\boldsymbol{v}})+\delta_{1}\left(\alpha\tilde{q},\operatorname{div}\boldsymbol{w}\right)_{\Omega}+\left(s_{0}\tilde{q},\tilde{q}\right)_{\Omega}+S(\tilde{q},\tilde{q}).$

By (21) with $\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}=1$ and the assumption
that $\alpha$ has a uniformly positive lower bound, we have

$\displaystyle\left(\alpha\tilde{q},\operatorname{div}\boldsymbol{w}\right)_{\Omega}\geq
C_{1}^{\prime\prime}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}-C_{2}^{\prime\prime}S(\tilde{q},\tilde{q})$

with $C_{1}^{\prime\prime},C_{2}^{\prime\prime}$ depending on
$C_{1}^{\prime}$, $C_{2}^{\prime}$ and the lower bound of $\alpha$. Note that

$\displaystyle\delta_{1}|a_{\boldsymbol{u}}(\boldsymbol{w},\tilde{\boldsymbol{v}})|\leq\frac{1}{2}\|\tilde{\boldsymbol{v}}\|_{\boldsymbol{V}_{h}}^{2}+\frac{\delta_{1}^{2}}{2}\|\boldsymbol{w}\|_{\boldsymbol{V}_{h}}^{2}=\frac{1}{2}\|\tilde{\boldsymbol{v}}\|_{\boldsymbol{V}_{h}}^{2}+\frac{\delta_{1}^{2}}{2}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}.$

If we use these two inequalities to (23), then we get

$\displaystyle B((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle\geq
a_{\boldsymbol{u}}(\tilde{\boldsymbol{v}},\tilde{\boldsymbol{v}})-\delta_{1}|a_{\boldsymbol{u}}(\boldsymbol{w},\tilde{\boldsymbol{v}})|+\delta_{1}\left(C_{1}^{\prime\prime}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}-C_{2}^{\prime\prime}S(\tilde{q},\tilde{q})\right)$
$\displaystyle\quad+\left(s_{0}\tilde{q},\tilde{q}\right)_{\Omega}+S(\tilde{q},\tilde{q})$
$\displaystyle\geq\frac{1}{2}a_{\boldsymbol{u}}(\tilde{\boldsymbol{v}},\tilde{\boldsymbol{v}})-\frac{\delta_{1}^{2}}{2}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}+\delta_{1}\left(C_{1}^{\prime\prime}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}-C_{2}^{\prime\prime}S(\tilde{q},\tilde{q})\right)$
$\displaystyle\quad+\left(s_{0}\tilde{q},\tilde{q}\right)_{\Omega}+S(\tilde{q},\tilde{q}).$

If we choose a sufficiently small $\delta_{1}$ satisfying
$\delta_{1}C_{1}^{\prime\prime}\geq\delta_{1}^{2}$ and
$\delta_{1}C_{2}^{\prime\prime}\leq 1/2$, then there exists $C_{4}>0$
depending only on $C_{1}^{\prime\prime}$, $C_{2}^{\prime\prime}$ such that

$\displaystyle B((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle\geq
C_{4}\|\delta_{1}(\tilde{\boldsymbol{v}},\tilde{q})\|_{\boldsymbol{V}_{h}\times
Q_{h}}^{2}.$

By the definition of $(\boldsymbol{v},q)$ it is not difficult to see that
$\|(\boldsymbol{v},q)\|_{\boldsymbol{V}_{h}\times Q_{h}}\leq
C_{3}\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\boldsymbol{V}_{h}\times Q_{h}}$
with a constant $C_{3}$ depending only on $\delta_{1}$. Therefore, the
assertion is proved. ∎

Before we prove the a priori error analysis, we show that the proposed method
needs considerably smaller number of degrees of freedom (DOFs). For the
discretization of $\boldsymbol{u}$ the MTW element has minimal inter-element
continuity for discrete Korn inequality ([24]), so we only need to compare the
degrees of freedom for $p$ when $p$ is discretized by the dual form of mixed
method with the lowest order Raviart-Thomas element.

######  Remark 3.8.

For easy computation of DOFs we consider a triangulation of the unit cube
domain $\Omega$ such that $\Omega$ is subdivided into $N\times N\times N$
subcubes and each subcube is decomposed to 6 tetrahedra. Since each subcube
has 6 internal faces, the numbers of vertices, faces, and tetrahedra are
$(N+1)^{3}$, $6N^{2}(N+1)+6N^{3}$, $6N^{3}$. The total number of DOFs of mixed
method discretization of $p$ with the lowest order Raviart-Thomas element, is
$6N^{2}(N+1)+6N^{3}+6N^{3}\sim 18N^{3}$. After hybridization (see, e.g., [1]),
only one face DOF is needed, so this can be reduced to $\sim 12N^{3}$. The
lowest order enriched Galerkin method needs only one DOF for vertices and one
DOF for tetrahedra, so the total number of DOFs is $(N+1)^{3}+6N^{3}\sim
7N^{3}$ which still is considerably smaller than the one of hybridized mixed
method. For hexahedral finite elements the numbers of vertices, faces, cubes
are $(N+1)^{3}$, $3N^{2}(N+1)$, $N^{3}$. Therefore, the hybridized method with
the lowest order Raviart-Thomas element needs approximately $3N^{3}$ DOFs
whereas the IOS-EG method needs approximately $2N^{3}$ DOFs.

###  3.3. Semidiscrete solutions and well-posedness

According to the variational formulation (3) we define the semidiscrete
problem as follows: For compatible numerical initial data
$(\boldsymbol{u}_{h}(0),p_{h}(0))\in\boldsymbol{V}_{h}\times Q_{h}$ satisfying

(24) $\displaystyle
a_{\boldsymbol{u}}\left(\boldsymbol{u}_{h}(0),\boldsymbol{v}\right)-\left(\alpha
p_{h}(0),\operatorname{div}\boldsymbol{v}\right)_{\Omega}=\left(\boldsymbol{f}(0),\boldsymbol{v}\right)_{\Omega}\qquad\forall\boldsymbol{v}\in\boldsymbol{V}_{h}$

we find a semidiscrete solution
$(\boldsymbol{u}_{h},p_{h}):C^{1}(0,\infty;\boldsymbol{V}_{h}\times Q_{h})$
such that

(25a) $\displaystyle
a_{\boldsymbol{u}}\left(\boldsymbol{u}_{h}(t),\boldsymbol{v}\right)-\left(\alpha
p_{h}(t),\operatorname{div}\boldsymbol{v}\right)_{\Omega}$
$\displaystyle=\left(\boldsymbol{f}(t),\boldsymbol{v}\right)_{\Omega},$ (25b)
$\displaystyle-\left(\alpha\operatorname{div}\dot{\boldsymbol{u}}_{h}(t),q\right)_{\Omega}-S\left(\dot{p}_{h}(t),q\right)-\left(s_{0}\dot{p}_{h}(t),q\right)_{\Omega}-a_{p}\left(p_{h}(t),q\right)$
$\displaystyle=-\left(g(t),q\right)_{\Omega}$

for all $(\boldsymbol{v},q)\in\boldsymbol{V}_{h}\times Q_{h}$ and for all
$t\geq 0$.

######  Theorem 3.9.

Suppose that $(\boldsymbol{u}_{h}(0),p_{h}(0))\in\boldsymbol{V}_{h}\times
Q_{h}$ satisfies (24). Then, (25) has a unique solution.

###### Proof.

For well-posedness of (25) we consider a matrix realization of (25). We
identify $\boldsymbol{V}_{h}$ and $Q_{h}$ to
${\mathbb{R}}^{N_{\boldsymbol{u}}}$ and ${\mathbb{R}}^{N_{p}}$ with certain
global bases of $\boldsymbol{V}_{h}$ and $Q_{h}$, and denote the vector
realizations of $\boldsymbol{u}_{h}$, $p_{h}$ and the $L^{2}$ projections of
$\boldsymbol{f},g$ to $\boldsymbol{V}_{h},Q_{h}$ by
$\vec{\boldsymbol{u}}_{h}$, $\vec{p}_{h}$, $\vec{\boldsymbol{f}}$, $\vec{g}$.
Let

$\displaystyle\mathbb{A}_{\boldsymbol{u}},\quad\mathbb{B},\quad\mathbb{S},\quad\mathbb{A}_{p}$

be the matrices generated by the bilinear forms

$\displaystyle
a_{\boldsymbol{u}}(\boldsymbol{v},\tilde{\boldsymbol{v}}),\quad-\left(\alpha\operatorname{div}\boldsymbol{v},q\right)_{\Omega},\quad
S(q,\tilde{q})+\left(s_{0}q,\tilde{q}\right)_{\Omega},\quad
a_{p}(q,\tilde{q})$

for $\boldsymbol{v},\tilde{\boldsymbol{v}}\in\boldsymbol{V}_{h}$,
$q,\tilde{q}\in Q_{h}$ with the identifications
$\boldsymbol{V}_{h}\leftrightarrow{\mathbb{R}}^{N_{\boldsymbol{u}}}$ and
$Q_{h}\leftrightarrow{\mathbb{R}}^{N_{p}}$. Then, (25) in matrix form is

$\displaystyle\mathbb{A}_{\boldsymbol{u}}{\vec{\boldsymbol{u}}}_{h}(t)+\mathbb{B}^{T}\vec{p}_{h}(t)$
$\displaystyle=\vec{\boldsymbol{f}(t)},$
$\displaystyle\mathbb{B}\dot{\vec{\boldsymbol{u}}}_{h}(t)-\mathbb{S}\dot{\vec{p}}_{h}(t)-\mathbb{A}_{p}\vec{p}_{h}(t)$
$\displaystyle=-\vec{g}(t).$

Since $\mathbb{A}_{\boldsymbol{u}}$ is invertible
$\vec{\boldsymbol{u}}_{h}=-\mathbb{A}_{\boldsymbol{u}}^{-1}\mathbb{B}^{T}\vec{p}_{h}+\mathbb{A}_{\boldsymbol{u}}^{-1}\vec{\boldsymbol{f}}$
from the the first equation. If we use this to the second equation, then we
get

$\displaystyle-\left(\mathbb{B}\mathbb{A}_{\boldsymbol{u}}^{-1}\mathbb{B}+\mathbb{S}\right)\dot{\vec{p}}_{h}(t)-\mathbb{A}_{p}\vec{p}_{h}(t)=-\mathbb{B}\mathbb{A}_{\boldsymbol{u}}^{-1}\dot{\vec{\boldsymbol{f}}}(t)-\vec{g}(t).$

This is a standard form of system of linear ordinary differential equations,
and this equation has a unique solution if
$\mathbb{B}\mathbb{A}_{\boldsymbol{u}}^{-1}\mathbb{B}+\mathbb{S}$ is
invertible. By Lemma 3.7 the matrix

$\displaystyle\begin{pmatrix}\mathbb{A}_{\boldsymbol{u}}&\mathbb{B}^{T}\\\
\mathbb{B}&-\mathbb{S}\end{pmatrix}$

is invertible, so its Schur complement
$-(\mathbb{B}\mathbb{A}_{\boldsymbol{u}}^{-1}\mathbb{B}+\mathbb{S})$ is
invertible as well. Therefore, (25) has a unique solution. ∎

##  4\. Error analysis of semidiscrete solutions

For error analysis we define
$\Pi_{h}^{\boldsymbol{V}}:\boldsymbol{V}\to\boldsymbol{V}_{h}$ and
$\Pi_{h}^{Q}:Q\to Q_{h}$ by

$\displaystyle
a_{\boldsymbol{u}}\left(\Pi_{h}^{\boldsymbol{V}}\boldsymbol{v},\tilde{\boldsymbol{v}}\right)$
$\displaystyle=a_{\boldsymbol{u}}\left(\boldsymbol{v},\tilde{\boldsymbol{v}}\right)$
$\displaystyle\forall\tilde{\boldsymbol{v}}\in\boldsymbol{V}_{h},$
$\displaystyle a_{p}\left(\Pi_{h}^{Q}q,\tilde{q}\right)$
$\displaystyle=a_{p}\left(q,\tilde{q}\right)$
$\displaystyle\forall\tilde{q}\in Q_{h}.$

It is known that

(26a)
$\displaystyle\|\boldsymbol{v}-\Pi_{h}^{\boldsymbol{V}}\boldsymbol{v}\|_{a_{\boldsymbol{u}}}\leq
Ch^{m}\|\boldsymbol{v}\|_{m+1},$ (26b)
$\displaystyle\|q-\Pi_{h}^{Q}q\|_{0}\leq C\|q-\Pi_{h}^{Q}q\|_{h}\leq
Ch^{m}\|q\|_{m+1},$ (26c)
$\displaystyle\sum_{e\in\mathcal{E}_{h}^{0}}h_{e}^{-1-\beta}\int_{e}\left[\\!\left[{\Pi_{h}^{Q}q}\right]\\!\right]^{2}\,ds\leq
Ch^{2m}\|q\|_{m+1}^{2}$

in [23, 37] and [18]. Moreover, if $\Omega$ satisfies the full elliptic
regularity assumption, then

$\displaystyle\|q-\Pi_{h}^{Q}q\|_{0}\leq Ch^{m+1}\|q\|_{m+1}.$

######  Theorem 4.1.

Suppose that $(\boldsymbol{u},p)$ and $(\boldsymbol{u}_{h},p_{h})$ are the
solutions of (3) and (25) with compatible numerical initial data
$(\boldsymbol{u}_{h}(0),p_{h}(0))\in\boldsymbol{V}_{h}\times Q_{h}$ satisfying
(24) and

(27a) $\displaystyle
2\mu\|\boldsymbol{u}(0)-\boldsymbol{u}_{h}(0)\|_{0}^{2}+\lambda\|\operatorname{div}(\boldsymbol{u}(0)-\boldsymbol{u}_{h}(0))\|_{0}^{2}$
$\displaystyle\leq Ch^{m}\|\boldsymbol{u}(0)\|_{m+1},$ (27b)
$\displaystyle\|p(0)-p_{h}(0)\|_{0},\|p(0)-p_{h}(0)\|_{1,h}$
$\displaystyle\leq Ch^{m}\|p(0)\|_{m+1}.$

Then,

$\displaystyle\|\boldsymbol{u}(t)-\boldsymbol{u}_{h}(t)\|_{\boldsymbol{V}_{h}}+\|p(t)-p_{h}(t)\|_{0,s_{0}}$
(28) $\displaystyle\leq
Ch^{m}\left(\|\boldsymbol{u}(0)\|_{m+1}+\|p(0)\|_{m+1}\right)+C\lambda^{-\frac{1}{2}}h^{m}\|p\|_{W^{1,1}(0,t;H^{m+1})}$
$\displaystyle\quad+Ch^{m}\|\boldsymbol{u}\|_{W^{1,1}(0,t;H^{m+1})}$

and

(29) $\|p(t)-p_{h}(t)\|_{1,h}\\\ \leq
Ch^{m}\left(\|p(0)\|_{m+1}+\lambda^{-\frac{1}{2}}\|p\|_{W^{1,2}(0,t;H^{m+1})}+\|\boldsymbol{u}\|_{W^{1,2}(0,t;H^{m+1})}\right).$

###### Proof.

Let $e_{\boldsymbol{u}}(t):=\boldsymbol{u}(t)-\boldsymbol{u}_{h}(t)$ and
$e_{p}(t):=p(t)-p_{h}(t)$. Now, subtracting (25) from (3), we can derive our
error equations as

(30a) $\displaystyle
a_{\boldsymbol{u}}\left(e_{\boldsymbol{u}}(t),\boldsymbol{v}\right)-\left(\alpha
e_{p}(t),\operatorname{div}\boldsymbol{v}\right)_{\Omega}$ $\displaystyle=0,$
(30b)
$\displaystyle-\left(\alpha\operatorname{div}\dot{e}_{\boldsymbol{u}}(t),q\right)_{\Omega}-S\left(\dot{e}_{p}(t),q\right)-\left(s_{0}\dot{e}_{p}(t),q\right)_{\Omega}+a_{p}\left(e_{p}(t),q\right)$
$\displaystyle=0$

for all $\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$, $q\in Q_{h}$. Defining

$\displaystyle e_{\boldsymbol{u}}^{h}(t)$
$\displaystyle:=\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u}(t)-\boldsymbol{u}_{h}(t),$
$\displaystyle e_{\boldsymbol{u}}^{I}(t)$
$\displaystyle:=\boldsymbol{u}(t)-\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u}(t),$
$\displaystyle e_{p}^{h}(t)$ $\displaystyle:=\Pi_{h}^{Q}p(t)-p_{h}(t),$
$\displaystyle e_{p}^{I}(t)$ $\displaystyle:=p(t)-\Pi_{h}^{Q}p(t),$

one can see that
$e_{\boldsymbol{u}}(t)=e_{\boldsymbol{u}}^{I}(t)+e_{\boldsymbol{u}}^{h}(t)$,
and
$a_{\boldsymbol{u}}\left(e_{\boldsymbol{u}}^{I}(t),\boldsymbol{v}\right)=0$
for all $\boldsymbol{v}\in\boldsymbol{V}_{h}$ by the definition of
$\Pi_{h}^{\boldsymbol{V}}$. Thus, we can rewrite (30) as

(31a) $\displaystyle
a_{\boldsymbol{u}}\left(e_{\boldsymbol{u}}^{h}(t),\boldsymbol{v}\right)-\left(\alpha\left(e_{p}^{h}(t)+e_{p}^{I}(t)\right),\operatorname{div}\boldsymbol{v}\right)_{\Omega}$
$\displaystyle=0,$
$\displaystyle-\left(\alpha\operatorname{div}(\dot{e}_{\boldsymbol{u}}^{h}(t)+\dot{e}_{\boldsymbol{u}}^{I}(t)),q\right)_{\Omega}-S(\dot{e}_{p}^{h}(t)+\dot{e}_{p}^{I}(t),q)$
(31b)
$\displaystyle-\left(s_{0}(\dot{e}_{p}^{h}(t)+\dot{e}_{p}^{I}(t)),q\right)_{\Omega}-a_{p}\left(e_{p}^{h}(t),q\right)$
$\displaystyle=0$

for $(\boldsymbol{v},q)\in\boldsymbol{V}_{h}\times Q_{h}$. For our error
analysis take $\boldsymbol{v}=\dot{e}_{\boldsymbol{u}}^{h}(t)$ and
$q=-e_{p}^{h}(t)$ and add (31a) and (31b) to obtain the following:

(32)
$a_{\boldsymbol{u}}\left(e_{\boldsymbol{u}}^{h}(t),\dot{e}_{\boldsymbol{u}}^{h}(t)\right)-\left(\alpha
e_{p}^{I}(t),\operatorname{div}\dot{e}_{\boldsymbol{u}}^{h}(t)\right)_{\Omega}\\\
+\left(\alpha\operatorname{div}\dot{e}_{\boldsymbol{u}}^{I}(t),e_{p}^{h}(t)\right)_{\Omega}+S\left(\dot{e}_{p}^{h}(t)+\dot{e}_{p}^{I}(t),e_{p}^{h}(t)\right)\\\
+\left(s_{0}(\dot{e}_{p}^{h}(t)+\dot{e}_{p}^{I}(t)),e_{p}^{h}(t)\right)_{\Omega}+a_{p}\left(e_{p}^{h}(t),e_{p}^{h}(t)\right)=0.$

We can rewrite this as

(33)
$\frac{1}{2}\frac{d}{dt}a_{\boldsymbol{u}}\left(e_{\boldsymbol{u}}^{h}(t),e_{\boldsymbol{u}}^{h}(t)\right)+\frac{1}{2}\frac{d}{dt}\left(s_{0}e_{p}^{h}(t),e_{p}^{h}(t)\right)_{\Omega}\\\
+a_{p}\left(e_{p}^{h}(t),e_{p}^{h}(t)\right)+\frac{1}{2}\frac{d}{dt}S\left(e_{p}^{h}(t),e_{p}^{h}(t)\right)\\\
=\left(\alpha
e_{p}^{I}(t),\operatorname{div}\dot{e}_{\boldsymbol{u}}^{h}(t)\right)_{\Omega}-\left(\alpha\operatorname{div}\dot{e}_{\boldsymbol{u}}^{I}(t),e_{p}^{h}(t)\right)_{\Omega}-S\left(\dot{e}_{p}^{I}(t),e_{p}^{h}(t)\right).$

Let us define $X(t)\geq 0$ by

$\displaystyle
X(t)^{2}:=a_{\boldsymbol{u}}\left(e_{\boldsymbol{u}}^{h}(t),e_{\boldsymbol{u}}^{h}(t)\right)+(s_{0}e_{p}^{h}(t),e_{p}^{h}(t))_{\Omega}+S\left(e_{p}^{h}(t),e_{p}^{h}(t)\right).$

For the proof of (28) we first prove

(34) $X(t)\leq X(0)+C\lambda^{-\frac{1}{2}}\max_{0\leq s\leq
t}\|e_{p}^{I}(s)\|_{0}\\\
+C\left(\lambda^{-\frac{1}{2}}\int_{0}^{t}\|\dot{e}_{p}^{I}(s)\|_{0}\,ds+\int_{0}^{t}\|\operatorname{div}\dot{e}_{\boldsymbol{u}}^{I}(s)\|_{0}\,ds\right)$

for any $t>0$.

To prove it, we integrate (33) from 0 to $t$ and get

(35)
$\frac{1}{2}X(t)^{2}-\frac{1}{2}X(0)^{2}+\int_{0}^{t}a_{p}\left(e_{p}^{h}(s),e_{p}^{h}(s)\right)\,ds\\\
=\int_{0}^{t}\left(\alpha
e_{p}^{I}(s),\operatorname{div}\dot{e}_{\boldsymbol{u}}^{h}(s)\right)_{\Omega}\,ds-\int_{0}^{t}\left(\alpha\operatorname{div}\dot{e}_{p}^{I}(s),e_{p}^{h}(s)\right)_{\Omega}\,ds-S\left(\dot{e}_{p}^{I}(t),e_{p}^{h}(t)\right).$

By the integration by parts in $t$,

$\displaystyle\int_{0}^{t}\left(\alpha
e_{p}^{I}(s),\operatorname{div}\dot{e}_{\boldsymbol{u}}^{h}(s)\right)_{\Omega}\,ds$
$\displaystyle=\left(\alpha
e_{p}^{I}(t),\operatorname{div}e_{\boldsymbol{u}}^{h}(t)\right)_{\Omega}-\left(\alpha
e_{p}^{I}(0),\operatorname{div}e_{\boldsymbol{u}}^{h}(0)\right)_{\Omega}$
$\displaystyle\quad-\int_{0}^{t}\left(\alpha\dot{e}_{p}^{I}(s),\operatorname{div}e_{\boldsymbol{u}}^{h}(s)\right)_{\Omega}\,ds.$

Plugging this into (35) yields

$\displaystyle\frac{1}{2}X(t)^{2}+\int_{0}^{T}a_{p}\left(e_{p}^{h}(s),e_{p}^{h}(s)\right)\,ds$
(36) $\displaystyle=\frac{1}{2}X(0)^{2}+\left(\alpha
e_{p}^{I}(t),\operatorname{div}e_{\boldsymbol{u}}^{h}(t)\right)_{\Omega}-\left(\alpha
e_{p}^{I}(0),\operatorname{div}e_{\boldsymbol{u}}^{h}(0)\right)_{\Omega}$
$\displaystyle\quad-\int_{0}^{T}\left(\alpha\dot{e}_{p}^{I}(s),\operatorname{div}e_{\boldsymbol{u}}^{h}(s)\right)_{\Omega}\,ds-\int_{0}^{T}\left(\alpha\operatorname{div}\dot{e}_{p}^{I}(s),e_{p}^{h}(s)\right)_{\Omega}\,ds$
$\displaystyle\quad-S\left(\dot{e}_{p}^{I}(t),e_{p}^{h}(t)\right).$

We now prove (34) assuming that

(37) $\displaystyle X(t)=\max_{0\leq s\leq t}X(s).$

With this assumption we can bound some terms in (36) as follows:

(38) $\displaystyle\left|\left(\alpha
e_{p}^{I}(t),\operatorname{div}e_{\boldsymbol{u}}^{h}(t)\right)_{\Omega}\right|$
$\displaystyle\leq C\lambda^{-\frac{1}{2}}\|e_{p}^{I}(t)\|_{0}X(t),$ (39)
$\displaystyle\left|\left(\alpha
e_{p}^{I}(0),\operatorname{div}e_{\boldsymbol{u}}^{h}(0)\right)_{\Omega}\right|$
$\displaystyle\leq C\lambda^{-\frac{1}{2}}\|e_{p}^{I}(0)\|_{0}X(0),$ (40)
$\displaystyle\left|\int_{0}^{t}\left(\alpha\dot{e}_{p}^{I}(s),\operatorname{div}e_{\boldsymbol{u}}^{h}(s)\right)_{\Omega}\,ds\right|$
$\displaystyle\leq
C\lambda^{-\frac{1}{2}}\int_{0}^{t}\|\dot{e}_{p}^{I}(s)\|_{0}\,dsX(t),$ (41)
$\displaystyle\left|\int_{0}^{t}\left(\alpha\operatorname{div}\dot{e}_{p}^{I}(s),e_{p}^{h}(s)\right)_{\Omega}\,ds\right|$
$\displaystyle\leq
C\int_{0}^{T}\|\operatorname{div}\dot{e}_{\boldsymbol{u}}^{I}(s)\|_{0}\,dsX(t).$

Plugging these into (36), we derive

$X(t)^{2}+2\int_{0}^{t}a_{p}\left(e_{p}^{h}(s),e_{p}^{h}(s)\right)\,ds\leq
X(0)^{2}+C\lambda^{-\frac{1}{2}}\|e_{p}^{I}(0)\|_{0}X(0)\\\
+C\left(\lambda^{-\frac{1}{2}}\|e_{p}^{I}(t)\|_{0}+\lambda^{-\frac{1}{2}}\int_{0}^{t}\|\dot{e}_{p}^{I}(s)\|_{0}\,ds+\int_{0}^{T}\|\operatorname{div}\dot{e}_{\boldsymbol{u}}^{I}(s)\|_{0}\,ds\right)X(t).$

If we divide by $X(t)$, then we can obtain (34). For general $t$ such that
(37) is not true, note that there exists $0\leq t_{M}<t$ such that

$\displaystyle X(t_{M})=\max_{0\leq s\leq t_{M}}X(s),\qquad X(t)<X(t_{M}).$

By the same argument as before for $t_{M}$ we can obtain

(42) $X(t_{M})\leq X(0)+C\lambda^{-\frac{1}{2}}\|e_{p}^{I}(0)\|_{0}\\\
+C\left(\lambda^{-\frac{1}{2}}\|e_{p}^{I}(t_{M})\|_{0}+\lambda^{-\frac{1}{2}}\int_{0}^{t_{M}}\|\dot{e}_{p}^{I}(s)\|_{0}\,ds+\int_{0}^{t_{M}}\|\operatorname{div}\dot{e}_{\boldsymbol{u}}^{I}(s)\|_{0}\,ds\right).$

We can now derive (34) by

$\displaystyle X(t)<X(t_{M}),\qquad\max_{0\leq s\leq
t_{M}}\|e_{p}^{I}(s)\|_{0}\leq\max_{0\leq s\leq t}\|e_{p}^{I}(t)\|_{0},$

and

$\lambda^{-1}\int_{0}^{t_{M}}\|\dot{e}_{p}^{I}(s)\|_{0}+\int_{0}^{t_{M}}\|\operatorname{div}\dot{e}_{\boldsymbol{u}}^{I}(s)\|_{0}\,ds\\\
\leq\lambda^{-1}\int_{0}^{t}\|\dot{e}_{p}^{I}(s)\|_{0}+\int_{0}^{t}\|\operatorname{div}\dot{e}_{\boldsymbol{u}}^{I}(s)\|_{0}\,ds.$

Then, (28) follows by the triangle inequality, (34), (27), and the
interpolation error bounds (26). ∎

##  5\. Preconditioners for fast solver algorithms

In this section we present an abstract form of preconditioners for our
discretization of poroelasticity equations by adopting the operator
preconditioning framework in [25].

After the Crank-Nicolson time discretization we need to solve the static
equation

$\displaystyle\mathcal{B}((\boldsymbol{u}^{n+1},p^{n+1}),({\boldsymbol{v}},{q}))=\mathcal{F}^{n+1}({\boldsymbol{v}},{q}),\qquad({\boldsymbol{v}},{q})\in\boldsymbol{V}_{h}\times
Q_{h}$

where $\mathcal{F}^{n+1}$ is determined by the previous time step solution
$(\boldsymbol{u}^{n},p^{n})$ and the right-hand side terms, and
$\mathcal{B}:(\boldsymbol{V}_{h}\times Q_{h})\times(\boldsymbol{V}_{h}\times
Q_{h})\rightarrow{\mathbb{R}}$ is defined by

$\displaystyle\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle=a_{\boldsymbol{u}}(\boldsymbol{v},\tilde{\boldsymbol{v}})-\left(\alpha
q,\operatorname{div}\tilde{\boldsymbol{v}}\right)_{\Omega}-\left(\alpha\tilde{q},\operatorname{div}\boldsymbol{v}\right)_{\Omega}-\left(s_{0}q,\tilde{q}\right)_{\Omega}-S\left(q,\tilde{q}\right)$
$\displaystyle\quad-\frac{\Delta t}{2}a_{p}(q,\tilde{q}).$

Let $\tilde{Q}_{h}$ be the space $Q_{h}$ with the norm

$\displaystyle\|q\|_{\tilde{Q}_{h}}^{2}=\|q^{c}\|_{0,\lambda^{-1}}^{2}+\|q^{c}\|_{0,s_{0}}^{2}+\|q^{0}\|_{0,\lambda^{-1}}^{2}+\|q^{0}\|_{0,s_{0}}^{2}+S(q^{0},\tilde{q}^{0})+\frac{\Delta
t}{2}\|q\|_{h}^{2}.$

Where $q^{0},\tilde{q}^{0}$ are mean-value zero on $\Omega$. For simplicity
let us define $\|(\boldsymbol{v},q)\|_{\mathcal{X}}$ by

$\displaystyle\|(\boldsymbol{v},q)\|_{\mathcal{X}}^{2}=\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}^{2}+\|q\|_{\tilde{Q}_{h}}^{2}.$

We will prove that the linear operator
$\mathcal{L}_{\mathcal{B}}:\boldsymbol{V}_{h}\times\tilde{Q}_{h}\rightarrow(\boldsymbol{V}_{h}\times\tilde{Q}_{h})^{*}$
defined by

(43)
$\displaystyle\left<\mathcal{L}_{\mathcal{B}}(\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q})\right>_{\left((\boldsymbol{V}_{h}\times\tilde{Q}_{h})^{*},\boldsymbol{V}_{h}\times\tilde{Q}_{h}\right)}:=\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q})),$

is a bounded isomorphism from $\boldsymbol{V}_{h}\times\tilde{Q}_{h}$ to the
dual space $(\boldsymbol{V}_{h}\times\tilde{Q}_{h})^{*}$. Then, a parameter-
robust preconditioner can be constructed as the inverse of the block diagonal
operator

(44)
$\displaystyle\mathcal{P}=\begin{pmatrix}\mathcal{P}_{\boldsymbol{V}_{h}}&0\\\
0&\mathcal{P}_{\tilde{Q}_{h}}\end{pmatrix}$

such that
$\mathcal{P}_{\boldsymbol{V}_{h}}:\boldsymbol{V}_{h}\rightarrow\boldsymbol{V}_{h}^{*}$
and $\mathcal{P}_{\tilde{Q}_{h}}:\tilde{Q}_{h}\rightarrow\tilde{Q}_{h}^{*}$
are the bounded linear operators defined by

$\displaystyle\left<\mathcal{P}_{\boldsymbol{V}_{h}}\boldsymbol{v},\tilde{\boldsymbol{v}}\right>_{(\boldsymbol{V}_{h}^{*},\boldsymbol{V}_{h})}$
$\displaystyle:=a_{\boldsymbol{u}}(\boldsymbol{v},\tilde{\boldsymbol{v}}),$
$\displaystyle\left<\mathcal{P}_{\tilde{Q}_{h}}q,\tilde{q}\right>_{(\tilde{Q}_{h}^{*},\tilde{Q}_{h})}$
$\displaystyle:=\left(s_{0}q^{c},\tilde{q}^{c}\right)_{\Omega}+\left(\lambda^{-1}q^{c},\tilde{q}^{c}\right)_{\Omega}$
$\displaystyle\quad+\frac{\Delta
t}{2}\left(\int_{\Omega}\underline{\boldsymbol{\kappa}}\nabla q^{c}\cdot\nabla
q^{c}\,dx+\sum_{e\in\mathcal{E}_{h}^{\partial}}\gamma
h_{e}^{-1}\left<\left[\\!\left[{q^{c}}\right]\\!\right],\left[\\!\left[{\tilde{q}^{c}}\right]\\!\right]\right>_{e}\right)$
$\displaystyle\quad+\left(s_{0}q^{0},\tilde{q}^{0}\right)_{\Omega}+\left(\lambda^{-1}q^{0},\tilde{q}^{0}\right)_{\Omega}+S(q^{0},\tilde{q}^{0})$
$\displaystyle\quad+\frac{\Delta
t}{2}\left(\sum_{e\in\mathcal{E}_{h}^{0}}\gamma
h_{e}^{-1-\beta}\left<\left[\\!\left[{q^{0}}\right]\\!\right],\left[\\!\left[{\tilde{q}^{0}}\right]\\!\right]\right>_{e}+\sum_{e\in\mathcal{E}_{h}^{\partial}}\gamma
h_{e}^{-1}\left<\left[\\!\left[{q^{0}}\right]\\!\right],\left[\\!\left[{\tilde{q}^{0}}\right]\\!\right]\right>_{e}\right).$

If
$\mathcal{L}_{\mathcal{B}}:\boldsymbol{V}_{h}\times\tilde{Q}_{h}\rightarrow(\boldsymbol{V}_{h}\times\tilde{Q}_{h})^{*}$
is an isomorphism, then a preconditioner of the form (44) is expected to give
parameter-robust preconditioners based on norm equivalence (cf. [25]).

By the definition (43),
$\mathcal{L}_{\mathcal{B}}:\boldsymbol{V}_{h}\times\tilde{Q}_{h}\rightarrow(\boldsymbol{V}_{h}\times\tilde{Q}_{h})^{*}$
is a bounded isomorphism independent of $h,\Delta t,\lambda$ if $\mathcal{B}$
is a bounded bilinear form on $\boldsymbol{V}_{h}\times\tilde{Q}_{h}$, and the
following inf-sup condition holds: There exists $C>0$ independent of $h$,
$\lambda$, $\Delta t$ such that

(45)
$\displaystyle\inf_{(\tilde{\boldsymbol{v}},\tilde{q})\in\boldsymbol{V}_{h}\times
Q_{h}}\sup_{(\boldsymbol{v},q)\in\boldsymbol{V}_{h}\times
Q_{h}}\frac{\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))}{\|(\boldsymbol{v},q)\|_{\mathcal{X}}\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\mathcal{X}}}\geq
C>0.$

The main result in this section is proving the inf-sup condition (45). As a
consequence, $\mathcal{L}_{\mathcal{B}}$ is an isomorphism, so preconditioners
of the form (44) give robust preconditioners.

Before we prove this we state a necessary result which is proved in [18].

######  Lemma 5.1.

For given $\tilde{q}=\tilde{q}^{c}+\tilde{q}^{0}$, one can find $r^{0}\in
Q_{h}^{0}$ satisfying

(46)
$\displaystyle\sum_{e\in\mathcal{E}_{h}^{0}}\left<h_{e}^{-1}\left[\\!\left[{r^{0}}\right]\\!\right],\left[\\!\left[{r^{0}}\right]\\!\right]\right>_{e}\leq
C_{r}\sum_{e\in\mathcal{E}_{h}^{0}}\left<h_{e}^{-1}\left[\\!\left[{\tilde{q}^{0}}\right]\\!\right],\left[\\!\left[{\tilde{q}^{0}}\right]\\!\right]\right>_{e},$
(47)
$\displaystyle\sum_{e\in\mathcal{E}_{h}^{\partial}}\left<h_{e}^{-1}\left[\\!\left[{r^{0}}\right]\\!\right],\left[\\!\left[{r^{0}}\right]\\!\right]\right>_{e}+\sum_{e\in\mathcal{E}_{h}^{0}}\left<h_{e}^{-1-\beta}\left[\\!\left[{r^{0}}\right]\\!\right],\left[\\!\left[{r^{0}}\right]\\!\right]\right>_{e}$
$\displaystyle\quad\leq
C_{r}\sum_{e\in\mathcal{E}_{h}^{0}}\left<h_{e}^{-1-\beta}\left[\\!\left[{\tilde{q}^{0}}\right]\\!\right],\left[\\!\left[{\tilde{q}^{0}}\right]\\!\right]\right>_{e},$
(48) $\displaystyle\|r^{0}\|_{0}\leq C_{r}\|\tilde{q}^{0}\|_{0},$

for some $C_{r}>0$ independent of $h$, and $q=\tilde{q}+\delta r^{0}$ with
sufficiently small (but independent of $h$) $\delta>0$ satisfies

(49) $\displaystyle a_{p}(q,\tilde{q})\geq C_{EG}\|q\|_{h}^{2}$

with $C_{EG}$ depending on the choice of $\delta$. In other words, the inf-sup
condition

$\displaystyle\inf_{\tilde{q}\in Q_{h}}\sup_{q\in
Q_{h}}\frac{a_{p}(q,\tilde{q})}{\|q\|_{h}\|\tilde{q}\|_{h}}\geq
C_{EG}^{\prime}>0$

holds with $C_{EG}^{\prime}$ independent of $h$.

We can now prove the main result.

######  Theorem 5.2.

The inf-sup condition (45) holds.

###### Proof of Theorem 5.2.

For given $(\tilde{\boldsymbol{v}},\tilde{q})\in\boldsymbol{V}_{h}\times
Q_{h}$ it suffices to show that there exists
$(\boldsymbol{v},q)\in\boldsymbol{V}_{h}\times Q_{h}$ such that

$\displaystyle\|(\boldsymbol{v},q)\|_{\mathcal{X}}$ $\displaystyle\leq
C_{5}\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\mathcal{X}},$
$\displaystyle\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle\geq
C_{6}\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\mathcal{X}}^{2}$

for some constants $C_{5},C_{6}>0$ independent of $h$, $\lambda$, $\Delta t$.

Suppose that
$(0,0)\not=(\tilde{\boldsymbol{v}},\tilde{q})\in\boldsymbol{V}_{h}\times
Q_{h}$ is given. By Lemma 5.1, there exists $r^{0}\in Q_{h}^{0}$ satisfying
(47), (48) such that $a_{p}(\tilde{q}+\delta_{0}r^{0},\tilde{q})\geq
C_{3}\|\tilde{q}\|_{h}^{2}$ and $\|\tilde{q}+\delta_{0}r^{0}\|_{h}\leq
C_{4}\|\tilde{q}\|_{h}$ with sufficiently small $\delta_{0}>0$ for
$C_{3},C_{4}>0$ independent of $h$, $\lambda$, $\Delta t$. By the same
argument in the proof of Lemma 3.7, there exists
$\boldsymbol{w}_{1}\in\boldsymbol{V}_{h}$ such that
$\|\boldsymbol{w}_{1}\|_{\boldsymbol{V}_{h}}=\|\tilde{q}\|_{0,\lambda^{-1}}$
and

(50)
$\displaystyle{\left(\operatorname{div}\boldsymbol{w}_{1},\tilde{q}\right)_{\Omega}}\geq
C_{1}^{\prime}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}-C_{2}^{\prime}S(\tilde{q},\tilde{q})$

for all $\tilde{q}\in Q_{h}$. Since
$\operatorname{div}\boldsymbol{V}_{h}=Q_{h}^{0}$, there exists
$\boldsymbol{w}_{2}\in\boldsymbol{V}_{h}$ such that
$\operatorname{div}\boldsymbol{w}_{2}=\lambda^{-1}\tilde{q}^{0}$ and
$\|\boldsymbol{w}_{2}\|_{\boldsymbol{V}_{h}}\leq
C\|\tilde{q}^{0}\|_{0,\lambda^{-1}}$.

Take $(\boldsymbol{v},q)$ as
$\boldsymbol{v}=\tilde{\boldsymbol{v}}+\delta_{1}\boldsymbol{w}_{1}$ and
$q=-\tilde{q}-\delta_{0}r^{0}$ with $\delta_{0},\delta_{1}>0$ which will be
determined later. Then,

(51)
$\displaystyle\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle=\|\tilde{\boldsymbol{v}}\|_{\boldsymbol{V}_{h}}^{2}+\delta_{1}a_{\boldsymbol{u}}(\boldsymbol{w}_{1},\tilde{\boldsymbol{v}})+\delta_{1}\left(\alpha\tilde{q},\operatorname{div}\boldsymbol{w}_{1}\right)_{\Omega}$
(52) $\displaystyle\quad+\delta_{0}\left(\alpha
r^{0},\operatorname{div}\tilde{\boldsymbol{v}}\right)_{\Omega}+\|\tilde{q}\|_{0,s_{0}}^{2}+\delta_{0}\left(s_{0}r^{0},\tilde{q}\right)_{\Omega}$
(53)
$\displaystyle\quad+S(\tilde{q},\tilde{q})+\delta_{0}S(r^{0},\tilde{q})+\frac{\Delta
t}{2}a_{p}\left(\tilde{q}+\delta_{0}r^{0},q\right).$

By (50), the parallelogram law,

(54) $\displaystyle
2\left(\|\tilde{q}\|_{0,\lambda^{-1}}^{2}+\|\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}\right)=\|\tilde{q}+\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}+\|\tilde{q}-\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}\geq\|\tilde{q}^{c}\|_{0,\lambda^{-1}}^{2}$

and Lemma 3.4, we have

$\displaystyle\delta_{1}\left(\alpha\tilde{q},\operatorname{div}\boldsymbol{w}_{1}\right)_{\Omega}$
$\displaystyle\geq\delta_{1}C_{1}^{\prime\prime}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}-\delta_{1}C_{2}^{\prime\prime}S(\tilde{q},\tilde{q})$
$\displaystyle\geq\frac{\delta_{1}C_{1}^{\prime\prime}}{2}\|\tilde{q}^{c}\|_{0,\lambda^{-1}}^{2}-\delta_{1}C_{1}^{\prime\prime}\|\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}-\delta_{1}C_{2}^{\prime\prime}S(\tilde{q},\tilde{q})$
(55)
$\displaystyle\geq\frac{\delta_{1}C_{1}^{\prime\prime}}{2}\|\tilde{q}^{c}\|_{0,\lambda^{-1}}^{2}-\delta_{1}\left(C_{P}C_{1}^{\prime\prime}+C_{2}^{\prime\prime}\right)S(\tilde{q},\tilde{q})$

with $C_{1}^{\prime\prime},C_{2}^{\prime\prime}$ which depend on the lower
bound of $\alpha$ but are independent of $h$, $\lambda$, $\Delta t$. By
Young’s inequality,

(56)
$\displaystyle\delta_{1}|a_{\boldsymbol{u}}(\boldsymbol{w}_{1},\tilde{\boldsymbol{v}})|$
$\displaystyle\leq\frac{1}{4}\|\tilde{\boldsymbol{v}}\|_{\boldsymbol{V}_{h}}^{2}+\delta_{1}^{2}\|\boldsymbol{w}_{1}\|_{\boldsymbol{V}_{h}}^{2}=\frac{1}{4}\|\tilde{\boldsymbol{v}}\|_{\boldsymbol{V}_{h}}^{2}+\delta_{1}^{2}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}.$

Moreover, we can obtain

$\displaystyle\delta_{0}|\left(\alpha
r^{0},\operatorname{div}\tilde{\boldsymbol{v}}\right)_{\Omega}|$
$\displaystyle\leq\delta_{0}\|r^{0}\|_{0,\lambda^{-1}}\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}\leq
C_{r}\delta_{0}\|\tilde{q}^{0}\|_{0,\lambda^{-1}}\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}$
(57) $\displaystyle\leq
C_{r}^{2}\delta_{0}^{2}\|\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}+\frac{1}{4}\|\boldsymbol{v}\|_{\boldsymbol{V}_{h}}^{2},$
$\displaystyle\delta_{0}|\left(s_{0}r^{0},\tilde{q}\right)_{\Omega}|$
$\displaystyle\leq
C_{r}^{2}\delta_{0}^{2}\|\tilde{q}^{0}\|_{0,s_{0}}^{2}+\frac{1}{2}\|\tilde{q}\|_{0,s_{0}}^{2}$
(58) $\displaystyle\leq
C_{r}^{2}C_{P}^{2}\delta_{0}^{2}S(\tilde{q},\tilde{q})+\frac{1}{2}\|\tilde{q}\|_{0,s_{0}}^{2},$
(59) $\displaystyle\delta_{0}|S(r^{0},\tilde{q})|$ $\displaystyle\leq
C_{r}\delta_{0}S(\tilde{q},\tilde{q}).$

If we use (55)–(59) to (51), then

$\displaystyle\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle\geq\frac{1}{2}\|\tilde{\boldsymbol{v}}\|_{\boldsymbol{V}_{h}}^{2}-\delta_{1}^{2}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}+\frac{\delta_{1}C_{1}^{\prime\prime}}{2}\|\tilde{q}^{c}\|_{0,\lambda^{-1}}^{2}$
$\displaystyle\quad-\delta_{1}\left(C_{P}C_{1}^{\prime\prime}+C_{2}^{\prime\prime}\right)S(\tilde{q},\tilde{q})$
$\displaystyle\quad-
C_{r}^{2}\delta_{0}^{2}\|\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}-C_{r}^{2}C_{P}^{2}\delta_{0}^{2}\|\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}+\frac{1}{2}\|\tilde{q}\|_{0,s_{0}}^{2}$
$\displaystyle\quad+S(\tilde{q},\tilde{q})-C_{r}\delta_{0}S(\tilde{q},\tilde{q})+\frac{C_{3}\Delta
t}{2}\|\tilde{q}\|_{h}^{2}.$

Since
$\|\tilde{q}^{c}\|_{0,s_{0}}^{2}=\|\tilde{q}+(-\tilde{q}^{0})\|_{0,s_{0}}^{2}\leq
2(\|\tilde{q}\|_{0,s_{0}}^{2}+\|\tilde{q}^{0}\|_{0,s_{0}}^{2})$, we have

$\displaystyle\frac{1}{2}\|\tilde{q}\|_{0,s_{0}}^{2}\geq\delta_{0}\|\tilde{q}\|_{0,s_{0}}^{2}\geq\frac{\delta_{0}}{2}\|\tilde{q}^{c}\|_{0,s_{0}}^{2}-\delta_{0}\|\tilde{q}^{0}\|_{0,s_{0}}^{2}$

if $\delta_{0}\leq\frac{1}{2}$. If we apply this to the previous inequality
assuming that $\delta_{0}\leq\frac{1}{2}$, then we have

$\displaystyle\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))$
$\displaystyle\geq\frac{1}{2}\|\tilde{\boldsymbol{v}}\|_{\boldsymbol{V}_{h}}^{2}-\delta_{1}^{2}\|\tilde{q}\|_{0,\lambda^{-1}}^{2}+\frac{\delta_{1}C_{1}^{\prime\prime}}{2}\|\tilde{q}^{c}\|_{0,\lambda^{-1}}^{2}$
$\displaystyle\quad-\delta_{1}\left(\frac{C_{P}C_{1}^{\prime\prime}}{2}+C_{2}^{\prime\prime}\right)S(\tilde{q},\tilde{q})$
$\displaystyle\quad-
C_{r}^{2}\delta_{0}^{2}\|\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}-C_{r}^{2}C_{P}^{2}\delta_{0}^{2}\|\tilde{q}^{0}\|_{0,\lambda^{-1}}^{2}+\frac{\delta_{0}}{2}\|\tilde{q}^{c}\|_{0,s_{0}}^{2}-\delta_{0}\|\tilde{q}^{0}\|_{0,s_{0}}^{2}$
$\displaystyle\quad+S(\tilde{q},\tilde{q})-C_{r}\delta_{0}S(\tilde{q},\tilde{q})+\frac{C_{3}\Delta
t}{2}\|\tilde{q}\|_{h}^{2}.$

If we choose a sufficiently small $\delta_{0}$, $\delta_{1}$ satisfying

$\displaystyle\delta_{1}\left(C_{P}C_{1}^{\prime\prime}+C_{2}^{\prime\prime}\right)+C_{r}\delta_{0}\leq\frac{1}{4},$
$\displaystyle\max\\{\delta_{1}^{2}+C_{r}^{2}\delta_{0}^{2}+C_{r}^{2}C_{P}^{2}\delta_{0}^{2},\delta_{0}\\}C_{P}^{2}\leq\frac{1}{4},$

we obtain

$\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))\\\
\geq\frac{1}{2}\|\tilde{\boldsymbol{v}}\|_{\boldsymbol{V}_{h}}^{2}+\frac{\delta_{1}C_{1}^{\prime\prime}}{2}\|\tilde{q}^{c}\|_{0,\lambda^{-1}}^{2}+\frac{\delta_{0}}{2}\|\tilde{q}^{c}\|_{0,s_{0}}^{2}+\frac{1}{2}S(\tilde{q},\tilde{q})+\frac{C_{3}\Delta
t}{2}\|\tilde{q}\|_{h}^{2}.$

Since $\|\tilde{q}^{0}\|_{0,\lambda^{-1}}+\|\tilde{q}^{0}\|_{0,s_{0}}\leq
C_{P}S(\tilde{q}^{0},\tilde{q}^{0})^{1/2}$, the above inequality gives

$\displaystyle\mathcal{B}((\boldsymbol{v},q),(\tilde{\boldsymbol{v}},\tilde{q}))\geq
C\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\mathcal{X}}^{2}$

with $C>0$ depending on $C_{1}^{\prime\prime}$, $C_{2}$, $C_{3}$,
$\delta_{0}$, $\delta_{1}$.

On the other hand, by the definition of $(\boldsymbol{v},q)$, it is easy to
check that $\|(\boldsymbol{v},q)\|_{\mathcal{X}}\leq
C\|(\tilde{\boldsymbol{v}},\tilde{q})\|_{\mathcal{X}}$ with $C>0$ depending on
$\delta_{0}$, $\delta_{1}$. Therefore, the conclusion follows. ∎

| $\|u-u_{h}\|_{a_{\boldsymbol{u}}}/\|u\|_{a_{\boldsymbol{u}}}$ | $\|u-u_{h}\|_{1}/\|u\|_{1}$ | $\|u-u_{h}\|_{0}/\|u\|_{0}$ | $\|p-p_{h}\|_{1}/\|p\|_{1}$ | $\|p-p_{h}\|_{0}/\|p\|_{0}$
---|---|---|---|---|---
$\beta$ | $\nu$ | $N$ | error | rate | error | rate | error | rate | error | rate | error | rate
1 | 0.3 | $8$ | 1.1861e-01 | – | 1.0361e-01 | – | 1.0997e-02 | – | 1.0178e-01 | – | 9.9130e-03 | –
$16$ | 5.9533e-02 | 0.99 | 5.1796e-02 | 1.00 | 2.8000e-03 | 1.97 | 5.1074e-02 | 0.99 | 2.5883e-03 | 1.94
$32$ | 2.9788e-02 | 1.00 | 2.5882e-02 | 1.00 | 7.0630e-04 | 1.99 | 2.5545e-02 | 1.0 | 6.9010e-04 | 1.91
$64$ | 1.4895e-02 | 1.00 | 1.2935e-02 | 1.00 | 1.7740e-04 | 1.99 | 1.2767e-02 | 1.0 | 1.7800e-04 | 1.95
0.499 | $8$ | 1.2029e-01 | – | 1.0223e-01 | – | 1.1052e-02 | – | 1.0180e-01 | – | 1.0211e-02 | –
$16$ | 6.0466e-02 | 0.99 | 5.1133e-02 | 1.00 | 2.8143e-03 | 1.97 | 5.1080e-02 | 0.99 | 2.6653e-03 | 1.94
$32$ | 3.0272e-02 | 1.00 | 2.5555e-02 | 1.00 | 7.1040e-04 | 1.99 | 2.5542e-02 | 1.0 | 6.8090e-04 | 1.97
$64$ | 1.5141e-02 | 1.00 | 1.2772e-02 | 1.00 | 1.7850e-04 | 1.99 | 1.2767e-02 | 1.0 | 1.7460e-04 | 1.96
2 | 0.3 | $8$ | 1.1861e-01 | – | 1.0361e-01 | – | 1.0997e-02 | – | 1.0179e-01 | – | 9.9147e-03 | –
$16$ | 5.9533e-02 | 0.99 | 5.1796e-02 | 1.00 | 2.8001e-03 | 1.97 | 5.1082e-02 | 0.99 | 2.5884e-03 | 1.94
$32$ | 2.9788e-02 | 1.00 | 2.5882e-02 | 1.00 | 7.0630e-04 | 1.99 | 2.5546e-02 | 1.0 | 6.9010e-04 | 1.91
$64$ | 1.4895e-02 | 1.00 | 1.2935e-02 | 1.00 | 1.7740e-04 | 1.99 | 1.2767e-02 | 1.0 | 1.7800e-04 | 1.95
0.499 | $8$ | 1.2029e-01 | – | 1.0223e-01 | – | 1.1052e-02 | – | 1.0181e-01 | – | 1.0213e-02 | –
$16$ | 6.0466e-02 | 0.99 | 5.1133e-02 | 1.00 | 2.8143e-03 | 1.97 | 5.1080e-02 | 0.99 | 2.6654e-03 | 1.94
$32$ | 3.0272e-02 | 1.00 | 2.5545e-02 | 1.00 | 7.1040e-04 | 1.99 | 2.5543e-02 | 1.0 | 6.8090e-04 | 1.97
$64$ | 1.5141e-02 | 1.00 | 1.2772e-02 | 1.00 | 1.7850e-04 | 1.99 | 1.2767e-02 | 1.0 | 1.7460e-04 | 1.96
Table 1.  The convergence results of errors with the manufactured solution (60) for $\beta=1,2$ and the Crank–Nicolson scheme with $\Delta t=1/N$.  | $N~{}(DOFs)$
---|---
$\beta$ | $\Delta t$ | $\nu$ | $8~{}(3201)$ | $16~{}(12545)$ | $32~{}(49665)$ | $64~{}(197633)$
$0$ | $10^{-1}$ | $0.3$ | $148$ | $195$ | $243$ | $307$
$0.499$ | $185$ | $192$ | $185$ | $163$
$10^{-2}$ | $0.3$ | $152$ | $210$ | $276$ | $343$
$0.499$ | $214$ | $254$ | $252$ | $231$
$10^{-3}$ | $0.3$ | $152$ | $218$ | $295$ | $391$
$0.499$ | $228$ | $313$ | $333$ | $336$
$1$ | $10^{-1}$ | $0.3$ | $58$ | $57$ | $57$ | $57$
$0.499$ | $91$ | $92$ | $93$ | $98$
$10^{-2}$ | $0.3$ | $60$ | $60$ | $60$ | $58$
$0.499$ | $94$ | $97$ | $102$ | $102$
$10^{-3}$ | $0.3$ | $61$ | $62$ | $59$ | $61$
$0.499$ | $97$ | $100$ | $105$ | $106$
$2$ | $10^{-1}$ | $0.3$ | $39$ | $36$ | $35$ | $33$
$0.499$ | $62$ | $61$ | $63$ | $62$
$10^{-2}$ | $0.3$ | $39$ | $38$ | $36$ | $35$
$0.499$ | $63$ | $63$ | $64$ | $64$
$10^{-3}$ | $0.3$ | $40$ | $39$ | $38$ | $36$
$0.499$ | $65$ | $64$ | $66$ | $65$
Table 2. The numbers of iteration for $\beta=0,1,2$ of the MinRes algorithm
with abstract preconditioners of the form (44). See Section 6 for more details
on implementation.

##  6\. Numerical results

In this section we present numerical results which support our theoretical
findings on error analysis and preconditioning. All numerical experiments are
done with the Firedrake package [30].

In our numerical experiments we set $\Omega=[0,1]\times[0,1]$ with
$\Gamma_{d}=\\{0\\}\times[0,1]$ and $\Gamma_{p}=\partial\Omega$. We use the
structured meshes obtained by dividing $\Omega$ into $N\times N$ subsquares
and the dividing each subsquare into two triangles. For convergence rate tests
we refine the meshes with $N=8,16,32,64$.

To check convergence rates of the errors of numerical and exact solutions, we
consider a manufactured solution with

(60a) $\displaystyle\boldsymbol{u}$ $\displaystyle=\left(-\pi
x^{2}\left(1-x\right)^{2}\sin^{2}\left(\pi y\right)\cos\left(2t\right),-\pi
y^{2}\left(1-y\right)^{2}\sin^{2}\left(\pi
x\right)\cos\left(2t\right)\right),$ (60b) $\displaystyle p$
$\displaystyle=x\left(1-x\right)y\left(1-y\right)\cos\left(t\right).$

Note that $\operatorname{div}\boldsymbol{u}\not=0$, so the $L^{2}$ norm of
$\boldsymbol{f}$ may grow unboundedly as $\lambda\rightarrow+\infty$. In order
to illustrate convergence rates of the errors independent of $\lambda$, we
compute the relative errors measured by

$\displaystyle\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{a_{\boldsymbol{u}}}/\|\boldsymbol{u}\|_{a_{\boldsymbol{u}}},\quad\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1}/\|\boldsymbol{u}\|_{1},\quad\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0}/\|\boldsymbol{u}\|_{0},$
$\displaystyle\|p-p_{h}\|_{1}/\|p\|_{1},\quad\|p-p_{h}\|_{0}/\|p\|_{0}.$

In Table 1 we present the results for $\beta=1,2$ cases with the
Crank–Nicolson time discretization with time step size $\Delta t=1/N$. For all
relative error quantities and for $\nu=0.3,0.499$ we observe that optimal
convergence rates are obtained.

For preconditioning we also present preconditioning test results. We compare
the results for all combinations of $\beta=0,1,2$, $\nu=0.3,0.499$, $\Delta
t=10^{-1},10^{-2},10^{-3}$. We use preconditioners which have the abstract
form (44). For implementation of $\mathcal{P}_{\boldsymbol{V}_{h}}$ we use the
geometric multigrid preconditioner developed in [2]. For implementation of
$\mathcal{P}_{\tilde{Q}_{h}}$ we use the preconditioners proposed in [21, 18]
combining algebraic multigrid preconditioners for $Q_{h}^{c}$ and the block
Jacobi preconditioner for $Q_{h}^{0}$. Numbers of iterations are given in
Table LABEL:table:scalar-kappa-precond. One can observe that $\beta=0$ case,
the original enriched Galerkin method for $p$, shows strong dependence on mesh
refinement and $\nu\approx 0.5$. For $\beta=1$, the dependence on mesh
refinement disappears but the numbers of iteration still show some dependence
for $\nu\approx 0.5$. For $\beta=2$, the numbers of iteration are $\approx
65$. We remark that these numbers are usually regarded as optimal iteration
numbers in parameter-robust preconditioning with MinRes algorithm [20]. The
iteration numbers become much smaller in practice if GmRes algorithm is used
but the norm equivalence discussion in Section 5 is not valid for GmRes, so we
do not include GmRes in this paper.

##  7\. Conclusion

In this paper we developed a new discretization for poroelasticity. We use the
primal formulations for the elastic displacement and for the pore pressure
minimizing the number of variables. In our discretization with nonconforming
and the interior over-stabilized enriched Galerkin methods, numerical
solutions are robust for nearly incompressible elasticity parameters, are
locally mass conservative. Furthermore, we showed a theoretical discussion how
to construct robust preconditioinors and the robust preconditioner
construction method was verified by numerical experiments.

## Statements and Declarations

Funding Jeonghun J. Lee is partially supported by Baylor’s Undergraduate
Research Council grant, Baylor’s start-up grant, and National Science
Foundation (DMS-2110781)

Conflict of Interest The authors declared that they have no conflict of
interest.

## References

  * [1] D. N. Arnold and F. Brezzi, _Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates_ , RAIRO Modél. Math. Anal. Numér. 19 (1985), no. 1, 7–32.
  * [2] Francis R. A. Aznaran, Patrick E. Farrell, and Robert C. Kirby, _Transformations for piola-mapped elements_ , arXiv preprint arXiv:2110.13224 (2021).
  * [3] M. A. Biot, _Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range_ , J. Acoust. Soc. Amer. 28 (1956), 168–178. MR 134056
  * [4] by same author, _Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range_ , J. Acoust. Soc. Amer. 28 (1956), 179–191. MR 134057
  * [5] Daniele Boffi, Michele Botti, and Daniele A. Di Pietro, _A nonconforming high-order method for the Biot problem on general meshes_ , SIAM J. Sci. Comput. 38 (2016), no. 3, A1508–A1537. MR 3504993
  * [6] Lorenzo Botti, Michele Botti, and Daniele A. Di Pietro, _An abstract analysis framework for monolithic discretisations of poroelasticity with application to hybrid high-order methods_ , Comput. Math. Appl. 91 (2021), 150–175. MR 4253882
  * [7] Susanne C. Brenner, _Poincaré-Friedrichs inequalities for piecewise $H^{1}$ functions_, SIAM J. Numer. Anal. 41 (2003), no. 1, 306–324. MR 1974504
  * [8] Annalisa Buffa and Christoph Ortner, _Compact embeddings of broken Sobolev spaces and applications_ , IMA J. Numer. Anal. 29 (2009), no. 4, 827–855. MR 2557047
  * [9] X. Feng, Z. Ge, and Y. Li, _Analysis of a multiphysics finite element method for a poroelasticity model_ , IMA J. Numer. Anal. 38 (2018), no. 1, 330–359.
  * [10] G. Fu, _A high-order HDG method for the Biot’s consolidation model_ , Comput. Math. Appl. (2018).
  * [11] Qingguo Hong and Johannes Kraus, _Parameter-robust stability of classical three-field formulation of Biot’s consolidation model_ , Electron. Trans. Numer. Anal. 48 (2018), 202–226. MR 3820123
  * [12] Xiaozhe Hu, Carmen Rodrigo, Francisco J. Gaspar, and Ludmil T. Zikatanov, _A nonconforming finite element method for the Biot’s consolidation model in poroelasticity_ , J. Comput. Appl. Math. 310 (2017), 143–154. MR 3544596
  * [13] G. Kanschat and B. Rivière, _A finite element method with strong mass conservation for Biot’s linear consolidation model_ , J. Sci. Comput. (2018).
  * [14] Johannes Kraus, Philip L. Lederer, Maria Lymbery, and Joachim Schöberl, _Uniformly well-posed hybridized discontinuous Galerkin/hybrid mixed discretizations for Biot’s consolidation model_ , Comput. Methods Appl. Mech. Engrg. 384 (2021), Paper No. 113991, 23. MR 4274934
  * [15] Bishnu P. Lamichhane, _A nonconforming finite element method for the Stokes equations using the Crouzeix-Raviart element for the velocity and the standard linear element for the pressure_ , Internat. J. Numer. Methods Fluids 74 (2014), no. 3, 222–228. MR 3148255
  * [16] Jeonghun J. Lee, _Robust error analysis of coupled mixed methods for Biot’s consolidation model_ , Journal of Scientific Computing 69 (2016), no. 2, 610–632.
  * [17] Jeonghun J. Lee, _Robust three-field finite element methods for Biot’s consolidation model in poroelasticity_ , BIT (2017).
  * [18] Jeonghun J. Lee and Omar Ghattas, _Interior over-stabilized enriched Galerkin methods for second order elliptic equations_ , arXiv preprint arXiv:1706.05169 (2022).
  * [19] Jeonghun J. Lee, Kent-Andre Mardal, and Ragnar Winther, _Parameter-robust discretization and preconditioning of Biot’s consolidation model_ , SIAM Journal on Scientific Computing 39 (2017), no. 1, A1–A24.
  * [20] by same author, _Parameter-robust discretization and preconditioning of Biot’s consolidation model_ , SIAM J. Sci. Comp. 39 (2017), no. 1, A1–A24.
  * [21] Sanghyun Lee, Young-Ju Lee, and Mary F. Wheeler, _A locally conservative enriched Galerkin approximation and efficient solver for elliptic and parabolic problems_ , SIAM J. Sci. Comput. 38 (2016), no. 3, A1404–A1429. MR 3499558
  * [22] Sanghyun Lee and Son-Young Yi, _Locking-free and locally-conservative enriched Galerkin method for poroelasticity_ , J. Sci. Comput. 94 (2023), no. 1, Paper No. 26, 23. MR 4527558
  * [23] Kent Andre Mardal, Xue-Cheng Tai, and Ragnar Winther, _A robust finite element method for Darcy-Stokes flow_ , SIAM J. Numer. Anal. 40 (2002), no. 5, 1605–1631. MR 1950614
  * [24] Kent-Andre Mardal and Ragnar Winther, _An observation on Korn’s inequality for nonconforming finite element methods_ , Math. Comp. 75 (2006), no. 253, 1–6. MR 2176387
  * [25] by same author, _Preconditioning discretizations of systems of partial differential equations_ , Numer. Linear Algebra Appl. 18 (2011), no. 1, 1–40. MR 2769031
  * [26] Márcio A Murad, Vidar Thomée, and Abimael FD Loula, _Asymptotic behavior of semidiscrete finite-element approximations of Biot’s consolidation problem_ , SIAM J. Numer. Anal. 33 (1996), no. 3, 1065–1083.
  * [27] R. Oyarzúa and R. Ruiz-Baier, _Locking-free finite element methods for poroelasticity_ , SIAM J. Numer. Anal. 54 (2016), no. 5, 2951–2973.
  * [28] P. J. Phillips and M. F. Wheeler, _Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach_ , Computat. Geosci. 13 (2009), 5–12.
  * [29] Phillip Joseph Phillips and Mary F. Wheeler, _A coupling of mixed and discontinuous Galerkin finite-element methods for poroelasticity_ , Comput. Geosci. 12 (2008), no. 4, 417–435. MR 2461315
  * [30] Florian Rathgeber, David A. Ham, Lawrence Mitchell, Michael Lange, Fabio Luporini, Andrew T. T. McRae, Gheorghe-Teodor Bercea, Graham R. Markall, and Paul H. J. Kelly, _Firedrake: automating the finite element method by composing abstractions_ , ACM Trans. Math. Software 43 (2017), no. 3, Art. 24, 27. MR 3615280
  * [31] B. Rivière, Jun Tan, and Travis Thompson, _Error analysis of primal discontinuous Galerkin methods for a mixed formulation of the Biot equations_ , Comput. Math. Appl. 73 (2017), no. 4, 666 – 683.
  * [32] C. Rodrigo, F. J. Gaspar, X. Hu, and L. T. Zikatanov, _Stability and monotonicity for some discretizations of the Biot’s consolidation model_ , Comput. Methods Appl. Mech. Engrg. 298 (2016), 183–204. MR 3427711
  * [33] Juan Enrique Santos, _Elastic wave propagation in fluid-saturated porous media. I. The existence and uniqueness theorems_ , RAIRO Modél. Math. Anal. Numér. 20 (1986), no. 1, 113–128. MR 844519
  * [34] Juan Enrique Santos and Ernesto Jorge Oreña, _Elastic wave propagation in fluid-saturated porous media. II. The Galerkin procedures_ , RAIRO Modél. Math. Anal. Numér. 20 (1986), no. 1, 129–139. MR 844520
  * [35] R. E. Showalter, _Diffusion in poro-elastic media_ , J. Math. Anal. Appl. 251 (2000), no. 1, 310–340. MR 1790411
  * [36] by same author, _Diffusion in poro-elastic media_ , J. Math. Anal. Appl. 251 (2000), no. 1, 310–340. MR 1790411
  * [37] Xue-Cheng Tai and Ragnar Winther, _A discrete de Rham complex with enhanced smoothness_ , Calcolo 43 (2006), no. 4, 287–306. MR 2283095
  * [38] Alexander Ženíšek, _The existence and uniqueness theorem in Biot’s consolidation theory_ , Apl. Mat. 29 (1984), no. 3, 194–211. MR 747212
  * [39] Herbert F. Wang, _Theory of linear poroelasticity with applications to geomechanics and hydrogeology_ , Princeton Series in Geophysics, Princeton University Press, 2000.
  * [40] S.-Y. Yi, _A coupling of nonconforming and mixed finite element methods for Biot’s consolidation model_ , Numer. Meth. Part. D. E. 29 (2013), no. 5, 1749–1777.
  * [41] Son-Young Yi, _Convergence analysis of a new mixed finite element method for Biot’s consolidation model_ , Numer. Methods Partial Differential Equations 30 (2014), no. 4, 1189–1210. MR 3200272
  * [42] by same author, _A study of two modes of locking in poroelasticity_ , SIAM J. Numer. Anal. 55 (2017), no. 4, 1915–1936.
  * [43] Son-Young Yi, Sanghyun Lee, and Ludmil Zikatanov, _Locking-free enriched Galerkin method for linear elasticity_ , SIAM J. Numer. Anal. 60 (2022), no. 1, 52–75. MR 4358467
  * [44] Kosaku Yosida, _Functional analysis_ , 6th ed., Springer Classics in Mathematics, Springer-Verlag, 1980.
  * [45] O. C. Zienkiewicz and T. Shiomi, _Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution_ , International Journal for Numerical and Analytical Methods in Geomechanics 8 (1984), no. 1, 71–96.