---
abstract: |
    In this paper, we design and analyze a Hybrid High-Order discretization method for the steady motion of non-Newtonian, incompressible fluids in the Stokes approximation of small velocities. The proposed method has several appealing features including the support of general meshes and high-order, unconditional inf-sup stability, and orders of convergence that match those obtained for Leray–Lions scalar problems. A complete well-posedness and convergence analysis of the method is carried out under new, general assumptions on the strain rate-shear stress law, which encompass several common examples such as the power-law and Carreau–Yasuda models. Numerical examples complete the exposition.\
    **Keywords:** Hybrid High-Order methods, non-Newtonian fluids, power-law, Carreau–Yasuda law, discrete Korn inequality
author:
- 'Michele Botti[^1]'
- 'Daniel Castanon Quiroz [^2]'
- 'Daniele A. Di Pietro [^3]'
- 'André Harnist [^4]'
bibliography:
- 'pstokes.bib'
title: 'A Hybrid High-Order method for creeping flows of non-Newtonian fluids'
---

Introduction {#sec:introduction}
============

In this paper, we design and analyze a Hybrid High-Order (HHO) discretization method for the steady motion of a non-Newtonian, incompressible fluid in the Stokes approximation of small velocities. Notable applications include ice sheet dynamics [@Isaac.Stadler.ea:15], mantle convection [@Schubert.Turcotte.ea:01], chemical engineering [@Ko.Pustejovska.ea:18], and biological fluids rheology [@Lai.Kuei.ea:78; @Galdi.Rannacher.ea:18]. We focus on fluids with shear-rate-dependent viscosity, whose behavior is characterized by a nonlinear strain rate-shear stress function. Physical interpretations and discussions of non-Newtonian fluid models can be found, e.g., in [@Bird.Armstrong.ea:87; @Malek.Rajagopal.ea:95]. Typical examples that are frequently used in the applications include the power-law and Carreau–Yasuda model.

The earliest investigations of fluids with shear-dependent viscosities date back to the pioneering work of Ladyzhenskaya [@Ladyzhenskaya:69]. For a detailed mathematical study of the well-posedness and regularity of the continuous problem, see also [@Malek.Rajagopal:05; @Ruzicka.Diening:07; @Diening.Ettwein:08; @Beirao-da-Veiga:09; @Berselli.Ruzicka:20] and references therein. Early results on the numerical analysis of non-Newtonian fluid flow problems were given in [@Sandri:93; @Barrett.Liu:94; @Glowinski.Rappaz:03]. Later, these results were improved in [@Belenki.Berselli.ea:12] and [@Hirn:13] by proving error estimates that are optimal for fluids with shear thinning behavior (described by a power law exponent ${r}\le 2$). In [@Belenki.Berselli.ea:12], the authors considered a conforming inf-sup stable finite element discretization, while in [@Hirn:13] a low-order scheme with local projection stabilization was proposed. In both works, the use of Orlicz functions is instrumental to unify the treatment of the shear thinning and shear thickening cases (also called pseudoplastic and dilatant, respectively; cf. Example \[ex:Carreau–Yasuda\]). More recently, a finite element method based on a four-field formulation of the nonlinear Stokes equations has been analyzed in [@Sandri:14]. Other notable contributions on the numerical approximation of generalized Stokes problems include [@Diening.Kreuzer.ea:13; @Isaac.Stadler.ea:15; @Kreuzer.Suli:16; @Ko.Suli:18].

The main issues to be accounted for in the numerical solution of non-Newtonian fluid flow problems are the presence of local features emerging from the nonlinear strain rate-shear stress relation, the incompressibility condition leading to indefinite systems, the roughly varying model coefficients, and, possibly, complex geometries requiring unstructured and highly-adapted meshes. The HHO method provides several advantages to deal with the complex nature of the problem, such as the support of general polygonal or polyhedral meshes, the possibility to select the approximation order, and unconditional inf-sup stability. Moreover, HHO schemes can be efficiently implemented thanks to the possibility of statically condensing a large subset of the unknowns for linearized versions of the problem encountered, e.g., when solving the nonlinear system by the Newton method. Hybrid High-Order methods have been successfully applied to the simulation of incompressible flows of Newtonian fluids governed by the Stokes [@Aghili.Boyaval.ea:15] and Navier–Stokes equations [@Di-Pietro.Krell:18; @Botti.Di-Pietro.ea:19*1], possibly driven by large irrotational volumetric forces [@Di-Pietro.Ern.ea:16; @Castanon-Quiroz.Di-Pietro:20]. Works related to the problem of creeping flows of non-Newtonian fluids are [@Botti.Di-Pietro.ea:17] and [@Di-Pietro.Droniou:17; @Di-Pietro.Droniou:17*1], respectively dealing with nonlinear elasticity and Leray–Lions problems. Going from nonlinear coercive elliptic equations to the nonlinear Stokes system involves additional difficulties arising from the pressure and the divergence constraint. Finally, we mention, in passing, that HHO methods are members of a wider family of polytopal methods that also includes, e.g., Virtual Element methods; cf., e.g., [@Beirao-da-Veiga.Lovadina.ea:17; @Beirao-da-Veiga.Lovadina.ea:18] for their application to Newtonian incompressible flows.

The HHO discretization presented in this paper is inspired by the previously mentioned works. It hinges on discontinuous polynomial unknowns on the mesh and on its skeleton, from which discrete differential operators are reconstructed. The reconstruction operators are then used to define a consistency term inspired by the weak formulation of the creeping flow problem and a cleverly designed stabilization term penalizing boundary residuals. We carry out a complete analysis of the proposed method. In particular, under general assumptions on the strain rate-shear stress function, we derive error estimates for the velocity and pressure approximations. The energy-norm error estimate for the velocity given in Theorem \[thm:error.estimate\] is optimal in the sense that it yields the same convergence orders established in [@Di-Pietro.Droniou:17*1 Theorem 7] for the scalar Leray–Lions elliptic problem. A key tool in our analysis is provided by Lemma \[lem:discrete.korn.inequality\], in which we prove a generalization of the discrete Korn inequality of [@Botti.Di-Pietro.ea:19 Lemma 1] to the non-Hilbertian case. The other main contributions are a novel formulation of the requirements on the strain rate-shear stress function allowing a unified treatment of pseudoplastic and dilatant fluids and the identification of a set of general assumptions on the nonlinear stabilization function ensuring the desired consistency properties along with the well-posedness of the discrete problem.

The rest of the paper is organized as follows. In Section \[sec:continuous.setting\] we introduce the strong and weak formulations of the nonlinear Stokes problem and present the assumptions on the strain rate-shear stress function. The construction of the HHO discretization is given in Section \[sec:discrete.setting\] by defining the discrete counterparts of the viscous and coupling terms. Section \[sec:discrete.setting\] also contains the proof of the discrete Korn inequality and the discussion on the nonlinear stabilization function. Section \[sec:well-posedness\] establishes the well-posedness of the discrete problem by proving the Hölder continuity and the strong-monotonicity of the viscous term, as well as the inf-sup stability of the pressure-velocity coupling. In Section \[sec:error.estimate\], we show the consistency of the discrete viscous function and coupling bilinear form. These results are then used to prove the error estimate. In Section \[sec:num.res\], we investigate the performance of the method by performing a convergence test with analytical solution on various families of refined meshes. Finally, in Appendix \[sec:properties.stress\] we provide a sufficient condition for the strain rate-shear stress law to fulfil the assumptions presented in Section \[sec:continuous.setting\].

Continuous setting {#sec:continuous.setting}
==================

Let $\Omega \subset {\mathbb{R}}^d$, $d\in\{2,3\}$, denote a bounded, connected, polyhedral open set with Lipschitz boundary $\partial\Omega$. We consider a possibly non-Newtonian fluid occupying $\Omega$ and subjected to a volumetric force field $\b f : \Omega \to {\mathbb{R}}^d$. Its flow is governed by the generalized Stokes problem, which consists in finding the velocity field $\b u : \Omega \to {\mathbb{R}}^d$ and the pressure field $p : \Omega \to {\mathbb{R}}$ such that

\[eq:stokes.continuous\] $$\begin{aligned}
{2} 
    -{\b\nabla{\cdot}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b u) + {\b\nabla}p &= \b f &\qquad& \mbox{  in  } \Omega, \label{eq:stokes.continuous:momentum} \\
    {\b\nabla{\cdot}}\b u &= 0 &\qquad& \mbox{ in }  \Omega, \label{eq:stokes.continuous:mass} \\
    \b u &= \b 0 &\qquad& \mbox{ on } \partial \Omega, \label{eq:stokes.continuous:bc} \\
    \int_\Omega p(\b x){\,\mathrm{d}}\b x &= 0, \label{eq:stokes.continuous:closure}
  \end{aligned}$$

where ${\b\nabla{\cdot}}$ denotes the divergence operator applied to vector fields, ${\b{\nabla}_{\mathrm{s}}}$ is the symmetric part of the gradient operator ${\b\nabla}$ applied to vector fields, and, denoting by ${{\mathbb{R}}^{d \times d}_{\mathrm{s}}}$ the set of square, symmetric, real-valued $d\times d$ matrices, ${\b\sigma}: \Omega \times {{\mathbb{R}}^{d \times d}_{\mathrm{s}}} \to {{\mathbb{R}}^{d \times d}_{\mathrm{s}}}$ is the strain rate-shear stress law. In what follows, we formulate assumptions on ${\b\sigma}$ that encompass common models for non-Newtonian fluids and state a weak formulation for problem that will be used as a starting point for its discretization.

Strain rate-shear stress law {#sec:strain.rate.shear.stress.law}
----------------------------

We define the Frobenius inner product such that, for all $\b\tau= (\tau_{ij})_{1 \le i,j \le d}$ and $\b\eta= (\eta_{ij})_{1 \le i,j \le d}$ in ${{\mathbb{R}}^{d \times d}}$, $\b\tau : \b\eta \coloneqq \sum_{i,j=1}^d \tau_{ij}\eta_{ij}$, and we denote by $|\b\tau|_{d \times d}\coloneqq \sqrt{\b\tau : \b\tau}$ the corresponding norm.

\[ass:stress\] Let a real number ${r}\in (1,+\infty)$ be fixed, denote by ${r}' \coloneqq \frac{{r}}{{r}-1} \in (1,+\infty)$ the conjugate exponent of ${r}$, and define the singular exponent of ${r}$ by $$\label{eq:sing}
    {r^{\circ}}\coloneq \min({r},2) \in (1,2].$$ The strain rate-shear stress law satisfies

\[eq:ass:sigma\] $$\begin{gathered}
      {\b\sigma}(\b x,\b 0) = \b 0 \text{ for almost every } \b x \in \Omega,\label{eq:ass-stress:0}
      \\
      {\b\sigma}: \Omega \times {{\mathbb{R}}^{d \times d}_{\mathrm{s}}} \to {{\mathbb{R}}^{d \times d}_{\mathrm{s}}} \text{ is measurable}.\label{eq:ass-stress:power-framed} 
    \end{gathered}$$ Moreover, there exist real numbers $\sigma_{\mathrm{de}} \in [0,+\infty)$ and $\sigma_{\mathrm{hc}},\sigma_{\mathrm{sm}} \in (0,+\infty)$ such that, for all $\b\tau,\b\eta \in {{\mathbb{R}}^{d \times d}_{\mathrm{s}}}$ and almost every $\b x \in \Omega$, we have the Hölder continuity property \[eq:power-framed:s.holder.continuity.strong.monotonicity\] $$\begin{aligned}
        \left|
        {\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta)
        \right|_{d\times d} &\le \sigma_{\mathrm{hc}} \left(\sigma_{\mathrm{de}}^{r}+|\b\tau|_{d\times d}^{r}+|\b\eta|_{d\times d}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}| \b\tau-\b\eta |_{d\times d}^{{r^{\circ}}-1},\label{eq:power-framed:s.holder.continuity} 
      \end{aligned}$$ and the strong monotonicity property $$\begin{aligned}
        \left({\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta)\right):(\b\tau-\b\eta) \left(\sigma_{\mathrm{de}}^{r}+|\b\tau|_{d\times d}^{r}+|\b\eta|_{d\times d}^{r}\right)^\frac{2-{r^{\circ}}}{{r}} \ge \sigma_{\mathrm{sm}}|\b\tau-\b\eta|_{d\times d}^{{r}+2-{r^{\circ}}}.\label{eq:power-framed:s.strong.monotonicity}
      \end{aligned}$$

Some remarks are in order.

Assumption can be relaxed by taking ${\b\sigma}(\cdot,\b 0)  \in L^{{r}'}(\Omega,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}})$. This modification requires only minor changes in the analysis, not detailed for the sake of conciseness.

Inequalities – can be proved starting from the following assumptions, which correspond to the conditions characterizing an ${r}$-power-framed function: For all $\b\tau,\b\eta \in {{\mathbb{R}}^{d \times d}_{\mathrm{s}}}$ with $\b\tau \neq \b\eta$ and almost every $\b x \in \Omega$, $$\begin{aligned} 
      |{\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta)|_{d \times d} &\le \sigma_{\mathrm{hc}} \left(\sigma_{\mathrm{de}}^{r}+|\b\tau|_{d \times d}^{r}+|\b\eta|_{d \times d}^{r}\right)^\frac{{r}-2}{{r}}| \b\tau-\b\eta |_{d \times d},
      \\
      \left({\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta)\right):\left(\b\tau-\b\eta\right) &\ge \sigma_{\mathrm{sm}}\left(\sigma_{\mathrm{de}}^{r}+|\b\tau|_{d \times d}^{r}+|\b\eta|_{d \times d}^{r}\right)^\frac{{r}-2}{{r}}|\b\tau-\b\eta|_{d \times d}^{2}.
    \end{aligned}$$ These relations are reminiscent of the ones used in [@Di-Pietro.Droniou:17*1] in the context of scalar Leray–Lions problems. The advantage of assumptions -, expressed in terms of the singular index ${r^{\circ}}$, is that they enable a unified treatment of the cases ${r}< 2$ and ${r}\ge 2$ in the proofs of Lemma \[lem:ah:holder.continuity.strong.monotonicity\], Theorem \[thm:well-posedness\], Lemma \[lem:consistency:ah\], and Theorem \[thm:error.estimate\] below.

Inequalities and give $$\label{eq:power-framed:constants.bound}
    \sigma_{\mathrm{sm}} \leq \sigma_{\mathrm{hc}}.$$ Indeed, let $\b\tau \in {{\mathbb{R}}^{d \times d}_{\mathrm{s}}}$ be such that $|\b\tau|_{d\times d} > 0$. Using the strong monotonicity (with $\b \eta = \b 0$), the Cauchy–Schwarz inequality, and the Hölder continuity (again with $\b \eta = \b 0$), we infer that $$\begin{aligned}
    \sigma_{\mathrm{sm}}\left(\sigma_{\mathrm{de}}^{r}+|\b\tau|_{d\times d}^{r}\right)^\frac{{r^{\circ}}-2}{{r}}|\b\tau|_{d\times d}^{{r}+2-{r^{\circ}}} &\leq  {\b\sigma}(\cdot,\b\tau):\b\tau\leq |{\b\sigma}(\cdot,\b\tau)|_{d\times d}|\b\tau|_{d\times d}\leq \sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{de}}^{r}+|\b\tau|_{d\times d}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}| \b\tau|_{d\times d}^{{r^{\circ}}}
  \end{aligned}$$ almost everywhere in $\Omega$. Hence, $\frac{\sigma_{\mathrm{sm}}}{\sigma_{\mathrm{hc}}} \le \left(\frac{\sigma_{\mathrm{de}}^{r}+|\b\tau|_{d\times d}^{r}}{| \b\tau|_{d\times d}^{r}}\right)^\frac{|{r}-2|}{{r}}$. Letting $|\b\tau|_{d\times d} \to +\infty$ gives .

\[ex:Carreau–Yasuda\] $(\mu,\delta,a,{r})$-Carreau–Yasuda fluids, introduced in [@Yasuda.Armstrong.Cohen:81] and later generalized in [@Hirn:13 Eq. (1.2)], are fluids for which it holds, for almost every $\b x\in\Omega$ and all $\b\tau \in {{\mathbb{R}}^{d \times d}_{\mathrm{s}}}$, $$\label{eq:Carreau--Yasuda}
    {\b\sigma}(\b x,\b\tau) = \mu(\b x)\left(\delta^{a(\b x)}+|\b\tau|_{d \times d}^{a(\b x)}\right)^\frac{{r}-2}{a(\b x)}\b\tau,$$ where $\mu : \Omega \to [\mu_-,\mu_+]$ is a measurable function with $\mu_-,\mu_+ \in (0,+\infty)$ corresponding to the local flow consistency index, $\delta \in [0,+\infty)$ is the degeneracy parameter, $a : \Omega \to [a_-,a_+]$ is a measurable function with $a_-,a_+ \in (0,+\infty)$ expressing the local transition flow behavior index, and ${r}\in (1,+\infty)$ is the flow behavior index. The Carreau–Yasuda law is a generalization of the Carreau law (corresponding to $a_- = a_+ = 2$) that takes into account the different local levels of flow behavior in the fluid. The degenerate case $\delta=0$ corresponds to the power-law model. Non-Newtonian fluids described by constitutive laws with a $(\mu,\delta,a,{r})$-structure exhibit a different behavior according to the value of ${r}$. If ${r}> 2$, then the fluid shows shear thickening behavior and is called *dilatant*. Examples of dilatant fluids are wet sand and oobleck. The case ${r}< 2$, on the other hand, corresponds to *pseudoplastic* fluids having shear thinning behavior, such as blood. Finally, if ${r}= 2$, then the fluid is Newtonian and becomes the classical (linear) Stokes problem. We show in Appendix \[sec:properties.stress\] that the strain rate-shear stress law is an ${r}$-power-framed function with $\sigma_{\mathrm{de}} = \delta$, $$\sigma_{\mathrm{hc}} = \begin{cases}
      \frac{\mu_+}{{r}-1}2^{\left[-\left(\frac{1}{a_+}-\frac{1}{{r}}\right)^\ominus-1\right]({r}-2)+\frac{1}{r}} & \text{if } {r}< 2, 
      \\
      \mu_+({r}-1)2^{\left(\frac{1}{a_-}-\frac{1}{{r}}\right)^\oplus({r}-2)} & \text{if } {r}\ge 2,
    \end{cases}
    \quad\text{and}\quad
    \sigma_{\mathrm{sm}} = \begin{cases}
      \mu_-({r}-1)2^{\left(\frac{1}{a_-}-\frac{1}{{r}}\right)^\oplus({r}-2)} & \text{if } {r}\le 2, 
      \\
      \frac{\mu_-}{{r}-1}2^{\left[-\left(\frac{1}{a_+}-\frac{1}{{r}}\right)^\ominus-1\right]({r}-2)-1} & \text{if } {r}> 2,
    \end{cases}$$ where $\xi^\oplus\coloneq\max(0,\xi)$ and $\xi^\ominus\coloneq-\min(0,\xi)$ denote, respectively, the positive and negative parts of a real number $\xi$. As a consequence, it matches Assumption \[ass:stress\].

Weak formulation {#sec:weak.formulation}
----------------

From this point on, we omit both the integration variable and the measure from integrals, as they can be in all cases inferred from the context. We define the following velocity and pressure spaces embedding, respectively, the homogeneous boundary condition for the velocity and the zero-average constraint for the pressure: $$\b U \coloneqq \left\{\b v \in W^{1,{r}}(\Omega,{\mathbb{R}}^d)\ : \ \b v{\ \!\!_{|_{\partial\Omega}}} = \b 0 \right\},
\qquad
P \coloneqq L^{{r}'}_0(\Omega,{\mathbb{R}}) \coloneqq \left\{q \in L^{{r}'}(\Omega,{\mathbb{R}})\ : \ \textstyle\int_\Omega q = 0 \right\}.$$ Assuming $\b f \in L^{{r}'}(\Omega,{\mathbb{R}}^d)$, the weak formulation of problem reads: Find $(\b u,p) \in \b U \times P$ such that

\[eq:stokes.weak\] $$\begin{aligned}
{2}
     a(\b u,\b v)+b(\b v,p) &= \displaystyle\int_\Omega \b f \cdot \b v &\qquad \forall \b v \in \b U,\label{eq:stokes.weak:momentum} \\
     -b(\b u,q) &= 0 &\qquad \forall q \in P,  \label{eq:stokes.weak:mass} 
  \end{aligned}$$

where the function $a : \b U \times \b U \to {\mathbb{R}}$ and the bilinear form $b : \b U \times L^{{r}'}(\Omega,{\mathbb{R}}) \to {\mathbb{R}}$ are defined such that, for all $\b v,\b w \in \b U$ and all $q \in L^{{r}'}(\Omega,{\mathbb{R}})$, $$\label{eq:a.b}
  a(\b w,\b v) \coloneqq \displaystyle\int_\Omega {\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w) : {\b{\nabla}_{\mathrm{s}}}\b v,\qquad
  b(\b v,q)  \coloneqq -\displaystyle\int_\Omega ({\b\nabla{\cdot}}\b v) q.$$

The test space in can be extended to $L^{{r}'}(\Omega,{\mathbb{R}})$ since, for all $\b v \in \b U$, the divergence theorem and the fact that $\b v{\ \!\!_{|_{\partial\Omega}}} = \b 0$ yield $b(\b v,1) = -\int_\Omega {\b\nabla{\cdot}}\b v = - \int_{\partial\Omega} \b v \cdot \b n_{\partial \Omega} = 0$, with $\b n_{\partial \Omega}$ denoting the unit vector normal to $\partial\Omega$ and pointing out of $\Omega$.

\[rem:a-priori\] It can be checked that, under Assumption \[ass:stress\], the continuous problem admits a unique solution $(\b u,p) \in \b U \times P$; see, e.g., [@Hirn:13 Section 2.4], where slightly stronger assumptions are considered. For future use, we also note the following a priori bound on the velocity: $$\label{eq:continuous.solution:bounds:uh}
    |\b u|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)} \lesssim \left(\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\right)^\frac{1}{{r}-1}+\left(\sigma_{\mathrm{de}}^{2-{r^{\circ}}}\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\right)^\frac{1}{{r}+1-{r^{\circ}}}.$$ To prove , use the strong-monotonicity of ${\b\sigma}$, sum written for $\b v=\b u$ to written for $q=p$, and use the Hölder and Korn inequalities to write $$\begin{aligned}
    \sigma_{\mathrm{sm}}\left(
    \sigma_{\mathrm{de}}^{r}+ \|{\b{\nabla}_{\mathrm{s}}}\b u\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}\right)^\frac{{r^{\circ}}-2}{{r}} \|{\b{\nabla}_{\mathrm{s}}}\b u\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{{r}+2-{r^{\circ}}}
    &\lesssim  a(\b u,\b u)
    = \displaystyle\int_\Omega \b f \cdot \b u
    \lesssim \| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\|{\b{\nabla}_{\mathrm{s}}}\b u\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})},
  \end{aligned}$$ that is, $$\label{eq:continuous.solution:bounds:uh:1}
    \mathcal{N}\coloneqq \left(
    \sigma_{\mathrm{de}}^{r}+ \|{\b{\nabla}_{\mathrm{s}}}\b u\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}\right)^\frac{{r^{\circ}}-2}{{r}} \|{\b{\nabla}_{\mathrm{s}}}\b u\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{{r}+1-{r^{\circ}}} \lesssim \sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}.$$ Observing that $\|{\b{\nabla}_{\mathrm{s}}}\b u\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{{r}+1-{r^{\circ}}} \lesssim \max\left(\|{\b{\nabla}_{\mathrm{s}}}\b u\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})},\sigma_{\mathrm{de}}\right)^{2-{r^{\circ}}}\mathcal{N} $, we obtain, enumerating the cases for the maximum and summing the corresponding bounds, $\|{\b{\nabla}_{\mathrm{s}}}\b u\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})} \lesssim \mathcal{N}^\frac{1}{{r}-1}+(\sigma_{\mathrm{de}}^{2-{r^{\circ}}}\mathcal{N})^\frac{1}{{r}+1-{r^{\circ}}} $. Combining this inequality with gives .

