2023 1

CARABALLO et al. Biot–Brinkman equations in vorticity form

Corresponding author Ricardo Ruiz-Baier.

Chilean Research and Technology Council through the ANID program for
international research visits; by the Australian Government through the
National Computational Infrastructure (NCI) under the ANU Merit Allocation
Scheme (ANUMAS); and by the Australian Research Council through the Future
Fellowship grant FT220100496 and Discovery Project grant DP22010316.

# Robust finite element methods and solvers for the Biot–Brinkman equations in
vorticity form

Ruben Caraballo  Chansophea Wathanak In  Alberto F. Martín  Ricardo Ruiz-Baier
GIMNAP, Departamento de Matemática, Universidad del Bío-Bío, Concepción, Chile
School of Mathematics, Monash University, 9 Rainforest Walk, Melbourne VIC
3800, Australia  School of Computing, Australian National University, Acton
ACT 2601, Australia  Universidad Adventista de Chile, Casilla 7-D, Chillán,
Chile<EMAIL_ADDRESS>Caraballo R  In CW  Martín AF  Ruiz-Baier
R

(Date Month Year; Date Month Year)

###### Abstract

[Abstract]In this paper, we propose a new formulation and a suitable finite
element method for the steady coupling of viscous flow in deformable porous
media using divergence-conforming filtration fluxes. The proposed method is
based on the use of parameter-weighted spaces, which allows for a more
accurate and robust analysis of the continuous and discrete problems.
Furthermore, we conduct a solvability analysis of the proposed method and
derive optimal error estimates in appropriate norms. These error estimates are
shown to be robust in a variety of regimes, including the case of large Lamé
parameters and small permeability and storativity coefficients. To illustrate
the effectiveness of the proposed method, we provide a few representative
numerical examples, including convergence verification and testing of
robustness of block-diagonal preconditioners with respect to model parameters.

, , , and . Robust finite element methods and solvers for the Biot–Brinkman
equations in vorticity form. Numer Methods PDEs. 2023;00(00):1–18.

###### keywords:

Biot–Brinkman coupled problem, deformable porous media, vorticity-based
formulation, mixed finite element methods.
Mathematics Subject Classifications (2000): 65N30, 65N15, 76S05, 35Q74.

††articletype: Article Type††journal: Numerical Methods for Partial
Differential Equations††volume: 00

##  1 Introduction

We address the analysis of the Biot–Brinkman equations, which serve as a model
for filtration of viscous flow in deformable porous media 3, 22, 54, 55. The
system has been recently analysed in 35 for the case of multiple network
poroelasticity, using
$\mathbf{H}(\mathop{\mathrm{div}}\nolimits,\Omega)$-conforming displacements
and filtration fluxes (or seepage velocities) for each compartment, also
designing robust preconditioners. Here we propose a reformulation for only one
fluid compartment but using the vorticity field (defined as the curl of the
filtration velocity) as an additional unknown in the system, and we also
include the total pressure, following 50, 42. Such an approach enables us to
avoid the notorious problem of locking or non-physical pressure oscillations
when approximating poroelastic models and it has led to a number of
developments including extensions to multiple network models, interfacial
free-flow and poromechanics coupling, nonlinear interaction with species
transport, reformulations into four and more fields systems, and using other
discretisations such as discontinuous Galerkin, nonconforming FEM, weak
Galerkin, and virtual elements. See, e.g., 16, 20, 36, 38, 41, 43, 51, 52, 56,
58, 59.

The formulation of viscous flow equations using vorticity, velocity and
pressure has been used and analysed (in terms of solvability of the continuous
and discrete formulations and deriving error estimates) extensively in, e.g.,
1, 2, 7, 8, 6, 9, 14, 23, 27, 28, 57. Methods based on vorticity formulations
are useful for visualisation of rotational flows and they are convenient when
dealing with rotation-based boundary conditions. The coupling with other
effects such as mass and energy transport has also been addressed, see for
example 4, 15, 44. These contributions include fully mixed finite elements,
augmented forms, spectral methods, and Galerkin least-squares stabilised types
of discretisations. At the continuous level, one appealing property of some of
these vorticity formulations is that the $\mathbf{H}^{1}$-conformity of the
filtration flux is relaxed and velocity is sought in
$\mathbf{H}(\mathop{\mathrm{div}}\nolimits,\Omega)$ and the vorticity is
sought in either $\mathbf{H}(\mathop{\mathbf{curl}}\nolimits,\Omega)$ or
$\mathbf{L}^{2}(\Omega)$. Then, using simply a _conforming_ method, resulting
discretisations are readily mass conservative. Also, mixed methods that look
for both vorticity and velocity will typically deliver approximate vorticity
_with the same accuracy_ as velocity (as opposed to schemes needing numerical
differentiation to get vorticity from approximate velocity). Another advantage
of using a vorticity-based method is that an exactness sequence exists between
the spaces for velocity, vorticity, and pressure, yielding a framework that
can be straightforwardly analysed with finite element exterior calculus 34. We
also stress that solving directly for vorticity enables a more thorough
investigation of vortical patterns, interactions at boundary layers, and the
formation of turbulent eddies, all of which are crucial in a diverse array of
applications.

Note however that, differently from the aforementioned works, in the case of
the Biot–Brinkman problem, the divergence of the fluid velocity is not zero
(or a prescribed fluid source), but it depends on the velocity of the solid
and on the rate of change of fluid pressure. Also, an additional term of
grad–div type appears in the momentum balance for the fluid.

In this paper, we prove the well-posedness of the continuous and discrete
formulations for the coupling of mechanics and fluid flow in fluid-saturated
deformable porous media using Banach–Nečas–Babuška theory in parameter-
weighted Hilbert spaces. The appropriate choice of weighting parameters yields
automatically a framework for robust operator preconditioning in the Biot
equations, following the approach from 42. This operator scaling yields
robustness with respect to the elastic parameters, storativity, Biot–Willis
coefficient, and with respect to permeability. For the Brinkman component, our
present formulation is such that the filtration flux terms have a different
weight in their $\mathbf{L}^{2}(\Omega)$ and
$\mathbf{H}(\mathop{\mathrm{div}}\nolimits,\Omega)$ contributions, which
requires a different treatment for the analysis of the pore pressure terms. To
address this issue it suffices to appeal to the recent theory in 13 (see also
46), which was developed for Darcy equations using non-standard sum spaces,
and we appropriately modify the scalings in the momentum equation. This
modification entails the use of (discontinuous) Laplacian operators in the
fluid pressure preconditioning.

Our proposed approach also offers a novel contribution to the field of
operator preconditioning for the interaction of mechanics and fluid flow in
fluid-saturated deformable porous media, which are challenging to solve, and
the design of efficient preconditioners is highly problem-dependent 29.
Previous works have explored the use of block-diagonal preconditioners, Schur
complements, and pressure-correction methods, which have improved the
convergence rate and computational efficiency of numerical solutions for
poromechanics problems 24, 32, 37, 45. In this work, we also derive parameter-
robust solvers, but following 35 and also 38, 43, 51. Our results confirm
that, additionally to the Laplacian contribution needed in the Riesz
preconditioner associated with the fluid pressure mentioned above, we also
require off-diagonal contributions in the total pressure and fluid pressure
coupling terms (as employed in 16). The overall parameter scalings that we
propose are motivated by the stability analysis and we verify computationally
that robustness holds for this particular choice. Note that we only discuss
one type of boundary conditions, but the extension to other forms can be
adapted accordingly.

Research on preconditioning techniques for advanced discretisations of block
multiphysics systems has also crystallised in a number of high quality open
source software packages. One of the earliest efforts was the BKPIT C++
package 25 which, following an object-oriented approach, provides an
extensible framework for the implementation of algebraic block
preconditioners, such as block Jacobi or block Gauss–Seidel. More recently, as
the field of physics-based and discretisation tailored preconditioners has
evolved with breakthrough inventions in approximate block factorisation and
Schur-complement methods towards ever faster and scalable iterative solvers
for large-scale systems, more sophisticated block preconditioning software is
available. The widely-used PETSc package offers the PCFIELDSPLIT subsystem 19
to design and compose complex block preconditioners. The Firedrike library
extends PETSc block preconditioning capabilities and algebraic composability
further 40. Another tightly-related effort is the Teko Trilinos package 26,
which provides a high-level interface to compose block preconditioners using a
functional programming style in C++. The authors in 10 present a generic
software framework in object-oriented Fortran to build block recursive
algebraic factorisation preconditioners for double saddle-point systems, as
those arising in MagnetoHydroDynamics (MHD).

In this paper, we build upon the high momentum gained in the last years by the
Julia programming language for scientific and numerical computing. In
particular, the realisation of the numerical discretisation and
preconditioning algorithms is conducted with the Gridap finite element
software package 11. We leverage the flexibility of this framework, and its
composability with others in the Julia package ecosystem, such as
LinearOperators 49, to prototype natural Riesz map preconditioners in the sum
spaces described above, leading to complex multiphysics coupling solvers that
are robust with respect to physical parameters variations and mesh resolution.
For the sake of reproducibility, the Julia software used in this paper is
available publicly/openly at 21.

The remainder of this article is organised as follows. The presentation of the
new form of Biot–Brinkman equations and its weak formulation are given in
Section 2. The modification of the functional structure to include parameter
weights and the unique solvability analysis for the continuous problem are
addressed in Section 3. The definition of the finite element discretisation
and the specification of the well-posedness theory for the discrete problem is
carried out in Section 4. The error analysis (tailored for a specific family
of finite elements but applicable to other combinations of discrete spaces as
well) is detailed also in that section. Numerical experiments are collected in
Section 5 and we close in Section 6 with a brief summary and a discussion on
possible extensions.

##  2 Model problem and its weak formulation

###  2.1 Preliminaries

Let us consider a simply connected bounded and Lipschitz domain
$\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$ occupied by a poroelastic
domain with one incompressible fluid network incorporating viscosity. The
domain boundary is denoted as $\Gamma:=\partial\Omega$. Throughout the text,
given a normed space $S$, by $\mathbf{S}$ and $\mathbb{S}$ we will denote the
vector and tensor extensions, $S^{d}$ and $S^{d\times d}$, respectively. In
addition, by $\mathrm{L}^{2}(\Omega)$ we will denote the usual Lebesgue space
of square integrable functions and $\mathrm{H}^{m}(\Omega)$ denotes the usual
Sobolev space with weak derivatives of order up to $m\geq 0$ in
$\mathrm{L}^{2}(\Omega)$, and use the convention that
$\mathrm{H}^{0}(\Omega)=\mathrm{L}^{2}(\Omega)$.

Next, we recall the definition of the following Hilbert spaces

$\displaystyle\mathbf{H}^{1+m}(\Omega):=\\{\boldsymbol{\gamma}\in\mathbf{H}^{m}(\Omega):\boldsymbol{\nabla}\boldsymbol{\gamma}\in\mathbb{H}^{m}(\Omega)\\},$
$\displaystyle\mathbf{H}^{m}(\mathop{\mathrm{div}}\nolimits,\Omega)\\!:=\\!\\{\boldsymbol{\zeta}\in\mathbf{H}^{m}(\Omega)\\!:\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\in\mathrm{H}^{m}(\Omega)\\},\quad\mathbf{H}^{m}(\mathop{\mathbf{curl}}\nolimits,\Omega)\\!:=\\!\\{\boldsymbol{\theta}\in\mathbf{H}^{m}(\Omega)\\!:\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\in\mathbf{H}^{m}(\Omega)\\},$

and in the case that $m=0$ the latter two spaces are denoted
$\mathbf{H}(\mathop{\mathrm{div}}\nolimits,\Omega)$ and
$\mathbf{H}(\mathop{\mathbf{curl}}\nolimits,\Omega)$, and we use the following
notation for the typical norms associated with such spaces

$\displaystyle\left\|\boldsymbol{u}\right\|_{1,\Omega}^{2}:=\left\|\boldsymbol{u}\right\|^{2}_{0,\Omega}+\left\|\boldsymbol{\nabla}\boldsymbol{u}\right\|^{2}_{0,\Omega},\qquad\left\|\boldsymbol{v}\right\|_{\mathop{\mathrm{div}}\nolimits,\Omega}^{2}:=\left\|\boldsymbol{v}\right\|^{2}_{0,\Omega}+\left\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\right\|^{2}_{0,\Omega},\qquad\left\|\boldsymbol{\zeta}\right\|_{\mathop{\mathbf{curl}}\nolimits,\Omega}^{2}:=\left\|\boldsymbol{\zeta}\right\|^{2}_{0,\Omega}+\left\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\zeta}\right\|^{2}_{0,\Omega},$

respectively. Furthermore, in view of the boundary conditions on $\Gamma$ (to
be made precise below), we also use the following notation for relevant
subspaces

$\displaystyle\mathbf{H}^{1}_{0}(\Omega)$
$\displaystyle:=\\{\boldsymbol{\gamma}\in\mathbf{H}^{1}(\Omega):\boldsymbol{\gamma}=\boldsymbol{0}\text{
on }\Gamma\\},$
$\displaystyle\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega)$
$\displaystyle:=\\{\boldsymbol{\zeta}\in\mathbf{H}(\mathop{\mathrm{div}}\nolimits,\Omega):\boldsymbol{\zeta}\cdot\boldsymbol{n}=0\text{
on }\Gamma\\},$ (2.1)
$\displaystyle\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega)$
$\displaystyle:=\\{\boldsymbol{\theta}\in\mathbf{H}(\mathop{\mathbf{curl}}\nolimits,\Omega):\boldsymbol{\theta}\times\boldsymbol{n}=\boldsymbol{0}\text{
on }\Gamma\\},$ $\displaystyle\mathrm{L}^{2}_{0}(\Omega)$
$\displaystyle:=\\{\zeta\in\mathrm{L}^{2}(\Omega):\int_{\Omega}\zeta=0\\}.$

For a generic functional space $X$ and a scalar $\eta>0$, the weighted space
$\eta X$ refers to the same $X$ but endowed with the norm $\eta\|\cdot\|_{X}$.
In addition, we recall the definition of the norm of intersection $X\cap Y$
and sum $X+Y$ of Hilbert spaces $X,Y$

$\|z\|^{2}_{X\cap
Y}=\|z\|^{2}_{X}+\|z\|^{2}_{Y},\qquad\|z\|^{2}_{X+Y}=\inf_{\small\begin{array}[]{c}z=x+y\\\
x\in X,y\in Y\end{array}}\|x\|^{2}_{X}+\|y\|^{2}_{Y}.$ (2.2)

###  2.2 The governing equations

The viscous filtration flow through the deformed porous skeleton can be
described by the following form of the Biot–Brinkman equations in steady form,
representing the mixture momentum, fluid momentum, and mixture mass balance,
respectively

$\displaystyle-\mathop{\mathbf{div}}\nolimits(2\mu\boldsymbol{\varepsilon}(\boldsymbol{u})+[\lambda\mathop{\mathrm{div}}\nolimits\boldsymbol{u}-\alpha
p]\mathbf{I})$ $\displaystyle=\boldsymbol{b},$ (2.3a)
$\displaystyle{\frac{1}{\kappa}}\boldsymbol{v}-{\frac{\nu}{\kappa}}\boldsymbol{\Delta}\boldsymbol{v}+{\nabla
p}$ $\displaystyle=\widehat{\boldsymbol{f}},$ (2.3b) $\displaystyle-
c_{0}p-\alpha\mathop{\mathrm{div}}\nolimits\boldsymbol{u}-\mathop{\mathrm{div}}\nolimits\boldsymbol{v}$
$\displaystyle=g,$ (2.3c)

where $\boldsymbol{u}$ is the displacement of the skeleton,
$\boldsymbol{\varepsilon}(\boldsymbol{u})=\frac{1}{2}(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{\tt
t})$ is the tensor of infinitesimal strains, $\boldsymbol{v}$ is the
filtration flux, and $p$ is the pressure head. The parameters are the external
body load $\boldsymbol{b}$, the external force applied on the fluid
$\widehat{\boldsymbol{f}}$, the kinematic viscosity of the interstitial fluid
$\nu$, the hydraulic conductance (permeability field, here assumed a positive
constant) $\kappa$, the Lamé coefficients of the solid structure
$\lambda,\mu$, the storativity $c_{0}$, and the Biot–Willis modulus $\alpha$.
Equations (2.3a) are equipped with boundary conditions of clamped boundary for
the solid phase and slip filtration velocity

