# Analysis of a semi-augmented mixed finite element method for double-
diffusive natural convection in porous media ††thanks: This work was partially
supported by Vicerrectoría de Investigación, through the project C00-89, Sede
de Occidente, Universidad de Costa Rica; by CONICYT-Chile through the project
1190241; and by Universidad Nacional, Costa Rica, through the project 0140-20.

Mario Alvarez Eligio Colmenares Filánder A. Sequeira Sección de Matemática,
Sede de Occidente, Universidad de Costa Rica, San Ramón, Alajuela, Costa Rica,
email<EMAIL_ADDRESS>de Ciencias Básicas,
Facultad de Ciencias, Universidad del Bío-Bío, Campus Fernando May, Chillán,
Chile, email<EMAIL_ADDRESS>de Matemática, Universidad
Nacional, Campus Omar Dengo, Heredia, Costa Rica, email:
<EMAIL_ADDRESS>
###### Abstract

In this paper we study a stationary double-diffusive natural convection
problem in porous media given by a Navier-Stokes/Darcy type system, for
describing the velocity and the pressure, coupled to a vector advection-
diffusion equation describing the heat and substance concentration, of a
viscous fluid in a porous media with physical boundary conditions. The model
problem is rewritten in terms of a first-order system, without the pressure,
based on the introduction of the strain tensor and a nonlinear pseudo-stress
tensor in the fluid equations. After a variational approach, the resulting
weak model is then augmented using appropriate redundant penalization terms
for the fluid equations along with a standard primal formulation for the heat
and substance concentration. Then, it is rewritten as an equivalent fixed-
point problem. Well-posedness and uniqueness results for both the continuous
and the discrete schemes are stated, as well as the respective convergence
result under certain regularity assumptions combined with the Lax-Milgram
theorem, and the Banach and Brouwer fixed-point theorems. In particular,
Raviart-Thomas elements of order $k$ are used for approximating the pseudo-
stress tensor, piecewise polynomials of degree $\leq k$ and $\leq k+1$ are
utilized for approximating the strain tensor and the velocity, respectively,
and the heat and substance concentration are approximated by means of Lagrange
finite elements of order $\leq k+1$. Optimal a priori error estimates are
derived and confirmed through some numerical examples that illustrate the
performance of the proposed semi-augmented mixed-primal scheme.

Keywords: Double-diffusive natural convection, Oberbeck-Boussinesq model,
augmented formulation, mixed-primal finite element method, fixed point theory,
a priori error analysis.

Mathematical subject classifications (2020): 65N30, 65N12, 65N15, 76D05,
76R10, 76M10, 80A19.

##  1 Introduction

In nature and several technological applications, transport phenomena widely
occur as a result of a combined heat and mass transfer that are driven by
buoyancy effects due to both temperature and concentration variations (see,
e.g., [15, 37, 39]). Such processes, also known as thermosolutal or double-
diffusive natural convection, involving fluid circulation in a porous media,
are frequently found in astrophysics, oceanology, metallurgy, electrophysics
and geophysics, but also appear in several engineering applications such as
filtration processes, geothermal energy exploitation, spreading on porous
substrates, bio-film growth, gasification of biomass, to name a few.

From the mathematical point of view the Darcy-Oberbeck-Boussinesq model allows
to adequately describe and quantify this complex flow by means of a nonlinear
partial differential equations system. More precisely, the momentum and
conservation of fluid mass give rise to a Navier-Stokes/Darcy type system for
describing the fluid flow in the porous media which, in turn, is coupled via
buoyancy forces and convective mass and heat transfer to a vector advection-
diffusion equation for describing the substance concentration and the
temperature, as a result of an energy and mass transfer balance (see, e.g.,
[37, 39]).

Many computational techniques have been developed so far in order to
numerically solve and simulate this problem and related ones (see [1, 3, 10,
17, 18, 29, 30, 34, 38, 44], and the references therein). Particularly, the
contributions [1, 3, 38, 29, 30, 34, 44] deal with double-diffusive convection
in a cavity, whereas in [10, 17, 18] the authors consider the phenomenon in a
porous media.

In [38], a finite volume method is proposed and applied to solve agro-food
processes, whereas some methods based on finite elements for this problem are
[1, 3]. In [1], the authors proposed stabilized finite element formulations
based on the SUPG (Streamline-Upwind/Petrov-Galerkin) and PSPG (Pressure-
Stabilization/Petrov-Galerkin) methods to solve the problem in unsteady state.
Numerical simulations in two and three dimensions illustrate the accuracy and
performance of this technique. However, the theoretical analyses of the
associate continuous and the discrete variational problems as well as the
convergence of the method are not carried out there, and the method only
allows to carry out low-order approximations of the main unknowns. On the
other hand, in [3] the problem is considered in steady state and analyzed by
using the boundary control theory. The authors formulate and prove solvability
results for the corresponding boundary control problem, state local uniqueness
and stability of optimal solutions.

Focusing in double-diffusive viscous flow in porous media, [17] proposes a
technique consisting of a projection-based stabilization method in the
unsteady state. There, the convergence of the velocity, temperature and
concentration in the semi-discrete case are derived. In addition, some
numerical experiments are reported to confirm optimal order error estimates
and to compare their results with previous ones. In [10] the authors construct
a divergence-conforming primal scheme and establish the existence and
uniqueness results for the continuous problem and the discrete scheme as well
as convergence properties. On the other hand, in [18] the authors propose a
high-order fully-mixed method based on the introduction of the velocity
gradient, the temperature gradient and the concentration gradient as new
unknowns into the problem. The resulting formulation has a saddle-point
structure on reflexive Banach spaces for both the Navier-Stokes/Darcy and the
thermal energy conservation equations. There, it is particularly shown that
the discrete scheme is well-posed and an a priori error estimate for the
Galerkin scheme is also derived under sufficiently small data. However,
feasible finite element subspaces must be constructed over meshes with a
macro-element structure in order to satisfy an inf-sup compatibility
condition, which in turn significantly increased the computational cost,
especially in the three-dimensional case.

According to the above discussion and in order to contribute with the design
and analysis of new mixed finite schemes for simulating double-diffusive
convection in porous media. The main goal of this work is precisely to propose
a new semi-augmented mixed finite element method in which the strain tensor
and a pseudo-stress tensor are introduced as primary unknowns of interest in
the fluid equations and the pressure is eliminated from the system by its own
definition. To avoid any inf-sup restriction, to guarantee greater flexibility
regarding finite element spaces and lower computational cost than [18], the
respective variational formulation is augmented by using appropriate redundant
penalization Galerkin type-terms based on the constitutive and equilibrium
equations combined with a primal formulation for the heat and substance
concentration equations in standard Hilbert spaces. In this way, the
aforementioned strain tensor, the nonlinear pseudostress, the velocity and a
vector field whose components are the temperature of the fluid and the
substance concentration are the main unknowns of the resulting coupled system.
Moreover, physical boundary conditions are considered. Indeed, a no-slip
condition (that is, homogeneous Dirichlet condition) for the fluid velocity, a
prescribed temperature and substance concentration on a Dirichlet boundary and
no heat/mass flow across an isolated surface/homogeneous Neumann condition.

Concerning the solvability analysis, we proceed similarly to [4, 19] by
introducing an equivalent fixed-point setting. According to this, combining
the Lax-Milgram theorem with the classical Banach and Brouwer fixed-point
theorems, we establish the respective solvability of the continuous problem
and the associated Galerkin scheme, under suitable regularity assumptions (see
[7], for further details), a feasible choice of parameters and, in the
discrete case, for any family of finite element subspaces. To handle the non-
homogeneous Dirichlet condition for the temperature and concentration, we
carry out a rigorous treatment of the boundary data throughout the analysis by
means of appropriate extensions involving the Scott-Zhang interpolator (in the
discrete case), which allows us to establish the well-posedness of our scheme,
along with its convergence result and the respective a priori error bounds. Up
to the best of our knowledge, because of the difficulties that can arise in
the analysis, the physically relevant non-homogeneous Dirichlet condition case
is usually either omitted or not considered; this is what motivates us to
contribute in this direction as well.

A Strang-type lemma, enables us to derive the corresponding Céa estimate and
to provide optimal a priori error bounds for the Galerkin solution. In turn,
the pressure can be recovered by a post-processed of the discrete solutions,
preserving the same rate of convergence. Finally, some numerical experiments
are presented to illustrate the performance of the technique, including the
unsteady case with unknown solution, to confirm the expected orders.

The contents of this paper is presented as follows. At the end of this
section, we introduce some standard notations to be used throughout the
manuscript. In Section 2, we introduce the model problem and the data. We also
derive an equivalent first-order equations in terms of additional variables.
Then, in Section 3, we derive the semi-augmented mixed-primal variational
formulation and establish its well-posedness. The associated Galerkin scheme
is introduced and analyzed in Section 4. In Section 5, we derive the
corresponding Céa estimate and, finally, in Section 6 we present some
numerical examples illustrating the performance of our semi-augmented mixed-
primal finite element method.

###  1.1 Notations

Let us denote by $\Omega\subseteq\mathrm{R}^{n}$, with $n\in\\{2,3\\}$, a
given bounded domain with polygonal/polyhedral boundary $\Gamma$ with outward
unit normal vector ${\boldsymbol{n}}$ and let
$\Gamma_{\mathrm{D}}\,,\Gamma_{\mathrm{N}}\,\subseteq\,\Gamma$ be such that
$\Gamma_{\mathrm{D}}\,\cap\,\Gamma_{\mathrm{N}}\,=\,\varnothing$,
$|\Gamma_{\mathrm{D}}|\,>\,0$ and
$\Gamma\,=\,\overline{\Gamma}_{\mathrm{D}}\,\cup\,\overline{\Gamma}_{\mathrm{N}}$.
Standard notation will be adopted for Lebesgue spaces $\mathrm{L}^{p}(\Omega)$
with norm $\|\cdot\|_{0,p,\Omega}$ or $\|\cdot\|_{0,\Omega}$ if $p=2$, and
Sobolev spaces $\mathrm{H}^{s}(\Omega)$ with norm $\|\cdot\|_{s,\Omega}$, and
semi-norm $|\cdot|_{s,\Omega}$. In particular, when $\mathrm{A}$ denotes a
generic scalar functional space, then we will denote its vectorial and tensor
counterparts by $\mathbf{A}$ and $\mathbb{A}$, respectively.

For vector fields $\boldsymbol{v}=(v_{i})_{1\leq i\leq n}$ and
$\boldsymbol{w}=(w_{i})_{1\leq i\leq n}$, we set the gradient, divergence and
dyadic product operators, as

$\nabla\boldsymbol{v}:=\left(\frac{\partial v_{i}}{\partial
x_{j}}\right)_{1\leq i,j\leq
n},\quad\mathrm{div}\,\boldsymbol{v}:=\sum_{j=1}^{n}\frac{\partial
v_{j}}{\partial
x_{j}}\,,\quad\mathrm{and}\quad\boldsymbol{v}\otimes\boldsymbol{w}:=(v_{i}w_{j})_{1\leq
i,j\leq n}\,,$

respectively. Furthermore, given the tensor fields
$\boldsymbol{\tau}=(\tau_{ij})_{1\leq i,j\leq n}$ and
$\boldsymbol{\zeta}=(\zeta_{ij})_{1\leq i,j\leq n}$, we let
$\mathbf{div}\,\boldsymbol{\tau}$ be the divergence operator $\mathrm{div}$
acting along the rows of $\boldsymbol{\tau}$, and define the transpose, the
trace, the tensor inner product, and the deviatoric tensor, respectively, as

$\boldsymbol{\tau}^{\mathrm{t}}\ :=\ (\tau_{ji})_{1\leq i,j\leq
n},\quad\mathrm{tr}(\boldsymbol{\tau}):=\sum_{j=1}^{n}\tau_{ii},\quad\boldsymbol{\tau}:\boldsymbol{\zeta}:=\sum_{j=1}^{n}\tau_{ij}\zeta_{ij}\,,\quad\mathrm{and}\quad\boldsymbol{\tau}^{\mathrm{d}}:=\boldsymbol{\tau}-\frac{1}{n}\mathrm{tr}(\boldsymbol{\tau})\mathbb{I},$

where $\mathbb{I}$ stands for the identity tensor in $\mathrm{R}^{n\times n}$.
We recall that the space

$\mathbb{H}(\mathbf{div};\Omega)\ :=\
\Big{\\{}\boldsymbol{\tau}\in\mathbb{L}^{2}(\Omega):\quad\mathbf{div}\,\boldsymbol{\tau}\in\mathbf{L}^{2}(\Omega)\Big{\\}}\,,$

equipped with the usual norm

$\|\boldsymbol{\tau}\|^{2}_{\mathbf{div},\Omega}:=\|\boldsymbol{\tau}\|^{2}_{0,\Omega}+\|\mathbf{div}\,\boldsymbol{\tau}\|^{2}_{0,\Omega},$

is a Hilbert space. Finally, we employ $\boldsymbol{0}$ to denote a generic
null vector and use $c$ or $C$, with or without subscripts, bars, tildes or
hats, to mean generic positive constants independent of the discretization
parameters, which may take different values at different places.

##  2 The model problem

This section introduces the mathematical model we address in the present work
as well as the auxiliary unknowns that are introduced and considered in the
subsequent variational formulation. Under the Oberbeck-Boussinesq
approximation framework, double-diffusive natural convection phenomenon in a
porous media is described in terms of a Navier-Stokes-Brinkman model coupled
to a system of diffusion-advection equations. In the stationary state, the
problem consists of finding the velocity $\boldsymbol{u}$, the pressure $p$,
and the vector $\boldsymbol{\varphi}=(\varphi_{1},\varphi_{2})^{\mathrm{t}}$
when $n=2$ and $\boldsymbol{\varphi}=(\varphi_{1},\varphi_{2},0)^{\mathrm{t}}$
when $n=3$, where $\varphi_{1}$ and $\varphi_{2}$ are the temperature and the
concentration fields, respectively, of a fluid in a confined porous enclosure
$\Omega$, satisfying the set of equations:

$\begin{array}[]{rcl}\gamma\,\boldsymbol{u}-2\,\mathbf{div}\left(\nu(\boldsymbol{\varphi})\mathbf{e}(\boldsymbol{u})\right)+\left(\boldsymbol{u}\cdot\nabla\right)\boldsymbol{u}+\nabla\,p&=&\left(\boldsymbol{\alpha}\cdot\boldsymbol{\varphi}\right){\boldsymbol{g}}{\quad\hbox{in}\quad}\Omega\,,\\\\[4.30554pt]
\mathrm{div}\,\boldsymbol{u}&=&0{\quad\hbox{in}\quad}\Omega\,,\\\\[4.30554pt]
-{\mathbf{div}}(\mathbb{K}\nabla\boldsymbol{\varphi})\,+\,(\nabla\boldsymbol{\varphi})\boldsymbol{u}&=&\mathbf{0}{\quad\hbox{in}\quad}\Omega\,,\end{array}$
(2.1)

where $\mathbf{e}(\boldsymbol{u})$ stands for the symmetric part of the
velocity gradient, i.e.,
$\mathbf{e}(\boldsymbol{u}):=\frac{1}{2}(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{\mathrm{t}})$
. The data are the gravitational force $\boldsymbol{g}$, the positive constant
$\gamma$ corresponding to the reciprocal of the Darcy number, the (thermal and
solute) expansion coefficients vector $\boldsymbol{\alpha}$, and the diagonal
matrix of (thermal and mass) diffusion constants
$\mathbb{K}=\mathrm{diag}(k_{i})_{1\leq i\leq
n}\in\mathbb{L}^{\infty}(\Omega)$, with $k_{i}=0$ when $i=3$, which is assumed
to be a uniformly positive definite tensor, which means that there exists a
positive constant $k_{0}$ such that

$(\mathbb{K}\,\boldsymbol{x})\cdot\boldsymbol{x}\,\geq\,k_{0}\,|\boldsymbol{x}|^{2}\,\qquad\forall\,\boldsymbol{x}\in\mathrm{R}^{n}\,,\quad\mbox{with}\quad
k_{0}=\min\\{k_{1},k_{2},\dots,k_{n}\\}\,.$ (2.2)

In turn, the kinematic viscosity
$\nu:\mathrm{R}^{+}\times\mathrm{R}^{+}\longrightarrow\mathrm{R}^{+}$ is
assumed to be a bounded and Lipschitz continuous function that might depend on
both the temperature and the mass concentration. That is, we assume the
existence of positive constants $\nu_{1},\,\nu_{2},$ and $L_{\nu}$ that
satisfy

$\nu_{1}\,\leq\,\nu(\boldsymbol{\phi})\,\leq\,\nu_{2}{\quad\hbox{and}\quad}|\nu(\boldsymbol{\phi})-\nu(\boldsymbol{\psi})|\,\leq\,L_{\nu}\,|\boldsymbol{\phi}-\boldsymbol{\psi}|\qquad\forall\,\boldsymbol{\phi},\,\boldsymbol{\psi}\in\mathrm{R}\times\mathrm{R}^{+}\,.$
(2.3)

The system (2.1) is finally completed with a non-slip condition on the whole
boundary for the velocity and physical boundary conditions for both the
temperature and the concentration fields, that is

$\boldsymbol{u}\,=\,\mathbf{0}\,\hskip 5.69046pt\mbox{on
$\Gamma$\,,}\qquad\displaystyle\boldsymbol{\varphi}\,=\,\boldsymbol{\varphi}_{{\mathrm{D}}}\,\,\,\mbox{on
${\Gamma_{\mathrm{D}}}$}\qquad\mbox{and}\qquad\displaystyle(\nabla\boldsymbol{\varphi})\,{\boldsymbol{n}}=\mathbf{0}\,\,\,\mbox{on}\,\,{\Gamma_{\mathrm{N}}}\,,$
(2.4)

where
$\boldsymbol{\varphi}_{\mathrm{D}}\in\mathbf{H}^{1/2}(\Gamma_{\mathrm{D}})$ is
a known trace function on ${\Gamma_{\mathrm{D}}}.$

Next, proceeding similarly to [13], we introduce the strain tensor
$\boldsymbol{t}:=\mathbf{e}(\boldsymbol{u})$ as an auxiliary variable and then
define

$\boldsymbol{\sigma}\ :=\
2\nu(\boldsymbol{\varphi})\boldsymbol{t}-(\boldsymbol{u}\otimes\boldsymbol{u})-\big{(}p+c_{\boldsymbol{u}}\big{)}\mathbb{I}\quad\mathrm{in}\quad\Omega\,,$
(2.5)

as an additional tensorial unknown, where $c_{\boldsymbol{u}}$ is a constant
to be suitably defined next (see equation (2.7) below). Thus, noting that the
incompressibility condition of the fluid $({\mathrm{div}}\,\boldsymbol{u}=0)$
implies that
${\mathbf{div}}(\boldsymbol{u}\otimes\boldsymbol{u})=(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}$,
we get from the first relation of (2.1) the equilibrium equation

$\gamma\boldsymbol{u}-{\mathbf{div}}\,\boldsymbol{\sigma}\ =\
(\boldsymbol{\alpha}\cdot\boldsymbol{\varphi})\,{\boldsymbol{g}}\quad\mathrm{in}\quad\Omega\,,$

and taking deviatoric part in (2.5), we find that the constitutive relation
defined by (2.5) can be written as

$\boldsymbol{\sigma}^{\mathrm{d}}\ =\
2\nu(\boldsymbol{\varphi})\boldsymbol{t}-(\boldsymbol{u}\otimes\boldsymbol{u})^{\mathrm{d}}\,{\quad\hbox{in}\quad}\,\,\Omega\,$

where $\boldsymbol{t}^{\mathrm{d}}=\boldsymbol{t}$ since
${\mathrm{div}}\,\boldsymbol{u}=0.$ Thus, the pressure is eliminated from the
system, however taking trace in (2.5), we find that it can be easily recovered
according to the formula

$p\,=\,-\dfrac{1}{n}\Big{[}\mathrm{tr}(\boldsymbol{\sigma}+(\boldsymbol{u}\otimes\boldsymbol{u}))\Big{]}\,-\,c_{\boldsymbol{u}}\,,\qquad\mbox{in}\qquad\Omega\,.$
(2.6)

Now, since $p\in\mathrm{L}^{2}_{0}(\Omega)$ (which as usual is clearly
required for uniqueness of an eventual pressure solution to (2.1)), the
equation (2.6) suggests to define

$c_{\boldsymbol{u}}\,=\,-\dfrac{1}{n|\Omega|}\int_{\Omega}\mathrm{tr}(\boldsymbol{u}\otimes\boldsymbol{u})\,=\,-\dfrac{1}{n|\Omega|}\|\boldsymbol{u}\|_{0,\Omega}^{2}\,,$
(2.7)

