# A Hybrid High-Order method for the steady incompressible Navier-Stokes problem

Daniele Di Pietro, Stella Krell

► To cite this version:

Daniele Di Pietro, Stella Krell. A Hybrid High-Order method for the steady incompressible Navier-Stokes problem. Journal of Scientific Computing, 2018, 74 (3), pp.1677-1705. 10.1007/s10915-017-0512-x. hal-01349519v2

HAL Id: hal-01349519

https://hal.science/hal-01349519v2

Submitted on 13 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Hybrid High-Order method for the steady incompressible Navier-Stokes problem*

Daniele A. Di Pietro†<fnref n="1" /> and Stella Krell‡<fnref n="2" />

January 13, 2017

## Abstract

In this work we introduce and analyze a novel Hybrid High-Order method for the steady incompressible Navier-Stokes equations. The proposed method is inf-sup stable on general polyhedral meshes, supports arbitrary approximation orders, and is (relatively) inexpensive thanks to the possibility of statically condensing a subset of the unknowns at each nonlinear iteration. We show under general assumptions the existence of a discrete solution, which is also unique provided a data smallness condition is verified. Using a compactness argument, we prove convergence of the sequence of discrete solutions to minimal regularity exact solutions for general data. For more regular solutions, we prove optimal convergence rates for the energy-norm of the velocity and the L<fnref n="2" />-norm of the pressure under a standard data smallness assumption. More precisely, when polynomials of degree $k \ge 0$ at mesh elements and faces are used, both quantities are proved to converge as $h^{k+1}$ (with h denoting the meshsize).

2010 Mathematics Subject Classification: 65N08, 65N30, 65N12, 35Q30, 76D05

Keywords: Hybrid High-Order, incompressible Navier-Stokes, polyhedral meshes, compactness, error estimates

## 1 Introduction

In this work we introduce and analyze a novel Hybrid High-Order (HHO) method for the steady incompressible Navier-Stokes equations. The proposed method is inf-sup stable on general meshes including polyhedral elements and nonmatching interfaces, it supports arbitrary approximation order, and has a reduced computational cost thanks to the possibility of statically condensing a subset of both velocity and pressure degrees of freedom (DOFs) at each nonlinear iteration. A complete analysis covering general data is provided.

Let $\Omega \subset \mathbb{R}^d$, $d \in \{2, 3\}$, denote a bounded connected open polyhedral domain. We additionally assume that $\Omega$ does not have cracks, i.e., it lies on one side of its boundary. The incompressible Navier-Stokes problem consists in finding the velocity field $u: \Omega \to \mathbb{R}^d$ and the pressure field $p: \Omega \to \mathbb{R}$ such that

$$-\nu \Delta u + \nabla u \cdot \mathbf{u} + \nabla p = \mathbf{f} \quad \text{in } \Omega, \quad (1a)$$

$$\operatorname{div} u = 0 \quad \text{in } \Omega, \quad (1b)$$

$$\mathbf{u} = \mathbf{0} \quad \text{on } \partial\Omega, \quad (1c)$$

$$\int_{\Omega} p = 0, \quad (1d)$$

<footnote n="1">¹Université de Montpellier, Institut Montpelliérain Alexander Grothendieck, 34095 Montpellier, France</footnote>

<footnote n="2">²Université de Nice, Laboratoire J. A. Dieudonné, avenue Valrose, 06000 Nice, France</footnote>

<footnote>*The work of D. A. Di Pietro was partially supported by Agence Nationale de la Recherche project HHOMM (ANR-15-CE40-0005).</footnote>

<footnote>†daniele.di-pietro@umontpellier.fr</footnote>

<footnote>‡stella.krell@unice.fr</footnote>
where $\nu > 0$ denotes the (constant) kinematic viscosity and $f \in L^2(\Omega)^d$ a volumetric body force. For the convective term, here and in what follows we will use the matrix-vector product notation with vector quantities treated as column vectors, so that, e.g., $\nabla u u = (\sum_{j=1}^d (\partial_j u_i) u_j)_{1 \le i \le d}$ with $u_i$ denoting the $i$th component of $u$ and $\partial_j$ the derivative along the $j$th direction. Let $U := H_0^1(\Omega)^d$ and $P := L_0^2(\Omega)$ with $L_0^2(\Omega) := \{q \in L^2(\Omega) | \int_{\Omega} q = 0\}$. A classical weak formulation of problem (1) reads: Find $(u, p) \in U \times P$ such that

$$ \nu a(u, v) + t(u, u, v) + b(v, p) = \int_{\Omega} f \cdot v \quad \forall v \in U, \quad (2a) $$

$$ -b(u, q) = 0 \quad \forall q \in P, \quad (2b) $$

with bilinear forms $a: U \times U \to \mathbb{R}$ and $b: U \times P \to \mathbb{R}$ and trilinear form $t: U \times U \times U \to \mathbb{R}$ such that

$$ a(u, v) := \int_{\Omega} \nabla u : \nabla v, \quad b(v, q) := - \int_{\Omega} (\operatorname{div} v) q, \quad t(w, u, v) := \frac{1}{2} \int_{\Omega} v^T \nabla u w - \frac{1}{2} \int_{\Omega} u^T \nabla v w, \quad (3) $$

where $^T$ denotes the transpose operator. Notice that all the discussion below can be extended to more general boundary conditions, but we stick to the homogeneous Dirichlet case for simplicity of exposition. Also, the modifications to handle variable kinematic viscosities are briefly discussed in Remark 5.

The literature on the numerical approximation of problem (2) is vast, and giving a detailed account lies out of the scope of the present work. We therefore mention here only those numerical methods which share similar features with our approach. The possibility to increase the approximation order and, possibly, to use general meshes and analogues to the ones considered here are supported by discontinuous Galerkin (DG) methods. Their application to the incompressible Navier-Stokes problem has been considered in several works starting from the early 00's; a non exhaustive list of references includes [3-5,15,19,33,35,42,43]; cf. also [20, Chapter 6] for a pedagogical introduction. The Hybridizable discontinuous Galerkin (HDG) method of [10, 14] has also been applied to the discretization of the incompressible Navier-Stokes equations in several recent works [11, 31, 39, 41, 44]. Albeit this is not explicitly pointed out in all of the above references, also HDG methods often support general meshes as well as the possibility to increase the approximation order. The relation between HDG and HHO methods (originally introduced in [21] in the context of quasi-incompressible linear elasticity) has been recently explored in [13] for a pure diffusion problem. High-order discretizations on general meshes are also provided by the Virtual Element Method; see, e.g., [6], where a divergence-free virtual space for the Stokes problem is proposed.

The HHO scheme studied in this work hinges on DOFs located at mesh elements and faces that are discontinuous polynomials of degree $k \ge 0$. Based on these DOFs, by solving local problems inside each element we obtain reconstructions of the velocity and of its gradient that are used to formulate the diffusive and convective terms in the momentum balance equation, as well as a reconstruction of the divergence used in the velocity-pressure coupling terms. More precisely, the discretization of the viscous term stems from a variation of the diffusive bilinear form originally introduced in [22]; for the convective term, we use a skew-symmetric formulation designed so as not to contribute to the kinetic energy balance; the velocity-pressure coupling is, on the other hand, analogous to that of [2, 23].

The resulting method has several advantageous features: (i) it supports general meshes possibly including polyhedral elements and nonmatching interfaces (resulting, e.g., from nonconforming mesh refinement); (ii) it allows one to increase the spatial approximation order to accelerate convergence in the presence of (locally) regular solutions; (iii) it is (relatively) inexpensive thanks to the possibility of statically condensing all element-based velocity DOFs and all but one pressure DOF per element at each non-linear iteration (cf. Remark 10 below for further details). Additionally, thanks to the underlying fully discontinuous polynomial spaces, the proposed method can be expected to accommodate abrupt spatial variations of the unknowns. Existence of a discrete solution is proved in Theorem 12 below for general data resorting to classical arguments in nonlinear analysis [16]. Uniqueness, on the other hand, is shown in Theorem 13 below under a standard smallness assumption on the volumetric body force.

A complete convergence analysis of the method is carried out. First, using a compactness argument inspired by the recent literature on finite volume methods (cf., e.g., [12, 28-30]), we show in Theorem 14 below that the sequence of discrete solutions on a refined mesh family converges (up to a subsequence)
to the continuous one for general data and without assuming more regularity for the exact solution than required by the weak formulation. Convergence extends to the whole sequence when the continuous solution is unique. The use of compactness techniques in the context of high-order methods is quite original; see, e.g., [19, 20, 38] concerning the incompressible Navier-Stokes problem. Key technical results required to prove convergence by compactness are the discrete Sobolev embeddings and Rellich-Kondrachov compactness results recently proved in [17] in the context of nonlinear Leray-Lions problems.

Then, in Theorem 15 below, we prove error estimates for regular exact solutions under a suitable data smallness assumption. When polynomials of degree $k \ge 0$ are used, we show that both the energy-norm of the velocity and the $L^2$-norm of the pressure converge as $h^{k+1}$ ($h$ denotes here the meshsize). These convergence rates are similar to the ones recently derived in [41] for a HDG method with pressure and velocity spaces chosen as in [26, 45]. A major difference with respect to [41] is that we obtain them here using polynomials $k$ instead of $(k+1)$ inside mesh elements (this is precisely one of the major outcomes of the HHO technology identified in [13]). Another difference with respect to [11, 41] is that our trilinear form is expressed in terms of a discrete gradient reconstruction and designed so that it does not contribute to the kinetic energy balance, a feature which simplifies several arguments in the analysis; cf. Remark 8 for further details. We also show numerically that the $L^2$-norm of the error on the velocity converges as $h^{k+2}$.

The rest of this paper is organized as follows. In Section 2 we introduce mesh-related notations and recall a few basic results on broken functional spaces. In Section 3 we define the local reconstructions, formulate the discretizations of the various terms appearing in (2), and state the discrete problem. In Section 4 we discuss the existence and uniqueness of a discrete solution, prove convergence to minimal regularity exact solutions for general data, and estimate the convergence rate for smooth exact solutions and small data. The theoretical results are illustrated on a numerical example including a comparison with a HDG-inspired trilinear form. In Section 5 we give proof of the properties of the viscous and velocity-pressure coupling bilinear forms and of the convective trilinear form used in the analysis.

## 2 Mesh and basic results

Let $\mathcal{H} \subset \mathbb{R}_*^+$ denote a countable set of meshsizes having 0 as its unique accumulation point. We consider refined mesh sequences $(\mathcal{T}_h)_{h \in \mathcal{H}}$ where, for all $h \in \mathcal{H}$, $\mathcal{T}_h = \{T\}$ is a finite collection of nonempty disjoint open polyhedral elements such that $\bar{\Omega} = \cup_{T \in \mathcal{T}_h} \bar{T}$ and $h = \max_{T \in \mathcal{T}_h} h_T$ ($h_T$ stands for the diameter of $T$). For all $h \in \mathcal{H}$, the set of mesh faces $\mathcal{F}_h$ is a finite collection of disjoint subsets $F$ of $\bar{\Omega}$ such that, for any $F \in \mathcal{F}_h$, $F$ is an open subset of a hyperplane of $\mathbb{R}^d$ such that $|F|_{d-1} > 0$ and $|\bar{F} \setminus F|_{d-1} = 0$ with $|\cdot|_{d-1}$ denoting the $(d-1)$-dimensional Hausdorff measure. Additionally, it holds that $\cup_{T \in \mathcal{T}_h} \partial T = \cup_{F \in \mathcal{F}_h} \bar{F}$ and, for all $F \in \mathcal{F}_h$, (i) either there exist $T_1, T_2 \in \mathcal{T}_h$ such that $F = \partial T_1 \cap \partial T_2$ (and $\bar{F}$ is an interface) or (ii) there exists $T \in \mathcal{T}_h$ such that $F = \partial T \cap \partial \Omega$ (and $\bar{F}$ is a boundary face). The set of interfaces is denoted by $\mathcal{F}_h^i$, the set of boundary faces by $\mathcal{F}_h^b$, and it holds that $\mathcal{F}_h = \mathcal{F}_h^i \cup \mathcal{F}_h^b$. For all $T \in \mathcal{T}_h$, the set $\mathcal{F}_T := \{F \in \mathcal{F}_h \mid F \subset \partial T\}$ collects the faces lying on the boundary of $T$ and, for all $F \in \mathcal{F}_T$, we denote by $n_{TF}$ the normal to $F$ pointing out of $T$. A normal vector $n_F$ is associated to each internal face by fixing once and for all an (arbitrary) orientation, whereas for boundary faces $n_F$ points out of $\Omega$.

We assume that $(\mathcal{T}_h)_{h \in \mathcal{H}}$ is admissible in the sense of [20, Chapter 1], i.e., for all $h \in \mathcal{H}$, $\mathcal{T}_h$ admits a matching simplicial submesh $\mathcal{T}_h$ and there exists a real number $\rho > 0$ (the mesh regularity parameter) independent of $h$ such that the following conditions hold: (i) for all $h \in \mathcal{H}$ and every simplex $S \in \mathcal{T}_h$ of diameter $h_S$ and inradius $r_S$, $\rho h_S \le r_S$; (ii) for all $h \in \mathcal{H}$, all $T \in \mathcal{T}_h$, and all $S \in \mathcal{T}_h$ such that $S \subset T$, $\rho h_T \le h_S$. We refer to [20, Chapter 1] and [17, 18] for a set of geometric and functional analytic results valid on admissible meshes. We recall, in particular, that, under these regularity assumptions, the number of faces of each element is uniformly bounded.

Let $X$ be a subset of $\mathbb{R}^d$ and, for an integer $l \ge 0$, denote by $\mathbb{P}^l(X)$ the space spanned by the restrictions to $X$ of polynomials in the space variables of total degree $l$. In what follows, the set $X$ will represent a mesh element or face. We denote by $\pi_X^l: \mathbb{P}^l(X) \to \mathbb{P}^l(X)$ the $L^2$-orthogonal projector such that, for all
$v \in L^1(X),$

$$ \int_X (v - \pi_X^l v)w = 0 \quad \forall w \in \mathbb{P}^l(X). \qquad (4) $$

The vector- and matrix-valued $L^2$-orthogonal projectors, both denoted by $\pi_X^l$, are obtained applying $\pi_X^l$ component-wise. The following optimal $W^{s,p}$-approximation properties are proved in [17, Appendix A.2] using the classical theory of [25] (cf. also [9, Chapter 4]): There is $C > 0$ such that, for all $l \ge 0$, all $h \in \mathcal{H}$, all $T \in \mathcal{T}_h$, all $s \in \{1, \dots, l+1\}$, all $p \in [1, +\infty]$, all $v \in W^{s,p}(T)$, and all $m \in \{0, \dots, s-1\}$,

$$ |v - \pi_T^l v|_{W^{m,p}(T)} + h_T^{\frac{1}{p}} |v - \pi_T^l v|_{W^{m,p}(\mathcal{T}_T)} \le C h_T^{s-m} |v|_{W^{s,p}(T)}, \qquad (5) $$

