# A posteriori error estimates of Darcy flows with Robin-type jump interface
conditions

Jeonghun J. Lee  Department of Mathematics, Baylor University, Waco, Texas,
USA<EMAIL_ADDRESS>
(Date: April, 2023)

###### Abstract.

In this work we develop an a posteriori error estimator for mixed finite
element methods of Darcy flow problems with Robin-type jump interface
conditions. We construct an energy-norm type a posteriori error estimator
using the Stenberg post-processing. The reliability of the estimator is proved
using an interface-adapted Helmholtz-type decomposition and an interface-
adapted Scott–Zhang type interpolation operator. A local efficiency and the
reliability of post-processed pressure are also proved. Numerical results
illustrating adaptivity algorithms using our estimator are included.

###### Key words and phrases:

mixed finite element methods, a posteriori error estimates, Robin boundary
conditions

###### 2000 Mathematics Subject Classification:

Primary: 65N30, 65N15

##  1\. Introduction

Fluid flow in porous media appears in various fields of science and
engineering applications. Therefore, mathematical modeling and numerical
methods for finding accurate numerical solutions of porous media flow have
been important problems in computational mathematics. Recently, mathematical
models in which porous media domains have low-dimensional fault (or fracture)
structures are considered for accurate descriptions of more realistic porous
media flow. In [23], some porous media flow models with fault/fracture
structures were proposed in which fluid flow on fractures and in surrounding
porous media are governed by separate partial differential equations with
coupling conditions. In [19], a reduced model was derived under the assumption
that there is no fluid flow along fault/fracture structures because of very
low permeability on fault/fracture. In the reduced models, the fluid flow and
the pressure jump on faults are related by a Robin-type interface condition.
We remark that similar models with nonlinear extensions are used for porous
media flows with semi-permeable membrane structures in consideration of their
applications to chemical processes in biological tissues (see, e.g., [11,
12]).

The purpose of this paper is to obtain a posteriori error estimators for the
model derived in [19] with the dual mixed form of finite element methods. We
remark that a posteriori error estimate results for the more complex models in
[23, 2] (see [13, 16, 31] for a posteriori error estimates), do not imply a
posteriori error estimate results for the model that we are interested in.
This is because a less number of error terms makes local efficiency of a
posteriori error estimators more difficult.

We also remark that the problem in this paper can be viewed as a
generalization of mixed finite element methods for Poisson equations with
Robin boundary conditions which was studied in [18]. A priori and a posteriori
error estimates are done in [18] using the mesh-dependent norm approach (cf.
[5, 22, 29]) but the saturation assumption is necessary for the reliability of
the estimator. The analysis in this paper does not need such an assumption for
reliability, and it also gives a new reliability estimate for post-processed
pressure.

The paper is organized as follows. In Section 2 we present background notions
on function spaces, governing equations, finite element discretization. We
define our a posteriori error estimator and prove its reliability and local
efficiency in Section 3. In Section 4 and 5, we present numerical experiment
results which show performance of our a posteriori error estimator, and
conclusion with summary. Finally, some calculus identities which are used in
our analysis are explained in Section 6 as appendix.

##  2\. Preliminaries for governing equations

###  2.1. Notation and definitions

For a bounded domain $D\subset\mathbb{R}^{n}$ ($n=2,3$) with positive
$n$-dimensional Lebesgue measure, we use the convention that
$(u,v)_{D}=\int_{D}uv\,dx$ for a subdomain $D\subset\overline{\Omega}$ which
has positive $n$-dimensional Lebesgue measure. Similarly,
$\left<u,v\right>_{D}=\int_{D}uv\,dS$ for a subdomain
$D\subset\overline{\Omega}$ which has positive $(n-1)$ or $(n-2)$-dimensional
Lebesgue measure up to context.

###  2.2. Governing equations and variational formulation

Let $\Omega\subset\mathbb{R}^{n}$, $n=2,3$ be a homologically trivial bounded
domain with polygonal/polyhedral boundary. We assume that a fault $\Gamma$ is
a union of disjoint ($n-1$)-dimensional piecewise linear submanifolds in
$\Omega$. Each connected component of $\Gamma$ is a union of linear segments
(if $n=2$) or as a union of planar domains such that the boundary of each
planar domain is a union of linear segments. We also assume that there are two
open subdomains $\Omega_{+},\Omega_{-}\subset(\Omega\setminus\Gamma)$ with
Lipschitz boundaries such that

$\displaystyle\overline{\Omega}=\overline{\Omega_{+}}\cup\overline{\Omega_{-}},\qquad\Gamma\subset\partial\Omega_{+}\cap\partial\Omega_{-},$

and only one side of $\Gamma$ is in contact with $\Omega_{+}$ or $\Omega_{-}$.
Let $\boldsymbol{n}_{+}$ and $\boldsymbol{n}_{-}$ be the two unit normal
vector fields on $\Gamma$ with opposite directions
($\boldsymbol{n}_{+}=-\boldsymbol{n}_{-}$) such that $\boldsymbol{n}_{\pm}$
correspond to the unit outward normal vector fields from $\Omega{\pm}$ (see
Figure 1).

Suppose that $\Gamma_{D}$, $\Gamma_{N}$ are disjoint ($n-1$)-dimensional open
submanifolds in $\partial\Omega$ such that
$\overline{\Gamma_{D}}\cup\overline{\Gamma_{N}}=\partial\Omega$. In this paper
we assume the following:

$\displaystyle\text{ For any }q\in L^{2}(\Omega)\text{ there exists
}\boldsymbol{w}\in H^{1}(\Omega;{\mathbb{R}}^{n})$ (1) $\displaystyle\text{
such that
}\boldsymbol{w}|_{\Gamma\cup\Gamma_{N}}=0,\operatorname{div}\boldsymbol{w}=q\text{
and }\|\boldsymbol{w}\|_{1}\leq C\|q\|_{0}$ $\displaystyle\text{ with a
constant }C>0\text{ depending on }\Omega,\Gamma,\Gamma_{N}.$

The assumption (1) is a weak assumption. For example, if both of
$\partial\Omega_{+}\cap\partial\Omega\cap\Gamma_{N}$ and
$\partial\Omega_{-}\cap\partial\Omega\cap\Gamma_{N}$ have positive
($n-1$)-dimensional Lebesgue measures, then (1) holds.

For any $q\in L^{2}(\Omega)$ with sufficient regularity, we use
$q_{+}|_{\Gamma}$ and $q_{-}|_{\Gamma}$ to denote the traces of
$q|_{\Omega_{+}}$ and $q|_{\Omega_{-}}$ on $\Gamma$. For simplicity we use
$\left[\\!\left[{q}\right]\\!\right]|_{\Gamma}:=q_{+}|_{\Gamma}-q_{-}|_{\Gamma}$.
Note that the continuity of $q$ on $\Gamma$ is not assumed in general, so
$\left[\\!\left[{q}\right]\\!\right]|_{\Gamma}\not=0$. For a vector-valued
function $\boldsymbol{v}$ on $\Omega$ with enough regularity (e.g.,
$\boldsymbol{v}\in H^{s}(\Omega\setminus\Gamma;\mathbb{R}^{n})$ with $s>1/2$),
$\boldsymbol{v}_{+}|_{\Gamma}$ and $\boldsymbol{v}_{-}|_{\Gamma}$ are well-
defined as the traces of $\boldsymbol{v}$ from $\Omega_{+}$ and $\Omega_{-}$.
We say that $\boldsymbol{v}$ satisfies normal continuity on $\Gamma$ if
$\boldsymbol{v}_{+}|_{\Gamma}\cdot\boldsymbol{n}_{+}=-\boldsymbol{v}_{-}|_{\Gamma}\cdot\boldsymbol{n}_{-}$
on $\Gamma$.

$\Gamma$$n_{+}$$n_{-}$$\Omega_{+}$$\Omega_{-}$ Figure 1. A model domain with
interface $\Gamma$

For governing equations assume that $\kappa$ is a symmetric positive definite
tensor on $\Omega$. In Darcy flow problems, the pressure $p$ and fluid flow
$\boldsymbol{u}$ satisfy Darcy’s law $\boldsymbol{u}=-\kappa\nabla p$ in
$\Omega$. Conservation of mass gives $\operatorname{div}\boldsymbol{u}=f$ for
given source/sink function $f$ on $\Omega$. The pressure and flux boundary
conditions are given by

$\displaystyle p=g_{D}\text{ on
}\Gamma_{D},\quad{\boldsymbol{u}}\cdot\boldsymbol{n}=g_{N}\text{ on
}\Gamma_{N},$

and the interface condition on $\Gamma$ is
$\alpha\boldsymbol{u}_{+}\cdot\boldsymbol{n}_{+}-\left[\\!\left[{p}\right]\\!\right]|_{\Gamma}=0$
with $\alpha>0$. Summarizing these equations and conditions, a strong form of
equations with a dual mixed formulation of the Darcy flow equation in domain
$\Omega$ with fault $\Gamma$ reads:

(2) $\displaystyle\kappa^{-1}\boldsymbol{u}+\nabla p=0\text{ in
}\Omega,\qquad\operatorname{div}\boldsymbol{u}=f\text{ in }\Omega,$ (3)
$\displaystyle\boldsymbol{u}\cdot\boldsymbol{n}=g_{N}\text{ on
}\Gamma_{N},\quad p=g_{D}\text{ on }\Gamma_{D},$ (4)
$\displaystyle\alpha\boldsymbol{u}_{+}\cdot\boldsymbol{n}_{+}-\left[\\!\left[{p}\right]\\!\right]=0\text{
on }\Gamma.$

Throughout this paper we assume that $\alpha$ is constant on $\Gamma$ and

(5) $\displaystyle 1\lesssim\alpha<\infty$

with a uniform lower bound, and we do not consider the limit cases $\alpha\to
0^{+}$ or $\alpha\to+\infty$. The limit case $\alpha=0$ becomes the classical
Darcy flow problems without fault which does not need the interface condition
(4). The $\alpha=+\infty$ case corresponds to the problem that no fluid flows
across $\Gamma$ which needs $\boldsymbol{u}\cdot\boldsymbol{n}|_{\Gamma}=0$ as
an interface condition. This case needs a completely different way to
implement the interface condition
$\boldsymbol{u}\cdot\boldsymbol{n}|_{\Gamma}=0$ with the dual mixed finite
element methods, the numerical method that we use in this paper. Therefore,
$\alpha=+\infty$ case cannot be covered by the work in this paper. However,
our analysis does not need a uniform upper bound of $\alpha$, so the results
in the paper are valid for nearly impermeable $\Gamma$, i.e., for arbitrarily
large but finite $\alpha$.

Hereafter, we assume $\partial\Omega=\Gamma_{D}$, $g_{D}=0$, $\kappa=1$ for
simplicity of discussions. Let $Q=L^{2}(\Omega)$, and
$H(\operatorname{div},\Omega)$ be the space of ${\mathbb{R}}^{n}$-valued
$L^{2}$ functions on $\Omega$ such that its distributional divergence is in
$L^{2}(\Omega)$. We define

$\displaystyle\boldsymbol{V}:=\\{\boldsymbol{v}\in
H(\operatorname{div},\Omega):\boldsymbol{v}_{+}\cdot\boldsymbol{n}_{+}|_{\Gamma}=-\boldsymbol{v}_{-}\cdot\boldsymbol{n}_{-}|_{\Gamma}\in
L^{2}(\Gamma)\\}$

with two norms

(6)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{v}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}$
$\displaystyle=\left(\int_{\Omega}\kappa^{-1}\boldsymbol{v}\cdot\boldsymbol{v}\,dx+\sum_{F\subset\Gamma}\int_{F}\alpha(\boldsymbol{v}\cdot\boldsymbol{n})(\boldsymbol{v}\cdot\boldsymbol{n})\,ds\right)^{\frac{1}{2}},$
(7) $\displaystyle\|\boldsymbol{v}\|_{\boldsymbol{V}}$
$\displaystyle=\left(\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{v}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}+\|\operatorname{div}\boldsymbol{v}\|_{0}^{2}\right)^{\frac{1}{2}}.$

By multiplying $\boldsymbol{v}\in\boldsymbol{V}$ to the first equation in (2)
and taking the integration by parts with $g_{D}=0$, $\kappa=1$,

$\displaystyle\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}\,dx+\int_{\Omega}\nabla
p\cdot\boldsymbol{v}\,dx$
$\displaystyle=\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}\,dx+\int_{\Gamma}p_{+}\boldsymbol{v}_{+}\cdot\boldsymbol{n}_{+}\,ds+\int_{\Gamma}p_{-}\boldsymbol{v}_{-}\cdot\boldsymbol{n}_{-}\,ds$
$\displaystyle\quad-\int_{\Omega}p\operatorname{div}\boldsymbol{v}\,dx.$

After using
$\boldsymbol{v}_{+}\cdot\boldsymbol{n}_{+}=-\boldsymbol{v}_{-}\cdot\boldsymbol{n}_{-}$
on $\Gamma$, and the interface condition (4), we obtain

