# The Biot-Stokes coupling using total pressure: formulation, analysis and
application to interfacial flow in the eye

Ricardo Ruiz-Baier<EMAIL_ADDRESS>School of Mathematics, Monash
University, 9 Rainforest Walk, Clayton, Victoria 3800, Australia.  Institute
of Computer Science and Mathematical Modelling, Sechenov University, Moscow,
Russian Federation.  Universidad Adventista de Chile, Casilla 7-D, Chillán,
Chile.  Matteo Taffetani<EMAIL_ADDRESS>Department of Engineering
Mathematics, University of Bristol, BS8 1TR Bristol, United Kingdom.  Hans D.
Westermeyer<EMAIL_ADDRESS>Department of Clinical Sciences, College of
Veterinary Medicine, North Carolina State University, 1060 William Moore
Drive, Raleigh, NC 27607, USA.  Ivan Yotov<EMAIL_ADDRESS>Department of
Mathematics, University of Pittsburgh, Pittsburgh, PA 15260, USA.

###### Abstract

We consider a multiphysics model for the flow of Newtonian fluid coupled with
Biot consolidation equations through an interface, and incorporating total
pressure as an unknown in the poroelastic region. A new mixed-primal finite
element scheme is proposed solving for the pairs fluid velocity - pressure and
displacement - total poroelastic pressure using Stokes-stable elements, and
where the formulation does not require Lagrange multipliers to set up the
usual transmission conditions on the interface. The stability and well-
posedness of the continuous and semi-discrete problems are analysed in detail.
Our numerical study is framed in the context of applicative problems
pertaining to heterogeneous geophysical flows and to eye poromechanics. For
the latter, we investigate different interfacial flow regimes in Cartesian and
axisymmetric coordinates that could eventually help describe early morphologic
changes associated with glaucoma development in canine species.

###### keywords:

Porous media flow , Biot consolidation , total pressure , transmission problem
, mixed finite element methods , eye fluid poromechanics.

###### MSC:

65M60 , 65M12 , 76S05 , 74F10 , 92C35.

††journal: Computer Methods in Applied Mechanics and Engineering

##  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 the spatial growth of biofilm in
fluids. When the interaction with a free fluid is considered, the mechanics of
the fluid and poroelastic domains are coupled through balance of forces and
continuity conditions that adopt diverse forms depending on the expected
behaviour in the specific application (see, e.g., [26, 49, 50, 51] and the
references therein). The particular problem we consider in this paper as
motivation for the design of the finite element formulation is the interfacial
flow of aqueous humour between the anterior chamber and the trabecular
meshwork (which is a deformable porous structure) in the eye, and how such
phenomenon relates to early stages of glaucoma.

Glaucoma encompasses a group of mechanisms that lead to decreased retinal
function, impaired visual fields and blindness. The main risk factor for
glaucoma in canines is an abnormal increase in the intra-ocular pressure
(which under physiologically normal conditions is balanced between aqueous
humour production and outflow to the venous drainage system [32]). We are
interested in modelling the flow behaviour of aqueous humour within the
anterior chamber and its interaction with the poroelastic properties of
particular compartments in the drainage outlet located between the base of the
iris and the limbus, which, in the dog eye and most other non-primate species
consists of an array of thin tissue columns (pectinate ligaments) [46] which
mark the boundary of the trabecular meshwork with the anterior chamber.
Sketches of the regions of interest are depicted in Figure 1.1. Our focus is
on how the physical changes associated with pectinate ligament dysplasia, a
change seen in all dogs with primary angle closure glaucoma, affect aqueous
humour flow through this boundary. We stress that the ciliary cleft anatomy of
all carnivorous mammals is fairly similar to that of the dog. In fact, this
anatomy is preserved across much of herbivorous mammals as well (see, e.g.,
[43]). Therefore, even though the dog is probably the most studied due to its
status as a companion animal, this work likely applies to most carnivorous and
herbivorous mammals.

The flow within the anterior chamber will be modelled by Navier-Stokes and
Stokes’ law for Newtonian fluids, whereas the filtration of aqueous humour
through the deformable trabecular meshwork and towards the angular aqueous
plexus will be described by Darcy’s law. Pressure differences are generated by
production (from the ciliary muscle) and drainage (to angular aqueous plexus
and then linked to the veins at the surface of the sclera through collecting
channels) of aqueous humour.

Other effects that could contribute to modification of the flow patterns and
that we do not consider here, are thermal properties (buoyancy mechanisms due
to temperature gradients from inner to outer cornea) [11], cross-link
interaction between fibrils in the cornea [30], pressure changes due to
phacodonesis (vibration of the lens while the head or eye itself moves) and
Rapid Eye Movement during sleep [29], and nonlinear flow conditions in the
filtration region (incorporated in [37] through Darcy-Forchheimer models).

In contrast with [25, 28, 42, 52], here we consider that the coalescing of the
pectinate ligaments results in marked changes in porosity properties of the
anterior chamber - trabecular meshwork interface, which could eventually lead
to progressive collapse of the ciliary cleft. We further postulate that these
modifications of the tissue’s microstructure could be induced by forces
exerted by the flow that concentrate at the interface between the dysplastic
pectinate ligament and the anterior chamber, and which occur over a timescale
much larger than that of the ocular pulsating flow. In fact, evidence of the
compliance of the trabecular meshwork can be found in, e.g., [35]. One of the
earliest modelling works including a coupling between aqueous humour in the
anterior chamber with complying structures is presented in [33], where
mechanical properties of the bovine iris were employed to set an elastic
interface to represent blinking. Other fluid-structure interaction models have
been recently developed in [54], suggesting that flow conditions in the
trabecular meshwork and the outlets could be largely affected by the changes
of permeability in microstructure, and [3], where poroelastic properties of
the choroid and viscoelastic response of the vitreous body are used to set up
a more complete 3D model of larger scale that discards a dedicated
physiological description of the trabecular meshwork and considers instead a
windkessel model.

Figure 1.1: Histology sketches of regions of interest (left, centre), and
ultrasound image taken from data (right).

In the general context of single phase fluid / poromechanical coupling, there
are already a variety of finite element formulations starting from the work
[44], which focuses on the effects of secondary consolidation. More recently,
partitioned finite element formulations using domain decomposition and or
Nitsche’s approach for single and double poroelastic layers in contact with a
single phase fluid can be found in [9, 16, 17, 39]. Monolithic couplings have
been analysed in [7] for a mixed Darcy formulation using a Lagrange multiplier
to impose flux continuity (see also [5] for the extension to the case of non-
Newtonian fluids), as well in [19, 20] for a primal Darcy formulation. Ghost
penalty methods have been employed for cut FEM methods valid in the regime of
large deformations in [2].

Here, and drawing inspiration from the formulation in [7] and [27, 40, 45], we
rewrite the poroelasticity equations using three fields (displacement, fluid
pressure and total pressure). Compared to previous works [7, 5, 9, 16, 17, 19,
20, 39], which employ the classical displacement formulation, an advantage of
the present approach, inherited from [45], is that the formulation is free of
poroelastic locking, meaning that it is robust with respect to the Lamé
parameters of the poroelastic structure. This is of particular importance when
we test variations of the flow response to changes in the material properties
of the skeleton and when the solid approaches the incompressibility limit. The
present work also stands as an extension of the formulation recently employed
in [50] (where only the case of intrinsic incompressible constituent in the
poroelastic region were considered) to obtain approximate solutions for
heterogeneous poroelasticity coupled with Stokes flow in channels (and using
also heterogeneous elastic moduli); while the PDE analysis, numerical aspects,
and applicability of the formalism to more realistic scenarios have not yet
been addressed. In this work, under adequate assumptions, the analysis of the
weak formulation is carried out, using a (time continuous) semi-discrete
Galerkin approximation and a weak compactness argument. The well-posedness of
the semi-discrete formulation is established using the theory of differential
algebraic equations (see, e.g., [15]) and using similar results to those
obtained in [7, 53]. A conforming mixed finite element scheme of general order
is used. Furthermore, a fully discrete scheme based on backward Euler’s time
discretisation is considered, and the unique solvability and convergence for
the fully discrete scheme are established.

We have organised the contents of this paper in the following manner. Section
2 outlines the model problem, motivating each term in the balance equations
and stating the interfacial and boundary conditions. Section 3 states the weak
form of the governing equations in Cartesian and axisymmetric coordinates.
Then, in Section 4, we address the construction of the finite element scheme,
the well-posedness of the continuous and discrete problems, the stability of
the fully discrete system in matrix form. Section 5 states the fully-discrete
scheme and presents the error estimates. In Section 6 we collect computational
results consisting in verification of spatio-temporal convergence and analysis
of different cases on simplified and more physiologically accurate geometries,
including also a typical application in reservoir modelling. One of the
examples involves large displacements near the interface, in which case a
harmonic extension operator is used to deform the fluid domain. We close with
a summary, some remarks and a discussion on model generalisations in Section
7.

##  2 Governing equations

Let us consider a spatial domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$
disjointly split into $\Omega_{F}$ and $\Omega_{P}$ representing,
respectively, the regions where a chamber filled with incompressible fluid and
the deformable porous structure are located. We will denote by
$\boldsymbol{n}$ the unit normal vector on the boundary $\partial\Omega$, and
by $\Sigma=\Omega_{F}\cap\Omega_{P}$ the interface between the two subdomains.
We also define the boundaries $\Gamma_{F}=\partial\Omega_{F}\setminus\Sigma$
and $\Gamma_{P}=\partial\Omega_{P}\setminus\Sigma$, and adopt the convention
that on $\Sigma$ the normal vector points from $\Omega_{F}$ to $\Omega_{P}$.
See a rough sketch in Figure 2.1, that represents the geometry of the anterior
segment in the eye distinguishing between the anterior chamber $\Omega_{F}$
and the trabecular meshwork $\Omega_{P}$. The domain is sketched as an
axisymmetric region, for which more specific properties will be listed later
on.

In presenting the set of governing equations for the coupled fluid -
poroelastic system we first focus on the fluid domain, then on the poroelastic
domain and, finally, on the initial, boundary and interfacial conditions.

Figure 2.1: Schematic diagram of multidomain configuration on a segmented and
meshed geometry, including the location of boundaries and interface. The inlet
region $\Gamma^{\mathrm{in}}$ and the symmetry axis
$\Gamma^{\mathrm{axisymm}}$ are part of the boundary
$\Gamma^{\boldsymbol{u}}_{F}$, whereas the outlet region
$\Gamma^{\mathrm{out}}$ is part of the boundary $\Gamma^{\boldsymbol{d}}_{P}$.

###  2.1 Fluid domain

In the fluid domain $\Omega_{F}$, the problem is governed by the momentum and
mass conservation equations. Defining the fluid velocity $\boldsymbol{u}$ and
the fluid pressure $p_{F}$, the resulting system is written as

$\displaystyle\rho_{f}(\partial_{t}\boldsymbol{u}+\boldsymbol{u}\cdot\nabla\boldsymbol{u})-\mathop{\mathbf{div}}\nolimits[2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I}]$
$\displaystyle=\rho_{f}\boldsymbol{g}$ $\displaystyle\text{in
$\Omega_{F}\times(0,T]$},$ (2.1a)
$\displaystyle\mathop{\mathrm{div}}\nolimits\boldsymbol{u}$ $\displaystyle=0$
$\displaystyle\text{in $\Omega_{F}\times(0,T]$},$ (2.1b)

where $\rho_{f},\mu_{f}$ are the density and dynamic viscosity of the fluid
(e.g. the aqueous humour if we think to the application for the filtration in
the eye and thus $\Omega_{F}$ is the anterior chamber), $\boldsymbol{g}$ is
the gravity acceleration,
$\boldsymbol{\epsilon}(\boldsymbol{u})=\frac{1}{2}(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\tt
t})$ is the strain rate tensor and $\partial_{t}$ indicates derivatives with
respect to time. Note that, in typical flow conditions of the eye anterior
chamber, the Reynolds number is low (approximately 1.25, while the reduced
Reynolds number is between 0.06-0.08. See, e.g., [18, 29, 35, 42]). Therefore,
for the stability analysis and for some of the numerical tests later on, we
will restrict the fluid model to Stokes’ equations.

###  2.2 Poroelastic domain

The poroelastic domain $\Omega_{P}$ is a biphasic material constituted by a
linear elastic solid phase (potentially intrinsically compressible) and an
intrinsically incompressible fluid phase. In the context of the eye
poromechanics, the trabecular meshwork region occupying $\Omega_{P}$ is
constituted by three distinctive tissues, the uveal meshwork, the
juxtacanalicular meshwork, and the corneoscleral network; they have different
micromechanical properties that, from our modelling perspective, can be
regarded as a single poroelastic domain with heterogeneous porosity
distribution that, in turn, means possible heterogeneity in the material
properties. In addition, we anticipate that, although the fluid viscosity is
relevant at the scale of the pore, we assume that the fluid can be treated as
inviscid at the macroscale. Calling $p_{P}$ the fluid pressure and
$\boldsymbol{d}$ the solid displacement, here we introduce the formulation
that assumes pressure and displacement as primary variables as presented, for
example, in [24]

$\displaystyle\partial_{t}\left(C_{0}p_{P}+\alpha\mathop{\mathrm{div}}\nolimits\boldsymbol{d}\right)-\mathop{\mathrm{div}}\nolimits\biggl{(}{\frac{\kappa}{\mu_{f}}(\nabla
p_{P}-\rho_{f}\boldsymbol{g})}\biggr{)}$ $\displaystyle=0$
$\displaystyle\text{in $\Omega_{P}\times(0,T]$},$ (2.2a)
$\displaystyle-\mathop{\mathbf{div}}\nolimits[2\mu_{s}\boldsymbol{\epsilon}(\boldsymbol{d})+\lambda\left(\mathop{\mathrm{div}}\nolimits\boldsymbol{d}\right)\mathbf{I}-\alpha
p_{P}\mathbf{I}]$ $\displaystyle=\rho_{m}\boldsymbol{f}$
$\displaystyle\text{in $\Omega_{P}\times(0,T]$}.$ (2.2b)

The first equation can be derived from the conservation of mass for the fluid
phase once employing the Darcy’s law and the relation between fluid content -
pressure - hydrostatic deformation of the solid phase. The storage capacity
$C_{0}$ is related to the intrinsic compressibility of the solid phase, while
$\kappa$ is the permeability (assumed isotropic but heterogeneous). The second
equation is the conservation of the momentum for the mixture, where
$\boldsymbol{f}$ is a (possibly fluid pressure-dependent) body load,
$\boldsymbol{\epsilon}(\boldsymbol{d})=\frac{1}{2}(\nabla\boldsymbol{d}+\nabla\boldsymbol{d}^{\tt
t})$ is the infinitesimal strain tensor, $\rho_{s}$ is the density of the
porous matrix and $\rho_{m}$ is the average density of the poroelastic body,
$\lambda,\mu_{s}$ are the Lamé constants of the solid; the term in the
divergence in the left-hand side of (2.2b) is known as the effective stress or
Terzaghi stress and the parameter $\alpha$, also known as Biot-Willis
poroelastic coefficient, depends on the intrinsic compressibility of the solid
phase ($\alpha=1$ when the solid phase is intrinsically incompressible). We
refer an interested reader to [24] or [23] for further details.

For sake of robustness of the formulation with respect to $\lambda$, we
introduce the _total pressure_ $\varphi:=\alpha
p_{P}-\lambda\mathop{\mathrm{div}}\nolimits\boldsymbol{d}$, as an additional
unknown in the system (following [40, 45]), and rewrite the pressure-
displacement formulation in (2.2) in term of the solid displacement
$\boldsymbol{d}$, the fluid pressure $p_{P}$, and the total pressure
$\varphi$, as

$\displaystyle-\mathop{\mathbf{div}}\nolimits[2\mu_{s}\boldsymbol{\epsilon}(\boldsymbol{d})-\varphi\mathbf{I}]$
$\displaystyle=\rho_{m}\boldsymbol{f}$ $\displaystyle\text{in
$\Omega_{P}\times(0,T]$},$ (2.3a) $\displaystyle\varphi-\alpha
p_{P}+\lambda\mathop{\mathrm{div}}\nolimits\boldsymbol{d}$ $\displaystyle=0$
$\displaystyle\text{in $\Omega_{P}\times(0,T]$},$ (2.3b)
$\displaystyle\bigl{(}C_{0}+\frac{\alpha^{2}}{\lambda}\bigr{)}\partial_{t}p_{P}-\frac{\alpha}{\lambda}\partial_{t}\varphi-\mathop{\mathrm{div}}\nolimits\biggl{(}{\frac{\kappa}{\mu_{f}}(\nabla
p_{P}-\rho_{f}\boldsymbol{g})}\biggr{)}$ $\displaystyle=0$
$\displaystyle\text{in $\Omega_{P}\times(0,T]$}.$ (2.3c)

In contrast with the formulations in [16, 7], here we do not employ the fluid
velocity in the porous domain as a separate unknown. Consistently with [41],
in the framework of linear poroelasticity, the permeability depends on the
pressure and the displacement fields only at higher orders; however, it can be
heterogeneous, with a spatial distribution dependent, for example, on an
initial porosity distribution. We thus simply write

$\kappa=\kappa(\boldsymbol{x}).$ (2.4)

Likewise, it is also possible to assume heterogeneity of the Lamé constants,
as in [50]. In such a case, we need to assume that there exist constants
$\lambda_{\min}$ and $\lambda_{\max}$ such that
$0<\lambda_{\min}\leq\lambda\leq\lambda_{\max}$. In the analysis we show that
the estimates are independent of $\lambda_{\max}$, i.e., the results are
uniform in the almost incompressible limit, even for the heterogeneous case.
Heterogeneous permeability and/or heterogeneous Lamé parameters will be used
in some of the numerical examples of Section 6.

###  2.3 Initial, boundary and transmission conditions

To close the system composed by (2.1a), (2.1b), (2.2b), (2.2a), we need to
provide suitable initial data, boundary conditions, and adequate transmission
assumptions. Without losing generality, we suppose that

$\boldsymbol{u}(0)=\boldsymbol{u}_{0},\quad p_{P}(0)=p_{P,0}\quad\text{in
$\Omega\times\\{0\\}$,}$ (2.5)

while for the rest of the variables we will construct compatible initial data.
In particular, and for sake of the energy estimates to be addressed in Section
4.2, we will require initial displacement and an initial total pressure (which
in turn is computed from the initial displacement and the initial fluid
pressure).

On the boundary $\Gamma_{F}$ of the fluid domain we can apply conditions on
either the velocity, or the stress tensor; we thus decomposed it between
$\Gamma_{F}^{\boldsymbol{u}}$ and $\Gamma_{F}^{\boldsymbol{\sigma}}$, with
$|\Gamma_{F}^{\boldsymbol{u}}|>0$, where we impose, respectively, no slip
velocities and zero normal total stresses as

$\displaystyle\boldsymbol{u}$ $\displaystyle=\boldsymbol{0}$
$\displaystyle\text{on $\Gamma_{F}^{\boldsymbol{u}}\times(0,T]$},$ (2.6a)
$\displaystyle[2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I}]\boldsymbol{n}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\text{on
$\Gamma_{F}^{\boldsymbol{\sigma}}\times(0,T]$}.$ (2.6b)

Similarly, on the boundary $\Gamma_{P}$ of the poroelastic domain we can
prescribe conditions on either the displacement or the traction and either the
pressure or the fluid flux $\mathbf{q}=-\frac{\kappa}{\mu_{f}}(\nabla
p_{P}-\rho_{f}\boldsymbol{g})$. We thus divide the boundary into
$\Gamma_{P}^{p_{P}}$ and $\Gamma_{P}^{\boldsymbol{d}}$, with
$|\Gamma_{P}^{p_{P}}|>0$ and $|\Gamma_{P}^{\boldsymbol{d}}|>0$, where we
apply, respectively,

$\displaystyle\boldsymbol{d}=\boldsymbol{0}\quad\text{and}\quad\frac{\kappa}{\mu_{f}}(\nabla
p_{P}-\rho_{f}\boldsymbol{g})\cdot\boldsymbol{n}$ $\displaystyle=0$
$\displaystyle\text{on $\Gamma_{P}^{\boldsymbol{d}}\times(0,T]$},$ (2.7a)
$\displaystyle[2\mu_{s}\boldsymbol{\epsilon}(\boldsymbol{d})-\varphi\mathbf{I}]\boldsymbol{n}=\boldsymbol{0}\quad\text{and}\quad
p_{P}$ $\displaystyle={0}$ $\displaystyle\text{on
$\Gamma_{P}^{p_{P}}\times(0,T]$}.$ (2.7b)

Next, and following [9, 16, 44, 49], we consider transmission conditions on
$\Sigma$ accounting for the continuity of normal fluxes, momentum
conservation, balance of fluid normal stresses, and the so-called Beavers-
Joseph-Saffman (BJS) condition for tangential fluid forces

