# An Eulerian finite element method for the linearized Navier–Stokes problem
in an evolving domain ††thanks: Partially supported by the NSF through the
grants DMS-2309425 (Neilan) and DMS-2309197 (Olshanskii)

Michael Neilan Department of Mathematics, University of Pittsburgh,
Pittsburgh, PA 15260, USA. Email<EMAIL_ADDRESS>sites.pitt.edu/~neilan
Maxim Olshanskii Department of Mathematics, University of Houston, Houston,
Texas 77204-3008, USA. Email<EMAIL_ADDRESS>www.math.uh.edu/~molshan

###### Abstract

The paper addresses an error analysis of an Eulerian finite element method
used for solving a linearized Navier–Stokes problem in a time-dependent
domain. In this study, the domain’s evolution is assumed to be known and
independent of the solution to the problem at hand. The numerical method
employed in the study combines a standard Backward Differentiation Formula
(BDF)-type time-stepping procedure with a geometrically unfitted finite
element discretization technique. Additionally, Nitsche’s method is utilized
to enforce the boundary conditions. The paper presents a convergence estimate
for several velocity–pressure elements that are inf-sup stable. The estimate
demonstrates optimal order convergence in the energy norm for the velocity
component and a scaled $L^{2}(H^{1})$-type norm for the pressure component.

###### keywords:

Interface Stokes problem, evolving interface, cutFEM

###### AMS:

65M12, 65M15, 65M60

##  1 Introduction

Fluid equations formulated in time-dependent domains are essential components
of mathematical models used in a wide range of applications, including
cardiovascular research and aerospace engineering [15, 2]. The analysis of
such equations is a classical topic in mathematical fluid mechanics [34, 28,
35, 29]. Moreover, a significant body of literature addresses the development
of computational methods aimed at numerically solving these problems. Well-
established numerical techniques include immersed boundary methods, fictitious
domain methods, methods based on Lagrangian and arbitrary Lagrangian–Eulerian
formulations, space–time finite element formulations, level-set methods, and
extended finite element methods; see, e.g., [21, 31, 37, 27, 17, 14, 13, 19].

In this paper, we focus on an Eulerian finite element method that utilizes a
time-independent triangulation of $\mathbb{R}^{3}$ to solve a system of
governing equations within a volume $\Omega(t)$ that smoothly evolves through
the background mesh, a typical configuration for unfitted finite element
methods. Specifically, we consider the CutFEM unfitted finite element method
[8] that incorporates Nitsche’s method for boundary condition imposition and
employs a ghost-penalty stabilization [7] to handle instabilities arising from
arbitrary small “cuts” made by $\Omega(t)$ within the background simplices.

For time stepping, we adopt an Eulerian procedure suggested in [23] that
relies on the implicit extension of the solution from $\Omega(t)$ to its
neighborhood of $\mathcal{O}(\Delta t)$. This combination of the CutFEM method
and implicit extension-based time stepping was initially applied to two-phase
flow problems in[11], demonstrating its efficacy when used in conjunction with
the level-set method for interface capturing. Recent studies in [9] and [38]
have addressed the analysis of this method, considering equal-order stabilized
and Taylor-Hood elements, respectively. Both of these analyses identified a
major challenge: the lack of a weak divergence-free property of the time
difference of the finite element solutions
$(\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1})/\Delta t$ with respect to the
discrete pressure space at time $t^{n}$. The absence of this property makes it
challenging to bound this term in a suitable norm and precluding optimal-order
estimates for the pressure. This observation has also been made in the
literature on adaptive-in-time finite element methods, where the pressure
space varies in time due to mesh adaptation [3, 4]. The use of equal-order
finite elements and pressure stabilization in [9] allows the authors to
establish the optimal error estimate for velocity. However, for inf-sup stable
Taylor-Hood elements, the coupling between pressure and velocity appears
stronger, and the sub-optimality in pressure also hindered the authors of [38]
from obtaining the optimal order estimation for the velocity error. It is
worth noting that [38] also quantified the error resulting from an approximate
reconstruction of the evolving “exact” domain, $\Omega(t)$.

Despite the aforementioned theoretical challenges, numerical experiments have
demonstrated optimal order convergence rates [38]. This raises the question of
whether the analysis can be enhanced to provide support for the observed
numerical evidence. This is the question addressed in the present paper. The
setup of the problem and the methods here is similar to [38], but we consider
general inf-sup stable unfitted finite element pairs, essentially those
covered in the analysis by Guzmán et al. [20].

The main result established in this paper can be summarized as follows:
Optimal convergence rates are proven for the energy norm of velocity and a
scaled $L^{2}(H^{1})$-norm of the pressure under the constraint
$h^{2}\lesssim\Delta t\lesssim h$, where $h$ represents the mesh size and
$\Delta t$ denotes the time step. This bridges the gap in the analysis up to
the selection of the pressure norm. Notably, the use of a non-standard
pressure norm is vital in mitigating the lack of divergence-free property in
the discrete time derivative. This argument aligns with the analysis in a
recent study [30], which analyzed a finite element method for the Navier-
Stokes equations posed on time-dependent surfaces.

In general, there is a scarcity of literature addressing error bounds for
fully discrete solutions of fluid equations in evolving domains. However,
under the simplifying assumption that the motion of the domain is given and
decoupled from the flow solution, error bounds for the Arbitrary Lagrangian-
Eulerian (ALE) and quasi-Lagrangian finite element methods for Stokes, Navier-
Stokes and coupled Stokes–parabolic equations in moving domains can be found
in [33, 24, 22]. Similarly, error bounds for the unfitted characteristic
finite element method within the same setup are provided in [25].

The remainder of the paper is organized in five sections and an appendix.
Section 2 formulates the linearized Navier-Stokes problem in evolving domains
and introduces suitable extension operators utilized in the analysis. In
particular, the numerical analysis relies on existence of a sufficiently
regular divergence-free extension of the fluid velocity field in a
neighborhood of $\Omega(t)$. The fully discrete numerical method based on a
Nitsche-based CutFEM formulation is given in Section 3. Here, we present the
scheme for general finite element Stokes pairs satisfying certain assumptions.
Stability and convergence analysis is the subject of Section 4. In Section 5,
we list three standard finite element pairs satisfying the assumptions.
Finally, a proof of a ‘discrete’ trace estimate is found in Appendix A.

##  2 Problem formulation

We consider a time-dependent domain $\Omega(t)\subset\mathbb{R}^{3}$ with
boundary $\Gamma(t):=\partial\Omega(t)$ whose motion is assumed to be known a
priori. In particular, we assume a smooth solenoidal vector field
$\mathbf{w}:\mathbb{R}^{3}\times[0,T]\to\mathbb{R}^{3}$, for some final time
$T>0$ such that the normal velocity of the boundary is specified via

$V_{\Gamma}=\mathbf{w}\cdot{\mathbf{n}}_{\Gamma}\quad\text{on}~{}\Gamma(t),$
(1)

where $\mathbf{n}_{\Gamma}$ denotes the outward unit normal of $\Gamma(t)$. We
then consider the Oseen problem in the moving volume $\Omega(t)$:

$\begin{split}\mathbf{u}_{t}+(\mathbf{w}\cdot\nabla)\mathbf{u}-\Delta\mathbf{u}+\nabla
p&=\mathbf{f}~{}\quad\text{in }\Omega(t),\\\ \operatorname{\rm
div}\mathbf{u}&=0~{}\quad\text{in }\Omega(t),\\\
\mathbf{u}&=\mathbf{w}\quad\text{on }\Gamma(t),\end{split}$ (2)

with initial condition $\mathbf{u}|_{t=0}=\mathbf{u}_{0}$ in
$\Omega_{0}:=\Omega(0)$. As mentioned in the introduction, unfitted finite
element methods for (2) were recently addressed in [9, 38] with suboptimal
error bounds. We note that the previous studies [9, 38] ignore the advection
term $(\mathbf{w}\cdot\nabla)\mathbf{u}$ in (2). While this term does not lead
to any additional difficulties in the analysis, we believe it is mechanically
relevant to include it in this simplified model. By a standard argument, we
can re-write the above problem for

$\mathbf{u}=0\quad\text{on}~{}~{}\Gamma(t).$ (3)

We assume the smooth velocity field
$\mathbf{w}\,:\,\mathbb{R}^{3}\times[0,T]\to\mathbb{R}^{3}$ is such that it
defines the flow map $\Phi_{t}:\,\Omega(0)\to\Omega(t)$ as the material
evolution of the fluid volume: For $\mathbf{z}\in\Omega_{0}$, the trajectory
$\mathbf{x}(t,\mathbf{z})=\Phi_{t}(\mathbf{z})$ solves

$\displaystyle\begin{cases}\mathbf{x}(0,\mathbf{z})=\mathbf{z},\\\
\frac{d}{dt}\mathbf{x}(t,\mathbf{z})=\mathbf{w}(t,\mathbf{x}(\mathbf{z},t))\qquad
t\in(0,T]\end{cases}$ (4)

for some final time $T>0$. Equation (4) defines a smooth bijection between
$\Omega_{0}$ and $\Omega(t)$ for every $t\in[0,T]$. If $\partial\Omega_{0}\in
C^{p}$ and $\mathbf{w}\in\mathbf{C}^{p}(\mathbb{R}^{3})$, then $\Gamma(t)\in
C^{p}$; the flow map $\Phi_{t}$ also preserves the connectivity of
$\Omega(t)$.

Summarizing, we are interested in the analysis of a finite element method for
solving (2) with $\Omega(t)=\Phi_{t}(\Omega(0))$ and homogeneous Dirichlet
boundary conditions (3).

###  2.1 Extensions

Let $\Omega(t)\subset\widehat{\Omega}$ for all $t\in[0,T]$, for a bounded
polyhedral domain $\widehat{\Omega}\subset\mathbb{R}^{3}$. We define the two
space-time domains $\mathcal{Q}$ and $\hat{\mathcal{Q}}$ as follows:

$\mathcal{Q}:=\bigcup\limits_{t\in[0,T]}\Omega(t)\times\\{t\\}\subset\widehat{\mathcal{Q}}:=\widehat{\Omega}\times[0,T]\subset\mathbb{R}^{4}.$

For a domain $D\subset\mathbb{R}^{3}$ and some $\delta>0$ we use the notation
$\mathcal{O}_{\delta}(D)$ for the $\delta$-neighborhood of $D$:

$\mathcal{O}_{\delta}(D)=\\{x\in\mathbb{R}^{3}:\ {\rm
dist}(x,D)\leq\delta\\}.$

Denoting by $\mathbf{V}(t)=\\{\mathbf{v}\in\mathbf{H}^{1}_{0}(\Omega(t)):\
{\rm div}\,\mathbf{v}=0\\}$, the subspace of divergence-free functions in
$\mathbf{H}^{1}_{0}(\Omega(t))$, our goal now is to define an extension
operator $\mathcal{E}\,:\,\mathbf{V}(t)\to\mathbf{H}^{1}(\widehat{\Omega})$
that preserves the divergence–free condition. To this end, we note that since
$\operatorname{\rm div}\mathbf{u}=0$, we can write
$\mathbf{u}=\nabla\times{\bm{\psi}}$ in $\Omega(t)$ with a stream function
that satisfies ${\bm{\psi}}\in{\bf W}^{k+1,p}(\Omega(t))$ and

$\|{\bm{\psi}}\|_{W^{k+1,p}(\Omega(t))}\lesssim\|\mathbf{u}\|_{W^{k,p}(\Omega(t))}\quad\text{for}~{}\mathbf{u}\in
W^{k,p}(\Omega(t)),$ (5)

$k\geq 0$, $1<p<\infty$; see [16, 12].

######  Remark 2.1.

Here, the statement $A\lesssim B$ (resp., $A\gtrsim B$) to mean $A\leq cB$
(resp., $A\geq cB$) for some constant $c>0$ independent of the spatial and
temporal discretization parameters $h$ and $\Delta t$ introduced below and
time $t$. The statement $A\simeq B$ means $A\lesssim B$ and $A\gtrsim B$.

For ${\bm{\psi}}_{0}={\bm{\psi}}\circ\Phi_{t}$ we consider Stein’s extension:
Since the boundary of $\Omega_{0}$ is smooth, there is a continuous linear
extension operator $\mathcal{E}_{0}:L^{2}(\Omega_{0})\to
L^{2}(\mathbb{R}^{3})$, ($\mathcal{E}_{0}{\bm{\psi}}_{0}={\bm{\psi}}_{0}$ in
$\Omega_{0}$), with the following properties [36, Section VI.3.1]:

$\|{\mathcal{E}_{0}}{\bm{\psi}}_{0}\|_{W^{k,p}(\mathbb{R}^{3})}\leq
C_{\Omega_{0}}\|{\bm{\psi}}_{0}\|_{W^{k,p}(\Omega_{0})},\quad\text{for}~{}{\bm{\psi}}_{0}\in
W^{k,p}(\Omega_{0}),~{}~{}k=0,\dots,m+1,~{}~{}1\leq p\leq\infty,$ (6)