where $W^{m,p}(\mathcal{T}_T)$ is spanned by functions that are in $W^{m,p}(F)$ for all $F \in \mathcal{T}_T$. At the global level, the space of broken polynomial functions on $\mathcal{T}_h$ of degree $l$ is denoted by $\mathbb{P}^l(\mathcal{T}_h)$, and $\pi_h^l$ is the corresponding $L^2$-orthogonal projector. The broken gradient operator on $\mathcal{T}_h$ is denoted by $\nabla_h$.

Let $p \in [1, +\infty]$. We recall the following continuous trace inequality: There is $C > 0$ such that, for all $h \in \mathcal{H}$ and all $T \in \mathcal{T}_h$ it holds for all $v \in W^{1,p}(T)$,

$$ h_T^{\frac{1}{p}} \|v\|_{L^p(\partial T)} \le C (\|v\|_{L^p(T)} + h_T \|\nabla v\|_{L^p(T)^d}). \qquad (6) $$

Let an integer $l \ge 0$ be fixed. Using (6) followed by the discrete inverse inequality

$$ \|\nabla v\|_{L^p(T)^d} \le C h_T^{-1} \|v\|_{L^p(T)}, \qquad (7) $$

valid for all $T \in \mathcal{T}_h$ and $v \in \mathbb{P}^l(T)$ with $C > 0$ independent of $h$ and of $T$, we obtain the following discrete trace inequality: There is $C > 0$ such that, for all $h \in \mathcal{H}$ and all $T \in \mathcal{T}_h$ it holds for all $v \in \mathbb{P}^l(T)$,

$$ h_T^{\frac{1}{p}} \|v\|_{L^p(\partial T)} \le C \|v\|_{L^p(T)}. \qquad (8) $$

Throughout the paper, we often write $a \lesssim b$ (resp. $a \gtrsim b$) to mean $a \le Cb$ (resp. $a \ge Cb$) with real number $C > 0$ independent of the meshsize $h$ and of the kinematic viscosity $\nu$. Constants are named when needed in the discussion.

# 3 Discretization

In this section we define the discrete counterparts of the various terms appearing in (2), state their properties, and formulate the discrete problem.

## 3.1 Discrete spaces

Let a polynomial degree $k \ge 0$ be fixed. We define the following hybrid space containing element-based and face-based velocity DOFs:

$$ \underline{U}_h^k := \left( \times_{T \in \mathcal{T}_h} \mathbb{P}^k(T)^d \right) \times \left( \times_{F \in \mathcal{F}_h} \mathbb{P}^k(F)^d \right). \qquad (9) $$

For the elements of $\underline{U}_h^k$ we use the underlined notation $\underline{v}_h = ((v_T)_{T \in \mathcal{T}_h}, (v_F)_{F \in \mathcal{F}_h})$. We define the global interpolator $\underline{I}_h^k : H^1(\Omega)^d \to \underline{U}_h^k$ such that, for all $v \in H^1(\Omega)^d$,

$$ \underline{I}_h^k v := ((\underline{\pi}_T^k v)_{T \in \mathcal{T}_h}, (\underline{\pi}_F^k v)_{F \in \mathcal{F}_h}). $$

For every mesh element $T \in \mathcal{T}_h$, we denote by $\underline{U}_h^k$ and $\underline{I}_h^k$ the restrictions to $T$ of $\underline{U}_h^k$ and $\underline{I}_h^k$, respectively. Similarly, $\underline{v}_T = (v_T, (v_F)_{F \in \mathcal{F}_T})$ denotes the restriction to $T$ of a generic vector $\underline{v}_h \in \underline{U}_h^k$. Also, for an
element $\underline{v}_h \in \underline{U}_h^k$, we denote by $v_h$ (no underline) the broken polynomial function in $\mathbb{P}^k(\mathcal{T}_h)^d$ such that $v_{h|T} = v_T$ for all $T \in \mathcal{T}_h$. Finally, we define on $\underline{U}_h^k$ the following seminorm:

$$ \|\underline{v}_h\|_{1,h}^2 := \sum_{T \in \mathcal{T}_h} \|\underline{v}_T\|_{1,T}^2, \qquad (10) $$

where, for all $T \in \mathcal{T}_h$,

$$ \|\underline{v}_T\|_{1,T}^2 := \|\nabla v_T\|_T^2 + |\underline{v}_T|_{1,\partial T}^2, \quad |\underline{v}_T|_{1,\partial T}^2 := \sum_{F \in \mathcal{F}_T} h_F^{-1} \|v_F - v_T\|_F^2. \qquad (11) $$

The following boundedness property holds for the global interpolator $\underline{I}_h^k$: For all $v \in H^1(\Omega)^d$,

$$ \|\underline{I}_h^k v\|_{1,h} \le C_I \|v\|_{H^1(\Omega)^d}, \qquad (12) $$

with real number $C_I > 0$ independent of $h$. Its proof relies on the stability properties of the $L^2$-projectors on elements and faces proved in [17].

The following velocity and pressure spaces embed the homogeneous boundary conditions for the velocity and the zero-average constraint for the pressure, respectively:

$$ \underline{U}_{h,0}^k := \{ \underline{v}_h \in \underline{U}_h^k \mid v_F = 0 \quad \forall F \in \mathcal{F}_h^b \}, \quad P_h^k := \{ q_h \in \mathbb{P}^k(\mathcal{T}_h) \mid \int_{\Omega} q_h = 0 \}. \qquad (13) $$

It is a simple matter to check that the map $\|\cdot\|_{1,h}$ defines a norm on $\underline{U}_{h,0}^k$. We also note the following discrete Sobolev embeddings, a consequence of [17, Proposition 5.4]: For all $p \in [1, +\infty)$ if $d = 2, p \in [1, 6]$ if $d = 3$, it holds for all $\underline{v}_h \in \underline{U}_{h,0}^k$,

$$ \|\underline{v}_h\|_{L^p(\Omega)^d} \le C_s \|\underline{v}_h\|_{1,h}, \qquad (14) $$

with real number $C_s > 0$ independent of $h$.

## 3.2 Reconstructions of differential operators

Let an element $T \in \mathcal{T}_h$ be fixed. For any polynomial degree $l \ge 0$, we define the local gradient reconstruction operator $G_T^l: \underline{U}_T^k \to \mathbb{P}^l(T)^{d \times d}$ such that, for all $\underline{v}_T \in \underline{U}_T^k$ and all $\tau \in \mathbb{P}^l(T)^{d \times d}$,

$$ \int_T G_T^l \underline{v}_T : \tau = - \int_T v_T \cdot (\mathrm{div} \, \tau) + \sum_{F \in \mathcal{F}_T} \int_F v_F \cdot (\tau \, n_{TF}) \qquad (15a) $$

$$ = \int_T \nabla v_T : \tau + \sum_{F \in \mathcal{F}_T} \int_F (v_F - v_T) \cdot (\tau \, n_{TF}), \qquad (15b) $$

where we have used integration by parts to pass to the second line. In (15a), the right-hand is designed to resemble an integration by parts formula where the role of the function in volumetric and boundary integrals is played by element-based and face-based DOFs, respectively.

For the discretization of the viscous term, we will need the local velocity reconstruction operator $\underline{r}_T^{k+1}: \underline{U}_T^k \to \mathbb{P}^{k+1}(T)^d$ obtained in a similar spirit as $G_T^k$ and such that, for all $\underline{v}_T \in \underline{U}_T^k$,

$$ \int_T \nabla \underline{r}_T^{k+1} \underline{v}_T : \nabla w = - \int_T v_T \cdot \Delta w + \sum_{F \in \mathcal{F}_T} \int_F v_F \cdot (\nabla w \cdot n_{TF}) \quad \forall w \in \mathbb{P}^{k+1}(T)^d, \qquad (16) $$

with closure condition $\int_T (\underline{r}_T^{k+1} \underline{v}_T - v_T) = 0$.

Finally, we define the discrete divergence operator $D_T^k: \underline{U}_T^k \to \mathbb{P}^k(T)$ such that, for all $\underline{v}_T \in \underline{U}_T^k$ and all $q \in \mathbb{P}^k(T)$,

$$ \int_T D_T^k \underline{v}_T q = - \int_T v_T \cdot \nabla q + \sum_{F \in \mathcal{F}_T} \int_F (v_F \cdot n_{TF}) q \qquad (17a) $$

$$ = \int_T (\mathrm{div} \, v_T) q + \sum_{F \in \mathcal{F}_T} \int_F (v_F - v_T) \cdot n_{TF} q, \qquad (17b) $$
where we have used integration by parts to pass to the second line. By definition, we have

$$D_T^k = \mathrm{tr}(G_T^k). \qquad (18)$$

We also define global versions of the gradient, velocity reconstruction, and divergence operators letting $G_h^l : \underline{U}_h^k \to \mathbb{P}^l(\mathcal{T}_h)^{d \times d}$, $r_h^{k+1} : \underline{U}_h^k \to \mathbb{P}^{k+1}(\mathcal{T}_h)^d$, and $D_h^k : \underline{U}_h^k \to \mathbb{P}^k(\mathcal{T}_h)$ be such that, for all $\underline{v}_h \in \underline{U}_h^k$ and all $T \in \mathcal{T}_h$,

$$(G_h^l \underline{v}_h)|_T := G_T^l \underline{v}_T, \quad (r_h^{k+1} \underline{v}_h)|_T := r_T^{k+1} \underline{v}_T, \quad (D_h^k \underline{v}_h)|_T := D_T^k \underline{v}_T.$$

**Proposition 1** (Properties of $G_h^l$). The global discrete gradient operator $G_h^l$ satisfies the following properties:

1) Boundedness. For all $l \ge 0$ and all $\underline{v}_h \in \underline{U}_h^k$, it holds

$$\|G_h^l \underline{v}_h\|_{L^2(\Omega)^{d \times d}} \lesssim \|\underline{v}_h\|_{1,h}. \qquad (19)$$

2) Consistency. For all $l \ge 0$ and all $\underline{v} \in H^m(\Omega)^d$ with $m = l + 2$ if $l \le k$, $m = k + 1$ otherwise,

$$\|G_h^l \underline{I}_h^k \underline{v} - \nabla \underline{v}\|_{L^2(\Omega)^{d \times d}} + \left( \sum_{T \in \mathcal{T}_h} h_T \|G_T^l \underline{I}_T^k \underline{v} - \nabla \underline{v}\|_{L^2(\partial T)^{d \times d}}^2 \right)^{\frac{1}{2}} \lesssim h^{m-1} \|\underline{v}\|_{H^m(\Omega)^d}. \qquad (20)$$

As a consequence, for all $\underline{\Phi} \in C_c^\infty(\Omega)^d$ and all $l, k$ such that $m > 1$ (i.e., provided $l = 0$ if $k = 0$), $G_h^l \underline{I}_h^k \underline{\Phi} \to \nabla \underline{\Phi}$ strongly in $L^2(\Omega)^{d \times d}$.

3) Sequential consistency. Let $(\underline{v}_h)_{h \in \mathcal{H}}$ denote a sequence in $(\underline{U}_{h,0}^k)_{h \in \mathcal{H}}$ bounded in the $\|\cdot\|_{1,h}$-norm. Then, there is $\underline{v} \in \underline{U}$ such that, as $h \to 0$, up to a subsequence,

* $\underline{v}_h \to \underline{v}$ strongly in $L^p(\Omega)^d$ for all $p \in [1, +\infty)$ if $d = 2$, $p \in [1, 6)$ if $d = 3$;

* $G_h^l \underline{v}_h \to \nabla \underline{v}$ weakly in $L^2(\Omega)^{d \times d}$ for all $l \ge 0$.

*Proof.* 1) Boundedness. Let an element $T \in \mathcal{T}_h$ be fixed, make $\underline{\tau} = G_T^l \underline{v}_T$ in (15b) and use the Cauchy-Schwarz inequality followed by the discrete trace inequality (8) with $p = 2$ to infer, for all $\underline{v}_T \in \underline{U}_T^k$,

$$\|G_T^l \underline{v}_T\|_{L^2(T)^{d \times d}} \lesssim \|\underline{v}_T\|_{1,T}.$$

Squaring the above inequality and summing over $T \in \mathcal{T}_h$, (19) follows.

2) Consistency. Let $\underline{v} \in H^m(\Omega)^d$. For all $T \in \mathcal{T}_h$, plugging the definition of $\underline{I}_T^k \underline{v}$ into (15a), we get, for all $\underline{\tau} \in \mathbb{P}^l(T)^{d \times d}$,

$$\int_T (G_T^l \underline{I}_T^k \underline{v} - \nabla \underline{v}) : \underline{\tau} = - \int_T (\underline{\pi}_T^k \underline{v} - \underline{v}) \cdot (\mathrm{div} \, \underline{\tau}) + \sum_{F \in \mathcal{F}_T} \int_F (\underline{\pi}_F^k \underline{v} - \underline{v}) \cdot (\underline{\tau} \, \underline{n}_{TF}). \qquad (21)$$

Recalling the definition (4) of $\underline{\pi}_T^k$ and $\underline{\pi}_F^k$, we get from the previous expression that

$$\forall n \le k, \quad \int_T (G_T^l \underline{I}_T^k \underline{v} - \nabla \underline{v}) : \underline{\tau} = 0 \quad \forall \underline{\tau} \in \mathbb{P}^n(T)^{d \times d}, \qquad (22)$$

since $(\mathrm{div} \, \underline{\tau}) \in \mathbb{P}^{n-1}(T)^d \subset \mathbb{P}^k(T)^d$ and $\underline{\tau}|_F \underline{n}_{TF} \in \mathbb{P}^n(F)^d \subset \mathbb{P}^k(F)^d$. If $l \le k$, this shows in particular that for all $T \in \mathcal{T}_h$ it holds

$$G_T^l \underline{I}_T^k \underline{v} = \underline{\pi}_T^l \nabla \underline{v}, \qquad (23)$$

and (20) is an immediate consequence of the approximation properties (5) of the $L^2$-orthogonal projector.

On the other hand, if $l > k$, making $\underline{\tau} = G_T^l \underline{I}_T^k \underline{v} - \underline{\pi}_h^l \nabla \underline{v}$ in (21) and using the Cauchy-Schwarz, discrete inverse (7) and trace (8) inequalities (both with $p = 2$) to bound the right-hand side, we infer

$$\|G_T^l \underline{I}_T^k \underline{v} - \underline{\pi}_h^l \nabla \underline{v}\|_{L^2(T)^{d \times d}} \lesssim h^k \|\underline{v}\|_{H^{k+1}(T)^d}. \qquad (24)$$
Hence, using the triangle inequality, we arrive at

$$ \|G_T^l \underline{I}_T^k \boldsymbol{v} - \nabla \boldsymbol{v}\|_{L^2(T)^{d \times d}} \lesssim h_T^k \|\boldsymbol{v}\|_{H^{k+1}(T)^d}, $$

and the bound for the first term in (20) follows squaring and summing over $T \in \mathcal{T}_h$. To prove the bound for the second term we observe that it holds, for all $T \in \mathcal{T}_h$,