Discrete setting {#sec:discrete.setting}
================

In this section, we recall the notion of polyhedral mesh along with the definitions and properties of $L^2$-orthogonal projectors on local and broken polynomial spaces. Then, after introducing the spaces of discrete unknowns for the velocity and the pressure, we prove a discrete Korn inequality, define the discrete counterparts of the function $a$ and of the bilinear form $b$, and formulate the HHO scheme.

Mesh and notation for inequalities up to a multiplicative constant
------------------------------------------------------------------

We define a mesh as a couple $\mathcal M_h\coloneq({\mathcal{T}}_h,{\mathcal{F}}_h)$, where ${\mathcal{T}}_h$ is a finite collection of polyhedral elements $T$ such that $h=\max_{T\in{\mathcal{T}}_h}h_T$ with $h_T$ denoting the diameter of $T$, while ${\mathcal{F}}_h$ is a finite collection of planar faces $F$ with diameter $h_F$. Notice that, here and in what follows, we use the three-dimensional nomenclature also when $d=2$, i.e., we speak of polyhedra and faces rather than polygons and edges. It is assumed henceforth that the mesh $\mathcal M_h$ matches the geometrical requirements detailed in [@Di-Pietro.Droniou:20 Definition 1.7]. In order to have the boundedness property for the interpolator, we additionally assume that the mesh elements are star-shaped with respect to every point of a ball of radius uniformly comparable to the element diameter; see [@Di-Pietro.Droniou:20 Lemma 7.12] for the Hilbertian case. Boundary faces lying on $\partial\Omega$ and internal faces contained in $\Omega$ are collected in the sets ${\mathcal{F}}_h^{\rm b}$ and ${\mathcal{F}}_h^{\rm i}$, respectively. For every mesh element $T\in{\mathcal{T}}_h$, we denote by ${\mathcal{F}}_T$ the subset of ${\mathcal{F}}_h$ containing the faces that lie on the boundary $\partial T$ of $T$. For every face $F \in {\mathcal{F}}_h$, we denote by ${\mathcal{T}}_F$ the subset of ${\mathcal{T}}_h$ containing the one (if $F\in{\mathcal{F}}_h^{\rm b}$) or two (if $F\in{\mathcal{F}}_h^{\rm i}$) elements on whose boundary $F$ lies. For each mesh element $T\in{\mathcal{T}}_h$ and face $F\in{\mathcal{F}}_T$, $\b n_{TF}$ denotes the (constant) unit vector normal to $F$ pointing out of $T$.

Our focus is on the $h$-convergence analysis, so we consider a sequence of refined meshes that is regular in the sense of [@Di-Pietro.Droniou:20 Definition 1.9] with regularity parameter uniformly bounded away from zero. The mesh regularity assumption implies, in particular, that the diameter of a mesh element and those of its faces are comparable uniformly in $h$ and that the number of faces of one element is bounded above by an integer independent of $h$.

To avoid the proliferation of generic constants, we write henceforth $a\lesssim b$ (resp., $a\gtrsim b$) for the inequality $a\le Cb$ (resp., $a\ge Cb$) with real number $C>0$ independent of $h$, of the constants $\sigma_{\mathrm{de}},\sigma_{\mathrm{hc}},\sigma_{\mathrm{sm}}$ in Assumption \[ass:stress\], and, for local inequalities, of the mesh element or face on which the inequality holds. We also write $a\simeq b$ to mean $a\lesssim b$ and $b\lesssim a$. The dependencies of the hidden constants are further specified when needed.

Projectors and broken spaces
----------------------------

Given $X \in {\mathcal{T}}_h \cup {\mathcal{F}}_h$ and $l \in {\mathbb{N}}$, we denote by ${\mathbb{P}}^l(X,{\mathbb{R}})$ the space spanned by the restriction to $X$ of scalar-valued, $d$-variate polynomials of total degree $\le l$. The local $L^2$-orthogonal projector ${\pi_{X}^{l}} : L^{1}(X,{\mathbb{R}}) \to {\mathbb{P}}^l(X,{\mathbb{R}})$ is defined such that, for all $v \in L^{1}(X,{\mathbb{R}})$, $$\label{eq:proj}
  \displaystyle\int_X ({\pi_{X}^{l}} v-v) w = 0 \qquad \forall w \in  {\mathbb{P}}^{l}(X,{\mathbb{R}}).$$ When applied to vector-valued fields in $L^1(X,{\mathbb{R}}^d)$ (resp., tensor-valued fields in $L^1(X,{{\mathbb{R}}^{d \times d}})$), the $L^2$-orthogonal projector mapping on ${\mathbb{P}}^l(X,{\mathbb{R}}^d)$ (resp., ${\mathbb{P}}^l(X,{{\mathbb{R}}^{d \times d}})$) acts component-wise and is denoted in boldface font. Let $T\in{\mathcal{T}}_h$, $n\in[0,l+1]$ and $m\in[0,n]$. The following $(n,{r},m)$-approximation properties of ${\pi_{T}^{l}}$ hold: For any $v\in W^{n,{r}}(T,{\mathbb{R}})$,

\[eq:proj:app\] $$\label{eq:proj:app:T}
  |v-{\pi_{T}^{l}}v|_{W^{m,{r}}(T,{\mathbb{R}})} \lesssim h_T^{n-m}|v|_{W^{n,{r}}(T,{\mathbb{R}})}.$$ The above property will also be used in what follows with ${r}$ replaced by its conjugate exponent ${r}'$. If, additionally, $n\ge 1$, we have the following $(n,{r}')$-trace approximation property: $$\label{eq:proj:app:F}
    \|v-{\pi_{T}^{l}}v\|_{L^{{r}'}(\partial T,{\mathbb{R}})}\lesssim h_T^{n-\frac{1}{{r}'}}|v|_{W^{n,{r}'}(T,{\mathbb{R}})}.$$

The hidden constants in are independent of $h$ and $T$, but possibly depend on $d$, the mesh regularity parameter, $l$, $n$, and ${r}$. The approximation properties are proved for integer $n$ and $m$ in [@Di-Pietro.Droniou:17 Appendix A.2] (see also [@Di-Pietro.Droniou:20 Theorem 1.45]), and can be extended to non-integer vales using standard interpolation techniques (see, e.g., [@Lions.Magenes:72 Theorem 5.1]).

At the global level, for a given integer $l\ge 0$, we define the broken polynomial space ${\mathbb{P}}^l({\mathcal{T}}_h,{\mathbb{R}})$ spanned by functions in $L^1(\Omega,{\mathbb{R}})$ whose restriction to each mesh element $T\in{\mathcal{T}}_h$ lies in ${\mathbb{P}}^l(T,{\mathbb{R}})$, and we define the global $L^2$-orthogonal projector ${\pi_{h}^{l}} : L^{1}(\Omega,{\mathbb{R}}) \to {\mathbb{P}}^l({\mathcal{T}}_h,{\mathbb{R}})$ such that, for all $v \in L^{1}(\Omega,{\mathbb{R}})$ and all $T \in {\mathcal{T}}_h$, $$({\pi_{h}^{l}} v){\ \!\!_{|_{T}}} \coloneq {\pi_{T}^{l}} v{\ \!\!_{|_{T}}}.$$ Broken polynomial spaces are subspaces of the broken Sobolev spaces $$W^{n,{r}}({\mathcal{T}}_h,{\mathbb{R}})\coloneq\left\{ v\in L^{r}(\Omega,{\mathbb{R}})\ : \ v{\ \!\!_{|_{T}}}\in W^{n,{r}}(T,{\mathbb{R}})\quad\forall T\in{\mathcal{T}}_h\right\}.$$ We define the broken gradient operator ${\b\nabla_h}: W^{1,1}({\mathcal{T}}_h,{\mathbb{R}}) \rightarrow L^1(\Omega,{\mathbb{R}}^d)$ such that, for all $v \in W^{1,1}({\mathcal{T}}_h,{\mathbb{R}})$ and all $T \in {\mathcal{T}}_h$, $({\b\nabla_h}v){\ \!\!_{|_{T}}} \coloneq {\b\nabla}v{\ \!\!_{|_{T}}}$. We define similarly the broken gradient acting on vector fields along with its symmetric part ${\b{\nabla}_{{\mathrm{s}},h}}$, as well as the broken divergence operator ${\b\nabla_h \cdot}$ acting on tensor fields. The global $L^2$-orthogonal projector ${\b{\pi}_{h}^{l}}$ mapping vector-valued fields in $L^1(\Omega,{\mathbb{R}}^d)$ (resp., tensor-valued fields in $L^1(\Omega,{{\mathbb{R}}^{d \times d}})$) on ${\mathbb{P}}^l({\mathcal{T}}_h,{\mathbb{R}}^d)$ (resp., ${\mathbb{P}}^l({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}})$) is obtained applying ${\pi_{h}^{l}}$ component-wise.

Discrete spaces and norms
-------------------------

Let an integer $k\ge 1$ be fixed. The HHO space of discrete velocity unknowns is $${{{\und{\b{U}}}}_{h}^{k}} \coloneqq \left\{
{\und{\b{v}}}_h \coloneq ((\b v_T)_{T \in {\mathcal{T}}_h},(\b v_F)_{F\in {\mathcal{F}}_h}) \ : \ \b v_T \in {\mathbb{P}}^k(T,{\mathbb{R}}^d)\ \ \forall T \in {\mathcal{T}}_h\ \mbox{ and }\ \b v_F \in {\mathbb{P}}^k(F,{\mathbb{R}}^d)\ \ \forall F \in {\mathcal{F}}_h \right\}.$$ The interpolation operator ${{\und{\b{I}}}_{h}^{k}} : W^{1,1}(\Omega,{\mathbb{R}}^d) \to  {{{\und{\b{U}}}}_{h}^{k}} $ maps a function $\b v \in W^{1,1}(\Omega,{\mathbb{R}}^d)$ on the vector of discrete unknowns ${{\und{\b{I}}}_{h}^{k}}\b v$ defined as follows: $${{\und{\b{I}}}_{h}^{k}} \b v \coloneqq (({\b{\pi}_{T}^{k}} \b v{\ \!\!_{|_{T}}})_{T \in {\mathcal{T}}_h},({\b{\pi}_{F}^{k}} \b v{\ \!\!_{|_{F}}})_{F \in {\mathcal{F}}_h}).$$ For all $T \in {\mathcal{T}}_h$, we denote by ${{{\und{\b{U}}}}_{T}^{k}}$ and ${{\und{\b{I}}}_{T}^{k}}$ the restrictions of ${{\und{\b{I}}}_{h}^{k}}$ and ${{{\und{\b{U}}}}_{h}^{k}}$ to $T$, respectively and, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$, we let ${\und{\b{v}}}_T \coloneqq (\b v_T,(\b v_F)_{F\in {\mathcal{F}}_T}) \in {{{\und{\b{U}}}}_{T}^{k}}$ denote the vector collecting the discrete unknowns attached to $T$ and its faces. Furthermore, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$, we define the broken polynomial field $\b v_h\in{\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}}^d)$ obtained patching element unknowns, that is, $$(\b v_h){\ \!\!_{|_{T}}} \coloneqq \b v_T\qquad\forall T \in {\mathcal{T}}_h.$$

We define on ${{{\und{\b{U}}}}_{h}^{k}}$ the $W^{1,{r}}(\Omega,{\mathbb{R}}^d)$-like seminorm $\| {\cdot} \|_{{\boldsymbol{\varepsilon}},{r},h}$ such that, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$,

\[eq:norm.epsilon.r\] $$\begin{gathered}
\label{eq:norm.epsilon.r.h}
    \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} \coloneqq \left(\displaystyle\sum_{T \in {\mathcal{T}}_h}\| {\und{\b{v}}}_T \|_{{\boldsymbol{\varepsilon}},{r},T}^{r}\right)^\frac{1}{{r}}
    \\\label{eq:norm.epsilon.r.T}
    \text{with 
      $\| {\und{\b{v}}}_T \|_{{\boldsymbol{\varepsilon}},{r},T} \coloneqq \left(\| {\b{\nabla}_{\mathrm{s}}}\b v_T \|^{r}_{L^{r}(T,{{\mathbb{R}}^{d \times d}})} + \displaystyle\sum_{F \in {\mathcal{F}}_T} h_F^{1-{r}} \| \b v_F - \b v_T\|^{r}_{L^{r}(F,{\mathbb{R}}^d)}\right)^\frac{1}{{r}}$
      for all $T \in {\mathcal{T}}_h$.
    }
  \end{gathered}$$

The following boundedness property for ${{\und{\b{I}}}_{T}^{k}}$ can be proved adapting the arguments of [@Di-Pietro.Droniou:20 Proposition 6.24] and requires the star-shaped assumption on the mesh elements: For all $T \in {\mathcal{T}}_h$ and all $\b v \in W^{1,{r}}(T,{\mathbb{R}}^d)$, $$\label{eq:I:boundedness}
  \|{{\und{\b{I}}}_{T}^{k}} \b v\|_{{\boldsymbol{\varepsilon}},{r},T} \lesssim | \b v |_{W^{1,{r}}(T,{\mathbb{R}}^d)},$$ where the hidden constant depends only on $d$, the mesh regularity parameter, ${r}$, and $k$.

The discrete velocity and pressure are sought in the following spaces, which embed, respectively, the homogeneous boundary condition for the velocity and the zero-average constraint for the pressure: $${{{\und{\b{U}}}}_{h,0}^{k}} \coloneqq \left\{ {\und{\b{v}}}_h = ((\b v_T)_{T \in {\mathcal{T}}_h},(\b v_F)_{F\in {\mathcal{F}}_h}) \in {{{\und{\b{U}}}}_{h}^{k}} \ : \ \b v_F = \b 0 \quad \forall F \in {{\mathcal{F}}_h^{\mathrm{b}}}\right\},\quad
{P_{h}^{k}} \coloneqq {\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}})\cap P.$$ By the discrete Korn inequality proved in Lemma \[lem:discrete.korn.inequality\] below, $\| {\cdot} \|_{{\boldsymbol{\varepsilon}},{r},h}$ is a norm on ${{{\und{\b{U}}}}_{h,0}^{k}}$ (the proof is obtained reasoning as in [@Di-Pietro.Droniou:20 Corollary 2.16]).

