Poroelasticity with stress-dependent diffusionGómez-Vargas, Mardal, Ruiz-
Baier, Vinje

# Twofold saddle-point formulation of Biot poroelasticity
with stress-dependent diffusion††thanks: Updated: .BGV acknowledges support
from Vicerrectoría de Investigación, through the project 540-C0-202, Sede de
Occidente, Universidad de Costa Rica. KAM acknowledges support from the
Research Council of Norway, grants 300305 and 301013\. RRB acknowledges
support from the Monash Mathematics Research Fund S05802-3951284 and from the
Ministry of Science and Higher Education of the Russian Federation within the
framework of state support for the creation and development of World-Class
Research Centers “Digital biodesign and personalised healthcare” No.
075-15-2020-926.

Bryan Gómez-Vargas Sección de Matemática, Sede de Occidente, Universidad de
Costa Rica, San Ramón, Alajuela, Costa Rica<EMAIL_ADDRESS>Kent-André Mardal Department of Mathematics, Division of Mechanics, University
of Oslo, Norway; and Simula Research Laboratory, 1325 Lysaker, Norway ().
<EMAIL_ADDRESS>Ricardo Ruiz-Baier School of Mathematics, Monash
University, 9 Rainforest Walk, Melbourne 3800 VIC, Australia; and Institute of
Computer Science and Mathematical Modelling, Sechenov University, Moscow,
Russian Federation; and Universidad Adventista de Chile, Casilla 7-D, Chillán,
Chile<EMAIL_ADDRESS>Vegard Vinje Department for Numerical
Analysis and Scientific Computing, Simula Research Laboratory, Lysaker, Norway
<EMAIL_ADDRESS>
###### Abstract

We introduce a stress/total-pressure formulation for poroelasticity that
includes the coupling with steady nonlinear diffusion modified by stress. The
nonlinear problem is written in mixed-primal form, coupling a perturbed
twofold saddle-point system with an elliptic problem. The continuous
formulation is analysed in the framework of abstract fixed-point theory and
Fredholm alternative for compact operators. A mixed finite element method is
proposed and its stability and convergence analysis are carried out. We also
include a few illustrative numerical tests. The resulting model can be used to
study waste removal in the brain parenchyma, where diffusion of a tracer alone
or a combination of advection and diffusion are not sufficient to explain the
alterations in rates of filtration observed in porous media samples.

###### keywords:

Stress-altered diffusion; Poroelasticity; Perturbed saddle-point; Mixed finite
elements; Brain multiphysics.

65N60, 92C50.

##  1 Introduction

Poroelastic structures are found in many applications of industrial and
scientific relevance. Examples include the interaction between soft permeable
tissue and blood flow, or the study of biofilm growth and distribution near
fluids [39]. We are concerned with one particular application involving the
transport of cerebrospinal fluid (CSF) and tracer within the brain parenchyma,
and how this can contribute to better explain the mechanisms that permit
solutes from the brain interstitial space to move. These processes are key in
eliminating waste from the brain and in maintaining a healthy function of the
nervous system (see, for instance, the review [36]), and which address
important questions in the context of the overall glymphatic function [27].

Detailed maps of permeability and simulations using realistic hydrostatic
pressure gradients and image-based reconstruction of intersitial space,
suggest that diffusive effects in combination with advection along vasculature
dominate over bulk flow [25]. However it appears that (even heterogeneous)
Fickian diffusion and convective effects only are not sufficient to explain
the main features of the physiological propagation of solutes along the
various elements of the brain structure [11, 42]. We thus reconsider how other
changes related to the mechanical state, such as the generation and nonuniform
localisation of stresses in the brain, could impact the filtration properties
at the spatio-temporal scales of interest for the application at hand.
Pressure, strain, and stress changes can affect the tissue microstructure by
rearranging the local distribution of compliant arrays, opening and closing
pores, as well as generating preferential paths for the accumulation of large
solvent molecules [45]. The combination of such effects can result in
variations of averaged coefficients usually employed in classical models, such
as diffusion (but also other parameters such as bulk moduli, relative
permeability, and Biot-Willis coefficient [34]). For this effect we propose a
simplified phenomenological model for the steady interaction between
poroelastic stress and diffusion of solutes (such as gadobutrol). The coupling
of poromechanics and diffusion is a fundamental aspect not only for
investigating the mechanical behaviour of the brain function, but also in the
description and design of other complex systems such as tissue engineering
scaffolds [10].

Linear models of poroelasticity seem to be sufficient to capture the most
important strain patterns when coupled to, e.g., immune system equations
modelling the blood-brain barrier in the brain (see the comparison in [30]).
In the present model the linear poroelasticity equations are formulated using
the stress as an unknown. This is of particular importance as this field will
be central to the coupling with the diffusion process. The recent literature
contains various formulations for Biot equations using poroelastic stress. For
instance, including weakly symmetric four-field schemes [4], adaptive mixed
formulations [15], multipoint stress-flux mixed methods [1], and domain
decomposition approaches [28]. In contrast, here we employ a double saddle-
point formulation designed to solve for total poroelastic stress,
displacement, rotation, total pressure (drawing inspiration from [31, 35] and
generalising these formulations to include poroelastic and active stress), and
fluid pressure. The Biot equations are then coupled with a nonlinear diffusion
equation for the CSF tracer, and the discretisation proceeds with classical
mixed-primal finite elements. Mixed-primal and fully mixed methods for stress-
modified diffusion problems have been investigated in [19, 20], and their
analysis hinges on techniques using a separation of elastostatics and
diffusion equations by means of a fixed-point scheme that requires additional
assumptions on regularity. This is precisely the approach followed here. These
models, as well as modification of transport properties in other soft tissues
such as cardiac muscle [38] allow to recover enhanced anisotropy and
heterogeneity of stress-assisted diffusion. However, here the coupling with
total poroelastic stress is used to _reduce_ or hinder diffusion patterns
which have the potential to produce a much better match with experimental
studies that indicate a decrease of effective diffusion by at least an order
of magnitude with respect to the base state observed during sleep. One of the
main challenges of the work lies in deriving a variational formulation
involving total poroelastic stress, and whose analysis does not need
additional conditions on the Lamé parameters. The resulting form has the
structure of a double saddle-point problem, where terms that explicitly
involve $\lambda_{s}$, appear. Exploiting the compactness of these terms, an
applying the Fredholm alternative together with a new extended version of the
Babuška-Brezzi theory, allow us to assert unique solvability robustly with
respect to $\lambda_{s}$.

The content of this paper has been laid out as follows. Section 2 provides
details of the model, describing the components of the balance equations and
stating boundary conditions. Section 3 is devoted to the analysis of
solvability of the weak form, using fixed-point arguments in conjunction with
the Fredholm alternative and the Babuška-Brezzi theory. We address in Section
4 the solvability and stability analysis of the discrete problem. A priori
error estimates are derived in Section 5, and in Section 6 we collect
computational results, consisting in verification of convergence and
simulation of different cases on simple geometries.

##  2 Governing equations

Let us consider a Lipschitz polyhedral domain $\Omega\subset\mathbb{R}^{d}$,
$d=2,3$ occupied by the brain parenchyma, here modelled as a mixture
consisting of an elastic structure and an interconnected porous space fully
saturated by interstitial fluid. We will denote by $\boldsymbol{n}$ the
outward unit normal vector on the boundary $\partial\Omega$. We also define
the boundaries $\Gamma$ and $\Sigma$ where different types of boundary
conditions will be applied. We assume that the parenchyma is isotropic, that
the deformation process is stationary, and that its constitutive stress-strain
relation is linear. As usual in the framework of linear poroelasticity, the
solid-fluid mixture is described in terms of the solid displacement
$\boldsymbol{u}:\Omega\to\mathbb{R}^{d}$ accounting for the movement of the
material from its initial undeformed configuration, and the pore fluid
pressure $p:\Omega\to\mathbb{R}$ (which, through Darcy’s law, also defines the
specific discharge or filtration velocity). The Biot system states momentum
and mass balances

(2.1)
$-\operatorname*{div}(\boldsymbol{\sigma})=\rho_{s}\boldsymbol{f}\quad\mathrm{and}\quad
c_{0}p+\alpha\operatorname*{div}(\boldsymbol{u})-\operatorname*{div}\left(\frac{\mathbf{\kappa}}{\mu_{f}}\nabla
p-\rho_{f}\boldsymbol{g}\right)=m\qquad\text{in $\Omega$,}$

respectively, where $\boldsymbol{f}:\Omega\to\mathbb{R}^{d}$ is the vector
field of body loads and $m:\Omega\to\mathbb{R}$ is a sources/sink of fluid. We
also consider the presence of a CSF tracer within the poroelastic parenchymal
domain. We denote its concentration by $\omega:\Omega\to\mathbb{R}$ and its
movement in the parenchyma is governed by

(2.2)
$\phi\omega-\operatorname*{div}(D(\boldsymbol{\sigma})\nabla\omega)=\phi\ell\qquad\text{in
$\Omega$},$

where
$\boldsymbol{\sigma}=2\mu_{s}\boldsymbol{\varepsilon}(\boldsymbol{u})+(\lambda_{s}\mathrm{div}(\boldsymbol{u})-\alpha
p-\beta\omega)\mathbb{I}$ is the poroelastic Cauchy stress tensor with
$\boldsymbol{\varepsilon}(\boldsymbol{u}):=\frac{1}{2}(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\tt
t})$ denoting the infinitesimal strain tensor, the Lamé constants
$\mu_{s},\lambda_{s}>0$ for homogeneous and isotropic materials, the Biot-
Willis constant $0<\alpha\leq 1$, the $d\times d$ identity matrix
$\mathbb{I}$, and $\beta$ a scaling of active stress that indicates a two-way
coupling between diffusion and motion. The right-hand side
$\ell:\Omega\to\mathbb{R}$ is a non-negative drainage or source coefficient
that encodes absorption from the lymph nodes and/or from capillaries [12], and
$\phi$ is porosity. In (2.1) $c_{0}$ is the mass storativity coefficient, and
$\kappa$ is the permeability tensor (symmetric, uniformly positive definite,
and bounded). The tensor function $D:\mathbb{R}^{d\times
d}\rightarrow\mathbb{R}^{d\times d}$ is a stress-dependent diffusivity
accounting for an altered diffusion acting in the poroelastic domain. This
term may assume the simple form

(2.3)
$D(\boldsymbol{\sigma})=\eta_{0}D_{0}+\eta_{1}\exp(-\eta_{1}\operatorname*{tr}\boldsymbol{\sigma}),$

with $D_{0}$ being the effective diffusion of typical solutes such as
potassium ions, and $\eta_{0},\eta_{1}$ modulating the intensity of the
stress-altered diffusion. Other constitutive equations (including power-law,
polynomial, and anisotropic forms) can be found in, e.g., [19], whereas
pressure-dependent diffusion models can be found in [24]. Alternative
constitutive relations could be given by a dependence on the volumetric part
of the effective poroelastic stress, or on the porosity (which in the linear
regime is the total amount of fluid) such as
$D=D_{0}-\eta(c_{0}p+\alpha\operatorname*{div}\boldsymbol{u})$. On the
boundary $\Gamma$ we prescribe tracer concentration, displacement and normal
specific discharge; whereas on $\Sigma$ we fix zero normal total stress, zero
flux for the tracer, and prescribed fluid pressure.

Besides the main variables we also employ the total pressure (the volumetric
contributions to the poroelastic Cauchy stress) $\widetilde{p}:=\alpha
p-\lambda_{s}\,\mathrm{div}\,\boldsymbol{u}$, the poroelastic stress tensor
$\boldsymbol{\sigma}=2\mu_{s}\boldsymbol{\varepsilon}(\boldsymbol{u})-\widetilde{p}\mathbb{I}-\beta\omega\mathbb{I}$,
and the tensor $\boldsymbol{\gamma}$ (rotation Lagrange multiplier assisting
the weakly imposition of the symmetry of $\boldsymbol{\sigma}$), and the
strain
$\boldsymbol{t}:=\boldsymbol{\varepsilon}(\boldsymbol{u})=\nabla\boldsymbol{u}-\boldsymbol{\gamma}$.
Therefore the governing equations consist in finding $\boldsymbol{\sigma}$,
$\boldsymbol{\rho}$, $\boldsymbol{t}$, $\tilde{p}$, $p$, and $\omega$,
satisfying

(2.4a)
$\displaystyle\boldsymbol{\sigma}=2\mu_{s}\boldsymbol{t}-\tilde{p}\mathbb{I}-\beta\omega\mathbb{I}\quad\text{in
}\Omega,\quad\boldsymbol{\sigma}=\boldsymbol{\sigma}^{\tt t}\quad\text{in
}\Omega,\quad-\operatorname*{\mathbf{div}}\boldsymbol{\sigma}=\rho_{s}\boldsymbol{f}\quad\text{in
}\Omega,\quad\boldsymbol{t}=\boldsymbol{\varepsilon}(\boldsymbol{u})\quad\text{in
}\Omega,$ (2.4b)
$\displaystyle\quad\boldsymbol{\gamma}=\nabla\boldsymbol{u}-\boldsymbol{t}\quad\text{in
}\Omega,\quad\left(c_{0}+\frac{\alpha^{2}}{\lambda_{s}}\right)p-\frac{\alpha}{\lambda_{s}}\tilde{p}-\operatorname*{div}\left(\frac{\mathbf{\kappa}}{\mu_{f}}\nabla
p-\rho_{f}\boldsymbol{g}\right)=m\quad\text{in }\Omega,$ (2.4c)
$\displaystyle\quad\frac{1}{\lambda_{s}}\tilde{p}+\operatorname*{tr}\boldsymbol{t}-\frac{\alpha}{\lambda_{s}}p=0\quad\text{in
}\Omega,\quad\boldsymbol{u}=\boldsymbol{0},\quad\frac{\kappa}{\mu_{f}}\nabla
p\cdot\boldsymbol{n}=0\quad\text{on
}\Gamma,\quad\boldsymbol{\sigma}\boldsymbol{n}=\boldsymbol{0},\quad
p=0\quad\text{on }\Sigma,$ (2.4d)
$\displaystyle\phi\omega-\operatorname*{div}(D(\boldsymbol{\sigma})\nabla\omega)=\phi\ell\quad\text{in
}\Omega,\quad\omega=0\quad\text{on }\Gamma,\quad\text{and}\quad
D(\boldsymbol{\sigma})\nabla\omega\cdot\boldsymbol{n}=0\quad\text{on }\Sigma.$

The second and third equations in (2.4a) state, respectively, the conservation
of angular and linear momentum, while the last equation in (2.4a) and the
first equality in (2.4c) encode other constitutive relations. Homogeneous
boundary data have been assumed only for sake of conciseness, but they can be
generalised to the non-homogeneous case by using classical lifting arguments.
An advantage of using this modified form of the system is that all equations
are robust with respect to the Lamé parameters.

##  3 Continuous mixed-primal formulation and well-posedness analysis

Let us define the tensor functional spaces for stress and rotations

$\mathbb{H}_{\Sigma}(\mathbf{div},\Omega):=\\{\boldsymbol{\tau}\in\mathbb{H}({\mathbf{div}},\Omega):\boldsymbol{\tau}\boldsymbol{n}=\boldsymbol{0}\
\mathrm{on}\,\Sigma\\},\
\mathbb{L}^{2}_{\mathrm{skew}}(\Omega):=\\{\boldsymbol{\eta}\in\mathbb{L}^{2}(\Omega):\boldsymbol{\eta}+\boldsymbol{\eta}^{\mathrm{t}}\\!=\boldsymbol{0}\\}.$

A weak formulation for the nonlinearly coupled problem (2.4a)-(2.4d) can be
derived as usual, and it consists in finding
$\bigl{(}\boldsymbol{t},(\boldsymbol{\sigma},\tilde{p}),(\boldsymbol{u},\boldsymbol{\gamma}),p,\omega\bigr{)}\in\mathbb{L}^{2}(\Omega)\times[\mathbb{H}_{\Sigma}(\mathbf{div},\Omega)\times\mathrm{L}^{2}(\Omega)]\times[{\mathbf{L}}^{2}(\Omega)\times\mathbb{L}^{2}_{\mathrm{skew}}(\Omega)]\times\mathrm{H}^{1}_{\Sigma}(\Omega)\times\mathrm{H}^{1}_{\Gamma}(\Omega)$
such that

(3.1a)
$\displaystyle-\bigl{(}c_{0}+\frac{\alpha^{2}}{\lambda_{s}}\bigr{)}\int_{\Omega}pq+\frac{\alpha}{\lambda_{s}}\int_{\Omega}\widetilde{p}q-\frac{1}{\mu_{f}}\int_{\Omega}\boldsymbol{\kappa}\nabla
p\cdot\nabla
q=\rho_{f}\int_{\Omega}\boldsymbol{\kappa}\boldsymbol{g}\cdot\nabla
q-\int_{\Omega}mq\quad\forall\,q\in\mathrm{H}^{1}_{\Sigma}(\Omega),$ (3.1b)
$\displaystyle-\int_{\Omega}\boldsymbol{t}\colon\boldsymbol{\tau}-\int_{\Omega}\boldsymbol{u}\cdot\mathbf{div}\,\boldsymbol{\tau}-\int_{\Omega}\boldsymbol{\gamma}\colon\boldsymbol{\tau}=0\;\;\;\forall\boldsymbol{\tau}\\!\in\mathbb{H}_{\Sigma}(\mathbf{div};\Omega),\quad-\int_{\Omega}\boldsymbol{v}\cdot\mathbf{div}\,\boldsymbol{\sigma}=\rho_{s}\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}\quad\forall\,\boldsymbol{v}\in\mathbf{L}^{2}(\Omega),$
(3.1c)
$\displaystyle-\frac{1}{\lambda_{s}}\int_{\Omega}\widetilde{p}\,\widetilde{q}-\int_{\Omega}\mathrm{tr}(\boldsymbol{t})\widetilde{q}+\frac{\alpha}{\lambda_{s}}\int_{\Omega}p\widetilde{q}=0\quad\forall\,\widetilde{q}\in\mathrm{L}^{2}(\Omega),\quad-\int_{\Omega}\boldsymbol{\sigma}\colon\boldsymbol{\eta}=0\quad\forall\,\boldsymbol{\eta}\in\mathbb{L}_{\mathrm{skew}}^{2}(\Omega),$
(3.1d) $\displaystyle
2\mu_{s}\int_{\Omega}\boldsymbol{t}\colon\boldsymbol{r}-\int_{\Omega}\boldsymbol{\sigma}\colon\boldsymbol{r}-\int_{\Omega}\widetilde{p}\,\mathrm{tr}(\boldsymbol{r})=\beta\int_{\Omega}\omega\mathrm{tr}(\boldsymbol{r})\quad\forall\,\boldsymbol{r}\in\mathbb{L}^{2}(\Omega),$
(3.1e)
$\displaystyle\int_{\Omega}\phi\omega\theta+\int_{\Omega}D(\boldsymbol{\sigma})\nabla\omega\cdot\nabla\theta=\int_{\Omega}\ell\phi\theta\quad\forall\,\theta\in\mathrm{H}^{1}_{\Gamma}(\Omega).$

We further assume that $\boldsymbol{\kappa}$ is uniformly bounded and positive
definite, and that $D(\boldsymbol{\sigma})$ is of class $C^{1}$, uniformly
positive definite, Lipschitz continuous and that satisfies a particular
bounding property: There exist positive constants $\kappa_{1},\kappa_{2}$, and
$d_{1},d_{2},d_{3}$, such that

(3.2)
$\kappa_{1}|\boldsymbol{v}|^{2}\leq\boldsymbol{v}^{\mathrm{t}}\kappa(\cdot)\boldsymbol{v}\leq\kappa_{2}|\boldsymbol{v}|^{2},\quad\mathrm{and}\quad
d_{1}|\boldsymbol{v}|^{2}\leq\boldsymbol{v}^{\mathrm{t}}D(\cdot)\boldsymbol{v}\leq
d_{2}|\boldsymbol{v}|^{2},\quad|D(\boldsymbol{\sigma})-D(\boldsymbol{\zeta})|\leq
d_{3}|\boldsymbol{\sigma}-\boldsymbol{\zeta}|,$

for all $\boldsymbol{v}\in\mathbb{R}^{d}$, and for all
$\boldsymbol{\sigma},\boldsymbol{\zeta}\in\mathbb{R}^{d\times d}$. Next, we
regroup unknowns and spaces in (3.1a)-(3.1d) as follows:

$\vec{\boldsymbol{\sigma}}:=(\boldsymbol{\sigma},\widetilde{p}),\;\;\vec{\boldsymbol{\tau}}:=(\boldsymbol{\tau},\widetilde{q}),\;\;\vec{\boldsymbol{u}}:=(\boldsymbol{u},\boldsymbol{\gamma}),\quad\vec{\boldsymbol{v}}:=(\boldsymbol{v},\boldsymbol{\eta}),\;\;\mathrm{and}\quad\mathbf{X}:=\mathbb{H}_{\Sigma}(\mathbf{div};\Omega)\times\mathrm{L}^{2}(\Omega),\;\;\mathbf{M}:=\mathbf{L}^{2}(\Omega)\times\mathbb{L}^{2}_{\mathrm{skew}}(\Omega).$