with any fixed $m\geq 0$. Here, the extension operator is performed component-
wise, i.e.,
$(\mathcal{E}_{0}{\bm{\psi}}_{0})_{i}=\mathcal{E}_{0}({\bm{\psi}}_{0})_{i}$
for $i=1,2,3$. For the extension
$\mathcal{E}_{\psi}{\bm{\psi}}:=(\mathcal{E}_{0}{\bm{\psi}}_{0})\circ\Phi^{-1}_{t}$
of ${\bm{\psi}}$, the following estimates follow from the analysis in [23]:

$\displaystyle\|\mathcal{E}_{\psi}{\bm{\psi}}\|_{H^{k}(\widehat{\Omega})}$
$\displaystyle\lesssim\|{\bm{\psi}}\|_{H^{k}(\Omega(t))},\quad{\small
k=0,\dots,m+1},\quad\|\mathcal{E}_{\psi}{\bm{\psi}}\|_{W^{4,5}(\widehat{\mathcal{Q}})}\lesssim\|{\bm{\psi}}\|_{W^{4,5}(\mathcal{Q})},$
(7)
$\displaystyle\|(\mathcal{E}_{\psi}{\bm{\psi}})_{t}\|_{H^{m}(\widehat{\Omega})}$
$\displaystyle\lesssim(\|{\bm{\psi}}\|_{H^{m+1}(\Omega(t))}+\|{\bm{\psi}}_{t}\|_{H^{m}(\Omega(t))}).$

We now define the velocity extension as follows

$\mathcal{E}{}\mathbf{u}(t):=\nabla\times(\mathcal{E}_{\psi}{\bm{\psi}}),\quad\text{for
each}~{}t\in[0,T].$ (8)

By construction there holds

$\operatorname{\rm
div}\mathcal{E}\mathbf{u}=0\quad\text{in}\,\widehat{\Omega}.$

For $\mathbf{u}\in L^{\infty}(0,T;\mathbf{H}^{m}(\Omega(t)))\cap{\bf
W}^{3,5}(\mathcal{Q})$ such that $\operatorname{\rm div}\mathbf{u}=0$ in
$\Omega(t)$ for all $t\in(0,T)$ and any fixed integer $m\geq 0$, the following
estimates follow from (5), (7), Poincare-Friedrich’s inequality, and the
embedding $\ W^{3,5}(\widehat{\mathcal{Q}})\subset
W^{2,\infty}(\widehat{\mathcal{Q}})$:

$\displaystyle\|\mathcal{E}\mathbf{u}\|_{H^{k}(\widehat{\Omega})}$
$\displaystyle\lesssim\|\mathbf{u}\|_{H^{k}(\Omega(t))},\quad{\small
k=0,\dots,m},$ (9a)
$\displaystyle\|\nabla(\mathcal{E}\mathbf{u})\|_{\widehat{\Omega}}$
$\displaystyle\lesssim\|\nabla\mathbf{u}\|_{\Omega(t)},$ (9b)
$\displaystyle\|\mathcal{E}\mathbf{u}\|_{W^{2,\infty}(\widehat{\mathcal{Q}})}$
$\displaystyle\lesssim\|\mathbf{u}\|_{W^{3,5}(\mathcal{Q})},$ (9c)

Here, we use the standard notation for the $L^{2}$-norm
$\|\cdot\|_{D}=\|\cdot\|_{L^{2}(D)}$ for some domain $D$. Furthermore, for
$\mathbf{u}\in L^{\infty}(0,T;\mathbf{H}^{m}(\Omega(t)))$ such that
$\mathbf{u}_{t}\in L^{\infty}(0,T;\mathbf{H}^{m-1}(\Omega(t)))$ it holds

$\|(\mathcal{E}\mathbf{u})_{t}\|_{H^{m-1}(\widehat{\Omega})}\lesssim(\|\mathbf{u}\|_{H^{m}(\Omega(t))}+\|\mathbf{u}_{t}\|_{H^{m-1}(\Omega(t))}).$
(10)

With an abuse of notation, we define the extension of the pressure as

$\displaystyle\mathcal{E}p(t)=(\mathcal{E}_{0}(p\circ\Phi_{t}))\circ\Phi_{t}^{-1},\quad\text{for
each}~{}t\in[0,T].$ (11)

Then estimates (9a),(9c), with $\mathcal{E}\mathbf{u}$ and $\mathbf{u}$
replaced by $\mathcal{E}p$ and $p$, respectively, are satisfied (cf. [23,
Lemma 3.3]). For the analysis, we need $\mathcal{E}\mathbf{u}$ and
$\mathcal{E}p$ defined in
$\mathcal{O}_{\delta}(\Omega(t))\subset\widehat{\Omega}$, a
$\delta$-neighborhood of $\Omega(t)$ with $\delta\simeq\Delta t$.

##  3 The Fully Discrete Finite Element Method

We adopt the basic framework in [23, 38, 9] to build a Nitsche-based CutFEM
spatial discretization of the Stokes problem on an evolving domain.

###  3.1 Approximate geometries

Recall that $\widehat{\Omega}\subset\mathbb{R}^{3}$ is a polyhedral domain
with $\Omega(t)\subset\widehat{\Omega}$ for all $t\in[0,T]$. For simplicity,
we consider a time discretization with a uniform time-step $\Delta t=T/N$ for
some $N\in\mathbb{N}$. We set $t_{n}=n\Delta t$, $\Omega^{n}=\Omega(t_{n})$,
$\Gamma^{n}=\Gamma(t_{n})$, and
$(\mathbf{u}^{n},p^{n})=(\mathbf{u}(t_{n}),p(t_{n}))$. We further set
$\mathbf{w}^{n}_{\infty}=\|\mathbf{w}(t_{n})\cdot\mathbf{n}_{\Gamma}\|_{L^{\infty}(\Gamma^{n})}$.
For practical purposes such as numerical integration, and similar to [23, 38,
9], we assume that the domains $\Omega^{n}$ are given by their approximations
$\Omega_{h}^{n}$ (cf. (12)–(13) below). The boundary of $\Omega_{h}^{n}$ is
denoted by $\Gamma_{h}^{n}$. To ensure that discrete solutions are well
defined at subsequent time-steps, we extend the computational domain by a
layer of thickness $\delta_{h}$ with
$c_{\delta_{h}}\mathbf{w}^{n}_{\infty}\Delta t\leq\delta_{h}$ with constant
$1\leq c_{\delta_{h}}=O(1)$ such that ${\rm
dist}(\Omega_{h}^{n},\Omega_{h}^{n+1})\leq\delta_{h}$ for all $n$.

We assume there is a bijective, Lipschitz continuous map
$\Psi_{n}:\mathcal{O}_{\delta_{h}}(\Omega_{h}^{n})\to\mathcal{O}_{\delta_{h}}(\Omega^{n})$
that connects the approximate and exact domains at each time step. In
particular, we assume $\Psi_{n}$ satisfies
$\mathcal{O}_{\delta_{h}}(\Omega^{n})=\Psi_{n}(\mathcal{O}_{\delta_{h}}(\Omega_{h}^{n}))$,
$\Omega^{n}=\Psi_{n}(\Omega_{h}^{n})$, $\Gamma^{n}=\Psi_{n}(\Gamma_{h}^{n})$,
and the existence of a positive integer $q$ such that

$\displaystyle\|\Psi_{n}-{\rm
id}\|_{W^{j,\infty}(\mathcal{O}_{\delta_{h}}(\Omega_{h}^{n}))}\lesssim
h^{q+1-j}\qquad j=0,1.$ (12)

We refer to $q$ as the geometric order of approximation. Such a mapping has
been constructed in [18] based on isoprametric mappings of geometries defined
via level sets. Note that (12) implies

${\rm dist}(\Omega^{n},\Omega_{h}^{n})\lesssim h^{q+1}.$ (13)

###  3.2 Triangulations

We let $\mathcal{T}_{h}$ denote a shape-regular and quasi-uniform simplicial
triangulation of the background domain $\widehat{\Omega}$ with
$h=\max_{T\in\mathcal{T}_{h}}{\rm diam}(T)$. Note the quasi-uniformity implies
a constant $c>0$ such that $h\leq c\,{\rm diam}(T)=:h_{T}$ for all
$T\in\mathcal{T}_{h}$.

We then define, for each time step $n$, the active triangulation and
corresponding domain induced by the background triangulation (cf. Figure 1):

$\displaystyle\mathcal{T}_{h,e}^{n}=\\{T\in\mathcal{T}_{h}:\ {\rm
dist}(\mathbf{x},\Omega_{h}^{n})\leq\delta_{h}\
\exists\,\mathbf{x}\in\bar{T}\\},\qquad\Omega_{h,e}^{n}={\rm
int}\left(\bigcup_{T\in\mathcal{T}_{h,e}^{n}}\bar{T}\right).$

We further definite the set of interior elements for $\Omega_{h}^{n}$ and
associated domain at time step $n$:

$\mathcal{T}_{h,i}^{n}=\\{T\in\mathcal{T}_{h,e}^{n}:\ {\rm
int}(T)\subset\Omega_{h}^{n}\\},\qquad\Omega_{h,i}^{n}={\rm
int}\left(\bigcup_{T\in\mathcal{T}_{h,i}^{n}}\bar{T}\right),$

and denote by $\mathcal{F}_{h,i}^{n}$ (resp., $\mathcal{F}_{h,e}^{n}$) the set
of interior faces of $\mathcal{T}_{h,i}^{n}$ (resp., $\mathcal{T}_{h,e}^{n}$),
i.e.,

$\mathcal{F}_{h,*}^{n}=\\{F=\partial T_{1}\cap\partial T_{2}:\
T_{1},T_{2}\in\mathcal{T}_{h,*}^{n},\ T_{1}\neq T_{2}\\}\qquad*\in\\{i,e\\}.$

We further set $h_{F}={\rm diam}(F)$ for all $F\in\mathcal{F}_{h,e}$.
Following [23, 38], we define the elements in a strip around $\Gamma_{h}^{n}$:

$\mathcal{T}_{\Gamma_{h}}^{n}:=\\{T\in\mathcal{T}_{h,e}^{n}:\ {\rm
dist}(x,\Gamma_{h}^{n})\leq\delta_{h}\ \exists\,x\in\bar{T}\\},$

and define the set of faces in this strip:

$\mathcal{F}_{\Gamma_{h}}^{n}:=\\{F=\partial T_{1}\cap\partial T_{2}:\
T_{1}\in\mathcal{T}_{h,e}^{n},\ T_{2}\in\mathcal{T}_{\Gamma_{h}}^{n},\
T_{1}\neq T_{2},\ |\partial T_{1}\cap\partial T_{2}|>0\\}.$

For any sub-triangulation $\mathcal{S}_{h}\subset\mathcal{T}_{h}$ and
$m\in\mathbb{N}$, we set $H^{m}(\mathcal{S}_{h})$ to be the piecewise Sobolev
space with respect to $\mathcal{S}_{h}$, i.e., $q\in H^{m}(\mathcal{S}_{h})$
implies $q$ is an $L^{2}$ function on the domain induced by $\mathcal{S}_{h}$
and $q|_{T}\in H^{m}(T)$ for all $T\in\mathcal{S}_{h}$. Analogous vector-
valued spaces are denoted in boldface.

Fig. 1: Left: A depiction of the mesh in two dimensions, with the interior
triangulation $\mathcal{T}_{h,i}^{n}$ in red and
$\mathcal{T}_{h,e}^{n}\backslash\mathcal{T}_{h,i}^{n}$ in blue. Right: The
triangulation $\mathcal{T}_{\Gamma_{h}}^{n}$ (blue) around a
$\delta_{h}$-neighborhood tubular region of $\Gamma_{h}$.

###  3.3 Finite Element Spaces and Assumptions

We denote by $\EuScript{P}_{m}(\mathcal{T}_{h})$ the space of piecewise
polynomials of degree $m$ with respect to $\mathcal{T}_{h}$, and set
$\EuScript{P}_{m}^{c}(\mathcal{T}_{h})=\EuScript{P}_{m}(\mathcal{T}_{h})\cap
H^{1}(\widehat{\Omega})$ to be its subspace of continuous, piecewise
polynomials. Analogous vector-valued spaces are denoted in boldface. We
consider a Stokes finite element pair $\mathbf{V}_{h}\times
Q_{h}\subset\mathbf{H}^{1}(\widehat{\Omega})\times L^{2}(\widehat{\Omega})$,
consisting of piecewise polynomial spaces with respect to $\mathcal{T}_{h}$,
and assume the following inclusions

$\displaystyle\bm{\EuScript{P}}_{\underline{m}_{v}}^{c}(\mathcal{T}_{h})$
$\displaystyle\subset\mathbf{V}_{h}\subset\bm{\EuScript{P}}_{\overline{m}_{v}}^{c}(\mathcal{T}_{h}),$
(14)

for some integers $1\leq\underline{m}_{v}\leq\overline{m}_{v}$. We further
assume there exists $m_{q}\in\mathbb{N}_{0}$ such that

$\displaystyle Q_{h}=\EuScript{P}_{m_{q}}(\mathcal{T}_{h})\quad\text{or}\quad
Q_{h}=\EuScript{P}_{m_{q}}^{c}(\mathcal{T}_{h}).$ (15)

