# POD-ROMs for incompressible flows including snapshots of the temporal
derivative of the full order solution

Bosco García-Archilla Departamento de Matemática Aplicada II, Universidad de
Sevilla, Sevilla, Spain. Research is supported by Spanish MCINYU under grants
PGC2018-096265-B-I00 and PID2019-104141GB-I00<EMAIL_ADDRESS>Volker John
Weierstrass Institute for Applied Analysis and Stochastics, Leibniz Institute
in Forschungsverbund Berlin e. V. (WIAS), Mohrenstr. 39, 10117 Berlin,
Germany. Freie Universität of Berlin, Department of Mathematics and Computer
Science, Arnimallee 6, 14195 Berlin, Germany. Julia Novo Departamento de
Matemáticas, Universidad Autónoma de Madrid, Spain. Research is supported by
Spanish MINECO under grants PID2019-104141GB-I00 and VA169P20
<EMAIL_ADDRESS>
###### Abstract

In this paper we study the influence of including snapshots that approach the
velocity time derivative in the numerical approximation of the incompressible
Navier-Stokes equations by means of proper orthogonal decomposition (POD)
methods. Our set of snapshots includes the velocity approximation at the
initial time from a full order mixed finite element method (FOM) together with
approximations to the time derivative at different times. The approximation at
the initial velocity can be replaced by the mean value of the velocities at
the different times so that implementing the method to the fluctuations, as
done mostly in practice, only approximations to the time derivatives are
included in the set of snapshots. For the POD method we study the differences
between projecting onto $L^{2}$ and $H^{1}$. In both cases pointwise in time
error bounds can be proved. Including grad-div stabilization both in the FOM
and POD methods error bounds with constants independent on inverse powers of
the viscosity can be obtained.

AMS subject classifications. 65M12, 65M15, 65M60.
Keywords. incompressible Navier–Stokes equations, proper orthogonal
decomposition (POD), reduced order models (ROMs), snapshots of the temporal
derivative of the full order solution, grad-div stabilization, robust
pointwise in time estimates

##  1 Introduction

It is well known that the computational cost of direct numerical simulations,
also called full order methods (FOMs), can be reduced by using reduced order
models (ROMs). In this paper, we study ROMs based on proper orthogonal
decomposition (POD) methods, so-called POD-ROMs. The computation of the
reduced basis uses solutions of a FOM, so-called snapshots.

We study incompressible flow problems that are modeled by means of the
incompressible Navier–Stokes equations

$\begin{array}[]{rcll}\boldsymbol{u}_{t}-\nu\Delta\boldsymbol{u}+(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}+\nabla
p&=&\boldsymbol{f}&\text{in }\ (0,T]\times\Omega,\\\
\nabla\cdot\boldsymbol{u}&=&0&\text{in }\ (0,T]\times\Omega,\end{array}$ (1)

in a bounded domain $\Omega\subset{\mathbb{R}}^{d}$, $d\in\\{2,3\\}$, with
initial condition $\boldsymbol{u}(0)=\boldsymbol{u}^{0}$. In (1),
$\boldsymbol{u}$ is the velocity field, $p$ the kinematic pressure, $\nu>0$
the kinematic viscosity coefficient, and $\boldsymbol{f}$ represents the
accelerations due to external body forces acting on the fluid. The
Navier–Stokes equations (1) have to be complemented with boundary conditions.
For simplicity, we only consider homogeneous Dirichlet boundary conditions
$\boldsymbol{u}=\boldsymbol{0}$ on $[0,T]\times\partial\Omega$.