$\displaystyle\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}\,dx-\int_{\Omega}p\operatorname{div}\boldsymbol{v}\,dx+\int_{\Gamma}\alpha(\boldsymbol{u}_{+}\cdot\boldsymbol{n}_{+})(\boldsymbol{v}_{+}\cdot\boldsymbol{n}_{+})\,ds=0$

which can be written as
$\left(\boldsymbol{u},\boldsymbol{v}\right)_{\Omega}-\left(p,\operatorname{div}\boldsymbol{v}\right)_{\Omega}+\left<\alpha\boldsymbol{u}_{+}\cdot\boldsymbol{n}_{+},\boldsymbol{v}_{+}\cdot\boldsymbol{n}_{-}\right>_{\Gamma}=0$.
In the following, we use
$\left<\alpha\boldsymbol{u}\cdot\boldsymbol{n},\boldsymbol{v}\cdot\boldsymbol{n}\right>_{\Gamma}$
to denote
$\left<\alpha\boldsymbol{u}_{+}\cdot\boldsymbol{n}_{+},\boldsymbol{v}_{+}\cdot\boldsymbol{n}_{+}\right>_{\Gamma}$.
From this and an immediate variational form of the second equation in (2), we
have a variational problem to find
$({\boldsymbol{u}},p)\in\boldsymbol{V}\times Q$ such that

(8a)
$\displaystyle\left(\boldsymbol{u},\boldsymbol{v}\right)_{\Omega}+\left<\alpha\boldsymbol{u}\cdot\boldsymbol{n},\boldsymbol{v}\cdot\boldsymbol{n}\right>_{\Gamma}-\left(p,\operatorname{div}\boldsymbol{v}\right)_{\Omega}$
$\displaystyle=0$ $\displaystyle\forall\boldsymbol{v}\in\boldsymbol{V},$ (8b)
$\displaystyle\left(\operatorname{div}\boldsymbol{u},q\right)_{\Omega}$
$\displaystyle=\left(f,q\right)_{\Omega}$ $\displaystyle\forall q\in Q.$

The stability of this system with an inf-sup condition

(9) $\displaystyle\inf_{\boldsymbol{v}\in\boldsymbol{V}}\sup_{q\in
Q}\frac{\left(q,\operatorname{div}\boldsymbol{v}\right)_{\Omega}}{\|\boldsymbol{v}\|_{\boldsymbol{V}}\|q\|_{0}}\geq
C>0$

was studied in [20].

###  2.3. Discretization with finite elements

We introduce finite element spaces for discretization. In the rest of the
paper we assume that $k\geq 1$ is a fixed integer. For a $d$-dimensional
simplex $D\subset{\mathbb{R}}^{n}$ ($d=n,n-1,n-2$), $\mathcal{P}_{k}(D)$ is
the space of polynomials on $D$ of degree $\leq k$. Similarly,
$\mathcal{P}_{k}(D;\mathbb{R}^{d})$ is the space of $\mathbb{R}^{d}$-valued
polynomials of degree $\leq k$ on the $d$-dimensional simplex $D$. For given
$k\geq 1$ and $n$-dimensional simplex $T\in\mathcal{T}_{h}$ let us define

(10) $\displaystyle{\boldsymbol{V}}_{k-1}^{RTN}(T)$
$\displaystyle=\mathcal{P}_{k-1}(T;\mathbb{R}^{n})+\begin{pmatrix}x_{1}\\\
\vdots\\\ x_{n}\end{pmatrix}\mathcal{P}_{k-1}(T),$ (11)
$\displaystyle\boldsymbol{V}_{k}^{BDM}(T)$
$\displaystyle=\mathcal{P}_{k}(T;\mathbb{R}^{n}).$

The Raviart–Thomas element (for $n=2$, [26]) and the first kind of Nedelec
$H(\operatorname{div})$ element (for $n=3$, [24]) are defined by

(12) $\displaystyle{\boldsymbol{V}}_{h}^{RTN}$
$\displaystyle=\\{{\boldsymbol{v}}\in\boldsymbol{V}\,:\,{\boldsymbol{v}}|_{T}\in{\boldsymbol{V}}_{k-1}^{RTN}(T),\quad\forall
T\in\mathcal{T}_{h}\\}.$

The Brezzi–Douglas–Marini element (for $n=2$, [10]) and the second kind of
Nedelec $H(\operatorname{div})$ element (for $n=3$, [25]) are defined by

(13) $\displaystyle\boldsymbol{V}_{h}^{BDM}$
$\displaystyle=\\{{\boldsymbol{v}}\in\boldsymbol{V}\,:\,{\boldsymbol{v}}|_{T}\in\boldsymbol{V}_{k}^{BDM}(T),\quad\forall
T\in\mathcal{T}_{h}\\}.$

The finite element spaces $Q_{h}^{0}$, $Q_{h}$, $Q_{h}^{*}$ are defined by

$\displaystyle Q_{h}^{0}$ $\displaystyle=\\{q\in
Q\,:\,q|_{T}\in\mathcal{P}_{0}(T)\quad\forall T\in\mathcal{T}_{h}\\},$
$\displaystyle Q_{h}$ $\displaystyle=\\{q\in
Q\,:\,q|_{T}\in\mathcal{P}_{k-1}(T)\quad\forall T\in\mathcal{T}_{h}\\},$ (14)
$\displaystyle Q_{h}^{*}$ $\displaystyle=\\{q\in
Q\,:\,q|_{T}\in\mathcal{P}_{k+1}(T)\quad\forall T\in\mathcal{T}_{h}\\},$

and $P_{h},P_{h}^{0}$ are the $L^{2}$ orthogonal projections to $Q_{h}$,
$Q_{h}^{0}$. For face $F$ or edge $E$ with integer $m\geq 0$, $P_{F}^{m}$ and
$P_{E}^{m}$ are the $L^{2}$ orthogonal projections to $\mathcal{P}_{m}(F)$ and
$\mathcal{P}_{m}(E)$.

It is well-known that $\operatorname{div}\boldsymbol{V}_{h}=Q_{h}$ for
$\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$ or
$\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{BDM}$ (cf. [9]). In addition, it is
proved in [20] with the assumption (1) that the pair
$(\boldsymbol{V}_{h},Q_{h})$ satisfies

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

with $C>0$ independent of mesh sizes. In the rest of this paper our
discussions are common for $\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$ or
$\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{BDM}$ unless we specify
$\boldsymbol{V}_{h}$ in our statements.

The discrete problem of (8) with $\boldsymbol{V}_{h}\times Q_{h}$ is to find
$(\boldsymbol{u}_{h},p_{h})\in\boldsymbol{V}_{h}\times Q_{h}$ such that

(16a)
$\displaystyle\left(\boldsymbol{u}_{h},\boldsymbol{v}\right)_{\Omega}+\left<\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n},\boldsymbol{v}\cdot\boldsymbol{n}\right>_{\Gamma}-\left(p_{h},\operatorname{div}\boldsymbol{v}\right)_{\Omega}$
$\displaystyle=0$ $\displaystyle\forall\boldsymbol{v}\in\boldsymbol{V}_{h},$
(16b)
$\displaystyle\left(\operatorname{div}\boldsymbol{u}_{h},q\right)_{\Omega}$
$\displaystyle=\left(f,q\right)_{\Omega}$ $\displaystyle\forall q\in Q_{h}.$

An a priori error analysis for $(\boldsymbol{u}_{h},p_{h})$ is proved in [20].

We finish this section by introducing a post-processed numerical solution of
$p$, a variant of the post-processing in [29]. Suppose that
$(\boldsymbol{u}_{h},p_{h})\in\boldsymbol{V}_{h}\times Q_{h}$ is a solution of
(16). Following the idea in [17, 14], a post-processed solution $p_{h}^{*}\in
Q_{h}^{*}$ is defined by

(17a) $\displaystyle\left(\nabla p_{h}^{*},\nabla q\right)_{T}$
$\displaystyle=-\left({\boldsymbol{u}}_{h},\nabla q\right)_{T}$
$\displaystyle\forall q\in Q_{h}^{*}(T),$ (17b)
$\displaystyle\left(p_{h}^{*},q\right)_{T}$
$\displaystyle=\left(P_{h}^{0}p_{h},q\right)_{T}$ $\displaystyle\forall q\in
Q_{h}^{0}(T)$

for $T\in\mathcal{T}_{h}$.

Here we remark that $Q_{h}^{*}$ in (14) contains all piecewise quadratic
polynomials because $k\geq 1$. This will be used in the proof of Lemma 3.6
which allows an elegant local efficiency proof of a posteriori error
estimator.

##  3\. A posteriori error estimate

In this section we define a posteriori error estimator and prove its
reliability and local efficiency. Let $(\boldsymbol{u}_{h},p_{h})$ be a
solution of (16) and $p_{h}^{*}$ be the post-processed pressure defined in
(17). Let

(18) $\displaystyle m=k-1\quad\text{ if
}\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}\quad\text{ and }\quad
m=k\quad\text{ if }\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{BDM}.$

Then, a posteriori error estimator $\eta$ for $n=3$ is defined by

(19)
$\displaystyle\eta=\left(\sum_{T\in\mathcal{T}_{h}}\eta_{T}^{2}+\sum_{F\in\mathcal{F}_{h}}\eta_{F}^{2}\right)^{\frac{1}{2}}$

where

(20) $\displaystyle\eta_{T}$ $\displaystyle:=\|\boldsymbol{u}_{h}+\nabla
p_{h}^{*}\|_{0,T},$ (21) $\displaystyle\eta_{F}$
$\displaystyle:=\begin{cases}h_{F}^{-1/2}\|\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\|_{0,F}&\text{
if }F\in\mathcal{F}_{h}^{0}\\\
\alpha^{-1/2}\|(I-P_{F}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\|_{0,F}&\text{
if }F\in\mathcal{F}_{h}^{\Gamma}\end{cases}.$

If $n=2$, $\eta$ is similarly defined by replacing $\eta_{F}$ by $\eta_{E}$
for edges in $\mathcal{E}_{h}^{0}$ and $\mathcal{E}_{h}^{\Gamma}$ with the
same formula in (21), so we omit the detailed definition.

###  3.1. Reliability estimate

We prove the reliability of the a posteriori error estimator $\eta$ in (19).
The main result is the following.

######  Theorem 3.1.

Suppose that $({\boldsymbol{u}},p)$, $(\boldsymbol{u}_{h},p_{h})$ are
solutions of (8), (16), and $\eta$ is defined by (19). Then, there exists
$C>0$ independent of mesh sizes and $\alpha$ in (5) such that

(22)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}-\boldsymbol{u}_{h}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}\leq
C\eta+\frac{1}{\pi}\left(\sum_{T\in\mathcal{T}_{h}}\operatorname{osc}(f,T)^{2}\right)^{\frac{1}{2}},\quad\operatorname{osc}(f,T):={h_{T}}\|f-P_{h}f\|_{0,T}.$

To prove this theorem, we split $\boldsymbol{u}-\boldsymbol{u}_{h}$ into two
components which are orthogonal with an inner product given by the bilinear
form of $\boldsymbol{u}$ and $\boldsymbol{v}$ in (8a). In the lemma below, we
first show that there is an orthogonal decomposition of
$\boldsymbol{u}-\boldsymbol{u}_{h}$.

######  Lemma 3.1.

Suppose that $(\tilde{{\boldsymbol{u}}},\tilde{p})$ is the solution of

(23a)
$\displaystyle\left(\tilde{\boldsymbol{u}},\boldsymbol{v}\right)_{\Omega}-\left(\tilde{p},\operatorname{div}\boldsymbol{v}\right)_{\Omega}+\left<\alpha\tilde{\boldsymbol{u}}\cdot\boldsymbol{n},\boldsymbol{v}\cdot\boldsymbol{n}\right>_{\Gamma}$
$\displaystyle=0$ $\displaystyle\forall\boldsymbol{v}\in\boldsymbol{V},$ (23b)
$\displaystyle\left(\operatorname{div}\tilde{\boldsymbol{u}},q\right)_{\Omega}$
$\displaystyle=\left(P_{h}f,q\right)_{\Omega}$ $\displaystyle\forall q\in Q.$

Then,

(24)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}-\boldsymbol{u}_{h}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}=\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}+\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\tilde{\boldsymbol{u}}-\boldsymbol{u}_{h}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}.$

###### Proof.

Then, it is easy to see that

(25a)
$\displaystyle\left(\boldsymbol{u}-\tilde{\boldsymbol{u}},\boldsymbol{v}\right)_{\Omega}-\left(p-\tilde{p},\operatorname{div}\boldsymbol{v}\right)_{\Omega}+\left<\alpha(\boldsymbol{u}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n},\boldsymbol{v}\cdot\boldsymbol{n}\right>_{\Gamma}$
$\displaystyle=0,$ (25b)
$\displaystyle\left(\operatorname{div}(\boldsymbol{u}-\tilde{\boldsymbol{u}}),q\right)_{\Omega}$
$\displaystyle=\left(f-P_{h}f,q\right)_{\Omega}$