We set $\mathbf{V}_{h}^{n}\subset\mathbf{H}^{1}(\Omega_{h,e}^{n})$ to be the
restriction of $\mathbf{V}_{h}$ to $\Omega_{h,e}^{n}$, and let $Q_{h}^{n}$ be
the restriction of $Q_{h}$ to $\Omega_{h,e}^{n}$ with a zero-mean constraint
on $\Omega_{h,i}^{n}$, i.e.,

$Q_{h}^{n}=\\{q|_{\Omega_{h,e}^{n}}:\ \exists\,q\in Q_{h}\text{ such that
}\int_{\Omega_{h,i}^{n}}q\,dx=0\\}.$

Note that, by construction, $\Omega_{h}^{n}\subset\Omega_{h,e}^{n-1}$, and
therefore functions in $\mathbf{V}_{h}^{n-1}\times Q_{h}^{n-1}$ are well
defined on $\Omega_{h}^{n}$.

We define the Nitsche-type norms on
$\mathbf{H}^{1}(\Omega_{h}^{n})\cap\mathbf{H}^{2}(\mathcal{T}_{h,e}^{n})\big{|}_{\Omega_{h}^{n}}$:

$|\\!|\\!|\mathbf{v}|\\!|\\!|_{n}^{2}:=\|\nabla\mathbf{v}\|_{\Omega_{h}^{n}}^{2}+h^{-1}\|\mathbf{v}\|_{\Gamma^{n}_{h}}^{2}+h\|\nabla\mathbf{v}\|_{\Gamma_{h}^{n}}^{2},\qquad$

and further define the norm for piecewise smooth functions on the extended
domains:

$|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}^{2}:=\|\nabla\mathbf{v}\|_{\Omega_{h,e}^{n}}^{2}+|\\!|\\!|\mathbf{v}|\\!|\\!|_{n}^{2}.\qquad$

Likewise, we define the weighted $H^{1}$-seminorm with respect to the interior
mesh $\mathcal{T}_{h,i}^{n}$:

$|\\!|\\!|q|\\!|\\!|_{n,i}^{2}:=\sum_{T\in\mathcal{T}_{h,i}^{n}}h_{T}^{2}\|\nabla
q\|_{T}^{2}+\sum_{F\in\mathcal{F}_{h,i}^{n}}h_{F}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2},$

where $\bigl{[}\bigl{[}\cdot\bigr{]}\bigr{]}$ denotes the jump operator across
an interior face. Note that $|\\!|\\!|\cdot|\\!|\\!|_{n,i}$ is a norm on
$Q_{h}^{n}|_{\Omega_{h,i}^{n}}$. Similarly, we define weighted seminorm over
the extended domain $\Omega_{h,e}^{n}$:

$|\\!|\\!|q|\\!|\\!|_{n,e}^{2}:=\sum_{T\in\mathcal{T}_{h,e}^{n}}h_{T}^{2}\|\nabla
q\|_{T}^{2}+\sum_{F\in\mathcal{F}_{h,e}^{n}}h_{F}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}.$

Note that $|\\!|\\!|q|\\!|\\!|_{n,e}$ is a norm on $Q_{h}^{n}$, and _it will
be our main pressure norm for stability and error analysis._

In addition to the inclusions (14)–(15) , we make the following assumptions to
ensure stability of the discretization presented below.

######  Assumption 3.1.

Assume that, given $q\in Q_{h}^{n}$, there exists
$\mathbf{v}_{h}\in\mathbf{V}_{h}^{n}$ that satisfies

$\displaystyle|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}$
$\displaystyle\lesssim|\\!|\\!|q|\\!|\\!|_{n,i},$ (16a)
$\displaystyle|\\!|\\!|q|\\!|\\!|_{n,i}^{2}$ $\displaystyle\leq
b^{n}_{h}(\mathbf{v},q):=\int_{\Omega_{h}^{n}}({\rm
div}\,\mathbf{v})q\,dx-\int_{\Gamma_{h}^{n}}(\mathbf{v}\cdot\mathbf{n})q\,ds,$
(16b) $\displaystyle\|\mathbf{v}\|_{\Omega_{h}^{n}}$ $\displaystyle\lesssim
h|\\!|\\!|q|\\!|\\!|_{n,i}.$ (16c)

######  Remark 3.1.

The first two statements (16a)–(16b) are assumptions related to discrete inf-
sup stability, but where the $L^{2}$ norm of the pressure function is replaced
with the weighted $H^{1}$-norm. A variation of these conditions is shown to
hold in the context of CutFEM for many standard stable Stokes pairs in [20].
Using a Verfürth-type trick, it is shown in this reference that, if
(16a)–(16b) is satisfied then the discrete inf-sup condition with $L^{2}$
pressure norm holds:

$\theta\|q\|_{\Omega_{h}^{n}}\leq\sup_{\mathbf{v}\in\mathbf{V}_{h}^{n}}\frac{b_{h}^{n}(\mathbf{v},q)}{|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}}+|q|_{J^{n}_{h}}\qquad\forall
q\in{Q}_{h}^{n},$

where $\theta>0$ is independent of $h$ and how $\Gamma_{h}^{n}$ cuts through
the triangulation $\mathcal{T}_{h}$, and $|\cdot|_{J_{h}^{n}}$ is given by
(21) below. We show below in Section 5 that the third condition (16c) is
satisfied for several canonical pairs as well.

######  Remark 3.2.

Assumption 3.1 can be modified and slightly weakened by replacing
$\Omega_{h,i}^{n}$ and $\mathcal{T}_{h,i}^{n}$ by a smaller domain and mesh,
respectively, provided the pressure ghost-penalty compensates for the smaller
domain. In particular, let
$\tilde{\mathcal{T}}_{h,i}^{n}\subset\mathcal{T}_{h,i}^{n}$ be a sub-mesh with
corresponding domain $\tilde{\Omega}_{h,i}^{n}={\rm
int}\left(\bigcup_{T\in\tilde{\mathcal{T}}_{h,i}^{n}}\bar{T}\right)$. Then if

$\|q\|_{\Omega_{h,e}^{n}}\lesssim\|q\|_{\tilde{\Omega}_{h,i}^{n}}+|q|_{J^{n}_{h}}\qquad\forall
q\in Q_{h}^{n},$ (17)

then we can replace $\mathcal{T}_{h,i}^{n}$ by $\tilde{\mathcal{T}}_{h,i}^{n}$
and $\Omega_{h,i}^{n}$ by $\tilde{\Omega}_{h,i}^{n}$ in Assumption 3.1. This
modified assumption is used in the case $\mathbf{V}_{h}\times Q_{h}$ is the
Taylor-Hood pair.

###  3.4 The CutFEM Discretization

The finite element method based on the backward Euler temporal discretization
seeks, at each time step, the pair
$(\mathbf{u}_{h}^{n},p_{h}^{n})\in\mathbf{V}_{h}^{n}\times Q_{h}^{n}$ such
that

$\int_{\Omega_{h}^{n}}\Big{(}\frac{\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}}{\Delta
t}\Big{)}\cdot\mathbf{v}\,dx+a^{n}_{h}(\mathbf{u}_{h}^{n},\mathbf{v})-b_{h}^{n}(\mathbf{v},p_{h}^{n})+b_{h}^{n}(\mathbf{u}_{h}^{n},q)+\gamma_{J}J_{h}^{n}(p_{h}^{n},q)=F^{n}(\mathbf{v},q),$
(18)

for all $\mathbf{v}\in\mathbf{V}_{h}^{n},\,q\in Q_{h}^{n}$. Here,
$b_{h}^{n}(\cdot,\cdot)$ is given by (16b), and the bilinear form
$a_{h}^{n}(\cdot,\cdot)$ is defined as

$\displaystyle\widehat{a}_{h}^{n}(\mathbf{u},\mathbf{v})$
$\displaystyle=\int_{\Omega_{h}^{n}}\nabla\mathbf{u}:\nabla\mathbf{v}\,dx+\int_{\Omega_{h}^{n}}(\mathbf{w}\cdot\nabla\mathbf{u})\cdot\mathbf{v}\,dx-\int_{\Gamma_{h}^{n}}\left([(\nabla\mathbf{u})\mathbf{n}]\cdot\mathbf{v}+[(\nabla\mathbf{v})\mathbf{n}]\cdot\mathbf{u}-\frac{\eta}{h}\mathbf{u}\cdot\mathbf{v}\right)\,ds,$
$\displaystyle a_{h}^{n}(\mathbf{u},\mathbf{v})$
$\displaystyle=\widehat{a}_{h}^{n}(\mathbf{u},\mathbf{v})+\gamma_{s}s^{n}_{h}(\mathbf{u},\mathbf{v}),$

where $\gamma_{s},\gamma_{J},\eta\geq 1$ are user-defined constants. The
bilinear forms $s_{h}^{n}(\cdot,\cdot)$ and $J_{h}^{n}(\cdot,\cdot)$ consist
of ghost-penalty terms acting on $\mathbf{V}^{n}_{h}\times\mathbf{V}_{h}^{n}$
and $Q_{h}^{n}\times Q_{h}^{n}$, respectively, defined on an $O(\delta_{h})$
neighborhood of $\Gamma_{h}^{n}$:

$\begin{split}s_{h}^{n}(\mathbf{u},\mathbf{v})&=\sum_{F\in\mathcal{F}_{\Gamma_{h}}^{n}}\sum_{k=1}^{\overline{m}_{v}}h^{2k-1}\int_{F}\bigl{[}\bigl{[}\partial^{k}_{\mathbf{n}_{F}}\mathbf{u}\bigr{]}\bigr{]}\bigl{[}\bigl{[}\partial^{k}_{\mathbf{n}_{F}}\mathbf{v}\bigr{]}\bigr{]}\,ds,\\\
J_{h}^{n}(p,q)&=\sum_{F\in\mathcal{F}_{\Gamma_{h}}^{n}}\sum_{k=0}^{m_{q}}h^{2k+1}\int_{F}\bigl{[}\bigl{[}\partial^{k}_{\mathbf{n}_{F}}p\bigr{]}\bigr{]}\bigl{[}\bigl{[}\partial^{k}_{\mathbf{n}_{F}}q\bigr{]}\bigr{]}\,ds,\end{split}$
(19)

and $\partial^{k}_{\mathbf{n}_{F}}$ denotes the $k$th-order directional
derivative with respect to the normal of the face $F$. Here,
$\overline{m}_{v}$ and $m_{q}$ are the integers in (14)–(15). Finally,
$F^{n}(\mathbf{v},q)$ is a bounded linear functional on
$\mathbf{V}^{n}_{h}\times Q_{h}^{n}$ with

$\|F^{n}\|_{*}:=\sup_{(\mathbf{v},q)\in\mathbf{V}_{h}^{n}\times
Q_{h}^{n}}\frac{F^{n}(\mathbf{v},q)}{(|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}^{2}+|\\!|\\!|q|\\!|\\!|_{n,e}^{2})^{\frac{1}{2}}}<\infty.$
(20)

In (18) it is given by

$F^{n}(\mathbf{v},q)=\int_{\Omega^{n}_{h}}{\bf f}^{n}\cdot\mathbf{v}\,dx,$

but later we will consider a more general $F^{n}$ for the purpose of analysis.

######  Remark 3.3.

The ghost-penalty bilinear forms (19) both stabilize the solution of problem
(18) due to irregular cuts as well as yield implicit extensions to
$\Omega_{h,e}^{n}$. These terms also aid in algebraic stabilization, as the
resulting condition number of the system is insensitive to how
$\Gamma_{h}^{n}$ intersects $\mathcal{T}_{h}$. The pressure ghost-
stabilization form $J_{h}^{n}(\cdot,\cdot)$ ensures numerical stability as it
provides an inf-sup-type stability condition of the pair
$\mathbf{V}_{h}^{n}\times Q_{h}^{n}$ (cf. Remark 3.1).

There are now several types of ghost-penalty stabilization besides the
“derivative jump version” used in (19). These include the “direct version”
[32] as well as the “local projection stabilization version” [7]. In
principle, we can replace (19) with any choice of these types of ghost penalty
versions, and the stability and convergence analysis presented below carries
through with only superficial modifications. However, for clarity of
presentation, we only focus on the derivative jump version in detail below.

Finally, we remark that the extension of the discrete pressure approximation
to all of $\Omega_{h,e}^{n}$ is not required, in particular, the pressure
ghost-penalty stabilization $J_{h}^{n}(\cdot,\cdot)$ only needs to be defined
on a single layer of elements around $\Gamma_{h}^{n}$ to ensure stability.
However, we use the set of faces $\mathcal{F}_{\Gamma_{h}}^{n}$ for both terms
in (19) to simplify the presentation.

##  4 Stability and Convergence Analysis

We denote by

$|\mathbf{v}|_{s_{h}^{n}}=\sqrt{s_{h}^{n}(\mathbf{v},\mathbf{v})}\quad\text{and}\quad|q|_{J^{n}_{h}}=\sqrt{J_{h}^{n}(q,q)}$
(21)

the semi-norms induced by the bilinear forms $s_{h}^{n}(\cdot,\cdot)$ and
$J_{h}^{n}(\cdot,\cdot)$, respectively. We assume that the Nitsche penalty
parameter $\eta$ is chosen sufficiently large (but independent of $h$ and the
mesh-interface cut) such that $a_{h}^{n}(\cdot,\cdot)$ is coercive on
$\mathbf{V}_{h}^{n}$ (cf. [10]). In particular, we assume $\eta>0$ is chosen
such that