Discrete Korn inequality {#sec:discrete.korn.inequality}
------------------------

We prove in this section a discrete counterpart of the following Korn inequality (see [@Geymonat.Suquet:86 Theorem 1]): For all $\b v \in\b U$. $$\label{eq:Korn}
\|\b v\|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)} \lesssim \|{\b{\nabla}_{\mathrm{s}}}\b v\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}.$$ We start by recalling a few preliminary results. The first concerns inequalities between sums of powers, and will be often used in what follows without necessarily recalling this fact explicitly each time. Let an integer $n\ge 1$ and a real number $m \in (0,+\infty)$ be given. Then, for all $a_1,\ldots,a_n \in (0,+\infty) $, we have $$\label{eq:sum-power}
  n^{-(m-1)^\ominus}\sum_{i=1}^n a_i^m \le \left(\sum_{i=1}^n a_i\right)^m  \le n^{(m-1)^\oplus}\sum_{i=1}^n a_i^m.$$ If $m=1$, then holds with the equal sign. If $m < 1$, [@Ursell:59 Eqs. (5) and (3)] with $\alpha = 1$ and $\beta = m$ give $
n^{m-1}\sum_{i=1}^n a_i^m \le \left(\sum_{i=1}^n a_i\right)^m  \le \sum_{i=1}^n a_i^m.
$ If, on the other hand, $m > 1$, [@Ursell:59 Eqs. (3) and (5)] with $\alpha = m$ and $\beta = 1$ give $
\sum_{i=1}^n a_i^m \le \left(\sum_{i=1}^n a_i\right)^m  \le n^{m-1}\sum_{i=1}^n a_i^m.
$ Gathering the above cases yields .

The second preliminary result concerns the node-averaging interpolator. Let $\mathfrak T_h$ be a matching simplicial submesh of $\mathcal M_h$ in the sense of [@Di-Pietro.Droniou:20 Definition 1.8]. The node-averaging operator ${\b{I}_{{\mathrm{av}},h}^{k}} : {\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}}^d) \to {\mathbb{P}}^k(\mathfrak T_h,{\mathbb{R}}^d) \cap W^{1,{r}}(\Omega,{\mathbb{R}}^d)$ is such that, for all $\b v_h \in {\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}}^d)$ and all Lagrange node $V$ of $\mathfrak T_h$, denoting by $\mathfrak T_V$ the set of simplices sharing $V$, $$({\b{I}_{{\mathrm{av}},h}^{k}}\b v_h)(V) \coloneq
\begin{cases}
  \frac{1}{\rm{card}(\mathfrak T_V)}\sum_{\b\tau\in\mathfrak T_V} \b v_h{\ \!\!_{|_{\b\tau}}}(V) & \text{if }\ V \in \Omega,
  \\ \b 0 & \text{if }\ V \in \partial\Omega.
\end{cases}$$ For all $F \in {{\mathcal{F}}_h^{\mathrm{i}}}$, denote by $T_1,T_2\in{\mathcal{T}}_h$ the elements sharing $F$, taken in an arbitrary but fixed order. We define the jump operator such that, for any function $\b v\in W^{1,1}({\mathcal{T}}_h,{\mathbb{R}}^d)$, $[\b v]_F \coloneq (\b v{\ \!\!_{|_{T_1}}}){\ \!\!_{|_{F}}}-(\b v{\ \!\!_{|_{T_2}}}){\ \!\!_{|_{F}}}$. This definition is extended to boundary faces $F\in{{\mathcal{F}}_h^{\mathrm{b}}}$ by setting $[\b v]_F \coloneq \b v{\ \!\!_{|_{F}}}$.

For all $\b v_h \in {\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}}^d)$, it holds $$\label{eq:Iav:bound}
    |\b v_h-{\b{I}_{{\mathrm{av}},h}^{k}}\b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}\lesssim \sum_{F\in{\mathcal{F}}_h} h_F^{1-{r}} \|[\b v_h]_F\|^{r}_{L^{r}(F,{\mathbb{R}}^d)}.$$

Combining [@Di-Pietro.Droniou:20 Eq. (4.13)] (which corresponds to for ${r}=2$) with the local Lebesgue embeddings of [@Di-Pietro.Droniou:20 Lemma 1.25] (see also [@Di-Pietro.Droniou:17 Lemma 5.1]) gives, for any $T\in{\mathcal{T}}_h$, $$\label{eq:Iav:bound:0}
    \|\b v_h-{\b{I}_{{\mathrm{av}},h}^{k}}\b v_h\|_{L^{r}(T,{\mathbb{R}}^d)}^{r}\lesssim \sum_{F\in{\mathcal{F}}_{\mathcal V,T}}h_F \|[\b v_h]_F\|^{r}_{L^{r}(F,{\mathbb{R}}^d)},$$ where ${\mathcal{F}}_{\mathcal V,T}$ collects the faces whose closure has non-empty intersection with $\overline{T}$. Using the local inverse inequality of [@Di-Pietro.Droniou:20 Lemma 1.28] (see also [@Di-Pietro.Droniou:17 Eq. (A.1)]) we can write $$\begin{aligned}
    |\b v_h-{\b{I}_{{\mathrm{av}},h}^{k}}\b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}&\lesssim \sum_{T\in{\mathcal{T}}_h} h_T^{-{r}}\|\b v_h-{\b{I}_{{\mathrm{av}},h}^{k}}\b v_h\|_{L^{r}(T,{\mathbb{R}}^d)}^{r}\\
    &\lesssim \sum_{T\in{\mathcal{T}}_h}\sum_{F\in{\mathcal{F}}_{\mathcal V,T}} h_F^{1-{r}} \|[\b v_h]_F\|^{r}_{L^{r}(F,{\mathbb{R}}^d)}
    \\
    &\lesssim \sum_{F\in{\mathcal{F}}_h} \sum_{T\in{\mathcal{T}}_{\mathcal V,F}} h_F^{1-{r}} \|[\b v_h]_F\|^{r}_{L^{r}(F,{\mathbb{R}}^d)}
    \\
    &\le\max_{F\in{\mathcal{F}}_h}\operatorname{card}({\mathcal{T}}_{\mathcal V,F}) \sum_{F\in{\mathcal{F}}_h} h_F^{1-{r}} \|[\b v_h]_F\|^{r}_{L^{r}(F,{\mathbb{R}}^d)},
  \end{aligned}$$ where we have used the fact that $h_T^{-{r}} \le h_F^{-{r}}$ along with inequality to pass to the second line, and we have exchanged the sums after setting ${\mathcal{T}}_{\mathcal V,F} \coloneqq \big\{T \in {\mathcal{T}}_h : \overline F \cap \overline T \neq \emptyset\big\}$ for all $F \in {\mathcal{F}}_h$ to pass to the third line. Observing that $\max_{F\in{\mathcal{F}}_h}\operatorname{card}({\mathcal{T}}_{\mathcal V,F}) \lesssim 1$ (since, for any $F\in{\mathcal{F}}_h$, $\operatorname{card}({\mathcal{T}}_{\mathcal V,F})$ is bounded by the left-hand side of [@Di-Pietro.Droniou:20 Eq. (4.23)] written for any $T\in{\mathcal{T}}_h$ to which $F$ belongs), follows.

(Discrete Korn inequality)\[lem:discrete.korn.inequality\] We have, for all ${\und{\b{v}}}_h  \in {{{\und{\b{U}}}}_{h,0}^{k}}$, $$\label{eq:discrete.Korn}
    \|  \b v_h \|_{L^{r}(\Omega,{\mathbb{R}}^d)}^{r}+|\b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}\lesssim \|   {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}.$$

Let ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}$. Using a triangle inequality followed by , we can write $$\begin{aligned}
    |\b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}&\lesssim |{\b{I}_{{\mathrm{av}},h}^{k}} \b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}+|\b v_h-{\b{I}_{{\mathrm{av}},h}^{k}}\b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}\\
    &\lesssim \|{\b{\nabla}_{\mathrm{s}}}({\b{I}_{{\mathrm{av}},h}^{k}} \b v_h) \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}+|\b v_h-{\b{I}_{{\mathrm{av}},h}^{k}}\b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}\\
    &\lesssim \|{\b{\nabla}_{{\mathrm{s}},h}}\b v_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}+|\b v_h-{\b{I}_{{\mathrm{av}},h}^{k}}\b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}\\
    &\lesssim \|{\b{\nabla}_{{\mathrm{s}},h}}\b v_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}+\sum_{F\in{\mathcal{F}}_h} h_F^{1-{r}} \|[\b v_h]_F\|^{r}_{L^{r}(F,{\mathbb{R}}^d)},
  \end{aligned}$$ where we have used the continuous Korn inequality to pass to the second line, we have inserted $\pm{\b{\nabla}_{{\mathrm{s}},h}}\b v_h$ into the first norm and used a triangle inequality followed by to pass to the third line, and we have invoked the bound to conclude. Observing that, for any $F\in{\mathcal{F}}_h$, $|[\b v_h]_F| \leq \sum_{T\in{\mathcal{T}}_F}|\b v_F-\b v_T|$ by a triangle inequality, and using , we can continue writing $$|\b v_h|_{W^{1,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{r}\lesssim \|{\b{\nabla}_{{\mathrm{s}},h}}\b v_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}+\sum_{F\in{\mathcal{F}}_h}\sum_{T\in{\mathcal{T}}_F}h_F^{1-{r}} \| \b v_F - \b v_T\|^{r}_{L^{r}(F,{\mathbb{R}}^d)}
  = \|{\und{\b{v}}}_h\|_{{\boldsymbol{\varepsilon}},{r},h}^{r},$$ where we have exchanged the sums over faces and elements and recalled definition to conclude. This proves the bound for the second term in the left-hand side of . Combining this result with the global discrete Sobolev embeddings of [@Di-Pietro.Droniou:17 Proposition 5.4] yields the bound for the first term in .

Viscous term
------------

### Local symmetric gradient reconstruction

For all $T \in {\mathcal{T}}_h$, we define the local symmetric gradient reconstruction ${\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},T}} : {{{\und{\b{U}}}}_{T}^{k}} \to {\mathbb{P}}^{k}(T,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}})$ such that, for all ${\und{\b{v}}}_T \in {{{\und{\b{U}}}}_{T}^{k}}$, $$\label{eq:G}
  \displaystyle\int_T {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},T}} {\und{\b{v}}}_T : \b\tau = \int_T {\b{\nabla}_{\mathrm{s}}}\b v_T : \b\tau + \sum_{F \in {\mathcal{F}}_T} \int_F (\b v_F-\b v_T)\cdot (\b\tau \b n_{TF})\qquad \forall \b\tau \in  {\mathbb{P}}^{k}(T,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}}).$$ This symmetric gradient reconstruction, originally introduced in [@Botti.Di-Pietro.ea:17 Section 4.2], is designed so that the following relation holds (see, e.g., [@Botti.Di-Pietro.ea:17*1 Proposition 5] or [@Di-Pietro.Droniou:20 Section 7.2.5]): For all $\b v\in W^{1,1}(T,{\mathbb{R}}^d)$, $$\label{eq:G:proj}
  {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},T}} ({{\und{\b{I}}}_{T}^{k}} \b v) = {\b{\pi}_{T}^{k}}({\b{\nabla}_{\mathrm{s}}}\b v).$$ The global symmetric gradient reconstruction ${\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} : {{{\und{\b{U}}}}_{h}^{k}} \to {\mathbb{P}}^{k}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}})$ is obtained patching the local contributions, that is, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$, $$\label{eq:Gh}
  ({\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{v}}}_h){\ \!\!_{|_{T}}} \coloneq {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},T}} {\und{\b{v}}}_T \qquad \forall T\in{\mathcal{T}}_h.$$

### Discrete viscous function

The discrete counterpart of the function $a$ defined by is the function ${\mathrm{a}}_h : {{{\und{\b{U}}}}_{h}^{k}} \times {{{\und{\b{U}}}}_{h}^{k}} \to {\mathbb{R}}$ such that, for all ${\und{\b{v}}}_h,{\und{\b{w}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$, $$\label{eq:ah}
  {\mathrm{a}}_h({\und{\b{w}}}_h, {\und{\b{v}}}_h) \coloneqq \displaystyle\int_\Omega {\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{w}}}_h): {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{v}}}_h +\gamma {\mathrm{s}}_h({\und{\b{w}}}_h,{\und{\b{v}}}_h).$$ In the above definition, recalling , $\gamma$ is a stabilization parameter such that $$\label{eq:gamma}
  \gamma \in [\sigma_{\mathrm{sm}},\sigma_{\mathrm{hc}}],$$ while the stabilization function ${\mathrm{s}}_h : {{{\und{\b{U}}}}_{h}^{k}} \times {{{\und{\b{U}}}}_{h}^{k}} \to {\mathbb{R}}$ is such that, for all ${\und{\b{v}}}_h,{\und{\b{w}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$, $$\label{eq:sh}
  {\mathrm{s}}_h({\und{\b{w}}}_h,{\und{\b{v}}}_h) \coloneqq \displaystyle\sum_{T \in {\mathcal{T}}_h}{\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{v}}}_T),$$ where the local contributions are assumed to satisfy the following assumption.

\[ass:sT\] For all $T \in {\mathcal{T}}_h$, the local stabilization function ${\mathrm{s}}_T:{{{\und{\b{U}}}}_{T}^{k}}\times{{{\und{\b{U}}}}_{T}^{k}}\to{\mathbb{R}}$ is linear in its second argument and satisfies the following properties, with hidden constants independent of both $h$ and $T$:

\[eq:ass:sT\]

1.  \[ass:sT:stability.boundedness\] *Stability and boundedness.* Recalling the definition of the local $\|{\cdot}\|_{{\boldsymbol{\varepsilon}},{r},T}$-seminorm, for all ${\und{\b{v}}}_T\in{{{\und{\b{U}}}}_{T}^{k}}$ it holds: $$\label{eq:sT:stability.boundedness}
          \| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},T}}{\und{\b{v}}}_T \|_{L^{r}(T,{{\mathbb{R}}^{d \times d}})}^{r}+ {\mathrm{s}}_T({\und{\b{v}}}_T,{\und{\b{v}}}_T)
          \simeq \|{\und{\b{v}}}_T\|_{{\boldsymbol{\varepsilon}},{r},T}^{r}.$$

2.  *Polynomial consistency.*\[ass:sT:polynomial.consistency\] For all $\b w\in{\mathbb{P}}^{k+1}(T,{\mathbb{R}}^d)$ and all ${\und{\b{v}}}_T\in{{{\und{\b{U}}}}_{T}^{k}}$, $$\label{eq:sT:polynomial.consistency}
          {\mathrm{s}}_T({{\und{\b{I}}}_{T}^{k}} \b w,{\und{\b{v}}}_T) = 0.$$

3.  *Hölder continuity.*\[ass:sT:holder-continuity\] For all ${\und{\b{u}}}_T, {\und{\b{v}}}_T, {\und{\b{w}}}_T \in {{{\und{\b{U}}}}_{T}^{k}}$, it holds, setting ${\und{\b{e}}}_T\coloneq{\und{\b{u}}}_T - {\und{\b{w}}}_T$, $$\label{eq:sT:holder-continuity}   
         \hspace{-0.5cm} \left|{\mathrm{s}}_T({\und{\b{u}}}_T,{\und{\b{v}}}_T)-{\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{v}}}_T)\right|
          \lesssim 
          \left({\mathrm{s}}_T({\und{\b{u}}}_T,{\und{\b{u}}}_T)+{\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{w}}}_T)\right)^\frac{{r}-{r^{\circ}}}{{r}}{\mathrm{s}}_T({\und{\b{e}}}_T,{\und{\b{e}}}_T)^\frac{{r^{\circ}}-1}{{r}}{\mathrm{s}}_T({\und{\b{v}}}_T,{\und{\b{v}}}_T)^\frac{1}{{r}}.$$

4.  *Strong monotonicity.*\[ass:sT:strong-monotonicity\] For all ${\und{\b{u}}}_T, {\und{\b{w}}}_T \in {{{\und{\b{U}}}}_{T}^{k}}$ , it holds, setting again ${\und{\b{e}}}_T\coloneq{\und{\b{u}}}_T - {\und{\b{w}}}_T$, $$\label{eq:sT:strong-monotonicity} 
          \left({\mathrm{s}}_T({\und{\b{u}}}_T,{\und{\b{e}}}_T) - {\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{e}}}_T)\right)\left( {\mathrm{s}}_T({\und{\b{u}}}_T,{\und{\b{u}}}_T)+{\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{w}}}_T)\right)^\frac{2-{r^{\circ}}}{{r}} \gtrsim {\mathrm{s}}_T({\und{\b{e}}}_T,{\und{\b{e}}}_T)^\frac{{r}+2-{r^{\circ}}}{{r}}.$$

If ${r}=2$, ${\mathrm{s}}_T$ can be any symmetric bilinear form satisfying [[(S\[ass:sT:stability.boundedness\])]{}]{}–[[(S\[ass:sT:polynomial.consistency\])]{}]{}. Indeed, property [[(S\[ass:sT:holder-continuity\])]{}]{} coincides in this case with the Cauchy–Schwarz inequality, while, by linearity of ${\mathrm{s}}_T$, property [[(S\[ass:sT:strong-monotonicity\])]{}]{} holds with the equal sign.

\[lem:sT:consist\] For any $T\in{\mathcal{T}}_h$ and any ${\mathrm{s}}_T$ satisfying Assumption \[ass:sT\], it holds, for all $\b w \in W^{k+2,{r}}(T,{\mathbb{R}}^d)$ and all ${\und{\b{v}}}_T \in {{{\und{\b{U}}}}_{T}^{k}}$, $$\label{eq:sT:consist}
    |{\mathrm{s}}_T({{\und{\b{I}}}_{T}^{k}} \b w,{\und{\b{v}}}_T)| \lesssim h_T^{(k+1)({r^{\circ}}-1)}| \b w |_{W^{1,{r}}(T,{\mathbb{R}}^d)}^{{r}-{r^{\circ}}}|\b w|_{W^{k+2,{r}}(T,{\mathbb{R}}^d)}^{{r^{\circ}}-1}\| {\und{\b{v}}}_T \|_{{\boldsymbol{\varepsilon}},{r},T},$$ where the hidden constant is independent of $h$, $T$, and $\b w$.