$$ \begin{aligned} h_T^{\frac{1}{2}} \|G_T^l \underline{I}_T^k \boldsymbol{v} - \nabla \boldsymbol{v}\|_{L^2(\partial T)^{d \times d}} &\lesssim h_T^{\frac{1}{2}} \|G_T^l \underline{I}_T^k \boldsymbol{v} - \pi_T^l \nabla \boldsymbol{v}\|_{L^2(\partial T)^{d \times d}} + h_T^{\frac{1}{2}} \|\pi_T^l \nabla \boldsymbol{v} - \nabla \boldsymbol{v}\|_{L^2(\partial T)^{d \times d}} \\ &\lesssim \|G_T^l \underline{I}_T^k \boldsymbol{v} - \pi_T^l \nabla \boldsymbol{v}\|_{L^2(T)^{d \times d}} + h_T^k \|\boldsymbol{v}\|_{H^{k+1}(T)^d} \end{aligned} $$

where we have inserted $\pm \pi_T^l \nabla \boldsymbol{v}$ inside the norm and used the triangle inequality in the first line, while we have used the discrete trace inequality (8) with $p = 2$, and the approximation properties of $\pi_T^l$ to pass to the second line. Using (24) to estimate the first term in the right-hand side of the above inequality, the bound for the second term in (20) is readily proved squaring and summing over $T \in \mathcal{T}_h$.

3) *Sequential consistency.* The proof for $l = k$ in the scalar case is given in [17, Proposition 5.6]. A close inspection shows that the arguments still stand when $l \neq k$ provided that we replace $\pi_T^k \Phi$ by $\pi_T^0 \Phi$. $\square$

*Remark 2 (Commuting property for $D_h^k$).* Combining (23) with (18), it is a simple matter to infer the following commuting property for $D_h^k$: For all $\boldsymbol{v} \in \mathcal{U}$,

$$ D_h^k \underline{I}_h^k \boldsymbol{v} = \pi_h^k (\operatorname{div} \boldsymbol{v}). \qquad (25) $$

This property is crucial to prove the inf-sup condition of Proposition 9 below using classical arguments in the analysis of saddle-point problems (cf., e.g., the reference textbook [8]).

## 3.3 Viscous term

The viscous term is discretized by means of the bilinear form $a_h$ on $\underline{U}_h^k \times \underline{U}_h^k$ such that, for all $\underline{u}_h, \underline{v}_h \in \underline{U}_h^k$,

$$ a_h(\underline{u}_h, \underline{v}_h) := \int_{\Omega} G_h^k \underline{u}_h : G_h^k \underline{v}_h + s_h(\underline{u}_h, \underline{v}_h), \qquad (26) $$

with stabilization bilinear form $s_h$ defined as follows:

$$ s_h(\underline{u}_h, \underline{v}_h) := \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_T} h_F^{-1} \int_F \delta_{TF}^k \underline{u}_T \cdot \delta_{TF}^k \underline{v}_T, $$

where, for all $T \in \mathcal{T}_h$ and all $F \in \mathcal{F}_T$, we have introduced the face-based residual operator $\delta_{TF}^k : \underline{U}_T^k \to \mathbb{P}^k(F)^d$ such that

$$ \delta_{TF}^k \underline{v}_T := \pi_F^k (\underline{v}_F - \underline{r}_T^{k+1} \underline{v}_T - \pi_T^k (\underline{v}_T - \underline{r}_T^{k+1} \underline{v}_T)). \qquad (27) $$

This specific form of the penalized residual ensures the following consistency property (cf. [22, Remark 6] for further insight): For all $\boldsymbol{v} \in H^{k+2}(\Omega)^d$,

$$ s_h(\underline{I}_h^k \boldsymbol{v}, \underline{I}_h^k \boldsymbol{v})^{\frac{1}{2}} \lesssim h^{k+1} \|\boldsymbol{v}\|_{H^{k+2}(\Omega)}. \qquad (28) $$

The proof of the following result is postponed to Section 5.1.

**Proposition 3 (Properties of $a_h$).** The bilinear form $a_h$ has the following properties:

1) Stability and boundedness. It holds, for all $\underline{v}_h \in \underline{U}_h^k$,

$$ C_a^{-1} \|\underline{v}_h\|_{1,h}^2 \lesssim \|\underline{v}_h\|_{a,h}^2 := a_h(\underline{v}_h, \underline{v}_h) \lesssim C_a \|\underline{v}_h\|_{1,h}^2, \qquad (29) $$

with real number $C_a > 0$ independent of $h$. Consequently, the map $\|\cdot\|_{a,h}$ defines a norm on $\underline{U}_{h,0}^k$ uniformly equivalent to $\|\cdot\|_{1,h}$.
2) Consistency. For all $\mathbf{v} \in \mathcal{U} \cap H^{k+2}(\Omega)^d$, it holds

$$ \sup_{\underline{\mathbf{w}}_h \in \underline{\mathcal{U}}_{h,0}^k, \|\underline{\mathbf{w}}_h\|_{1,h}=1} \left| \int_{\Omega} \Delta \mathbf{v} \cdot \mathbf{w}_h + a_h(\underline{\mathbf{I}}_h^k \mathbf{v}, \underline{\mathbf{w}}_h) \right| \lesssim h^{k+1} \|\mathbf{v}\|_{H^{k+2}(\Omega)^d}. \quad (30) $$

3) Sequential consistency. Let $(\underline{\mathbf{v}}_h)_{h \in \mathcal{H}}$ denote a sequence in $(\underline{\mathcal{U}}_{h,0}^k)_{h \in \mathcal{H}}$ bounded in the $\|\cdot\|_{1,h}$-norm with limit $\mathbf{v} \in \mathcal{U}$ (cf. point 3) in Proposition 1). Then, it holds for all $\mathbf{\Phi} \in C_c^\infty(\Omega)^d$, as $h \to 0$, up to a subsequence,

$$ a_h(\underline{\mathbf{v}}_h, \underline{\mathbf{I}}_h^k \mathbf{\Phi}) \to a(\mathbf{v}, \mathbf{\Phi}). $$

Some remarks are of order.

Remark 4 (Alternative viscous bilinear form). An alternative choice corresponding to the original HHO bilinear form of [22] is

$$ a_h(\underline{\mathbf{u}}_h, \underline{\mathbf{v}}_h) := \int_{\Omega} \nabla_h \mathbf{r}_h^{k+1} \underline{\mathbf{u}}_h : \nabla_h \mathbf{r}_h^{k+1} \underline{\mathbf{v}}_h + s_h(\underline{\mathbf{u}}_h, \underline{\mathbf{v}}_h), $$

where the difference with respect to (26) lies in the fact that $\nabla_h \mathbf{r}_h^{k+1}$ replaces $\mathbf{G}_h^k$ in the consistency term. Properties 1)–2) in Proposition 3 are straightforward consequences of [22, Lemma 4 and Theorem 8], respectively. Property 3), on the other hand, would require proving for $\nabla_h \mathbf{r}_h^{k+1}$ sequential consistency as in point 3) of Proposition 1.

Remark 5 (Variable kinematic viscosity). A more general form for the viscous term in (1a) accommodating variable kinematic viscosities $\nu: \Omega \to \mathbb{R}$ is

$$ -\operatorname{div} \boldsymbol{\sigma}(\mathbf{u}), \quad \boldsymbol{\sigma}(\mathbf{u}) = 2\nu \nabla_s \mathbf{u}, $$

where $\nabla_s$ denotes the symmetric gradient operator. Our discretization can be modified to accommodate this case adapting the ideas developed in [21] in the framework of linear elasticity. Assume, for the sake of simplicity, that $\nu$ is piecewise constant on a partition of $\Omega$, and that for all $h \in \mathcal{H}$ the mesh $\mathcal{T}_h$ is compliant with the partition (so that jumps of $\nu$ only occur at interfaces). For all $T \in \mathcal{T}_h$, we define the discrete symmetric gradient operator $\mathbf{G}_{s,T}^k := \frac{1}{2}(\mathbf{G}_T^k + (\mathbf{G}_T^k)^\mathrm{T})$ (with $\mathbf{G}_T^k$ defined by (15a)) and we use instead of (16) the velocity reconstruction such that, for all $\underline{\mathbf{v}}_T \in \underline{\mathcal{U}}_T^k$,

$$ \int_T \nabla_s \mathbf{r}_T^{k+1} \underline{\mathbf{v}}_T : \nabla_s \mathbf{w} = - \int_T \mathbf{v}_T \cdot \operatorname{div}(\nabla_s \mathbf{w}) + \sum_{F \in \mathcal{F}_T} \int_F \mathbf{v}_F \cdot (\nabla_s \mathbf{w} \cdot \mathbf{n}_{TF}) \quad \forall \mathbf{w} \in \mathbb{P}^{k+1}(T)^d \quad (31a) $$

and

$$ \int_T (\mathbf{r}_T^{k+1} \underline{\mathbf{v}}_T - \mathbf{v}_T) = 0, \quad \int_T \nabla_s \mathbf{r}_T^{k+1} \underline{\mathbf{v}}_T = \frac{1}{2} \sum_{F \in \mathcal{F}_T} \int_F (\mathbf{n}_{TF} \otimes \mathbf{v}_F - \mathbf{v}_F \otimes \mathbf{n}_{TF}). \quad (31b) $$

Letting $\mathbf{G}_{s,h}^k : \underline{\mathcal{U}}_h^k \to \mathbb{P}^{k+1}(\mathcal{T}_h)^{d \times d}$ be such that, for all $\underline{\mathbf{v}}_h \in \underline{\mathcal{U}}_h^k$, $(\mathbf{G}_{s,h}^k \underline{\mathbf{v}}_h)|_T = \mathbf{G}_{s,T}^k \underline{\mathbf{v}}_T$, the viscous term in (44a) below is discretized by means of the bilinear form

$$ a_{\nu,h}(\underline{\mathbf{u}}_h, \underline{\mathbf{v}}_h) := \int_{\Omega} 2\nu \mathbf{G}_{s,h}^k \underline{\mathbf{u}}_h : \mathbf{G}_{s,h}^k \underline{\mathbf{v}}_h + s_{\nu,h}(\underline{\mathbf{u}}_h, \underline{\mathbf{v}}_h), $$

with stabilization bilinear form

$$ s_{\nu,h}(\underline{\mathbf{u}}_h, \underline{\mathbf{v}}_h) := \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_T} \frac{2\nu_T}{h_F} \int_F \delta_{TF}^k \underline{\mathbf{u}}_T \cdot \delta_{TF}^k \underline{\mathbf{v}}_T, $$

where $\nu_T := \nu|_T \in \mathbb{P}^0(T)$ and $\delta_{TF}^k$ is formally defined as in (27) but using the velocity reconstruction operator defined by (31). In the analysis, the main difference with respect to constant kinematic viscosities is that the polynomial degree $k$ should be taken $\ge 1$ in order to ensure coercivity by a discrete Korn inequality (cf. in particular [21, Lemma 4] for insight into this point).
## 3.4 Convective term

For the discretization of the convective term, we consider here the following trilinear form on $\underline{U}_h^k \times \underline{U}_h^k \times \underline{U}_h^k$ expressed in terms of the discrete gradient operator $G_h^{2k}$:

$$t_h(\underline{w}_h, \underline{u}_h, \underline{v}_h) := \frac{1}{2} \int_{\Omega} \underline{v}_h^T G_h^{2k} \underline{u}_h \underline{w}_h - \frac{1}{2} \int_{\Omega} \underline{u}_h^T G_h^{2k} \underline{v}_h \underline{w}_h. \quad (32)$$

This expression mimicks the continuous one given in (3) with $G_h^{2k}$ replacing the continuous gradient operator. Notice that, in the practical implementation, one does not need to actually compute $G_h^{2k}$ to evaluate $t_h$. Instead, the following expression can be used, obtained by applying (15a) twice to expand the terms involving $G_h^{2k}$:

$$t_h(\underline{w}_h, \underline{u}_h, \underline{v}_h) = \sum_{T \in \mathcal{T}_h} t_T(\underline{w}_T, \underline{u}_T, \underline{v}_T),$$

where, for all $T \in \mathcal{T}_h$,

$$t_T(\underline{w}_T, \underline{u}_T, \underline{v}_T) := -\frac{1}{2} \int_T \underline{u}_T^T \nabla \underline{v}_T \underline{w}_T + \frac{1}{2} \int_T \underline{v}_T^T \nabla \underline{u}_T \underline{w}_T \\ + \frac{1}{2} \sum_{F \in \mathcal{F}_T} \int_F (\underline{u}_F \cdot \underline{v}_T)(\underline{w}_T \cdot \underline{n}_{TF}) - \frac{1}{2} \sum_{F \in \mathcal{F}_T} \int_F (\underline{v}_F \cdot \underline{u}_T)(\underline{w}_T \cdot \underline{n}_{TF}). \quad (33)$$

The proof of the following result is postponed to Section 5.2.

**Proposition 6 (Properties of $t_h$).** The trilinear form $t_h$ has the following properties:

1) Skew-symmetry. For all $\underline{v}_h, \underline{w}_h \in \underline{U}_{h,0}^k$, it holds

$$t_h(\underline{w}_h, \underline{v}_h, \underline{v}_h) = 0. \quad (34)$$

2) Boundedness. For all $\underline{u}_h, \underline{v}_h, \underline{w}_h \in \underline{U}_{h,0}^k$, it holds

$$|t_h(\underline{w}_h, \underline{u}_h, \underline{v}_h)| \le C_t \|\underline{w}_h\|_{1,h} \|\underline{u}_h\|_{1,h} \|\underline{v}_h\|_{1,h}, \quad (35)$$

with real number $C_t > 0$ independent of $h$.

3) Consistency. For all $\underline{v} \in \underline{U} \cap H^{k+2}(\Omega)^d$ such that $\mathrm{div} \, \underline{v} = 0$, it holds

$$\sup_{\underline{w}_h \in \underline{U}_{h,0}^k, \|\underline{w}_h\|_{1,h}=1} \left| \int_{\Omega} \underline{w}_h^T \nabla \underline{v} \, \underline{v} - t_h(\underline{I}_h^k \underline{v}, \underline{I}_h^k \underline{v}, \underline{w}_h) \right| \lesssim h^{k+1} \|\underline{v}\|_{H^2(\Omega)^d} \|\underline{v}\|_{H^{k+2}(\Omega)^d}. \quad (36)$$

4) Sequential consistency. Let $(\underline{v}_h)_{h \in \mathcal{H}}$ denote a sequence in $(\underline{U}_{h,0}^k)_{h \in \mathcal{H}}$ bounded in the $\|\cdot\|_{1,h}$-norm with limit $\underline{v} \in \underline{U}$ (cf. point 3) in Proposition 1). Then, for all $\underline{\Phi} \in C_c^\infty(\Omega)^d$ it holds, as $h \to 0$, up to a subsequence,

$$t_h(\underline{v}_h, \underline{v}_h, \underline{I}_h^k \underline{\Phi}) \to t(\underline{v}, \underline{v}, \underline{\Phi}). \quad (37)$$

Some remarks are of order.