Then, given
$(H_{\omega},F,G,J)\in\mathbb{L}^{2}(\Omega)^{\prime}\times\mathbf{M}^{\prime}\times\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}\times\mathrm{H}^{1}_{\Gamma}(\Omega)^{\prime}$,
we seek
$(\boldsymbol{t},\vec{\boldsymbol{\sigma}},\vec{\boldsymbol{u}},p,\omega)\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)\times\mathrm{H}^{1}_{\Gamma}(\Omega)$
such that

(3.3a) $\displaystyle[A(\boldsymbol{t}),\boldsymbol{r}]$ $\displaystyle+\quad$
$\displaystyle[B^{*}(\vec{\boldsymbol{\sigma}}),\boldsymbol{r}]$
$\displaystyle=$ $\displaystyle\;\;[H_{\omega},\boldsymbol{r}],$ (3.3b)
$\displaystyle[B(\boldsymbol{t}),\vec{\boldsymbol{\tau}}]\;\;\,$
$\displaystyle-\quad$
$\displaystyle[C(\vec{\boldsymbol{\sigma}}),\vec{\boldsymbol{\tau}}]\;\;\,$
$\displaystyle+$
$\displaystyle\;\;[B^{*}_{1}(\vec{\boldsymbol{u}}),\vec{\boldsymbol{\tau}}]\;\;$
$\displaystyle+$ $\displaystyle\;\;[B_{2}^{*}(p),\vec{\boldsymbol{\tau}}]\;$
$\displaystyle=$ $\displaystyle\;\;\mathrm{O},$ (3.3c)
$\displaystyle[B_{1}(\vec{\boldsymbol{\sigma}}),\vec{\boldsymbol{v}}]$
$\displaystyle=$ $\displaystyle\;\;[F,\vec{\boldsymbol{v}}],$ (3.3d)
$\displaystyle[B_{2}(\vec{\boldsymbol{\sigma}}),q]$ $\displaystyle-$
$\displaystyle\;\;[D(p),q]$ $\displaystyle=$ $\displaystyle\;\;[G,q],$ (3.3e)
$\displaystyle[\mathscr{A}_{\boldsymbol{\sigma}}(\omega),\theta]$
$\displaystyle=$ $\displaystyle\;\;[J,\theta],$

for all
$(\boldsymbol{r},\vec{\boldsymbol{\tau}},\vec{\boldsymbol{v}},q,\theta)\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)\times\mathrm{H}^{1}_{\Gamma}(\Omega)$,
where the linear bounded operators
$A:\mathbb{L}^{2}(\Omega)\rightarrow\mathbb{L}^{2}(\Omega)^{\prime}$,
$C:\mathbf{X}\rightarrow\mathbf{X}^{\prime}$,
$B:\mathbb{L}^{2}(\Omega)\rightarrow\mathbf{X}^{\prime}$,
$B_{1}:\mathbf{X}\rightarrow\mathbf{M}^{\prime}$,
$B_{2}:\mathbf{X}\rightarrow\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime},$
$D:\mathrm{H}^{1}_{\Sigma}(\Omega)\rightarrow\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}$,
$B^{*}:\mathbf{X}\rightarrow\mathbb{L}^{2}(\Omega)^{\prime}$,
$B_{1}^{*}:\mathbf{M}\rightarrow\mathbf{X}^{\prime}$,
$B_{2}^{*}:\mathrm{H}^{1}_{\Sigma}(\Omega)\rightarrow\mathbf{X}^{\prime}$, and
$\mathscr{A}_{\boldsymbol{\sigma}}:\mathrm{H}^{1}_{\Gamma}(\Omega)\rightarrow\mathrm{H}^{1}_{\Gamma}(\Omega)^{\prime}$
are defined as

(3.4a)
$\displaystyle[A({\boldsymbol{t}}),{\boldsymbol{r}}]:=2\mu_{s}\int_{\Omega}\boldsymbol{t}\colon\boldsymbol{r},\quad$
$\displaystyle[B({\boldsymbol{r}}),\vec{\boldsymbol{\tau}}]:=-\int_{\Omega}\boldsymbol{\tau}\colon\boldsymbol{r}-\int_{\Omega}\widetilde{q}\,\mathrm{tr}(\boldsymbol{r}),\quad[B_{1}(\vec{\boldsymbol{\tau}}),\vec{\boldsymbol{v}}]:=-\int_{\Omega}\boldsymbol{v}\cdot\mathbf{div}\,\boldsymbol{\tau}-\int_{\Omega}\boldsymbol{\eta}\colon\boldsymbol{\tau},$
(3.4b)
$\displaystyle[C(\vec{\boldsymbol{\sigma}}),\vec{\boldsymbol{\tau}}]:=\frac{1}{\lambda_{s}}\int_{\Omega}\widetilde{p}\,\widetilde{q},$
$\displaystyle\quad[B_{2}(\vec{\boldsymbol{\tau}}),q]:=\frac{\alpha}{\lambda_{s}}\int_{\Omega}\widetilde{q}q,\quad[D(p),q]:=\left(c_{0}+\frac{\alpha^{2}}{\lambda_{s}}\right)\int_{\Omega}pq+\frac{1}{\mu_{f}}\int_{\Omega}\boldsymbol{\kappa}\nabla
p\cdot\nabla q,$ (3.4c)
$\displaystyle[\mathscr{A}_{\boldsymbol{\sigma}}(\omega),\theta]:=$
$\displaystyle\int_{\Omega}\phi\omega\theta+\int_{\Omega}D(\boldsymbol{\sigma})\nabla\omega\cdot\nabla\theta,\quad[H_{\omega},\boldsymbol{r}]:=\beta\int_{\Omega}\omega\mathrm{tr}(\boldsymbol{r}),$
(3.4d) $\displaystyle[F,\vec{\boldsymbol{v}}]:=$
$\displaystyle\rho_{s}\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v},\quad[G,q]:=-\int_{\Omega}mq+\rho_{f}\int_{\Omega}\boldsymbol{\kappa}\boldsymbol{g}\cdot\nabla
q,\quad[J,\theta]:=\int_{\Omega}\ell\phi\theta,$

for each
${\boldsymbol{t}},{\boldsymbol{r}}\in\mathbb{L}^{2}(\Omega),\;\vec{\boldsymbol{\sigma}},\vec{\boldsymbol{\tau}}\in\mathbf{X}$,
$\vec{\boldsymbol{u}},\vec{\boldsymbol{v}}\in\mathbf{M},$
$p,q\in\mathrm{H}^{1}_{\Sigma}(\Omega)$ and
$\omega,\theta\in\mathrm{H}^{1}_{\Gamma}(\Omega).$ The symbol $\mathrm{O}$
stands here for the null functional, and $[\cdot,\cdot]$ denotes the duality
pairing induced by the operators and functionals.

###  3.1 Fixed-point-strategy

In this section, we utilise a fixed-point strategy to prove that (3.3a)-(3.3e)
is well-posed. Let
$\textbf{S}:\mathrm{H}_{\Gamma}^{1}(\Omega)\to\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$
be the operator defined by

$\textbf{S}(\omega):=(\textbf{S}_{1}(\omega),(\textbf{S}_{2}(\omega),\textbf{S}_{3}(\omega)),(\textbf{S}_{4}(\omega),\textbf{S}_{5}(\omega)),\textbf{S}_{6}(\omega))=(\boldsymbol{t},(\boldsymbol{\sigma},\widetilde{p}),(\boldsymbol{u},\boldsymbol{\gamma}),p)\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega),$

where
$(\boldsymbol{t},(\boldsymbol{\sigma},\widetilde{p}),(\boldsymbol{u},\boldsymbol{\gamma}),p)$
is the unique solution of (3.3a)-(3.3d) with $\omega$ given. In turn, we let
${\widetilde{\mathbf{S}}}:\mathbb{H}_{\Sigma}(\mathbf{div},\Omega)\to\mathrm{H}^{1}_{\Gamma}(\Omega)$
be the operator

${\widetilde{\mathbf{S}}}(\boldsymbol{\sigma})\,:=\,\omega\quad\forall\,\boldsymbol{\sigma}\,\in\,\mathbb{H}_{\Sigma}(\mathbf{div},\Omega),$

where $\omega$ is the unique solution of (3.3e) with $\boldsymbol{\sigma}$
given. Then, we define the operator
$\mbox{{{T}}}:\mathrm{H}^{1}_{\Gamma}(\Omega)\to\mathrm{H}^{1}_{\Gamma}(\Omega)$
as

$\mbox{{{T}}}(\omega):=\widetilde{\mathbf{S}}(\mathbf{S}_{2}(\omega))\quad\forall\,\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega),$

and realise that solving (3.3a)-(3.3e) is equivalent to finding
$\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$ such that

(3.5) $\mbox{{{T}}}(\omega)=\omega.$

###  3.2 Well-posedness of the uncoupled problems

In this section, we analyse the solvability of problems defining $\mathbf{S}$
and $\widetilde{\mathbf{S}}$. In view of the analysis of Section 5 and Lemma
3.15 below, the well-posedness and stability of the uncoupled problem related
with the operator S for a given $\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$,
will be developed considering generic functionals
$\widetilde{H}_{\omega},\widetilde{F}_{1},\widetilde{F}$ and $\widetilde{G}$,
instead of $H_{\omega},\mathrm{O},F$ and $G$, respectively, in (3.3a)-(3.3d).
We begin by investigating the well-posedness of this uncoupled problem. First,
we recall from [7] the decomposition

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

That is, for each $\boldsymbol{\tau}\in\mathbb{H}(\mathbf{div},\Omega)$ there
exist unique

(3.6)
$\boldsymbol{\tau}_{0}:=\boldsymbol{\tau}-{\displaystyle\left\\{\frac{1}{d|\Omega|}\int_{\Omega}\mathrm{tr}(\boldsymbol{\tau})\right\\}}\mathbb{I}\in\mathbb{H}_{0}(\mathbf{div},\Omega)\quad\text{and}\quad\hat{d}:=\frac{1}{d|\Omega|}\int_{\Omega}\mathrm{tr}(\boldsymbol{\tau})\in\mathbb{R},$

such that $\boldsymbol{\tau}=\boldsymbol{\tau}_{0}+\hat{d}\mathbb{I}$.
Moreover, we recall the following results which will be useful in the
forthcoming analysis. We remit to [18, Lemma 2.3], and [18, Lemma 2.4],
respectively, for further details.

######  Lemma 3.1.

There exists $c_{1}>0$, depending only on $\Omega$, such that

(3.7)
$c_{1}\|\boldsymbol{\tau}_{0}\|^{2}_{\mathbb{L}^{2}(\Omega)}\leq\|\boldsymbol{\tau}^{\mathrm{d}}\|^{2}_{\mathbb{L}^{2}(\Omega)}+\|\mathbf{div}(\boldsymbol{\tau})\|^{2}_{\mathbb{L}^{2}(\Omega)}\quad\forall\,\boldsymbol{\tau}\in\mathbb{H}(\mathbf{div},\Omega).$

######  Lemma 3.2.

There exists $c_{2}>0$, depending only on $\Sigma$ and $\Omega$, such that

(3.8)
$c_{2}\|\boldsymbol{\tau}\|^{2}_{\mathbb{H}(\mathbf{div},\Omega)}\leq\|\boldsymbol{\tau}_{0}\|^{2}_{\mathbb{H}(\mathbf{div},\Omega)}\quad\forall\,\boldsymbol{\tau}\in\mathbb{H}_{\Sigma}(\mathbf{div},\Omega).$

On the other hand, we can notice that the kernel of ${B}_{1}$ is given by

$\mathbb{V}:=\left\\{\vec{\boldsymbol{\tau}}\in\mathbf{X}:[B_{1}(\vec{\boldsymbol{\tau}}),\vec{\boldsymbol{v}}]=0\
\forall\,\vec{\boldsymbol{v}}\in\mathbf{M}\right\\}=\left\\{\vec{\boldsymbol{\tau}}\in\mathbf{X}:\mathbf{div}\,{\boldsymbol{\tau}}=0\
\mathrm{and}\ {\boldsymbol{\tau}}={\boldsymbol{\tau}}^{\mathrm{t}}\
\mathrm{in}\ \Omega\right\\}.$

The following two results will serve to establish the inf-sup conditions for
$B$ and $B_{1}$ (cf. second and third equations in (3.4a)).

######  Lemma 3.3.

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

$\sup_{\vec{\boldsymbol{\tau}}\in\mathbf{X}\setminus\\{\boldsymbol{0}\\}}\frac{[B_{1}(\vec{\boldsymbol{\tau}}),\vec{\boldsymbol{v}}]}{\left\lVert\vec{\boldsymbol{\tau}}\right\rVert_{\mathbf{X}}}\geq\hat{\beta}_{1}\left\lVert\vec{\boldsymbol{v}}\right\rVert_{\mathbf{M}}\quad\forall\,\vec{\boldsymbol{v}}\in\mathbf{M}.$

######  Proof 3.4.

The (equivalent) surjectivity of $B_{1}$ follows as a slight modification to
the proof of [17, Section 2.4.3].

######  Lemma 3.5.

There exists ${\beta}>0,$ such that

(3.9)
$\quad\sup_{{\boldsymbol{r}}\in\mathbb{L}^{2}(\Omega)\setminus\\{\boldsymbol{0}\\}}\frac{[B({\boldsymbol{r}}),\vec{\boldsymbol{\tau}}]}{\left\lVert{\boldsymbol{r}}\right\rVert_{\mathbb{L}^{2}(\Omega)}}\geq\hat{\beta}\left\lVert\vec{\boldsymbol{\tau}}\right\rVert_{\mathbf{X}}\quad\forall\,\vec{\boldsymbol{\tau}}\in\mathbb{V}.$

######  Proof 3.6.

It follows by using (3.7) and (3.8). For more details see [18, Lemma 2.5].

Before analysing (3.3a)-(3.3d), we study the reduced problem: For each
$\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$, find
$(\hat{\boldsymbol{t}},\hat{\vec{\boldsymbol{\sigma}}},\hat{p})\in\mathbb{L}^{2}(\Omega)\times\mathbb{V}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$
such that

(3.10a) $\displaystyle[A(\hat{\boldsymbol{t}}),\boldsymbol{r}]$
$\displaystyle+\quad$
$\displaystyle[B^{*}(\hat{\vec{\boldsymbol{\sigma}}}),\boldsymbol{r}]$
$\displaystyle=$ $\displaystyle\;\;[\widetilde{H}_{\omega},\boldsymbol{r}],$
(3.10b) $\displaystyle[B(\hat{\boldsymbol{t}}),\vec{\boldsymbol{\tau}}]\;\;\,$
$\displaystyle-\quad$
$\displaystyle[C(\hat{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{\tau}}]\;\;\,$
$\displaystyle+$
$\displaystyle\;\;[B_{2}^{*}(\hat{p}),\vec{\boldsymbol{\tau}}]\;$
$\displaystyle=$
$\displaystyle\;\;[\widetilde{F}_{1},\vec{\boldsymbol{\tau}}],$ (3.10c)
$\displaystyle[B_{2}(\hat{\vec{\boldsymbol{\sigma}}}),q]$ $\displaystyle-$
$\displaystyle\;\;[D(\hat{p}),q]$ $\displaystyle=$
$\displaystyle\;\;[\widetilde{G},q],$

for all
$(\boldsymbol{r},\vec{\boldsymbol{\tau}},q)\in\mathbb{L}^{2}(\Omega)\times\mathbb{V}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$,
where
$\hat{\vec{\boldsymbol{\sigma}}}:=(\hat{\boldsymbol{\sigma}},\hat{\widetilde{p}})$,
and $\widetilde{H}_{\omega},\widetilde{F}_{1}$ and $\widetilde{G}$ are the
generic functionals mentioned at the beginning of Section 3.2. The unique
solvability of (3.10a)-(3.10c) proceeds using Fredholm’s alternative. Let us
recast (3.10a)-(3.10c) with $\omega$ given, as the following equivalent
operator problem: find
$\vec{\boldsymbol{t}}:=(\hat{\boldsymbol{t}},\hat{\vec{\boldsymbol{\sigma}}},\hat{p})\in\mathbf{H}:=\mathbb{L}^{2}(\Omega)\times\mathbb{V}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$
such that

(3.11) $(\mathcal{S}+\mathcal{T})\vec{\boldsymbol{t}}=\mathcal{F},$

where the linear operators
$\mathcal{S},\mathcal{T}:\mathbf{H}\to\mathbf{H}^{*}$, and
$\mathcal{F}\in\mathbf{H}^{*}$ are defined, for all
$\vec{\boldsymbol{r}}:=(\boldsymbol{r},\vec{\boldsymbol{\tau}},q)\in\mathbf{H}$,
as

$\displaystyle[\mathcal{S}(\vec{\boldsymbol{t}}),\vec{\boldsymbol{r}}]:=[A(\hat{\boldsymbol{t}}),\boldsymbol{r}]+[B^{*}(\hat{\vec{\boldsymbol{\sigma}}}),\boldsymbol{r}]-[B(\hat{\boldsymbol{t}}),\vec{\boldsymbol{\tau}}]+[C(\hat{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{\tau}}]+[D(\hat{p}),q],$
$\displaystyle[\mathcal{T}(\vec{\boldsymbol{t}}),\vec{\boldsymbol{r}}]:=-[B_{2}^{*}(\hat{p}),\vec{\boldsymbol{\tau}}]-[B_{2}(\hat{\vec{\boldsymbol{\sigma}}}),q],\qquad[\mathcal{F},\vec{\boldsymbol{r}}]:=[\widetilde{H}_{\omega},r]-[\widetilde{F}_{1},\vec{\boldsymbol{\tau}}]-[\widetilde{G},q].$

The three upcoming lemmas establish that $\mathcal{S}$ is invertible,
$\mathcal{T}$ is compact, and $\mathcal{S}+\mathcal{T}$ is injective.
Fredholm’s theory will yield the well-posedness of (3.11) and, equivalently,
that of (3.10a)-(3.10c), with $\omega\in\mathrm{H}_{\Gamma}^{1}(\Omega)$
given.

######  Lemma 3.7.

The operator $\mathcal{S}:\mathbf{H}\to\mathbf{H}^{*}$ is invertible.

Asserting the invertibility of $\mathcal{S}$ is equivalent to proving the
unique solvability of the following uncoupled problems: find
$(\hat{\boldsymbol{t}},\hat{\vec{\boldsymbol{\sigma}}})\in\mathbb{L}^{2}(\Omega)\times\mathbb{V}$
such that

(3.12a) $\displaystyle[A(\hat{\boldsymbol{t}}),\boldsymbol{r}]$
$\displaystyle+\quad$
$\displaystyle[B^{*}(\hat{\vec{\boldsymbol{\sigma}}}),\boldsymbol{r}]$
$\displaystyle=$ $\displaystyle\;\;[F_{\omega},\boldsymbol{r}]$
$\displaystyle\quad\forall\,\boldsymbol{r}\in\mathbb{L}^{2}(\Omega),$ (3.12b)
$\displaystyle[B(\hat{\boldsymbol{t}}),\vec{\boldsymbol{\tau}}]\;\;\,$
$\displaystyle-\quad$
$\displaystyle[C(\hat{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{\tau}}]\;\;\,$
$\displaystyle=$ $\displaystyle\;\;[F_{\mathbf{X}},\vec{\boldsymbol{\tau}}]$
$\displaystyle\forall\,\vec{\boldsymbol{\tau}}\in\mathbb{V},$

and: find $\hat{p}\in\mathrm{H}^{1}_{\Sigma}(\Omega)$ such that

(3.13)
$[D(\hat{p}),q]=[F_{H},q]\quad\forall\,q\in\mathrm{H}^{1}_{\Sigma}(\Omega),$

where, defining
$\mathcal{F}:=(\mathcal{F}_{\omega},\mathcal{F}_{\mathbf{X}},\mathcal{F}_{H})$,
the functionals $F_{\omega},F_{\mathbf{X}}$, and $F_{H}$ are the ones induced
by the operators
$\mathcal{F}_{\omega},\mathcal{F}_{\mathbf{X}},\mathcal{F}_{H}$, respectively.
The unique solution of (3.13) simply follows by Lax-Milgram’s lemma. In turn,
for (3.12a)-(3.12b), we verify the conditions given by the linear version of
[21, Lemma 2.1]. We first note that

(3.14)
$[A(\boldsymbol{r}),\boldsymbol{r}]=2\mu_{s}\|\boldsymbol{r}\|^{2}_{\mathbb{L}^{2}(\Omega)}\quad\forall\,\boldsymbol{r}\in\mathbb{L}^{2}(\Omega),$

which means that the operator $A$ is elliptic on $\mathbb{L}^{2}(\Omega)$.
Moreover, we observe that

(3.15)
$[C(\vec{\boldsymbol{\tau}}),\vec{\boldsymbol{\tau}}]=\frac{1}{\lambda_{s}}\|\widetilde{q}\|^{2}_{\mathrm{L}^{2}(\Omega)}\geq
0\quad\forall\,\vec{\boldsymbol{\tau}}\in\mathbb{V},$

which shows that the operator $C$ is semi positive definite on $\mathbb{V}$.
In this way, having in mind (3.14), (3.15), and the inf-sup condition for $B$
given by Lemma 3.5, we simply apply [21, Lemma 2.1], and obtain the desired
result.

######  Lemma 3.8.

The operator $\mathcal{T}:\mathbf{H}\to\mathbf{H}^{*}$ is compact.

######  Proof 3.9.

Let us to define the operator
$\mathbb{B}:\mathrm{L}^{2}(\Omega)\to\mathrm{H}^{1}_{\Sigma}(\Omega)$ as

$\langle\mathbb{B}(\widetilde{q}),q\rangle_{\mathrm{L}^{2}(\Omega)}:=\frac{\alpha}{\lambda_{s}}\int_{\Omega}\widetilde{q}q\quad\forall\,\widetilde{q}\in\mathrm{L}^{2}(\Omega),\;\forall\,q\in\mathrm{H}^{1}_{\Sigma}(\Omega),$

where $\langle\cdot,\cdot\rangle_{\mathrm{L}^{2}(\Omega)}$ stands for the
${\mathrm{L}^{2}(\Omega)}$-inner product. Then, thanks to the compactness of
the adjoint operator $\mathbb{B}^{*}$ (see, [35, Lemma 2.2]), we deduce that
the following operator is also compact

$\mathcal{T}(\vec{\boldsymbol{t}})=(\mathbf{0},(\boldsymbol{0},\mathbb{B}(\hat{\widetilde{p}})),\mathbb{B}^{*}(\hat{p})).$

######  Lemma 3.10.

The operator $(\mathcal{S}+\mathcal{T}):\mathbf{H}\to\mathbf{H}^{*}$ is
injective.

######  Proof 3.11.

It is sufficient to show that the unique solution of the homogeneous problem
of (3.10a)-(3.10c) is the null-operator in the space $\mathbf{H}.$ Thus, we
begin by considering
$\widetilde{H}_{\omega}=\widetilde{F}_{1}=\widetilde{G}=0$ in (3.10a)-(3.10c),
and then, taking $\boldsymbol{r}=\hat{\boldsymbol{t}}$, and
$\vec{\boldsymbol{\tau}}=\hat{\vec{\boldsymbol{\sigma}}}$, in (3.10a), and
(3.10b) respectively, and applying some computations, we obtain

(3.16)
$2\mu_{s}\|\hat{\boldsymbol{t}}\|^{2}_{\mathbb{L}^{2}(\Omega)}+\frac{1}{\lambda_{s}}\|\hat{\widetilde{p}}\|^{2}_{\mathrm{L}^{2}(\Omega)}-\frac{\alpha}{\lambda_{s}}\|\hat{\widetilde{p}}\|_{\mathrm{L}^{2}(\Omega)}\|\hat{p}\|_{\mathrm{H}^{1}(\Omega)}\leq
0.$

In turn, choosing $q=\hat{p}$ in (3.10c), we have

(3.17)
$\mathrm{min}\left\\{c_{0},\frac{\kappa_{1}}{\mu_{f}}\right\\}\|\hat{p}\|^{2}_{\mathrm{H}^{1}(\Omega)}+\frac{\alpha^{2}}{\lambda_{s}}\|\hat{p}\|^{2}_{\mathrm{L}^{2}(\Omega)}-\frac{\alpha}{\lambda_{s}}\|\hat{\widetilde{p}}\|_{\mathrm{L}^{2}(\Omega)}\|\hat{p}\|_{\mathrm{H}^{1}(\Omega)}\leq
0.$

Thus, by adding (3.16) and (3.17), and applying Young’s inequality, we readily
obtain
$\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}=\|\hat{p}\|_{\mathrm{H}^{1}(\Omega)}=0$.
On the other hand, from Lemma 3.5, and using (3.10a), we get

$\hat{\beta}\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}\leq\sup_{\mbox{\scriptsize$\boldsymbol{r}\in\mathbb{L}^{2}(\Omega)\atop\boldsymbol{r}\neq
0$}}\frac{[B(\boldsymbol{r}),\hat{\vec{\boldsymbol{\sigma}}}]}{\|{\boldsymbol{r}\|_{\mathbb{L}^{2}(\Omega)}}}=\sup_{\mbox{\scriptsize$\boldsymbol{r}\in\mathbb{L}^{2}(\Omega)\atop\boldsymbol{r}\neq
0$}}\frac{-[A(\hat{\boldsymbol{t}}),\boldsymbol{r}]}{\|{\boldsymbol{r}\|_{\mathbb{L}^{2}(\Omega)}}}\leq
2\mu_{s}\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)},$

