$\frac{1}{s^{\rho}}-\frac{1}{t^{\rho}}\leq c(t-s)^{\theta}s^{-(\rho+\theta)}.$

Let $\rho\geq 0$ and $\theta,\theta^{\prime}\in[0,1]$. There exists a constant
such that

$\frac{1}{(\tau^{\prime}-u)^{\rho}}-\frac{1}{(\tau^{\prime}-s)^{\rho}}-\frac{1}{(\tau-u)^{\rho}}+\frac{1}{(\tau-s)^{\rho}}\leq
c(\tau-\tau^{\prime})^{\theta}(u-s)^{\theta^{\prime}}(\tau^{\prime}-u)^{-(\rho+\theta+\theta^{\prime})}$

###### Proof.

Let us remark that

$\frac{1}{s^{\rho}}-\frac{1}{t^{\rho}}=(1+\rho)\int_{s}^{t}r^{-(\rho+1)}\mathop{}\\!\mathrm{d}r\leq(1+\rho)(t-s)s^{-(\rho+1)}.$

The result follows by a standard interpolation between the two inequalities.
For the second inequality, set

$B=\frac{1}{(\tau^{\prime}-u)^{\rho}}-\frac{1}{(\tau^{\prime}-s)^{\rho}}-\frac{1}{(\tau-u)^{\rho}}+\frac{1}{(\tau-s)^{\rho}}.$

We have, thanks to the same integral representation,

$B\lesssim(\tau^{\prime}-u)^{-\rho},$
$B\lesssim(\tau-\tau^{\prime})(\tau^{\prime}-u)^{-(\rho+1)},$
$B\lesssim(u-s)(\tau^{\prime}-u)^{-(\rho+1)},$

and

$B\lesssim(u-s)(\tau-\tau^{\prime})(\tau^{\prime}-u)^{-(\rho+2)}.$

Let us suppose, without loss of generality that $\theta^{\prime}\leq\theta$
that is $\theta^{\prime}=\alpha\theta$ for some $0\leq\alpha<1$. We have,
using the first and the last inequalities,

$B\lesssim(\tau-\tau^{\prime})^{\theta}(u-s)^{\theta}(\tau^{\prime}-u)^{-(\rho+2\theta)}.$

Using the first and the second one, we have

$B\lesssim(\tau-\tau^{\prime})^{\theta}(\tau^{\prime}-u)^{-(\rho+2\theta)}.$

We interpolate those last two inequalities to have

$B\lesssim(\tau-\tau^{\prime})^{\theta}(u-s)^{\alpha\theta}(\tau^{\prime}-u)^{-(\rho+2\alpha\theta+(1-\alpha)\theta)}=(\tau-\tau^{\prime})^{\theta}(u-s)^{\theta^{\prime}}(\tau^{\prime}-u)^{-(\rho+\theta+\theta^{\prime})}.$

∎

######  Corollary 45.

Let $w$ be an admissible weight and $\alpha\in(0,2]$. Let
$\kappa\in\mathbb{R}$ and let $1\leq p,q\leq+\infty$. Let $0\leq t$, then for
all $\rho\geq 0$,

$\|P^{\frac{\alpha}{2}}_{t}f\|_{B^{\kappa+\alpha\rho}_{p,q}(w)}\lesssim
t^{-\rho}\|f\|_{B^{\kappa}_{p,q}(w)}.$

For all $\theta\in[0,1]$ and all $\rho>0$,

$\|(P^{\frac{\alpha}{2}}_{t}-P^{\frac{\alpha}{2}}_{s})f\|_{B^{\kappa+\alpha\rho}_{p,q}(w)}\lesssim|(t-s)^{\theta}s^{-(\rho+\theta)}\|f\|_{B^{\kappa}_{p,q}(w)}.$

Furthermore, for any $\theta,\theta^{\prime}\in[0,1]$ and all $\rho\geq 0$,
$s\leq u<\tau^{\prime}\leq\tau$,

$\left\|\big{(}(P^{\frac{\alpha}{2}}_{\tau-s}-P^{\frac{\alpha}{2}}_{\tau-u})-(P^{\frac{\alpha}{2}}_{\tau^{\prime}-s}-P^{\frac{\alpha}{2}}_{\tau^{\prime}-u})\big{)}f\right\|_{B^{\kappa+\alpha\rho}_{p,q}(w)}\lesssim(\tau-\tau^{\prime})^{\theta}(u-s)^{\theta^{\prime}}(\tau^{\prime}-u)^{-(\rho+\theta+\theta^{\prime})}\|f\|_{B^{\kappa}_{p,q}(w)}.$

###### Proof.

The first bound is a direct consequence of the first bound of Proposition 43.
Indeed, for all $j\geq 0$,

$\|\Delta_{j}P^{\frac{\alpha}{2}}_{t}f\|_{L^{P}(w)}=\|P_{t}(\Delta_{j}f)\|_{L^{P}(w)}\lesssim
e^{-c2^{\alpha j}t}\|\Delta_{j}f\|_{L^{P}(w)}\lesssim 2^{-\alpha\rho
j}t^{-\rho}\|\Delta_{j}f\|_{L^{P}(w)},$

and the result follows via the definition of weighted Besov spaces. For the
second bound, let us remark that we have for all $0<s\leq t$, and thanks to
Proposition 43,

$\|\Delta_{j}(P^{\frac{\alpha}{2}}_{t}-P^{\frac{\alpha}{2}}_{s})f\|_{L^{P}(w)}=\lesssim\left(\frac{1}{s^{\rho}}-\frac{1}{s^{\rho}}\right)2^{-\alpha\rho
j}\|\Delta_{j}f\|_{L^{P}(w)}\lesssim 2^{-\alpha\rho
j}(t-s)^{\theta}s^{-(\rho+\theta)}\|\Delta_{j}f\|_{L^{p}(w)}.$

where the last inequality comes from the previous lemma. This allows us to
derive the result by using the definition of $B$ spaces. Finally, note that we
have thanks to Proposition 43 and Lemma 44,

$\Big{\|}\Delta_{j}\big{(}(P^{\frac{\alpha}{2}}_{\tau-u}-P^{\frac{\alpha}{2}}_{\tau-s})-(P^{\frac{\alpha}{2}}_{\tau^{\prime}-u}-P^{\frac{\alpha}{2}}_{\tau^{\prime}-s})\big{)}f\Big{\|}_{L^{p}(w)}\lesssim
2^{-\alpha\rho
j}(\tau-\tau^{\prime})^{\theta}(u-s)^{\theta^{\prime}}(\tau^{\prime}-u)^{\rho+\theta+\theta^{\prime}}\|\Delta_{j}f\|_{L^{p}(w)}.$

And again one can conclude with the definition of the weighted Besov spaces. ∎

##  Appendix B Cauchy-Lipschitz theorem for mSHE in standard case.

We give a short proof of local well-posedness of the mSHE in a simple context.
For more on this, one can consult [30] and the generalization for more
irregular noise with more involve techniques [19] and [21]. Note that

######  Theorem 46.

Let $\alpha\in(0,2]$. Let $\vartheta\in(0,\min(1,\alpha))$, let
$\vartheta<\beta<\alpha-\vartheta$. Let $g\in\mathcal{C}^{2}$ and let
$u_{0}\in\mathcal{C}^{\beta}$ and $\xi\in\mathcal{C}^{-\vartheta}$ There
exists a unique local solution of the equation

$\partial_{t}u=\Delta^{\frac{\alpha}{2}}u+g(u)\xi$

in the mild form

$u(t,\cdot)=P^{\frac{\alpha}{2}}_{t}u_{0}+\int_{s}^{t}P^{\frac{\alpha}{2}}_{t-s}\xi
g(u(s,\cdot))\mathop{}\\!\mathrm{d}s.$

###### Proof.

First, let us take $u\in\mathcal{C}^{\beta}$. We have