**Remark 7 (Design guidelines).** The trilinear form $t_h$ appears in the analysis carried out in Section 4 only through its properties detailed in Proposition 6, with the sole exception of **Step 4** in the proof of Theorem 14 (strong convergence of the pressure), which requires a more intimate use of its expression. Such properties can therefore be intended as design guidelines.

**Remark 8 (Comparison with a HDG discretization).** A discretization of the convective term inspired by the recent HDG literature is obtained replacing $t_h$ by

$$t_h^{\mathrm{HDG}}(\underline{w}_h, \underline{u}_h, \underline{v}_h) := \sum_{T \in \mathcal{T}_h} t_T^{\mathrm{HDG}}(\underline{w}_T, \underline{u}_T, \underline{v}_T), \quad (38)$$
where, for all $T \in \mathcal{T}_h$,

$$t_T^{\text{HDG}}(\underline{w}_T, \underline{u}_T, \underline{v}_T) := -\frac{1}{2} \int_T \underline{u}_T^\text{T} \nabla \underline{v}_T \underline{w}_T + \frac{1}{2} \int_T \underline{v}_T^\text{T} \nabla \underline{u}_T \underline{w}_T \\ + \frac{1}{2} \sum_{F \in \mathcal{F}_T} \int_F (\underline{u}_F \cdot \underline{v}_T)(\underline{w}_F \cdot \underline{n}_{TF}) - \frac{1}{2} \sum_{F \in \mathcal{F}_T} \int_F (\underline{v}_F \cdot \underline{u}_T)(\underline{w}_F \cdot \underline{n}_{TF}) \\ + \frac{\eta}{2} \sum_{F \in \mathcal{F}_T} \int_F (\underline{u}_F - \underline{u}_T) \cdot (\underline{v}_F - \underline{v}_T) |\underline{w}_F \cdot \underline{n}_{TF}|.$$

This discretization has been recently proposed in [41] (cf. Definition 3.3 therein and also [11]), where a HDG method is considered with element-based DOFs that are polynomials of degree ($k+1$) (recall that here we use polynomials of degree $k$, cf. (9)) and the viscous term is discretized as in [37, 40] in order to improve the convergence rates to match the ones of HHO methods (cf. [13] for further details, in particular Remark 2.2). Comparing the above expression of $t_T^{\text{HDG}}$ with (33), we observe the following differences: (i) $w_F$ replaces $w_T$ in both terms in the second line; (ii) a nonnegative stabilization corresponding to the term in the third line is added, including a user-dependent parameter $\eta \ge 0$ (taken equal to 1 in [41]). It can be proved that all the properties listed in Proposition 6 hold for $t_h^{\text{HDG}}$ with $\eta = 0$.

## 3.5 Pressure-velocity coupling

The pressure-velocity coupling is realized by means of the bilinear form $b_h$ on $\underline{U}_h^k \times P_h^k$ such that, for all $(\underline{v}_h, q_h) \in \underline{U}_h^k \times P_h^k$,

$$b_h(\underline{v}_h, q_h) := - \int_{\Omega} D_h^k \underline{v}_h q_h. \qquad (39)$$

The proof of the following result is postponed to Section 5.3.

**Proposition 9 (Properties of $b_h$).** The bilinear form $b_h$ has the following properties:

1) Inf-sup stability. For all $q_h \in P_h^k$ (with $P_h^k$ defined by (13)), it holds

$$\|q_h\|_{L^2(\Omega)} \lesssim \sup_{\underline{v}_h \in \underline{U}_{h,0}^k, \|\underline{v}_h\|_{1,h}=1} b_h(\underline{v}_h, q_h). \qquad (40)$$

2) Consistency. For all $q \in H^{k+1}(\Omega)$ it holds

$$\sup_{\underline{v}_h \in \underline{U}_{h,0}^k, \|\underline{v}_h\|_{1,h}=1} \left| \int_{\Omega} \nabla q \cdot \underline{v}_h - b_h(\underline{v}_h, \pi_h^k q) \right| \lesssim h^{k+1} \|q\|_{H^{k+1}(\Omega)}. \qquad (41)$$

3) Sequential consistency. We have sequential consistency for $b_h$ in the following sense:

* Let $(q_h)_{h \in \mathcal{H}}$ denote a sequence in $(P_h^k)_{h \in \mathcal{H}}$ bounded in the $L^2(\Omega)$-norm and weakly converging to $q \in P$ as $h \to 0$, up to a subsequence. Then, for all $\Phi \in C_c^\infty(\Omega)^d$ it holds

$$b_h(\underline{I}_h^k \Phi, q_h) \to b(\Phi, q). \qquad (42)$$

* Let $(\underline{v}_h)_{h \in \mathcal{H}}$ denote a sequence in $(\underline{U}_{h,0}^k)_{h \in \mathcal{H}}$ bounded in the $\|\cdot\|_{1,h}$-norm with limit $\underline{v} \in \underline{U}$ (cf. point 3) in Proposition 1). Then, for all $\varphi \in C_c^\infty(\Omega)$ it holds

$$b_h(\underline{v}_h, \pi_h^k \varphi) \to b(\underline{v}, \varphi). \qquad (43)$$
## 3.6 Discrete problem

The discrete problem reads: Find $(\underline{u}_h, p_h) \in \underline{U}_{h,0}^k \times P_h^k$ such that

$$ \nu a_h(\underline{u}_h, \underline{v}_h) + t_h(\underline{u}_h, \underline{u}_h, \underline{v}_h) + b_h(\underline{v}_h, p_h) = \int_{\Omega} \mathbf{f} \cdot \mathbf{v}_h \quad \forall \underline{v}_h \in \underline{U}_{h,0}^k, \quad (44a) $$

$$ -b_h(\underline{u}_h, q_h) = 0 \quad \forall q_h \in P_h^k. \quad (44b) $$

Remark 10 (Efficient implementation). When solving the system of nonlinear algebraic equations corresponding to (44) by a first-order (Newton-like) algorithm, all element-based velocity DOFs and all but one pressure DOF per element can be locally eliminated at each iteration by computing the corresponding Schur complement element-wise. As all the computations are local, this static condensation procedure is a trivially parallel task which can fully benefit from multi-thread and multi-processor architectures. For the details, we refer to [23, Section 6.2], where the Stokes problem is considered (the only variation here is that also the linearized convective term appears in the matrices therein denoted by $A_T$). As a result, after the elimination to boundary DOFs corresponding to Dirichlet boundary conditions, we end up solving at each iteration a linear system of size

$$ d \operatorname{card}(\mathcal{F}_h^i) \binom{k+d-1}{d-1} + \operatorname{card}(\mathcal{T}_h). $$

# 4 Analysis of the method

In this section we study the existence and uniqueness of the solution to the HHO scheme (44), prove convergence to the exact solution for general data, and derive convergence rates under a standard data smallness assumption.

## 4.1 Existence and uniqueness

The existence of a solution to problem (44) can be proved using the following topological degree lemma (cf., e.g., [16]), as originally proposed in [27] in the context of finite volumes for nonlinear hyperbolic problems; see also [19, 30] for the Navier-Stokes equations.

**Lemma 11** (Topological degree). Let $W$ be a finite-dimensional functional space equipped with a norm $\|\cdot\|_W$, and let the function $\Psi: W \times [0, 1] \to W$ satisfy the following assumptions:

1) $\Psi$ is continuous;

2) There exists $\mu > 0$ such that, for any $(w, \rho) \in W \times [0, 1]$, $\Psi(w, \rho) = 0$ implies $\|w\|_W \neq \mu$;

3) $\Psi(\cdot, 0)$ is an affine function and the equation $\Psi(w, 0) = 0$ has a solution $w \in W$ such that $\|w\|_W < \mu$. Then, there exists $w \in V$ such that $\Psi(w, 1) = 0$ and $\|w\|_W < \mu$.

**Theorem 12** (Existence and a priori bounds). There exists a solution $(\underline{u}_h, p_h) \in \underline{U}_{h,0}^k \times P_h^k$ to (44), which satisfies the a priori bounds

$$ \|\underline{u}_h\|_{1,h} \le C_a C_s \nu^{-1} \|\mathbf{f}\|_{L^2(\Omega)^d}, \quad \|p_h\|_{L^2(\Omega)} \le C \left( \|\mathbf{f}\|_{L^2(\Omega)^d} + \nu^{-2} \|\mathbf{f}\|_{L^2(\Omega)^d}^2 \right), \quad (45) $$

with $C_a$ and $C_s$ as in (29) and (14), respectively, and $C > 0$ real number independent of both $h$ and $\nu$.

*Proof.* We consider the finite-dimensional space $W_h^k := \underline{U}_{h,0}^k \times P_h^k$ equipped with the norm

$$ \| (\underline{u}_h, r_h) \|_{W,h} := \| \underline{w}_h \|_{1,h} + \nu^{-1} \| r_h \|_{L^2(\Omega)} $$
and the function $\Psi: W_h^k \times [0, 1] \to W_h^k$ such that, for given $(\underline{w}_h, r_h) \in W_h^k$ and $\rho \in [0, 1]$, $(\xi_h, \zeta_h) = \Psi((\underline{w}_h, r_h), \rho)$ is defined as the unique solution of

$$ (\xi_h, \underline{v}_h)_{0,h} = \nu a_h(\underline{w}_h, \underline{v}_h) + \rho t_h(\underline{w}_h, \underline{w}_h, \underline{v}_h) + b_h(\underline{v}_h, r_h) - \int_{\Omega} \mathbf{f} \cdot \mathbf{v}_h \quad \forall \underline{v}_h \in \underline{U}_{h,0}^k, \quad (46a) $$

$$ \int_{\Omega} \zeta_h q_h = -b_h(\underline{w}_h, q_h) \qquad \forall q_h \in P_h^k, \quad (46b) $$

where $(\cdot, \cdot)_{0,h}$ is the $L^2$-like scalar product on $W_h^k$ defined by

$$ (\underline{w}_h, \underline{v}_h)_{0,h} := \int_{\Omega} \underline{w}_h \cdot \underline{v}_h + \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_T} h_F \int_F (\underline{w}_F - \underline{w}_T) \cdot (\underline{v}_F - \underline{v}_T). $$

We next check the assumptions of the topological degree lemma.

1) Since $W_h^k$ is a finite-dimensional space, the bilinear forms $a_h$ and $b_h$, the trilinear form $t_h$, and the scalar products are continuous, and so is the case for the function $\Psi$.

2) Let $(\underline{w}_h, r_h) \in W_h^k$ be such that $\Psi((\underline{w}_h, r_h), \rho) = (\underline{0}, 0)$ for some $\rho \in [0, 1]$. We next show that

$$ \| (\underline{w}_h, r_h) \|_{W,h} \le (C_a C_s + C) \nu^{-1} \| \mathbf{f} \|_{L^2(\Omega)^d} + C \nu^{-3} \| \mathbf{f} \|_{L^2(\Omega)^d}^2, $$

and point 2) in Lemma 11 is verified for

$$ \mu = (C_a C_s + C) \nu^{-1} \| \mathbf{f} \|_{L^2(\Omega)^d} + C \nu^{-3} \| \mathbf{f} \|_{L^2(\Omega)^d}^2 + \epsilon $$

with $\epsilon > 0$. Recalling the coercivity of $a_h$ expressed by the first inequality in (29), making $\underline{v}_h = \underline{w}_h$ in (46a) and observing that $t_h(\underline{w}_h, \underline{w}_h, \underline{w}_h) = 0$ owing to skew-symmetry (34) and that $b_h(\underline{w}_h, r_h) = 0$ owing to (46b) with $q_h = r_h$, we have

$$ \nu C_a^{-1} \| \underline{w}_h \|_{1,h}^2 \le \nu \| \underline{w}_h \|_{a,h}^2 = \int_{\Omega} \mathbf{f} \cdot \mathbf{w}_h \le \| \mathbf{f} \|_{L^2(\Omega)^d} \| \underline{w}_h \|_{L^2(\Omega)^d} \le C_s \| \mathbf{f} \|_{L^2(\Omega)^d} \| \underline{w}_h \|_{1,h}, $$

where we have used the discrete Poincaré inequality (14) with $p = 2$ to conclude. The bound on $\underline{w}_h$ follows. To prove the bound on $r_h$, we proceed as follows:

$$ \begin{aligned} \| r_h \|_{L^2(\Omega)} &\lesssim \sup_{\underline{v}_h \in \underline{U}_{h,0}^k, \| \underline{v}_h \|_{1,h} = 1} b_h(\underline{v}_h, r_h) \\ &= \sup_{\underline{v}_h \in \underline{U}_{h,0}^k, \| \underline{v}_h \|_{1,h} = 1} \left( \int_{\Omega} \mathbf{f} \cdot \mathbf{v}_h - \nu a_h(\underline{w}_h, \underline{v}_h) - \rho t_h(\underline{w}_h, \underline{w}_h, \underline{v}_h) \right) \\ &\lesssim \| \mathbf{f} \|_{L^2(\Omega)^d} + \nu \| \underline{w}_h \|_{1,h} + \rho \| \underline{w}_h \|_{1,h}^2 \lesssim \| \mathbf{f} \|_{L^2(\Omega)^d} + \nu^{-2} \| \mathbf{f} \|_{L^2(\Omega)^d}^2, \end{aligned} $$

where we have used the inf-sup condition (40) on $b_h$ in the first line and (46a) to pass to the second line; the Cauchy–Schwarz and the discrete Poincaré inequalities together with the boundedness of $a_h$ and $t_h$ expressed by the second inequality in (29) and by (35), respectively, are used to pass to the third line; the bound on the velocity and the fact that $\rho \le 1$ allow to conclude.

3) $\Psi(\cdot, 0)$ is an affine function from $W_h^k$ to $W_h^k$. The fact that $\Psi(\cdot, 0)$ is invertible corresponds to the well-posedness of the HHO scheme for the Stokes problem, and can therefore be proved using the arguments of [23, Lemma 3] (which classically rely on the coercivity of $a_h$ expressed by the first inequality in (29) and the inf-sup condition (40) for $b_h$). Additionally, the unique solution $(\underline{w}_h, r_h) \in W_h^k$ to the equation $\Psi((\underline{w}_h, r_h), 0) = 0$ satisfies $\| (\underline{w}_h, r_h) \|_{W,h} < \mu$ as a consequence of point 2).

The existence of a solution to (44) is an immediate consequence of Lemma 11. Observing that, if $(\underline{w}_h, p_h) \in W_h^k$ solves (44), then $\Psi((\underline{w}_h, p_h), 1) = (\underline{0}, 0)$, the bounds (45) follow from point 2) above. $\square$

We next consider uniqueness, which can be classically proved under a data smallness condition.
**Theorem 13** (Uniqueness of the discrete solution). Assume that the right-hand side verifies

$$ \|\boldsymbol{f}\|_{L^2(\Omega)^d} \le \frac{\nu^2}{2C_a^2 C_t C_s} \qquad (47) $$

with $C_a, C_t$ and $C_s$ as in (29), (35), (14) respectively. Then, the solution $(\underline{\boldsymbol{u}}_h, p_h) \in \underline{U}_{h,0}^k \times P_h^k$ of (44) is unique.