The proof adapts the arguments of [@Di-Pietro.Droniou:20 Propositon 2.14]. Using the polynomial consistency property [[(S\[ass:sT:polynomial.consistency\])]{}]{}, we can write $$\begin{aligned}
      |{\mathrm{s}}_T({{\und{\b{I}}}_{T}^{k}} \b w,{\und{\b{v}}}_T)|
      &= |{\mathrm{s}}_T({{\und{\b{I}}}_{T}^{k}} \b w,{\und{\b{v}}}_T)-{\mathrm{s}}_T({{\und{\b{I}}}_{T}^{k}}({\b{\pi}_{T}^{k+1}} \b w),{\und{\b{v}}}_T)|
      \\ &\lesssim {\mathrm{s}}_T({{\und{\b{I}}}_{T}^{k}} \b w,{{\und{\b{I}}}_{T}^{k}} \b w)^\frac{{r}-{r^{\circ}}}{{r}}{\mathrm{s}}_T({{\und{\b{I}}}_{T}^{k}}(\b w-{\b{\pi}_{T}^{k+1}} \b w),{{\und{\b{I}}}_{T}^{k}}(\b w-{\b{\pi}_{T}^{k+1}} \b w))^\frac{{r^{\circ}}-1}{{r}}{\mathrm{s}}_T({\und{\b{v}}}_T,{\und{\b{v}}}_T)^\frac{1}{{r}} \\
      &\lesssim \|{{\und{\b{I}}}_{T}^{k}} \b w\|_{{\boldsymbol{\varepsilon}},{r},T}^{{r}-{r^{\circ}}}\|{{\und{\b{I}}}_{T}^{k}}(\b w-{\b{\pi}_{T}^{k+1}} \b w)\|_{{\boldsymbol{\varepsilon}},{r},T}^{{r^{\circ}}-1}\|{\und{\b{v}}}_T\|_{{\boldsymbol{\varepsilon}},{r},T}\\
      &\lesssim | \b w |_{W^{1,{r}}(T,{\mathbb{R}}^d)}^{{r}-{r^{\circ}}}|\b w-{\b{\pi}_{T}^{k+1}} \b w|_{W^{1,{r}}(T,{\mathbb{R}}^d)}^{{r^{\circ}}-1}\|{\und{\b{v}}}_T\|_{{\boldsymbol{\varepsilon}},{r},T}\\
      &\lesssim h_T^{(k+1)({r^{\circ}}-1)}| \b w |_{W^{1,{r}}(T,{\mathbb{R}}^d)}^{{r}-{r^{\circ}}}|\b w|_{W^{k+2,{r}}(T,{\mathbb{R}}^d)}^{{r^{\circ}}-1}\|{\und{\b{v}}}_T\|_{{\boldsymbol{\varepsilon}},{r},T},
    \end{aligned}$$ where we have used the Hölder continuity [[(S\[ass:sT:holder-continuity\])]{}]{} and observed that, by the consistency property [[(S\[ass:sT:polynomial.consistency\])]{}]{}, ${\mathrm{s}}_T({{\und{\b{I}}}_{T}^{k}}({\b{\pi}_{T}^{k+1}} \b w),{{\und{\b{I}}}_{T}^{k}}({\b{\pi}_{T}^{k+1}} \b w))=0$ to pass to the second line, we have used the boundedness property [[(S\[ass:sT:stability.boundedness\])]{}]{} to pass to the third line, the boundedness of ${{\und{\b{I}}}_{T}^{k}}$ to pass to the fourth line, and the $(k+2,{r},1)$-approximation property of ${\b{\pi}_{T}^{k+1}}$ to conclude.

In what follows, we will need generalized versions of the continuous and discrete Hölder inequalities, recalled hereafter for the sake of convenience. Let $X \subset {\mathbb{R}}^d$ be measurable, $n \in {\mathbb{N}}^*$, and let $t,p_1,\ldots,p_n \in (0,+\infty\rbrack$ be such that $\sum_{i=1}^n\frac{1}{p_i} = \frac{1}{t}$. The continuous $(t;p_1,\ldots,p_n)$-Hölder inequality reads: For any $(f_1,\ldots,f_n) \in \bigtimes_{i=1}^n L^{p_i}(X,{\mathbb{R}})$, $$\label{eq:holder}
  \left\| \prod_{i=1}^n f_i \right\|_{L^t(X,{\mathbb{R}})} \le\ \prod_{i=1}^n\| f_i \|_{L^{p_i}(X,{\mathbb{R}})}.$$ Let $m \in {\mathbb{N}}^*$. For all $f : \{1,\ldots,m\} \to {\mathbb{R}}$ and all $q \in [1,+\infty)$, setting $\|f\|_q \coloneqq \left(\sum_{i=1}^m |f(i)|^q\right)^\frac{1}{q}$, and $\|f\|_\infty \coloneqq \max_{1 \le i \le m}|f(i)|$, the discrete $(t;p_1,\ldots,p_n)$-Hölder inequality reads: For any $f_1,\ldots,f_n : \{1,\ldots,m\} \to {\mathbb{R}}$, $$\label{eq:discrete.holder}
  \left\| \prod_{i=1}^n f_i \right\|_{t} \le\ \prod_{i=1}^n\| f_i \|_{p_i}.$$

Let ${\mathrm{s}}_h$ be given by with, for all $T\in{\mathcal{T}}_h$, ${\mathrm{s}}_T$ satisfying Assumption \[ass:sT\]. Then it holds, for all ${\und{\b{v}}}_h \in{{{\und{\b{U}}}}_{h}^{k}}$,

\[eq:sh:properties\] $$\label{eq:sh:stability.boundedness}
      \| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}+ {\mathrm{s}}_h({\und{\b{v}}}_h,{\und{\b{v}}}_h)
      \simeq \|{\und{\b{v}}}_h\|_{{\boldsymbol{\varepsilon}},{r},h}^{r}.$$ Furthermore, for all ${\und{\b{u}}}_h, {\und{\b{v}}}_h,{\und{\b{w}}}_h\in{{{\und{\b{U}}}}_{h}^{k}}$ it holds, setting ${\und{\b{e}}}_h\coloneq{\und{\b{u}}}_h - {\und{\b{w}}}_h$, $$\begin{gathered}
\label{eq:sh:holder-continuity}
      \left|{\mathrm{s}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h)-{\mathrm{s}}_h({\und{\b{w}}}_h,{\und{\b{v}}}_h)\right|
      \lesssim
      \left({\mathrm{s}}_h({\und{\b{u}}}_h,{\und{\b{u}}}_h)+{\mathrm{s}}_h({\und{\b{w}}}_h,{\und{\b{w}}}_h)\right)^\frac{{r}-{r^{\circ}}}{{r}}{\mathrm{s}}_h({\und{\b{e}}}_h,{\und{\b{e}}}_h)^\frac{{r^{\circ}}-1}{{r}}{\mathrm{s}}_h({\und{\b{v}}}_h,{\und{\b{v}}}_h)^\frac{1}{{r}},
      \\\label{eq:sh:strong-monotonicity}
      \left(
      {\mathrm{s}}_h({\und{\b{u}}}_h,{\und{\b{e}}}_h) - {\mathrm{s}}_h({\und{\b{w}}}_h,{\und{\b{e}}}_h)
      \right)\left(
      {\mathrm{s}}_h({\und{\b{u}}}_h,{\und{\b{u}}}_h)+{\mathrm{s}}_h({\und{\b{w}}}_h,{\und{\b{w}}}_h)
      \right)^\frac{2-{r^{\circ}}}{{r}} \gtrsim {\mathrm{s}}_h({\und{\b{e}}}_h,{\und{\b{e}}}_h)^\frac{{r}+2-{r^{\circ}}}{{r}} .
    \end{gathered}$$

Finally, for any $\b w \in \b U\cap W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)$, it holds $$\label{eq:sh:consist}
    \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}},\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1}{\mathrm{s}}_h({{\und{\b{I}}}_{h}^{k}} \b w,{\und{\b{v}}}_h)  \lesssim h^{(k+1)({r^{\circ}}-1)}| \b w |_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)}^{{r}-{r^{\circ}}}|\b w|_{W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{{r^{\circ}}-1}.$$ Above, the hidden constants are independent of $h$ and of the arguments of ${\mathrm{s}}_h$.

For the sake of conciseness, we only sketch the proof and leave the details to the reader. Summing over $T \in {\mathcal{T}}_h$ immediately yields . The Hölder continuity property follows applying to the quantity in the left-hand side triangle inequalities, using , and concluding with a discrete $(1;\frac{{r}}{{r}-{r^{\circ}}},\frac{{r}}{{r^{\circ}}-1},{r})$-Hölder inequality. Moving to , starting from $|{\mathrm{s}}_h({\und{\b{e}}}_h,{\und{\b{e}}}_h)|$, we use and apply a discrete $(1;\frac{{r}+2-{r^{\circ}}}{2-{r^{\circ}}},\frac{{r}+2-{r^{\circ}}}{{r}})$-Hölder inequality to conclude. Finally, to prove we start from ${\mathrm{s}}_h({\und{\b{I}}}_h^k\b w,{\und{\b{v}}}_h)$, expand this quantity according to , use, for all $T \in {\mathcal{T}}_h$, the local consistency property together with $h_T \le h$, invoke the discrete $(1;\frac{{r}}{{r}-{r^{\circ}}},\frac{{r}}{{r^{\circ}}-1},{r})$-Hölder inequality, and pass to the supremum to conclude.

### An example of viscous stabilization function

Taking inspiration from the scalar case (cf., e.g., [@Di-Pietro.Droniou:17 Eq. (4.11c)]), a local stabilization function that matches Assumption \[ass:sT\] can be obtained setting, for all ${\und{\b{v}}}_T,{\und{\b{w}}}_T \in {{{\und{\b{U}}}}_{T}^{k}}$, $$\label{eq:sT}
  {\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{v}}}_T) \coloneqq \int_{\partial T} |{\b{\Delta}^{k}_{\partial T}} {\und{\b{w}}}_T|^{{r}-2}{\b{\Delta}^{k}_{\partial T}} {\und{\b{w}}}_T \cdot {\b{\Delta}^{k}_{\partial T}} {\und{\b{v}}}_T,$$ where, denoting by ${\mathbb{P}}^k({\mathcal{F}}_T,{\mathbb{R}}^d)$ the space of vector-valued broken polynomials of total degree $\le k$ on ${\mathcal{F}}_T$, the boundary residual operator ${\b{\Delta}^{k}_{\partial T}} : {{{\und{\b{U}}}}_{T}^{k}} \to {\mathbb{P}}^k({\mathcal{F}}_T,{\mathbb{R}}^d)$ is such that, for all ${\und{\b{v}}}_T \in {{{\und{\b{U}}}}_{T}^{k}}$, $$({\b{\Delta}^{k}_{\partial T}} {\und{\b{v}}}_T){\ \!\!_{|_{F}}} \coloneq
h_F^{-\frac1{{r}'}}\left(
  {\b{\pi}_{F}^{k}}({\b{{\mathrm{r}}}^{k+1}_{T}} {\und{\b{v}}}_T-\b v_F)-{\b{\pi}_{T}^{k}}({\b{{\mathrm{r}}}^{k+1}_{T}} {\und{\b{v}}}_T-\b v_T)
  \right)
\qquad\forall F\in{\mathcal{F}}_T,$$ with velocity reconstruction ${\b{{\mathrm{r}}}^{k+1}_{T}} : {{{\und{\b{U}}}}_{T}^{k}} \to {\mathbb{P}}^{k+1}(T,{\mathbb{R}}^d)$ such that $$\begin{gathered}
  \displaystyle\int_T ({\b{\nabla}_{\mathrm{s}}}{\b{{\mathrm{r}}}^{k+1}_{T}} {\und{\b{v}}}_T - {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},T}} {\und{\b{v}}}_T) : {\b{\nabla}_{\mathrm{s}}}\b w = 0\qquad\forall\b w \in  {\mathbb{P}}^{k+1}(T,{\mathbb{R}}^d), \\
  \text{
    $\displaystyle\int_T {\b{{\mathrm{r}}}^{k+1}_{T}} {\und{\b{v}}}_T = \int_T \b v_T$, and
    $\int_T {\b{\nabla}_{\mathrm{ss}}}{\b{{\mathrm{r}}}^{k+1}_{T}} {\und{\b{v}}}_T = \frac{1}{2}\sum_{F \in {\mathcal{F}}_T} \int_F ({\b v_F \otimes \b n_{TF}} - {\b n_{TF} \otimes \b v_F})$.
  }\end{gathered}$$ Above, ${\b{\nabla}_{\mathrm{ss}}}$ denotes the skew-symmetric part of the gradient operator ${\b\nabla}$ applied to vector fields and ${ \otimes }$ is the tensor product such that, for all $\b x = (x_i)_{1 \le i \le d}$ and $\b y = (y_i)_{1 \le i \le d}$ in ${\mathbb{R}}^d$, ${\b x \otimes \b y} \coloneq (x_i y_j)_{1 \le i,j \le d} \in {{\mathbb{R}}^{d \times d}}$.

\[lem:sT\] The local stabilization function defined by satisfies Assumption \[ass:sT\].

The proof of [[(S\[ass:sT:stability.boundedness\])]{}]{} for ${r}=2$ is given in [@Botti.Di-Pietro.ea:17 Eq. (25)]. The result can be generalized to ${r}\neq 2$ using the same arguments of [@Di-Pietro.Droniou:17 Lemma 5.2]. Property [[(S\[ass:sT:polynomial.consistency\])]{}]{} is an immediate consequence of the fact that ${\b{\Delta}^{k}_{\partial T}}({{\und{\b{I}}}_{T}^{k}}\b w) = \b 0$ for any $\b w\in{\mathbb{P}}^{k+1}(T,{\mathbb{R}}^d)$, which can be proved reasoning as in [@Di-Pietro.Droniou:20 Proposition 2.6].

Let us prove [[(S\[ass:sT:holder-continuity\])]{}]{}. First, we remark that, since the function $\alpha \mapsto \alpha^{{r}-2}$ verifies the conditions in , we can apply Theorem \[thm:1d.power-framed\] to infer that the function ${\mathbb{R}}^d\ni\b x \mapsto |\b x|^{{r}-2}\b x$ satisfies for all $\b x,\b y \in {\mathbb{R}}^d$,

\[eq:rho.holder.strong\] $$\begin{gathered}
      \big|
      |\b x|^{{r}-2}\b x-|\b y|^{{r}-2}\b y
      \big| \lesssim \big(
      |\b x|^{r}+|\b y|^{r}\big)^\frac{{r}-{r^{\circ}}}{{r}}|\b x-\b y|^{{r^{\circ}}-1},\label{eq:rho:s.holder.continuity}
      \\
      \big(
      |\b x|^{{r}-2}\b x-|\b y|^{{r}-2}\b y
      \big) \cdot (\b x-\b y) \big(
      |\b x|^{r}+|\b y|^{r}\big)^\frac{2-{r^{\circ}}}{{r}} \gtrsim |\b x-\b y|^{{r}+2-{r^{\circ}}}.\label{eq:rho:s.strong.monotonicity}
    \end{gathered}$$

Recalling , we can write $$\begin{aligned}
    \left|{\mathrm{s}}_T({\und{\b{u}}}_T,{\und{\b{v}}}_T)-{\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{v}}}_T)\right| &\leq  \int_{\partial T} \left||{\b{\Delta}^{k}_{\partial T}} {\und{\b{u}}}_T|^{{r}-2}{\b{\Delta}^{k}_{\partial T}} {\und{\b{u}}}_T-|{\b{\Delta}^{k}_{\partial T}} {\und{\b{w}}}_T|^{{r}-2}{\b{\Delta}^{k}_{\partial T}} {\und{\b{w}}}_T\right|| {\b{\Delta}^{k}_{\partial T}} {\und{\b{v}}}_T|\\
    &\lesssim  \int_{\partial T} \left(|{\b{\Delta}^{k}_{\partial T}} {\und{\b{u}}}_T|^{r}+|{\b{\Delta}^{k}_{\partial T}} {\und{\b{w}}}_T|^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}| {\b{\Delta}^{k}_{\partial T}} {\und{\b{e}}}_T|^{{r^{\circ}}-1}| {\b{\Delta}^{k}_{\partial T}} {\und{\b{v}}}_T|\\
    &\le \left({\mathrm{s}}_T({\und{\b{u}}}_T,{\und{\b{u}}}_T)+{\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{w}}}_T)\right)^\frac{{r}-{r^{\circ}}}{{r}}{\mathrm{s}}_T({\und{\b{e}}}_T,{\und{\b{e}}}_T)^\frac{{r^{\circ}}-1}{{r}}{\mathrm{s}}_T({\und{\b{v}}}_T,{\und{\b{v}}}_T)^\frac{1}{{r}},
  \end{aligned}$$ where we have used to pass to the second line and the $(1;\frac{{r}}{{r}-{r^{\circ}}},\frac{{r}}{{r^{\circ}}-1},{r})$-Hölder inequality to conclude.

Moving to [[(S\[ass:sT:strong-monotonicity\])]{}]{}, and the $(1;\frac{{r}+2-{r^{\circ}}}{2-{r^{\circ}}},\frac{{r}+2-{r^{\circ}}}{{r}})$-Hölder inequality yield $$\begin{aligned}
    &{\mathrm{s}}_T({\und{\b{e}}}_T,{\und{\b{e}}}_T) \\
    &\quad = \int_{\partial T} |{\b{\Delta}^{k}_{\partial T}}{\und{\b{u}}}_T-{\b{\Delta}^{k}_{\partial T}}{\und{\b{w}}}_T|^{r}\\
    &\quad
    \lesssim \int_{\partial T} \left(|{\b{\Delta}^{k}_{\partial T}} {\und{\b{u}}}_T|^{r}+|{\b{\Delta}^{k}_{\partial T}} {\und{\b{w}}}_T|^{r}\right)^\frac{2-{r^{\circ}}}{{r}+2-{r^{\circ}}}\left[
      \left(|{\b{\Delta}^{k}_{\partial T}}{\und{\b{u}}}_T|^{{r}-2}{\b{\Delta}^{k}_{\partial T}}{\und{\b{u}}}_T-|{\b{\Delta}^{k}_{\partial T}}{\und{\b{w}}}_T|^{{r}-2}{\b{\Delta}^{k}_{\partial T}}{\und{\b{w}}}_T\right)\cdot {\b{\Delta}^{k}_{\partial T}}{\und{\b{e}}}_T
      \right]^\frac{{r}}{{r}+2-{r^{\circ}}}\\
    &\quad
    \le \left( {\mathrm{s}}_T({\und{\b{u}}}_T,{\und{\b{u}}}_T)+{\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{w}}}_T)\right)^\frac{2-{r^{\circ}}}{{r}+2-{r^{\circ}}}\left({\mathrm{s}}_T({\und{\b{u}}}_T,{\und{\b{e}}}_T)-{\mathrm{s}}_T({\und{\b{w}}}_T,{\und{\b{e}}}_T)\right)^\frac{{r}}{{r}+2-{r^{\circ}}}.\qedhere
  \end{aligned}$$

Pressure-velocity coupling
--------------------------