$\boldsymbol{u}=\boldsymbol{0}\quad\text{and}\quad\boldsymbol{v}\cdot\boldsymbol{n}=0\quad\text{on}\quad\Gamma.$

Next we introduce the rescaled filtration vorticity vector

$\boldsymbol{\omega}:={\sqrt{\frac{\nu}{\kappa}}}\mathop{\mathbf{curl}}\nolimits\boldsymbol{v},$
(2.4)

which has a different weight than that used in 9 for Brinkman flows. We also
define, following the developments in 50, 43, the additional total pressure
field

$\varphi:=-\lambda\mathop{\mathrm{div}}\nolimits\boldsymbol{u}+\alpha p.$

In order to rewrite the fluid momentum balance in terms of the rescaled
filtration vorticity (2.4), we employ the following vector identity, valid for
a generic vector field $\boldsymbol{v}$:

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

These steps, together with a rescaling of the external force
$\boldsymbol{f}={\nu}\widehat{\boldsymbol{f}}$, lead to the following
equations (mixture momentum, constitutive equation for total pressure, fluid
momentum, constitutive equation for filtration vorticity, and mixture mass
balance)

$\displaystyle-\mathop{\mathbf{div}}\nolimits(2\mu\boldsymbol{\varepsilon}(\boldsymbol{u})-\varphi\mathbf{I})$
$\displaystyle=\boldsymbol{b},$ (2.5a)
$\displaystyle-\frac{1}{\lambda}\varphi+\frac{\alpha}{\lambda}p-\mathop{\mathrm{div}}\nolimits\boldsymbol{u}$
$\displaystyle=0,$ (2.5b)
$\displaystyle{\frac{1}{\kappa}}\boldsymbol{v}+{\sqrt{\frac{\nu}{\kappa}}}\mathop{\mathbf{curl}}\nolimits\boldsymbol{\omega}-{\frac{\nu}{\kappa}}\nabla(\mathop{\mathrm{div}}\nolimits\boldsymbol{v})+\nabla
p$ $\displaystyle=\boldsymbol{f},$ (2.5c)
$\displaystyle-\boldsymbol{\omega}+{\sqrt{\frac{\nu}{\kappa}}}\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}$
$\displaystyle=\boldsymbol{0},$ (2.5d)
$\displaystyle-(c_{0}+\frac{\alpha^{2}}{\lambda})p+\frac{\alpha}{\lambda}\varphi-\mathop{\mathrm{div}}\nolimits\boldsymbol{v}$
$\displaystyle=g;$ (2.5e) and the boundary conditions now read
$\boldsymbol{u}=\boldsymbol{0}\quad\text{and}\quad\boldsymbol{v}\cdot\boldsymbol{n}=0\quad\text{and}\quad\boldsymbol{\omega}\times\boldsymbol{n}=\boldsymbol{0}\qquad\text{on}\quad\Gamma.$
(2.5f)

###  2.3 Weak formulation

Let us assume that $\boldsymbol{b},\boldsymbol{f}\in\mathbf{L}^{2}(\Omega)$,
$g\in\mathrm{L}^{2}(\Omega)$ and that all other model coefficients are
positive constants. We proceed to multiply the governing equations by suitable
test functions and to integrate by parts over the domain. Note that for the
divergence-based terms we appeal to the usual form of the Gauss formula,
whereas for curl-based terms we use the following result from, e.g., 33:

$\int_{\Omega}\mathop{\mathbf{curl}}\nolimits\boldsymbol{\gamma}\cdot\boldsymbol{\theta}=\int_{\Omega}\boldsymbol{\gamma}\cdot\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}+\int_{\partial\Omega}(\boldsymbol{\theta}\times\boldsymbol{n})\cdot\boldsymbol{\gamma}\qquad\forall\boldsymbol{\theta}\in\mathbf{H}(\mathop{\mathbf{curl}}\nolimits,\Omega),\
\boldsymbol{\gamma}\in\mathbf{H}^{1}(\Omega),$

together with the scalar triple product identity

$(\boldsymbol{\theta}\times\boldsymbol{n})\cdot\boldsymbol{v}=\boldsymbol{\theta}\cdot(\boldsymbol{n}\times\boldsymbol{v}).$

We observe in advance that equations (2.5b),(2.5e) suggest that, in the limit
of $c_{0}\to 0,\lambda\to\infty$, both pressures are not uniquely defined and
so we require to add a constraint on their mean value (to be zero, for
example). We then arrive at the following weak formulation for (2.5a): Find
$(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)\in\mathbf{H}^{1}_{0}(\Omega)\times\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega)\times\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega)\times\mathrm{L}^{2}_{0}(\Omega)\times\mathrm{L}^{2}_{0}(\Omega)$
such that

$\displaystyle
2\mu\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{\gamma})-\int_{\Omega}\varphi\mathop{\mathrm{div}}\nolimits\boldsymbol{\gamma}$
$\displaystyle=\int_{\Omega}\boldsymbol{b}\cdot\boldsymbol{\gamma}\quad\forall\boldsymbol{\gamma}\in\mathbf{H}^{1}_{0}(\Omega),$
(2.6a)
$\displaystyle{\frac{1}{\kappa}}\int_{\Omega}\\!\boldsymbol{v}\cdot\boldsymbol{\zeta}+{\sqrt{\frac{\nu}{\kappa}}}\int_{\Omega}\\!\\!\mathop{\mathbf{curl}}\nolimits\boldsymbol{\omega}\cdot\boldsymbol{\zeta}+{\frac{\nu}{\kappa}}\int_{\Omega}\\!\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}-\int_{\Omega}p\,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}$
$\displaystyle=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{\zeta}\quad\forall\boldsymbol{\zeta}\in\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega),$
(2.6b)
$\displaystyle-\int_{\Omega}\boldsymbol{\omega}\cdot\boldsymbol{\theta}+{\sqrt{\frac{\nu}{\kappa}}}\int_{\Omega}\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\cdot\boldsymbol{v}$
$\displaystyle=0\quad\forall\boldsymbol{\theta}\in\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega),$
(2.6c)
$\displaystyle-\frac{1}{\lambda}\int_{\Omega}\varphi\,\psi+\frac{\alpha}{\lambda}\int_{\Omega}p\psi-\int_{\Omega}\psi\,\mathop{\mathrm{div}}\nolimits\boldsymbol{u}$
$\displaystyle=0\quad\forall\psi\in\mathrm{L}_{0}^{2}(\Omega),$ (2.6d)
$\displaystyle-(c_{0}+\frac{\alpha^{2}}{\lambda})\int_{\Omega}p\,q+\frac{\alpha}{\lambda}\int_{\Omega}\varphi\,q-\int_{\Omega}q\,\mathop{\mathrm{div}}\nolimits\boldsymbol{v}$
$\displaystyle=\int_{\Omega}g\,q\quad\forall q\in\mathrm{L}_{0}^{2}(\Omega),$
(2.6e)

where we have also used the boundary conditions (2.5f).

Let us now define the following bilinear forms and linear functionals

$\displaystyle
a_{1}(\boldsymbol{u},\boldsymbol{\gamma}):=2\mu\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{\gamma}),\qquad
b_{1}(\boldsymbol{\gamma},\psi):=-\int_{\Omega}\psi\mathop{\mathrm{div}}\nolimits\boldsymbol{\gamma},\qquad\hat{b}_{1}(\boldsymbol{\zeta},q):=-\int_{\Omega}q\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta},$
$\displaystyle
a_{2}(\boldsymbol{v},\boldsymbol{\zeta}):={\frac{1}{\kappa}}\int_{\Omega}\boldsymbol{v}\cdot\boldsymbol{\zeta}+{\frac{\nu}{\kappa}}\int_{\Omega}\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta},\qquad
b_{2}(\boldsymbol{\theta},\boldsymbol{\zeta}):={\sqrt{\frac{\nu}{\kappa}}}\int_{\Omega}\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\cdot\boldsymbol{\zeta},$
$\displaystyle
a_{3}(\boldsymbol{\omega},\boldsymbol{\theta}):=\int_{\Omega}\boldsymbol{\omega}\cdot\boldsymbol{\theta},\qquad
a_{4}(\varphi,\psi):=\frac{1}{\lambda}\int_{\Omega}\varphi\,\psi,\qquad
b_{3}(q,\psi):=\frac{\alpha}{\lambda}\int_{\Omega}q\,\psi,\qquad
a_{5}(p,q):=\biggl{(}c_{0}+\frac{\alpha^{2}}{\lambda}\biggr{)}\int_{\Omega}p\,q,$
$\displaystyle
B(\boldsymbol{\gamma}):=\int_{\Omega}\boldsymbol{b}\cdot\boldsymbol{\gamma},\qquad
F(\boldsymbol{\zeta}):=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{\zeta},\qquad
G(q):=\int_{\Omega}g\,q,$

with which (2.6) is rewritten as follows: find
$(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)\in\mathbf{H}^{1}_{0}(\Omega)\times\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega)\times\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega)\times\mathrm{L}_{0}^{2}(\Omega)\times\mathrm{L}_{0}^{2}(\Omega)$
such that

$\displaystyle a_{1}(\boldsymbol{u},\boldsymbol{\gamma})$
$\displaystyle+\;b_{1}(\boldsymbol{\gamma},\varphi)$ $\displaystyle=$
$\displaystyle\;B(\boldsymbol{\gamma})$
$\displaystyle\qquad\forall\boldsymbol{\gamma}\in\mathbf{H}^{1}_{0}(\Omega),$
(2.7a) $\displaystyle a_{2}(\boldsymbol{v},\boldsymbol{\zeta})$
$\displaystyle+\;b_{2}(\boldsymbol{\omega},\boldsymbol{\zeta})$
$\displaystyle+\;\hat{b}_{1}(\boldsymbol{\zeta},p)$ $\displaystyle=$
$\displaystyle\;F(\boldsymbol{\zeta})$
$\displaystyle\qquad\forall\boldsymbol{\zeta}\in\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega),$
(2.7b) $\displaystyle b_{2}(\boldsymbol{\theta},\boldsymbol{v})$
$\displaystyle-\;a_{3}(\boldsymbol{\omega},\boldsymbol{\theta})$
$\displaystyle=$ $\displaystyle\;0$
$\displaystyle\qquad\forall\boldsymbol{\theta}\in\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega),$
(2.7c) $\displaystyle b_{1}(\boldsymbol{u},\psi)$
$\displaystyle-\;a_{4}(\varphi,\psi)$ $\displaystyle+\;b_{3}(p,\psi)$
$\displaystyle=$ $\displaystyle\;0$
$\displaystyle\qquad\forall\psi\in\mathrm{L}_{0}^{2}(\Omega),$ (2.7d)
$\displaystyle\hat{b}_{1}(\boldsymbol{v},q)$
$\displaystyle+\;b_{3}(q,\varphi)$ $\displaystyle-\;a_{5}(p,q)$
$\displaystyle=$ $\displaystyle\;G(q)$ $\displaystyle\qquad\forall
q\in\mathrm{L}_{0}^{2}(\Omega).$ (2.7e)

##  3 Solvability analysis

###  3.1 Preliminaries

The well-posedness analysis for (2.7) will be put in the framework of the
abstract Banach–Nečas–Babuška theory, which we state next (see, e.g., 30).

######  Theorem 3.1.

Let $(E_{1},\|\cdot\|_{E_{1}})$ be a reflexive Banach space,
$(E_{2},\|\cdot\|_{E_{2}})$ a Banach space, and $T:E_{1}\rightarrow
E_{2}^{\prime}$ a bounded, linear form satisfying the followings conditions:

  1. (BNB1) For each $y\in E_{2}\setminus\\{0\\}$, there exists $x\in E_{1}$ such that

$\langle T(x),y\rangle_{E_{2}^{\prime},E_{2}}\neq 0.$ (3.1)

  2. (BNB2) There exists $c>0$ such that

$\|T(x)\|_{E_{2}^{\prime}}\geq c\|x\|_{E_{1}}\text{ for all }x\in E_{1}.$
(3.2)

Then, for every $x^{*}\in E_{2}^{\prime}$ there exists a unique $x\in E_{1}$
such that

$T(x)=x^{*}.$

Let us first consider the product space

$\displaystyle\mathbf{X}$
$\displaystyle:=\mathbf{H}^{1}_{0}(\Omega)\times\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega)\times\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega)\times\mathrm{L}_{0}^{2}(\Omega)\times\mathrm{L}_{0}^{2}(\Omega),$

and, using the notation
$\vec{\boldsymbol{x}}:=(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)\in\mathbf{X}$,
we proceed to equip this space with the norm

$\displaystyle\|\vec{\boldsymbol{x}}\|_{\mathbf{X}}^{2}$
$\displaystyle:=2\mu\|\boldsymbol{\varepsilon}(\boldsymbol{u})\|^{2}_{0,\Omega}+{\frac{1}{\kappa}}\|\boldsymbol{v}\|^{2}_{0,\Omega}+{\frac{\nu}{\kappa}}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|^{2}_{0,\Omega}+\|\boldsymbol{\omega}\|^{2}_{0,\Omega}+\nu\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\omega}\|^{2}_{0,\Omega}+{c_{0}\|p\|^{2}_{0,\Omega}+\frac{1}{\lambda}\|\varphi+\alpha
p\|^{2}_{0,\Omega}.}$ (3.3)

######  Remark 3.2.

Note that the scaling of the vorticity seminorm with the viscosity (fifth term
in the right-hand side of (3.3) is not readily apparent from the weak
formulation (2.7), but it is suggested by the stability analysis in Lemma 3.3,
below.

Let us also introduce the bilinear form

$\displaystyle\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}\rangle$
$\displaystyle:=a_{1}(\boldsymbol{u},\boldsymbol{\gamma})+b_{1}(\boldsymbol{\gamma},\varphi)+a_{2}(\boldsymbol{v},\boldsymbol{\zeta})+b_{2}(\boldsymbol{\omega},\boldsymbol{\zeta})+\hat{b}_{1}(\boldsymbol{\zeta},p)+b_{2}(\boldsymbol{\theta},\boldsymbol{v})-a_{3}(\boldsymbol{\omega},\boldsymbol{\theta})$
$\displaystyle\qquad+b_{1}(\boldsymbol{u},\psi)-a_{4}(\varphi,\psi)+b_{3}(p,\psi)+\hat{b}_{1}(\boldsymbol{v},q)+b_{3}(q,\varphi)-a_{5}(p,q),$
(3.4)

induced by the operator
$\mathcal{A}_{\epsilon}:\mathbf{X}\rightarrow\mathbf{X}^{\prime}$ (where the
subscript $\epsilon$ indicates dependence with respect to the model parameters
$\kappa,\alpha,\mu,\nu,c_{0},\lambda$), and again we emphasise that we have
different scalings than those used in 5, 35.

From the Cauchy–Schwarz inequality we readily have the following bounds for
the bilinear forms in (2.7)

$\displaystyle a_{1}(\boldsymbol{u},\boldsymbol{\gamma})\leq
2\mu|\boldsymbol{u}|_{1,\Omega}|\boldsymbol{\gamma}|_{1,\Omega},\qquad
b_{1}(\boldsymbol{\gamma},\varphi)\leq\|\varphi\|_{0,\Omega}|\boldsymbol{\gamma}|_{1,\Omega},\qquad\hat{b}_{1}(\boldsymbol{\zeta},q)\leq\|q\|_{0,\Omega}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega},$
$\displaystyle
a_{2}(\boldsymbol{v},\boldsymbol{\zeta})\leq\dfrac{1}{\kappa}\|\boldsymbol{v}\|_{0,\Omega}\|\boldsymbol{\zeta}\|_{0,\Omega}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega},\qquad
b_{2}(\boldsymbol{\theta},\boldsymbol{\zeta})\leq\sqrt{\frac{\nu}{\kappa}}\|\mathop{\mathrm{curl}}\nolimits\boldsymbol{\theta}\|_{0,\Omega}\|\boldsymbol{\zeta}\|_{0,\Omega},$
$\displaystyle
a_{3}(\boldsymbol{\omega},\boldsymbol{\theta})\leq\|\boldsymbol{\omega}\|_{0,\Omega}\|\boldsymbol{\theta}\|_{0,\Omega},\qquad
a_{4}(\varphi,\psi)\leq\dfrac{1}{\lambda}\|\varphi\|_{0,\Omega}\|\psi\|_{0,\Omega},$
$\displaystyle
b_{3}(q,\psi)\leq\frac{\alpha}{\lambda}\|q\|_{0,\Omega}\|\psi\|_{0,\Omega},\qquad
a_{5}(p,q)\leq\left(c_{0}+\dfrac{\alpha^{2}}{\lambda}\right)\|p\|_{0,\Omega}\|q\|_{0,\Omega},$
(3.5)