so as to get the equivalence

$\int_{\Omega}p\,=\,0\qquad\mbox{if, and only
if,}\qquad\int_{\Omega}\mathrm{tr}\,{\boldsymbol{\sigma}}\,=0\,.$ (2.8)

According to the above discusion, the system (2.1) and (2.4) equivalently
reads: Find $\boldsymbol{t}$, ${\boldsymbol{\sigma}}$ and
$\boldsymbol{\varphi}$ such that

$\begin{array}[]{c}\displaystyle
2\nu(\boldsymbol{\varphi})\boldsymbol{t}-\boldsymbol{\sigma}^{\mathrm{d}}-(\boldsymbol{u}\otimes\boldsymbol{u})^{\mathrm{d}}\
=\
\mathbf{0}{\quad\hbox{in}\quad}\Omega\,,\quad\displaystyle\gamma\boldsymbol{u}-{\mathbf{div}}\,\boldsymbol{\sigma}\
=\
(\boldsymbol{\alpha}\cdot\boldsymbol{\varphi})\,{\boldsymbol{g}}{\quad\hbox{in}\quad}\Omega\,,\\\\[8.61108pt]
\boldsymbol{t}+\boldsymbol{\omega}(\boldsymbol{u})\ =\
\nabla\boldsymbol{u}{\quad\hbox{in}\quad}\Omega\,,\quad-{\mathbf{div}}(\mathbb{K}\nabla\boldsymbol{\varphi})\,+\,(\nabla\boldsymbol{\varphi})\boldsymbol{u}\
=\ \mathbf{0}{\quad\hbox{in}\quad}\Omega\,,\\\\[8.61108pt]
\displaystyle\boldsymbol{u}\ =\
\boldsymbol{0}{\quad\hbox{on}\quad}\Gamma\,,\quad\boldsymbol{\varphi}\ =\
\boldsymbol{\varphi}_{{\mathrm{D}}}{\quad\hbox{on}\quad}\Gamma_{\mathrm{D}}\,,\quad(\nabla\boldsymbol{\varphi})\,{\boldsymbol{n}}\
=\
\mathbf{0}{\quad\hbox{on}\quad}\Gamma_{\mathrm{N}}\,,{\quad\hbox{and}\quad}\int_{\Omega}\mathrm{tr}\,{\boldsymbol{\sigma}}\
=\ 0\,,\end{array}$ (2.9)

where
$\boldsymbol{\omega}(\boldsymbol{u}):=\frac{1}{2}(\nabla\boldsymbol{u}-(\nabla\boldsymbol{u})^{\mathrm{t}})$
is the skew-symmetric part of the velocity gradient. We emphasize that the
introduction of the variables $\boldsymbol{t}$ and ${\boldsymbol{\sigma}}$ as
new unknowns in the system allows us to equivalently rewrite the Navier-
Stokes-Brinkman model (first row of (2.1)) in terms of a first-order set of
equations. Also, observe that according to (2.6)-(2.8), the zero mean value
restriction of the pressure on the domain is imposed via the respective
restriction on $\mathrm{tr}\,{\boldsymbol{\sigma}}$ in the last equation of
(2.9). Note further that the incompressibility condition of the fluid is
implicitly present via the equilibrium relation defined by
${\boldsymbol{\sigma}}$ according to the second equation in the first row of
(2.9).

##  3 The continuous formulation

In this section we introduce and analyze the weak formulation proposed for the
system described by (2.9). In Section 3.1 we firstly deduce an augmented mixed
variational formulation of (2.9) and then we rewrite it as a fixed-point
problem in Section 3.2, whose analysis is addressed through Sections 3.3 and
3.4.

###  3.1 The semi-augmented mixed-primal variational formulation

To begin with, the fact that the trace of the tensor solution
${\boldsymbol{\sigma}}$ of system (2.9) has zero mean value in $\Omega$ (see
last equation of (2.9)) suggests to introduce the space

$\mathbb{H}_{0}(\mathbf{div};\Omega)\ :=\
\left\\{\boldsymbol{\tau}\in\mathbb{H}(\mathbf{div};\Omega):\quad\int_{\Omega}\mathrm{tr}\,\boldsymbol{\tau}=0\right\\}.$

In addition, by their own definitions, we introduce the following space for
the strain tensor $\boldsymbol{t}$, as

$\mathbb{L}_{\mathrm{tr}}^{2}(\Omega)\ :=\
\Big{\\{}\boldsymbol{r}\in\mathbb{L}^{2}(\Omega):\quad\boldsymbol{r}^{\mathrm{t}}=\boldsymbol{r}\quad\mathrm{and}\quad\mathrm{tr}\,\boldsymbol{r}=0\Big{\\}}.$

In turn, because of the boundary conditions for the temperature and
concentrations (see second and third equations of last row in (2.9)) we
consider the closed subspace

$\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)\ :=\
\Big{\\{}\boldsymbol{\psi}\in\mathbf{H}^{1}(\Omega):\quad\boldsymbol{\psi}\big{|}_{{\Gamma_{\mathrm{D}}}}=\mathbf{0}\,\Big{\\}}\,,$

for which, from the generalized Poincaré inequality, we know that there exists
$c_{\text{gp}}>0$, depending only on $\Omega$ and ${\Gamma_{\mathrm{D}}}$,
such that

$\|\boldsymbol{\psi}\|_{1,\Omega}\;\leq\;c_{\text{gp}}|\boldsymbol{\psi}|_{1,\Omega}\qquad\forall\
\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)\,.$ (3.1)

Now, testing the first equation from first row in (2.9) with
$\boldsymbol{r}\in\mathbb{L}_{\mathrm{tr}}^{2}(\Omega)$, we obtain

$2\int_{\Omega}\nu(\boldsymbol{\varphi})\boldsymbol{t}:\boldsymbol{r}-\int_{\Omega}\boldsymbol{\sigma}^{\mathrm{d}}:\boldsymbol{r}-\int_{\Omega}(\boldsymbol{u}\otimes\boldsymbol{u})^{\mathrm{d}}:\boldsymbol{r}\
=\ 0\qquad\forall\ \boldsymbol{r}\in\mathbb{L}_{\mathrm{tr}}^{2}(\Omega)\,.$
(3.2)

At this point, we readily note that in order to bound the third term on the
left hand side of (3.2), and thanks to the continuous injection
$\boldsymbol{i}:\mathbf{H}^{1}(\Omega)\rightarrow\mathbf{L}^{4}(\Omega)$ (see
e.g. [2] or [41]), we requiere the unknown $\boldsymbol{u}$ to live in
$\mathbf{H}^{1}_{0}(\Omega):=\\{\boldsymbol{v}\in\mathbf{H}^{1}(\Omega):\;\boldsymbol{v}|_{\Gamma}=\bf{0}\\}$.
Indeed, by applying Cauchy-Schwarz and Hölder inequalities, we deduce that
there exists a positive constant $c_{1}(\Omega):=\|\boldsymbol{i}\|^{2}$, such
that

$\left|\int_{\Omega}(\boldsymbol{u}\otimes\boldsymbol{w})^{\mathrm{d}}:\boldsymbol{r}\right|\,\leq
c_{1}(\Omega)\|\boldsymbol{u}\|_{1,\Omega}\|\boldsymbol{w}\|_{1,\Omega}\|\boldsymbol{r}\|_{0,\Omega}\,\quad\forall\,\boldsymbol{u},\boldsymbol{w}\in\mathbf{H}^{1}_{0}(\Omega),\qquad\forall\,\boldsymbol{r}\in\mathbb{L}^{2}(\Omega)\,.$
(3.3)

Next, multiplying the first equation from second row in (2.9) by a test
function $\boldsymbol{\tau}\in\mathbb{H}_{0}(\mathbf{div};\Omega)$,
integrating by parts, using the Dirichlet condition for $\boldsymbol{u}$, and
the identity
$\boldsymbol{t}:\boldsymbol{\tau}=\boldsymbol{t}:\boldsymbol{\tau}^{\mathrm{d}}$
(since $\boldsymbol{t}$ is trace-free), we get

$\int_{\Omega}\boldsymbol{t}:\boldsymbol{\tau}^{\mathrm{d}}+\int_{\Omega}\boldsymbol{\omega}(\boldsymbol{u}):\boldsymbol{\tau}+\int_{\Omega}\boldsymbol{u}\cdot\mathbf{div}\,\boldsymbol{\tau}\
=\ 0\qquad\forall\,\boldsymbol{\tau}\in\mathbb{H}_{0}(\mathbf{div};\Omega).$

Likewise, the equilibrium relation associated to ${\boldsymbol{\sigma}}$
(second equation from first row in (2.9)) is written as

$\gamma\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}-\int_{\Omega}\boldsymbol{v}\cdot\mathbf{div}\,\boldsymbol{\sigma}\
=\
\int_{\Omega}(\boldsymbol{\alpha}\cdot\boldsymbol{\varphi}){\boldsymbol{g}}\cdot\boldsymbol{v}\qquad\forall\
\boldsymbol{v}\in\mathbf{H}^{1}_{0}(\Omega)\,,$

whereas, the symmetry of the pseudo-stress tensor is imposed in an ultra-weak
sense (see e.g [8]) through the identity

$\displaystyle-\int_{\Omega}\boldsymbol{\sigma}:\boldsymbol{\omega}(\boldsymbol{v})\
=\ 0\,\qquad\forall\,\boldsymbol{v}\in\mathbf{H}^{1}_{0}(\Omega)\,.$

As for the equation associated to the temperature and concentration (second
equation from second row in (2.9)), we simply multiply it by a test function
$\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)$ and,
after integrating by parts, and employing the Neumann boundary contidions on
${\Gamma_{\mathrm{N}}}$, we find

$\int_{\Omega}\mathbb{K}\nabla\boldsymbol{\varphi}:\nabla\boldsymbol{\psi}\ +\
\int_{\Omega}(\nabla\boldsymbol{\varphi})\boldsymbol{u}\cdot\boldsymbol{\psi}\
=\
0\qquad\forall\,\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)\,.$

In this way, a preliminary weak formulation for the coupled problem (2.9)
takes the form: Find
$(\boldsymbol{t},\boldsymbol{\sigma},\boldsymbol{u},\boldsymbol{\varphi})\in\mathbb{L}_{\mathrm{tr}}^{2}(\Omega)\times\mathbb{H}_{0}(\mathbf{div};\Omega)\times\mathbf{H}^{1}_{0}(\Omega)\times\mathbf{H}^{1}(\Omega)$,
with
$\boldsymbol{\varphi}\big{|}_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{\mathrm{D}}$,
such that

$\begin{array}[]{rcll}\displaystyle{2\int_{\Omega}\nu(\boldsymbol{\varphi})\boldsymbol{t}:\boldsymbol{r}-\int_{\Omega}\boldsymbol{\sigma}^{\mathrm{d}}:\boldsymbol{r}-\int_{\Omega}(\boldsymbol{u}\otimes\boldsymbol{u})^{\mathrm{d}}:\boldsymbol{r}}&=&0&\forall\,\boldsymbol{r}\in\mathbb{L}_{\mathrm{tr}}^{2}(\Omega)\,,\\\\[12.91663pt]
\displaystyle{\int_{\Omega}\boldsymbol{t}:\boldsymbol{\tau}^{\mathrm{d}}+\int_{\Omega}\boldsymbol{\omega}(\boldsymbol{u}):\boldsymbol{\tau}+\int_{\Omega}\boldsymbol{u}\cdot\mathbf{div}\,\boldsymbol{\tau}}&=&0&\forall\,\boldsymbol{\tau}\in\mathbb{H}_{0}(\mathbf{div};\Omega)\,,\\\\[12.91663pt]
\displaystyle-\int_{\Omega}\boldsymbol{v}\cdot\mathbf{div}\,\boldsymbol{\sigma}-\int_{\Omega}\boldsymbol{\omega}(\boldsymbol{v}):\boldsymbol{\sigma}+\gamma\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}&=&\displaystyle\int_{\Omega}(\boldsymbol{\alpha}\cdot\boldsymbol{\varphi}){\boldsymbol{g}}\cdot\boldsymbol{v}&\forall\,\boldsymbol{v}\in\mathbf{H}^{1}_{0}(\Omega)\,,\\\\[12.91663pt]
\displaystyle\int_{\Omega}\mathbb{K}\nabla\boldsymbol{\varphi}:\nabla\boldsymbol{\psi}+\int_{\Omega}(\nabla\boldsymbol{\varphi})\boldsymbol{u}\cdot\boldsymbol{\psi}&=&0&\forall\,\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)\,.\end{array}$
(3.4)

In order to analyze the variational formulation (3.4), and similarly as in
[13, Section 2] (see also [4, 8, 21]), we additionally augment (3.4) by
incorporating the following residual Galerkin type-terms coming from the
constitutive and equilibrium equations for the fluid,

$\begin{array}[]{rcll}\displaystyle\kappa_{1}\int_{\Omega}(\mathbf{e}(\boldsymbol{u})-\boldsymbol{t}):\mathbf{e}(\boldsymbol{v})&=&0&\forall\,\boldsymbol{v}\in\mathbf{H}^{1}_{0}(\Omega)\,,\\\\[17.22217pt]
\displaystyle-\kappa_{2}\gamma\int_{\Omega}\boldsymbol{u}\cdot\mathbf{div}\,\boldsymbol{\tau}+\kappa_{2}\int_{\Omega}\mathbf{div}\,\boldsymbol{\sigma}\cdot\mathbf{div}\,\boldsymbol{\tau}&=&\displaystyle-\kappa_{2}\int_{\Omega}(\boldsymbol{\alpha}\cdot\boldsymbol{\varphi}){\boldsymbol{g}}\cdot\mathbf{div}\,\boldsymbol{\tau}&\forall\,\boldsymbol{\tau}\in\mathbb{H}_{0}(\mathbf{div};\Omega)\,,\\\\[17.22217pt]
\displaystyle\quad-\kappa_{3}\int_{\Omega}\left\\{2\nu(\boldsymbol{\varphi})\boldsymbol{t}-\boldsymbol{\sigma}^{\mathrm{d}}-(\boldsymbol{u}\otimes\boldsymbol{u})^{\mathrm{d}}\right\\}:\boldsymbol{\tau}^{\mathrm{d}}&=&0&\forall\,\boldsymbol{\tau}\in\mathbb{H}_{0}(\mathbf{div};\Omega)\,,\end{array}$
(3.5)

where $(\kappa_{1},\kappa_{2},\kappa_{3})$ is a vector of positive parameters
to be specified later in Section 3.3.

Hence, letting

$\vec{\boldsymbol{\xi}}\ :=\
(\boldsymbol{t},\boldsymbol{\sigma},\boldsymbol{u})\in\mathbb{H}:=\mathbb{L}_{\mathrm{tr}}^{2}(\Omega)\times\mathbb{H}_{0}(\mathbf{div};\Omega)\times\mathbf{H}^{1}_{0}(\Omega),$

where $\mathbb{H}$ is endowed with the natural norm

$\|\vec{\boldsymbol{\eta}}\|_{\mathbb{H}}:=\Big{\\{}\|\boldsymbol{r}\|_{0,\Omega}^{2}+\|\boldsymbol{\tau}\|_{\mathbf{div},\Omega}^{2}+\|\boldsymbol{v}\|_{1,\Omega}^{2}\Big{\\}}^{1/2},\quad\forall\,\vec{\boldsymbol{\eta}}:=(\boldsymbol{r},\boldsymbol{\tau},\boldsymbol{v})\in\mathbb{H}\,,$

and adding up (3.4) with (3.5), we arrive at the following semi-augmented
mixed-primal formulation for the double-diffusive natural convection problem
in porous media: Find
$(\vec{\boldsymbol{\xi}},\boldsymbol{\varphi})\in\mathbb{H}\times\mathbf{H}^{1}(\Omega)$
, with
$\boldsymbol{\varphi}\big{|}_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{\mathrm{D}}$,
such that

$\displaystyle\mathbf{A}_{\boldsymbol{\varphi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})+\mathbf{B}_{\boldsymbol{u}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})$
$\displaystyle=$ $\displaystyle
F_{\boldsymbol{\varphi}}(\vec{\boldsymbol{\eta}})\qquad\quad\forall\,\vec{\boldsymbol{\eta}}\in\mathbb{H},$
(3.6a) $\displaystyle\mathbf{C}(\boldsymbol{\varphi},\boldsymbol{\psi})$
$\displaystyle=$ $\displaystyle
G_{\boldsymbol{u},\boldsymbol{\varphi}}(\boldsymbol{\psi})\qquad\forall\,\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega),$
(3.6b)

where, given $\boldsymbol{\phi}\in\mathbf{H}^{1}(\Omega)$ and
$\boldsymbol{w}\in\mathbf{H}_{0}^{1}(\Omega)$,
$\mathbf{A}_{\boldsymbol{\phi}}$, $\mathbf{B}_{\boldsymbol{w}}$ and
$\mathbf{C}$ are the bilinear forms defined, respectively, as

$\begin{array}[]{rcl}\mathbf{A}_{\boldsymbol{\phi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})&:=&\displaystyle
2\int_{\Omega}\nu(\boldsymbol{\phi})\boldsymbol{t}:\big{\\{}\,\boldsymbol{r}-\kappa_{3}\boldsymbol{\tau}^{\mathrm{d}}\,\big{\\}}\,+\,\int_{\Omega}\boldsymbol{t}:\\{\boldsymbol{\tau}^{\mathrm{d}}-\kappa_{1}\mathbf{e}(\boldsymbol{v})\\}\,-\,\int_{\Omega}\boldsymbol{\sigma}^{\mathrm{d}}:\big{\\{}\boldsymbol{r}-\kappa_{3}\boldsymbol{\tau}^{\mathrm{d}}\big{\\}}\\\\[12.91663pt]
&&\displaystyle+\
(1-\kappa_{2}\gamma)\int_{\Omega}\boldsymbol{u}\cdot{\mathbf{div}}\,\boldsymbol{\tau}\,-\,\int_{\Omega}\boldsymbol{v}\cdot\mathbf{div}\,\boldsymbol{\sigma}\,+\,\int_{\Omega}\boldsymbol{\omega}(\boldsymbol{u}):\boldsymbol{\tau}-\int_{\Omega}\boldsymbol{\sigma}:\boldsymbol{\omega}(\boldsymbol{v})\\\\[12.91663pt]
&&\displaystyle+\
\gamma\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}\,+\,\kappa_{2}\int_{\Omega}\mathbf{div}\,\boldsymbol{\sigma}\cdot\mathbf{div}\,\boldsymbol{\tau}\,+\,\kappa_{1}\int_{\Omega}\mathbf{e}(\boldsymbol{u}):\mathbf{e}(\boldsymbol{v})\,,\end{array}$
(3.7)

and

$\mathbf{B}_{\boldsymbol{w}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})\
:=\
\int_{\Omega}(\boldsymbol{u}\otimes\boldsymbol{w})^{\mathrm{d}}:\big{\\{}\kappa_{3}\boldsymbol{\tau}^{\mathrm{d}}-\boldsymbol{r}\big{\\}}\,,$
(3.8)

for all $\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}\in\mathbb{H}$ and

$\mathbf{C}(\boldsymbol{\chi},\boldsymbol{\psi})\ :\
=\,\int_{\Omega}\mathbb{K}\nabla\boldsymbol{\chi}:\nabla\boldsymbol{\psi}\qquad\forall\,\boldsymbol{\chi},\boldsymbol{\psi}\in\mathbf{H}^{1}(\Omega)\,.$
(3.9)

In turn, $F_{\boldsymbol{\phi}}$ and $G_{\boldsymbol{w},\boldsymbol{\phi}}$
are linear functionals defined as

$F_{\boldsymbol{\phi}}(\vec{\boldsymbol{\eta}})\ :=\
\int_{\Omega}(\boldsymbol{\alpha}\cdot\boldsymbol{\phi}){\boldsymbol{g}}\cdot\\{\boldsymbol{v}-\kappa_{2}\mathbf{div}\,\boldsymbol{\tau}\big{\\}}\,\qquad\forall\,\vec{\boldsymbol{\eta}}\in\mathbb{H}\,,$
(3.10)

and

$G_{\boldsymbol{w},\boldsymbol{\phi}}(\boldsymbol{\psi})\ :=\
-\int_{\Omega}(\nabla\boldsymbol{\phi})\boldsymbol{w}\cdot\boldsymbol{\psi}\,\qquad\forall\,\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)\,.$
(3.11)

###  3.2 The fixed point approach