For all $T \in {\mathcal{T}}_h$, we define the local divergence reconstruction ${{\mathrm{D}}^{k}_{T}} : {{{\und{\b{U}}}}_{T}^{k}} \to {\mathbb{P}}^{k}(T,{\mathbb{R}})$ by setting, for all ${\und{\b{v}}}_T \in {{{\und{\b{U}}}}_{T}^{k}}$, ${{\mathrm{D}}^{k}_{T}}{\und{\b{v}}}_T \coloneq {\mathrm{tr}}({\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},T}}{\und{\b{v}}}_T)$. We have the following characterization of ${{\mathrm{D}}^{k}_{T}}$: For all ${\und{\b{v}}}_T \in {{{\und{\b{U}}}}_{T}^{k}}$, $$\label{eq:D}
  \int_T {{\mathrm{D}}^{k}_{T}} {\und{\b{v}}}_T~ q = \int_T ({\b\nabla{\cdot}}\b v_T)~q
  + \sum_{F \in {\mathcal{F}}_T} \int_F (\b v_F-\b v_T)\cdot \b n_{TF}~q \qquad \forall q \in  {\mathbb{P}}^{k}(T,{\mathbb{R}}),$$ as can be checked writing for $\b\tau = q\mathrm{I}_d$. Taking the trace of , it is inferred that, for all $T\in{\mathcal{T}}_h$ and all $\b v\in W^{1,1}(T,{\mathbb{R}}^d)$, ${{\mathrm{D}}^{k}_{T}} ({{\und{\b{I}}}_{T}^{k}} \b v) = {\pi_{T}^{k}}({\b\nabla{\cdot}}\b v)$. The pressure-velocity coupling is realized by the bilinear form ${\mathrm{b}}_h : {{{\und{\b{U}}}}_{h}^{k}} \times {\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}}) \to {\mathbb{R}}$ such that, for all $({\und{\b{v}}}_h,q_h) \in {{{\und{\b{U}}}}_{h}^{k}} \times {\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}})$, setting $q_T\coloneq (q_h){\ \!\!_{|_{T}}}$ for all $T\in{\mathcal{T}}_h$, $$\label{eq:bh}
  {\mathrm{b}}_h({\und{\b{v}}}_h,q_h) \coloneqq -\sum_{T \in {\mathcal{T}}_h}\int_T  {{\mathrm{D}}^{k}_{T}} {\und{\b{v}}}_T~ q_T.$$

Discrete problem
----------------

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

\[eq:stokes.discrete\] $$\begin{aligned}
{2}
    \label{eq:stokes.discrete:momentum} {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h) + {\mathrm{b}}_h({\und{\b{v}}}_h,p_h) &= \displaystyle\int_\Omega \b f \cdot \b v_h &\qquad& \forall {\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}, \\
    \label{eq:stokes.discrete:mass} -{\mathrm{b}}_h({\und{\b{u}}}_h,q_h) &= 0 &\qquad& \forall q_h \in {P_{h}^{k}}. 
  \end{aligned}$$

Before proceding, some remarks are in order.

The space of test functions in can be extended to ${\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}})$ since, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}$, the divergence theorem together with the fact that $\b v_F = \b 0$ for all $F \in {{\mathcal{F}}_h^{\mathrm{b}}}$ and $\sum_{T\in {\mathcal{T}}_F} \int_F \b v_F \cdot \b n_{TF} = 0$ for all $F \in {{\mathcal{F}}_h^{\mathrm{i}}}$, yield $${\mathrm{b}}_h({\und{\b{v}}}_h,1) = -\sum_{T \in {\mathcal{T}}_h} \sum_{F \in {\mathcal{F}}_T} \int_F \b v_F \cdot \b n_{TF} = -\sum_{F \in {{\mathcal{F}}_h^{\mathrm{i}}}}\sum_{T \in {\mathcal{T}}_F}  \int_F \b v_F \cdot \b n_{TF} = 0.$$

When solving the system of nonlinear algebraic equations corresponding to by a first-order (e.g., Newton) algorithm, all element-based velocity unknowns and all but one pressure unknown per element can be locally eliminated at each iteration by computing the corresponding Schur complement element-wise. As all the computations are local, this procedure is an embarrassingly parallel task which can fully benefit from multi-thread and multi-processor architectures. This implementation strategy has been described for the linear Stokes problem in [@Di-Pietro.Ern.ea:16 Section 6.2]. After further eliminating the boundary unknowns by strongly enforcing the boundary condition , we end up solving, at each iteration of the nonlinear solver, a linear system of size $d{\mathrm{card}}({{\mathcal{F}}_h^{\mathrm{i}}}){k+d-1\choose d-1} + {\mathrm{card}}({\mathcal{T}}_h)$.

Well-posedness {#sec:well-posedness}
==============

In this section, after studying the stability properties of the viscous function ${\mathrm{a}}_h$ and of the velocity-pressure coupling bilinear form ${\mathrm{b}}_h$, we prove the well-posedness of problem .

Hölder continuity and strong monotonicity of the viscous function
-----------------------------------------------------------------

\[lem:ah:holder.continuity.strong.monotonicity\] For all ${\und{\b{u}}}_h, {\und{\b{v}}}_h, {\und{\b{w}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$, setting ${\und{\b{e}}}_h\coloneq{\und{\b{u}}}_h - {\und{\b{w}}}_h$, it holds

\[eq:ah:holder.continuity.strong.monotonicity\] $$\begin{gathered}
\label{eq:ah:holder.continuity}
      \left|
      {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h)-{\mathrm{a}}_h({\und{\b{w}}}_h,{\und{\b{v}}}_h)
      \right| \lesssim \sigma_{\mathrm{hc}}\left( \sigma_{\mathrm{de}}^{r}+  \| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\| {\und{\b{w}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r^{\circ}}-1}\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h},
      \\\label{eq:ah:strong.monotonicity}
      \left({\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{e}}}_h)-{\mathrm{a}}_h({\und{\b{w}}}_h,{\und{\b{e}}}_h)\right)\left( \sigma_{\mathrm{de}}^{r}+\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\| {\und{\b{w}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{2-{r^{\circ}}}{{r}} \gtrsim \sigma_{\mathrm{sm}}\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r}+2-{r^{\circ}}}.
    \end{gathered}$$

\(i) *Hölder continuity.* Denote by $|\Omega|_d$ the measure of $\Omega$. Using a Cauchy–Schwarz inequality followed by the Hölder continuity of ${\b\sigma}$, we can write $$\label{C2:rge2:1}
    \begin{aligned}
      &\left|\int_\Omega
      \big(
      {\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{u}}}_h)-{\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{w}}}_h)
      \big):{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h\right| \\
      &\quad \le \sigma_{\mathrm{hc}} \int_\Omega  \left(\sigma_{\mathrm{de}}^{r}+|{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{u}}}_h|_{d\times d}^{r}+|{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{w}}}_h|_{d\times d}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{e}}}_h|_{d\times d}^{{r^{\circ}}-1} |{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h|_{d \times d}  \\
      &\quad \lesssim \sigma_{\mathrm{hc}}\left(|\Omega|_d\sigma_{\mathrm{de}}^{r}+\|  {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{u}}}_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}+\|  {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{w}}}_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}} \\
      & \qquad \times \|  {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{e}}}_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{{r^{\circ}}-1}\| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}\\
      &\quad \lesssim \sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{de}}^{r}+ \| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\| {\und{\b{w}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r^{\circ}}-1}\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h},
    \end{aligned}$$ where we have used the $(1;\frac{{r}}{{r}-{r^{\circ}}},\frac{{r}}{{r^{\circ}}-1},{r})$-Hölder inequality in the second bound and the global seminorm equivalence together with the fact that $|\Omega|_d\lesssim 1$ (since $\Omega$ is bounded) to conclude. For the stabilization term, combining the Hölder continuity of ${\mathrm{s}}_h$ and the seminorm equivalence readily gives $$\label{C2:rge2:2}
    \left|
    {\mathrm{s}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h)-{\mathrm{s}}_h({\und{\b{w}}}_h,{\und{\b{v}}}_h)
    \right|
    \lesssim \left(\sigma_{\mathrm{de}}^{r}+  \| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\| {\und{\b{w}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r^{\circ}}-1}\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h},$$ where we have additionally noticed that $\sigma_{\mathrm{de}}^{r}\ge 0$ to add this term to the quantity inside parentheses. Using the definition of ${\mathrm{a}}_h$, a triangle inequality followed by and , and recalling that $\gamma \le \sigma_{\mathrm{hc}}$ (cf. ), follows.\
(ii) *Strong monotonicity.* Using the strong monotonicity of ${\b\sigma}$ and the $(1;\frac{{r}+2-{r^{\circ}}}{2-{r^{\circ}}},\frac{{r}+2-{r^{\circ}}}{{r}})$-Hölder inequality , we get $$\label{eq:ah:sm:1}
    \hspace{-0.2cm}\begin{aligned}
      &\sigma_{\mathrm{sm}}^\frac{{r}}{{r}+2-{r^{\circ}}}\| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{e}}}_h\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}\\
      &\leq \int_\Omega \left(\sigma_{\mathrm{de}}^{r}+|{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{u}}}_h|_{d\times d}^{r}+|{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{w}}}_h|_{d\times d}^{r}\right)^{\frac{2-{r^{\circ}}}{{r}+2-{r^{\circ}}}}\left)
      \left({\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{u}}}_h)-{\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{w}}}_h)\right):{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{e}}}_h
      \right)^\frac{{r}}{{r}+2-{r^{\circ}}}\\
      &\lesssim \left(\sigma_{\mathrm{de}}^{r}+\| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{u}}}_h\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}+\| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{w}}}_h\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}\right)^\frac{2-{r^{\circ}}}{{r}+2-{r^{\circ}}} \\
      &\qquad \times \left(
        \int_\Omega \left(
        {\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{u}}}_h)-{\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{w}}}_h)
        \right):{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{e}}}_h
        \right)^\frac{{r}}{{r}+2-{r^{\circ}}}\\
      &\lesssim \left(\sigma_{\mathrm{de}}^{r}+\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\| {\und{\b{w}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{2-{r^{\circ}}}{{r}+2-{r^{\circ}}} \left(
        \int_\Omega \left(
        {\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{u}}}_h)-{\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{w}}}_h)
        \right):{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{e}}}_h
        \right)^\frac{{r}}{{r}+2-{r^{\circ}}},
    \end{aligned}$$ where the conclusion follows from the global seminorm equivalence . Additionally, using the strong monotonicity of ${\mathrm{s}}_h$ together with the fact that $\sigma_{\mathrm{sm}} \le \gamma$ (cf. ) and invoking again the seminorm equivalence , we readily obtain $$\label{eq:ah:sm:2}
    \sigma_{\mathrm{sm}}^\frac{{r}}{{r}+2-{r^{\circ}}}{\mathrm{s}}_h({\und{\b{e}}}_h,{\und{\b{e}}}_h)
    \lesssim \left(\sigma_{\mathrm{de}}^{r}+\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\| {\und{\b{w}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{2-{r^{\circ}}}{{r}+2-{r^{\circ}}}\left(
      \gamma{\mathrm{s}}_h({\und{\b{u}}}_h,{\und{\b{e}}}_h)-\gamma{\mathrm{s}}_h({\und{\b{w}}}_h,{\und{\b{e}}}_h)
      \right)^\frac{{r}}{{r}+2-{r^{\circ}}}.$$ Finally, combining again the norm equivalence with and , and using yields $$\begin{aligned}
    \sigma_{\mathrm{sm}}^\frac{{r}}{{r}+2-{r^{\circ}}}\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\lesssim &\left( \sigma_{\mathrm{de}}^{r}+\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\| {\und{\b{w}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{2-{r^{\circ}}}{{r}+2-{r^{\circ}}}\left(
      {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{e}}}_h)-{\mathrm{a}}_h({\und{\b{w}}}_h,{\und{\b{e}}}_h)
      \right)^\frac{{r}}{{r}+2-{r^{\circ}}}.
  \end{aligned}$$ Raising this inequality to the power $\frac{{r}-2-{r^{\circ}}}{{r}}$ yields .

Stability of the pressure-velocity coupling
-------------------------------------------

\[lem:bh:inf-sup\] It holds, for all $q_h \in {P_{h}^{k}}$, $$\label{eq:bh:inf-sup}
    \| q_h \|_{L^{{r}'}(\Omega,{\mathbb{R}})} \lesssim \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}, \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1} {\mathrm{b}}_h({\und{\b{v}}}_h,q_h),$$ with hidden constant depending only on $d$, $k$, ${r}$, $\Omega$, and the mesh regularity parameter.

The proof follows the classical Fortin argument (cf., e.g., [@Boffi.Brezzi.ea:13 Section 8.4]), adapted here to the non-Hilbertian setting: we first prove that ${{\und{\b{I}}}_{h}^{k}}$ is a Fortin operator, then combine this fact with the continuous inf-sup condition.\
(i) *Fortin operator.* We need to prove that the following properties hold for any $\b v\in W^{1,{r}}(\Omega,{\mathbb{R}}^d)$:

\[eq:fortin\] $$\begin{gathered}
      \| {{\und{\b{I}}}_{h}^{k}} \b v\|_{{\boldsymbol{\varepsilon}},{r},h} \lesssim | \b v |_{W^{1,{r}}(\Omega,{\mathbb{R}}^d),}\label{eq:fortin:boundedness}
      \\ 
      {\mathrm{b}}_h({{\und{\b{I}}}_{h}^{k}} \b v,q_h) = b(\b v,q_h)\qquad\forall q_h\in{\mathbb{P}}^k({\mathcal{T}}_h,{\mathbb{R}}).\label{eq:fortin:consistency}
    \end{gathered}$$

Property is obtained by raising both sides of to the power ${r}$, summing over $T \in {\mathcal{T}}_h$, then taking the $r$th root of the resulting inequality. The proof of is given, e.g., in [@Di-Pietro.Droniou:20 Lemma 8.12].\
(ii) *Inf-sup condition on ${\mathrm{b}}_h$.* Let $q_h \in {P_{h}^{k}}$ and set $c_h \coloneqq \int_\Omega |q_h|^{{r}'-2}q_h$. Using a triangle inequality, the Hölder inequality, and the fact that $|\Omega|_d\lesssim 1$, we get $$\label{eq:inf-sup:qT}
  \| |q_h|^{{r}'-2}q_h-c_h \|_{L^{{r}}(\Omega,{\mathbb{R}})}
  \le \| q_h\|_{L^{{r}'}(\Omega,{\mathbb{R}})}^{{r}'-1}+|c_h||\Omega|_d^\frac{1}{{r}}
  \le \left(1+|\Omega|_d\right)\| q_h\|_{L^{{r}'}(\Omega,{\mathbb{R}})}^{{r}'-1} \lesssim \| q_h\|_{L^{{r}'}(\Omega,{\mathbb{R}})}^{{r}'-1},$$ where we have used the fact that $|c_h|\le\| q_h\|_{L^{{r}'}(\Omega,{\mathbb{R}})}^{{r}'-1} |\Omega|_d^{\frac1{{r}'}}$ along with $\frac1{r}+\frac1{{r}'}=1$ in the second bound and the fact that $|\Omega|_d\lesssim 1$ to conclude. Since $q_h \in L^{{r}'}(\Omega,{\mathbb{R}})$, bound implies that $|q_h|^{{r}'-2}q_h-c_h \in L^{{r}}_0(\Omega,{\mathbb{R}}) \coloneq \left\{ q \in L^{r}(\Omega,{\mathbb{R}}) : \int_\Omega q = 0 \right\}$ by construction. Thus, using the surjectivity of the continuous divergence operator ${\b\nabla{\cdot}}: \b U \to L^{{r}}_0(\Omega,{\mathbb{R}})$, (c.f. [@Duran.Muschietti.ea:10] and also [@Bogovski:79 Theorem 1]), we infer that there exists $\b v_{q_h} \in \b U$ such that $$\label{eq:velocity.lifting}
    -{\b\nabla{\cdot}}\b v_{q_h} = |q_h|^{{r}'-2}q_h-c_h \quad {\mathrm{and}} \quad | \b v_{q_h}|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)} \lesssim \| |q_h|^{{r}'-2}q_h-c_h \|_{L^{r}(\Omega,{\mathbb{R}})}.$$ Denote by $\$$ the supremum in . Using the fact that $q_h$ has zero mean value over $\Omega$, the equality in together with the definition of $b$, and the second Fortin property , we have $$\|q_h \|_{L^{{r}'}(\Omega,{\mathbb{R}})}^{{r}'}
  {=} \int_\Omega \big(|q_h|^{{r}'-2}q_h-c_h\big) q_h
  = b(\b v_{q_h},q_h)
  = {\mathrm{b}}_h({{\und{\b{I}}}_{h}^{k}}\b v_{q_h},q_h) 
  \le \$  \| {{\und{\b{I}}}_{h}^{k}} \b v_{q_h} \|_{{\boldsymbol{\varepsilon}},{r},h}
  \lesssim \$ \| q_h \|_{L^{{r}'}(\Omega,{\mathbb{R}})}^{{r}'-1},$$ where, to conclude, we have used followed by and . Simplifying yields .

Well-posedness {#well-posedness}
--------------

We are now ready to prove the main result of this section.

\[thm:well-posedness\] There exists a unique solution $({\und{\b{u}}}_h,p_h) \in {{{\und{\b{U}}}}_{h,0}^{k}} \times {P_{h}^{k}}$ to the discrete problem . Additionally, the following a priori bounds hold:

\[eq:discrete.solution:bounds\] $$\begin{aligned}
      \| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} &\lesssim \left(\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\right)^\frac{1}{{r}-1}+\left(\sigma_{\mathrm{de}}^{2-{r^{\circ}}}\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\right)^\frac{1}{{r}+1-{r^{\circ}}}, \label{eq:discrete.solution:bounds:uh}\\
      \| p_h \|_{L^{{r}'}(\Omega,{\mathbb{R}})} &\lesssim \sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}+\sigma_{\mathrm{de}}^{|{r}-2|({r^{\circ}}-1)}\left(\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\right)^\frac{{r^{\circ}}-1}{{r}+1-{r^{\circ}}}\right). \label{eq:discrete.solution:bounds:ph}
    \end{aligned}$$

\(i) *Existence.* Denote by ${P_{h}^{k,*}}$ the dual space of ${P_{h}^{k}}$ and let $B_h : {{{\und{\b{U}}}}_{h,0}^{k}} \to {P_{h}^{k,*}}$ be such that, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}$, $$\langle B_h{\und{\b{v}}}_h, q_h \rangle \coloneqq -{\mathrm{b}}_h({\und{\b{v}}}_h,q_h) \qquad \forall q_h \in {P_{h}^{k}}.$$ Here and in what follows, $\langle{\cdot},{\cdot}\rangle$ denotes the appropriate duality pairing as inferred from its arguments. Define the following subspace of ${{{\und{\b{U}}}}_{h,0}^{k}}$ spanned by vectors of discrete unknowns with zero discrete divergence: $$\label{eq:Whk}
    {\und{\b{W}}}_h^k \coloneq \operatorname{Ker}(B_h) = \left\{
    {\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}} : {\mathrm{b}}_h({\und{\b{v}}}_h,q_h) = 0 \quad \forall q_h \in {P_{h}^{k}}
    \right\},$$ and consider the following problem: Find ${\und{\b{u}}}_h\in{\und{\b{W}}}_h^k$ such that $$\label{eq:existence:auxiliary}
    {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h) = \int_\Omega \b f\cdot \b v_h\qquad\forall {\und{\b{v}}}_h\in{\und{\b{W}}}_h^k.$$ Existence of a solution to this problem for a fixed $h$ can be proved adapting the arguments of [@Di-Pietro.Droniou:17 Theorem 4.5]. Specifically, equip ${\und{\b{W}}}_h^k$ with an inner product $(\cdot,\cdot)_{\b W,h}$ (which need not be further specified), denote by $\|{\cdot}\|_{\b W,h}$ the induced norm, and let $\b\Phi_h:{\und{\b{W}}}_h^k\to {\und{\b{W}}}_h^k$ be such that, for all ${\und{\b{w}}}_h\in{\und{\b{W}}}_h^k$, $(\b\Phi_h({\und{\b{w}}}_h),{\und{\b{v}}}_h)_{\b W,h} = {\mathrm{a}}_h({\und{\b{w}}}_h,{\und{\b{v}}}_h)$ for all ${\und{\b{v}}}_h\in{\und{\b{W}}}_h^k$. The strong monotonicity of ${\mathrm{a}}_h$ yields, for any ${\und{\b{v}}}_h\in{\und{\b{W}}}_h^k$ such that $\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} \ge \sigma_{\mathrm{de}}$, $$(\b\Phi_h({\und{\b{v}}}_h),{\und{\b{v}}}_h)_{\b W,h}\ge \sigma_{\mathrm{sm}} (\sigma_{\mathrm{de}}^{r}+ \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r})^\frac{{r^{\circ}}-2}{{r}} \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r}+2-{r^{\circ}}} \gtrsim \sigma_{\mathrm{sm}}\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\ge C^{r}\sigma_{\mathrm{sm}}\| {\und{\b{v}}}_h\|_{\b W,h}^{r},$$ where $C$ denotes the constant (possibly depending on $h$) in the equivalence of the norms $\|{\cdot}\|_{{\boldsymbol{\varepsilon}},{r},h}$ and $\|{\cdot}\|_{\b W,h}$ (which holds since ${\und{\b{W}}}_h^k$ is finite-dimensional). This shows that $\b\Phi_h$ is coercive hence, by [@Deimling:85 Theorem 3.3], surjective. Let now ${\und{\b{w}}}_h\in{\und{\b{W}}}_h^k$ be such that $({\und{\b{w}}}_h,{\und{\b{v}}}_h)_{\b W,h}=\int_\Omega \b f\cdot \b v_h$ for all ${\und{\b{v}}}_h\in{\und{\b{W}}}_h^k$. By the surjectivity of $\b\Phi_h$, there exists ${\und{\b{u}}}_h\in {\und{\b{W}}}_h^k$ such that $\b\Phi_h({\und{\b{u}}}_h)={\und{\b{w}}}_h$ which, by definition of ${\und{\b{w}}}_h$ and $\b\Phi_h$, is a solution to the discrete problem .

The proof of existence now continues as in the linear case; see, e.g., [@Boffi.Brezzi.ea:13 Theorem 4.2.1]. Denote by ${{{\und{\b{U}}}}_{h,0}^{k,*}}$ the dual space of ${{{\und{\b{U}}}}_{h,0}^{k}}$ and consider the linear mapping $\ell_h\in{{{\und{\b{U}}}}_{h,0}^{k,*}}$ such that, for all ${\und{\b{v}}}_h\in{{{\und{\b{U}}}}_{h,0}^{k}}$, $$\langle\ell_h,{\und{\b{v}}}_h\rangle\coloneq \int_\Omega \b f\cdot \b v_h - {\mathrm{a}}_h({\und{\b{u}}}_h, {\und{\b{v}}}_h).$$ Thanks to , $\ell_h$ vanishes identically for every ${\und{\b{v}}}_h\in{\und{\b{W}}}_h^k$, that is to say, $\ell_h$ lies in the polar space of ${\und{\b{W}}}_h^k$ which, denoting by $B_h^*:{P_{h}^{k}}\to{{{\und{\b{U}}}}_{h,0}^{k,*}}$ the adjoint operator of $B_h$, coincides in our case with $\operatorname{Im}(B_h^*)$ (see, e.g., [@Boffi.Brezzi.ea:13 Theorem 4.14]). Hence, $\ell_h\in\operatorname{Im}(B_h^*)$, and there exists therefore a $p_h\in {P_{h}^{k}}$ such that $B_h^* p_h = \ell_h$. This means that, for all ${\und{\b{v}}}_h\in{{{\und{\b{U}}}}_{h,0}^{k}}$, $${\mathrm{b}}_h({\und{\b{v}}}_h,p_h)
  = \langle B_h^* p_h,{\und{\b{v}}}_h\rangle
  = \langle\ell_h,{\und{\b{v}}}_h\rangle
  = \int_\Omega \b f\cdot \b v_h - {\mathrm{a}}_h({\und{\b{u}}}_h, {\und{\b{v}}}_h),$$ i.e., the $({\und{\b{u}}}_h,p_h)$ satisfies the discrete momentum equation . On the other hand, since ${\und{\b{u}}}_h\in{\und{\b{W}}}_h^k$, we also have, by the definition of ${\und{\b{W}}}_h^k$, $  {\mathrm{b}}_h({\und{\b{u}}}_h,q_h) = 0$ for all $q_h\in {P_{h}^{k}}$, which shows that the discrete mass equation is also verified. In conclusion, $({\und{\b{u}}}_h,p_h)\in{{{\und{\b{U}}}}_{h,0}^{k}}\times {P_{h}^{k}}$ solves .\
(ii) *Uniqueness.* We start by proving uniqueness for the velocity. Let $({\und{\b{u}}}_h,p_h),({\und{\b{u}}}'_h,p'_h) \in {{{\und{\b{U}}}}_{h,0}^{k}} \times {P_{h}^{k}}$ be two solutions of . Making ${\und{\b{v}}}_h = {\und{\b{u}}}_h - {\und{\b{u}}}'_h$ in written first for $({\und{\b{u}}}_h, p_h)$ then for $({\und{\b{u}}}'_h, p'_h)$, then taking the difference and observing that ${\mathrm{b}}_h({\und{\b{u}}}_h-{\und{\b{u}}}'_h,p_h)={\mathrm{b}}_h({\und{\b{u}}}_h - {\und{\b{u}}}'_h,p'_h)=0$ by , we infer that $${\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{u}}}_h - {\und{\b{u}}}'_h)-{\mathrm{a}}_h({\und{\b{u}}}'_h,{\und{\b{u}}}_h - {\und{\b{u}}}'_h) = 0.$$ Thus, the strong monotonicity of ${\mathrm{a}}_h$ yields $\| {\und{\b{u}}}_h - {\und{\b{u}}}'_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 0$, which implies ${\und{\b{u}}}_h = {\und{\b{u}}}'_h$ since $\|{\cdot}\|_{{\boldsymbol{\varepsilon}},{r},h}$ is a norm on ${{{\und{\b{U}}}}_{h,0}^{k}}$. Moreover, using the inf-sup stability of ${\mathrm{b}}_h$ and written first for ${\und{\b{u}}}_h$ then for ${\und{\b{u}}}_h'$, we get $$\begin{aligned}
    \| p_h-p'_h \|_{L^{{r}'}(\Omega,{\mathbb{R}})}
    &\lesssim \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}},\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1} {\mathrm{b}}_h({\und{\b{v}}}_h,p_h-p'_h)
    \\
    &=  \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}},\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1} \left({\mathrm{a}}_h({\und{\b{u}}}'_h,{\und{\b{v}}}_h)-{\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h)\right) = 0,
   \end{aligned}$$ hence $p_h=p'_h$.\
(iii) *A priori estimates.* Using the strong monotonicity of ${\mathrm{a}}_h$ (with ${\und{\b{w}}}_h = {\und{\b{0}}}$), equation together with , and the Hölder inequality together with the discrete Korn inequality , we obtain $$\label{eq:well-posedness:0}
    \begin{aligned}
      \sigma_{\mathrm{sm}}\big(
      \sigma_{\mathrm{de}}^{r}+ \| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\big)^\frac{{r^{\circ}}-2}{{r}} \| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r}+2-{r^{\circ}}}
      &\lesssim  {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{u}}}_h)
      = \displaystyle\int_\Omega \b f \cdot \b u_h
      \lesssim \| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}.
    \end{aligned}$$ We then conclude as in the continuous case to infer (see Remark \[rem:a-priori\]). To prove the bound on the pressure, we use the inf-sup stability of ${\mathrm{b}}_h$ to write $$\begin{aligned}
    \| p_h \|_{L^{{r}'}(\Omega,{\mathbb{R}})}
    &\lesssim \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}, \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1} {\mathrm{b}}_h({\und{\b{v}}}_h,p_h)
    \\
    &= \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}, \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}
     = 1} \left(\displaystyle\int_\Omega \b f \cdot \b v_h - {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h) \right) \\
    &\lesssim \| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)} +\sigma_{\mathrm{hc}}(\sigma_{\mathrm{de}}^{r}+\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r})^\frac{{r}-{r^{\circ}}}{{r}}\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r^{\circ}}-1} \\   
    &\lesssim \sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}+\sigma_{\mathrm{de}}^{|{r}-2|({r^{\circ}}-1)}\left(\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\right)^\frac{{r^{\circ}}-1}{{r}+1-{r^{\circ}}}\right),
  \end{aligned}$$ where we have used the discrete momentum equation to pass to the second line, the Hölder and discrete Korn inequalities together with the Hölder continuity of ${\mathrm{a}}_h$ to pass to the third line, and the a priori bound on the velocity together with $\frac{\sigma_{\mathrm{hc}}}{\sigma_{\mathrm{sm}}}\ge 1$ (see ) to conclude.