for all $\boldsymbol{u},\boldsymbol{\gamma}\in\mathbf{H}^{1}(\Omega)$,
$\boldsymbol{v},\boldsymbol{\zeta}\in\mathbf{H}(\mathop{\mathrm{div}}\nolimits,\Omega)$,
$\boldsymbol{\omega},\boldsymbol{\theta}\in\mathbf{H}(\mathop{\mathbf{curl}}\nolimits,\Omega)$,
$\varphi,\psi,p,q\in\mathrm{L}^{2}(\Omega)$.

######  Lemma 3.3.

Consider the bilinear form defined in (3.4). For all
$\vec{\boldsymbol{x}}\in\mathbf{X}$, there exists
$\vec{\boldsymbol{y}}\in\mathbf{X}$ such that

$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}\rangle\gtrsim\|\vec{\boldsymbol{x}}\|^{2}_{\mathbf{X}}\quad\text{
and
}\quad\|\vec{\boldsymbol{y}}\|_{\mathbf{X}}\leq\sqrt{2}\|\vec{\boldsymbol{x}}\|_{\mathbf{X}}.$
(3.6)

######  Proof 3.4.

Using first and second Young’s inequalities it is not difficult to assert that

$(\boldsymbol{\zeta},\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta})_{0,\Omega}\leq{\dfrac{1}{2}}\|\boldsymbol{\zeta}\|_{0,\Omega}^{2}+{\dfrac{1}{2}}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\|_{0,\Omega}^{2}\quad\text{
and
}\quad(q,\psi)_{0,\Omega}\leq\dfrac{\varepsilon}{2}\|q\|_{0,\Omega}^{2}+\dfrac{1}{2\varepsilon}\|\psi\|_{0,\Omega}^{2}.$
(3.7)

Next, for a given
$\vec{\boldsymbol{z}}:=(\boldsymbol{\gamma},\boldsymbol{\zeta},\boldsymbol{\theta},\psi,q)$
we can construct
$\vec{\boldsymbol{y}}:=\bigl{(}\boldsymbol{\gamma},\boldsymbol{\zeta}+\frac{1}{2}\sqrt{\kappa\nu}\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta},-\boldsymbol{\theta},-\psi,-q\bigr{)}$.
We then invoke (3.7), so that we can ensure that

$\displaystyle\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{z}}),\vec{\boldsymbol{y}}\rangle$
$\displaystyle=2\mu|\boldsymbol{\gamma}|_{1,\Omega}^{2}+\dfrac{1}{\kappa}\|\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\dfrac{\sqrt{\nu}}{2\sqrt{\kappa}}(\boldsymbol{\zeta},\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta})_{0,\Omega}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\dfrac{\nu}{2}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\|_{0,\Omega}^{2}$
$\displaystyle\qquad+\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\dfrac{1}{\lambda}\|\psi\|_{0,\Omega}^{2}-\dfrac{\alpha}{\lambda}(q,\psi)_{0,\Omega}-\dfrac{\alpha}{\lambda}(q,\psi)+\left(c_{0}+\dfrac{\alpha^{2}}{\lambda}\right)\|q\|_{0,\Omega}^{2}$
$\displaystyle\geq
2\mu|\boldsymbol{\gamma}|_{1,\Omega}^{2}+\dfrac{1}{\kappa}\|\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\dfrac{\sqrt{\nu}}{2\sqrt{\kappa}}(\boldsymbol{\zeta},\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta})_{0,\Omega}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\dfrac{\nu}{2}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\|_{0,\Omega}^{2}+\|\boldsymbol{\theta}\|_{0,\Omega}^{2}$
$\displaystyle\qquad+\left(c_{0}+\dfrac{\alpha^{2}}{\lambda}-\dfrac{\varepsilon\alpha}{2\lambda}\right)\|q\|_{0,\Omega}^{2}+\left(\dfrac{1}{\lambda}-\dfrac{\alpha}{2\lambda\varepsilon}\right)\|\psi\|_{0,\Omega}^{2}.$

Now, taking $\varepsilon:=\dfrac{5\alpha}{4}+\dfrac{\lambda c_{0}}{\alpha}$ we
can deduce that

$\displaystyle
c_{0}+\dfrac{\alpha^{2}}{\lambda}-\dfrac{\varepsilon\alpha}{2\lambda}$
$\displaystyle\geq\dfrac{3}{8}\left(c_{0}+\dfrac{\alpha^{2}}{\lambda}\right),$
(3.8a)
$\displaystyle\dfrac{1}{\lambda}-\dfrac{\alpha}{2\lambda\varepsilon}=\dfrac{\dfrac{1}{4}+\dfrac{\lambda
c_{0}}{\alpha^{2}}}{\dfrac{5}{4}+\dfrac{\lambda
c_{0}}{\alpha^{2}}}\dfrac{1}{\lambda}$
$\displaystyle\geq\dfrac{1}{5}\dfrac{1}{\lambda},$ (3.8b)

and using (3.8a) and (3.8b) we readily obtain the bound

$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{z}}),\vec{\boldsymbol{y}}\rangle\geq\dfrac{1}{10}\|\vec{\boldsymbol{z}}\|_{\mathbf{X}}^{2}.$

Finally, the definition of the preliminary $\mathbf{X}$-norm (3.3) and
triangle inequality yield the estimates

$\displaystyle\|\vec{\boldsymbol{y}}\|_{\mathbf{X}}^{2}$
$\displaystyle=2\mu|\boldsymbol{\gamma}|_{1,\Omega}^{2}+\frac{1}{\kappa}\left\|\boldsymbol{\zeta}+\dfrac{\sqrt{\kappa\nu}}{2}\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\right\|_{0,\Omega}^{2}+{\dfrac{\nu}{\kappa}}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\nu\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\|_{0,\Omega}^{2}+\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\left(c_{0}+\dfrac{\alpha^{2}}{\lambda}\right)\|q\|_{0,\Omega}^{2}+\dfrac{1}{\lambda}\|\psi\|_{0,\Omega}^{2}$
$\displaystyle\leq
2\mu|\boldsymbol{\gamma}|_{1,\Omega}^{2}+\dfrac{2\nu}{\kappa}\|\boldsymbol{\zeta}\|_{0,\Omega}^{2}+{\dfrac{\nu}{\kappa}}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\dfrac{3\nu}{2}\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\theta}\|_{0,\Omega}^{2}+\|\boldsymbol{\theta}\|_{0,\Omega}^{2}+\left(c_{0}+\dfrac{\alpha^{2}}{\lambda}\right)\|q\|_{0,\Omega}^{2}+\dfrac{1}{\lambda}\|\psi\|_{0,\Omega}^{2}$
$\displaystyle\leq 2\|\vec{\boldsymbol{z}}\|_{\mathbf{X}}^{2},$

which completes the proof.

Using the Banach–Nečas–Babuška result, from (3.1) and (3.6) we immediately
conclude that problem (2.7) has a unique solution
$(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)$ in the space
$\mathbf{X}$. However we can note that the bound on
$(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)$ in the
$\mathbf{X}$-norm will degenerate with $\varepsilon$ (tending either to zero
or to infinity).

As a preliminary result required in the sequel, we recall the following
classical inf-sup condition for the bilinear form $\hat{b}_{1}$ (which
coincides with that of the Stokes problem, see, e.g., 33).

######  Lemma 3.5.

There exists $\beta_{1}>0$ such that

$\sup_{\boldsymbol{\gamma}\in\mathbf{H}^{1}_{0}(\Omega)\setminus\\{\boldsymbol{0}\\}}\dfrac{\hat{b}_{1}(\psi,\boldsymbol{\gamma})}{|\boldsymbol{\gamma}|_{1,\Omega}}=\sup_{\boldsymbol{\gamma}\in\mathbf{H}^{1}_{0}(\Omega)\setminus\\{\boldsymbol{0}\\}}\dfrac{-(\psi,\mathop{\mathrm{div}}\nolimits\boldsymbol{\gamma})_{0,\Omega}}{|\boldsymbol{\gamma}|_{1,\Omega}}\geq\beta_{1}\|\psi\|_{0,\Omega}\qquad\text{
for all }\psi\in\mathrm{L}_{0}^{2}(\Omega).$

###  3.2 Parameter-robust well-posedness

Before addressing the well-posedness of (2.7) robustly with respect to
$\varepsilon$, we first note that the bilinear forms defining the solution
operator suggest to modify the metric in $\mathbf{X}$ and include the
following particular parameter-weighting of the functional spaces

$\displaystyle\mathbf{X}_{\epsilon}$
$\displaystyle:={\sqrt{2\mu}}\mathbf{H}^{1}_{0}(\Omega)\times\biggl{[}{\sqrt{\frac{1}{\kappa}}}\mathbf{L}^{2}(\Omega)\cap\sqrt{\frac{\nu}{\kappa}}\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega)\biggr{]}\times\biggl{[}\mathbf{L}^{2}(\Omega)\cap\sqrt{\nu}\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega)\biggr{]}$
$\displaystyle\qquad\times\biggl{[}{\sqrt{\frac{1}{\lambda}}}\mathrm{L}^{2}(\Omega)\cap{\sqrt{\frac{1}{2\mu}}}\mathrm{L}_{0}^{2}(\Omega)\biggr{]}\times\biggl{[}{\sqrt{c_{0}+\frac{\alpha^{2}}{\lambda}}}\mathrm{L}^{2}(\Omega)\cap\left(\sqrt{\kappa}\mathrm{H}^{1}_{0}(\Omega)+\sqrt{\frac{\kappa}{\nu}}\mathrm{L}_{0}^{2}(\Omega)\right)\biggr{]}.$

An important observation is that $\mathbf{X}_{\epsilon}$ contains the same
vectors $\vec{\boldsymbol{x}}$ that are in $\mathbf{X}$, but which are bounded
in the norm $\|\cdot\|_{\epsilon}$ to be defined later.

Note also that, proceeding similarly as in 13 (see also 46), we have
decomposed the space for fluid pressure as the sum

$\mathrm{L}_{0}^{2}(\Omega)={\sqrt{\dfrac{\kappa}{\nu}}}\mathrm{L}_{0}^{2}(\Omega)+{\sqrt{\kappa}}{\mathrm{H}^{1}_{0}(\Omega)}\cap\mathrm{L}_{0}^{2}(\Omega),$

and have endowed it with the norm $\|\cdot\|_{r}$ defined, thanks to (2.2), as
follows

$\|q\|^{2}_{r}:=\inf_{\displaystyle{s\in\sqrt{\kappa}\mathrm{H}^{1}_{0}(\Omega)\cap\mathrm{L}_{0}^{2}(\Omega)}}\left\\{\left(\dfrac{\kappa}{\nu}+c_{0}\right)\|q-s\|_{0,\Omega}^{2}+c_{0}\|s\|_{0,\Omega}^{2}+\|\sqrt{\kappa}\nabla
s\|_{0,\Omega}^{2}\right\\}.$ (3.9)

With this norm we see that, for example, the boundedness of the bilinear form
$\hat{b}_{1}(\cdot,\cdot)$ can be written as follows

$\displaystyle(q,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta})_{0,\Omega}$
$\displaystyle=(q-s+s,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta})_{0,\Omega}=-(\nabla
s,\boldsymbol{\zeta})_{0,\Omega}+(q-s,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta})_{0,\Omega}$
$\displaystyle\leq\left\\{\kappa|s|_{1,\Omega}^{2}+\dfrac{\kappa}{\nu}\|q-s\|_{0,\Omega}^{2}\right\\}^{1/2}\left\\{\dfrac{1}{\kappa}\|\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega}^{2}\right\\}^{1/2},$

for all $s\in\mathrm{H}_{0}^{1}(\Omega)$, so from (3.9) we have

$(q,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta})_{0,\Omega}\leq\|q\|_{r}\left\\{\dfrac{1}{\kappa}\|\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega}^{2}\right\\}^{1/2}.$

On the other hand, we have the following robust-in-$\epsilon$ inf-sup
condition for $\hat{b}_{1}(\cdot,\cdot)$.

######  Lemma 3.6.

There exists $\beta_{0}>0$ independent of the parameters in $\epsilon$, such
that

$\sup_{\boldsymbol{\zeta}\in\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega)\setminus\\{\boldsymbol{0}\\}}\dfrac{-(q,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta})}{\left\\{\dfrac{1}{\kappa}\|\boldsymbol{\zeta}\|_{0,\Omega}^{2}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}\|_{0,\Omega}^{2}\right\\}^{1/2}}\geq\beta_{0}\|q\|_{r}\qquad\text{
for all }q\in\mathrm{L}_{0}^{2}(\Omega).$ (3.10)

######  Proof 3.7.

Owing to 33, we know that there exists $\beta_{1}>0$ (independent of the
parameters) such that

$\|\nabla q\|_{-1,\Omega}\geq\beta_{1}\|q\|_{0,\Omega}\text{ for all
}q\in\mathrm{L}_{0}^{2}(\Omega).$

Then, for the operator
$\nabla^{-1}:\nabla\mathrm{L}_{0}^{2}(\Omega)\rightarrow\mathrm{L}_{0}^{2}(\Omega)$
(where $\nabla\mathrm{L}_{0}^{2}(\Omega)$ is a closed subspace of
$\mathbf{H}^{-1}(\Omega)$), we can deduce that

$\|\nabla^{-1}\|_{\mathcal{L}(\nabla\mathrm{L}_{0}^{2}(\Omega),\mathrm{L}_{0}^{2}(\Omega))}\leq\beta_{1}^{-1}.$

Using the Poincaré inequality we can find a positive constant $c:=c(\Omega)$
such that

$\|q\|_{0,\Omega}\leq c|q|_{1,\Omega}\text{ for all
}q\in{\mathrm{H}^{1}(\Omega)}\cap\mathrm{L}_{0}^{2}(\Omega),$

or, equivalently,

$\|\nabla^{-1}\|_{\displaystyle\mathcal{L}(\nabla(\mathrm{H}^{1}(\Omega)\cap\mathrm{L}_{0}^{2}(\Omega)),\mathrm{H}^{1}(\Omega)\cap\mathrm{L}_{0}^{2}(\Omega))}\leq
c.$ (3.11)

Then, we have that

$\|\nabla
q\|_{\mathbf{L}^{2}(\Omega)+\nu^{-1/2}\mathbf{H}^{-1}(\Omega)}\geq\max\\{c,\beta_{1}^{-1}\\}\displaystyle\inf_{{\begin{subarray}{c}q=q_{1}+q_{2}\\\
q_{1}\in\mathrm{H}_{0}^{1}(\Omega)\cap\mathrm{L}_{0}^{2}(\Omega),\\\
q_{2}\in\mathrm{L}_{0}^{2}(\Omega)\end{subarray}}}\left\\{|q_{1}|_{1,\Omega}^{2}+\dfrac{1}{\nu}\|q_{2}\|_{0,\Omega}^{2}\right\\}^{1/2}.$
(3.12)

Multiplying (3.12) by $\sqrt{\kappa}$ and applying algebraic manipulations, we
can conclude that the inf-sup condition (3.10) holds with
$\beta_{0}:=\max\\{c,\beta_{1}^{-1}\\}$.

As done in the previous subsection, the unique solvability analysis will also
follow from the Banach–Nečas–Babuška theory, but now using the space
$\mathbf{X}_{\epsilon}$ endowed with the new norm

$\|\vec{\boldsymbol{x}}\|_{\epsilon}^{2}:=\|\vec{\boldsymbol{x}}\|_{\mathbf{X}}^{2}+\|p\|_{r}^{2}+\dfrac{1}{2\mu}\|\varphi\|_{0,\Omega}^{2}.$
(3.13)