In this Section we describe a fixed point strategy that allow us to solve the
coupled problem given by (3.6). We start by denoting
$\mathbf{H}:=\mathbf{H}_{0}^{1}(\Omega)\times\mathbf{H}^{1}(\Omega)$, and
defining the operator
$\text{\Large$\mathpzc{A}$}:\mathbf{H}\rightarrow\mathbb{H}$ by

$\text{\Large$\mathpzc{A}$}(\boldsymbol{w},\boldsymbol{\phi})\,=\,(\text{\Large$\mathpzc{A}$}_{1}(\boldsymbol{w},\boldsymbol{\phi}),\,\text{\Large$\mathpzc{A}$}_{2}(\boldsymbol{w},\boldsymbol{\phi}),\,\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi}))\,:=\,\vec{\boldsymbol{\xi}}\quad\quad\forall\,(\boldsymbol{w},\boldsymbol{\phi})\in\mathbf{H}$

where
$\vec{\boldsymbol{\xi}}=(\boldsymbol{t},\boldsymbol{\sigma},\boldsymbol{u})\in\mathbb{H}$
is the unique solution of the augmented mixed formulation given by (3.6a),
with $(\boldsymbol{w},\boldsymbol{\phi})$ instead
$(\boldsymbol{u},\boldsymbol{\varphi})$, that is:

$\begin{array}[]{rlll}\mathbf{A}_{\boldsymbol{\phi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})\,+\,\mathbf{B}_{\boldsymbol{w}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})&=&F_{\boldsymbol{\phi}}(\vec{\boldsymbol{\eta}})&\quad\forall\,\vec{\boldsymbol{\eta}}\in\mathbb{H}\,,\end{array}$
(3.12)

where the bilinear forms $\mathbf{A}_{\boldsymbol{\phi}}$,
$\mathbf{B}_{\boldsymbol{w}}$ and the functional $F_{\boldsymbol{\phi}}$ are
defined exactly as in (3.7), (3.8), and (3.10), respectively. In addition, we
also introduce the operator
$\text{\Large$\mathpzc{B}$}:\mathbf{H}\rightarrow\mathbf{H}^{1}(\Omega)$
defined as

$\text{\Large$\mathpzc{B}$}(\boldsymbol{w},\boldsymbol{\phi})\,:=\,\boldsymbol{\varphi}\quad\quad\forall\,(\boldsymbol{w},\boldsymbol{\phi})\in\mathbf{H},$

where $\boldsymbol{\varphi}$ is the unique solution of the problem (3.6b),
with $(\boldsymbol{w},\boldsymbol{\phi})$ instead
$(\boldsymbol{u},\boldsymbol{\varphi})$, that is:

$\mathbf{C}(\boldsymbol{\varphi},\boldsymbol{\psi})\ =\
G_{\boldsymbol{w},\boldsymbol{\phi}}(\boldsymbol{\psi})\qquad\forall\,\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega),$
(3.13)

where the bilinear form $\mathbf{C}$, and the functional
$G_{\boldsymbol{w},\boldsymbol{\phi}}$ are defined by (3.9) and (3.11),
respectively.

In this way, by introducing the operator
$\text{\large$\mathscr{L}$}:\mathbf{H}\rightarrow\mathbf{H}$ as

$\text{\large$\mathscr{L}$}(\boldsymbol{w},\boldsymbol{\phi})\,:=\,(\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi}),\,\text{\Large$\mathpzc{B}$}(\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi}),\boldsymbol{\phi}))\quad\quad\forall\,(\boldsymbol{w},\boldsymbol{\phi})\in\mathbf{H}$
(3.14)

we readily realize that solving (3.6) is equivalent to seeking a fixed point
of $\mathscr{L}$, that is: Find
$(\boldsymbol{u},\boldsymbol{\varphi})\in\mathbf{H}$ such that

$\text{\large$\mathscr{L}$}(\boldsymbol{u},\boldsymbol{\varphi})\,=\,(\boldsymbol{u},\boldsymbol{\varphi})\,.$
(3.15)

The following two sections establish the well-posedness of (3.15).

###  3.3 Well-posedness of the uncoupled problems

We begin by recalling the following lemmas which are useful to prove the
ellipticity of the bilinear form
$\mathbf{A}_{\boldsymbol{\phi}}+\mathbf{B}_{\boldsymbol{w}}$.

###### Lemma 3.1

There exists $c_{2}(\Omega)>0$ such that

$c_{2}(\Omega)\|\boldsymbol{\tau}_{0}\|^{2}_{0,\Omega}\ \leq\
\|\boldsymbol{\tau}^{\mathrm{d}}\|^{2}_{0,\Omega}\,+\,\|\mathbf{div}(\boldsymbol{\tau})\|^{2}_{0,\Omega}\quad\quad\forall\,\boldsymbol{\tau}=\boldsymbol{\tau}_{0}+c\,\mathbb{I}\in\mathbb{H}(\mathbf{div};\Omega)\,,$

with $\boldsymbol{\tau}_{0}\in\mathbb{H}_{0}(\mathbf{div};\Omega)$ and
$c\in\mathrm{R}$.

Proof. See [9, Proposition 3.1]. $\square$

###### Lemma 3.2

There holds

$\frac{1}{2}|\boldsymbol{v}|_{1,\Omega}^{2}\;\leq\;\|\mathbf{e}(\boldsymbol{v})\|_{0,\Omega}^{2}\quad\quad\forall\,\boldsymbol{v}\in\mathbf{H}^{1}_{0}(\Omega)\,.$

Proof. See [36, Theorem 10.1]. $\square$

We now provide sufficient conditions under which the uncoupled problems (3.12)
and (3.13) are indeed uniquely solvable.

###### Lemma 3.3

Assume that
$\kappa_{1}\in\big{(}0,2\delta_{3}\big{(}2\nu_{1}-\frac{\kappa_{3}\nu_{2}}{\delta_{1}}\big{)}\big{)}$,
$\kappa_{2}\in(0,2\delta_{2})$,
$\kappa_{3}\in\big{(}0,\frac{2\nu_{1}\delta_{1}}{\nu_{2}}\big{)}$ with
$\delta_{1}\in\big{(}0,\frac{1}{\nu_{2}}\big{)}$,
$\delta_{2}\in\big{(}0,\frac{2}{\gamma}\big{)}$, $\delta_{3}\in(0,1)$. Then,
there exists $r_{0}>0$ such that for each $r\in(0,r_{0})$, problem (3.12) has
a unique solution
$\text{\Large$\mathpzc{A}$}(\boldsymbol{w},\boldsymbol{\phi}):=\vec{\boldsymbol{\xi}}\in\mathbb{H}$,
for each $(\boldsymbol{w},\boldsymbol{\phi})\in\mathbf{H}$ with
$\|\boldsymbol{w}\|_{1,\Omega}\leq r$. Moreover, there exists
$C_{\\!\text{\footnotesize$\mathpzc{A}$}}$, independent of
$(\boldsymbol{w},\boldsymbol{\phi})$, such that

$\|\text{\Large$\mathpzc{A}$}(\boldsymbol{w},\boldsymbol{\phi})\|_{\mathbb{H}}\;=\;\|\vec{\boldsymbol{\xi}}\|_{\mathbb{H}}\;\leq\;C_{\\!\text{\footnotesize$\mathpzc{A}$}}\;|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\phi}\|_{0,\Omega}\quad\forall\,(\boldsymbol{w},\boldsymbol{\phi})\in\mathbf{H}\,.$
(3.16)

Proof. We begin by deriving the continuity of the bilinear forms
$\mathbf{A}_{\boldsymbol{\phi}}$ and $\mathbf{B}_{\boldsymbol{w}}$ (sf. (3.7)
and (3.8), respectively). Indeed, applying Cauchy-Schwarz’s inequality, the
assumptions (2.3), and the fact that
$\|\mathbf{e}(\boldsymbol{v})\|_{0,\Omega}\leq|\boldsymbol{v}|_{1,\Omega}$ and
$\|\boldsymbol{\omega}(\boldsymbol{v})\|_{0,\Omega}\leq|\boldsymbol{v}|_{1,\Omega}\;\;\forall\,\boldsymbol{v}\in\mathbf{H}^{1}(\Omega)$,
we deduce that there exists a positive constant
$C_{\mathbf{A}_{\boldsymbol{\phi}}}$, depending on
$\kappa_{1},\kappa_{2},\kappa_{3},\nu_{2},\gamma$, such that

$|\mathbf{A}_{\boldsymbol{\phi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})|\;\leq\;C_{\mathbf{A}_{\boldsymbol{\phi}}}\|\vec{\boldsymbol{\xi}}\|_{\mathbb{H}}\|\vec{\boldsymbol{\eta}}\|_{\mathbb{H}}\qquad\forall\,\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}\in\mathbb{H}\,,$
(3.17)

where
$C_{\mathbf{A}_{\boldsymbol{\phi}}}:=3\max\\{2,\gamma+\kappa_{1},\kappa_{2}+\kappa_{3},2\nu_{2},1+2\nu_{2}\kappa_{3},1+|1-\gamma\kappa_{2}|\\}$.
In turn, by applying Hölder’s inequality and the continuous injection
$\mathbf{i}:\mathbf{H}^{1}(\Omega)\to\mathbf{L}^{4}(\Omega)$, with constant
$c_{1}(\Omega)$, we obtain that

$|\mathbf{B}_{\boldsymbol{w}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})|\;\leq\;c_{1}(\Omega)(2+\kappa_{3}^{2})^{1/2}\|\boldsymbol{w}\|_{1,\Omega}\,\|\vec{\boldsymbol{\xi}}\|_{\mathbb{H}}\,\|\vec{\boldsymbol{\eta}}\|_{\mathbb{H}}\qquad\forall\,\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}\in\mathbb{H}\,.$
(3.18)

Then, it follows from (3.17) and (3.18), that there exists a positive constant
denote by $\|\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{w}}\|$,
that depends on
$\kappa_{1},\kappa_{2},\kappa_{3},\nu_{2},\gamma,c_{1}(\Omega)$ and
$\|\boldsymbol{w}\|_{1,\Omega}$, such that

$|(\mathbf{A}_{\boldsymbol{\phi}}+\mathbf{B}_{\boldsymbol{w}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})|\;\leq\;\|\mathbf{A}_{\boldsymbol{\phi}}+\mathbf{B}_{\boldsymbol{w}}\|\,\|\vec{\boldsymbol{\xi}}\|_{\mathbb{H}}\,\|\vec{\boldsymbol{\eta}}\|_{\mathbb{H}}\qquad\forall\,\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}\in\mathbb{H}\,.$
(3.19)

We now address the proof of the ellipticity for the bilinear form
$\mathbf{A}_{\boldsymbol{\phi}}+\mathbf{B}_{\boldsymbol{w}}$. To this end, we
first derive this property for the bilinear form
$\mathbf{A}_{\boldsymbol{\phi}}$. In fact, from (3.7), by applying (2.3) and
the Cauchy-Schwarz inequality, together with Lemma 3.2, it followss that

$\displaystyle\mathbf{A}_{\boldsymbol{\phi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}})$
$\displaystyle=$ $\displaystyle
2\int_{\Omega}\nu(\boldsymbol{\phi})\boldsymbol{t}:\boldsymbol{t}\,-\,2\kappa_{3}\int_{\Omega}\nu(\boldsymbol{\phi})\boldsymbol{t}:\boldsymbol{\sigma}^{\mathrm{d}}\,-\,\kappa_{1}\int_{\Omega}\boldsymbol{t}:\mathbf{e}(\boldsymbol{u})\,+\,\kappa_{3}\|\boldsymbol{\sigma}^{\mathrm{d}}\|_{0,\Omega}^{2}$
$\displaystyle-\
\kappa_{2}\gamma\int_{\Omega}\boldsymbol{u}\cdot{\mathbf{div}}\,\boldsymbol{\sigma}\,+\,\gamma\|\boldsymbol{u}\|_{0,\Omega}^{2}\,+\,\kappa_{2}\|\mathbf{div}\,\boldsymbol{\sigma}\|_{0,\Omega}^{2}\,+\,\kappa_{1}\|\mathbf{e}(\boldsymbol{u})\|_{0,\Omega}^{2}$
$\displaystyle\geq$ $\displaystyle
2\nu_{1}\|\boldsymbol{t}\|_{0,\Omega}^{2}\,-\,2\kappa_{3}\nu_{2}\|\boldsymbol{t}\|_{0,\Omega}\|\boldsymbol{\sigma}^{\mathrm{d}}\|_{0,\Omega}\,-\,\kappa_{1}\|\boldsymbol{t}\|_{0,\Omega}\|\mathbf{e}(\boldsymbol{u})\|_{0,\Omega}\,+\,\kappa_{3}\|\boldsymbol{\sigma}^{\mathrm{d}}\|_{0,\Omega}^{2}$
$\displaystyle-\
\kappa_{2}\gamma\|\boldsymbol{u}\|_{0,\Omega}\|{\mathbf{div}}\,\boldsymbol{\sigma}\|_{0,\Omega}\,+\,\gamma\|\boldsymbol{u}\|_{0,\Omega}^{2}\,+\,\kappa_{2}\|\mathbf{div}\,\boldsymbol{\sigma}\|_{0,\Omega}^{2}\,+\,\kappa_{1}\|\mathbf{e}(\boldsymbol{u})\|_{0,\Omega}^{2}$
$\displaystyle\geq$ $\displaystyle
2\nu_{1}\|\boldsymbol{t}\|_{0,\Omega}^{2}\,-\,2\kappa_{3}\nu_{2}\|\boldsymbol{t}\|_{0,\Omega}\|\boldsymbol{\sigma}^{\mathrm{d}}\|_{0,\Omega}\,-\,\kappa_{1}\|\boldsymbol{t}\|_{0,\Omega}|\boldsymbol{u}|_{1,\Omega}\,+\,\kappa_{3}\|\boldsymbol{\sigma}^{\mathrm{d}}\|_{0,\Omega}^{2}$
$\displaystyle-\
\kappa_{2}\gamma\|\boldsymbol{u}\|_{0,\Omega}\|{\mathbf{div}}\,\boldsymbol{\sigma}\|_{0,\Omega}\,+\,\gamma\|\boldsymbol{u}\|_{0,\Omega}^{2}\,+\,\kappa_{2}\|\mathbf{div}\,\boldsymbol{\sigma}\|_{0,\Omega}^{2}\,+\,\frac{\kappa_{1}}{2}|\boldsymbol{u}|_{1,\Omega}^{2}$

Next, employing the Young inequality and gathering similar terms, we obtain

$\displaystyle\mathbf{A}_{\boldsymbol{\phi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}})$
$\displaystyle\geq$
$\displaystyle\left(2\nu_{1}-\frac{\kappa_{3}\nu_{2}}{\delta_{1}}-\frac{\kappa_{1}}{2\delta_{3}}\right)\|\boldsymbol{t}\|_{0,\Omega}^{2}\,+\,\kappa_{3}\big{(}1-\nu_{2}\delta_{1}\big{)}\|\boldsymbol{\sigma}^{\mathrm{d}}\|_{0,\Omega}^{2}\,+\,\kappa_{2}\left(1-\frac{\gamma\delta_{2}}{2}\right)\|{\mathbf{div}}\,\boldsymbol{\sigma}\|_{0,\Omega}^{2}$
$\displaystyle+\
\gamma\left(1-\frac{\kappa_{2}}{2\delta_{2}}\right)\|\boldsymbol{u}\|_{0,\Omega}^{2}\,+\,\frac{\kappa_{1}}{2}\big{(}1-\delta_{3}\big{)}|\boldsymbol{u}|_{1,\Omega}^{2}$

and then, choosing $\kappa_{1}$, $\kappa_{2}$, $\kappa_{3}$, $\delta_{1}$,
$\delta_{2}$ and $\delta_{3}$ in the ranges specified of the statement of the
present Lemma, we deduce that there exists a positive constant
$\alpha(\Omega)$, independent of $(\boldsymbol{w},\boldsymbol{\varphi})$, such
that

$\mathbf{A}_{\boldsymbol{\phi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}})\
\geq\
\alpha(\Omega)\,\|\vec{\boldsymbol{\xi}}\|_{\mathbb{H}}^{2}\qquad\forall\,\vec{\boldsymbol{\xi}}\in\mathbb{H}\,.$
(3.20)

Next, by combining (3.18) and (3.20), we find that

$(\mathbf{A}_{\boldsymbol{\phi}}+\mathbf{B}_{\boldsymbol{w}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}})\
\geq\
\Big{\\{}\alpha(\Omega)-c_{1}(\Omega)(2+\kappa_{3}^{2})^{1/2}\|\boldsymbol{w}\|_{1,\Omega}\Big{\\}}\|\vec{\boldsymbol{\xi}}\|_{\mathbb{H}}^{2}\qquad\forall\,\vec{\boldsymbol{\xi}}\in\mathbb{H}\,,$

from which, we deduce that

$(\mathbf{A}_{\boldsymbol{\phi}}+\mathbf{B}_{\boldsymbol{w}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}})\;\geq\;\frac{\alpha(\Omega)}{2}\|\vec{\boldsymbol{\xi}}\|_{\mathbb{H}}^{2}\qquad\forall\,\vec{\boldsymbol{\xi}}\in\mathbb{H}\,,$
(3.21)

provided that $\|\boldsymbol{w}\|_{1,\Omega}\leq r_{0}$, with

$r_{0}\ :=\ \frac{\alpha(\Omega)}{2c_{1}(\Omega)(2+\kappa_{3}^{2})^{1/2}}\,,$
(3.22)

which confirms the ellipticity of
$\mathbf{A}_{\boldsymbol{\phi}}+\mathbf{B}_{\boldsymbol{w}}$. On the other
hand, by applying the Cauchy-Schwarz inequality, we deduce that
$F_{\boldsymbol{\phi}}\in\mathbb{H}^{\prime}$ (cf. (3.10)) with

$\|F_{\boldsymbol{\phi}}\|\ \leq\
\sqrt{2}\,(2+\kappa_{2}^{2})^{1/2}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\phi}\|_{0,\Omega}\,.$
(3.23)

Consequently, a straightforward application of the Lax-Milgram lemma implies
that there exists a unique solution $\vec{\boldsymbol{\xi}}\in\mathbb{H}$ of
(3.12). Finally, from (3.21) and (3.23), and performing simple algebraic
manipulations, we derive (3.16), with $\displaystyle
C_{\\!\text{\footnotesize$\mathpzc{A}$}}:=\frac{2\,\sqrt{2}\,(2+\kappa_{2}^{2})^{1/2}}{\alpha(\Omega)}>0$,
independent of $(\boldsymbol{w},\boldsymbol{\phi})$. $\square$

###### Remark 3.1

At this point, we remark that for computational purposes, the constant
$\alpha(\Omega)$ defined in Lemma 3.3, can be maximized by choosing the
parameters $\delta_{1}$, $\delta_{2}$, $\delta_{3}$, $\kappa_{1}$,
$\kappa_{2}$ and $\kappa_{3}$ at the middle points of their feasible ranges.
Thus, adequate choices for these parameters are
$\delta_{1}:=\frac{1}{2\nu_{2}}$, $\delta_{2}:=\frac{1}{\gamma}$ and
$\delta_{3}:=\frac{1}{2}$, which establish that

$\kappa_{1}\ =\ \frac{\nu_{1}}{2}\,,\qquad\kappa_{2}\ =\
\frac{1}{\gamma}\qquad\text{and}\qquad\kappa_{3}\ =\
\frac{\nu_{1}}{2\nu_{2}^{2}}\,.$ (3.24)

###### Remark 3.2

In order to establish the solvability of the problem (3.13), associated to the
operator $\mathpzc{B}$, we first point out that
$\boldsymbol{\varphi}_{{\mathrm{D}}}\in\mathbf{H}^{1/2}({\Gamma_{\mathrm{D}}})$
can be continuously extended in the trace sense to $\mathbf{H}^{1/2}(\Gamma)$.
Indeed, the set

$\mathcal{H}(\boldsymbol{\varphi}_{{\mathrm{D}}})\ =\
(\gamma_{0}|_{{\Gamma_{\mathrm{D}}}})^{-1}(\\{\boldsymbol{\varphi}_{{\mathrm{D}}}\\})\
=\
\big{\\{}\boldsymbol{\varphi}\in\mathbf{H}^{1}(\Omega):\quad\gamma_{0}(\boldsymbol{\varphi})|_{{\Gamma_{\mathrm{D}}}}\;=\;\boldsymbol{\varphi}_{{\mathrm{D}}}\big{\\}}$