This paper studies the impact of including approximations of the temporal
derivative of the velocity in the set of snapshots. The idea consists in
taking, in addition to the mixed finite element approximation to the velocity
at the initial time, $\left\\{\boldsymbol{u}_{h}(t_{0})\right\\}$, the time
derivatives of the mixed finite element approximations,
$\left\\{\boldsymbol{u}_{h,t}(t_{j})\right\\}$. These temporal derivatives can
be easily computed using the right-hand side of the mixed finite element
Galerkin equation. In the present paper we follow an idea presented in the
recent paper [22] in which it is shown that there is no need to include in the
set of snapshots other than one approximation to the velocity at a fixed time,
instead of the full set $\left\\{\boldsymbol{u}_{h}(t_{j})\right\\}_{j=0}^{M}$
as it is usually done in the literature. The numerical analysis in [22] is
carried out for the heat equation and with difference quotients,
$\left\\{(\boldsymbol{u}_{h}(t_{j})-\boldsymbol{u}_{h}(t_{j-1})/(t_{j}-t_{j-1})\right\\}$,
approaching the time derivative. In the present paper we consider instead the
Galerkin time derivatives although the analysis for the difference quotients
case is essentially the same (or even simpler). Actually, in practice, any
approximation to the time derivative can work equally well. We also prove that
the snapshot at the initial value can be replaced by the mean value
$\overline{\boldsymbol{u}}_{h}=(M+1)^{-1}\sum_{j=0}^{M}\boldsymbol{u}_{h}(t_{j})$,
which can be more efficient in the numerical simulations. It is standard to
apply the POD method to the fluctuations
$\left\\{\boldsymbol{y}_{h}(t_{j})\right\\}_{j=0}^{M}=\left\\{\boldsymbol{u}_{h}(t_{j})-\overline{\boldsymbol{u}}_{h}\right\\}_{j=0}^{M}$
that have no mean by definition. Then, with the method we propose only
approximations to the derivatives are needed in this case.

Several works in the literature have already studied the subject of increasing
the set of snapshots with approximations to the time derivatives. However,
apart from [22], all of them include as snapshots the set of the
approximations at different times, instead of only one snapshot. Also,
starting with the pioneering paper [21], most of these papers include a
different type of approximation than considered in this paper, namely
difference quotients. In particular, to the best of our knowledge, this is the
first paper that studies the inclusion of the temporal derivatives of the
mixed finite element velocity approximations in order to generate the reduced
order basis for the incompressible Navier–Stokes equations. With this more
general setting we can deduce that in practice any approximation to the time
derivative produces essentially the same results.

The initial motivation for investigating a different approximation than
difference quotients in the set of snapshots is that the results for
difference quotients are ambivalent. On the one hand, from the theoretical
point of view, the inclusion of the difference quotients possesses some
advantages. First of all, it allows to prove optimal error bounds for the POD-
ROM when the POD basis functions are based on the projection onto the Hilbert
space $X=H_{0}^{1}(\Omega)^{d}$, see [21, 15, 27]. In this way, the standard
finite element error analysis is mimicked, in which the Ritz or Stokes
projection is used to split the error in a projection error and a discrete
remainder. It was observed that if the POD basis functions are based on the
projection onto the Hilbert space $X=L^{2}(\Omega)^{d}$, the difference
quotients are not needed to prove optimal error bounds in certain norms, see
[4, 15, 27, 23]. However, as pointed out in [19], even in this case, the
inclusion of the difference quotients allows to get pointwise estimates in
time that generally cannot be proved if there are no difference quotients in
the set of snapshots. On the other hand, from the numerical point of view, it
is not clear that the difference quotients should be included in the actual
simulations with the POD-ROM. In fact, it is reported in [18, 17] that the
POD-ROM without the difference quotients performs considerably better than
with the difference quotients.

Trying to keep the theoretical advantages of including approximations of the
temporal derivative of the velocity in the set of snapshots but relaxing their
drawbacks in practical simulations, we study in this paper, both theoretically
and numerically, the inclusion of time derivatives of the discrete velocity.
As in [22], our approach has half the number of snapshots as in the standard
POD finite quotient approach. Also, as in the difference quotient case, we are
able to get pointwise in time estimates using as projecting space
$X=H_{0}^{1}(\Omega)^{d}$ and $L^{2}(\Omega)^{d}$. To this end, we follow [22,
Lemma 3.3] (see also [19, Lemma 3.6]) to prove that the $X$-norm of a function
at any point in time (let’s say $t_{j}$) is bounded in terms of the $X$-norm
of the value at $t_{0}$ plus the mean values of its time derivative taken in a
time interval (let’s say $\left\\{t_{0},t_{1},\ldots,t_{M}\right\\}$, up to an
error that tends to zero with the length of the time step. From the numerical
point of view, including the snapshots of the time derivatives avoids the
potential problem of performing badly conditioned operations in the
computation of the snapshots like
$\left\\{(\boldsymbol{u}_{h}(t_{j})-\boldsymbol{u}_{h}(t_{j-1}))/(t_{j}-t_{j-1})\right\\}$
because both numerator and denominator may suffer from numerical cancellation.
Finally, we mimic model reduction ideas coming from dynamical systems, where
using snapshots from the time derivative is a more common approach, see for
example [20].

For the error analysis in the present paper, instead of considering a concrete
fully discrete scheme from which the values
$\left\\{\boldsymbol{u}_{h}(t_{j})\right\\}$ are taken, we consider a
continuous-in-time method, which has some advantages. In practice, one always
computes the snapshots with a fully discrete method but the error analysis
based on the continuous-in-time method holds for any time integrator used in
the FOM. With this approach one can use different time steps for the FOM,
which produces the snapshots, and for the POD-ROM method. Our error analysis
takes into account the temporal error coming from the POD-ROM. Finally,
following [23], we analyze the case in which stabilized approximations are
computed both for the FOM and the POD-ROM. More precisely, the considered
finite element method is based on a Galerkin discretization plus grad-div
stabilization with pairs of inf-sup stable elements. For the POD-ROM we also
use grad-div stabilization. In this way, the constants in the error bounds for
the snapshots do not depend explicitly on inverse powers of the viscosity,
i.e., they do not blow up for small viscosity coefficients, see [8]. Adapting
the results from [23], the same holds for the error bounds of the POD-ROM. The
importance of such so-called robust methods is stated in the survey [9]: In
the case of small viscosity coefficients and coarse grids, only robust
estimates provide useful information about the behavior of a numerical method
on coarse grids if the analytic solution is smooth.

In the numerical studies, we compare the different approaches obtained by
taking $X=H_{0}^{1}(\Omega)^{d}$ and $X=L^{2}(\Omega)^{d}$ in combination with
one of the following sets: the set of snapshots at different times, the set of
difference quotients, and the set of Galerkin time derivatives. We cannot
deduce from these studies that any of the approaches is much better than the
other ones and the necessary comprehensive numerical studies are outside the
scope of this paper, which is a rigorous numerical analysis from which
interesting properties and sharp bounds for the different methods can be
deduced.

The paper is organized as follows. In Section 2 we state some preliminaries
and notations. The POD method and some a priori bounds for the projection of
the FOM approximation onto the POD space are shown in Section 3. The analysis
of the POD method is included in Section 4 with a first subsection for the
case $X=H_{0}^{1}(\Omega)^{d}$ and a second one for the case
$X=L^{2}(\Omega)^{d}$. As stated above, Section 5 is devoted to study the
performance of the methods with some numerical experiments. Finally, we have
included an appendix in which we get robust bounds for the second time
derivative of the FOM approximation in some norms, since this time derivative
appears in our bounds in the a priori error analysis.

##  2 Preliminaries and notations

Standard symbols will be used for Lebesgue and Sobolev spaces, with the usual
convention that $W^{s,2}(\Omega)=H^{s}(\Omega)$, $s\geq 1$. The inner product
in $L^{2}(\Omega)^{d}$, $d\geq 1$, is denoted by $(\cdot,\cdot)$.

The following Sobolev imbeddings [1] will be used in the analysis: For
$q\in[1,\infty)$, there exists a constant $C=C(\Omega,q)$ such that

$\|v\|_{L^{q^{\prime}}}\leq
C\|v\|_{W^{s,q}},\,\,\quad\frac{1}{q^{\prime}}\geq\frac{1}{q}-\frac{s}{d}>0,\quad
q<\infty,\quad v\in W^{s,q}(\Omega)^{d}.$ (2)

We will denote by $C_{p}$ the constant in the Poincaré inequality

$\|\boldsymbol{v}\|_{0}\leq
C_{p}\|\nabla\boldsymbol{v}\|_{0},\quad\boldsymbol{v}\in
H_{0}^{1}(\Omega)^{d}.$ (3)

The following inequality can be found in [16, Remark 3.35]

$\|\nabla\cdot\boldsymbol{v}\|_{0}\leq\|\nabla\boldsymbol{v}\|_{0},\quad\boldsymbol{v}\in
H_{0}^{1}(\Omega)^{d}.$ (4)

Let us denote by ${\boldsymbol{V}}=H_{0}^{1}(\Omega)^{d}$ and
$Q=L_{0}^{2}(\Omega)=\\{q\in L^{2}(\Omega)\ \mid\ (q,1)=0\\}$.

Let $\mathcal{T}_{h}=(\sigma_{j}^{h},\phi_{j}^{h})_{j\in J_{h}}$, $h>0$, be a
family of partitions of $\overline{\Omega}$, where $h$ denotes the maximum
diameter of the mesh cells $\sigma_{j}^{h}\in\mathcal{T}_{h}$, and
$\phi_{j}^{h}$ are the mappings from the reference simplex $\sigma_{0}$ onto
$\sigma_{j}^{h}$. We shall assume that the family of partitions is shape-
regular and quasi-uniform. On these partitions, we define the following finite
element spaces

$\displaystyle Y_{h}^{l}$ $\displaystyle=$ $\displaystyle\left\\{v_{h}\in
C^{0}(\overline{\Omega})\ \mid\ {v_{h}}_{\mid_{K}}\in{\mathbb{P}}_{l}(K),\
\forall\ K\in\mathcal{T}_{h}\right\\},\ l\geq
1,\quad{\boldsymbol{Y}}_{h}^{l}=\left(Y_{h}^{l}\right)^{d},$
$\displaystyle{\boldsymbol{X}}_{h}^{l}$ $\displaystyle=$
$\displaystyle{\boldsymbol{Y}}_{h}^{l}\cap H_{0}^{1}(\Omega)^{d},\quad
Q_{h}^{l}=Y_{h}^{l}\cap L_{0}^{2}(\Omega),$
$\displaystyle{\boldsymbol{V}}_{h,l}$ $\displaystyle=$
$\displaystyle{\boldsymbol{X}}_{h}^{l}\cap\left\\{{\boldsymbol{v}}_{h}\in
H_{0}^{1}(\Omega)^{d}\ \mid\ (q_{h},\nabla\cdot{\boldsymbol{v}}_{h})=0\
\forall\ q_{h}\in Q_{h}^{l-1}\right\\},\quad l\geq 2.$ (5)

The space ${\boldsymbol{V}}_{h,l}$ is the space of discretely divergence-free
functions.

Since the family of partitions is quasi-uniform, the following inverse
inequality holds for each $\boldsymbol{v}_{h}\in Y_{h}^{l}$, e.g., see [6,
Theorem 3.2.6],

$|\boldsymbol{v}_{h}|_{W^{m,p}(K)}\leq
c_{\mathrm{inv}}h_{K}^{n-m-d\left(\frac{1}{q}-\frac{1}{p}\right)}|\boldsymbol{v}_{h}|_{W^{n,q}(K)},$
(6)

where $0\leq n\leq m\leq 1$, $1\leq q\leq p\leq\infty$, and $h_{K}$ is the
diameter of $K\in\mathcal{T}_{h}$.

The analysis uses a modified Stokes projection $\boldsymbol{s}_{h}^{m}\ :\
{\boldsymbol{V}}\rightarrow{\boldsymbol{V}}_{h,l}$ that was introduced in [7]
and that is defined by

$(\nabla\boldsymbol{s}_{h}^{m},\nabla\boldsymbol{\varphi}_{h})=(\nabla\boldsymbol{u},\nabla\boldsymbol{\varphi}_{h}),\quad\forall\
\boldsymbol{\varphi}_{h}\in\boldsymbol{V}_{h,l}.$ (7)

This projection satisfies the following error bound, see [7],

$\|\boldsymbol{u}-\boldsymbol{s}_{h}^{m}\|_{0}+h\|\boldsymbol{u}-\boldsymbol{s}_{h}^{m}\|_{1}\leq
C\|\boldsymbol{u}\|_{j}h^{j},\quad 1\leq j\leq l+1.$ (8)

From [5], we also have

$\|\nabla\boldsymbol{s}_{h}^{m}\|_{\infty}\leq
C\|\nabla\boldsymbol{u}\|_{\infty},$ (9)

and from [12, Lemma 3.8]

$\displaystyle\|\boldsymbol{s}_{h}^{m}\|_{\infty}$ $\displaystyle\leq$
$\displaystyle C(\|\boldsymbol{u}\|_{d-2}\|\boldsymbol{u}\|_{2})^{1/2},$ (10)
$\displaystyle\|\nabla\boldsymbol{s}_{h}^{m}\|_{L^{2d/(d-1)}}$
$\displaystyle\leq$ $\displaystyle
C\bigl{(}\|\boldsymbol{u}\|_{1}\|\boldsymbol{u}\|_{2}\bigr{)}^{1/2},$ (11)

where all constants $C$ in (9) – (11) are independent of $\nu$.

We consider the mixed finite element pair known as Hood–Taylor elements [2,
28] $({\boldsymbol{X}}_{h}^{l},Q_{h}^{l-1})$, $l\geq 2$. For these elements a
uniform inf-sup condition is satisfied (see [2]), that is, there exists a
constant $\beta_{\rm is}>0$ independent of the mesh size $h$ such that

$\inf_{q_{h}\in
Q_{h}^{l-1}}\sup_{\boldsymbol{v}_{h}\in{\boldsymbol{X}}_{h}^{l}}\frac{(q_{h},\nabla\cdot\boldsymbol{v}_{h})}{\|\boldsymbol{v}_{h}\|_{1}\|q_{h}\|_{L^{2}/{\mathbb{R}}}}\geq\beta_{\rm{is}}.$
(12)

As a direct method, or full order method, we consider a Galerkin method with
grad-div stabilization. The semi-discrete method reads as follows: Find
$(\boldsymbol{u}_{h},p_{h})\in{\boldsymbol{X}}_{h}^{l}\times Q_{h}^{l-1}$ such
that

$\displaystyle\left(\boldsymbol{u}_{h,t},\boldsymbol{v}_{h}\right)+\nu(\nabla\boldsymbol{u}_{h},\nabla\boldsymbol{v}_{h})+b(\boldsymbol{u}_{h},\boldsymbol{u}_{h},\boldsymbol{v}_{h})$
(13)
$\displaystyle-(p_{h},\nabla\cdot\boldsymbol{v}_{h})+\mu(\nabla\cdot\boldsymbol{u}_{h},\nabla\cdot\boldsymbol{v}_{h})$
$\displaystyle=$
$\displaystyle({\boldsymbol{f}},\boldsymbol{v}_{h})\quad\forall\
\boldsymbol{v}_{h}\in{\boldsymbol{X}}_{h}^{l},$
$\displaystyle(\nabla\cdot\boldsymbol{u}_{h},q_{h})$ $\displaystyle=$
$\displaystyle 0\quad\forall\ q_{h}\in Q_{h}^{l-1},$

where $\mu$ is the positive grad-div stabilization parameter.

It is well-known that considering the discretely divergence-free space
$\boldsymbol{V}_{h,l}$, we can remove the pressure from (13) since
$\boldsymbol{u}_{h}\in\boldsymbol{V}_{h,l}$ satisfies

$\left(\boldsymbol{u}_{h,t},\boldsymbol{v}_{h}\right)+\nu(\nabla\boldsymbol{u}_{h},\nabla\boldsymbol{v}_{h})+b(\boldsymbol{u}_{h},\boldsymbol{u}_{h},\boldsymbol{v}_{h})+\mu(\nabla\cdot\boldsymbol{u}_{h},\nabla\cdot\boldsymbol{v}_{h})=({\boldsymbol{f}},\boldsymbol{v}_{h}),\quad\forall\
\boldsymbol{v}_{h}\in{\boldsymbol{V}}_{h,l}.$ (14)

For this method the following bound holds, see [8],

$\|\boldsymbol{u}(\cdot,t)-\boldsymbol{u}_{h}(\cdot,t)\|_{0}+h\|\boldsymbol{u}(\cdot,t)-\boldsymbol{u}_{h}(\cdot,t)\|_{1}\leq
C(\boldsymbol{u},p,l+1)h^{l},\quad t\in(0,T].$ (15)

where the constant $C(\boldsymbol{u},p,l+1)$ does not explicitly depend on
inverse powers of $\nu$. Actually, only the first term on the left-hand side
of (15) is considered in [8] but the estimate for the second term follows then
from (8) and the inverse inequality (6). Numerical studies presented in [9]
show that the estimate from [8] is sharp. From the error analysis performed in
[8], it can be seen that an optimal order for the error of the velocity
gradient, in $L^{2}(0,T;L^{2}(\Omega)^{d})$, is obtained with a constant that
depends on $\nu^{-1}$, i.e., this estimate is not robust.

##  3 Proper orthogonal decomposition

We consider a POD method. Let us fix $T>0$ and $M>0$ and take $\Delta t=T/M$.
For $N=M+1$ we define the following space

${\cal\boldsymbol{U}}=\mbox{span}\left\\{\boldsymbol{y}_{h}^{1},\boldsymbol{y}_{h}^{2},\ldots,\boldsymbol{y}_{h}^{N}\right\\}=\mbox{span}\left\\{\sqrt{N}\boldsymbol{u}_{h}^{0},\tau\boldsymbol{u}_{h,t}^{1},\ldots,\tau\boldsymbol{u}_{h,t}^{M}\right\\},$

where we use the notation
$\boldsymbol{u}_{h}^{j}=\boldsymbol{u}_{h}(\cdot,t_{j})$ for the
approximations at time instance $t_{j}=j\Delta t$ and
$\boldsymbol{u}_{h,t}^{j}=\boldsymbol{u}_{h,t}(\cdot,t_{j})$ are the snapshots
of the temporal derivatives. The factor $\tau$ in front of the temporal
derivatives is a time scale and it makes the snapshots dimensionally correct,
i.e., all snapshots possess the same physical unit. Let $d_{v}$ be the
dimension of $\cal\boldsymbol{U}$.

The correlation matrix corresponding to the snapshots is given by
$K_{v}=((k_{i,j}^{v}))\in{\mathbb{R}}^{N\times N}$, with the entries

$k_{i,j}^{v}=\frac{1}{N}\left(\boldsymbol{y}_{h}^{i},\boldsymbol{y}_{h}^{j}\right)_{X},\quad
i,j=1,\ldots,N,$

and $(\cdot,\cdot)_{X}$ is the inner product in $X$, which is either
$L^{2}(\Omega)^{d}$ or $H_{0}^{1}(\Omega)^{d}$. Following [21], we denote by
$\lambda_{1}\geq\lambda_{2},\ldots\geq\lambda_{d_{v}}>0$ the positive
eigenvalues of $K_{v}$ and by
$\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{d_{v}}\in{\mathbb{R}}^{N}$
associated eigenvectors of Euclidean norm $1$. Then, the (orthonormal) POD
basis functions of $\cal\boldsymbol{U}$ are given by

$\boldsymbol{\varphi}_{k}=\frac{1}{\sqrt{N}}\frac{1}{\sqrt{\lambda_{k}}}\sum_{j=1}^{N}v_{k}^{j}\boldsymbol{y}_{h}^{j},$
(16)

where $v_{k}^{j}$ is the $j$-th component of the eigenvector
$\boldsymbol{v}_{k}$. The following error estimate is known from [21,
Proposition 1]

$\frac{1}{N}\sum_{j=1}^{N}\left\|\boldsymbol{y}_{h}^{j}-\sum_{k=1}^{r}(\boldsymbol{y}_{h}^{j},\boldsymbol{\varphi}_{k})_{X}\boldsymbol{\varphi}_{k}\right\|_{X}^{2}=\sum_{k=r+1}^{d_{v}}\lambda_{k},$

from which one can deduce

$\left\|\boldsymbol{u}_{h}^{0}-\sum_{k=1}^{r}(\boldsymbol{u}_{h}^{0},\boldsymbol{\varphi}_{k})_{X}\boldsymbol{\varphi}_{k}\right\|_{X}^{2}+\frac{\tau^{2}}{M+1}\sum_{j=1}^{M}\left\|\boldsymbol{u}_{h,t}^{j}-\sum_{k=1}^{r}(\boldsymbol{u}_{h,t}^{j},\boldsymbol{\varphi}_{k})_{X}\boldsymbol{\varphi}_{k}\right\|_{X}^{2}=\sum_{k=r+1}^{d_{v}}\lambda_{k}.$
(17)

We will denote the mass matrix of the POD basis by
$M^{v}=((m_{i,j}^{v}))\in{\mathbb{R}}^{d_{v}\times d_{v}}$, where
$m_{i,j}^{v}=(\boldsymbol{\varphi}_{j},\boldsymbol{\varphi}_{i})_{X}$. In the
case $X=H_{0}^{1}(\Omega)^{d}$, for any
$\boldsymbol{v}\in{\cal\boldsymbol{U}}$, the following inverse inequality
holds, see [21, Lemma 2],

$\|\nabla\boldsymbol{v}\|_{0}\leq\sqrt{\|(M^{v})^{-1}\|_{2}}\|\boldsymbol{v}\|_{0}.$
(18)

The stiffness matrix of the POD basis is given by
$S^{v}=((s_{i,j}^{v}))\in{\mathbb{R}}^{d_{v}\times d_{v}}$, with the entries
$s_{i,j}^{v}=(\nabla\boldsymbol{\varphi}_{j},\nabla\boldsymbol{\varphi}_{i})_{X}$.
If $X=L^{2}(\Omega)^{d}$, the following inequality holds for for all
$\boldsymbol{v}\in{\cal\boldsymbol{U}}$, see [21, Lemma 2],

$\|\nabla\boldsymbol{v}\|_{0}\leq\sqrt{\|S^{v}\|_{2}}\|\boldsymbol{v}\|_{0}.$
(19)

The following lemma will be the basis for proving pointwise in time estimates.
We follow [22, Lemma 3.3] (see also [19, Lemma 3.6]).

###### Lemma 1

Let $T>0$, $\Delta t=T/M$, $t^{n}=n\Delta t$, $n=0,1,\ldots M$, let $X$ be a
normed space, $\boldsymbol{z}(t,\boldsymbol{x})$ be a function defined for
$t\in[0,T]$, $\boldsymbol{x}\in\Omega$ with
$\boldsymbol{z}^{n}=\boldsymbol{z}(t_{n},\cdot)\in X$,
$\boldsymbol{z}_{t}^{n}=\boldsymbol{z}_{t}(t_{n},\cdot)\in X$ and
$\boldsymbol{z}_{tt}(t,\cdot)\in X$ for any $t\in[0,T]$. Then, the following
estimate holds

$\max_{0\leq k\leq
N}\|\boldsymbol{z}^{k}\|_{X}^{2}\leq{3}\|\boldsymbol{z}^{0}\|_{X}^{2}+\frac{3T^{2}}{M}\sum_{n=1}^{M}\|\boldsymbol{z}_{t}^{n}\|_{X}^{2}+\frac{4T}{3}(\Delta
t)^{2}\int_{0}^{T}\|\boldsymbol{z}_{tt}(s)\|_{X}^{2}\ ds.$ (20)

###### Proof:

For each $k$ we have

$\boldsymbol{z}^{k}=\boldsymbol{z}^{0}+\int_{t_{0}}^{t_{k}}\boldsymbol{z}_{t}\
ds.$ (21)

Adding and subtracting terms leads to

$\boldsymbol{z}^{k}=\boldsymbol{z}^{0}+\int_{t_{0}}^{t_{k}}\boldsymbol{z}_{t}\
ds=\boldsymbol{z}^{0}+\sum_{n=0}^{k-1}\Delta
t\boldsymbol{z}_{t}^{n}+\sum_{n=0}^{k-1}\int_{t_{n}}^{t_{n+1}}\left(\boldsymbol{z}_{t}(s)-\boldsymbol{z}_{t}(t_{n})\right)\
ds.$ (22)

To bound the last term on the right-hand side, we first notice that

$\boldsymbol{z}_{t}(s)-\boldsymbol{z}_{t}(t_{n})=\int_{t_{n}}^{s}\boldsymbol{z}_{tt}(\sigma)\
d\sigma.$

With the Cauchy–Schwarz inequality, it follows that

$\displaystyle\left\|\int_{t_{n}}^{t_{n+1}}(\boldsymbol{z}_{t}(s)-\boldsymbol{z}_{t}(t_{n}))\
ds\right\|_{X}$ $\displaystyle\leq$
$\displaystyle\int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{z}_{t}(s)-\boldsymbol{z}_{t}(t_{n})\right\|_{X}\
ds$ $\displaystyle\leq$
$\displaystyle\int_{t_{n}}^{t_{n+1}}(s-t_{n})^{1/2}\left(\int_{t_{n}}^{s}\left\|\boldsymbol{z}_{tt}(\sigma)\right\|_{X}^{2}\
d\sigma\right)^{1/2}\ ds$ $\displaystyle\leq$ $\displaystyle\frac{2}{3}(\Delta
t)^{3/2}\left(\int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{z}_{tt}(s)\right\|_{X}^{2}\
ds\right)^{1/2}.$

Consequently, for the last term on the right-hand side of (22), we obtain,
using the Cauchy–Schwarz inequality for sums,

$\displaystyle\left\|\sum_{n=0}^{k-1}\left(\int_{t_{n}}^{t_{n+1}}(\boldsymbol{z}_{t}(s)-\boldsymbol{z}_{t}(t_{n}))\
ds\right)\right\|_{X}$ $\displaystyle\leq$
$\displaystyle\sum_{n=0}^{k-1}\left\|\int_{t_{n}}^{t_{n+1}}(\boldsymbol{z}_{t}(s)-\boldsymbol{z}_{t}(t_{n}))\
ds\right\|_{X}$ $\displaystyle\leq$
$\displaystyle\frac{2}{3}\sum_{n=0}^{k-1}\left[(\Delta
t)^{3/2}\left(\int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{z}_{tt}(s)\right\|_{X}^{2}\
ds\right)^{1/2}\right]$ $\displaystyle\leq$
$\displaystyle\frac{2}{3}\left(\sum_{n=0}^{k-1}(\Delta
t)^{3}\right)^{1/2}\left(\int_{t_{0}}^{t_{k}}\left\|\boldsymbol{z}_{tt}(t)\right\|_{X}^{2}\
ds\right)^{1/2}$ $\displaystyle\leq$ $\displaystyle\frac{2}{3}T^{1/2}\Delta
t\left(\int_{t_{0}}^{t_{k}}\left\|\boldsymbol{z}_{tt}(t)\right\|_{X}^{2}\
ds\right)^{1/2}.$

Taking norms in (22) gives the estimate

$\displaystyle\|\boldsymbol{z}^{k}\|_{X}$ $\displaystyle\leq$
$\displaystyle\|\boldsymbol{z}^{0}\|_{X}+\sum_{n=l}^{k-1}\Delta
t\|\boldsymbol{z}_{t}^{n}\|_{X}+\frac{2}{3}T^{1/2}\Delta
t\left(\int_{t_{l}}^{t_{k}}\left\|\boldsymbol{z}_{tt}(t)\right\|_{X}^{2}\
ds\right)^{1/2}$ $\displaystyle\leq$
$\displaystyle\|\boldsymbol{z}^{0}\|_{X}+T^{1/2}(\Delta
t)^{1/2}\left(\sum_{n=1}^{M}\|\boldsymbol{z}_{t}^{n}\|_{X}^{2}\right)^{1/2}+\frac{2}{3}T^{1/2}\Delta
t\left(\int_{0}^{T}\left\|\boldsymbol{z}_{tt}(t)\right\|_{X}^{2}\
ds\right)^{1/2},$

from which we conclude (20). $\Box$

The above lemma also holds true changing the initial value by the mean value.

###### Lemma 2

With the assumptions of Lemma 1, it holds

$\max_{0\leq k\leq
N}\|\boldsymbol{z}^{k}\|_{X}^{2}\leq{3}\|\overline{\boldsymbol{z}}\|_{X}^{2}+\frac{12T^{2}}{M}\sum_{n=1}^{M}\|\boldsymbol{z}_{t}^{n}\|_{X}^{2}+\frac{16T}{3}(\Delta
t)^{2}\int_{0}^{T}\|\boldsymbol{z}_{tt}(s)\|_{X}^{2}\ ds,$ (23)

where
$\overline{\boldsymbol{z}}=\frac{1}{M+1}\sum_{j=0}^{M}\boldsymbol{z}^{j}$.

###### Proof:

We first observe that (21) gives

$\overline{\boldsymbol{z}}=\boldsymbol{z}^{0}+\frac{1}{M+1}\left\\{\int_{t_{0}}^{t_{1}}\boldsymbol{z}_{t}\
ds+\int_{t_{0}}^{t_{2}}\boldsymbol{z}_{t}\
ds+\ldots+\int_{t_{0}}^{t_{M}}\boldsymbol{z}_{t}\ ds\right\\},$

so that

$\|\boldsymbol{z}^{0}\|_{X}\leq\|\overline{\boldsymbol{z}}\|_{X}+\frac{1}{M+1}\sum_{j=1}^{M}\left\|\int_{t_{0}}^{t_{j}}\boldsymbol{z}_{t}\
ds\right\|_{X}.$

Since, we have obtained in the proof of Lemma 1

$\left\|\int_{t_{0}}^{t_{j}}\boldsymbol{z}_{t}\right\|_{X}\leq T^{1/2}(\Delta
t)^{1/2}\left(\sum_{n=1}^{M}\|\boldsymbol{z}_{t}^{n}\|_{X}^{2}\right)^{1/2}+\frac{2}{3}T^{1/2}\Delta
t\left(\int_{0}^{T}\left\|\boldsymbol{z}_{tt}(t)\right\|_{X}^{2}\
ds\right)^{1/2},$ (24)

it follows that

$\|\boldsymbol{z}^{0}\|_{X}\leq\|\overline{\boldsymbol{z}}\|_{X}+T^{1/2}(\Delta
t)^{1/2}\left(\sum_{n=1}^{M}\|\boldsymbol{z}_{t}^{n}\|_{X}^{2}\right)^{1/2}+\frac{2}{3}T^{1/2}\Delta
t\left(\int_{0}^{T}\left\|\boldsymbol{z}_{tt}(t)\right\|_{X}^{2}\
ds\right)^{1/2}.$ (25)

Now, taking into account (21), (24), and (25), we obtain for any $k$

$\|\boldsymbol{z}^{k}\|_{X}\leq\|\overline{\boldsymbol{z}}\|_{X}+2T^{1/2}(\Delta
t)^{1/2}\left(\sum_{n=1}^{M}\|\boldsymbol{z}_{t}^{n}\|_{X}^{2}\right)^{1/2}+\frac{4}{3}T^{1/2}\Delta
t\left(\int_{0}^{T}\left\|\boldsymbol{z}_{tt}(t)\right\|_{X}^{2}\
ds\right)^{1/2},$

from which we conclude (23). $\Box$

###### Remark 1

In the sequel we will apply Lemma 1 and assume that we have in the set of
snapshots $\boldsymbol{y}_{h}^{1}=\sqrt{N}\boldsymbol{u}_{h}^{0}$. However,
applying Lemma 2 instead of Lemma 1, we can substitute the first snapshot by
$\boldsymbol{y}_{h}^{1}=\sqrt{N}\overline{\boldsymbol{u}}_{h}$, where
$\overline{\boldsymbol{u}}_{h}$ is the mean value
$\overline{\boldsymbol{u}}_{h}=\frac{1}{M+1}\sum_{j=0}^{M}\boldsymbol{u}_{h}^{j}$,
to obtain the same results. The set of snapshots in this case is

${\cal\boldsymbol{U}}=\mbox{span}\left\\{\boldsymbol{y}_{h}^{1},\boldsymbol{y}_{h}^{2},\ldots,\boldsymbol{y}_{h}^{N}\right\\}=\mbox{span}\left\\{\sqrt{N}\overline{\boldsymbol{u}}_{h},\tau\boldsymbol{u}_{h,t}^{1},\ldots,\tau\boldsymbol{u}_{h,t}^{M}\right\\},$

It is standard in numerical simulations to subtract the mean and apply the POD
method to the fluctuations, see [3], which by definition have zero mean. Then,
the new approach we propose, in which the first snapshot is a weighted mean
and the rest are weighted time derivatives, has the advantage that applying
the POD method to the fluctuations only includes snapshots for the
approximations to the time derivatives (since the mean is zero). This
procedure was applied in our numerical studies.

In the sequel, we will denote by
${\cal\boldsymbol{U}}^{r}=\mbox{span}\\{\boldsymbol{\varphi}_{1},\boldsymbol{\varphi}_{2},\ldots,\boldsymbol{\varphi}_{r}\\}$,
$1\leq r\leq d_{v},$ and by $P_{r}^{v}\ :\
{\boldsymbol{X}}_{h}^{l}\to{\cal\boldsymbol{U}}^{r}$, the $X$-orthogonal
projection onto ${\cal\boldsymbol{U}}^{r}$.

Taking $\boldsymbol{z}=P_{r}^{v}\boldsymbol{u}_{h}-\boldsymbol{u}_{h}$ in (20)
and applying (17) yields

$\displaystyle\max_{0\leq n\leq
M}\|P_{r}^{v}\boldsymbol{u}_{h}^{n}-\boldsymbol{u}_{h}^{n}\|_{X}^{2}$
$\displaystyle\leq$
$\displaystyle\left(3+6\frac{T^{2}}{\tau^{2}}\right)\sum_{k={r+1}}^{d_{v}}\lambda_{k}$
$\displaystyle+\frac{4T}{3}(\Delta
t)^{2}\int_{0}^{T}\|P_{r}^{v}\boldsymbol{u}_{h,tt}(s)-\boldsymbol{u}_{h,tt}(s)\|_{X}^{2}\
ds.$

Let us bound the second term on the right-hand side above. Using the triangle
inequality and taking into account that
$\|P_{r}^{v}\boldsymbol{w}\|_{X}^{2}\leq\|\boldsymbol{w}\|_{X}^{2}$, we
conclude that

$\|P_{r}^{v}\boldsymbol{u}_{h,tt}(s)-\boldsymbol{u}_{h,tt}(s)\|_{X}^{2}\leq
2\|P_{r}^{v}\boldsymbol{u}_{h,tt}(s)\|_{X}^{2}+2\|\boldsymbol{u}_{h,tt}(s)\|_{X}^{2}\leq
4\|\boldsymbol{u}_{h,tt}(s)\|_{X}^{2},$

so that

$\max_{0\leq n\leq
M}\|P_{r}^{v}\boldsymbol{u}_{h}^{n}-\boldsymbol{u}_{h}^{n}\|_{X}^{2}\leq
C_{X}^{2}:=\left(3+6\frac{T^{2}}{\tau^{2}}\right)\sum_{k={r+1}}^{d_{v}}\lambda_{k}+\frac{16T}{3}(\Delta
t)^{2}\int_{0}^{T}\|\boldsymbol{u}_{h,tt}(s)\|_{X}^{2}\ ds.$ (26)

Finally, from (26) we obtain

$\displaystyle\frac{1}{M}\sum_{j=1}^{M}\left\|\boldsymbol{u}_{h}^{j}-\sum_{k=1}^{r}(\boldsymbol{u}_{h}^{j},\boldsymbol{\varphi}_{k})_{X}\boldsymbol{\varphi}_{k}\right\|_{X}^{2}\leq
C_{X}^{2}.$ (27)

The integral term in (26) and (27) is bounded in the appendix. In case one
uses difference quotients instead of Galerkin time derivatives in the set of
snapshots this term does not appear.

The remainder of this section is devoted to proving a priori bounds for
$P_{r}^{v}\boldsymbol{u}_{h}^{j}$, $j=0,\ldots,M$. The a priori bounds for the
orthogonal projection $P_{r}^{v}\boldsymbol{u}_{h}^{j}$, $j=0,\cdots,M$, are
obtained from a priori bounds for the Galerkin approximation
$\boldsymbol{u}_{h}^{j}$, $j=0,\cdots,M$. We argue as in [23]. Then, we start
getting a priori bounds for the stabilized approximation
$\boldsymbol{u}_{h}^{n}$. We follow the same arguments we introduced in [11].

We start with the $L^{\infty}$ norm. Using (6), (10), (15), and (8), we get

$\displaystyle\max_{0\leq j\leq M}\|\boldsymbol{u}_{h}^{j}\|_{\infty}$
$\displaystyle\leq$
$\displaystyle\|\boldsymbol{u}_{h}^{j}-\boldsymbol{s}_{h}^{m}(\cdot,t_{j})\|_{\infty}+\|\boldsymbol{s}_{h}^{m}(\cdot,t_{j})\|_{\infty}$
$\displaystyle\leq$ $\displaystyle
Ch^{-d/2}\|\boldsymbol{u}_{h}^{j}-\boldsymbol{s}_{h}^{m}(\cdot,t_{j})\|_{0}+C(\|\boldsymbol{u}^{j}\|_{d-2}\|\boldsymbol{u}^{j}\|_{2})^{1/2}$
$\displaystyle\leq$ $\displaystyle
Ch^{-d/2}C(\boldsymbol{u},p,3)h^{2}+C(\|\boldsymbol{u}^{j}\|_{d-2}\|\boldsymbol{u}^{j}\|_{2})^{1/2}$
$\displaystyle\leq$ $\displaystyle C_{\boldsymbol{u},{\rm
inf}}:=C\left(C(\boldsymbol{u},p,3)+\left(\|\boldsymbol{u}\|_{L^{\infty}(H^{d-2})}\|\boldsymbol{u}\|_{L^{\infty}(H^{2})}\right)^{1/2}\right).$

The $L^{\infty}$ norm of the gradient is bounded in a similar way. Using (6),
(9), (15), and (8), we obtain

$\displaystyle\max_{0\leq j\leq M}\|\nabla\boldsymbol{u}_{h}^{j}\|_{\infty}$
$\displaystyle\leq$
$\displaystyle\|\nabla\boldsymbol{u}_{h}^{j}-\nabla\boldsymbol{s}_{h}^{m}(\cdot,t_{j})\|_{\infty}+\|\nabla\boldsymbol{s}_{h}^{m}(\cdot,t_{j})\|_{\infty}$
(29) $\displaystyle\leq$ $\displaystyle
Ch^{-d/2}\|\boldsymbol{u}_{h}^{j}-\boldsymbol{s}_{h}^{m}(\cdot,t_{j})\|_{1}+C\|\nabla\boldsymbol{u}^{j}\|_{\infty}$
$\displaystyle\leq$ $\displaystyle
Ch^{-d/2}C(\boldsymbol{u},p,d+1)h^{d-1}+C\|\nabla\boldsymbol{u}^{j}\|_{\infty}$
$\displaystyle\leq$ $\displaystyle C_{\boldsymbol{u},1,{\rm
inf}}:=C\left(C(\boldsymbol{u},p,d+1)+\|\nabla\boldsymbol{u}\|_{L^{\infty}(L^{\infty})}\right).$

Note that the estimate for $d=3$ requires the use cubic elements for the
velocity. Finally, applying (6), (11), (15), and (8) leads to the following
bound of the $L^{2d/(d-1)}$ norm of the velocity gradient

$\displaystyle\max_{0\leq j\leq
M}\|\nabla\boldsymbol{u}_{h}^{j}\|_{L^{2d/(d-1)}}$ $\displaystyle\leq$
$\displaystyle\|\nabla(\boldsymbol{u}_{h}^{j}-\boldsymbol{s}_{h}^{m}(\cdot,t_{j}))\|_{L^{2d/(d-1)}}+\|\nabla\boldsymbol{s}_{h}^{m}(\cdot,t_{j})\|_{L^{2d/(d-1)}}$
$\displaystyle\leq$ $\displaystyle
Ch^{-1/2}\|\boldsymbol{u}_{h}^{j}-\boldsymbol{s}_{h}^{m}(\cdot,t_{j})\|_{1}+C\bigl{(}\|\boldsymbol{u}\|_{1}\|\boldsymbol{u}\|_{2}\bigr{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle
Ch^{-1/2}C(\boldsymbol{u},p,3)h+C\bigl{(}\|\boldsymbol{u}^{j}\|_{1}\|\boldsymbol{u}^{j}\|_{2}\bigr{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle C_{\boldsymbol{u},{\rm
ld}}:=C\left(C(\boldsymbol{u},p,3)+\left(\|\boldsymbol{u}\|_{L^{\infty}(H^{1})}\|\boldsymbol{u}\|_{L^{\infty}(H^{2})}\right)^{1/2}\right).$

Now, we prove a priori bounds in the same norms for

$P_{r}^{v}\boldsymbol{u}_{h}^{j}=(P_{r}^{v}\boldsymbol{u}_{h}^{j}-\boldsymbol{u}_{h}^{j})+\boldsymbol{u}_{h}^{j}.$

Since we have already proved error bounds for the second term on the right-
hand side, we only need to bound the first one. First, we consider the case
$X=L^{2}(\Omega)^{d}$. From (26) we get for $j=0,\ldots,M$,

$\displaystyle\|\boldsymbol{u}_{h}^{j}-P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{0}\leq
C_{L^{2}}.$ (31)

Applying the inverse inequality (6), (3), and (31) gives

$\|P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}\leq\|\boldsymbol{u}_{h}^{j}\|_{\infty}+c_{\rm
inv}h^{-d/2}\|\boldsymbol{u}_{h}^{j}-P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{0}.$

Then

$\displaystyle C_{\rm inf}:=\max_{0\leq j\leq
M}\|P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}\leq C_{\boldsymbol{u},{\rm
inf}}+c_{\rm inv}h^{-d/2}C_{L^{2}},$ (32)

and

$K_{\rm inf}:=\Delta
t\sum_{j=0}^{M}\|P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}^{2}\leq TC_{\rm
inf}^{2}.$ (33)

Now, we observe that from (19) and (31), we get

$\|\nabla(\boldsymbol{u}_{h}^{j}-P_{r}^{v}\boldsymbol{u}_{h}^{j})\|_{0}\leq\|S^{v}\|_{2}^{1/2}C_{L^{2}}.$

Applying this inequality together with the inverse inequality (6) and (29), we
obtain

$C_{1,\rm inf}:=\max_{0\leq j\leq M}\|\nabla
P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}\leq C_{\boldsymbol{u},1,{\rm
inf}}+c_{\rm inv}h^{-d/2}\|S^{v}\|_{2}^{1/2}C_{L^{2}},$ (34)

and then, as before,

$\displaystyle K_{1,{\rm inf}}:=\Delta t\sum_{j=0}^{M}\|\nabla
P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}\leq TC_{1,\rm inf}.$ (35)

Finally, arguing in the same way but applying (3) instead of (29), we find

$C_{\rm ld}:=\max_{0\leq j\leq M}\|\nabla
P_{r}^{v}\boldsymbol{u}^{j}\|_{L^{2d/(d-1)}}\leq C_{\boldsymbol{u},{\rm
ld}}+c_{\rm inv}h^{-1/2}\|S^{v}\|_{2}^{1/2}C_{L^{2}}.$ (36)

The case $X=H_{0}^{1}(\Omega)^{d}$ is simpler. As before, from (26), we get
for $j=0,\ldots,M$,

$\|\nabla(\boldsymbol{u}_{h}^{j}-P_{r}^{v}\boldsymbol{u}_{h}^{j})\|_{0}\leq
C_{H^{1}}.$

Applying Poincaré’s inequality yields

$\|\boldsymbol{u}_{h}^{j}-P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{0}\leq
C_{p}C_{H^{1}}.$

Arguing as before, we obtain

$\displaystyle C_{\rm inf}:=\max_{0\leq j\leq
M}\|P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}$ $\displaystyle\leq$
$\displaystyle C_{\boldsymbol{u},{\rm inf}}+c_{\rm
inv}C_{p}h^{-d/2}C_{H^{1}},$ (37) $\displaystyle C_{1,\rm inf}:=\max_{0\leq
j\leq M}\|\nabla P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}$
$\displaystyle\leq$ $\displaystyle C_{\boldsymbol{u},1,{\rm inf}}+c_{\rm
inv}h^{-d/2}C_{H^{1}},$ (38) $\displaystyle C_{\rm ld}:=\max_{0\leq j\leq
M}\|\nabla P_{r}^{v}\boldsymbol{u}^{j}\|_{L^{2d/(d-1)}}$ $\displaystyle\leq$
$\displaystyle C_{\boldsymbol{u},{\rm ld}}+c_{\rm inv}h^{-1/2}C_{H^{1}}.$ (39)