Problem (2.7) is written as: Find
$\vec{\boldsymbol{x}}\in\mathbf{X}_{\epsilon}$ such that

$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}\rangle=\mathcal{F}(\vec{\boldsymbol{y}})\qquad\text{for
all }\vec{\boldsymbol{y}}\in\mathbf{X}_{\epsilon},$

or in operator form as follows

$\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}})=\mathcal{F}\qquad\text{in
$\mathbf{X}^{\prime}_{\epsilon}$},$ (3.14)

and the norm of the solution operator is defined as

$\|\mathcal{A}_{\epsilon}\|_{\mathcal{L}(\mathbf{X}_{\epsilon},\mathbf{X}_{\epsilon}^{\prime})}:=\sup_{\vec{\boldsymbol{x}},\vec{\boldsymbol{y}}\in\mathbf{X}_{\epsilon}\setminus\\{\vec{\boldsymbol{0}}\\}}\frac{|\langle\mathcal{A}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}\rangle|}{\|\vec{\boldsymbol{x}}\|_{\epsilon}\|\vec{\boldsymbol{y}}\|_{\epsilon}}.$

Translating Theorem 3.1 to the present context, we aim to prove that the
operator $\mathcal{A}_{\epsilon}$ is continuous, that is

$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}\rangle\lesssim\|\vec{\boldsymbol{x}}\|_{\epsilon}\|\vec{\boldsymbol{y}}\|_{\epsilon}\qquad\forall\vec{\boldsymbol{x}},\vec{\boldsymbol{y}}\in\mathbf{X}_{\epsilon},$
(3.15)

and that the following global inf-sup condition is satisfied

$\sup_{\vec{\boldsymbol{y}}\in\mathbf{X}_{\epsilon}\setminus\\{\vec{\boldsymbol{0}}\\}}\frac{\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}\rangle}{\|\vec{\boldsymbol{y}}\|_{\epsilon}}\gtrsim\|\vec{\boldsymbol{x}}\|_{\epsilon}\qquad\forall\vec{\boldsymbol{x}}\in\mathbf{X}_{\epsilon}.$
(3.16)

######  Theorem 3.8.

Let $\|\cdot\|_{\epsilon}$ be defined as in (3.13). Then the bilinear form
induced by $\mathcal{A}_{\epsilon}$ (cf. (3.4)) is continuous and inf-sup
stable under the norm $\|\cdot\|_{\epsilon}$, i.e., the conditions (3.15) and
(3.16) are satisfied.

######  Proof 3.9.

For the continuity of $\mathcal{A}_{\epsilon}$, it suffices to use (3.1), the
norm definition (3.13), Cauchy–Schwarz inequality, the definition of
$\mathcal{A}_{\epsilon}$, and (3.10), to arrive at

$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}\rangle\leq
2\|\vec{\boldsymbol{x}}\|_{\epsilon}\|\vec{\boldsymbol{y}}\|_{\epsilon}.$

For the global inf-sup, we take a given
$\vec{\boldsymbol{x}}\in\mathbf{X}_{\epsilon}$, and for lemma 3.3, there
exists $\vec{\boldsymbol{y}}_{1}\in\mathbf{X}$ such that

$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}_{1}\rangle\geq\dfrac{1}{10}\|\vec{\boldsymbol{x}}\|^{2}_{\mathbf{X}}\quad\text{and}\quad\|\vec{\boldsymbol{y}}_{1}\|_{\mathbf{X}}\leq\sqrt{2}\|\vec{\boldsymbol{x}}\|_{\mathbf{X}}.$

Using Lemma 3.6 and Lemma 3.5 we can find $\boldsymbol{\zeta}_{2}$,
$\boldsymbol{\gamma}_{3}$ and constants $C_{1}$, $C_{2}$, $\hat{C}_{1}$,
$\hat{C}_{2}$ such that

$-(p,\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}_{2})\geq
C_{1}\|p\|_{r}^{2}\quad\text{and}\quad\left\\{\dfrac{1}{\kappa}\|\boldsymbol{\zeta}_{2}\|_{0,\Omega}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}_{2}\|_{0,\Omega}^{2}\right\\}^{1/2}\leq
C_{2}\|p\|_{r}.$

and

$-(\varphi,\mathop{\mathrm{div}}\nolimits\boldsymbol{\gamma})\geq\hat{C}_{1}\frac{1}{2\mu}\|\varphi\|_{0,\Omega}^{2}\quad\text{and}\quad\sqrt{2\mu}|\boldsymbol{\gamma}_{3}|_{1,\Omega}\leq\hat{C}_{2}\frac{1}{\sqrt{2\mu}}\|\varphi\|_{0,\Omega}.$

Taking
$\vec{\boldsymbol{y}}_{2}:=(\boldsymbol{0},\boldsymbol{\zeta}_{2},\boldsymbol{0},0,0)$,
$\vec{\boldsymbol{y}}_{3}:=(\boldsymbol{\gamma}_{3},\boldsymbol{0},\boldsymbol{0},0,0)$,
$\delta>0$, and $\hat{\delta}>0$, we have that

$\displaystyle\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),10\vec{\boldsymbol{y}}_{1}+\delta\vec{\boldsymbol{y}}_{2}+\hat{\delta}\vec{\boldsymbol{y}}_{3}\rangle$
$\displaystyle\geq\|\vec{\boldsymbol{x}}\|^{2}_{\mathbf{X}}+\delta(a_{2}(\boldsymbol{v},\boldsymbol{\zeta}_{2})+b_{2}(\boldsymbol{\omega},\boldsymbol{\zeta}_{2})+\hat{b}_{1}(\boldsymbol{\zeta}_{2},p))+\hat{\delta}(a_{1}(\boldsymbol{u},\boldsymbol{\gamma}_{3})+b_{1}(\boldsymbol{\gamma}_{3},\varphi))$
$\displaystyle\geq\|\vec{\boldsymbol{x}}\|^{2}_{\mathbf{X}}-\delta\left\\{\dfrac{1}{\kappa}\|\boldsymbol{v}\|_{0,\Omega}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}\right\\}^{1/2}\left\\{\dfrac{1}{\kappa}\|\boldsymbol{\zeta}_{2}\|_{0,\Omega}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}_{2}\|_{0,\Omega}^{2}\right\\}^{1/2}$
$\displaystyle\qquad-\delta\dfrac{\nu}{\sqrt{\kappa}}\|\mathop{\mathrm{curl}}\nolimits\boldsymbol{\omega}\|_{0,\Omega}\|\boldsymbol{\zeta}_{2}\|_{0,\Omega}+\delta
C_{1}\|p\|_{r}^{2}-2\mu\hat{\delta}|\boldsymbol{u}|_{1,\Omega}|\boldsymbol{\gamma}_{3}|_{1,\Omega}+\hat{\delta}\frac{\hat{C}_{1}}{2\mu}\|\varphi\|_{0,\Omega}^{2}$
$\displaystyle\geq\|\vec{\boldsymbol{x}}\|^{2}_{\mathbf{X}}-\dfrac{1}{2}\left(\dfrac{1}{\kappa}\|\boldsymbol{v}\|_{0,\Omega}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}\right)-\dfrac{\nu}{2}\|\mathop{\mathrm{curl}}\nolimits\boldsymbol{\omega}\|_{0,\Omega}^{2}-\mu|\boldsymbol{u}|_{1,\Omega}^{2}+\delta
C_{1}\|p\|_{r}^{2}$
$\displaystyle\qquad+\hat{\delta}\frac{\hat{C}_{1}}{2\mu}\|\varphi\|_{0,\Omega}^{2}-\delta^{2}\dfrac{1}{\kappa}\|\boldsymbol{\zeta}_{2}\|_{0,\Omega}-\delta^{2}\dfrac{\nu}{2\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}_{2}\|_{0,\Omega}^{2}-2\mu\hat{\delta}^{2}|\boldsymbol{\gamma}_{3}|_{1,\Omega}^{2}$
$\displaystyle\geq\dfrac{1}{2}\|\vec{\boldsymbol{x}}\|^{2}_{\mathbf{X}}+\delta(C_{1}-\delta
C_{2}^{2})\|p\|_{r}^{2}+\hat{\delta}\frac{1}{2\mu}\left(\hat{C}_{1}-\hat{\delta}\hat{C_{2}}^{2}\right)\|\varphi\|_{0,\Omega}^{2}.$

Then, choosing $\delta:=\dfrac{C_{1}}{2C_{2}^{2}}$ and
$\hat{\delta}:=\dfrac{\hat{C}_{1}}{2\hat{C}_{2}^{2}}$, we can deduce the
estimates

$\displaystyle\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),10\vec{\boldsymbol{y}}_{1}+\delta\vec{\boldsymbol{y}}_{2}+\hat{\delta}\vec{\boldsymbol{y}}_{3}\rangle$
$\displaystyle\geq\dfrac{1}{2}\|\vec{\boldsymbol{x}}\|^{2}_{\mathbf{X}}+\dfrac{C_{1}^{2}}{4C_{2}^{2}}\|p\|_{r}^{2}+\dfrac{\hat{C}_{1}^{2}}{4\hat{C}_{2}^{2}}\frac{1}{2\mu}\|\varphi\|_{0,\Omega}^{2}\geq\dfrac{1}{2}\min\left\\{1,\dfrac{C_{1}^{2}}{2C_{2}^{2}},\dfrac{\hat{C}_{1}^{2}}{2\hat{C}_{2}^{2}}\right\\}\|\vec{\boldsymbol{x}}\|^{2}_{\epsilon},\qquad\text{and}$
$\displaystyle\|10\vec{\boldsymbol{y}}_{1}+\delta\vec{\boldsymbol{y}}_{2}+\hat{\delta}\vec{\boldsymbol{y}}_{3}\|_{\epsilon}$
$\displaystyle\leq 10\sqrt{2}\|\vec{\boldsymbol{x}}\|_{\epsilon}+\delta
C_{2}\|p\|_{r}+\hat{\delta}\hat{C}_{2}\frac{1}{\sqrt{2\mu}}\|\varphi\|_{0,\Omega}\leq\max\left\\{10\sqrt{2},\dfrac{C_{1}}{2C_{2}},\dfrac{\hat{C}_{1}}{2\hat{C}_{2}}\right\\}\|\vec{\boldsymbol{x}}\|_{\epsilon}.$

And from these relations we can conclude that:

$\sup_{\vec{\boldsymbol{y}}\in\mathbf{X}_{\epsilon}\setminus\\{\vec{\boldsymbol{0}}\\}}\frac{\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}\rangle}{\|\vec{\boldsymbol{y}}\|_{\epsilon}}\geq\frac{\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),10\vec{\boldsymbol{y}}_{1}+\delta\vec{\boldsymbol{y}}_{2}+\hat{\delta}\vec{\boldsymbol{y}}_{3}\rangle}{\|10\vec{\boldsymbol{y}}_{1}+\delta\vec{\boldsymbol{y}}_{2}+\hat{\delta}\vec{\boldsymbol{y}}_{3}\|_{\epsilon}}\geq\dfrac{1}{2}\dfrac{\min\left\\{1,\dfrac{C_{1}^{2}}{2C_{2}^{2}},\dfrac{\hat{C}_{1}^{2}}{2\hat{C}_{2}^{2}}\right\\}\|\vec{\boldsymbol{x}}\|^{2}_{\epsilon}}{\max\left\\{10\sqrt{2},\dfrac{C_{1}}{2C_{2}},\dfrac{\hat{C}_{1}}{2\hat{C}_{2}}\right\\}\|\vec{\boldsymbol{x}}\|_{\epsilon}}\gtrsim\|\vec{\boldsymbol{x}}\|_{\epsilon}.$

###  3.3 Operator preconditioning

We recall from, e.g., 39, 46, that since $\mathcal{A}_{\epsilon}$ maps
$\mathbf{X}_{\epsilon}$ to its dual, when solving the discrete version of
(3.14), iterative methods are not directly applicable unless a modified
problem is considered, for example

$\mathcal{B}\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}})=\mathcal{B}\mathcal{F}\qquad\text{in
$\mathbf{X}_{\epsilon}$,}$

where $\mathcal{B}:\mathbf{X}_{\epsilon}^{\prime}\to\mathbf{X}_{\epsilon}$ is
an appropriately defined isomorphism. As usual, one can take $\mathcal{B}$ as
the Riesz map (self-adjoint and positive definite) whose inverse defines a
scalar product $(\cdot,\cdot)_{\mathbf{X}_{\epsilon}}$ on
$\mathbf{X}_{\epsilon}$, and the operator $\mathcal{B}\mathcal{A}_{\epsilon}$
is also self-adjoint with respect to this inner product. Then

$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\boldsymbol{y}\rangle=(\mathcal{B}\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}})_{\mathbf{X}_{\epsilon}}\quad\text{and}\quad\|\mathcal{B}\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}})\|_{\mathbf{X}_{\epsilon}}=\|\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}})\|_{\mathbf{X}_{\epsilon}^{\prime}},$

and therefore, using the definition of the operator norms, it is readily
deduced that

$\|\mathcal{B}\mathcal{A}_{\epsilon}\|_{\mathcal{L}(\mathbf{X}_{\epsilon},\mathbf{X}_{\epsilon})}=\|\mathcal{A}_{\epsilon}\|_{\mathcal{L}(\mathbf{X}_{\epsilon},\mathbf{X}_{\epsilon}^{\prime})}\quad\text{and}\quad\|(\mathcal{B}\mathcal{A}_{\epsilon})^{-1}\|_{\mathcal{L}(\mathbf{X}_{\epsilon},\mathbf{X}_{\epsilon})}=\|\mathcal{A}_{\epsilon}^{-1}\|_{\mathcal{L}(\mathbf{X}_{\epsilon},\mathbf{X}_{\epsilon}^{\prime})}.$

Then, if an appropriate metric is chosen such that the norms of
$\mathcal{A}_{\epsilon}$ and of $\mathcal{A}_{\epsilon}^{-1}$ are bounded by
constants independent of the model parameters $\epsilon$, then the condition
number of the preconditioned system will also be independent of the model
parameters.

Proceeding then similarly as in, e.g., 43, 13, 16, we consider the following
block-diagonal preconditioners focusing on the case of mixed boundary
conditions (and therefore not including contributions related to the zero-mean
value of fluid and total pressure)

$\displaystyle\mathcal{B}_{1}$
$\displaystyle=\begin{pmatrix}{(-\mathop{\mathbf{div}}\nolimits(2\mu\boldsymbol{\varepsilon}))}^{-1}&0&0&0&0\\\
0&\\!\\!\bigl{(}{\kappa^{-1}\mathbf{I}-\frac{\nu}{\kappa}\nabla}\mathop{\mathrm{div}}\nolimits\bigr{)}^{-1}&0&0&0\\\
0&0&\\!\\!(\mathbf{I}+{\nu}\mathop{\mathbf{curl}}\nolimits)^{-1}&0&0\\\
0&0&0&\\!\\!\bigl{(}\bigl{(}\frac{1}{\lambda}+\frac{1}{2\mu}\bigr{)}I\bigr{)}^{-1}&0\\\
0&0&0&0&\\!\\!\bigl{(}\bigl{(}c_{0}+\frac{\alpha^{2}}{\lambda}+{{\kappa}}\bigr{)}I\bigr{)}^{-1}\end{pmatrix},$
(3.17a) $\displaystyle\mathcal{B}_{2}$
$\displaystyle=\begin{pmatrix}{(-\mathop{\mathbf{div}}\nolimits(2\mu\boldsymbol{\varepsilon}))}^{-1}&0&0&0&0\\\
0&\\!\\!\bigl{(}{\kappa^{-1}}\mathbf{I}-\frac{\nu}{\kappa}\nabla\mathop{\mathrm{div}}\nolimits\bigr{)}^{-1}&0&0&0\\\
0&0&\\!\\!(\mathbf{I}+{\nu}\mathop{\mathbf{curl}}\nolimits)^{-1}&0&0\\\
0&0&0&\\!\\!\bigl{(}\bigl{(}\frac{1}{\lambda}+\frac{1}{2\mu}\bigr{)}I\bigr{)}^{-1}&0\\\
0&0&0&0&\\!\\!\bigl{(}(c_{0}+\frac{\alpha^{2}}{\lambda})I-{\kappa}\Delta\bigr{)}^{-1}\end{pmatrix},$
(3.17b) $\displaystyle\mathcal{B}_{3}$
$\displaystyle=\begin{pmatrix}{(-\mathop{\mathbf{div}}\nolimits(2\mu\boldsymbol{\varepsilon}))}^{-1}&0&0&0&0\\\
0&\bigl{(}{\kappa^{-1}}\mathbf{I}-{\bigl{(}1+\frac{\nu}{\kappa}\bigr{)}\nabla}\mathop{\mathrm{div}}\nolimits\bigr{)}^{-1}&0&0&0\\\
0&0&(\mathbf{I}+{\nu}\mathop{\mathbf{curl}}\nolimits)^{-1}&0&0\\\
0&0&0&\lx@intercol\hfil\hbox{\multirowsetup{$\widehat{\mathcal{B}}$}}\hfil\lx@intercol\\\
0&0&0&\end{pmatrix},$ (3.17c)