$a_{h}^{n}(\mathbf{v},\mathbf{v})\geq\frac{1}{2}|\\!|\\!|\mathbf{v}|\\!|\\!|_{n}^{2}+\gamma_{s}|\mathbf{v}|_{s_{h}^{n}}^{2}\qquad\forall\mathbf{v}\in\mathbf{V}_{h}^{n}.$
(22)

Similar to [23, 38], we assume that elements in the strip
$\mathcal{T}_{h,e}^{n}\backslash\mathcal{T}_{h,i}^{n}$ can be reached from an
uncut element in $\mathcal{T}_{h,i}^{n}$ by a path that crosses at most $L$
faces with $L\lesssim(1+\frac{\delta_{h}}{h})$; we refer to [23, 38] to see
why this is a reasonable assumption and how it relates to the shape-regularity
of the triangulation $\mathcal{T}_{h}$. We consider the setting where $L$ is
uniformly bounded with respect to the discretization parameters, i.e., when
$\delta_{h}\lesssim h$. Recalling that
$c_{\delta_{h}}\mathbf{w}_{\infty}^{n}\Delta t\leq\delta_{h}$ with $1\leq
c_{\delta_{h}}=O(1)$, this brings us to the time-step restriction:

$\Delta t\lesssim h.$ (23)

The condition (23) and $\|\mathbf{w}\|_{L^{\infty}(\mathcal{Q})}\lesssim 1$
implies

$L\lesssim 1.$ (24)

Thanks to (24) and standard properties of the stabilization terms (see, e.g.,
[26, Lemma 5.1]), we have the following norm equivalences for all
$\mathbf{v}\in\mathbf{V}_{h}^{n}$ and $q\in Q_{h}^{n}$:

$\begin{split}\|\mathbf{v}\|_{\Omega_{h,e}^{n}}^{2}&\simeq\|\mathbf{v}\|_{\Omega_{h}^{n}}^{2}+h^{2}|\mathbf{v}|_{s_{h}^{n}}^{2},\\\
|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}^{2}&\simeq|\\!|\\!|\mathbf{v}|\\!|\\!|_{n}^{2}+|\mathbf{v}|_{s_{h}^{n}}^{2},\\\
|\\!|\\!|q|\\!|\\!|_{n,e}^{2}&\simeq|\\!|\\!|q|\\!|\\!|_{n,i}^{2}+|q|_{J_{h}^{n}}^{2},\\\
\|q\|_{\Omega_{h,e}^{n}}^{2}&\simeq\|q\|_{\Omega_{h,i}^{n}}^{2}+|q|_{J_{h}^{n}}^{2}.\end{split}$
(25)

###  4.1 Preliminary Results

In this section, we collect some preliminary results used in the stability and
the convergence analysis of the finite element method (18).

######  Lemma 1.

For $h$ sufficiently small, there holds for all
$\mathbf{v}\in\mathbf{V}_{h}^{n-1}$,

$\|\mathbf{v}\|_{\Omega_{h}^{n}}^{2}\leq\|\mathbf{v}\|_{\Omega_{h,e}^{n-1}}^{2}\leq(1+c_{1}\Delta
t)\|\mathbf{v}\|^{2}_{\Omega_{h}^{n-1}}+\frac{\Delta
t}{4}|\\!|\\!|\mathbf{v}|\\!|\\!|_{n-1}^{2}+\Delta
tL|\mathbf{v}|_{s_{h}^{n-1}}^{2}$ (26)

for a constant $c_{1}>0$ independent of $h$, $\Delta t$ and how the boundary
cuts through the triangulation.

###### Proof.

From [23, Lemma 5.7], we have

$\displaystyle\|\mathbf{v}\|_{\Omega_{h,e}^{n-1}}^{2}$
$\displaystyle\leq(1+c_{1}(\epsilon)\Delta
t)\|\mathbf{v}\|_{\Omega_{h}^{n-1}}^{2}+c_{2}(\epsilon)\Delta
t\|\nabla\mathbf{v}\|_{\Omega_{h}^{n-1}}^{2}+c_{3}(\epsilon,h)\Delta
tL|\mathbf{v}|_{s_{h}^{n-1}}^{2},$

with

$\displaystyle
c_{1}(\epsilon)=c^{\prime}c_{\delta_{h}}\mathbf{w}_{\infty}^{n}(1+\epsilon^{-1}),\qquad$
$\displaystyle
c_{2}(\epsilon)=c^{\prime}c_{\delta_{h}}\mathbf{w}_{\infty}^{n}\epsilon,$
$\displaystyle c_{3}(\epsilon,h)=c_{2}(\epsilon)+c_{4}(\epsilon,h),\qquad$
$\displaystyle
c_{4}(\epsilon,h)=h^{2}c^{\prime}c_{\delta_{h}}\mathbf{w}_{\infty}^{n}(1+\epsilon^{-1}),$

$c^{\prime}>0$ is a generic constant, and $\epsilon>0$ is arbitrary. The
result (26) follows from the inequality
$\|\nabla\mathbf{v}\|_{\Omega_{h}^{n-1}}\leq|\\!|\\!|\mathbf{v}|\\!|\\!|_{n-1}$
and by taking $\epsilon$ such that $c_{2}(\epsilon)=\frac{1}{4}$ and $h$
sufficiently small such that $c_{4}(\epsilon,h)\leq 1$. ∎

######  Lemma 2.

There holds the following discrete Poincare inequality

$\|\mathbf{v}\|_{\Omega_{h}^{n}}\leq
c_{P}|\\!|\\!|\mathbf{v}|\\!|\\!|_{n}\qquad\forall\mathbf{v}\in\mathbf{V}_{h}^{n}.$
(27)

###### Proof.

See [26, Lemma 7.2]. ∎

The following continuity estimate for the bilinear form
$a^{n}_{h}(\cdot,\cdot)$ is essentially given in [10] (also see [26, 38]) and
follows from the Cauchy-Schwarz inequality, so the proof is omitted.

######  Lemma 3.

There holds

$\displaystyle a^{n}_{h}(\mathbf{u},\mathbf{v})$
$\displaystyle\lesssim|\\!|\\!|\mathbf{u}|\\!|\\!|_{n}|\\!|\\!|\mathbf{v}|\\!|\\!|_{n}+\gamma_{s}|\mathbf{u}|_{s_{h}^{n}}|\mathbf{v}|_{s_{h}^{n}}\quad$
$\displaystyle\forall\mathbf{u},\mathbf{v}\in\mathbf{H}^{\overline{m}_{v}+1}(\mathcal{T}_{h,e}^{n})\cap\mathbf{H}^{1}(\Omega_{h,e}^{n}).$
(28)

The next result states a discrete trace inequality for discontinuous piecewise
polynomial functions. Its proof is given in Appendix A (also see [6, Theorem
4.4])).

######  Lemma 4.

There holds

$\|q\|_{\Gamma_{h}^{n}}\lesssim h^{-1}|\\!|\\!|q|\\!|\\!|_{n,e}\qquad\forall
q\in Q_{h}^{n}.$ (29)

###  4.2 Stability analysis

In this section, we derive stability results for the finite element method
(18). First, we state the energy estimate for the finite element velocity in
the following lemma. This result is essentially given in [38, Theorem 5.9],
but we provide a proof for completeness.

######  Lemma 5.

There holds for $h$ sufficiently small, any $\varepsilon>0$ and $k=1,2,\dots$

$\begin{split}&\|\mathbf{u}_{h}^{k}\|^{2}_{\Omega^{k}_{h}}+\sum_{n=1}^{k}\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|_{\Omega^{n}_{h}}^{2}+\Delta
t\sum_{n=1}^{k}\left(\frac{1}{4}|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+(2\gamma_{s}-L-\frac{1}{2})|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2}+2\gamma_{J}|p_{h}^{n}|_{J_{h}^{n}}^{2}\right)\\\
&\hskip 54.2025pt\leq{\rm
exp}(ct_{k})\left(\|\mathbf{u}_{h}^{0}\|_{\Omega^{0}_{h}}^{2}+\frac{\Delta
t}{4}|\\!|\\!|\mathbf{u}^{0}_{h}|\\!|\\!|_{0}^{2}+\Delta
tL|\mathbf{u}^{0}_{h}|_{s_{h}^{0}}^{2}\right.\\\ &\hskip
126.47249pt\left.+\Delta
t(c_{e}+\varepsilon^{-1})\sum_{n=0}^{k}\|F^{n}\|_{*}^{2}+\Delta
t{\varepsilon}\sum_{n=0}^{k}|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2}\right).\end{split}$
(30)

with constants $c$ and $c_{e}$ independent of the discretization parameters.

###### Proof.

Taking $\mathbf{v}=\mathbf{u}_{h}^{n-1}$ and $q=p_{h}^{n}$ in (18), adding the
two statements, applying (22), and using the algebraic identity
$(a-b)a=\frac{1}{2}(a^{2}-b^{2})+\frac{1}{2}(a-b)^{2}$ yields

$\frac{1}{2}\|\mathbf{u}_{h}^{n}\|_{\Omega_{h}^{n}}^{2}-\frac{1}{2}\|\mathbf{u}_{h}^{n-1}\|_{\Omega_{h}^{n}}^{2}+\frac{1}{2}\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|_{\Omega_{h}^{n}}^{2}+\Delta
t\left(\frac{1}{2}|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+\gamma_{s}|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2}+\gamma_{J}|p_{h}^{n}|_{J^{n}_{h}}^{2}\right)\leq\Delta
tF^{n}(\mathbf{u}_{h}^{n},p_{h}^{n})$

Using (20) and the Cauchy-Schwarz inequality we estimate the right-hand side
as follows

$\begin{split}F^{n}(\mathbf{u}_{h}^{n},p_{h}^{n})&\leq\|F^{n}\|_{*}(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n,e}^{2}+|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2})^{\frac{1}{2}}\leq\|F^{n}\|_{*}(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n,e}+|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e})\\\
&\leq\sqrt{c_{e}/2}\|F^{n}\|_{*}(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2})^{\frac{1}{2}}+\|F^{n}\|_{*}|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}\\\
&\leq\frac{1}{2}\Big{(}c_{e}+\epsilon^{-1}\Big{)}\|F^{n}\|_{*}^{2}+\frac{1}{4}(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2})+{\frac{\varepsilon}{2}}|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2},\end{split}$

where $c_{e}\geq 1$ satisfies
$|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n,e}^{2}\leq\frac{c_{e}}{2}(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2})$
(cf. (25)). This yields

$\|\mathbf{u}_{h}^{n}\|_{\Omega_{h}^{n}}^{2}-\|\mathbf{u}_{h}^{n-1}\|_{\Omega_{h}^{n}}^{2}+\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|_{\Omega_{h}^{n}}^{2}+{\Delta
t}\left(\frac{1}{2}|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+(2\gamma_{s}-\frac{1}{2})|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2}+2\gamma_{J}|p_{h}^{n}|_{J^{n}_{h}}^{2}\right)\\\
\leq{\Delta
t}\Big{(}\big{(}c_{e}+\varepsilon^{-1}\big{)}\|F^{n}\|_{*}^{2}+\varepsilon|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2}\Big{)}.$
(31)

Applying the estimate (26) (with $\mathbf{v}=\mathbf{u}_{h}^{n-1}$) into (31)
and summing the result over $n=1,\ldots,k$ yields

$\|\mathbf{u}_{h}^{k}\|_{\Omega_{h}^{k}}^{2}+\sum_{n=1}^{k}\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|_{\Omega_{h}^{n}}^{2}+{\Delta
t}\sum_{n=1}^{k}\left(\frac{1}{4}|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+(2\gamma_{s}-L-\frac{1}{2})|\mathbf{u}^{n}_{h}|_{s_{h}^{n}}^{2}+2\gamma_{J}|p_{h}^{n}|_{J_{h}^{n}}^{2}\right)\\\
\leq\|\mathbf{u}_{h}^{0}\|_{\Omega_{h}^{0}}+\frac{\Delta
t}{4}|\\!|\\!|\mathbf{u}_{h}^{0}|\\!|\\!|_{0}^{2}+\Delta
tL|\mathbf{u}_{h}^{0}|_{s_{h}^{0}}^{2}+c_{1}\Delta
t\sum_{n=0}^{k-1}\|\mathbf{u}_{h}^{n}\|_{\Omega_{h}^{n}}^{2}\\\ +\Delta
t\sum_{n=1}^{k}\Big{(}\big{(}c_{e}+\varepsilon^{-1}\big{)}\|F^{n}\|_{*}^{2}+\varepsilon|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2}\Big{)}.$

The estimate (30) now follows from a discrete Gronwall inequality. ∎

For the complete stability result we need to estimate the pressure term on the
right-hand side of (30). The estimate is given in the next lemma.

######  Lemma 6.

Assume $h^{2}\lesssim\Delta t$. Then

$\Delta t\sum_{n=1}^{k}|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2}\lesssim{\rm
exp}(ct_{k})\left(\|\mathbf{u}_{h}^{0}\|_{\Omega^{0}_{h}}^{2}+{\Delta
t}(|\\!|\\!|\mathbf{u}^{0}_{h}|\\!|\\!|_{0}^{2}+|\mathbf{u}^{0}_{h}|_{s_{h}^{0}}^{2})+\Delta
t\sum_{n=0}^{k}\|F^{n}\|_{*}^{2}\right).$