From (37), it follows that

$K_{\rm inf}:=\Delta
t\sum_{j=0}^{M}\|P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}^{2}\leq TC_{\rm
inf}^{2}$ (40)

and from (38) that

$\displaystyle K_{1,{\rm inf}}$ $\displaystyle:=$ $\displaystyle\Delta
t\sum_{j=0}^{M}\|\nabla P_{r}^{v}\boldsymbol{u}_{h}^{j}\|_{\infty}\leq
TC_{1,\rm inf}.$ (41)

###### Remark 2

Let us observe that the factor $\|S^{v}\|_{2}^{1/2}$ appearing in (34), (35),
and (36) does not appear in (38), (39), and (41). Hence, for a comparable
value of $\sqrt{\lambda_{r+1}}$ in the first and the second bounds, the second
bounds, i.e., those corresponding to the case $X=H_{0}^{1}(\Omega)^{d}$, are
smaller. $\Box$

##  4 The POD-ROM method

###  4.1 The case $X=H_{0}^{1}(\Omega)^{d}$

We now consider the grad-div POD-ROM model. For the sake of simplifying the
error analysis, we use the implicit Euler method as time integrator: For
$n\geq 1$, find $\boldsymbol{u}_{r}^{n}\in{\cal\boldsymbol{U}}^{r}$ such that
for all $\boldsymbol{\varphi}\in{\cal\boldsymbol{U}}^{r}$