from which, we deduce that
$\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}=0$, concluding the proof.

As announced above, thanks to Lemmas 3.7, 3.8 and 3.10, and the Fredholm’s
alternative, we deduce the existence of a unique solution to (3.10a)-(3.10c),
again, with a given $\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$. We complement
the above result with the stability of (3.10a)-(3.10c), for a given
$\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$.

######  Lemma 3.12.

For each $\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$, there exists a constant
$\hat{C}>0$ independent of $\lambda_{s}$, such that

$\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}+\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}+\|\hat{p}\|_{\mathrm{H}^{1}(\Omega)}\leq\hat{C}\left(\|\widetilde{H}_{\omega}\|_{\mathbb{L}^{2}(\Omega)^{\prime}}+\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}+\|\widetilde{G}\|_{\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}}\right).$

######  Proof 3.13.

We proceed similarly as in [29, Theorem 2.1] (see also [5, Theorem 4.3.1]), to
obtain estimates and then use linearity. Firstly, we assume that
$\widetilde{F}_{1}=0$ and bound the solution in terms of
$\widetilde{H}_{\omega}$ and $\widetilde{G}$, and secondly, we assume that
$\widetilde{H}_{\omega}=0,\,\widetilde{G}=0$ and deduce an estimate for the
solution in terms of $\widetilde{F}_{1}$.

Step 1. ($\widetilde{F}_{1}=0$). Proceeding exactly as in the proof of Lemma
3.10, it can be deduced the bounds

(3.18a)
$\displaystyle\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}+\|\hat{p}\|_{\mathrm{H}^{1}(\Omega)}$
$\displaystyle\leq
C_{1}\left(\|\widetilde{H}_{\omega}\|_{\mathbb{L}^{2}(\Omega)^{\prime}}+\|\widetilde{G}\|_{\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}}\right),$
(3.18b)
$\displaystyle\hat{\beta}\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}$
$\displaystyle\leq\sup_{\mbox{\scriptsize$\boldsymbol{r}\in\mathbb{L}^{2}(\Omega)\atop\boldsymbol{r}\neq
0$}}\frac{[B(\boldsymbol{r}),\hat{\vec{\boldsymbol{\sigma}}}]}{\|{\boldsymbol{r}\|_{\mathbb{L}^{2}(\Omega)}}}=\sup_{\mbox{\scriptsize$\boldsymbol{r}\in\mathbb{L}^{2}(\Omega)\atop\boldsymbol{r}\neq
0$}}\frac{[\widetilde{H}_{\omega},\boldsymbol{r}]-[A(\hat{\boldsymbol{t}}),\boldsymbol{r}]}{\|{\boldsymbol{r}\|_{\mathbb{L}^{2}(\Omega)}}}\leq\|\widetilde{H}_{\omega}\|_{\mathbb{L}^{2}(\Omega)^{\prime}}+2\mu_{s}\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)},$

where $C_{1}$ is a constant independent of $\lambda_{s}$. Therefore, the
desired estimate then follows from (3.18a) and (3.18b).

Step 2. ($\widetilde{H}_{\omega}=0$ and $\widetilde{G}=0$). We start by
defining

$[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),({\vec{\boldsymbol{\tau}}},q)]:=[C(\hat{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{\tau}}]-[B_{2}^{*}(\hat{p}),\vec{\boldsymbol{\tau}}]-[B_{2}(\hat{\vec{\boldsymbol{\sigma}}}),q]+D[(\hat{p}),q]=\frac{1}{\lambda_{s}}\int_{\Omega}(\hat{\widetilde{p}}-\alpha\hat{p})(\hat{\widetilde{q}}-\alpha
q)+c_{0}\int_{\Omega}\hat{p}q,$

and noticing, for all
$\hat{\vec{\boldsymbol{\sigma}}},\vec{\boldsymbol{\tau}}\in\mathbb{V}$ and
$\hat{p},q\in\mathrm{H}^{1}_{\Sigma}(\Omega)$, that

$\displaystyle|[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),(\vec{\boldsymbol{\tau}},q)]|$
$\displaystyle\leq[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),(\hat{\vec{\boldsymbol{\sigma}}},\hat{p})]^{1/2}[S(\vec{\boldsymbol{\tau}},q),(\vec{\boldsymbol{\tau}},q)]^{1/2}$
(3.19)
$\displaystyle=\left(\frac{1}{\lambda_{s}}\|\hat{\widetilde{p}}-\alpha\hat{p}\|^{2}_{\mathrm{L}^{2}(\Omega)}+c_{0}\|\hat{p}\|^{2}_{\mathrm{L}^{2}(\Omega)}\right)^{1/2}\left(\frac{1}{\lambda_{s}}\|\hat{\widetilde{q}}-\alpha
q\|^{2}_{\mathrm{L}^{2}(\Omega)}+c_{0}\|q\|^{2}_{\mathrm{L}^{2}(\Omega)}\right)^{1/2}.$

Additionally, from the inf-sup condition (3.9) it can be deduced that

(3.20)
$\sup_{\vec{\boldsymbol{\tau}}\in\mathbb{V}\setminus\\{\boldsymbol{0}\\}}\frac{[B({\boldsymbol{r}^{\perp}}),\vec{\boldsymbol{\tau}}]}{\|\vec{\boldsymbol{\tau}}\|_{\mathbf{X}}}\geq\hat{\beta}\left\lVert\boldsymbol{r}^{\perp}\right\rVert_{\mathbb{L}^{2}(\Omega)}\quad\forall\,\boldsymbol{r}^{\perp}\in\widetilde{\mathbb{V}}^{\perp},$

where
$\widetilde{\mathbb{V}}:=\left\\{\hat{\boldsymbol{r}}\in\mathbb{L}^{2}(\Omega):[B(\hat{\boldsymbol{r}}),\vec{\boldsymbol{\tau}}]=0\
\forall\,\vec{\boldsymbol{\tau}}\in\mathbb{V}\right\\}.$ Now, we let
$\hat{\boldsymbol{t}}_{0}\in\widetilde{\mathbb{V}}$ and
$\hat{\boldsymbol{t}}^{\perp}\in\widetilde{\mathbb{V}}^{\perp}$ be such that
$\hat{\boldsymbol{t}}=\hat{\boldsymbol{t}}_{0}+\hat{\boldsymbol{t}}^{\perp}$,
and observe from (3.10a)-(3.10c) that there hold

(3.21a)
$\displaystyle[A(\hat{\boldsymbol{t}}),\hat{\boldsymbol{t}}]+[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),(\hat{\vec{\boldsymbol{\sigma}}},\hat{p})]$
$\displaystyle=-\widetilde{F}_{1}(\hat{\vec{\boldsymbol{\sigma}}}),$ (3.21b)
$\displaystyle[B(\hat{\boldsymbol{t}}^{\perp}),\vec{\boldsymbol{\tau}}]-[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),(\vec{\boldsymbol{\tau}},q)]$
$\displaystyle=\widetilde{F}_{1}(\vec{\boldsymbol{\tau}})\quad\forall\,\vec{\boldsymbol{\tau}}\in\mathbb{V},\quad\forall\,q\in\mathrm{H}^{1}_{\Sigma}(\Omega).$

Thus, from (3.20) with $\boldsymbol{r}^{\perp}=\hat{\boldsymbol{t}}^{\perp}$,
estimate (3.21b), the continuity of $\widetilde{F}_{1}$ and the bound (3.13),
we get

$\displaystyle\hat{\beta}\|\hat{\boldsymbol{t}}^{\perp}\|_{\mathbb{L}^{2}(\Omega)}$
$\displaystyle\leq\sup_{\vec{\boldsymbol{\tau}}\in\mathbb{V}\setminus\\{\boldsymbol{0}\\}}\frac{[B({\hat{\boldsymbol{t}}^{\perp}}),\vec{\boldsymbol{\tau}}]}{\|\vec{\boldsymbol{\tau}}\|_{\mathbf{X}}}=\sup_{\vec{\boldsymbol{\tau}}\in\mathbb{V}\setminus\\{\boldsymbol{0}\\}}\frac{\widetilde{F}_{1}(\vec{\boldsymbol{\tau}})+[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),(\vec{\boldsymbol{\tau}},q)]}{\|\vec{\boldsymbol{\tau}}\|_{\mathbf{X}}}$
(3.22)
$\displaystyle\leq\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}+[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),(\hat{\vec{\boldsymbol{\sigma}}},\hat{p})]^{1/2}\sup_{\vec{\boldsymbol{\tau}}\in\mathbb{V}\setminus\\{\boldsymbol{0}\\}}\frac{\left(\frac{1}{\lambda_{s}}\|\hat{\widetilde{q}}-\alpha
q\|^{2}_{\mathrm{L}^{2}(\Omega)}+c_{0}\|q\|^{2}_{\mathrm{L}^{2}(\Omega)}\right)^{1/2}}{\|\vec{\boldsymbol{\tau}}\|_{\mathbf{X}}}.$

Then, defining

$C_{S}:=\sup_{\vec{\boldsymbol{\tau}}\in\mathbb{V}\setminus\\{\boldsymbol{0}\\}}\frac{\left(\frac{1}{\lambda_{s}}\|\hat{\widetilde{q}}-\alpha
q\|^{2}_{\mathrm{L}^{2}(\Omega)}+c_{0}\|q\|^{2}_{\mathrm{L}^{2}(\Omega)}\right)^{1/2}}{\|\vec{\boldsymbol{\tau}}\|_{\mathbf{X}}},$

which can be seen as a constant independent of $\lambda_{s}$ if
$\lambda_{s}\rightarrow\infty$, and noticing from (3.21a) that

$[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),(\hat{\vec{\boldsymbol{\sigma}}},\hat{p})]\leq-\widetilde{F}_{1}(\hat{\vec{\boldsymbol{\sigma}}}),$

it can be deduced from (3.13) that

(3.23)
$\hat{\beta}\|\hat{\boldsymbol{t}}^{\perp}\|_{\mathbb{L}^{2}(\Omega)}\leq\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}+C_{S}[S(\hat{\vec{\boldsymbol{\sigma}}},\hat{p}),(\hat{\vec{\boldsymbol{\sigma}}},\hat{p})]^{1/2}\leq\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}+C_{S}\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}^{1/2}\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}^{1/2}.$

Furthermore, taking $\hat{\boldsymbol{t}}_{0}$ in (3.10a), applying the
ellipticity and continuity of $A$, and recalling (3.6), we easily get

(3.24)
$\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}\leq(1+\widetilde{C}_{1})\|\hat{\boldsymbol{t}}^{\perp}\|_{\mathbb{L}^{2}(\Omega)}.$

In this way, from (3.24), the inf-sup condition (3.9) and equation (3.10a), it
readily follows that

(3.25)
$\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}\leq\hat{\beta}^{-1}\sup_{\mbox{\scriptsize$\boldsymbol{r}\in\mathbb{L}^{2}(\Omega)\atop\boldsymbol{r}\neq
0$}}\frac{[B(\boldsymbol{r}),\hat{\vec{\boldsymbol{\sigma}}}]}{\|{\boldsymbol{r}\|_{\mathbb{L}^{2}(\Omega)}}}=\hat{\beta}^{-1}\sup_{\mbox{\scriptsize$\boldsymbol{r}\in\mathbb{L}^{2}(\Omega)\atop\boldsymbol{r}\neq
0$}}\frac{-[A(\hat{\boldsymbol{t}}),\boldsymbol{r}]}{\|{\boldsymbol{r}\|_{\mathbb{L}^{2}(\Omega)}}}\leq
2\mu_{s}\hat{\beta}^{-1}\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}\leq\widetilde{C}_{2}\|\hat{\boldsymbol{t}}^{\perp}\|_{\mathbb{L}^{2}(\Omega)},$

which combined with (3.23) and the Young’s inequality yields

$\|\hat{\boldsymbol{t}}^{\perp}\|_{\mathbb{L}^{2}(\Omega)}\leq\widetilde{C}_{3}\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}.$

Therefore, from the latter inequality, (3.24) and (3.25), we derive that

(3.26)
$\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}\leq\widetilde{C}_{4}\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}\quad\mathrm{and}\quad\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}\leq\widetilde{C}_{5}\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}},$

where $\widetilde{C}_{4}>0$ and $\widetilde{C}_{5}>0$ are constants
independent of $\lambda_{s}$. Finally, the ellipticity of $D$ and (3.10c)
yield

(3.27)
$\|\hat{p}\|^{2}_{\mathrm{H}^{1}(\Omega)}\leq\widetilde{C}_{6}\|\hat{p}\|_{\mathrm{H}^{1}(\Omega)}\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}.$

In this way, the result follows from (3.26) and (3.27), and after apply some
computations.

We summarise the foregoing results in the following lemma.

######  Lemma 3.14.

For each $\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$ the problem
(3.10a)-(3.10c) has a unique solution
$(\hat{\boldsymbol{t}},\hat{\vec{\boldsymbol{\sigma}}},\hat{p})\in\mathbb{L}^{2}(\Omega)\times\mathbb{V}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$.
Moreover, there exists $\hat{C}>0$ independent of $\omega$ and $\lambda_{s}$,
such that

$\|\hat{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}+\|\hat{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}+\|\hat{p}\|_{\mathrm{H}^{1}(\Omega)}\leq\hat{C}\left(\|\widetilde{H}_{\omega}\|_{\mathbb{L}^{2}(\Omega)^{\prime}}+\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}+\|\widetilde{G}\|_{\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}}\right).$

Finally, the main result is given next.

######  Lemma 3.15.

For each $\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$ the problem (3.3a)-(3.3d)
has a unique solution
$(\boldsymbol{t},\vec{\boldsymbol{\sigma}},\vec{\boldsymbol{u}},p)\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$.
Moreover, there exists $C>0$ independent of $\omega$ and $\lambda_{s}$, such
that

(3.28)
$\|\boldsymbol{t}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}\|_{\mathbf{X}}+\|\vec{\boldsymbol{u}}\|_{\mathbf{M}}+\|p\|_{\mathrm{H}^{1}(\Omega)}\leq
C\left(\|\widetilde{H}_{\omega}\|_{\mathbb{L}^{2}(\Omega)^{\prime}}+\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}+\|\widetilde{F}\|_{\mathbf{M}^{\prime}}+\|\widetilde{G}\|_{\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}}\right).$

######  Proof 3.16.

We proceed similarly as in [21, Theorem 2.1]. Recalling from [23, Lemma 4.1]
that the result given in Lemma 3.3 implies that
$B_{1}:\mathbb{V}^{\perp}\rightarrow\mathbf{M}^{\prime}$ and
$B_{1}^{*}:\mathbf{M}\rightarrow\mathbb{V}^{\circ}$ are isomorphisms with

(3.29) $\|B_{1}\|,\quad\|(B_{1}^{*})^{-1}\|\leq\frac{1}{\hat{\beta_{1}}},$

where $\mathbb{V}^{\circ}$ stands for the set of functionals in $\mathbf{X}$
that vanish on the elements of $\mathbb{V}$. Now, let
$\vec{\boldsymbol{\sigma}}_{0}:=B_{1}^{-1}(\widetilde{F})\in\mathbb{V}^{\perp}$,
and notice that by using (3.29) it can be deduce that

(3.30)
$\|\vec{\boldsymbol{\sigma}}_{0}\|_{\mathbf{X}}\leq\frac{1}{\beta_{1}}\|\widetilde{F}\|_{\mathbf{M}^{\prime}}.$

With this in mind, we define the functionals
$\underline{H}_{\omega}:=\widetilde{H}_{\omega}-B^{*}(\vec{\boldsymbol{\sigma}}_{0})$,
$\underline{F}_{1}:=\widetilde{F}_{1}+C(\vec{\boldsymbol{\sigma}}_{0})$ and
$\underline{G}:=\widetilde{G}-B_{2}(\vec{\boldsymbol{\sigma}}_{0})$, and
consider the following problem: Find
$(\underline{\boldsymbol{t}},\underline{\vec{\boldsymbol{\sigma}}},\underline{p})\in\mathbb{L}^{2}(\Omega)\times\mathbb{V}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$
such that

$\displaystyle[A(\underline{\boldsymbol{t}}),\boldsymbol{r}]$
$\displaystyle+\quad$
$\displaystyle[B^{*}(\underline{\vec{\boldsymbol{\sigma}}}),\boldsymbol{r}]$
$\displaystyle=$ $\displaystyle\;\;[\underline{H}_{\omega},\boldsymbol{r}],$
(3.31)
$\displaystyle[B(\underline{\boldsymbol{t}}),\vec{\boldsymbol{\tau}}]\;\;\,$
$\displaystyle-\quad$
$\displaystyle[C(\underline{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{\tau}}]\;\;\,$
$\displaystyle+$
$\displaystyle\;\;[B_{2}^{*}(\underline{p}),\vec{\boldsymbol{\tau}}]\;$
$\displaystyle=$ $\displaystyle\;\;[\underline{F},\vec{\boldsymbol{\tau}}],$
$\displaystyle[B_{2}(\underline{\vec{\boldsymbol{\sigma}}}),q]$
$\displaystyle-$ $\displaystyle\;\;[D(\underline{p}),q]$ $\displaystyle=$
$\displaystyle\;\;[\underline{G},q],$

for all
$(\boldsymbol{r},\vec{\boldsymbol{\tau}},q)\in\mathbb{L}^{2}(\Omega)\times\mathbb{V}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$,
where the involved operators are exactly the ones defining (3.10a)-(3.10c).
Therefore, by noticing that
$\underline{H}_{\omega}\in\mathbb{L}^{2}(\Omega)^{\prime}$,
$\underline{F}\in\mathbb{V}^{\circ}$ and
$\underline{G}\in\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}$, we can simply take
$\widetilde{H}_{\omega}$, $\widetilde{F}_{1}$ and $\widetilde{G}$ in
(3.10a)-(3.10c) as $\underline{H}_{\omega}$, $\underline{F}$ and
$\underline{G}$, respectively, and apply Lemma 3.12 to assert the existence
and uniqueness of a solution
$(\underline{\boldsymbol{t}},\underline{\vec{\boldsymbol{\sigma}}},\underline{p})\in\mathbb{L}^{2}(\Omega)\times\mathbb{V}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$
to problem (3.31).