is a closed and convex since the usual trace operator
$\gamma_{0}:\mathbf{H}^{1}(\Omega)\to\mathbf{H}^{1/2}(\Gamma)$ is linear and
continuous. Then, from The Best Approximation Theorem [24, Theorem 7] there
exists a unique
$\widetilde{\boldsymbol{\varphi}}:=E(\boldsymbol{\varphi}_{{\mathrm{D}}})\in\mathbf{H}^{1}(\Omega)$,
where $E$ denotes the extension of $\boldsymbol{\varphi}_{{\mathrm{D}}}$, such
that
$\gamma_{0}(\widetilde{\boldsymbol{\varphi}})|_{{\Gamma_{\mathrm{D}}}}\;=\;\boldsymbol{\varphi}_{\mathrm{D}}$
with

$\|\widetilde{\boldsymbol{\varphi}}\|_{1,\Omega}\ =\
\inf_{\boldsymbol{\varphi}\in\mathcal{H}(\boldsymbol{\varphi}_{{\mathrm{D}}})}\|\boldsymbol{\varphi}\|_{1,\Omega}\
=\ \text{\rm dist}({\bf 0},\mathcal{H}(\boldsymbol{\varphi}_{{\mathrm{D}}}))\
=\ \|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}.$
(3.25)

In this way, the suitable extension of
$\boldsymbol{\varphi}_{{\mathrm{D}}}\in\mathbf{H}^{1/2}({\Gamma_{\mathrm{D}}})$
is not but the corresponding one of
$\gamma_{0}(\widetilde{\boldsymbol{\varphi}})\in\mathbf{H}^{1/2}(\Gamma)$ to
$\mathbf{H}^{1}(\Omega)$.

###### Lemma 3.4

For each $(\boldsymbol{w},\boldsymbol{\phi})\in\mathbf{H}$, problem (3.13) has
a unique solution $\boldsymbol{\varphi}\in\mathbf{H}^{1}(\Omega)$, with
$\boldsymbol{\varphi}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}}}$.
Moreover, there exists a constant $C_{\\!\text{\footnotesize$\mathpzc{B}$}}>0$
independent of $(\boldsymbol{w},\boldsymbol{\phi})$, such that

$\|\text{\Large$\mathpzc{B}$}(\boldsymbol{w},\boldsymbol{\phi})\|_{1,\Omega}\;=\;\|\boldsymbol{\varphi}\|_{1,\Omega}\;\leq\;C_{\\!\text{\footnotesize$\mathpzc{B}$}}\big{\\{}\|\boldsymbol{w}\|_{1,\Omega}\|\boldsymbol{\phi}\|_{1,\Omega}\,+\,\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}\big{\\}}\,.$
(3.26)

Proof. We begin by nothing, according to the aforementioned in Remark 3.2,
that given
$\boldsymbol{\varphi}_{{\mathrm{D}}}\in\mathbf{H}^{1/2}({\Gamma_{\mathrm{D}}})$
there exists $\boldsymbol{\varphi}_{1}\in\mathbf{H}^{1}(\Omega)$ such that

$\boldsymbol{\varphi}_{1}|_{{\Gamma_{\mathrm{D}}}}\ =\
\boldsymbol{\varphi}_{{\mathrm{D}}}\qquad\text{and}\qquad\|\boldsymbol{\varphi}_{1}\|_{1,\Omega}\
=\ \|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}\,.$
(3.27)

Then, we consider the auxiliary linear problem: Find
$\boldsymbol{\varphi}_{0}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)$
such that

$\mathbf{C}(\boldsymbol{\varphi}_{0},\boldsymbol{\psi})\ =\
G_{\boldsymbol{w},\boldsymbol{\phi}}(\boldsymbol{\psi})\,-\,\mathbf{C}(\boldsymbol{\varphi}_{1},\boldsymbol{\psi})\qquad\forall\,\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)$
(3.28)

where $\mathbf{C}$ and $G_{\boldsymbol{w},\boldsymbol{\phi}}$ are defined in
(3.9) and (3.11), respectively. In addition, from (3.9) and the Cauchy-
Schwarz’s inequality, we deduce that

$|\mathbf{C}(\boldsymbol{\chi},\boldsymbol{\psi})|\ \leq\
\|\mathbb{K}\|_{\infty,\Omega}\|\boldsymbol{\chi}\|_{1,\Omega}\|\boldsymbol{\psi}\|_{1,\Omega}\qquad\forall\,\boldsymbol{\chi},\boldsymbol{\psi}\in\mathbf{H}^{1}(\Omega)\,,$
(3.29)

which, in particular, confirms the boundedness of $\mathbf{C}$. Then, from
(3.9), using (2.2) and the Poncairé inequality (3.1), we get

$\mathbf{C}(\boldsymbol{\psi},\boldsymbol{\psi})\ \geq\
\widetilde{\alpha}\,\|\boldsymbol{\psi}\|_{1,\Omega}^{2}\qquad\forall\,\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)\,,$
(3.30)

which proves that $\mathbf{C}$ is
$\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)$-elliptic with constant

$\widetilde{\alpha}\ :=\ k_{0}\,c_{\text{gp}}^{-1}\,.$ (3.31)

Next, applying the Cauchy-Schwarz’s inequality, the boundedness of the
continuous injection
$\boldsymbol{i}:\mathbf{H}^{1}(\Omega)\to\mathbf{L}^{4}(\Omega)$, with
constant $c_{1}(\Omega):=\|\boldsymbol{i}\|^{2}$, (3.29), and the second
identity in (3.27), we easily deduce that

$|G_{\boldsymbol{w},\boldsymbol{\phi}}(\boldsymbol{\psi})-\mathbf{C}(\boldsymbol{\varphi}_{1},\boldsymbol{\psi})|\
\leq\
\Big{\\{}c_{1}(\Omega)\|\boldsymbol{w}\|_{1,\Omega}\|\boldsymbol{\phi}\|_{1,\Omega}+\|\mathbb{K}\|_{\infty,\Omega}\|\boldsymbol{\varphi}_{\mathrm{D}}\|_{1/2,{\Gamma_{\mathrm{D}}}}\Big{\\}}\|\boldsymbol{\psi}\|_{1,\Omega}\quad\forall\,\boldsymbol{\psi}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)\,,$
(3.32)

which establishes the boundedness of the right-hand side of (3.28).
Consequently, a direct application of the Lax-Milgram lemma implies that there
exists a unique
$\boldsymbol{\varphi}_{0}\in\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)$
that satisfies (3.28) with

$\|\boldsymbol{\varphi}_{0}\|_{1,\Omega}\ \leq\
\frac{1}{\widetilde{\alpha}}\,\big{\\{}c_{1}(\Omega)\|\boldsymbol{w}\|_{1,\Omega}\|\boldsymbol{\phi}\|_{1,\Omega}\,+\,\|\mathbb{K}\|_{\infty,\Omega}\|\boldsymbol{\varphi}_{\mathrm{D}}\|_{1/2,{\Gamma_{\mathrm{D}}}}\big{\\}}\,.$

On the other hand, we now set
$\boldsymbol{\varphi}:=\boldsymbol{\varphi}_{0}+\boldsymbol{\varphi}_{1}$,
which satisfies that
$\boldsymbol{\varphi}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{0}|_{{\Gamma_{\mathrm{D}}}}+\boldsymbol{\varphi}_{1}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}}}$.
In addition, it is easy to check that $\boldsymbol{\varphi}$ is a solution of
problem (3.13), where the uniqueness comes from (3.30). Finally,
$\boldsymbol{\varphi}$ verifies the estimate (3.26), with constant
$C_{\\!\text{\footnotesize$\mathpzc{B}$}}:=\max\\{\widetilde{\alpha}^{-1}c_{1}(\Omega),\widetilde{\alpha}^{-1}\|\mathbb{K}\|_{\infty,\Omega}+1\\}$.
$\square$

We end this section by introducing a further regularity hypotheses on the
problem definig $\mathpzc{A}$, which will be employed to derive a Lipschitz-
continuity property for the operator $\mathscr{L}$. More precisely, we assume
that for each $(\boldsymbol{w},\boldsymbol{\phi})\in\mathbf{H}$ with
$\|\boldsymbol{w}\|_{1,\Omega}+\|\boldsymbol{\phi}\|_{1,\Omega}\leq r$, $r>0$
given, there holds
$(\boldsymbol{r},\boldsymbol{\tau},\boldsymbol{v})=\text{\Large$\mathpzc{A}$}(\boldsymbol{w},\boldsymbol{\phi})\in\mathbb{L}^{2}_{\mathrm{tr}}(\Omega)\cap\mathbb{H}^{\varepsilon}(\Omega)\times\mathbb{H}_{0}({\mathbf{div}};\Omega)\cap\mathbb{H}^{\varepsilon}(\Omega)\times\mathbf{H}^{1+\varepsilon}(\Omega)$,
for some $\varepsilon\in(0,1)$ (when $n=2$) or $\varepsilon\in[\frac{1}{2},1)$
(when $n=3$), with

$\|\boldsymbol{r}\|_{\varepsilon,\Omega}\,+\,\|\boldsymbol{\tau}\|_{\varepsilon,\Omega}\,+\,\|\boldsymbol{v}\|_{\varepsilon,\Omega}\
\leq\
\widehat{C}(r)\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\phi}\|_{0,\Omega}\,,$
(3.33)

where $\widehat{C}(r)$ is a positive constant independent of
$(\boldsymbol{w},\boldsymbol{\phi})$ but depending on the upper bound $r$ of
its norm. The reason of the estipulated ranges for $\varepsilon$ will be
clarified in the forthcomming analysis (specifically in the proofs of Lemmas
3.6 and 3.8 below). Also, we pay attention to the fact the while the estimate
(3.33) will be employed only to bound
$\|\boldsymbol{r}\|_{\varepsilon,\Omega}$, we have stated it including the
terms $\|\boldsymbol{\tau}\|_{\varepsilon,\Omega}$ and
$\|\boldsymbol{v}\|_{1+\varepsilon,\Omega}$ as well, since due to the first
and fourth equations of (2.9), the regularities of
$\boldsymbol{r},\boldsymbol{\tau}$ and $\boldsymbol{v}$ will most likely be
connected.

###  3.4 Solvability analysis of the fixed point equation

We begin by emphasizing that the well-posedness of the uncoupled problems
(3.12) and (3.13) confirms that the operators $\mathpzc{A}$, $\mathpzc{B}$,
and $\mathscr{L}$ (cf. (3.14)) are well defined, and hence now we can address
the solvability analysis of the fixed point problem presented in (3.15). To
this end, we will verify below the hypotheses of the Banach fixed-point
theorem.

###### Lemma 3.5

Given $r\in(0,r_{0})$, with $r_{0}$ given by (3.22), we let

$W_{r}\ :=\
\big{\\{}(\boldsymbol{w},\boldsymbol{\phi})\in\mathbf{H}:\quad\|(\boldsymbol{w},\boldsymbol{\phi})\|_{\mathbf{H}}\,\leq\,r\big{\\}}\,,$

and assume that the data satisfy

$c(r)\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,+\,C_{\\!\text{\footnotesize$\mathpzc{B}$}}\,\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}\
\leq\ r$ (3.34)

where
$c(r):=rC_{\\!\text{\footnotesize$\mathpzc{A}$}}(1+rC_{\\!\text{\footnotesize$\mathpzc{B}$}})$,
with $C_{\\!\text{\footnotesize$\mathpzc{A}$}}$ and
$C_{\\!\text{\footnotesize$\mathpzc{B}$}}$ given by (3.16) and (3.26),
respectively. Then, there holds $\text{\large$\mathscr{L}$}(W_{r})\subseteq
W_{r}$.

Proof. It follows similar as in [20, Lemma 3.5] $\square$

###### Lemma 3.6

Let $r\in(0,r_{0})$ with $r_{0}$ given by (3.22). Then, there exists a
constant $L_{\\!\text{\footnotesize$\mathpzc{A}$}}>0$, depending on the
stabilization parameters $\kappa_{2},\kappa_{3}$, and the constants
$L_{\nu},c_{1}(\Omega),\alpha(\Omega),C_{\varepsilon}$ (cf. (2.3), (3.3),
(3.20) and (3.41), respectively), such that for all
$(\boldsymbol{w},\boldsymbol{\phi}),(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\in\mathbf{H}$,
with
$\|\boldsymbol{w}\|_{1,\Omega},\|\widetilde{\boldsymbol{w}}\|_{1,\Omega}\,\leq\;r$,
there holds

$\begin{array}[]{l}\|\text{\Large$\mathpzc{A}$}(\boldsymbol{w},\boldsymbol{\phi})-\text{\Large$\mathpzc{A}$}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{\mathbb{H}}\
\leq\
L_{\\!\text{\footnotesize$\mathpzc{A}$}}\Big{\\{}\|\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi})\|_{1,\Omega}\,\|\boldsymbol{w}-\widetilde{\boldsymbol{w}}\|_{1,\Omega}\\\\[4.30554pt]
\qquad\qquad\qquad\qquad+\
\|\text{\Large$\mathpzc{A}$}_{1}(\boldsymbol{w},\boldsymbol{\phi})\|_{\varepsilon,\Omega}\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{0,n/\varepsilon,\Omega}\
+\
|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}\Big{\\}}\,.\end{array}$
(3.35)

Proof. Given
$(\boldsymbol{w},\boldsymbol{\phi}),(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})$
as stated, we let
$\vec{\boldsymbol{\xi}}=(\boldsymbol{t},\boldsymbol{\sigma},\boldsymbol{u})=\text{\Large$\mathpzc{A}$}(\boldsymbol{w},\boldsymbol{\phi})$
and
$\vec{\boldsymbol{\zeta}}=(\widetilde{\boldsymbol{t}},\widetilde{\boldsymbol{\sigma}},\widetilde{\boldsymbol{u}})=\text{\Large$\mathpzc{A}$}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})$,
which according to the definition of operator $\mathpzc{A}$ (cf. (3.12)),
means that:

$\mathbf{A}_{\boldsymbol{\phi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})+\mathbf{B}_{\boldsymbol{w}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})\
=\
F_{\boldsymbol{\phi}}(\vec{\boldsymbol{\eta}})\qquad\text{and}\qquad\mathbf{A}_{\widetilde{\boldsymbol{\phi}}}(\vec{\boldsymbol{\zeta}},\vec{\boldsymbol{\eta}})+\mathbf{B}_{\widetilde{\boldsymbol{w}}}(\vec{\boldsymbol{\zeta}},\vec{\boldsymbol{\eta}})\
=\
F_{\widetilde{\boldsymbol{\phi}}}(\vec{\boldsymbol{\eta}})\qquad\forall\,\vec{\boldsymbol{\eta}}\in\mathbb{H}\,.$

Then, subtracting both identities, replacing
$\vec{\boldsymbol{\zeta}}=\vec{\boldsymbol{\zeta}}-\vec{\boldsymbol{\xi}}+\vec{\boldsymbol{\xi}}$,
and using the bilinearity of
$\mathbf{A}_{\boldsymbol{\phi}}+\mathbf{B}_{\boldsymbol{w}}$ for any
$\boldsymbol{\phi}$ and $\boldsymbol{w}$, it follows from (3.12) that:

$\displaystyle(\mathbf{A}_{\widetilde{\boldsymbol{\phi}}}+\mathbf{B}_{\widetilde{\boldsymbol{w}}})(\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}},\vec{\boldsymbol{\eta}})\
=\
(F_{\boldsymbol{\phi}}-F_{\widetilde{\boldsymbol{\phi}}})(\vec{\boldsymbol{\eta}})\,+\,(\mathbf{A}_{\widetilde{\boldsymbol{\phi}}}-\mathbf{A}_{\boldsymbol{\phi}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})\,+\,\mathbf{B}_{\widetilde{\boldsymbol{w}}-\boldsymbol{w}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})\qquad\forall\,\vec{\boldsymbol{\eta}}\in\mathbb{H}\,.$
(3.36)

Moreover, applying the ellipticity of
$\mathbf{A}_{\widetilde{\boldsymbol{\phi}}}+\mathbf{B}_{\widetilde{\boldsymbol{w}}}$
(cf. (3.21)), and then employing (3.36) with
$\vec{\boldsymbol{\eta}}:=\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}$, we
find that

$\displaystyle\frac{\alpha(\Omega)}{2}\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}\|_{\mathbb{H}}^{2}$
$\displaystyle\leq$
$\displaystyle(\mathbf{A}_{\widetilde{\boldsymbol{\phi}}}+\mathbf{B}_{\widetilde{\boldsymbol{w}}})(\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}},\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}})$
(3.37) $\displaystyle=$
$\displaystyle(F_{\boldsymbol{\phi}}-F_{\widetilde{\boldsymbol{\phi}}})(\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}})\,+\,(\mathbf{A}_{\widetilde{\boldsymbol{\phi}}}-\mathbf{A}_{\boldsymbol{\phi}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}})\,+\,\mathbf{B}_{\widetilde{\boldsymbol{w}}-\boldsymbol{w}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}})\,.$

Then, for the first and third terms on the right-hand side of (3.37), we
employ the Cauchy-Schwarz and Hölder inequalities, together with the
continuous injection
$\mathbf{i}:\mathbf{H}^{1}(\Omega)\to\mathbf{L}^{4}(\Omega)$, similar as in
(3.23) and (3.18), in order to obtain that

$\begin{array}[]{rcl}|(F_{\boldsymbol{\phi}}-F_{\widetilde{\boldsymbol{\phi}}})(\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}})|&=&\displaystyle\Bigg{|}\int_{\Omega}(\boldsymbol{\alpha}(\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}})){\boldsymbol{g}}\cdot\big{\\{}(\boldsymbol{u}-\widetilde{\boldsymbol{u}})-\kappa_{2}\,{\mathbf{div}}(\boldsymbol{\sigma}-\widetilde{\boldsymbol{\sigma}})\big{\\}}\Bigg{|}\\\\[12.91663pt]
&\leq&\sqrt{2}\,(2+\kappa_{2}^{2})^{1/2}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}\,\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}\|_{\mathbb{H}}\end{array}$
(3.38)

and

$\begin{array}[]{rcl}|\mathbf{B}_{\boldsymbol{w}-\widetilde{\boldsymbol{w}}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}})|&=&\displaystyle\Bigg{|}\int_{\Omega}(\boldsymbol{u}\otimes(\widetilde{\boldsymbol{w}}-\boldsymbol{w}))^{\tt{d}}:\big{\\{}\kappa_{3}(\boldsymbol{\sigma}-\widetilde{\boldsymbol{\sigma}})^{\tt{d}}-(\boldsymbol{t}-\widetilde{\boldsymbol{t}})\big{\\}}\Bigg{|}\\\\[12.91663pt]
&\leq&c_{1}(\Omega)\,(2+\kappa_{3}^{2})^{1/2}\,\|\boldsymbol{u}\|_{1,\Omega}\,\|\boldsymbol{w}-\widetilde{\boldsymbol{w}}\|_{1,\Omega}\,\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}\|_{\mathbb{H}}\end{array}$
(3.39)

On the other hand, for the second term in the right-hand side of (3.37), we
apply the Lipschitz continuity property for $\nu$ given in (2.3), the Cauchy-
Schwarz and Hölder inequalities, and the definition of
$\mathbf{A}_{\boldsymbol{\varphi}}$ (cf. (3.7)), to obtain that

$\begin{array}[]{l}\displaystyle|(\mathbf{A}_{\widetilde{\boldsymbol{\phi}}}-\mathbf{A}_{\boldsymbol{\phi}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}})|\;=\;\Bigg{|}\int_{\Omega}2(\nu(\boldsymbol{\phi})-\nu(\widetilde{\boldsymbol{\phi}}))\boldsymbol{t}:\\{(\boldsymbol{t}-\widetilde{\boldsymbol{t}})-\kappa_{3}(\boldsymbol{\sigma}-\widetilde{\boldsymbol{\sigma}})^{\tt{d}}\\}\Bigg{|}\\\\[17.22217pt]
\hskip
28.45274pt\quad\quad\;\leq\;2L_{\nu}(2+\kappa_{3}^{2})^{1/2}\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{0,2q,\Omega}\|\boldsymbol{t}\|_{0,2p,\Omega}\;\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}\|_{\mathbb{H}},\end{array}$
(3.40)

with $p,q\in(1,+\infty)$ such that $\displaystyle\frac{1}{p}+\frac{1}{q}=1.$
We now proceed as in [7, Lemma 3.9]. In fact, given the further
$\varepsilon$-regularity assumed in (3.33), we reall that the Sobolev
embedding theorem (see e.g [2, Theorem 4.12]) establishes the continuous
injection
$i_{\varepsilon}:\mathrm{H}^{\varepsilon}(\Omega)\to\mathrm{L}^{2p}(\Omega)$
with boundedness constant $C_{\varepsilon}$, where