**Proof.** Let $(\underline{\boldsymbol{u}}_{1,h}, p_{1,h}) \in \underline{U}_{h,0}^k \times P_h^k$ and $(\underline{\boldsymbol{u}}_{2,h}, p_{2,h}) \in \underline{U}_{h,0}^k \times P_h^k$ solve (44), and let

$$ \underline{\boldsymbol{w}}_h := \underline{\boldsymbol{u}}_{1,h} - \underline{\boldsymbol{u}}_{2,h} \quad \text{and} \quad r_h := p_{1,h} - p_{2,h}. $$

Taking the difference of the discrete momentum balance equation (44a) written for $(\underline{\boldsymbol{u}}_h, p_h) = (\underline{\boldsymbol{u}}_{1,h}, p_{1,h})$ and $(\underline{\boldsymbol{u}}_h, p_h) = (\underline{\boldsymbol{u}}_{2,h}, p_{2,h})$, we infer that it holds for all $\underline{\boldsymbol{v}}_h \in \underline{U}_{h,0}^k$,

$$ \nu a_h(\underline{\boldsymbol{w}}_h, \underline{\boldsymbol{v}}_h) + t_h(\underline{\boldsymbol{u}}_{1,h}, \underline{\boldsymbol{w}}_h, \underline{\boldsymbol{v}}_h) + t_h(\underline{\boldsymbol{w}}_h, \underline{\boldsymbol{u}}_{2,h}, \underline{\boldsymbol{v}}_h) + b_h(\underline{\boldsymbol{v}}_h, r_h) = 0. \qquad (48) $$

Making $\underline{\boldsymbol{v}}_h = \underline{\boldsymbol{w}}_h$ in the above equation, observing that $t_h(\underline{\boldsymbol{u}}_{1,h}, \underline{\boldsymbol{w}}_h, \underline{\boldsymbol{w}}_h) = 0$ owing to skew-symmetry (cf. point 1) in Proposition 6), that $b_h(\underline{\boldsymbol{w}}_h, r_h) = 0$ (this is a consequence of the discrete mass balance equation (44b) written for $\underline{\boldsymbol{u}}_{1,h}$ and $\underline{\boldsymbol{u}}_{2,h}$ with $q_h = r_h$), and using the boundedness (29) of $a_h$ and (35) of $t_h$, we obtain

$$ (\nu C_a^{-1} - C_t \|\underline{\boldsymbol{u}}_{2,h}\|_{1,h}) \|\underline{\boldsymbol{w}}_h\|_{1,h}^2 \le 0. $$

By the a priori bounds (45) and the assumption (47) on $\boldsymbol{f}$, the first factor in the left-hand side is > 0. As a result, we infer $\underline{\boldsymbol{w}}_h = \underline{\boldsymbol{0}}$, thus proving uniqueness for the velocity. Plugging this result into (48), it is inferred that for all $\underline{\boldsymbol{v}}_h \in \underline{U}_{h,0}^k$ it holds $b_h(\underline{\boldsymbol{v}}_h, r_h) = 0$. The inf-sup stability (cf. point 1) in Proposition 9) then gives

$$ \|r_h\|_{L^2(\Omega)} \lesssim \sup_{\underline{\boldsymbol{v}}_h \in \underline{U}_{h,0}^k, \|\underline{\boldsymbol{v}}_h\|_{1,h}=1} b_h(\underline{\boldsymbol{v}}_h, r_h) = 0, $$

which proves uniqueness for the pressure. $\square$

## 4.2 Convergence to minimal regularity solutions

**Theorem 14** (Convergence to minimal regularity solutions). Let $(\mathcal{T}_h)_{h \in \mathcal{H}}$ denote and admissible mesh sequence as in Section 2, and let $((\underline{\boldsymbol{u}}_h, p_h))_{h \in \mathcal{H}}$ be such that, for all $h \in \mathcal{H}$, $(\underline{\boldsymbol{u}}_h, p_h) \in \underline{U}_{h,0}^k \times P_h^k$ solves (44). Then, it holds up to a subsequence with $(\underline{\boldsymbol{u}}, p) \in \mathcal{U} \times \mathcal{P}$ solution of the continuous problem (2),

1) $\underline{\boldsymbol{u}}_h \to \boldsymbol{u}$ strongly in $L^p(\Omega)^d$ for all $p \in [1, +\infty)$ if $d = 2$, $p \in [1, 6)$ if $d = 3$;

2) $\mathbf{G}_h^k \underline{\boldsymbol{u}}_h \to \nabla \boldsymbol{u}$ strongly in $L^2(\Omega)^{d \times d}$;

3) $s_h(\underline{\boldsymbol{u}}_h, \underline{\boldsymbol{u}}_h) \to 0$;

4) $p_h \to p$ strongly in $L^2(\Omega)$.

If, in addition, the solution to (2) is unique (which is the case if the smallness condition detailed in [32, Eq. (2.12), Chapter IV] holds for $\boldsymbol{f}$), convergence extends to the whole sequence.

**Proof.** The proof proceeds in several steps. In **Step 1** we prove the existence of a limit for the sequence of discrete solutions. In **Step 2** we show that this limit is indeed a solution of the continuous problem (2). In **Step 3** we prove the strong convergence of the velocity gradient and of the jumps, and in **Step 4** the strong convergence of the pressure.

**Step 1. Existence of a limit.** Since, for all $h \in \mathcal{H}$, $(\underline{\boldsymbol{u}}_h, p_h) \in \underline{U}_{h,0}^k \times P_h^k$ solves (44), we infer combining the a priori bound (45) and point 3) in Proposition 1 that there exists $(\boldsymbol{u}, p) \in \mathcal{U} \times \mathcal{P}$ such that, up to a subsequence as $h \to 0$,

(i) $\underline{\boldsymbol{u}}_h \to \boldsymbol{u}$ strongly in $L^p(\Omega)^d$ for all $p \in [1, +\infty)$ if $d = 2$, $p \in [1, 6)$ if $d = 3$;

(ii) $\mathbf{G}_h^l \underline{\boldsymbol{u}}_h \to \nabla \boldsymbol{u}$ weakly in $L^2(\Omega)^{d \times d}$ for all $l \ge 0$;
(iii) $p_h \to p$ weakly in $L^2(\Omega)$.

**Step 2.** *Identification of the limit.* We next prove that $(\mathbf{u}, p) \in U \times P$ is a solution of (2). Let $\mathbf{\Phi} \in C_c^\infty(\Omega)^d$. We apply the sequential consistency of the viscous, convective and pressure terms (respectively expressed by point 3) in Proposition 3, point 4) in Proposition 6 and point 3) in Proposition 9) to infer

$$ \nu a_h(\underline{\mathbf{u}}_h, \underline{\mathbf{I}}_h^k \mathbf{\Phi}) + t_h(\underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h, \underline{\mathbf{I}}_h^k \mathbf{\Phi}) + b_h(\underline{\mathbf{I}}_h^k \mathbf{\Phi}, p_h) \to \nu a(\mathbf{u}, \mathbf{\Phi}) + t(\mathbf{u}, \mathbf{u}, \mathbf{\Phi}) + b(\mathbf{\Phi}, p). $$

Furthermore, we have $\pi_h^k \mathbf{\Phi} \to \mathbf{\Phi}$ strongly in $L^2(\Omega)^d$, which implies

$$ \int_{\Omega} f \cdot \pi_h^k \mathbf{\Phi} \to \int_{\Omega} f \cdot \mathbf{\Phi}. $$

Finally, point 3) of Proposition 9 gives for all $\varphi \in C_c^\infty(\Omega)$

$$ b_h(\underline{\mathbf{u}}_h, \pi_h^k \varphi) \to b(\mathbf{u}, \varphi). $$

As a result, we can conclude by density that $(\mathbf{u}, p) \in U \times P$ is a solution of (2) and point 1) is proved.

**Step 3.** *Strong convergence of the velocity gradient and of the jumps.* Making $\underline{\mathbf{v}}_h = \underline{\mathbf{u}}_h$ in (44a) and observing that $t_h(\underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h) = 0$ owing to skew-symmetry (34) and that $b_h(\underline{\mathbf{u}}_h, p_h) = 0$ owing to (44b) with $q_h = p_h$, we have

$$ \nu \| \mathbf{G}_h^k \mathbf{u}_h \|_{L^2(\Omega)^{d \times d}}^2 \le \nu a_h(\underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h) = \int_{\Omega} f \cdot \mathbf{u}_h. $$

Since $\mathbf{u}_h$ converges to $\mathbf{u}$ strongly in $L^2(\Omega)^d$ and $\mathbf{u}$ is a solution of (2), we have

$$ \nu \limsup \| \mathbf{G}_h^k \mathbf{u}_h \|_{L^2(\Omega)^{d \times d}}^2 \le \limsup \int_{\Omega} f \cdot \mathbf{u}_h = \int_{\Omega} f \cdot \mathbf{u} = \nu \| \nabla \mathbf{u} \|_{L^2(\Omega)^{d \times d}}^2. $$

This estimate combined with the weak convergence of $\mathbf{G}_h^k \mathbf{u}_h$ to $\nabla \mathbf{u}$ implies the strong convergence of the velocity gradient $\mathbf{G}_h^k \underline{\mathbf{u}}_h$ to $\nabla \mathbf{u}$ in $L^2(\Omega)^{d \times d}$. On the other hand, we also obtain that $a_h(\underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h)$ converges to $\| \nabla \mathbf{u} \|_{L^2(\Omega)^{d \times d}}^2$, and finally we get

$$ s_h(\underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h) = a_h(\underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h) - \int_{\Omega} \mathbf{G}_h^k \underline{\mathbf{u}}_h : \mathbf{G}_h^k \underline{\mathbf{u}}_h \to 0. \quad (49) $$

This proves points 2) and 3).

**Step 4.** *Strong convergence of the pressure.* Observing that $p_h \in P$, from the surjectivity of the continuous divergence operator from $U$ to $P$ we infer the existence of $\mathbf{v}_{p_h} \in U$ such that

$$ \mathrm{div} \, \mathbf{v}_{p_h} = p_h \quad \text{and} \quad \| \mathbf{v}_{p_h} \|_{H^1(\Omega)^d} \lesssim \| p_h \|_{L^2(\Omega)}. \quad (50) $$

We let, for all $h \in \mathcal{H}$, $\hat{\mathbf{v}}_{p_h, h} := \underline{\mathbf{I}}_h^k \mathbf{v}_{p_h}$, and study the properties of the sequence $(\hat{\mathbf{v}}_{p_h, h})_{h \in \mathcal{H}}$. For all $h \in \mathcal{H}$, it holds

$$ \| \hat{\mathbf{v}}_{p_h, h} \|_{1, h} \lesssim \| \mathbf{v}_{p_h} \|_{H^1(\Omega)^d} \lesssim \| p_h \|_{L^2(\Omega)} \lesssim \| \mathbf{f} \|_{L^2(\Omega)^d} + \nu^{-2} \| \mathbf{f} \|_{L^2(\Omega)^d}^2, \quad (51) $$

where we have used the boundedness (12) of $\underline{\mathbf{I}}_h^k$ in the first inequality, (50) in the second, and the a priori bound (45) on the pressure to conclude. Then, by point 3) in Proposition 1, there exists $\mathbf{v}_p \in U$ such that $\hat{\mathbf{v}}_{p_h, h} \to \mathbf{v}_p$ strongly in $L^p(\Omega)^d$ for all $p \in [1, 4]$ and $\mathbf{G}_h^l \hat{\mathbf{v}}_{p_h, h} \to \nabla \mathbf{v}_p$ weakly in $L^2(\Omega)^{d \times d}$ for all $l \ge 0$. Moreover, by uniqueness of the limit in the distribution sense, it holds that

$$ \mathrm{div} \, \mathbf{v}_p = p \quad (52) $$

Making $\underline{\mathbf{v}}_h = \hat{\mathbf{v}}_{p_h, h}$ in the discrete momentum balance equation (44a) and recalling the commuting property (25), we have

$$ \| p_h \|_{L^2(\Omega)}^2 = -b_h(\hat{\mathbf{v}}_{p_h, h}, p_h) = \nu a_h(\underline{\mathbf{u}}_h, \hat{\mathbf{v}}_{p_h, h}) + t_h(\underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h, \hat{\mathbf{v}}_{p_h, h}) - \int_{\Omega} f \cdot \hat{\mathbf{v}}_{p_h, h}. \quad (53) $$
We study the limit of the three terms on the right of (53) using the convergence properties for the discrete solution proved in the previous points. Combining the strong converge of $G_h^k \hat{u}_h$ with the weak convergence of $G_h^k \hat{v}_{p_h,h}$ gives

$$ \int_{\Omega} G_h^k \underline{u}_h : G_h^k \hat{v}_{p_h,h} \to \int_{\Omega} \nabla \underline{u} : \nabla v_p. $$

Moreover, the convergence (49) of the jumps of $\underline{u}_h$ and the uniform bound (51) imply

$$ s_h(\underline{u}_h, \hat{v}_{p_h,h}) \le s_h(\underline{u}_h, \underline{u}_h)^{\frac{1}{2}} \times s_h(\hat{v}_{p_h,h}, \hat{v}_{p_h,h})^{\frac{1}{2}} \to 0, $$

so that, in conclusion, we have for the viscous term

$$ a_h(\underline{u}_h, \hat{v}_{p_h,h}) \to a(\underline{u}, v_p). $$

Observing that the convergence properties of the sequences $(\underline{u}_h)_{h \in \mathcal{H}}$ and $(\hat{v}_{p_h,h})_{h \in \mathcal{H}}$ are respectively analogous to those of the sequences $(\underline{I}_h^k \Phi)_{h \in \mathcal{H}}$ and $(\underline{v}_h)_{h \in \mathcal{H}}$ in point 4) of Proposition 6, we can prove proceeding in a similar way that

$$ t_h(\underline{u}_h, \underline{u}_h, \hat{v}_{p_h,h}) \to t(\underline{u}, \underline{u}, v_p). $$

Finally, by strong convergence of $\hat{v}_{p_h,h}$ to $v_p$ in $L^2(\Omega)^d$, we readily infer for the source term

$$ \int_{\Omega} f \cdot \hat{v}_{p_h,h} \to \int_{\Omega} f \cdot v_p. $$

Collecting the above convergence results and using the momentum balance equation (2a) together with (52) leads to

$$ \limsup \|p_h\|_{L^2(\Omega)}^2 \le \nu a(\underline{u}, v_p) + t(\underline{u}, \underline{u}, v_p) - \int_{\Omega} f \cdot v_p = -b(v_p, p) = \|p\|_{L^2(\Omega)}^2, $$

and the strong convergence of the pressure in $L^2(\Omega)$ stated in point 4) follows. $\square$

## 4.3 Convergence rates for small data

**Theorem 15** (Convergence rates for small data). Let $(\underline{u}, p) \in \mathcal{U} \times P$ and $(\underline{u}_h, p_h) \in \underline{\mathcal{U}}_{h,0}^k \times P_h^k$ solve problems (2) and (44), respectively, and assume uniqueness (which holds, in particular, if both smallness conditions [32, Eq. (2.12), Chapter IV] and (47) are verified). Assume, moreover, the additional regularity $(\underline{u}, p) \in H^{k+2}(\Omega)^d \times H^{k+1}(\Omega)$, as well as