On the other hand, since
$(\underline{F}-B(\underline{\boldsymbol{t}})+C(\underline{\vec{\boldsymbol{\sigma}}})-B^{*}_{2}(\underline{p}))\in\mathbb{V}^{\circ}$,
we can take
$\underline{\vec{\boldsymbol{u}}}:=(B_{1}^{*})^{-1}(\underline{F}-B(\underline{\boldsymbol{t}})+C(\underline{\vec{\boldsymbol{\sigma}}})-B^{*}_{2}(\underline{p}))\in\mathbf{M}$,
from which, it is clear that

(3.32)
$[B_{1}^{*}(\underline{\vec{\boldsymbol{u}}}),{\vec{\boldsymbol{\tau}}}]=[\underline{F}-B(\underline{\boldsymbol{t}})+C(\underline{\vec{\boldsymbol{\sigma}}})-B^{*}_{2}(\underline{p}),{\vec{\boldsymbol{\tau}}}]\quad\forall\,{\vec{\boldsymbol{\tau}}}\in\mathbf{X},$

and therefore,

(3.33)
$\|{\underline{\vec{\boldsymbol{u}}}}\|_{\mathbf{M}}\leq\frac{1}{\beta_{1}}\|\underline{F}-B(\underline{\boldsymbol{t}})+C(\underline{\vec{\boldsymbol{\sigma}}})-B^{*}_{2}(\underline{p})\|.$

In this manner, noting also that
$B_{1}(\underline{\vec{\boldsymbol{\sigma}}}+\vec{\boldsymbol{\sigma}}_{0})=B_{1}(\vec{\boldsymbol{\sigma}}_{0})=\widetilde{F}$,
we finally conclude for each $\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$, that
$(\underline{\boldsymbol{t}},\underline{\vec{\boldsymbol{\sigma}}}+\vec{\boldsymbol{\sigma}}_{0},\underline{\vec{\boldsymbol{u}}},\underline{p})\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$
is solution of (3.3a)-(3.3d).

Now, for the uniqueness, consider another solution
$({\boldsymbol{t}},{\vec{\boldsymbol{\sigma}}},{\vec{\boldsymbol{u}}},{p})\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$
of (3.3a)-(3.3d). It is fairly directly seen, from (3.3a)-(3.3d), that
$({\boldsymbol{t}},{\vec{\boldsymbol{\sigma}}}-\vec{\boldsymbol{\sigma}}_{0},{p})\in\mathbb{L}^{2}(\Omega)\times\mathbb{V}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$
is also a solution of (3.31), and therefore
$({\boldsymbol{t}},{\vec{\boldsymbol{\sigma}}}-\vec{\boldsymbol{\sigma}}_{0},{p})=(\underline{\boldsymbol{t}},\underline{\vec{\boldsymbol{\sigma}}},\underline{p})$,
which together with (3.32) gives
$({\boldsymbol{t}},{\vec{\boldsymbol{\sigma}}},\vec{\boldsymbol{u}},{p})=(\underline{\boldsymbol{t}},\underline{\vec{\boldsymbol{\sigma}}}+\vec{\boldsymbol{\sigma}}_{0},\underline{\vec{\boldsymbol{u}}},\underline{p})$.

On the other hand, in order to deduce the stability estimate (3.28) we begin
by using (3.33), to obtain

$\|{{\vec{\boldsymbol{u}}}}\|_{\mathbf{M}}\leq\frac{1}{\beta_{1}}\|-B({\boldsymbol{t}})+C({\vec{\boldsymbol{\sigma}}_{0}})-B^{*}_{2}({p})\|,$

from which, by using the first and second inequalities in (3.4b), the inf-sup
condition (3.9), and the bound (3.30), we deduce that there exists
$\underline{C}_{1}>0$ independent of $\lambda_{s}$ such that

(3.34)
$\|{{\vec{\boldsymbol{u}}}}\|_{\mathbf{M}}\leq\underline{C}_{1}\left\\{\|\boldsymbol{t}\|_{\mathbb{L}^{2}(\Omega)}+\frac{1}{\lambda_{s}}\|\widetilde{F}\|_{{}_{\mathbf{M}}^{\prime}}+\frac{\alpha}{\lambda_{s}}\|p\|_{\mathrm{H}^{1}(\Omega)}\right\\}.$

Moreover, applying Lemma 3.14 to (3.31), we deduce the existence of
$\underline{C}>0$, independent of $\lambda_{s}$, such that

$\|\boldsymbol{t}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}-\vec{\boldsymbol{\sigma}}_{0}\|_{\mathbf{X}}+\|p\|_{\mathrm{H}^{1}(\Omega)}\leq\underline{C}\left\\{\|\widetilde{H}_{\omega}-B^{*}(\vec{\boldsymbol{\sigma}}_{0})\|+\|\widetilde{F}_{1}+C(\vec{\boldsymbol{\sigma}}_{0})\|+\|\widetilde{G}-B_{2}(\vec{\boldsymbol{\sigma}}_{0})\|\right\\}.$

Thus, from the definitions in (3.4c) and (3.4d), we proceed as before to get

(3.35)
$\|\boldsymbol{t}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}\|_{\mathbf{X}}+\|p\|_{\mathrm{H}^{1}(\Omega)}\\\
\leq\underline{C}_{2}\left\\{\|\widetilde{H}_{\omega}\|_{\mathbb{L}^{2}(\Omega)^{\prime}}+\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}+\left(1+\frac{1}{\lambda_{s}}+\frac{\alpha}{\lambda_{s}}\right)\|\widetilde{F}\|_{\mathbf{M}^{\prime}}+\|\widetilde{G}\|_{\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}}\right\\},$

where $\underline{C}_{2}>0$ is independent of $\lambda_{s}$. In this way,
replacing (3.35) back into (3.34), we deduce that

$\displaystyle\|\boldsymbol{t}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}\|_{\mathbf{X}}+\|\vec{\boldsymbol{u}}\|_{\mathbf{M}}+\|p\|_{\mathrm{H}^{1}(\Omega)}$
$\displaystyle\leq\underline{C}_{3}\left(1+\frac{\alpha}{\lambda_{s}}\right)\left\\{\|\widetilde{H}_{\omega}\|_{\mathbb{L}^{2}(\Omega)^{\prime}}+\|\widetilde{F}_{1}\|_{\mathbf{X}^{\prime}}+\left(1+\frac{2}{\lambda_{s}}+\frac{\alpha}{\lambda_{s}}\right)\|\widetilde{F}\|_{\mathbf{M}^{\prime}}+\|\widetilde{G}\|_{\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}}\right\\}.$

Therefore, noting that the constants $1+\frac{\alpha}{\lambda_{s}}$ and
$1+\frac{2}{\lambda_{s}}+\frac{\alpha}{\lambda_{s}}$ are understood as
independent of $\lambda_{s}$ in the limit $\lambda_{s}\to\infty$, we can apply
elementary algebraic computations to deduce the existence of a constant $C$
independent of $\lambda_{s}$, such that (3.28) holds.

We establish now the well-posedness of the operator $\widetilde{\mathbf{S}}$.
The following result asserts the unique solvability of (3.3e) for a given
$\boldsymbol{\sigma}\in\mathbb{H}_{\Sigma}(\mathbf{div},\Omega)$.

######  Lemma 3.17.

For each $\boldsymbol{\sigma}\in\mathbb{H}_{\Sigma}(\mathbf{div},\Omega)$, the
problem (3.3e) has a unique solution
$\omega:=\widetilde{\mathbf{S}}(\boldsymbol{\sigma})\in\mathrm{H}^{1}_{\Gamma}(\Omega)$,
and there holds

(3.36)
$\|\widetilde{\mathbf{S}}(\boldsymbol{\sigma})\|:=\|\omega\|_{\mathrm{H}^{1}(\Omega)}\leq\widetilde{\alpha}^{-1}\|J\|_{\mathrm{H}^{1}_{\Gamma}(\Omega)^{\prime}}.$

######  Proof 3.18.

We begin by applying the uniformly positive definiteness of $D$ (cf. (3.2)) to
obtain

(3.37)
$\mathscr{A}_{\boldsymbol{\sigma}}(\omega,\omega)\geq\phi\|\omega\|^{2}_{\mathrm{L}^{2}(\Omega)}+D_{0}|\omega|^{2}_{\mathrm{H}^{1}(\Omega)}\geq\widetilde{\alpha}\|\omega\|^{2}_{\mathrm{H}^{1}(\Omega)},$

which implies the ellipticity of
$\mathscr{A}_{\boldsymbol{\sigma}}(\cdot,\cdot)$, with
$\widetilde{\alpha}:=\mathrm{min}\left\\{D_{0},\phi\right\\}$. Thus, a
straightforward application of Lax-Milgram’s lemma (see, e.g., [17, Thm.
1.1]), proves that for each
$\boldsymbol{\sigma}\in\mathbb{H}_{\Sigma}(\mathbf{div},\Omega)$, problem
(3.3e) has a unique solution
$\omega=\widetilde{\mathbf{S}}(\boldsymbol{\sigma})\in\mathrm{H}^{1}_{\Gamma}(\Omega)$.
Moreover, the corresponding continuous dependence on the data is formulated as

(3.38)
$\|\omega\|_{\mathrm{H}^{1}(\Omega)}\leq\widetilde{\alpha}^{-1}\|J\|_{\mathrm{H}^{1}_{\Gamma}(\Omega)^{\prime}}.$

As $\mathbf{S}$ and $\widetilde{\mathbf{S}}$ are well-defined, we can
guarantee the well-posedness of the operator $\mathbf{T}$.

###  3.3 Solvability of the fixed-point problem

With the aim to use the Schauder fixed-point theorem on $\mathbf{T}$, in what
follows we establish sufficient conditions under which $\mathbf{T}$ maps a
closed ball of $\mathrm{H}^{1}_{\Gamma}(\Omega)$ into itself. Indeed, from now
on we let

(3.39)
$\mathbf{W}:=\left\\{\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega):\;\;\left\lVert\omega\right\rVert_{\mathrm{H}^{1}(\Omega)}\leq
r:=\widetilde{\alpha}^{-1}\|J\|_{\mathrm{H}^{1}_{\Gamma}(\Omega)^{\prime}}\right\\}.$

######  Lemma 3.19.

For the closed ball $\mathbf{W}$, it holds that
$\mathbf{T}(\mathbf{W})\subseteq\mathbf{W}$.

######  Proof 3.20.

It suffices to recall the definition of $\mathbf{T}$ (3.5), and then apply the
estimate (3.36).

We now verify the hypotheses of the Schauder fixed-point theorem. Before
starting the result to be proved, we show the Lipschitz continuity property
for the operators $\mathbf{S}$ and $\widetilde{\mathbf{S}}$, and then for
$\mathbf{T}$.

######  Lemma 3.21.

There exists a constant $C_{\mathbf{S}}>0$, independent of $\lambda_{s}$ such
that

$\displaystyle\|\mathbf{S}(\omega_{1})-\mathbf{S}(\omega_{2})\|$
$\displaystyle\leq
C_{\mathbf{S}}\beta\|\omega_{1}-\omega_{2}\|_{\mathrm{L}^{2}(\Omega)}\quad\forall\,\omega_{1},\omega_{2}\in\mathrm{H}^{1}_{\Gamma}(\Omega),$
$\displaystyle\|\widetilde{\mathbf{S}}(\boldsymbol{\sigma})-\widetilde{\mathbf{S}}(\widetilde{\boldsymbol{\sigma}})\|$
$\displaystyle\leq\frac{d_{3}}{\widetilde{\alpha}}\|\widetilde{\boldsymbol{\sigma}}-\boldsymbol{\sigma}\|_{\mathbb{L}^{2}(\Omega)}\|\widetilde{\mathbf{S}}(\widetilde{\boldsymbol{\sigma}})\|_{W^{1,\infty}(\Omega)}\quad\forall\,{\boldsymbol{\sigma}},\widetilde{\boldsymbol{\sigma}}\in\mathbb{H}_{\Sigma}(\mathbf{div},\Omega),$

where $\beta,\widetilde{\alpha}$ and $d_{3}$ are the constants given in
(2.4a), (3.37) and (3.2), respectively.

######  Proof 3.22.

We begin with the proof for $\mathbf{S}$. Let
$(\boldsymbol{t}_{1},\vec{\boldsymbol{\sigma}}_{1},\vec{\boldsymbol{u}}_{1},p_{1})$
and
$(\boldsymbol{t}_{2},\vec{\boldsymbol{\sigma}}_{2},\vec{\boldsymbol{u}}_{2},p_{2})$
be two solutions of the problem (3.3a)-(3.3d), with $\omega_{1}$ and
$\omega_{2}\in\mathrm{H}^{1}_{\Gamma}(\Omega)$ given, such that
$\mathbf{S}(\omega_{1}):=(\boldsymbol{t}_{1},\vec{\boldsymbol{\sigma}}_{1},\vec{\boldsymbol{u}}_{1},p_{1})$
and
$\mathbf{S}(\omega_{2}):=(\boldsymbol{t}_{2},\vec{\boldsymbol{\sigma}}_{2},\vec{\boldsymbol{u}}_{2},p_{2})$.
Then, by applying the linearity of the involved bilinear forms and
functionals, we obtain

$\displaystyle A(\boldsymbol{t}_{1}-\boldsymbol{t}_{2},\boldsymbol{r})$
$\displaystyle+\quad$ $\displaystyle
B(\boldsymbol{r},\vec{\boldsymbol{\sigma}}_{1}-\vec{\boldsymbol{\sigma}}_{2})$
$\displaystyle=$ $\displaystyle\;H_{\omega_{1}-\omega_{2}}(\boldsymbol{r}),$
$\displaystyle
B(\boldsymbol{t}_{1}-\boldsymbol{t}_{2},\vec{\boldsymbol{\tau}})\;\;\,$
$\displaystyle-\quad$ $\displaystyle
C(\vec{\boldsymbol{\tau}},\vec{\boldsymbol{\sigma}}_{1}-\vec{\boldsymbol{\sigma}}_{2})\;\;\,$
$\displaystyle+$
$\displaystyle\;\;B_{1}(\vec{\boldsymbol{\tau}},\vec{\boldsymbol{u}}-\vec{\boldsymbol{u}_{2}})\;\;$
$\displaystyle+$
$\displaystyle\;\;B_{2}(\vec{\boldsymbol{\tau}},p_{1}-p_{2})\;$
$\displaystyle=$ $\displaystyle\;0,$ $\displaystyle
B_{1}(\vec{\boldsymbol{\sigma}}_{1}-\vec{\boldsymbol{\sigma}}_{2},\vec{\boldsymbol{v}})$
$\displaystyle=$ $\displaystyle\;0,$ $\displaystyle
B_{2}(\vec{\boldsymbol{\sigma}},q)$ $\displaystyle-$ $\displaystyle\;\;D(p,q)$
$\displaystyle=$ $\displaystyle\;0,$

for each
$(\boldsymbol{r},\vec{\boldsymbol{\tau}},\vec{\boldsymbol{v}},p)\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)$.
Proceeding as in the proof of Lemma 3.14, we deduce that

(3.40)
$\|\mathbf{S}(\omega_{1})-\mathbf{S}(\omega_{2})\|=\|(\boldsymbol{t}_{1},\vec{\boldsymbol{\sigma}}_{1},\vec{\boldsymbol{u}}_{1},p_{1})-(\boldsymbol{t}_{2},\vec{\boldsymbol{\sigma}}_{2},\vec{\boldsymbol{u}}_{2},p_{2})\|\leq
C_{\mathbf{S}}\beta\|\omega_{1}-\omega_{2}\|_{\mathrm{L}^{2}(\Omega)}.$

In turn, for the problem defined by $\widetilde{\mathbf{S}}$, given
$\boldsymbol{\sigma}$ and
$\widetilde{\boldsymbol{\sigma}}\in\mathbb{H}_{\Sigma}(\mathbf{div},\Omega)$,
we let $\omega$ and $\widetilde{\omega}$ be two solutions of problem (3.3e),
such that $\omega:={\widetilde{\mathbf{S}}}(\boldsymbol{\sigma})$ and
$\widetilde{\omega}:=\widetilde{\mathbf{S}}(\widetilde{\boldsymbol{\sigma}})$.
Then, adding and subtracting appropriate terms, and utilising the ellipticity
of $\mathscr{A}_{\boldsymbol{\sigma}}$, we find that

$\displaystyle\widetilde{\alpha}\|\omega-\widetilde{\omega}\|^{2}_{\mathrm{H}^{1}(\Omega)}$
$\displaystyle\leq\mathscr{A}_{\boldsymbol{\sigma}}(\omega,\omega-\widetilde{\omega})-\mathscr{A}_{\boldsymbol{\sigma}}(\widetilde{\omega},\omega-\widetilde{\omega})=(\mathscr{A}_{\widetilde{\boldsymbol{\sigma}}}-\mathscr{A}_{{\boldsymbol{\sigma}}})(\widetilde{\omega},\omega-\widetilde{\omega})$
(3.41)
$\displaystyle=\int_{\Omega}(D(\widetilde{\boldsymbol{\sigma}})-D({\boldsymbol{\sigma}}))\nabla\widetilde{\omega}\cdot\nabla(\omega-\widetilde{\omega})\leq
d_{3}\|\widetilde{\boldsymbol{\sigma}}-\boldsymbol{\sigma}\|_{\mathbb{L}^{2}(\Omega)}\|\widetilde{\omega}\|_{W^{1,\infty}(\Omega)}\|\omega-\widetilde{\omega}\|_{\mathrm{H}^{1}(\Omega)}.$

Thus, the proof ends after noticing that from (3.22) we obtain

(3.42)
$\|\widetilde{\mathbf{S}}(\boldsymbol{\sigma})-\widetilde{\mathbf{S}}(\widetilde{\boldsymbol{\sigma}})\|_{\mathrm{H}^{1}(\Omega)}:=\|\omega-\widetilde{\omega}\|^{2}_{\mathrm{H}^{1}(\Omega)}\leq\frac{d_{3}}{\widetilde{\alpha}}\|\widetilde{\boldsymbol{\sigma}}-\boldsymbol{\sigma}\|_{\mathbb{L}^{2}(\Omega)}\|\widetilde{\mathbf{S}}(\widetilde{\boldsymbol{\sigma}})\|_{W^{1,\infty}(\Omega)}.$

Now, we are able to show the announced property of the operator $\mathbf{T}$.

######  Lemma 3.23.

There exists a constant $C_{\mathbf{T}}>0$, independent of $\lambda_{s}$ such
that

(3.43)
$\|\mathbf{T}(\omega_{1})-\mathbf{T}(\omega_{2})\|_{\mathrm{H}^{1}(\Omega)}\leq
C_{\mathbf{T}}d_{3}\beta\|\mathbf{T}(\omega_{2})\|_{W^{1,\infty}(\Omega)}\|\omega_{1}-\omega_{2}\|_{\mathrm{L}^{2}(\Omega)}\quad\forall\,\omega_{1},\omega_{2}\in\mathrm{H}^{1}_{\Gamma}(\Omega).$

######  Proof 3.24.

It suffices to recall from Section 3.1 that
$\mathbf{T}(\omega)=\widetilde{\mathbf{S}}(\mathbf{S}_{2}(\omega))$ for all
$\omega\in\mathrm{H}^{1}_{\Gamma}(\Omega)$, and apply Lemma 3.21.

The next lemma establishes the continuity and compactness of $\mathbf{T}$ on
$\mathbf{W}$.

######  Lemma 3.25.

Let $\mathbf{W}$ be as in (3.39). Then, $\mathbf{T}:\mathbf{W}\to\mathbf{W}$
is continuous and $\overline{\mathbf{T}(\mathbf{W})}$ is compact.

######  Proof 3.26.

The continuity of
$i_{c}:\mathrm{H}^{1}(\Omega)\rightarrow\mathrm{L}^{2}(\Omega)$ in combination
with the estimate (3.43), give the continuity of the operator $\mathbf{T}$.
Finally, thanks to the compactness of $i_{c}$ and the fact that every bounded
sequence in a Hilbert space has a weakly convergent subsequence, we can assert
the compactness of $\overline{\mathbf{T}(\mathbf{W})}$, concluding the proof.

We can now establish the existence of a solution to (3.3a)-(3.3e), which
follows straightforwardly from Schauder’s theorem together with Lemmas 3.23
and 3.25.