$\left(\frac{\boldsymbol{u}_{r}^{n}-\boldsymbol{u}_{r}^{n-1}}{\Delta
t},\boldsymbol{\varphi}\right)+\nu(\nabla\boldsymbol{u}_{r}^{n},\nabla\boldsymbol{\varphi})+b_{h}(\boldsymbol{u}_{r}^{n},\boldsymbol{u}_{r}^{n},\boldsymbol{\varphi})+\mu(\nabla\cdot\boldsymbol{u}_{r}^{n},\nabla\cdot\boldsymbol{\varphi})=(\boldsymbol{f}^{n},\boldsymbol{\varphi}).$
(42)

Taking $t=t_{j}$ in (13) and considering the second equation we get
$\boldsymbol{u}_{h}^{j}\in\boldsymbol{V}_{h,l}$. Differentiating the second
equation in (13) with respect to $t$ and taking again $t=t_{j}$, we also
obtain $\boldsymbol{u}_{h,t}^{j}\in\boldsymbol{V}_{h,l}$. As a consequence, we
observe that ${\cal\boldsymbol{U}}^{r}\subset\boldsymbol{V}_{h,l}$ so that
$\boldsymbol{u}_{r}^{n}$ belongs to the space $\boldsymbol{V}_{h,l}$ of
discretely divergence-free functions.

Choosing $t=t_{n}$ in (14) yields for all
$\boldsymbol{v}_{h}\in{\boldsymbol{V}}_{h,l}$

$\left(\boldsymbol{u}_{h,t}^{n},\boldsymbol{v}_{h}\right)+\nu(\nabla\boldsymbol{u}_{h}^{n},\nabla\boldsymbol{v}_{h})+b(\boldsymbol{u}_{h}^{n},\boldsymbol{u}_{h}^{n},\boldsymbol{v}_{h})+\mu(\nabla\cdot\boldsymbol{u}_{h}^{n},\nabla\cdot\boldsymbol{v}_{h})=({\boldsymbol{f}^{n}},\boldsymbol{v}_{h}).$
(43)

Using the notation
$\boldsymbol{\eta}_{h}^{n}=P_{r}^{v}\boldsymbol{u}_{h}^{n}-\boldsymbol{u}_{h}^{n}$,
a straightforward calculation gives

$\displaystyle\left(\frac{P_{r}^{v}\boldsymbol{u}^{n}_{h}-P_{r}^{v}\boldsymbol{u}^{n-1}_{h}}{\Delta
t},\boldsymbol{\varphi}\right)+\nu(\nabla
P_{r}^{v}\boldsymbol{u}^{n}_{h},\nabla\boldsymbol{\varphi})+b_{h}(P_{r}^{v}\boldsymbol{u}^{n}_{h},P_{r}^{v}\boldsymbol{u}^{n}_{h},\boldsymbol{\varphi})+\mu(\nabla\cdot
P_{r}^{v}\boldsymbol{u}^{n}_{h},\nabla\cdot\boldsymbol{\varphi})$
$\displaystyle=$
$\displaystyle(\boldsymbol{f}^{n},\boldsymbol{\varphi})+\left(\frac{P_{r}^{v}\boldsymbol{u}^{n}_{h}-P_{r}^{v}\boldsymbol{u}^{n-1}_{h}}{\Delta
t}-\boldsymbol{u}_{h,t}^{n},\boldsymbol{\varphi}\right)+\mu(\nabla\cdot\boldsymbol{\eta}_{h}^{n},\nabla\cdot\boldsymbol{\varphi})$
$\displaystyle+b_{h}(P_{r}^{v}\boldsymbol{u}^{n}_{h},P_{r}^{v}\boldsymbol{u}^{n}_{h},\boldsymbol{\varphi})-b_{h}(\boldsymbol{u}^{n}_{h},\boldsymbol{u}^{n}_{h},\boldsymbol{\varphi})\quad\forall\
\boldsymbol{\varphi}\in{\cal\boldsymbol{U}}^{r}.$

Subtracting (4.1) from (42) and denoting by
$\boldsymbol{e}_{r}^{n}=\boldsymbol{u}_{r}^{n}-P_{r}^{v}\boldsymbol{u}_{h}^{n}\in{\cal\boldsymbol{U}}^{r}$
leads to

$\displaystyle\left(\frac{\boldsymbol{e}_{r}^{n}-\boldsymbol{e}_{r}^{n-1}}{\Delta
t},\boldsymbol{\varphi}\right)+\nu(\nabla\boldsymbol{e}_{r}^{n},\nabla\boldsymbol{\varphi})+\mu(\nabla\cdot\boldsymbol{e}_{r}^{n},\nabla\cdot\boldsymbol{\varphi})$
(45)
$\displaystyle+b_{h}(\boldsymbol{u}_{r}^{n},\boldsymbol{u}_{r}^{n},\boldsymbol{\varphi})-b_{h}(P_{r}^{v}\boldsymbol{u}^{n}_{h},P_{r}^{v}\boldsymbol{u}^{n}_{h},\boldsymbol{\varphi})$
$\displaystyle=$
$\displaystyle\left(\boldsymbol{u}_{h,t}^{n}-\frac{P_{r}^{v}\boldsymbol{u}^{n}_{h}-P_{r}^{v}\boldsymbol{u}^{n-1}_{h}}{\Delta
t},\boldsymbol{\varphi}\right)-\mu(\nabla\cdot\boldsymbol{\eta}_{h}^{n},\nabla\cdot\boldsymbol{\varphi})$
$\displaystyle+b_{h}(\boldsymbol{u}^{n}_{h},\boldsymbol{u}^{n}_{h},\boldsymbol{\varphi})-b_{h}(P_{r}^{v}\boldsymbol{u}_{r}^{n},P_{r}^{v}\boldsymbol{u}_{r}^{n},\boldsymbol{\varphi}),\quad\forall\
\boldsymbol{\varphi}\in{\cal\boldsymbol{U}}^{r}.$

Taking now $\boldsymbol{\varphi}=\boldsymbol{e}_{r}^{n}$ yields

$\displaystyle\frac{1}{2\Delta
t}\left(\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}-\|\boldsymbol{e}_{r}^{n-1}\|_{0}^{2}\right)+\nu\|\nabla\boldsymbol{e}_{r}^{n}\|_{0}^{2}+\mu\|\nabla\cdot\boldsymbol{e}_{r}^{n}\|_{0}^{2}$
(46) $\displaystyle\leq$
$\displaystyle\left(\boldsymbol{u}_{h,t}^{n}-\frac{P_{r}^{v}\boldsymbol{u}^{n}_{h}-P_{r}^{v}\boldsymbol{u}^{n-1}_{h}}{\Delta
t},\boldsymbol{e}_{r}^{n}\right)+\big{(}b_{h}(P_{r}^{v}\boldsymbol{u}^{n}_{h},P_{r}^{v}\boldsymbol{u}^{n}_{h},\boldsymbol{e}_{r}^{n})-b_{h}(\boldsymbol{u}_{r}^{n},\boldsymbol{u}_{r}^{n},\boldsymbol{e}_{r}^{n})\big{)}$
$\displaystyle-\mu(\nabla\cdot\boldsymbol{\eta}_{h}^{n},\nabla\cdot\boldsymbol{e}_{r}^{n})+\big{(}b_{h}(\boldsymbol{u}^{n}_{h},\boldsymbol{u}^{n}_{h},\boldsymbol{e}_{r}^{n})-b_{h}(P_{r}^{v}\boldsymbol{u}_{r}^{n},P_{r}^{v}\boldsymbol{u}_{r}^{n},\boldsymbol{e}_{r}^{n})\big{)}$
$\displaystyle=$ $\displaystyle I+II+III+IV.$

The first term on the right-hand side of (46) is estimated by using the
Cauchy–Schwarz and Young inequalities

$\displaystyle|I|\leq\frac{T}{2}\left\|\boldsymbol{u}_{h,t}^{n}-\frac{P_{r}^{v}\boldsymbol{u}^{n}_{h}-P_{r}^{v}\boldsymbol{u}^{n-1}_{h}}{\Delta
t}\right\|_{0}^{2}+\frac{1}{2T}\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}.$ (47)

To bound the second term on the right-hand side of (46), we use, following
[23, Eq. (45)], the skew-symmetry of the trilinear term, Hölder’s inequality,
(37), (38), and Young’s inequality to obtain

$\displaystyle|II|=|b_{h}(\boldsymbol{e}_{r}^{n},P_{r}^{v}\boldsymbol{u}_{h}^{n},\boldsymbol{e}_{r}^{n})|$
$\displaystyle\leq$ $\displaystyle\|\nabla
P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}+\frac{1}{2}\|\nabla\cdot\boldsymbol{e}_{r}^{n}\|_{0}\|P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{r}^{n}\|_{0}$
(48) $\displaystyle\leq$ $\displaystyle\left(\|\nabla
P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}+\frac{\|P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}^{2}}{4\mu}\right)\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}+\frac{\mu}{4}\|\nabla\cdot\boldsymbol{e}_{r}^{n}\|_{0}^{2}.$

For the third term, the application of the Cauchy–Schwarz and Young
inequalities leads to