###### Proof.

Let $\mathbf{v}\in\mathbf{V}_{h}^{n}$ satisfy (16) with $q=p_{h}^{n}$. Then
using the identity (18) and bounds in (16), (20), (28), and (25), we have

$\displaystyle|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,i}^{2}$ $\displaystyle\leq
b^{n}_{h}(\mathbf{v},p_{h}^{n})$
$\displaystyle=\int_{\Omega^{n}_{h}}\frac{\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}}{\Delta
t}\cdot\mathbf{v}\,dx+a_{h}^{n}(\mathbf{u}^{n}_{h},\mathbf{v})-F^{n}(\mathbf{v},0)$
$\displaystyle\lesssim\left(\frac{1}{\Delta
t}\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|_{\Omega_{h}^{n}}\|\mathbf{v}\|_{\Omega_{h}^{n}}+|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}|\\!|\\!|\mathbf{v}|\\!|\\!|_{n}+\gamma_{s}|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}|\mathbf{v}|_{s_{h}^{n}}+\|F^{n}\|_{*}|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}\right)$
$\displaystyle\lesssim\left(\frac{h}{\Delta
t}\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|_{\Omega^{n}_{h}}+|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}+\gamma_{s}|\mathbf{u}_{h}^{n}|^{2}_{s_{h}^{n}}+\|F^{n}\|_{*}\right)|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,i}.$

Thus, we have

$\displaystyle\Delta t|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,i}^{2}$
$\displaystyle\lesssim\frac{h^{2}}{\Delta
t}\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|^{2}_{\Omega^{n}_{h}}+\Delta
t(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2}+\|F^{n}\|_{*}^{2}).$
(32)

Combining this with (25) leads to the estimate of the pressure norm in the
extended domain:

$\displaystyle\Delta t|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2}$
$\displaystyle\lesssim\frac{h^{2}}{\Delta
t}\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|^{2}_{\Omega^{n}_{h}}+\Delta
t(|p_{h}^{n}|_{J_{h}^{n}}^{2}+|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2}+\|F^{n}\|_{\ast}^{2}).$

Summing inequality over $n=1,\ldots,k$, and using $h^{2}\lesssim\Delta t$ and
(30) gets

$\displaystyle\Delta t\sum_{n=1}^{k}|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2}$
$\displaystyle\lesssim\sum_{n=1}^{k}\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|^{2}_{\Omega^{n}_{h}}+\Delta
t\sum_{n=1}^{k}(|p_{h}^{n}|_{J_{h}^{n}}^{2}+|\\!|\\!|\mathbf{u}^{n}_{h}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2})+\Delta
t\sum_{n=1}^{k}\|F^{n}\|_{*}^{2}.$

All terms on the right-hand side of the last inequality are estimated in (30).
Thus, by applying (30) with $\varepsilon$ sufficiently small but independent
of the discretization parameters proves the lemma. ∎

######  Remark 4.1.

The corresponding BDF2 scheme is analogous to (18), but where the discrete
time derivative is replaced by
$\frac{3\mathbf{u}^{n+1}_{h}-4\mathbf{u}_{h}^{n-1}+\mathbf{u}_{h}^{n-2}}{2\Delta
t}$, and the computational mesh is enlarged. In particular, $\delta_{h}$ is
replaced by $2\delta_{h}$ in the definition of $\mathcal{T}_{h,e}^{n}$ so that
functions in $\mathbf{V}_{h}^{n-2}$ are well defined in $\Omega_{h}^{n}$. In
this setting, a stability result holds for the discrete velocity similar to
(30), but where
$\|\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}\|_{\Omega_{h}^{n}}^{2}$ is replaced
by
$\|\mathbf{u}_{h}^{n}-2\mathbf{u}_{h}^{n-1}+\mathbf{u}_{h}^{n-2}\|_{\Omega_{h}^{n}}^{2}$.
The proof of this result is similar to that of Lemma 5 but using a different
polarization identity, and so we omit the details.

The stability of the discrete pressure solution in the BDF2 scheme is more
subtle and requires a different argument than Lemma 6. Analogous to (32),
there holds

$\displaystyle\Delta t|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,i}^{2}$
$\displaystyle\lesssim\frac{h^{2}}{\Delta
t}\|3\mathbf{u}_{h}^{n}-4\mathbf{u}_{h}^{n-1}+\mathbf{u}_{h}^{n-2}\|^{2}_{\Omega^{n}_{h}}+\Delta
t(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2}+\|F^{n}\|_{*}^{2}),$

and therefore

$\displaystyle\sum_{n=1}^{k}\Delta t|\\!|\\!|p_{h}^{n}|\\!|\\!|_{n,e}^{2}$
$\displaystyle\lesssim\frac{h^{2}}{\Delta
t}\sum_{n=0}^{k}\|\mathbf{u}_{h}^{n}\|^{2}_{\Omega^{n}_{h}}+\Delta
t\sum_{k=1}^{n}(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2}+\Delta
t|p_{h}^{n}|_{J_{h}^{n}}^{2})+\Delta t\sum_{n=1}^{k}\|F^{n}\|_{*}^{2}$
$\displaystyle\leq\frac{Th^{2}}{\Delta t^{2}}\max_{0\leq n\leq
N}\|\mathbf{u}_{h}^{n}\|^{2}_{\Omega^{n}_{h}}+\Delta
t\sum_{k=1}^{n}(|\\!|\\!|\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\mathbf{u}_{h}^{n}|_{s_{h}^{n}}^{2}+\Delta
t|p_{h}^{n}|_{J_{h}^{n}}^{2})+\Delta t\sum_{n=1}^{k}\|F^{n}\|_{*}^{2}.$

Thus, for $h\lesssim\Delta t$, the terms in the right-hand side of this
expression are uniformly bounded, hence obtaining a stability estimate for the
discrete pressure solution. Note that when combined with (23), we have the
relation $\Delta t\simeq h$ in the case of BDF2.

###  4.3 Consistency

The consistency of the scheme (18) largely follows the arguments in [38, Lemma
5.14]. First, we identify the extensions of the smooth exact solution
$\mathcal{E}\mathbf{u}$ and $\mathcal{E}p$ with $\mathbf{u}$ and $p$,
respectively, both of which satisfy (9). Recall that for $\mathbf{u}$, we
consider the divergence-free extension from (8). We then set
$\mathbb{U}^{n}=\mathbf{u}^{n}-\mathbf{u}_{h}^{n}$ and
$\mathbb{P}^{n}=p^{n}-p_{h}^{n}$ to denote the errors at $t_{n}$.

######  Lemma 7.

There holds for all $(\mathbf{v},q)\in\mathbf{V}^{n}_{h}\times Q_{h}^{n}$,

$\displaystyle\int_{\Omega_{h}^{n}}\frac{\mathbb{U}^{n}-\mathbb{U}^{n-1}}{\Delta
t}\cdot\mathbf{v}\,dx+a_{h}^{n}(\mathbb{U}^{n},\mathbf{v})-b_{h}^{n}(\mathbf{v},\mathbb{P}^{n})+b_{h}^{n}(\mathbb{U}^{n},q)+\gamma_{J}J_{h}^{n}(\mathbb{P}^{n},q)=\mathfrak{C}^{n}_{c}(\mathbf{v},q),$

where the consistency error $\mathfrak{C}^{n}_{c}(\mathbf{v},q)$ satisfies

$\displaystyle|\mathfrak{C}^{n}_{c}(\mathbf{v},q)|\lesssim
h^{q}\|\mathbf{u}^{n}\|_{H^{2}(\Omega^{n})}|\\!|\\!|q|\\!|\\!|_{n,e}$ (33)
$\displaystyle\quad+(\Delta t+h^{q}+h^{m_{1}}+h^{m_{2}})\left(\|{\bf
f}^{n}\|_{H^{1}(\Omega^{n})}+\|\mathbf{u}\|_{W^{3,5}(\mathcal{Q})}+\|p^{n}\|_{H^{m_{2}}(\Omega^{n})}+\|\mathbf{u}^{n}\|_{H^{m_{1}+1}(\Omega^{n})}\right)|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e},$

for any integers $m_{1},m_{2}$ satisfying $m_{1}\geq\overline{m}_{v}$ and
$m_{2}\geq m_{q}+1$.

###### Proof.

Recall that
$\Psi_{n}:\mathcal{O}_{\delta_{h}}(\Omega_{h}^{n})\to\mathcal{O}_{\delta_{h}}(\Omega_{h})$
is the mapping that connects the approximate and exact domains and satisfies
(12). Testing (2) with $\mathbf{v}^{\ell}:=\mathbf{v}\circ\Psi^{-1}_{n}$,
$\mathbf{v}\in\mathbf{V}_{h}^{n}$, and $q^{\ell}:=q\circ\Psi^{-1}_{n}$, $q\in
Q_{h}^{n}$, and integrating by parts we arrive at the identity

$\int_{\Omega^{n}}\frac{\partial\mathbf{u}^{n}}{\partial
t}\mathbf{v}^{\ell}\,dx+\int_{\Omega^{n}}\nabla\mathbf{u}^{n}:\nabla\mathbf{v}^{\ell}\,dx-\int_{\Gamma^{n}}[(\nabla\mathbf{u}^{n})\mathbf{n}]\cdot\mathbf{v}^{\ell}\,ds+\int_{\Omega^{n}}(\mathbf{w}\cdot\nabla\mathbf{u}^{n})\cdot\mathbf{v}^{\ell})\,dx\\\
-\int_{\Omega^{n}}p^{n}\operatorname{\rm
div}\mathbf{v}^{\ell}\,dx+\int_{\Gamma^{n}}p^{n}(\mathbf{v}^{\ell}\cdot\mathbf{n})\,ds-\int_{\Omega^{n}}q^{\ell}\operatorname{\rm
div}\mathbf{u}^{n}\,dx=\int_{\Omega^{n}}{\bf
f}^{n}\cdot\mathbf{v}^{\ell}\,dx.$

Subtracting this identity from (18) gives the consistency term:

$\displaystyle\mathfrak{C}^{n}_{c}(\mathbf{v},q)$
$\displaystyle=\underbrace{\int_{\Omega^{n}}{\bf
f}^{n}\cdot\mathbf{v}^{\ell}\,dx-\int_{\Omega^{n}_{h}}{\bf
f}^{n}\cdot\mathbf{v}\,dx}_{=:\mathfrak{T}_{1}}+\underbrace{\int_{\Omega_{h}^{n}}\frac{{\mathbf{u}}^{n}-{\mathbf{u}}^{n-1}}{\Delta
t}\cdot\mathbf{v}\,dx-\int_{\Omega^{n}}\frac{\partial\mathbf{u}^{n}}{\partial
t}\mathbf{v}^{\ell}\,dx}_{=:\mathfrak{T}_{2}}$
$\displaystyle+\underbrace{\widehat{a}_{h}^{n}(\mathbf{u}^{n},\mathbf{v})-\int_{\Omega^{n}}\nabla\mathbf{u}^{n}:\nabla\mathbf{v}^{\ell}\,dx+\int_{\Gamma^{n}}[(\nabla\mathbf{u}^{n})\mathbf{n}]\cdot\mathbf{v}^{\ell}\,ds-\int_{\Omega^{n}}(\mathbf{w}\cdot\nabla\mathbf{u}^{n})\cdot\mathbf{v}^{\ell}\,dx}_{=:\mathfrak{T}_{3}}$
$\displaystyle+\underbrace{\int_{\Omega^{n}}p^{n}\operatorname{\rm
div}\mathbf{v}^{\ell}\,dx-\int_{\Gamma^{n}}p^{n}(\mathbf{v}^{\ell}\cdot\mathbf{n})\,ds-
b_{h}^{n}(\mathbf{v},p^{n})}_{=:\mathfrak{T}_{4}}$
$\displaystyle-\underbrace{b_{h}^{n}(\mathbf{u}^{n},q)}_{=:\mathfrak{T}_{5}}+\underbrace{\gamma_{J}J_{h}^{n}(p^{n},q)+\gamma_{s}s_{h}^{n}(\mathbf{u}^{n},\mathbf{v})}_{=:\mathfrak{T}_{6}}.$

Estimates for $\mathfrak{T}_{1}$ and $\mathfrak{T}_{4}$ are exactly the same
as in [38, Lemma 5.14]:

$\displaystyle|\mathfrak{T}_{1}|$ $\displaystyle\lesssim h^{q}\|{\bf
f}^{n}\|_{H^{1}(\Omega^{n})}\|\mathbf{v}\|_{\Omega_{h}^{n}},\qquad$
$\displaystyle|\mathfrak{T}_{4}|\lesssim(h^{q}\|p^{n}\|_{H^{1}(\Omega^{n})}+h^{m_{2}}\|p^{n}\|_{H^{m_{2}}(\Omega^{n})})|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}$
(34)

for any $m_{2}\geq 1$. Likewise, the arguments in [23, Lemma 5.6] and [38,
Lemma 5.14] show