######  Theorem 3.27.

Let $\mathbf{W}$ be as in (3.39). Then, there exists $\overline{C}>0$,
independent of $\lambda_{s}$, such that $\mathrm{(\ref{VF-
poro_1})-\eqref{VF_diff}}$ has at least one solution
$(\boldsymbol{t},\vec{\boldsymbol{\sigma}},\vec{\boldsymbol{u}},p,\omega)\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)\times\mathrm{H}^{1}_{\Gamma}(\Omega)$
with $\omega\in\mathbf{W}$, satisfying

(3.44)
$\|\boldsymbol{t}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}\|_{\mathbf{X}}+\|\vec{\boldsymbol{u}}\|_{\mathbf{M}}+\|p\|_{\mathrm{H}^{1}(\Omega)}+\|\omega\|_{\mathrm{H}^{1}(\Omega)}\leq\overline{C}\left(\|F\|_{\mathbf{M}^{\prime}}+\|G\|_{\mathrm{H}^{1}_{\Sigma}(\Omega)^{\prime}}+\|J\|_{\mathrm{H}^{1}_{\Gamma}(\Omega)^{\prime}}\right).$

######  Proof 3.28.

It only remains to prove (3.44). For that, it is suffices to apply the result
(3.28) with
$\widetilde{H}_{\omega}=H_{\omega},\widetilde{F}_{1}=\mathrm{O},\widetilde{F}=F$
and $\widetilde{G}=G$, the estimate (3.36), the corresponding bound for
$\|H_{\omega}\|_{\mathbb{L}^{2}(\Omega)^{\prime}}$ and the fact that
$\omega\in\mathbf{W}$.

Finally, we establish the uniqueness of solution to (3.3a)-(3.3e) based on an
additional smallness assumption.

######  Theorem 3.29.

Let
$(\boldsymbol{t},\vec{\boldsymbol{\sigma}},\vec{\boldsymbol{u}},p,\omega)\in\mathbb{L}^{2}(\Omega)\times\mathbf{X}\times\mathbf{M}\times\mathrm{H}^{1}_{\Sigma}(\Omega)\times\mathrm{H}^{1}_{\Gamma}(\Omega)\cap
W^{1,\infty}(\Omega)$ be a solution to problem (3.3a)-(3.3e), and assume that
there exists $M>0$, such that

(3.45) $\|\omega\|_{W^{1,\infty}(\Omega)}\leq
M<\frac{\widetilde{\alpha}}{d_{3}\beta C_{\mathbf{S}}}.$

Then, the solution of problem (3.3a)-(3.3e) is unique.

######  Proof 3.30.

It follows the same steps as in the proof of Lemma 3.21 to obtain the bounds
(3.40) and (3.42), and therefore, the result is merely an application of the
assumption (3.45).

##  4 Finite element method and solvability of the discrete problem

In this section we introduce and analyse the Galerkin scheme associated with
(3.3a)-(3.3e). We consider generic finite dimensional subspaces

$\displaystyle\mathbb{H}^{\boldsymbol{t}}_{h}\subseteq\mathbb{L}^{2}(\Omega),\quad\mathbb{H}^{\boldsymbol{\sigma}}_{h}\subseteq\mathbb{H}_{\Sigma}(\mathbf{div},\Omega),\quad\mathrm{H}^{\widetilde{p}}_{h}(\Omega)\subseteq\mathrm{L}^{2}(\Omega),\quad\mathbf{H}^{\boldsymbol{u}}_{h}\subseteq\mathbf{L}^{2}(\Omega),$
(4.1)
$\displaystyle\mathbb{H}^{\gamma}_{h}(\Omega)\subseteq\mathbb{L}^{2}_{\mathrm{skew}}(\Omega),\quad\mathrm{H}^{p}_{h}\subseteq\mathrm{H}^{1}_{\Sigma}(\Omega),\quad\mathrm{and}\quad\mathrm{H}^{\omega}_{h}\subseteq\mathrm{H}^{1}_{\Gamma}(\Omega),$

which will be specified later on. Then, we define

$\displaystyle\vec{\boldsymbol{\sigma}}_{h}:=(\boldsymbol{\sigma}_{h},\widetilde{p}_{h}),\quad\vec{\boldsymbol{\tau}}_{h}:=(\boldsymbol{\tau}_{h},\widetilde{q}_{h}),\quad\vec{\boldsymbol{u}}_{h}:=(\boldsymbol{u}_{h},\boldsymbol{\gamma}_{h}),\quad\vec{\boldsymbol{v}}_{h}:=(\boldsymbol{v}_{h},\boldsymbol{\eta}_{h}),\quad\mathbf{X}_{h}:=\mathbb{H}_{h}^{\boldsymbol{\sigma}}\times\mathrm{H}^{\widetilde{p}}_{h},\quad\mathbf{M}_{h}:=\mathbf{H}^{\boldsymbol{u}}_{h}\times\mathbb{H}^{\gamma}_{h},$

and a Galerkin scheme for (3.3a)-(3.3e) reads: Find
$(\boldsymbol{t}_{h},\vec{\boldsymbol{\sigma}},\vec{\boldsymbol{u}}_{h},$
$p_{h},\omega_{h})\in\mathbb{H}^{\boldsymbol{t}}_{h}\times\mathbf{X}_{h}\times\mathbf{M}_{h}\times\mathrm{H}^{p}_{h}\times\mathrm{H}^{\omega}_{h}$
such that

(4.2a) $\displaystyle[A(\boldsymbol{t}_{h}),\boldsymbol{r}_{h}]$
$\displaystyle\,+\,$
$\displaystyle[B^{*}(\vec{\boldsymbol{\sigma}}_{h}),\boldsymbol{r}_{h}]$
$\displaystyle\;=\;$ $\displaystyle[H_{\omega_{h}},\boldsymbol{r}_{h}]$
$\displaystyle\quad\forall\,\boldsymbol{r}_{h}\in\mathbb{H}^{\boldsymbol{t}}_{h},$
(4.2b) $\displaystyle[B(\boldsymbol{t}_{h}),\vec{\boldsymbol{\tau}}_{h}]$
$\displaystyle\,-\,$
$\displaystyle[C(\vec{\boldsymbol{\sigma}}_{h}),\vec{\boldsymbol{\tau}}_{h}]$
$\displaystyle\,+\,$
$\displaystyle[B_{1}^{*}(\vec{\boldsymbol{u}}_{h}),\vec{\boldsymbol{\tau}}_{h}]$
$\displaystyle\,+\,$
$\displaystyle[B_{2}^{*}(p_{h}),\vec{\boldsymbol{\tau}}_{h}]$
$\displaystyle\;=\;$ $\displaystyle 0$
$\displaystyle\forall\,\vec{\boldsymbol{\tau}}_{h}\in\mathbf{X}_{h},$ (4.2c)
$\displaystyle[B_{1}(\vec{\boldsymbol{\sigma}}_{h}),\vec{\boldsymbol{v}}_{h}]$
$\displaystyle\;=\;$ $\displaystyle[F,\vec{\boldsymbol{v}}_{h}]$
$\displaystyle\forall\,\vec{\boldsymbol{v}}_{h}\in\mathbf{M}_{h},$ (4.2d)
$\displaystyle[B_{2}(\vec{\boldsymbol{\sigma}}_{h}),q_{h}]$
$\displaystyle\,-\,$ $\displaystyle[D(p_{h}),q_{h}]$ $\displaystyle\;=\;$
$\displaystyle[G,q_{h}]$ $\displaystyle\forall\,p_{h}\in\mathrm{H}^{p}_{h},$
(4.2e)
$\displaystyle[\mathscr{A}_{\boldsymbol{\sigma}_{h}}(\omega_{h}),\theta_{h}]$
$\displaystyle\;=\;$ $\displaystyle[J,\theta_{h}]$
$\displaystyle\forall\,\theta_{h}\in\mathrm{H}^{\omega}_{h},$

In order to address the well-posedness of (4.2a)-(4.2e), we use again a fixed-
point strategy. Let us define
$\textbf{S}_{h}:\mathrm{H}_{h}^{\omega}\to\mathbb{H}^{\boldsymbol{t}}_{h}\times\mathbf{X}_{h}\times\mathbf{M}_{h}\times\mathrm{H}^{p}_{h}$
as

$\displaystyle\textbf{S}_{h}(\omega_{h})$
$\displaystyle:=(\textbf{S}_{1,h}(\omega_{h}),(\textbf{S}_{2,h}(\omega_{h}),\textbf{S}_{3,h}(\omega_{h})),(\textbf{S}_{4,h}(\omega_{h}),\textbf{S}_{5,h}(\omega_{h})),\textbf{S}_{6,h}(\omega_{h}))$
$\displaystyle\;=(\boldsymbol{t}_{h},(\boldsymbol{\sigma}_{h},\widetilde{p}_{h}),(\boldsymbol{u}_{h},\boldsymbol{\gamma}_{h}),p_{h})\in\mathbb{H}^{\boldsymbol{t}}_{h}\times\mathbf{X}_{h}\times\mathbf{M}_{h}\times\mathrm{H}^{p}_{h},$

where
$(\boldsymbol{t}_{h},(\boldsymbol{\sigma}_{h},\widetilde{p}_{h}),(\boldsymbol{u}_{h},\boldsymbol{\gamma}_{h}),p_{h})$,
is the unique solution of (4.2a)-(4.2d) with $\omega_{h}$ given. In turn, let
${\widetilde{\mathbf{S}}}_{h}:\mathbb{H}_{h}^{\boldsymbol{\sigma}}\to\mathrm{H}^{\omega}_{h}$
be the operator defined by

${\widetilde{\mathbf{S}}}_{h}(\boldsymbol{\sigma}_{h})\,:=\,\omega_{h}\quad\forall\,\boldsymbol{\sigma}_{h}\,\in\,\mathbb{H}_{h}^{\boldsymbol{\sigma}},$

where $\omega_{h}$ is the unique solution of (4.2e) with
$\boldsymbol{\sigma}_{h}$ given. Finally, by introducing the operator
$\mbox{{{T}}}_{h}:\mathrm{H}^{\omega}_{h}\to\mathrm{H}^{\omega}_{h}$ as

$\mbox{{{T}}}_{h}(\omega_{h}):=\widetilde{\mathbf{S}}_{h}(\mathbf{S}_{2,h}(\omega_{h}))\quad\forall\,\omega_{h}\in\mathrm{H}^{\omega}_{h},$

we see that solving (4.2a)-(4.2e) is equivalent to seeking a fixed point of
$\mbox{{{T}}}_{h}$, that is: find $\omega_{h}\in\mathrm{H}^{\omega}_{h}$ such
that

(4.3) $\mbox{{{T}}}_{h}(\omega_{h})=\omega_{h}.$

###  4.1 Well-definedness of the operator $\mathbf{T}_{h}$

Here we establish the solvability of (4.2a)-(4.2e) by studying the equivalent
fixed-point problem (4.3). We begin by introducing some needed notations and
preliminary results, as well as specific finite element subspaces satisfying
(4.1).

Let us denote by $\mathcal{T}_{h}$ a regular partition of $\overline{\Omega}$
into triangles (or tetrahedra in 3D) $K$ of diameter $h_{K}$, where
$h:=\max\left\\{h_{K}:\ K\in\mathcal{T}_{h}\right\\}$ is the meshsize. 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 local Raviart-Thomas space of order $k$ as
$\mathbf{RT}_{k}(K):=\mathbf{P}_{k}(K)\oplus\mathrm{P}_{k}(K)\,\boldsymbol{x}$,
where $\mathbf{P}_{k}(K)=[\mathrm{P}_{k}(K)]^{d},$ and $\boldsymbol{x}$ is a
generic vector in $\mathbb{R}^{d}$. Now, let $b_{K}$ be the element bubble
function defined as the unique polynomial in $\mathrm{P}_{d+1}(K)$ vanishing
on $\partial K$ with $\int_{K}b_{K}=1$. Then, for each $K\in\mathcal{T}_{h}$
we consider the bubble space of order $k$, defined as

$\mathbf{B}_{k}(K):=\begin{cases}\mathbf{curl}^{\mathrm{t}}(b_{K}\mathrm{P}_{k}(K))&\textrm{in}\quad\mathbb{R}^{2},\\\
\nabla\times(b_{K}\mathbf{P}_{k}(K))&\textrm{in}\quad\mathbb{R}^{3}.\end{cases}.$

On the other hand, we observe thanks to Lemma 3.7, that the ellipticity of
$A:\mathbb{H}^{\boldsymbol{t}}_{h}\times\mathbb{H}^{\boldsymbol{t}}_{h}\rightarrow\mathbb{R}$
is satisfied for any finite dimensional subspace
$\mathbb{H}^{\boldsymbol{t}}_{h}$, and with the same constant from (3.14).
Therefore $\mathbb{H}^{\boldsymbol{t}}_{h}$ is chosen such that the discrete
inf-sup condition for $B$ holds. Moreover, thanks to the discrete analogue of
$B_{1}$ (cf. first equation in (3.4a)), an inf-sup condition can be ensured by
using the classical PEERSk elements introduced in [2], that is

$\displaystyle\mathbb{H}^{\boldsymbol{\sigma}}_{h}$
$\displaystyle:=\left\\{{\boldsymbol{\tau}}_{h}\in\mathbb{H}_{\Sigma}(\mathbf{div},\Omega):\quad{\boldsymbol{\tau}}_{h}|_{K}\in[\mathbf{RT}_{k}(K)]^{d}\oplus[\mathbf{B}_{k}(K)]^{d}\quad\forall\,K\in\mathcal{T}_{h}\right\\},$
(4.4) $\displaystyle\mathbf{H}^{\boldsymbol{u}}_{h}$
$\displaystyle:=\left\\{\boldsymbol{v}_{h}\in{\mathbf{L}}^{2}(\Omega):\quad\boldsymbol{v}_{h}|_{K}\in\textbf{P}_{k}(K)\quad\forall\,K\in\mathcal{T}_{h}\right\\},$
$\displaystyle\mathbb{H}^{\boldsymbol{\gamma}}_{h}$
$\displaystyle:=\left\\{\boldsymbol{\eta}_{h}\in\mathbb{L}^{2}_{\mathrm{skew}}(\Omega):\quad\boldsymbol{\eta}_{h}\in\mathbf{C}(\overline{\Omega})\quad\mathrm{and}\quad\boldsymbol{\eta}_{h}|_{K}\in\mathbb{P}_{k+1}(K)\quad\forall\,K\in\mathcal{T}_{h}\right\\},$

or by employing the well-known Arnold-Falk-Winther (AFWk, [3]) family of order
$k\geq 0$, that is

(4.5) $\displaystyle\mathbb{H}^{\boldsymbol{\sigma}}_{h}$
$\displaystyle:=\left\\{{\boldsymbol{\tau}}_{h}\in\mathbb{H}_{\Sigma}(\mathbf{div},\Omega):\
{\boldsymbol{\tau}}_{h}|_{K}\in\mathbf{BDM}_{k+1}(K)\quad\forall\,K\in\mathcal{T}_{h}\right\\},$
$\displaystyle\mathbf{H}^{\boldsymbol{u}}_{h}$
$\displaystyle:=\left\\{\boldsymbol{v}_{h}\in{\mathbf{L}}^{2}(\Omega):\
\boldsymbol{v}_{h}|_{K}\in\textbf{P}_{k}(K)\quad\forall\,K\in\mathcal{T}_{h}\right\\},\
\mathbb{H}^{\boldsymbol{\rho}}_{h}:=\left\\{\boldsymbol{\eta}_{h}\in\mathbb{L}^{2}_{\mathrm{skew}}(\Omega):\
\boldsymbol{\eta}_{h}|_{K}\in\mathbb{P}_{k}(K)\quad\forall\,K\in\mathcal{T}_{h}\right\\}.$

Also, we notice that the kernel of ${B}_{1}$ is given by

(4.6)
$\mathbb{V}_{h}:=\left\\{\vec{\boldsymbol{\tau}}_{h}\in\mathbf{X}_{h}:\;\;[B_{1}(\vec{\boldsymbol{\tau}}_{h}),\vec{\boldsymbol{u}_{h}}]=0\quad\forall\,\vec{\boldsymbol{u}}_{h}\in\mathbf{M}_{h}\right\\}.$

Then, by using the definition given for
$\mathbb{H}^{\boldsymbol{\sigma}}_{h}$, $\mathbf{H}^{\boldsymbol{u}}_{h},$ and
$\mathbb{H}^{\boldsymbol{\gamma}}_{h}$, (4.6) becomes
$\mathbb{V}_{h}:={\mathbb{V}}^{\boldsymbol{\sigma}}_{h}\times\mathrm{H}^{p}_{h}$,
where

(4.7)
${\mathbb{V}}^{\boldsymbol{\sigma}}_{h}:=\left\\{{\boldsymbol{\tau}}_{h}\in\mathbb{H}^{\boldsymbol{\sigma}}_{h}:\;\;\mathbf{div}\,\boldsymbol{\tau}_{h}=0\quad\mathrm{in}\quad\Omega\quad\mathrm{and}\quad\int_{\Omega}\boldsymbol{\eta}_{h}\colon\boldsymbol{\tau}_{h}=0\quad\boldsymbol{\forall}\,\boldsymbol{\eta}_{h}\in\mathbb{H}^{\boldsymbol{\rho}}_{h}\right\\},$

and therefore, if we use the finite elements (4.1), we can proceed as in [22,
Section 2.4], to extend the results given in [18] to the case $k\geq 1$, and
define an appropriate space for $\mathbb{H}^{\boldsymbol{t}}_{h}$ as

(4.8) $\displaystyle\mathbb{H}^{\boldsymbol{t}}_{h}$
$\displaystyle:=\left\\{{\boldsymbol{r}}_{h}\in\mathbb{L}^{2}(\Omega):\;\;{\boldsymbol{r}}_{h}|_{K}\in\mathbb{P}_{k}(K)\oplus[\mathbf{B}_{k}(K)]^{d}\oplus([\mathbf{B}_{k}(K)]^{d})^{\mathrm{d}}\quad\forall\,K\in\mathcal{T}_{h}\right\\},$

where $([\mathbf{B}_{k}(K)]^{d})^{\mathrm{d}}$ stands for the deviatoric
tensor of $[\mathbf{B}_{k}(K)]^{d}$. In turn, if (4.5) is employed, we simply
take the part of the AFWk element that approximates $\boldsymbol{\sigma}$
without requiring $\mathbb{H}(\mathbf{div},\Omega)$\- conformity, that is

(4.9) $\displaystyle\mathbb{H}^{\boldsymbol{t}}_{h}$
$\displaystyle:=\left\\{{\boldsymbol{r}}_{h}\in\mathbb{L}^{2}(\Omega):\;\;{\boldsymbol{r}}_{h}|_{K}\in\mathbf{BDM}_{k+1}(K)\quad\forall\,K\in\mathcal{T}_{h}\right\\}.$

Additionally, note that $\mathrm{H}^{\widetilde{p}}_{h}$ does not require any
specific condition. It is therefore simply chosen as

(4.10)
$\mathrm{H}^{\widetilde{p}}_{h}:=\left\\{\widetilde{q}_{h}\in\mathrm{L}^{2}(\Omega):\;\;\widetilde{q}_{h}|_{K}\in\mathrm{P}_{k}(K)\quad\forall\,K\in\mathcal{T}_{h}\right\\}.$

Finally, for pressure and solute concentration we consider Lagrange finite
elements of degree $\leq k+1$, namely

(4.11a)
$\mathrm{H}^{p}_{h}:=\left\\{q_{h}\in\mathrm{C}(\overline{\Omega})\,\cap\,\mathrm{H}^{1}_{\Sigma}(\Omega)\quad
q_{h}|_{K}\in\textrm{P}_{k+1}(K)\;\;\quad\forall\,K\in\mathcal{T}_{h}\right\\},$
(4.11b)
$\mathrm{H}^{\omega}_{h}:=\left\\{\theta_{h}\in\mathrm{C}(\overline{\Omega})\,\cap\,\mathrm{H}^{1}_{\Gamma}(\Omega)\quad\theta_{h}|_{K}\in\textrm{P}_{k+1}(K)\;\;\quad\forall\,K\in\mathcal{T}_{h}\right\\}.$

We now require some preliminary results to establish the well-posedness of the
operator $\mathbf{S}_{h}$.

######  Lemma 4.1.

There exists $\hat{\beta}_{1_{d}}>0$, independent of $h$, such that

(4.12)
$\quad\sup_{\vec{\boldsymbol{\tau}}_{h}\in\mathbf{X}_{h}\setminus\\{\boldsymbol{0}\\}}\frac{[B_{1}(\vec{\boldsymbol{\tau}}_{h}),\vec{\boldsymbol{v}}_{h}]}{\left\lVert\vec{\boldsymbol{\tau}}_{h}\right\rVert_{\mathbf{X}}}\geq\hat{\beta}_{1_{d}}\left\lVert\vec{\boldsymbol{v}}_{h}\right\rVert_{\mathbf{M}}\quad\forall\,\vec{\boldsymbol{v}}_{h}\in\mathbf{M}_{h}.$

######  Lemma 4.2.

There exists $c_{3}>0$, independent of $h$, such that