$\displaystyle|III|\leq\mu\|\nabla\boldsymbol{\eta}_{h}^{n}\|_{0}^{2}+\frac{\mu}{4}\|\nabla\cdot\boldsymbol{e}_{r}^{n}\|_{0}^{2}.$
(49)

For estimating the fourth term, we follow [23, Eq. (50)]. Using Hölder’s
inequality, (3), (3), The Sobolev imbedding (2) with $s=1$ and $q=2$, (4), and
Young’s inequality leads to

$\displaystyle|IV|$ $\displaystyle\leq$
$\displaystyle|b_{h}(P_{r}^{v}\boldsymbol{u}_{h}^{n},\boldsymbol{\eta}_{h}^{n},\boldsymbol{e}_{r}^{n})|+|b_{h}(\boldsymbol{\eta}_{h}^{n},\boldsymbol{u}_{h}^{n},\boldsymbol{e}_{r}^{n})|$
(50) $\displaystyle\leq$
$\displaystyle\|P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}\|\nabla\boldsymbol{\eta}_{h}^{n}\|_{0}\|\boldsymbol{e}_{r}^{n}\|_{0}+\frac{1}{2}\|\nabla\cdot
P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{L^{2d/(d-1)}}\|\boldsymbol{\eta}_{h}^{n}\|_{L^{2d}}\|\boldsymbol{e}_{r}^{n}\|_{0}$
$\displaystyle+\|\boldsymbol{\eta}_{h}^{n}\|_{L^{2d}}\|\nabla\boldsymbol{u}_{h}^{n}\|_{L^{2d/(d-1)}}\|\boldsymbol{e}_{r}^{n}\|_{0}+\frac{1}{2}\|\nabla\cdot\boldsymbol{\eta}_{h}^{n}\|_{0}\|\boldsymbol{u}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{r}^{n}\|_{0}$
$\displaystyle\leq$ $\displaystyle C_{\rm
inf}\|\nabla\boldsymbol{\eta}_{h}^{n}\|_{0}\|\boldsymbol{e}_{r}^{n}\|_{0}+CC_{\rm
ld}\|\nabla\boldsymbol{\eta}_{h}^{n}\|_{0}\|\boldsymbol{e}_{r}^{n}\|_{0}$
$\displaystyle+C\|\nabla\boldsymbol{\eta}_{h}^{n}\|_{0}C_{\boldsymbol{u},\rm
ld}\|\boldsymbol{e}_{r}^{n}\|_{0}+\frac{1}{2}\|\nabla\boldsymbol{\eta}_{h}^{n}\|_{0}C_{\boldsymbol{u},\rm
inf}\|\boldsymbol{e}_{r}^{n}\|_{0}$ $\displaystyle\leq$ $\displaystyle
C_{m}^{2}T\|\nabla\boldsymbol{\eta}_{h}^{n}\|_{0}^{2}+\frac{1}{2T}\|\boldsymbol{e}_{r}^{n}\|_{0}^{2},$

where

$C_{m}=C(C_{\inf}+C_{\boldsymbol{u},\rm inf}+C_{\rm ld}+C_{\boldsymbol{u},\rm
ld}).$ (51)

Inserting (47), (48), (49), and (50) into (46) and adding over the time
instances, we get

$\displaystyle\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}+2\nu\sum_{j=1}^{n}\Delta
t\|\nabla\boldsymbol{e}_{r}^{j}\|_{0}^{2}+\mu\sum_{j=1}^{n}\Delta
t\|\nabla\cdot\boldsymbol{e}_{r}^{j}\|_{0}^{2}$ (52) $\displaystyle\leq$
$\displaystyle\|\boldsymbol{e}_{r}^{0}\|_{0}^{2}+\sum_{j=1}^{n}\Delta
t\left(2\|\nabla
P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}+\frac{\|P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}^{2}}{2\mu}+\frac{2}{T}\right)\|\boldsymbol{e}_{r}^{j}\|_{0}^{2}$
$\displaystyle+2(\mu+C_{m}^{2}T)\sum_{j=1}^{n}\Delta
t\|\nabla\boldsymbol{\eta}_{h}^{j}\|_{0}^{2}+T\sum_{j=1}^{n}\Delta
t\left\|\boldsymbol{u}_{h,t}^{j}-\frac{P_{r}^{v}\boldsymbol{u}^{j}_{h}-P_{r}^{v}\boldsymbol{u}^{j-1}_{h}}{\Delta
t}\right\|_{0}^{2}.$

Adding $\pm P_{r}^{v}\boldsymbol{u}_{h,t}^{j}$ for bounding the last term on
the right-hand side of (52) gives

$\boldsymbol{u}_{h,t}^{j}-\frac{P_{r}^{v}\boldsymbol{u}^{j}_{h}-P_{r}^{v}\boldsymbol{u}^{j-1}_{h}}{\Delta
t}=(\boldsymbol{u}_{h,t}^{j}-P_{r}^{v}\boldsymbol{u}_{h,t}^{j})+\left(P_{r}^{v}\boldsymbol{u}_{h,t}^{j}-\frac{P_{r}^{v}\boldsymbol{u}^{j}_{h}-P_{r}^{v}\boldsymbol{u}^{j-1}_{h}}{\Delta
t}\right),$ (53)

from what follows

$\displaystyle\sum_{j=1}^{n}\Delta
t\left\|\boldsymbol{u}_{h,t}^{j}-\frac{P_{r}^{v}\boldsymbol{u}^{j}_{h}-P_{r}^{v}\boldsymbol{u}^{j-1}_{h}}{\Delta
t}\right\|_{0}^{2}$ $\displaystyle\leq$ $\displaystyle 2\sum_{j=1}^{n}\Delta
t\left\|P_{r}^{v}\boldsymbol{u}_{h,t}^{j}-\boldsymbol{u}_{h,t}^{j}\right\|_{0}^{2}+2\sum_{j=1}^{n}\Delta
t\left\|P_{r}^{v}\left(\boldsymbol{u}_{h,t}^{j}-\frac{\boldsymbol{u}^{j}_{h}-\boldsymbol{u}^{j-1}_{h}}{\Delta
t}\right)\right\|_{0}^{2}.$

To bound the first term on the right-hand side of (4.1), we apply Poincaré’s
inequality (3) and (17), taking into account that $(M+1)/M\leq 2$, to get

$\sum_{j=1}^{n}\Delta
t\left\|P_{r}^{v}\boldsymbol{u}_{h,t}^{j}-\boldsymbol{u}_{h,t}^{j}\right\|_{0}^{2}\leq
C_{p}^{2}\sum_{j=1}^{n}\Delta
t\left\|\nabla\left(P_{r}^{v}\boldsymbol{u}_{h,t}^{j}-\boldsymbol{u}_{h,t}^{j}\right)\right\|_{0}^{2}\leq\frac{2C_{p}^{2}T}{\tau^{2}}\sum_{k=r+1}^{d_{v}}\lambda_{k}.$
(55)

For the second one, we also use Poincaré’s inequality (3) and then notice
that, utilizing a property of a projection in Hilbert spaces, $\|\nabla
P_{r}^{v}\boldsymbol{w}\|_{0}^{2}\leq\|\nabla\boldsymbol{w}\|_{0}^{2}$ for
every $\boldsymbol{w}\in H_{0}^{1}(\Omega)^{d}$, so that

$\sum_{j=1}^{n}\Delta
t\left\|P_{r}^{v}\left(\boldsymbol{u}_{h,t}^{j}-\frac{\boldsymbol{u}^{j}_{h}-\boldsymbol{u}^{j-1}_{h}}{\Delta
t}\right)\right\|_{0}^{2}\leq C_{p}^{2}\sum_{j=1}^{n}\Delta
t\left\|\nabla\left(\boldsymbol{u}_{h,t}^{j}-\frac{\boldsymbol{u}_{h}^{j}-\boldsymbol{u}_{h}^{j-1}}{\Delta
t}\right)\right\|_{0}^{2}.$ (56)

Inserting (55) and (56) into (4.1), we obtain

$\displaystyle\sum_{j=1}^{n}\Delta
t\left\|\boldsymbol{u}_{h,t}^{j}-\frac{P_{r}^{v}\boldsymbol{u}^{j}_{h}-P_{r}^{v}\boldsymbol{u}^{j-1}_{h}}{\Delta
t}\right\|_{0}^{2}$ (57) $\displaystyle\leq$
$\displaystyle\frac{2C_{p}^{2}T}{\tau^{2}}\sum_{k=r+1}^{d_{v}}\lambda_{k}+C_{p}^{2}\sum_{j=1}^{n}\Delta
t\left\|\nabla\left(\boldsymbol{u}_{h,t}^{j}-\frac{\boldsymbol{u}^{j}_{h}-\boldsymbol{u}^{j-1}_{h}}{\Delta
t}\right)\right\|_{0}^{2}.$

For the second term on the right-hand side of (57), a standard argument gives

$C_{p}^{2}\sum_{j=1}^{n}\Delta
t\left\|\nabla\left(\boldsymbol{u}_{h,t}^{j}-\frac{\boldsymbol{u}^{j}_{h}-\boldsymbol{u}^{j-1}_{h}}{\Delta
t}\right)\right\|_{0}^{2}\leq CC_{p}^{2}(\Delta
t)^{2}\int_{0}^{T}\|\nabla(\boldsymbol{u}_{h,tt})\|_{0}^{2}\ ds.$ (58)

Arguing as in [14, Proposition 3.2], one can get an a priori bound for

$\int_{0}^{T}\|\nabla(\boldsymbol{u}_{h,tt})\|_{0}^{2}$ (59)

However, the error analysis in [14] is not valid for high Reynolds numbers. In
the appendix we get an error bound for this term with constants independent of
inverse powers of the viscosity.

Assuming that

$\Delta t\left(2C_{1,\rm inf}+\frac{C_{\rm
inf}^{2}}{2\mu}+\frac{2}{T}\right)\leq\frac{1}{2},$ (60)

denoting by

$C_{u}=2K_{1,{\rm inf}}+\frac{K_{\rm inf}^{2}}{2\mu}+2,$ (61)

applying Gronwall’s Lemma [14, Lemma 5.1], (27), and taking into account that
$(M+1)/M\leq 2$, we obtain from (52)

$\displaystyle\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}+2\nu\sum_{j=1}^{n}\Delta
t\|\nabla\boldsymbol{e}_{r}^{j}\|_{0}^{2}+\mu\sum_{j=1}^{n}\Delta
t\|\nabla\cdot\boldsymbol{e}_{r}^{j}\|_{0}^{2}$ $\displaystyle\leq$
$\displaystyle
e^{2C_{u}}\biggl{(}\|\boldsymbol{e}_{r}^{0}\|_{0}^{2}+\left(2T(\mu+C_{m}^{2}T)(3+6(T/\tau)^{2})+4C_{p}^{2}(T/\tau)^{2}\right)\sum_{k=r+1}^{d_{v}}\lambda_{k}\
$
$\displaystyle{}+\left(CC_{p}^{2}T+2T(\mu+C_{m}^{2}T)\frac{16T}{3}\right)(\Delta
t)^{2}\int_{0}^{T}\|\nabla(\boldsymbol{u}_{h,tt})\|_{0}^{2}\ ds\biggr{)}.$

###### Theorem 4.1

Let $\boldsymbol{u}$ be the velocity in the Navier–Stokes equations (1), let
$\boldsymbol{u}_{r}$ be the grad-div POD stabilized approximation defined in
(42), assume that the solution $(\boldsymbol{u},p)$ of (1) is regular enough
and that (60) holds. Then, the following bound is valid

$\displaystyle\sum_{j=1}^{n}\Delta
t\|\boldsymbol{u}_{r}^{j}-\boldsymbol{u}^{j}\|_{0}^{2}$ $\displaystyle\leq$
$\displaystyle
3Te^{2C_{u}}\left[\|\boldsymbol{e}_{r}^{0}\|_{0}^{2}+\left(2T(\mu+C_{m}^{2}T)(3+6(T/\tau)^{2})+4C_{p}^{2}(T/\tau)^{2}\right)\sum_{k=r+1}^{d_{v}}\lambda_{k}\right.$
$\displaystyle\left.{}+\left(CC_{p}^{2}T+2T(\mu+C_{m}^{2}T)\frac{16T}{3}+16T^{2}C_{p}^{2}\right)\int_{0}^{T}\|\nabla(\boldsymbol{u}_{h,tt})\|_{0}^{2}\
ds\right]$
$\displaystyle+3TC(\boldsymbol{u},p,l)^{2}h^{2l}+3TC_{p}^{2}\left(3+6(T/\tau)^{2}\right)\sum_{k=r+1}^{d_{v}}\lambda_{k}.$

###### Proof:

We have $\sum_{j=1}^{n}\Delta t\|\boldsymbol{e}_{r}^{j}\|_{0}^{2}\leq
T\max_{1\leq j\leq n}\|\boldsymbol{e}_{r}^{j}\|_{0}^{2}$ and

$\sum_{j=1}^{n}\Delta
t\|\boldsymbol{u}_{r}^{j}-\boldsymbol{u}^{j}\|_{0}^{2}\leq
3\left(\sum_{j=1}^{n}\Delta
t\|\boldsymbol{e}_{r}^{j}\|_{0}^{2}+\sum_{j=1}^{n}\Delta
t\|P_{r}^{v}\boldsymbol{u}_{h}^{j}-\boldsymbol{u}_{h}^{j}\|_{0}^{2}+\sum_{j=1}^{n}\Delta
t\|\boldsymbol{u}_{h}^{j}-\boldsymbol{u}^{j}\|_{0}^{2}\right).$

Inserting the estimates (4.1), (27), and (15) leads directly to (4.1). $\Box$

A robust estimate for (59), in the case $l\geq 3$, is derived in the appendix.
Hence, for $l\geq 3$, there is no explicit appearance of inverse powers of the
viscosity coefficient in the error bound (4.1). The technical reason for not
obtaining a robust estimate for $l=2$ with the technique from the appendix is
the gradient in front of $\boldsymbol{u}_{h,tt}$. This gradient was introduced
with the transition from the $L^{2}(\Omega)^{d}$ norm to the corresponding
norm of the gradient in (56) in order to be able to apply the Hilbert space
argument. Note that with the approach presented in the appendix the
boundedness can be shown also for $l=2$, but not the robustness of the bound.

###### Remark 3

With the error decomposition in the proof of Theorem 4.1 and applying (18) or
the inverse inequality (6) to (4.1), one can also prove a robust error bound
for $\sum_{j=1}^{n}\Delta
t\|\nabla(\boldsymbol{u}_{r}^{j}-\boldsymbol{u}^{j})\|_{0}^{2}$. $\Box$

We can apply Lemma 1 to get pointwise estimates with respect to time both in
$L^{2}(\Omega)$ and $H^{1}(\Omega)$. Let us prove pointwise estimates in
$L^{2}(\Omega)$, the argument for proving bounds in $H^{1}(\Omega)$ is the
same. Since