$\displaystyle\boldsymbol{u}\cdot\boldsymbol{n}$
$\displaystyle=\left(\partial_{t}\boldsymbol{d}-\frac{\kappa}{\mu_{f}}(\nabla
p_{P}-\rho_{f}\boldsymbol{g})\right)\cdot\boldsymbol{n}$
$\displaystyle\text{on $\Sigma\times(0,T]$},$ (2.8a)
$\displaystyle(2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I})\boldsymbol{n}$
$\displaystyle=(2\mu_{s}\boldsymbol{\epsilon}(\boldsymbol{d})-\varphi\mathbf{I})\boldsymbol{n}$
$\displaystyle\text{on $\Sigma\times(0,T]$},$ (2.8b)
$\displaystyle-\boldsymbol{n}\cdot(2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I})\boldsymbol{n}$
$\displaystyle={p_{P}}$ $\displaystyle\text{on $\Sigma\times(0,T]$},$ (2.8c)
$\displaystyle-{\boldsymbol{t}^{j}}\cdot(2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I})\boldsymbol{n}$
$\displaystyle=\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\cdot{\boldsymbol{t}^{j},\quad
1\leq j\leq d-1}$ $\displaystyle\text{on $\Sigma\times(0,T]$},$ (2.8d)

where $\gamma>0$ is the slip rate coefficient (or tangential resistance
parameter), and we recall that the normal $\boldsymbol{n}$ on the interface is
understood as pointing from the fluid domain $\Omega_{F}$ towards the porous
structure $\Omega_{P}$, while $\boldsymbol{t}^{1}$ stands for the tangent
vector on $\Sigma$ (for the case of $d=2$, while for 3D
$\boldsymbol{t}^{1},\boldsymbol{t}^{2}$ represent the two tangent vectors on
the interface, normal to $\boldsymbol{n}$).

##  3 Weak formulation

Apart from the nomenclature introduced at the beginning of the section,
conventional notation will be adopted throughout the paper. For Lipschitz
domains $\Xi$ in $\mathbb{R}^{d}$, and for $s\in\mathbb{N}$,
$k\in\mathbb{N}\cup\infty$ we denote by $W^{k,s}(\Xi)$ the space of all
$L^{s}(\Xi)$ integrable functions with weak derivatives up to order $s$ being
also $L^{s}(\Xi)$ integrable. As usual, for the special case of $s=2$ we write
$H^{s}(\Xi):=W^{k,2}(\Xi)$ and use boldfaces to refer to vector-valued
functions and function spaces, e.g., $\mathbf{H}^{s}(\Xi):=[H^{s}(\Xi)]^{d}$.
We will further utilise the Bochner space-time norms, for a separable Banach
space $\mathrm{V}$ and $f:(0,T)\to\mathrm{V}$,
$\|f\|^{2}_{L^{2}(0,T;\mathrm{V})}:=\int^{T}_{0}\|f(t)\|^{2}_{\mathrm{V}}\,\,\mathrm{d}t$
and
$\|f\|_{L^{\infty}(0,T;\mathrm{V})}:=\mathop{\mathrm{ess\,sup}}\limits_{t\in(0,T)}\|f(t)\|_{\mathrm{V}}$.
By $C$ we will denote generic constants that are independent of the mesh size.

###  3.1 Cartesian case

The initial step in deriving the finite element scheme consists in stating a
weak form for (2.1a)-(2.3c). We proceed to test these equations against
suitable smooth functions and to integrate over the corresponding subdomain.
After applying integration by parts wherever adequate, we formally end up with
the following remainder on the interface

$I_{\Sigma}=-\langle(2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I})\boldsymbol{n},\boldsymbol{v}\rangle_{\Sigma}+\langle(2\mu_{s}\boldsymbol{\epsilon}(\boldsymbol{d})-\varphi\mathbf{I})\boldsymbol{n},\boldsymbol{w}\rangle_{\Sigma}+\langle\frac{\kappa}{\mu_{f}}\nabla
p_{P}\cdot\boldsymbol{n},q_{P}\rangle_{\Sigma},$

where $\langle\cdot,\cdot\rangle_{\Sigma}$ denotes the duality pairing between
the trace functional space space $H^{1/2}(\Sigma)$ and its dual
$H^{-1/2}(\Sigma)$. Then, as in, e.g., [36], we proceed to use each of the
transmission conditions (2.8a)-(2.8d), yielding the expression

$I_{\Sigma}=\langle{p_{P}},(\boldsymbol{v}-\boldsymbol{w})\cdot\boldsymbol{n}\rangle_{\Sigma}+{\sum_{j=1}^{d-1}}\langle\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\cdot{\boldsymbol{t}^{j}},(\boldsymbol{v}-\boldsymbol{w})\cdot{\boldsymbol{t}^{j}}\rangle_{\Sigma}-\langle(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\cdot\boldsymbol{n},q_{P}\rangle_{\Sigma}.$

This interfacial term is well-defined because of the regularity of the
entities involved, and this implies that we do not require additional Lagrange
multipliers to realise the coupling conditions. Also, in view of the boundary
conditions we define the Hilbert spaces

$\displaystyle\mathbf{H}^{1}_{\star}(\Omega_{F})=\\{\boldsymbol{v}\in\mathbf{H}^{1}(\Omega_{F}):\boldsymbol{v}|_{\Gamma_{F}^{\boldsymbol{u}}}=\boldsymbol{0}\\},\quad\mathbf{H}^{1}_{\star}(\Omega_{P})=\\{\boldsymbol{w}\in\mathbf{H}^{1}(\Omega_{P}):\boldsymbol{w}|_{\Gamma_{P}^{\boldsymbol{d}}}=\boldsymbol{0}\\},$
$\displaystyle H^{1}_{\star}(\Omega_{P})=\\{q_{P}\in
H^{1}(\Omega_{P}):q_{P}|_{\Gamma_{P}^{p_{P}}}=0\\},$

associated with the classical norms in
$\mathbf{H}^{1}(\Omega_{F}),\mathbf{H}^{1}(\Omega_{P})$, and
$H^{1}(\Omega_{P})$, respectively. Consequently we have the following mixed
weak form: For $t\in[0,T]$, find
$\boldsymbol{u}\in\mathbf{H}^{1}_{\star}(\Omega_{F})$, $p_{F}\in
L^{2}(\Omega_{F})$, $\boldsymbol{d}\in\mathbf{H}^{1}_{\star}(\Omega_{P})$,
$p_{P}\in H^{1}_{\star}(\Omega_{P})$, $\varphi\in L^{2}(\Omega_{P})$, such
that

$\displaystyle
a_{1}^{F}(\partial_{t}\boldsymbol{u},\boldsymbol{v})+a_{2}^{F}(\boldsymbol{u},\boldsymbol{v})+c^{F}(\boldsymbol{u},\boldsymbol{u};\boldsymbol{v})+b^{F}_{1}(\boldsymbol{v},p_{F})$
$\displaystyle\qquad\quad+{b^{\Sigma}_{2}(\boldsymbol{v},p_{P})}+b_{3}^{\Sigma}(\boldsymbol{v},\partial_{t}\boldsymbol{d})=F^{F}(\boldsymbol{v})$
$\displaystyle\forall\boldsymbol{v}\in\mathbf{H}^{1}_{\star}(\Omega_{F}),$
(3.1a) $\displaystyle-b^{F}_{1}(\boldsymbol{u},q_{F})=0$ $\displaystyle\forall
q_{F}\in L^{2}(\Omega_{F}),$ (3.1b) $\displaystyle
b_{3}^{\Sigma}(\boldsymbol{u},\boldsymbol{w})+{b_{4}^{\Sigma}(\boldsymbol{w},p_{P})}+a_{1}^{P}(\boldsymbol{d},\boldsymbol{w})$
$\displaystyle\qquad\quad+a_{2}^{\Sigma}(\partial_{t}\boldsymbol{d},\boldsymbol{w})+b_{1}^{P}(\boldsymbol{w},\varphi)=F^{P}(\boldsymbol{w})$
$\displaystyle\forall\boldsymbol{w}\in\mathbf{H}^{1}_{\star}(\Omega_{P}),$
(3.1c) $\displaystyle-
b_{2}^{\Sigma}(\boldsymbol{u},q_{P})-b_{4}^{\Sigma}(\partial_{t}\boldsymbol{d},q_{P})+a_{3}^{P}(\partial_{t}p_{P},q_{P})$
$\displaystyle\qquad\quad+a_{4}^{P}(p_{P},q_{P})-b_{2}^{P}(\partial_{t}\varphi,q_{P})=G(q_{P})$
$\displaystyle\forall q_{P}\in H^{1}_{\star}(\Omega_{P}),$ (3.1d)
$\displaystyle-
b_{1}^{P}(\boldsymbol{d},\psi)-b_{2}^{P}(\psi,p_{P})+a_{5}^{P}(\varphi,\psi)=0$
$\displaystyle\forall\psi\in L^{2}(\Omega_{P}),$ (3.1e)

where the bilinear and trilinear forms and linear functionals are defined as

$\displaystyle
a_{1}^{F}\\!(\boldsymbol{u},\boldsymbol{v})=\rho_{f}\\!\int_{\Omega_{F}}\\!\\!\\!\boldsymbol{u}\cdot\boldsymbol{v},\quad
a_{2}^{F}\\!(\boldsymbol{u},\boldsymbol{v})=2\mu_{f}\\!\int_{\Omega_{F}}\\!\\!\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v})+{\sum_{j=1}^{d-1}}\langle\frac{\gamma\mu_{f}}{\sqrt{\kappa}}\boldsymbol{u}\cdot{\boldsymbol{t}^{j}},\boldsymbol{v}\cdot{\boldsymbol{t}^{j}}\rangle_{\Sigma},$
$\displaystyle
c^{F}\\!(\boldsymbol{u},\boldsymbol{w};\boldsymbol{v})=\rho_{f}\\!\\!\int_{\Omega_{F}}\\!\\!\\!(\boldsymbol{u}\cdot\nabla\boldsymbol{w})\cdot\boldsymbol{v},\quad
b^{F}_{1}(\boldsymbol{v},q_{F})=-\int_{\Omega_{F}}q_{F}\mathop{\mathrm{div}}\nolimits\boldsymbol{v},\quad
b_{1}^{P}(\boldsymbol{w},\psi)=-\int_{\Omega_{P}}\psi\mathop{\mathrm{div}}\nolimits\boldsymbol{w},$
$\displaystyle b_{2}^{\Sigma}(\boldsymbol{v},q_{P})=\langle
q_{P},\boldsymbol{v}\cdot\boldsymbol{n}\rangle_{\Sigma},\quad
b_{3}^{\Sigma}(\boldsymbol{v},\boldsymbol{w})=-{\sum_{j=1}^{d-1}}\langle\frac{\gamma\mu_{f}}{\sqrt{\kappa}}\boldsymbol{v}\cdot{\boldsymbol{t}^{j}},\boldsymbol{w}\cdot{\boldsymbol{t}^{j}}\rangle_{\Sigma},\quad
b_{4}^{\Sigma}(\boldsymbol{w},q_{P})=-\langle
q_{P},\boldsymbol{w}\cdot\boldsymbol{n}\rangle_{\Sigma},$ $\displaystyle
a_{1}^{P}(\boldsymbol{d},\boldsymbol{w})=2\mu_{s}\int_{\Omega_{P}}\boldsymbol{\epsilon}(\boldsymbol{d}):\boldsymbol{\epsilon}(\boldsymbol{w}),\quad
a_{2}^{\Sigma}(\boldsymbol{d},\boldsymbol{w})={\sum_{j=1}^{d-1}}\langle\frac{\gamma\mu_{f}}{\sqrt{\kappa}}\boldsymbol{d}\cdot{\boldsymbol{t}^{j}},\boldsymbol{w}\cdot{\boldsymbol{t}^{j}}\rangle_{\Sigma},$
(3.2) $\displaystyle
a_{3}^{P}(p_{P},q_{P})=\bigl{(}C_{0}+\frac{\alpha^{2}}{\lambda}\bigr{)}\int_{\Omega_{P}}p_{P}q_{P},\quad
a_{4}^{P}(p_{P},q_{P})=\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\nabla
p_{P}\cdot\nabla q_{P},$ $\displaystyle
b_{2}^{P}(\psi,q_{P})=\frac{\alpha}{\lambda}\int_{\Omega_{P}}\psi q_{P},\quad
a_{5}^{P}(\varphi,\psi)=\frac{1}{\lambda}\int_{\Omega_{P}}\varphi\psi,\quad
F^{F}(\boldsymbol{v})=\rho_{f}\int_{\Omega_{F}}\boldsymbol{g}\cdot\boldsymbol{v},$
$\displaystyle
F^{P}(\boldsymbol{w})=\rho_{s}\int_{\Omega_{P}}\boldsymbol{f}\cdot\boldsymbol{w},\quad
G(q_{P})=\int_{\Omega_{P}}\rho_{f}\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\nabla
q_{P}-\langle\rho_{f}\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\boldsymbol{n},q_{P}\rangle_{\Sigma}.$

###  3.2 Axisymmetric case

For the specific application of interfacial flow in the eye, the radial
symmetry of the domain and of the flow conditions could be better represented
using axisymmetric formulations as in [8, 25, 42]. Then, the domain as well as
the expected flow properties are all symmetric with respect to the axis of
symmetry $\Gamma^{\mathrm{axisymm}}$. The model equations can be written in
the meridional domain $\Omega$ (making abuse of notation, and referring to
Figure 2.1). In such a setting the fluid velocity and solid displacement only
possess radial and vertical components and we recall that the divergence
operator of the generic vector field $\boldsymbol{v}$ in axisymmetric
coordinates (in radial and height variables $r,z$) is

$\mathop{\mathrm{div}_{\mathrm{a}}}\nolimits\boldsymbol{v}:=\partial_{z}v_{z}+\frac{1}{r}\partial_{r}(rv_{r}),$

while the notation of the gradient coincides with that in Cartesian
coordinates. The weak formulation (3.1) adopts the following modifications
(again making abuse of notation, the unknowns are denoted the same as in the
Cartesian case): Find $\boldsymbol{u}\in\widehat{\mathbf{V}}$,
$p_{F}\in\widehat{\mathrm{Q}}^{F}$, $\boldsymbol{d}\in\widehat{\mathbf{W}}$,
$p_{P}\in\widehat{\mathrm{Q}}^{P}$, $\varphi\in\widehat{\mathrm{Z}}$, such
that