$2p\ =\
\left\\{\begin{array}[]{cl}\displaystyle\frac{2}{1-\varepsilon}&\text{if
}n=2\\\\[8.61108pt] \displaystyle\frac{6}{3-2\varepsilon}&\text{if
}n=3\end{array}\right.$

and $2q=\frac{n}{\varepsilon}$, and therefore, there holds

$\|\boldsymbol{t}\|_{0,2p,\Omega}\;\leq\;C_{\varepsilon}\,\|\boldsymbol{t}\|_{\varepsilon,\Omega}\qquad\forall\,\boldsymbol{t}\in\mathbb{H}^{\varepsilon}(\Omega)\,.$
(3.41)

In this way, denoting

$L_{\\!\text{\footnotesize$\mathpzc{A}$}}\ :=\
\frac{2}{\alpha(\Omega)}\max\Big{\\{}c_{1}(\Omega)\,(2+\kappa_{3}^{2})^{1/2},\
2\,L_{\nu}(2+\kappa_{3}^{2})^{1/2}\,C_{\varepsilon},\
\sqrt{2}\,(2+\kappa_{2}^{2})^{1/2}\Big{\\}}$

from the inequalities (3.37), (3.38), (3.39), (3.40), and recalling that
$\boldsymbol{t}=\text{\Large$\mathpzc{A}$}_{1}(\boldsymbol{w},\boldsymbol{\phi})$
and
$\boldsymbol{u}=\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi})$,
yields (3.35) and concludes the proof. $\square$

###### Lemma 3.7

There exists $L_{\\!\text{\footnotesize$\mathpzc{B}$}}>0$, depending on
$\widetilde{\alpha}$ and $c_{1}(\Omega)$ (cf. (3.30) and (3.32), respectively)
such that for all
$(\boldsymbol{w},\boldsymbol{\phi}),(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\in
W_{r}$ (cf. Lemma 3.5)there holds

$\|\text{\Large$\mathpzc{B}$}(\boldsymbol{w},\boldsymbol{\phi})-\text{\Large$\mathpzc{B}$}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{1,\Omega}\
\leq\
L_{\\!\text{\footnotesize$\mathpzc{B}$}}\Big{\\{}\|\boldsymbol{w}-\widetilde{\boldsymbol{w}}\|_{1,\Omega}\,+\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}\Big{\\}}\,.$
(3.42)

Proof. Given
$(\boldsymbol{w},\boldsymbol{\phi}),(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\in
W_{r}$, we let
$\boldsymbol{\varphi},\widetilde{\boldsymbol{\varphi}}\in\mathbf{H}^{1}(\Omega)$
be the corresponding solutions of (3.13), that is
$\boldsymbol{\varphi}:=\text{\Large$\mathpzc{B}$}(\boldsymbol{w},\boldsymbol{\phi})$
and
$\widetilde{\boldsymbol{\varphi}}:=\text{\Large$\mathpzc{B}$}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})$.
Then, since
$\boldsymbol{\varphi}|_{{\Gamma_{\mathrm{D}}}}=\widetilde{\boldsymbol{\varphi}}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}}}$,
we realize that $\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}}$
belongs to $\mathbf{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)$. In this way,
applying the ellipticity of $\mathbf{C}$ (cf. (3.30)), using (3.13) and
(3.11), adding and subtracting $\widetilde{\boldsymbol{\phi}}$, and employing
the Hölder inequality, the continuous injection
$\mathbf{i}:\mathbf{H}^{1}(\Omega)\to\mathbf{L}^{4}(\Omega)$, and the
definition of $W_{r}$ (cf. Lemma 3.5), we readily deduce that

$\displaystyle\widetilde{\alpha}\;\|\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}}\|_{1,\Omega}^{2}$
$\displaystyle\leq$
$\displaystyle\mathbf{C}(\boldsymbol{\varphi},\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}})-\mathbf{C}(\widetilde{\boldsymbol{\varphi}},\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}})$
$\displaystyle=$ $\displaystyle
G_{\boldsymbol{w},\boldsymbol{\phi}}(\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}})-G_{\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}}}(\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}})$
$\displaystyle=$
$\displaystyle-\int_{\Omega}\Big{\\{}\big{(}\nabla(\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}})\big{)}\boldsymbol{w}+(\nabla\widetilde{\boldsymbol{\phi}})(\boldsymbol{w}-\widetilde{\boldsymbol{w}})\Big{\\}}\cdot(\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}})$
$\displaystyle\leq$ $\displaystyle
c_{1}(\Omega)\Big{\\{}\|\boldsymbol{w}\|_{1,\Omega}|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}|_{1,\Omega}+|\widetilde{\boldsymbol{\phi}}|_{1,\Omega}\|\boldsymbol{w}-\widetilde{\boldsymbol{w}}\|_{1,\Omega}\Big{\\}}\|\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}}\|_{1,\Omega}$
$\displaystyle\leq$ $\displaystyle
rc_{1}(\Omega)\,\Big{\\{}\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}+\|\boldsymbol{w}-\widetilde{\boldsymbol{w}}\|_{1,\Omega}\Big{\\}}\|\boldsymbol{\varphi}-\widetilde{\boldsymbol{\varphi}}\|_{1,\Omega}$

which give (3.42) with
$L_{\\!\text{\footnotesize$\mathpzc{B}$}}:=\frac{rc_{1}(\Omega)}{\widetilde{\alpha}}$.
$\square$

###### Lemma 3.8

Let $r$ and $W_{r}$ as in Lemma 3.5. Then, there exists a positive constant
$L_{\,\tt{Lip}}$, depending on $r$, $|\boldsymbol{\alpha}|$,
$|{\boldsymbol{g}}|$, and the constants
$C_{\\!\text{\footnotesize$\mathpzc{A}$}}$,
$L_{\\!\text{\footnotesize$\mathpzc{A}$}}$, and
$L_{\\!\text{\footnotesize$\mathpzc{B}$}}$ ((3.16), (3.35), and (3.42),
respectively), such that for all
$(\boldsymbol{w},\boldsymbol{\phi}),(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\in
W_{r}$ there holds

$\begin{array}[]{l}\|\text{\large$\mathscr{L}$}(\boldsymbol{w},\boldsymbol{\phi})-\text{\large$\mathscr{L}$}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{\mathbf{H}}\
\leq\
L_{\,\tt{Lip}}\,\|(\boldsymbol{w},\boldsymbol{\phi})-(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{\mathbf{H}}\,.\end{array}$
(3.44)

Proof. Given $r\in(0,r_{0})$ and
$(\boldsymbol{w},\boldsymbol{\phi}),(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\in
W_{r}$, we first observe, according to the definition of $\mathscr{L}$ (cf.
(3.14)), and the Lipschitz-continuity of $\mathpzc{B}$ (cf. (3.42)), that

$\displaystyle\|\text{\large$\mathscr{L}$}(\boldsymbol{w},\boldsymbol{\phi})-\text{\large$\mathscr{L}$}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{\mathbf{H}}$
$\displaystyle=$
$\displaystyle\|\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi})-\text{\Large$\mathpzc{A}$}_{3}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{1,\Omega}\,+\,\|\text{\Large$\mathpzc{B}$}(\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi}),\boldsymbol{\phi})-\text{\Large$\mathpzc{B}$}(\text{\Large$\mathpzc{A}$}_{3}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}}),\widetilde{\boldsymbol{\phi}})\|_{1,\Omega}$
$\displaystyle\leq$
$\displaystyle(1+L_{\\!\text{\footnotesize$\mathpzc{B}$}})\,\|\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi})-\text{\Large$\mathpzc{A}$}_{3}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{1,\Omega}\,+\,L_{\\!\text{\footnotesize$\mathpzc{B}$}}\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}$

from which, employing the Lipschitz-continuity of $\mathpzc{A}$ (cf. (3.6)),
yields

$\begin{array}[]{l}\|\text{\large$\mathscr{L}$}(\boldsymbol{w},\boldsymbol{\phi})-\text{\large$\mathscr{L}$}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{\mathbf{H}}\
\leq\
(1+L_{\\!\text{\footnotesize$\mathpzc{B}$}})L_{\\!\text{\footnotesize$\mathpzc{A}$}}\Big{\\{}\|\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi})\|_{1,\Omega}\,\|\boldsymbol{w}-\widetilde{\boldsymbol{w}}\|_{1,\Omega}\\\\[8.61108pt]
\qquad\qquad+\
\|\text{\Large$\mathpzc{A}$}_{1}(\boldsymbol{w},\boldsymbol{\phi})\|_{\varepsilon,\Omega}\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{0,n/\varepsilon,\Omega}\,+\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}\Big{\\}}\,+\,L_{\\!\text{\footnotesize$\mathpzc{B}$}}\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}\end{array}$
(3.45)

Then, applying the bound (3.16) to estimate the term
$\|\text{\Large$\mathpzc{A}$}_{3}(\boldsymbol{w},\boldsymbol{\phi})\|_{1,\Omega}$,
employing the continuous injection of $\mathbf{H}^{1}(\Omega)$ into
$\mathbf{L}^{n/\varepsilon}(\Omega)$ with boundedness constant
$\widetilde{C}_{\varepsilon}$, using the estimate (3.33) to estimate the term
$\|\text{\Large$\mathpzc{A}$}_{1}(\boldsymbol{w},\boldsymbol{\phi})\|_{\varepsilon,\Omega}$,
noting that $\|\boldsymbol{\phi}\|_{1,\Omega}\leq r$, and performing some
algebraic manipulations, we get from (3.45) that

$\displaystyle\|\text{\large$\mathscr{L}$}(\boldsymbol{w},\boldsymbol{\phi})-\text{\large$\mathscr{L}$}(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{\mathbf{H}}$
$\displaystyle\leq$
$\displaystyle(1+L_{\\!\text{\footnotesize$\mathpzc{B}$}})L_{\\!\text{\footnotesize$\mathpzc{A}$}}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\Big{\\{}rC_{\\!\text{\footnotesize$\mathpzc{A}$}}\|\boldsymbol{w}-\widetilde{\boldsymbol{w}}\|_{1,\Omega}\,+\,(r\widetilde{C}_{\varepsilon}\widehat{C}(r)+1)\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}\Big{\\}}$
$\displaystyle+\,L_{\\!\text{\footnotesize$\mathpzc{B}$}}\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{1,\Omega}$
$\displaystyle\leq$ $\displaystyle
L_{\,\tt{Lip}}\|(\boldsymbol{w},\boldsymbol{\phi})-(\widetilde{\boldsymbol{w}},\widetilde{\boldsymbol{\phi}})\|_{\mathbf{H}}$

In this way, (3.44) follows from the foregoing inequality by defining

$L_{\,\tt{Lip}}\ :=\
rC_{\\!\text{\footnotesize$\mathpzc{A}$}}(1+L_{\\!\text{\footnotesize$\mathpzc{B}$}})L_{\\!\text{\footnotesize$\mathpzc{A}$}}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,+\,(1+L_{\\!\text{\footnotesize$\mathpzc{B}$}})L_{\\!\text{\footnotesize$\mathpzc{A}$}}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|(r\widetilde{C}_{\varepsilon}\widehat{C}(r)+1)\,+\,L_{\\!\text{\footnotesize$\mathpzc{B}$}}$
(3.46)

in order to complete the proof. $\square$

###### Theorem 3.9

Suppose that the parameters $\kappa_{1},\kappa_{2}$ and $\kappa_{3}$ satisfy
the conditions required by Lemma 3.3. Let $r$ and $W_{r}$ as in Lemma 3.5, and
assume that the data satisfy (3.34) and allow to have

$L_{\,\tt{Lip}}\ <\ 1\,,$ (3.47)

with $L_{\,\tt{Lip}}$ defined in (3.46). Then, problem (3.6) has a unique
solution
$(\vec{\boldsymbol{\xi}},\boldsymbol{\varphi})\in\mathbb{H}\times\mathbf{H}^{1}(\Omega)$
such that
$\boldsymbol{\varphi}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}}}$,
with $(\boldsymbol{u},\boldsymbol{\varphi})\in W_{r}$ (cf. Lemma 3.5), and
there holds

$\|\vec{\boldsymbol{\xi}}\|_{\mathbb{H}}\ \leq\
C_{\\!\text{\footnotesize$\mathpzc{A}$}}\,r\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,,$
(3.48)

and

$\|\boldsymbol{\varphi}\|_{1,\Omega}\ \leq\
C_{\\!\text{\footnotesize$\mathpzc{B}$}}\,\big{\\{}r\,\|\boldsymbol{u}\|_{1,\Omega}\,+\,\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}\big{\\}}\,.$

Proof. It follows as a combination of Lemmas 3.5 and 3.8, the assumption
(3.47), the Banach fixed-point theorem, and the a priori estimates (3.16) and
(3.26). We omit further details. $\square$

##  4 Galerkin scheme

In this section we introduce and analyze the Galerkin scheme of the semi-
augmented mixed-primal problem (3.6). To this end, we now let
$\mathcal{T}_{h}$ be a regular triangulation of $\Omega$ by triangles $K$ (in
$\mathrm{R}^{2}$) or tetrahedra $K$ ( in $\mathrm{R}^{3}$) of diameter
$h_{K}$, and define the meshsize $h:=\max\\{h_{K}:\;K\in\mathcal{T}_{h}\\}$.
In addition, given an integer $k\geq 0$, for each $K\in\mathcal{T}_{h}$ we let
$\mathrm{P}_{k}(K)$ be the space of polynomial functions on $K$ of degree
$\leq\,k$, and define the corresponding local Raviart-Thomas space of order
$k$ as

$\mathbf{RT}_{k}(K)\,:=\,\mathbf{P}_{k}(K)\,\oplus\,\mathrm{P}_{k}(K)\,{\boldsymbol{x}}\,,$

where, according to the notations described in Section 1.1,
$\mathbf{P}_{k}(K)\,=\,[\mathrm{P}_{k}(K)]^{n}$, and ${\boldsymbol{x}}$ is the
generic vector in $\mathrm{R}^{n}$. Then, we consider piecewise polynomials of
degree $\leq k$ for approximating entries of the strain rate $\boldsymbol{t}$,
the global Raviart-Thomas space of order $k$ to approximate rows of the
pseudostress $\boldsymbol{\sigma}$, and the Lagrange space given by the
continuous piecewise polynomial vectors of degree $\leq k+1$ for the velocity
$\boldsymbol{u}$, respectively, that is

$\displaystyle\mathbb{H}^{\boldsymbol{t}}_{h}$ $\displaystyle:=$
$\displaystyle\Big{\\{}\boldsymbol{r}_{h}\in\mathbb{L}^{2}_{\tt{tr}}(\Omega):\quad\boldsymbol{r}_{h}|_{K}\in\mathbb{P}_{k}(K)\quad\forall\,K\in\mathcal{T}_{h}\Big{\\}}\,,$
(4.1) $\displaystyle\mathbf{H}^{\boldsymbol{\sigma}}_{h}$ $\displaystyle:=$
$\displaystyle\Big{\\{}\boldsymbol{\tau}_{h}\in\mathbb{H}_{0}({\mathbf{div}};\Omega):\quad\mathbf{c}^{\tt
t}\,\boldsymbol{\tau}_{h}|_{K}\in\mathbf{RT}_{k}(K)\quad\forall\,\mathbf{c}\in\mathrm{R}^{n}\,,\quad\forall\,K\in\mathcal{T}_{h}\Big{\\}}\,,$
(4.2) $\displaystyle\mathbf{H}^{\boldsymbol{u}}_{h}$ $\displaystyle:=$
$\displaystyle\Big{\\{}\boldsymbol{v}_{h}\in\mathbf{C}(\overline{\Omega}):\quad\boldsymbol{v}_{h}|_{K}\in\mathbf{P}_{k+1}(K)\quad\forall\,K\in\mathcal{T}_{h},\quad\boldsymbol{v}_{h}={\bf
0}\quad\mbox{on}\quad\Gamma\Big{\\}}\,.$ (4.3)

For the unknown $\boldsymbol{\varphi}$ containing the tempeture and the
concentration into its coordinates, we let
$\mathbf{H}_{h}^{\boldsymbol{\varphi}}\subset\mathbf{H}^{1}(\Omega)$ denote
the Lagrange space of degree $\leq k+1$ with respect to $\mathcal{T}_{h}$
(similar to $\mathbf{H}^{\boldsymbol{u}}_{h}$), and set

$\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}\;:=\;\Big{\\{}\boldsymbol{\psi}_{h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h}:\quad\boldsymbol{\psi}_{h}|_{{\Gamma_{\mathrm{D}}}}=\mathbf{0}\Big{\\}}$
(4.4)

to be the analogous space with homogeneous Dirichlet boundary conditions. We
define

${\color[rgb]{0,0,0}\boldsymbol{\varphi}_{{\mathrm{D}},h}\ :=\
\mathcal{I}_{h}^{SZ}(E(\boldsymbol{\varphi}_{{\mathrm{D}}}))\big{|}_{{\Gamma_{\mathrm{D}}}}}$
(4.5)

to be the approximate Dirichlet boundary data, where
$\mathcal{I}_{h}^{SZ}:\mathbf{H}^{1}(\Omega)\to\mathbf{H}^{\boldsymbol{\varphi}}_{h}$
denotes the Scott-Zhang interpolant operator of degree $k+1$, which satisfies
the following stability and approximation properties, respectively, (see [25,
Lemma 1.130]).

###### Lemma 4.1

Let $p$ and $\ell$ satisfy $1\leq p<+\infty$ and $\ell\geq 1$ if $p=1$, and
$1/p<\ell$ otherwise. Then, there exists a positive constant $c$, independent
of $h$, such that the following properties hold:

  * $(i)$

For all $0\leq m\leq\mbox{min}\\{1,\ell\\}$,

$\forall h,\;\;\forall v\in
W^{\ell,p}(\Omega),\quad\quad\|\mathcal{I}_{h}^{SZ}(v)\|_{\,m,p,\Omega}\leq\
c\,\|v\|_{\,\ell,p,\Omega}.$ (4.6)

  * $(ii)$

Provided $\ell\leq k+1$, for all $0\leq m\leq\ell$,

$\forall h,\;\;\forall K\in\mathcal{T}_{h},\;\forall v\in
W^{\ell,p}(\Delta_{K}),\quad\quad\|v-\mathcal{I}_{h}^{SZ}(v)\|_{\,m,p,K}\leq\
c\,h^{\ell-m}_{K}\,|v|_{\,\ell,p,\Delta_{K}}$ (4.7)

where $\Delta_{K}$ denotes the set of elements in $\mathcal{T}_{h}$, sharing
at least one vertex with $K$ .

Hence, $\mathbf{\boldsymbol{\varphi}}_{{\mathrm{D}},h}$ belongs to the
discrete trace space on ${\Gamma_{\mathrm{D}}}$ given by

$\mathbf{H}^{1/2}_{h}({\Gamma_{\mathrm{D}}})\;:=\;\Big{\\{}\boldsymbol{\psi}_{{\mathrm{D}},h}\in\mathbf{C}({\Gamma_{\mathrm{D}}}):\quad\boldsymbol{\psi}_{{\mathrm{D}},h}|_{e}\in\mathbf{P}_{k+1}(e)\quad\forall\,e\in\mathcal{E}_{{\Gamma_{\mathrm{D}}}}\Big{\\}},$

where $\mathcal{E}_{{\Gamma_{\mathrm{D}}}}$ stands for the set of edges/faces
on ${\Gamma_{\mathrm{D}}}$.

In this way, defining
$\mathbb{H}_{h}:=\mathbb{H}^{\boldsymbol{t}}_{h}\times\mathbf{H}^{\boldsymbol{\sigma}}_{h}\times\mathbf{H}^{\boldsymbol{u}}_{h}$
and denoting
$\vec{\boldsymbol{\xi}}_{h}:=(\boldsymbol{t}_{h},\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})$,
the underlying Galerkin scheme given by the discrete counterpart of (3.6),
reads: Find
$(\vec{\boldsymbol{\xi}}_{h},\boldsymbol{\varphi}_{h})\in\mathbb{H}_{h}\times\mathbf{H}^{\boldsymbol{\varphi}}_{h}$,
with
$\boldsymbol{\varphi}_{h}\big{|}_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}},h}$,
such that

$\displaystyle\mathbf{A}_{\boldsymbol{\varphi}_{h}}(\vec{\boldsymbol{\xi}}_{h},\vec{\boldsymbol{\eta}}_{h})+\mathbf{B}_{\boldsymbol{u}_{h}}(\vec{\boldsymbol{\xi}}_{h},\vec{\boldsymbol{\eta}}_{h})$
$\displaystyle=$ $\displaystyle
F_{\boldsymbol{\varphi}_{h}}(\vec{\boldsymbol{\eta}}_{h})\qquad\quad\;\;\forall\,\vec{\boldsymbol{\eta}}_{h}\in\mathbb{H}_{h},$
(4.8a)
$\displaystyle\mathbf{C}(\boldsymbol{\varphi}_{h},\boldsymbol{\psi}_{h})$
$\displaystyle=$ $\displaystyle
G_{\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h}}(\boldsymbol{\psi}_{h})\qquad\forall\,\boldsymbol{\psi}_{h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}.$
(4.8b)