with

${\widehat{\mathcal{B}}=\begin{pmatrix}\biggl{(}\displaystyle{\frac{{\displaystyle{1}}}{{\displaystyle{\lambda}}}}+\displaystyle{\frac{{\displaystyle{1}}}{{\displaystyle{2\mu}}}}\biggr{)}I&\displaystyle{\frac{{\displaystyle{\alpha}}}{{\displaystyle{\lambda}}}}I\\\
\displaystyle{\frac{{\displaystyle{\alpha}}}{{\displaystyle{\lambda}}}}I&\biggl{(}1+c_{0}+\displaystyle{\frac{{\displaystyle{\alpha^{2}}}}{{\displaystyle{\lambda}}}}\biggr{)}I\end{pmatrix}^{-1}+\begin{pmatrix}\biggl{(}\displaystyle{\frac{{\displaystyle{1}}}{{\displaystyle{\lambda}}}}+\displaystyle{\frac{{\displaystyle{1}}}{{\displaystyle{2\mu}}}}\biggr{)}I&\displaystyle{\frac{{\displaystyle{\alpha}}}{{\displaystyle{\lambda}}}}I\\\
\displaystyle{\frac{{\displaystyle{\alpha}}}{{\displaystyle{\lambda}}}}I&\biggl{(}c_{0}+\displaystyle{\frac{{\displaystyle{\alpha^{2}}}}{{\displaystyle{\lambda}}}}\biggr{)}I-{\kappa}\Delta\end{pmatrix}^{-1}}.$

Note that only $\mathcal{B}_{3}$ results from the Riesz map corresponding to
$\mathbf{X}_{\epsilon}$ with the complete norm as in (3.13), while
$\mathcal{B}_{1},\mathcal{B}_{2}$ are approximations of $\mathcal{B}_{3}$. In
particular, $\mathcal{B}_{1}$ simply considers the parameter weighting
suggested by the weak formulation (2.6) (but taking into account the scaling
of the vorticity seminorm according to Remark 3.2) combined with the Riesz map
associated with the natural regularity of that formulation, and
$\mathcal{B}_{2}$ includes also the sum of spaces leading to the pressure
Laplacian forms which are key in achieving robustness for Darcy-type problems
13. The full form $\mathcal{B}_{3}$ also includes the non-standard
Brezzi–Braess type of block $\widehat{\mathcal{B}}$ for total and fluid
pressures, which is needed in perturbed saddle-point problems with penalty as
proposed in 16 Section 4 (see also 17).

##  4 Analysis of a finite element method

Let $\mathcal{T}_{h}$ denote a family of tetrahedral meshes (triangular in 2D)
on $\Omega$ and denote by $\mathcal{E}_{h}$ the set of all facets (edges in
2D) in the mesh. By $h_{K}$ we denote the diameter of the element $K$ and by
$h_{F}$ we denote the length/area of the facet $F$. As usual, by $h$ we denote
the maximum of the diameters of elements in $\mathcal{T}_{h}$. For all meshes
we assume that they are sufficiently regular (there exists a uniform positive
constant $\eta_{1}$ such that each element $K$ is star-shaped with respect to
a ball of radius greater than $\eta_{1}h_{K}$). It is also assumed that there
exists $\eta_{2}>0$ such that for each element and every facet $F\in\partial
K$, we have that $h_{F}\geq\eta_{2}h_{K}$, see, e.g., 53, 30. By
$\mathrm{P}_{k}(\Theta)$ we will denote be the space of polynomials of total
degree at most $k$ defined locally on the domain $\Theta$, and denote by
$\mathbf{P}_{k}(\Theta)$ and $\mathbb{P}_{k}(\Theta)$ their vector- and
tensor-valued counterparts, respectively. By $\mathcal{E}_{h}$ we will denote
the set of all facets and will distinguish between facets lying on the
interior and the two sub-boundaries
$\mathcal{E}_{h}=\mathcal{E}_{h}^{\mathrm{int}}\cup\mathcal{E}_{h}^{\Gamma}$.

For smooth scalar fields $w$ defined on $\mathcal{T}_{h}$, the symbol
$w^{\pm}$ denotes the traces of $\boldsymbol{w}$ on $e$ that are the
extensions from the interior of the two elements $K^{+}$ and $K^{-}$ sharing
the facet $e$. The symbols
$\left\\{\\!\\!\left\\{\cdot\right\\}\\!\\!\right\\}$ and
$\left[\\!\left[\cdot\right]\\!\right]$ denote, respectively, the average and
jump operators defined as
$\left\\{\\!\\!\left\\{w\right\\}\\!\\!\right\\}:=\frac{1}{2}(w^{-}+w^{+})$,
$\left[\\!\left[w\right]\\!\right]:=(w^{-}-w^{+})$. The element-wise action of
a differential operator is denoted with a subindex $h$, for example,
$\nabla_{h}$ will denote the broken gradient operator.

The discrete spaces that we consider herein correspond, for $k\geq 0$, to the
generalised Taylor–Hood element pair ($\mathbf{P}_{k+2}-\mathrm{P}_{k+1}$) for
the displacement / total pressure approximation, the H(div)-conforming
Raviart–Thomas elements of degree $k$ (denoted $\mathbf{RT}_{k}$) for velocity
approximation, the H(curl)-conforming Nédélec elements of the first kind and
order $k+1$ (denoted $\mathbf{ND}_{k+1}$) for filtration vorticity (see, e.g.,
18, 31 for precise definitions of these families of spaces), and piecewise
polynomials of degree $k$ for the approximation of interstitial pressure

$\displaystyle\mathbf{U}_{h}$
$\displaystyle:=\\{\boldsymbol{u}_{h}\in\mathbf{H}^{1}_{0}(\Omega):\boldsymbol{u}_{h}|_{K}\in\mathbf{P}_{k+2}(K),\quad\forall
K\in\mathcal{T}_{h}\\},$ $\displaystyle\mathbf{V}_{h}$
$\displaystyle:=\\{\boldsymbol{v}_{h}\in\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega):\boldsymbol{v}_{h}|_{K}\in\mathbf{RT}_{k}(K),\quad\forall
K\in\mathcal{T}_{h}\\},$ $\displaystyle\mathbf{W}_{h}$
$\displaystyle:=\\{\boldsymbol{\omega}_{h}\in\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega):\boldsymbol{\omega}_{h}|_{K}\in\mathbf{ND}_{k+1}(K),\quad\forall
K\in\mathcal{T}_{h}\\},$ $\displaystyle{\mathbf{M}^{m}_{h}}$
$\displaystyle:={\\{\boldsymbol{u}_{h}\in\mathbf{L}^{2}(\Omega):\boldsymbol{u}_{h}|_{K}\in\mathbf{P}_{m}(K),\quad\forall
K\in\mathcal{T}_{h}\\}},$ (4.1) $\displaystyle\mathrm{U}_{h}$
$\displaystyle:=\\{u_{h}\in\mathrm{H}^{1}_{0}(\Omega):u_{h}|_{K}\in\mathrm{P}_{k}(K),\quad\forall
K\in\mathcal{T}_{h}\\},$ $\displaystyle\mathrm{Z}_{h}$
$\displaystyle:=\\{\psi_{h}\in{C^{0}(\Omega)}:\psi_{h}|_{K}\in{\mathrm{P}_{k+1}(K)},\quad\forall
K\in\mathcal{T}_{h}\\},$ $\displaystyle\mathrm{Q}_{h}$
$\displaystyle:=\\{q_{h}\in\mathrm{L}^{2}(\Omega):q_{h}|_{K}\in\mathrm{P}_{k}(K),\quad\forall
K\in\mathcal{T}_{h}\\}.$

Note that other combinations of finite element families are feasible as well
(as long as appropriate discrete inf-sup conditions are satisfied). In 2D we
will consider the same space $\mathrm{Q}_{h}$ for both total and fluid
pressures and we will use the $\mathbf{P}_{2}-\mathrm{P}_{0}$ pair for
displacement and total pressure approximation.

The (conforming) finite element scheme associated with (2.7) reads: Find

$(\boldsymbol{u}_{h},\boldsymbol{v}_{h},\boldsymbol{\omega}_{h},\varphi_{h},p_{h})\in\mathbf{X}_{h}:=\mathbf{U}_{h}\times\mathbf{V}_{h}\times\mathbf{W}_{h}\times\mathrm{Z}_{h}\times\mathrm{Q}_{h},$

such that

$\displaystyle a_{1}(\boldsymbol{u}_{h},\boldsymbol{\gamma}_{h})$
$\displaystyle+\;b_{1}(\boldsymbol{\gamma}_{h},\varphi_{h})$ $\displaystyle=$
$\displaystyle\;B(\boldsymbol{\gamma}_{h})$
$\displaystyle\qquad\forall\boldsymbol{\gamma}_{h}\in\mathbf{U}_{h},$ (4.2a)
$\displaystyle a_{2}(\boldsymbol{v}_{h},\boldsymbol{\zeta}_{h})$
$\displaystyle+\;b_{2}(\boldsymbol{\omega}_{h},\boldsymbol{\zeta}_{h})$
$\displaystyle+\;\hat{b}_{1}(\boldsymbol{\zeta}_{h},p_{h})$ $\displaystyle=$
$\displaystyle\;F(\boldsymbol{\zeta}_{h})$
$\displaystyle\qquad\forall\boldsymbol{\zeta}_{h}\in\mathbf{V}_{h},$ (4.2b)
$\displaystyle b_{2}(\boldsymbol{\theta}_{h},\boldsymbol{v}_{h})$
$\displaystyle-\;a_{3}(\boldsymbol{\omega}_{h},\boldsymbol{\theta}_{h})$
$\displaystyle=$ $\displaystyle\;0$
$\displaystyle\qquad\forall\boldsymbol{\theta}_{h}\in\mathbf{W}_{h},$ (4.2c)
$\displaystyle b_{1}(\boldsymbol{u}_{h},\psi_{h})$
$\displaystyle-\;a_{4}(\varphi_{h},\psi_{h})$
$\displaystyle+\;b_{3}(p_{h},\psi_{h})$ $\displaystyle=$ $\displaystyle\;0$
$\displaystyle\qquad\forall\psi_{h}\in\mathrm{Z}_{h},$ (4.2d)
$\displaystyle\hat{b}_{1}(\boldsymbol{v}_{h},q_{h})$
$\displaystyle+\;b_{3}(q_{h},\varphi_{h})$
$\displaystyle-\;a_{5}(p_{h},q_{h})$ $\displaystyle=$
$\displaystyle\;G(q_{h})$ $\displaystyle\qquad\forall q_{h}\in\mathrm{Q}_{h}.$
(4.2e)

Similarly as in the continuous case, we define

$\displaystyle\mathbf{X}_{\epsilon,h}$
$\displaystyle:={2\mu}\mathbf{U}_{h}\times\left[\sqrt{\frac{1}{\kappa}}\mathbf{M}_{h}^{k+1}\cap\sqrt{\dfrac{\nu}{\kappa}}\mathbf{V}_{h}\right]\times[\mathbf{M}_{h}^{k+1}\cap\sqrt{\nu}\mathbf{W}_{h}]\times\left[{\dfrac{1}{\lambda}}\mathrm{Q}_{h}\cap{\dfrac{1}{2\mu}}\mathrm{Z}_{h}\right]\times\left[\sqrt{c_{0}}\mathrm{Q}_{h}\cap\left(\sqrt{\dfrac{\kappa}{\nu}}\mathrm{Q}_{h}+\sqrt{\kappa}\mathrm{U}_{h}\right)\right],$

and the associated discrete norm is

$\displaystyle\|\vec{\boldsymbol{x}}_{h}\|_{\varepsilon,h}^{2}$
$\displaystyle:=2\mu\|\boldsymbol{\varepsilon}(\boldsymbol{u}_{h})\|^{2}_{0,\Omega}+\frac{1}{\kappa}\|\boldsymbol{v}_{h}\|^{2}_{0,\Omega}+\frac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}_{h}\|^{2}_{0,\Omega}+\|\boldsymbol{\omega}_{h}\|^{2}_{0,\Omega}+\nu\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\omega}_{h}\|^{2}_{0,\Omega}+\frac{1}{2\mu}\|\varphi_{h}\|_{0,\Omega}^{2}$
$\displaystyle\quad+\inf_{s_{h}\in\mathrm{U}_{h}}\biggl{[}\left(c_{0}+\frac{\kappa}{\nu}\right)\|p_{h}-s_{h}\|_{0,\Omega}^{2}+c_{0}\|s_{h}\|^{2}+\bigl{\|}\sqrt{\kappa}\nabla_{h}s_{h}\bigr{\|}^{2}_{0,\Omega}\biggr{]}+\frac{1}{\lambda}\|\varphi_{h}+\alpha
p_{h}\|^{2}_{0,\Omega}+c_{0}\|p_{h}\|^{2}_{0,\Omega}.$ (4.3)

As in the continuous case, now problem (4.2) is written as: Find
$\vec{\boldsymbol{x}}_{h}\in\mathbf{X}_{\epsilon,h}$ such that

$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}_{h}),\vec{\boldsymbol{y}}_{h}\rangle=\mathcal{F}_{h}(\vec{\boldsymbol{y}}_{h})\qquad\text{for
all }\vec{\boldsymbol{y}}_{h}\in\mathbf{X}_{\epsilon,h},$

or in operator form as follows

$\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}_{h})=\mathcal{F}_{h}\qquad\text{in
$\mathbf{X}^{\prime}_{\epsilon,h}$}.$ (4.4)

######  Lemma 4.1.

There exists $\beta_{1}>0$ independent of the parameters in $\epsilon$ and
$h$, such that

$\sup_{\boldsymbol{\zeta}_{h}\in\mathbf{V}_{h}\setminus\\{\boldsymbol{0}\\}}\dfrac{-(q_{h},\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}_{h})_{0,\Omega}}{\left\\{\dfrac{1}{\kappa}\|\boldsymbol{\zeta}_{h}\|_{0,\Omega}^{2}+\dfrac{\nu}{\kappa}\|\mathop{\mathrm{div}}\nolimits\boldsymbol{\zeta}_{h}\|_{0,\Omega}^{2}\right\\}^{1/2}}\geq\beta_{1}\|q_{h}\|_{r}\qquad\text{
for all }q_{h}\in\mathrm{Q}_{h}.$ (4.5)

######  Proof 4.2.

The proof requires to assume that the continuous inf-sup condition holds.
Then, similarly to the proof of that result (Lemma 3.6), the first part (steps
until (3.11)) is a consequence of the fact that the spaces $\mathbf{V}_{h}$
and $\mathrm{Q}_{h}$ satisfy the usual discrete inf-sup condition for the
Stokes problem. Then, it suffices to follow the scaling argument in (3.12),
also valid at the discrete level, to complete the desired condition.

Analogously to the continuous case, we need the following conditions to be
satisfied to guarantee existence and uniqueness to problem (4.4):