for all $\boldsymbol{v}\in\boldsymbol{V}$ and $q\in Q$. If
$\boldsymbol{v}=\tilde{\boldsymbol{u}}-\boldsymbol{u}_{h}\in\boldsymbol{V}$,
then

$\displaystyle\left(\boldsymbol{u}-\tilde{\boldsymbol{u}},\tilde{\boldsymbol{u}}-\boldsymbol{u}_{h}\right)_{\Omega}+\left<\alpha(\boldsymbol{u}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n},(\tilde{\boldsymbol{u}}-\boldsymbol{u}_{h})\cdot\boldsymbol{n}\right>_{\Gamma}=0$

because $\operatorname{div}(\tilde{\boldsymbol{u}}-\boldsymbol{u}_{h})=0$.
Then, (24) follows from this orthogonality. ∎

As a consequence of (24), it suffices to estimate
$\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
and
$\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
by the right-hand side terms in (22). We first estimate
$\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}$.

######  Lemma 3.2.

Suppose that $(\tilde{{\boldsymbol{u}}},\tilde{p})$ is defined as in Lemma
3.1. Then,

(26)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}\leq\frac{1}{\pi}\left(\sum_{T\in\mathcal{T}_{h}}\operatorname{osc}(f,T)^{2}\right)^{\frac{1}{2}}.$

###### Proof.

Taking arbitrary test function vanishing near $\Gamma$ in (25a), we obtain

$\displaystyle\boldsymbol{u}-\tilde{\boldsymbol{u}}=\nabla(p-\tilde{p})\quad\text{
in }{L^{2}(\Omega_{0};\mathbb{R}^{n})}$

for every open set $\Omega_{0}$ such that
$\overline{\Omega_{0}}\subset\Omega\setminus\Gamma$. By the dominated
convergence theorem,
$\boldsymbol{u}-\tilde{\boldsymbol{u}}=\nabla(p-\tilde{p})$ on every
$T\in\mathcal{T}_{h}$. If
$\boldsymbol{v}=\boldsymbol{u}-\tilde{\boldsymbol{u}}$ in (23a), then

(27)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}=\left(p-\tilde{p},f-P_{h}f\right)_{\Omega}$

by (6) and (25b). By the Cauchy–Schwarz and element-wise Poincare inequalities
for mean-value zero functions (cf. [6]),

$\displaystyle\left(p-\tilde{p},f-P_{h}f\right)_{\Omega}$
$\displaystyle=\left(p-\tilde{p}-P_{h}(p-\tilde{p}),f-P_{h}f\right)_{\Omega}$
$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}\frac{h_{T}}{\pi}\|\nabla(p-\tilde{p})\|_{0,T}\|f-P_{h}f\|_{0,T}$
$\displaystyle\leq\frac{1}{\pi}\|\nabla(p-\tilde{p})\|_{0}\left(\sum_{T\in\mathcal{T}_{h}}{h_{T}^{2}}\|f-P_{h}f\|_{0,T}^{2}\right)^{\frac{1}{2}}$
(28)
$\displaystyle\leq\frac{1}{\pi}\|\boldsymbol{u}-\tilde{\boldsymbol{u}}\|_{0}\left(\sum_{T\in\mathcal{T}_{h}}\operatorname{osc}(f,T)^{2}\right)^{\frac{1}{2}}.$

Combining (27) and (28), we can obtain (26). ∎

To estimate
$\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}$
by the a posteriori error estimator $\eta$, we need an auxiliary finite
element space $S_{h}$. We choose different $S_{h}$ for
$\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$ and
$\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{BDM}$ and for $n=2,3$.

We first define $S_{h}$ for $\boldsymbol{V}_{h}^{BDM}$. Note that
$\boldsymbol{V}_{h}^{BDM}=\mathcal{P}_{k}\Lambda^{n-1}(\mathcal{T}_{h})$ in
the language of the finite element exterior calculus ([3, 4]). If
$\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{BDM}$, then

(29) $\displaystyle S_{h}=\mathcal{P}_{k+1}\Lambda^{n-2}(\mathcal{T}_{h}),$

which is the Lagrange finite element of degree $k+1$ if $n=2$ and is the
Nedelec edge element of the 2nd kind ([25]) with degree $k+1$ if $n=3$.

If $\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$ and $n=2$, then

(30) $\displaystyle S_{h}=\mathcal{P}_{k}\Lambda^{0}(\mathcal{T}_{h}),$

the Lagrange finite element of degree $k$.

If $\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$ and $n=3$, then we define a
new finite element space $S_{h}$ obtained by enriching
$\mathcal{P}_{k}\Lambda^{1}(\mathcal{T}_{h})$ with curl-free edge bubble
functions which will be described below.

For an edge $E$ in a triangulation $\mathcal{T}_{h}$ for $n=3$, define $M_{E}$
by

(31) $\displaystyle M_{E}=\bigcup_{T\in\mathcal{T}_{h},E\subset\partial T}T.$

Note that each face $F\subset\partial M_{E}$ does not contain $E$. Denoting
the barycentric coordinate which vanishes on $F$ by $\lambda_{F}$, we define
$b_{E}$ by

(32) $\displaystyle b_{E}=\prod_{F\subset\partial M_{E}}\lambda_{F}.$

Since every tetrahedron in $M_{E}$ does not have more than two distinct faces
which do not contain $E$, $b_{E}|_{T}\in\mathcal{P}_{2}(T)$ for every
tetrahedron $T\subset M_{E}$. In the discussion below,
$\mathcal{P}_{k+1}^{\perp}(E)$ is the space of polynomials with degree $(k+1)$
on $E$ which are orthogonal to all polynomials with degree $k$ on $E$, and
$\bar{\mathcal{P}}_{k+1}^{\perp}(E)$ is the space of polynomials on $M_{E}$
which are constant on every plane perpendicular to $E$ and the restriction of
the polynomials on $E$ are in $\mathcal{P}_{k+1}^{\perp}(E)$.

In the following lemma, for an edge $E$, $\boldsymbol{t}_{E}$ is a unit
tangential vector of $E$ and $\frac{\partial}{\partial\boldsymbol{t}_{E}}$ is
the derivative along the direction of $\boldsymbol{t}_{E}$.

######  Lemma 3.3.

For an edge $E$,

(33)
$\displaystyle\mathcal{B}(E)=\operatorname{span}\left\\{\nabla\left(\frac{\partial
q}{\partial\boldsymbol{t}_{E}}b_{E}\right)\,:\,q\in\bar{\mathcal{P}}_{k+1}^{\perp}(E)\right\\}.$

For a tetrahedron $T$ and a fixed face $\tilde{F}\subset\partial T$ let
$S(T,\tilde{F})$ be

$\displaystyle
S(T,\tilde{F})=\mathcal{P}_{k}(T;\mathbb{R}^{3})+\oplus_{E\subset\partial\tilde{F}}\mathcal{B}(E).$

and a set of local degrees of freedom for $\tau\in S(T,\tilde{F})$ is given by

(34) $\displaystyle\tau$
$\displaystyle\longmapsto\int_{E}\tau\cdot\boldsymbol{t}_{E}\tilde{q}\,dl$
$\displaystyle\forall\tilde{q}\in\mathcal{P}_{k}(E)\text{ if
}E\not\subset\partial\tilde{F},$ (35) $\displaystyle\tau$
$\displaystyle\longmapsto\int_{E}\tau\cdot\boldsymbol{t}_{E}\tilde{q}\,dl$
$\displaystyle\forall\tilde{q}\in\mathcal{P}_{k+1}(E)\text{ if
}E\subset\partial\tilde{F},$ (36) $\displaystyle\tau$
$\displaystyle\longmapsto\int_{F}(\tau\times\boldsymbol{n}_{F})\cdot\boldsymbol{q}\,ds$
$\displaystyle\forall\boldsymbol{q}\in\boldsymbol{V}_{k-1}^{RTN}(F)\text{ if
}k\geq 2,$ (37) $\displaystyle\tau$
$\displaystyle\longmapsto\int_{T}\tau\cdot\xi\,dx$
$\displaystyle\forall\xi\in\boldsymbol{V}_{k-2}^{RTN}(T)\text{ if }k\geq 3.$

Then, $\tau\in S(T,\tilde{F})$ is uniquely determined by (34), (35), (36),
(37).

###### Proof.

Suppose that $\tau\in S(T,\tilde{F})$ and all DOFs of $\tau$ given by (34),
(35), (36), (37) vanish. Let $E_{i}$, $i=1,2,3$ be the edges of $\tilde{F}$.
We can write $\tau$ as $\tau=\sum_{i=0}^{3}\tau_{i}$ with
$\tau_{0}\in\mathcal{P}_{k}(T;\mathbb{R}^{3})$ and

(38) $\displaystyle\tau_{i}=\nabla\left(\frac{\partial
q_{i}}{\partial\boldsymbol{t}_{E_{i}}}b_{E_{i}}\right)\in\mathcal{B}(E_{i}),\quad
q_{i}\in\bar{\mathcal{P}}_{k+1}^{\perp}(E_{i}),i=1,2,3.$

By the vanishing DOFs assumption, for $E_{i}\subset\partial\tilde{F}$,

$\displaystyle\int_{E_{i}}\tau\cdot\boldsymbol{t}_{E_{i}}\tilde{q}\,dl=\int_{E_{i}}\tau_{0}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}\,dl+\sum_{j=1}^{3}\int_{E_{i}}\tau_{j}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}\,dl=0$

for all $\tilde{q}\in\mathcal{P}_{k+1}(E_{i})$. By the definition (33),
$\tau_{i}\cdot\boldsymbol{t}_{E_{j}}|_{E_{j}}=0$ if $j\not=i$, so we get

(39)
$\displaystyle\int_{E_{i}}\tau_{0}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}\,dl+\int_{E_{i}}\tau_{i}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}\,dl=0.$

Consider the decomposition
$\tilde{q}=\tilde{q}_{0}+\tilde{q}_{1}\in\mathcal{P}_{k}(E_{i})\oplus\mathcal{P}_{k+1}^{\perp}(E_{i})$.
Then,

(40)
$\displaystyle\int_{E_{i}}\tau_{0}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}_{1}\,dl=0.$

Taking the integration by parts twice gives

$\displaystyle\int_{E_{i}}\tau_{i}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}_{0}\,dl$
$\displaystyle=\int_{E_{i}}\frac{\partial}{\partial\boldsymbol{t}_{E_{i}}}\left(\frac{\partial
q_{i}}{\partial\boldsymbol{t}_{E_{i}}}b_{E_{i}}\right)\tilde{q}_{0}\,dl$ (41)
$\displaystyle=-\int_{E_{i}}\left(\frac{\partial
q_{i}}{\partial\boldsymbol{t}_{E_{i}}}b_{E_{i}}\right)\frac{\partial\tilde{q}_{0}}{\partial\boldsymbol{t}_{E_{i}}}\,dl$
$\displaystyle=\int_{E_{i}}q_{i}\frac{\partial}{\partial\boldsymbol{t}_{E_{i}}}\left(b_{E_{i}}\frac{\partial\tilde{q}_{0}}{\partial\boldsymbol{t}_{E_{i}}}\right)\,dl$
$\displaystyle=0$

where the last identity follows from
$\tilde{q}_{0}\in\mathcal{P}_{k}({E_{i}})$,
$b_{E_{i}}|_{E_{i}}\in\mathcal{P}_{2}({E_{i}})$, and
$q_{i}|_{E_{i}}\in\mathcal{P}_{k+1}^{\perp}({E_{i}})$. Therefore, (39) is
reduced to

$\displaystyle\int_{E_{i}}\tau_{0}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}_{0}\,dl+\int_{E_{i}}\tau_{i}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}_{1}\,dl=0.$

If $\tilde{q}_{0}=\tau_{0}\cdot\boldsymbol{t}_{E_{i}}$, $\tilde{q}_{1}=-q_{i}$
in this formula, and use the identity

$\displaystyle\int_{E_{i}}\tau_{i}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}_{1}\,dl$
$\displaystyle=-\int_{E_{i}}\frac{\partial}{\partial\boldsymbol{t}_{E_{i}}}\left(\frac{\partial
q_{i}}{\partial\boldsymbol{t}_{E_{i}}}b_{E_{i}}\right){q}_{i}\,dl=\int_{E_{i}}\left(\frac{\partial
q_{i}}{\partial\boldsymbol{t}_{E_{i}}}\right)^{2}b_{E_{i}}\,dl,$

then

$\displaystyle\int_{E_{i}}\tau_{0}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}_{0}\,dl+\int_{E_{i}}\tau_{i}\cdot\boldsymbol{t}_{E_{i}}\tilde{q}_{1}\,dl$
$\displaystyle=\int_{E_{i}}(\tau_{0}\cdot\boldsymbol{t}_{E_{i}})^{2}\,dl+\int_{E_{i}}\left(\frac{\partial
q_{i}}{\partial\boldsymbol{t}_{E_{i}}}\right)^{2}b_{E_{i}}\,dl$
$\displaystyle=0.$