$\|g(u(t,\cdot)\|_{\mathcal{C}^{\beta}}\leq\|g\|_{\mathcal{C}^{2}}\|u(t,.)\|_{\mathcal{C}^{\beta}}.$

For $u\in C([0,T];\mathcal{C}^{\beta})$, let us remark that thanks to standard
Bony estimates in Besov-Hölder spaces (Proposition 40)

$\|\xi
g(u(t,\cdot)\|_{\mathcal{C}^{-\vartheta}}\lesssim\|g\|_{\mathcal{C}^{2}}\|\xi\|_{\mathcal{C}^{-\vartheta}}\|u(t,.)\|_{\mathcal{C}^{\beta}}.$

Finally for $0\leq s<t\leq T$,

$\|P_{t-s}\xi
g(u(t,\cdot))\|_{\mathcal{C}^{\beta}}\lesssim\frac{1}{(t-s)^{\frac{\beta+\rho}{\alpha}}}\|g\|_{\mathcal{C}^{2}}\|\xi\|_{\mathcal{C}^{-\vartheta}}\|u(t,.)\|_{\mathcal{C}^{\beta}}.$

Hence, the application $\Gamma$

$\Gamma(u)(t,x)=P_{t}u_{0}+\int_{0}^{t}P_{t-s}\xi
g(u(s,\cdot))\mathop{}\\!\mathrm{d}s$

is well-defined from $C^{0}([0,T];\mathcal{C}^{\beta})$ to itself.
Furthermore, Let us remark that for $u,v\in\mathcal{C}^{\beta}$,

$\|g(u)-g(v)\|_{\mathcal{C}^{\beta}}\lesssim\|g\|_{\mathcal{C}^{2}}\|u-v\|_{\mathcal{C}^{\beta}},$

hence for $u,v\in C^{0}([0,T];\mathcal{C}^{\beta})$,

$\sup_{t\in[0,T]}\|\Gamma(u)(t,\cdot)-\Gamma(v)(t,\cdot)\|_{\beta}\lesssim\int_{0}^{t}\frac{1}{(t-s)^{\frac{\beta+\rho}{\alpha}}}\|g\|_{\mathcal{C}^{2}}\|\xi\|_{\mathcal{C}^{-\vartheta}}\|u(t,.)-v(t,.)\|_{\mathcal{C}^{\beta}}\mathop{}\\!\mathrm{d}s\\\
\lesssim
T^{1-\frac{\beta+\rho}{\alpha}}\sup_{t\in[0,T]}\|u(t,.)-v(t,.)\|_{\mathcal{C}^{\beta}}.$

Hence, be standard Schauder fixed point, for $T$ small enough, there is a
unique $u\in C^{0}([0,T];\mathcal{C}^{\beta})$. ∎

## References

  * [1] S. Athreya, O. Butkovsky, K. Lê, and L. Mytnik.  Well-posedness of stochastic heat equation with distributional drift and skew stochastic heat equation.  arXiv:2011.13498 [math], November 2020.  arXiv: 2011.13498.
  * [2] H. Bahouri, J-Y. Chemin, and R. Danchin.  Fourier analysis and nonlinear partial differential equations, volume 343 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer, Heidelberg, 2011.
  * [3] I. Bailleul, A. Debussche, and M. Hofmanová.  Quasilinear generalized parabolic Anderson model equation.  Stoch. Partial Differ. Equ. Anal. Comput., 7(1):40–63, 2019.
  * [4] C. Bellingeri, P. Friz, and M. Gerencsér.  Singular paths spaces and applications, 2020.
  * [5] J. Benedikt, V. Bobkov, P. Girg, L. Kotrla, and P. Takac.  Nonuniqueness of solutions of initial-value problems for parabolic p-Laplacian.  Electron. J. Differ. Equ, 38:2015, 2015.
  * [6] O. Butkovsky and L. Mytnik.  Regularization by noise and flows of solutions for a stochastic heat equation.  Ann. Probab., 47(1):165–212, 01 2019.
  * [7] R. Catellier and M. Gubinelli.  Averaging along irregular curves and regularisation of ODEs.  Stochastic Process. Appl., 126(8):2323–2366, 2016.
  * [8] L. Coutin, R. Duboscq, and A. Réveillac.  The Itô-Tanaka Trick: a non-semimartingale approach, July 2019\.
  * [9] A. M. Davie.  Differential equations driven by rough paths: an approach via discrete approximation.  Appl. Math. Res. Express. AMRX, (2):Art. ID abm009, 40, 2007.
  * [10] A. Deya and S. Tindel.  Rough Volterra equations. I. The algebraic integration setting.  Stoch. Dyn., 9(3):437–477, 2009.
  * [11] P. Friz.  Mini-course on Rough Paths, 2009.
  * [12] P. Friz and M. Hairer.  A course on rough paths.  Universitext. Springer, Cham, 2014.  With an introduction to regularity structures.
  * [13] H. Fujita and S. Watanabe.  On the uniqueness and non-uniqueness of solutions of initial value problems for some quasi-linear parabolic equations.  Communications on Pure and Applied Mathematics, 21(6):631–652, 1968\.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.3160210609.
  * [14] L. Galeati.  Nonlinear young differential equations: a review, 2020.
  * [15] L. Galeati and M. Gubinelli.  Noiseless regularisation by noise, 2020.
  * [16] L. Galeati and M. Gubinelli.  Prevalence of rho-irregularity and related properties, 2020.
  * [17] L. Galeati and F. A. Harang.  Regularization of multiplicative sdes through additive noise, 2020. Arxiv.
  * [18] M. Gubinelli, P. Imkeller, and N. Perkowski.  Paracontrolled distributions and singular PDEs.  Forum of Mathematics, Pi, 3(e6), 2015.
  * [19] M. Gubinelli, P. Imkeller, and N. Perkowski.  A Fourier analytic approach to pathwise stochastic integration.  Electron. J. Probab., 21, 2016.
  * [20] M. Hairer and X. M. Li.  Averaging dynamics driven by fractional Brownian motion.  Ann. Probab., 48(4):1826–1860, 07 2020.
  * [21] M. Hairer and N. S. Pillai.  Regularity of laws and ergodicity of hypoelliptic SDEs driven by rough paths.  Ann. Probab., 41(4):2544–2598, 2013.
  * [22] F. A. Harang and C. Ling.  Regularity of local times associated to volterra-lévy processes and path-wise regularization of stochastic differential equations, 2020. Arxiv.
  * [23] F. A. Harang and N. Perkowski.  C-infinity regularization of odes perturbed by noise, 2020. Arxiv.
  * [24] F. A. Harang and S. Tindel.  Volterra equations driven by rough signals, 2019. Arxiv.
  * [25] Y. Hu, D. Nualart, and J. Song.  A nonlinear stochastic heat equation: Hölder continuity and smoothness of the density of the solution.  Stochastic Processes and their Applications, 123(3):1083 – 1103, 2013.
  * [26] T. Hytönen, J. Neerven, M. Veraar, and L. Weis.  Analysis in Banach Spaces: Volume I: Martingales and Littlewood-Paley Theory.  Springer, 1st ed. 2016 edition, November 2016.
  * [27] M. Kwasnicki.  Ten equivalent definitions of the fractional laplace operator.  Fractional Calculus and Applied Analysis, 20(1):7–51, February 2017\.  Publisher: De Gruyter Section: Fractional Calculus and Applied Analysis.
  * [28] K. Lê.  A stochastic sewing lemma and applications.  arXiv:1810.10500 [math], March 2020.  arXiv: 1810.10500.
  * [29] M. Lifshits and T. Simon.  Small deviations for fractional stable processes.  Annales de l’Institut Henri Poincare (B) Probability and Statistics, 41(4):725–752, July 2005.
  * [30] W. Liu and M. Röckner.  Stochastic Partial Differential Equations: An Introduction.  Universitext. Springer International Publishing, 2015.
  * [31] J. C. Mourrat and H. Weber.  Global well-posedness of the dynamic $\phi^{4}$ model in the plane.  Ann. Probab., 45(4):2398–2476, 2017.
  * [32] C. Mueller, L. Mytnik, and E. Perkins.  Nonuniqueness for a parabolic SPDE with $\frac{3}{4}-\varepsilon $-Hölder diffusion coefficients.  Annals of Probability, 42(5):2032–2112, September 2014.  Publisher: Institute of Mathematical Statistics.
  * [33] E. Neuman.  Pathwise uniqueness of the stochastic heat equation with spatially inhomogeneous white noise.  Annals of Probability, 46(6):3090–3187, November 2018.  Publisher: Institute of Mathematical Statistics.
  * [34] D. Nualart and Y. Ouknine.  Regularization of quasilinear heat equations by a fractional noise.  Stoch. Dyn., 4(2):201–221, 2004.
  * [35] I. Pinelis.  Optimum Bounds for the Distributions of Martingales in Banach Spaces.  Annals of Probability, 22(4):1679–1706, October 1994.  Publisher: Institute of Mathematical Statistics.
  * [36] G. Samorodnitsky and M. S. Taqqu.  Stable non-Gaussian random processes.  Stochastic Modeling. Chapman & Hall, New York, 1994.  Stochastic models with infinite variance.
  * [37] H. Triebel.  Theory of function spaces. III, volume 100 of Monographs in Mathematics.  Birkhäuser Verlag, Basel, 2006.
  * [38] X. Yang and X. Zhou.  Pathwise uniqueness for an SPDE with Hölder continuous coefficient driven by $\alpha $-stable noise.  Electronic Journal of Probability, 22, 2017.  Publisher: The Institute of Mathematical Statistics and the Bernoulli Society.


remarkRemark Vorticity formulation for the Oseen problem with variable
viscosityV. Anaya, R. Caraballo, B. Gómez-Vargas, D. Mora, R. Ruiz-Baier

# Velocity-vorticity-pressure formulation for the Oseen problem with variable
viscosity††thanks: Updated: . This work has been partially supported by DIUBB
through projects 2020127 IF/R and 194608 GI/C, by the National Agency for
Research and Development, ANID-Chile, through project Centro de Modelamiento
Matemático (AFB170001) of the PIA Program: Concurso Apoyo a Centros
Científicos y Tecnológicos de Excelencia con Financiamiento Basal, by the
Monash Mathematics Research Fund S05802-3951284, and by the HPC-Europa3
Transnational Access programme through grant HPC175QA9K.

V. Anaya GIMNAP, Departamento de Matemática, Universidad del Bío-Bío,
Concepción, Chile and CI2MA, Universidad de Concepción, Concepción, Chile ().
<EMAIL_ADDRESS>R. Caraballo GIMNAP, Departamento de Matemática, Universidad
del Bío-Bío, Concepción, Chile<EMAIL_ADDRESS>B.
Gómez-Vargas Sección de Matemática, Sede de Occidente, Universidad de Costa
Rica, San Ramón de Alajuela, Costa Rica<EMAIL_ADDRESS>D.
Mora GIMNAP, Departamento de Matemática, Universidad del Bío-Bío, Concepción,
Chile and CI2MA, Universidad de Concepción, Concepción, Chile ().
<EMAIL_ADDRESS>R. Ruiz-Baier School of Mathematical Sciences, Monash
University, 9 Rainforest Walk, Melbourne VIC 3800, Australia; and Universidad
Adventista de Chile, Casilla 7-D, Chillán, Chile ().
<EMAIL_ADDRESS>
###### Abstract

We propose and analyse an augmented mixed finite element method for the Oseen
equations written in terms of velocity, vorticity, and pressure with non-
constant viscosity and homogeneous Dirichlet boundary condition for the
velocity. The weak formulation includes least-squares terms arising from the
constitutive equation and from the incompressibility condition, and we show
that it satisfies the hypotheses of the Babuška-Brezzi theory. Repeating the
arguments of the continuous analysis, the stability and solvability of the
discrete problem are established. The method is suited for any Stokes inf-sup
stable finite element pair for velocity and pressure, while for vorticity any
generic discrete space (of arbitrary order) can be used. A priori and a
posteriori error estimates are derived using two specific families of discrete
subspaces. Finally, we provide a set of numerical tests illustrating the
behaviour of the scheme, verifying the theoretical convergence rates, and
showing the performance of the adaptive algorithm guided by residual a
posteriori error estimation.

###### keywords:

Oseen equations, velocity-vorticity-pressure formulation, mixed finite element
methods, variable viscosity, a priori and a posteriori error analysis,
adaptive mesh refinement.

65N30, 65N12, 76D07, 65N15.

##  1 Introduction

Using vorticity as additional field in the formulation of incompressible flow
equations can be advantageous in a number of applicative problems [47].
Starting from the seminal works [25, 26] that focused on Stokes equations and
where vorticity was sought in
$\mathrm{H}(\mathop{\mathbf{curl}}\nolimits,\Omega)$, several different
problems including Brinkman, Navier-Stokes, and related flow problems written
in terms of vorticity have been studied from the viewpoint of numerical
analysis of finite volume and mixed finite element methods exhibiting diverse
properties and specific features. Some of these contributions include [3, 2,
4, 7, 23, 12, 45, 46].

The starting point is the Oseen equations in the case of variable viscosity,
and written in terms of velocity $\boldsymbol{u}$ and pressure $p$, as follows
(see [35]):

(1.1a)
$\displaystyle\sigma\boldsymbol{u}-2\boldsymbol{\mathop{\mathrm{div}}\nolimits}(\nu\boldsymbol{\varepsilon}(\boldsymbol{u}))+(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}+\nabla
p$ $\displaystyle=\boldsymbol{f}$ $\displaystyle\mbox{ in }\Omega,$ (1.1b)
$\displaystyle\mathop{\mathrm{div}}\nolimits\boldsymbol{u}$ $\displaystyle=0$
$\displaystyle\mbox{ in }\Omega,$ (1.1c) $\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\mbox{ on }\Gamma,$ (1.1d)
$\displaystyle(p,1)_{0,\Omega}$ $\displaystyle=0,$

where $\sigma>0$ is inversely proportional to the time-step,
$\boldsymbol{f}\in\mathrm{L}^{2}(\Omega)^{d}$ is a force density,
$\boldsymbol{\beta}\in\mathrm{H}^{1}(\Omega)^{d}$ is the convecting velocity
field (not necessarily divergence-free), and $\nu\in W^{1,\infty}(\Omega)$ is
the kinematic viscosity of the fluid, satisfying

(1.2) $0<\nu_{0}\leq\nu\leq\nu_{1}.$

Such a set of equations will appear, for instance, in the linearisation of
non-Newtonian flow problems, as well as in applications where viscosity may
depend on temperature, concentration or volume fractions, or other fields
where the fluid flow patterns depend on marked spatial distributions of
viscosity [37, 41, 42, 44]. The specific literature related to the analysis of
numerical schemes for the Oseen equations in terms of vorticity includes the
non-conforming exponentially accurate least-squares spectral method proposed
in [40], least-squares methods proposed in [48] for Oseen and Navier-Stokes
equations with velocity boundary conditions, the family of vorticity-based
first-order Oseen-type systems studied in [21], the enhanced accuracy
formulation in terms of velocity-vorticity-helicity investigated in [11], and
the recent mixed (exactly divergence-free) and DG discretisations for Oseen’s
problem in velocity-vorticity-pressure form given in [5]. However, in most of
these references, the derivation of the variational formulations depends on
the viscosity being constant. This is attributed to the fact that the usual
vorticity-based weak formulation results from exploiting the following
identity

(1.3)
$\mathop{\mathbf{curl}}\nolimits(\mathop{\mathbf{curl}}\nolimits\boldsymbol{v})=-\Delta\boldsymbol{v}+\nabla(\mathop{\mathrm{div}}\nolimits\boldsymbol{v}),$

applied to the viscous term. However for a more general friction term of the
form
$-\boldsymbol{\mathop{\mathrm{div}}\nolimits}(\nu\boldsymbol{\varepsilon}(\boldsymbol{u}))$,
where $\boldsymbol{\varepsilon}(\boldsymbol{u})$ is the strain rate tensor,
the decomposition gives other additional terms that do not permit the direct
recasting of the coupled system as done in the cited references above.

Extensions to cover the case of variable viscosity do exist in the literature.
For instance, [27] addresses the well-posedness of the vorticity–velocity
formulation of the Stokes problem with varying density and viscosity, and the
equivalence of the vorticity–velocity and velocity–pressure formulations in
appropriate functional spaces is proved. More recently, in [6] we have taken a
different approach and employed an augmented vorticity-velocity-pressure
formulation for Brinkman equations with variable viscosity. Here we extend
that analysis to the generalised Oseen equations with variable viscosity, and
address in particular how to deal with the additional challenges posed by the
presence of the convective term that did not appear in the Brinkman momentum
equation.

We will employ the so-called augmented formulations (also known as Galerkin
least-squares methods), which can be regarded as a stabilisation technique
where some terms are added to the variational formulation. Augmented finite
elements have been considered in several works with applications in fluid
mechanics (see, e.g., [8, 9, 13, 20, 17, 18, 34, 43] and the references
therein). These methods enjoy appealing advantages as those described in
length in, e.g., [14, 16], and reformulations of the set of equations
following this approach are also of great importance in the design of block
preconditioners (see [10, 31] for an application in Oseen and Navier-Stokes
equations in primal form, [30] for stress-velocity-pressure formulations for
non-Newtonian flows, or [19, 29] for stress-displacement-pressure mixed
formulations for hyperelasticity). In the particular context of our mixed
formulation for Oseen equations, the augmentation assists us in deriving the
Babuška-Brezzi property of ellipticity on the kernel needed for the top-left
diagonal block.

The formulation that we employ is non-symmetric, and the augmentation terms
appear from least-squares contributions associated with the constitutive
equation and the incompressibility constraint. The mixed variational
formulation is shown to be well-posed under a condition on the viscosity
bounds (a generalisation of the usual condition needed in Oseen equations,
(cf. Theorem 2.5 and Remark 2.8). Then we establish the well-posedness of the
discrete problem for generic inf-sup stable finite elements (for velocity and
pressure) in combination with a generic space for vorticity approximation. We
obtain error estimates for two stable families of finite elements. We also
derive a reliable and efficient residual-based a posteriori error estimator
for the mixed problem, which can be fully computed locally. In summary, the
advantages of the proposed method are the possibility to obtain directly the
vorticity field with optimal accuracy and without the need of postprocessing;
moreover, different from many existing finite element methods with vorticity
field as unknown, the present contribution supports variable viscosity and no-
slip boundary condition in a natural way.

The contents of the paper have been structured as follows. Functional spaces
and recurrent notation is collected in the remainder of this section. Section
2 presents the governing equations in terms of velocity, vorticity and
pressure; we state an augmented formulation, and we perform the solvability
analysis invoking the Babuška–Brezzi theory. The finite element discretisation
is introduced in Section 3, where we also derive the stability analysis and
optimal error estimates for two families of stable elements. In Section 4, we
develop the a posteriori error analysis. Several numerical tests illustrating
the convergence of the proposed method under different scenarios are reported
in Section 5.

Preliminaries. Let $\Omega$ be a bounded domain of $\mathbb{R}^{d}$, $d=2,3$,
with Lipschitz boundary $\Gamma=\partial\Omega$. For any $s\geq 0$, the
notation $\left\|\cdot\right\|_{s,\Omega}$ stands for the norm of the
Hilbertian Sobolev spaces $\mathrm{H}^{s}(\Omega)$ or
$\mathrm{H}^{s}(\Omega)^{d}$, with the usual convention
$\mathrm{H}^{0}(\Omega):=\mathrm{L}^{2}(\Omega)$.

Moreover, $c$ and $C$, with or without subscripts, tildes, or hats, will
represent a generic constant independent of the mesh parameter $h$, assuming
different values in different occurrences. In addition, for any vector field
$\boldsymbol{v}=(v_{i})_{i=1}^{3}$ and any scalar field $q$ we recall the
notation:

$\mathop{\mathrm{div}}\nolimits\boldsymbol{v}=\sum_{i=1}^{3}\partial_{i}v_{i},\
\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}=\begin{pmatrix}\partial_{2}v_{3}-\partial_{3}v_{2}\\\
\partial_{3}v_{1}-\partial_{1}v_{3}\\\
\partial_{1}v_{2}-\partial_{2}v_{1}\end{pmatrix},\ \nabla
q=\begin{pmatrix}\partial_{1}q\\\ \partial_{2}q\\\
\partial_{3}q\end{pmatrix},$

whereas for dimension $d=2$, the curl of a vector $\boldsymbol{v}$ and a
scalar $q$ are scalar function $\partial_{1}v_{2}-\partial_{2}v_{1}$ and the
vector $\mathop{\mathbf{curl}}\nolimits q=(\partial_{2}q,\partial_{1}q)^{t}$,
respectively.

Recall that, according to [35, Theorem 2.11], for a generic domain
$\Omega\subseteq\mathbb{R}^{3}$, the relevant integration by parts formula
corresponds to

$\int_{\Omega}\mathop{\mathbf{curl}}\nolimits\boldsymbol{\omega}\cdot\boldsymbol{v}=\int_{\Omega}\boldsymbol{\omega}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}+\langle\boldsymbol{\omega}\times\boldsymbol{n},\boldsymbol{v}\rangle_{\Gamma},$

which in 2D reads as

(1.4)
$\int_{\Omega}\mathop{\mathbf{curl}}\nolimits\omega\cdot\boldsymbol{v}=\int_{\Omega}\omega\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}-\langle\boldsymbol{v}\cdot\boldsymbol{t},\omega\rangle_{\Gamma}.$

##  2 Vorticity-based formulation

With the aim of proposing a vorticity-based formulation for (1.1), we consider
the following identities