$\|\boldsymbol{u}_{r}^{n}-\boldsymbol{u}^{n}\|_{0}^{2}\leq
3\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}+3\|P_{r}^{v}\boldsymbol{u}_{h}^{n}-\boldsymbol{u}_{h}^{n}\|_{0}^{2}+3\|\boldsymbol{u}_{h}^{n}-\boldsymbol{u}^{n}\|_{0}^{2},$

we can utilize (4.1) and (15) to bound the first and third term on the right-
hand side. For bounding the second term, (26) is utilized. More precisely,
taking into account we are analyzing the case $X=H_{0}^{1}(\Omega)^{d}$ we
apply Poincaré’s inequality (3) to bound
$\|P_{r}^{v}\boldsymbol{u}_{h}^{n}-\boldsymbol{u}_{h}^{n}\|_{0}^{2}\leq
C_{p}^{2}\|\nabla(P_{r}^{v}\boldsymbol{u}_{h}^{n}-\boldsymbol{u}_{h}^{n})\|_{0}^{2}$
and then (26). Collecting the estimates (4.1), (26), (15) proves the following
theorem.

###### Theorem 4.2

Let the assumptions of Theorem 4.1 be satisfied, then the following bound is
valid

$\displaystyle\max_{0\leq n\leq
M}\|\boldsymbol{u}_{r}^{n}-\boldsymbol{u}^{n}\|_{0}^{2}$ $\displaystyle\leq$
$\displaystyle
3e^{2C_{u}}\left[\|\boldsymbol{e}_{r}^{0}\|_{0}^{2}+\left(2T(\mu+C_{m}^{2}T)(3+6(T/\tau)^{2})+4C_{p}^{2}(T/\tau)^{2}\right)\sum_{k=r+1}^{d_{v}}\lambda_{k}\right.$
$\displaystyle+\left.\left(CC_{p}^{2}T+2T(\mu+C_{m}^{2}T)\frac{16T}{3}+16TC_{p}^{2}\right)(\Delta
t)^{2}\int_{0}^{T}\|\nabla(\boldsymbol{u}_{h,tt})\|_{0}^{2}\ ds\right]$
$\displaystyle+3C_{p}^{2}(3+6(T/\tau)^{2})\sum_{k=r+1}^{d_{v}}\lambda_{k}+C(\boldsymbol{u},p,l+1)h^{l},$

where the constants for $l\geq 3$ do not blow up for small viscosity
coefficients since (59) is bounded.

###  4.2 The case $X=L^{2}(\Omega)^{d}$

For the sake of brevity, we are going to mention only the differences with
respect to the analysis for the case $X=H_{0}^{1}(\Omega)^{d}$.

Observing that the orthogonality property of $P_{r}^{v}$ affects now the term
with the approximation of the temporal derivative (instead of the viscous
term) we obtain, instead of (45), the following error equation

$\displaystyle\left(\frac{\boldsymbol{e}_{r}^{n}-\boldsymbol{e}_{r}^{n-1}}{\Delta
t},\boldsymbol{\varphi}\right)+\nu(\nabla\boldsymbol{e}_{r}^{n},\nabla\boldsymbol{\varphi})+\mu(\nabla\cdot\boldsymbol{e}_{r}^{n},\nabla\cdot\boldsymbol{\varphi})$
(64)
$\displaystyle+b_{h}(\boldsymbol{u}_{r}^{n},\boldsymbol{u}_{r}^{n},\boldsymbol{\varphi})-b_{h}(P_{r}^{v}\boldsymbol{u}^{n}_{h},P_{r}^{v}\boldsymbol{u}^{n}_{h},\boldsymbol{\varphi})$
$\displaystyle=$
$\displaystyle\left(\boldsymbol{u}_{h,t}^{n}-\frac{\boldsymbol{u}_{h}^{n}-\boldsymbol{u}_{h}^{n-1}}{\Delta
t},\boldsymbol{\varphi}\right)-\nu(\nabla\boldsymbol{\eta}_{h}^{n},\nabla\boldsymbol{\varphi})-\mu(\nabla\cdot\boldsymbol{\eta}_{h}^{n},\nabla\cdot\boldsymbol{\varphi})$
$\displaystyle+b_{h}(\boldsymbol{u}^{n}_{h},\boldsymbol{u}^{n}_{h},\boldsymbol{\varphi})-b_{h}(P_{r}^{v}\boldsymbol{u}^{n}_{h},P_{r}^{v}\boldsymbol{u}^{n}_{h},\boldsymbol{\varphi})\quad\forall\
\boldsymbol{\varphi}\in{\cal\boldsymbol{U}}^{r}.$

Using the same techniques as for the other case, we infer, instead of (52),
that

$\displaystyle\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}+2\nu\sum_{j=1}^{n}\Delta
t\|\nabla\boldsymbol{e}_{r}^{j}\|_{0}^{2}+\mu\sum_{j=1}^{n}\Delta
t\|\nabla\cdot\boldsymbol{e}_{r}^{j}\|_{0}^{2}$ (65) $\displaystyle\leq$
$\displaystyle\|\boldsymbol{e}_{r}^{0}\|_{0}^{2}+\sum_{j=1}^{n}\Delta
t\left(2\|\nabla
P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}+\frac{\|P_{r}^{v}\boldsymbol{u}_{h}^{n}\|_{\infty}^{2}}{2\mu}+\frac{2}{T}\right)\|\boldsymbol{e}_{r}^{j}\|_{0}^{2}$
$\displaystyle+2(\nu+\mu+C_{m}^{2}T)\sum_{j=1}^{n}\Delta
t\|\nabla\boldsymbol{\eta}_{h}^{j}\|_{0}^{2}+CT(\Delta
t)^{2}\int_{0}^{T}\|\boldsymbol{u}_{h,tt}(s)\|_{0}^{2}\ ds.$

From (65) we continue as before, compare [23, Theorem 5.3], and take into
account that applying (19) and (27), it is

$\sum_{j=1}^{n}\Delta t\|\nabla\boldsymbol{\eta}_{h}^{j}\|_{0}^{2}\leq
T\|S^{v}\|_{2}\left((3+6(T/\tau)^{2})\sum_{k=r+1}^{d_{v}}\lambda_{k}+\frac{16T^{2}}{3}(\Delta
t)^{2}\int_{0}^{T}\|\boldsymbol{u}_{h,tt}(s)\|_{0}^{2}\ ds\right).$

Then, instead of (4.1) we conclude

$\displaystyle\|\boldsymbol{e}_{r}^{n}\|_{0}^{2}+2\nu\sum_{j=1}^{n}\Delta
t\|\nabla\boldsymbol{e}_{r}^{j}\|_{0}^{2}+\mu\sum_{j=1}^{n}\Delta
t\|\nabla\cdot\boldsymbol{e}_{r}^{j}\|_{0}^{2}$ (66) $\displaystyle\leq$
$\displaystyle
e^{2C_{u}}\biggl{(}\|\boldsymbol{e}_{r}^{0}\|_{0}^{2}+2T(\nu+\mu+C_{m}^{2}T)\|S^{v}\|_{2}(3+6(T/\tau)^{2})\sum_{k=r+1}^{d_{v}}\lambda_{k}\
$ $\displaystyle{}+(\Delta
t)^{2}\left(CT+2T(\nu+\mu+C_{m}^{2}T)\|S^{v}\|_{2}(16T^{2})/3\right)\int_{0}^{T}\|\boldsymbol{u}_{h,tt}(s)\|_{0}^{2}\
ds\biggr{)}.$

A robust estimate for the following term and $l\geq 2$ is proved in the
appendix

$\int_{0}^{T}\|\boldsymbol{u}_{h,tt}\|_{0}^{2}\ ds.$ (67)

In [23, Section 5], there is an error analysis for the case in which a fully
discrete Galerkin method with the same implicit Euler method is utilized for
the FOM.

###### Remark 4

Let us observe that at this point it seems that the snapshots from the time
derivative are not helpful in the $L^{2}(\Omega)^{d}$ case. However, it turns
out that these snapshots are needed for obtaining pointwise estimates in time.
In fact, we can repeat the arguments of the previous section to get such
estimates for the $L^{2}(\Omega)^{d}$ error, and also for the
$H^{1}(\Omega)^{d}$ error applying the inverse inequality (19), that cannot be
obtained, at least with the same arguments, without adding those snapshots.
$\Box$

###### Remark 5

Comparing the $L^{2}(\Omega)^{d}$ and the $H_{0}^{1}(\Omega)^{d}$ cases, we
can observe that in the $L^{2}(\Omega)^{d}$ case the eigenvalues are
multiplied by $\|S^{v}\|_{2}$ in the error bound (66), which increases the
size of this term, compared with (4.1). In addition, the factor
$\|S^{v}\|_{2}^{1/2}$ appears in the definition (61) of $C_{u}$, which gives a
bigger constant in the exponential term of the error bound, and also in the
assumption for the time steps (60), i.e., the condition for applying
Gronwall’s lemma leads to a severe time step restriction if $\|S^{v}\|_{2}$ is
very large. On the other hand, using $L^{2}(\Omega)^{d}$ as projection space
leads to an integral term in the error bound, (67), that can be bounded in a
robust way for $l\geq 2$, whereas $l\geq 3$ is necessary if the projection
space is $H_{0}^{1}(\Omega)^{d}$, to bound the term (59), at least with the
approach from the appendix. Thus, both approaches possess advantages with
respect to certain aspects of the error analysis. $\Box$

###### Remark 6

A popular pair of finite element spaces that leads to weakly divergence-free
discrete velocity fields is the Scott–Vogelius pair
$({\boldsymbol{X}}_{h}^{l},Q_{h,\mathrm{disc}}^{l-1})$. It is inf-sup stable
for $l\geq d$ on so-called barycentric-refined grids [29]. A favorable feature
of the Scott–Vogelius pair is that it leads to so-called pressure-robust
velocity estimates, i.e., the velocity error bounds do not depend on the
pressure. In particular, an estimate of form (15) was derived in [26] with
$C(\boldsymbol{u},l+1)$. In estimating the integral term as performed in the
appendix, the term $(\sigma_{2},\nabla\cdot\boldsymbol{\varphi}_{h})$ vanishes
if $\boldsymbol{\varphi}_{h}$ is weakly divergence-free, such that the bound
of the integral term does not contain the pressure. Finally, the snapshots are
only from the velocity field, which can be (formally) computed by solving an
equation of type (14), which does not contain the pressure. Hence all terms in
estimates (27) and (17) are independent of the pressure. Applying the same
analysis as for the Taylor–Hood pair of spaces, even with some
simplifications, gives for the Scott–Vogelius pair pressure-robust error
bounds of the same type, under the conditions on the inf-sup stability
mentioned above. $\Box$

##  5 Numerical studies

We now present some results for the well-known benchmark problem defined in
[25]. The domain is given by

$\Omega=(0,2.2)\times(0,0.41)/\left\\{(x,y)\mid(x-0.2)^{2}+(y-0.2)^{2}\leq
0.0025\right\\}.$

On the inflow boundary the velocity is prescribed by

$\boldsymbol{u}(0,y)=\frac{6}{0.41^{2}}\sin\left(\frac{\pi
t}{8}\right)\left(\begin{array}[]{c}y(0.41-y)\\\ 0\end{array}\right)$

and on the outflow boundary $x=2.2$ we set the so-called “do nothing” boundary
condition. On the rest of the boundary the velocity is set
$\boldsymbol{u}={\bf 0}$. It is well known that for $\nu=0.001$ there is a
stable periodic orbit.

We use piecewise quadratic and linear elements for velocity and pressure,
respectively, on the same mesh as in [10], resulting in 27168 degrees of
freedom for the velocity and 3480 for the pressure. We call this grid here the
main grid. The grid obtained from this one by regular refinement (107328
degrees of freedom for the velocity and 13584 for the pressure), called
refined grid henceforth, was used only to compute a reference solution and the
errors of the snapshots. On both grids we computed the periodic orbit by
finding the fixed point of the return map to a Poincaré section (see e.g.,
[24]). The periods, computed with a relative error below $10^{-6}$ were
$T=0.331761$ and $T_{r}=0.331338$ on the main and refined grid, respectively.
We computed snapshots over one period with a spacing of $T/1024$ and
$T_{r}/1024$. Taking as exact the results on the refined grid, the relative
errors of the snapshots are shown in Fig. 1, both in velocity and in
acceleration, and both in the $L^{2}$ and $H^{1}$ norms. It can be seen that
the errors in velocity are around 0.003 in the $L^{2}$ norm and around 0.03 in
the $H^{1}$ norm, and that the acceleration errors are, approximately, eight
times larger than those of the velocity in the $L^{2}$ norm, and four times
larger in the $H^{1}$ norm.

Figure 1: Snapshot errors in the $L^{2}$ norm (left) and in the $H^{1}$ norm
(right).

Fig. 2 shows the first 256 singular values $\sigma_{k}=\sqrt{\lambda_{k}}$
relative to their Euclidean norm, that is,

$\sigma_{k}/\biggl{(}\sum_{j=1}^{N}\sigma_{k}^{2}\biggr{)}^{1/2},$ (68)

where $N$ is the number of snapshots, both when the inner product is that of
$L^{2}$ and $H^{1}$, in the cases where the elements in the data sets are
$\boldsymbol{y}_{h}^{j}=\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$,
$\boldsymbol{y}_{h}^{j}=T\boldsymbol{u}_{u,t}^{j}$, and
$\boldsymbol{y}_{h}^{j}=T\delta_{t}\boldsymbol{u}_{h}^{j}=T(\boldsymbol{u}_{h}^{j}-\boldsymbol{u}_{h}^{j-1})/\Delta
t$.

Figure 2: First 256 singular values of data set relative to their Euclidean
norm (68), for the $L^{2}$ inner product (left) and the $H^{1}$ inner product
(right).

We see that with both inner products, the singular values are slightly larger
when the data set is the time derivatives or their approximation, but already
for $k=50$ their value is considerable smaller than the approximation errors
in Fig. 1. In fact, for the singular values corresponding to the data set
given by the fluctuations
$\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$ (blue line), there are
only $16$ values above $10^{-3}$ in Fig. 2 (left) and $14$ above $10^{-2}$ in
the right plot. Note that $10^{-3}$ and $10^{-2}$ are about one third of the
average errors in Fig. 1 for the $L^{2}$ and the $H^{1}$ norm, so it
reasonable to use a POD basis with no more elements, since although with a
larger basis we may better approximate the elements in the data set, these
have already larger approximation errors. Consequently, in the sequel, we use
a POD basis with 16 elements if the $L^{2}$ inner product is utilized and $14$
in case of the $H^{1}$ inner product.

Next we present in Fig. 3 the projection errors of the snapshots,

$\left\|\boldsymbol{u}_{h}^{j}-P_{r}^{v}\boldsymbol{u}_{h}^{j}\right\|_{i}\bigg{/}\left\|\boldsymbol{u}_{h}^{j}\right\|_{i},\quad
j=1,\ldots,N,\quad i\in\\{0,1\\},$ (69)