Since $b_{E_{i}}>0$ on $E_{i}$, $\frac{\partial
q_{i}}{\partial\boldsymbol{t}_{E_{i}}}|_{E_{i}}=0$, so $q_{i}|_{E_{i}}$ is
constant. Furthermore, $q_{i}|_{E_{i}}\in\mathcal{P}_{k+1}^{\perp}(E_{i})$, so
$q_{i}=0$. As a consequence, $\tau_{i}=0$ for $i=1,2,3$. Then, $\tau=0$
follows by a standard unisolvency proof of the Nedelec edge elements of the
2nd kind. ∎

We now define $S_{h}$ for $n=3$ and
$\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$. The enriched
$H(\operatorname{curl})$ element $S_{h}$ with the shape functions

(42) $\displaystyle
S_{h}=\mathcal{P}_{k}(\mathcal{T}_{h};\mathbb{R}^{3})+\bigoplus_{E\in\mathcal{E}_{h}^{\Gamma}}\mathcal{B}(E)$

and the global degrees of freedom (34) for
$E\in\mathcal{E}_{h}\setminus\mathcal{E}_{h}^{\Gamma}$, (35) for
$E\in\mathcal{E}_{h}^{\Gamma}$, (36) for $F\in\mathcal{F}_{h}$, (37) for
$E\in\mathcal{T}_{h}$.

For $S_{h}$ defined by (29), (30), (42) depending on $n$ and
$\boldsymbol{V}_{h}$, we can check that

(43) $\displaystyle\operatorname{curl}S_{h}\subset\boldsymbol{V}_{h}.$

Here we show existence of an appropriate interpolation.

######  Lemma 3.4.

Suppose that $\Psi\in H^{1}(\Omega)$ and
$\operatorname{curl}\Psi\cdot\boldsymbol{n}|_{\Gamma}\in L^{2}(\Gamma)$. Let
$S_{h}$ be defined by (29), (30), (42) depending on $n$ and
$\boldsymbol{V}_{h}$, and recall $l$ defined in (18). Then, there exists
$I_{h}\Psi\in S_{h}$ such that

(44)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\operatorname{curl}I_{h}\Psi\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}$
$\displaystyle\leq
C(\|\Psi\|_{1}+\|\alpha^{1/2}\operatorname{curl}\Psi\cdot\boldsymbol{n}\|_{0,\Gamma}),$
(45) $\displaystyle\|\Psi-I_{h}\Psi\|_{0}$ $\displaystyle\leq
C(\|\Psi\|_{1}+\|\alpha\operatorname{curl}\Psi\cdot\boldsymbol{n}\|_{0,\Gamma}),$
(46) $\displaystyle\operatorname{curl}I_{h}\Psi\cdot\boldsymbol{n}|_{E}$
$\displaystyle=P_{E}^{m}(\operatorname{curl}\Psi\cdot\boldsymbol{n}|_{E}),\quad\forall
E\in\mathcal{E}_{h}^{\Gamma},\text{ if }n=2,$ (47)
$\displaystyle\operatorname{curl}I_{h}\Psi\cdot\boldsymbol{n}|_{F}$
$\displaystyle=P_{F}^{m}(\operatorname{curl}\Psi\cdot\boldsymbol{n}|_{F}),\quad\forall
F\in\mathcal{F}_{h}^{\Gamma},\text{ if }n=3,$

and for $E\in\mathcal{E}_{h}^{\Gamma}$,

(48) $\displaystyle\begin{cases}(I_{h}\Psi-\Psi)(v)&=0\quad\text{ if
}n=2,v\in\partial E,\\\
\int_{E}(I_{h}\Psi-\Psi)\cdot\boldsymbol{t}_{E}\tilde{q}\,dl&=0\quad\text{ if
}n=3,\tilde{q}\in\mathcal{P}_{k+1}(E).\end{cases}$

###### Proof.

Suppose that $n=2$. In this case,
$\operatorname{curl}\Psi\cdot\boldsymbol{n}_{\Gamma}\in L^{2}(\Gamma)$ means
that the tangential derivative of $\Psi$ along $\Gamma$ is in $L^{2}(\Gamma)$
(cf. (95) in Appendix). Since $\Psi|_{\Gamma}\in L^{2}(\Gamma)$, we have
$\Psi|_{\Gamma}\in H^{1}(\Gamma)$. Let $\mathcal{N}_{h}$ be the set of vertex
nodes in $\overline{\Omega}$ which determines the degrees of freedom of
$S_{h}$, a Lagrange finite element. By the Sobolev embedding on the
1-dimensional submanifold $\Gamma$, vertex evaluation of $\Psi$ on the nodes
on $\Gamma$ is well-defined. On $E\in\mathcal{E}_{h}^{\Gamma}$, we define
$I_{h}^{\Gamma}\Psi|_{E}$ by

(49) $\displaystyle I_{h}^{\Gamma}\Psi(v)$ $\displaystyle=\Psi(v)$
$\displaystyle v\in\partial E,$ (50)
$\displaystyle\int_{E}I_{h}^{\Gamma}\Psi\frac{\partial}{\partial\boldsymbol{t}_{E}}\tilde{q}\,dl$
$\displaystyle=\int_{E}\Psi\frac{\partial}{\partial\boldsymbol{t}_{E}}\tilde{q}\,dl$
$\displaystyle\forall\tilde{q}\in\mathcal{P}_{m}(E),\tilde{q}\perp\mathcal{P}_{0}(E)$

for $m$ defined in (18), and $I_{h}^{\Gamma}\Psi(v)=0$ for
$v\in\mathcal{N}_{h}\setminus\Gamma$. Then, by a standard scaling argument,

$\displaystyle\|I_{h}^{\Gamma}\Psi\|_{0}+\|\operatorname{curl}I_{h}^{\Gamma}\Psi\|_{0}\leq
C(\|\Psi\|_{1}+\|\operatorname{curl}\Psi\cdot\boldsymbol{n}\|_{0,\Gamma})$

with $C>0$ independent of mesh sizes. Let $I_{h}^{SZ}$ be a Scott–Zhang
interpolation (cf. [28]) which takes $\Gamma$ as a vanishing interface and
satisfies

$\displaystyle\|I_{h}^{SZ}\Phi\|_{0}+\|\operatorname{curl}I_{h}^{SZ}\Phi\|_{0}\leq
C\|\Phi\|_{1}$

for $\Psi\in H^{1}(\Omega)$. If we define $I_{h}\Psi$ by

(51) $\displaystyle I_{h}\Psi=I_{h}^{SZ}(\Psi-
I_{h}^{\Gamma}\Psi)+I_{h}^{\Gamma}\Psi,$

then Therefore, $I_{h}\Psi$ is bounded by
$\|\Psi\|_{1}+\|\operatorname{curl}\Psi\cdot\boldsymbol{n}\|_{0,\Gamma}$.

We now check (44), (45), (46), (48). First, (48) is a consequence of (49) and
(51). By (49), (50),
$\operatorname{curl}I_{h}\Psi\cdot\boldsymbol{n}|_{\Gamma}=I_{h}^{\Gamma}\Psi\cdot\boldsymbol{n}|_{\Gamma}$,
and the integration by parts on every $E\in\mathcal{E}_{h}^{\Gamma}$, (46)
follows.

Furthermore, if $\Psi\in S_{h}$, then $I_{h}^{\Gamma}\Psi=\Psi|_{\Gamma}$, so
$I_{h}$ is the identity map on $S_{h}$ because $I_{h}^{SZ}$ is the identity
map for the elements in $S_{h}$ which vanish on $\Gamma$. By the
Bramble–Hilbert lemma, $\|\Psi-I_{h}\Psi\|_{0}\leq
Ch(\|\Psi\|_{1}+\|\operatorname{curl}\Psi\cdot\boldsymbol{n}\|_{0,\Gamma})$.

Suppose that $n=3$. First, $\Psi\times\boldsymbol{n}_{\Gamma}\in
L^{2}(\Gamma;\mathbb{R}^{2})$ and
$\operatorname{curl}\Psi\cdot\boldsymbol{n}_{\Gamma}\in L^{2}(\Gamma)$ impliy
that the tangential component of $\Psi$ on $\Gamma$ is in the rotated
$H(\operatorname{div})$ space on $\Gamma$ (cf. (96) in Appendix). Since
$\Psi|_{\Gamma}\in H^{s}(\Gamma)$ with $s>0$ as a trace of
$H^{1}(\Omega;\mathbb{R}^{3})$, $\Psi\times\boldsymbol{n}_{\Gamma}\in
L^{r}(\Gamma)$ for $r>2$ by Sobolev embedding, so

(52)
$\displaystyle\int_{E}I_{h}^{\Gamma}\Psi\cdot\boldsymbol{t}_{E}\tilde{q}\,dl$
$\displaystyle=\int_{E}\Psi\cdot\boldsymbol{t}_{E}\tilde{q}\,dl,$
$\displaystyle E\in\mathcal{E}_{h}^{\Gamma},\tilde{q}\in\mathcal{P}_{k+1}(E),$
(53)
$\displaystyle\int_{F}I_{h}^{\Gamma}\Psi\times\boldsymbol{n}_{F}\cdot\xi\,ds$
$\displaystyle=\int_{F}(\Psi\times\boldsymbol{n}_{F})\cdot\xi\,ds,$
$\displaystyle
F\in\mathcal{F}_{h}^{\Gamma},\xi\in\mathcal{P}_{m-2}(F;\mathbb{R}^{2})$

are well-defined (cf. [9]). $I_{h}\Psi\in S_{h}$ is defined by taking zeros
for all other degrees of freedom, and $\|I_{h}^{\Gamma}\Psi\|_{0}$ is bounded
by $\|\Psi\|_{1}+\|\operatorname{curl}\Psi\cdot\boldsymbol{n}\|_{0,\Gamma}$.
There exists a Scott–Zhang type interpolation $I_{h}^{SZ}$ for
$H(\operatorname{curl})$ elements (see [15]) with vanishing interface
$\Gamma$, so define $I_{h}$ as in (51). By an argument similar to the proof
for $n=2$, (45) can be obtained, and (48) follows from (52) and (51). Finally,
(47) follows from (52), (53), and the integration by parts on every face
$F\in\mathcal{F}_{h}^{\Gamma}$ (see (97) in Appendix for details). ∎

We now prove a reliability estimate of
$\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}$.

######  Theorem 3.2.

Suppose that $(\tilde{{\boldsymbol{u}}},\tilde{p})$ is defined as in Lemma
3.1. Then, there exists $C>0$ independent of mesh sizes such that

(54)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}\leq
C\eta.$

###### Proof.

We now estimate
$\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$.
First, since $\Omega$ is homologically trivial, there exists $\Psi\in
H^{1}(\Omega,\mathbb{R}^{n})$ such that
$\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}=\operatorname{curl}\Psi$ and

$\displaystyle\|\Psi\|_{1}$ $\displaystyle\leq
C\|\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\|_{0},$
$\displaystyle\left<\alpha\operatorname{curl}\Psi\cdot\boldsymbol{n},\operatorname{curl}\Psi\cdot\boldsymbol{n}\right>_{\Gamma}$
$\displaystyle=\left<\alpha(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n},(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n}\right>_{\Gamma}$

because $\operatorname{div}(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}})=0$.
Thus,

(55)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}=\left(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}},\operatorname{curl}\Psi\right)_{\Omega}+\left<\alpha(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n},\operatorname{curl}\Psi\cdot\boldsymbol{n}\right>_{\Gamma}.$

Since

(56a)
$\displaystyle\left(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}},\boldsymbol{v}\right)_{\Omega}-\left(p_{h}-\tilde{p},\operatorname{div}\boldsymbol{v}\right)_{\Omega}+\left<\alpha(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n},\boldsymbol{v}\cdot\boldsymbol{n}\right>_{\Gamma}$
$\displaystyle=0$

for all $\boldsymbol{v}\in\boldsymbol{V}_{h}$, we have

(57)
$\displaystyle\left(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}},\operatorname{curl}I_{h}\Psi\right)_{\Omega}+\left<\alpha(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n},\operatorname{curl}I_{h}\Psi\cdot\boldsymbol{n}\right>_{\Gamma}=0$

for $I_{h}$ in Lemma 3.4 because of $I_{h}\Psi\in S_{h}$ and (43). Applying
(57) to (55),

(58)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
$\displaystyle=\left(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}},\operatorname{curl}(\Psi-
I_{h}\Psi)\right)_{\Omega}$
$\displaystyle\quad+\left<\alpha(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n},\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{\Gamma}.$

Since $\tilde{\boldsymbol{u}}=-\nabla\tilde{p}$ and
$\alpha\tilde{\boldsymbol{u}}\cdot\boldsymbol{n}=\left[\\!\left[{\tilde{p}}\right]\\!\right]$
on $\Gamma$, we can further obtain