Error estimate {#sec:error.estimate}
==============

In this section, after studying the consistency of the viscous and pressure-velocity coupling terms, we prove an energy error estimate.

Consistency of the viscous function {#sec:consistency:viscous.function}
-----------------------------------

\[lem:consistency:ah\] Let $\b w \in \b U \cap W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)$ be such that ${\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w) \in W^{1,{r}'}(\Omega,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}}) \cap W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}})$. Define the viscous consistency error linear form $\mathcal E_{{\mathrm{a}},h}(\b w;\cdot) : {{{\und{\b{U}}}}_{h}^{k}} \to {\mathbb{R}}$ such that, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$, $$\label{eq:Eah}
    \mathcal E_{{\mathrm{a}},h}(\b w;{\und{\b{v}}}_h) \coloneqq \int_\Omega ({\b\nabla{\cdot}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)) \cdot \b v_h + {\mathrm{a}}_h({{\und{\b{I}}}_{h}^{k}} \b w,{\und{\b{v}}}_h).$$ Then, under Assumptions \[ass:stress\] and \[ass:sT\], we have $$\begin{gathered}
\label{eq:consistency:ah}
    \sup_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}},\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1} \mathcal E_{{\mathrm{a}},h}(\b w;{\und{\b{v}}}_h)
    \lesssim
    h^{(k+1)({r^{\circ}}-1)}\bigg[
      \sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{de}}^{r}+ |\b w|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}|\b w|_{W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{{r^{\circ}}-1}
      \\
      +|{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)|_{W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}})}        
      \bigg].
  \end{gathered}$$

Let ${\und{\b{\hat w}}}_h \coloneqq {{\und{\b{I}}}_{h}^{k}} \b w$ and ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}$. Expanding ${\mathrm{a}}_h$ according to its definition in the expression of $\mathcal E_{{\mathrm{a}},h}$, inserting $\pm\left(
  \int_\Omega{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w): {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h
  +  \int_\Omega{\b{\pi}_{h}^{k}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w) : {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h
  \right)$, and rearranging, we obtain $$\begin{gathered}
\label{eq:consistency:ah:EJ}
    \mathcal E_{{\mathrm{a}},h}(\b w;{\und{\b{v}}}_h)
    =\!\!
    \underbrace{      \int_\Omega ({\b\nabla{\cdot}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)) \cdot \b v_h
      {+} \int_\Omega {\b{\pi}_{h}^{k}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)\! :\! {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h
    }_{\mathcal T_1}
    +\! \underbrace{      \int_\Omega\! \left( {\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w) - {\b{\pi}_{h}^{k}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)\right)\! :\! {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h
    }_{\mathcal T_2}
    \\
    {+} \underbrace{      \int_\Omega\!\left( {\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{\hat w}}}_h) - {\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)\right)\!:\! {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h
    }_{\mathcal T_3}
    + \underbrace{\vphantom{\int_\Omega}\gamma {\mathrm{s}}_h({\und{\b{\hat w}}}_h,{\und{\b{v}}}_h)}_{\mathcal T_4}.
  \end{gathered}$$ We proceed to estimate the terms in the right-hand side. For the first term, we start by noticing that $$\label{eq:consistency:ah:null}
    \sum_{T \in {\mathcal{T}}_h}\sum_{F \in {\mathcal{F}}_T} \int_F  \b v_F \cdot\left({\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)\b n_{TF}\right) = 0$$ as a consequence of the continuity of the normal trace of ${\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)$ together with the single-valuedness of $\b v_F$ across each interface $F\in{{\mathcal{F}}_h^{\mathrm{i}}}$ and of the fact that $\b v_F=\b 0$ for every boundary face $F\in{{\mathcal{F}}_h^{\mathrm{b}}}$. Using an element by element integration by parts on the first term of $\mathcal T_1$ along with the definitions of ${\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}$ and of ${\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},T}}$, we can write $$\begin{aligned}
    \mathcal T_1
    &= \cancel{\int_\Omega \left({\b{\pi}_{h}^{k}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)- {\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)\right) : {\b{\nabla}_{{\mathrm{s}},h}}\b v_h} \\
    &\qquad + \sum_{T \in {\mathcal{T}}_h}\sum_{F \in {\mathcal{F}}_T} \left(\int_F  (\b v_F-\b v_T)\cdot({\b{\pi}_{T}^{k}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w))\b n_{TF}+\int_F\b v_T\cdot\left({\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)\b n_{TF}\right) \right) \\
    &=  \sum_{T \in {\mathcal{T}}_h}\sum_{F \in {\mathcal{F}}_T} \int_F  (\b v_F-\b v_T)\cdot\left({\b{\pi}_{T}^{k}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)-{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)\right)\b n_{TF},
  \end{aligned}$$ where we have used the definition of ${\b{\pi}_{h}^{k}}$ together with the fact that ${\b{\nabla}_{{\mathrm{s}},h}}\b v_h \in {\mathbb{P}}^{k-1}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}}) \subset {\mathbb{P}}^{k}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}})$ to cancel the term in the first line, and we have inserted and rearranged to conclude. Therefore, applying the Hölder inequality together with the bound $h_F \le h_T$, we infer $$\label{eq:consistency:ah:T1}
    \begin{aligned}
      \left|\mathcal T_1\right|
      &\le \left(\displaystyle\sum_{T \in {\mathcal{T}}_h}h_T \|{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)- {\b{\pi}_{T}^{k}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w) \|_{L^{{r}'}(\partial T,{{\mathbb{R}}^{d \times d}})}^{{r}'} \right)^\frac{1}{{r}'}\left(\sum_{T \in {\mathcal{T}}_h}\sum_{F \in {\mathcal{F}}_T}h_F^{1-{r}}\| \b v_F-\b v_T\|_{L^{r}(F,{\mathbb{R}}^d)}^{r}\right)^\frac{1}{{r}}
      \\
      &\lesssim h^{(k+1)({r^{\circ}}-1)} |{\b\sigma}(\cdot, {\b{\nabla}_{\mathrm{s}}}\b w)|_{W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}})}\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h},
    \end{aligned}$$ where the conclusion follows using the $((k+1)({r^{\circ}}-1),{r}')$-trace approximation properties of ${\b{\pi}_{T}^{k}}$ along with $h_T \le h$ for the first factor and the definition of the $\|{\cdot}\|_{{\boldsymbol{\varepsilon}},{r},h}$-norm for the second.

For the second term, we use the Hölder inequality and the seminorm equivalence to write $$\label{eq:consistency:ah:T2}
    \begin{aligned}
      \left|\mathcal T_2\right|
      &= \| {\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w)-{\b{\pi}_{h}^{k}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w) \|_{L^{{r}'}(\Omega,{{\mathbb{R}}^{d \times d}})}
      \| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}{\und{\b{v}}}_h\|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}
      \\      
      &\lesssim h^{(k+1)({r^{\circ}}-1)}|{\b\sigma}(\cdot, {\b{\nabla}_{\mathrm{s}}}\b w)|_{W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}})} \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h},
    \end{aligned}$$ where the conclusion follows from the $((k+1)({r^{\circ}}-1),{r}',0)$-approximation properties of ${\b{\pi}_{T}^{k}}$ along with $h_T\le h$ for the first factor and the global norm equivalence for the second.

For the third term, using the Hölder inequality and again , we get $$\label{eq:consistency:ah:T3:0}
    \left|\mathcal T_3\right|
    \le \|{\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{\hat w}}}_h)- {\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w) \|_{L^{{r}'}(\Omega,{{\mathbb{R}}^{d \times d}})}\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}.$$ We estimate the first factor as follows: $$\begin{aligned}
      &\|{\b\sigma}(\cdot,{\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{\hat w}}}_h)- {\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b w) \|_{L^{{r}'}(\Omega,{{\mathbb{R}}^{d \times d}})}
      \\
      &\quad
      \le \sigma_{\mathrm{hc}} \left\| \left(\sigma_{\mathrm{de}}^{r}+| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{\hat w}}}_h |_{d \times d}^{r}+  | {\b{\nabla}_{\mathrm{s}}}\b w |_{d \times d}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{\hat w}}}_h - {\b{\nabla}_{\mathrm{s}}}\b w |_{d \times d}^{{r^{\circ}}-1}\right\|_{L^{{r}'}(\Omega,{\mathbb{R}})} \\
      &\quad
      \lesssim \sigma_{\mathrm{hc}} \left(\sigma_{\mathrm{de}}^{r}+\| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{\hat w}}}_h \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}+\|  {\b{\nabla}_{\mathrm{s}}}\b w \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}\| {\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}} {\und{\b{\hat w}}}_h - {\b{\nabla}_{\mathrm{s}}}\b w \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{{r^{\circ}}-1} \\
      &\quad
      \lesssim \sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{de}}^{r}+\|{\und{\b{\hat w}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+  |\b w|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}\| {\b{\pi}_{h}^{k}}({\b{\nabla}_{\mathrm{s}}}\b w) - {\b{\nabla}_{\mathrm{s}}}\b w \|_{L^{r}(\Omega,{{\mathbb{R}}^{d \times d}})}^{{r^{\circ}}-1} \\
      &\quad
      \lesssim h^{(k+1)({r^{\circ}}-1)}\sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{de}}^{r}+|\b w|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}|\b w|_{W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{{r^{\circ}}-1},
    \end{aligned}$$ where we have used the Hölder continuity of ${\b\sigma}$ in the first bound, the $({r}';\frac{{r}}{{r}-{r^{\circ}}},\frac{{r}}{{r^{\circ}}-1})$-Hölder inequality in the second, the boundedness of $\Omega$ along with and the commutation property of ${\b{{\mathrm{G}}}^{k}_{{\mathrm{s}},h}}$ in the third, and we have concluded invoking the $(k+1,{r},0)$-approximation property of ${\b{\pi}_{T}^{k}}$. Plugging this estimate into , we get $$\label{eq:consistency:ah:T3}
    \left|\mathcal T_3\right|
    \lesssim h^{(k+1)({r^{\circ}}-1)}\sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{de}}^{r}+|\b w|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}|\b w|_{W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{{r^{\circ}}-1}
    \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}.$$