$\displaystyle\rho_{f}\int_{\Omega_{F}}\\!\partial_{t}\boldsymbol{u}\cdot\boldsymbol{v}\,r\,\mathrm{d}r\,\mathrm{d}z+2\mu_{f}\int_{\Omega_{F}}\\!\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v})r\,\mathrm{d}r\,\mathrm{d}z+2\mu_{f}\int_{\Omega_{F}}\\!\frac{1}{r}u_{r}v_{r}\,\mathrm{d}r\,\mathrm{d}z-\int_{\Omega_{F}}p_{F}\mathop{\mathrm{div}_{\mathrm{a}}}\nolimits\boldsymbol{v}\,r\,\mathrm{d}r\,\mathrm{d}z$
$\displaystyle\quad+\int_{\Sigma}\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\\!\cdot\boldsymbol{t}\boldsymbol{v}\cdot\boldsymbol{t}r\,\mathrm{d}r\,\mathrm{d}z+\\!\\!{\int_{\Sigma}p_{F}\boldsymbol{v}\cdot\boldsymbol{n}r\,\mathrm{d}r\,\mathrm{d}z}=\rho_{f}\\!\\!\int_{\Omega_{F}}\\!\boldsymbol{g}\cdot\boldsymbol{v}r\,\mathrm{d}r\,\mathrm{d}z\quad\forall\boldsymbol{v}\in\widehat{\mathbf{V}},$
(3.3a)
$\displaystyle-\int_{\Omega_{F}}q_{F}\mathop{\mathrm{div}_{\mathrm{a}}}\nolimits\boldsymbol{u}\,r\,\mathrm{d}r\,\mathrm{d}z=0\qquad\forall
q_{F}\in\widehat{\mathrm{Q}}^{F},$ (3.3b) $\displaystyle
2\mu_{s}\int_{\Omega_{P}}\boldsymbol{\epsilon}(\boldsymbol{d}):\boldsymbol{\epsilon}(\boldsymbol{w})\,r\,\mathrm{d}r\,\mathrm{d}z+2\mu_{s}\int_{\Omega_{P}}\\!\frac{1}{r}d_{r}w_{r}\,\mathrm{d}r\,\mathrm{d}z-\int_{\Sigma}\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\cdot\boldsymbol{t}\boldsymbol{w}\cdot\boldsymbol{t}r\,\mathrm{d}r\,\mathrm{d}z$
$\displaystyle\qquad-\int_{\Omega_{P}}\varphi\mathop{\mathrm{div}_{\mathrm{a}}}\nolimits\boldsymbol{w}\,r\,\mathrm{d}r\,\mathrm{d}z-{\int_{\Sigma}p_{P}\boldsymbol{w}\cdot\boldsymbol{n}\,\mathrm{d}r\,\mathrm{d}z}=\rho_{s}\int_{\Omega_{P}}\boldsymbol{f}\cdot\boldsymbol{w}\,r\,\mathrm{d}r\,\mathrm{d}z\qquad\forall\boldsymbol{w}\in\widehat{\mathbf{W}},$
(3.3c)
$\displaystyle\bigl{(}C_{0}+\frac{\alpha^{2}}{\lambda}\bigr{)}\int_{\Omega_{P}}\partial_{t}p_{P}q_{P}\,r\,\mathrm{d}r\,\mathrm{d}z+\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\nabla
p_{P}\cdot\nabla
q_{P}\,r\,\mathrm{d}r\,\mathrm{d}z-\frac{\alpha}{\lambda}\int_{\Omega_{P}}\partial_{t}\varphi
q_{P}\,r\,\mathrm{d}r\,\mathrm{d}z$
$\displaystyle\quad+\int_{\Sigma}q_{P}(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\cdot\boldsymbol{n}r\,\mathrm{d}r\,\mathrm{d}z=\int_{\Omega_{P}}\\!\\!\rho_{f}\boldsymbol{g}\cdot\nabla
q_{P}r\,\mathrm{d}r\,\mathrm{d}z-\int_{\Sigma}\\!\rho_{f}\boldsymbol{g}\cdot\boldsymbol{n}q_{P}r\,\mathrm{d}r\,\mathrm{d}z\quad\forall
q_{P}\in\widehat{\mathrm{Q}}^{P},$ (3.3d)
$\displaystyle-\int_{\Omega_{F}}\psi\mathop{\mathrm{div}_{\mathrm{a}}}\nolimits\boldsymbol{v}\,r\,\mathrm{d}r\,\mathrm{d}z+\frac{\alpha}{\lambda}\int_{\Omega_{P}}p_{P}\psi\,r\,\mathrm{d}r\,\mathrm{d}z-\frac{1}{\lambda}\int_{\Omega_{P}}\varphi\psi\,r\,\mathrm{d}r\,\mathrm{d}z=0\qquad\forall\psi\in\widehat{\mathrm{Z}}.$
(3.3e)

Here the functional spaces are now defined as

$\displaystyle\widehat{\mathbf{V}}:=\\{\boldsymbol{v}\in
V_{1}^{1}(\Omega_{F})\times\mathrm{H}_{1}^{1}(\Omega_{F}):\boldsymbol{v}|_{\Gamma_{F}^{\boldsymbol{u}}}=\boldsymbol{0}\
\text{and}\
\boldsymbol{v}\cdot\boldsymbol{n}|_{\Gamma^{\text{axisymm}}}=0\\},\quad\widehat{\mathrm{Q}}^{F}:=L^{2}_{1}(\Omega_{F}),$
$\displaystyle\widehat{\mathbf{W}}:=\\{\boldsymbol{w}\in
V_{1}^{1}(\Omega_{P})\times\mathrm{H}_{1}^{1}(\Omega_{P}):\boldsymbol{w}|_{\Gamma_{F}^{\boldsymbol{d}}}=\boldsymbol{0}\\},\
\widehat{\mathrm{Q}}^{P}:=\\{q_{P}\in\mathrm{H}_{1}^{1}(\Omega_{P}):q_{P}|_{\Gamma_{P}^{p_{P}}}=0\\},\
\widehat{\mathrm{Z}}:=L^{2}_{1}(\Omega_{P}),$

where, for $m\in\mathbb{R}$, $1\leq p<\infty$ the weighted functional spaces
adopt the specification

$L_{m}^{p}(\Omega_{i})=\\{v:\|v\|^{p}_{m,p,\Omega_{i}}:=\int_{\Omega_{i}}|v|^{p}r^{m}\,\mathrm{d}r\,\mathrm{d}z<\infty\\},$

and

$\mathrm{H}_{1}^{1}(\Omega_{i}):=\\{v\in L_{1}^{2}(\Omega_{i}):\nabla
v\in\mathbf{L}_{1}^{2}(\Omega_{i})\\},\quad
V_{1}^{1}(\Omega_{i}):=\mathrm{H}_{1}^{1}(\Omega_{i})\cap
L_{-1}^{2}(\Omega_{i}).$

##  4 Well-posedness of the weak formulation

The following analysis is confined to the Cartesian case. Furthermore, we
focus on the quasi-static Biot-Stokes model, i.e., we neglect the terms
$a_{1}^{F}$ and $c^{F}$ in (3.1a), as this is the typical flow regime for the
application of interest. We also restrict our attention to the case
$\tilde{\alpha}=1$. The solvability analysis is based on a Galerkin argument,
where one considers the semi-discrete continuous in time formulation with a
discretisation parameter $h$. We establish that it has a unique solution and
derive stability bounds. Then, owing to a weak compactness argument, we pass
to the limit $h\to 0$ and obtain existence and uniqueness of a weak solution.

###  4.1 Semi-discrete mixed finite element formulation

In addition to the assumptions stated before on the domain geometry, to avoid
additional technicalities, we operate under the condition that $\Omega$ is a
polytope. We denote by $\\{\mathcal{T}_{h}\\}_{h>0}$ a shape-regular family of
finite element partitions of $\bar{\Omega}$, conformed by tetrahedra (or
triangles in 2D) $K$ of diameter $h_{K}$, with mesh size
$h:=\max\\{h_{K}:K\in\mathcal{T}_{h}\\}$. The finite-dimensional subspaces for
fluid velocity, fluid pressure, porous displacement, porous fluid pressure,
and total pressure are denoted
$\mathbf{V}_{h},\mathrm{Q}^{F}_{h},\mathbf{W}_{h},\mathrm{Q}^{P}_{h},\mathrm{Z}_{h}$,
respectively. It is sufficient for the following analysis that the pairs
$(\mathbf{V}_{h},\mathrm{Q}^{F}_{h})$ and $(\mathbf{W}_{h},\mathrm{Z}_{h})$
are Stokes inf-sup stable, in the sense that there exist positive constants
$\beta_{F}$ and $\beta_{P}$ independent of $h$ such that

$\displaystyle\forall\,q_{F,h}\in\mathrm{Q}^{F}_{h},\quad\sup_{\boldsymbol{0}\neq\boldsymbol{v}_{h}\in\mathbf{V}_{h}}\frac{b_{1}^{F}(\boldsymbol{v}_{h},q_{F,h})}{\|\boldsymbol{v}_{h}\|_{\mathbf{H}^{1}(\Omega_{F})}}\geq\beta_{F}\|q_{F,h}\|_{L^{2}(\Omega_{F})},$
(4.1a)
$\displaystyle\forall\,\psi_{h}\in\mathrm{Z}_{h},\quad\sup_{\boldsymbol{0}\neq\boldsymbol{w}_{h}\in\mathbf{W}_{h}}\frac{b_{1}^{P}(\boldsymbol{w}_{h},\psi_{h})}{\|\boldsymbol{w}_{h}\|_{\mathbf{H}^{1}(\Omega_{P})}}\geq\beta_{P}\|\psi_{h}\|_{L^{2}(\Omega_{P})}.$
(4.1b)

Feasible choices are Taylor-Hood, the MINI element, Crouzeix-Raviart, Scott-
Vogelius, Guzmán-Neilan, Bernardi-Raugel, equal-order stabilised methods
(including total pressure projection stabilisation), divergence-conforming
stabilised methods, and many others. If taking, for example, generalised
Taylor-Hood elements of degree $(k,k-1)$ for
$\mathbf{V}_{h}\times\mathrm{Q}_{h}^{F}$ and
$\mathbf{W}_{h}\times\mathrm{Z}_{h}$, alongside a piecewise continuous and
polynomial space of degree $k$ for $\mathrm{Q}_{h}^{P}$, then their
interpolation properties will yield a method of overall order $k$ in space
when the displacement and velocity errors are measured in the
$\mathbf{H}^{1}$-norm, the total pressure and Stokes fluid pressure in the
$L^{2}$-norm, and the Biot fluid pressure in the $H^{1}$-norm. A more general
case, with different polynomial degrees for the different finite element
spaces, is addressed in Section 5.1, below.

We look for
$(\boldsymbol{u}_{h},p_{F,h},\boldsymbol{d}_{h},p_{P,h},\varphi_{h}):[0,T]\to\mathbf{V}_{h}\times\mathrm{Q}^{F}_{h}\times\mathbf{W}_{h}\times\mathrm{Q}^{P}_{h}\times\mathrm{Z}_{h}=:\mathbf{H}_{h}$
such that for a.e. $t\in(0,T]$,

$\displaystyle
a_{2}^{F}(\boldsymbol{u}_{h},\boldsymbol{v}_{h})+b^{F}_{1}(\boldsymbol{v}_{h},p_{F,h})+b^{\Sigma}_{2}(\boldsymbol{v}_{h},p_{P,h})+b_{3}^{\Sigma}(\boldsymbol{v}_{h},\partial_{t}\boldsymbol{d}_{h})$
$\displaystyle=F^{F}(\boldsymbol{v}_{h})$
$\displaystyle\\!\forall\boldsymbol{v}_{h}\in\mathbf{V}_{h},$ (4.2a)
$\displaystyle-b^{F}_{1}(\boldsymbol{u}_{h},q_{F,h})$ $\displaystyle=0$
$\displaystyle\\!\forall q_{F,h}\in\mathrm{Q}_{h}^{F},$ (4.2b) $\displaystyle
b_{3}^{\Sigma}(\boldsymbol{u}_{h},\boldsymbol{w}_{h})+b_{4}^{\Sigma}(\boldsymbol{w}_{h},p_{P,h})+a_{1}^{P}(\boldsymbol{d}_{h},\boldsymbol{w}_{h})$
$\displaystyle+a_{2}^{\Sigma}(\partial_{t}\boldsymbol{d}_{h},\boldsymbol{w}_{h})+b_{1}^{P}(\boldsymbol{w}_{h},\varphi_{h})$
$\displaystyle=F^{P}(\boldsymbol{w}_{h})$
$\displaystyle\\!\forall\boldsymbol{w}_{h}\in\mathbf{W}_{h},$ (4.2c)
$\displaystyle-
b_{2}^{\Sigma}(\boldsymbol{u}_{h},q_{P,h})-b_{4}^{\Sigma}(\partial_{t}\boldsymbol{d}_{h},q_{P,h})+a_{3}^{P}(\partial_{t}p_{P,h},q_{P,h})$
$\displaystyle+a_{4}^{P}(p_{P,h},q_{P,h})-b_{2}^{P}(\partial_{t}\varphi_{h},q_{P,h})$
$\displaystyle=G(q_{P,h})$ $\displaystyle\\!\forall
q_{P,h}\in\mathrm{Q}_{h}^{P},$ (4.2d) $\displaystyle-
b_{1}^{P}(\boldsymbol{d}_{h},\psi_{h})-b_{2}^{P}(\psi_{h},p_{P,h})+a_{5}^{P}(\varphi_{h},\psi_{h})$
$\displaystyle=0$ $\displaystyle\\!\forall\psi_{h}\in\mathrm{Z}_{h}.$ (4.2e)

This is a system of differential-algebraic equations (DAE) that can be written
in an operator form as

$\partial_{t}\mathcal{N}\underline{\boldsymbol{u}}_{h}(t)+\mathcal{M}\underline{\boldsymbol{u}}_{h}(t)=\mathcal{F}(t),$
(4.3)

with $\underline{\boldsymbol{u}}_{h}:=[\boldsymbol{u}_{h}\ p_{F,h}\
\boldsymbol{d}_{h}\ p_{P,h}\ \varphi_{h}]^{\tt t}$,

$\displaystyle\mathcal{N}=\left[\begin{array}[]{cc:ccc}\boldsymbol{0}&\boldsymbol{0}&(\mathcal{B}_{3}^{\Sigma})^{\prime}&\boldsymbol{0}&\boldsymbol{0}\\\\[4.30554pt]
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\hdashline\boldsymbol{0}&\boldsymbol{0}&\mathcal{A}_{2}^{\Sigma}&\boldsymbol{0}&\boldsymbol{0}\vphantom{\int^{X^{X}}}\\\\[4.30554pt]
\boldsymbol{0}&\boldsymbol{0}&-\mathcal{B}_{4}^{\Sigma}&\mathcal{A}_{3}^{P}&-(\mathcal{B}_{2}^{P})^{\prime}\\\\[4.30554pt]
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\end{array}\right]\\!,\quad\mathcal{M}=\left[\begin{array}[]{cc:ccc}\mathcal{A}_{2}^{F}&(\mathcal{B}_{1}^{F})^{\prime}&\boldsymbol{0}&(\mathcal{B}_{2}^{\Sigma})^{\prime}&\boldsymbol{0}\\\\[4.30554pt]
-\mathcal{B}_{1}^{F}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\hdashline\mathcal{B}_{3}^{\Sigma}&\boldsymbol{0}&\mathcal{A}_{1}^{P}&(\mathcal{B}_{4}^{\Sigma})^{\prime}&(\mathcal{B}_{1}^{P})^{\prime}\vphantom{\int^{X^{X}}}\\\\[4.30554pt]
-\mathcal{B}_{2}^{\Sigma}&\boldsymbol{0}&\boldsymbol{0}&\mathcal{A}_{4}^{P}&\boldsymbol{0}\\\\[4.30554pt]
\boldsymbol{0}&\boldsymbol{0}&-\mathcal{B}_{1}^{P}&-\mathcal{B}_{2}^{P}&\mathcal{A}_{5}^{P}\end{array}\right]\\!,\quad\mathcal{F}=\left[\begin{array}[]{c}\mathcal{F}^{F}\\\
\boldsymbol{0}\\\ \mathcal{F}^{P}\\\ \mathcal{G}\\\
\boldsymbol{0}\end{array}\right],$

where the operators in calligraphic letters are induced by the forms in (3.2).
The symbol $(\cdot)^{\prime}$ denotes the dual operator resulting in the
transpose block matrix of a given elementary block.

We next discuss the initial condition for (4.3). Recall from (2.5) that we are
given initial data $p_{P,0}$ and note that we do not use the initial data
$\boldsymbol{u}_{0}$, since we are considering the quasi-static Stokes model.
We assume that $p_{P,0}\in H^{1}_{\star}(\Omega_{p})$ and take
$p_{P,h}(0)=p_{P,h,0}$, where $p_{P,h,0}\in\mathrm{Q}_{h}^{P}$ is the
$L^{2}$-projection of $p_{P,0}$. In addition, we need initial data
$\boldsymbol{d}_{h,0}\in\mathbf{W}_{h}$ and $\varphi_{h,0}\in Z_{h}$ such that
$\boldsymbol{d}_{h}(0)=\boldsymbol{d}_{h,0}$ and
$\varphi_{h}(0)=\varphi_{h,0}$. The construction of such data will be
discussed in Theorem 4.1.

###  4.2 Well-posedness of the semi-discrete problem

The well-posedness of the semi-discrete formulation could be analysed by
recasting the system as a general parabolic operator with degenerate time
derivative, as recently proposed in [5] for the interaction of non-Newtonian
fluids and poroelastic media. This would require, however, to include the
solid velocity as a new unknown. Instead, we study it using the theory of DAE,
following the approach from [7].

In the forthcoming analysis we will appeal to the Poincaré inequality

$\|\nabla q_{P}\|_{L^{2}(\Omega_{P})}\geq
c_{P}\|q_{P}\|_{H^{1}(\Omega_{P})},\quad\forall\,q_{P}\in
H^{1}_{\star}(\Omega_{P}),$ (4.4)

Korn’s inequality

$\|\boldsymbol{\epsilon}(\boldsymbol{v})\|_{\mathbf{L}^{2}(\Omega_{F})}\geq
c_{K}^{F}\|\boldsymbol{v}\|_{\mathbf{H}^{1}(\Omega_{F})},\qquad\|\boldsymbol{\epsilon}(\boldsymbol{w})\|_{\mathbf{L}^{2}(\Omega_{P})}\geq
c_{K}^{P}\|\boldsymbol{w}\|_{\mathbf{H}^{1}(\Omega_{P})},$ (4.5)

for all $\boldsymbol{v}\in\mathbf{H}^{1}_{\star}(\Omega_{F})$,
$\boldsymbol{w}\in\mathbf{H}^{1}_{\star}(\Omega_{P})$, and the trace
inequality

$\|\boldsymbol{v}\|_{\mathbf{L}^{2}(\Sigma)}\leq
C_{\Sigma}^{F}\|\boldsymbol{v}\|_{\mathbf{H}^{1}(\Omega_{F})},\qquad\|q_{P}\|_{L^{2}(\Sigma)}\leq
C_{\Sigma}^{P}\|q_{P}\|_{H^{1}(\Omega_{P})},$ (4.6)

for all $\boldsymbol{v}\in\mathbf{H}^{1}(\Omega_{F})$, $q_{P}\in
H^{1}(\Omega_{P})$. We also assume that there exist constants
$0<\kappa_{1}<\kappa_{2}<\infty$ such that

$\kappa_{1}\leq\kappa(\boldsymbol{x})\leq\kappa_{2}\quad\text{for a.e.
}\boldsymbol{x}\in\Omega_{P}.$ (4.7)

###### Lemma 4.1

The bilinear forms that appear in (4.2) are continuous in the spaces
$\mathbf{V}$, $\mathrm{Q}^{F}$, $\mathbf{W}$, $\mathrm{Q}^{P}$, and
$\mathrm{Z}$. If $\boldsymbol{f}:[0,T]\to\mathbf{L}^{2}(\Omega_{P})$, then the
linear functionals on the right-hand sides are also continuous.

Proof. The statement of the lemma follows from the use of the Cauchy-Schwarz
inequality and the trace inequality (4.6). $\square$

###### Theorem 4.1

For each $\boldsymbol{f}\in H^{1}(0,T;\mathbf{L}^{2}(\Omega_{P}))$ and
$p_{P,0}\in H^{1}_{\star}(\Omega_{p})$, there exist initial data
$\boldsymbol{u}_{h,0}\in\mathbf{V}_{h}$, $p_{F,h,0}\in\mathrm{Q}^{F}_{h}$,
$\boldsymbol{d}_{h,0}\in\mathbf{W}_{h}$, and $\varphi_{h,0}\in Z_{h}$ such
that the semi-discrete problem (4.3) with initial conditions
$p_{P,h}(0)=p_{P,h,0}$, $\boldsymbol{d}_{h}(0)=\boldsymbol{d}_{h,0}$, and
$\varphi_{h}(0)=\varphi_{h,0}$ has a unique solution satisfying

$\displaystyle\|\boldsymbol{u}_{h}\|_{L^{2}(0,T;\mathbf{H}^{1}(\Omega_{F}))}+\|p_{F,h}\|_{L^{2}(0,T;L^{2}(\Omega_{F}))}+\|\boldsymbol{d}_{h}\|_{L^{\infty}(0,T;\mathbf{H}^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\sum_{j=1}^{d-1}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{L^{2}(0,T;L^{2}(\Sigma))}+\|p_{P,h}\|_{L^{\infty}(0,T;L^{2}(\Omega_{P}))}+\|p_{P,h}\|_{L^{2}(0,T;H^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\|\varphi_{h}\|_{L^{2}(0,T;L^{2}(\Omega_{P}))}+\frac{1}{\sqrt{\lambda}}\|\alpha\,p_{P,h}-\varphi_{h}\|_{L^{\infty}(0,T;L^{2}(\Omega_{P}))}$
$\displaystyle\leq
C\left(\|\boldsymbol{g}\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega))}+\|\boldsymbol{g}\cdot\boldsymbol{n}\|_{L^{2}(0,T;L^{2}(\Sigma))}+\|\boldsymbol{f}\|_{H^{1}(0,T;\mathbf{L}^{2}(\Omega_{P}))}+\|p_{P,0}\|_{H^{1}(\Omega_{P})}\right),$
(4.8a) and
$\displaystyle\|\boldsymbol{u}_{h}\|_{L^{\infty}(0,T;\mathbf{H}^{1}(\Omega_{F}))}+\|p_{F,h}\|_{L^{\infty}(0,T;L^{2}(\Omega_{F}))}+\|\partial_{t}\boldsymbol{d}_{h}\|_{L^{2}(0,T;\mathbf{H}^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\sum_{j=1}^{d-1}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{L^{\infty}(0,T;L^{2}(\Sigma))}+\|\partial_{t}p_{P,h}\|_{L^{2}(0,T;L^{2}(\Omega_{P}))}+\|p_{P,h}\|_{L^{\infty}(0,T;H^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\|\varphi_{h}\|_{L^{\infty}(0,T;L^{2}(\Omega_{P}))}+\frac{1}{\sqrt{\lambda}}\|\alpha\,\partial_{t}p_{P,h}-\partial_{t}\varphi_{h}\|_{L^{2}(0,T;L^{2}(\Omega_{P}))}$
$\displaystyle\leq
C\left(\|\boldsymbol{g}\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega))}+\|\boldsymbol{g}\cdot\boldsymbol{n}\|_{L^{2}(0,T;L^{2}(\Sigma))}+\|\boldsymbol{f}\|_{H^{1}(0,T;\mathbf{L}^{2}(\Omega_{P}))}+\|p_{P,0}\|_{H^{1}(\Omega_{P})}\right),$
(4.8b)

with a constant $C$ independent of $\lambda_{\max}$ and $h$.

Proof. To establish existence, we employ [15, Th. 2.3.1], which asserts that
(4.3) has a solution if the matrix pencil $s\mathcal{N}+\mathcal{M}$ is
nonsingular for some $s\neq 0$. The solvability of the associated initial
value problem requires initial data that is consistent with the DAE system. To
deal with this issue, we first consider a related DAE system by including new
variables $\theta_{h}^{j}\in\mathbf{W}_{h}\cdot\boldsymbol{t}^{j}$,
$j=1,\ldots d-1$ and equations

$\langle\theta_{h}^{j},\boldsymbol{w}_{h}\cdot\boldsymbol{t}^{j}\rangle_{\Sigma}=\langle\partial_{t}\boldsymbol{d}_{h}\cdot\boldsymbol{t}^{j},\boldsymbol{w}_{h}\cdot\boldsymbol{t}^{j}\rangle_{\Sigma}\quad\forall\boldsymbol{w}_{h}\in\mathbf{W}_{h},\quad
j=1,\ldots d-1,$ (4.9)

and replacing $\partial_{t}\boldsymbol{d}_{h}\cdot\boldsymbol{t}^{j}$ with
$\theta_{h}^{j}$ in (4.2a) and (4.2c). Let
$\tilde{}\underline{\boldsymbol{u}}_{h}:=[\boldsymbol{u}_{h}\ p_{F,h}\
\boldsymbol{d}_{h}\ p_{P,h}\ \varphi_{h}\
\theta_{h}^{1}\cdots\,\theta_{h}^{d-1}]^{\tt t}$ and denote the extended DAE
system by

$\partial_{t}\tilde{\mathcal{N}}\tilde{}\underline{\boldsymbol{u}}_{h}(t)+\tilde{\mathcal{M}}\tilde{}\underline{\boldsymbol{u}}_{h}(t)=\tilde{\mathcal{F}}(t).$
(4.10)

Clearly any solution of (4.10) also solves (4.3). We will apply [15, Th.
2.3.1] to (4.10). We will show that the matrix
$\tilde{\mathcal{N}}+\tilde{\mathcal{M}}$ is nonsingular by proving that the
system
$(\tilde{\mathcal{N}}+\tilde{\mathcal{M}})\tilde{}\underline{\boldsymbol{u}}_{h}=\boldsymbol{0}$
has only the zero solution. By eliminating $\theta_{h}^{j}$, this system
results in
$(\mathcal{N}+\mathcal{M})\underline{\boldsymbol{u}}_{h}=\boldsymbol{0}$.
Using that

$\mathcal{N}+\mathcal{M}=\left[\begin{array}[]{cc:ccc}\mathcal{A}_{2}^{F}&(\mathcal{B}_{1}^{F})^{\prime}&(\mathcal{B}_{3}^{\Sigma})^{\prime}&(\mathcal{B}_{2}^{\Sigma})^{\prime}&\boldsymbol{0}\\\\[4.30554pt]
-\mathcal{B}_{1}^{F}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\hdashline\mathcal{B}_{3}^{\Sigma}&\boldsymbol{0}&\mathcal{A}_{1}^{P}+\mathcal{A}_{2}^{\Sigma}&(\mathcal{B}_{4}^{\Sigma})^{\prime}&(\mathcal{B}_{1}^{P})^{\prime}\vphantom{\int^{X^{X}}}\\\\[4.30554pt]
-\mathcal{B}_{2}^{\Sigma}&\boldsymbol{0}&-\mathcal{B}_{4}^{\Sigma}&\mathcal{A}_{3}^{P}+\mathcal{A}_{4}^{P}&-(\mathcal{B}_{2}^{P})^{\prime}\\\\[4.30554pt]
\boldsymbol{0}&\boldsymbol{0}&-\mathcal{B}_{1}^{P}&-\mathcal{B}_{2}^{P}&\mathcal{A}_{5}^{P}\end{array}\right],$

the equation $\underline{\boldsymbol{u}}_{h}^{\tt
t}(\mathcal{N}+\mathcal{M})\underline{\boldsymbol{u}}_{h}=0$ gives

$\displaystyle
2\mu_{f}\\!\int_{\Omega_{F}}\\!\\!\boldsymbol{\epsilon}(\boldsymbol{u}_{h}):\boldsymbol{\epsilon}(\boldsymbol{u}_{h})+\sum_{j=1}^{d-1}\langle\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(\boldsymbol{u}_{h}-\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\,,(\boldsymbol{u}_{h}-\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\rangle_{\Sigma}+2\mu_{s}\int_{\Omega_{P}}\boldsymbol{\epsilon}(\boldsymbol{d}_{h}):\boldsymbol{\epsilon}(\boldsymbol{d}_{h})$
$\displaystyle\qquad+C_{0}\int_{\Omega_{P}}p_{P,h}^{2}+\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\nabla
p_{P,h}\cdot\nabla
p_{P,h}+\frac{1}{\lambda}\int_{\Omega_{P}}(\alpha\,p_{P,h}-\varphi_{h})^{2}=0,$

which implies that $\boldsymbol{u}_{h}=\boldsymbol{0}$,
$\boldsymbol{d}_{h}=\boldsymbol{0}$, $p_{P,h}=0$, and $\varphi_{h}=0$.
Equations (4.9) imply that $\theta_{h}^{j}=0$. The inf-sup condition (4.1a)
with $q_{F,h}=p_{F,h}$ and (4.2a) give $p_{F,h}=0$. Therefore (4.10) has a
solution.

We proceed with the construction of the initial data. We first note that there
exists a solution to (4.10) satisfying $p_{P,h}(0)=p_{P,h,0}$, since this
initial condition is associated with the differential equation (4.2d). We need
to find initial values for the rest of the variables that are consistent with
the DAE system. Let us set $\theta_{h}^{j}(0)=\theta^{j}_{h,0}=0$ and consider
the Stokes problem (4.2a)–(4.2b) for $\boldsymbol{u}_{h,0}$ and $p_{F,h,0}$
with data $p_{P,h,0}$ and $\theta^{j}_{h,0}$, which is now decoupled and well-
posed from the Stokes finite element theory. Finally, let
$\boldsymbol{d}_{h,0}$ and $\varphi_{h,0}$ solve the problem coupling (4.2c)
and (4.2e) with data $p_{P,h,0}$, $\theta^{j}_{h,0}$, and
$\boldsymbol{u}_{h,0}$. The well-posedness of this problem follows from the
theory of saddle-point problems [12], due to the inf-sup condition (4.1b), see
also [45]. We further note that taking $t\to 0$ in (4.9) implies that
$\partial_{t}\boldsymbol{d}_{h}(0)\cdot\boldsymbol{t}^{j}|_{\Sigma}=\theta^{j}_{h,0}=0$.
Now, taking $t\to 0$ in (4.2a)–(4.2b), (4.2c), and (4.2e) and using the above
construction of the initial data, we conclude that
$\boldsymbol{u}_{h}(0)=\boldsymbol{u}_{h,0}$, $p_{F,h}(0)=p_{F,h,0}$,
$\boldsymbol{d}_{h}(0)=\boldsymbol{d}_{h,0}$, and
$\varphi_{h}(0)=\varphi_{h,0}$.

We proceed with the stability bound (4.8a). Differentiating (4.2e) in time and
taking
$(\boldsymbol{v}_{h},q_{F,h},\linebreak\boldsymbol{w}_{h},q_{P,h},\psi_{h})=(\boldsymbol{u}_{h},p_{F,h},\partial_{t}\boldsymbol{d}_{h},p_{P,h},\varphi_{h})$
in (4.2) gives

$\displaystyle
2\mu_{f}\\!\int_{\Omega_{F}}\\!\\!\boldsymbol{\epsilon}(\boldsymbol{u}_{h}):\boldsymbol{\epsilon}(\boldsymbol{u}_{h})+\sum_{j=1}^{d-1}\langle\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\,,(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\rangle_{\Sigma}+\frac{1}{2}\frac{\partial}{\partial
t}2\mu_{s}\int_{\Omega_{P}}\boldsymbol{\epsilon}(\boldsymbol{d}_{h}):\boldsymbol{\epsilon}(\boldsymbol{d}_{h})$
$\displaystyle\qquad+\frac{1}{2}\frac{\partial}{\partial
t}C_{0}\int_{\Omega_{P}}p_{P,h}^{2}+\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\nabla
p_{P,h}\cdot\nabla p_{P,h}+\frac{1}{2}\frac{\partial}{\partial
t}\frac{1}{\lambda}\int_{\Omega_{P}}(\alpha\,p_{P,h}-\varphi_{h})^{2}$
$\displaystyle\qquad\quad=\rho_{f}\int_{\Omega_{F}}\boldsymbol{g}\cdot\boldsymbol{u}_{h}+\rho_{s}\int_{\Omega_{P}}\boldsymbol{f}\cdot\partial_{t}\boldsymbol{d}_{h}+\rho_{f}\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\nabla
p_{P,h}-\rho_{f}\langle\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\boldsymbol{n},p_{P,h}\rangle_{\Sigma}.$
(4.11)

Integrating from $0$ to $t\in(0,T]$, we obtain

$\displaystyle
2\mu_{f}\int_{0}^{t}\|\boldsymbol{\epsilon}(\boldsymbol{u}_{h})\|_{\mathbf{L}^{2}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\gamma\mu_{f}\int_{0}^{t}\|\kappa^{-1/4}(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\mu_{s}\|\boldsymbol{\epsilon}(\boldsymbol{d}_{h})(t)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}$
$\displaystyle\quad\qquad\qquad+\frac{C_{0}}{2}\|p_{P,h}(t)\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{\mu_{f}}\int_{0}^{t}\|\kappa^{1/2}\nabla
p_{P,h}\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{2\lambda}\|(\alpha\,p_{P,h}-\varphi_{h})(t)\|^{2}_{L^{2}(\Omega_{P})}$
$\displaystyle\qquad\quad=\mu_{s}\|\boldsymbol{\epsilon}(\boldsymbol{d}_{h})(0)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{C_{0}}{2}\|p_{P,h}(0)\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{2\lambda}\|(\alpha\,p_{P,h}-\varphi_{h})(0)\|^{2}_{L^{2}(\Omega_{P})}$
$\displaystyle\qquad\quad\qquad+\rho_{f}\int_{0}^{t}\int_{\Omega_{F}}\boldsymbol{g}\cdot\boldsymbol{u}_{h}-\rho_{s}\int_{0}^{t}\int_{\Omega_{P}}\partial_{t}\boldsymbol{f}\cdot\boldsymbol{d}_{h}+\rho_{s}\int_{\Omega_{P}}\boldsymbol{f}(t)\cdot\boldsymbol{d}_{h}(t)-\rho_{s}\int_{\Omega_{P}}\boldsymbol{f}(0)\cdot\boldsymbol{d}_{h}(0)$
$\displaystyle\qquad\quad\qquad+\rho_{f}\int_{0}^{t}\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\nabla
p_{P,h}-\rho_{f}\int_{0}^{t}\langle\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\boldsymbol{n},p_{P,h}\rangle_{\Sigma},$
(4.12)

where we have integrated by parts in time the second term on the right-hand
side in (4.11). Then, on the left-hand side we use Korn’s inequality (4.5),
the Poincaré inequality (4.4), and the permeability bound (4.7), whereas on
the right-hand side we use the Cauchy-Schwarz inequality, the trace inequality
(4.6), and Young’s inequality, obtaining

$\displaystyle
2\mu_{f}(c_{K}^{F})^{2}\int_{0}^{t}\|\boldsymbol{u}_{h}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\frac{\gamma\mu_{f}}{\sqrt{\kappa}_{2}}\int_{0}^{t}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\mu_{s}(c_{K}^{P})^{2}\|\boldsymbol{d}_{h}(t)\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}$
$\displaystyle\qquad\qquad+\frac{C_{0}}{2}\|p_{P,h}(t)\|_{L^{2}(\Omega_{P})}^{2}+\frac{\kappa_{1}c_{P}^{2}}{\mu_{f}}\int_{0}^{t}\|p_{P,h}\|_{H^{1}(\Omega_{P})}^{2}+\frac{1}{2\lambda}\|(\alpha\,p_{P,h}-\varphi_{h})(t)\|^{2}_{L^{2}(\Omega_{P})}$
$\displaystyle\qquad\leq\mu_{s}\|\boldsymbol{\epsilon}(\boldsymbol{d}_{h})(0)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{C_{0}}{2}\|p_{P,h}(0)\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{2\lambda}\|(\alpha\,p_{P,h}-\varphi_{h})(0)\|^{2}_{L^{2}(\Omega_{P})}$
$\displaystyle\qquad\quad+\frac{\epsilon}{2}\left(\rho_{f}\int_{0}^{t}\|\boldsymbol{u}_{h}\|_{\mathbf{L}^{2}(\Omega_{F})}^{2}+\rho_{s}\|\boldsymbol{d}_{h}(t)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{\rho_{f}\kappa_{2}}{\mu_{f}}((C_{\Sigma}^{P})^{2}+1)\int_{0}^{t}\|p_{P,h}\|^{2}_{H^{1}(\Omega_{P})}\right)$
$\displaystyle\qquad\quad+\frac{1}{2\epsilon}\left(\rho_{f}\int_{0}^{t}\|\boldsymbol{g}\|_{\mathbf{L}^{2}(\Omega_{F})}^{2}+\rho_{s}\|\boldsymbol{f}(t)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{\rho_{f}\kappa_{2}}{\mu_{f}}\int_{0}^{t}(\|\boldsymbol{g}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\|\boldsymbol{g}\cdot\boldsymbol{n}\|_{L^{2}(\Sigma)}^{2})\right)$
$\displaystyle\qquad\quad+\frac{\rho_{s}}{2}\int_{0}^{t}\|\boldsymbol{d}_{h}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{\rho_{s}}{2}\int_{0}^{t}\|\partial_{t}\boldsymbol{f}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{\rho_{s}}{2}\|\boldsymbol{d}_{h}(0)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{\rho_{s}}{2}\|\boldsymbol{f}(0)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}.$
(4.13)

Taking $\epsilon$ sufficiently small and employing Gronwall’s inequality for
the term
$\displaystyle\frac{\rho_{s}}{2}\int_{0}^{t}\|\boldsymbol{d}_{h}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}$,
we obtain

$\displaystyle\quad\int_{0}^{t}\|\boldsymbol{u}_{h}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\int_{0}^{t}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\|\boldsymbol{d}_{h}(t)\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}$
$\displaystyle\qquad\qquad\qquad+\|p_{P,h}(t)\|_{L^{2}(\Omega_{P})}^{2}+\int_{0}^{t}\|p_{P,h}\|_{H^{1}(\Omega_{P})}^{2}+\frac{1}{\lambda}\|(\alpha\,p_{P,h}-\varphi_{h})(t)\|^{2}_{L^{2}(\Omega_{P})}$
(4.14) $\displaystyle\leq
C\biggl{(}\int_{0}^{t}(\|\boldsymbol{g}\|_{\mathbf{L}^{2}(\Omega)}^{2}+\|\boldsymbol{g}\cdot\boldsymbol{n}\|_{L^{2}(\Sigma)}^{2})+\|\boldsymbol{f}(t)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}$
$\displaystyle\qquad\qquad+\|\boldsymbol{f}(0)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\int_{0}^{t}\|\partial_{t}\boldsymbol{f}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\|p_{P,0}\|_{H^{1}(\Omega_{P})}^{2}\biggr{)},$

with a constant $C$ independent of $\lambda_{\max}$. In the above inequality
we have bounded the initial data terms by
$C\|p_{P,0}\|_{H^{1}(\Omega_{P})}^{2}$. This bound follows from the classical
stability bound for the Stokes problem (4.2a)–(4.2b), which allows to obtain
$\|\boldsymbol{u}_{h,0}\|_{\mathbf{H}^{1}(\Omega_{F})}\leq
C\|p_{P,h,0}\|_{H^{1}(\Omega_{P})}$; a stability bound for the saddle-point
problem (4.2c), (4.2e) to obtain

$\|\boldsymbol{d}_{h,0}\|_{\mathbf{H}^{1}(\Omega_{P})}+\|\varphi_{h,0}\|_{L^{2}(\Omega_{p})}\leq
C(\|\boldsymbol{u}_{h,0}\|_{\mathbf{H}^{1}(\Omega_{F})}+\|p_{P,h,0}\|_{H^{1}(\Omega_{P})}),$

(cf. [45]), and the $H^{1}$-stability of the $L^{2}$-projection
$\|p_{P,h,0}\|_{H^{1}(\Omega_{P})}\leq C\|p_{P,0}\|_{H^{1}(\Omega_{P})}$ (see,
e.g., [14]).

Next we proceed with bounding $p_{F,h}$ and $\varphi_{h}$. The inf-sup
condition (4.1a) together with (4.2a) give

$\displaystyle\beta_{F}\|p_{F,h}\|_{L^{2}(\Omega_{F})}$
$\displaystyle\leq\sup_{\boldsymbol{0}\neq\boldsymbol{v}_{h}\in\mathbf{V}_{h}^{\Sigma}}\frac{b_{1}^{F}(\boldsymbol{v}_{h},p_{F,h})}{\|\boldsymbol{v}_{h}\|_{\mathbf{H}^{1}(\Omega_{F})}}$
$\displaystyle=\sup_{\boldsymbol{0}\neq\boldsymbol{v}_{h}\in\mathbf{V}_{h}}\frac{-a_{2}^{F}(\boldsymbol{u}_{h},\boldsymbol{v}_{h})-b^{\Sigma}_{2}(\boldsymbol{v}_{h},p_{P,h})-b_{3}^{\Sigma}(\boldsymbol{v}_{h},\partial_{t}\boldsymbol{d}_{h})+F^{F}(\boldsymbol{v}_{h})}{\|\boldsymbol{v}_{h}\|_{\mathbf{H}^{1}(\Omega_{F})}}$
$\displaystyle\leq
2\mu_{f}\|\boldsymbol{u}_{h}\|_{\mathbf{H}^{1}(\Omega_{F})}+\sum_{j=1}^{d-1}\frac{\gamma\mu_{f}C_{\Sigma}^{F}}{\sqrt{\kappa}_{1}}\\!\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{\mathbf{L}^{2}(\Sigma)}$
$\displaystyle\qquad+C_{\Sigma}^{F}C_{\Sigma}^{P}\|p_{P,h}\|_{H^{1}(\Omega_{P})}\\!+\rho_{f}\|\boldsymbol{g}\|_{\mathbf{L}^{2}(\Omega_{F})},$

implying

$\int_{0}^{t}\|p_{F,h}\|_{L^{2}(\Omega_{F})}^{2}\leq
C\int_{0}^{t}\Big{(}\|\boldsymbol{u}_{h}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\|p_{P,h}\|_{H^{1}(\Omega_{P})}^{2}+\|\boldsymbol{g}\|_{\mathbf{L}^{2}(\Omega_{F})}^{2}\Big{)}.$
(4.15)

Finally, using the inf-sup condition (4.1b) and (4.2c), we obtain

$\displaystyle\beta_{P}\|\varphi_{h}\|_{L^{2}(\Omega_{P})}$
$\displaystyle\leq\sup_{\boldsymbol{0}\neq\boldsymbol{w}_{h}\in\mathbf{W}_{h}}\frac{b_{1}^{P}(\boldsymbol{w}_{h},\varphi_{h})}{\|\boldsymbol{w}_{h}\|_{\mathbf{H}^{1}(\Omega_{P})}}$
$\displaystyle=\sup_{\boldsymbol{0}\neq\boldsymbol{w}_{h}\in\mathbf{W}_{h}}\\!\\!\frac{-b_{3}^{\Sigma}(\boldsymbol{u}_{h},\boldsymbol{w}_{h})-b_{4}^{\Sigma}(\boldsymbol{w}_{h},p_{P,h})-a_{1}^{P}(\boldsymbol{d}_{h},\boldsymbol{w}_{h})-a_{2}^{\Sigma}(\partial_{t}\boldsymbol{d}_{h},\boldsymbol{w}_{h})+F^{P}(\boldsymbol{w}_{h})}{\|\boldsymbol{w}_{h}\|_{\mathbf{H}^{1}(\Omega_{P})}}$
$\displaystyle\leq(C_{\Sigma}^{P})^{2}\|p_{P,h}\|_{H^{1}(\Omega_{P})}+\sum_{j=1}^{d-1}\frac{\gamma\mu_{f}C_{\Sigma}^{P}}{\sqrt{\kappa}_{1}}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{\mathbf{L}^{2}(\Sigma)}$
$\displaystyle\qquad+2\mu_{s}\|\boldsymbol{d}_{h}\|_{\mathbf{H}^{1}(\Omega_{P})}+\rho_{s}\|\boldsymbol{f}\|_{\mathbf{L}^{2}(\Omega_{P})},$

yielding

$\int_{0}^{t}\|\varphi_{h}\|_{L^{2}(\Omega_{P})}^{2}\leq
C\int_{0}^{t}\Big{(}\|p_{P,h}\|_{H^{1}(\Omega_{P})}^{2}+\sum_{j=1}^{d-1}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\|\boldsymbol{d}_{h}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\boldsymbol{f}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}\Big{)}.$
(4.16)

Combining (4.14)–(4.16) and using Gronwall’s inequality for the third term on
the right-hand side in (4.16), we obtain (4.8a).

The above argument implies that the solution of (4.3) under the initial
conditions $p_{P,h}(0)=p_{P,h,0}$,
$\boldsymbol{d}_{h}(0)=\boldsymbol{d}_{h,0}$, and
$\varphi_{h}(0)=\varphi_{h,0}$ is unique. In particular, taking $p_{P,h,0}=0$,
$\boldsymbol{d}_{h,0}=\boldsymbol{0}$, $\varphi_{h,0}=0$,
$\boldsymbol{g}=\boldsymbol{0}$, and $\boldsymbol{f}=\boldsymbol{0}$, (4.13)
implies that (4.14) holds with right-hand side zero. Together with (4.15) and
(4.16), this gives that all components of the solution are zero, therefore the
solution is unique.

We next prove the higher regularity stability bound (4.8b). To that end, we
differentiate in time (4.2a), (4.2c), and (4.2e) and take
$(\boldsymbol{v}_{h},q_{F,h},\boldsymbol{w}_{h},q_{P,h},\psi_{h})=(\boldsymbol{u}_{h},\partial_{t}p_{F,h},\partial_{t}\boldsymbol{d}_{h},\partial_{t}p_{P,h},\partial_{t}\varphi_{h})$
in (4.2), obtaining

$\displaystyle\frac{1}{2}\frac{\partial}{\partial
t}2\mu_{f}\\!\int_{\Omega_{F}}\\!\\!\\!\boldsymbol{\epsilon}(\boldsymbol{u}_{h}):\boldsymbol{\epsilon}(\boldsymbol{u}_{h})+\sum_{j=1}^{d-1}\frac{1}{2}\frac{\partial}{\partial
t}\langle\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\\!\cdot\boldsymbol{t}^{j},(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\\!\cdot\boldsymbol{t}^{j}\rangle_{\Sigma}$
$\displaystyle\quad+2\mu_{s}\int_{\Omega_{P}}\\!\\!\boldsymbol{\epsilon}(\partial_{t}\boldsymbol{d}_{h}):\boldsymbol{\epsilon}(\partial_{t}\boldsymbol{d}_{h})+C_{0}\int_{\Omega_{P}}(\partial_{t}p_{P,h})^{2}+\frac{1}{2}\frac{\partial}{\partial
t}\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\nabla p_{P,h}\cdot\nabla p_{P,h}$
$\displaystyle\quad+\frac{1}{\lambda}\int_{\Omega_{P}}(\alpha\,\partial_{t}p_{P,h}-\partial_{t}\varphi_{h})^{2}=\rho_{s}\int_{\Omega_{P}}\partial_{t}\boldsymbol{f}\cdot\partial_{t}\boldsymbol{d}_{h}+\rho_{f}\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\partial_{t}\nabla
p_{P,h}-\rho_{f}\langle\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\boldsymbol{n},\partial_{t}p_{P,h}\rangle_{\Sigma},$
(4.17)

where we have used that $\partial_{t}\boldsymbol{g}=\boldsymbol{0}$.
Integration from $0$ to $t\in(0,T]$ gives

$\displaystyle\mu_{f}\|\boldsymbol{\epsilon}(\boldsymbol{u}_{h})(t)\|_{\mathbf{L}^{2}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\frac{1}{2}\gamma\mu_{f}\|\kappa^{-1/4}(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}(t)\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\mu_{s}\int_{0}^{t}\|\boldsymbol{\epsilon}(\partial_{t}\boldsymbol{d}_{h})\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}$
$\displaystyle\qquad+C_{0}\int_{0}^{t}\|\partial_{t}p_{P,h}(t)\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{2\mu_{f}}\|\kappa^{1/2}\nabla
p_{P,h}(t)\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{\lambda}\int_{0}^{t}\|\alpha\,\partial_{t}p_{P,h}-\partial_{t}\varphi_{h}\|^{2}_{L^{2}(\Omega_{P})}$
$\displaystyle=\mu_{f}\|\boldsymbol{\epsilon}(\boldsymbol{u}_{h})(0)\|_{\mathbf{L}^{2}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\frac{1}{2}\gamma\mu_{f}\|\kappa^{-1/4}(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}(0)\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\frac{1}{2\mu_{f}}\|\kappa^{1/2}\nabla
p_{P,h}(0)\|_{L^{2}(\Omega_{P})}^{2}$ $\displaystyle\
+\rho_{s}\int_{0}^{t}\int_{\Omega_{P}}\partial_{t}\boldsymbol{f}\cdot\partial_{t}\boldsymbol{d}_{h}+\rho_{f}\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot(\nabla
p_{P,h}(t)-\nabla
p_{P,h}(0))-\rho_{f}\langle\frac{\kappa}{\mu_{f}}\boldsymbol{g}\cdot\boldsymbol{n},p_{P,h}(t)-p_{P,h}(0)\rangle_{\Sigma},$

where we have integrated by parts the last two terms in (4.17) and used that
$\partial_{t}\boldsymbol{g}=\boldsymbol{0}$. Next, on the left-hand side we
use Korn’s inequality (4.5), the Poincaré inequality (4.4), and the
permeability bound (4.7), while on the right-hand side we invoke Cauchy-
Schwarz inequality, the trace inequality (4.6), and Young’s inequality,
yielding

$\displaystyle\mu_{f}(c_{K}^{F})^{2}\|\boldsymbol{u}_{h}(t)\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\frac{1}{2}\frac{\gamma\mu_{f}}{\sqrt{\kappa}_{2}}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}(t)\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\mu_{s}(c_{K}^{P})^{2}\int_{0}^{t}\|\partial_{t}\boldsymbol{d}_{h}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}$
$\displaystyle\qquad\quad\qquad+C_{0}\int_{0}^{t}\|\partial_{t}p_{P,h}(t)\|_{L^{2}(\Omega_{P})}^{2}+\frac{\kappa_{1}c_{P}^{2}}{2\mu_{f}}\|p_{P,h}(t)\|_{H^{1}(\Omega_{P})}^{2}+\frac{1}{\lambda}\int_{0}^{t}\|\alpha\,\partial_{t}p_{P,h}-\partial_{t}\varphi_{h}\|^{2}_{L^{2}(\Omega_{P})}$
$\displaystyle\qquad\quad\leq\mu_{f}\|\boldsymbol{\epsilon}(\boldsymbol{u}_{h})(0)\|_{\mathbf{L}^{2}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\frac{1}{2}\frac{\gamma\mu_{f}}{\sqrt{\kappa}_{1}}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}(0)\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\frac{\kappa_{2}}{2\mu_{f}}\|\nabla
p_{P,h}(0)\|_{L^{2}(\Omega_{P})}^{2}$
$\displaystyle\qquad\quad\qquad+\frac{\epsilon}{2}\left(\rho_{s}\int_{0}^{t}\|\partial_{t}\boldsymbol{d}_{h}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{\rho_{f}\kappa_{2}}{\mu_{f}}((C_{\Sigma}^{P})^{2}+1)(\|p_{P,h}(0)\|^{2}_{H^{1}(\Omega_{P})}+\|p_{P,h}(t)\|^{2}_{H^{1}(\Omega_{P})})\right)$
$\displaystyle\qquad\quad\qquad+\frac{1}{2\epsilon}\left(\rho_{s}\int_{0}^{t}\|\partial_{t}\boldsymbol{f}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\frac{\rho_{f}\kappa_{2}}{\mu_{f}}(\|\boldsymbol{g}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\|\boldsymbol{g}\cdot\boldsymbol{n}\|_{L^{2}(\Sigma)}^{2})\right).$

In addition, bounding the initial data terms as in (4.14) and taking
$\epsilon$ sufficiently small, we can assert that

$\displaystyle\|\boldsymbol{u}_{h}(t)\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}(t)\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\int_{0}^{t}\|\partial_{t}\boldsymbol{d}_{h}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}$
$\displaystyle\qquad\qquad\quad+\int_{0}^{t}\|\partial_{t}p_{P,h}(t)\|_{L^{2}(\Omega_{P})}^{2}+\|p_{P,h}(t)\|_{H^{1}(\Omega_{P})}^{2}+\frac{1}{\lambda}\int_{0}^{t}\|\alpha\,\partial_{t}p_{P,h}-\partial_{t}\varphi_{h}\|^{2}_{L^{2}(\Omega_{P})}$
$\displaystyle\qquad\quad\leq
C\left(\int_{0}^{t}\|\partial_{t}\boldsymbol{f}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\|\boldsymbol{g}\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}+\|\boldsymbol{g}\cdot\boldsymbol{n}\|_{L^{2}(\Sigma)}^{2}+\|p_{P,0}\|_{H^{1}(\Omega_{P})}^{2}\right),$
(4.18)