$|\mathfrak{T}_{2}|\lesssim(\Delta
t+h^{q})\|\mathbf{u}\|_{W^{2,\infty}(\widehat{\mathcal{Q}})}\|\mathbf{v}\|_{\Omega_{h}^{n}}\lesssim(\Delta
t+h^{q})\|\mathbf{u}\|_{W^{3,5}(\mathcal{Q})}\|\mathbf{v}\|_{\Omega_{h}^{n}},$
(35)

where we used (9c) in the last inequality.

Unlike the problem considered in [38], the bilinear form
$\widehat{a}_{h}^{n}(\cdot,\cdot)$ includes convective terms. Nonetheless, the
same arguments in [38, Lemma 5.14] are valid, yielding the following estimate:

$\begin{split}|\mathfrak{T}_{3}|&\lesssim(h^{q}\|\mathbf{u}\|_{W^{2,\infty}(\widehat{\mathcal{Q}})}+h^{m_{1}}\|\mathbf{u}^{n}\|_{H^{m_{1}+1}(\Omega^{n})})|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}\\\
&\lesssim(h^{q}\|\mathbf{u}\|_{W^{3,5}(\mathcal{Q})}+h^{m_{1}}\|\mathbf{u}^{n}\|_{H^{m_{1}+1}(\Omega^{n})})|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e},\end{split}$
(36)

where $m_{1}\geq 1$ is only dictated by the regularity of $\mathbf{u}^{n}$,
and we have again used (9c).

On the other hand, the estimate of
$\mathfrak{T}_{5}=b_{h}^{n}(\mathbf{u}^{n},q)$ should involve the elementwise
scaled $H^{1}$-norm for the pressure (which is nonstandard and not provided in
[38]). Since the extension of $\mathbf{u}$ is divergence–free, the estimate of
$\mathfrak{T}_{5}$ reduces to estimating the boundary term:

$\displaystyle\mathfrak{T}_{5}$
$\displaystyle=-\int_{\Gamma_{h}^{n}}(\mathbf{u}^{n}\cdot\mathbf{n})q\,ds.$

Since $\Psi_{n}(\Gamma_{h}^{n})=\Gamma^{n}$, there holds

$\mathbf{u}^{n}\circ\Psi_{n}=0\qquad\text{on $\Gamma_{h}^{n}$}.$

Using the estimate
$\|\mathbf{u}^{n}-\mathbf{u}^{n}\circ\Psi_{n}\|_{\Gamma_{h}^{n}}\lesssim
h^{q+1}\|\mathbf{u}^{n}\|_{H^{2}(\Omega^{n})}$ (cf. [18, Lemma 7.3]) and the
discrete trace inequality in Lemma 4, we have

$\begin{split}|\mathfrak{T}_{5}|&=\left|\int_{\Gamma_{h}^{n}}(\mathbf{u}^{n}-\mathbf{u}^{n}\circ\Psi_{n})\cdot\mathbf{n}\,q\,ds\right|\leq\|\mathbf{u}^{n}-\mathbf{u}^{n}\circ\Psi_{n}\|_{\Gamma^{n}_{h}}\|q\|_{\Gamma^{n}_{h}}\lesssim
h^{q+1}\|\mathbf{u}^{n}\|_{H^{2}(\Omega^{n})}\|q\|_{\Gamma^{n}_{h}}\\\
&\lesssim
h^{q}\|\mathbf{u}^{n}\|_{H^{2}(\Omega^{n})}|\\!|\\!|q|\\!|\\!|_{n,e}.\end{split}$
(37)

Finally, the consistency term involving ghost stabilization $\mathfrak{T}_{6}$
vanishes provided
$\mathbf{u}^{n}\in\mathbf{H}^{\overline{m}_{v}+1}(\Omega_{h,e}^{n})$ and
$p^{n}\in H^{m_{q}+1}(\Omega_{h,e}^{n})$. The estimate (33) then follows from
(34)–(37) and the discrete Poincare inequality (27). ∎

###  4.4 Error Estimates

In this section, we combine the stability and consistency estimates to obtain
error estimates for the finite element method (18). As a first step, let
$(\mathbf{u}^{n}_{I},p^{n}_{I})\in\mathbf{V}_{h}^{n}\times Q_{h}^{n}$ be
approximations to the exact solution satisfying

$\begin{split}|\\!|\\!|\mathbf{u}^{n}-\mathbf{u}^{n}_{I}|\\!|\\!|_{n,e}+|\mathbf{u}^{n}-\mathbf{u}^{n}_{I}|_{s_{h}^{n}}&\lesssim
h^{\underline{m}_{v}}\|\mathbf{u}^{n}\|_{H^{\overline{m}_{v}+1}(\Omega_{h,e}^{n})}\lesssim
h^{\underline{m}_{v}}\|\mathbf{u}^{n}\|_{H^{\overline{m}_{v}+1}(\Omega^{n})},\\\
|\\!|\\!|p^{n}-p^{n}_{I}|\\!|\\!|_{n,e}+|p^{n}-p^{n}_{I}|_{J_{h}^{n}}&\lesssim
h^{{m}_{q}+1}\|p^{n}\|_{H^{{m}_{q}+1}(\Omega_{h,e}^{n})}\lesssim
h^{{m}_{q}+1}\|p^{n}\|_{H^{{m}_{q}+1}(\Omega^{n})},\end{split}$ (38a) and
$\begin{split}h^{-1}\|\mathbf{u}^{n}-\mathbf{u}^{n}_{I}\|_{\Omega_{h}^{n}}&\lesssim
h^{\underline{m}_{v}}\|\mathbf{u}\|_{H^{\overline{m}_{v}+1}(\Omega_{h,e}^{n})}\lesssim
h^{\underline{m}_{v}}\|\mathbf{u}\|_{H^{\overline{m}_{v}+1}(\Omega^{n})},\\\
\|p^{n}-p^{n}_{I}\|_{\Omega_{h}^{n}}+h^{1/2}\|p^{n}-p_{I}^{n}\|_{\Gamma_{h}^{n}}&\lesssim
h^{{m}_{q}+1}\|p^{n}\|_{H^{{m}_{q}+1}(\Omega_{h,e}^{n})}\lesssim
h^{{m}_{q}+1}\|p^{n}\|_{H^{{m}_{q}+1}(\Omega^{n})}.\end{split}$ (38b)

The existence of such $\mathbf{u}_{I}^{n}$ and $p_{I}^{n}$ satisfying (38)
follows from the inclusions (14)–(15) and standard scaling and interpolation
arguments. We also assume the initial condition of the finite element method
(18) is $\mathbf{u}_{h}^{0}=\mathbf{u}_{I}^{0}$.

We then split the error into its interpolation and discretization errors:

$\displaystyle\mathbb{U}^{n}$
$\displaystyle=\underbrace{(\mathbf{u}^{n}-\mathbf{u}^{n}_{I})}_{=:{\bm{\eta}^{n}}}+\underbrace{(\mathbf{u}_{I}^{n}-\mathbf{u}_{h}^{n})}_{=:\mathbf{e}^{n}_{h}\in\mathbf{V}_{h}^{n}},\qquad\mathbb{P}^{n}=\underbrace{(p^{n}-p^{n}_{I})}_{=:\zeta^{n}}+\underbrace{(p_{I}^{n}-p_{h}^{n})}_{=:{\rm
d}_{h}^{n}\in Q_{h}^{n}}.$

Then the pair $(\mathbf{e}^{n}_{h},{\rm d}_{h}^{n})\in\mathbf{V}_{h}^{n}\times
Q_{h}^{n}$ satisfies

$\displaystyle\int_{\Omega_{h}^{n}}\frac{\mathbf{e}^{n}_{h}-\mathbf{e}_{h}^{n-1}}{\Delta
t}\cdot\mathbf{v}\,dx+a_{h}^{n}(\mathbf{e}_{h}^{n},\mathbf{v})-b_{h}^{n}(\mathbf{v},{\rm
d}_{h}^{n})+b_{h}^{n}(\mathbf{e}_{h}^{n},q)+\gamma_{J}J_{h}({\rm
d}_{h}^{n},q)=\mathfrak{C}^{n}_{c}(\mathbf{v},q)+\mathfrak{C}^{n}_{I}(\mathbf{v},q),$
(39)

for all $(\mathbf{v},q)\in\mathbf{V}_{h}^{n}\times Q_{h}^{n}$, where
$\mathfrak{C}^{n}_{c}(\mathbf{v},q)$ is given in Lemma 7 and

$\displaystyle\mathfrak{C}^{n}_{I}(\mathbf{v},q)$
$\displaystyle=-\underbrace{\int_{\Omega_{h}^{n}}\frac{{\bm{\eta}}^{n}-{\bm{\eta}}^{n-1}}{\Delta
t}\cdot\mathbf{v}\,dx}_{=:\mathfrak{T}_{7}}-\underbrace{a_{h}^{n}({\bm{\eta}}^{n},\mathbf{v})}_{=:\mathfrak{T}_{8}}+\underbrace{b_{h}^{n}(\mathbf{v},\zeta^{n})}_{=:\mathfrak{T}_{9}}-\underbrace{b_{h}^{n}({\bm{\eta}}^{n},q)}_{=:\mathfrak{T}_{10}}-\underbrace{\gamma_{J}J_{h}^{n}(\zeta^{n},q)}_{=:\mathfrak{T}_{11}}.$

We now bound the terms in $\mathfrak{C}^{n}_{I}(\mathbf{v},q)$ individually.

First, by continuity estimates and the approximation properties (38), we have

$\displaystyle|\mathfrak{T}_{i}|\lesssim(h^{\underline{m}_{v}}\|\mathbf{u}^{n}\|_{H^{\overline{m}_{v}+1}(\Omega^{n})}+h^{m_{q}+1}\|p^{n}\|_{H^{m_{q}+1}(\Omega^{n})})|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}\quad
i=8,9,11.$ (40)

For the temporal interpolation error there holds by [23, Lemma 5.7] and the
discrete Poincare inequality (27),

$\begin{split}|\mathfrak{T}_{7}|&\lesssim
h^{\underline{m}_{v}}\sup_{t\in[0,T]}\left(\|\mathbf{u}\|_{H^{\overline{m}_{v}+1}(\Omega(t))}+\|\mathbf{u}_{t}\|_{H^{\overline{m}_{v}}(\Omega(t))}\right)\|\mathbf{v}\|_{\Omega_{h}^{n}}\\\
&\lesssim
h^{\underline{m}_{v}}\sup_{t\in[0,T]}\left(\|\mathbf{u}\|_{H^{\overline{m}_{v}+1}(\Omega(t))}+\|\mathbf{u}_{t}\|_{H^{\overline{m}_{v}}(\Omega(t))}\right)|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}.\end{split}$
(41)

For $\mathfrak{T}_{10}$, we integrate by parts to obtain

$\mathfrak{T}_{10}=\int_{\Omega_{h}^{n}}(\operatorname{\rm
div}{\bm{\eta}}^{n})q\,dx-\int_{\Gamma_{h}^{n}}({\bm{\eta}}^{n}\cdot\mathbf{n})q\,ds=\int_{\Omega_{h}^{n}}{\bm{\eta}}^{n}\cdot\nabla
q\,dx+\sum_{F\in\mathcal{F}_{h,e}^{n}}\int_{F\cap\Omega_{h}^{n}}{\bm{\eta}}^{n}\cdot\mathbf{n}\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\,ds.$

Consequently by an elementwise trace inequality and (38), there holds

$\begin{split}|\mathfrak{T}_{10}|&\leq\left(\sum_{T\in\mathcal{T}_{h,e}^{n}}h_{T}^{-2}\|{\bm{\eta}}^{n}\|_{T}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{h,e}^{n}}h_{T}^{2}\|\nabla
q\|_{T}^{2}\right)^{\frac{1}{2}}\\\
&\qquad+\left(\sum_{F\in\mathcal{F}_{h,e}^{n}}h_{F}^{-1}\|{\bm{\eta}}^{n}\|_{F}^{2}\right)^{\frac{1}{2}}\left(\sum_{F\in\mathcal{F}_{h,e}^{n}}h_{F}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}\right)^{\frac{1}{2}}\\\
&\lesssim\left(\sum_{T\in\mathcal{T}_{h,e}^{n}}(h_{T}^{-2}\|{\bm{\eta}}^{n}\|_{T}^{2}+\|\nabla{\bm{\eta}}^{n}\|_{T}^{2})\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{h,e}^{n}}h_{T}^{2}\|\nabla
q\|_{T}^{2}+\sum_{F\in\mathcal{F}_{h,e}^{n}}h_{F}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}\right)^{\frac{1}{2}}\\\
&\lesssim
h^{\underline{m}_{v}}\|\mathbf{u}^{n}\|_{H^{\overline{m}_{v}+1}(\Omega^{n})}|\\!|\\!|q|\\!|\\!|_{n,e}.\end{split}$
(42)

Summarizing (40)–(42), we proved the bound

$\begin{split}|\mathfrak{C}_{I}^{n}(\mathbf{v},q)|\lesssim\Big{(}h^{\underline{m}_{v}}\sup_{t\in[0,T]}(\|\mathbf{u}\|_{H^{\overline{m}_{v}+1}(\Omega(t))}+\|\mathbf{u}_{t}\|_{H^{\overline{m}_{v}}(\Omega(t))})+h^{m_{q}+1}&\|p^{n}\|_{H^{m_{q}+1}(\Omega^{n})}\Big{)}\\\
&\times\big{(}|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}+|\\!|\\!|q|\\!|\\!|_{n,e}\big{)}.\end{split}$
(43)