$\displaystyle-2\boldsymbol{\mathop{\mathrm{div}}\nolimits}(\nu\boldsymbol{\varepsilon}(\boldsymbol{u}))$
$\displaystyle=-2\nu\boldsymbol{\mathop{\mathrm{div}}\nolimits}(\boldsymbol{\varepsilon}(\boldsymbol{u}))-2\boldsymbol{\varepsilon}(\boldsymbol{u})\nabla\nu=-\nu\Delta\boldsymbol{u}-2\boldsymbol{\varepsilon}(\boldsymbol{u})\nabla\nu$
$\displaystyle=\nu\mathop{\mathbf{curl}}\nolimits(\mathop{\mathbf{curl}}\nolimits\boldsymbol{u})-\nu\nabla(\mathop{\mathrm{div}}\nolimits\boldsymbol{u})-2\boldsymbol{\varepsilon}(\boldsymbol{u})\nabla\nu.$

Therefore, problem (1.1) rewrites as

(2.5a)
$\displaystyle\sigma\boldsymbol{u}+\nu\mathop{\mathbf{curl}}\nolimits\boldsymbol{\omega}-2\boldsymbol{\varepsilon}(\boldsymbol{u})\nabla\nu+(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}+\nabla
p$ $\displaystyle=\boldsymbol{f}$ $\displaystyle\mbox{ in }\Omega,$ (2.5b)
$\displaystyle\boldsymbol{\omega}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\mbox{ in }\Omega,$ (2.5c)
$\displaystyle\mathop{\mathrm{div}}\nolimits\boldsymbol{u}$ $\displaystyle=0$
$\displaystyle\mbox{ in }\Omega,$ (2.5d) $\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\mbox{ on }\Gamma,$ (2.5e)
$\displaystyle(p,1)_{0,\Omega}$ $\displaystyle=0,$

where we have considered the definition of the vorticity and have applied the
incompressibility condition. The equations state, respectively, the momentum
conservation, the constitutive relation, the mass balance, the no-slip
boundary condition, and the pressure closure condition.

###  2.1 Variational formulation for the Oseen equations with non-constant
viscosity

In this section, we propose a mixed variational formulation of system
(2.5a)-(2.5e). First, we endow the space $\mathrm{H}_{0}^{1}(\Omega)^{d}$ with
the following norm:

$|\\!|\\!|\boldsymbol{v}|\\!|\\!|_{1,\Omega}^{2}:=\|\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}{\color[rgb]{1,0,0},}$

and note that for $\mathrm{H}_{0}^{1}(\Omega)^{d}$ the above norm is
equivalent to the usual norm. In particular, we have that there exists a
positive constant $C_{pf}$ such that:

$\|\boldsymbol{v}\|_{1,\Omega}^{2}\leq
C_{pf}(\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2})\qquad\forall\boldsymbol{v}\in\mathrm{H}_{0}^{1}(\Omega)^{d},$

where the above inequality is a consequence of the identity

(2.6)
$\|\nabla\boldsymbol{v}\|_{0,\Omega}^{2}=\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2},$

which follows from (1.3) and the Poincaré inequality. Moreover, in order to
establish a weak formulation for (2.5), we will use the following identity:

(2.7)
$\mathop{\mathbf{curl}}\nolimits(\phi\boldsymbol{v})=\nabla\phi\times\boldsymbol{v}+\phi\mathop{\mathbf{curl}}\nolimits\boldsymbol{v},$

valid for any vector field $\boldsymbol{v}$ and any scalar field $\phi$.

After testing each equation of (2.5a)-(2.5d) against adequate functions, using
(2.7), and imposing the boundary conditions, we end up with the following
system:

$\displaystyle\int_{\Omega}(\sigma\boldsymbol{u}+(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u})\cdot\boldsymbol{v}-2\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{u})\nabla\nu\cdot\boldsymbol{v}+\int_{\Omega}\nu\boldsymbol{\omega}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}+\int_{\Omega}\boldsymbol{\omega}\cdot(\nabla\nu\times\boldsymbol{v})-\int_{\Omega}p\mathop{\mathrm{div}}\nolimits\boldsymbol{v}$
$\displaystyle=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v},$
$\displaystyle\int_{\Omega}\nu\boldsymbol{\theta}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}-\int_{\Omega}\nu\boldsymbol{\omega}\cdot\boldsymbol{\theta}$
$\displaystyle=\,0,$
$\displaystyle-\int_{\Omega}q\mathop{\mathrm{div}}\nolimits\boldsymbol{u}$
$\displaystyle=\,0,$

for all
$(\boldsymbol{v},\boldsymbol{\theta},q)\in\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2}\times\mathrm{L}_{0}^{2}(\Omega)$,
where $\mathrm{L}_{0}^{2}(\Omega)\,:=\,\\{q\in\mathrm{L}^{2}(\Omega):\
(q,1)_{0,\Omega}\,=\,0\\}$.

Contrary to what is usually found in the the standard velocity-pressure mixed
formulation, the ellipticity on the kernel condition for the Babuška-Brezzi
theory is not straightforward in the above mixed formulation. Here is where
the augmentation contributes to simplify the analysis. We introduce the
following residual terms arising from equations (2.5b) and (2.5c):

$\kappa_{1}\int_{\Omega}(\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}-\boldsymbol{\omega})\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}=0,\qquad\kappa_{2}\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}=0\qquad\forall\boldsymbol{v}\in\mathrm{H}_{0}^{1}(\Omega)^{d},$

where $\kappa_{1}$ and $\kappa_{2}$ are positive parameters to be specified
later on.

In this way, we propose the following augmented variational formulation for
(2.5):

Find
$((\boldsymbol{u},\boldsymbol{\omega}),p)\in(\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2})\times\mathrm{L}_{0}^{2}(\Omega)$
such that

(2.8a) $\displaystyle
A((\boldsymbol{u},\boldsymbol{\omega}),(\boldsymbol{v},\boldsymbol{\theta}))+B((\boldsymbol{v},\boldsymbol{\theta}),p)=$
$\displaystyle\;F(\boldsymbol{v},\boldsymbol{\theta})$
$\displaystyle\qquad\forall(\boldsymbol{v},\boldsymbol{\theta})\in\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2},$
(2.8b) $\displaystyle B((\boldsymbol{u},\boldsymbol{\omega}),q)=$
$\displaystyle\;0$ $\displaystyle\qquad\forall
q\in\mathrm{L}_{0}^{2}(\Omega),$

where the bilinear forms and the linear functional are defined by

$\displaystyle
A((\boldsymbol{u},\boldsymbol{\omega}),(\boldsymbol{v},\boldsymbol{\theta}))$
$\displaystyle:=\int_{\Omega}(\sigma\boldsymbol{u}+(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u})\cdot\boldsymbol{v}+\int_{\Omega}\nu\boldsymbol{\omega}\cdot\boldsymbol{\theta}+\int_{\Omega}\nu\boldsymbol{\omega}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}-\int_{\Omega}\nu\boldsymbol{\theta}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}$
$\displaystyle\quad+\kappa_{1}\int_{\Omega}\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}+\kappa_{2}\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}-\kappa_{1}\int_{\Omega}\boldsymbol{\omega}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}$
(2.9a)
$\displaystyle\quad-2\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{u})\nabla\nu\cdot\boldsymbol{v}+\int_{\Omega}\boldsymbol{\omega}\cdot(\nabla\nu\times\boldsymbol{v}),$
(2.9b) $\displaystyle B((\boldsymbol{v},\boldsymbol{\theta}),q)$
$\displaystyle:=-\int_{\Omega}q\mathop{\mathrm{div}}\nolimits\boldsymbol{v},$
(2.9c) $\displaystyle F(\boldsymbol{v},\boldsymbol{\theta})$
$\displaystyle:=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v},$

for all
$(\boldsymbol{u},\boldsymbol{\omega}),(\boldsymbol{v},\boldsymbol{\theta})\in\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2}$,
and $q\in\mathrm{L}_{0}^{2}(\Omega)$.

As we will address in full detail in the next section, the augmented mixed
formulation will permit us to analyse the problem directly under the classical
Babuška-Brezzi theory [16].

###  2.2 Well-posedness analysis

In this section, we will address the well-posedness of the proposed weak
formulation (2.8).

In our analysis, we will need to invoke the following inequality, which is a
consequence of the Sobolev embedding
$\mathrm{H}^{1}(\Omega)\hookrightarrow\mathrm{L}^{4}(\Omega)$

(2.10)
$\left|\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}(\boldsymbol{u}\cdot\boldsymbol{v})\right|\leq\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega}|\\!|\\!|\boldsymbol{u}|\\!|\\!|_{1,\Omega}|\\!|\\!|\boldsymbol{v}|\\!|\\!|_{1,\Omega}.$

We will also make use of the following identity (cf. [35, Lemma 2.2])

(2.11)
$\int_{\Omega}[(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}]\cdot\boldsymbol{v}+\int_{\Omega}[(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{v}]\cdot\boldsymbol{u}=-\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}(\boldsymbol{u}\cdot\boldsymbol{v}).$

The continuity of the bilinear forms and the linear functional (cf.
(2.9)-(2.9c)), will be a consequence of the following lemma, whose proof
follows standard arguments in combination with (1.2).

######  Lemma 2.1.

The following estimates hold

$\displaystyle\left|\sigma\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}\right|\leq\sigma\|\boldsymbol{u}\|_{0,\Omega}\|\boldsymbol{v}\|_{0,\Omega},\qquad\left|\int_{\Omega}\nu\boldsymbol{\omega}\cdot\boldsymbol{\theta}\right|\leq\nu_{1}\|\boldsymbol{\omega}\|_{0,\Omega}\|\boldsymbol{\theta}\|_{0,\Omega},$
$\displaystyle\left|\int_{\Omega}[(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}]\cdot\boldsymbol{v}\right|\leq\widehat{C}|\\!|\\!|\boldsymbol{\beta}|\\!|\\!|_{1,\Omega}\|\nabla\boldsymbol{u}\|_{0,\Omega}|\\!|\\!|\boldsymbol{v}|\\!|\\!|_{1,\Omega},$
$\displaystyle\left|\int_{\Omega}\nu\boldsymbol{\theta}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\right|\leq\nu_{1}\|\boldsymbol{\theta}\|_{0,\Omega}|\\!|\\!|\boldsymbol{v}|\\!|\\!|_{1,\Omega},\qquad\left|\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{u})\nabla\nu\cdot\boldsymbol{v}\right|\leq\|\nabla\nu\|_{\infty,\Omega}\|\boldsymbol{\varepsilon}(\boldsymbol{u})\|_{0,\Omega}\|\boldsymbol{v}\|_{0,\Omega},$
$\displaystyle\left|\int_{\Omega}\boldsymbol{\theta}\cdot(\nabla\nu\times\boldsymbol{v})\right|\leq
2\|\nabla\nu\|_{\infty,\Omega}\|\boldsymbol{v}\|_{0,\Omega}\|\boldsymbol{\theta}\|_{0,\Omega},\qquad|F(\boldsymbol{v},\boldsymbol{\theta})|\leq\|\boldsymbol{f}\|_{0,\Omega}\|\boldsymbol{v}\|_{0,\Omega}.$

As a consequence of the above lemma, there exist constants
$C_{1},C_{2},C_{3}>0$ such that

$\displaystyle|A((\boldsymbol{u},\boldsymbol{\omega}),(\boldsymbol{v},\boldsymbol{\theta}))|\leq
C_{1}\|(\boldsymbol{u},\boldsymbol{\omega})\|\|(\boldsymbol{v},\boldsymbol{\theta})\|,\qquad|B((\boldsymbol{v},\boldsymbol{\theta}),q)|\leq
C_{2}\|(\boldsymbol{v},\boldsymbol{\theta})\|\|q\|_{0,\Omega},$
$\displaystyle|F(\boldsymbol{v},\boldsymbol{\theta})|\leq
C_{3}\|(\boldsymbol{v},\boldsymbol{\theta})\|,$

with the product space norm defined as

$\|(\boldsymbol{v},\boldsymbol{\theta})\|^{2}:=|\\!|\\!|\boldsymbol{v}|\\!|\\!|_{1,\Omega}^{2}+\|\boldsymbol{\theta}\|_{0,\Omega}^{2}.$

The following lemma states the ellipticity of the bilinear form
$A(\cdot,\cdot)$.

######  Lemma 2.2.

Assume that

(2.12) $\sigma>\frac{9\|\nabla\nu\|_{\infty,\Omega}^{2}}{\nu_{0}}\quad\mbox{
and
}\quad\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega}<\min\left\\{\sigma-\frac{9\|\nabla\nu\|_{\infty,\Omega}^{2}}{\nu_{0}},\dfrac{\nu_{0}}{12}\right\\}.$

Then, if we choose $\kappa_{1}=\frac{2}{3}\nu_{0}$ and
$\kappa_{2}>\dfrac{\nu_{0}}{3}$, there exists a constant $\alpha>0$ such that

$A((\boldsymbol{v},\boldsymbol{\theta}),(\boldsymbol{v},\boldsymbol{\theta}))\geq\alpha\|(\boldsymbol{v},\boldsymbol{\theta})\|^{2}\qquad\forall(\boldsymbol{v},\boldsymbol{\theta})\in\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2}.$

######  Proof 2.3.

Let
$(\boldsymbol{v},\boldsymbol{\theta})\in\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2}$.
As a consequence of Lemma 2.1, we have that

$\displaystyle\left|2\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{v})\nabla\nu\cdot\boldsymbol{v}\right|$
$\displaystyle\leq
2\|\nabla\nu\|_{\infty,\Omega}\left(\dfrac{\nu_{0}}{12\|\nabla\nu\|_{\infty,\Omega}}\|\nabla\boldsymbol{v}\|_{0,\Omega}^{2}+\dfrac{3\|\nabla\nu\|_{\infty,\Omega}}{\nu_{0}}\|\boldsymbol{v}\|_{0,\Omega}^{2}\right)$
(2.13)
$\displaystyle=\dfrac{\nu_{0}}{6}(\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2})+\dfrac{6\|\nabla\nu\|_{\infty,\Omega}^{2}}{\nu_{0}}\|\boldsymbol{v}\|_{0,\Omega}^{2},$

where we have used (2.6). Moreover, using that
$\|(\nabla\nu\times\boldsymbol{v})\|_{0,\Omega}\leq
2\|\nabla\nu\|_{\infty,\Omega}\|\boldsymbol{v}\|_{0,\Omega}$, we get