$\displaystyle\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}_{h}),\vec{\boldsymbol{y}}_{h}\rangle$
$\displaystyle\lesssim\|\vec{\boldsymbol{x}}_{h}\|_{\epsilon}\|\vec{\boldsymbol{y}}_{h}\|_{\epsilon,h}\qquad\forall\vec{\boldsymbol{x}}_{h},\vec{\boldsymbol{y}}_{h}\in\mathbf{X}_{\epsilon,h},$
(4.6a)
$\displaystyle\sup_{\vec{\boldsymbol{y}}_{h}\in\mathbf{X}_{\epsilon,h}\setminus\\{\vec{\boldsymbol{0}}\\}}\frac{\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}_{h}),\vec{\boldsymbol{y}}_{h}\rangle}{\|\vec{\boldsymbol{y}}_{h}\|_{\epsilon,h}}$
$\displaystyle\gtrsim\|\vec{\boldsymbol{x}}_{h}\|_{\epsilon,h}\qquad\forall\vec{\boldsymbol{x}}_{h}\in\mathbf{X}_{\epsilon,h}.$
(4.6b)

######  Theorem 4.3.

Let $\|\cdot\|_{\epsilon,h}$ be defined as in (4.3). Then the bilinear form
induced by $\mathcal{A}_{\epsilon}$ (cf. (3.4)) is continuous and inf-sup
stable under the norm $\|\cdot\|_{\epsilon,h}$, i.e., the conditions (4.6a)
and (4.6b) are satisfied.

######  Proof 4.4.

We use again that the pairs $(\mathbf{U}_{h},\mathrm{Z}_{h})$ and
$(\mathbf{V}_{h},\mathrm{Q}_{h})$ are inf-sup stable spaces for the usual
bilinear form in Stokes problem, which ensures that we can find discrete
versions $\vec{\boldsymbol{y}}_{h,2}$, $\vec{\boldsymbol{y}}_{h,3}$ of the
tuples $\vec{\boldsymbol{y}}_{2}$ and $\vec{\boldsymbol{y}}_{3}$,
respectively, constructed in Theorem 3.8. We also note that
$\mathop{\mathbf{curl}}\nolimits\mathbf{W}_{h}\subset\mathbf{V}_{h}$, and
therefore we can prove the discrete version of Lemma 3.3. Then the desired
result is a consequence of repeating the arguments used in the proof of
Theorem 3.8.

######  Theorem 4.5.

For given $\boldsymbol{b},\boldsymbol{f}\in\mathbf{L}^{2}(\Omega)$ and $g\in
L^{2}(\Omega)$, problem (4.2) has a unique solution
$(\boldsymbol{u}_{h},\boldsymbol{v}_{h},\boldsymbol{\omega}_{h},\varphi_{h},q_{h})\in\mathbf{X}_{\epsilon,h}$.
In addition, the solution satisfies the following continuous dependence on
data

$\|(\boldsymbol{u}_{h},\boldsymbol{v}_{h},\boldsymbol{\omega}_{h},\varphi_{h},p_{h})\|_{\epsilon}\lesssim(\|\boldsymbol{b}\|_{0,\Omega}+\|\boldsymbol{f}\|_{0,\Omega}+\|g\|_{0,\Omega}),$

and the following Céa estimate

$\displaystyle\|(\boldsymbol{u}-\boldsymbol{u}_{h},\boldsymbol{v}-\boldsymbol{v}_{h},\boldsymbol{\omega}-\boldsymbol{\omega}_{h},\varphi-\varphi_{h},p-p_{h})\|_{\epsilon}\leq\bigl{(}1+\alpha^{-1}\|\mathcal{A}_{\epsilon}\|_{\mathcal{L}(\mathbf{X}_{\epsilon},\mathbf{X}^{\prime}_{\epsilon})}\bigr{)}\|(\boldsymbol{u}-\boldsymbol{\gamma}_{h},\boldsymbol{v}-\boldsymbol{\zeta}_{h},\boldsymbol{\omega}-\boldsymbol{\theta}_{h},\varphi-\psi_{h},p-q_{h})\|_{\epsilon},$

for all
$(\boldsymbol{\gamma}_{h},\boldsymbol{\zeta}_{h},\boldsymbol{\theta}_{h},\psi_{h},q_{h})\in\mathbf{X}_{\epsilon,h}$,
where $\alpha$ is the positive constant associated with (4.6b).

######  Proof 4.6.

The existence and uniqueness of the solution is obtained in a similar way to
its counterpart in the continuous level. For the corresponding Céa estimate,
we proceed to denote as
$\vec{\boldsymbol{x}}:=(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)$,
$\vec{\boldsymbol{x_{h}}}:=(\boldsymbol{u}_{h},\boldsymbol{v}_{h},\boldsymbol{\omega}_{h},\varphi_{h},p_{h})$
and
$\vec{\boldsymbol{y}_{h}}:=(\boldsymbol{\gamma}_{h},\boldsymbol{\zeta}_{h},\boldsymbol{\theta}_{h},\psi_{h},q_{h})$.
From (4.6b), we can infer that there exists a positive constant $\alpha$
independent of the parameters such that

$\sup_{\vec{\boldsymbol{z}}_{h}\in\mathbf{X}_{\epsilon,h}\setminus\\{\vec{\boldsymbol{0}}\\}}\frac{\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}_{h}),\vec{\boldsymbol{z}}_{h}\rangle}{\|\vec{\boldsymbol{z}}_{h}\|_{\epsilon}}\geq\alpha\|\vec{\boldsymbol{x}}_{h}\|_{\epsilon}\qquad\forall\vec{\boldsymbol{x}}_{h}\in\mathbf{X}_{\epsilon,h}.$

Using the error equation, we readily obtain that
$\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}),\vec{\boldsymbol{y}}_{h}\rangle=\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{x}}_{h}),\vec{\boldsymbol{y}}_{h}\rangle$.
Furthermore, since
$\vec{\boldsymbol{y}}_{h}-\vec{\boldsymbol{x}}_{h}\in\mathbf{X}_{\epsilon,h}$
we can deduce that

$\displaystyle\|\vec{\boldsymbol{x}}-\vec{\boldsymbol{x}}_{h}\|_{\epsilon}$
$\displaystyle\leq\|\vec{\boldsymbol{x}}-\vec{\boldsymbol{y}}_{h}\|_{\epsilon}+\|\vec{\boldsymbol{y}}_{h}-\vec{\boldsymbol{x}}_{h}\|_{\epsilon}$
$\displaystyle\leq\|\vec{\boldsymbol{x}}-\vec{\boldsymbol{y}}_{h}\|_{\epsilon}+\alpha^{-1}\sup_{\vec{\boldsymbol{z}}_{h}\in\mathbf{X}_{\epsilon,h}\setminus\\{\vec{\boldsymbol{0}}\\}}\frac{\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{y}}_{h}-\vec{\boldsymbol{x}}_{h}),\vec{\boldsymbol{z}}_{h}\rangle}{\|\vec{\boldsymbol{z}}_{h}\|_{\epsilon}}$
$\displaystyle\leq\|\vec{\boldsymbol{x}}-\vec{\boldsymbol{y}}_{h}\|_{\epsilon}+\alpha^{-1}\sup_{\vec{\boldsymbol{z}}_{h}\in\mathbf{X}_{\epsilon,h}\setminus\\{\vec{\boldsymbol{0}}\\}}\frac{\langle\mathcal{A}_{\epsilon}(\vec{\boldsymbol{y}}_{h}-\vec{\boldsymbol{x}}),\vec{\boldsymbol{z}}_{h}\rangle}{\|\vec{\boldsymbol{z}}_{h}\|_{\epsilon}}$
$\displaystyle\leq\|\vec{\boldsymbol{x}}-\vec{\boldsymbol{y}}_{h}\|_{\epsilon}+\alpha^{-1}\sup_{\vec{\boldsymbol{z}}_{h}\in\mathbf{X}_{\epsilon,h}\setminus\\{\vec{\boldsymbol{0}}\\}}\frac{\|\mathcal{A}_{\epsilon}\|_{\mathcal{L}(\mathbf{X}_{\epsilon},\mathbf{X}^{\prime}_{\epsilon})}\|\vec{\boldsymbol{y}}_{h}-\vec{\boldsymbol{x}}\|_{\epsilon}\|\vec{\boldsymbol{z}}_{h}\|_{\epsilon}}{\|\vec{\boldsymbol{z}}_{h}\|_{\epsilon}}$
$\displaystyle\leq(1+\alpha^{-1}\|\mathcal{A}_{\epsilon}\|_{\mathcal{L}(\mathbf{X}_{\epsilon},\mathbf{X}^{\prime}_{\epsilon})})\|\vec{\boldsymbol{x}}-\vec{\boldsymbol{y}}_{h}\|_{\epsilon},$

which finishes the proof.

Let us recall, from, e.g., 30 Chapters 16, 17 and 22, the following
approximation properties of the finite element subspaces (4.1), which are
obtained using the classical interpolation theory. Assume that
$(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)\in\mathbf{H}^{1+s}(\Omega)\times{\mathbf{H}^{s}(\mathop{\mathrm{div}}\nolimits,\Omega)\times\mathbf{H}^{s}(\mathop{\mathbf{curl}}\nolimits,\Omega)}\times{H}^{s}(\Omega)\times{H}^{s}(\Omega)$,
for some $s\in(1/2,k+1]$. Then there exists $C>0$, independent of $h$, such
that

$\displaystyle\|\boldsymbol{\varepsilon}(\boldsymbol{u}-\mathcal{I}_{h}(\boldsymbol{u}))\|_{0,\Omega}$
$\displaystyle\leq Ch^{s}|\boldsymbol{u}|_{1+s,\Omega},$ (4.7a)
$\displaystyle\|p-\Pi_{h}(p)\|_{0,\Omega}$ $\displaystyle\leq
Ch^{s}|p|_{s,\Omega},$ (4.7b)
$\displaystyle\|\boldsymbol{v}-\mathcal{I}_{h}^{RT}(\boldsymbol{v})\|_{0,\Omega}$
$\displaystyle\leq Ch^{s}|\boldsymbol{v}|_{s,\Omega},$ (4.7c)
$\displaystyle\|\boldsymbol{\omega}-\mathcal{I}_{h}^{N}(\boldsymbol{\omega})\|_{0,\Omega}$
$\displaystyle\leq Ch^{s}|\boldsymbol{\omega}|_{s,\Omega},$ (4.7d)

where
$\Pi_{h}:\mathrm{L}_{0}^{2}\rightarrow\mathrm{Z}_{h}\subset\mathrm{Q}_{h}$ is
the $L^{2}$-projection and
$\mathcal{I}_{h}:\mathbf{H}^{1}(\Omega)\rightarrow\mathbf{U}_{h}$,
$\mathcal{I}_{h}^{RT}:\mathbf{H}_{0}(\mathop{\mathrm{div}}\nolimits,\Omega)\rightarrow\mathbf{V}_{h}$,
$\mathcal{I}_{h}^{N}:\mathbf{H}_{0}(\mathop{\mathbf{curl}}\nolimits,\Omega)\rightarrow\mathbf{W}_{h}$
are the Lagrange, Raviart–Thomas and Nédelec interpolators respectively.

######  Theorem 4.7.

Let
$(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)\in\mathbf{X}_{\epsilon}$
and
$(\boldsymbol{u}_{h},\boldsymbol{v}_{h},\boldsymbol{\omega}_{h},\varphi_{h},p_{h})\in\mathbf{X}_{\epsilon,h}$
be the unique solutions to the continuous and discrete problems (2.7) and
(4.2), respectively. Assume that
$(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)\in\mathbf{H}^{1+s}(\Omega)\times{\mathbf{H}^{s}(\mathop{\mathrm{div}}\nolimits,\Omega)\times\mathbf{H}^{s}(\mathop{\mathbf{curl}}\nolimits,\Omega)}\times{H}^{s}(\Omega)\times{H}^{s}(\Omega)$,
for some $s\in(1/2,{k+1}]$. Then, there exists $C>0$, independent of the mesh
size $h$ and of the model parameters $\epsilon$, such that

$\displaystyle\|(\boldsymbol{u}-\boldsymbol{u}_{h},\boldsymbol{v}-\boldsymbol{v}_{h},\boldsymbol{\omega}-\boldsymbol{\omega}_{h},\varphi-\varphi_{h},p-p_{h})\|_{\epsilon}\leq
Ch^{s}\|(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)\|_{s,\epsilon},$

where

$\displaystyle\|(\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p)\|_{s,\epsilon}^{2}:=$
$\displaystyle
2\mu|\boldsymbol{u}|_{1+s,\Omega}^{2}+\frac{1}{\kappa}|\boldsymbol{v}|_{s,\Omega}^{2}+\frac{\nu}{\kappa}|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}|_{s,\Omega}^{2}+|\boldsymbol{\omega}|_{s,\Omega}^{2}+\nu|\mathop{\mathbf{curl}}\nolimits\boldsymbol{\omega}|_{s,\Omega}^{2}+\frac{1}{2\mu}|\varphi|_{s,\Omega}^{2}+\left(c_{0}+\frac{\kappa}{\nu}\right)|p|_{s,\Omega}^{2}+\frac{1}{\lambda}|\varphi+\alpha
p|_{s,\Omega}^{2}.$

######  Proof 4.8.

This result follows immediately after choosing the tuple

$(\boldsymbol{\gamma}_{h},\boldsymbol{\zeta}_{h},\boldsymbol{\theta}_{h},\psi_{h},q_{h}):=(\mathcal{I}_{h}(\boldsymbol{u}),\mathcal{I}_{h}^{RT}(\boldsymbol{v}),\mathcal{I}_{h}^{N}(\boldsymbol{\omega}),\Pi_{h}(\varphi),\Pi_{h}(p)),$

in Theorem 4.5, and then using the estimates (4.7a)-(4.7d) together with the
following commuting properties of the Raviart–Thomas and Nédélec operators

$\mathop{\mathrm{div}}\nolimits\mathcal{I}_{h}^{RT}(\boldsymbol{v})=\Pi_{h}(\mathop{\mathrm{div}}\nolimits\boldsymbol{v}),\qquad\mathop{\mathbf{curl}}\nolimits\mathcal{I}_{h}^{N}(\boldsymbol{\omega})=\Pi_{h}(\mathop{\mathbf{curl}}\nolimits\boldsymbol{\omega}).$

Finally, it suffices to invoke the fact that
$\|q\|_{r}\leq\left(\frac{\kappa}{\nu}+c_{0}\right)\|q\|_{0,\Omega}$.

######  Remark 4.9.

The error estimate above holds true also in the limit cases of near
incompressibility ($\lambda\to\infty$) and near impermeability ($\kappa\to
0$). This robustness in these cases (along with variations in other
parameters) is explored in Section 5.2. We also stress that in the non-viscous
regime ($\nu=0$) we do not have a velocity Laplacian in the fluid momentum
equation (2.3b), and no vorticity is required. In that case the formulation in
(2.7) recovers the 4-field Biot formulation from 16. Optimal convergence for
the Biot limit is confirmed numerically in Section 5.1, below.

##  5 Numerical verification

The aim of this section is to experimentally validate the theoretical results
presented in Section 4. We certify by error convergence verification the
proposed finite element method, and then we apply the formulation in the
simulation of a representative viscous flow in a poroelastic channel. We also
evaluate the robustness of the preconditioners in (3.17). The numerical
implementation uses the open-source finite element framework Gridap 11, and is
available in the public domain 21. For the solution of the linear systems in
the accuracy verification tests we employ the sparse direct method MUMPS,
while in the preconditioner tests we use the preconditioner MINRES iterative
solver; see the corresponding section below for the full details underlying
the set up of the preconditioner tests.

Table 5.1: Accuracy test in 2D. Error history (errors for each field variable
on a sequence of successively refined grids) using non-weighted spaces,
convergence rates, and $\ell^{\infty}$-norm of the loss of mass (the residual
of the mass conservation equation at the discrete level, eqn. (5.1)).