From (33) and (43) it follows that the functionals $\mathfrak{C}_{c}^{n}$ and
$\mathfrak{C}_{I}^{n}$ are bounded as

$\displaystyle|\mathfrak{C}_{c}^{n}(\mathbf{v},q)|+|\mathfrak{C}_{I}^{n}(\mathbf{v},q)|\leq
C(\Delta
t+h^{q}+h^{\underline{m}_{v}}+h^{{m}_{q}})(|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}+|\\!|\\!|q|\\!|\\!|_{n,e}),$

where $C>0$ depends on Sobolev norms of the exact solution and the source
function.

Note that $\mathbf{e}^{n}_{h}$ and $d_{h}^{n}$ satisfy the same FE formulation
(18) as $\mathbf{u}^{n}_{h}$ and $p_{h}^{n}$ but with the zero initial
condition and the right-hand side functional given by
$F^{n}(\mathbf{v},q)=\mathfrak{C}_{c}^{n}(\mathbf{v},q)+\mathfrak{C}_{I}^{n}(\mathbf{v},q)$.
Therefore, we can apply the stability results from Lemmas 5 and 6 to estimate
$\mathbf{e}^{n}_{h}$ and $d_{h}^{n}$:

$\|\mathbf{e}_{h}^{k}\|^{2}_{\Omega^{k}_{h}}+\Delta
t\sum_{n=1}^{k}\left(|\\!|\\!|\mathbf{e}_{h}^{n}|\\!|\\!|_{n}^{2}+|\\!|\\!|d_{h}^{n}|\\!|\\!|_{n,e}^{2}\right)\leq
C(\Delta t+h^{q}+h^{\underline{m}_{v}}+h^{{m}_{q}+1})^{2}.$

Applying the triangle inequality and the estimates (38) one more time leads to
our final result.

######  Theorem 8.

Assume the solution $(\mathbf{u},p)$ to (2) is sufficiently smooth and let
$\mathbf{u}^{n}_{h},p^{n}_{h}$ be the solution to (18). Assume that the
discretization parameters satisfy $h^{2}\lesssim\Delta t\lesssim h$. The
following error estimate holds

$\max_{1\leq n\leq
N}\|\mathbf{u}(t^{n})-\mathbf{u}_{h}^{n}\|^{2}_{\Omega^{n}_{h}}+\Delta
t\sum_{n=1}^{N}\left(|\\!|\\!|\mathbf{u}(t^{n})-\mathbf{u}_{h}^{n}|\\!|\\!|_{n}^{2}+|\\!|\\!|p(t^{n})-p_{h}^{n}|\\!|\\!|_{n,e}^{2}\right)\\\
\leq C(\Delta t+h^{q}+h^{\underline{m}_{v}}+h^{{m}_{q}+1})^{2},$ (44)

with a constant $C$ independent of discretization parameters but dependent on
the solution $(\mathbf{u},p)$ and final time $T$.

##  5 Examples of finite element pairs satisfying Assumption 3.1

In this section, we show that several canonical finite element pairs for the
Stokes problem satisfy the three inequalities (16) in Assumption 3.1.

###  5.1 The Mini element

For a tetrahedron $T\in\mathcal{T}_{h}$, let $b_{T}\in\EuScript{P}_{4}(T)$
denote the standard quartic bubble function, i.e., the product of the
barycentric coordinates of $T$. The lowest-order Mini pair with respect to
$\mathcal{T}_{h}$ is given by [1]

$\displaystyle\mathbf{V}_{h}$
$\displaystyle=\\{\mathbf{v}\in\mathbf{H}^{1}(\widehat{\Omega}):\
\mathbf{v}|_{T}\in\bm{\EuScript{P}}_{1}(T)+b_{T}\bm{\EuScript{P}}_{0}(T)\
\forall T\in\mathcal{T}_{h}\\},$ $\displaystyle Q_{h}$ $\displaystyle=\\{q\in
H^{1}(\widehat{\Omega}):\ q|_{T}\in\EuScript{P}_{1}(T)\ \forall
T\in\mathcal{T}_{h}\\}.$

In this setting we can take $\underline{m}_{v}=1$, $\overline{m}_{v}=4$ and
$m_{q}=1$.

We now verify conditions (16). Given $q\in Q_{h}^{n}$, we set
$\mathbf{v}\in\mathbf{V}_{h}^{n}$ so that
$\mathbf{v}|_{T}=h_{T}^{2}b_{T}\nabla q|_{T}$ for all
$T\in\mathcal{T}_{h}^{n,i}$. The function $\mathbf{v}$ is extended to
$\Omega^{n}_{e}$ by zero. The results in [20, Section 6.5] show that
(16a)–(16b) is satisfied. We also have by a simple scaling argument

$\displaystyle\|\mathbf{v}\|_{\Omega^{n}}^{2}$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}^{n,i}}\|\mathbf{v}\|_{T}^{2}=\sum_{T\in\mathcal{T}_{h}^{n,i}}h_{T}^{4}\|b_{T}\nabla
q\|_{T}^{2}\simeq\sum_{T\in\mathcal{T}_{h}^{n,i}}h_{T}^{4}\|\nabla
q\|_{T}^{2}\lesssim h^{2}|\\!|\\!|q|\\!|\\!|_{n,i}^{2}.$

Thus, (16c) is satisfied as well.

###  5.2 Taylor-Hood

The (generalized) Taylor–Hood finite element pair is given by

$\displaystyle\mathbf{V}_{h}$
$\displaystyle=\\{\mathbf{v}\in\mathbf{H}^{1}(\widehat{\Omega}):\
\mathbf{v}|_{T}\in\bm{\EuScript{P}}_{m}(T)\ \forall T\in\mathcal{T}_{h}\\},$
$\displaystyle Q_{h}$ $\displaystyle=\\{q\in H^{1}(\widehat{\Omega}):\
q|_{T}\in\EuScript{P}_{m-1}(T)\ \forall T\in\mathcal{T}_{h}\\},$

where $m\geq 2$. Thus, in this case $\overline{m}_{v}=\underline{m}_{v}=m$ and
$m_{q}=m-1$ in (14)–(15). Denote by $\mathcal{E}_{h}^{n,i}$ the set of
interior one-dimensional edges of the the interior triangulation
$\mathcal{T}_{h}^{n,i}$. We then denote by $\tilde{\mathcal{T}}_{h}^{n,i}$ the
members in $\mathcal{T}_{h}^{n,i}$ that have at least three edges in
$\mathcal{E}^{n,i}_{h}$ (cf. Remark 3.2). We assume that the domain of
pressure ghost-stabilization is chosen such that (17) is satisfied. This is
the case provided $c_{\delta_{h}}$ is sufficiently large (but still $O(1)$).

We denote the set of interior edges of $\tilde{\mathcal{T}}_{h}^{n,i}$ by
$\tilde{\mathcal{E}}_{h}^{n,i}$. Then for $e\in\tilde{\mathcal{E}}_{h}^{n,i}$,
we let $\phi_{e}$ denote the quadratic bubble function associated with $e$,
and let ${\bf t}_{e}$ be a unit tangent vector of $e$. Note that $\phi_{e}$
has support on the tetrahedra that have $e$ as an edge, and the number of such
tetrahedra is uniformly bounded due to the shape-regularity of
$\tilde{\mathcal{T}}_{h,i}^{n}$.

For a given $q\in Q_{h}^{n}$, we define

$\mathbf{v}=\sum_{e\in\tilde{\mathcal{E}}_{h}^{n,i}}h_{e}^{2}\phi_{e}(\nabla
q\cdot{\bf t}_{e}){\bf t}_{e}.$

Because $q$ is continuous, we see that $\nabla q\cdot{\bf t}_{e}$ is single-
valued on $e$, and thus $\mathbf{v}$ is continuous and a piecewise polynomial
of degree $m$; hence, $\mathbf{v}\in\mathbf{V}_{h}^{n}$.

It is shown in [20, Section 6.1] that (16a)–(16b) is satisfied, thus is it
remains to show (16c). This follows from the identity
$\|\phi_{e}\|_{\infty}=1$ and the shape-regularity of the triangulation:

$\displaystyle\|\mathbf{v}\|_{\tilde{\Omega}_{h}^{n}}\lesssim\sum_{T\in\tilde{\mathcal{T}}_{h,i}^{n}}h_{T}^{2}\|\nabla
q\|_{T}^{2}.$

###  5.3 $\EuScript{P}_{3}-\EuScript{P}_{0}$

As our final example, we consider the $\bm{\EuScript{P}}_{3}-\EuScript{P}_{0}$
pair. In particular, the discrete velocity space is the cubic Lagrange space,
and the discrete pressure space consists of piecewise constants:

$\displaystyle\mathbf{V}_{h}$
$\displaystyle=\bm{\EuScript{P}}_{3}^{c}(\mathcal{T}_{h})=\\{\mathbf{v}\in\mathbf{H}^{1}(\widehat{\Omega}):\
\mathbf{v}|_{T}\in\bm{\EuScript{P}}_{3}(T)\ \forall T\in\mathcal{T}_{h}\\},$
$\displaystyle Q_{h}$ $\displaystyle=\EuScript{P}_{0}(\mathcal{T}_{h})=\\{q\in
L^{2}(\widehat{\Omega}):\ q|_{T}\in\EuScript{P}_{0}(T)\ \forall
T\in\mathcal{T}_{h}\\}.$

For each interior face $F\in\mathcal{F}_{h,i}^{n}$ with $F=\partial
T_{1}\cap\partial T_{2}$, we denote by $\mathbf{n}_{j}$ the outward unit
normal of $\partial T_{j}$ restricted to $F$. Then for given $q\in Q^{n}_{h}$,
we define $\mathbf{v}\in\mathbf{V}^{n}_{h}$ such that for all
$F\in\mathcal{F}_{h,i}^{n}$,

$\int_{F}\mathbf{v}\cdot\mathbf{n}_{1}\,ds=-h_{F}\int_{F}(q_{1}\mathbf{n}_{1}+q_{2}\mathbf{n}_{2})\cdot\mathbf{n}_{1}\,ds=h_{F}\int_{F}\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\cdot\mathbf{n}_{1}\,ds,$

where $q_{j}=q|_{T_{j}}$. Note that this condition implies
$\int_{F}\mathbf{v}\cdot\mathbf{n}_{F}\,ds=-h_{F}\int_{F}\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\cdot\mathbf{n}_{F}\,ds$
for any unit normal of $F\in\mathcal{F}_{h,i}^{n}$. We further specify that
$\mathbf{v}=0$ on all vertices and edges in $\mathcal{T}_{h,i}^{n}$,
$\mathbf{v}\times\mathbf{n}_{F}=0$ on all faces $F\in\mathcal{F}_{h,i}^{n}$,
and $\mathbf{v}=0$ on the boundary of $\Omega_{h,i}^{n}$. We extend
$\mathbf{v}$ to $\Omega_{h,e}^{n}$ by zero.

By the divergence theorem, and using that $q$ is piecewise constant, we have

$\displaystyle b_{h}^{n}(\mathbf{v},q)=\int_{\Omega_{h,i}^{n}}({\rm
div}\,\mathbf{v})q\,dx$
$\displaystyle=-\sum_{T\in\mathcal{T}_{h,i}^{n}}\int_{\partial
T}q(\mathbf{v}\cdot\mathbf{n}_{\partial
T})\,ds=\sum_{F\in\mathcal{F}_{h,i}^{n}}h_{F}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}=|\\!|\\!|q|\\!|\\!|_{n,i}^{2}.$

Thus, (16b) is satisfied. A scaling argument also yields on each
$T\in\mathcal{T}_{h,i}^{n}$,

$\displaystyle|\mathbf{v}|_{H^{m}(T)}$ $\displaystyle\lesssim
h_{T}^{2-2m}\sum_{\mathcal{F}_{h,i}^{n}\ni F\subset\partial
T}h_{F}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}.$

Consequently, by another scaling argument,

$\displaystyle|\\!|\\!|\mathbf{v}|\\!|\\!|_{n,e}^{2}$
$\displaystyle\lesssim\|\nabla\mathbf{v}\|_{\Omega_{h,i}^{n}}^{2}+h^{-2}\|\mathbf{v}\|_{\Omega_{h,i}^{n}}^{2}\lesssim\sum_{F\in\mathcal{F}_{h,i}^{n}}h_{F}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}=|\\!|\\!|q|\\!|\\!|_{n,i}^{2},$
$\displaystyle\|\mathbf{v}\|_{\Omega_{h}^{n}}^{2}$
$\displaystyle=\|\mathbf{v}\|_{\Omega_{h,i}^{n}}^{2}\lesssim
h^{2}\sum_{F\in\mathcal{F}_{h,i}^{n,i}}h_{F}\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\|_{F}^{2}=h^{2}|\\!|\\!|q|\\!|\\!|_{n,i}^{2},$

and therefore (16a) and (16c) are satisfied as well.

##  Appendix A Proof of Lemma 4