with a constant $C$ independent of $\lambda_{\max}$. Next, using the inf-sup
conditions (4.1a) and (4.1b), and proceeding similarly to the derivations of
(4.15) and (4.16), we obtain

$\|p_{F,h}(t)\|_{L^{2}(\Omega_{F})}^{2}\\!\leq
C\\!\Big{(}\|\boldsymbol{u}_{h}(t)\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\|(\boldsymbol{u}_{h}\\!-\\!\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}(t)\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\|p_{P,h}(t)\|_{H^{1}(\Omega_{P})}^{2}+\|\boldsymbol{g}\|_{\mathbf{L}^{2}(\Omega_{F})}^{2}\Big{)},$

and

$\|\varphi_{h}(t)\|_{L^{2}(\Omega_{P})}^{2}\leq\\!C\\!\Big{(}\|p_{P,h}(t)\|_{H^{1}(\Omega_{P})}^{2}+\sum_{j=1}^{d-1}\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}(t)\|_{\mathbf{L}^{2}(\Sigma)}^{2}\\!+\\!\|\boldsymbol{d}_{h}(t)\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\boldsymbol{f}(t)\|_{\mathbf{L}^{2}(\Omega_{P})}^{2}\Big{)}.$
(4.19)

Finally, combining (4.18)–(4.19) and employing (4.8a) for the control of
$\|\boldsymbol{d}_{h}(t)\|_{\mathbf{H}^{1}(\Omega_{P})}$, we obtain the second
bound (4.8b). $\square$