DoF | $h$ | $e_{1}(\boldsymbol{u})\\!$ | r | $e_{\mathop{\mathrm{div}}\nolimits}(\boldsymbol{v})\\!$ | r | $e_{\mathop{\mathbf{curl}}\nolimits}(\boldsymbol{\omega})$ | r | $e_{0}(\varphi)$ | r | $e_{0}(p)$ | r | $\mathrm{loss}_{h}$
---|---|---|---|---|---|---|---|---|---|---|---|---
Discretisation with $k=0$
93 | 0.7071 | 2.36e+00 | $\star$ | 1.98e+00 | $\star$ | 8.43e+00 | $\star$ | 3.10e+00 | $\star$ | 4.77e-01 | $\star$ | 2.66e-15
309 | 0.3536 | 1.04e+00 | 1.18 | 1.36e+00 | 0.53 | 6.74e+00 | 0.32 | 1.75e+00 | 0.83 | 2.59e-01 | 0.88 | 8.88e-15
1125 | 0.1768 | 4.68e-01 | 1.16 | 6.94e-01 | 0.98 | 3.47e+00 | 0.96 | 8.98e-01 | 0.96 | 9.13e-02 | 1.50 | 1.66e-14
4293 | 0.0884 | 2.27e-01 | 1.04 | 3.49e-01 | 0.99 | 1.74e+00 | 0.99 | 4.52e-01 | 0.99 | 3.87e-02 | 1.24 | 3.68e-14
16773 | 0.0442 | 1.13e-01 | 1.01 | 1.74e-01 | 1.00 | 8.73e-01 | 1.00 | 2.26e-01 | 1.00 | 1.84e-02 | 1.08 | 1.20e-13
66309 | 0.0221 | 5.65e-02 | 1.00 | 8.73e-02 | 1.00 | 4.37e-01 | 1.00 | 1.13e-01 | 1.00 | 9.05e-03 | 1.02 | 1.99e-13
Discretisation with $k=1$
221 | 0.7071 | 8.54e-01 | $\star$ | 1.27e+00 | $\star$ | 5.40e+00 | $\star$ | 1.26e+00 | $\star$ | 2.04e-01 | $\star$ | 8.95e-15
789 | 0.3536 | 1.90e-01 | 2.17 | 3.01e-01 | 2.07 | 1.47e+00 | 1.87 | 3.92e-01 | 1.69 | 3.72e-02 | 2.46 | 2.30e-14
2981 | 0.1768 | 4.80e-02 | 1.99 | 7.75e-02 | 1.96 | 3.86e-01 | 1.93 | 1.01e-01 | 1.96 | 7.69e-03 | 2.27 | 5.56e-14
11589 | 0.0884 | 1.20e-02 | 2.00 | 1.95e-02 | 1.99 | 9.84e-02 | 1.97 | 2.54e-02 | 1.99 | 1.82e-03 | 2.08 | 1.15e-13
45701 | 0.0442 | 3.00e-03 | 2.00 | 4.89e-03 | 2.00 | 2.48e-02 | 1.99 | 6.36e-03 | 2.00 | 4.50e-04 | 2.02 | 2.40e-13
181509 | 0.0221 | 7.50e-04 | 2.00 | 1.22e-03 | 2.00 | 6.23e-03 | 1.99 | 1.59e-03 | 2.00 | 1.12e-04 | 2.00 | 5.00e-13

The symbol $\star$ denotes that no experimental convergence rate has been
computed at the first refinement level.

###  5.1 Accuracy tests

Let us consider the unit square domain $\Omega=(0,1)^{2}$ together with the
manufactured solutions

$\displaystyle\boldsymbol{u}(x,y)=\begin{pmatrix}\sin(\pi[x+y])\\\
\cos(\pi[x^{2}+y^{2}])\end{pmatrix},\qquad p(x,y)=\sin(\pi x+y)\sin(\pi y),$
$\displaystyle\boldsymbol{v}(x,y)=\begin{pmatrix}\sin(\pi x)\sin(\pi y)\\\
\cos(\pi x)\cos(2\pi
y)\end{pmatrix},\qquad\omega(x,y)={\sqrt{\frac{\nu}{\kappa}}}\mathrm{curl}\,\boldsymbol{v},\qquad\varphi(x,y)=-\lambda\mathop{\mathrm{div}}\nolimits\boldsymbol{u}+\alpha
p.$

The model parameters assume the arbitrary values $\nu=1$, $\lambda=1$,
$\mu=1$, $\kappa=1$, $c_{0}=1$, $\alpha=1$; and the loading and source terms,
together with essential boundary conditions are computed from the manufactured
solutions above. The mean value of fluid and total pressure is prescribed to
coincide with the mean values of the manufactured pressures, which is
implemented by means of a real Lagrange multiplier. Sequences of successively
refined uniform tetrahedral meshes are used to compute approximate solutions
and to generate the error history (error decay with the mesh size $h$ and
experimental convergence rates, using norms in non-weighted spaces for each
individual unknown and at each refinement level). We display the error history
in Table 5.1, where the method confirms an asymptotic optimal convergence of
$O(h^{k+1})$ for each variable and for both polynomial degrees. We remark that
for this test we use continuous and piecewise polynomials of degree $k+2$ for
displacement and discontinuous piecewise polynomials of degree $k$ for total
pressure. Note that in 2D, the filtration vorticity is a scalar field
$\omega=\sqrt{\frac{\nu}{\kappa}}\mathrm{curl}\,\boldsymbol{v}$ and the
appropriate functional space is $\mathrm{H}^{1}_{0}(\Omega)$. In the discrete
setting we then select

$\mathrm{W}_{h}:=\\{\omega_{h}\in{\mathrm{H}^{1}_{0}(\Omega)}:\omega_{h}|_{K}\in\mathrm{P}_{1}(K),\
\forall K\in\mathcal{T}_{h}\\}.$

Moreover, to strengthen the numerical evidence, in the last column of the
table we report the ($L^{2}$-projection onto $\mathrm{Z}_{h}$ of the) loss of
mass

$\mathrm{loss}_{h}:=-(c_{0}+\frac{\alpha^{2}}{\lambda})p_{h}+\frac{\alpha}{\lambda}\varphi_{h}-\mathop{\mathrm{div}}\nolimits(\boldsymbol{v}_{h})-g,$
(5.1)

for each refinement level, confirming a satisfaction of the mass conservation
equation on the order of machine accuracy. This test has been performed using
pure Dirichlet boundary conditions.

In regards to Remark 4.9, we perform exactly the same error history as
reported in Table 5.1, now in the non-viscous regime (setting $\nu=0$). Of
course, then vorticity is zero and optimal convergence is attained for all
other fields also in this case, as shown in Table 5.2.

Table 5.2: Accuracy test in 2D. Error history for the non-viscous regime with
$\nu=0$ (errors for each field variable except vorticity on a sequence of
successively refined grids) using non-weighted spaces, convergence rates, and
$\ell^{\infty}$-norm of the loss of mass (the residual of the mass
conservation equation at the discrete level, eqn. (5.1)).

DoF | $h$ | $e_{1}(\boldsymbol{u})\\!$ | r | $e_{\mathop{\mathrm{div}}\nolimits}(\boldsymbol{v})\\!$ | r | $e_{0}(\varphi)$ | r | $e_{0}(p)$ | r | $\mathrm{loss}_{h}$
---|---|---|---|---|---|---|---|---|---|---
Discretisation with $k=0$
93 | 0.7071 | 2.36e+00 | $\star$ | 1.93e+00 | $\star$ | 3.08e+00 | $\star$ | 2.61e-01 | $\star$ | 2.61e-15
309 | 0.3536 | 1.04e+00 | 1.18 | 1.34e+00 | 0.52 | 1.74e+00 | 0.83 | 1.40e-01 | 0.90 | 1.04e-14
1’125 | 0.1768 | 4.68e-01 | 1.15 | 6.91e-01 | 0.96 | 8.97e-01 | 0.95 | 7.15e-02 | 0.97 | 8.17e-14
4’293 | 0.0884 | 2.27e-01 | 1.04 | 3.48e-01 | 0.99 | 4.52e-01 | 0.99 | 3.59e-02 | 0.99 | 3.37e-13
16’773 | 0.0442 | 1.13e-01 | 1.01 | 1.74e-01 | 1.00 | 2.26e-01 | 1.00 | 1.80e-02 | 1.00 | 3.64e-12
66’309 | 0.0221 | 5.65e-02 | 1.00 | 8.73e-02 | 1.00 | 1.13e-01 | 1.00 | 9.00e-03 | 1.00 | 2.43e-11
Discretisation with $k=1$
221 | 0.7071 | 8.54e-01 | $\star$ | 1.26e+00 | $\star$ | 1.25e+00 | $\star$ | 9.43e-02 | $\star$ | 7.75e-15
789 | 0.3536 | 1.90e-01 | 2.17 | 3.01e-01 | 2.06 | 3.90e-01 | 1.68 | 2.54e-02 | 1.89 | 2.37e-14
2’981 | 0.1768 | 4.81e-02 | 1.99 | 7.76e-02 | 1.95 | 1.01e-01 | 1.96 | 6.48e-03 | 1.97 | 5.83e-14
11’589 | 0.0884 | 1.20e-02 | 2.00 | 1.96e-02 | 1.99 | 2.53e-02 | 1.99 | 1.63e-03 | 1.99 | 1.25e-13
45’701 | 0.0442 | 3.00e-03 | 2.00 | 4.91e-03 | 2.00 | 6.35e-03 | 2.00 | 4.07e-04 | 2.00 | 4.19e-13
181’509 | 0.0221 | 7.51e-04 | 2.00 | 1.23e-03 | 2.00 | 1.59e-03 | 2.00 | 1.02e-04 | 2.00 | 7.12e-12

The symbol $\star$ denotes that no experimental convergence rate has been
computed at the first refinement level.

###  5.2 Robustness of the convergence rates

Table 5.3: Accuracy test in 3D. Error history (total error in the weighted
norm that leads to the definition of $\mathcal{B}_{3}$ on a sequence of
successively refined grids).

DoF | $h$ | $e((\boldsymbol{u},\boldsymbol{v},\boldsymbol{\omega},\varphi,p))$ | r
---|---|---|---
Discretisation with $k=0$
4’675 | 0.6124 | 1.02e+01 | $\star$
19’939 | 0.3062 | 3.87e+00 | 1.402
11’0083 | 0.1531 | 1.18e+00 | 1.718
718’339 | 0.0765 | 3.24e-01 | 1.862
5’162’755 | 0.0383 | 8.51e-02 | 1.928
Discretisation with $k=1$
12’766 | 0.6124 | 3.03e+00 | $\star$
55’514 | 0.3062 | 7.32e-01 | 2.048
310’450 | 0.1531 | 1.34e-01 | 2.454
2’041’058 | 0.0765 | 2.03e-02 | 2.715

The symbol $\star$ denotes that no experimental convergence rate has been
computed at the first refinement level.

Next we conduct a second test of convergence, now in the unit cube domain
$\Omega=(0,1)^{3}$ and with mixed boundary conditions (displacement, normal
velocity and tangential vorticity are prescribed on the sides $x=0$, $y=0$,
$z=0$ and known normal stress, flux, and pressures are imposed on the
remainder of the boundary) considering closed-form solutions to the vorticity-
based Biot–Brinkman equations as follows

$\displaystyle\boldsymbol{u}(x,y,z)=\frac{1}{10}\begin{pmatrix}\sin(\pi[x+y+z])\\\
\cos(\pi[x^{2}+y^{2}+z^{2}])\\\
\sin(\pi[x+y+z])\cos(\pi[x+y+z])\end{pmatrix},\qquad p(x,y,z)=\sin(\pi
x)\cos(\pi y)\sin(\pi z),$
$\displaystyle\boldsymbol{v}(x,y,z)=\begin{pmatrix}\sin^{2}(\pi x)\sin(\pi
y)\sin(2\pi z)\\\ \sin(\pi x)\sin^{2}(\pi y)\sin(2\pi z)\\\ -[\sin(2\pi
x)\sin(\pi y)+\sin(\pi x)\sin(2\pi y)]\sin^{2}(\pi
z)\end{pmatrix},\qquad\boldsymbol{\omega}(x,y,z)={\sqrt{\frac{\nu}{\kappa}}}\mathop{\mathbf{curl}}\nolimits\boldsymbol{v},\qquad\varphi(x,y,z)=-\lambda\mathop{\mathrm{div}}\nolimits\boldsymbol{u}+\alpha
p.$

The 3D example uses generalised Taylor–Hood elements for the
displacement/total pressure pair, but differently from (4.1), we take one
polynomial degree higher for the approximation of velocity, vorticity and
fluid pressure. This to maintain an overall $O(h^{k+2})$ convergence. This
time we consider the parameter values $\nu=0.1$, $\lambda=100$, $\mu=10$,
$\kappa=10^{-3}$, $c_{0}=0.1$, $\alpha=0.1$, and in Table 5.3 we only show the
error decay measured in the total weighted norm that leads to the definition
of $\mathcal{B}_{3}$, that is
$\|\vec{\boldsymbol{x}}-\vec{\boldsymbol{x}_{h}}\|_{\varepsilon}=\sqrt{(\vec{\boldsymbol{x}}-\vec{\boldsymbol{x}_{h}})\cdot\mathcal{B}_{3}^{-1}(\vec{\boldsymbol{x}}-\vec{\boldsymbol{x}_{h}})}$.
This shows an asymptotic $h^{k+2}$ order of convergence. For sake of
illustration, we depict in Figure 5.1 the approximate solutions obtained after
4 levels of uniform refinement.

Figure 5.1: Accuracy test in 3D. Approximate solutions of the Biot–Brinkman
equations on a relatively coarse mesh. Displacements on the deformed
configuration, velocity streamlines, vorticity vectors, total and fluid
pressures computed with the lowest-order method.

We also explore the dependence of the order of convergence upon variation of
the base-line parameter values used in Table 5.3. We consider several sample
values of model parameters $\mu\in[1,10^{8}]$, $\lambda\in[1,10^{8}]$,
$\nu\in[10^{-6},1]$, $\kappa\in[10^{-8},1]$, $\alpha\in[0,1]$, and
$c_{0}\in[10^{-8},1]$. These ranges are encountered in typical applications of
poromechanics of subsurface flows and of linear Biot consolidation of soft
tissues 16, 24, 43, 47, 56. For the sake of brevity, however, we only report a
subset of representative results we obtained with $\mu=1$, $\alpha=1$,
$\lambda=\\{1,10^{8}\\}$, $\nu=\\{10^{-8},1\\}$, $\kappa=\\{10^{-8},1\\}$, and
$c_{0}=\\{10^{-8},1\\}$. In any case, the software in 21 is written such that
the reader might also run additional tests with other parameter values as per-
required. We consider the unit cube domain discretised into uniform
tetrahedral meshes. In particular, we consider four mesh resolutions, from one
up to four levels of uniform refinement of the unit cube discretised with two
tetrahedra. We use mixed boundary conditions, with $\Gamma$ being conformed by
the faces $x=0$, $y=0$, $z=0$ and $\Sigma$ the remainder of the boundary.

The results are shown in Figure 5.2 for $k=0$. Overall, the different
convergence curves confirm $O(h^{k+2})$ convergence in the weighted norm that
leads to the definition of $\mathcal{B}_{3}$ regardless of the parameter
values, as predicted by the theoretical analysis.

Figure 5.2: Error convergence curves in 3D for $k=0$ and different
combinations of the physical parameters values. The errors are measured in the
weighted norm that leads to the definition of $\mathcal{B}_{3}$.

###  5.3 Preconditioner robustness with respect to model parameters

Finally, we proceed to study the robustness of $\mathcal{B}_{j}^{-1}$ in
(3.17) with respect to varying model parameters and increasing mesh
resolution. We use the same combination of parameter values as in the previous
section, namely $\mu=1$, $\alpha=1$, $\lambda=\\{1,10^{8}\\}$,
$\nu=\\{10^{-8},1\\}$, $\kappa=\\{10^{-8},1\\}$, and $c_{0}=\\{10^{-8},1\\}$,
and five levels of uniform refinement of the unit cube discretised with two
tetrahedra. We used the same boundary conditions as in the previous section.
We used the definition of FE spaces in (4.1) for $k=0$. The action of the
different inverses arising in the diagonal blocks of $\mathcal{B}_{j}^{-1}$ in
(3.17) is implemented with the UMFPACK direct solver. For each combination of
parameter values and mesh resolution, we used these preconditioners to
accelerate the convergence of the MINRES iterative solver. Convergence is
claimed whenever the Euclidean norm of the (unpreconditioned) residual of the
whole system is reduced by a factor of $10^{6}$, and stopped otherwise if the
number of iterations reaches an upper bound of 500 iterations. The (discrete)
Laplacian operator required for $\mathcal{B}_{2},\mathcal{B}_{3}$ acts on
discontinuous pore pressure approximations so we use (for a given piecewise-
defined field $\eta$) the following form (see 13)