$\displaystyle\left|\int_{\Omega}\boldsymbol{\theta}\cdot(\nabla\nu\times\boldsymbol{v})\right|$
$\displaystyle\leq
2\|\nabla\nu\|_{\infty,\Omega}\left(\dfrac{\nu_{0}}{6\|\nabla\nu\|_{\infty,\Omega}}\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\dfrac{3\|\nabla\nu\|_{\infty,\Omega}}{2\nu_{0}}\|\boldsymbol{v}\|_{0,\Omega}^{2}\right)$
$\displaystyle=\dfrac{\nu_{0}}{3}\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\dfrac{3\|\nabla\nu\|_{\infty,\Omega}^{2}}{\nu_{0}}\|\boldsymbol{v}\|_{0,\Omega}^{2},$
$\displaystyle\left|\kappa_{1}\int_{\Omega}\boldsymbol{\theta}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\right|$
$\displaystyle\leq\kappa_{1}\left(\dfrac{\nu_{0}}{3\kappa_{1}}\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\dfrac{3\kappa_{1}}{4\nu_{0}}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}\right)$
(2.14)
$\displaystyle=\dfrac{\nu_{0}}{3}\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\dfrac{3\kappa_{1}^{2}}{4\nu_{0}}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}.$

Thus, using the Cauchy-Schwarz inequality, (2.3)-(2.14), (2.11) and (2.10), we
obtain

$\displaystyle
A((\boldsymbol{v},\boldsymbol{\theta}),(\boldsymbol{v},\boldsymbol{\theta}))\geq$
$\displaystyle\sigma\|\boldsymbol{v}\|_{0,\Omega}^{2}+\int_{\Omega}[(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{v}]\cdot\boldsymbol{v}+\int_{\Omega}\nu|\boldsymbol{\theta}|^{2}+\kappa_{1}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\kappa_{2}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}$
$\displaystyle-\kappa_{1}\int_{\Omega}\boldsymbol{\theta}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}-2\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{v})\nabla\nu\cdot\boldsymbol{v}+\int_{\Omega}\boldsymbol{\theta}\cdot(\nabla\nu\times\boldsymbol{v})$
$\displaystyle\geq$
$\displaystyle\sigma\|\boldsymbol{v}\|_{0,\Omega}^{2}-\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega}{|\\!|\\!|\boldsymbol{v}|\\!|\\!|_{1,\Omega}^{2}}+\nu_{0}\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\kappa_{1}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\kappa_{2}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}$
$\displaystyle-\dfrac{\nu_{0}}{3}\|\boldsymbol{\theta}\|_{0,\Omega}^{2}-\dfrac{3\kappa_{1}^{2}}{4\nu_{0}}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}-\dfrac{\nu_{0}}{6}(\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2})$
$\displaystyle-\dfrac{6\|\nabla\nu\|_{\infty,\Omega}^{2}}{\nu_{0}}\|\boldsymbol{v}\|_{0,\Omega}^{2}-\dfrac{\nu_{0}}{3}\|\boldsymbol{\theta}\|_{0,\Omega}^{2}-\dfrac{3\|\nabla\nu\|_{\infty,\Omega}^{2}}{\nu_{0}}\|\boldsymbol{v}\|_{0,\Omega}^{2}$
$\displaystyle=$
$\displaystyle\dfrac{\nu_{0}}{3}\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\left(\dfrac{\nu_{0}}{6}-\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega}\right)\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}$
$\displaystyle+\left(\kappa_{2}-\dfrac{\nu_{0}}{6}-\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega}\right)\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}$
$\displaystyle+\left(\sigma-\dfrac{9\|\nabla\nu\|_{\infty,\Omega}^{2}}{\nu_{0}}-\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega}\right)\|\boldsymbol{v}\|_{0,\Omega}^{2}.$

Now, using assumption (2.12), we have

$A((\boldsymbol{v},\boldsymbol{\theta}),(\boldsymbol{v},\boldsymbol{\theta}))\geq\alpha\|(\boldsymbol{v},\boldsymbol{\theta})\|^{2},$

where

$\alpha:=\min\left\\{{\dfrac{\nu_{0}}{3}},\dfrac{\nu_{0}}{6}-\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega},\kappa_{2}-\dfrac{\nu_{0}}{6}-\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega},\sigma-\dfrac{9\|\nabla\nu\|_{0,\Omega}^{2}}{\nu_{0}}-\widehat{C}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}\|_{0,\Omega}\right\\},$

which is clearly positive according to (2.12) and the assumptions on
$\kappa_{1}$ and $\kappa_{2}$.

Now we recall the following result related to the inf-sup condition: There
exists $C>0$, depending only on $\Omega$, such that (cf. [33])

$\sup_{0\neq\boldsymbol{v}\in\mathrm{H}_{0}^{1}(\Omega)^{d}}\frac{\left|\displaystyle\int_{\Omega}q\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\right|}{\|\boldsymbol{v}\|_{1,\Omega}}\geq
C\|q\|_{0,\Omega}\quad\forall q\in\mathrm{L}_{0}^{2}(\Omega).$

As a consequence, we immediately have the following lemma.

######  Lemma 2.4.

There exists $\gamma>0$, independent of $\nu$, such that

$\sup_{0\neq(\boldsymbol{v},\boldsymbol{\theta})\in\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2}}\frac{|B((\boldsymbol{v},\boldsymbol{\theta}),q)|}{\|(\boldsymbol{v},\boldsymbol{\theta})\|}\geq\gamma\|q\|_{0,\Omega}\quad\forall
q\in\mathrm{L}_{0}^{2}(\Omega).$

We state the well-posedness of problem (2.8) in the next theorem.

######  Theorem 2.5.

Assume that the hypotheses of Lemma 2.2 hold true. Then, there exists a unique
solution
$((\boldsymbol{u},\boldsymbol{\omega}),p)\in(\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2})\times\mathrm{L}_{0}^{2}(\Omega)$
to problem (2.8). Moreover, there exists $C>0$ such that

$\|(\boldsymbol{u},\boldsymbol{\omega})\|+\|p\|_{0,\Omega}\leq
C\|\boldsymbol{f}\|_{0,\Omega}.$

######  Proof 2.6.

The proof follows from Lemmas 2.2 and 2.4, and a direct consequence of the
Babuška-Brezzi Theorem ([16, Theorem II.1.1]).

######  Remark 2.7.

The unique solution of problem (2.8) also solves (2.5a)-(2.5e). The
equivalence follows essentially from applying integration by parts backwardly
in (2.8) and using suitable test functions. This is employed in Section 4 to
prove the efficiency of the a posteriori error estimator.

######  Remark 2.8.

If the convective velocity $\boldsymbol{\beta}\in\mathrm{H}^{1}(\Omega)^{d}$
is solenoidal (i.e., $\mathop{\mathrm{div}}\nolimits\boldsymbol{\beta}=0$ in
$\Omega$), then problem (2.8) is well-posed after choosing
$\kappa_{1}=\frac{2}{3}\nu_{0}$, $\kappa_{2}>\dfrac{\nu_{0}}{3}$, and assuming

(2.15) $\sigma\nu_{0}>{9\|\nabla\nu\|_{\infty,\Omega}^{2}}.$

##  3 Numerical discretisation

Let $\\{\mathcal{T}_{h}(\Omega)\\}_{h>0}$ be a shape-regular family of
partitions of the polygonal/polyhedral region $\bar{\Omega}$, by
triangles/tetrahedrons $T$ of diameter $h_{T}$, with the meshsize defined as
$h:=\max\\{h_{T}:\,T\in\mathcal{T}_{h}(\Omega)\\}$. In what follows, given an
integer $k\geq 0$ and a subset $S$ of $\mathbb{R}^{d}$, $\mathbb{P}_{k}(S)$
denotes the space of polynomial functions defined on $S$ and being of degree
$\leq$ $k$.

Now, we consider generic finite dimensional subspaces
$\mathbf{V}_{h}\subseteq\mathrm{H}_{0}^{1}(\Omega)^{d}$,
$\mathbf{W}_{h}\subseteq\mathrm{L}^{2}(\Omega)^{d(d-1)/2}$ and
$Q_{h}\subseteq\mathrm{L}_{0}^{2}(\Omega)$ such that the following discrete
inf-sup holds

(3.16)
$\sup_{0\neq(\boldsymbol{v}_{h},\boldsymbol{\theta}_{h})\in\mathbf{V}_{h}\times\mathbf{W}_{h}}\frac{|B((\boldsymbol{v}_{h},\boldsymbol{\theta}_{h}),q_{h})|}{\|(\boldsymbol{v}_{h},\boldsymbol{\theta}_{h})\|}\geq\gamma_{0}\|q_{h}\|_{0,\Omega}\quad\forall
q_{h}\in Q_{h},$

where $\gamma_{0}>0$ is independent of $h$.

In this way, the above inf-sup condition can be obtained if
($\mathbf{V}_{h},Q_{h}$) is an inf-sup stable pair for the classical Stokes
problem. Moreover, the discrete space
$\mathbf{W}_{h}\subseteq\mathrm{L}^{2}(\Omega)^{d(d-1)/2}$ for the vorticity
can be taken as continuous or discontinuous polynomial space. Here we will
consider both options.

Now, we are in a position to introduce the finite element scheme related to
problem (2.8): Find
$((\boldsymbol{u}_{h},\boldsymbol{\omega}_{h}),p_{h})\in(\mathbf{V}_{h}\times\mathbf{W}_{h})\times
Q_{h}$ such that

(3.17)
$\begin{split}A((\boldsymbol{u}_{h},\boldsymbol{\omega}_{h}),(\boldsymbol{v}_{h},\boldsymbol{\theta}_{h}))+B((\boldsymbol{v}_{h},\boldsymbol{\theta}_{h}),p_{h})=&\;F(\boldsymbol{v}_{h},\boldsymbol{\theta}_{h})\qquad\forall(\boldsymbol{v}_{h},\boldsymbol{\theta}_{h})\in\mathbf{V}_{h}\times\mathbf{W}_{h},\\\
B((\boldsymbol{u}_{h},\boldsymbol{\omega}_{h}),q_{h})=&\;0\qquad\qquad\quad\forall
q\in Q_{h}.\end{split}$

The next step is to establish the unique solvability and convergence of the
discrete problem (3.17).

######  Theorem 3.1.

Assume that the hypotheses of Lemma 2.2 hold true. Let
$\mathbf{V}_{h}\subseteq\mathrm{H}_{0}^{1}(\Omega)^{d}$,
$\mathbf{W}_{h}\subseteq\mathrm{L}^{2}(\Omega)^{d(d-1)/2}$ and
$Q_{h}\subseteq\mathrm{L}_{0}^{2}(\Omega)$ satisfy (3.16). Then, there exists
a unique
$((\boldsymbol{u}_{h},\boldsymbol{\omega}_{h}),p_{h})\in(\mathbf{V}_{h}\times\mathbf{W}_{h})\times
Q_{h}$ solution to (3.17). Moreover, there exist
$\hat{C}_{1},\,\hat{C}_{2}>0$, independent of $h$, such that

$\|(\boldsymbol{u}_{h},\boldsymbol{\omega}_{h})\|+\|p_{h}\|_{0,\Omega}\leq\hat{C}_{1}\|\boldsymbol{f}\|_{0,\Omega},$

and

(3.18)
$\begin{split}&\|(\boldsymbol{u},\boldsymbol{\omega})-(\boldsymbol{u}_{h},\boldsymbol{\omega}_{h})\|+\|p-p_{h}\|_{0,\Omega}\\\
&\qquad\qquad\qquad\leq\hat{C}_{2}\inf_{(\boldsymbol{v}_{h},\boldsymbol{\theta}_{h},q_{h})\in\mathbf{V}_{h}\times\mathbf{W}_{h}\times\mathrm{Q}_{h}}(|\\!|\\!|\boldsymbol{u}-\boldsymbol{v}_{h}|\\!|\\!|_{1,\Omega}+\|\boldsymbol{\omega}-\boldsymbol{\theta}_{h}\|_{0,\Omega}+\|p-q_{h}\|_{0,\Omega}),\end{split}$

where
$((\boldsymbol{u},\boldsymbol{\omega}),p)\in(\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2})\times\mathrm{L}_{0}^{2}(\Omega)$
is the unique solution of (2.8).

###  3.1 Discrete subspaces and error estimates

In this section, we will define explicit families of finite element subspaces
yielding the unique solvability of the discrete scheme (3.17). In addition, we
derive the corresponding rate of convergence for each family.

####  3.1.1 Taylor-Hood-$\mathbb{P}_{k}$

We start by introducing a family based on Taylor-Hood [36] finite elements for
velocity and pressure, and continuous or discontinuous piecewise polynomial
spaces for vorticity. More precisely, for any $k\geq 1$, we consider:

(3.19) $\begin{split}\mathbf{V}_{h}:&=\\{\boldsymbol{v}_{h}\in
C(\overline{\Omega})^{d}:\boldsymbol{v}_{h}|_{K}\in\mathbb{P}_{k+1}(K)^{d}\quad\forall
K\in\mathcal{T}_{h}\\}\cap\mathrm{H}_{0}^{1}(\Omega)^{d},\\\
Q_{h}:&=\\{q_{h}\in
C(\overline{\Omega}):q_{h}|_{K}\in\mathbb{P}_{k}(K)\quad\forall
K\in\mathcal{T}_{h}\\}\cap\mathrm{L}_{0}^{2}(\Omega),\\\
\mathbf{W}_{h}^{1}:&=\\{\boldsymbol{\theta}_{h}\in
C(\overline{\Omega})^{d(d-1)/2}:\boldsymbol{\theta}_{h}|_{K}\in\mathbb{P}_{k}(K)^{d(d-1)/2}\quad\forall
K\in\mathcal{T}_{h}\\},\\\
\mathbf{W}_{h}^{2}:&=\\{\boldsymbol{\theta}_{h}\in\mathrm{L}^{2}(\Omega)^{d(d-1)/2}:\boldsymbol{\theta}_{h}|_{K}\in\mathbb{P}_{k}(K)^{d(d-1)/2}\quad\forall
K\in\mathcal{T}_{h}\\}.\end{split}$

It is well known that $(\mathbf{V}_{h},Q_{h})$ satisfies the inf-sup condition
(3.16) [15]. In addition, we will consider continuous ($\mathbf{W}_{h}^{1}$)
and discontinuous ($\mathbf{W}_{h}^{2}$) polynomial approximations for
vorticity.

Now, we recall the approximation properties of the spaces specified in (3.19).
Assume that $\boldsymbol{u}\in\mathrm{H}^{1+s}(\Omega)^{d}$,
$p\in\mathrm{H}^{s}(\Omega)$ and
$\boldsymbol{\omega}\in\mathrm{H}^{s}(\Omega)^{d(d-1)/2}$, for some
$s\in(1/2,k+1]$. Then there exists $C>0$, independent of $h$, such that