###### Remark 4.1

We emphasise that, even though initial data was constructed for all variables,
the initial value problem for (4.3) involves initial conditions only for
$p_{P,h}$, $\boldsymbol{d}_{h}$, and $\varphi_{h}$.

###  4.3 Existence, uniqueness, and stability of the weak solution

###### Theorem 4.2

For each $\boldsymbol{f}\in H^{1}(0,T;\mathbf{L}^{2}(\Omega_{P}))$ and
$p_{P,0}\in H^{1}_{\star}(\Omega_{p})$, there exist initial data
$\boldsymbol{u}_{0}\in\mathbf{H}^{1}_{\star}(\Omega_{F})$, $p_{F,0}\in
L^{2}(\Omega_{F})$, $\boldsymbol{d}_{0}\in\mathbf{H}^{1}_{\star}(\Omega_{P})$,
and $\varphi_{0}\in L^{2}(\Omega_{P})$ such that the weak formulation (3.1)
with $a_{1}^{F}=0$, $c^{F}=0$, and $\tilde{\alpha}=1$, complemented with the
initial conditions $p_{P}(0)=p_{P,0}$, $\boldsymbol{d}(0)=\boldsymbol{d}_{0}$,
and $\varphi(0)=\varphi_{0}$, has a unique solution satisfying