$$ \|f\|_{L^2(\Omega)^d} \le \frac{\nu^2}{2C_I C_a C_t (1 + C_P^2)}, \quad (54) $$

with $C_I, C_a$ and $C_t$ as in (12), (29) and (35), respectively, and $C_P$ Poincaré constant only depending on $\Omega$ such that, for all $\underline{v} \in \mathcal{U}$, $\|\underline{v}\|_{L^2(\Omega)^d} \le C_P \|\nabla \underline{v}\|_{L^2(\Omega)^{d \times d}}$. Let

$$ \underline{u}_h := \underline{I}_h^k \underline{u}, \qquad \hat{p}_h := \pi_h^k p. $$

Then, there is $C > 0$ independent of both $h$ and $\nu$ such that

$$ \|\underline{u}_h - \hat{\underline{u}}_h\|_{1,h} + \nu^{-1} \|p_h - \hat{p}_h\|_{L^2(\Omega)} \le C h^{k+1} \left( (1 + \nu^{-1} \|\underline{u}\|_{H^2(\Omega)^d}) \|\underline{u}\|_{H^{k+2}(\Omega)^d} + \nu^{-1} \|p\|_{H^{k+1}(\Omega)} \right). \quad (55) $$

**Corollary 16** (Convergence rates for small data). Under the above assumptions, it holds

$$ \|G_h^k \underline{u}_h - \nabla \underline{u}\|_{L^2(\Omega)^{d \times d}} + s_h(\underline{u}_h, \underline{u}_h)^{\frac{1}{2}} + \nu^{-1} \|p_h - p\|_{L^2(\Omega)} \lesssim h^{k+1} \left( (1 + \nu^{-1} \|\underline{u}\|_{H^2(\Omega)^d}) \|\underline{u}\|_{H^{k+2}(\Omega)^d} + \nu^{-1} \|p\|_{H^{k+1}(\Omega)} \right), $$

where the second term in the left-hand side accounts for the jumps of the discrete solution.
Proof. Using the triangle inequality, we infer

$$
\begin{aligned}
& \| \mathbf{G}_h^k \underline{\mathbf{u}}_h - \nabla \mathbf{u} \|_{L^2(\Omega)^{d \times d}} + s_h(\underline{\mathbf{u}}_h, \underline{\mathbf{u}}_h)^{\frac{1}{2}} + \nu^{-1} \| p_h - p \|_{L^2(\Omega)} \le \\
& \qquad \| \nabla \mathbf{u} - \mathbf{G}_h^k \hat{\underline{\mathbf{u}}}_h \|_{L^2(\Omega)^{d \times d}} + s_h(\hat{\underline{\mathbf{u}}}_h, \hat{\underline{\mathbf{u}}}_h)^{\frac{1}{2}} + \nu^{-1} \| p - \hat{p}_h \|_{L^2(\Omega)} \\
& \qquad + \| \underline{\mathbf{u}}_h - \hat{\underline{\mathbf{u}}}_h \|_{a,h} + \nu^{-1} \| p_h - \hat{p}_h \|_{L^2(\Omega)}.
\end{aligned}
$$

The terms in the second line can be estimated recalling the consistency property (20) for the gradient reconstruction and using the approximation properties (5) of the L²-orthogonal projector and the consistency properties (28) of $s_h$. For the terms in the third line, recall the norm equivalence (29) and use (55). □

Remark 17 (Extension to other hybrid discretizations). The following proof extends without modifications to any bilinear forms $a_h$ and $b_h$ and trilinear form $t_h$ that match, respectively, the properties 1) and 2) in Proposition 3, 1)–3) in Proposition 9, and 1) and 2) in Proposition 6, respectively. Such properties can therefore be intended as design guidelines.

*Proof of Theorem 15.* Let, for the sake of brevity, $(\underline{e}_h, \epsilon_h) := (\underline{u}_h - \hat{\underline{u}}_h, p_h - \hat{p}_h)$. The proof proceeds in three steps: in **Step 1**, we identify the consistency error and derive a lower bound in terms of $\|\underline{e}_h\|_{1,h}$ using the data smallness assumption, in **Step 2** we estimate the error on the velocity and in **Step 3** the error on the pressure.

**Step 1. Consistency error and lower bound.** It is readily inferred from the discrete momentum balance equation (44a) that it holds, for all $\underline{v}_h \in \underline{U}_{h,0}^k$,

$$ \nu a_h(\underline{e}_h, \underline{v}_h) + t_h(\underline{u}_h, \underline{u}_h, \underline{v}_h) - t_h(\hat{\underline{u}}_h, \hat{\underline{u}}_h, \underline{v}_h) + b_h(\underline{v}_h, \epsilon_h) = \mathcal{E}_h(\underline{v}_h), \quad (56) $$

with consistency error

$$ \mathcal{E}_h(\underline{v}_h) := (f, v_h) - a_h(\hat{\underline{u}}_h, \underline{v}_h) - t_h(\hat{\underline{u}}_h, \hat{\underline{u}}_h, \underline{v}_h) - b_h(\underline{v}_h, \hat{p}_h). $$

Making $\underline{v}_h = \underline{e}_h$ in (56), and observing that $t_h(\underline{u}_h, \underline{u}_h, \underline{e}_h) = t_h(\underline{u}_h, \hat{\underline{u}}_h, \underline{e}_h)$ owing to the skew-symmetry property (34), and that $b_h(\underline{e}_h, \epsilon_h) = b_h(\underline{u}_h, \epsilon_h) - b_h(\hat{\underline{u}}_h, \epsilon_h) = 0$ owing to (44b) and since $D_h^k \hat{\underline{u}}_h = \pi_h^k(\text{div } \mathbf{u}) = 0$ (cf. (25) and (1b)), we infer

$$
\begin{aligned}
\mathcal{E}_h(\underline{e}_h) &= \nu \|\underline{e}_h\|_{a,h}^2 + t_h(\underline{e}_h, \hat{\underline{u}}_h, \underline{e}_h) \\
&\geq \nu C_a^{-1} \|\underline{e}_h\|_{1,h}^2 - C_t \|\hat{\underline{u}}_h\|_{1,h} \|\underline{e}_h\|_{1,h}^2 \quad (57) \\
&\geq (\nu C_a^{-1} - C_t C_I (1 + C_P^2) \nu^{-1}) \|\mathbf{f}\|_{H^{-1}(\Omega)^d} \|\underline{e}_h\|_{1,h}^2 \geq \nu \|\underline{e}_h\|_{1,h}^2,
\end{aligned}
$$

where we have used the coercivity of $a_h$ expressed by the first inequality in (29) together with the boundedness (35) of $t_h$ to pass to the second line, the boundedness (12) of $I_h^k$ with the standard a priori estimate $\|\mathbf{u}\|_{H^1(\Omega)^d} \le (1 + C_P^2)\nu^{-1}\|\mathbf{f}\|_{H^{-1}(\Omega)^d}$ on the exact solution to infer

$$ \|\hat{\underline{u}}_h\|_{1,h} \le C_I \|\mathbf{u}\|_{H^1(\Omega)^d} \le C_I (1 + C_P^2) \nu^{-1} \|\mathbf{f}\|_{H^{-1}(\Omega)^d}, \quad (58) $$

and the data smallness assumption (54) to conclude.

**Step 2. Estimate on the velocity.** Observing that $\mathbf{f} = -\nu\Delta\mathbf{u} + \nabla\mathbf{u}\mathbf{u} + \nabla p$ a.e. in $\Omega$ (cf. (1a)), it holds for all $\underline{v}_h \in \underline{U}_{h,0}^k$,

$$ |\mathcal{E}_h(\underline{v}_h)| \le \nu \left| \int_{\Omega} (\Delta \mathbf{u}) \cdot \mathbf{v}_h + a_h(\hat{\underline{u}}_h, \underline{v}_h) \right| + \left| \int_{\Omega} \mathbf{v}_h^\mathrm{T} \nabla \mathbf{u} \mathbf{u} - t_h(\hat{\underline{u}}_h, \hat{\underline{u}}_h, \underline{v}_h) \right| + \left| \int_{\Omega} \nabla p \cdot \mathbf{v}_h - b_h(\underline{v}_h, \hat{p}_h) \right|. $$

Using (30), (36) and (41), respectively, to estimate the three terms in the right-hand side, it is readily inferred that

$$ S := \sup_{\underline{v}_h \in \underline{U}_{h,0}^k, \|\underline{v}_h\|_{1,h}=1} |\mathcal{E}_h(\underline{v}_h)| \lesssim \nu h^{k+1} (1 + \nu^{-1} \|\mathbf{u}\|_{H^2(\Omega)^d}) \|\mathbf{u}\|_{H^{k+2}(\Omega)^d} + h^{k+1} \|p\|_{H^{k+1}(\Omega)}, \quad (59) $$

so that, in particular,

$$ |\mathcal{E}_h(\underline{e}_h)| \lesssim S \|\underline{e}_h\|_{1,h} \lesssim [\nu h^{k+1} (1 + \nu^{-1} \|\mathbf{u}\|_{H^2(\Omega)^d}) \|\mathbf{u}\|_{H^{k+2}(\Omega)^d} + h^{k+1} \|p\|_{H^{k+1}(\Omega)}] \|\underline{e}_h\|_{1,h}. \quad (60) $$

Combining (57) with (60), the estimate on the velocity in (55) follows.
**Step 3.** *Estimate on the pressure.* Let us now estimate the error on the pressure. We have

$$
\begin{align}
\|\epsilon_h\|_{L^2(\Omega)} &\lesssim \sup_{\underline{v}_h \in \underline{U}_{h,0}^k, \|\underline{v}_h\|_{1,h}=1} b_h(\underline{v}_h, \epsilon_h) \nonumber \\
&= \sup_{\underline{v}_h \in \underline{U}_{h,0}^k, \|\underline{v}_h\|_{1,h}=1} [\mathcal{E}_h(\underline{v}_h) - \nu a_h(\underline{e}_h, \underline{v}_h) - t_h(\underline{u}_h, \underline{u}_h, \underline{v}_h) + t_h(\hat{\underline{u}}_h, \hat{\underline{u}}_h, \underline{v}_h)] \nonumber \\
&= \sup_{\underline{v}_h \in \underline{U}_{h,0}^k, \|\underline{v}_h\|_{1,h}=1} [\mathcal{E}_h(\underline{v}_h) - \nu a_h(\underline{e}_h, \underline{v}_h) - t_h(\underline{e}_h, \underline{u}_h, \underline{v}_h) - t_h(\hat{\underline{u}}_h, \underline{e}_h, \underline{v}_h)] \tag{61} \\
&\lesssim S + \nu (1 + \nu^{-1}\|\underline{u}_h\|_{1,h} + \nu^{-1}\|\hat{\underline{u}}_h\|_{1,h}) \|\underline{e}_h\|_{1,h} \nonumber \\
&\lesssim S + \nu (1 + \nu^{-2}\|\underline{f}\|_{L^2(\Omega)^d} + \nu^{-2}\|\underline{f}\|_{H^{-1}(\Omega)^d}) \|\underline{e}_h\|_{1,h} \lesssim S + \nu \|\underline{e}_h\|_{1,h}, \nonumber
\end{align}
$$

In (61), we have used the inf-sup inequality (40) on $b_h$ in the first line and the error equation (56) to pass to the second line; to pass to the third line, we have inserted $\pm t_h(\hat{\underline{u}}_h, \underline{u}_h, \underline{v}_h)$ and used the linearity of $t_h$ in its first and second arguments; to pass to the fourth line, we have used the boundedness (29) of $a_h$ and (35) of $t_h$; to pass to the fifth line, we have used the a priori bounds (45) on $\|\underline{u}_h\|_{1,h}$ and (58) on $\|\hat{\underline{u}}_h\|_{1,h}$; the data smallness assumption (54) gives the conclusion. The estimate on the pressure then follows using (59) and (55), respectively, to further bound the addends in the right-hand side of (61). $\square$

## 4.4 Numerical example

To close this section, we provide a numerical example that demonstrates the convergence properties of our method. We solve on the two-dimensional square domain $\Omega := (-0.5, 1.5) \times (0, 2)$ the Dirichlet problem corresponding to the exact solution $(\underline{u}, p)$ of [36] with $\underline{u} = (u_1, u_2)$ such that, introducing the Reynolds number $Re := (2\nu)^{-1}$ and letting $\lambda := Re - (\text{Re}^2 + 4\pi^2)^{\frac{1}{2}}$,

$$u_1(\underline{x}) := 1 - \exp(\lambda x_1) \cos(2\pi x_2), \quad u_2(\underline{x}) := \frac{\lambda}{2\pi} \exp(\lambda x_1) \sin(2\pi x_2),$$

and pressure given by

$$p(\underline{x}) := -\frac{1}{2} \exp(2\lambda x_1) + \frac{\lambda}{2} (\exp(4\lambda) - 1).$$

We take here $\nu = 1$ and consider two sequences of refined meshes obtained by linearly mapping on $\Omega$ the mesh family 2 of [34] and the (predominantly) hexagonal mesh family of [24] (both meshes were originally defined on the unit square). The implementation uses the static condensation procedure discussed in Remark 10. The convergence results for $k = 2$ and $k = 3$ are reported in Figures 1 and 2, respectively. Using the notation of Theorem 15, we separately plot the $H^1$-error on the velocity $\|\underline{u}_h - \hat{\underline{u}}_h\|_{1,h}$, the $L^2$-error on the pressure $\|p_h - \hat{p}_h\|_{L^2(\Omega)}$, and the $L^2$-error on the velocity $\|\underline{u}_h - \hat{\underline{u}}_h\|_{L^2(\Omega)^d}$. For the sake of completeness, we consider both the trilinear forms $t_h$ given by (32) and $t_h^{\text{HDG}}$ given by (38) (with $\eta = 0$). In both cases, we obtain similar results, and the $H^1$-error on the velocity as well as the $L^2$-error on the pressure converge as $h^{k+1}$ as expected. The $L^2$-error on the velocity, on the other hand, converges as $h^{k+2}$. Notice also that this means that the error $\|\underline{r}_h^{k+1}\underline{u}_h - \underline{u}\|_{L^2(\Omega)^d}$ can be proved to converge as $h^{k+2}$ following a similar reasoning as in [2, Corollary 4.6]. The details are omitted for the sake of brevity.

# 5 Properties of the discrete bilinear and trilinear forms

We gather in this section the proofs of Propositions 3, 6 and 9.

## 5.1 Viscous bilinear form

*Proof of Proposition 3.* We only sketch the proof and provide references for the details.

1) This norm equivalence follows taking $p = 2$ in [17, Lemma 5.2], where the scalar case is considered; cf. also [22, Lemma 4], where a slightly different expression for $a_T$ is studied (cf. Remark 4).
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>(a) $t_h, k = 2$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) $t_h, k = 3$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) $t_h^{\text{HDG}}, \eta = 0, k = 2$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) $t_h^{\text{HDG}}, \eta = 0, k = 3$</figcaption></figure>