(3.20a)
$\displaystyle\inf_{\boldsymbol{v}_{h}\in\mathbf{V}_{h}}|\\!|\\!|\boldsymbol{u}-\boldsymbol{v}_{h}|\\!|\\!|_{1,\Omega}$
$\displaystyle\leq Ch^{s}\|\boldsymbol{u}\|_{\mathrm{H}^{1+s}(\Omega)^{d}},$
(3.20b) $\displaystyle\inf_{q_{h}\in Q_{h}}\|p-q_{h}\|_{0,\Omega}$
$\displaystyle\leq Ch^{s}\|p\|_{\mathrm{H}^{s}(\Omega)},$ (3.20c)
$\displaystyle\inf_{\boldsymbol{\theta}_{h}\in\mathbf{W}_{h}^{1}}\|\boldsymbol{\omega}-\boldsymbol{\theta}_{h}\|_{0,\Omega}$
$\displaystyle\leq
Ch^{s}\|\boldsymbol{\omega}\|_{\mathrm{H}^{s}(\Omega)^{d(d-1)/2}},$ (3.20d)
$\displaystyle\inf_{\boldsymbol{\theta}_{h}\in\mathbf{W}_{h}^{2}}\|\boldsymbol{\omega}-\boldsymbol{\theta}_{h}\|_{0,\Omega}$
$\displaystyle\leq
Ch^{s}\|\boldsymbol{\omega}\|_{\mathrm{H}^{s}(\Omega)^{d(d-1)/2}}.$

The following theorem provides the rate of convergence of the augmented mixed
scheme (3.17).

######  Theorem 3.2.

Let $k\geq 1$ be an integer, and let $\mathbf{V}_{h},Q_{h}$ and $W^{i}_{h}$,
$i=1,2$ be specified by (3.19). Let
$(\boldsymbol{u},\boldsymbol{\omega},p)\in\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2}\times\mathrm{L}_{0}^{2}(\Omega)$
and
$(\boldsymbol{u}_{h},\boldsymbol{\omega}_{h},p_{h})\in\mathbf{V}_{h}\times\mathbf{W}_{h}^{i}\times
Q_{h}$ be the unique solutions to the continuous and discrete problems (2.8)
and (3.17), respectively. Assume that
$\boldsymbol{u}\in\mathrm{H}^{1+s}(\Omega)^{d}$,
$\boldsymbol{\omega}\in\mathrm{H}^{s}(\Omega)^{d(d-1)/2}$ and
$p\in\mathrm{H}^{s}(\Omega)$, for some $s\in(1/2,k+1]$. Then, there exists
$\hat{C}>0$, independent of $h$, such that

$\|(\boldsymbol{u},\boldsymbol{\omega})-(\boldsymbol{u}_{h},\boldsymbol{\omega}_{h})\|+\|p-p_{h}\|_{0,\Omega}\leq\hat{C}h^{s}(\|\boldsymbol{u}\|_{\mathrm{H}^{1+s}(\Omega)^{d}}+\|\boldsymbol{\omega}\|_{\mathrm{H}^{s}(\Omega)^{d(d-1)/2}}+\|p\|_{\mathrm{H}^{s}(\Omega)}).$

######  Proof 3.3.

The proof follows from (3.18) and the approximation properties
(3.20a)-(3.20d).

####  3.1.2 MINI-element-$\mathbb{P}_{k}$

The second finite element family uses the so-called MINI-element for velocity
and pressure, and continuous or discontinuous piecewise polynomials for
vorticity. Let us introduce the following spaces (see [16, Sections 8.6 and
8.7], for further details):

$\displaystyle\mathbf{U}_{h}:$ $\displaystyle=\\{\boldsymbol{v}_{h}\in
C(\overline{\Omega})^{d}:\boldsymbol{v}_{h}|_{K}\in\mathbb{P}_{k}(K)^{d}\quad\forall
K\in\mathcal{T}_{h}\\},$ $\displaystyle\mathbb{B}(b_{K}\nabla H_{h}):$
$\displaystyle=\\{\boldsymbol{v}_{hb}\in\mathrm{H}^{1}(\Omega)^{d}:\boldsymbol{v}_{hb}|_{K}=b_{K}\nabla(q_{h})|_{K}\,\text{for
some}\,q_{h}\in H_{h}\\},$

where $b_{K}$ is the standard (cubic or quartic) bubble function
$\lambda_{1}\cdots\lambda_{d+1}\in\mathbb{P}_{d+1}(K)$, and let us define the
following finite element subspaces:

(3.21) $\begin{split}Q_{h}:&=\\{q_{h}\in
C(\overline{\Omega}):q_{h}|_{K}\in\mathbb{P}_{k}(K)\quad\forall
K\in\mathcal{T}_{h}\\}\cap\mathrm{L}_{0}^{2}(\Omega),\\\
\mathbf{V}_{h}:&=\mathbf{U}_{h}\oplus\mathbb{B}(b_{K}\nabla
Q_{h})\cap\mathrm{H}_{0}^{1}(\Omega)^{d},\\\
\mathbf{W}_{h}^{1}:&=\\{\boldsymbol{\theta}_{h}\in
C(\overline{\Omega})^{d(d-1)/2}:\boldsymbol{\theta}_{h}|_{K}\in\mathbb{P}_{k}(K)^{d(d-1)/2}\quad\forall
K\in\mathcal{T}_{h}\\},\\\
\mathbf{W}_{h}^{2}:&=\\{\boldsymbol{\theta}_{h}\in\mathrm{L}^{2}(\Omega)^{d(d-1)/2}:\boldsymbol{\theta}_{h}|_{K}\in\mathbb{P}_{k}(K)^{d(d-1)/2}\quad\forall
K\in\mathcal{T}_{h}\\}.\end{split}$

The rate of convergence of our augmented mixed finite element scheme
considering the above discrete spaces (3.21) is as follows.

######  Theorem 3.4.

Let $k\geq 1$ be an integer, and let $\mathbf{V}_{h},Q_{h}$ and $W^{i}_{h}$,
$i=1,2$ be given by (3.21). Let
$(\boldsymbol{u},\boldsymbol{\omega},p)\in\mathrm{H}_{0}^{1}(\Omega)^{d}\times\mathrm{L}^{2}(\Omega)^{d(d-1)/2}\times\mathrm{L}_{0}^{2}(\Omega)$
and
$(\boldsymbol{u}_{h},\boldsymbol{\omega}_{h},p_{h})\in\mathbf{V}_{h}\times\mathbf{W}_{h}^{i}\times
Q_{h}$ be the unique solutions to the continuous and discrete problems (2.8)
and (3.17), respectively. Assume that
$\boldsymbol{u}\in\mathrm{H}^{1+s}(\Omega)^{d}$,
$\boldsymbol{\omega}\in\mathrm{H}^{s}(\Omega)^{d(d-1)/2}$ and
$p\in\mathrm{H}^{s}(\Omega)$, for some $s\in(1/2,k]$. Then, there exists
$\hat{C}>0$, independent of $h$, such that

$\|(\boldsymbol{u},\boldsymbol{\omega})-(\boldsymbol{u}_{h},\boldsymbol{\omega}_{h})\|+\|p-p_{h}\|_{0,\Omega}\leq\hat{C}h^{s}(\|\boldsymbol{u}\|_{\mathrm{H}^{1+s}(\Omega)^{d}}+\|\boldsymbol{\omega}\|_{\mathrm{H}^{s}(\Omega)^{d(d-1)/2}}+\|p\|_{\mathrm{H}^{s}(\Omega)}).$

##  4 A posteriori error estimator

In this section, we propose a residual-based a posteriori error estimator and
prove its reliability and efficiency. The analysis restricts to the two-
dimensional case and using continuous finite element approximations for
vorticity. Nevertheless, the extension to 3D and to discontinuous vorticity
follows straightforwardly.

For each $T\in{\mathcal{T}}_{h}$ we let ${\mathcal{E}}(T)$ be the set of edges
of $T$, and we denote by ${\mathcal{E}}_{h}$ the set of all edges in
${\mathcal{T}}_{h}$, that is

${\mathcal{E}}_{h}={\mathcal{E}}_{h}(\Omega)\cup{\mathcal{E}}_{h}(\Gamma),$

where ${\mathcal{E}}_{h}(\Omega):=\\{e\in{\mathcal{E}}_{h}:e\subset\Omega\\}$,
and ${\mathcal{E}}_{h}(\Gamma):=\\{e\in{\mathcal{E}}_{h}:e\subset\Gamma\\}$.
In what follows, $h_{e}$ stands for the diameter of a given edge
$e\in{\mathcal{E}}_{h}$, $\boldsymbol{t}_{e}=(-n_{2},n_{1})$, where
$\boldsymbol{n}_{e}=(n_{1},n_{2})$ is a fix unit normal vector of $e$. Now,
let $q\in\mathrm{L}^{2}(\Omega)$ such that $q|_{T}\in C(T)$ for each
$T\in{\mathcal{T}}_{h}$, then given $e\in{\mathcal{E}}_{h}(\Omega)$, we denote
by $[q]$ the jump of $q$ across $e$, that is
$[q]:=(q|_{T^{\prime}})|_{e}-(q|_{T^{\prime\prime}})|_{e}$, where $T^{\prime}$
and $T^{\prime\prime}$ are the triangles of ${\mathcal{T}}_{h}$ sharing the
edge $e$. Moreover, let $\boldsymbol{v}\in\mathrm{L}^{2}(\Omega)^{2}$ such
that $\boldsymbol{v}|_{T}\in C(T)^{2}$ for each $T\in{\mathcal{T}}_{h}$. Then,
given $e\in{\mathcal{E}}_{h}(\Omega)$, we denote by
$[\boldsymbol{v}\cdot\boldsymbol{t}]$ the tangential jump of $\boldsymbol{v}$
across $e$, that is,
$[\boldsymbol{v}\cdot\boldsymbol{t}]:=\left((\boldsymbol{v}|_{T^{\prime}})|_{e}-(\boldsymbol{v}|_{T^{\prime\prime}})|_{e}\right)\cdot\boldsymbol{t}_{e}$,
where $T^{\prime}$ and $T^{\prime\prime}$ are the triangles of
${\mathcal{T}}_{h}$ sharing the edge $e$.

Next, let $k\geq 1$ be an integer, and let $\mathbf{V}_{h},Q_{h}$ and
$\mathbf{W}_{h}^{1}$ be given as in (3.19) or (3.21). Let
$(\boldsymbol{u},\omega,p)\in\mathrm{H}_{0}^{1}(\Omega)^{2}\times\mathrm{L}^{2}(\Omega)\times\mathrm{L}_{0}^{2}(\Omega)$
and
$(\boldsymbol{u}_{h},\omega_{h},p_{h})\in\mathbf{V}_{h}\times\mathbf{W}_{h}^{1}\times
Q_{h}$ be the unique solutions to the continuous and discrete problems (2.8)
and (3.17), respectively. We introduce for each $T\in{\mathcal{T}}_{h}$ the
local a posteriori error indicator and its global counterpart as

$\displaystyle\Theta_{T}^{2}:=$ $\displaystyle
h_{T}^{2}\|\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h}\|_{0,T}^{2}$ (4.22)
$\displaystyle+\|\omega_{h}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h}\|_{0,T}^{2}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}\|_{0,T}^{2},\qquad\Theta^{2}:=\sum_{T\in{\mathcal{T}}_{h}}\Theta_{T}^{2}.$

Let us now establish reliability and efficiency of (4.22).

###  4.1 Reliability

We begin by recalling that the continuous dependence result given in Theorem
2.5 is equivalent to the global inf–sup condition for the continuous
formulation (2.8). Then, applying this estimate to the error
$(\boldsymbol{u}-\boldsymbol{u}_{h},\omega-\omega_{h},p-p_{h})$, we obtain

(4.23)
$\|(\boldsymbol{u},\omega)-(\boldsymbol{u}_{h},\omega_{h})\|+\|p-p_{h}\|_{0,\Omega}\leq
C_{glob}\sup_{(\boldsymbol{v}_{h},\theta_{h},q_{h})\in}\frac{\mathcal{R}(\boldsymbol{v},\theta,q)}{\|(\boldsymbol{v},\theta,q)\|},$

where the residual functional $\mathcal{R}$ is defined by

(4.24)
$\mathcal{R}(\boldsymbol{v},\theta,q)=A((\boldsymbol{u}-\boldsymbol{u}_{h},\omega-\omega_{h}),(\boldsymbol{v},\theta))+B((\boldsymbol{v},\theta),{p-p_{h}})+B((\boldsymbol{u}-\boldsymbol{u}_{h},\omega-\omega_{h}),q),$

for all
$(\boldsymbol{v},\theta,q)\in\mathrm{H}_{0}^{1}(\Omega)^{2}\times\mathrm{L}^{2}(\Omega)\times\mathrm{L}_{0}^{2}(\Omega)$.

Some technical results are provided beforehand. Let us first recall the
Clément-type interpolation operator
$\mathcal{I}_{h}:\mathrm{H}^{1}_{0}(\Omega)\to Y_{h}$, where
$Y_{h}:=\\{v_{h}\in
C(\overline{\Omega})\cap\mathrm{H}^{1}_{0}(\Omega):v_{h}\Big{|}_{T}\in\mathbb{P}_{1}(T),\forall
T\in\mathcal{T}_{h}\\}.$ This operator satisfies the following local
approximation properties (cf. [22]).

######  Lemma 4.1.

There exist positive constants $C_{1}$ and $C_{2}$ such that for all
$v\in\mathrm{H}^{1}_{0}(\Omega)$ there hold

(4.25a) $\displaystyle\|v-\mathcal{I}_{h}v\|_{0,T}$ $\displaystyle\leq
C_{1}h_{T}|v|_{1,w_{T}}\quad\forall T\in\mathcal{T}_{h},$ (4.25b)
$\displaystyle\|v-\mathcal{I}_{h}v\|_{0,e}$ $\displaystyle\leq
C_{2}h^{1/2}_{e}|v|_{1,w_{e}}\quad\forall e\in{\mathcal{E}}_{h}(\Omega),$

where $w_{T}:=\bigcup\\{T^{\prime}\in\mathcal{T}_{h}:T^{\prime}\cap
T\neq\emptyset\\}$ and
$w_{e}:=\bigcup\\{T^{\prime}\in\mathcal{T}_{h}:T^{\prime}\cap
e\neq\emptyset\\}$.

The main result of this section is stated as follows.

######  Theorem 4.2.

There exists a positive constant $C_{\mathrm{rel}}$, independent of $h$, such
that

(4.26)
$\|(\boldsymbol{u},\omega)-(\boldsymbol{u}_{h},\omega_{h})\|+\|p-p_{h}\|_{0,\Omega}\leq{C}_{\mathrm{rel}}\;\Theta.$

######  Proof 4.3.

From (4.24) and the continuous problem (2.8), we have that,