(59)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
$\displaystyle=\left(\boldsymbol{u}_{h},\operatorname{curl}(\Psi-
I_{h}\Psi)\right)_{\Omega}$
$\displaystyle\quad+\left<\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n},\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{\Gamma}$

by the integration by parts. For $m$ defined in (18),
$(\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n})|_{E}\in\mathcal{P}_{m}(E)$ for
$E\in\mathcal{E}_{h}^{\Gamma}$ if $n=2$ and
$(\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n})|_{F}\in\mathcal{P}_{m}(F)$ for
$F\in\mathcal{F}_{h}^{\Gamma}$ if $n=3$. Therefore, (59) is reduced to

(60)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
$\displaystyle=\left(\boldsymbol{u}_{h},\operatorname{curl}(\Psi-
I_{h}\Psi)\right)_{\Omega}$

by (46) and (47).

If $n=2$, a simple algebra and triangle-wise integration by parts give

$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
$\displaystyle=\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\operatorname{curl}(\Psi-I_{h}\Psi)\right)_{\Omega}-\left(\nabla
p_{h}^{*},\operatorname{curl}(\Psi-I_{h}\Psi)\right)_{\Omega}$ (61)
$\displaystyle=\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\operatorname{curl}(\Psi-
I_{h}\Psi)\right)_{\Omega}-\sum_{T\in\mathcal{T}_{h}}\left<\nabla
p_{h}^{*}\cdot\boldsymbol{t}_{\partial T},\Psi-I_{h}\Psi\right>_{\partial T}.$

By edge-wise integration by parts using $(\Psi-I_{h}\Psi)(v)=0$ for every
endpoint $v$ of edges $E\in\mathcal{E}_{h}^{\Gamma}$,

$\displaystyle\sum_{T\in\mathcal{T}_{h}}\left<\nabla
p_{h}^{*}\cdot\boldsymbol{t}_{\partial T},\Psi-I_{h}\Psi\right>_{\partial T}$
$\displaystyle\quad=\pm\sum_{E\in\mathcal{E}_{h}}\left<\nabla\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\cdot\boldsymbol{t}_{E},\Psi-
I_{h}\Psi\right>_{E}$ (62)
$\displaystyle\quad=\pm\sum_{E\in\mathcal{E}_{h}^{\Gamma}}\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{E}\pm\sum_{E\in\mathcal{E}_{h}\setminus\mathcal{E}_{h}^{\Gamma}}\left<\nabla\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\cdot\boldsymbol{t}_{E},\Psi-
I_{h}\Psi\right>_{E}$
$\displaystyle\quad=\pm\sum_{E\in\mathcal{E}_{h}^{\Gamma}}\left<(I-P_{E}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{E}$
$\displaystyle\qquad\pm\sum_{E\in\mathcal{E}_{h}\setminus\mathcal{E}_{h}^{\Gamma}}\left<\nabla\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\cdot\boldsymbol{t}_{E},\Psi-
I_{h}\Psi\right>_{E}.$

Here we use $\pm$ due to sign ambiguity of the definitions of
$\left[\\!\left[{p_{h}^{*}}\right]\\!\right]$ and $\boldsymbol{t}_{E}$.
However, we will use the Cauchy–Schwarz inequality to estimate the terms that
this ambiguous sign is involved, so the exact sign is not important in the
rest of discussions. By this and (61),

$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
$\displaystyle=\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\operatorname{curl}(\Psi-
I_{h}\Psi)\right)_{\Omega}\pm\sum_{E\in\mathcal{E}_{h}\setminus\mathcal{E}_{h}^{\Gamma}}\left<\nabla\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\cdot\boldsymbol{t}_{E},\Psi-
I_{h}\Psi\right>_{E}$ (63)
$\displaystyle\quad\pm\sum_{E\in\mathcal{E}_{h}^{\Gamma}}\left<(I-P_{E}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{E}$
$\displaystyle=:I_{2,a}+I_{2,b}+I_{2,c}.$

By the Cauchy–Schwarz inequality

$\displaystyle|I_{2,a}|$ $\displaystyle\leq\|\boldsymbol{u}_{h}+\nabla
p_{h}^{*}\|_{0}\|\operatorname{curl}(\Psi-I_{h}\Psi)\|_{0}$ (64)
$\displaystyle\leq
C\left(\sum_{T\in\mathcal{T}_{h}}\eta_{T}^{2}\right)^{\frac{1}{2}}\left(\|\Psi\|_{1}+\|\operatorname{curl}\Psi\cdot\boldsymbol{n}\|_{L^{2}(\Gamma)}\right)$
$\displaystyle\leq
C\left(\sum_{T\in\mathcal{T}_{h}}\eta_{T}^{2}\right)^{\frac{1}{2}}\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}.$

By element-wise inverse inequality and an approximation property of $\Psi-
I_{h}\Psi$,

$\displaystyle|I_{2,b}|$
$\displaystyle\leq\sum_{E\in\mathcal{E}_{h}^{0}\setminus\mathcal{E}_{h}^{\Gamma}}\left|\left<\nabla\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\cdot\boldsymbol{t}_{E},\Psi-
I_{h}\Psi\right>_{E}\right|$ (65) $\displaystyle\leq
C\left(\sum_{E\in\mathcal{E}_{h}\setminus\mathcal{E}_{h}^{\Gamma}}h_{E}^{-1}\|\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\|_{0,E}^{2}\right)^{\frac{1}{2}}\|\Psi\|_{1}$
$\displaystyle\leq
C\left(\sum_{E\in\mathcal{E}_{h}\setminus\mathcal{E}_{h}^{\Gamma}}\eta_{E}^{2}\right)^{\frac{1}{2}}\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}.$

For $I_{2,c}$,

$\displaystyle\left|\left<(I-P_{l}^{E})\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{E}\right|$ (66)
$\displaystyle\quad\leq\alpha^{-1/2}\|(I-P_{E}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\|_{0,E}\|\alpha^{1/2}\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\|_{0,E}$ $\displaystyle\quad\leq
2\alpha^{-1/2}\|(I-P_{E}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\|_{0,E}\|\alpha^{1/2}\operatorname{curl}\Psi\cdot\boldsymbol{n}\|_{0,E}.$

Combining (63), (64), (65), (66), we obtain

$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}\leq
C\eta.$

If $n=3$, then

$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
$\displaystyle=\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\operatorname{curl}(\Psi-I_{h}\Psi)\right)_{\Omega}-\left(\nabla
p_{h}^{*},\operatorname{curl}(\Psi-I_{h}\Psi)\right)_{\Omega}$ (67)
$\displaystyle=\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\operatorname{curl}(\Psi-
I_{h}\Psi)\right)_{\Omega}-\sum_{T\in\mathcal{T}_{h}}\left<\nabla
p_{h}^{*},\boldsymbol{n}\times(\Psi-I_{h}\Psi)\right>_{\partial T}$

by tetrahedron-wise integration by parts. By face-wise integration by parts
(97) and by the property (48),

$\displaystyle\sum_{T\in\mathcal{T}_{h}}\left<\nabla
p_{h}^{*},\boldsymbol{n}\times(\Psi-I_{h}\Psi)\right>_{\partial T\cap\Gamma}$
$\displaystyle=\pm\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\left<\nabla\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\boldsymbol{n}\times(\Psi-
I_{h}\Psi)\right>_{F}$ (68)
$\displaystyle=\pm\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{F}$
$\displaystyle=\pm\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\left<(I-P_{F}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{F}.$

We can proceed using this and (67) to obtain

$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}$
$\displaystyle=\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\operatorname{curl}(\Psi-
I_{h}\Psi)\right)_{\Omega}\pm\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}\left<\nabla\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\times\boldsymbol{n},(\Psi-
I_{h}\Psi)\right>_{F}$ (69)
$\displaystyle\quad\pm\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\left<(I-P_{F}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\operatorname{curl}(\Psi-
I_{h}\Psi)\cdot\boldsymbol{n}\right>_{F}$
$\displaystyle=:I_{3,a}+I_{3,b}+I_{3,c}.$

By the arguments which are completely similar to the ones for $n=2$, we can
obtain

(70) $\displaystyle|I_{3,a}|$ $\displaystyle\leq
C\left(\sum_{T\in\mathcal{T}_{h}}\|\boldsymbol{u}_{h}+\nabla
p_{h}^{*}\|_{0,T}^{2}\right)^{\frac{1}{2}}\|\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\|_{0},$
(71) $\displaystyle|I_{3,b}|$ $\displaystyle\leq
C\left(\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}h_{F}^{-1}\|\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\|_{0,F}^{2}\right)^{\frac{1}{2}}\|\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}}\|_{0},$
(72) $\displaystyle|I_{3,c}|$ $\displaystyle\leq
C\left(\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\alpha^{-1}\|(I-P_{F}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\|_{0,F}^{2}\right)^{\frac{1}{2}}\|\alpha^{1/2}(\boldsymbol{u}_{h}-\tilde{\boldsymbol{u}})\cdot\boldsymbol{n}\|_{0,\Gamma}.$

Applying these estimates to (69), we can obtain (54). ∎

######  Remark 3.5.

The new $S_{h}$ space in (42) for $n=3$,
$\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$, is necessary for (68). More
precisely, the new $S_{h}$ allows an interpolation $I_{h}$ satisfying (48) for
$n=3$, $\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$, which is necessary for
the first equality in (68).

###  3.2. Local Efficiency

In this subsection we show local efficiency of the a posteriori error
estimator. We give a detailed proof for $n=3$ because the two-dimensional case
is almost same.

We prove a lemma employing the techniques in [30, 17].

######  Lemma 3.6.

For $F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}$,

(73) $\displaystyle\int_{F}\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\,ds$
$\displaystyle=0\quad\text{ if
}F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma},$ (74)
$\displaystyle\int_{F}(\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n}-\left[\\!\left[{p_{h}^{*}}\right]\\!\right])\,ds$
$\displaystyle=0\quad\text{ if }F\in\mathcal{F}_{h}^{\Gamma}.$

###### Proof.

For $F\in\mathcal{F}_{h}$ let $\boldsymbol{v}_{F}$ be a test function in the
lowest order Raviart–Thomas finite element such that

$\displaystyle\boldsymbol{v}_{F}\cdot\boldsymbol{n}|_{F^{\prime}}=\begin{cases}1&\text{
if }F^{\prime}=F\\\ 0&\text{ if }F^{\prime}\not=F\end{cases}\quad\text{ for
}F^{\prime}\in\mathcal{F}_{h}.$

If we take this $\boldsymbol{v}_{F}$ in (16a) and
$F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}$, then

$\displaystyle\left(\boldsymbol{u}_{h},\boldsymbol{v}_{F}\right)_{\Omega}-\left(p_{h},\operatorname{div}\boldsymbol{v}_{F}\right)_{\Omega}$
$\displaystyle=0.$

Since $\operatorname{div}\boldsymbol{v}_{F}$ is a piecewise constant function
and the mean-values of $p_{h}^{*}$, $p_{h}$ on every tetrahedron are same, we
have

$\displaystyle\left(\boldsymbol{u}_{h},\boldsymbol{v}_{F}\right)_{\Omega}-\left(p_{h}^{*},\operatorname{div}\boldsymbol{v}_{F}\right)_{\Omega}=0,$

and the element-wise integration by parts gives

$\displaystyle\left(\boldsymbol{u}_{h},\boldsymbol{v}_{F}\right)_{\Omega}-\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],1\right>_{F}+\left(\nabla
p_{h}^{*},\boldsymbol{v}_{F}\right)_{\Omega}=0.$

By the forms of shape functions of the lowest order Raviart–Thomas elements
((10) with $k=1$), $\operatorname{curl}\boldsymbol{v}_{F}=0$ on every
$T\in\mathcal{T}_{h}$, so there exists $\phi_{T}\in\mathcal{P}_{2}(T)$ for
every $T\subset\operatorname{supp}\boldsymbol{v}_{F}$ such that
$\boldsymbol{v}_{F}|_{T}=\nabla\phi_{T}$. By this, (17a), and the above
equation, we have

$\displaystyle-\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],1\right>_{F}=0,$

so (73) follows.

If $F\in\mathcal{F}_{h}^{\Gamma}$, then

$\displaystyle\left(\boldsymbol{u}_{h},\boldsymbol{v}_{F}\right)_{\Omega}+\left<\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n},\boldsymbol{v}_{F}\cdot\boldsymbol{n}\right>_{F}-\left(p_{h},\operatorname{div}\boldsymbol{v}_{F}\right)_{\Omega}$
$\displaystyle=0.$

The integration by parts gives

$\displaystyle\left(\boldsymbol{u}_{h},\boldsymbol{v}_{F}\right)_{\Omega}+\left<\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n}-\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\boldsymbol{v}_{F}\cdot\boldsymbol{n}\right>_{F}-\left(\nabla
p_{h}^{*},\boldsymbol{v}_{F}\right)_{\Omega}$ $\displaystyle=0,$

and we can conclude (74) by a same argument for the proof of (73). ∎

######  Theorem 3.3.