Figure 1: Cartesian mesh family, errors versus $h$. The triangles indicate reference slopes. The trilinear forms $t_h$ and $t_h^{\text{HDG}}$ are defined by (32) and (38), respectively.

2) We adapt the arguments of [22, Theorem 8]. For the sake of brevity, we let $\hat{v}_h := I_h^k v$ in what follows. Integrating by parts element-by-element, and using the single-valuedness of $\nabla v \cdot n_F$ at interfaces and the fact that $w_F = 0$ on boundary faces to insert $w_F$ into the second term, we have

$$ \int_{\Omega} \Delta v \cdot w_h = - \sum_{T \in \mathcal{T}_h} \left( \int_T \nabla v : \nabla w_T + \sum_{F \in \mathcal{F}_T} \int_F (w_F - w_T)^{\mathrm{T}} \nabla v \cdot n_{TF} \right). \quad (62) $$

On the other hand, using on each element $T \in \mathcal{T}_h$ the definition (15b) of $G_T^k$ (with $v_T = w_T$ and $\tau = G_T^k \hat{v}_T$), we have

$$ a_h(\hat{v}_h, \underline{w}_h) = \sum_{T \in \mathcal{T}_h} \left( \int_T G_T^k \hat{v}_T : \nabla w_T + \sum_{F \in \mathcal{F}_T} \int_F (w_F - w_T)^{\mathrm{T}} G_T^k \hat{v}_T n_{TF} \right) + s_h(\hat{v}_h, \underline{w}_h). \quad (63) $$

Summing (62) and (63), observing that the first terms in parentheses cancel out as a result of the Euler equation (22) for $G_T^k$, and using the Cauchy–Schwarz inequality followed by the trace approximation properties (20) of $G_T^k$, the consistency properties (28) of $s_h$, and the norm equivalence (29), we get

$$ \begin{aligned} & \left| \int_{\Omega} \Delta v \cdot w_h + a_h(\hat{v}_h, \underline{w}_h) \right| \\ & \lesssim \left( \sum_{T \in \mathcal{T}_h} h_T \| G_T^k \hat{v}_T - \nabla v \|_{L^2(\partial T)^{d \times d}}^2 + s_h(\hat{v}_h, \underline{v}_h) \right)^{\frac{1}{2}} \times \left( \sum_{T \in \mathcal{T}_h} |\underline{w}_T|^2_{1, \partial T} + s_h(\underline{w}_h, \underline{w}_h) \right)^{\frac{1}{2}} \\ & \lesssim h^{k+1} \|v\|_{H^{k+2}(\Omega)^d} \|\underline{w}_h\|_{1,h}, \end{aligned} $$
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>(a) $t_h, k = 2$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) $t_h, k = 3$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) $t_h^{\text{HDG}}, \eta = 0, k = 2$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(d) $t_h^{\text{HDG}}, \eta = 0, k = 3$</figcaption></figure>

Figure 2: Hexagonal mesh family, errors versus $h$. The triangles indicate reference slopes. The trilinear forms $t_h$ and $t_h^{\text{HDG}}$ are defined by (32) and (38), respectively.

which concludes the proof of (30).

3) The sequential consistency can be proved following steps 1) and 2) of [17, Theorem 4.6], where the scalar case is considered. $\square$

## 5.2 Convective trilinear form

*Proof of Proposition 6.* 1) Skew-symmetry. This property is straightforward from the definition of $t_h$.

2) Boundedness. For all $\underline{w}_h, \underline{u}_h, \underline{v}_h \in \underline{U}_{h,0}^k$, using Hölder inequalities, we have

$$ \begin{aligned} |t_h(\underline{w}_h, \underline{u}_h, \underline{v}_h)| & \lesssim \|\underline{v}_h\|_{L^4(\Omega)^d} \|\underline{G}_h^{2k} \underline{u}_h\|_{L^2(\Omega)^{d \times d}} \|\underline{w}_h\|_{L^4(\Omega)^d} + \|\underline{u}_h\|_{L^4(\Omega)^d} \|\underline{G}_h^{2k} \underline{v}_h\|_{L^2(\Omega)^{d \times d}} \|\underline{w}_h\|_{L^4(\Omega)^d} \\ & \lesssim \|\underline{v}_h\|_{1,h} \|\underline{u}_h\|_{1,h} \|\underline{w}_h\|_{1,h}, \end{aligned} $$

where the conclusion follows using several times the discrete Sobolev embedding (14) with $p = 4$ and the boundedness (19) of $\underline{G}_h^{2k}$.

3) Consistency. Set, for the sake of brevity, $\hat{\underline{v}}_h := \underline{I}_h^k \underline{v}$. Integrating by parts element-by-element, recalling that $\text{div} \underline{v} = 0$, and using the single-valuedness of $(\underline{v} \cdot \underline{n}_F)\underline{v}$ at interfaces together with the fact that $\underline{w}_F = \underline{0}$ on boundary faces to insert $\underline{w}_F$ into the third term, we have

$$ \int_{\Omega} \underline{w}_h^\mathrm{T} \nabla \underline{v} \underline{v} = \frac{1}{2} \left( \int_{\Omega} \underline{w}_h^\mathrm{T} \nabla \underline{v} \underline{v} - \int_{\Omega} \underline{v}^\mathrm{T} \nabla_h \underline{w}_h \underline{v} - \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_T} \int_F (\underline{v} \cdot \underline{n}_{TF})(\underline{w}_F - \underline{w}_T) \cdot \underline{v} \right). \quad (64) $$
On the other hand, using on each element $T \in \mathcal{T}_h$ the definition (15b) of $G_T^{2k}$ (with $\underline{v}_T = \underline{w}_T$ and $\tau = \hat{v}_T \otimes \hat{v}_T$), we have

$$t_h(\hat{\underline{v}}_h, \hat{\underline{v}}_h, \underline{w}_h) = \frac{1}{2} \left( \int_{\Omega} \underline{w}_h^T G_h^{2k} \hat{\underline{v}}_h \hat{\underline{v}}_h - \int_{\Omega} \hat{\underline{v}}_h^T \nabla_h \underline{w}_h \hat{\underline{v}}_h - \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_T} \int_F (\hat{\underline{v}}_T \cdot n_{TF})(\underline{w}_F - \underline{w}_T) \cdot \hat{\underline{v}}_T \right). \quad (65)$$

Subtracting (65) from (64) and inserting in the right-hand side of the resulting expression the quantity

$$\pm \frac{1}{2} \left( \int_{\Omega} \underline{w}_h^T G_h^{2k} \hat{\underline{v}}_h \underline{v} + \int_{\Omega} \underline{v}^T \nabla_h \underline{w}_h \hat{\underline{v}}_h + \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_T} \int_F (\underline{v} \cdot n_{TF})(\underline{w}_F - \underline{w}_T) \cdot \hat{\underline{v}}_T \right),$$

we arrive at

$$\int_{\Omega} \underline{w}_h^T \nabla \underline{v} \underline{v} - t_h(\underline{I}_h^k \underline{v}, \underline{I}_h^k \underline{v}, \underline{w}_h) = \frac{1}{2} \int_{\Omega} \underline{w}_h^T (\nabla \underline{v} - G_h^{2k} \hat{\underline{v}}_h) \underline{v} + \frac{1}{2} \int_{\Omega} \underline{w}_h^T G_h^{2k} \hat{\underline{v}}_h (\underline{v} - \hat{\underline{v}}_h) \\ + \frac{1}{2} \int_{\Omega} (\hat{\underline{v}}_h - \underline{v})^T \nabla_h \underline{w}_h \hat{\underline{v}}_h + \frac{1}{2} \int_{\Omega} \underline{v}^T \nabla_h \underline{w}_h (\hat{\underline{v}}_h - \underline{v}) \\ + \frac{1}{2} \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_T} \int_F ((\hat{\underline{v}}_T - \underline{v}) \cdot n_{TF})(\underline{w}_F - \underline{w}_T) \cdot \hat{\underline{v}}_T \\ + \frac{1}{2} \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_T} \int_F (\underline{v} \cdot n_{TF})(\underline{w}_F - \underline{w}_T) \cdot (\hat{\underline{v}}_T - \underline{v}). \quad (66)$$

Denote by $\mathcal{T}_1, \dots, \mathcal{T}_6$ the addends in the right-hand side of the above expression. For the first term, using for all $T \in \mathcal{T}_h$ the Euler equation (22) with $\tau = \underline{w}_T \otimes \underline{\pi}_T^0 \underline{v} \in \mathbb{P}^k(T)^{d \times d}$, we infer

$$\mathcal{T}_1 = \frac{1}{2} \int_{\Omega} \underline{w}_h^T (\nabla \underline{v} - G_h^{2k} \hat{\underline{v}}_h)(\underline{v} - \underline{\pi}_h^0 \underline{v}).$$

Hence, using the Hölder inequality followed by the approximation properties (20) of $G_h^{2k}$ and (5) of $\underline{\pi}_h^0$ (with $m=0, p=4, s=1$), we obtain

$$|\mathcal{T}_1| \lesssim h^{k+1} \|\underline{w}_h\|_{L^4(\Omega)^d} \|\underline{v}\|_{H^{k+1}(\Omega)^a} \|\underline{v}\|_{W^{1,4}(\Omega)^d} \lesssim h^{k+1} \|\underline{w}_h\|_{1,h} \|\underline{v}\|_{H^{k+1}(\Omega)^d} \|\underline{v}\|_{H^2(\Omega)^d}, \quad (67)$$

where the conclusion follows using the discrete Sobolev embedding (14) with $p=4$ to bound the first factor and the continuous injection $H^2(\Omega) \to W^{1,4}(\Omega)$ valid in $d \in \{2,3\}$ on domains satisfying the cone condition to bound the third (cf. [1, Theorem 4.12]).

Using again the Hölder inequality, the boundedness (19) of $G_h^{2k}$ and (12) of $\underline{I}_h^k$ to infer $\|G_h^{2k} \hat{\underline{v}}_h\|_{L^2(\Omega)^{d \times d}} \lesssim \|\underline{v}\|_{H^1(\Omega)^d}$, and the approximation properties (5) of $\underline{\pi}_h^k$ (with $m=0, p=4$, and $s=k+1$), we infer

$$|\mathcal{T}_2| \lesssim h^{k+1} \|\underline{w}_h\|_{L^4(\Omega)^d} \|\underline{v}\|_{H^1(\Omega)^d} \|\underline{v}\|_{W^{k+1,4}(\Omega)^d} \lesssim h^{k+1} \|\underline{w}_h\|_{1,h} \|\underline{v}\|_{H^1(\Omega)^d} \|\underline{v}\|_{H^{k+2}(\Omega)^d}, \quad (68)$$

where the conclusion follows from the discrete Sobolev embedding (14) with $p=4$ together with the continuous injection $H^{k+2}(\Omega) \to W^{k+1,4}(\Omega)$ valid for all $k \ge 0$ and $d \in \{2,3\}$ on domains satisfying the cone condition (cf. [1, Theorem 4.12]).

Proceeding similarly, we have for the third and fourth terms

$$|\mathcal{T}_3| + |\mathcal{T}_4| \lesssim h^{k+1} \|\nabla_h \underline{w}_h\|_{L^2(\Omega)^{d \times d}} (\|\hat{\underline{v}}_h\|_{L^4(\Omega)^d} + \|\underline{v}\|_{L^4(\Omega)^d}) \|\underline{v}\|_{W^{k+1,4}(\Omega)^d} \\ \lesssim h^{k+1} \|\underline{w}_h\|_{1,h} \|\underline{v}\|_{H^1(\Omega)^d} \|\underline{v}\|_{H^{k+2}(\Omega)^d}, \quad (69)$$

where, to pass to the second line, we have used the definition (10) of the $\|\cdot\|_{1,h}$-norm to bound the first factor, the discrete and continuous Sobolev embeddings to estimate the $L^4(\Omega)^d$-norms in the second factor, the boundedness (12) of $\underline{I}_h^k$ to further bound $\|\hat{\underline{v}}_h\|_{1,h} \lesssim \|\underline{v}\|_{H^1(\Omega)^d}$, and the continuous injection $H^{k+2}(\Omega) \to W^{k+1,4}(\Omega)$ to conclude.
Finally, for the fifth and sixth term, using Hölder inequalities and the trace approximation properties (5) of the $L^2$-orthogonal projector (with $m = 0, p = 4$, and $s = k + 1$), we obtain

$$ \begin{aligned} |\mathcal{T}_5| + |\mathcal{T}_6| &\lesssim h^{k+1} \|\mathbf{v}\|_{W^{k+1,4}(\Omega)^d} \left( \sum_{T \in \mathcal{T}_h} |\underline{\mathbf{w}}|_{1, \partial T}^2 \right)^{\frac{1}{2}} \times \left( \sum_{T \in \mathcal{T}_h} h_T (\|\mathbf{v}\|_{L^4(\partial T)^d}^4 + \|\hat{\mathbf{v}}\|_{L^4(\partial T)^d}^4) \right)^{\frac{1}{4}} \\ &\lesssim h^{k+1} \|\mathbf{v}\|_{H^{k+2}(\Omega)^d} \|\underline{\mathbf{w}}_h\|_{1,h} \|\mathbf{v}\|_{H^2(\Omega)^d}, \end{aligned} \quad (70) $$

where, to pass to the second line, we have used the continuous injection $H^{k+2}(\Omega) \to W^{k+1,4}(\Omega)$ for the first factor, the definition (10) of the $\|\cdot\|_{1,h}$-norm for the second factor, and the continuous (6) and discrete (8) trace inequalities with $p = 4$ followed by the continuous injection $H^2(\Omega) \to W^{1,4}(\Omega)$ for the third factor. Taking absolute values in (66), and using (67)-(70) to bound the right-hand side, (36) follows.

4) Sequential consistency. We have, letting for the sake of brevity $\hat{\mathbf{\Phi}}_h := \underline{\mathbf{I}}_h^k \mathbf{\Phi}$,

$$ t_h(\underline{\mathbf{v}}_h, \underline{\mathbf{v}}_h, \hat{\mathbf{\Phi}}_h) = \frac{1}{2} \int_{\Omega} \hat{\mathbf{\Phi}}_h^\mathrm{T} \mathbf{G}_h^{2k} \underline{\mathbf{v}}_h \mathbf{v}_h - \frac{1}{2} \int_{\Omega} \mathbf{v}_h^\mathrm{T} \mathbf{G}_h^{2k} \hat{\mathbf{\Phi}}_h \mathbf{v}_h := \mathcal{T}_1 + \mathcal{T}_2. $$

Since $\mathbf{v}_h \to \mathbf{v}$ and $\hat{\mathbf{\Phi}}_h \to \mathbf{\Phi}$ strongly in $L^4(\Omega)^d$, $\hat{\mathbf{\Phi}}_h \otimes \mathbf{v}_h \to \mathbf{\Phi} \otimes \mathbf{v}$ strongly in $L^2(\Omega)^{d \times d}$. Hence, recalling that $\mathbf{G}_h^{2k} \underline{\mathbf{v}}_h \to \nabla \mathbf{v}$ weakly in $L^2(\Omega)^{d \times d}$ owing to point 3) in Proposition 1, we infer that