$c_{3}\|\boldsymbol{\tau}_{h}\|^{2}_{\mathbb{H}(\mathbf{div},\Omega)}\leq\|\boldsymbol{\tau}_{0,h}\|^{2}_{\mathbb{H}(\mathbf{div},\Omega)}\quad\forall\,\boldsymbol{\tau}_{h}\in\mathbb{H}_{h}^{\boldsymbol{\sigma}}.$

######  Proof 4.3.

It follows the same arguments from [18, Lemma 2.4], but now using that
$\boldsymbol{\tau}_{h}\boldsymbol{n}=\boldsymbol{0}$ for all
$\boldsymbol{\tau}_{h}\in\mathbb{H}^{\boldsymbol{\sigma}}_{h}$.

######  Lemma 4.4.

There exists $\hat{\beta}_{d}>0,$ independent of $h$, such that for all

(4.13)
$\quad\sup_{{\boldsymbol{r}}_{h}\in\mathbb{H}^{\boldsymbol{t}}_{h}\setminus\\{\boldsymbol{0}\\}}\frac{[B({\boldsymbol{r}}_{h}),\vec{\boldsymbol{\sigma}}_{h}]}{\left\lVert{\boldsymbol{r}}_{h}\right\rVert_{\mathbb{L}^{2}(\Omega)}}\geq\hat{\beta}_{d}\left\lVert\vec{\boldsymbol{\sigma}}_{h}\right\rVert_{\mathbf{X}}\quad\forall\,\vec{\boldsymbol{\sigma}}_{h}\in\mathbb{V}_{h}.$

######  Proof 4.5.

The proof follows applying Lemmas 3.1 and 4.2, and noticing that
$\boldsymbol{\tau}^{\mathrm{d}}_{h}$ and
$(\widetilde{q}_{h}\mathbb{I}+\boldsymbol{\tau}_{h})\in\mathbb{H}^{\boldsymbol{t}}_{h}$
for each $\vec{\boldsymbol{\tau}}_{h}\in\mathbb{V}_{h}$. We refer to [18,
Lemma 2.5] for further details.

Now, we are in a position to establish the well-definedness of
$\mathbf{T}_{h}$. First we guarantee that the discrete problems defined by
$\mathbf{S}_{h}$ and $\widetilde{\mathbf{S}}_{h}$ are well-posed. In what
follows, we show this result for $\mathbf{S}_{h}$, considering as in the
continuous case, generic functionals in (4.2a)-(4.2d), namely
$\widetilde{H}_{w_{h}},\widetilde{F}_{1,h},\widetilde{F}_{h}$, and
$\widetilde{G}_{h}$ instead of ${H}_{w_{h}},\mathrm{O},{F}$, and ${G}$,
respectively.

######  Lemma 4.6.

For each $\omega_{h}\in\mathrm{H}^{\omega}_{h}$ the problem (4.2a)-(4.2d) has
a unique solution
$(\boldsymbol{t}_{h},\vec{\boldsymbol{\sigma}}_{h},\vec{\boldsymbol{u}}_{h},p_{h})\\\
\in\mathbb{H}^{\boldsymbol{t}}_{h}\times\mathbf{X}_{h}\times\mathbf{M}_{h}\times\mathrm{H}^{p}_{h}$.
Moreover, there exists $C_{d}>0$ independent of $h$ and $\lambda_{s}$, such
that

$\displaystyle\|\mathbf{S}_{h}(\omega_{h})\|$
$\displaystyle=\|\boldsymbol{t}_{h}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}_{h}\|_{\mathbf{X}}+\|\vec{\boldsymbol{u}}_{h}\|_{\mathbf{M}}+\|p_{h}\|_{\mathrm{H}^{1}(\Omega)}$
(4.14) $\displaystyle\leq
C_{d}\left(\|\widetilde{H}_{\omega_{h}}\|_{(\mathbb{H}^{\boldsymbol{t}}_{h})^{\prime}}+\|\widetilde{F}_{1,h}\|_{\mathbf{X}^{\prime}}+\|\widetilde{F}_{h}\|_{\mathbf{M}_{h}^{\prime}}+\|\widetilde{G}_{h}\|_{(\mathrm{H}^{p}_{h})^{\prime}}\right).$

######  Proof 4.7.

It follows similar to the corresponding proof for the continuos case. In fact,
for the well-posedness of the discrete version of problem (3.31), we take
advantage that for finite dimension, it suffices to prove that the homogeneous
problem only has the trivial solution, which follows analogously to the proof
of Lemma 3.10, whereas for the corresponding stability result (4.6), we
proceed exactly as in Lemma 3.12. Therefore, by noticing that $B$ satisfies
the inf-sup condition (4.12), the proof follows the same steps as in Lemma
3.15 but now using the arguments given by [21, Theorem 3.1].

Now we establish the unique solvability of the nonlinear problem (4.2e).

######  Lemma 4.8.

For each $\boldsymbol{\sigma}_{h}\in\mathbb{H}^{\boldsymbol{\sigma}}_{h}$, the
problem (4.2e) has a unique solution
$\omega_{h}:=\widetilde{\mathbf{S}}_{h}(\boldsymbol{\sigma}_{h})\in\mathrm{H}^{\omega}_{h}$.
Moreover, with the same constant $\widetilde{\alpha}$ given in Lemma 3.17,
there holds

(4.15)
$\|\widetilde{\mathbf{S}}_{h}(\boldsymbol{\sigma}_{h})\|:=\|\omega_{h}\|_{\mathrm{H}^{1}(\Omega)}\leq\widetilde{\alpha}^{-1}\|J\|_{(\mathrm{H}^{\omega}_{h})^{\prime}}.$

######  Proof 4.9.

It follows identically as in the proof of Lemma 3.17 but now with elements
living in $\mathrm{H}^{\omega}_{h}$ (cf. (4.11b)).

###  4.2 Discrete solvability analysis

In this section we address the solvability of (4.3). We verify the hypotheses
of the Brouwer fixed-point theorem to prove that (4.3) has at least one fixed
point. Most of the details can be omitted since they follow straightforwardly
by adapting the results given in Section 3.3.

Let
$\mathbf{W}_{h}:=\left\\{\omega_{h}\in\mathrm{H}^{\omega}_{h}:\;\;\left\lVert\omega_{h}\right\rVert_{\mathrm{H}^{1}(\Omega)}\leq
r\right\\}$, be a compact and convex subset of $\mathrm{H}^{\omega}_{h}$, with
$r>0$ defined as in Lemma 3.17. We now provide the discrete analogues of Lemma
3.19.

######  Lemma 4.10.

For the closed ball $\mathbf{W}_{h}$, it holds that
$\mathbf{T}_{h}(\mathbf{W}_{h})\subseteq\mathbf{W}_{h}$.

On the other hand, it is not difficult to prove that for the corresponding
discrete case, the operators $\mathbf{S}$ and $\mathbf{S}_{h}$ satisfy the
bounds given by Lemma 3.21 with a constant $C_{\mathbf{S}_{d}}>0$ instead of
$C_{\mathbf{S}}$, and therefore we have

######  Lemma 4.11.

There exists a constant $C_{\mathbf{T}_{d}}>0$, independent of $h$ and
$\lambda_{s}$ such that

$\|\mathbf{T}_{h}(\omega_{1,h})-\mathbf{T}_{h}(\omega_{2,h})\|_{\mathrm{H}^{1}(\Omega)}\leq
C_{\mathbf{T}_{d}}d_{3}\beta\|\mathbf{T}_{h}(\omega_{2,h})\|_{W^{1,\infty}(\Omega)}\|\omega_{1,h}-\omega_{2,h}\|_{\mathrm{L}^{2}(\Omega)}\quad\forall\,\omega_{1,h},\omega_{2,h}\in\mathrm{H}^{\omega}_{h}.$

Finally, thanks to Lemmas 4.10 and 4.11, and the continuity of $i_{c}$ (cf.
proof of Lemma 3.25) a straightforward application of the Brouwer fixed-point
theorem implies the main result of this section, stated as follows.

######  Theorem 4.12.

Let $\mathbf{W}_{h}$ be defined as at the beginning of this section. Then,
there exists $\overline{C}_{d}>0$ independent of $\lambda_{s}$ and $h$ such
that the coupled problem $\mathrm{(\ref{VFh-poro_1})-\eqref{VFh_diff}}$ has at
least one solution
$(\boldsymbol{t}_{h},\vec{\boldsymbol{\sigma}}_{h},\vec{\boldsymbol{u}}_{h},p_{h},\omega_{h})$
$\in\mathbb{H}^{\boldsymbol{t}}_{h}\times\mathbf{X}_{h}\times\mathbf{M}_{h}\times\mathrm{H}^{p}_{h}\times\mathrm{H}^{\omega}_{h}$
with $\omega_{h}\in\mathbf{W}_{h}$, satisfying the bound

$\|\boldsymbol{t}_{h}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}_{h}\|_{\mathbf{X}}+\|\vec{\boldsymbol{u}}_{h}\|_{\mathbf{M}}+\|p_{h}\|_{\mathrm{H}^{1}(\Omega)}+\|\omega_{h}\|_{\mathrm{H}^{1}(\Omega)}\leq\overline{C}_{d}\left(\|F\|_{\mathbf{M}_{h}^{\prime}}+\|G\|_{(\mathrm{H}^{p}_{h})^{\prime}}+\|J\|_{(\mathrm{H}^{\omega}_{h})^{\prime}}\right).$

######  Proof 4.13.

We proceed similarly as in the continuous case. In fact, from (4.6) and
(4.15), we get

$\displaystyle\|\boldsymbol{t}_{h}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}_{h}\|_{\mathbf{X}}+\|\vec{\boldsymbol{u}}_{h}\|_{\mathbf{M}}+\|p_{h}\|_{\mathrm{H}^{1}(\Omega)}+\|\omega_{h}\|_{\mathrm{H}^{1}(\Omega)}$
$\displaystyle\qquad\qquad\leq\overline{C}_{1,d}\left(\|\widetilde{H}_{\omega_{h}}\|_{(\mathbb{H}^{\boldsymbol{t}}_{h})^{\prime}}+\|\widetilde{F}_{1,h}\|_{\mathbf{X}_{h}^{\prime}}+\|\widetilde{F}_{h}\|_{\mathbf{M}_{h}^{\prime}}+\|\widetilde{G}_{h}\|_{(\mathrm{H}^{p}_{h})^{\prime}}+\|J\|_{\mathrm{H}^{1}_{\Gamma}(\Omega)^{\prime}}\right).$

Then, the result follows by setting
$\widetilde{H}_{\omega_{h}}=H_{\omega_{h}},\widetilde{F}_{1,h}=\mathrm{O},\widetilde{F}_{h}=F$
and $\widetilde{G}_{h}=G$, noticing that
$\|H_{\omega_{h}}\|_{(\mathbb{H}^{\boldsymbol{t}}_{h})^{\prime}}\leq\sqrt{d}\beta\|\omega_{h}\|_{\mathrm{H}^{1}(\Omega)}$,
and applying the fact that $\omega_{h}\in\mathbf{W}_{h}$.

##  5 Error analysis

In this section, we derive the optimal a priori error estimate. For this
purpose, we first establish a Céa estimate formulated in the following
theorem.

######  Theorem 5.1.

Let us
$(\boldsymbol{t},\vec{\boldsymbol{\sigma}},\vec{\boldsymbol{u}},p,\omega)$ and
$(\boldsymbol{t}_{h},\vec{\boldsymbol{\sigma}}_{h},\vec{\boldsymbol{u}}_{h},p_{h},\omega_{h})$
be the solutions of problems (3.3a)-(3.3e) and (4.2a)-(4.2e), respectively,
and assume that

(5.1)
$\hat{C}_{1,d}\beta<\frac{1}{2},\quad\textrm{and}\quad\|\omega\|_{\mathrm{W}^{1,\infty}(\Omega)}\leq\hat{M}:=\min\\{M,\widetilde{M}\\},$

with $M$, $\hat{C}_{1,d}$ and $\widetilde{M}$ specified in (3.45), (5.2) and
(5.10), respectively. Then, there exists $C_{C\textrm{\'{e}}a}>0$ such that

$\displaystyle\|\boldsymbol{t}-\boldsymbol{t}_{h}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}-\vec{\boldsymbol{\sigma}}_{h}\|_{\mathbf{X}}+\|\vec{\boldsymbol{u}}-\vec{\boldsymbol{u}}_{h}\|_{\mathbf{M}}+\|p-p_{h}\|_{\mathrm{H}^{1}(\Omega)}+\|\omega-\omega_{h}\|_{\mathrm{H}^{1}(\Omega)}$
(5.2) $\displaystyle\quad\leq
C_{C\textrm{\'{e}}a}\left(\mathrm{dist}(\boldsymbol{t},\mathbb{H}^{\boldsymbol{t}}_{h})+\mathrm{dist}(\vec{\boldsymbol{\sigma}},\mathbb{H}^{\vec{\boldsymbol{\sigma}}}_{h})+\mathrm{dist}(\vec{\boldsymbol{u}},\mathbf{H}^{\vec{\boldsymbol{u}}}_{h})+\mathrm{dist}(p,\mathrm{H}^{p}_{h})+\mathrm{dist}(\omega,\mathrm{H}^{\omega}_{h})\right).$

######  Proof 5.2.

For sake of notational convenience we define
$\mathbf{e}_{\boldsymbol{t}}=\boldsymbol{t}-\boldsymbol{t}_{h},\,\mathbf{e}_{\vec{\boldsymbol{\sigma}}}=\vec{\boldsymbol{\sigma}}-\vec{\boldsymbol{\sigma}}_{h},\,\mathbf{e}_{\vec{\boldsymbol{u}}}=\vec{\boldsymbol{u}}-\vec{\boldsymbol{u}}_{h}$,
$\mathrm{e}_{p}=p-p_{h},$ $\mathrm{e}_{\omega}=\omega-\omega_{h}$. As usual,
for arbitrary
$(\hat{\boldsymbol{r}_{h}},\hat{\vec{\boldsymbol{\tau}}}_{h},\hat{\vec{\boldsymbol{v}}}_{h},\hat{q}_{h},\hat{\theta}_{h})\in\mathbb{H}^{\boldsymbol{t}}_{h}\times\mathbf{X}_{h}\times\mathbf{M}_{h}\times\mathrm{H}^{p}_{h}\times\mathrm{H}^{\omega}_{h}$,
they are decomposed as

(5.3)
$\mathbf{e}_{{\boldsymbol{t}}}=\boldsymbol{\xi}_{{\boldsymbol{t}}}+\boldsymbol{\chi}_{{\boldsymbol{t}}},\quad\mathbf{e}_{\vec{\boldsymbol{\sigma}}}=\boldsymbol{\xi}_{\vec{\boldsymbol{\sigma}}}+\boldsymbol{\chi}_{\vec{\boldsymbol{\sigma}}},\quad\mathbf{e}_{\vec{\boldsymbol{u}}}=\boldsymbol{\xi}_{\vec{\boldsymbol{u}}}+\boldsymbol{\chi}_{\vec{\boldsymbol{u}}},\quad\mathrm{e}_{{p}}=\xi_{{p}}+\chi_{{p}},\quad\textrm{and}\quad\mathrm{e}_{{\omega}}=\xi_{{\omega}}+\chi_{{\omega}},$

where

$\displaystyle\boldsymbol{\xi}_{{\boldsymbol{t}}}:=\boldsymbol{t}-\hat{\boldsymbol{r}}_{h}\in\mathbb{L}^{2}(\Omega),\
\boldsymbol{\chi}_{{\boldsymbol{t}}}:=\hat{\boldsymbol{r}}_{h}-\boldsymbol{t}_{h}\in\mathbb{H}^{\boldsymbol{t}}_{h},\
\boldsymbol{\xi}_{\vec{\boldsymbol{\sigma}}}:=\vec{\boldsymbol{\sigma}}-\hat{\vec{\boldsymbol{\tau}}}_{h}\in\mathbf{X},\
\boldsymbol{\chi}_{\vec{\boldsymbol{\sigma}}}:=\hat{\vec{\boldsymbol{\tau}}}_{h}-\vec{\boldsymbol{\sigma}}_{h}\in\mathbf{X}_{h},\
\boldsymbol{\chi}_{\vec{\boldsymbol{u}}}:=\hat{\vec{\boldsymbol{v}}}_{h}-\vec{\boldsymbol{u}}_{h}\in\mathbf{M}_{h},$
$\displaystyle\xi_{{p}}:=p-\hat{q}_{h}\in\mathrm{H}^{1}_{\Sigma}(\Omega),\quad\chi_{{p}}:=\hat{q}_{h}-p_{h}\in\mathrm{H}^{p}_{h},\quad\xi_{\omega}:=\omega-\hat{\theta}_{h}\in\mathrm{H}^{1}_{\Gamma}(\Omega),\quad\chi_{\omega}:=\hat{\theta}_{h}-\omega_{h}\in\mathrm{H}^{\omega}_{h}.$

Therefore, by subtracting (4.2a)-(4.2e) and (3.3a)-(3.3e), we easily get the
classical Galerkin orthogonality, which together with the decompositions
(5.3), implies that

$\displaystyle[A({\boldsymbol{\chi}}_{\boldsymbol{t}}),\boldsymbol{r}_{h}]$
$\displaystyle\,+\,$
$\displaystyle[B^{*}({\boldsymbol{\chi}}_{\vec{\boldsymbol{\sigma}}}),\boldsymbol{r}_{h}]$
$\displaystyle\;=\;$
$\displaystyle[H^{h}+H_{\omega-\omega_{h}},\boldsymbol{r}_{h}],$
$\displaystyle[B({\boldsymbol{\chi}}_{\boldsymbol{t}}),\vec{\boldsymbol{\tau}}_{h}]$
$\displaystyle\,-\,$
$\displaystyle[C({\boldsymbol{\chi}}_{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{\tau}}_{h}]$
$\displaystyle\,+\,$
$\displaystyle[B_{1}^{*}({\boldsymbol{\chi}}_{\vec{\boldsymbol{u}}}),\vec{\boldsymbol{\tau}}_{h}]$
$\displaystyle\,+\,$
$\displaystyle[B_{2}^{*}({\chi}_{p}),\vec{\boldsymbol{\tau}}_{h}]$
$\displaystyle\;=\;$
$\displaystyle[\widetilde{F}^{h}_{1},\vec{\boldsymbol{\tau}}_{h}],$ (5.4)
$\displaystyle[B_{1}({\boldsymbol{\chi}}_{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{v}}_{h}]$
$\displaystyle\;=\;$ $\displaystyle[F^{h},\vec{\boldsymbol{v}}_{h}],$
$\displaystyle[B_{2}({\boldsymbol{\chi}}_{\vec{\boldsymbol{\sigma}}}),q_{h}]$
$\displaystyle\,-\,$ $\displaystyle[D(\chi_{p}),q_{h}]$ $\displaystyle\;=\;$
$\displaystyle[G^{h},p_{h}],$
$\displaystyle[\mathscr{A}_{\boldsymbol{\sigma}_{h}}(\chi_{\omega}),\theta_{h}]$
$\displaystyle\;=\;$
$\displaystyle[J^{h}+\mathscr{A}_{\boldsymbol{\sigma}_{h}}(\omega)-\mathscr{A}_{\boldsymbol{\sigma}}(\omega),\theta_{h}],$

for all
$\boldsymbol{r}_{h}\in\mathbb{H}^{\boldsymbol{t}}_{h},\vec{\boldsymbol{\tau}}_{h}\in\mathbf{X}_{h},\vec{\boldsymbol{v}}_{h}\in\mathbf{M}_{h},p_{h}\in\mathrm{H}^{p}_{h},$
and $\theta_{h}\in\mathrm{H}^{\omega}_{h},$ and where

$\displaystyle[H^{h},\boldsymbol{r}_{h}]:=-[A({\boldsymbol{\xi}}_{\boldsymbol{t}}),\boldsymbol{r}_{h}]-[B^{*}({\boldsymbol{\xi}}_{\vec{\boldsymbol{\sigma}}}),\boldsymbol{r}_{h}],\quad[{F}^{h},\vec{\boldsymbol{v}}]:=-[B_{1}({\boldsymbol{\xi}}_{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{v}}_{h}],\quad[G^{h},p_{h}]:=-[B_{2}({\boldsymbol{\xi}}_{\vec{\boldsymbol{\sigma}}}),q_{h}]+[D(\xi_{p}),q_{h}],$
$\displaystyle[J^{h},\theta_{h}]:=-[\mathscr{A}_{\boldsymbol{\sigma}_{h}}(\xi_{p}),\theta_{h}],\quad[\widetilde{F}^{h}_{1},\vec{\boldsymbol{\tau}}]:=-[B({\boldsymbol{\xi}}_{\boldsymbol{t}}),\vec{\boldsymbol{\tau}}_{h}]+[C({\boldsymbol{\xi}}_{\vec{\boldsymbol{\sigma}}}),\vec{\boldsymbol{\tau}}_{h}]-[B_{1}^{*}({\boldsymbol{\xi}}_{\vec{\boldsymbol{u}}}),\vec{\boldsymbol{\tau}}_{h}]-[B_{2}^{*}({\xi}_{p}),\vec{\boldsymbol{\tau}}_{h}].$

Then, proceeding as in the proof of Theorem 4.12 but now for problem (5.4), we
deduce that there exists $\hat{C}_{1,d}>0$, independent of $\lambda_{s}$ and
$h$, such that

$\displaystyle\|\boldsymbol{\chi}_{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}+\|\boldsymbol{\chi}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}+\|\boldsymbol{\chi}_{\vec{\boldsymbol{u}}}\|_{\mathbf{M}}+\|\chi_{p}\|_{\mathrm{H}^{1}(\Omega)}+\|\chi_{\omega}\|_{\mathrm{H}^{1}(\Omega)}$
(5.5)
$\displaystyle\;\;\leq\hat{C}_{1,d}\left(\|H^{h}+H_{\omega-\omega_{h}}\|_{(\mathbb{H}^{\boldsymbol{t}})^{\prime}}+\|\widetilde{F}^{h}_{1}\|_{\mathbf{X}_{h}^{\prime}}+\|F^{h}\|_{\mathbf{M}_{h}^{\prime}}+\|G^{h}\|_{(\mathrm{H}^{p}_{h})^{\prime}}+\|J^{h}+\mathscr{A}_{\boldsymbol{\sigma}_{h}}(\omega)-\mathscr{A}_{\boldsymbol{\sigma}}(\omega)\|_{(\mathrm{H}^{\omega}_{h})^{\prime}}\right).$

In this way, we now proceed to bound the terms on the right-hand side of
(5.2). We begin by noticing thanks to the continuity of the operators $A$ and
$B^{*}$, and the definition of $H_{\omega-\omega_{h}}$ (cf. (3.4c)) that there
holds

(5.6)
$\|H^{h}+H_{\omega-\omega_{h}}\|_{(\mathbb{H}^{\boldsymbol{t}})^{\prime}}\leq\hat{C}_{1}\left(\|\boldsymbol{\xi}_{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}+\|\boldsymbol{\xi}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}\right)+\beta\|\mathrm{e}_{\omega}\|_{\mathrm{H}^{1}(\Omega)},$

with $\hat{C}_{1}$ independent of $\lambda_{s}$ and $h$. Proceeding in a
similar way as before, it can be deduced that

$\displaystyle\|F^{h}\|_{\mathbf{M}_{h}^{\prime}}$
$\displaystyle\leq\hat{C}_{2}\|\boldsymbol{\xi}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}},\qquad\|G^{h}\|_{(\mathrm{H}^{p}_{h})^{\prime}}\leq\hat{C}_{3}\left(\frac{\alpha}{\lambda_{s}}\|\boldsymbol{\xi}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}+\|\xi_{p}\|_{\mathrm{H}^{1}(\Omega)}\right),$
(5.7) $\displaystyle\|\widetilde{F}^{h}_{1}\|_{\mathbf{X}_{h}^{\prime}}$
$\displaystyle\leq\hat{C}_{4}\left(\|\boldsymbol{\xi}_{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}+\frac{1}{\lambda_{s}}\|\boldsymbol{\xi}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}+\|\boldsymbol{\xi}_{\vec{\boldsymbol{u}}}\|_{\mathbf{M}}+\frac{\alpha}{\lambda_{s}}\|\xi_{p}\|_{\mathrm{H}^{1}(\Omega)}\right),$