For $\alpha$ satisfying (5) with a uniform lower bound, there exist $C>0$
independent of mesh sizes and $\alpha$ such that the following local
efficiency holds:

(75) $\displaystyle\eta_{T}$ $\displaystyle\leq
C\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0,T},$ (76) $\displaystyle\eta_{F}$
$\displaystyle\leq C\sum_{F\subset\partial
T}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0,T},\quad\text{ if
}F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma},$ (77)
$\displaystyle\eta_{F}$ $\displaystyle\leq C\left(\sum_{F\subset\partial
T}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0,T}+\alpha^{1/2}\|(\boldsymbol{u}-\boldsymbol{u}_{h})\cdot\boldsymbol{n}\|_{0,F}\right),\quad\text{
if }F\in\mathcal{F}_{h}^{\Gamma}.$

###### Proof.

We first show (75). By the triangle inequality $\|\boldsymbol{u}_{h}+\nabla
p_{h}^{*}\|_{0,T}\leq\|(\boldsymbol{u}_{h}-\boldsymbol{u})\|_{0,T}+\|\boldsymbol{u}+\nabla
p_{h}^{*}\|_{0,T}$. By Lemma 3.7 in [14], there exists $C>0$ such that

$\displaystyle\|\boldsymbol{u}_{h}+\nabla p_{h}^{*}\|_{0,T}\leq
C\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{0,T},$

so we have

(78) $\displaystyle\eta_{T}=\|\boldsymbol{u}_{h}+\nabla p_{h}^{*}\|_{0,T}\leq
C\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{0,T}.$

We now prove (76). For $F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}$,
using (73) and $\left[\\!\left[{p}\right]\\!\right]|_{F}=0$, we have (see [1]
or [14, Lemma 3.5])

$\displaystyle\eta_{F}^{2}$
$\displaystyle=\int_{F}h_{F}^{-1}((I-P_{F}^{0})\left[\\!\left[{p_{h}^{*}-p}\right]\\!\right])^{2}\,ds$
$\displaystyle\leq 2\sum_{T\in\mathcal{T}_{h},F\subset\partial
T}h_{F}^{-1}\int_{\partial T\cap F}((I-P_{F}^{0})(p-p_{h}^{*}))^{2}\,ds$ (79)
$\displaystyle\leq C\sum_{T\in\mathcal{T}_{h},F\subset\partial
T}\|\nabla(p-p_{h}^{*})\|_{0,T}^{2}$
$\displaystyle=C\sum_{T\in\mathcal{T}_{h},F\subset\partial
T}\|\boldsymbol{u}+\nabla p_{h}^{*}\|_{0,T}^{2}.$

By the triangle inequality $\|\boldsymbol{u}+\nabla
p_{h}^{*}\|_{0,T}\leq\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0,T}+\|\boldsymbol{u}_{h}+\nabla
p_{h}^{*}\|_{0,T}$ and (78),

(80) $\displaystyle\eta_{F}^{2}$ $\displaystyle\leq
C\sum_{T\in\mathcal{T}_{h},F\subset\partial
T}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0,T}^{2}.$

For (77), using
$\alpha\boldsymbol{u}\cdot\boldsymbol{n}=\left[\\!\left[{p}\right]\\!\right]$
on $\Gamma$ and
$(I-P_{F}^{m})(\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n})=0$,

$\displaystyle\eta_{F}^{2}$
$\displaystyle=\int_{F}\alpha^{-1}\left((I-P_{F}^{m})\left[\\!\left[{p_{h}^{*}}\right]\\!\right]\right)^{2}\,ds$
$\displaystyle=\int_{F}\alpha^{-1}\left((I-P_{F}^{m})(\left[\\!\left[{p_{h}^{*}}\right]\\!\right]-\left[\\!\left[{p}\right]\\!\right]+\alpha\boldsymbol{u}\cdot\boldsymbol{n}-\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n})\right)^{2}\,ds.$

Using (5), we can obtain

(81) $\displaystyle\eta_{F}^{2}$ $\displaystyle\leq
2\sum_{T\in\mathcal{T}_{h},F\subset\partial T}\int_{\partial T\cap
F}\left((I-P_{F}^{m})(p_{h}^{*}-p)\right)^{2}\,ds$
$\displaystyle\quad+2\int_{F}\alpha(\boldsymbol{u}\cdot\boldsymbol{n}-\boldsymbol{u}_{h}\cdot\boldsymbol{n})^{2}\,ds.$

Then, the first term on the right-hand side can be estimated as in (79) with a
natural assumption $h_{F}\lesssim 1$, so (77) follows. ∎

###  3.3. Reliability of post-processed pressure

In this subsection we prove that $\eta$ is a realiable a posteriori error
estimator for the error of post-processed pressure $p_{h}^{*}$. We show a
proof for $n=3$ because the same argument works for $n=2$.

######  Theorem 3.4.

For $p_{h}^{*}$ in (17) and $\eta$ in (19), there exists $C>0$ independent of
mesh sizes such that

(82) $\displaystyle\|p-p_{h}^{*}\|_{0}\leq
C\eta+\frac{C}{\pi}\left(\sum_{T\in\mathcal{T}_{h}}{h_{T}^{2}}\|f-P_{h}f\|_{0,T}^{2}\right)^{\frac{1}{2}}.$

###### Proof.

By the assumption (1), there exists $\boldsymbol{w}\in
H^{1}(\Omega;\mathbb{R}^{3})$ such that

(83)
$\displaystyle\operatorname{div}\boldsymbol{w}=p-p_{h}^{*},\quad\|\boldsymbol{w}\|_{1}\leq
C\|p-p_{h}^{*}\|_{0},\quad\boldsymbol{w}\cdot\boldsymbol{n}|_{\Gamma}=0.$

Using this $\boldsymbol{w}$, the equation $-\nabla p=\boldsymbol{u}$ in
$\Omega$, continuity of $p$ on $\Omega\setminus\Gamma$, and $p=0$ on
$\partial\Omega$,

$\displaystyle\|p-p_{h}^{*}\|_{0}^{2}$
$\displaystyle=\left(p-p_{h}^{*},\operatorname{div}\boldsymbol{w}\right)_{\Omega}$
$\displaystyle=-\left(\nabla(p-p_{h}^{*}),\boldsymbol{w}\right)_{\Omega}\pm\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}\left<\left[\\!\left[{p-p_{h}^{*}}\right]\\!\right],\boldsymbol{w}\cdot\boldsymbol{n}\right>_{F}$
$\displaystyle=\left(\boldsymbol{u}+\nabla
p_{h}^{*},\boldsymbol{w}\right)_{\Omega}\pm\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\boldsymbol{w}\cdot\boldsymbol{n}\right>_{F}$
$\displaystyle=\left(\boldsymbol{u}-\boldsymbol{u}_{h},\boldsymbol{w}\right)_{\Omega}+\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\boldsymbol{w}\right)_{\Omega}\pm\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\boldsymbol{w}\cdot\boldsymbol{n}\right>_{F}.$

Then, (82) follows by

$\displaystyle|\left(\boldsymbol{u}-\boldsymbol{u}_{h},\boldsymbol{w}\right)_{\Omega}|$
$\displaystyle\leq\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0}\|\boldsymbol{w}\|_{0},$
$\displaystyle|\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\boldsymbol{w}\right)_{\Omega}|$
$\displaystyle\leq\left(\sum_{T\in\mathcal{T}_{h}}\eta_{T}^{2}\right)^{\frac{1}{2}}\|\boldsymbol{w}\|_{0},$
$\displaystyle\left|\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\boldsymbol{w}\cdot\boldsymbol{n}\right>_{F}\right|$
$\displaystyle\leq\left(\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}\eta_{F}^{2}\right)^{\frac{1}{2}}\|\boldsymbol{w}\|_{1},$

and (22). ∎

With an assumption of (partial) elliptic regularity, we show that an improved
estimate of $\|p-p_{h}^{*}\|_{0}$ is obtained. Consider the dual problem to
find $(\bar{\boldsymbol{u}},\bar{p})\in\boldsymbol{V}\times Q$ satisfying

(84a)
$\displaystyle\left(\bar{\boldsymbol{u}},\boldsymbol{v}\right)_{\Omega}-\left(\bar{p},\operatorname{div}\boldsymbol{v}\right)_{\Omega}+\left<\alpha\bar{\boldsymbol{u}}\cdot\boldsymbol{n},\boldsymbol{v}\cdot\boldsymbol{n}\right>_{\Gamma}$
$\displaystyle=0$ $\displaystyle\forall\boldsymbol{v}\in\boldsymbol{V},$ (84b)
$\displaystyle\left(\operatorname{div}\bar{\boldsymbol{u}},q\right)_{\Omega}$
$\displaystyle=\left(p-p_{h}^{*},q\right)_{\Omega}$ $\displaystyle\forall q\in
Q,$

and assume that

(85) $\displaystyle\bar{\boldsymbol{u}}\in
H^{\beta}(\Omega;\mathbb{R}^{n}),\quad\frac{1}{2}<\beta\leq
1,\qquad\|\boldsymbol{u}\|_{\beta}\leq C\|p-p_{h}^{*}\|_{0}.$

The boundary condition of this problem is $\bar{p}=0$ on $\partial\Omega$. By
the inf-sup condition (9), there exists $C>0$ such that $\|\bar{p}\|_{0}\leq
C\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\bar{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}$.
Furthermore, by taking $\boldsymbol{v}=\bar{\boldsymbol{u}}$, $q=\bar{p}$, we
can obtain

$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\bar{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}^{2}\leq\left(p-p_{h}^{*},\bar{p}\right)_{\Omega}\leq
C\|p-p_{h}^{*}\|_{0}\|\bar{p}\|_{0},$

so,

(86)
$\displaystyle\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\bar{\boldsymbol{u}}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}\leq
C\|p-p_{h}^{*}\|_{0}.$

######  Corollary 3.7.

Suppose that (85) holds and $\boldsymbol{V}_{h}=\boldsymbol{V}_{h}^{RTN}$ with
$k=1$. Then, there exists $C>0$ independent of mesh sizes such that

$\displaystyle\|p-p_{h}^{*}\|_{0}$
$\displaystyle\leq\left(\sum_{T\in\mathcal{T}_{h}}h_{T}^{2\beta}\eta_{T}^{2}\right)^{\frac{1}{2}}+\left(\sum_{T\in\mathcal{T}_{h}}\operatorname{osc}(f,T)^{2}\right)^{\frac{1}{2}}$
$\displaystyle\quad+\left(\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\min\\{\eta_{F},\alpha^{1/2}h_{F}^{\beta}\eta_{F}\\}^{2}+\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}h_{F}^{2\beta-1}\eta_{F}^{2}\right)^{\frac{1}{2}}.$

###### Proof.

For $\bar{\boldsymbol{u}}$ in (84) and $\tilde{\boldsymbol{u}}$ in (23)

$\displaystyle\|p-p_{h}^{*}\|_{0}^{2}$
$\displaystyle=\left(p-p_{h}^{*},\operatorname{div}\bar{\boldsymbol{u}}\right)_{\Omega}$
$\displaystyle=-\left(\nabla(p-p_{h}^{*}),\bar{\boldsymbol{u}}\right)_{\Omega}+\sum_{T\in\mathcal{T}_{h}}\left<{p-p_{h}^{*}},\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{\partial
T}$ $\displaystyle=\left(\boldsymbol{u}+\nabla
p_{h}^{*},\bar{\boldsymbol{u}}\right)_{\Omega}+\sum_{T\in\mathcal{T}_{h}}\left<{p-p_{h}^{*}},\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{\partial
T}.$ (87)
$\displaystyle=\left(\boldsymbol{u}-\boldsymbol{u}_{h},\bar{\boldsymbol{u}}\right)_{\Omega}+\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\bar{\boldsymbol{u}}\right)_{\Omega}+\sum_{T\in\mathcal{T}_{h}}\left<{p-p_{h}^{*}},\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{\partial
T}.$

Taking $\boldsymbol{v}=\boldsymbol{u}-\boldsymbol{u}_{h}$ in (84a),

$\displaystyle\left(\boldsymbol{u}-\boldsymbol{u}_{h},\bar{\boldsymbol{u}}\right)_{\Omega}$
$\displaystyle=-\left<\alpha\bar{\boldsymbol{u}}\cdot\boldsymbol{n},(\boldsymbol{u}-\boldsymbol{u}_{h})\cdot\boldsymbol{n}\right>_{\Gamma}+\left(\bar{p},\operatorname{div}(\boldsymbol{u}-\boldsymbol{u}_{h})\right)_{\Omega}$
(88)
$\displaystyle=-\left<\alpha\bar{\boldsymbol{u}}\cdot\boldsymbol{n},(\boldsymbol{u}-\boldsymbol{u}_{h})\cdot\boldsymbol{n}\right>_{\Gamma}+\left(\bar{p},f-P_{h}f\right)_{\Omega}.$

By (87), (88), (4), and $\left[\\!\left[{p}\right]\\!\right]|_{F}=0$ for
$F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}$,