Finally, using the fact that $\gamma \le \sigma_{\mathrm{hc}}$ together with the consistency of ${\mathrm{s}}_h$ and the norm equivalence , we obtain for the fourth term $$\label{eq:consistency:ah:T4}
    \left|\mathcal T_4\right|
    \lesssim h^{(k+1)({r^{\circ}}-1)}\sigma_{\mathrm{hc}}|\b w|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)}^{{r}-{r^{\circ}}}|\b w|_{W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{{r^{\circ}}-1}
    \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}.$$

Plug the bounds , , , and into and pass to the supremum to conclude.

Consistency of the pressure-velocity coupling bilinear form
-----------------------------------------------------------

Let $q \in W^{1,{r}'}(\Omega,{\mathbb{R}}) \cap W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{\mathbb{R}})$. Let $\mathcal E_{{\mathrm{b}},h}(q;\cdot) : {{{\und{\b{U}}}}_{h}^{k}} \to {\mathbb{R}}$ be the pressure consistency error linear form such that, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h}^{k}}$, $$\label{eq:Ebh}
    \mathcal E_{{\mathrm{b}},h}(q;{\und{\b{v}}}_h) \coloneqq \int_\Omega {\b\nabla}q \cdot \b v_h - {\mathrm{b}}_h({\und{\b{v}}}_h,{\pi_{h}^{k}} q).$$ Then, we have that $$\label{eq:consistency:bh}
    \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}},\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1} \mathcal E_{{\mathrm{b}},h}(q;{\und{\b{v}}}_h) \lesssim h^{(k+1)({r^{\circ}}-1)} |q|_{W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{\mathbb{R}})}.$$

Let ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}$. Integrating by parts element by element, we can reformulate the first term in the right-hand side of as follows: $$\label{eq:consistency:bh:1}
    \displaystyle\int_\Omega {\b\nabla}q\cdot \b v_h = - \sum_{T \in {\mathcal{T}}_h} \left( \int_T q({\b\nabla{\cdot}}\b v_T) +\sum_{F \in {\mathcal{F}}_T} \int_F q(\b v_F - \b v_T)\cdot \b n_{TF} \right),$$ where the introduction of $\b v_F$ in the boundary term is justified by the fact that the jumps of $q$ vanish across interfaces by the assumed regularity and that $\b v_F=\b 0$ on every boundary face $F \in {{\mathcal{F}}_h^{\mathrm{b}}}$. On the other hand, expanding, for each $T \in {\mathcal{T}}_h$, ${{\mathrm{D}}^{k}_{T}}$ according to its definition , we get $$\label{eq:consistency:bh:2}
    -{\mathrm{b}}_h({\und{\b{v}}}_h, {\pi_{h}^{k}} q) = \sum_{T \in {\mathcal{T}}_h} \left( \int_T {\pi_{T}^{k}} q~({\b\nabla{\cdot}}\b v_T) +\sum_{F \in {\mathcal{F}}_T} \int_F {\pi_{T}^{k}} q~(\b v_F - \b v_T)\cdot \b n_{TF} \right).$$ Summing and and observing that the first terms in parentheses cancel out by the definition of ${\pi_{T}^{k}}$ since ${\b\nabla{\cdot}}\b v_T \in {\mathbb{P}}^{k-1}(T,{\mathbb{R}}) \subset {\mathbb{P}}^k(T,{\mathbb{R}})$ for all $T \in{\mathcal{T}}_h$, we can write $$\begin{aligned}
    \mathcal E_{{\mathrm{b}},h}(q;{\und{\b{v}}}_h) &= \sum_{T \in {\mathcal{T}}_h} \left( \cancel{\int_T ({\pi_{T}^{k}} q-q) ({\b\nabla{\cdot}}\b v_T)} +\sum_{F \in {\mathcal{F}}_T} \int_F ({\pi_{T}^{k}} q-q) (\b v_F - \b v_T)\cdot \b n_{TF} \right)   \\
    &\le \left(\sum_{T \in {\mathcal{T}}_h} h_T\| {\pi_{T}^{k}} q-q \|_{L^{{r}'}(\partial T,{\mathbb{R}})}^{{r}'} \right)^\frac{1}{{r}'} \left(\sum_{T \in {\mathcal{T}}_h}\sum_{F \in {\mathcal{F}}_T} h_F^{1-{r}}\| \b v_F-\b v_T \|_{L^{{r}}(F,{\mathbb{R}}^d)}^{{r}} \right)^\frac{1}{{r}} \\
    &\lesssim h^{(k+1)({r^{\circ}}-1)} |q|_{W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{\mathbb{R}})} \| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h},
  \end{aligned}$$ where we have used the Hölder inequality along with $h_F\ge h_T$ whenever $F \in {\mathcal{F}}_T$ in the second line and the $((k+1)({r^{\circ}}-1),{r}')$-trace approximation property of ${\pi_{T}^{k}}$ together with the bound $h_F \le h$ and the definition of the $\|{\cdot}\|_{{\boldsymbol{\varepsilon}},{r},h}$-norm to conclude. Passing to the supremum yields .

Error estimate {#error-estimate}
--------------

\[thm:error.estimate\] Let $(\b u,p) \in \b U \times P$ and $({\und{\b{u}}}_h,p_h) \in {{{\und{\b{U}}}}_{h,0}^{k}} \times {P_{h}^{k}}$ solve and , respectively. Assume the additional regularity $\b u \in W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)$, ${\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b u) \in W^{1,{r}'}(\Omega,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}}) \cap W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}})$, and $p \in W^{1,{r}'}(\Omega,{\mathbb{R}}) \cap W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{\mathbb{R}})$. Then, under Assumptions \[ass:stress\] and \[ass:sT\],

\[eq:error.estimate\] $$\begin{aligned}
\label{eq:error.estimate:velocity}
      \| {\und{\b{u}}}_h - {{\und{\b{I}}}_{h}^{k}} \b u \|_{{\boldsymbol{\varepsilon}},{r},h}
      &\lesssim 
      h^\frac{(k+1)({r^{\circ}}-1)}{{r}+1-{r^{\circ}}}      
      \left(\sigma_{\mathrm{sm}}^{-1}\mathcal N_{\b f}^{2-{r^{\circ}}}\mathcal N_{{\b\sigma},\b u,p}\right)^\frac{1}{{r}+1-{r^{\circ}}},
      \\
      \label{eq:error.estimate:pressure}
      \| p_h - {\pi_{h}^{k}} p \|_{L^{{r}'}(\Omega,{\mathbb{R}})}
      &\lesssim
        h^{(k+1)({r^{\circ}}-1)}\mathcal N_{{\b\sigma},\b u,p}
        + 
        h^{\frac{(k+1)({r^{\circ}}-1)^2}{{r}+1-{r^{\circ}}}} \sigma_{\mathrm{hc}}\mathcal N_{\b f}^{|{r}-2|({r^{\circ}}-1)}
        \left(\sigma_{\mathrm{sm}}^{-1}\mathcal N_{{\b\sigma},\b u,p}\right)^\frac{{r^{\circ}}-1}{{r}+1-{r^{\circ}}},
    \end{aligned}$$

where we have set, for the sake of brevity, $$\begin{aligned}
    \mathcal N_{{\b\sigma},\b u,p}
    &\coloneq
    \sigma_{\mathrm{hc}}\left(\sigma_{\mathrm{de}}^{r}+|\b u|_{W^{1,{r}}(\Omega,{\mathbb{R}}^d)}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}|\b u|_{W^{k+2,{r}}({\mathcal{T}}_h,{\mathbb{R}}^d)}^{{r^{\circ}}-1}
    \\
    &\quad
    + |{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b u)|_{W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{{\mathbb{R}}^{d \times d}})}
    + |p|_{W^{(k+1)({r^{\circ}}-1),{r}'}({\mathcal{T}}_h,{\mathbb{R}})},
    \\
    \mathcal N_{\b f} &\coloneqq \sigma_{\mathrm{de}}+\left(\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\right)^\frac{1}{{r}-1}+\left(\sigma_{\mathrm{de}}^{2-{r^{\circ}}}\sigma_{\mathrm{sm}}^{-1}\| \b f \|_{L^{{r}'}(\Omega,{\mathbb{R}}^d)}\right)^\frac{1}{{r}+1-{r^{\circ}}}.
  \end{aligned}$$

\[rem:ocv\] From , neglecting higher-order terms, we infer asymptotic convergence rates of $\mathcal O_{\mathrm{vel}}^k \coloneqq \frac{(k+1)({r^{\circ}}-1)}{{r}+1-{r^{\circ}}}$ for the velocity and $\mathcal O_{\mathrm{pre}}^k \coloneqq \frac{(k+1)({r^{\circ}}-1)^2}{{r}+1-{r^{\circ}}}$ for the pressure, that is, $$\label{eq:asymptotic.order}
    \mathcal O_{\mathrm{vel}}^k
    = \begin{cases}
      (k+1)(r-1) & \text{if $r<2$}, \\
      \tfrac{k+1}{r-1} & \text{if $r\ge 2$,}
    \end{cases}\quad \text{and} \quad 
    \mathcal O_{\mathrm{pre}}^k
      = \begin{cases}
        (k+1)(r-1)^2 & \text{if $r<2$}, \\
        \tfrac{k+1}{r-1} & \text{if $r\ge 2$.}
      \end{cases}$$ Notice that, owing to the presence of higher-order terms in the right-hand sides of , higher convergence rates may be observed before attaining the asymptotic ones; see Section \[sec:num.res\].

Let $({\und{\b{e}}}_h, \epsilon_h) \coloneqq ({\und{\b{u}}}_h - {\und{\b{\hat u}}}_h,p_h - \hat p_h) \in {{{\und{\b{U}}}}_{h,0}^{k}} \times {P_{h}^{k}}$ where ${\und{\b{\hat u}}}_h \coloneqq{{\und{\b{I}}}_{h}^{k}} \b u \in {{{\und{\b{U}}}}_{h,0}^{k}}$ and $\hat p_h \coloneqq {\pi_{h}^{k}} p \in {P_{h}^{k}}$.\
**Step 1.** *Consistency error.* Let $\mathcal E_h : {{{\und{\b{U}}}}_{h,0}^{k}} \to {\mathbb{R}}$ be the consistency error linear form such that, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}$, $$\label{eq:Eh}
    \mathcal E_h({\und{\b{v}}}_h) \coloneqq \int_\Omega \b f \cdot \b v_h - {\mathrm{a}}_h({\und{\b{\hat u}}}_h,{\und{\b{v}}}_h)-{\mathrm{b}}_h({\und{\b{v}}}_h,{\hat p}_h).$$ Using in the above expression the fact that $\b f = -{\b\nabla{\cdot}}{\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b u)+{\b\nabla}p$ almost everywhere in $\Omega$ to write $\mathcal E_h({\und{\b{v}}}_h)=\mathcal E_{{\mathrm{a}},h}(\b u;{\und{\b{v}}}_h) + \mathcal E_{{\mathrm{b}},h}(p;{\und{\b{v}}}_h)$, and invoking the consistency properties of ${\mathrm{a}}_h$ and of ${\mathrm{b}}_h$, we obtain $$\label{eq:error.estimate:step1:eh0}
    \$\coloneq
    \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}},\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1}\mathcal E_h({\und{\b{v}}}_h)
    \lesssim
    h^{(k+1)({r^{\circ}}-1)} \mathcal N_{{\b\sigma},\b u,p}.$$\
**Step 2.** *Error estimate for the velocity.* Using the strong monotonicity of ${\mathrm{a}}_h$, we get $$\label{eq:error.estimate:step2:eh0}
    \begin{aligned}
      \| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r}+2-{r^{\circ}}}
      &\lesssim \sigma_{\mathrm{sm}}^{-1}\left(
      \sigma_{\mathrm{de}}^{r}+\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\|{\und{\b{\hat u}}}_h\|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{2-{r^{\circ}}}{{r}}\left(
      {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{e}}}_h)-{\mathrm{a}}_h({\und{\b{\hat u}}}_h,{\und{\b{e}}}_h)
      \right) \\
      &\lesssim \sigma_{\mathrm{sm}}^{-1}\mathcal N_{\b f}^{2-{r^{\circ}}}\left(
      {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{e}}}_h)-{\mathrm{a}}_h({\und{\b{\hat u}}}_h,{\und{\b{e}}}_h)
      \right),
    \end{aligned}$$ where we have used the a priori bound on the discrete solution along with the boundedness of the global interpolator and the a priori bound on the continuous solution to conclude. Using then the discrete mass equation along with (written for $\b v=\b u$) and the continuous mass equation to write ${\mathrm{b}}_h({{\und{\b{I}}}_{h}^{k}} \b u,q_h) = b(\b u,q_h) = 0$, we get ${\mathrm{b}}_h({\und{\b{e}}}_h,q_h) = 0$ for all $q_h \in P^k_h$. Hence, combining this result with and the discrete momentum equation (with ${\und{\b{v}}}_h = {\und{\b{e}}}_h$), we obtain $$\label{eq:error.estimate:step1:eh1}
    {\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{e}}}_h)-{\mathrm{a}}_h({\und{\b{\hat u}}}_h,{\und{\b{e}}}_h)
    = \int_\Omega \b f \cdot \b e_h - {\mathrm{a}}_h({\und{\b{\hat u}}}_h,{\und{\b{e}}}_h)-\cancel{{\mathrm{b}}_h({\und{\b{e}}}_h,p_h)}
    = \mathcal E_h({\und{\b{e}}}_h).$$ Plugging into , we get $$\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r}+2-{r^{\circ}}}
  \le \sigma_{\mathrm{sm}}^{-1}\mathcal N_{\b f}^{2-{r^{\circ}}}\$\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}.$$ Simplifying, using , and taking the $({r}+1-{r^{\circ}})$th root of the resulting inequality yields .\