where $\hat{C}_{2},\hat{C}_{3},\hat{C}_{4}>0$ are independent of $\lambda_{s}$
and $h$. On the other hand, in order to bound the last term on the right-hand
side of (5.2), we proceed as in (3.22) to get

$\|\mathscr{A}_{\boldsymbol{\sigma}_{h}}(\omega)-\mathscr{A}_{\boldsymbol{\sigma}}(\omega)\|_{(\mathrm{H}^{\omega}_{h})^{\prime}}\leq
d_{3}\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{\mathbb{L}^{2}(\Omega)}\|\omega\|_{\mathrm{W}^{1,\infty}(\Omega)},$

and then, from the latter, we conclude that

(5.8)
$\|J^{h}+\mathscr{A}_{\boldsymbol{\sigma}_{h}}(\omega)-\mathscr{A}_{\boldsymbol{\sigma}}(\omega)\|_{(\mathrm{H}^{\omega}_{h})^{\prime}}\leq\hat{C}_{5}\|\xi_{\omega}\|_{\mathrm{H}^{1}(\Omega)}+d_{3}\|\mathbf{e}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}\|\omega\|_{\mathrm{W}^{1,\infty}(\Omega)},$

with $\hat{C}_{5}$ independent of $\lambda_{s}$ and $h$. Therefore, from
(5.6), (5.2) and $\eqref{error_J}$, we obtain that there exists
$\hat{C}_{2,d}>0$ independent of $\lambda_{s}$ and $h$, such that

$\displaystyle\|\boldsymbol{\chi}_{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}+\|\boldsymbol{\chi}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}+\|\boldsymbol{\chi}_{\vec{\boldsymbol{u}}}\|_{\mathbf{M}}+\|\chi_{p}\|_{\mathrm{H}^{1}(\Omega)}+\|\chi_{\omega}\|_{\mathrm{H}^{1}(\Omega)}$
$\displaystyle\;\;\leq\hat{C}_{2,d}\left(\|\boldsymbol{\xi}_{\boldsymbol{t}}\|_{\mathbb{L}^{2}(\Omega)}+\left(1+\frac{1}{\lambda_{s}}+\frac{\alpha}{\lambda_{s}}\right)\|\boldsymbol{\xi}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}+\|\boldsymbol{\xi}_{\vec{\boldsymbol{u}}}\|_{\mathbf{M}}+\left(1+\frac{\alpha}{\lambda_{s}}\right)\|\xi_{p}\|_{\mathrm{H}^{1}(\Omega)}+\|\xi_{\omega}\|_{\mathrm{H}^{1}(\Omega)}\right)$
(5.9)
$\displaystyle\;\;\quad+d_{3}\|\mathbf{e}_{\vec{\boldsymbol{\sigma}}}\|_{\mathbf{X}}\widetilde{M}+\hat{C}_{1,d}\beta\|\mathrm{e}_{\omega}\|_{\mathrm{H}^{1}(\Omega)},$

with $\|\omega\|_{\mathrm{W}^{1,\infty}(\Omega)}\leq\widetilde{M}$, and where
$\left(1+\frac{1}{\lambda_{s}}+\frac{\alpha}{\lambda_{s}}\right)$ and
$\left(1+\frac{\alpha}{\lambda_{s}}\right)$ can be seen as constants
independent of $\lambda_{s}$ if $\lambda_{s}\rightarrow\infty$. In this way,
choosing $\widetilde{M}$ such that

(5.10) $\widetilde{M}\leq\frac{1}{2d_{3}},$

the desired result follows simply by triangle’s inequality in (5.3), the
estimation (5.2), and assumption (5.1).

The main result of this section is given by the following lemma.

######  Theorem 5.3.

In addition to the hypotheses of Theorems 3.29, 4.12 and 5.1, assume that
there exist $s>0$ and $l>1/2$ such that
$\boldsymbol{\sigma}\in\mathbb{H}^{s}(\Omega)$,
$\mathbf{div}\,\boldsymbol{\sigma}\in\mathbf{H}^{s}(\Omega),\,\boldsymbol{u}\in\mathbf{H}^{s}(\Omega),\,\boldsymbol{\gamma}\in\mathbb{H}^{s}(\Omega),\,\widetilde{p}\in\mathrm{H}^{s}(\Omega),\,p\in\mathrm{H}^{1+s}(\Omega),\,\omega\in\mathrm{H}^{1+s}(\Omega),$
and $\boldsymbol{t}\in\mathbb{H}^{l}(\Omega)$. Then, there exist
$\widetilde{C}_{1},\widetilde{C}_{2}>0$, independent of $h$, such that, with
the finite element subspaces defined by $\mathrm{(\ref{discrete-spaces-
elasticity-AFW})}$, $\mathrm{(\ref{FEM_T1})}$, $\mathrm{(\ref{FEM_PT1})}$,
$\mathrm{(\ref{FE_pressure})},$ and $\mathrm{(\ref{FE_concentration})},$ there
holds

(5.11)
$\displaystyle\begin{split}&\|\boldsymbol{t}-\boldsymbol{t}_{h}\|_{\mathbb{L}^{2}(\Omega)}+\|\vec{\boldsymbol{\sigma}}-\vec{\boldsymbol{\sigma}}_{h}\|_{\mathbf{X}}+\|\vec{\boldsymbol{u}}-\vec{\boldsymbol{u}}_{h}\|_{\mathbf{M}}+\|p-p_{h}\|_{\mathrm{H}^{1}(\Omega)}+\|\omega-\omega_{h}\|_{\mathrm{H}^{1}(\Omega)}\\\
&\qquad\leq\widetilde{C}_{1}h^{\mathrm{min}\\{s,k+1\\}}\,\left\\{\left\lVert\boldsymbol{\sigma}\right\rVert_{\mathbb{H}^{s}(\Omega)}+\,\,\left\lVert\mathbf{div}\,\boldsymbol{\sigma}\right\rVert_{s,\Omega}\,+\,\left\lVert\boldsymbol{u}\right\rVert_{s,\Omega}+\left\lVert\boldsymbol{\gamma}\right\rVert_{\mathbb{H}^{s}(\Omega)}+\|\widetilde{p}\|_{\mathrm{H}^{s}(\Omega)}\right.\\\
&\qquad\quad\left.+\,\|p\|_{\mathrm{H}^{s}(\Omega)}+\left\lVert\omega\right\rVert_{\mathrm{H}^{s+1}(\Omega)}\right\\}+\widetilde{C}_{2}h^{\mathrm{min}\\{l,k+1\\}}\left\lVert\boldsymbol{t}\right\rVert_{\mathbb{H}^{l}(\Omega)}.\end{split}$

######  Proof 5.4.

It follows as a combination of the approximation properties of the spaces
(4.5), (4.9), (4.10), (4.11a), and (4.11b) (see, e.g., [5, 17]), the Céa
estimate (5.1) and the smallness assumption.

######  Remark 5.5.

Similar results are obtained when (4.1) and (4.8) are employed instead of
(4.5) and (4.9), respectively. For the lowest-order case, we notice from the
classical approximation properties of the finite element subspaces, that it
suffices to consider $l>0$, and therefore, the last term in (LABEL:rates)
should be changed by
$\widetilde{C}_{1}h^{\mathrm{min}\\{s,k+1\\}}\sum_{K\in\mathcal{T}_{h}}\left\lVert\boldsymbol{t}\right\rVert_{\mathbb{H}^{s}(K)}$
(see [18, Theorem 3.2]), whereas the other terms remain unchanged. Even if
approximation results for (4.8) in the high-order case do not seem to be
available in the literature, the numerical results in the next section
demonstrate an asymptotic $\mathcal{O}(h^{k+1})$ convergence (see Table 6.1).

##  6 Numerical tests

###  6.1 Example 1: Convergence against manufactured solutions

The accuracy of the discretisation is verified through two convergence tests
using manufactured solutions, and implemented using the FEniCS library [33].
The model parameters are simply taken as
$\mu_{s}=\lambda_{s}=1=c_{0}=\kappa=\alpha=\ell=\rho_{f}=\rho_{s}=\mu_{f}=\beta=1$,
and $\phi=0.5$, and $\boldsymbol{g}=\boldsymbol{0}$. Variations (of several
orders of magnitude) in $\lambda_{s},\kappa,\alpha,\mu_{f}$ will be also
considered to study the robustness of the formulation with respect to these
parameters. The stress-altered diffusion is (2.3) with $D_{0}=0.01$ and
$\eta=0.01$. We use the following closed-form smooth solutions to (2.4)

$\displaystyle\omega=e^{x_{1}}+\cos(\pi x_{1})\cos(\pi x_{2}),\
\boldsymbol{u}=\frac{1}{10}\begin{pmatrix}-\cos(x_{1})\sin(x_{2})+\frac{x_{1}^{2}}{\lambda_{s}}\\\
\sin(x_{1})\cos(x_{2})+\frac{x_{2}^{2}}{\lambda_{s}}\end{pmatrix},\ p=\sin(\pi
x_{1})\sin(\pi x_{2}),\ \widetilde{p}=\alpha
p-\lambda_{s}\operatorname*{div}\boldsymbol{u},$

which are used to produce non-homogeneous forcing and source terms. For sake
of simplicity only one type of boundary conditions is considered, i.e.,
$\partial\Omega=\Gamma$.The remaining boundary conditions in (2.4c) (for
displacement and fluid flux) are imposed naturally. This also implies that the
stress is not uniquely defined and we constraint its trace, using a real
Lagrange multiplier. We generate seven successively refined meshes (congruent
right-angled triangular partitions) for the domain $\Omega=(0,1)^{2}$ and
calculate errors between approximate and exact solutions $e(\cdot)$ (measured
in $\mathrm{H}^{1}$-norm for tracer and fluid pressure, in the tensor
$\mathbf{H}(\operatorname*{div})$-norm for stress, and in the tensorial,
vectorial and scalar $\mathrm{L}^{2}$-norm for strain, rotation, displacement
and total pressure). The mixed finite element methods are defined by the
PEERSk and enriched piecewise polynomial spaces specified in
(4.1),(4.8),(4.11a).

DoF | $e_{0}(\boldsymbol{t})$ | rate | $e_{\operatorname*{\mathbf{div}}}(\boldsymbol{\sigma})$ | rate | $e_{1}(\tilde{p})$ | rate | $e_{0}(\boldsymbol{u})$ | rate | $e_{0}(\boldsymbol{\gamma})$ | rate | $e_{1}(p)$ | rate | $e_{1}(\omega)$ | rate
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
scheme with $k=0$
148 | 1.9e-01 | – | 1.2e+0 | – | 2.9e-01 | – | 2.8e-02 | – | 3.3e-01 | – | 1.1e+0 | – | 1.6e+0 | –
540 | 8.7e-02 | 1.15 | 6.2e-01 | 0.99 | 1.4e-01 | 1.03 | 9.0e-03 | 1.64 | 9.8e-02 | 1.76 | 7.4e-01 | 0.62 | 8.9e-01 | 0.87
2068 | 4.2e-02 | 1.05 | 3.1e-01 | 1.00 | 6.7e-02 | 1.06 | 3.9e-03 | 1.19 | 3.9e-02 | 1.41 | 4.1e-01 | 0.85 | 4.6e-01 | 0.97
8100 | 2.1e-02 | 1.03 | 1.6e-01 | 1.00 | 3.3e-02 | 1.03 | 1.9e-03 | 1.06 | 1.8e-02 | 1.35 | 2.1e-01 | 0.95 | 2.3e-01 | 0.99
32068 | 1.0e-02 | 1.02 | 7.8e-02 | 1.00 | 1.6e-02 | 1.01 | 9.4e-04 | 1.01 | 7.8e-03 | 1.18 | 1.1e-01 | 0.98 | 1.2e-01 | 1.00
127620 | 5.0e-03 | 1.01 | 3.9e-02 | 1.00 | 8.2e-03 | 1.00 | 4.7e-04 | 1.00 | 3.7e-03 | 1.10 | 5.4e-02 | 1.00 | 5.8e-02 | 1.00
509188 | 2.5e-03 | 1.01 | 1.9e-02 | 1.00 | 4.1e-03 | 1.00 | 2.3e-04 | 1.00 | 1.5e-03 | 1.20 | 2.7e-02 | 1.00 | 2.9e-02 | 1.00
scheme with $k=1$
436 | 1.3e-02 | – | 2.3e-01 | – | 7.7e-02 | – | 1.8e-03 | – | 6.9e-03 | – | 3.8e-01 | – | 4.8e-01 | –
1652 | 4.6e-03 | 1.51 | 5.8e-02 | 1.97 | 2.0e-02 | 1.95 | 4.4e-04 | 2.02 | 2.4e-03 | 1.35 | 1.2e-01 | 1.70 | 1.3e-01 | 1.88
6436 | 1.3e-03 | 1.76 | 1.5e-02 | 1.99 | 5.0e-03 | 1.99 | 1.0e-04 | 2.07 | 8.7e-04 | 1.44 | 3.2e-02 | 1.88 | 3.3e-02 | 1.96
25412 | 3.6e-04 | 1.91 | 3.7e-03 | 2.00 | 1.2e-03 | 2.00 | 2.5e-05 | 2.04 | 2.5e-04 | 1.78 | 8.2e-03 | 1.95 | 8.4e-03 | 1.99
100996 | 9.2e-05 | 1.96 | 9.2e-04 | 2.00 | 3.1e-04 | 2.00 | 6.3e-06 | 2.01 | 6.7e-05 | 1.91 | 2.1e-03 | 1.98 | 2.1e-03 | 2.00
402692 | 2.3e-05 | 1.98 | 2.3e-04 | 2.00 | 7.8e-05 | 2.00 | 1.6e-06 | 2.00 | 1.7e-05 | 1.96 | 5.2e-04 | 1.99 | 5.3e-04 | 2.00
1608196 | 5.9e-06 | 1.99 | 5.7e-05 | 2.00 | 1.9e-05 | 2.00 | 3.9e-07 | 2.00 | 4.4e-06 | 1.98 | 1.3e-04 | 2.00 | 1.3e-04 | 2.00
Table 6.1: Verification of convergence for the method with PEERSk and enriched
piecewise polynomials with degrees $k=0$ and $k=1$, and using unity
parameters. Errors and convergence rates are tabulated for strain, stress,
total pressure, displacement, rotation Lagrange multiplier, fluid pressure,
and tracer concentration.

Table 6.1 shows such error decay for polynomial orders $k=0$ and $k=1$, and we
can clearly observe a convergence of $O(h^{k+1})$ for all field variables in
their natural norms, which is consistent with the theoretical error bounds.
Table 6.2 shows the error history associated with the extension to 3D, using
the exact solutions

$\displaystyle\omega=\cos(\pi x_{1})\cos(\pi x_{2})\cos(\pi x_{3}),\
\boldsymbol{u}=\frac{1}{10}\begin{pmatrix}\sin(x_{1})\cos(x_{2})\cos(x_{3})+\frac{x_{1}^{2}}{\lambda_{s}}\\\
-2\cos(x_{1})\sin(x_{2})\cos(x_{3})+\frac{x_{2}^{2}}{\lambda_{s}}\\\
\cos(x_{1})\cos(x_{2})\sin(x_{3})+\frac{x_{3}^{2}}{\lambda_{s}}\end{pmatrix},\
p=\sin(\pi x_{1})\sin(\pi x_{2})\sin(\pi x_{3}),$

and $\widetilde{p}=\alpha p-\lambda_{s}\operatorname*{div}\boldsymbol{u}$, on
$\Omega=(0,1)^{3}$, with the same unit parameters, and focusing on the lowest-
order scheme. Optimal, first-order convergence is also achieved in this case.
For completeness, we also verify the error decay of the scheme resulting from
using AFWk elements (4.5) and (4.9), which also delivers optimal convergence.

DoF | $e_{0}(\boldsymbol{t})$ | rate | $e_{\operatorname*{\mathbf{div}}}(\boldsymbol{\sigma})$ | rate | $e_{1}(\tilde{p})$ | rate | $e_{0}(\boldsymbol{u})$ | rate | $e_{0}(\boldsymbol{\gamma})$ | rate | $e_{1}(p)$ | rate | $e_{1}(\omega)$ | rate
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
scheme with PEERSk and enriched piecewise polynomials
1984 | 1.7e-01 | – | 1.5e+0 | 0.00 | 2.0e-01 | – | 2.0e-02 | – | 1.4e-01 | – | 1.4e+0 | – | 1.6e+0 | –
6477 | 9.8e-02 | 0.72 | 1.0e+0 | 0.87 | 1.5e-01 | 0.74 | 1.3e-02 | 0.99 | 8.3e-02 | 1.29 | 9.4e-01 | 0.63 | 1.2e+0 | 0.73
29281 | 6.4e-02 | 0.83 | 6.4e-01 | 0.95 | 8.9e-02 | 0.98 | 7.7e-03 | 1.08 | 2.7e-02 | 2.16 | 6.7e-01 | 0.76 | 7.7e-01 | 0.90
168297 | 3.5e-02 | 0.92 | 3.6e-01 | 0.99 | 4.7e-02 | 1.09 | 4.1e-03 | 1.07 | 1.0e-02 | 1.71 | 4.1e-01 | 0.84 | 4.3e-01 | 0.98
1125049 | 2.1e-02 | 0.95 | 1.9e-01 | 1.00 | 2.4e-02 | 1.06 | 2.1e-03 | 1.04 | 3.7e-03 | 1.59 | 2.2e-01 | 0.96 | 2.3e-01 | 1.00
8194137 | 1.2e-02 | 0.96 | 9.6e-02 | 1.01 | 1.3e-02 | 1.02 | 1.1e-03 | 1.00 | 1.6e-03 | 1.49 | 1.2e-01 | 0.97 | 1.2e-01 | 0.98
scheme with AFWk elements
2551 | 7.1e-02 | – | 1.4e+0 | – | 1.9e-01 | – | 1.9e-02 | – | 3.8e-02 | – | 1.1e+0 | – | 1.6e+0 | –
8067 | 4.3e-02 | 1.24 | 1.0e+0 | 0.90 | 1.4e-01 | 0.78 | 1.2e-02 | 1.04 | 2.2e-02 | 1.38 | 9.4e-01 | 0.43 | 1.2e+0 | 0.74
35383 | 2.3e-02 | 1.23 | 6.1e-01 | 0.96 | 8.5e-02 | 0.99 | 7.4e-03 | 1.02 | 1.1e-02 | 1.35 | 6.7e-01 | 0.66 | 7.7e-01 | 0.91
198831 | 1.0e-02 | 1.34 | 3.4e-01 | 0.99 | 4.6e-02 | 1.06 | 4.0e-03 | 1.02 | 5.0e-03 | 1.32 | 4.1e-01 | 0.84 | 4.3e-01 | 0.98
1310431 | 4.6e-03 | 1.31 | 1.8e-01 | 1.00 | 2.4e-02 | 1.03 | 2.1e-03 | 1.01 | 2.3e-03 | 1.22 | 2.3e-01 | 0.94 | 2.3e-01 | 0.99
Table 6.2: Verification of space convergence in 3D, with polynomial degree
$k=0$, and using unity parameters. Errors history for strain, stress, total
pressure, displacement, rotation Lagrange multiplier, fluid pressure, and
tracer concentration.