$\bigl{(}-\eta\Delta_{h}p_{h},q_{h}\bigr{)}_{0,\Omega}=\sum_{K\in\mathcal{T}_{h}}(\eta\nabla_{h}p_{h},\nabla_{h}p_{h})_{0,K}{+\sum_{e\in\mathcal{E}^{\mathrm{int}}_{h}}\bigl{\langle}\frac{\left\\{\\!\\!\left\\{\eta\right\\}\\!\\!\right\\}}{\left\\{\\!\\!\left\\{h_{e}\right\\}\\!\\!\right\\}}\left[\\!\left[p_{h}\right]\\!\right],\left[\\!\left[q_{h}\right]\\!\right]\bigr{\rangle}_{e}}+\sum_{e\in\mathcal{E}^{\Gamma}_{h}}\langle\frac{\eta}{h_{e}}{p_{h}},{q_{h}}\rangle_{e},$
(5.2)

while for the lowest-order case we only keep the second and third terms on the
right-hand side of (5.2).

In Figures 5.3-5.4 we show the number of preconditioned MINRES iterations
versus number of DoFs for the $\mathcal{B}_{1}$, $\mathcal{B}_{2}$, and
$\mathcal{B}_{3}$ preconditioners with the particular parameter value
combinations mentioned above. Each of the plots in these figures contain 4
curves corresponding each to one of the four combinations of
$\kappa=\\{10^{-8},1\\}$, and $c_{0}=\\{10^{-8},1\\}$. To facilitate the
comparison among $\mathcal{B}_{1}$, $\mathcal{B}_{2}$, and $\mathcal{B}_{3}$,
the 6 plots in each figure are grouped into two groups of three horizontally
adjacent plots each. Each group correspond to a combination of
$\lambda=\\{1,10^{8}\\}$, $\nu=\\{10^{-8},1\\}$; the particular combination
corresponding to a group is indicated in the title of the three plots in the
group.

Figure 5.3: Comparison of parameter robustness for preconditioners
$\mathcal{B}_{1}$, $\mathcal{B}_{2}$, and $\mathcal{B}_{3}$.

Figure 5.4: Comparison of parameter robustness for preconditioners
$\mathcal{B}_{1}$, $\mathcal{B}_{2}$, and $\mathcal{B}_{3}$.

From Figures 5.3-5.4, we observe that the number of MINRES iterations with
$\mathcal{B}_{3}$ reaches an asymptotically constant regime with increasing
mesh resolution for all combinations of parameter values tested. This is also
the case of $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ in the majority of cases,
except for some combinations of parameter values in which $\nu=10^{-8}$ (see,
e.g., $\mathcal{B}_{1}$ and $\mu=1$, $\lambda=1$, $\alpha=1$, $\nu=10^{-8}$,
$\kappa=1$, and $c_{0}=\\{10^{-8},1\\}$), where preconditioner efficiency
(number of MINRES iterations) significantly degrades (increases) with mesh
resolution. (We recall from Remark 4.9 that in the limit $\nu\rightarrow 0$
the system at hand becomes a four-field Biot system 16.) For these
combinations of parameters, the coupling among the two pressures in the last
leading block of $\mathcal{B}_{3}$ is essential to retain mesh independence
convergence. This observation agrees with the experiments in 16 for four-field
formulations of Biot poroelasticity and simpler problems (including Herrmann
elasticity and reaction-diffusion equation), where comparisons were performed
against sub-optimal preconditioners. On the other hand, we observe a
relatively low sensitivity of the number of MINRES iterations (and thus
robustness) with respect to the value of the model of parameters for all
preconditioners (leaving aside the aforementioned combination of parameter
values). For example, for $\mu=1$, $\lambda=10^{8}$, $\alpha=1$,
$\nu=10^{-8}$, $\mathcal{B}_{3}$, and finest mesh resolution, the number of
iterations varies between 40 and 50 despite the disparity of scales in the
values of $c_{0}$ and $\kappa$. It is worth noting that $\mathcal{B}_{1}$ and
$\mathcal{B}_{2}$ lead to a similar or even lower number of iterations than
$\mathcal{B}_{3}$ in most cases, despite it being a computationally cheaper
preconditioner. This is the case, e.g., for and $\mathcal{B}_{1}$, $\mu=1$,
$\lambda=1$, $\alpha=1$, $\nu=1$, and all combinations of $\kappa$ and $c_{0}$
tested.

##  6 Concluding remarks

We have presented a new formulation for the Biot–Brinkman problem using
rescaled vorticity and total pressure, and have carried out the stability and
solvability analysis of the continuous and discrete problems using parameter-
weighted norms. We have derived theoretical error estimates and have
subsequently confirmed them numerically; and we have constructed suitable
preconditioners that achieve mesh-robustness iteration numbers when varying
elastic and porous media flow model parameters. In order to be able to apply
these algorithms to more realistic problems one needs to efficiently exploit
the vast amount of hardware parallelism available in high-end supercomputers.
As future work, we plan to address the parallelisation of the algorithms
proposed using the GridapDistributed.jl package 12.

Parts of the theoretical framework advanced herein (in particular, the use of
a vorticity-based formulation for the filtration equation) extend naturally to
more complex setups, for example to the multiple network generalised
Biot–Brinkman model from 35. Further improvements to the model and to the
theory include the rigorous treatment of different types of boundary
conditions, the interfacial coupling with free-flow or with elasticity, and
robust a posteriori error estimates.

*Financial disclosure

None reported.

*Conflict of interest

The authors declare no potential conflict of interests.

## References

  * 1 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.
  * 2 M. Amara, E. Chacón Vera, and D. Trujillo, Vorticity–velocity–pressure formulation for Stokes problem. Math. Comp., 73(248) (2004) 1673–1697.
  * 3 I. Ambartsumyan, E. Khattatov, I. Yotov, and P. Zunino, Simulation of flow in fractured poroelastic media: A comparison of different discretization approaches. In: I. Dimov, I. Faragó, and L. Vulkov (eds). Finite Difference Methods,Theory and Applications. FDM 2014. Lecture Notes in Computer Science, vol 9045. Springer, Cham (2015) 3–14.
  * 4 V. Anaya, M. Bendahmane, D. Mora, and R. Ruiz-Baier, On a vorticity-based formulation for reaction-diffusion-Brinkman systems. Netw. Heterog. Media, 13(1) (2018) 69–94.
  * 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, R. Caraballo, B. Gómez-Vargas, D. Mora, and R. Ruiz-Baier, Velocity-vorticity-pressure formulation for the Oseen problem with variable viscosity. Calcolo, 58(4) (2021) e44(1–25).
  * 7 V. Anaya, R. Caraballo, S. Caucao, L.F. Gatica, R. Ruiz-Baier, and I. Yotov, A vorticity-based mixed formulation for the unsteady Brinkman–Forchheimer equations. Comput. Methods Appl. Mech. Engrg., 404 (2023) e115829(1–30).
  * 8 V. Anaya, G.N. Gatica, D. Mora, and R. Ruiz-Baier, An augmented velocity-vorticity-pressure formulation for the Brinkman equations. Int. J. Numer. Methods Fluids, 79(3) (2015) 109–137.
  * 9 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.
  * 10 S. Badia, A. F. Martín, and R. Planas, Block recursive LU preconditioners for the thermally coupled incompressible inductionless MHD problem. J. Comput. Phys., 274 (2014), 562–591.
  * 11 S. Badia and F. Verdugo, Gridap: An extensible finite element toolbox in julia. J. Open Source Softw., 5 (2020), 2520.
  * 12 S. Badia, A. F. Martín, and F. Verdugo, GridapDistributed: a massively parallel finite element toolbox in Julia, J. Open Source Softw., 7(74) (2022), 4157.
  * 13 T. Bærland, M. Kuchta, K.-A. Mardal, and T. Thompson, An observation on the uniform preconditioners for the mixed Darcy problem, Numer. Methods PDEs, 36(6) (2020), 1718–1734.
  * 14 C. Bernardi and N. Chorfi, Spectral discretization of the vorticity, velocity, and pressure formulation of the Stokes problem. SIAM J. Numer. Anal., 44(2) (2007) 826–850.
  * 15 C. Bernardi, S. Dib, V. Girault, F. Hecht, F. Murat, and T. Sayah, Finite element methods for Darcy’s problem coupled with the heat equation. Numer. Math., 139 (2018) 315–348.
  * 16 W. Boon, M. Kuchta, K.-A. Mardal, and R. Ruiz-Baier, Robust preconditioners and stability analysis for perturbed saddle-point problems – application to conservative discretizations of Biot’s equations utilizing total pressure, SIAM J. Sci. Comput., 43 (2021) B961–B983.
  * 17 D. Braess, Stability of saddle point problems with penalty. RAIRO Modél. Math. Anal. Numér., 30 (1996) 731–742.
  * 18 F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods. Springer-Verlag, 1991.
  * 19 J. Brown, M. G. Knepley, D. A. May, L. C. McInnes, and B. Smith, Composable linear solvers for multiphysics, In 2012 11th International Symposium on Parallel and Distributed Computing, (2012) 55–62
  * 20 R. Bürger, S. Kumar, D. Mora, R. Ruiz-Baier, and N. Verma, Virtual element methods for the three-field formulation of time-dependent linear poroelasticity. Adv. Comput. Math., 47 (2021) e2.
  * 21 R. Caraballo, C. W. In, A. F. Martín, R. Ruiz-Baier, Software used in “Robust finite element methods and solvers for the Biot-Brinkman equations in vorticity form”, DOI: https://doi.org/10.5281/zenodo.10078398 (2023)
  * 22 F.J. Carrillo and I.C. Bourg, Capillary and viscous fracturing during drainage in porous media. Phys. Rev. E, 103 (2021) e063106.
  * 23 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.
  * 24 S. Chen, Q. Hong, J. Xu, and K. Yang, Robust block preconditioners for poroelasticity. Comput. Methods Appl. Mech. Engrg., 369 (2020) e113229.
  * 25 E. Chow, and M. A. Heroux, An object-oriented framework for block preconditioning. ACM Trans. Math. Softw., 24(2) (1998) 159–183.
  * 26 E.C Cyr, J. N. Shadid, and R. S. Tuminaro, Teko: A Block Preconditioning Capability with Concrete Example Applications in Navier–Stokes and MHD. SIAM J. Sci. Comput., 38(5) (2016) S307–S331.
  * 27 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.
  * 28 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.
  * 29 Y. Efendiev, J. Galvis, and Y. Vassilevski, Preconditioning of coupled systems and applications. Numer. Linear Alg. Appl., 16(7–8) (2009) 899–924.
  * 30 A. Ern and J.-L. Guermond, Finite Elements I: Approximation and Interpolation. Vol. 72 of Texts in Applied Mathematics. Springer, Cham (2021).
  * 31 G.N. Gatica, A Simple Introduction to the Mixed Finite Element Method. Theory and Applications. SpringerBriefs in Mathematics. Springer, Cham, 2014.
  * 32 L. Giraud, C. Geuzaine, and J. Dominguez, An efficient preconditioner for elasto-poroelasticity based on the pressure-correction method. Int. J. Numer. Methods Engrg., 89(10) (2011) 1139–1164.
  * 33 V. Girault and P.A. Raviart, Finite Element Methods for Navier–Stokes Equations. Theory and Algorithms. Springer, Berlin (1986).
  * 34 M.-L. Hanot, An arbitrary order and pointwise divergence-free finite element scheme for the incompressible 3D Navier–Stokes equations. SIAM J. Numer. Anal., 61(2) (2023), 784–811.
  * 35 Q. Hong, J. Kraus, M. Kuchta, M. Lymbery, K.-A. Mardal and M. E. Rognes, Robust approximation of generalized Biot-Brinkman problems, J. Sci. Comput., 93 (2022) e77.
  * 36 Q. Hong, J. Kraus, M. Lymbery, and F. Philo, Conservative discretizations and parameter-robust preconditioners for Biot and multiple-network flux-based poroelasticity models, Numer. Linear Algebra Appl., 26 (2019) e2242.
  * 37 S.K. Jha, Y Efendiev, J. Galvis, and Y. Vassilevski, Block-diagonal preconditioning for coupled systems in subsurface flow simulation. J. Comput. Phys., 299 (2015) 203–224.
  * 38 G. Ju, M. Cai, J. Li, and J. Tian, Parameter-robust multiphysics algorithms for Biot model with application in brain edema simulation. Math. Comput. Simul., 177 (2020) 385–403.
  * 39 R.C. Kirby, From functional analysis to iterative methods. SIAM Rev., 52(2) (2010) 269–293.
  * 40 R.C. Kirby and L. Mitchell, Solver Composition Across the PDE/Linear Algebra Barrier. SIAM J. Sci. Comput. 40(1) (2018) C76–C98.
  * 41 S. Kumar, R. Oyarzúa, R. Ruiz-Baier, and R. Sandilya, Conservative discontinuous finite volume and mixed schemes for a new four-field formulation in poroelasticity, ESAIM: Math. Model. Numer. Anal., 54 (2020) 273–299.
  * 42 J. Lee, K.-A. Mardal, and R. Winther, Parameter-robust discretization and preconditioning of Biot’s consolidation model. SIAM J. Sci. Comput., 39 (2017) A1–A24.
  * 43 J. J. Lee, E. Piersanti, K.-A. Mardal, and M. E. Rognes, A mixed finite element method for nearly incompressible multiple-network poroelasticity. SIAM J. Sci. Comput. 41 (2019) A722–A747.
  * 44 P. Lenarda, M. Paggi, and R. Ruiz-Baier, Partitioned coupling of advection-diffusion-reaction systems and Brinkman flows. J. Comput. Phys., 344 (2017) 281–302.
  * 45 X. Liu, Y. Zhang, and Z. Chen, A block-diagonal preconditioner for the solution of coupled PDEs in poromechanics. Numer. Methods PDEs, 36(4) (2020) 1407–1425.
  * 46 K.-A. Mardal and R. Winther, Preconditioning discretizations of systems of partial differential equations, Numer. Linear Alg. Appl., 18 (2011) 1–40.
  * 47 A. Mikelic, M.F. Wheeler, and T. Wick, Phase-field modeling of a fluid-driven fracture in a poroelastic medium. Comput. Geosci., 19(6) (2015) 1171–1195.
  * 48 P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford University Press, New York, 2003.
  * 49 D. Orban, and A. S. Siqueira, LinearOperators.jl, DOI: https://doi.org/10.5281/zenodo.2559294 (2020)
  * 50 R. Oyarzúa and R. Ruiz-Baier, Locking-free finite element methods for poroelasticity. SIAM J. Numer. Anal., 54 (2016) 2951–2973.
  * 51 E. Piersanti, J.J. Lee, T. Thompson, K.-A. Mardal, and M.E. Rognes, Parameter robust preconditioning by congruence for multiple-network poroelasticity. SIAM J. S. Comput., 43 (2021) B984–B1007.
  * 52 W. Qi, P. Seshaiyer, and J. Wang, Finite element method with the total stress variable for Biot’s consolidation model. Numer. Methods PDEs, 37(3) (2021) 2409–2428.
  * 53 A. Quarteroni, Numerical Models for Differential Problems. Springer-Verlag Milano (2009).
  * 54 K.R. Rajagopal, On a hierarchy of approximate models for flows of incompressible fluids through porous solids. Math. Model. Methods Appl. Sci., 17 (2007) 215–252.
  * 55 E. Rohan, J. Turjanicová, and V. Lukeš, The Biot–Darcy–Brinkman model of flow in deformable double porous media; homogenization and numerical modelling. Comput. Math. Appl., 78 (2019) 3044–3066.
  * 56 R. Ruiz-Baier, M. Taffetani, H.D. Westermeyer, and I. Yotov, The Biot–Stokes coupling using total pressure: formulation, analysis and application to interfacial flow in the eye. Comput. Methods Appl. Mech. Engrg., 389 (2022) e114384.
  * 57 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.
  * 58 N. Verma, B. Gómez-Vargas, L. M. De Oliveira Vilaca, S. Kumar, and R. Ruiz-Baier, Well-posedness and discrete analysis for advection-diffusion-reaction in poroelastic media, Appl. Anal., 101(14) (2022) 4914–4941.
  * 59 J. Zhang, C. Zhou, Y. Cao, and A.J. Meir, A locking free numerical approximation for quasilinear poroelasticity problems. Comput. Math. Appl., 80(16) (2020) 1538–1554.