for data sets
$\boldsymbol{y}_{h}^{j}=\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$,
$\boldsymbol{y}_{h}^{j}=T\boldsymbol{u}_{u,t}^{j}$, and
$\boldsymbol{y}_{h}^{j}=T\delta_{t}\boldsymbol{u}_{h}^{j}=(\boldsymbol{u}_{h}^{j}-\boldsymbol{u}_{h}^{j-1})/\Delta
t$, when $r=16$ and $r=14$ for POD basis based on $L^{2}$ and $H^{1}$ inner
products, respectively.

Figure 3: Snapshot projection errors (69) for the correlation matrix based on
$L^{2}$ product (left) and $H^{1}$ products (right). Results for data sets
$\boldsymbol{y}_{h}^{j}=\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$
(blue), $\boldsymbol{y}_{h}^{j}=T\boldsymbol{u}_{u,t}^{j}$ (magenta), and
$\boldsymbol{y}_{h}^{j}=T\delta_{t}\boldsymbol{u}_{h}^{j}=(\boldsymbol{u}_{h}^{j}-\boldsymbol{u}_{h}^{j-1})/\Delta
t$ (black).

We notice that the results on the left plot are quite independent of the data
set used, whereas those on the right plot show some advantage for the data set
based on the fluctuations,
$\boldsymbol{y}_{h}^{j}=\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$,
although in all cases the projection errors are well below the upper bounds

$\biggl{(}\sum_{k=r+1}^{N}\lambda_{k}\biggr{)}^{1/2}\bigg{/}\biggl{(}\sum_{k=1}^{N}\lambda_{k}\biggr{)}^{1/2}$

which, recall, are $0.001$ and $0.01$ for POD basis based on $L^{2}$ and
$H^{1}$, respectively.

Fig. 4 depicts the POD-ROM errors

$\left\|\boldsymbol{u}_{r}^{j}-\boldsymbol{u}_{h}^{j}\right\|_{i}\bigg{/}{\left\|\boldsymbol{u}_{h}^{j}\right\|_{i}},\quad
j=1,\ldots,N,\quad i\in\\{0,1\\},$ (70)

with the same color convention as in Fig. 3.

Figure 4: POD-ROM errors (70) for correlation matrix based on $L^{2}$ product
(left) and $H^{1}$ products (right). Results for data sets
$\boldsymbol{y}_{h}^{j}=\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$
(blue), $\boldsymbol{y}_{h}^{j}=T\boldsymbol{u}_{u,t}^{j}$ (magenta) and
$\boldsymbol{y}_{h}^{j}=T\delta_{t}\boldsymbol{u}_{h}^{j}=(\boldsymbol{u}_{h}^{j}-\boldsymbol{u}_{h}^{j-1})/\Delta
t$ (black).

For the time integration, instead of the backward Euler method in (42), we
used the two step BDF formula with fixed step size $\Delta t=T/N$, $N=1024$,
except for the first step, where the backward Euler method was applied. The
parameter for the grad-div term was as in the computation of the snapshots
$\mu=0.01$. We see that the results are very independent of the data set used,
for both inner products, and that the POD-ROM errors are slightly larger than
the corresponding projection errors, except for data sets based on snapshots,
$\boldsymbol{y}_{h}^{j}=\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$
(blue line) and POD basis based on $H^{1}$ inner products (right plot) which
are approximately twice the size of the corresponding errors in Fig. 3.

Finally, lift and drag coefficients are studied. For the FOM, they were
computed as indicated in [16] (see also [10]). For the POD-ROM method they
were computed as in [11]. In Fig. 5 we see the evolution of the drag and lift
coefficients along 12 periods for the reference solution (computed on the
refined grid), the FOM, and the POD-ROM approximation with the fluctions
$\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$ of the snapshots
computed in the first period as data set with the $L^{2}$ inner product (all
other POD-ROM approximations gave similar results).

Figure 5: Drag (left) and lift (right) coefficients for the reference solution
(black), the FOM (magenta), and the POD-ROM approximation with snapshots
$\boldsymbol{u}_{h}^{j}$ as data set and $L^{2}$ inner product.

At first sight, the agreement is excellent. Fig. 6 compares the first and
twelfth periods. There is a phase difference between the black line and the
other two, since, as commented above, the “exact” solution computed on the
refined grid has a slightly different period than that computed on the main
grid. We notice however, that the difference between the other two lines is
hardly altered from the first period to the last one.

Figure 6: Drag (left) and lift (right) coefficients for the reference solution
(black), the FOM (magenta), and the POD-ROM approximation with snapshots
$\boldsymbol{u}_{h}^{j}$ as data set and $L^{2}$ inner product.

In addition, the results for the lift coefficient corresponding to the FOM and
the POD-ROM approximation are on top of each other (in fact the maximum
relative error of this coefficient with respect to that of the FOM in the 12
periods is below $1.7\times 10^{-4}$). Some differences can be seen in the
drag coefficient. However, the maximum relative error with respect the
coefficient computed with the FOM is below 0.0011.

The summary of the results is that there are no significant differences in the
POD-ROM approximation, i.e., they are quite independent of the data set used.
Only some minor differences on the projection errors depending on the data
set, when the $H^{1}$ inner product is used, could be observed in the
computation of the POD basis.

##  6 Conclusions

We analyzed reduced order models for the incompressible Navier–Stokes
equations based on proper orthogonal decomposition methods. The influence of
including approximations to the time derivative in the set of snapshots was
studied. Our set of snapshots is constituted by the approximation to the
velocity at the initial time together with approximations to the time
derivative at different times. The approximation to the velocity at the
initial time can be replaced by the mean value and then only approximations to
the time derivatives are required for applying the POD method to the
fluctuations. The Galerkin time derivative can be replaced by any other
approximation as the standard difference quotient.

We studied the differences between projecting onto $L^{2}$ and $H^{1}$. We
proved that including the Galerkin time derivatives (or the difference
quotients) leads to pointwise estimates for both projections. In the $L^{2}$
case, error bounds can be proved in the no-time-derivatives case (with only
snapshots for the velocity) as shown in the literature, e.g., see [23].
However, the time derivatives approach is useful also in this case to get
pointwise estimates. In the numerical analysis, we utilized the projection of
the continuous-in-time Galerkin approximation since this allows to use in
practice any time integrator for computing the snapshots. It is easy to
include the corresponding error in time in the bounds of the present paper.
Also, different times can be considered to compute the set of snapshots and
the fully discrete POD approximations. Finally, as in one of the methods in
[23], we added grad-div stabilization to the approximations computed with the
FOM and POD methods to be able to prove error bounds in which the constants do
not depend on inverse powers of the viscosity.

In the numerical studies, we compared three different sets of snapshots, for
both inner products, where the elements in the data sets were
$\boldsymbol{y}_{h}^{j}=\boldsymbol{u}_{h}^{j}-\overline{\boldsymbol{u}}_{h}$,
$\boldsymbol{y}_{h}^{j}=T\boldsymbol{u}_{u,t}^{j}$, and
$\boldsymbol{y}_{h}^{j}=T\delta_{t}\boldsymbol{u}_{h}^{j}=T(\boldsymbol{u}_{h}^{j}-\boldsymbol{u}_{h}^{j-1})/\Delta
t$. While the snapshots were computed using only one period, we showed that
very good approximations are obtained to the lift and drag coefficients in a
time interval of 12 periods. In our numerical studies there were no
significant differences between the different procedures. More comprehensive
studies have to further investigate this topic.

Altogether, we think that the rigorous numerical analysis presented in this
paper shows interesting properties and sharp bounds for the different methods.
It also supports the idea of the recent paper [22] which shows that in the
case of the heat equation there is no need to include in the set of snapshots
other than one approximation to the solution at a fixed time. We prove that
for the approximation to the incompressible Navier–Stokes equations the same
situation holds. Moreover, in case of applying the POD method to fluctuations,
as it is standard, only snapshots of the time derivative are needed.

##  Appendix A Robust bounds for the terms with the second temporal derivative

In this section we derive bounds for the terms (67) and (59). As a
consequence, the terms on the right-hand sides of (4.1) and (66) that contain
(59) and (67), respectively, are bounded.

The constants in the bounds below do not blow up as $\nu\to 0$. It will be
assumed that all functions are sufficiently smooth such that the performed
operations are well defined. It is

$\|\boldsymbol{u}_{h,tt}\|_{0}\leq\|\boldsymbol{u}_{h,tt}-\boldsymbol{s}_{h,tt}\|_{0}+\|\boldsymbol{u}_{tt}\|_{0}+\|\boldsymbol{s}_{h,tt}-\boldsymbol{u}_{tt}\|_{0},$
(71)

and likewise for the norm of the gradient. The first term on the right-hand
side will be bounded below. The second term is bounded because the solution is
sufficiently smooth, and then the last term can be bounded with (8), where the
notation $\boldsymbol{s}_{h,tt}$ has to be understood that the modified Stokes
projection is applied to $\boldsymbol{u}_{tt}$.

Denote $\boldsymbol{e}_{h}=\boldsymbol{u}_{h}-\boldsymbol{s}_{h}$. The error
equation for (13) is given by

$\displaystyle\left(\boldsymbol{e}_{h,t},\boldsymbol{\varphi}_{h}\right)+\nu(\nabla\boldsymbol{e}_{h},\nabla\boldsymbol{\varphi}_{h})+b(\boldsymbol{u}_{h},\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})-b(\boldsymbol{u},\boldsymbol{u},\boldsymbol{\varphi}_{h})+\mu(\nabla\cdot\boldsymbol{e}_{h},\nabla\cdot\boldsymbol{\varphi}_{h})$
(72) $\displaystyle=$
$\displaystyle(\sigma_{1},\boldsymbol{\varphi}_{h})+(\sigma_{2},\nabla\cdot\boldsymbol{\varphi}_{h})\quad\forall\
\boldsymbol{\varphi}_{h}\in{\boldsymbol{V}}_{h}^{l},$

with $\sigma_{1}=\boldsymbol{u}_{t}-\boldsymbol{s}_{h,t}$,
$\sigma_{2}=(p-P_{Q_{h}}p)+\mu(\nabla\cdot(\boldsymbol{u}-\boldsymbol{s}_{h})$,
and $P_{Q_{h}}p$ is the best approximation of $p$ in $Q_{h}^{l}$. Note that if
a (discrete) velocity function is in $\boldsymbol{V}$ (or
${\boldsymbol{V}}_{h}^{l}$), then also the temporal derivatives of this
function are in the same space. Hence, taking
$\boldsymbol{\varphi}_{h}=\boldsymbol{e}_{h,t}\in{\boldsymbol{V}}_{h}^{l}$ in
(72) and using the inverse inequality (6) yields

$\displaystyle\|\boldsymbol{e}_{h,t}\|_{0}^{2}+\frac{d}{dt}\frac{1}{2}\nu\|\nabla\boldsymbol{e}_{h}\|_{0}^{2}+\frac{d}{dt}\frac{1}{2}\mu\|\nabla\cdot\boldsymbol{e}_{h}\|_{0}^{2}\leq|b(\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{u},\boldsymbol{e}_{h,t})|$
(73)
$\displaystyle+|b(\boldsymbol{u}_{h},\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{e}_{h,t})|+\|\sigma_{1}\|_{0}\|\boldsymbol{e}_{h,t}\|_{0}+c_{\rm
inv}\|\sigma_{2}\|_{0}h^{-1}\|\boldsymbol{e}_{h,t}\|_{0}.$

The first term on the right-hand side is estimated with Hölder’s inequality,
(15), and Young’s inequality

$\displaystyle|b(\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{u},\boldsymbol{e}_{h,t})|$
$\displaystyle\leq$
$\displaystyle\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{0}\|\nabla\boldsymbol{u}\|_{\infty}\|\boldsymbol{e}_{h,t}\|_{0}+\frac{1}{2}\|\nabla\cdot(\boldsymbol{u}_{h}-\boldsymbol{u})\|_{0}\|\boldsymbol{u}\|_{\infty}\|\boldsymbol{e}_{h,t}\|_{0}$
$\displaystyle\leq$ $\displaystyle
Ch^{2(l-1)}+\frac{1}{8}\|\boldsymbol{e}_{h,t}\|_{0}^{2}.$

For the second term, in addition an estimate from [16, Lemma 6.11], the
inverse inequality (6), and the condition $l\geq 2$ are utilized

$\displaystyle|b(\boldsymbol{u}_{h},\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{e}_{h,t})|$
$\displaystyle\leq$
$\displaystyle|b(\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{e}_{h,t})|+|b(\boldsymbol{u},\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{e}_{h,t})|$
$\displaystyle\leq$ $\displaystyle
C\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{1}^{2}\|\boldsymbol{e}_{h,t}\|_{1}+\|\boldsymbol{u}\|_{\infty}\|\nabla(\boldsymbol{u}_{h}-\boldsymbol{u})\|_{0}\|\boldsymbol{e}_{h,t}\|_{0}$
$\displaystyle\leq$ $\displaystyle
C\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{1}^{2}c_{\rm
inv}h^{-1}\|\boldsymbol{e}_{h,t}\|_{0}+C\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{1}\|\boldsymbol{e}_{h,t}\|_{0}$
$\displaystyle\leq$ $\displaystyle
C\left(h^{4l-6}+h^{2l-2}\right)+\frac{1}{8}\|\boldsymbol{e}_{h,t}\|_{0}^{2}\leq
Ch^{2(l-1)}+\frac{1}{8}\|\boldsymbol{e}_{h,t}\|_{0}^{2}.$

For the last two terms, we obtain, with (8), the $L^{2}(\Omega)$ best
approximation error for the pressure, and Young’s inequality, the following
bound

$\displaystyle\|\sigma_{1}\|_{0}\|\boldsymbol{e}_{h,t}\|_{0}+c_{\rm
inv}\|\sigma_{2}\|_{0}h^{-1}\|\boldsymbol{e}_{h,t}\|_{0}$ $\displaystyle\leq$
$\displaystyle 2\|\sigma_{1}\|_{0}^{2}+c_{\rm
inv}^{2}\|\sigma_{2}\|_{0}^{2}h^{-2}+\frac{1}{4}\|\boldsymbol{e}_{h,t}\|_{0}^{2}$
$\displaystyle\leq$ $\displaystyle
Ch^{2(l-1)}+\frac{1}{4}\|\boldsymbol{e}_{h,t}\|_{0}^{2}.$

Absorbing terms in the left-hand side of (73) and collecting the other terms
on the right-hand side leads to

$\frac{1}{2}\left(\|\boldsymbol{e}_{h,t}\|_{0}^{2}+\frac{d}{dt}\nu\|\nabla\boldsymbol{e}_{h}\|_{0}^{2}+\frac{d}{dt}\mu\|\nabla\cdot\boldsymbol{e}_{h}\|_{0}^{2}\right)\leq
Ch^{2(l-1)},$ (74)

such that, assuming for simplicity
$\boldsymbol{u}_{h}(0)=\boldsymbol{s}_{h}(0)$, we derive the estimate

$\int_{0}^{t}\|\boldsymbol{e}_{h,t}(s)\|_{0}^{2}\
ds+\nu\|\nabla\boldsymbol{e}_{h}(t)\|_{0}^{2}+\mu\|\nabla\cdot\boldsymbol{e}_{h}(t)\|_{0}^{2}\leq
Ch^{2(l-1)}.$ (75)

Considering $t\to 0$ in (72), the viscous term and the grad-div stabilization
term vanish, since $\boldsymbol{e}_{h}(0)=\boldsymbol{0}$. A pressure at $t=0$
can be defined as proposed in [13]. Then, one can perform the same analysis as
above and obtains, instead of (74),

$\frac{1}{2}\|\boldsymbol{e}_{h,t}(0)\|_{0}^{2}\leq Ch^{2(l-1)}.$

Applying the inverse inequality yields

$\|(\nabla\boldsymbol{e}_{h,t})(0)\|_{0}^{2}\leq
Ch^{2(l-2)},\quad\|(\nabla\cdot\boldsymbol{e}_{h,t})(0)\|_{0}^{2}\leq
Ch^{2(l-2)},$ (76)

such that the norms on the left-hand sides are bounded for $l\geq 2$.

Taking the time derivative of (72) gives

$\displaystyle\left(\boldsymbol{e}_{h,tt},\boldsymbol{\varphi}_{h}\right)+\nu(\nabla\boldsymbol{e}_{h,t},\nabla\boldsymbol{\varphi}_{h})+b(\boldsymbol{u}_{h,t},\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})-b(\boldsymbol{u}_{t},\boldsymbol{u},\boldsymbol{\varphi}_{h})$
$\displaystyle+b(\boldsymbol{u}_{h},\boldsymbol{u}_{h,t},\boldsymbol{\varphi}_{h})-b(\boldsymbol{u},\boldsymbol{u}_{t},\boldsymbol{\varphi}_{h})+\mu(\nabla\cdot\boldsymbol{e}_{h,t},\nabla\cdot\boldsymbol{\varphi}_{h})$
$\displaystyle=$
$\displaystyle(\sigma_{1,t},\boldsymbol{\varphi}_{h})+(\sigma_{2,t},\nabla\cdot\boldsymbol{\varphi}_{h})\quad\forall\
\boldsymbol{\varphi}_{h}\in{\boldsymbol{V}}_{h}^{l},$