$\displaystyle\|\boldsymbol{u}\|_{L^{2}(0,T;\mathbf{H}^{1}(\Omega_{F}))}+\|p_{F}\|_{L^{2}(0,T;L^{2}(\Omega_{F}))}+\|\boldsymbol{d}\|_{L^{\infty}(0,T;\mathbf{H}^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\sum_{j=1}^{d-1}\|(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\cdot\boldsymbol{t}^{j}\|_{L^{2}(0,T;L^{2}(\Sigma))}+\|p_{P}\|_{L^{\infty}(0,T;L^{2}(\Omega_{P}))}+\|p_{P}\|_{L^{2}(0,T;H^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\|\varphi\|_{L^{2}(0,T;L^{2}(\Omega_{P}))}+\frac{1}{\sqrt{\lambda}}\|\alpha\,p_{P}-\varphi\|_{L^{\infty}(0,T;L^{2}(\Omega_{P}))}$
$\displaystyle\leq
C\left(\|\boldsymbol{g}\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega))}+\|\boldsymbol{g}\cdot\boldsymbol{n}\|_{L^{2}(0,T;L^{2}(\Sigma))}+\|\boldsymbol{f}\|_{H^{1}(0,T;\mathbf{L}^{2}(\Omega_{P}))}+\|p_{P,0}\|_{H^{1}(\Omega_{P})}\right),$
(4.20a) and
$\displaystyle\|\boldsymbol{u}\|_{L^{\infty}(0,T;\mathbf{H}^{1}(\Omega_{F}))}+\|p_{F}\|_{L^{\infty}(0,T;L^{2}(\Omega_{F}))}+\|\partial_{t}\boldsymbol{d}\|_{L^{2}(0,T;\mathbf{H}^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\sum_{j=1}^{d-1}\|(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\cdot\boldsymbol{t}^{j}\|_{L^{\infty}(0,T;L^{2}(\Sigma))}+\|\partial_{t}p_{P}\|_{L^{2}(0,T;L^{2}(\Omega_{P}))}+\|p_{P}\|_{L^{\infty}(0,T;H^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\|\varphi\|_{L^{\infty}(0,T;L^{2}(\Omega_{P}))}+\frac{1}{\sqrt{\lambda}}\|\alpha\,\partial_{t}p_{P}-\partial_{t}\varphi\|_{L^{2}(0,T;L^{2}(\Omega_{P}))}$
$\displaystyle\leq
C\left(\|\boldsymbol{g}\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega))}+\|\boldsymbol{g}\cdot\boldsymbol{n}\|_{L^{2}(0,T;L^{2}(\Sigma))}+\|\boldsymbol{f}\|_{H^{1}(0,T;\mathbf{L}^{2}(\Omega_{P}))}+\|p_{P,0}\|_{H^{1}(\Omega_{P})}\right),$
(4.20b)

with a constant $C$ independent of $\lambda_{\max}$.

Proof. From Theorem 4.1, there exists an infinite sequence
$\\{(\boldsymbol{u}_{h},p_{F,h},\boldsymbol{d}_{h},p_{P,h},\varphi_{h})\\}_{h>0}$
satisfying (4.3) for each $h$ such that $\\{\boldsymbol{u}_{h}\\}_{h>0}$ is
bounded in $L^{2}(0,T;\mathbf{H}^{1}(\Omega_{F}))$, $\\{p_{F,h}\\}_{h>0}$ is
bounded in $L^{2}(0,T;\\-L^{2}(\Omega_{F}))$, $\\{\boldsymbol{d}_{h}\\}_{h>0}$
is bounded in $H^{1}(0,T;\mathbf{H}^{1}(\Omega_{P}))$, the sequence
$\\{p_{P,h}\\}_{h>0}$ is bounded in $H^{1}(0,T;L^{2}(\Omega_{P}))$, as well as
in $L^{2}(0,T;H^{1}(\Omega_{P}))$, and the sequence $\\{\varphi_{h}\\}_{h>0}$
is bounded in $H^{1}(0,T;L^{2}(\Omega_{P}))$. Therefore there exist weakly
convergent subsequences, denoted in the same way, such that

$\displaystyle\boldsymbol{u}_{h}\rightharpoonup\boldsymbol{u}\text{ in
}L^{2}(0,T;\mathbf{H}^{1}(\Omega_{F})),\quad p_{F,h}\rightharpoonup
p_{F}\text{ in }L^{2}(0,T;L^{2}(\Omega_{F})),$
$\displaystyle\boldsymbol{d}_{h}\rightharpoonup\boldsymbol{d}\text{ in
}H^{1}(0,T;\mathbf{H}^{1}(\Omega_{P})),$ $\displaystyle p_{P,h}\rightharpoonup
p_{P}\text{ in }H^{1}(0,T;L^{2}(\Omega_{P}))\cap
L^{2}(0,T;H^{1}(\Omega_{P})),\quad\varphi_{h}\rightharpoonup\varphi\text{ in
}H^{1}(0,T;L^{2}(\Omega_{P})).$

Next, we fix a set of test functions
$(\boldsymbol{v}_{h},q_{F,h},\boldsymbol{w}_{h},q_{P,h},\psi_{h})\in
C^{0}(0,T;\mathbf{V}_{h}\times\mathrm{Q}^{F}_{h}\times\mathbf{W}_{h}\times\mathrm{Q}^{P}_{h}\times\mathrm{Z}_{h})$
in (4.2), integrate it in time from 0 to $T$, and take $h\to 0$. Since all
bilinear forms and linear functionals are continuous, cf. Lemma 4.1, we
conclude that $(\boldsymbol{u},p_{F},\boldsymbol{d},p_{P},\varphi)$ satisfy
the time-integrated version of (3.1) with this choice of test functions.
Finally, since the space
$C^{0}(0,T;\mathbf{V}_{h}\times\mathrm{Q}^{F}_{h}\times\mathbf{W}_{h}\times\mathrm{Q}^{P}_{h}\times\mathrm{Z}_{h})$
is dense in $L^{2}(0,T;\mathbf{H}^{1}_{\star}(\Omega_{F})\times
L^{2}(\Omega_{F})\times\mathbf{H}^{1}_{\star}(\Omega_{P})\times
H^{1}_{\star}(\Omega_{P})\times L^{2}(\Omega_{P}))$, we conclude that (3.1)
holds for a.e. $t\in(0,T)$.

It remains to handle the initial conditions. First, taking $h\to 0$ in
$p_{P,h}(0)=p_{P,h,0}$ gives $p_{P}(0)=p_{P,0}$. We further note that the
control of the terms
$\|(\boldsymbol{u}_{h}-\partial_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{L^{\infty}(0,T;L^{2}(\Sigma))}$
and $\|\boldsymbol{u}_{h}\|_{L^{\infty}(0,T;\mathbf{H}^{1}(\Omega_{F}))}$ in
(4.8b) implies that for all $t\in[0,T]$,
$\partial_{t}\boldsymbol{d}_{h}(t)\cdot\boldsymbol{t}^{j}\rightharpoonup\partial_{t}\boldsymbol{d}(t)\cdot\boldsymbol{t}^{j}$
in $L^{2}(\Sigma)$. Taking $t=0$ and $h\to 0$ in (4.9) and using that
$\partial_{t}\boldsymbol{d}_{h}(0)\cdot\boldsymbol{t}^{j}=0$ on $\Sigma$, we
conclude that $\partial_{t}\boldsymbol{d}(0)\cdot\boldsymbol{t}^{j}=0$ on
$\Sigma$. Next, the stability of the Stokes and elasticity problems for the
initial data in the proof of Theorem 4.1 implies that there exist weakly
convergent subsequences such that

$\displaystyle\boldsymbol{u}_{h,0}\rightharpoonup\boldsymbol{u}_{0}\text{ in
}\mathbf{H}^{1}(\Omega_{F}),\ \ p_{F,h,0}\rightharpoonup p_{F,0}\text{ in
}L^{2}(\Omega_{F}),\ \
\boldsymbol{d}_{h,0}\rightharpoonup\boldsymbol{d}_{0}\text{ in
}\mathbf{H}^{1}(\Omega_{P}),\ \ \varphi_{h,0}\rightharpoonup\varphi_{0}\text{
in }L^{2}(\Omega_{P}).$

Then, taking $t\to 0$ in (3.1a)–(3.1b), (3.1c), and (3.1e) and using that the
initial data satisfies the same equations, we conclude that
$\boldsymbol{u}(0)=\boldsymbol{u}_{0}$, $p_{F}(0)=p_{F,0}$,
$\boldsymbol{d}(0)=\boldsymbol{d}_{0}$, and $\varphi(0)=\varphi_{0}$.

Finally, the uniqueness of the solution under the initial conditions
$p_{P}(0)=p_{P,0}$, $\boldsymbol{d}(0)=\boldsymbol{d}_{0}$,
$\varphi(0)=\varphi_{0}$, and the stability bounds (4.20a) and (4.20b), follow
in the same way as in the proof of Theorem 4.1. $\square$

##  5 Fully discrete scheme

We apply a time discretisation to (4.2) using backward Euler’s method with
fixed time step $\Delta t=T/N$. Let $t_{n}=n\Delta t$, $n=0,\ldots,N$, be the
discrete times. Starting from the discrete initial data constructed in the
proof of Theorem 4.1, at each time iteration $n=1,\ldots,N$ we look for
$(\boldsymbol{u}^{n}_{h},p^{n}_{F,h},\boldsymbol{d}^{n}_{h},p^{n}_{P,h},\varphi^{n}_{h})\in\mathbf{V}_{h}\times\mathrm{Q}^{F}_{h}\times\mathbf{W}_{h}\times\mathrm{Q}^{P}_{h}\times\mathrm{Z}_{h}=:\mathbf{H}_{h}$
such that

$\displaystyle
a_{2}^{F}(\boldsymbol{u}_{h}^{n},\boldsymbol{v}_{h})+b^{F}_{1}(\boldsymbol{v}_{h},p_{F,h}^{n})+b^{\Sigma}_{2}(\boldsymbol{v}_{h},p_{P,h}^{n})+b_{3}^{\Sigma}(\boldsymbol{v}_{h},\partial_{t}^{n}\boldsymbol{d}_{h})$
$\displaystyle=F^{F,n}(\boldsymbol{v}_{h})$
$\displaystyle\\!\forall\boldsymbol{v}_{h}\in\mathbf{V}_{h},$ (5.1a)
$\displaystyle-b^{F}_{1}(\boldsymbol{u}_{h}^{n},q_{F,h})$ $\displaystyle=0$
$\displaystyle\\!\forall q_{F,h}\in\mathrm{Q}_{h}^{F},$ (5.1b) $\displaystyle
b_{3}^{\Sigma}(\boldsymbol{u}_{h}^{n},\boldsymbol{w}_{h})+b_{4}^{\Sigma}(\boldsymbol{w}_{h},p_{P,h}^{n})+a_{1}^{P}(\boldsymbol{d}_{h}^{n},\boldsymbol{w}_{h})$
$\displaystyle+a_{2}^{\Sigma}(\partial_{t}^{n}\boldsymbol{d}_{h},\boldsymbol{w}_{h})+b_{1}^{P}(\boldsymbol{w}_{h},\varphi_{h}^{n})$
$\displaystyle=F^{P}(\boldsymbol{w}_{h})$
$\displaystyle\\!\forall\boldsymbol{w}_{h}\in\mathbf{W}_{h},$ (5.1c)
$\displaystyle-
b_{2}^{\Sigma}(\boldsymbol{u}_{h}^{n},q_{P,h})-b_{4}^{\Sigma}(\partial_{t}^{n}\boldsymbol{d}_{h},q_{P,h})+a_{3}^{P}(\partial_{t}^{n}p_{P,h},q_{P,h})$
$\displaystyle+a_{4}^{P}(p_{P,h}^{n},q_{P,h})-b_{2}^{P}(\partial_{t}^{n}\varphi_{h},q_{P,h})$
$\displaystyle=G^{n}(q_{P,h})$ $\displaystyle\\!\forall
q_{P,h}\in\mathrm{Q}_{h}^{P},$ (5.1d) $\displaystyle-
b_{1}^{P}(\boldsymbol{d}_{h}^{n},\psi_{h})-b_{2}^{P}(\psi_{h},p_{P,h}^{n})+a_{5}^{P}(\varphi_{h}^{n},\psi_{h})$
$\displaystyle=0$ $\displaystyle\\!\forall\psi_{h}\in\mathrm{Z}_{h},$ (5.1e)

where, for a generic scalar or vector field $f$, we set $f^{n}:=f(t_{n})$ and
$\partial_{t}^{n}f:=\frac{1}{\Delta t}(f^{n}-f^{n-1})$. For convenience we
also define the global discrete time derivative $\tilde{\partial}_{t}f$ such
that $(\tilde{\partial}_{t}f)^{n}:=\partial_{t}^{n}f$ for $n=1,\ldots,N$. The
method requires solving at each time step the algebraic system

$\left[\begin{array}[]{cc:ccc}\mathcal{A}_{2}^{F}&(\mathcal{B}_{1}^{F})^{\prime}&\frac{1}{\Delta
t}(\mathcal{B}_{3}^{\Sigma})^{\prime}&(\mathcal{B}_{2}^{\Sigma})^{\prime}&\boldsymbol{0}\\\\[4.30554pt]
-\mathcal{B}_{1}^{F}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\hdashline\mathcal{B}_{3}^{\Sigma}&\boldsymbol{0}&\mathcal{A}_{1}^{P}+\frac{1}{\Delta
t}\mathcal{A}_{2}^{\Sigma}&(\mathcal{B}_{4}^{\Sigma})^{\prime}&(\mathcal{B}_{1}^{P})^{\prime}\vphantom{\int^{X^{X}}}\\\\[4.30554pt]
-\mathcal{B}_{2}^{\Sigma}&\boldsymbol{0}&-\frac{1}{\Delta
t}\mathcal{B}_{4}^{\Sigma}&\frac{1}{\Delta
t}\mathcal{A}_{3}^{P}+\mathcal{A}_{4}^{P}&-\frac{1}{\Delta
t}(\mathcal{B}_{2}^{P})^{\prime}\\\\[4.30554pt]
\boldsymbol{0}&\boldsymbol{0}&-\mathcal{B}_{1}^{P}&-\mathcal{B}_{2}^{P}&\mathcal{A}_{5}^{P}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{u}^{n}_{h}\\\\[4.30554pt]
p^{n}_{F,h}\\\\[4.30554pt] \boldsymbol{d}^{n}_{h}\\\\[4.30554pt]
p^{n}_{P,h}\\\\[4.30554pt]
\varphi^{n}_{h}\end{array}\right]=\left[\begin{array}[]{c}\tilde{\mathcal{F}}^{F,n}\\\\[4.30554pt]
\boldsymbol{0}\\\\[4.30554pt] \tilde{\mathcal{F}}^{P,n}\\\\[4.30554pt]
\tilde{\mathcal{G}}^{n}\\\\[4.30554pt] \boldsymbol{0}\end{array}\right],$
(5.2)

where the tilde-notation on the right-hand side vectors indicate that they
also receive contributions from the backward Euler time-discretisation.

###### Theorem 5.1

The fully discrete method (5.2) has a unique solution.

Proof. Consider the matrix obtained from the matrix in (5.2) by scaling the
third and fifth rows by $\frac{1}{\Delta t}$. It has the same structure as the
matrix $\mathcal{N}+\mathcal{M}$ that appears in the proof of Theorem 4.1,
which is shown to be nonsingular. Therefore the scaled matrix is nonsingular,
and so is the matrix in (5.2). $\square$

###  5.1 Error estimates for the fully discrete scheme

We proceed with the error analysis for the fully discrete scheme. We will make
use of the discrete space-time norms for $f:t_{n}\to\mathrm{V}$,
$n=1,\ldots,N$,

$\|f\|^{2}_{l^{2}(0,T;\mathrm{V})}:=\sum_{n=1}^{N}\Delta
t\|f^{n}\|_{\mathrm{V}}^{2},\quad\|f\|_{l^{\infty}(0,T;\mathrm{V})}:=\max\limits_{n=1,\ldots,N}\|f^{n}\|_{\mathrm{V}},$

In addition, we consider finite-dimensional subspaces of continuous and
piecewise polynomials of generic degrees

$k_{\boldsymbol{v}},k_{\boldsymbol{w}},k_{q_{P}}\geq 1,\quad
k_{q_{F}},k_{z}\geq 0,\quad\text{ for the
spaces}\quad\mathbf{V}_{h},\mathbf{W}_{h},\mathrm{Q}_{h}^{P},\mathrm{Q}_{h}^{F},\mathrm{Z}_{h},$
(5.3)

respectively.

Let $I^{\mathbf{V}}:\mathbf{H}^{1}_{\star}(\Omega_{F})\to\mathbf{V}_{h}$,
$I^{\mathbf{W}}:\mathbf{H}^{1}_{\star}(\Omega_{P})\to\mathbf{W}_{h}$, and
$I^{\mathrm{Q}^{P}}:H^{1}_{\star}(\Omega_{P})\to\mathrm{Q}_{h}^{P}$ be the
Scott-Zhang interpolants [48]. In addition, let
$I^{\mathrm{Q}^{F}}:L^{2}(\Omega_{F})\to\mathrm{Q}_{h}^{F}$ and
$I^{\mathrm{Z}}:L^{2}(\Omega_{P})\to\mathrm{Z}_{h}$ be the $L^{2}$-orthogonal
projections. These operators, alongside the polynomial degrees (5.3) have the
approximation properties (see, e.g., [48, 22])

$\displaystyle\|\boldsymbol{v}-I^{\mathbf{V}}\boldsymbol{v}\|_{\mathbf{H}^{s}(\Omega_{F})}\leq
Ch^{r_{\boldsymbol{v}}-s}\|\boldsymbol{v}\|_{\mathbf{H}^{r_{\boldsymbol{v}}}(\Omega_{F})},$
$\displaystyle 1\leq r_{\boldsymbol{v}}\leq k_{\boldsymbol{v}}+1,\ 0\leq s\leq
1,$ (5.4a)
$\displaystyle\|q_{F}-I^{\mathrm{Q}^{F}}q_{F}\|_{L^{2}(\Omega_{P})}\leq
Ch^{r_{q_{F}}}\|q_{F}\|_{H^{r_{q_{F}}}(\Omega_{F})},$ $\displaystyle 0\leq
r_{q_{F}}\leq k_{q_{F}}+1,$ (5.4b)
$\displaystyle\|\boldsymbol{w}-I^{\mathbf{W}}\boldsymbol{w}\|_{\mathbf{H}^{s}(\Omega_{P})}\leq
Ch^{r_{\boldsymbol{w}}-s}\|\boldsymbol{w}\|_{\mathbf{H}^{r_{\boldsymbol{w}}}(\Omega_{P})},$
$\displaystyle 1\leq r_{\boldsymbol{w}}\leq k_{\boldsymbol{w}}+1,\ 0\leq s\leq
1,$ (5.4c)
$\displaystyle\|q_{P}-I^{\mathrm{Q}^{P}}q_{P}\|_{H^{s}(\Omega_{P})}\leq
Ch^{r_{q_{P}}-s}\|q_{P}\|_{H^{r_{q_{P}}}(\Omega_{P})},$ $\displaystyle 1\leq
r_{q_{P}}\leq k_{q_{P}}+1,\ 0\leq s\leq 1,$ (5.4d)
$\displaystyle\|\psi-I^{\mathrm{Z}}\psi\|_{L^{2}(\Omega_{P})}\leq
Ch^{r_{z}}\|\psi\|_{H^{r_{z}}(\Omega_{P})},$ $\displaystyle 0\leq r_{z}\leq
k_{z}+1.$ (5.4e)

###### Theorem 5.2

Assume that the weak solution of (3.1) is sufficiently smooth. Then, for the
fully discrete solution (5.1), there exists a constant $C$ independent of
$\lambda_{\max}$, $h$, and $\Delta t$, such that

$\displaystyle\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{l^{2}(0,T;\mathbf{H}^{1}(\Omega_{F}))}+\|p_{F}-p_{F,h}\|_{l^{2}(0,T;L^{2}(\Omega_{F}))}+\|\boldsymbol{d}-\boldsymbol{d}_{h}\|_{l^{\infty}(0,T;\mathbf{H}^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\sum_{j=1}^{d-1}\|(\boldsymbol{u}-\tilde{\partial}_{t}\boldsymbol{d}_{h})\cdot\boldsymbol{t}^{j}\|_{l^{2}(0,T;L^{2}(\Sigma))}+\|p_{P}-p_{P,h}\|_{l^{\infty}(0,T;L^{2}(\Omega_{P}))}+\|p_{P}-p_{P,h}\|_{l^{2}(0,T;H^{1}(\Omega_{P}))}$
$\displaystyle\qquad+\|\varphi-\varphi_{h}\|_{l^{2}(0,T;L^{2}(\Omega_{P}))}+\frac{1}{\sqrt{\lambda}}\|(\alpha\,p_{P}-\varphi)-(\alpha\,p_{P,h}-\varphi_{h})\|_{l^{\infty}(0,T;L^{2}(\Omega_{P}))}$
$\displaystyle\quad\leq
C\Big{(}h^{r_{\boldsymbol{v}}}\|\boldsymbol{u}\|_{H^{1}(0,T;\mathbf{H}^{r_{\boldsymbol{v}}+1}(\Omega_{F}))}+h^{r_{q_{F}}}\|p_{F}\|_{H^{1}(0,T;H^{r_{q_{F}}}(\Omega_{F}))}+h^{r_{\boldsymbol{w}}}\|\boldsymbol{d}\|_{W^{2,\infty}(0,T;\mathbf{H}^{r_{\boldsymbol{w}}+1}(\Omega_{P}))}$
$\displaystyle\qquad\qquad+h^{r_{q_{P}}}\|p_{P}\|_{H^{1}(0,T;H^{r_{q_{P}}+1}(\Omega_{P}))}+h^{r_{z}}\|\varphi\|_{H^{1}(0,T;H^{r_{z}}(\Omega_{P}))}$
$\displaystyle\qquad\qquad+\Delta
t\big{(}\|\boldsymbol{d}\|_{H^{3}(0,T;\mathbf{L}^{2}(\Sigma))}+\|p_{P}\|_{H^{2}(0,T;L^{2}(\Omega_{P})}+\|\varphi\|_{H^{2}(0,T;L^{2}(\Omega_{P})}\big{)}\Big{)},$
(5.5)

with $1\leq r_{\boldsymbol{v}}\leq k_{\boldsymbol{v}}$, $0\leq r_{q_{F}}\leq
k_{q_{F}}+1$, $1\leq r_{\boldsymbol{w}}\leq k_{\boldsymbol{w}}$, $1\leq
r_{q_{P}}\leq k_{q_{P}}$, and $0\leq r_{z}\leq k_{z}+1$.

Proof. We decompose the numerical errors into approximation and discretisation
components:

$\displaystyle\boldsymbol{u}-\boldsymbol{u}_{h}=(\boldsymbol{u}-I^{\mathbf{V}}\boldsymbol{u})+(I^{\mathbf{V}}\boldsymbol{u}-\boldsymbol{u}_{h})=:e_{\boldsymbol{u},I}+e_{\boldsymbol{u},h},$
$\displaystyle
p_{F}-p_{F,h}=(p_{F}-I^{\mathrm{Q}^{F}}p_{F})+(I^{\mathrm{Q}^{F}}p_{F}-p_{F,h})=:e_{p_{F},I}+e_{p_{F},h},$
$\displaystyle\boldsymbol{d}-\boldsymbol{d}_{h}=(\boldsymbol{d}-I^{\mathbf{W}}\boldsymbol{d})+(I^{\mathbf{W}}\boldsymbol{d}-\boldsymbol{d}_{h})=:e_{\boldsymbol{d},I}+e_{\boldsymbol{d},h},$
$\displaystyle
p_{P}-p_{P,h}=(p_{P}-I^{\mathrm{Q}^{P}}p_{P})+(I^{\mathrm{Q}^{P}}p_{P}-p_{P,h})=:e_{p_{P},I}+e_{p_{P},h},$
$\displaystyle\varphi-\varphi_{h}=(\varphi-I^{\mathrm{Z}}\varphi)+(I^{\mathrm{Z}}\varphi-\varphi_{h})=:e_{\varphi,I}+e_{\varphi,h}.$

Denote the time discretisation errors as
$r_{\phi}^{n}:=\phi(t_{n})-\partial_{t}^{n}\phi$, for
$\phi\in\\{\boldsymbol{d},p_{P},\varphi\\}$. Subtracting (5.1) from (3.1), we
obtain the error system

$\displaystyle
a_{2}^{F}(e_{\boldsymbol{u},h}^{n},\boldsymbol{v}_{h})+b^{F}_{1}(\boldsymbol{v}_{h},e_{p_{F},h}^{n})+b^{\Sigma}_{2}(\boldsymbol{v}_{h},e_{p_{P},h}^{n})+b_{3}^{\Sigma}(\boldsymbol{v}_{h},\partial_{t}^{n}e_{\boldsymbol{d},h})$
$\displaystyle\quad=-a_{2}^{F}(e_{\boldsymbol{u},I}^{n},\boldsymbol{v}_{h})-b^{F}_{1}(\boldsymbol{v}_{h},e_{p_{F},I}^{n})-b^{\Sigma}_{2}(\boldsymbol{v}_{h},e_{p_{P},I}^{n})-b_{3}^{\Sigma}(\boldsymbol{v}_{h},\partial_{t}^{n}e_{\boldsymbol{d},I})-b_{3}^{\Sigma}(\boldsymbol{v}_{h},r_{\boldsymbol{d}}^{n}),$
(5.6a)
$\displaystyle-b^{F}_{1}(e_{\boldsymbol{u},h}^{n},q_{F,h})=b^{F}_{1}(e_{\boldsymbol{u},I}^{n},q_{F,h}),$
(5.6b) $\displaystyle
b_{3}^{\Sigma}(e_{\boldsymbol{u},h}^{n},\boldsymbol{w}_{h})+b_{4}^{\Sigma}(\boldsymbol{w}_{h},e_{p_{P},h}^{n})+a_{1}^{P}(e_{\boldsymbol{d},h}^{n},\boldsymbol{w}_{h})+a_{2}^{\Sigma}(\partial_{t}^{n}e_{\boldsymbol{d},h},\boldsymbol{w}_{h})+b_{1}^{P}(\boldsymbol{w}_{h},e_{\varphi,h}^{n})=-b_{3}^{\Sigma}(e_{\boldsymbol{u},I}^{n},\boldsymbol{w}_{h})$
$\displaystyle\quad-
b_{4}^{\Sigma}(\boldsymbol{w}_{h},e_{p_{P},I}^{n})-a_{1}^{P}(e_{\boldsymbol{d},I}^{n},\boldsymbol{w}_{h})-a_{2}^{\Sigma}(\partial_{t}^{n}e_{\boldsymbol{d},I},\boldsymbol{w}_{h})-b_{1}^{P}(\boldsymbol{w}_{h},e_{\varphi,I}^{n})-a_{2}^{\Sigma}(r_{\boldsymbol{d}}^{n},\boldsymbol{w}_{h}),$
(5.6c) $\displaystyle-
b_{2}^{\Sigma}(e_{\boldsymbol{u},h}^{n},q_{P,h})-b_{4}^{\Sigma}(\partial_{t}^{n}e_{\boldsymbol{d},h},q_{P,h})+a_{3}^{P}(\partial_{t}^{n}e_{p_{P},h},q_{P,h})+a_{4}^{P}(e_{p_{P},h}^{n},q_{P,h})-b_{2}^{P}(\partial_{t}^{n}e_{\varphi,h},q_{P,h})$
$\displaystyle\quad=b_{2}^{\Sigma}(e_{\boldsymbol{u},I}^{n},q_{P,h})+b_{4}^{\Sigma}(\partial_{t}^{n}e_{\boldsymbol{d},I},q_{P,h})-a_{3}^{P}(\partial_{t}^{n}e_{p_{P},I},q_{P,h})-a_{4}^{P}(e_{p_{P},I}^{n},q_{P,h})$
$\displaystyle\quad+b_{2}^{P}(\partial_{t}^{n}e_{\varphi,I},q_{P,h})-a_{3}^{P}(r_{p_{P}}^{n},q_{P,h})-b_{2}^{P}(r_{\varphi}^{n},q_{P,h}),$
(5.6d) $\displaystyle-
b_{1}^{P}(\partial_{t}^{n}e_{\boldsymbol{d},h},\psi_{h})-b_{2}^{P}(\psi_{h},\partial_{t}^{n}e_{p_{P},h})+a_{5}^{P}(\partial_{t}^{n}e_{\varphi,h},\psi_{h})=b_{1}^{P}(\partial_{t}^{n}e_{\boldsymbol{d},I},\psi_{h})+b_{2}^{P}(\psi_{h},\partial_{t}^{n}e_{p_{P},I}).$
(5.6e)

Equation (5.6e) has been obtained by taking the divided difference of the
error equation at $t_{n}$ and $t_{n-1}$ for $n=1,\ldots,N$, using that it is
satisfied by the initial data. We also used the orthogonality property of
$I^{\mathrm{Z}}$ to conclude that
$a_{5}^{P}(\partial_{t}^{n}e_{\varphi,I},\psi_{h})=0$. Now, taking
$(\boldsymbol{v}_{h},q_{F,h},\boldsymbol{w}_{h},q_{P,h},\psi_{h})=(e_{\boldsymbol{u},h}^{n},e_{p_{F},h}^{n},\partial_{t}^{n}e_{\boldsymbol{d},h},e_{p_{P},h}^{n},e_{\varphi,h}^{n})$
in (5.6), summing the equations, and using the identity

$(\partial_{t}^{n}\phi)\phi^{n}=\frac{1}{2}\partial_{t}^{n}(\phi^{2})+\frac{1}{2}\Delta
t(\partial_{t}^{n}\phi)^{2},$

results in

$\begin{split}&2\mu_{f}\\!\int_{\Omega_{F}}\\!\\!\boldsymbol{\epsilon}(e_{\boldsymbol{u},h}^{n}):\boldsymbol{\epsilon}(e_{\boldsymbol{u},h}^{n})+\sum_{j=1}^{d-1}\langle\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(e_{\boldsymbol{u},h}^{n}-\partial_{t}^{n}e_{\boldsymbol{d},h})\cdot\boldsymbol{t}^{j}\,,(e_{\boldsymbol{u},h}^{n}-\partial_{t}^{n}e_{\boldsymbol{d},h})\cdot\boldsymbol{t}^{j}\rangle_{\Sigma}\\\
&\qquad+\frac{1}{2}\partial_{t}^{n}2\mu_{s}\int_{\Omega_{P}}\boldsymbol{\epsilon}(e_{\boldsymbol{d},h}):\boldsymbol{\epsilon}(e_{\boldsymbol{d},h})+\frac{1}{2}\Delta
t\,2\mu_{s}\int_{\Omega_{P}}\partial_{t}^{n}\boldsymbol{\epsilon}(e_{\boldsymbol{d},h}):\partial_{t}^{n}\boldsymbol{\epsilon}(e_{\boldsymbol{d},h})\\\
&\qquad+\frac{1}{2}\partial_{t}^{n}C_{0}\int_{\Omega_{P}}(e_{p_{P},h})^{2}+\frac{1}{2}\Delta
tC_{0}\int_{\Omega_{P}}(\partial_{t}^{n}e_{p_{P},h})^{2}+\int_{\Omega_{P}}\frac{\kappa}{\mu_{f}}\nabla
e_{p_{P},h}^{n}\cdot\nabla e_{p_{P},h}^{n}\\\
&\qquad+\frac{1}{2}\partial_{t}^{n}\frac{1}{\lambda}\int_{\Omega_{P}}(\alpha\,e_{p_{P},h}-e_{\varphi,h})^{2}+\frac{1}{2}\Delta
t\frac{1}{\lambda}\int_{\Omega_{P}}\left(\partial_{t}^{n}(\alpha\,e_{p_{P},h}-e_{\varphi,h})\right)^{2}=\mathcal{L}^{n},\end{split}$
(5.7)

where $\mathcal{L}^{n}$ is the collection of terms that appear on the right-
hand sides in (5.6). Using the continuity of the bilinear forms (cf. Lemma
4.1) in combination with Young’s inequality, we have

$\begin{split}\mathcal{L}^{n}&\leq\epsilon\big{(}\|e_{\boldsymbol{u},h}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e_{p_{F},h}^{n}\|_{L^{2}(\Omega_{F})}^{2}+\|e_{p_{P},h}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\|e_{\varphi,h}^{n}\|_{L^{2}(\Omega_{P})}^{2}\big{)}\\\
&\quad+C_{\epsilon}\big{(}\|e_{\boldsymbol{u},I}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e_{p_{F},I}^{n}\|_{L^{2}(\Omega_{F})}^{2}+\|e_{p_{P},I}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{p_{P},I}\|_{H^{1}(\Omega_{P})}^{2}\\\
&\quad+\|\partial_{t}^{n}e_{\varphi,I}\|_{L^{2}(\Omega_{P})}^{2}+\|r_{\boldsymbol{d}}^{n}\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\|r_{p_{P}}^{n}\|_{L^{2}(\Omega_{P})}^{2}+\|r_{\varphi}^{n}\|_{L^{2}(\Omega_{P})}^{2}\big{)}+\mathcal{H}(\partial_{t}^{n}e_{\boldsymbol{d},h}),\end{split}$
(5.8)

where $\mathcal{H}(\partial_{t}^{n}e_{\boldsymbol{d},h})$ is the collection of
terms on the right-hand side of (5.6c) with
$w_{h}=\partial_{t}^{n}e_{\boldsymbol{d},h}$. Multiplying (LABEL:discr-
error-1) by $\Delta t$, summing for $n$ from 1 to $k\in\\{1,\ldots,N\\}$, and
using (5.8), we obtain

$\displaystyle\Delta
t\sum_{n=1}^{k}\Big{(}\|e_{\boldsymbol{u},h}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\|(e_{\boldsymbol{u},h}^{n}-\partial_{t}^{n}e_{\boldsymbol{d},h})\cdot\boldsymbol{t}^{j}\|_{L^{2}(\Sigma)}^{2}+\|e_{p_{P},h}^{n}\|_{H^{1}(\Omega_{P})}^{2}\Big{)}$
$\displaystyle\quad+\|e_{\boldsymbol{d},h}^{k}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|e_{p_{P},h}^{k}\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{\lambda}\|\alpha
e_{p_{P},h}^{k}-e_{\varphi,h}^{k}\|_{L^{2}(\Omega_{P})}^{2}$
$\displaystyle\quad\leq C\Delta
t\sum_{n=1}^{k}\Big{(}\epsilon\big{(}\|e_{\boldsymbol{u},h}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e_{p_{F},h}^{n}\|_{L^{2}(\Omega_{F})}^{2}+\|e_{p_{P},h}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\|e_{\varphi,h}^{n}\|_{L^{2}(\Omega_{P})}^{2}\big{)}$
$\displaystyle\qquad+C_{\epsilon}\big{(}\|e_{\boldsymbol{u},I}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e_{p_{F},I}^{n}\|_{L^{2}(\Omega_{F})}^{2}+\|e_{p_{P},I}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{p_{P},I}\|_{H^{1}(\Omega_{P})}^{2}$
$\displaystyle\qquad+\|\partial_{t}^{n}e_{\varphi,I}\|_{L^{2}(\Omega_{P})}^{2}+\|r_{\boldsymbol{d}}^{n}\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\|r_{p_{P}}^{n}\|_{L^{2}(\Omega_{P})}^{2}+\|r_{\varphi}^{n}\|_{L^{2}(\Omega_{P})}^{2}\big{)}+\mathcal{H}(\partial_{t}^{n}e_{\boldsymbol{d},h})\Big{)}$
$\displaystyle\qquad+C\big{(}\|e_{\boldsymbol{d},h}^{0}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|e_{p_{P},h}^{0}\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{\lambda}\|\alpha
e_{p_{P},h}^{0}-e_{\varphi,h}^{0}\|_{L^{2}(\Omega_{P})}^{2}\big{)},$ (5.9)

where we also used Korn’s inequality (4.5), the Poincaré inequality (4.4), and
the permeability bound (4.7). Next, for each term in
$\mathcal{H}(\partial_{t}^{n}e_{\boldsymbol{d},h})$ we use summation by parts:

$\Delta
t\sum_{n=1}^{k}\phi^{n}\partial_{t}^{n}e_{\boldsymbol{d},h}=\phi^{k}e_{\boldsymbol{d},h}^{k}-\phi^{1}e_{\boldsymbol{d},h}^{0}-\Delta
t\sum_{n=2}^{k}\partial_{t}^{n}\phi\,e_{\boldsymbol{d},h}^{n-1},$

where $\phi$ stands for any of the functions that appear in
$\mathcal{H}(\partial_{t}^{n}e_{\boldsymbol{d},h})$. Then, for the first term
in $\mathcal{H}(\partial_{t}^{n}e_{\boldsymbol{d},h})$ we write, using Young’s
inequality,

$\begin{split}&\Delta
t\sum_{n=1}^{k}b_{3}^{\Sigma}(e_{\boldsymbol{u},I}^{n},\partial_{t}^{n}e_{\boldsymbol{d},h})=b_{3}^{\Sigma}(e_{\boldsymbol{u},I}^{k},e_{\boldsymbol{d},h}^{k})-b_{3}^{\Sigma}(e_{\boldsymbol{u},I}^{1},e_{\boldsymbol{d},h}^{0})-\Delta
t\sum_{n=2}^{k}b_{3}^{\Sigma}(\partial_{t}^{n}e_{\boldsymbol{u},I},e_{\boldsymbol{d},h}^{n-1})\\\
&\quad\leq\epsilon\|e_{\boldsymbol{d},h}^{k}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+C_{\epsilon}\|e_{\boldsymbol{u},I}^{k}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+C\big{(}\|e_{\boldsymbol{d},h}^{0}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|e_{\boldsymbol{u},I}^{1}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}\big{)}\\\
&\qquad+C\Delta
t\sum_{n=2}^{k}\big{(}\|e_{\boldsymbol{d},h}^{n-1}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\boldsymbol{u},I}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}\big{)}.\end{split}$
(5.10)

Now, combining (5.1)–(LABEL:discr-error-4), and using for the rest of the
terms in $\mathcal{H}(\partial_{t}^{n}e_{\boldsymbol{d},h})$ bounds that are
similar to (LABEL:discr-error-4), we arrive at

$\displaystyle\Delta
t\sum_{n=1}^{k}\Big{(}\|e_{\boldsymbol{u},h}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\|(e_{\boldsymbol{u},h}^{n}-\partial_{t}^{n}e_{\boldsymbol{d},h})\cdot\boldsymbol{t}^{j}\|_{L^{2}(\Sigma)}^{2}+\|e_{p_{P},h}^{n}\|_{H^{1}(\Omega_{P})}^{2}\Big{)}$
$\displaystyle\quad+\|e_{\boldsymbol{d},h}^{k}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|e_{p_{P},h}^{k}\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{\lambda}\|\alpha
e_{p_{P},h}^{k}-e_{\varphi,h}^{k}\|_{L^{2}(\Omega_{P})}^{2}$
$\displaystyle\quad\leq C\Delta
t\sum_{n=1}^{k}\Big{(}\epsilon\big{(}\|e_{\boldsymbol{u},h}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e_{p_{F},h}^{n}\|_{L^{2}(\Omega_{F})}^{2}+\|e_{p_{P},h}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\|e_{\varphi,h}^{n}\|_{L^{2}(\Omega_{P})}^{2}\big{)}$
$\displaystyle\qquad+\|e_{\boldsymbol{u},I}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e_{p_{F},I}^{n}\|_{L^{2}(\Omega_{F})}^{2}+\|e_{p_{P},I}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\|r_{\boldsymbol{d}}^{n}\|_{\mathbf{L}^{2}(\Sigma)}^{2}+\|r_{p_{P}}^{n}\|_{L^{2}(\Omega_{P})}^{2}+\|r_{\varphi}^{n}\|_{L^{2}(\Omega_{P})}^{2}$
$\displaystyle\qquad+\|\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{p_{P},I}\|_{H^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\varphi,I}\|_{L^{2}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\boldsymbol{u},I}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}\Big{)}$
$\displaystyle\qquad+C\Big{(}\Delta
t\sum_{n=2}^{k}\big{(}\|e_{\boldsymbol{d},h}^{n-1}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}r_{\boldsymbol{d}}\|_{\mathbf{L}^{2}(\Sigma)}^{2}\big{)}$
$\displaystyle\qquad\quad+\|e_{\boldsymbol{d},h}^{0}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|e_{p_{P},h}^{0}\|_{L^{2}(\Omega_{P})}^{2}+\frac{1}{\lambda}\|\alpha
e_{p_{P},h}^{0}-e_{\varphi,h}^{0}\|_{L^{2}(\Omega_{P})}^{2}+\|\partial_{t}^{k}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{1}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}$
$\displaystyle\qquad\quad+\|e^{k}_{\boldsymbol{u},I}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e^{k}_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|e^{k}_{p_{P},I}\|_{H^{1}(\Omega_{P})}^{2}+\|e^{k}_{\varphi,I}\|_{L^{2}(\Omega_{P})}^{2}+\|r^{k}_{\boldsymbol{d}}\|_{\mathbf{L}^{2}(\Sigma)}^{2}$
$\displaystyle\qquad\quad+\|e^{1}_{\boldsymbol{u},I}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e^{1}_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|e^{1}_{p_{P},I}\|_{H^{1}(\Omega_{P})}^{2}+\|e^{1}_{\varphi,I}\|_{L^{2}(\Omega_{P})}^{2}+\|r^{1}_{\boldsymbol{d}}\|_{\mathbf{L}^{2}(\Sigma)}^{2}\Big{)}.$
(5.11)

We continue with bounding $\|e_{p_{F},h}^{n}\|_{L^{2}(\Omega_{F})}^{2}$ and
$\|e_{\varphi,h}^{n}\|_{L^{2}(\Omega_{P})}^{2}$, which appear on the right-
hand side above. The inf-sup condition (4.1a) and (5.6a) imply

$\displaystyle\Delta
t\sum_{n=1}^{k}\|e_{p_{F},h}^{n}\|_{L^{2}(\Omega_{F})}^{2}\leq C\Delta
t\sum_{n=1}^{k}\Big{(}\|e_{\boldsymbol{u},h}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\sum_{j=1}^{d-1}\|(e_{\boldsymbol{u},h}^{n}-\partial_{t}^{n}e_{\boldsymbol{d},h})\cdot\boldsymbol{t}^{j}\|_{L^{2}(\Sigma)}^{2}+\|e_{p_{P},h}^{n}\|_{H^{1}(\Omega_{P})}^{2}$
$\displaystyle\quad+\|e_{\boldsymbol{u},I}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e_{p_{F},I}^{n}\|_{L^{2}(\Omega_{F})}^{2}+\|e_{p_{P},I}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|r_{\boldsymbol{d}}^{n}\|_{\mathbf{L}^{2}(\Sigma)}^{2}\Big{)}.$

On the other hand, the inf-sup condition (4.1b) and (5.6c) allow us to get

$\begin{split}&\Delta
t\sum_{n=1}^{k}\|e_{\varphi,h}^{n}\|_{L^{2}(\Omega_{P})}^{2}\leq C\Delta
t\sum_{n=1}^{k}\Big{(}\|e_{p_{P},h}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\sum_{j=1}^{d-1}\|(e_{\boldsymbol{u},h}^{n}-\partial_{t}^{n}e_{\boldsymbol{d},h})\cdot\boldsymbol{t}^{j}\|_{L^{2}(\Sigma)}^{2}+\|e_{\boldsymbol{d},h}^{n}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}\\\
&\
+\|e_{\boldsymbol{u},I}^{n}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}+\|e_{p_{P},I}^{n}\|_{H^{1}(\Omega_{P})}^{2}+\|e_{\boldsymbol{d},I}^{n}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|e_{\varphi,I}^{n}\|_{L^{2}(\Sigma)}^{2}+\|r_{\boldsymbol{d}}^{n}\|_{\mathbf{L}^{2}(\Sigma)}^{2}\Big{)}.\end{split}$
(5.12)

Before combining (5.1)–(LABEL:phi-discr), we note that the terms involving
$\partial_{t}^{n}$ on the right-hand sides require special treatment. In
particular, it holds for $\phi(t)$ that

$(\partial_{t}^{n}\phi)^{2}=\frac{1}{\Delta
t^{2}}\Big{(}\int_{t_{n-1}}^{t_{n}}\partial_{t}\phi\Big{)}^{2}\leq\frac{1}{\Delta
t}\int_{t_{n-1}}^{t_{n}}(\partial_{t}\phi)^{2},$

implying, for $\phi:[0,T]\to V$, where $V$ is a Banach space with norm
$\|\cdot\|_{V}$, that

$\Delta
t\sum_{n=1}^{k}\|\partial_{t}^{n}\phi\|_{V}^{2}\leq\int_{0}^{t_{k}}\|\partial_{t}\phi\|_{V}^{2}.$

We then have

$\begin{split}&\Delta
t\sum_{n=1}^{k}\big{(}\|\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{p_{P},I}\|_{H^{1}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\varphi,I}\|_{L^{2}(\Omega_{P})}^{2}+\|\partial_{t}^{n}e_{\boldsymbol{u},I}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}\big{)}\\\
&\quad\leq
C\int_{0}^{t_{k}}\big{(}\|\partial_{t}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}+\|\partial_{t}e_{p_{P},I}\|_{H^{1}(\Omega_{P})}^{2}+\|\partial_{t}e_{\varphi,I}\|_{L^{2}(\Omega_{P})}^{2}+\|\partial_{t}e_{\boldsymbol{u},I}\|_{\mathbf{H}^{1}(\Omega_{F})}^{2}\big{)}.\end{split}$
(5.13)