**Step 3.** *Error estimate for the pressure.* Using the Hölder continuity of ${\mathrm{a}}_h$, we have, for all ${\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}}$, $$\label{eq:error.estimate:step3:ah}
   \begin{aligned}
     \left|{\mathrm{a}}_h({\und{\b{\hat u}}}_h,{\und{\b{v}}}_h)-{\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h)\right|
     &\lesssim \sigma_{\mathrm{hc}}\left( \sigma_{\mathrm{de}}^{r}+  \| {\und{\b{\hat u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}+\| {\und{\b{u}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{r}\right)^\frac{{r}-{r^{\circ}}}{{r}}\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r^{\circ}}-1}\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} \\
  &\lesssim \sigma_{\mathrm{hc}}\mathcal N_{\b f}^{{r}-{r^{\circ}}}\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r^{\circ}}-1}\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h},
  \end{aligned}$$ where the first factor is estimated as in . Thus, using the inf-sup condition , we can write $$\label{eq:error.estimate:step3:epsilonh}
    \begin{aligned}
      \| \epsilon_h \|_{L^{{r}'}(\Omega,{\mathbb{R}})} &\lesssim \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}},\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1} {\mathrm{b}}_h({\und{\b{v}}}_h,\epsilon_h) \\
      &= \sup\limits_{{\und{\b{v}}}_h \in {{{\und{\b{U}}}}_{h,0}^{k}},\| {\und{\b{v}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h} = 1} \left(\mathcal E_h({\und{\b{v}}}_h)+{\mathrm{a}}_h({\und{\b{\hat u}}}_h,{\und{\b{v}}}_h)-{\mathrm{a}}_h({\und{\b{u}}}_h,{\und{\b{v}}}_h)\right)\\
      &\lesssim \$+\sigma_{\mathrm{hc}}\mathcal N_{\b f}^{{r}-{r^{\circ}}}\| {\und{\b{e}}}_h \|_{{\boldsymbol{\varepsilon}},{r},h}^{{r^{\circ}}-1} \\
      &\lesssim  h^{(k+1)({r^{\circ}}-1)}\mathcal N_{{\b\sigma},\b u,p}
      +h^{(k+1)({r^{\circ}}-1)^2}\sigma_{\mathrm{hc}}\mathcal N_{\b f}^{|{r}-2|({r^{\circ}}-1)}    
      \left(\sigma_{\mathrm{sm}}^{-1}\mathcal N_{{\b\sigma},\b u,p}\right)^\frac{{r^{\circ}}-1}{{r}+1-{r^{\circ}}},
    \end{aligned}$$ where we have used the definition of the consistency error together with equation to pass to the second line, to pass to the third line (recall that $\$$ denotes here the supremum in the left-hand side of ), and the bounds and (proved in Step 2) to conclude.

Numerical examples {#sec:num.res}
==================

We consider a manufactured solution to problem in order to assess the convergence of the method, which was implemented within the SpaFEDTe library (cf. <https://spafedte.github.io>). Specifically, we take $\Omega=(0,1)^{2}$ and consider the $(1,0,1,{r})$-Carreau–Yasuda law (corresponding to the power-law model) with Sobolev exponent ${r}\in\{1.5, 1.75, 2, 2.25, 2.5, 2.75\}$. The exact velocity $\b u$ and pressure $p$ are given by, respectively, $$\b u(x,y) = \left(\sin\left(\tfrac{\pi}{2}x\right)\cos\left(\tfrac{\pi}{2}y\right),-\cos\left(\tfrac{\pi}{2}x\right)\sin\left(\tfrac{\pi}{2}y\right)\right),\quad 
  p(x,y) = \sin\left(\tfrac{\pi}{2} x\right)\sin\left(\tfrac{\pi}{2} y\right)-\tfrac{4}{\pi^2}.$$ The volumetric load $\b f$ and the Dirichlet boundary conditions are inferred from the exact solution. This solution matches the assumptions required in Theorem \[thm:error.estimate\] for $k = 1$, except the case $r = 1.5$ for which ${\b\sigma}(\cdot,{\b{\nabla}_{\mathrm{s}}}\b u) \notin W^{1,{r}'}(\Omega,{{\mathbb{R}}^{d \times d}_{\mathrm{s}}})$. We consider the HHO scheme for $k=1$ on three mesh families, namely Cartesian orthogonal, distorted triangular, and distorted Cartesian; see Figure \[fig:meshes\]. Overall, the results are in agreement with the theoretical predictions, and in some cases the expected asymptotic orders of convergence are exceeded. Specifically, for ${r}\neq 2$, the convergence rates computed on the last refinement surpass in some cases the theoretical ones. As noticed in Remark \[rem:ocv\], this suggests that the asymptotic order is still not attained. A similar phenomenon has been observed on certain meshes for the $p$-Laplace problem; see [@Di-Pietro.Droniou:17*1 Section 3.5.2] and [@Di-Pietro.Droniou.ea:18 Section 3.7].

![Coarsest Cartesian, distorted triangular, and distorted Cartesian mesh used in the numerical tests of Section \[sec:num.res\].\[fig:meshes\]](square.pdf "fig:") ![Coarsest Cartesian, distorted triangular, and distorted Cartesian mesh used in the numerical tests of Section \[sec:num.res\].\[fig:meshes\]](disttri.pdf "fig:") ![Coarsest Cartesian, distorted triangular, and distorted Cartesian mesh used in the numerical tests of Section \[sec:num.res\].\[fig:meshes\]](distsquare.pdf "fig:")

![Numerical results for the test case of Section \[sec:num.res\]. The slopes indicate the expected order of convergence expected from Theorem \[thm:error.estimate\], i.e. $\mathcal O_{\mathrm{vel}}^1 = 2(r-1)$ and $\mathcal O_{\mathrm{pre}}^1  = 2(r-1)^2$ for ${r}\in \{1.5,1.75,2\}$. \[tab:num.res.1\]](conv1.pdf)

![Numerical results for the test case of Section \[sec:num.res\]. The slopes indicate the expected order of convergence expected from Theorem \[thm:error.estimate\], i.e. $\mathcal O_{\mathrm{vel}}^1 = \mathcal O_{\mathrm{pre}}^1 = \frac{2}{r-1}$ for ${r}\in \{2.25,2.5,2.75\}$. \[tab:num.res.2\]](conv2.pdf)

Power-framed functions {#sec:properties.stress}
======================

In the following theorem, we introduce the notion of power-framed function and discuss sufficient conditions for this property to hold.

\[thm:1d.power-framed\] Let $U$ be a measurable subset of ${\mathbb{R}}^n$ with $n\ge1$, $(W,(\cdot,\cdot)_W)$ an inner product space, and $\b \sigma : U \times W \to W$. Assume that there exists a Carathéodory function $\varsigma : U \times \lbrack0,+\infty) \to {\mathbb{R}}$ such that, for all $\b\tau \in W$ and almost every $\b x \in U$,

\[eq:1d.power-framed:stress\] $${\b\sigma}(\b x,\b\tau) = \varsigma(\b x,\|\b\tau\|_W)\b\tau,$$ where $\|{\cdot}\|_W$ is the norm induced by $(\cdot,\cdot)_W$. Additionally assume that, for almost every $\b x \in U$, $\varsigma(\b x,\cdot)$ is differentiable on $(0,+\infty)$ and there exist $\varsigma_{\mathrm{de}} \in \lbrack0,+\infty)$ and $\varsigma_{\mathrm{sm}},\varsigma_{\mathrm{hc}} \in (0,+\infty)$ independent of $\b x$ such that, for all $\alpha \in (0,+\infty)$, $$\begin{aligned}
      \varsigma_{\mathrm{sm}} (\varsigma_{\mathrm{de}}^{r}+\alpha^{r})^\frac{{r}-2}{{r}} \leq \frac{\partial(\alpha\varsigma(\b x,\alpha))}{\partial \alpha} \leq \varsigma_{\mathrm{hc}}(\varsigma_{\mathrm{de}}^{r}+\alpha^{r})^\frac{{r}-2}{{r}}. \label{eq:1d.power-framed:eta}
    \end{aligned}$$

Then, ${\b\sigma}$ is an *${r}$-power-framed function*, i.e., for all $(\b\tau,\b\eta) \in W^2$ with $\b\tau \neq \b\eta$ and almost every $\b x \in U$, the function ${\b\sigma}$ verifies the Hölder continuity property

\[eq:od:power-framed:holder.continuity.strong.monotonicity\] $$\|{\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta)\|_W \le \sigma_{\mathrm{hc}} \left(\sigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r}\right)^\frac{{r}-2}{{r}}\| \b\tau-\b\eta \|_W,\label{eq:od:power-framed:holder.continuity}$$ and the strong monotonicity property $$\left({\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta),\b\tau-\b\eta\right)_W \ge \sigma_{\mathrm{sm}}\left(\sigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r}\right)^\frac{{r}-2}{{r}}\|\b\tau-\b\eta\|_W^{2},\label{eq:od:power-framed:strong.monotonicity}$$

with $\sigma_{\mathrm{de}} \coloneqq \varsigma_{\mathrm{de}}$, $\sigma_{\mathrm{hc}} \coloneqq 2^{2-{r^{\circ}}+{r}^{-1}\left\lceil\hspace{0.02cm}2-{r^{\circ}}\right\rceil}({r^{\circ}}-1)^{-1} \varsigma_{\mathrm{hc}}$, and $\sigma_{\mathrm{sm}} \coloneqq 2^{{r^{\circ}}-{r}-\left\lceil{r}^{-1}({r}-{r^{\circ}})\right\rceil}({r}+1-{r^{\circ}})^{-1} \varsigma_{\mathrm{sm}}$, where ${r^{\circ}}$ is given by and $\lceil{\cdot}\rceil$ is the ceiling function.

The boldface notation for the elements of $W$ is reminescent of the fact that Theorem \[thm:1d.power-framed\] is used with $W = {{\mathbb{R}}^{d \times d}_{\mathrm{s}}}$ in Corollary \[cor:Carreau–Yasuda\] to characterize the Carreau-Yasuda law as an ${r}$-power-framed function and in Lemma \[lem:sT\] with $W = {\mathbb{R}}^d$ to study the local stabilization function ${\mathrm{s}}_T$.

Let $\b x \in U$ be such that holds, and $\b\tau,\b\eta \in W$. By symmetry of inequalities and the fact that ${\b\sigma}$ is continuous, we can assume, without loss of generality, that $\|\b\tau\|_W > \|\b\eta\|_W > 0$.\
(i) *Strong monotonicity.* Let $\beta \in (0,+\infty)$ and let $g : \lbrack\beta,+\infty) \to {\mathbb{R}}$ be such that, for all $\alpha \in \lbrack\beta,+\infty)$, $$g(\alpha) \coloneqq \alpha\varsigma(\b x,\alpha)-\beta\varsigma(\b x,\beta)-C_{\mathrm{sm}}(\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r})^\frac{{r}-2}{{r}}(\alpha-\beta),
  \;\text{ with }\ C_{\mathrm{sm}} \coloneqq \tfrac{2^{{r^{\circ}}-{r}}}{{r}+1-{r^{\circ}}}\varsigma_{\mathrm{sm}}.$$ Differentiating $g$ and using the first inequality in , we obtain, for all $\alpha \in \lbrack\beta,+\infty)$, $$\begin{aligned}
    \frac{\partial}{\partial \alpha}g(\alpha)
    &\geq \varsigma_{\mathrm{sm}}(\varsigma_{\mathrm{de}}^{r}+\alpha^{r})^\frac{{r}-2}{{r}}-C_{\mathrm{sm}}\left(({r}-2)(\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r})^{-\frac{2}{{r}}}(\alpha-\beta)\alpha^{{r}-1}+(\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r})^\frac{{r}-2}{{r}}\right) \\
    &\geq \varsigma_{\mathrm{sm}}(\varsigma_{\mathrm{de}}^{r}+\alpha^{r})^\frac{{r}-2}{{r}}-({r}+1-{r^{\circ}})C_{\mathrm{sm}} (\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r})^\frac{{r}-2}{{r}} \\
    &\geq \varsigma_{\mathrm{sm}}2^{{r^{\circ}}-{r}}(\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r})^\frac{{r}-2}{{r}}-({r}+1-{r^{\circ}})C_{\mathrm{sm}} (\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r})^\frac{{r}-2}{{r}} = 0,
  \end{aligned}$$ where, to pass to the second line, we have removed negative contributions if ${r}< 2$ and used the fact that $(\alpha-\beta)\alpha^{{r}-1} \le \varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r}$ if ${r}\ge 2$, to pass to the third line we have used the fact that $t \mapsto t^{{r}-2}$ is non-increasing if ${r}< 2$, and the fact that $\beta \le \alpha$ otherwise, while the conclusion follows from the definition of $C_{\mathrm{sm}}$. This shows that $g$ is non-decreasing. Hence, for all $\alpha\in\lbrack\beta,+\infty)$, $g(\alpha)\geq g(\beta)=0$, i.e. $$\label{eq:1d.power-framed:bound:1}
    \alpha\varsigma(\b x,\alpha)-\beta\varsigma(\b x,\beta) \geq C_{\mathrm{sm}}(\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r})^\frac{{r}-2}{{r}}(\alpha-\beta).$$ Moreover, for all $\alpha,\beta \in (0,+\infty)$, using (with $\beta = 0$) along with the fact that $t \mapsto t^{{r}-2}$ is decreasing if ${r}< 2$ and inequality if ${r}\ge 2$, we infer that $$\label{eq:1d.power-framed:bound:2}
    \begin{aligned}
      \varsigma(\b x,\alpha)+\varsigma(\b x,\beta) &\ge C_{\mathrm{sm}}\left((\varsigma_{\mathrm{de}}^{r}+\alpha^{r})^\frac{{r}-2}{{r}}+(\varsigma_{\mathrm{de}}^{r}+\beta^{r})^\frac{{r}-2}{{r}}\right) \ge C_{\mathrm{sm}}2^{1-\left\lceil\frac{{r}-{r^{\circ}}}{{r}}\right\rceil}(\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r})^\frac{{r}-2}{{r}}.
    \end{aligned}$$ We conclude that ${\b\sigma}$ verifies by using and with $\alpha = \|\b\tau\|_W$ and $\beta = \|\b\eta\|_W$ as follows: $$\begin{aligned}
    &({\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta),\b\tau-\b\eta)_W
    \\
    &\quad= (\b\tau\varsigma(\b x,\|\b\tau\|_W)-\b\eta\varsigma(\b x,\|\b\eta\|_W),\b\tau-\b\eta)_W \\
    &\quad = \|\b\tau\|_W^2\varsigma(\b x,\|\b\tau\|_W)+\|\b\eta\|_W^2\varsigma(\b x,\|\b\eta\|_W) 
    -(\b\tau,\b\eta)_W \left[\varsigma(\b x,\|\b\tau\|_W) + \varsigma(\b x,\|\b\eta\|_W)\right] \\
    &\quad = \left[
      \|\b\tau\|_W\varsigma(\b x,\|\b\tau\|_W)-\|\b\eta\|_W\varsigma(\b x,\|\b\eta\|_W)
      \right](\|\b\tau\|_W-\|\b\eta\|_W)
    \\
    &\qquad
    + \left[
      \varsigma(\b x,\|\b\tau\|_W)+\varsigma(\b x,\|\b\eta\|_W)
      \right](\|\b\tau\|_W\|\b\eta\|_W-(\b\tau,\b\eta)_W)
    \\
    &\quad\geq  C_{\mathrm{sm}}2^{-\left\lceil\frac{{r}-{r^{\circ}}}{{r}}\right\rceil}\left(
    \varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r}\right)^\frac{{r}-2}{{r}}\left[
      (\|\b\tau\|_W-\|\b\eta\|_W)^2+2(\|\b\tau\|_W\|\b\eta\|_W-(\b\tau,\b\eta)_W)
      \right] \\
    &\quad= C_{\mathrm{sm}}2^{-\left\lceil\frac{{r}-{r^{\circ}}}{{r}}\right\rceil}\left(
    \varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r}\right)^\frac{{r}-2}{{r}}\|\b\tau-\b\eta\|_W^2.
  \end{aligned}$$\
(ii) *Hölder continuity.* Now, setting $C_{\mathrm{hc}} \coloneqq \frac{\varsigma_{\mathrm{hc}}}{{r^{\circ}}-1}$ and reasoning in a similar way as for the proof of to leverage the second inequality in , we have, for all $\alpha \in \lbrack\beta,+\infty)$, $$\label{eq:1d.power-framed:bound:3}
    \alpha\varsigma(\b x,\alpha)-\beta\varsigma(\b x,\beta) \le C_{\mathrm{hc}}\left(
    \varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r}\right)^\frac{{r}-2}{{r}}(\alpha-\beta).$$ First, let ${r}\ge 2$. Using (with $\beta=0$) and the fact that $t \mapsto t^{{r}-2}$ is non-decreasing, we have, for all $\alpha,\beta \in (0,+\infty)$, $$\label{eq:1d.power-framed:bound:4}
      \varsigma(\b x,\alpha)\varsigma(\b x,\beta)
      \le C_{\mathrm{hc}}^2\left(
      \varsigma_{\mathrm{de}}^{r}+\alpha^{r}\right)^\frac{{r}-2}{{r}}\left(
      \varsigma_{\mathrm{de}}^{r}+\beta^{r}\right)^\frac{{r}-2}{{r}} \le \left[
        C_{\mathrm{hc}}\left(\varsigma_{\mathrm{de}}^{r}+\alpha^{r}+\beta^{r}\right)^\frac{{r}-2}{{r}}
        \right]^2.$$ Thus, using inequalities and with $\alpha = \|\b\tau\|_W$ and $\beta = \|\b\eta\|_W$, we infer $$\label{eq:1d.power-framed:bound:5}
    \begin{aligned}
      &\|{\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta)\|_W^2
      \\
      &\quad = \left(
      \b\tau\varsigma(\b x,\|\b\tau\|_W)-\b\eta\varsigma(\b x,\|\b\eta\|_W),
      \b\tau\varsigma(\b x,\|\b\tau\|_W)-\b\eta\varsigma(\b x,\|\b\eta\|_W)
      \right)_W
      \\
      &\quad = \left[\|\b\tau\|_W\varsigma(\b x,\|\b\tau\|_W)-\|\b\eta\|_W\varsigma(\b x,\|\b\eta\|_W)\right]^2
      \\
      &\quad \qquad + 2\varsigma(\b x,\|\b\tau\|_W)\varsigma(\b x,\|\b\eta\|_W)\left[
        \|\b\tau\|_W\|\b\eta\|_W-(\b\tau,\b\eta)_W
        \right]
      \\
      &\quad \le  \left[
        C_{\mathrm{hc}}\left(\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r}\right)^\frac{{r}-2}{{r}}
        \right]^2\left[
        (\|\b\tau\|_W-\|\b\eta\|_W)^2+2(\|\b\tau\|_W\|\b\eta\|_W-(\b\tau,\b\eta)_W)
        \right] \\
      &\quad = \left[
        C_{\mathrm{hc}}\left(\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r}\right)^\frac{{r}-2}{{r}}\|\b\tau-\b\eta\|_W
        \right]^2,
    \end{aligned}$$ hence ${\b\sigma}$ verifies for ${r}\ge 2$. Assume now ${r}< 2$. Using a triangle inequality followed by and the left inequality in , it is inferred that $$\begin{aligned}
    \|{\b\sigma}(\b x,\b\tau)-{\b\sigma}(\b x,\b\eta)\|_W 
    &\leq \varsigma(\b x,\|\b\tau\|_W)\|\b\tau\|_W+\varsigma(\b x,\|\b\eta\|_W)\|\b\eta\|_W 
    \\
    & \leq C_{\mathrm{hc}}\left((\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r})^\frac{{r}-1}{{r}}+(\varsigma_{\mathrm{de}}^{r}+\|\b\eta\|_W^{r})^\frac{{r}-1}{{r}}\right) 
    \\
    &\leq 2^\frac{1}{{r}}C_{\mathrm{hc}}(2\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r})^\frac{{r}-1}{{r}} 
    \\
    & =  2^\frac{1}{{r}} C_{\mathrm{hc}}(2\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r})^\frac{{r}-2}{{r}}
    (2\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r})^\frac1{r},
    \\
    & \leq  2^\frac{1}{{r}} C_{\mathrm{hc}}(\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r})^\frac{{r}-2}{{r}}
    (2\varsigma_{\mathrm{de}}+\|\b\tau\|_W+\|\b\eta\|_W),
  \end{aligned}$$ where the last line follows from the fact that $t \mapsto t^{{r}-2}$ is decreasing and again . If $2\varsigma_{\mathrm{de}}+\|\b\tau\|_W+\|\b\eta\|_W\le 2^{2-r}\|\b\tau-\b\eta\|_W$, from the previous bound we directly get the conclusion, i.e. with $\sigma_{\mathrm{hc}}=2^{2-r+\frac{1}{r}}C_{\mathrm{hc}}$. Otherwise, using and a triangle inequality yields $$\label{eq:est_else}
  \begin{aligned}
    (\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r})^\frac{1}{{r}}(\varsigma_{\mathrm{de}}^{r}+\|\b\eta\|_W^{r})^\frac{1}{{r}} 
    &\ge 2^{-\frac{2}{r'}}(\varsigma_{\mathrm{de}}+\|\b\tau\|_W)(\varsigma_{\mathrm{de}}+\|\b\eta\|_W) \\
    &= 2^{-2(\frac{1}{r'}+1)}\left[
      \left(2\varsigma_{\mathrm{de}}+\|\b\tau\|_W+\|\b\eta\|_W\right)^2-\left(\|\b\tau\|_W-\|\b\eta\|_W\right)^2
      \right]
    \\
    &\ge 2^{-2(\frac{1}{r'}+1)}\left[
      \left(2\varsigma_{\mathrm{de}}+\|\b\tau\|_W+\|\b\eta\|_W\right)^2-\|\b\tau-\b\eta\|_W^2
      \right]
    \\
    &\ge 2^{-2(\frac{1}{r'}+1)}(1-4^{r-2})\left(2\varsigma_{\mathrm{de}}+\|\b\tau\|_W+\|\b\eta\|_W\right)^2\\
    &\ge 2^{\frac{2}{(r-2)r}-2}\left(\varsigma_{\mathrm{de}}^{r}+\|\b\tau\|_W^{r}+\|\b\eta\|_W^{r}\right)^\frac{2}{{r}},
    \end{aligned}$$ where we concluded with together with the fact that $2^{-2(\frac{1}{r'}+1)}\left(1-4^{r-2}\right) \ge 2^{\frac{2}{(r-2)r}-2}$. Finally, raising both sides of to the power ${r}-2$, we get a relation analogous to . Hence, proceeding as in , we infer .

\[cor:Carreau–Yasuda\] The strain rate-shear stress law of the $(\mu,\delta,a,{r})$-Carreau–Yasuda fluid defined in Example \[ex:Carreau–Yasuda\] is an ${r}$-power-framed function.

Let $\b x \in \Omega$ and $g : (0,+\infty) \to {\mathbb{R}}$ be such that, for all $\alpha\in(0,+\infty)$, $$g(\alpha)
  \coloneqq \frac{\partial}{\partial \alpha}\left[
  \alpha\mu(\b x)\left(\delta^{a(\b x)}+\alpha^{a(\b x)}
  \right)^\frac{{r}-2}{a(\b x)}
  \right]
  = \mu(\b x)\left(
  \delta^{a(\b x)}+\alpha^{a(\b x)}
  \right)^{\frac{{r}-2}{a(\b x)}-1}\left(
  \delta^{a(\b x)}+({r}-1)\alpha^{a(\b x)}
  \right).$$ We have for all $\alpha \in (0,+\infty)$, $$\mu_- ({r^{\circ}}-1) \left(
  \delta^{a(\b x)}+\alpha^{a(\b x)}
  \right)^{\frac{{r}-2}{a(\b x)}} \le g(\alpha)
  \leq \mu_+ ({r}+1-{r^{\circ}})\left(
  \delta^{a(\b x)}+\alpha^{a(\b x)}
  \right)^{\frac{{r}-2}{a(\b x)}},$$ and we conclude using together with Theorem \[thm:1d.power-framed\].

[^1]: [<michele.botti@polimi.it>]{}

[^2]: [<danielcq.mathematics@gmail.com>]{}

[^3]: [<daniele.di-pietro@umontpellier.fr>]{}

[^4]: [<andre.harnist@umontpellier.fr>]{}, corresponding author