Throughout the rest of this section we adopt the discrete analogue of the
fixed point strategy introduced in Section 3.2. Indeed, denoting
$\mathbf{H}_{h}:=\mathbf{H}^{\boldsymbol{u}}_{h}\times\mathbf{H}^{\boldsymbol{\varphi}}_{h}$,
we define the operator
$\text{\Large$\mathpzc{A}$}_{h}:\mathbf{H}_{h}\rightarrow\mathbb{H}_{h}$ by

$\text{\Large$\mathpzc{A}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\ =\
(\text{\Large$\mathpzc{A}$}_{1,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}),\text{\Large$\mathpzc{A}$}_{2,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}),\text{\Large$\mathpzc{A}$}_{3,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}))\
:=\
\vec{\boldsymbol{\xi}}_{h}\qquad\forall\,(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\in\mathbf{H}_{h}$

where
$\vec{\boldsymbol{\xi}}_{h}=(\boldsymbol{t}_{h},\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})\in\mathbb{H}_{h}$
is the unique solution of the discrete problem (4.8a) with
$(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})$ instead pf
$(\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})$, that is

$\mathbf{A}_{\boldsymbol{\phi}_{h}}(\vec{\boldsymbol{\xi}}_{h},\vec{\boldsymbol{\eta}}_{h})+\mathbf{B}_{\boldsymbol{w}_{h}}(\vec{\boldsymbol{\xi}}_{h},\vec{\boldsymbol{\eta}}_{h})\
=\
F_{\boldsymbol{\phi}_{h}}(\vec{\boldsymbol{\eta}}_{h})\qquad\forall\,\vec{\boldsymbol{\eta}}_{h}\in\mathbb{H}\,,$
(4.9)

where the bilinear forms $\mathbf{A}_{\boldsymbol{\phi}_{h}}$,
$\mathbf{B}_{\boldsymbol{w}_{h}}$ and the functional
$F_{\boldsymbol{\phi}_{h}}$ are those corresponding to (3.7), (3.8), and
(3.10), respectively, with $\boldsymbol{w}=\boldsymbol{w}_{h}$ and
$\boldsymbol{\phi}=\boldsymbol{\phi}_{h}$.

In addition, we introduce the operator
$\text{\Large$\mathpzc{B}$}_{h}:\mathbf{H}_{h}\rightarrow\mathbf{H}^{\boldsymbol{\varphi}}_{h}$
defined as

$\text{\Large$\mathpzc{B}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\ :=\
\boldsymbol{\varphi}_{h}\qquad\forall\,(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\in\mathbf{H}_{h}\,,$

where $\boldsymbol{\varphi}_{h}$ is the unique solution of the discrete
problem (4.8b) with $(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})$ instead pf
$(\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})$, that is

$\mathbf{C}(\boldsymbol{\varphi}_{h},\boldsymbol{\psi}_{h})\ =\
G_{\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}}(\boldsymbol{\psi}_{h})\qquad\forall\,\boldsymbol{\psi}_{h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}\,,$
(4.10)

where the bilinear form $\mathbf{C}$, and the functional
$G_{\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}}$ are defined as in (3.9) and
(3.11), respectively, with $\boldsymbol{w}=\boldsymbol{w}_{h}$ and
$\boldsymbol{\phi}=\boldsymbol{\phi}_{h}$.

Therefore, by introducing the operator
$\text{\large$\mathscr{L}$}_{h}:\mathbf{H}_{h}\rightarrow\mathbf{H}_{h}$ as

$\text{\large$\mathscr{L}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\ :=\
(\text{\Large$\mathpzc{A}$}_{3,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}),\,\text{\Large$\mathpzc{B}$}_{h}(\text{\Large$\mathpzc{A}$}_{3,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}),\boldsymbol{\phi}_{h}))\qquad\forall\,(\boldsymbol{w}_{h},\phi_{h})\in\mathbf{H}_{h}$

we realize that solving (4.8) is equivalent to seeking a fixed point of
$\text{\large$\mathscr{L}$}_{h}$, that is: Find
$(\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})\in\mathbf{H}_{h}$ such that

$\text{\large$\mathscr{L}$}_{h}(\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})\
=\ (\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})\,.$ (4.11)

Certainly, all the above makes sense if we guarantee that the discrete
problems (4.9) and (4.10) are well-posed. This is precisely the purpose of the
next section.

###  4.1 Well-posedness of the uncoupled problems

In this section, we establish the well-posedness of both (4.9) and (4.10),
thus confirming that the operators $\text{\Large$\mathpzc{A}$}_{h}$,
$\text{\Large$\mathpzc{B}$}_{h}$, and hence $\text{\large$\mathscr{L}$}_{h}$,
are well-defined. We begin with the corresponding result for
$\text{\Large$\mathpzc{A}$}_{h}$, which actually follows almost verbatim to
that of its continuous counterpart $\mathpzc{A}$ (see Lemma 3.3), and the
proof can be omitted.

###### Lemma 4.2

Assume that
$\kappa_{1}\in\big{(}0,2\delta_{3}\big{(}2\nu_{1}-\frac{\kappa_{3}\nu_{2}}{\delta_{1}}\big{)}\big{)}$,
$\kappa_{2}\in(0,2\delta_{2})$,
$\kappa_{3}\in\big{(}0,\frac{2\nu_{1}\delta_{1}}{\nu_{2}}\big{)}$ with
$\delta_{1}\in\big{(}0,\frac{1}{\nu_{2}}\big{)}$,
$\delta_{2}\in\big{(}0,\frac{2}{\gamma}\big{)}$, $\delta_{3}\in(0,1)$. Then,
there exists $r_{0}>0$ (cf. (3.22)) such that for each $r\in(0,r_{0})$,
problem (4.9) has a unique solution
$\text{\Large$\mathpzc{A}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}):=\vec{\boldsymbol{\xi}}_{h}\in\mathbb{H}_{h}$,
for each $(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\in\mathbf{H}_{h}$ with
$\|\boldsymbol{w}_{h}\|_{1,\Omega}\leq r$. Moreover, there exists
$C_{\\!\text{\footnotesize$\mathpzc{A}$}}$, independent of
$(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})$, such that

$\|\text{\Large$\mathpzc{A}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\|_{\mathbb{H}_{h}}\
=\ \|\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}_{h}}\ \leq\
C_{\\!\text{\footnotesize$\mathpzc{A}$}}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\phi}_{h}\|_{0,\Omega}\qquad\forall\,(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\in\mathbf{H}_{h}\,.$

We now provide the discrete version of Lemma 3.4.

###### Lemma 4.3

For each $(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\in\mathbf{H}_{h}$,
problem (4.10) has a unique solution
$\boldsymbol{\varphi}_{h}\in\mathbf{H}_{h}^{\boldsymbol{\varphi}}$, with
$\boldsymbol{\varphi}_{h}|_{\Gamma_{{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}},h}$.
Moreover, there exists a constant
$\widetilde{C}_{\\!\text{\footnotesize$\mathpzc{B}$}}>0$ independent of
$(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})$, such that

$\|\text{\Large$\mathpzc{B}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\|_{1,\Omega}\
=\ \|\boldsymbol{\varphi}_{h}\|_{1,\Omega}\ \leq\
\widetilde{C}_{\\!\text{\footnotesize$\mathpzc{B}$}}\big{\\{}\|\boldsymbol{w}_{h}\|_{1,\Omega}\|\boldsymbol{\phi}_{h}\|_{1,\Omega}\,+\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,+\,\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}\big{\\}}\,.$
(4.12)

Proof. Let
${\color[rgb]{0,0,0}\boldsymbol{\varphi}_{1,h}:=\mathcal{I}_{h}^{SZ}(\boldsymbol{\varphi}_{1})}\in\mathbf{H}_{h}^{\boldsymbol{\varphi}}$
which satisfies
$\boldsymbol{\varphi}_{1,h}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}},h}$.
Then, similar to Lemma 3.4, we consider the auxiliary discrete problem: Find
$\boldsymbol{\varphi}_{0,h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}$
such that

$\mathbf{C}(\boldsymbol{\varphi}_{0,h},\boldsymbol{\psi}_{h})\ =\
\widetilde{G}_{\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}}(\boldsymbol{\psi}_{h})\qquad\forall\,\boldsymbol{\psi}_{h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}\,,$
(4.13)

where

$\widetilde{G}_{\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}}(\boldsymbol{\psi}_{h})\
:=\
-\int_{\Omega}(\nabla\boldsymbol{\phi}_{h})\boldsymbol{w}_{h}\cdot\boldsymbol{\psi}_{h}\,-\,\int_{\Omega}\mathbb{K}\nabla\boldsymbol{\varphi}_{1,h}:\nabla\boldsymbol{\psi}_{h}\qquad\forall\,\boldsymbol{\psi}_{h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}\,.$

Next, the boundedness and the ellipticity of $\mathbf{C}$ are obtained exactly
as in the proof of Lemma 3.4 with the same ellipticity constant
$\widetilde{\alpha}$ given by (3.31). On the other hand, reasoning as in the
proof of Lemma 3.4, and applying the stability property of
$\mathcal{I}_{h}^{SZ}$ (cf. (4.6)) and the equinormic property of
$\boldsymbol{\varphi}_{1}$ (cf. (3.27)), we easily deduce that

$|\widetilde{G}_{\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}}(\boldsymbol{\psi}_{h})|\
\leq\
\Big{\\{}c_{1}(\Omega)\|\boldsymbol{w}_{h}\|_{1,\Omega}\|\boldsymbol{\phi}_{h}\|_{1,\Omega}\,+{\color[rgb]{0,0,0}\,c\,\|\mathbb{K}\|_{\infty,\Omega}\,\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}}\Big{\\}}\|\boldsymbol{\psi}_{h}\|_{1,\Omega},$

$\forall\,\boldsymbol{\psi}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}$,
which says that
$\widetilde{G}_{\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}}\in[\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}]^{\prime}$
and

$\|\widetilde{G}_{\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}}\|\leq
c_{1}(\Omega)\|\boldsymbol{w}_{h}\|_{1,\Omega}\|\boldsymbol{\phi}_{h}\|_{1,\Omega}\,+{\color[rgb]{0,0,0}\,c\,\|\mathbb{K}\|_{\infty,\Omega}\,\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}}.$

Therefore, a direct application of the Lax-Milgram lemma implies that there
exists a unique
$\boldsymbol{\varphi}_{0,h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}$
that satisfies (4.13) with

$\|\boldsymbol{\varphi}_{0,h}\|_{1,\Omega}\ \leq\
\frac{1}{\widetilde{\alpha}}\,\Big{\\{}c_{1}(\Omega)\|\boldsymbol{w}_{h}\|_{1,\Omega}\|\boldsymbol{\phi}_{h}\|_{1,\Omega}\,+{\color[rgb]{0,0,0}\,c\,\|\mathbb{K}\|_{\infty,\Omega}\,\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}}\Big{\\}}.$

Then,
$\boldsymbol{\varphi}_{h}\;:=\;\boldsymbol{\varphi}_{0,h}+\boldsymbol{\varphi}_{1,h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h}$,
which in fact satisfies that
$\boldsymbol{\varphi}_{h}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}},h}$,
is the unique solution of (4.10). In addition, the estimate (4.12) holds with
${\color[rgb]{0,0,0}\widetilde{C}_{\\!\text{\footnotesize$\mathpzc{B}$}}:=\max\\{\widetilde{\alpha}^{-1}\,c_{1}(\Omega),\widetilde{c}\,(\widetilde{\alpha}^{-1}\,\|\mathbb{K}\|_{\infty,\Omega}+1)\\}}$.
$\square$

###  4.2 Solvability analysis of the fixed point equation

In this section we establish the solvability of the fixed point problem (4.11)
by applying the Brouwer fixed-point theorem [16, Theorem 9.9-2]. To this end,
we begin with the discrete version of lemma 3.5.

###### Lemma 4.4

Given $r\in(0,r_{0})$, with $r_{0}$ given by (3.22), we let

$W_{r}^{h}\ :=\
\big{\\{}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\in\mathbf{H}_{h}:\quad\|(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\|_{\mathbf{H}}\,\leq\,r\big{\\}}\,,$

and assume that the data satisfy

$c(r)\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,+\,\widetilde{C}_{\\!\text{\footnotesize$\mathpzc{B}$}}\,\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}\
\leq\ r$ (4.14)

where
${\color[rgb]{0,0,0}c(r):=rC_{\\!\text{\footnotesize$\mathpzc{A}$}}(1+r\widetilde{C}_{\\!\text{\footnotesize$\mathpzc{B}$}})}$,
with $C_{\\!\text{\footnotesize$\mathpzc{A}$}}$ and
$\widetilde{C}_{\\!\text{\footnotesize$\mathpzc{B}$}}$ as in (3.16) and
(4.12), respectively. Then, there holds
$\text{\large$\mathscr{L}$}_{h}(W_{r}^{h})\subseteq W_{r}^{h}$.

In order to provide the discrete analogue of Lemma 3.6, we notice in advance
that, instead of the regularity assumptions employed in the continuous case,
which are not applicable in the present case, we simple utilize an
$\mathrm{L}^{4}$-$\mathrm{L}^{4}$-$\mathrm{L}^{2}$ argument.

###### Lemma 4.5

Let $r\in(0,r_{0})$ with $r_{0}$ given by (3.22). Then, there exists a
constant $\widetilde{L}_{\\!\text{\footnotesize$\mathpzc{A}$}}>0$, independent
of $r$, such that for all
$(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}),(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})\in
W_{r}^{h}$, with
$\|\boldsymbol{w}_{h}\|_{1,\Omega},\|\widetilde{\boldsymbol{w}}_{h}\|_{1,\Omega}\,\leq\;r$,
there holds

$\begin{array}[]{l}\|\text{\Large$\mathpzc{A}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})-\text{\Large$\mathpzc{A}$}_{h}(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})\|_{\mathbb{H}}\
\leq\
\widetilde{L}_{\\!\text{\footnotesize$\mathpzc{A}$}}\Big{\\{}\|\text{\Large$\mathpzc{A}$}_{3,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\|_{1,\Omega}\;\|\boldsymbol{w}_{h}-\widetilde{\boldsymbol{w}}_{h}\|_{1,\Omega}\\\\[4.30554pt]
\qquad\qquad\qquad\qquad+\
\|\text{\Large$\mathpzc{A}$}_{1,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\|_{0,4,\Omega}\,\|\boldsymbol{\phi}-\widetilde{\boldsymbol{\phi}}\|_{0,4,\Omega}\,+\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\phi}_{h}-\widetilde{\boldsymbol{\phi}}_{h}\|_{1,\Omega}\Big{\\}}\,.\end{array}$
(4.15)

Proof. It proceeds exactly as in the proof of Lemma 3.6, except for the
derivation of the discrete analogue of (3.40), where, instead of choosing the
values of $p,q$ determined by the regularity parameter $\varepsilon$, it
suffices to take $p=q=2$, thus obtaining

$\begin{array}[]{l}\displaystyle|(\mathbf{A}_{\widetilde{\boldsymbol{\phi}}_{h}}-\mathbf{A}_{\boldsymbol{\phi}_{h}})(\vec{\boldsymbol{\xi}}_{h},\vec{\boldsymbol{\xi}}_{h}-\vec{\boldsymbol{\zeta}}_{h})|\
\leq\
2\,L_{\nu}(2+\kappa_{3}^{2})^{1/2}\,\|\boldsymbol{t}_{h}\|_{0,4,\Omega}\,\|\boldsymbol{\phi}_{h}-\widetilde{\boldsymbol{\phi}}_{h}\,\|_{0,4,\Omega}\,\|\vec{\boldsymbol{\xi}}_{h}-\vec{\boldsymbol{\zeta}}_{h}\|_{\mathbb{H}}\,,\end{array}$

for all
$(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}),(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})$,
with
$\vec{\boldsymbol{\xi}}_{h}=(\boldsymbol{t}_{h},\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h}):=\text{\Large$\mathpzc{A}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\in\mathbb{H}_{h}$
and
$\vec{\boldsymbol{\zeta}}_{h}=(\widetilde{\boldsymbol{t}}_{h},\widetilde{\boldsymbol{\sigma}}_{h},\widetilde{\boldsymbol{u}}_{h}):=\text{\Large$\mathpzc{A}$}_{h}(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})\in\mathbb{H}_{h}$.
Thus, since the elements of $\mathbb{H}_{h}^{\,\boldsymbol{t}}$ are piecewise
polynomials, we can guarantee that
$\|\boldsymbol{t}_{h}\|_{0,4,\Omega}<\infty$ for each
$\boldsymbol{t}_{h}\in\mathbb{H}_{h}^{\,\boldsymbol{t}}$. The proof concludes
with
$\widetilde{L}_{\\!\text{\footnotesize$\mathpzc{A}$}}:=\frac{2}{\alpha(\Omega)}\max\\{c_{1}(\Omega)\,(2+\kappa_{3}^{2})^{1/2},$
$2\,L_{\nu}(2+\kappa_{3}^{2})^{1/2},\sqrt{2}\,(2+\kappa_{2}^{2})^{1/2}\\}$.
$\square$

The discrete version of Lemma 3.7 is given as follows.

###### Lemma 4.6

Let $L_{\\!\text{\footnotesize$\mathpzc{B}$}}>0$ as in Lemma 3.7. Then, for
all
$(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}),(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})\in
W_{r}^{h}$ there holds

$\begin{array}[]{l}\|\text{\Large$\mathpzc{B}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})-\text{\Large$\mathpzc{B}$}_{h}(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})\|_{1,\Omega}\
\leq\
L_{\\!\text{\footnotesize$\mathpzc{B}$}}\Big{\\{}\|\boldsymbol{w}_{h}-\widetilde{\boldsymbol{w}}_{h}\|_{1,\Omega}\,+\,\|\boldsymbol{\phi}_{h}-\widetilde{\boldsymbol{\phi}}_{h}\|_{1,\Omega}\Big{\\}}\,.\end{array}$
(4.16)

Proof. It corresponds to an adaptation of the proof of Lemma 3.7 to the
discrete context. $\square$

Now, combining Lemmas 4.5 and 4.6, and employing the continuous injection of
$\mathbf{H}^{1}(\Omega)$ into $\mathbf{L}^{4}(\Omega)$, we can prove the
discrete version of Lemma 3.8.

###### Lemma 4.7

Let $r$ and $W_{r}^{h}$ as in Lemma 4.4. Then, there exits a positive constant
$C$, depending only on $\widetilde{L}_{\\!\text{\footnotesize$\mathpzc{A}$}}$
and $L_{\\!\text{\footnotesize$\mathpzc{B}$}}$ (cf. (4.15) and (4.16),
respectively), such that for all
$(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h}),(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})\in
W_{r}^{h}$ there holds

$\displaystyle\|\text{\large$\mathscr{L}$}_{h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})-\text{\large$\mathscr{L}$}_{h}(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})\|_{\mathbf{H}}$
$\displaystyle\leq$ $\displaystyle
C\,\Big{\\{}\|\text{\Large$\mathpzc{A}$}_{3,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\|_{1,\Omega}\,+\,c_{1}(\Omega)\,\|\text{\Large$\mathpzc{A}$}_{1,h}(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})\|_{0,4,\Omega}$
$\displaystyle\phantom{C\,\Big{\\{}}+\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,+\,L_{\\!\text{\footnotesize$\mathpzc{B}$}}\Big{\\}}\,\|(\boldsymbol{w}_{h},\boldsymbol{\phi}_{h})-(\widetilde{\boldsymbol{w}}_{h},\widetilde{\boldsymbol{\phi}}_{h})\|_{\mathbf{H}}\,.$

More precisely, there holds
$C=\max\\{(1+L_{\\!\text{\footnotesize$\mathpzc{B}$}})\widetilde{L}_{\\!\text{\footnotesize$\mathpzc{A}$}},1\\}$.

Consequently, since the foregoing lemma confirms the continuity of
$\text{\large$\mathscr{L}$}_{h}$, by a straightforward application of Brouwer
fixed point theorem (cf. [16, Theorem 9.9-2]) on the convex and compact set
$W_{r}^{h}\subseteq\mathbf{H}_{h}$, we can provide the main result of this
section.

###### Theorem 4.8

Suppose that the parameters $\kappa_{1},\kappa_{2}$ and $\kappa_{3}$ satisfy
the conditions required by Lemma 4.2. Let $r$ and $W_{r}^{h}$ as in Lemma 4.4,
and assume that the data satisfy (4.14). Then, problem (4.8) has at least one
solution
$(\vec{\boldsymbol{\xi}}_{h},\boldsymbol{\varphi}_{h})\in\mathbb{H}_{h}\times\mathbf{H}_{h}^{\boldsymbol{\varphi}}$
such that
$\boldsymbol{\varphi}_{h}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}},h}$,
with $(\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})\in W_{r}^{h}$, and there
holds