$\displaystyle\|p-p_{h}^{*}\|_{0}^{2}$
$\displaystyle=\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\bar{\boldsymbol{u}}\right)_{\Omega}+\left<\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n}-\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{\Gamma}$
$\displaystyle\quad\pm\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}+\left(\bar{p},f-P_{h}f\right)_{\Omega}$
$\displaystyle=:J_{1}+J_{2}+J_{3}+J_{4}.$

To estimate $J_{1}$, note that
$\nabla:\mathcal{P}_{1}(T)\to\mathcal{P}_{0}(T;\mathbb{R}^{n})$ is surjective,
so we get

$\displaystyle|\left(\boldsymbol{u}_{h}+\nabla p_{h}^{*},\bar{u}\right)_{T}|$
$\displaystyle=\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\bar{u}-\nabla\phi\right)_{T}\quad\forall\phi\in\mathcal{P}_{1}(T)$
$\displaystyle\leq C\|\boldsymbol{u}_{h}+\nabla
p_{h}^{*}\|_{0,T}h_{T}^{\beta}\|\bar{\boldsymbol{u}}\|_{H^{\beta}(T)}$
$\displaystyle\leq
Ch_{T}^{\beta}\eta_{T}\|\bar{\boldsymbol{u}}\|_{H^{\beta}(T)}$

where we used (17a), a Poincare inequality for fractional Sobolev spaces ([7,
Lemma 3.1]) with a standard scaling argument, and (20). Then, by the
Cauchy–Schwarz inequality and the integral form of fractional Sobolev norm [8,
Chapter 14]),

$\displaystyle|J_{1}|$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}}|\left(\boldsymbol{u}_{h}+\nabla
p_{h}^{*},\bar{u}\right)_{T}|$ $\displaystyle\leq
C\left(\sum_{T\in\mathcal{T}_{h}}h_{T}^{2\beta}\eta_{T}^{2}\right)^{\frac{1}{2}}\|\bar{\boldsymbol{u}}\|_{H^{\beta}(\Omega;\mathbb{R}^{n})}$
$\displaystyle\leq
C\left(\sum_{T\in\mathcal{T}_{h}}h_{T}^{2\beta}\eta_{T}^{2}\right)^{\frac{1}{2}}\|p-p_{h}^{*}\|_{0}.$

To estimate $J_{2}$ and $J_{3}$, we first note that one can obtain

(89)
$\displaystyle\|\bar{\boldsymbol{u}}\cdot\boldsymbol{n}-P_{F}^{0}\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\|_{0,F}\leq
Ch^{\beta}\|\bar{\boldsymbol{u}}\|_{H^{\beta}(T;\mathbb{R}^{n})},\quad
F\subset\partial T$

by a trace theorem of fractional Sobolev spaces (see, e.g., [27]), the
Poincare inequality for fractional Sobolev spaces, and a standard scaling
argument. Since
$\boldsymbol{u}_{h}\cdot\boldsymbol{n}|_{F}\in\mathcal{P}_{0}(F)$ and
$P_{F}^{0}(\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n}-\left[\\!\left[{p_{h}^{*}}\right]\\!\right])=0$
for $F\in\mathcal{F}_{h}^{\Gamma}$ by (74),

(90)
$\displaystyle\left<\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n}-\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}$
$\displaystyle=\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right]-P_{F}^{0}\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}$
(91)
$\displaystyle=\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right]-P_{F}^{0}\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}-P_{F}^{0}\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}$

for $F\in\mathcal{F}_{h}^{\Gamma}$. From (90) and (91), we can obtain either

$\displaystyle|\left<\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n}-\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}|$
$\displaystyle\leq\eta_{F}\|\alpha^{1/2}\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\|_{0,F}$

or

$\displaystyle|\left<\alpha\boldsymbol{u}_{h}\cdot\boldsymbol{n}-\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}|$
$\displaystyle\leq
C\alpha^{1/2}h_{F}^{\beta}\eta_{F}\|\boldsymbol{u}\|_{H^{\beta}(T;\mathbb{R}^{n})},\quad
F\subset\partial T.$

The mean-value zero property (73) gives

$\displaystyle\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}=\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}-P_{F}^{0}\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}\quad\forall
F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma},$

so we can obtain

$\displaystyle|\left<\left[\\!\left[{p_{h}^{*}}\right]\\!\right],\bar{\boldsymbol{u}}\cdot\boldsymbol{n}\right>_{F}|$
$\displaystyle\leq
C\eta_{F}h_{F}^{\beta-1/2}\|\boldsymbol{u}\|_{H^{\beta}(T;\mathbb{R}^{n})},\quad
F\subset\partial T$

for $F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}$. Combining all of
these estimates,

$\displaystyle|J_{2}+J_{3}|$ $\displaystyle\leq
C\left(\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\min\\{\eta_{F},\alpha^{1/2}h_{F}^{\beta}\eta_{F}\\}^{2}+\sum_{F\in\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{\Gamma}}h_{F}^{2\beta-1}\eta_{F}^{2}\right)^{\frac{1}{2}}$
$\displaystyle\quad\times\|p-p_{h}^{*}\|_{0}.$

Finally, $J_{4}$ can be estimated as

$\displaystyle|J_{4}|$
$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}\frac{h_{T}}{\pi}\|\nabla\bar{p}\|_{0,T}\|f-P_{h}f\|_{0,T}$
$\displaystyle\leq\|\bar{\boldsymbol{u}}\|_{0}\frac{1}{\pi}\left(\sum_{T\in\mathcal{T}_{h}}\operatorname{osc}(f,T)^{2}\right)^{\frac{1}{2}}$
$\displaystyle\leq
C\|p-p_{h}^{*}\|_{0}\frac{1}{\pi}\left(\sum_{T\in\mathcal{T}_{h}}\operatorname{osc}(f,T)^{2}\right)^{\frac{1}{2}}$

by (86), so the conclusion follows. ∎

######  Remark 3.8.

In Corollary 3.7, if $\alpha\ll h_{F}^{2\beta}$, then an improved error bound
with the a posteriori error estimator terms are obtained. However, the data
oscillation error term with $\operatorname{osc}(f,T)$ is same, so data
oscillation can be a dominant factor in error bounds. If $f\in H^{1}(\Omega)$,
then an improved bound
$\operatorname{osc}(f,T)\leq\frac{h_{T}}{\pi}\|f\|_{1,T}$ can be obtained for
each $\operatorname{osc}(f,T)$. This observation can explain convergence of
$\|p-p_{h}^{*}\|_{0}$ faster than the one of
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0}$.

##  4\. Numerical results

In this section we present results of numerical experiments. All experiments
are done with the finite element package FEniCS (version 2019.1.0 [21]). In
particular, the marked elements for refinement are refined using the built-in
adaptive mesh refinement algorithm in FEniCS.

In the first numerical experiment let $\Omega=[0,1]\times[0,1]$ with fault
$\Gamma=\\{1/2\\}\times[1/4,3/4]\subset\Omega$ (see the first figure in Figure
2). The manufactured solution (see the middle and right figures in Figure 2)
for this test case is given by

(92) $\displaystyle p(x,y)=\begin{cases}0&\text{ if }y<\frac{1}{4}\text{ or
}y>\frac{3}{4}\\\ \sin\frac{3\pi
x}{2}\cos^{2}\left(2\pi\left(y-\frac{1}{2}\right)\right)&\text{ if
}x<\frac{1}{2}\text{ and }\frac{1}{4}\leq y\leq\frac{3}{4}\\\
-\sin\frac{3\pi(1-x)}{2}\cos^{2}\left(2\pi\left(y-\frac{1}{2}\right)\right)&\text{
if }x>\frac{1}{2}\text{ and }\frac{1}{4}\leq y\leq\frac{3}{4}.\end{cases}$

Figure 2. The domain with a vertical fault in numerical experiments (left
figure) and the graphs of the pressure field in (92) (middle and right
figures).  Figure 3. Comparison of convergence of errors for uniform and
adaptive refinements. The errors are computed with the manufactured solution
(92).

Figure 4. The first 4 adaptive mesh refinement for $\alpha=0.1$

Figure 5. The first 4 adaptive mesh refinement for $\alpha=10$

Figure 6. The first 4 adaptive mesh refinement for $\alpha=100$

We can compute ${\boldsymbol{u}}=-\nabla p$ and
$f=\operatorname{div}{\boldsymbol{u}}$ on $\Omega\setminus\Gamma$. Note that
the manufactured solutions are not smooth on

$\displaystyle\\{y=3/4,0\leq x\leq 1\,:\,(x,y)\in\Omega\\},$
$\displaystyle\\{y=1/4,0\leq x\leq 1\,:\,(x,y)\in\Omega\\},$
$\displaystyle\\{x=1/2,1/4\leq y\leq 3/4\,:\,(x,y)\in\Omega\\}.$

If we take an initial mesh which includes these segments in the set of edges,
then the manufactured solution $({\boldsymbol{u}},p)$ are smooth on every
triangle. Therefore, the $L^{2}$ errors
$\|{\boldsymbol{u}}-{\boldsymbol{u}}_{h}\|_{0}$ and $\|p-p_{h}\|_{0}$ will
converge with optimal convergence rates for uniform mesh refinement. Since
$p_{h}^{*}$ is an approximation of $p$ which is better than or as good as
$p_{h}$, we only present $\|p-p_{h}^{*}\|_{0}$ in our numerical experiments.
In this experiment, we use the lowest order BDM element for
$\boldsymbol{V}_{h}$ and the piecewise constant element for $Q_{h}$. , so the
optimal convergence rates of $\|{\boldsymbol{u}}-{\boldsymbol{u}}_{h}\|_{0}$
and $\|p-p_{h}\|_{0}$ are 2 and 1, respectively. The $L^{2}$ errors
$\|{\boldsymbol{u}}-{\boldsymbol{u}}_{h}\|_{0}$ and $\|p-p_{h}^{*}\|_{0}$ up
to degrees of freedom are given in Figure 3 (black graph), and one can see
that $\|p-p_{h}^{*}\|_{0}$ shows superconvergence. The errors for adaptive
mesh refinement are also given in Figure 3 (red graph), and one can see that
adaptive mesh refinement gives more optimal convergence of errors up to the
numbers of degrees of freedom. The effectivity index is computed by

$\displaystyle\sqrt{\eta^{2}+\frac{1}{\pi^{2}}\sum_{T\in\mathcal{T}_{h}}\operatorname{osc}(f,T)^{2}}/\|{\boldsymbol{u}}-{\boldsymbol{u}}_{h}\|_{0},$

and is asymptotically $1.57$.

In the second set of experiments, we present mesh adaptivity for nonsmooth
solutions. Since it is difficult to construct nonsmooth manufactured solutions
with the fault structure, we show adaptive mesh refinement by our a posteriori
error estimator for numerical solutions with given boundary conditions,
$\alpha=0.1,10,100$, and $f=1$. Assuming $\Omega=[0,1]\times[0,1]$ with the
same $\Gamma$, zero flux boundary conditions are imposed on the top and bottom
boundary components $\\{(x,y)\in\Omega\,:\,0\leq x\leq 1,y=0\text{ or }y=1\\}$
of $\partial\Omega$, and $p=0$ on the left side, $p=-1$ on the right side, are
imposed.

In Figures 4–6, we present the first 4 adaptive mesh refinements for the three
cases. In all cases, most mesh refinements are done around the two endpoints
of $\Gamma$. One can see that mesh refinements are more concentrated at the
two endpoints of $\Gamma$ as $\alpha$ is larger. As $\alpha$ is larger,
$\Gamma$ is more resistant against fluid flow across it. This is observed in
the numerical experiment results because color contrast across $\Gamma$,
corresponding to pressure jump, gets stronger as $\alpha$ is larger. Then, it
is naturally expected that fluid flow is more singular near the two endpoints
as $\alpha$ is larger. This explains the differences of mesh refinement
concentration for different $\alpha$ values. One can also see that there are
not many mesh adaptation in the interior of $\Gamma$. Nonetheless a pressure
jump term is involved in solutions around the interior of $\Gamma$, these
numerical test results show that solutions around the interior of $\Gamma$ are
not so singular.

##  5\. Conclusion

In this work we studied a recovery type a posteriori error estimator of the
Darcy flow model with Robin-type interface conditions. The reliability and the
local efficiency of the estimator are proved. In contrast to the previous work
in [18], we developed a new $H(\operatorname{div})$-based proof using a
modified Helmholtz decomposition, a modified Scott–Zhang interpolation,
edge/face-wise integration by parts. Moreover, we proved that the post-
processed pressure is bounded by the estimator, and a superconvergent upper
bound can be obtained under a (partial) elliptic regularity assumption of the
dual problem. Numerical test results are included to illustrate the adaptivity
results of our estimator.

##  6\. Appendix: integration by parts identities

In this section we present identities from the integration by parts that we
used in the paper.