We first note that if $Q_{h}^{n}\subset H^{1}(\Omega_{h,e}^{n})$, then a
standard trace inequality and the definition of
$|\\!|\\!|\cdot|\\!|\\!|_{n,e}$ yields

$\|q\|_{\Gamma_{h}^{n}}\lesssim\|q\|_{H^{1}(\Omega_{h}^{n})}\lesssim
h^{-1}|\\!|\\!|q|\\!|\\!|_{n,e}+\|q\|_{\Omega_{h}^{n}}.$ (45)

To establish (29) in this case, we first apply a standard Poincare-Friedrich
inequality

$\|q\|_{\Omega_{h,i}^{n}}\lesssim\|\nabla q\|_{\Omega_{h,i}^{n}}\qquad\forall
q\in\mathring{L}^{2}(\Omega_{h,i}^{n})\cap H^{1}(\Omega_{h,i}^{n}),$

and (25) to conclude

$\displaystyle\|q\|_{\Omega_{h}^{n}}\lesssim\|q\|_{\Omega_{h,i}^{n}}+|q|_{J_{h}^{n}}\lesssim\|\nabla
q\|_{\Omega_{h,i}^{n}}+|q|_{J_{h}^{n}}\lesssim
h^{-1}\big{(}|\\!|\\!|q|\\!|\\!|_{n,i}+|q|_{J_{h}^{n}}\big{)}\lesssim
h^{-1}|\\!|\\!|q|\\!|\\!|_{n,e}\quad\forall q\in Q_{h}^{n}.$

The estimate (29) then follows from this inequality and (45).

Thus, it suffices to prove (29) in the case $Q_{h}^{n}$ consists of
discontinuous polynomials. To this end, we introduce an enriching operator
$E_{h}:Q_{h}^{n}\to Q_{h}^{n}\cap H^{1}(\Omega_{h,e}^{n})$ constructed by
averaging [5]. Let

$\mathcal{T}^{n}_{T}=\\{T^{\prime}\in\mathcal{T}_{h,e}^{n}:\
\bar{T}\cap\bar{T}^{\prime}\neq\emptyset\\},$

and let $\mathcal{F}_{T}^{n,I}$ denote the set of interior faces of
$\mathcal{T}^{n}_{T}$. Then there holds

$\displaystyle|q-E_{h}q|_{H^{\ell}(T)}^{2}$ $\displaystyle\lesssim
h_{T}^{2-2\ell}\sum_{F\in\mathcal{F}^{n,I}_{T}}h_{F}^{-1}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{L^{2}(F)}^{2}\qquad\ell=0,1.$
(46)

It then follows from (46) and the trace inequality

$\displaystyle\|q\|_{T\cap\Gamma_{h}^{n}}\lesssim
h_{T}^{-1/2}\|q\|_{T}+h_{T}^{1/2}\|\nabla q\|_{T}\qquad\forall q\in H^{1}(T)$

that

$\begin{split}\|q-E_{h}q\|_{\Gamma_{h}}^{2}&=\sum_{T\in\mathcal{T}_{h,e}^{n}}\|q-E_{h}q\|_{T\cap\Gamma_{h}^{n}}^{2}\\\
&\lesssim\sum_{T\in\mathcal{T}_{h,e}^{n}}\left(h_{T}^{-1}\|q-E_{h}q\|_{T}^{2}+h_{T}\|\nabla(q-E_{h}q)\|_{T}^{2}\right)\\\
&\lesssim
h^{-1}\sum_{F\in\mathcal{F}_{h,e}^{n}}h_{F}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}\lesssim
h^{-1}|\\!|\\!|q|\\!|\\!|_{n,e}^{2}.\end{split}$ (47)

Furthermore by a standard trace inequality and (46), we have

$\begin{split}\|E_{h}q\|_{\Gamma_{h}^{n}}^{2}&\lesssim\|E_{h}q\|_{H^{1}(\Omega_{h}^{n})}^{2}\leq\|E_{h}q\|_{H^{1}(\Omega_{h,e}^{n})}^{2}\\\
&\lesssim\sum_{T\in\mathcal{T}_{h,e}^{n}}\|q\|_{H^{1}(T)}^{2}+\sum_{F\in\mathcal{F}^{n}_{h,e}}h_{F}^{-1}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}\\\
&\lesssim
h^{-2}|\\!|\\!|q|\\!|\\!|_{n,e}^{2}+\|q\|_{\Omega_{h,e}^{n}}^{2}.\end{split}$
(48)

Combining (47)–(48) yields

$\displaystyle\|q\|_{\Gamma_{h}^{n}}\lesssim
h^{-1}|\\!|\\!|q|\\!|\\!|_{n,e}+\|q\|_{\Omega_{h,e}^{n}}.$ (49)

Finally, since $q|_{\Omega_{h,i}^{n}}\in\mathring{L}^{2}(\Omega_{h,i}^{n})$,
we apply the discrete Poincare-Friedrich inequality [5, Theorem 10.6.12]

$\displaystyle\|q\|_{\Omega_{h,i}^{n}}^{2}$
$\displaystyle\lesssim\sum_{T\in\mathcal{T}_{h,i}^{n}}\|\nabla
q\|_{T}^{2}+\sum_{F\in\mathcal{F}_{h,i}^{n}}h_{F}^{-1}\left\|\bigl{[}\bigl{[}q\bigr{]}\bigr{]}\right\|_{F}^{2}\lesssim
h^{-2}|\\!|\\!|q|\\!|\\!|_{n,i}^{2},$

and (25) to conclude

$\displaystyle\|q\|_{\Omega_{h,e}^{n}}\lesssim\|q\|_{\Omega_{h,i}^{n}}+|q|_{J_{h}^{n}}\lesssim
h^{-1}\big{(}|\\!|\\!|q|\\!|\\!|_{n,i}+|q|_{J_{h}^{n}}\big{)}\lesssim
h^{-1}|\\!|\\!|q|\\!|\\!|_{n,e}.$

Combined with (49), we obtain (29).

## References

  * [1] D. N. Arnold, F. Brezzi, and M. Fortin, A stable finite element for the Stokes equations, Calcolo, 21 (1984), pp. 337–344 (1985).
  * [2] Y. Bazilevs, K. Takizawa, and T. E. Tezduyar, Computational fluid-structure interaction: methods and applications, John Wiley & Sons, 2013\.
  * [3] M. Besier and W. Wollner, On the pressure approximation in nonstationary incompressible flow simulations on dynamically varying spatial meshes, Internat. J. Numer. Methods Fluids, 69 (2012), pp. 1045–1064.
  * [4] A. Brenner, E. Bänsch, and M. Bause, A priori error analysis for finite element approximations of the Stokes problem on dynamic meshes, IMA J. Numer. Anal., 34 (2014), pp. 123–146.
  * [5] S. C. Brenner and L. R. Scott, The mathematical theory of finite element methods, vol. 15 of Texts in Applied Mathematics, Springer, New York, third ed., 2008.
  * [6] A. Buffa and C. Ortner, Compact embeddings of broken Sobolev spaces and applications, IMA J. Numer. Anal., 29 (2009), pp. 827–855.
  * [7] E. Burman, Ghost penalty, C. R. Math. Acad. Sci. Paris, 348 (2010), pp. 1217–1220.
  * [8] E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing, CutFEM: discretizing geometry and partial differential equations, Internat. J. Numer. Methods Engrg., 104 (2015), pp. 472–501.
  * [9] E. Burman, S. Frei, and A. Massing, Eulerian time-stepping schemes for the non-stationary Stokes equations on time-dependent domains, Numer. Math., 150 (2022), pp. 423–478.
  * [10] E. Burman and P. Hansbo, Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method, Appl. Numer. Math., 62 (2012), pp. 328–341.
  * [11] S. Claus and P. Kerfriden, A CutFEM method for two-phase flow problems, Comput. Methods Appl. Mech. Engrg., 348 (2019), pp. 185–206.
  * [12] M. Costabel and A. McIntosh, On Bogovskiĭ and regularized Poincaré integral operators for de Rham complexes on Lipschitz domains, Mathematische Zeitschrift, 265 (2010), pp. 297–320.
  * [13] F. Duarte, R. Gormaz, and S. Natesan, Arbitrary Lagrangian-Eulerian method for Navier-Stokes equations with moving boundaries, Comput. Methods Appl. Mech. Engrg., 193 (2004), pp. 4819–4836.
  * [14] L. Formaggia and F. Nobile, A stability analysis for the arbitrary Lagrangian Eulerian formulation with finite elements, East-West J. Numer. Math., 7 (1999), pp. 105–131.
  * [15] L. Formaggia, A. Quarteroni, and A. Veneziani, Cardiovascular Mathematics: Modeling and simulation of the circulatory system, vol. 1, Springer Science & Business Media, 2010.
  * [16] V. Girault and P. A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer, Berlin, 1986.
  * [17] R. Glowinski, T.-W. Pan, T. I. Hesla, and D. D. Joseph, A distributed Lagrange multiplier/fictitious domain method for particulate flows, Int. J. Multiph. Flow, 25 (1999), pp. 755–794.
  * [18] S. Gross, M. A. Olshanskii, and A. Reusken, A trace finite element method for a class of coupled bulk-interface transport problems, ESAIM Math. Model. Numer. Anal., 49 (2015), pp. 1303–1330.
  * [19] S. Gross and A. Reusken, Numerical methods for two-phase incompressible flows, vol. 40, Springer Science & Business Media, 2011.
  * [20] J. Guzmán and M. Olshanskii, Inf-sup stability of geometrically unfitted Stokes finite elements, Math. Comp., 87 (2018), pp. 2091–2112.
  * [21] C. W. Hirt, A. A. Amsden, and J. L. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds [J. Comput. Phys. 14 (1974), no. 3, 227–253], J. Comput. Phys., 135 (1997), pp. 198–216.  With an introduction by L. G. Margolin, Commemoration of the 30th anniversary {of J. Comput. Phys.}.
  * [22] I. Kesler, R. Lan, and P. Sun, The arbitrary Lagrangian-Eulerian finite element method for a transient stokes/parabolic interface problem, Int. J. Numer. Anal. Model., 18 (2021), pp. 339–361.
  * [23] C. Lehrenfeld and M. Olshanskii, An Eulerian finite element method for PDEs in time-dependent domains, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 585–614.
  * [24] A. Lozovskiy, M. A. Olshanskii, and Y. V. Vassilevski, A quasi-Lagrangian finite element method for the Navier-Stokes equations in a time-dependent domain, Comput. Methods Appl. Mech. Engrg., 333 (2018), pp. 55–73.
  * [25] C. Ma, T. Tian, and W. Zheng, High-order unfitted characteristic finite element methods for moving interface problem of Oseen equations, J. Comput. Appl. Math., 425 (2023), pp. Paper No. 115028, 29.
  * [26] A. Massing, M. G. Larson, A. Logg, and M. E. Rognes, A stabilized Nitsche fictitious domain method for the Stokes problem, J. Sci. Comput., 61 (2014), pp. 604–628.
  * [27] A. Masud and T. J. R. Hughes, A space-time Galerkin/least-squares finite element formulation of the Navier-Stokes equations for moving domain problems, Comput. Methods Appl. Mech. Engrg., 146 (1997), pp. 91–126.
  * [28] T. Miyakawa and Y. Teramoto, Existence and periodicity of weak solutions of the Navier-Stokes equations in a time dependent domain, Hiroshima Math. J., 12 (1982), pp. 513–528.
  * [29] J. Neustupa, Existence of a weak solution to the Navier-Stokes equation in a general time-varying domain by the Rothe method, Math. Methods Appl. Sci., 32 (2009), pp. 653–683.
  * [30] M. A. Olshanskii, A. Reusken, and P. Schwering, An Eulerian finite element method for tangential Navier–Stokes equations on evolving surfaces, arXiv preprint arXiv:2302.00779, (2023).
  * [31] C. S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys., 25 (1977), pp. 220–252.
  * [32] J. Preuß, Higher order unfitted isoparametric space-time FEM on moving domains, master’s thesis, NAM, University of Göttingen, 2018.
  * [33] J. San Martín, L. Smaranda, and T. Takahashi, Convergence of a finite element/ALE method for the Stokes equations in a domain depending on time, J. Comput. Appl. Math., 230 (2009), pp. 521–545.
  * [34] V. A. Solonnikov, Solvability of a problem on the motion of a viscous incompressible fluid bounded by a free surface, Mathematics of the USSR-Izvestiya, 11 (1977), p. 1323.
  * [35] , On the transient motion of an isolated volume of viscous incompressible fluid, Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, 51 (1987), pp. 1065–1087.
  * [36] E. M. Stein, Singular integrals and differentiability properties of functions (PMS-30), vol. 30, Princeton university press, 2016.
  * [37] T. E. Tezduyar, M. Behr, S. Mittal, and J. Liou, A new strategy for finite element computations involving moving boundaries and interfaces—the deforming-spatial-domain/space-time procedure. II. Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders, Comput. Methods Appl. Mech. Engrg., 94 (1992), pp. 353–371.
  * [38] H. von Wahl, T. Richter, and C. Lehrenfeld, An unfitted Eulerian finite element method for the time-dependent Stokes problem on moving domains, IMA J. Numer. Anal., 42 (2022), pp. 2505–2544.