$\|\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}\ \leq\
C_{\\!\text{\footnotesize$\mathpzc{A}$}}\,r\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,,$

and

$\|\boldsymbol{\varphi}_{h}\|_{1,\Omega}\ \leq\
\widetilde{C}_{\\!\text{\footnotesize$\mathpzc{B}$}}\,\big{\\{}r\,\|\boldsymbol{u}_{h}\|_{1,\Omega}\,+\,{\color[rgb]{0,0,0}\|\boldsymbol{\varphi}_{{\mathrm{D}}}\|_{1/2,{\Gamma_{\mathrm{D}}}}}\big{\\}}\,.$

##  5 A priori error analysis

We now aim to derive the a priori error estimates for the Galerkin scheme
given by (4.8). To this end, given
$((\boldsymbol{t},\boldsymbol{\sigma},\boldsymbol{u}),\boldsymbol{\varphi}):=(\vec{\boldsymbol{\xi}},\boldsymbol{\varphi})\in\mathbb{H}\times\mathbf{H}^{1}(\Omega)$,
with $(\boldsymbol{u},\boldsymbol{\varphi})\in W_{r}$ and
$((\boldsymbol{t}_{h},\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h}),\boldsymbol{\varphi}_{h}):=(\vec{\boldsymbol{\xi}}_{h},\boldsymbol{\varphi}_{h})\in\mathbb{H}_{h}\times\mathbf{H}^{\boldsymbol{\varphi}}_{h}$,
with $(\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h})\in W_{r}^{h}$ solutions of
(3.6) and (4.8), respectively, we first observe that the above problems can be
rewritten as two pairs of corresponding continuous and discrete formulations,
namely

$\displaystyle\mathbf{A}_{\boldsymbol{\varphi}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})+\mathbf{B}_{\boldsymbol{u}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}})$
$\displaystyle=$ $\displaystyle
F_{\boldsymbol{\varphi}}(\vec{\boldsymbol{\eta}})\qquad\quad\forall\,\vec{\boldsymbol{\eta}}\in\mathbb{H}\,,$
(5.1a)
$\displaystyle\mathbf{A}_{\boldsymbol{\varphi}_{h}}(\vec{\boldsymbol{\xi}}_{h},\vec{\boldsymbol{\eta}}_{h})+\mathbf{B}_{\boldsymbol{u}_{h}}(\vec{\boldsymbol{\xi}}_{h},\vec{\boldsymbol{\eta}}_{h})$
$\displaystyle=$ $\displaystyle
F_{\boldsymbol{\varphi}_{h}}(\vec{\boldsymbol{\eta}}_{h})\qquad\forall\,\vec{\boldsymbol{\eta}}_{h}\in\mathbb{H}_{h}\,,$
(5.1b)

and

$\displaystyle\mathbf{C}(\boldsymbol{\varphi},\boldsymbol{\psi})$
$\displaystyle=$ $\displaystyle
G_{\boldsymbol{u},\boldsymbol{\varphi}}(\boldsymbol{\psi})\qquad\quad\,\forall\,\boldsymbol{\psi}\in\mathrm{H}^{1}_{{\Gamma_{\mathrm{D}}}}(\Omega)\,,$
(5.2a)
$\displaystyle\mathbf{C}(\boldsymbol{\varphi}_{h},\boldsymbol{\psi}_{h})$
$\displaystyle=$ $\displaystyle
G_{\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h}}(\boldsymbol{\psi}_{h})\qquad\forall\,\boldsymbol{\psi}_{h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h,{\Gamma_{\mathrm{D}}}}\,.$
(5.2b)

Our goal is to obtain an upper bound for the error
$\displaystyle\|(\vec{\boldsymbol{\xi}},\boldsymbol{\varphi})-(\vec{\boldsymbol{\xi}}_{h},\boldsymbol{\varphi}_{h})\|_{\mathbb{H}\times\mathbf{H}^{1}(\Omega)}$.
For this purpose, we first recall from [43, Theorem 11.1] an abstract result
that corresponds to the standar Strang Lemma for elliptic variational
problems, which will be straightforwardly applied to the pair (5.1a)-(5.1b).

###### Lemma 5.1

Let $H$ be a Hilbert space, $F\in H^{\prime}$, and $A:H\times H\to\text{\rm
R}$ be a bounded and $H$-elliptic bilinear form. In addition, let
$\\{H_{h}\\}_{h>0}$ be a sequence of finite dimensional subspaces of $H$, and
for each $h>0$ consider a bounded bilinear form $A_{h}:H_{h}\times
H_{h}\to\text{\rm R}$ and a functional $F_{h}\in H_{h}^{\prime}$. Assume that
the family $\\{A_{h}\\}_{h>0}$ is uniformly elliptic, that is, there exists a
constant $\beta>0$, independent of $h$ such that

$A_{h}(v_{h},v_{h})\ \geq\ \beta\,\|v_{h}\|^{2}_{H}\qquad\forall\,v_{h}\in
H_{h}\,,\quad\forall\,h>0\,.$

In turn, let $u\in H$ and $u_{h}\in H_{h}$ such that

$A(u,v)\ =\ F(v)\qquad\forall\,v\in H\qquad\text{and}\qquad
A_{h}(u_{h},v_{h})\ =\ F_{h}(v_{h})\qquad\forall\,v_{h}\in H_{h}\,.$

Then for each $h>0$ there holds

$\displaystyle\|u-u_{h}\|_{H}$ $\displaystyle\leq$ $\displaystyle
C\,\Bigg{\\{}\sup_{\begin{subarray}{c}{w_{h}\in H_{h}}\\\ {w_{h}\neq\bf
0}\end{subarray}}\frac{\big{|}\,F(w_{h})-F_{h}(w_{h})\,\big{|}}{\|w_{h}\|_{H}}\,+\,\inf_{\begin{subarray}{c}{v_{h}\in
H_{h}}\\\ {v_{h}\neq\bf 0}\end{subarray}}\Bigg{(}\|u-v_{h}\|_{H}$
$\displaystyle\phantom{C\,\Bigg{\\{}}+\,\sup_{\begin{subarray}{c}{w_{h}\in
H_{h}}\\\ {w_{h}\neq\bf
0}\end{subarray}}\frac{\big{|}\,A(v_{h},w_{h})-A_{h}(v_{h},w_{h})\,\big{|}}{\|w_{h}\|_{H}}\Bigg{)}\Bigg{\\}}\,,$

where $C:=\beta^{-1}\,\max\\{1,\|A\|\\}$.

In what follows, as usual, we denote

$\mathrm{dist}\big{(}\vec{\boldsymbol{\xi}},\mathbb{H}_{h}\big{)}\ :=\
\inf_{\vec{\boldsymbol{\eta}}_{h}\in\mathbb{H}_{h}}\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\eta}}_{h}\|_{\mathbb{H}}\qquad\text{and}\qquad\mathrm{dist}\big{(}\boldsymbol{\varphi},\mathbf{H}^{\varphi}_{h}\big{)}\
:=\
\inf_{\boldsymbol{\psi}_{h}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h}}\|\boldsymbol{\varphi}-\boldsymbol{\psi}_{h}\|_{1,\Omega}\,.$

We now derive the a preliminary estimate for the error
$\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}=\|(\boldsymbol{t},\boldsymbol{\sigma},\boldsymbol{u})-(\boldsymbol{t}_{h},\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})\|_{\mathbb{H}}$.

###### Lemma 5.2

There exists a constant $C_{\tt{ST}}>0$ independent of $h$, such that

$\displaystyle\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}$
$\displaystyle\leq$ $\displaystyle
C_{\tt{ST}}\Big{\\{}\mathrm{dist}\big{(}\vec{\boldsymbol{\xi}},\mathbb{H}_{h}\big{)}\,+\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\,+\,\|\boldsymbol{t}\|_{\varepsilon,\Omega}\,\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}$
(5.3)
$\displaystyle\phantom{C_{\tt{ST}}\Big{\\{}}+\,\|\boldsymbol{u}\|_{1,\Omega}\,\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1,\Omega}\Big{\\}}\,.$

Proof. From Lemma 3.3 we have that the bilinear forms
$\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}}$ and
$\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}}$ are
both bounded and uniformly elliptic, with ellipticity constant
$\frac{\alpha(\Omega)}{2}$ (cf. (3.21)). In turn, $F_{\boldsymbol{\varphi}}$
and $F_{\boldsymbol{\varphi}_{h}}$ are linear bounded functionals in
$\mathbb{H}$ and $\mathbb{H}_{h}$ respectively. Then, a straightforward
application of Lemma 5.1 to the context given by (5.1a)-(5.1b), gives

$\begin{array}[]{l}\displaystyle\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}\
\leq\
C\,\Bigg{\\{}\sup_{\begin{subarray}{c}{\vec{\boldsymbol{\eta}}_{h}\in\mathbb{H}_{h}}\\\
{\vec{\boldsymbol{\eta}}_{h}\neq\bf
0}\end{subarray}}\frac{\big{|}\,F_{\boldsymbol{\varphi}}(\vec{\boldsymbol{\eta}}_{h})-F_{\boldsymbol{\varphi}_{h}}(\vec{\boldsymbol{\eta}}_{h})\,\big{|}}{\|\vec{\boldsymbol{\eta}}_{h}\|_{\mathbb{H}}}\\\\[21.52771pt]
\displaystyle\qquad\qquad+\
\inf_{\begin{subarray}{c}{\vec{\boldsymbol{\zeta}}_{h}\in\mathbb{H}_{h}}\\\
{\vec{\boldsymbol{\zeta}}_{h}\neq\bf
0}\end{subarray}}\Bigg{(}\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}_{h}\|_{\mathbb{H}}\,+\,\sup_{\begin{subarray}{c}{\vec{\boldsymbol{\eta}}_{h}\in\mathbb{H}_{h}}\\\
{\vec{\boldsymbol{\eta}}_{h}\neq\bf
0}\end{subarray}}\frac{\big{|}\,(\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}})(\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})-(\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}})(\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})\,\big{|}}{\|\vec{\boldsymbol{\eta}}_{h}\|_{\mathbb{H}}}\Bigg{)}\Bigg{\\}}\,,\end{array}$
(5.4)

where
$C:=\frac{2}{\alpha(\Omega)}\max\\{1,\|\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}}\|\\}$.
It is important to recall here, from (3.19), that
$\|\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}}\|$ depends
only on $\kappa_{1}$, $\kappa_{2}$, $\kappa_{3}$, $\nu_{2}$, $\gamma$,
$c_{1}(\Omega)$ and $\|\boldsymbol{u}\|_{1,\Omega}$, where
$\|\boldsymbol{u}\|_{1,\Omega}\leq r$. Furthermore, we now proceed to estimate
each term appearing at the right-hand side of (5.4). In order to do that, we
first apply the same arguments employed to obtain (3.38), to find that

$|F_{\boldsymbol{\varphi}}(\vec{\boldsymbol{\eta}}_{h})-F_{\boldsymbol{\varphi}_{h}}(\vec{\boldsymbol{\eta}}_{h})|\
\leq\
\sqrt{2}\,(2+\kappa_{2}^{2})^{1/2}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\,\|\vec{\boldsymbol{\eta}}_{h}\|_{\mathbb{H}}\,.$
(5.5)

Next, in order to estimate the last supremum in (5.4), we add and subtract
$\vec{\boldsymbol{\xi}}:=(\boldsymbol{t},\boldsymbol{\sigma},\boldsymbol{u})$,
we note that

$\displaystyle(\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}})(\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})-(\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}})(\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})$
$\displaystyle=$
$\displaystyle(\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}_{h})-(\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}_{h})-(\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}})(\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})+(\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}})(\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})$
$\displaystyle=$
$\displaystyle(\mathbf{A}_{\boldsymbol{\varphi}}-\mathbf{A}_{\boldsymbol{\varphi}_{h}})(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}_{h})\,+\,\mathbf{B}_{\boldsymbol{u}-\boldsymbol{u}_{h}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}_{h})-(\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}})(\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})+(\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}})(\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})$

where, applying the same approach used in (3.40) and (3.39), together with
(3.41) and the continuous embedding
$\mathbf{H}^{1}(\Omega)\to\mathbf{L}^{n/\varepsilon}(\Omega)$ with constant
$\widetilde{C}_{\varepsilon}$, and the boundedness of the bilinear forms
$\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}}$ and
$\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}}$, it
follows that

$\begin{array}[]{r}|(\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}})(\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})-(\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}})(\vec{\boldsymbol{\zeta}}_{h},\vec{\boldsymbol{\eta}}_{h})|\
\leq\
\Big{\\{}2\,L_{\nu}C_{\varepsilon}\widetilde{C}_{\varepsilon}(2+\kappa_{3}^{2})^{1/2}\,\|\boldsymbol{t}\|_{\varepsilon,\Omega}\,\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\\\\[8.61108pt]
\qquad+\
c_{1}(\Omega)\,(2+\kappa_{3}^{2})^{1/2}\,\|\boldsymbol{u}\|_{1,\Omega}\,\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1,\Omega}\\\\[8.61108pt]
\qquad+\
\big{(}\|\mathbf{A}_{\boldsymbol{\varphi}}+\mathbf{B}_{\boldsymbol{u}}\|+\|\mathbf{A}_{\boldsymbol{\varphi}_{h}}+\mathbf{B}_{\boldsymbol{u}_{h}}\|\big{)}\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\zeta}}_{h}\|_{\mathbb{H}}\Big{\\}}\|\vec{\boldsymbol{\eta}}_{h}\|_{\mathbb{H}}\,.\end{array}$
(5.6)

In this way, by replacing (5.5) and (5.6) back into (5.4), we obtain (5.3)
with $C_{\tt{ST}}$ is a positive constant depending on $\alpha(\Omega)$,
$L_{\nu}$, $C_{\varepsilon}$, $\widetilde{C}_{\varepsilon}$, $c_{1}(\Omega)$,
$\kappa_{1}$, $\kappa_{2}$, $\kappa_{3}$, $\nu_{2}$, $\gamma$, and $r$.
$\square$

The following result present a estimate for the error
$\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}$.

###### Lemma 5.3

Assume that $r$ satisfy that

$r\,\frac{c_{1}(\Omega)}{\widetilde{\alpha}}\ \leq\ \frac{1}{2}\,,$ (5.7)

where $\widetilde{\alpha}$ is defined in (3.31), and
$c_{1}(\Omega):=\|\boldsymbol{i}\|^{2}$ is the boundedness constants of the
continuous injection
$\boldsymbol{i}:\mathbf{H}^{1}(\Omega)\to\mathbf{L}^{4}(\Omega)$. Then, there
exists a constant $\widetilde{C}_{\tt{ST}}>0$, independent of $h$, such that

$\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\ \leq\
\widetilde{C}_{\tt{ST}}\,\|\boldsymbol{\varphi}-\mathcal{I}^{SZ}_{h}(\boldsymbol{\varphi})\|_{1,\Omega}\,+\,\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1,\Omega}\,.$
(5.8)

where $\mathcal{I}_{h}^{SZ}$ denotes the Scott-Zhang interpolant operator
introduced in Section 4.

Proof. We proceed similarly as in the proof of [6, Lemma 5.3]. Indeed, by
applying the triangle inequality we have that

$\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\ \leq\
\|\boldsymbol{\varphi}-{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}\|_{1,\Omega}\,+\,\|{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h}\|_{1,\Omega},$
(5.9)

where
${\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}\;=\;\mathcal{I}^{SZ}_{h}(\boldsymbol{\varphi})}\in\mathbf{H}^{\boldsymbol{\varphi}}_{h}$.
Moreover, nothing that
$\boldsymbol{\varphi}|_{{\Gamma_{\mathrm{D}}}}=E(\boldsymbol{\varphi}_{{\mathrm{D}}})|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}}}$,
we can employ [45, eq. (2.17)] and (4.5), to obtain that
$\boldsymbol{\psi}_{h}|_{{\Gamma_{\mathrm{D}}}}=\boldsymbol{\varphi}_{{\mathrm{D}},h}$.
Now, utilizing the ellipticity of bilinear the form $\mathbf{C}(\cdot,\cdot)$
on $\mathbf{H}_{h,{\Gamma_{\mathrm{D}}}}^{\boldsymbol{\varphi}}$ (see the
proof of Lemma 4.3) with constant $\widetilde{\alpha}$, along with the fact
that
$\mathbf{C}(\boldsymbol{\varphi}_{h},{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})=G_{\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h}}({\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})$
and
$\mathbf{C}(\boldsymbol{\varphi},{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})=G_{\boldsymbol{u},\boldsymbol{\varphi}}({\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})$
(see (5.2)), we deduce that

$\displaystyle\widetilde{\alpha}\,\|{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}^{2}$
$\displaystyle\leq$
$\displaystyle\mathbf{C}(\boldsymbol{\varphi},{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})-\mathbf{C}(\boldsymbol{\varphi}-{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}},{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})-\mathbf{C}(\boldsymbol{\varphi}_{h},{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})$
(5.10) $\displaystyle\leq$
$\displaystyle|G_{\boldsymbol{u},\boldsymbol{\varphi}}({\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})-G_{\boldsymbol{u}_{h},\boldsymbol{\varphi}_{h}}({\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})|\,+\,|\mathbf{C}(\boldsymbol{\varphi}-{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}},{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h})|\,.$

Next, we apply the estimate (3.4) to bound the first term on the right-hand
side of (5.10), whereas for second term we use the boundedness of
$\mathbf{C}(\cdot,\cdot)$ (cf. (3.29)), and (5.7). Then, it follows that

$\displaystyle\|{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}$
$\displaystyle\leq$
$\displaystyle\frac{\|\mathbb{K}\|_{\infty,\Omega}}{\widetilde{\alpha}}\,\|\boldsymbol{\varphi}-{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}\|_{1,\Omega}\
+\
\frac{rc_{1}(\Omega)}{\widetilde{\alpha}}\,\Big{\\{}\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\,+\,\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1,\Omega}\Big{\\}}$
(5.11) $\displaystyle\leq$
$\displaystyle\frac{\|\mathbb{K}\|_{\infty,\Omega}}{\widetilde{\alpha}}\,\|\boldsymbol{\varphi}-{\color[rgb]{0,0,0}\boldsymbol{\psi}_{h}}\|_{1,\Omega}\
+\
\frac{1}{2}\,\Big{\\{}\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\,+\,\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1,\Omega}\Big{\\}}\,.$

Then, by replacing (5.11) back into (5.9), we conclude (5.8) with
$\widetilde{C}_{\tt{ST}}:=2(1+\widetilde{\alpha}^{-1}\|\mathbb{K}\|_{\infty,\Omega})$.
$\square$

We now proceed to combine Lemmas 5.2 and 5.3 to derive the Céa estimate for
the total error

$\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}\,+\,\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\,.$

In fact, by replacing the estimate for
$\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}$ given by (5.8)
into the right-hand side of (5.3), and using the fact that
$\|\boldsymbol{u}\|_{1,\Omega}\leq
r\,C_{\\!\text{\footnotesize$\mathpzc{A}$}}|\boldsymbol{\alpha}||{\boldsymbol{g}}|$
(cf. (3.48)) and
$\|\boldsymbol{t}\|_{\varepsilon,\Omega}\leq\widehat{C}(r)|\boldsymbol{\alpha}||{\boldsymbol{g}}|\|\boldsymbol{\varphi}\|_{0,\Omega}$
(cf. (3.33)) along with $\|\boldsymbol{\varphi}\|_{1,\Omega}\leq r$, we find
that

$\displaystyle\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}$
$\displaystyle\leq$ $\displaystyle
C_{\tt{ST}}\,\mathrm{dist}\big{(}\vec{\boldsymbol{\xi}},\mathbb{H}_{h}\big{)}\,+\,C_{\tt{ST}}\widetilde{C}_{\tt{ST}}\,\Big{\\{}|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|+\|\boldsymbol{t}\|_{\varepsilon,\Omega}\Big{\\}}\,\|\boldsymbol{\varphi}-\mathcal{I}^{SZ}_{h}(\boldsymbol{\varphi})\|$
$\displaystyle+\
C_{\tt{ST}}\,\Big{\\{}|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|+\|\boldsymbol{t}\|_{\varepsilon,\Omega}+\|\boldsymbol{u}\|_{1,\Omega}\Big{\\}}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1,\Omega}$
$\displaystyle\leq$ $\displaystyle
C_{\tt{ST}}\,\mathrm{dist}\big{(}\vec{\boldsymbol{\xi}},\mathbb{H}_{h}\big{)}\,+\,C_{0}\,\|\boldsymbol{\varphi}-\mathcal{I}^{SZ}_{h}(\boldsymbol{\varphi})\|_{1,\Omega}\,+\,C_{1}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\,\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1,\Omega}\,,$