Choosing $\boldsymbol{\varphi}_{h}=\boldsymbol{e}_{h,tt}$ and arguing as
before leads to

$\displaystyle\|\boldsymbol{e}_{h,tt}\|_{0}^{2}+\frac{d}{dt}\frac{1}{2}\nu\|\nabla\boldsymbol{e}_{h,t}\|_{0}^{2}+\frac{d}{dt}\frac{1}{2}\mu\|\nabla\cdot\boldsymbol{e}_{h,t}\|_{0}^{2}$
$\displaystyle\leq$ $\displaystyle\|\sigma_{1,t}\|_{0}^{2}+c_{\rm
inv}^{2}h^{-2}\|\sigma_{2,t}\|_{0}^{2}+\frac{1}{4}\|\boldsymbol{e}_{h,tt}\|_{0}^{2}+|b(\boldsymbol{u}_{h,t}-\boldsymbol{u}_{t},\boldsymbol{u},\boldsymbol{e}_{h,tt})|$
$\displaystyle+|b(\boldsymbol{u}_{h,t},\boldsymbol{u}-\boldsymbol{u}_{h},\boldsymbol{e}_{h,tt})|+|b(\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{u}_{t},\boldsymbol{e}_{h,tt})|+|b(\boldsymbol{u}_{h},\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t},\boldsymbol{e}_{h,tt})|.$

Let us observe that the bounds for $\|\sigma_{1,t}\|_{0}$ and
$\|\sigma_{2,t}\|_{0}$ depend on the regularity of $\boldsymbol{u}_{t}$,
$\boldsymbol{u}_{tt}$, and $p_{t}$, but all are bounded if the solution is
sufficiently regular. For the nonlinear terms, the same arguments as above are
used. For the first term, we obtain

$\displaystyle|b(\boldsymbol{u}_{h,t}-\boldsymbol{u}_{t},\boldsymbol{u},\boldsymbol{e}_{h,tt})|$
(79) $\displaystyle\leq$
$\displaystyle\|\boldsymbol{u}_{h,t}-\boldsymbol{u}_{t}\|_{0}\|\nabla\boldsymbol{u}\|_{\infty}\|\boldsymbol{e}_{h,tt}\|_{0}+\frac{1}{2}\|\nabla\cdot(\boldsymbol{u}_{h,t}-\boldsymbol{u}_{t})\|_{0}\|\boldsymbol{u}\|_{\infty}\|\boldsymbol{e}_{h,tt}\|_{0}$
$\displaystyle\leq$ $\displaystyle
C\|\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t}\|_{1}^{2}+\frac{1}{16}\|\boldsymbol{e}_{h,tt}\|_{0}^{2}.$

It follows from (8), (75), and the inverse inequality (6) that

$\int_{0}^{t}\|\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t}\|_{1}^{2}\leq
Ch^{2(l-2)},$

so that absorbing the second term on the right-hand side of (79) in the left-
hand side of (A) and integrating in time, the corresponding first term on the
right-hand side of (79) is bounded. For the second nonlinear term on the
right-hand side of (A), we obtain

$\displaystyle|b(\boldsymbol{u}_{h,t},\boldsymbol{u}-\boldsymbol{u}_{h},\boldsymbol{e}_{h,tt})|$
$\displaystyle\leq$ $\displaystyle
C\|\boldsymbol{u}_{h,t}\|_{1}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1}c_{\rm
inv}h^{-1}\|\boldsymbol{e}_{h,tt}\|_{0}$ $\displaystyle\leq$ $\displaystyle
Ch^{2(l-2)}\|\boldsymbol{u}_{h,t}\|_{1}^{2}+\frac{1}{16}\|\boldsymbol{e}_{h,tt}\|_{0}^{2}.$

Again, the last term on the right-hand side can be absorbed in the left-hand
side of (A) and the integral with respect to time of the first term is
bounded. For the third nonlinear term we argue as for the second one to get

$|b(\boldsymbol{u}_{h}-\boldsymbol{u},\boldsymbol{u}_{t},\boldsymbol{e}_{h,tt})|\leq\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{1}\|\boldsymbol{u}_{t}\|_{1}\|\boldsymbol{e}_{h,tt}\|_{1}\leq
Ch^{2(l-2)}\|\boldsymbol{u}_{t}\|_{1}^{2}+\frac{1}{16}\|\boldsymbol{e}_{h,tt}\|_{0}^{2}.$

Using for the fourth nonlinear term the inverse inequality (6) and (8) gives

$\displaystyle|b(\boldsymbol{u}_{h},\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t},\boldsymbol{e}_{h,tt})|\leq\|\boldsymbol{u}_{h}\|_{\infty}\|\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t}\|_{1}\|\boldsymbol{e}_{h,tt}\|_{0}$
$\displaystyle+\frac{1}{2}\|\nabla\cdot(\boldsymbol{u}_{h}-\boldsymbol{s}_{h})\|_{\infty}\|\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t}\|_{0}\|\boldsymbol{e}_{h,tt}\|_{0}+\frac{1}{2}\|\nabla\cdot\boldsymbol{s}_{h}\|_{\infty}\|\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t}\|_{0}\|\boldsymbol{e}_{h,tt}\|_{0}$
$\displaystyle\leq$ $\displaystyle
C\left(\|\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t}\|_{1}^{2}+c_{\rm
inv}^{2}h^{-d}h^{2(l-1)}\|\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t}\|_{0}^{2}\right)+\frac{1}{16}\|\boldsymbol{e}_{h,tt}\|_{0}^{2}.$

Again, the last term on the right-hand side above is absorbed into the left-
hand side of (A), the integral of the first term is bounded, and the integral
for the second is

$c_{\rm
inv}^{2}h^{-d}h^{2(l-1)}\int_{0}^{t}\|\boldsymbol{u}_{t}-\boldsymbol{u}_{h,t}(s)\|_{0}^{2}\
ds\leq Ch^{4l-4-d},$

which is bounded since $l\geq 2$. Collecting all error bounds and taking (76)
into account, we conclude that

$\int_{0}^{t}\|\boldsymbol{e}_{h,tt}(s)\|_{0}^{2}\ ds\leq Ch^{2(l-2)},$ (80)

where the constant does not depend explicitly on inverse powers of $\nu$.

Estimate (80) can be applied, in combination with (71), in (66). To bound
(59), which is the term appearing in (4.1), the inverse inequality gives

$\int_{0}^{t}\|\nabla(\boldsymbol{e}_{h,tt})(s)\|_{0}^{2}\ ds\leq
Ch^{2(l-3)},$

such that a robust estimate is proved only for pairs of finite element spaces
with $l\geq 3$.

## References

  * [1] R. A. Adams.  Sobolev spaces.  Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, 1975.  Pure and Applied Mathematics, Vol. 65.
  * [2] F. Brezzi and R. S. Falk.  Stability of higher-order Hood-Taylor methods.  SIAM J. Numer. Anal., 28(3):581–590, 1991.
  * [3] A. Caiazzo, T. Iliescu, V. John, and S. Schyschlowa.  A numerical investigation of velocity-pressure reduced order models for incompressible flows.  J. Comput. Phys., 259:598–616, 2014.
  * [4] D. Chapelle, A. Gariah, and J. Sainte-Marie.  Galerkin approximation with proper orthogonal decomposition: new error estimates and illustrative examples.  ESAIM Math. Model. Numer. Anal., 46(4):731–757, 2012.
  * [5] H. Chen.  Pointwise error estimates for finite element solutions of the Stokes problem.  SIAM J. Numer. Anal., 44(1):1–28, 2006.
  * [6] P. G. Ciarlet.  The finite element method for elliptic problems, volume 40 of Classics in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002.  Reprint of the 1978 original [North-Holland, Amsterdam].
  * [7] J. de Frutos, B. García-Archilla, V. John, and J. Novo.  Grad-div stabilization for the evolutionary Oseen problem with inf-sup stable finite elements.  J. Sci. Comput., 66(3):991–1024, 2016.
  * [8] J. de Frutos, B. García-Archilla, V. John, and J. Novo.  Analysis of the grad-div stabilization for the time-dependent Navier-Stokes equations with inf-sup stable finite elements.  Adv. Comput. Math., 44(1):195–225, 2018.
  * [9] B. García-Archilla, V. John, and J. Novo.  On the convergence order of the finite element error in the kinetic energy for high Reynolds number incompressible flows.  Comput. Methods Appl. Mech. Engrg., 385:Paper No. 114032, 54, 2021\.
  * [10] B. García-Archilla and J. Novo.  Robust error bounds for the Navier-Stokes equations using implicit-explicit second order BDF method with variable steps, 2021.
  * [11] B. García-Archilla, J. Novo, and S. Rubino.  Error analysis of proper orthogonal decomposition data assimilation schemes with grad-div stabilization for the Navier-Stokes equations.  J. Comput. Appl. Math., 411:Paper No. 114246, 30, 2022.
  * [12] B. García-Archilla, J. Novo, and E. S. Titi.  Uniform in time error estimates for a finite element method applied to a downscaling data assimilation algorithm for the Navier-Stokes equations.  SIAM J. Numer. Anal., 58(1):410–429, 2020.
  * [13] J. G. Heywood and R. Rannacher.  Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization.  SIAM J. Numer. Anal., 19(2):275–311, 1982.
  * [14] J. G. Heywood and R. Rannacher.  Finite-element approximation of the nonstationary Navier-Stokes problem. IV. Error analysis for second-order time discretization.  SIAM J. Numer. Anal., 27(2):353–384, 1990.
  * [15] T. Iliescu and Z. Wang.  Are the snapshot difference quotients needed in the proper orthogonal decomposition?  SIAM J. Sci. Comput., 36(3):A1221–A1250, 2014.
  * [16] V. John.  Finite element methods for incompressible flow problems, volume 51 of Springer Series in Computational Mathematics.  Springer, Cham, 2016.
  * [17] V. John, B. Moreau, and J. Novo.  Error analysis of a supg-stabilized pod-rom method for convection-diffusion-reaction equations.  submitted.
  * [18] K. Kean and M. Schneier.  Error analysis of supremizer pressure recovery for POD based reduced-order models of the time-dependent Navier-Stokes equations.  SIAM J. Numer. Anal., 58(4):2235–2264, 2020.
  * [19] B. Koc, S. Rubino, M. Schneier, J. Singler, and T. Iliescu.  On optimal pointwise in time error bounds and difference quotients for the proper orthogonal decomposition.  SIAM J. Numer. Anal., 59(4):2163–2196, 2021.
  * [20] T. Kostova-Vassilevska and G. M. Oxberry.  Model reduction of dynamical systems by proper orthogonal decomposition: error bounds and comparison of methods using snapshots from the solution and the time derivatives.  J. Comput. Appl. Math., 330:553–573, 2018.
  * [21] K. Kunisch and S. Volkwein.  Galerkin proper orthogonal decomposition methods for parabolic problems.  Numer. Math., 90(1):117–148, 2001.
  * [22] S. K. Locke and J. R. Singler.  A new approach to proper orthogonal decomposition with difference quotiens.  arXiv:2106.10224v1 [math.NA] 18 Jun 2021.
  * [23] J. Novo and S. Rubino.  Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows.  SIAM J. Numer. Anal., 59(1):334–369, 2021.
  * [24] J. Sánchez, M. Net, B. García-Archilla, and C. Simó.  Newton-Krylov continuation of periodic orbits for Navier-Stokes flows.  J. Comput. Phys., 201(1):13–33, 2004.
  * [25] M. Schäfer, S. Turek, F. Durst, E. Krause, and R. Rannacher.  Benchmark Computations of Laminar Flow Around a Cylinder, pages 547–566.  Vieweg+Teubner Verlag, Wiesbaden, 1996.
  * [26] P. W. Schroeder and G. Lube.  Pressure-robust analysis of divergence-free and conforming FEM for evolutionary incompressible Navier-Stokes flows.  J. Numer. Math., 25(4):249–276, 2017.
  * [27] J. R. Singler.  New POD error expressions, error bounds, and asymptotic results for reduced order models of parabolic PDEs.  SIAM J. Numer. Anal., 52(2):852–876, 2014.
  * [28] C. Taylor and P. Hood.  A numerical solution of the Navier-Stokes equations using the finite element technique.  Internat. J. Comput. & Fluids, 1(1):73–100, 1973.
  * [29] S. Zhang.  A new family of stable mixed finite elements for the 3D Stokes equations.  Math. Comp., 74(250):543–554, 2005.