$\displaystyle\mathcal{R}(\boldsymbol{v},\theta,q)$
$\displaystyle=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}-\Big{(}A((\boldsymbol{u}_{h},\omega_{h}),(\boldsymbol{v},\theta))+B((\boldsymbol{v},\theta),p_{h})+B((\boldsymbol{u}_{h},\omega_{h}),q)\Big{)}$
$\displaystyle=\int_{\Omega}\big{(}\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu\big{)}\cdot\boldsymbol{v}-\int_{\Omega}\nu(\omega_{h}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h})\theta$
$\displaystyle\quad-\kappa_{1}\int_{\Omega}(\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h}-\omega_{h})\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}-\kappa_{2}\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}$
$\displaystyle\quad-\left(\int_{\Omega}\nu\omega_{h}\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}+\int_{\Omega}\omega_{h}(\nabla\nu\times\boldsymbol{v})\right)+\int_{\Omega}p_{h}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}+\int_{\Omega}q\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}.$

Using the identity
$\mathop{\mathbf{curl}}\nolimits(\nu\boldsymbol{v})=\nabla\nu\times\boldsymbol{v}+\nu\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}$
and integration by parts on the above residual (cf. (1.4)), we obtain

$\displaystyle\mathcal{R}(\boldsymbol{v},\theta,q)$
$\displaystyle=\int_{\Omega}(\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu)\cdot\boldsymbol{v}-\int_{\Omega}\nu(\omega_{h}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h})\theta$
$\displaystyle-\kappa_{1}\int_{\Omega}(\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h}-\omega_{h})\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}-\kappa_{2}\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}+\int_{\Omega}q\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}$
$\displaystyle-\sum_{T\in\mathcal{T}_{h}}\left(\int_{T}\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}\cdot\boldsymbol{v}-\langle\boldsymbol{v}\cdot\boldsymbol{t},\nu\omega_{h}\rangle_{\partial
T}-\int_{T}\nabla
p_{h}\cdot\boldsymbol{v}+\langle\boldsymbol{v}\cdot\boldsymbol{n},p_{h}\rangle_{\partial
T}\right)$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}}\int_{T}(\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h})\cdot\boldsymbol{v}$
$\displaystyle-\int_{\Omega}\nu(\omega_{h}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h})\theta-\kappa_{1}\int_{\Omega}(\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h}-\omega_{h})\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}-\kappa_{2}\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}+\int_{\Omega}q\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h},$

where we have used the fact that $\omega_{h}$ and $p_{h}$ are piecewise
continuous functions. Hence, since from (4.24) we have
$\mathcal{R}(\boldsymbol{v}_{h},\theta_{h},q_{h})=0$, we obtain

$\displaystyle\mathcal{R}(\boldsymbol{v},\theta,q)$
$\displaystyle=\mathcal{R}(\boldsymbol{v}-\boldsymbol{v}_{h},\theta-\theta_{h},q-q_{h})$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}}\int_{T}\big{(}\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h}\big{)}\cdot(\boldsymbol{v}-\boldsymbol{v}_{h})$
$\displaystyle\quad-\int_{\Omega}\nu(\omega_{h}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h})(\theta-\theta_{h})-\kappa_{1}\int_{\Omega}(\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h}-\omega_{h})\mathop{\mathbf{curl}}\nolimits(\boldsymbol{v}-\boldsymbol{v}_{h})$
$\displaystyle\quad-\kappa_{2}\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}\mathop{\mathrm{div}}\nolimits(\boldsymbol{v}-\boldsymbol{v}_{h})+\int_{\Omega}(q-q_{h})\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}.$

Thus, it suffices to take
$\boldsymbol{v}_{h}:=\mathcal{I}_{h}(\boldsymbol{v})$ (cf. Lemma 4.1), and
$\theta_{h}:=\Pi(\theta)$ and $q_{h}:=\Pi(q)$ with $\Pi$ being the
$L^{2}$-projection onto piecewise constants. And then, using the Cauchy-
Schwarz inequality, triangle inequality, properties for $\mathcal{I}_{h}$
given by Lemma 4.1 and [28, Lemma 1.127], and approximation properties for
$\Pi$, we obtain

$\displaystyle\mathcal{R}(\boldsymbol{v},\theta,q)$ $\displaystyle\leq
C_{1}\sum_{T\in\mathcal{T}_{h}}h_{T}\|\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h}\|_{0,T}|\boldsymbol{v}|_{1,w_{T}}$
$\displaystyle\quad+\sum_{T\in\mathcal{T}_{h}}(\nu_{1}+\kappa_{1})\|\omega_{h}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h}\|_{0,T}(C_{3}\|\theta\|_{0,T}+|\boldsymbol{v}-\boldsymbol{v}_{h}|_{1,T})$
$\displaystyle\quad+\sum_{T\in\mathcal{T}_{h}}(\kappa_{2}+1)\|\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}\|_{0,T}(|\boldsymbol{v}-\boldsymbol{v}_{h}|_{1,T}+C_{4}\|q\|_{0,T})$
$\displaystyle\leq\widehat{C}_{1}\left(\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\|f-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h}\|_{0,T}^{2}\right)^{1/2}\\!\\!\|\boldsymbol{v}\|_{1,\Omega}$
$\displaystyle\quad+\widehat{C}_{2}\left(\sum_{T\in\mathcal{T}_{h}}\|\omega_{h}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h}\|_{0,T}^{2}\right)^{1/2}(\|\theta\|_{0,\Omega}+{\|\boldsymbol{v}\|_{1,\Omega}})$
$\displaystyle\quad+\widehat{C}_{3}\left(\sum_{T\in\mathcal{T}_{h}}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}\|_{0,T}^{2}\right)^{1/2}(\|\boldsymbol{v}\|_{1,\Omega}+\|q\|_{0,\Omega}).$

And the proof of (4.26) follows from (4.23) and the above estimate.

###  4.2 Efficiency

This subsection deals with the efficiency of the a posteriori error estimator.
For simplicity, we will assume that the given convective velocity
$\boldsymbol{\beta}$ and the viscosity $\nu$ are polynomial functions both of
degree $s$. The general case can be proved by repeating the same arguments and
requiring an additional regularity for the data.

A major role in the proof of efficiency is played by element and edge bubbles
(locally supported non-negative functions), whose definition we recall in what
follows. For $T\in\mathcal{T}_{h}(\Omega)$ and $e\in{\mathcal{E}}(T)$, let
$\psi_{T}$ and $\psi_{e}$, respectively, be the interior and edge bubble
functions defined as in, e.g., [1]. Let $\psi_{T}\in\mathbb{P}_{3}(T)$ with
$supp(\psi_{T})\subset T,$ $\psi_{T}=0$ on $\partial T$ and $0\leq\psi_{T}\leq
1$ in $T.$ Moreover, let $\psi_{e}|_{T}\in\mathbb{P}_{2}(T)$ with
$supp(\psi_{e})\subset\Omega_{e}:=\\{T^{\prime}\in\mathcal{T}_{h}(\Omega):e\in{\mathcal{E}}(T^{\prime})\\},$
$\psi_{e}=0$ on $\partial T\setminus e,$ and $0\leq\psi_{e}\leq 1$ in
$\Omega_{e}.$ Again, let us recall an extension operator $E:C^{0}(e)\mapsto
C^{0}(T)$ that satisfies $E(q)\in\mathbb{P}_{k}(T)$ and $E(q)|_{e}=q$ for all
$q\in\mathbb{P}_{k}(e)$ and for all $k\in\mathbb{N}\cup\\{0\\}.$

We now summarise the properties of $\psi_{T},\psi_{e}$ and $E$ in the
following lemma (see [1, 49]).

######  Lemma 4.4.

The following properties hold:

  * (i)

For $T\in\mathcal{T}_{h}$ and for $v\in\mathbb{P}_{k}(T)$, there is a positive
constant $C_{1}$ such that

$\displaystyle
C_{1}^{-1}\;\|v\|^{2}_{0,T}\leq\int_{T}\psi_{T}v^{2}\,\mbox{d}x\leq
C_{1}\|v\|^{2}_{0,T},\qquad C_{1}^{-1}\;\|v\|^{2}_{0,T}\leq\|\psi
v\|^{2}_{0,T}+h_{T}^{2}|\psi v|^{2}_{1,T}\leq C_{1}\|v\|^{2}_{0,T}.$

  * (ii)

For $e\in{\mathcal{E}}_{h}$ and $v\in\mathbb{P}_{k}(e)$, there exists a
positive constant, say $C_{1}$, such that

$C_{1}^{-1}\;\|v\|^{2}_{0,e}\leq\int_{e}\psi_{e}v^{2}ds\leq
C_{1}\|v\|^{2}_{0,e}.$

  * (iii)

For $T\in\mathcal{T}_{h}$, $e\in{\mathcal{E}}(T)$ and for
$v\in\mathbb{P}_{k}(e)$, there is a positive constant, again say $C_{1}$, such
that

$\|\psi_{e}^{1/2}\;E(v)\|^{2}_{0,T}\leq C_{1}h_{e}\;\|v\|^{2}_{0,e}.$

The following classical result which states an inverse estimate will also be
used.

######  Lemma 4.5.

Let $k,l,m\in\mathbb{N}\cup\\{0\\}$ such that $l\leq m$. Then, there exists
$\tilde{C}>0$, depending only on $k,l,m$ and the shape regularity of the
triangulations, such that for each triangle $T$ there holds

$|q|_{m,T}\leq\tilde{C}h_{T}^{l-m}|q|_{l,T}\quad\forall
q\in\mathbb{P}_{k}(T).$

In order to prove the efficiency of the a posteriori error estimator, we will
bound each term defining $\Theta_{T}$ in terms of local errors.

######  Theorem 4.6.

There is a positive constant $C_{\mathrm{eff}}$, independent of $h$, such that

$C_{\mathrm{eff}}\;\Theta\leq\;\|(\boldsymbol{u},\omega)-(\boldsymbol{u}_{h},\omega_{h})\|+\|p-p_{h}\|_{0,\Omega}+\mathrm{h.o.t.},$

where $\mathrm{h.o.t.}$ denotes higher-order terms.

######  Proof 4.7.

Using that $\omega-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}=0$ and
$\mathop{\mathrm{div}}\nolimits\boldsymbol{u}=0$ in $\Omega$ (see (2.5b) and
(2.5c), respectively), we immediately have that

$\|\omega_{h}-\mathop{\mathbf{curl}}\nolimits\boldsymbol{u}_{h}\|_{0,T}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{u}_{h}\|_{0,T}\leq\|\mathop{\mathbf{curl}}\nolimits(\boldsymbol{u}-\boldsymbol{u}_{h})\|_{0,T}+\|\mathop{\mathrm{div}}\nolimits(\boldsymbol{u}-\boldsymbol{u}_{h})\|_{0,T}+\|\omega-\omega_{h}\|_{0,T}.$

On the other hand, with the help of the $\mathrm{L}^{2}(T)^{2}$-orthogonal
projection $\mathcal{P}_{T}^{\ell}$ onto $\mathbb{P}_{\ell}(T)^{2}$, for
$\ell\geq(s+k+1),$ with respect to the weighted $\mathrm{L}^{2}$-inner product
$(\psi_{T}\boldsymbol{f},\boldsymbol{g})_{0,T}$, for
$\boldsymbol{f},\boldsymbol{g}\in\mathrm{L}^{2}(T)^{2},$ it now follows that

$\displaystyle\|\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-$
$\displaystyle\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h}\|_{0,T}^{2}$
$\displaystyle=\|\boldsymbol{f}-\mathcal{P}_{T}^{\ell}(\boldsymbol{f})+\mathcal{P}_{T}^{\ell}(\boldsymbol{f})-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h}\|_{0,T}^{2}$
$\displaystyle\leq\|\boldsymbol{f}-\mathcal{P}_{T}^{\ell}(\boldsymbol{f})\|_{0,T}^{2}+\|\mathcal{P}_{T}^{\ell}(\boldsymbol{f})-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h}\|_{0,T}^{2}$
$\displaystyle=\|\boldsymbol{f}-\mathcal{P}_{T}^{\ell}(\boldsymbol{f})\|_{0,T}^{2}+\|\mathcal{P}_{T}^{\ell}(\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h})\|_{0,T}^{2}.$

For the second term on the right-hand side, an application of Lemma 4.4 shows
that

$\displaystyle\|\mathcal{P}_{T}^{\ell}(\boldsymbol{f}-\sigma\boldsymbol{u}_{h}$
$\displaystyle-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h})\|_{0,T}^{2}$
$\displaystyle\leq\|\psi_{T}^{1/2}\mathcal{P}_{T}^{\ell}(\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h})\|_{0,T}^{2}$
$\displaystyle=\int_{T}\psi_{T}\mathcal{P}_{T}^{\ell}(\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h})$
$\displaystyle\qquad\times(\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h}),$

where we have used the fact that $\mathcal{P}_{T}^{\ell}$ is the
$L^{2}(T)^{2}$-orthogonal projection. Thus, from the above inequality, and
(2.5a) (cf. Remark 2.7), we can deduce that

$\displaystyle\|\mathcal{P}_{T}^{\ell}(\boldsymbol{f}$
$\displaystyle-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h})\|_{0,T}^{2}$
$\displaystyle{\leq}\int_{T}\psi_{T}\mathcal{P}_{T}^{\ell}(\boldsymbol{f}-\sigma\boldsymbol{u}_{h}-\nu\mathop{\mathbf{curl}}\nolimits\omega_{h}-(\boldsymbol{\beta}\cdot\nabla)\boldsymbol{u}_{h}+2\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\nabla\nu-\nabla
p_{h})$
$\displaystyle\quad\times(\sigma(\boldsymbol{u}-\boldsymbol{u}_{h})+\nu\mathop{\mathbf{curl}}\nolimits(\omega-\omega_{h})+(\boldsymbol{\beta}\cdot\nabla)(\boldsymbol{u}-\boldsymbol{u}_{h})-2\boldsymbol{\varepsilon}(\boldsymbol{u}-\boldsymbol{u}_{h})\nabla\nu+\nabla(p-p_{h})).$

Next, using that the viscosity is a polynomial function, the bound follows by
integration by parts on the terms
$\mathop{\mathbf{curl}}\nolimits(\omega-\omega_{h})$ and $\nabla(p-p_{h})$,
Cauchy-Schwarz inequality and an inverse inequality (cf. Lemma 4.5). We end
the proof by observing that the required efficiency bound follows
straightforwardly from the estimates above, and after assuming additional
regularity for $\boldsymbol{f}$.

##  5 Numerical results

In this section, we present some numerical experiments carried out with the
schemes proposed and analysed in Section 3. We also present two numerical
examples in $\mathbb{R}^{2}$, confirming the reliability and efficiency of the
a posteriori error estimator $\Theta$ derived in Section 4, and showing the
behaviour of the associated adaptive algorithm. The solution of all linear
systems is carried out with the multifrontal massively parallel sparse direct
solver MUMPS.