where
$C_{0}:=C_{\tt{ST}}\widetilde{C}_{\tt{ST}}(1+r\widehat{C}(r))\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|$
and
$C_{1}:=C_{\tt{ST}}(1+r\widehat{C}(r)+rC_{\\!\text{\footnotesize$\mathpzc{A}$}})$.
In this way, assuming that the data $\boldsymbol{\alpha}$ and
${\boldsymbol{g}}$ satisfy that

$C_{1}\,|\boldsymbol{\alpha}|\,|{\boldsymbol{g}}|\ \leq\ \frac{1}{2}\,,$
(5.12)

we can conclude that

$\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}\ \leq\
2\,C_{\tt{ST}}\,\mathrm{dist}\big{(}\vec{\boldsymbol{\xi}},\mathbb{H}_{h}\big{)}\,+\,2\,C_{0}\,\|\boldsymbol{\varphi}-\mathcal{I}^{SZ}_{h}(\boldsymbol{\varphi})\|_{1,\Omega}.$
(5.13)

Consequently, we now can establish the following main result.

###### Theorem 5.4

Assume that $r$ and the data $\boldsymbol{\alpha}$ and ${\boldsymbol{g}}$ are
sufficiently small so that (5.7) and (5.12) hold, respectively. Then, there
exists a positive constant $C^{*}$, independent of $h$, such that

$\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}\,+\,\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\
\leq\
C^{*}\,\Big{\\{}\mathrm{dist}\big{(}\vec{\boldsymbol{\xi}},\mathbb{H}_{h}\big{)}\;+\;\|\boldsymbol{\varphi}-\mathcal{I}^{SZ}_{h}(\boldsymbol{\varphi})\|_{1,\Omega}\Big{\\}}\,.$
(5.14)

Proof. It follows straightforwardly from the Céa estimates (5.13) and (5.8).
$\square$

In order to provided the result concerning to the theoretical rate of
convergence of (4.8), we recall from [27], the approximation properties of the
specific finite element subspaces introduced in Section 4.

($\mathbf{AP}_{h}^{\boldsymbol{t}}$) There exists $c>0$, independent of $h$,
such that for each $s\in(0,k+1]$, and for each
$\boldsymbol{t}\in\mathbb{H}^{s}(\Omega)\cap\mathbb{L}^{2}_{\tt{tr}}(\Omega)$,
there holds

$\mathrm{dist}\big{(}\boldsymbol{t},\mathbb{H}^{\boldsymbol{t}}_{h}\big{)}\;\leq\;c\,h^{s}\,\|\boldsymbol{t}\|_{s,\Omega}\,.$

($\mathbf{AP}_{h}^{\boldsymbol{\sigma}}$) There exists $c>0$, independent of
$h$, such that for each $s\in(0,k+1]$, and for each
$\boldsymbol{\sigma}\in\mathbb{H}^{s}(\Omega)\cap\mathbb{H}_{0}({\mathbf{div}};\Omega)$,
with ${\mathbf{div}}\boldsymbol{\sigma}\in\mathbf{H}^{s}(\Omega)$, there holds

$\mathrm{dist}\big{(}\boldsymbol{\sigma},\mathbb{H}^{\boldsymbol{\sigma}}_{h}\big{)}\;\leq\;c\,h^{s}\,\big{\\{}\|\boldsymbol{\sigma}\|_{s,\Omega}\,+\,\|{\mathbf{div}}\boldsymbol{\sigma}\|_{s,\Omega}\big{\\}}\,.$

($\mathbf{AP}_{h}^{\boldsymbol{u}}$) There exists $c>0$, independent of $h$,
such that for each $s\in(0,k+1]$, and for each
$\boldsymbol{u}\in\mathbf{H}^{s+1}(\Omega)$, there holds

$\mathrm{dist}\big{(}\boldsymbol{u},\mathbf{H}^{\boldsymbol{u}}_{h}\big{)}\;\leq\;c\,h^{s}\,\|\boldsymbol{u}\|_{s+1,\Omega}\,.$

Finally, thanks to the approximation property of $\mathcal{I}^{SZ}_{h}$ given
in Lemma 4.1, there exists $c>0$, independent of $h$, such that for each
$s\in(0,k+1]$, and for each $\boldsymbol{\varphi}\in\mathbf{H}^{s+1}(\Omega)$,
there holds

$\|\boldsymbol{\varphi}-\mathcal{I}^{SZ}_{h}(\boldsymbol{\varphi})\|_{1,\Omega}\;\leq\;c\,h^{s}\,\|\boldsymbol{\varphi}\|_{s+1,\Omega}\,.$

Therefore, from the Céa estimate (5.14), employing the aforementioned
approximation properties, we can establish the following result.

###### Theorem 5.5

In addition to the hypotheses of Theorems 3.9, 4.8 and 5.4, assume that there
exists $s>0$ such that
$\boldsymbol{t}\in\mathbb{H}^{s}(\Omega),\boldsymbol{\sigma}\in\mathbb{H}^{s}(\Omega),{\mathbf{div}}\,\boldsymbol{\sigma}\in\mathbf{H}^{s}(\Omega),\boldsymbol{u}\in\mathbf{H}^{1+s}(\Omega),\boldsymbol{\varphi}\in\mathbf{H}^{1+s}(\Omega)$.
Then, there exists a positive constant $C$, independent of $h$, such that with
the finite element subspaces defined by (4.1), (4.2), (4.3), (4.4), there
holds

$\|\vec{\boldsymbol{\xi}}-\vec{\boldsymbol{\xi}}_{h}\|_{\mathbb{H}}+\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\;\leq\;C\,h^{\min\\{s,k+1\\}}\Big{\\{}\|\boldsymbol{t}\|_{s,\Omega}+\|\boldsymbol{\sigma}\|_{s,\Omega}+\|{\mathbf{div}}\,\boldsymbol{\sigma}\|_{s,\Omega}+\|\boldsymbol{u}\|_{1+s,\Omega}+\|\boldsymbol{\varphi}\|_{1+s,\Omega}\Big{\\}}.$

##  6 Numerical results

In this section we present three numerical experiments illustrating the
performance of our semi-augmented mixed-primal finite element scheme (4.8),
and confirming the rates of convergence provided by Theorem 5.5. More
precisely, we take the stabilization parameters $\kappa_{1}$, $\kappa_{2}$ and
$\kappa_{3}$ as in (3.24), which satisfies the assumption of Lemma 4.2. In
addition, the zero integral mean condition for tensors in the space (4.2) is
imposed via a real Lagrange multiplier. In turn, the nonlinear algebraic
systems obtained are solved by the fixed-point method with a tolerance of
$10^{-6}$, along with the Newton method for approximate the solution of (4.8a)
in each fixed-point’s iteration. We take as initial guess the solution of a
similar linear problem (in particular, satisfying the boundary conditions for
$\boldsymbol{u}_{h}$ and $\boldsymbol{\varphi}_{h}$). The numerical results
presented below were obtained using a C++ code, where the corresponding linear
systems arising from (4.8a) are solved using the BiCGSTAB method, whereas for
(4.8b) we employ the Conjugate Gradient method as the main solver. Finally, in
all experiments we let ${\boldsymbol{g}}=(0,-1)^{\mathrm{t}}$ be the
gravitational force, and utilizing structure triangulations of the
corresponding domain in 2D. Furthermore, for the first two examples we
consider polynomial degrees $k\in\\{0,1,2\\}$, whereas we only use $k=0$ in
the last example.

We now introduce some additional notation. The individual errors are denoted
by:

$\mathtt{e}(\boldsymbol{t})\ :=\
\|\boldsymbol{t}-\boldsymbol{t}_{h}\|_{0,\Omega}\,,\qquad\mathtt{e}(\boldsymbol{\sigma})\
:=\
\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{\mathbf{div},\Omega}\,,\qquad\mathtt{e}(\boldsymbol{u})\
:=\ \|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1,\Omega}\,,$
$\mathtt{e}(\boldsymbol{\varphi})\ :=\
\|\boldsymbol{\varphi}-\boldsymbol{\varphi}_{h}\|_{1,\Omega}\qquad\text{and}\qquad\mathtt{e}(p)\
:=\ \|p-p_{h}\|_{1,\Omega}\,,$

where, according to (2.6) and (2.7), $p_{h}$ can be computed as:

$p_{h}\ =\
-\dfrac{1}{n}\Big{[}\mathrm{tr}(\boldsymbol{\sigma}_{h}\,+\,(\boldsymbol{u}_{h}\otimes\boldsymbol{u}_{h}))\Big{]}\,+\,\dfrac{1}{n|\Omega|}\|\boldsymbol{u}_{h}\|_{0,\Omega}^{2}\,.$

On the other hand as is usual, we let $\mathtt{r}(\cdot)$ be the experimental
rate of convergence given by

$\mathtt{r}(\cdot)\ :=\
\frac{\mathrm{log}(\mathtt{e}(\cdot)\,/\,\mathtt{e}^{\prime}(\cdot))}{\mathrm{log}(h\,/\,h^{\prime})}\,,$

where $\mathtt{e}$ and $\mathtt{e}^{\prime}$ denote errors computed on two
consecutive meshes of sizes $h$ and $h^{\prime}$, respectively. In addition,
$N$ stands for the total number of degrees of freedom (unknowns) of (4.8),
that is,

$\displaystyle N$ $\displaystyle:=$ $\displaystyle 4\times\\{\text{number of
nodes in
}\mathcal{T}_{h}\\}\,+\,\big{\\{}2(k+1)+4k\big{\\}}\times\\{\text{number of
edges in }\mathcal{T}_{h}\\}$ $\displaystyle+\
\big{\\{}(k+1)(k+2)+2k(k+1)+2k(k-1)\big{\\}}\times\\{\text{number of elements
in }\mathcal{T}_{h}\\}\ +\ 1\,.$

Example 1. We first consider the square $\Omega=(0,1)^{2}$, and set
${\Gamma_{\mathrm{D}}}=\left\\{(s,0),(s,1)\in\text{\rm R}^{2}\,:\,0\leq s\leq
1\right\\}$, ${\Gamma_{\mathrm{N}}}=\Gamma\setminus{\Gamma_{\mathrm{D}}}$,
$\gamma=0.1$, $\nu({\boldsymbol{x}})=(x_{1}^{2}+x_{2}^{2}+1)^{-1}$ (here
$\nu_{1}=1$ and $\nu_{2}=2$), $\boldsymbol{\alpha}=(0.5,1.5)^{\mathrm{t}}$,
the thermal conductivity $\mathbb{K}=2\,\mathbb{I}$, and adequately
manufacture the data so that the exact solution is given by the smooth
functions

$\boldsymbol{u}({\boldsymbol{x}})\ =\ \left(\begin{array}[]{r}-\sin^{2}(2\pi
x_{1})\sin(4\pi x_{2})\\\\[4.30554pt] \sin(4\pi x_{1})\sin^{2}(2\pi
x_{2})\end{array}\right)\,,\qquad p({\boldsymbol{x}})\ =\
\cos(x_{1})\cos(x_{2})\,-\,\sin^{2}(1)\,,$

and

$\boldsymbol{\varphi}({\boldsymbol{x}})\ =\
\left(\begin{array}[]{c}x_{1}x_{2}\\\\[4.30554pt]
\exp(x_{1}+x_{2})\end{array}\right)\,,$

for all ${\boldsymbol{x}}:=(x_{1},x_{2})^{\mathrm{t}}\in\Omega$. In Table 1,
we summarize the convergence history of the finite element scheme (4.8) as
applied to Example 1. We notice there that the rate of convergence
$O(h^{k+1})$ predicted by Theorem 5.5 is attained by all the unknowns.

$k$ | $h$ | $N$ | $\mathtt{e}(\boldsymbol{t})$ | $\mathtt{r}(\boldsymbol{t})$ | $\mathtt{e}(\boldsymbol{\sigma})$ | $\mathtt{r}(\boldsymbol{\sigma})$ | $\mathtt{e}(\boldsymbol{u})$ | $\mathtt{r}(\boldsymbol{u})$ | $\mathtt{e}(\boldsymbol{\varphi})$ | $\mathtt{r}(\boldsymbol{\varphi})$ | $\mathtt{e}(p)$ | $\mathtt{r}(p)$
---|---|---|---|---|---|---|---|---|---|---|---|---
| 0.0404 | 17575 | 7.98e-01 | $--$ | 2.68e+00 | $--$ | 1.37e+00 | $--$ | 4.13e-02 | $--$ | 1.09e-01 | $--$
| 0.0314 | 28895 | 6.05e-01 | 1.10 | 2.09e+00 | 1.00 | 1.05e+00 | 1.05 | 3.20e-02 | 1.02 | 8.53e-02 | 0.98
| 0.0257 | 43015 | 4.87e-01 | 1.09 | 1.71e+00 | 1.00 | 8.51e-01 | 1.05 | 2.61e-02 | 1.01 | 6.99e-02 | 0.99
| 0.0218 | 59935 | 4.07e-01 | 1.07 | 1.45e+00 | 1.00 | 7.16e-01 | 1.04 | 2.20e-02 | 1.01 | 5.92e-02 | 1.00
0 | 0.0189 | 79655 | 3.50e-01 | 1.06 | 1.25e+00 | 1.00 | 6.18e-01 | 1.03 | 1.91e-02 | 1.01 | 5.13e-02 | 1.00
| 0.0129 | 170725 | 2.35e-01 | 1.04 | 8.54e-01 | 1.00 | 4.18e-01 | 1.02 | 1.30e-02 | 1.00 | 3.49e-02 | 1.00
| 0.0094 | 316805 | 1.71e-01 | 1.02 | 6.27e-01 | 1.00 | 3.05e-01 | 1.01 | 9.52e-03 | 1.00 | 2.56e-02 | 1.00
| 0.0071 | 562405 | 1.28e-01 | 1.01 | 4.70e-01 | 1.00 | 2.29e-01 | 1.01 | 7.14e-03 | 1.00 | 1.92e-02 | 1.00
| 0.0057 | 878005 | 1.02e-01 | 1.00 | 3.76e-01 | 1.00 | 1.83e-01 | 1.01 | 5.71e-03 | 1.00 | 1.54e-02 | 1.00
| 0.0404 | 59645 | 5.37e-02 | $--$ | 1.77e-01 | $--$ | 8.81e-02 | $--$ | 1.45e-04 | $--$ | 7.38e-03 | $--$
| 0.0314 | 98285 | 3.20e-02 | 2.07 | 1.07e-01 | 1.99 | 5.28e-02 | 2.04 | 8.48e-05 | 2.12 | 4.53e-03 | 1.94
1 | 0.0257 | 146525 | 2.12e-02 | 2.05 | 7.18e-02 | 2.00 | 3.51e-02 | 2.03 | 5.61e-05 | 2.06 | 3.05e-03 | 1.96
| 0.0218 | 204365 | 1.51e-02 | 2.04 | 5.14e-02 | 2.00 | 2.51e-02 | 2.02 | 4.00e-05 | 2.03 | 2.20e-03 | 1.97
| 0.0189 | 271805 | 1.13e-02 | 2.03 | 3.86e-02 | 2.00 | 1.88e-02 | 2.02 | 3.00e-05 | 2.01 | 1.66e-03 | 1.98
| 0.0129 | 583445 | 5.19e-03 | 2.02 | 1.80e-02 | 2.00 | 8.70e-03 | 2.01 | 1.39e-05 | 2.01 | 7.75e-04 | 1.98
| 0.0404 | 126215 | 2.70e-03 | $--$ | 8.82e-03 | $--$ | 3.76e-03 | $--$ | 3.25e-07 | $--$ | 3.45e-04 | $--$
| 0.0314 | 208175 | 1.26e-03 | 3.04 | 4.15e-03 | 2.99 | 1.76e-03 | 3.03 | 1.47e-07 | 3.17 | 1.63e-04 | 2.99
2 | 0.0257 | 310535 | 6.87e-04 | 3.03 | 2.28e-03 | 3.00 | 9.58e-04 | 3.02 | 7.96e-08 | 3.05 | 8.92e-05 | 2.99
| 0.0218 | 433295 | 4.15e-04 | 3.01 | 1.36e-03 | 3.08 | 5.79e-04 | 3.02 | 4.79e-08 | 3.04 | 5.40e-05 | 3.01
| 0.0189 | 576455 | 2.70e-04 | 3.00 | 8.83e-04 | 3.03 | 3.76e-04 | 3.01 | 3.10e-08 | 3.03 | 3.51e-05 | 3.00
Table 1: History of convergence for Example 1.

Example 2. Next, we adapt [11, Example 3], and consider the $L$-shaped domain
$\Omega=(-1,1)^{2}\,\setminus\,[0,1]^{2}$, and set
${\Gamma_{\mathrm{N}}}=\left\\{(s,0),(0,s)\in\text{\rm R}^{2}\,:\,0\leq s\leq
1\right\\}$, ${\Gamma_{\mathrm{D}}}=\Gamma\setminus\bar{\Gamma}_{N}$,
$\gamma=10^{-3}$, $\nu({\boldsymbol{x}})=1+\exp(-x_{1}^{2})$ (once again
$\nu_{1}=1$ and $\nu_{2}=2$), $\boldsymbol{\alpha}=(1,0.5)^{\mathrm{t}}$,
$\mathbb{K}=\left(\begin{smallmatrix}1&\;0\\\\[2.15277pt]
0&\;2\end{smallmatrix}\right)$, and adequately manufacture the data so that
the exact solution is given by the smooth functions

$\boldsymbol{u}({\boldsymbol{x}})\ =\
\left(\begin{array}[]{r}x_{2}^{2}\\\\[4.30554pt]
-x_{1}^{2}\end{array}\right)\,,\qquad p({\boldsymbol{x}})\ =\
(x_{1}^{2}\,+\,x_{2}^{2})^{1/3}\,-\,p_{0}\,,\qquad\text{and}\qquad\boldsymbol{\varphi}({\boldsymbol{x}})\
=\ \left(\begin{array}[]{c}\exp(-x_{1}^{2}-x_{2}^{2})\\\\[4.30554pt]
\exp(-x_{1}x_{2})\end{array}\right)\,,$

for all ${\boldsymbol{x}}:=(x_{1},x_{2})^{\mathrm{t}}\in\Omega$, where
$p_{0}\in\mathrm{R}$ is such that $\int_{\Omega}p=0$ holds
($p_{0}\approx\text{8.211056552903396e-01}$). In addition, we remark here that
the partial derivatives of $p$, and hence, in particular
${\mathbf{div}}(\boldsymbol{\sigma})$, are singular at the origin. Indeed,
according to the power $1/3$, there holds
$\boldsymbol{\sigma}\in\mathbb{H}^{5/3-\varepsilon}(\Omega)$ and
${\mathbf{div}}(\boldsymbol{\sigma})\in\textbf{H}^{2/3-\varepsilon}(\Omega)$
for each $\varepsilon>0$. In fact, in Table 2 we present the corresponding
convergence history of Example 2, where, as predicted in advance, we note that
the orders $O(h^{\min\\{k+1,5/3\\}})$ and $O(h^{2/3})$ are attained by
$(\boldsymbol{t}_{h},\boldsymbol{u}_{h})$ and $\boldsymbol{\sigma}_{h}$,
respectively. Once again, the rate of convergence predicted by Theorem 5.5 is
attained by all the unknowns, except for the variable
$\boldsymbol{\varphi}_{h}$ that preserves $O(h^{k+1})$. The foregoing
phenomenon could be a special feature of this example. Furthermore, the
results in Example 2 suggest that our approach should certainly be
strengthened with the further incorporation of an adaptive strategy based on a
suitable a-posteriori error estimates. This issue will also be addressed in a
forthcoming paper.

$k$ | $h$ | $N$ | $\mathtt{e}(\boldsymbol{t})$ | $\mathtt{r}(\boldsymbol{t})$ | $\mathtt{e}(\boldsymbol{\sigma})$ | $\mathtt{r}(\boldsymbol{\sigma})$ | $\mathtt{e}(\boldsymbol{u})$ | $\mathtt{r}(\boldsymbol{u})$ | $\mathtt{e}(\boldsymbol{\varphi})$ | $\mathtt{r}(\boldsymbol{\varphi})$ | $\mathtt{e}(p)$ | $\mathtt{r}(p)$
---|---|---|---|---|---|---|---|---|---|---|---|---