In this section $(a_{1}\;\cdots\;a_{n})^{t}$ denotes the column vector with
entries $a_{1},\cdots,a_{n}$. For differentiable functions
$\phi:\mathbb{R}^{2}\to\mathbb{R}$, $\Psi:\mathbb{R}^{2}\to\mathbb{R}^{2}$
with $\Psi=(\Psi_{1}\;\Psi_{2})^{t}$, we define $\operatorname{curl}$ and
$\operatorname{rot}$ by

(93) $\displaystyle\operatorname{curl}\phi=\begin{pmatrix}-\partial_{y}\phi\\\
\partial_{x}\phi\end{pmatrix},\qquad\operatorname{rot}\Psi=-\partial_{y}\Psi_{1}+\partial_{x}\Psi_{2}.$

For a triangle $T\subset\mathbb{R}^{2}$, $\boldsymbol{n}_{\partial T}$ is the
outward unit normal vector field on $\partial T$ and $\boldsymbol{t}_{\partial
T}$ is the unit tangential vector field along the counterclock-wise direction
of $\partial T$. Denoting $\boldsymbol{n}_{\partial T}$ by
$\boldsymbol{n}_{\partial T}=(n_{1}\;n_{2})^{t}$, note that
$\boldsymbol{t}_{\partial T}=(-n_{2}\;n_{1})^{t}$. By the integration by
parts,

$\displaystyle\int_{T}\operatorname{curl}\phi\cdot{\Psi}\,ds$
$\displaystyle=\int_{T}(-\partial_{y}\phi\Psi_{1}+\partial_{x}\phi\Psi_{2})\,ds$
(94) $\displaystyle=\int_{T}\nabla\phi\cdot\begin{pmatrix}\Psi_{2}\\\
-\Psi_{1}\end{pmatrix}\,ds$ $\displaystyle=\int_{\partial
T}\phi\begin{pmatrix}n_{1}\\\
n_{2}\end{pmatrix}\cdot\begin{pmatrix}\Psi_{2}\\\
-\Psi_{1}\end{pmatrix}\,dl-\int_{T}\phi(\partial_{x}\Psi_{2}-\partial_{y}\Psi_{1})\,ds$
$\displaystyle=\int_{\partial T}\phi\begin{pmatrix}-n_{2}\\\
n_{1}\end{pmatrix}\cdot\begin{pmatrix}\Psi_{1}\\\
\Psi_{2}\end{pmatrix}\,dl-\int_{T}\phi\operatorname{rot}{\Psi}\,ds$
$\displaystyle=\int_{\partial T}\phi\boldsymbol{t}_{\partial
T}\cdot\Psi\,dl-\int_{T}\phi\operatorname{rot}\Psi\,ds$

and for $E\subset\partial T$,

$\displaystyle\int_{E}\operatorname{curl}\phi\cdot\boldsymbol{n}_{\partial
T}\,dl$ $\displaystyle=\int_{E}\begin{pmatrix}-\partial_{y}\phi\\\
\partial_{x}\phi\end{pmatrix}\cdot\begin{pmatrix}n_{1}\\\
n_{2}\end{pmatrix}\,dl$ (95)
$\displaystyle=\int_{E}\nabla\phi\cdot\begin{pmatrix}n_{2}\\\
-n_{1}\end{pmatrix}\,dl$
$\displaystyle=-\int_{E}\nabla\phi\cdot\boldsymbol{t}_{\partial T}\,dl.$

Let $F$ be a triangle in the $xy$-plane in $\mathbb{R}^{3}$ and
$\boldsymbol{n}_{\partial F}=(n_{1}\;n_{2}\;0)^{t}$ be the unit outward normal
vector field of $F$ in $\mathbb{R}^{3}$. The tangential vector field on
$\partial F$ in $\mathbb{R}^{3}$ is $\boldsymbol{t}_{\partial
F}=(-n_{2}\;n_{1}\;0)^{t}$. For differentiable functions
$\phi:\mathbb{R}^{3}\to\mathbb{R}$, $\Psi:\mathbb{R}^{3}\to\mathbb{R}^{3}$
with $\Psi=({\Psi_{1}\;\Psi_{2}\;\Psi_{3}})^{t}$,
$\boldsymbol{n}_{F}:=({0\;0\;1})^{t}$, we get

(96) $\displaystyle\boldsymbol{n}_{F}\times\Psi=\begin{pmatrix}-\Psi_{2}\\\
\Psi_{1}\\\
0\end{pmatrix},\quad\operatorname{curl}\Psi\cdot\boldsymbol{n}_{F}=\partial_{x}\Psi_{2}-\partial_{y}\Psi_{1}.$

By these identities and the 3rd equality in (94),

$\displaystyle\int_{F}\operatorname{curl}\Psi\cdot\boldsymbol{n}_{F}\phi\,ds$
$\displaystyle=\int_{F}(\partial_{x}\Psi_{2}-\partial_{y}\Psi_{1})\phi\,ds$
(97) $\displaystyle=-\int_{F}\begin{pmatrix}-\partial_{y}\phi\\\
\partial_{x}\phi\end{pmatrix}\cdot\begin{pmatrix}\Psi_{1}\\\
\Psi_{2}\end{pmatrix}\,ds+\int_{\partial F}\phi\begin{pmatrix}n_{1}\\\
n_{2}\end{pmatrix}\cdot\begin{pmatrix}\Psi_{2}\\\ -\Psi_{1}\end{pmatrix}\,dl$
$\displaystyle=\int_{F}\begin{pmatrix}\partial_{x}\phi\\\
\partial_{y}\phi\end{pmatrix}\cdot\begin{pmatrix}-\Psi_{2}\\\
\Psi_{1}\end{pmatrix}\,ds+\int_{\partial F}\phi\begin{pmatrix}-n_{2}\\\
n_{1}\end{pmatrix}\cdot\begin{pmatrix}\Psi_{1}\\\ \Psi_{2}\end{pmatrix}\,dl$
$\displaystyle=\int_{F}\nabla\phi\cdot(\boldsymbol{n}_{F}\times\Psi)\,ds-\int_{\partial
F}\phi\boldsymbol{t}_{\partial F}\cdot\Psi\,dl.$

a

Funding Jeonghun J. Lee gratefully acknowledge support from the National
Science Foundation (DMS-2110781).

## References

  * [1] Mark Ainsworth, _A posteriori error estimation for lowest order Raviart-Thomas mixed finite elements_ , SIAM J. Sci. Comput. 30 (2007/08), no. 1, 189–204. MR 2377438
  * [2] Philippe Angot, Franck Boyer, and Florence Hubert, _Asymptotic and numerical modelling of flows in fractured porous media_ , M2AN Math. Model. Numer. Anal. 43 (2009), no. 2, 239–275. MR 2512496
  * [3] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther, _Mixed finite element methods for linear elasticity with weakly imposed symmetry_ , Math. Comp. 76 (2007), no. 260, 1699–1723 (electronic). MR 2336264 (2008k:74057)
  * [4] by same author, _Finite element exterior calculus: from Hodge theory to numerical stability_ , Bull. Amer. Math. Soc. (N.S.) 47 (2010), no. 2, 281–354. MR 2594630 (2011f:58005)
  * [5] I. Babuška, J. Osborn, and J. Pitkäranta, _Analysis of mixed methods using mesh dependent norms_ , Math. Comp. 35 (1980), no. 152, 1039–1062. MR 583486
  * [6] M. Bebendorf, _A note on the Poincaré inequality for convex domains_ , Z. Anal. Anwendungen 22 (2003), no. 4, 751–756. MR 2036927 (2004k:26025)
  * [7] José C. Bellido and Carlos Mora-Corral, _Existence for nonlocal variational problems in peridynamics_ , SIAM J. Math. Anal. 46 (2014), no. 1, 890–916. MR 3166960
  * [8] Susanne C. Brenner and L. Ridgway Scott, _The mathematical theory of finite element methods_ , Third ed., Springer, 2008. MR 515228 (80k:35056)
  * [9] F. Brezzi and M. Fortin, _Mixed and hybrid finite element methods_ , Springer Series in computational Mathematics, vol. 15, Springer, 1992. MR MR2233925 (2008i:35211)
  * [10] Franco Brezzi, Jr. Jim Douglas, and L. D. Marini, _Two families of mixed finite elements for second order elliptic problems_ , Numer. Math. 47 (1985), no. 2, 217–235. MR 799685 (87g:65133)
  * [11] Andrea Cangiani, Emmanuil H. Georgoulis, and Max Jensen, _Discontinuous Galerkin methods for mass transfer through semipermeable membranes_ , SIAM J. Numer. Anal. 51 (2013), no. 5, 2911–2934. MR 3121762
  * [12] Andrea Cangiani, Emmanuil H. Georgoulis, and Younis A. Sabawi, _Adaptive discontinuous Galerkin methods for elliptic interface problems_ , Math. Comp. 87 (2018), no. 314, 2675–2707. MR 3834681
  * [13] Huangxin Chen and Shuyu Sun, _A residual-based a posteriori error estimator for single-phase Darcy flow in fractured porous media_ , Numer. Math. 136 (2017), no. 3, 805–839. MR 3660303
  * [14] Bernardo Cockburn and Wujun Zhang, _An a posteriori error estimate for the variable-degree Raviart-Thomas method_ , Math. Comp. 83 (2014), no. 287, 1063–1082. MR 3167450
  * [15] Evan Gawlik, Michael J. Holst, and Martin W. Licht, _Local finite element approximation of Sobolev differential forms_ , ESAIM Math. Model. Numer. Anal. 55 (2021), no. 5, 2075–2099. MR 4319601
  * [16] F. Hecht, Z. Mghazli, I. Naji, and J. E. Roberts, _A residual a posteriori error estimators for a model for flow in porous media with fractures_ , J. Sci. Comput. 79 (2019), no. 2, 935–968. MR 3968997
  * [17] Kwang-Yeon Kim, _Guaranteed a posteriori error estimator for mixed finite element methods of elliptic problems_ , Appl. Math. Comput. 218 (2012), no. 24, 11820–11831. MR 2945185
  * [18] Juho Könnö, Dominik Schötzau, and Rolf Stenberg, _Mixed finite element methods for problems with Robin boundary conditions_ , SIAM J. Numer. Anal. 49 (2011), no. 1, 285–308. MR 2783226
  * [19] Jeonghun J. Lee, Tan Bui-Thanh, Umberto Villa, and Omar Ghattas, _Forward and inverse modeling of fault transmissibility in subsurface flows_ , Comput. Math. Appl. 128 (2022), 354–367. MR 4512460
  * [20] by same author, _Forward and inverse modeling of fault transmissibility in subsurface flows_ , Comput. Math. Appl. 128 (2022), 354–367. MR 4512460
  * [21] Anders Logg, Kent-Andre Mardal, and Garth N. Wells (eds.), _Automated solution of differential equations by the finite element method_ , Lecture Notes in Computational Science and Engineering, vol. 84, Springer, Heidelberg, 2012, The FEniCS book. MR 3075806
  * [22] Carlo Lovadina and Rolf Stenberg, _Energy norm a posteriori error estimates for mixed finite element methods_ , Math. Comp. 75 (2006), no. 256, 1659–1674 (electronic). MR 2240629 (2007h:65129)
  * [23] Vincent Martin, Jérôme Jaffré, and Jean E. Roberts, _Modeling fractures and barriers as interfaces for flow in porous media_ , SIAM J. Sci. Comput. 26 (2005), no. 5, 1667–1691. MR 2142590
  * [24] J.-C. Nédélec, _Mixed finite elements in ${\bf R}^{3}$_, Numer. Math. 35 (1980), no. 3, 315–341. MR 592160 (81k:65125)
  * [25] by same author, _A new family of mixed finite elements in ${\bf R}^{3}$_, Numer. Math. 50 (1986), no. 1, 57–81. MR 864305 (88e:65145)
  * [26] P.-A. Raviart and J. M. Thomas, _A mixed finite element method for 2nd order elliptic problems_ , Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975), Springer, Berlin, 1977, pp. 292–315. Lecture Notes in Math., Vol. 606\. MR 0483555 (58 #3547)
  * [27] Michael Renardy and Robert C. Rogers, _An introduction to partial differential equations_ , second ed., Texts in Applied Mathematics, vol. 13, Springer-Verlag, New York, 2004. MR 2028503
  * [28] L. Ridgway Scott and Shangyou Zhang, _Finite element interpolation of nonsmooth functions satisfying boundary conditions_ , Math. Comp. 54 (1990), no. 190, 483–493. MR 1011446
  * [29] Rolf Stenberg, _Postprocessing schemes for some mixed finite elements_ , RAIRO Modél. Math. Anal. Numér. 25 (1991), no. 1, 151–167. MR 1086845 (92a:65303)
  * [30] Martin Vohralík, _A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations_ , SIAM J. Numer. Anal. 45 (2007), no. 4, 1570–1599. MR 2338400
  * [31] Lina Zhao and Eric Chung, _An adaptive discontinuous Galerkin method for the Darcy system in fractured porous media_ , Comput. Geosci. 26 (2022), no. 6, 1581–1596. MR 4510543