To bound the term
$\|\partial_{t}^{n}\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}$
in (5.1), for any $\phi(t)$ we have, using the integral mean value theorem and
the mean value theorem,

$(\partial_{t}^{n}\partial_{t}^{n}\phi)^{2}=\frac{1}{\Delta
t^{4}}\Big{(}\int_{t_{n-1}}^{t^{n}}\partial_{t}\phi-\int_{t_{n-2}}^{t^{n-1}}\partial_{t}\phi\Big{)}^{2}=\frac{1}{\Delta
t^{2}}\big{(}\partial_{t}\phi(\xi^{n})-\partial_{t}\phi(\xi^{n-1})\big{)}^{2}=\partial_{tt}\phi(\xi),\quad\xi\in[t_{n-2},t_{n}].$

Therefore it holds that

$\Delta
t\sum_{n=2}^{k}\|\partial_{t}^{n}\partial_{t}^{n}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}\leq
C\mathop{\mathrm{ess\,sup}}\limits_{t\in(0,t_{k})}\|\partial_{tt}e_{\boldsymbol{d},I}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}.$

Next, we need to bound the time discretisation error. Taylor’s expansion gives

$r_{\phi}^{n}=\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n-1}}\partial_{tt}\phi(t)(t_{n-1}-t)\,\mathrm{d}t,$

thus, for $\phi:[0,T]\to V$,

$\|r_{\phi}^{n}\|_{V}\leq C\Delta
t\mathop{\mathrm{ess\,sup}}\limits_{t\in(t_{n-1},t_{n})}\|\partial_{tt}\phi\|_{V}\quad\text{and}\quad\Delta
t\sum_{n=1}^{k}\|r_{\phi}^{n}\|_{V}^{2}\leq C\Delta
t^{2}\int_{0}^{t_{k}}\|\partial_{tt}\phi\|_{V}^{2}.$

Similarly,

$\Delta t\sum_{n=1}^{k}\|\partial_{t}^{n}r_{\phi}\|_{V}^{2}\leq C\Delta
t^{2}\int_{0}^{t_{k}}\|\partial_{ttt}\phi\|_{V}^{2}.$

Finally, we need a bound on the initial discretisation error. Recalling the
construction of the discrete initial data in the proof of Theorem 4.1 and the
definition of the continuous initial data in the proof of Theorem 4.2, we note
that $(\boldsymbol{u}_{h}^{0},p_{F,h}^{0})$ is the Stokes elliptic projection
of $(\boldsymbol{u}(0),p_{F,h}(0))$ based on (4.2a)–(4.2b) at $t=0$ with
$\partial_{t}\boldsymbol{d}_{h}=0$ and the term
$b_{2}^{\Sigma}(p_{P,0}-p_{P,h,0},\boldsymbol{v}_{h})$ on the right-hand side.
In addition, $(\boldsymbol{d}_{h}^{0},\varphi_{h}^{0})$ is the elliptic
projection of $(\boldsymbol{d}(0),\varphi(0))$ based on the stable problem
(4.2c)–(4.2e) at $t=0$ with $\partial_{t}\boldsymbol{d}_{h}=0$ and the terms
$b_{3}^{\Sigma}(\boldsymbol{u}(0)-\boldsymbol{u}_{h}^{0},\boldsymbol{w}_{h})$,
$b_{4}^{\Sigma}(\boldsymbol{w}_{h},p_{P,0}-p_{P,h,0})$, and
$b_{2}^{P}(\psi_{h},p_{P,0}-p_{P,h,0})$ on the right-hand side. Classical
finite element analysis for these two problems implies that

$\begin{split}&\|e_{\boldsymbol{u},h}^{0}\|_{\mathbf{H}^{1}(\Omega_{F})}+\|e_{p_{F},h}^{0}\|_{L^{2}(\Omega_{F})}+\|e_{\boldsymbol{d},h}^{0}\|_{\mathbf{H}^{1}(\Omega_{P})}+\|e_{p_{P},h}^{0}\|_{L^{2}(\Omega_{P})}+\|e_{\varphi,h}^{0}\|_{L^{2}(\Omega_{P})}\\\
&\quad\leq
C\big{(}\|e_{\boldsymbol{u},I}^{0}\|_{\mathbf{H}^{1}(\Omega_{F})}+\|e_{p_{F},I}^{0}\|_{L^{2}(\Omega_{F})}+\|e_{\boldsymbol{d},I}^{0}\|_{\mathbf{H}^{1}(\Omega_{P})}+\|e_{p_{P},I}^{0}\|_{L^{2}(\Omega_{P})}+\|e_{\varphi,I}^{0}\|_{L^{2}(\Omega_{P})}\big{)}.\end{split}$
(5.14)

The assertion of the theorem follows from combining (5.1)–(LABEL:init-error),
using the discrete Gronwall inequality [47, Lemma 1.4.2] for the term $\Delta
t\sum_{n=2}^{k}\|e_{\boldsymbol{d},h}^{n-1}\|_{\mathbf{H}^{1}(\Omega_{P})}^{2}$,
and applying the triangle inequality and the approximation properties
(5.4a)–(5.4e). $\square$

##  6 Representative computational results

All routines have been implemented using the open source finite element
library FEniCS [4], as well as the specialised module multiphenics [10] for
handling subdomain- and boundary- restricted terms that we require to impose
transmission conditions across interfaces. The solvers are monolithic and the
solution of all linear systems is performed with the distributed direct solver
MUMPS. We present four examples: convergence tests (example 1), channel flow
behaviour (example 2), a simulation of subsurface flow with heterogeneous
random permeability (example 3), and the solution of an axisymmetric problem
using parameters relevant to eye poromechanics (example 4). For examples 1,2
and 4 we use Taylor-Hood elements for the pairs velocity-pressure and
displacement - total pressure, plus continuous and piecewise quadratic
elements for Biot fluid pressure. For example 3, the inf-sup stable pair used
is the MINI element, and the Biot fluid pressure is approximated with
continuous and piecewise linear elements.

###  6.1 Convergence tests against manufactured solutions

The accuracy of the spatio-temporal discretisation is verified using the
following closed-form solutions defined on the domains
$\Omega_{F}=(-1,1)\times(0,2)$, $\Omega_{P}=(-1,1)\times(-2,0)$, separated by
the interface $\Sigma=(-1,1)\times\\{0\\}$

$\begin{split}\boldsymbol{u}=\sin(t)\begin{pmatrix}-\cos(\pi x)\sin(\pi y)\\\
\sin(\pi x)\cos(\pi y)\end{pmatrix},\quad p_{F}=\sin(t)\cos(\pi x)\cos(\pi
y),\\\ \boldsymbol{d}=\cos(t)\mathop{\mathrm{curl}}\nolimits(\sin(\pi
xy)),\quad p_{P}=\cos(t)\sin(\pi x)\sin(\pi y),\quad\varphi=\alpha
p_{P}-\lambda\mathop{\mathrm{div}}\nolimits\boldsymbol{d}.\end{split}$ (6.1)

DoF | $h$ | $\texttt{e}_{\boldsymbol{u}}$ | rate | $\texttt{e}_{p_{F}}$ | rate | $\texttt{e}_{\boldsymbol{d}}$ | rate | $\texttt{e}_{p_{P}}$ | rate | $\texttt{e}_{\varphi}$ | rate
---|---|---|---|---|---|---|---|---|---|---|---
144 | 1.414 | 4.70604 | – | 0.82152 | – | 21.4632 | – | 2.80468 | – | 19.3402 | –
456 | 0.7071 | 1.66701 | 1.497 | 0.29604 | 1.472 | 9.75813 | 1.137 | 0.84917 | 1.721 | 8.20939 | 1.236
1608 | 0.3536 | 0.40411 | 2.044 | 0.06081 | 2.260 | 2.78206 | 1.811 | 0.24942 | 1.769 | 1.21860 | 2.422
6024 | 0.1768 | 0.09746 | 2.051 | 0.01376 | 2.138 | 0.72730 | 1.936 | 0.06567 | 1.943 | 0.25337 | 2.240
23304 | 0.0884 | 0.02405 | 2.014 | 0.00328 | 2.052 | 0.18545 | 1.984 | 0.01661 | 1.985 | 0.05166 | 2.231
91656 | 0.0442 | 0.00597 | 2.005 | 0.00081 | 2.018 | 0.04683 | 1.991 | 0.00410 | 1.991 | 0.00472 | 2.098
363528 | 0.0221 | 0.00151 | 2.001 | 0.00022 | 2.008 | 0.01172 | 1.998 | 0.00100 | 1.995 | 0.00061 | 2.015
Table 6.1: Example 1. Experimental errors associated with the spatial
discretisation and convergence rates for the approximate solutions
$\boldsymbol{u}_{h}$, $p_{F,h}$, $\boldsymbol{d}_{h}$, $p_{P,h}$, and
$\varphi_{h}$ using
$\mathbb{P}^{2}_{2}-\mathbb{P}_{1}-\mathbb{P}^{2}_{2}-\mathbb{P}_{2}-\mathbb{P}_{1}$
elements and computed at the last time step.

We consider $\Gamma_{P}^{p}$ to be the bottom horizontal segment,
$\Gamma_{P}^{\boldsymbol{d}}$ to be the lower vertical walls,
$\Gamma_{F}^{\boldsymbol{u}}$ to be the top horizontal segment, and
$\Gamma_{F}^{\boldsymbol{\sigma}}$ to be the upper vertical walls. the
synthetic model parameters are taken as

$\displaystyle\lambda=1000,\quad\mu_{s}=1,\quad\mu_{f}=0.1,\quad\alpha=\gamma=1,\quad
c_{0}=0.01,\quad\kappa=0.001,\quad\rho_{s}=1.2,\quad\rho_{f}=1,$

all regarded adimensional and do not have physical relevance in this case, as
we will be simply testing the convergence of the finite element
approximations. The manufactured solutions (6.1) are used to prescribe initial
conditions, essential non-homogeneous velocity and displacement boundary
conditions, as well as natural non-homogeneous flux conditions for fluid
pressure. These functions do not necessarily fulfil the interface conditions,
so additional terms are required giving modified relations on $\Sigma$:

$\displaystyle\boldsymbol{u}\cdot\boldsymbol{n}-(\partial_{t}\boldsymbol{d}-\frac{\kappa}{\mu_{f}}\nabla
p_{P})\cdot\boldsymbol{n}$
$\displaystyle=m^{1}_{\Sigma,\text{ex}},\qquad(2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I})\boldsymbol{n}-(2\mu_{s}\boldsymbol{\epsilon}(\boldsymbol{d})-\varphi\mathbf{I})\boldsymbol{n}=\boldsymbol{m}^{2}_{\Sigma,\text{ex}},$
$\displaystyle\boldsymbol{n}\cdot(2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I})\boldsymbol{n}+{p_{P}}$
$\displaystyle=m^{3}_{\Sigma,\text{ex}},\qquad\boldsymbol{n}\cdot(2\mu_{f}\boldsymbol{\epsilon}(\boldsymbol{u})-p_{F}\mathbf{I})\boldsymbol{t}+\frac{\gamma\mu_{f}}{\sqrt{\kappa}}(\boldsymbol{u}-\partial_{t}\boldsymbol{d})\cdot\boldsymbol{t}=m^{4}_{\Sigma,\text{ex}},$

and the additional scalar and vector terms $m^{i}_{\Sigma,\text{ex}}$
(computed with the exact solutions (6.1)) entail the following changes in the
linear functionals