Figure 6.1 illustrates the robustness with respect to large variations in
physical parameters including nearly incompressible materials (taking
$\lambda_{s}=10^{8}$), for low permeability (with $\kappa=10^{-12}$), for weak
Biot-Willis coupling ($\alpha=10^{-6}$), and with small fluid viscosity (with
$\mu_{f}=10^{-4}$). Analogous results (not shown here) were obtained for zero
storage coefficient, and for different values of porosity. For all these tests
(both 2D and 3D), in average, four Newton-Raphson iterations were required to
reach convergence across all mesh refinements.

Fig. 6.1: Verification of space convergence and robustness with respect to
changes in model parameters $\lambda_{s},\kappa,\alpha,\mu_{f}$.

###  6.2 Example 2: Testing the effect of stress-hindered diffusion

In this test we compare the filtration effects on models using different
degrees of stress-modified diffusion where the diffusion due to stress is
lower than the constant base-line diffusion. This simple test emphasises the
significance of anisotropy and heterogeneity in the transport of tracer and it
is more illustrative to consider the time-dependent case (adding a time
derivative of the first two terms in the second equation of (2.1) and of the
first term in (2.2), which we discretise using backward Euler’s method with a
constant time step. For this we employ a slab of tissue of size 1 mm2, and the
model parameters are as follows (see, e.g., [25, 40, 8])

$E=800\,\text{Pa},\quad\nu=0.495,\quad c_{0}=2\times
10^{-8},\quad\kappa=10^{-8}\,\text{mm}^{2},\quad\alpha=1,\quad\ell=0,\quad\rho_{s}=10^{-3}\,\text{mg/mm}^{3},$
$\mu_{f}=0.7\,\text{Pa s},\quad\beta=0.45,\quad\phi=0.2,\quad D_{0}=5.3\times
10^{-5}\text{mm}^{2}/\text{s},$

where $D_{0}$ is made sufficiently large to compare with stress-hindered
effects. On the top of the slab we impose a traction
$\boldsymbol{\sigma}\boldsymbol{n}=-\alpha p_{\mathrm{top}}\boldsymbol{n}$
with $p_{\mathrm{in}}=0.5\,\mathrm{atan}(t/10)$ mmHg/mm, on the bottom we
clamp the tissue $\boldsymbol{u}=\boldsymbol{0}$, and on the vertical walls we
set zero traction. The fluid pressure $p_{\mathrm{in}}$ is prescribed on the
top segment and set to $p_{0}=9$ mmHg/mm on the vertical walls, whereas we set
zero normal flux on the bottom. The tracer concentration
$\omega_{\mathrm{in}}=1$ mmol is fixed on the top and zero diffusive flux is
considered elsewhere on the boundary. A coarse mesh with 16384 triangular
elements is employed and the time step is $\Delta t=50$ s.

Fig. 6.2: Example 2. Stress-assisted diffusion. Top: Filtration of the CSF
tracer into the slab of tissue for pure reaction-diffusion (left), and stress-
hindered diffusion according to (6.1) (centre), and the anisotropic form (6.2)
(right). Middle row: approximate displacement (left), total pressure (centre),
and fluid pressure (right). Bottom: tensor glyph representation on the
deformed configuration, of the total Cauchy stress (left), the isotropic
diffusion (centre), and anisotropic one (right). All snapshots are taken at
time $t=1800$ s.

Three main scenarios are considered: (a) Pure reaction-diffusion of the tracer
with $D_{0}$; (b) Adding the effect of isotropic stress-assisted diffusion
according to

(6.1)
${D}(\boldsymbol{\sigma})=D_{0}\mathbb{I}+D_{0}\exp(-\eta_{0}\operatorname*{tr}\boldsymbol{\sigma})\mathbb{I},$

with $\eta_{0}=5\times 10^{-5}$; and (c) Modifying the functional form of the
apparent diffusion to

(6.2)
$\tilde{D}(\boldsymbol{\sigma})=\eta_{0}D_{0}\mathbb{I}-\eta_{2}D_{0}\boldsymbol{\sigma}+\eta_{2}D_{0}\boldsymbol{\sigma}^{2},$

with $\eta_{1}=0.02$ and $\eta_{2}=10^{-5}$. Figure 6.2 shows the results from
these tests. The first row illustrates how with pure diffusion (without
stress-hinderance) the tracer penetrates faster than in the cases with stress-
dependent diffusion (this is shown at $t=1800$ s). The deformation of the
parenchymal slab using case (c) is shown in the second row, where we can
observe a localisation of total pressure accumulation near the bottom corners
of the domain. The bottom row indicates the degree of anisotropy of the total
poroelastic stress and of the stress-altered diffusivity tensors according to
(6.1) and (6.2), also at $t=1800$ s. In case (c), apart from a slower
diffusion than in case (a), a slight deviation from plane solute transport is
seen (in the top-right panel of the figure), explained by the anisotropy
differences exhibited in the bottom row of the figure.

###  6.3 Example 3: Stress-hindered diffusion and transport in the brain

Next we perform two application tests investigating the filtration properties
of parenchymal brain tissue and the evolution of tracer concentration in its
sleeping vs awake state. In the first example, we create a 2D slab (1 mm deep,
1 $\mu$m wide) of a mouse brain, and simulated tracer enrichment from the
cortical surface into the cortex. The diffusion coefficient was here assumed
to take the form

(6.3) $D=D_{0}-D_{0}\exp(-\eta|\operatorname*{tr}(\boldsymbol{\sigma})|).$

Fig. 6.3: Example 3. Tracer infiltration in the mouse cortex (top panel) and
tracer clearance from the human brain (bottom panel) using the stress-hindered
formulation as described by (6.3). In the top-left we compare the sleeping
versus awake state after 15 minutes of transport from the mouse cortex 300
$mu$m into the tissue. The concentration was held fixed at the cortex. To the
right, the shown time evolution 100 $\mu$m below the cortex indicates that
transport by stress-hindered diffusion is an order of magnitude lower than for
free diffusion within the cortex. In the bottom panel we show the tracer
distribution in a brain slice of a human brain. In the left panel is the
tracer distribution after 12 hours of diffusive transport. In the awake state,
stress hindered diffusion results in slightly slower clearance. In the right
bottom panel, the average concentration for the entire piece of brain tissue
is plotted over time.

At the cortex we set a given concentration of $\omega=1$, while other
boundaries were assigned homogeneous Neumann conditions and no initial tracer
was assumed within the tissue. The depth (1 mm) was assumed to cover the
entire width of the mouse cortex. Over time we measured the concentration (100
$\mu$m into the tissue) similarly to the results reported by [44]. In the
awake state, the porosity was set to $\phi=0.14$, while in the awake state we
used $\phi=0.23$ [44]. We assumed that the known change in volume fraction was
associated with a similar change in the width of the cerebral cortex, which is
known to change during sleep [14]. Simulating the (compressed) awake state, we
therefore added a force
$\boldsymbol{\sigma}\boldsymbol{n}=-p_{0}\arctan(\frac{t}{50s})\boldsymbol{n}$,
where $p_{0}=0.5$ mmHg to the top of the domain, while the bottom was fixed.
The arctan function was used to ensure a smooth loading. In the sleeping
state, no forces were applied, and tracers were thus allowed to diffuse freely
according to the base diffusion $D_{0}$. In both cases the sidewalls were
assigned Neumann conditions. The resulting displacement at the cortex surface
was 0.09 mm (data now shown). In Figure 6.3 (top panel) we show the resulting
tracer concentration. We compare tracer concentrations after 15 minutes in the
awake versus the sleeping state, from 0 to 300 $\mu$m into the cortex (left),
and show the time evolution of tracer concentration at a slice 100 $\mu$m
below the surface (right). The applied forces in the awake state significantly
slows down transport into mouse cortex. Other parameters were here set as
follows (see also [43, 8])

$E=800\,\mathrm{Pa},\quad\nu=0.495,\quad c_{0}=2\times
10^{-8}\,\text{Pa}^{-1},\quad\kappa=10^{-8}\,\text{mm}^{2},\quad\alpha=1,\quad\ell=0,\quad\eta=2\times
10^{-5},$ $\rho_{s}=10^{-3}\,\text{g/mm}^{3},\quad\mu_{f}=0.7\times
10^{-3}\,\text{Pa}/\text{s},\quad\beta=0.35,\quad\quad D_{0}=5.3\times
10^{-2}\text{mm}^{3}/\text{s}.$

In the second scenario, we simulate tracer clearance from the human brain
using a 3D mesh generated from the the right hemisphere of the built-in
subject Bert in FreeSurfer [16]. The mesh consisted of 400,000 cells, and the
resulting system had 15M degrees of freedom. Here we have used an initial
condition of $\omega=0$ everywhere in the brain, while we set $\omega=0$ on
the brain surface. Parameters were set exactly as in the previous example,
except that now $\eta=2\times 10^{-1}$ and the cortex pressure is $p_{0}=0.01$
mmHg. On the ventricles we imposed no-slip conditions. In Figure 6.3 (bottom
panel) we show the tracer concentration in the sleeping versus awake brain
after 12 hours of transport (left), and also show the time evolution of the
average tracer concentration within the slice shown (right). Transport rates
differ between the sleeping versus the awake state, but occurs on
approximately the same time scale [42] in contrast to the mouse test case. Our
results from the human brain test case thus reflects differences in clearance
in the awake versus sleeping state observed experimentally [13].

We close this section mentioning that next steps in the process of model
refinement include the study of interfacial flow between the poroelastic
structure (at the macroscale) coupled with the surrounding layer of CSF and
the specific waste disposal system, which is still not well understood. For
this we will also use multiple network models [26, 37], the preconditioning
framework for perturbed saddle-point problems from [6], and interface
couplings such as those in [9, 32, 41]. Other envisaged generalisations deal
with large deformations and kinematics through a more complete mixture theory
following, e.g., [30, 10].

## References

  * [1] I. Ambartsumyan, E. Khattatov, and I. Yotov, A coupled multipoint stress - multipoint flux mixed finite element method for the Biot system of poroelasticity, Computer Methods in Applied Mechanics and Engineering, 372 (2020), p. e113407.
  * [2] D. Arnold, F. Brezzi, and J. Douglas, PEERS: A new mixed finite element method for plane elasticity, Japan Journal of Applied Mathematics, 1 (1984), pp. 347–367.
  * [3] D. Arnold, R. Falk, and R. Winther, Mixed finite element methods for linear elasticity with weakly imposed symmetry, Mathematical of Computations, 76 (2007), pp. 1699–1723.
  * [4] T. Bærland, J. J. Lee, K.-A. Mardal, and R. Winther, Weakly imposed symmetry and robust preconditioners for Biot’s consolidation model, Computational Methods in Applied Mathematics, 17 (2017), pp. 377–396.
  * [5] D. Boffi, F. Brezzi, and M. Fortin, Mixed Finite Element Methods and Applications, vol. 44, Springer-Verlag, Berlin, 2013.
  * [6] W. Boon, M. Kuchta, K.-A. Mardal, and R. Ruiz-Baier, Robust preconditioners and stability analysis for perturbed saddle-point problems – application to conservative discretizations of Biot’s equations utilizing total pressure, SIAM Journal on Scientific Computing, 43 (2021), pp. B961–B983.
  * [7] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, Berlin, 1991.
  * [8] S. Budday, R. Nay, R. de Rooij, P. Steinmann, T. Wyrobek, T. C. Ovaert, and E. Kuhl, Mechanical properties of gray and white matter brain tissue by indentation, Journal of the mechanical behavior of biomedical materials, 46 (2015), pp. 318–330.
  * [9] M. Bukač, I. Yotov, R. Zakerzadeh, and P. Zunino, Partitioning strategies for the interaction of a fluid with a poroelastic material based on a Nitsche’s coupling approach, Computer Methods in Applied Mechanics and Engineering, 292 (2015), pp. 138–170.
  * [10] F. Costanzo and S. T. Miller, An arbitrary Lagrangian-Eulerian finite element formulation for a poroelasticity problem stemming from mixture theory, Computer Methods in Applied Mechanics and Engineering, 323 (2017), pp. 64–97.
  * [11] M. Croci, V. Vinje, and M. E. Rognes, Uncertainty quantification of parenchymal tracer distribution using random diffusion and convective velocity fields, Fluids and Barriers of the CNS, 16 (2019), pp. 1–21.
  * [12] , Fast uncertainty quantification of tracer distribution in the brain interstitial fluid with multilevel and quasi Monte Carlo, International Journal for Numerical Methods in Biomedical Engineering, 37 (2020), p. e3412.
  * [13] P. K. Eide, V. Vinje, A. H. Pripp, K.-A. Mardal, and G. Ringstad, Sleep deprivation impairs molecular clearance from the human brain, Brain, 144 (2021), pp. 863–874.
  * [14] T. Elvsåshagen, N. Zak, L. B. Norbom, P. Ø. Pedersen, S. H. Quraishi, A. Bjørnerud, D. Alnæs, N. T. Doan, U. F. Malt, I. R. Groote, et al., Evidence for cortical structural plasticity in humans after a day of waking and sleep deprivation, Neuroimage, 156 (2017), pp. 214–223.
  * [15] A. Elyes, F. Radu, and J. Nordbotten, Adaptive poromechanics computations based on a posteriori error estimates for fully mixed formulations of Biot’s consolidation model, Computer Methods in Applied Mechanics and Engineering, 347 (2019), pp. 264–294.
  * [16] B. Fischl, Freesurfer, Neuroimage, 62 (2012), pp. 774–781.
  * [17] G. N. Gatica, A Simple Introduction to the Mixed Finite Element Method. Theory and Applications, Springer-Verlag, Berlin, 2014\.
  * [18] G. N. Gatica and L. F. Gatica, On the a priori and a posteriori error analysis of a two-fold saddle-point approach for nonlinear incompressible elasticity, International Journal for Numerical Methods in Engineering, 68 (2006), pp. 861–892.
  * [19] G. N. Gatica, B. Gómez-Vargas, and R. Ruiz-Baier, Analysis and mixed-primal finite element discretisations for stress-assisted diffusion problems, Computer Methods in Applied Mechanics and Engineering, 337 (2018), pp. 411–438.
  * [20] , Formulation and analysis of fully-mixed finite element methods for stress-assisted diffusion problems, Computers & Mathematics with Applications, 77 (2019), pp. 1312–1330.
  * [21] G. N. Gatica, N. Heuer, and S. Meddahi, On the numerical analysis of nonlinear twofold saddle point problems, IMA J. Numer. Anal, 23 (2003), pp. 301–330.
  * [22] G. N. Gatica, A. Márquez, and W. Rudolph, A priori and a posteriori error analyses of augmented twofold saddle point formulations for nonlinear elasticity problems, Computer Methods in Applied Mechanics and Engineering, 264 (2013), pp. 23–48.
  * [23] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, Springer Publishing Company, 1st ed., 1986.
  * [24] P. Grigoreva, E. N. Vilchevskaya, and W. H. Müller, Stress and Diffusion Assisted Chemical Reaction Front Kinetics in Cylindrical Structures, Springer International Publishing, Cham, 2019, pp. 53–72.
  * [25] K. E. Holter, B. Kehlet, A. Devor, T. J. Sejnowski, A. M. Dale, S. W. Omholt, O. P. Ottersen, E. A. Nagelhus, K.-A. Mardal, and K. H. Pettersen, Interstitial solute transport in 3D reconstructed neuropil occurs by diffusion rather than bulk flow, Proceedings of the National Academy of Sciences, 114 (2017), pp. 9894–9899.
  * [26] Q. Hong, J. Kraus, M. Lymbery, and F. Philo, Conservative discretizations and parameter-robust preconditioners for Biot and multiple-network flux-based poroelasticity models, Numerical Linear Algebra with Applications, 26 (2019), p. e2242.
  * [27] J. J. Iliff, M. Wang, Y. Liao, B. A. Plogg, W. Peng, G. A. Gundersen, H. Benveniste, G. E. Vates, R. Deane, S. A. Goldman, E. A. Nagelhus, and M. Nedergaard, A paravascular pathway facilitates CSF flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid $\beta$, Science Translational Medicine, 4 (2012), p. e147ra111.
  * [28] M. Jayadharan, E. Khattatov, and I. Yotov, Domain decomposition and partitioning methods for mixed finite element discretizations of the Biot system of poroelasticity.  Submitted preprint. Available from https://arxiv.org/abs/2010.15353, 2021.
  * [29] S. Kumar, R. Oyarzúa, R. Ruiz-Baier, and R. Sandilya, Conservative discontinuous finite volume and mixed schemes for a new four-field formulation in poroelasticity, ESAIM: Mathematical Modelling and Numerical Analysis, 54 (2020), pp. 273–299.
  * [30] G. E. Lang, D. Vella, S. L. Waters, and A. Goriely, Mathematical modelling of blood-brain barrier failure and oedema, Mathematical Medicine and Biology: A Journal of the IMA, 34 (2016), pp. 391–414.
  * [31] J. J. Lee, K.-A. Mardal, and R. Winther, Parameter-robust discretization and preconditioning of Biot’s consolidation model, SIAM Journal on Scientific Computing, 39 (2017), pp. A1–A24.
  * [32] T. Li and I. Yotov, A mixed elasticity formulation for fluid-poroelastic structure interaction.  Submitted preprint. Available from https://arxiv.org/abs/2011.00132, 2020.
  * [33] A. Logg, K.-A. Mardal, and G. Wells, Automated Solution of Differential Equations by the Finite Element Method, Springer-Verlag, Berlin, Germany, 2012.
  * [34] X. Ma and M. D. Zoback, Laboratory experiments simulating poroelastic stress changes associated with depletion and injection in low-porosity sedimentary rocks, Journal of Geophysical Research: Solid Earth, 122 (2017), pp. 2478–2503.
  * [35] R. Oyarzúa and R. Ruiz-Baier, Locking-free finite element methods for poroelasticity, SIAM Journal on Numerical Analysis, 54 (2016), pp. 2951–2973.
  * [36] W. M. Pardridge, CSF, blood-brain barrier, and brain drug delivery, Expert Opinion on Drug Delivery, 13 (2016), pp. 963–975.
  * [37] E. Piersanti, J. J. Lee, T. Thompson, K.-A. Mardal, and M. E. Rognes, Parameter robust preconditioning by congruence for multiple-network poroelasticity, SIAM Journal on Scientific Computing, 43 (2021), pp. B984–B1007.
  * [38] A. Propp, A. Gizzi, F. Levrero-Florencio, and R. Ruiz-Baier, An orthotropic electro-viscoelastic model for the heart with stress-assisted diffusion, Biomechanics and Modeling in Mechanobiology, 19 (2020), pp. 633–659.
  * [39] R. E. Showalter, Poroelastic filtration coupled to Stokes flow, in Control Theory of Partial Differential Equations, O. Imanuvilov, G. Leugering, R. Triggiani, and B.-Y. Zhang, eds., vol. 242 of Lecture Notes in Pure and Applied Mathematics, Chapman & Hall, Boca Raton, 2005, pp. 229–241.
  * [40] E. Syková and C. Nicholson, Diffusion in brain extracellular space, Physiological Reviews, 88 (2008), pp. 1277–1340.
  * [41] M. Taffetani, R. Ruiz-Baier, and S. L. Waters, Coupling Stokes flow with inhomogeneous poroelasticity.  Submitted preprint. Available from https://arxiv.org/abs/2103.01529, 2021.
  * [42] L. M. Valnes, S. K. Mitusch, G. Ringstad, P. K. Eide, S. W. Funke, and K.-A. Mardal, Apparent diffusion coefficient estimates based on 24 hours tracer movement support glymphatic transport in human cerebral cortex, Scientific reports, 10 (2020), pp. 1–12.
  * [43] V. Vinje, A. Eklund, K.-A. Mardal, M. E. Rognes, and K.-H. Støverud, Intracranial pressure elevation alters csf clearance pathways, Fluids and Barriers of the CNS, 17 (2020), pp. 1–19.
  * [44] L. Xie, H. Kang, Q. Xu, M. J. Chen, Y. Liao, M. Thiyagarajan, J. O’Donnell, D. J. Christensen, C. Nicholson, J. J. Iliff, et al., Sleep drives metabolite clearance from the adult brain, science, 342 (2013), pp. 373–377.
  * [45] T. Zohdi and P. Wriggers, On the effects of microstress on macroscopic diffusion processes, Acta Mechanica, 136 (1999), pp. 91–107.