We construct a series of uniformly successively refined triangular meshes for
$\Omega$ and compute individual errors

$e(\boldsymbol{u})=|\\!|\\!|\boldsymbol{u}-\boldsymbol{u}_{h}|\\!|\\!|_{1,\Omega},\quad
e(\boldsymbol{\omega})=\|\boldsymbol{\omega}-\boldsymbol{\omega}_{h}\|_{0,\Omega},\quad
e(p)=\|p-p_{h}\|_{0,\Omega},$

and convergence rates

(5.27)
$r(\boldsymbol{u})=\dfrac{\log(e(\boldsymbol{u})/\widehat{e}(\boldsymbol{u}))}{\log(h/\hat{h})},\qquad
r(\boldsymbol{\omega})=\dfrac{\log(e(\boldsymbol{\omega})/\widehat{e}(\boldsymbol{\omega}))}{\log(h/\hat{h})},\qquad
r(p)=\dfrac{\log(e(p)/\widehat{e}(p))}{\log(h/\hat{h})},$

where $e,\widehat{e}$ denote errors generated on two consecutive meshes of
sizes $h,\hat{h}$, respectively.

###  5.1 Example 1: Convergence test using manufactured solutions

The first test consists of approximating closed-form solutions on a two-
dimensional domain $\Omega=(0,1)^{2}$. We construct the forcing term
$\boldsymbol{f}$ so that the exact solution to (2.5a)-(2.5c) is given by the
following smooth functions

$\displaystyle
p(x,y):=\left(\left(x-\dfrac{1}{2}\right)^{3}y^{2}+(1-x)^{3}\left(y-\dfrac{1}{2}\right)^{3}\right),$
$\displaystyle\boldsymbol{u}(x,y):=\mathop{\mathbf{curl}}\nolimits(1000x^{2}(1-x)^{4}y^{3}(1-y)^{2}),\qquad\omega(x,y):=\mathop{\mathbf{curl}}\nolimits\boldsymbol{u},$

which satisfy the incompressibility constraint as well as the boundary
conditions. In addition, we take $\boldsymbol{\beta}=\boldsymbol{u},$ and two
specifications for the variable viscosity are considered,

$\nu_{a}(x,y)=\nu_{0}+(\nu_{1}-\nu_{0})xy,\
\nu_{b}(x,y)=\nu_{0}+(\nu_{1}-\nu_{0})exp(-10^{13}((x-0{\text{.}}5)^{10}+(y-0{\text{.}}5)^{10})),$

with $\nu_{0}=0{\text{.}}001$, $\nu_{1}=1$, and taking
$\kappa_{1}=\frac{2}{3}\nu_{0}$, $\kappa_{2}=\frac{\nu_{0}}{2}$ and
$\sigma=100$. The error history of the method introduced in Section 3.1.1 with
discontinuous finite elements for vorticity $(\mathbf{W}_{h}^{2})$ for $k=1$
and for the two different viscosity functions is collected in Tables 5.1 and
5.2, respectively. These values indicate optimal accuracy $O(h^{2})$ for
$k=1$, and for $\nu_{a}$ and $\nu_{b}$, according to Theorem 3.2.

$h$ | $|\\!|\\!|\boldsymbol{u}-\boldsymbol{u}_{h}|\\!|\\!|_{1,\Omega}$ | $r(\boldsymbol{u})$ | $\|\omega-\omega_{h}\|_{0,\Omega}$ | $r(\omega)$ | $\|p-p_{h}\|_{0,\Omega}$ | $r(p)$
---|---|---|---|---|---|---
0.7071 | 10.86 | – | 9.1110 | – | 2.5470 | –
0.3536 | 4.4240 | 1.3 | 3.5500 | 1.4 | 1.5330 | 0.7
0.1768 | 1.2540 | 1.8 | 0.9854 | 1.9 | 0.3493 | 2.1
0.0883 | 0.3492 | 1.8 | 0.2470 | 2.0 | 0.0622 | 2.4
0.0441 | 0.1096 | 1.7 | 0.0613 | 2.0 | 0.0107 | 2.5
0.0221 | 0.0327 | 1.8 | 0.0151 | 2.0 | 0.0020 | 2.4
0.0110 | 0.0075 | 2.1 | 0.0037 | 2.0 | 0.0004 | 2.2
Table 5.1: Example 1: convergence tests against analytical solutions on a sequence of uniformly refined triangulations of the domain $\Omega$ and the viscosity function $\nu_{a}$. $h$ | $|\\!|\\!|\boldsymbol{u}-\boldsymbol{u}_{h}|\\!|\\!|_{1,\Omega}$ | $r(\boldsymbol{u})$ | $\|\omega-\omega_{h}\|_{0,\Omega}$ | $r(\omega)$ | $\|p-p_{h}\|_{0,\Omega}$ | $r(p)$
---|---|---|---|---|---|---
0.7071 | 10.91 | – | 9.1340 | – | 2.1190 | –
0.3536 | 4.489 | 1.3 | 3.6710 | 1.3 | 1.4580 | 0.5
0.1768 | 1.367 | 1.7 | 1.1200 | 1.7 | 0.2789 | 2.4
0.0883 | 0.366 | 1.9 | 0.2951 | 1.9 | 0.0482 | 2.5
0.0441 | 0.113 | 1.7 | 0.0864 | 1.8 | 0.0070 | 2.8
0.0221 | 0.036 | 1.6 | 0.0220 | 2.0 | 0.0014 | 2.3
0.0110 | 0.007 | 2.1 | 0.0046 | 2.2 | 0.0003 | 2.2
Table 5.2: Example 1: convergence tests against analytical solutions on a
sequence of uniformly refined triangulations of the domain $\Omega$ and the
viscosity function $\nu_{b}$.

###  5.2 Example 2: Convergence in 3D

The aim of this numerical test is to assess the accuracy of the method in the
3D case. With this end, we consider the domain $\Omega:=(0,1)^{3}$ and take
$\boldsymbol{f}$ so that the exact solution is given by

$\displaystyle
p(x,y,z):=1-x^{2}-y^{2}-z^{2},\quad\varphi(x,y,z):=x^{2}(1-x)^{2}y^{2}(1-y)^{2}z^{2}(1-z)^{2},$
$\displaystyle\boldsymbol{u}(x,y,z)=\mathop{\mathbf{curl}}\nolimits\varphi,\qquad\boldsymbol{\omega}(x,y,z)=\mathop{\mathbf{curl}}\nolimits\boldsymbol{u},$

and we consider $\boldsymbol{\beta}=\boldsymbol{u}$, and
$\nu_{c}(x,y,z)=\nu_{0}+(\nu_{1}-\nu_{0})x^{2}y^{2}z^{2}$. The remaining
constants are $\nu_{0}=0{\text{.}}1$, $\nu_{1}=1$,
$\kappa_{1}=\frac{2}{3}\nu_{0}$, $\kappa_{2}=\frac{\nu_{0}}{2}$, and
$\sigma=1000$. We observe that the hypothesis of Lemma 2.2 are satisfied.
Additionally, we employ finite elements with $k=1$, that is, $\mathbf{V}_{h}$
approximating the velocity, and piecewise linear and continuous elements for
vorticity and pressure.

In Table 5.3, we summarise the convergence history for a sequence of uniform
meshes. For velocity we observe the $O(h)$ convergence predicted by Theorem
3.4, whereas the approximation of vorticity and pressure seem to be
superconvergent. Figure 5.1 displays velocity and vorticity streamlines as
well as the approximate pressure distribution.

t

$h$ | $|\\!|\\!|\boldsymbol{u}-\boldsymbol{u}_{h}|\\!|\\!|_{1,\Omega}$ | $r(\boldsymbol{u})$ | $\|\boldsymbol{\omega}-\boldsymbol{\omega}_{h}\|_{0,\Omega}$ | $r(\boldsymbol{\omega})$ | $\|p-p_{h}\|_{0,\Omega}$ | $r(p)$
---|---|---|---|---|---|---
0.866 | 0.01021 | – | 0.00299 | – | 0.04732 | –
0.433 | 0.00858 | 0.3 | 0.00125 | 1.3 | 0.01399 | 1.8
0.288 | 0.00665 | 0.6 | 0.00067 | 1.5 | 0.00572 | 2.2
0.216 | 0.00513 | 0.9 | 0.00043 | 1.5 | 0.00290 | 2.4
0.173 | 0.00398 | 1.1 | 0.00030 | 1.5 | 0.00171 | 2.4
0.144 | 0.00313 | 1.3 | 0.00023 | 1.5 | 0.00112 | 2.3
0.123 | 0.00251 | 1.3 | 0.00018 | 1.5 | 0.00079 | 2.2

Table 5.3: Example 2: experimental convergence using homogeneous Dirichlet
boundary conditions on a 3D domain $\Omega$ and using the viscosity function
$\nu_{c}$.

Figure 5.1: Example 2: Approximate solutions computed using the MINI-element.
Velocity streamlines (left) vorticity streamlines (centre) and pressure
distribution (right).

###  5.3 Example 3: A posteriori error estimates and adaptive mesh refinement

In this numerical test, we test the efficiency of the a posteriori error
estimator (4.22) and applying mesh refinement according to the local value of
the indicator. In this case, the convergence rates are obtained by replacing
the expression $\log(h/\hat{h})$ appearing in the computation of (5.27) by
$-\frac{1}{2}\log(N/\hat{N})$, where $N$ and $\hat{N}$ denote the
corresponding degrees of freedom of each triangulation.

Figure 5.2: Example 3: Snapshots of four grids, $\mathcal{T}_{h}^{1}$,
$\mathcal{T}_{h}^{4}$, $\mathcal{T}_{h}^{6}$, $\mathcal{T}_{h}^{10}$,
adaptively refined according to the a posteriori error indicator defined in
(4.22).

Now, we recall the definition of the so-called effectivity index as the ratio
between the total error and the global error estimator, i.e.,

${\tt
e}(\boldsymbol{u},\omega,p):=\Big{\\{}[e(\boldsymbol{u})]^{2}+[e(\omega)]^{2}+[e(p)]^{2}\Big{\\}}^{1/2},\qquad{\tt
eff}(\Theta):=\frac{{\tt e}(\boldsymbol{u},\omega,p)}{\Theta}.$

We will employ the family of finite elements introduced in Section 3.1.1 for
$k=1$, namely piecewise quadratic and continuous elements for velocity and
piecewise linear and continuous elements for vorticity and pressure fields.

The computational domain is the nonconvex L-shaped domain
$\Omega=(-1,1)^{2}\setminus(0,1)^{2}$, where problem (2.5a)-(2.5c) admits the
following exact solution

$\displaystyle
p(x,y):=\dfrac{1-x^{2}-y^{2}}{(x-0.025)^{2}+(y-0.025)^{2}}-12.742942014/3,$
$\displaystyle\varphi(x,y)=x^{2}(1-x)^{2}y^{2}(1-y)^{2}\exp(-50((x-0.025)^{2}+(y-0.025)^{2})),\quad\boldsymbol{u}=\mathop{\mathbf{curl}}\nolimits\varphi,\quad\omega=\mathop{\mathbf{curl}}\nolimits\boldsymbol{u},$

which satisfy the incompressibility constraint as well as the boundary
conditions. Convective velocity, viscosity, and other parameters are taken as

$\displaystyle\boldsymbol{\beta}=\boldsymbol{u},\quad\nu_{d}(x,y)=\nu_{0}+\dfrac{721}{16}(\nu_{1}-\nu_{0})x^{2}(1-x)y^{2}(1-y),\quad\nu_{0}=0.1,\quad\nu_{1}=1,$
$\displaystyle\nu_{e}(x,y)=\nu_{0}+(\nu_{1}-\nu_{0})\exp(-10^{12}((x-0.5)^{10}+(y-0.5)^{10})),\quad\kappa_{1}=\frac{2}{3}\nu_{0},\quad\kappa_{2}=\frac{\nu_{0}}{2},\quad\sigma=10.$

Pressure is singular near the reentrant corner of the domain and so we expect
hindered convergence of the approximations when a uniform (or quasi-uniform)
mesh refinement is applied. In contrast, if we apply the following adaptive
mesh refinement procedure from [49]:

  1. 1)

Start with a coarse mesh $\mathcal{T}_{h}$.

  2. 2)

Solve the discrete problem (3.17) for the current mesh $\mathcal{T}_{h}$.

  3. 3)

Compute $\Theta_{T}:=\Theta$ for each triangle $T\,\in\,\mathcal{T}_{h}$.

  4. 4)

Check the stopping criterion and decide whether to finish or go to next step.

  5. 5)

Use blue-green refinement on those $T^{\prime}\in\mathcal{T}_{h}$ whose
indicator $\Theta_{T^{\prime}}$ satisfies

$\Theta_{T^{\prime}}\,\geq\,\frac{1}{2}\,\max_{T\in\mathcal{T}_{h}}\left\\{\Theta_{T}:\,\,T\in\mathcal{T}_{h}\,\right\\}\,.$

  6. 6)

Define resulting meshes as current meshes $\mathcal{T}_{h}$ and
$\mathcal{T}_{h}$, and go to step 2,

we expect a recovering of the optimal convergence rates. In fact, this can be
observed from the bottom rows of Tables 5.4 and 5.5, for both $\nu_{d}$ and
$\nu_{e}$, respectively. Moreover, the efficiency indexes are around 1 for
both viscosities. The resulting meshes after a few adaptation steps are
reported in Figure 5.2, showing the expected refinement near the reentrant
corner.