$\displaystyle
F^{F}(\boldsymbol{v})=\rho_{f}\int_{\Omega_{F}}\boldsymbol{g}\cdot\boldsymbol{v}+\langle
m^{3}_{\Sigma,\text{ex}},\boldsymbol{v}\cdot\boldsymbol{n}\rangle_{\Sigma}+\langle
m^{4}_{\Sigma,\text{ex}},\boldsymbol{v}\cdot\boldsymbol{t}\rangle_{\Sigma},$
$\displaystyle
F^{P}(\boldsymbol{w})=\rho_{s}\int_{\Omega_{P}}\boldsymbol{f}\cdot\boldsymbol{w}+\int_{\Omega_{P}}\boldsymbol{m}^{2}_{\Sigma,\text{ex}}\cdot\boldsymbol{w}+\langle
m^{3}_{\Sigma,\text{ex}},\boldsymbol{w}\cdot\boldsymbol{n}\rangle_{\Sigma}+\langle
m^{4}_{\Sigma,\text{ex}},\boldsymbol{w}\cdot\boldsymbol{t}\rangle_{\Sigma},$
$\displaystyle G(q_{P})=\int_{\Omega_{P}}\rho_{f}\boldsymbol{g}\cdot\nabla
q_{P}-\langle\rho_{f}\boldsymbol{g}\cdot\boldsymbol{n},q_{P}\rangle_{\Sigma}-\langle
m^{1}_{\Sigma,\text{ex}},q_{P}\rangle_{\Sigma}.$

$\Delta t$ | $\hat{\texttt{e}}_{\boldsymbol{u}}$ | rate | $\hat{\texttt{e}}_{p_{F}}$ | rate | $\hat{\texttt{e}}_{\boldsymbol{d}}$ | rate | $\hat{\texttt{e}}_{p_{P}}$ | rate | $\hat{\texttt{e}}_{\varphi}$ | rate
---|---|---|---|---|---|---|---|---|---|---
0.5 | 10.549 | – | 2.5844 | – | 43.764 | – | 8.6738 | – | 38.734 | –
0.25 | 5.1408 | 0.984 | 1.2710 | 1.163 | 21.840 | 1.211 | 4.3673 | 1.125 | 19.371 | 1.055
0.125 | 2.2689 | 1.205 | 0.6485 | 1.182 | 10.661 | 1.204 | 2.1690 | 1.032 | 9.6901 | 0.962
0.0625 | 1.1365 | 1.107 | 0.3216 | 1.095 | 5.4517 | 1.121 | 1.0795 | 1.025 | 4.8614 | 0.981
0.03125 | 0.6813 | 1.004 | 0.1615 | 1.071 | 2.7326 | 1.043 | 0.5527 | 0.987 | 2.4396 | 0.992
Table 6.2: Example 1. Experimental cumulative errors associated with the
temporal discretisation and convergence rates for the approximate solutions
$\boldsymbol{u}_{h}$, $p_{F,h}$, $\boldsymbol{d}_{h}$, $p_{P,h}$, and
$\varphi_{h}$ using a backward Euler scheme.

We generate successively refined simplicial grids and use a sufficiently small
(non dimensional) time step $\Delta t=0.01$ and simulate a relatively short
time horizon $t_{\mathrm{final}}=3\Delta t$, to guarantee that the error
produced by the time discretisation does not dominate. Errors between the
approximate and exact solutions are tabulated against the number of degrees of
freedom in Table 6.1. This error history confirms the optimal convergence of
the finite element scheme (in this case, second-order) for all variables in
their respective norms, where a slightly better rate is seen for the total
pressure. In the table caption, $\mathbb{P}_{k}$ denotes the space of
piecewise polynomial functions being of total degree up to $k$.

For this example we have reincorporated the acceleration and the nonlinear
convective term, and the Newton-Raphson algorithm takes, in average, three
iterations to reach the prescribed tolerance of $10^{-8}$ on the residuals.

The convergence in time achieved by the backward Euler method is verified by
partitioning the time interval $(0,1)$ into successively refined uniform
discretisations and computing accumulated errors

$\hat{\texttt{e}}_{s}=\biggl{(}\sum_{n=1}^{N}\Delta
t\|s(t_{n+1})-s_{h}^{n+1}\|^{2}_{\star}\biggr{)}^{1/2},$

where $\|\cdot\|_{\star}$ denotes the appropriate space norm for the generic
vector or scalar field $s$. For this test we use a fixed mesh involving 10K
DoFs. The results are shown in Table 6.2, confirming the optimal first-order
convergence.

###  6.2 Channel filtration and stress build-up on interface deformation

Figure 6.1: Example 2. Infiltration into a deformable porous medium. Velocity
in the fluid domain, fluid pressure, extended displacement of fluid and solid
domains, solid displacement, fluid pressure in the poroelastic domain, and
total pressure. All snapshots are taken at time $t=2$, and the black outer
line indicates the location of the undeformed domain.

Although the model stated in Section 2 holds in the limit of small strains, it
is possible to have large displacements, likely located near the interface
(and without violating the model assumptions). In this scenario, the
discretisation might no longer be suitable. A simple remedy consists in
smoothly moving the fluid domain and the fluid mesh to avoid distortions
generated near the interface. We use a standard harmonic extension (see e.g.
[16]) that is solved at each time step, just after (5.2): Find
$\boldsymbol{d}_{h}^{*}=\boldsymbol{d}_{h}+\widehat{\boldsymbol{d}}_{h}$ such
that

$-D\Delta\widehat{\boldsymbol{d}}_{h}=\boldsymbol{0}\quad\text{ in
$\Omega_{F}$},\qquad\widehat{\boldsymbol{d}}_{h}=\boldsymbol{d}_{h}\text{ on
$\Sigma$},\qquad\text{ and
}\quad\widehat{\boldsymbol{d}}_{h}=\boldsymbol{0}\text{ on
$\partial\Omega_{F}.$}$ (6.2)

And then we perform an $L^{2}-$projection of both $\boldsymbol{d}_{h}$ and
$\widehat{\boldsymbol{d}}_{h}$ into $\mathbf{W}_{h}+\mathbf{V}_{h}$ and add
them to obtain the global displacement $\boldsymbol{d}_{h}^{*}$.

We illustrate the effect of using (6.2) by looking at the behaviour of normal
filtration into a 2D deformable porous medium. The same domains as in the
accuracy tests are employed here (that is, the single phase fluid domain
located on top of the poroelastic domain), however the boundary treatment is
as follows, assuming that the flow is driven by pressure differences only. On
the top segment we impose the fluid pressure
$p_{F}^{\mathrm{in}}=p_{0}\sin^{2}(\pi t)$ with $p_{0}=2$, and on the outlet
(the bottom segment) the fluid pressure $p_{P}^{\mathrm{out}}=0$. On the
vertical walls of $\Omega_{F}$ we set $\boldsymbol{u}=\boldsymbol{0}$ while on
the vertical walls of $\Omega_{P}$ we set the slip condition
$\boldsymbol{d}\cdot\boldsymbol{n}=0$ (and the porous structure is free to
deform on the outlet boundary, i.e., zero traction imposed). The permeability
is $\kappa=0.02$ and the remaining parameters are

$\displaystyle\lambda=10,\quad\mu_{s}=5,\quad\rho_{s}=1.1,\quad\rho_{f}=1,\quad\alpha=0.6,\quad\gamma=0.1,\quad
C_{0}=0.01,$

and we assume that there are no body forces nor gravity acting on the system.
In contrast with the convergence tests, for this example we use piecewise
linear and continuous finite elements for the approximation of $p_{P}$.

Figure 6.2: Example 2. Zoom of the meshes on the interface at times $t=0$,
$t=1$ and $t=2$. Effect of using or not the harmonic extension to move the
fluid domain (bottom and top, respectively).

The numerical results are presented in Figure 6.1. The effect of the interface
can be clearly seen in the top left panel where recirculation vortices replace
the parabolic profile at the inlet; also, in the poroelastic domain, we see
that, close to the interface, the solid displacement and the fluid pressure
are heterogeneous in the horizontal direction, before recovering the expected
constant value (constant in the horizontal direction) expected in the far
field. We also plot the evolution of the mesh deformation near the interface.
From Figure 6.2 one can see that for large enough interfacial displacements,
the elements close to it exhibit a large distortion.

###  6.3 Simulation of subsurface fracture flow

Next we include a test case that illustrates the applicability of the
formulation in hydraulic fracturing. The problem setup follows [6, Section
5.2.2] (except that we do not model tracer transport), considering a
rectangular domain $\Omega=(0,3.048)\,[\text{m}]\times(0,6.096)$ [m] including
a relatively large fracture regarded as a macro void, or open channel
$\Omega_{F}$ filled with an incompressible fluid (see Figure 6.3, left), and
the Biot domain is $\Omega_{P}=\Omega\setminus\Omega_{F}$. The heterogeneous
(but isotropic in the $xy$-plane) permeability $\kappa(\boldsymbol{x})$ is the
non-smooth pattern taken from the Cartesian SPE10 benchmark data / model 2
(see, e.g., [1, 21]), which we rescale as in [7] and project onto a piecewise
constant field defined on an unstructured triangular mesh for the poroelastic
geometry. There are 85 distinct layers within two general categories. We
choose layer 80 from the dataset, which corresponds to the Upper Ness region
exhibiting a fluvial fan pattern (flux channels of higher permeability and
porosity). The Lamé parameters (in [KPa]) are highly heterogeneous and
determined from the Poisson ratio $\nu=0.2$ and the Young modulus
$E(\boldsymbol{x})=10^{7}(1-2\phi(\boldsymbol{x}))^{2.1}$, where
$\phi(\boldsymbol{x})$ is the porosity field also taken from layer 80 of the
benchmark dataset (see Figure 6.3, right). It is plotted in logarithmic scale
and the contrast is of about $10^{8}$. No gravity and no external loads are
considered, and the remaining parameters are $\alpha=\gamma=1$,
$c_{0}=6.89\cdot 10^{-2}\,[\text{KPa}^{-1}]$,
$\mu_{f}=10^{-6}\,[\text{KPa}\cdot\text{s}]$.

Figure 6.3: Example 3. Schematic representation of sub-domains, location of
the free-porous interface, and configuration of sub-boundaries (left panel);
and material properties (porosity $\phi(\boldsymbol{x})$, permeability
$\kappa(\boldsymbol{x})$, and Young modulus $E(\boldsymbol{x})$) from layer 80
of the SPE10 benchmark dataset for reservoir simulations, herein projected
onto a coarse mesh for the poroelastic sub-domain.

Similarly as in [7], we set the flow initially at rest
$\boldsymbol{d}(0)=\boldsymbol{0}$ and the initial Biot pressures are
$p_{P}(0)=\varphi(0)=1000$ [KPa] (since $\alpha=1$). On the inlet boundary
$\Gamma^{F}_{\mathrm{in}}$ (the vertical segment on the Stokes boundary) we
impose the inflow velocity $\boldsymbol{u}=(10,0)^{\tt t}\,$[m/s], on the
bottom, right and top sub-boundaries $\Gamma^{P}_{1}$ of the Biot domain we
prescribe sliding conditions $\boldsymbol{d}\cdot\boldsymbol{n}=0$ [m]
together with a compatible normal-tangential stress condition and a fixed Biot
pressure $p_{P}=1000$ [KPa], and on the right sub-boundaries $\Gamma^{P}_{2}$
of the Biot domain we set stress-free conditions. The time discretisation uses
the fixed time step $\Delta t=60$ [s] and we run the simulation of injection
over a period of $T=10$ [hours]. The unstructured triangular mesh has 1629
elements for the Stokes sub-mesh and 18897 elements for the Biot domain. For
this test we have used the MINI element, consisting of continuous and
piecewise linear elements with bubble enrichment for velocity and
displacement, and continuous and piecewise linear elements for all remaining
fields. The injected fluid imposes an increase of pressure on the interface
and from there the expected channel-like progressive filtration from the
Stokes to the Biot domain is clearly observed in the Biot pressure plot in the
left panel of Figure 6.4, showing higher fluxes near the tip of the fracture.
The remaining panels show, at the final time, snapshots of Stokes velocity,
Biot displacement, and post-processed Biot filtration velocity. The
deformation of the poroelastic structure closer to the fracture is relatively
small, which is why we do not consider in this test the harmonic extension of
the fluid domain addressed in the previous example.

Figure 6.4: Example 3. Snapshots of the approximate solutions (Biot fluid
pressure, post-processed Biot fluid velocity and line integral convolution,
Biot displacement magnitude and displacement arrows on the deformed domain,
Stokes fluid pressure, and Stokes velocity magnitude with line integral
convolution) for fluid injection into a fracture porous medium using the
SPE10-based benchmark test. Results were obtained after $t=10$ [hours].

###  6.4 Application to interfacial flow in the eye

Figure 6.5: Example 4. Axisymmetric interfacial flow in the eye for a linearly
varying permeability profile. Velocity magnitude and streamlines (top), fluid
pressure (second row), displacement magnitude (third row), porous total
pressure, and fluid pressures on both domains extruded to the 3D case (bottom
centre and bottom right). All solutions are shown at times $t=4.2$ [s] (left)
and $t=5$ [s] (right column).

To finalise this section we include a problem pertaining to the application of
the interaction between aqueous humour in the anterior chamber and the
trabecular meshwork. Therein, one of the main driving questions is whether one
can observe deformation of the porous skeleton (and in particular of the
interface) that could drive a rise in intra-ocular pressure. For this test we
use the axisymmetric formulation (3.3), we include the convective term in
(2.1a) and we discard gravity. A large amount of data is available to
specialise the geometry and the mechanical properties (both fluidic and
elastic) of the eye to different animal species [18, 35, 29, 42]. In our case
(and consistent with a uniform temperature of $37^{\circ}$ and a
characteristic length of $6.8\cdot 10^{-3}\,[\text{m}]$) we impose

$\displaystyle\rho_{f}=998.7\,[\text{Kg}\cdot\text{m}^{-3}],\quad\mu_{f}=7.5\cdot
10^{-4}\,[\text{Pa}\cdot\text{s}],\quad p_{0}=0\,[\text{Pa}],\quad
E=2700\,[\text{Pa}],\quad\nu=0.47,$
$\displaystyle\lambda=14388\,[\text{Pa}],\quad\mu_{s}=918.36\,[\text{Pa}],\quad\alpha=1,\quad\rho_{s}=1102\,[\text{Kg}\cdot\text{m}^{-3}],\quad
C_{0}=0,\quad\gamma=0.1.$

As discussed in [25], an increase in the Beavers-Joseph-Saffmann friction
parameter $\gamma$ leads to higher pressure differences between the trabecular
meshwork and the anterior chamber. Note that even if the choice of $\alpha=1$,
$C_{0}=0$ indicates that both constituents (fluid and solid) are assumed
intrinsically incompressible, a Poisson ratio smaller than 0.5 implies
compressibility of the poroelastic medium (due to the possible rearrangement
of the porosity field, i.e., the fluid escaping the medium).

With reference to Figure 2.1, the length of the interface $\Sigma$ between the
trabecular meshwork and the anterior chamber is $5.7\cdot
10^{-4}\,[\text{m}]$, and the length of the separation between the trabecular
meshwork and angular aqueous plexus ($\Gamma^{\text{out}}$) is $3.4\cdot
10^{-4}\,[\text{m}]$. A parabolic profile for inlet velocity with a pulsating
magnitude $u_{\text{in}}=4.89\cdot 10^{-7}\sin^{2}(\pi
t)\,[\text{m}\cdot\text{s}^{-1}]$ (that has approximately the same frequency
as the heartbeat) is imposed on $\Gamma^{\mathrm{in}}$ (the magnitude is
obtained from the ratio $\frac{C\mu_{s}\kappa}{L\mu_{f}}$ with $C=0.2$
integrated through the thickness and the condition is imposed through a
Nitsche approach with penalty parameter equal to 1 since the relevant sub-
boundary is not aligned with the axes), no-slip conditions are prescribed
essentially on the walls, and a slip condition is considered for the fluid
velocity on the symmetry axis (also imposed essentially). On the outlet
$\Gamma^{\text{out}}$ we impose zero fluid pressure $p_{P}=p_{0}$. We simulate
the interfacial flow until $t=5$ [s] and use a time step of $\Delta t=0.1$
[s]. In this case, Taylor-Hood elements are used for the pairs [velocity,
fluid pressure] and [displacement, total pressure], alongside continuous and
piecewise quadratic elements for Biot fluid pressure. Figure 6.5 depicts the
numerical solutions at two time instants, showing the distribution on the
meridional axisymmetric domain of displacement, pressures and velocity
(including the Darcy velocity in the trabecular meshwork) and illustrating the
modification of flow patterns throughout the last cycle of the computation,
and using (for this first run of simulations), a constant permeability
$\kappa_{0}=5.0\cdot 10^{-12}\,[\text{m}^{2}]$. For sake of visualisation, we
have also plotted the fluid pressure on both sub-domains rotationally extruded
to a cut of the 3D domain.

Figure 6.6: Example 4. Axisymmetric interfacial flow in the eye. Top: Coarse
unstructured mesh indicating the axis of symmetry. Second row: four
permeability distributions (in $[\text{m}^{2}]$) projected on piecewise
constants over $\Omega_{P}$. Linearly decreasing (left), linearly increasing
(centre-left), with a log-random distribution (centre-right), and with
randomly distributed spots of smaller permeability (right). In all cases the
average permeability is $\bar{\kappa}=2.0\cdot 10^{-11}$. Third-fifth rows:
comparisons of resulting profiles (velocity vectors, fluid pressure, and post-
processed strain) at $t=6$ [s].

In addition, we compare the behaviour produced by four heterogeneous
permeability profiles (2.4). The first case has a gradient going from
$\kappa_{1}^{\max}=2.88\cdot 10^{-11}\,[\text{m}^{2}]$ on the interface,
linearly down to $\kappa^{\min}=10^{-14}\,[\text{m}^{2}]$ on the outlet. A
second synthetic permeability profile will decrease from the value
$\kappa_{2}^{\max}=6.55\cdot 10^{-11}\,[\text{m}^{2}]$ on the interface down
to $\kappa^{\min}$ on the outlet. These two profiles are generated by solving
a Laplace problem with mixed boundary conditions, setting $\kappa_{i}^{\max}$
and $\kappa^{\min}$ essentially, and no-flux naturally on the remainder of
$\Gamma_{P}$. A third permeability distribution is generated by a uniform
random distribution bounded by $\kappa^{\min}$ and
$\kappa_{3}^{\max}=3.99\cdot 10^{-11}\,[\text{m}^{2}]$, and a fourth case is
constructed by placing random points in $\Omega_{P}$ having permeability
$\kappa^{\min}$, and $\kappa_{4}^{\max}=2.28\cdot 10^{-11}\,[\text{m}^{2}]$
elsewhere (see the top panels of Figure 6.6). The maximum values
$\kappa_{i}^{\max}$ were tuned so that the average permeability
$\bar{\kappa}_{i}=\frac{1}{|\Omega_{P}|}\int_{\Omega_{P}}\kappa_{i}(r,z)\mathrm{d}r\,\mathrm{d}z$
is equal to $2.0\cdot 10^{-11}$ in all four cases.

The effect of spatial variations in permeability are evaluated by imposing an
inlet velocity profile (that is, equivalently, controlling the flow) with
$\boldsymbol{u}_{\mathrm{in}}\cdot\boldsymbol{n}=-0.1t\,[\text{m}\cdot\text{s}^{-1}]$
while prescribing zero fluid Biot pressure on $\Gamma^{\text{out}}$. From the
third row of Figure 6.6 we see that the pressure difference (intra-ocular
fluid pressure minus the pressure at the angular aqueous plexus) remains
roughly the same in all the cases consistently with the fact that
$\bar{\kappa}$ is kept constant. Nonetheless we see that different
permeability profiles give rise to slightly different spatial distributions
for velocity, for the Biot fluid pressure, and for the strain (post-processed
from the Biot displacement and from the harmonic extension Stokes
displacement). Compared with the constant permeability case, that is the
underlying field in the fourth permeability profile (right column), where the
pressure decays linearly toward the outlet, in the decreasing permeability
case (left column), the pressure gradient looks larger close to the inlet
while the strain are, in general everywhere, lower; for the increasing
permeability case (centre-left column) the pressure gradient is larger at the
inlet and the strain are, in general everywhere, larger. Taking into
consideration that here we investigate a much complex geometry, these
speculations are consistent with the description provided in [50] and obtained
in the limit of a slender geometry. The description provided here allows to
speculate that a rearrangement of the microstructure could not drive a
macroscopic change in the mechanics of the system within a short time horizon,
since the flow and the pressure drop remains of the same order of magnitude;
however, higher pressure and strain gradients could be the cause of a local
remodelling of the tissue that, in the long term and in the nonlinear regime,
could lead to a collapse of part of the pectinate ligaments.

##  7 Concluding remarks

We have introduced a new formulation for the coupling of Biot’s poroelasticity
system using total pressure and free flow described by the Stokes and Navier-
Stokes equations, and which does not require Lagrange multipliers to set up
the interface conditions between the two subdomains. The well-posedness of the
continuous problem has been proved, we have provided a rigorous stability
analysis. A mixed finite element method is defined for the proposed
formulation together with a corresponding reduction to the axisymmetric case
and the proposed schemes are robust with respect to the first Lamé parameter.