$$ \mathcal{T}_1 = \frac{1}{2} \int_{\Omega} \mathbf{G}_h^{2k} \underline{\mathbf{v}}_h : (\hat{\mathbf{\Phi}}_h \otimes \mathbf{v}_h) \to \frac{1}{2} \int_{\Omega} \mathbf{\Phi}^\mathrm{T} \nabla \mathbf{v} \mathbf{v}. $$

For the second term, observing that $\mathbf{v}_h \otimes \mathbf{v}_h \to \mathbf{v} \otimes \mathbf{v}$ and $\mathbf{G}_h^{2k} \hat{\mathbf{\Phi}}_h \to \nabla \mathbf{\Phi}$ strongly in $L^2(\Omega)^{d \times d}$, we readily get

$$ \mathcal{T}_2 = \frac{1}{2} \int_{\Omega} \mathbf{G}_h^{2k} \hat{\mathbf{\Phi}}_h : (\mathbf{v}_h \otimes \mathbf{v}_h) \to \frac{1}{2} \int_{\Omega} \mathbf{v}_h^\mathrm{T} \nabla \mathbf{\Phi} \mathbf{v}_h. $$

The conclusion follows from the above results recalling the definition (3) of $t$. $\square$

## 5.3 Velocity-pressure coupling bilinear form

*Proof of Proposition 9.* 1) Inf-sup stability. We deploy similar arguments as in [7, Lemma 4] and [23, Lemma 3]. Let $q_h \in P_h^k$ and denote by $S$ the supremum in (40). Observing that $q_h \in P$, from the surjectivity of the continuous divergence operator from $U$ to $P$ we infer the existence of $\mathbf{v}_{q_h} \in U$ such that $\mathrm{div} \, \mathbf{v}_{q_h} = q_h$ and $\|\mathbf{v}_{q_h}\|_{H^1(\Omega)^d} \lesssim \|q_h\|_{L^2(\Omega)}$. Then, we have

$$ \|q_h\|_{L^2(\Omega)}^2 = -b_h(\underline{\mathbf{I}}_h^k \mathbf{v}_{q_h}, q_h) \le S \|\underline{\mathbf{I}}_h^k \mathbf{v}_{q_h}\|_{1,h} \le S \|\mathbf{v}_{q_h}\|_{H^1(\Omega)^d} \le S \|q_h\|_{L^2(\Omega)}, $$

where we have used the commuting property (25) for $D_h^k$, the definition of the supremum, the boundeness (12) of $\underline{\mathbf{I}}_h^k$, and $\|\mathbf{v}_{q_h}\|_{H^1(\Omega)^d} \lesssim \|q_h\|_{L^2(\Omega)}$ to conclude.

2) Consistency. Integrating by parts element-by-element, and using the fact that the jumps of $q$ vanish at interfaces by the assumed regularity and that $\mathbf{v}_F = \mathbf{0}$ on boundary faces to insert $\mathbf{v}_F$ into the second term, we have

$$ \int_{\Omega} \nabla q \cdot \mathbf{v}_h = - \sum_{T \in \mathcal{T}_h} \left( \int_T q(\mathrm{div} \, \mathbf{v}_T) + \sum_{F \in \mathcal{F}_T} \int_F q(\mathbf{v}_F - \mathbf{v}_T) \cdot \mathbf{n}_{TF} \right). \quad (71) $$

On the other hand, using (17b) on each element $T \in \mathcal{T}_h$ to express the right-hand side of (39), we have

$$ -b_h(\underline{\mathbf{v}}_h, \pi_h^k q) = \sum_{T \in \mathcal{T}_h} \left( \int_T \pi_T^k q(\mathrm{div} \, \mathbf{v}_T) + \sum_{F \in \mathcal{F}_T} \int_F \pi_T^k q(\mathbf{v}_F - \mathbf{v}_T) \cdot \mathbf{n}_{TF} \right). \quad (72) $$

Summing (71) and (72), observing that the first terms in parentheses cancel out by the definition (4) of $\pi_T^k$ since $(\mathrm{div} \, \mathbf{v}_T) \in \mathbb{P}^{k-1}(T) \subset \mathbb{P}^k(T)$, and using for the second terms the Cauchy-Schwarz inequality
followed by the trace approximation properties (5) of $\pi_T^k$, we infer that

$$ \left| \int_{\Omega} \nabla q \cdot \mathbf{v}_h - b_h(\mathbf{v}_h, \pi_h^k q) \right| \le \left( \sum_{T \in \mathcal{T}_h} h_T \|\pi_T^k q - q\|_{L^2(\partial T)}^2 \right)^{\frac{1}{2}} \left( \sum_{T \in \mathcal{T}_h} |\mathbf{v}_T|_{1, \partial T}^2 \right)^{\frac{1}{2}} \lesssim h^{k+1} \|q\|_{H^{k+1}(\Omega)} \|\mathbf{v}_h\|_{1,h}. $$

Passing to the supremum in the above expression, (41) follows.

3) Sequential consistency. Recalling (18), $D_h^k = \text{tr}(\mathbf{G}_h^k)$ and the sequential consistency (42) is a straightforward consequence of point 2) in Proposition 1 combined with a weak-strong convergence argument. Similarly, the sequential consistency (43) follows from the fact that $D_h^k \mathbf{v}_h \to \text{div} \mathbf{v}$ weakly in $L^2(\Omega)$ as a consequence of point 3) in Proposition 1 and $\pi_h^k \varphi \to \varphi$ strongly in $L^2(\Omega)$. □

## References

[1] R. A. Adams and J. J. F. Fournier. Sobolev spaces, volume 140 of Pure and Applied Mathematics (Amsterdam). Elsevier/Academic Press, Amsterdam, second edition, 2003.

[2] J. Aghili, S. Boyaval, and D. A. Di Pietro. Hybridization of mixed high-order methods on general meshes and application to the Stokes equations. Comput. Meth. Appl. Math., 15(2):111–134, 2015.

[3] F. Bassi, L. Botti, A. Colombo, and S. Rebay. Agglomeration based discontinuous Galerkin discretization of the Euler and Navier-Stokes equations. Comput. & Fluids, 61:77–85, 2012.

[4] F. Bassi, A. Crivellini, D. A. Di Pietro, and S. Rebay. An artificial compressibility flux for the discontinuous Galerkin solution of the incompressible Navier-Stokes equations. J. Comput. Phys., 218(2):794–815, 2006.

[5] F. Bassi, A. Crivellini, D. A. Di Pietro, and S. Rebay. An implicit high-order discontinuous Galerkin method for steady and unsteady incompressible flows. Comp. & Fl., 36(10):1529–1546, 2007.

[6] L. Beirão da Veiga, C. Lovadina, and G. Vacca. Divergence free Virtual Elements for the Stokes problem on polygonal meshes. ESAIM: Math. Model. Numer. Anal. (M2AN), 2016. Published online. DOI 10.1051/m2an/2016032.

[7] D. Boffi, M. Botti, and D. A. Di Pietro. A nonconforming high-order method for the Biot problem on general meshes. SIAM J. Sci. Comput., 38(3):A1508–A1537, 2016.

[8] D. Boffi, F. Brezzi, and M. Fortin. Mixed finite element methods and applications, volume 44 of Springer Series in Computational Mathematics. Springer, Heidelberg, 2013.

[9] S. C. Brenner and L. R. Scott. The mathematical theory of finite element methods, volume 15 of Texts in Applied Mathematics. Springer, New York, third edition, 2008.

[10] P. Castillo, B. Cockburn, I. Perugia, and D. Schötzau. An a priori error analysis of the local discontinuous Galerkin method for elliptic problems. SIAM J. Numer. Anal., 38:1676–1706, 2000.

[11] A. Çeşmelioğlu, B. Cockburn, and W. Qiu. Analysis of an HDG method for the incompressible Navier-Stokes equations. Math. Comp., 2016. Published online. DOI: 10.1090/mcom/3195.

[12] C. Chainais-Hillairet, S. Krell, and A. Mouton. Convergence analysis of a DDFV scheme for a system describing miscible fluid flows in porous media. Numer. Methods Partial Differential Equations, 31(3):723–760, 2015.

[13] B. Cockburn, D. A. Di Pietro, and A. Ern. Bridging the Hybrid High-Order and Hybridizable Discontinuous Galerkin methods. ESAIM: Math. Model. Numer. Anal. (M2AN), 50(3):635–650, 2016.

[14] B. Cockburn, J. Gopalakrishnan, and R. Lazarov. Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems. SIAM J. Numer. Anal., 47(2):1319–1365, 2009.

[15] B. Cockburn, G. Kanschat, and D. Schötzau. A locally conservative LDG method for the incompressible Navier-Stokes equations. Math. Comp., 74(251):1067–1095, 2005.

[16] K. Deimling. Nonlinear functional analysis. Springer-Verlag, Berlin, 1985.

[17] D. A. Di Pietro and J. Droniou. A Hybrid High-Order method for Leray–Lions elliptic equations on general meshes. Math. Comp., 2016. Published online. DOI: 10.1090/mcom/3180.

[18] D. A. Di Pietro and J. Droniou. $W^{s,p}$-approximation properties of elliptic projectors on polynomial spaces, with application to the error analysis of a Hybrid High-Order discretisation of Leray–Lions problems, June 2016. Submitted. Preprint arXiv:1606.02832 [math.NA].
[19] D. A. Di Pietro and A. Ern. Discrete functional analysis tools for discontinuous Galerkin methods with application to the incompressible Navier-Stokes equations. *Math. Comp.*, 79:1303–1330, 2010.

[20] D. A. Di Pietro and A. Ern. *Mathematical aspects of discontinuous Galerkin methods*, volume 69 of *Mathématiques & Applications*. Springer-Verlag, Berlin, 2012.

[21] D. A. Di Pietro and A. Ern. A hybrid high-order locking-free method for linear elasticity on general meshes. *Comput. Meth. Appl. Mech. Engrg.*, 283:1–21, 2015.

[22] D. A. Di Pietro, A. Ern, and S. Lemaire. An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators. *Comput. Meth. Appl. Math.*, 14(4):461–472, 2014.

[23] D. A. Di Pietro, A. Ern, A. Linke, and F. Schieweck. A discontinuous skeletal method for the viscosity-dependent Stokes problem. *Comput. Meth. Appl. Mech. Engrg.*, 306:175–195, 2016.

[24] D. A. Di Pietro and S. Lemaire. An extension of the Crouzeix-Raviart space to general meshes with application to quasi-incompressible linear elasticity and Stokes flow. *Math. Comp.*, 84(291):1–31, 2015.

[25] T. Dupont and R. Scott. Polynomial approximation of functions in Sobolev spaces. *Math. Comp.*, 34(150):441–463, 1980.

[26] H. Egger and C. Waluga. hp analysis of a hybrid DG method for Stokes flow. *IMA J. Numer. Anal.*, 33(2):687–721, 2013.

[27] R. Eymard, T. Gallouët, M. Ghilani, and R. Herbin. Error estimates for the approximate solutions of a nonlinear hyperbolic equation given by finite volume schemes. *IMA J. Numer. Anal.*, 18(4):563–594, 1998.

[28] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In *Handbook of numerical analysis, Vol. VII*, Handb. Numer. Anal., VII, pages 713–1020. North-Holland, Amsterdam, 2000.

[29] R. Eymard, T. Gallouët, and R. Herbin. Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces. *IMA J. Numer. Anal.*, 30(4):1009–1043, 2010.

[30] R. Eymard, R. Herbin, and J.-C. Latché. Convergence analysis of a colocated finite volume scheme for the incompressible Navier-Stokes equations on general 2D or 3D meshes. *SIAM J. Numer. Anal.*, 45(1):1–36, 2007.

[31] G. Giorgiani, S. Fernández-Méndez, and A. Huerta. Hybridizable Discontinuous Galerkin with degree adaptivity for the incompressible Navier-Stokes equations. *Computers & Fluids*, 98:196–208, 2014.

[32] V. Girault and P.-A. Raviart. *Finite element methods for Navier-Stokes equations*. Springer-Verlag, Berlin, 1986.

[33] V. Girault, B. Rivière, and M. F. Wheeler. A discontinuous Galerkin method with nonoverlapping domain decomposition for the Stokes and Navier-Stokes problems. *Math. Comp.*, 74(249):53–84, 2005.

[34] R. Herbin and F. Hubert. Benchmark on discretization schemes for anisotropic diffusion problems on general grids. In R. Eymard and J.-M. Hérard, editors, *Finite Volumes for Complex Applications V*, pages 659–692. John Wiley & Sons, 2008.

[35] O. Karakashian and T. Katsaounis. A discontinuous Galerkin method for the incompressible Navier-Stokes equations. In *Discontinuous Galerkin methods (Newport, RI, 1999)*, volume 11 of *Lect. Notes Comput. Sci. Eng.*, pages 157–166. Springer, Berlin, 2000.

[36] L. S. G. Kovasznay. Laminar flow behind a two-dimensional grid. *Proc. Camb. Philos. Soc.*, 44:58–62, 1948.

[37] C. Lehrenfeld. *Hybrid Discontinuous Galerkin methods for solving incompressible flow problems*. PhD thesis, Rheinisch-Westfälischen Technischen Hochschule Aachen, 2010.

[38] C. Liu and N. J. Walkington. Convergence of numerical approximations of the incompressible Navier-Stokes equations with variable density and viscosity. *SIAM J. Numer. Anal.*, 45(3):1287–1304, 2007.

[39] N.C. Nguyen, J. Peraire, and B. Cockburn. An implicit high-order hybridizable discontinuous Galerkin method for the incompressible Navier-Stokes equations. *J. Comput. Phys.*, 230:1147–1170, 2011.

[40] I. Oikawa. A hybridized discontinuous Galerkin method with reduced stabilization. *J. Sci. Comput.*, 65:327–340, 2015.

[41] W. Qiu and K. Shi. A superconvergent HDG method for the incompressible Navier-Stokes equations on general polyhedral meshes. *IMA J. Numer. Anal.*, 36(4):1943–1967, 2016.

[42] B. Rivière and S. Sardar. Penalty-free discontinuous Galerkin methods for incompressible Navier-Stokes equations. *Math. Models Methods Appl. Sci.*, 24(6):1217–1236, 2014.
[43] M. Tavelli and M. Dumbser. A staggered semi-implicit discontinuous Galerkin method for the two dimensional incompressible Navier-Stokes equations. *Appl. Math. Comput.*, 248:70–92, 2014.

[44] M. P. Ueckermann and P. F. J. Lermusiaux. Hybridizable discontinuous Galerkin projection methods for Navier-Stokes and Boussinesq equations. *J. Comput. Phys.*, 306:390–421, 2016.

[45] J. Wang and X. Ye. A weak Galerkin finite element method for the Stokes equations. *Adv. Comput. Math.*, 42:155–174, 2016.