$N$ | $|\\!|\\!|\boldsymbol{u}-\boldsymbol{u}_{h}|\\!|\\!|_{1,\Omega}$ | $r(\boldsymbol{u})$ | $\|\omega-\omega_{h}\|_{0,\Omega}$ | $r(\omega)$ | $\|p-p_{h}\|_{0,\Omega}$ | $r(p)$ | ${\tt eff}(\Theta)$
---|---|---|---|---|---|---|---
661 | 49.68 | – | 8.821 | – | 6.685 | – | 1.133
999 | 32.37 | 2.07 | 5.069 | 2.68 | 3.985 | 2.50 | 1.157
1241 | 15.46 | 6.81 | 2.104 | 8.10 | 1.846 | 7.09 | 1.144
1881 | 9.058 | 2.57 | 1.396 | 1.97 | 1.057 | 2.68 | 1.098
2103 | 7.178 | 4.17 | 0.907 | 7.72 | 0.828 | 4.36 | 1.135
2621 | 5.645 | 2.18 | 0.754 | 1.67 | 0.655 | 2.12 | 1.120
3851 | 3.647 | 2.27 | 0.454 | 2.63 | 0.418 | 2.33 | 1.168
4267 | 3.243 | 2.29 | 0.401 | 2.46 | 0.365 | 2.61 | 1.156
5271 | 2.687 | 1.77 | 0.298 | 2.76 | 0.294 | 2.03 | 1.143
7819 | 1.754 | 2.16 | 0.194 | 2.18 | 0.191 | 2.22 | 1.155
Table 5.4: Example 3: Convergence history and effectivity indexes for the method introduced in Section 3.1.1, computed on a sequence of adaptively refined triangulations of the L-shaped domain and using viscosity $\nu_{d}$. $N$ | $|\\!|\\!|\boldsymbol{u}-\boldsymbol{u}_{h}|\\!|\\!|_{1,\Omega}$ | $r(\boldsymbol{u})$ | $\|\omega-\omega_{h}\|_{0,\Omega}$ | $r(\omega)$ | $\|p-p_{h}\|_{0,\Omega}$ | $r(p)$ | ${\tt eff}(\Theta)$
---|---|---|---|---|---|---|---
661 | 49.73 | – | 8.842 | – | 6.681 | – | 1.132
999 | 32.39 | 2.07 | 5.081 | 2.68 | 3.980 | 2.50 | 1.155
1241 | 15.50 | 6.79 | 2.122 | 8.05 | 1.838 | 7.12 | 1.138
1881 | 9.087 | 2.56 | 1.401 | 1.99 | 1.039 | 2.74 | 1.085
2103 | 7.213 | 4.14 | 0.914 | 7.65 | 0.806 | 4.55 | 1.114
2589 | 5.683 | 2.29 | 0.759 | 1.78 | 0.633 | 2.32 | 1.112
3771 | 3.734 | 2.23 | 0.461 | 2.64 | 0.406 | 2.35 | 1.113
5161 | 2.674 | 2.12 | 0.307 | 2.58 | 0.287 | 2.20 | 1.108
6867 | 1.946 | 2.22 | 0.207 | 2.77 | 0.205 | 2.36 | 1.116
9887 | 1.346 | 2.02 | 0.128 | 2.60 | 0.138 | 2.16 | 1.119
Table 5.5: Example 3: Convergence history and effectivity indexes for the
method introduced in Section 3.1.1, computed on a sequence of adaptively
refined triangulations of the L-shaped domain and using viscosity $\nu_{e}$.

###  5.4 Example 4: Steady blood flow in aortic arch

We finalise the set of examples with a simple simulation of pseudo-stationary
blood flow in an aorta. The patient-specific geometry [38, 39] has one inlet
(a segment that connects with the pre-aortic root coming from the aortic valve
in the heart) and four outlets (the left common carotid artery, the left
subclavian artery, the innominate artery, and the larger descending aorta). On
the inlet we impose a Poiseuille profile of magnitude 4, on the vessel walls
we set no-slip conditions, and on the remaining boundaries we set zero normal
stresses (more physiologically relevant boundary conditions can be considered
following, e.g., [24, 32]). The initial unstructured mesh has 46352
tetrahedral elements. The synthetic variable viscosity field is a smooth
exponential function
$\nu=\nu_{0}+(\nu_{1}-\nu_{0})\exp(-10^{3}[(x-0.1)^{6}+y-0.5)^{6}+(z-0.5)^{6}])$
with $\nu_{0}=10^{-3},\nu_{1}=10$ that entails an average Reynolds number of
approximately 60 (computed using the inlet diameter and maximal inlet
velocity), while the convecting velocity is computed as the solution of a
preliminary Stokes problem (on the initial coarse mesh), and we prescribe
$\sigma=1000$ and $\boldsymbol{f}=\sigma\boldsymbol{\beta}$. Then we compute
numerical solutions of the Oseen problem and apply four steps of adaptive mesh
refinement using a 3D version of the estimator (4.22) and the algorithm
described in the previous example. The results are portrayed in Figure 5.3,
plotting pressure distribution, velocity streamlines, vorticity, and a sample
of the resulting adaptive mesh which shows more refinement near the boundaries
of the descending aorta. For this test we have used a conforming approximation
of vorticity.

Figure 5.3: Example 4: Simulation of stationary blood flow in an aortic arch.
Approximate velocity, vorticity, and pressure (top panels), and samples of
adaptive mesh after one, two and three refinement steps, and visualising a cut
focusing on the boundaries (bottom row).

## References

  * [1] M. Ainsworth and J.T. Oden, A posteriori error estimation in finite element analysis. Wiley, New York, 2000.
  * [2] M. Alvarez, G.N. Gatica, and R. Ruiz-Baier, A vorticity-based fully-mixed formulation for the 3D Brinkman-Darcy problem, Comput. Methods Appl. Mech. Engrg., 307 (2016) 68–95.
  * [3] M. Amara, D. Capatina-Papaghiuc, and D. Trujillo, Stabilized finite element method for Navier-Stokes equations with physical boundary conditions, Math. Comp., 76(259) (2007) 1195–1217.
  * [4] M. Amara, E. Chacón Vera, and D. Trujillo, Vorticity–velocity–pressure formulation for Stokes problem, Math. Comp., 73(248) (2004) 1673–1697.
  * [5] V. Anaya, A. Bouharguane, D. Mora, C. Reales, R. Ruiz-Baier, N. Seloula and H. Torres, Analysis and approximation of a vorticity-velocity-pressure formulation for the Oseen equations, J. Sci. Comput., 88(3) (2019) 1577–1606.
  * [6] V. Anaya, B. Gómez-Vargas, D. Mora and R. Ruiz-Baier, Incorporating variable viscosity in vorticity-based formulations for Brinkman equations, C. R. Math. Acad. Sci. Paris, 357(6) (2019) 552–560.
  * [7] V. Anaya, D. Mora, R. Oyarzúa, and R. Ruiz-Baier, A priori and a posteriori error analysis of a mixed scheme for the Brinkman problem, Numer. Math., 133(4) (2016) 781-817.
  * [8] V. Anaya, D. Mora, and R. Ruiz-Baier, An augmented mixed finite element method for the vorticity-velocity-pressure formulation of the Stokes equations, Comput. Methods Appl. Mech. Engrg., 267 (2013) 261–274.
  * [9] T.P. Barrios, J.M. Cascón and M. González, Augmented mixed finite element method for the Oseen problem: a priori and a posteriori error analyses, Comput. Methods Appl. Mech. Engrg., 313 (2017) 216–238.
  * [10] M. Benzi, and M.A. Olshanskii, An augmented Lagrangian-based approach to the Oseen problem, SIAM J. Sci. Comput., 28(6) (2006) 2095–2113.
  * [11] M. Benzi, M.A. Olshanskii, L.G. Rebholz and Z. Wang, Assessment of a vorticity based solver for the Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg., 247-248 (2012) 216–225.
  * [12] C. Bernardi, and N. Chorfi, Spectral discretization of the vorticity, velocity, and pressure formulation of the Stokes problem, SIAM J. Numer. Anal., 44(2) (2006) 826–850.
  * [13] P.V. Bochev, Analysis of least-squares finite element methods for the Navier–Stokes equations, SIAM J. Numer. Anal., 34(5) (1997) 1817–1844.
  * [14] P.V. Bochev, and M. Gunzburger, Least-Squares finite element methods. Volume 166 of Applied Mathematical Sciences. Springer Verlag (2009).
  * [15] D. Boffi, Stability of higher order triangular Hood–Taylor methods for stationary Stokes equations, Math. Models Methods Appl. Sci., 2(4) (1994) 223–235.
  * [16] D. Boffi, F. Brezzi, and M. Fortin, Mixed Finite Element Methods and applications. Springer Series in Computational Mathematics, 44. Springer, Heidelberg (2013).
  * [17] Z. Cai and B. Chen, Least-squares method for the Oseen equation, Numer. Methods Partial Differential Equations, 32 (2016) 1289–1303.
  * [18] J. Camaño, R. Oyarzúa, R. Ruiz-Baier, and G. Tierra, Error analysis of an augmented mixed method for the Navier-Stokes problem with mixed boundary conditions, IMA J. Numer. Anal., 38(3) (2018) 1452–1484.
  * [19] J.O. Campos, R. Weber dos Santos, J. Sundnes, and B. Martins Rocha, Preconditioned augmented Lagrangian formulation for nearly incompressible cardiac mechanics, Int. J. Numer. Methods Biomed. Engrg., 34(4) (2018) e2948.
  * [20] C.L. Chang, and B.-N. Jiang, An error analysis of least-squares finite element method of velocity-pressure-vorticity formulation for the Stokes problem, Comput. Methods Appl. Mech. Engrg., 84(3) (1990) 247–255.
  * [21] C.L. Chang and S.-Y. Yang, Analysis of the $[L^{2},L^{2},L^{2}]$ least-squares finite element method for incompressible Oseen-type problems, Int. J. Numer. Anal. Model., 4(3-4) (2007) 402–424.
  * [22] P. Clément, Approximation by finite element functions using local regularisation, RAIRO Modél. Math. Anal. Numer., 9 (1975) 77–84.
  * [23] B. Cockburn, and J. Cui, An analysis of HDG methods for the vorticity-velocity-pressure formulation of the Stokes problem in three dimensions, Math. Comp., 81(279) (2012) 1355–1368.
  * [24] T. Dobroserdova, F. Liang, G. Panasenko, and Y. Vassilevski, Multiscale models of blood flow in the compliant aortic bifurcation, Appl. Math. Letters, 93C (2019) 98–104.
  * [25] H.-Y. Duan, and G.-P. Liang, On the velocity-pressure-vorticity least-squares mixed finite element method for the 3D Stokes equations, SIAM J. Numer. Anal., 41(6) (2003) 2114–2130.
  * [26] F. Dubois, M. Salaün, and S. Salmon, First vorticity-velocity-pressure numerical scheme for the Stokes problem, Comput. Methods Appl. Mech. Engrg., 192(44–46) (2003) 4877–4907.
  * [27] A. Ern, Vorticity-velocity formulation of the Stokes problem with variable density and viscosity, Math. Models Methods Appl. Sci., 8(2) (1998) 203–218.
  * [28] A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements. Applied Mathematical Sciences. Springer-Verlag, (2004).
  * [29] P.E. Farrell, L.F. Gatica, B.P. Lamichhane, R. Oyarzúa, and R. Ruiz-Baier, Mixed Kirchhoff stress - displacement - pressure formulations for incompressible hyperelasticity, Comput. Methods Appl. Mech. Engrg., 374 (2021) e113562.
  * [30] P.E. Farrell, and P.A. Gazca-Orozco, An augmented Lagrangian preconditioner for implicitly constituted non-Newtonian incompressible flow, SIAM J. Sci. Comput., 42(6) (2020) B1329–B1349.
  * [31] P.E. Farrell, L. Mitchell, and F. Wechsung, An augmented Lagrangian preconditioner for the 3D stationary incompressible Navier-Stokes equations at high Reynolds number, SIAM J. Sci. Comput., 41(5) (2019) A3073–A3096.
  * [32] L Formaggia, J.F. Gerbeau, F. Nobile, and A. Quarteroni, Numerical treatment of defective boundary conditions for the Navier-Stokes equations, SIAM J. Numer. Anal., 40 (2002) 376–401.
  * [33] G.N. Gatica, A Simple Introduction to the Mixed Finite Element Method. Theory and Applications. Springer Briefs in Mathematics, Springer, Cham Heidelberg New York Dordrecht London (2014).
  * [34] G.N. Gatica, R. Oyarzúa, and N. Valenzuela, A five-field augmented fully-mixed finite element method for the Navier-Stokes/Darcy coupled problem, Comput. Math. Appl., 80(8) (2020) 1944–1963.
  * [35] V. Girault, and P.A. Raviart, Finite element methods for Navier-Stokes equations. Theory and Algorithms. Springer-Verlag, Berlin (1986).
  * [36] P. Hood, and C. Taylor, Numerical solution of the Navier–Stokes equations using the finite element technique, Comput. Fluids, 1 (1973) 1–28
  * [37] V. John, K. Kaiser, and J. Novo, Finite element methods for the incompressible Stokes equations with variable viscosity, ZAMM Z. Angew. Math. Mech., 96(2) (2016) 205–216.
  * [38] E. Marchandise, P. Crosetto, C. Geuzaine, J.-F. Remacle, and E. Sauvage, Quality open source mesh generation for cardiovascular flow simulation. In: D. Ambrosi, A. Quarteroni, and G. Rozza, editors. Modeling of Physiological Flows. Milano: Springer (2011) 395–414.
  * [39] E. Marchandise, E. Sauvage, and J.-F. Remacle, Quality open source mesh generation for biological flow simulations. In: Proceedings of the ECCOMAS Thematic International Conference on Simulation and Modeling of Biological Flows (SIMBIO 2011) September 21-23, 2011, VUB, Brussels, Belgium.
  * [40] S. Mohapatra, and S. Ganesan, A non-conforming least squares spectral element formulation for Oseen equations with applications to Navier-Stokes equations, Numer. Funct. Anal. Optim., 37(10) (2016) 295–1311.
  * [41] P.R. Patil, and G. Vaidyanathan, Effect of variable viscosity on thermohaline convection in a porous medium, J. Hydrology, 57(1-2) (1982) 147–161.
  * [42] L.E. Payne, J.C. Song, and B. Straughan, Continuous dependence and convergence results for Brinkman and Forchheimer models with variable viscosity, Proc. Royal Soc. London Series A, 455 (1986) 1–20.
  * [43] J.P. Pontaza, and J.N. Reddy, Spectral/hp least-squares finite element formulation for the Navier-Stokes equations, J. Comput. Phys., 190(2) (2003) 523–549.
  * [44] J. Rudi, G. Stadler, and O. Ghattas, Weighted BFBT preconditioner for Stokes flow problems with highly heterogeneous viscosity, SIAM J. Sci. Comput., 39(5) (2017) S272–S297.
  * [45] M. Salaün, and S. Salmon, Numerical stabilization of the Stokes problem in vorticity-velocity-pressure formulation, Comput. Methods Appl. Mech. Engrg., 196(9-12) (2007) 1767–1786.
  * [46] M. Salaün, and S. Salmon, Low-order finite element method for the well-posed bidimensional Stokes problem, IMA J. Numer. Anal., 35 (2015) 427–453.
  * [47] C.G. Speziale, On the advantages of the vorticity-velocity formulations of the equations of fluid dynamics, J. Comput. Phys., 73(2) (1987) 476–480.
  * [48] C.-C. Tsai, and S.-Y. Yang, On the velocity-vorticity-pressure least-squares finite element method for the stationary incompressible Oseen problem, J. Comput. Appl. Math., 182(1) (2005) 211–232.
  * [49] R. Verfürth, A review of a posteriori error estimation and adaptive-mesh-refinement techniques. Wiley-Teubner, Chichester (1996).