text:
\section{Introduction}

In this paper we develop and study a domain decomposition method for the quasistatic Biot system of poroelasticity \cite{biot1941general} using mixed finite element subdomain discretizations with non-matching grids along the interfaces. The Biot system models flow of viscous fluids through deformable porous media. The system has a wide rage of applications, including in the geosciences, such as earthquakes, landslides, groundwater cleanup, and hydraulic fracturing, as well as in biomedicine, such as arterial flows and biological tissues. The model consists of an equilibrium equation for the solid medium coupled with a mass balance equation for the fluid flow through the medium. Various numerical methods for the Biot system have been developed in the literature, considering two-field displacement--pressure formulations \cite{Gaspar-FD-Biot,Jan-SINUM-Biot,Zik-MINI}, three-field displacement--pressure--Darcy velocity formulations \cite{phillips2007coupling1,Yi-Biot-locking,Zik-nonconf,Lee-Biot-three-field,phillips-DG,Whe-Xue-Yot-Biot}, three-field displacement--pressure--total pressure formulations
\cite{Lee-Mardal-Winther,ORB}, and four-field stress--displacement--pressure--Darcy velocity mixed formulations \cite{Yi-Biot-mixed,Ahmed-apost-CMAME}. In this work we consider the five-field weakly symmetric stress--displacement--rotation--pressure--Darcy velocity formulation \cite{Lee-Biot-five-field,msfmfe-Biot}. The four-field and five-field formulations lead to mixed finite element (MFE) approximations, which exhibit local mass and momentum conservation, accurate normal-continuous Darcy velocity and solid stress, as well as robust and locking-free behavior for a wide range of physical parameters. An additional advantage of the five-field weakly symmetric MFE method is that it can be reduced to a positive definite cell-centered scheme for the pressure and displacement only, as it is done in the multipoint
stress--multipoint flux MFE method developed in \cite{msfmfe-Biot}, through the use of a vertex quadrature rule and local elimination of some of the variables. We note that our analysis for the weakly symmetric formulation can be carried over to the strongly symmetric formulation found in \cite{Yi-Biot-mixed}.

Numerical methods for the Biot system of poroelasticity usually lead to
large algebraic systems, due to the coupling of unknowns, as well as the size of the domain and the wide range of scales associated with practical applications. Domain decomposition methods \cite{Toselli-Widlund,QV} a commonly used for solving large systems resulting from discretizations of partial differential equations, as they lead to parallel and efficient solution algorithms. In this work we focus on non-overlapping domain decomposition methods, where the domain is split into non-overlapping subdomains and the continuity of the solution variables at the subdomain interfaces is enforced through a suitable interface Lagrange multiplier. The global problem can be reduced to solving iteratively an interface problem, involving the solution of smaller subdomain systems at each iteration, which can be performed in parallel. Despite the extensive studies of numerical methods for the Biot system of poroelasticity, there have been relatively few results on domain decomposition methods for this problem and they have been mostly based on the two-field displacement--pressure formulation \cite{girault2011domain,HoracioNurbs1,FlorezDD,GosseletMonoDD}. To the best of our knowledge, the only paper on domain decomposition for Biot with a mixed formulation is \cite{dd-biot}, which is based on the five-field mixed formulation with weak stress symmetry. Two types of methods are developed in \cite{dd-biot}. One is a monolithic domain decomposition method, which involves solving the Biot system on each subdomain. The second is a partitioned method, which splits the Biot system into solving separate elasticity and flow equations, and applies domain decomposition for each of the equations. The developments in \cite{dd-biot} are motivated by earlier works on non-overlapping domain decomposition methods
for MFE discretizations of Darcy flow \cite{glowinski1988domain,cowsar1995balancing,arbogast2000mixed}
and elasticity \cite{eldar_elastdd}.

The domain decomposition methods in \cite{dd-biot} are limited to subdomain grids that match at the interfaces. In this paper we generalize the work in \cite{dd-biot} to enable the use of non-matching subdomain grids through the use of mortar finite elements \cite{arbogast2000mixed,APWY,Fritz-mortar-elast,Kim-elast-BDDC,pencheva2003balancing,eldar_elastdd}. This generality provides the flexibility to use different grid resolution in different subdomains, as well as a coarser mortar space, resulting in a multiscale approximation. We refer to the method as a multiscale mortar mixed finite element (MMMFE) method. The MMMFE method has been studied for mixed formulations of scalar elliptic equations in \cite{arbogast2000mixed,APWY,ganis2009implementation} and for weakly-symmetric mixed elasticity in \cite{eldar_elastdd}. Following the monolithic domain decomposition method from \cite{dd-biot}, we utilize a physically heterogeneous Lagrange multiplier vector consisting of interface displacement and pressure variables to impose weakly the continuity of the normal components of stress and velocity, respectively. In contrast to \cite{dd-biot}, we choose the Lagrange multiplier vector from a space of mortar finite elements defined on a separate interface grid, which allows for handling non-matching subdomain grids through projections from and onto the mortar finite element space. This also allows for the mortar space to be on a coarser scale $H$, see \cite{PWY,APWY,ganis2009implementation}, compared to a finer subdomain
grid size $h$. The multiscale capability adds an extra layer of flexibility over the methods from \cite{dd-biot}.

The main contributions of this paper are as follows. We first consider the semidiscrete continuous-in-time formulation and establish existence, uniqueness, and stability of the MMMFE method for the Biot system, employing the theory of degenerate evolutionary systems of partial differential equations with monotone operators. For the solvability of the associated resolvent problem we utilize, under a condition on the richness of the mortar finite element space,
an inf-sup condition for the mortar space, as well as inf-sup conditions for the stress and velocity spaces with weak interface continuity of normal components. Next, we establish a priori error estimates for the stress, displacement, rotation, pressure, and Darcy velocity, as well as the displacement and pressure mortar variables in their natural norms. We then consider a fully-discrete method based on the backward Euler time discretization. We show that the solution of the algebraic system at each time step can be reduced to solving a positive definite interface problem for the composite displacement--pressure mortar variable. Motivated by the multiscale flux basis from \cite{ganis2009implementation} and the multiscale stress basis from \cite{eldar_elastdd}, we propose the construction and use of a multiscale stress---flux basis, which makes the number of subdomain solves independent of the number of iterations required for the interface problem. Moreover, since the basis can be reused at each time step, the total number of subdomain solves is independent of the number of time steps. This illustrates that the multiscale basis results in a significant reduction of computational cost in the case of time-dependent problems. Finally, we present the results of several numerical tests designed to illustrate the well-posedness, stability, and accuracy of the proposed MMMFE method. We also consider a test based on data from the Society of Petroleum Engineers SPE10 benchmark, illustrating the multiscale capabilities of the method and the advantages of using a multiscale basis.

The rest of the paper is organized as follows. Section \ref{sec:Math-Formulation-mortar}
introduces the model problem and its domain decomposition mixed finite element approximation. The well-posedness, stability, and error analysis for the semi-discrete formulation is presented in Section \ref{sec:Analysis-of-MMMFE}. In Section \ref{sec:Implementation-Mortar} we discuss the fully discrete method, the reduction to an interface problem, and the construction of the multiscale stress--flux basis. Numerical results are reported in Section \ref{sec:Numerical-results-mortar}.

\section{Formulation of the method\label{sec:Math-Formulation-mortar}}

In this section, we introduce the mathematical model of interest and its
mixed finite element approximation. We also develop the framework for the multiscale mortar mixed finite element domain decomposition method. Projection operators
critical in the analysis of the method and various bounds associated
with them are introduced. Finally, we introduce the weakly continuous
spaces of stress, $\mathbb{X}_{h,0}$, and velocity, $Z_{h,0}$, and reformulate
the MMMFE method in terms of these spaces.

\subsection{Mathematical formulation of model Problem}

Let $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ be a simply connected
domain. We use the notation
$\mathbb{M}$, $\mathbb{S}$ and $\mathbb{N}$ for the spaces of $\ensuremath{d\times d}$
matrices, symmetric matrices, and skew-symmetric matrices, respectively,
all over the field of real numbers. Let $I\in\mathbb{S}$ represents the $d\times d$
identity matrix. The partial derivative operator with respect to time, $\frac{\partial
}{\partial t}$, is often abbreviated to $\partial_t$. $C$ denotes a generic positive
constant that is independent of the discretization parameters
$h$ and $H$. Throughout the thesis,
the divergence operator is the usual divergence for vector fields,
which produces vector field when applied to matrix field by taking
the divergence of each row. The divergence operator is represented by either $\nabla\cdot$ or $\dvr$.

For a set $G\subset\mathbb{R}^{d},$ the $L^{2}(G)$ inner product and norm
are denoted by $ (\cdot,\cdot)_G$ and $\|\cdot\|_G$, respectively,
for scalar, vector, or tensor valued functions. For any fraction $r$,
$\|\cdot\|_{r,G}$ denotes the $H^{r}(G)$-norm, for example $\|\cdot\|_{0,G}=\|\cdot\|_{G}$.
We omit subscript $\ensuremath{G}$ if $\ensuremath{G=\Omega}.$ For
a section of the domain or element boundary $S\subset\mathbb{R}^{d-1}$, we
write $\gnp[\cdot]{\cdot}_{S}$ and $\|\cdot\|_{S}$ for the $L^{2}(S)$
inner product (or duality pairing) and norm, respectively. We will
also use the spaces
\begin{align*}
 & H(\dvr;\Omega)=\{\zeta\in L^{2}(\Omega,\mathbb{R}^{d}):\dvr\zeta\in L^{2}(\Omega)\},\\
 & H(\dvr;\Omega,\mathbb{M})=\{\tau\in L^{2}(\Omega,\mathbb{M}):\dvr\tau\in L^{2}(\Omega,\mathbb{R}^{d})\},
\end{align*}
with the norm $\|\tau\|_{\dvr}=\left(\|\tau\|^{2}+\|\dvr\tau\|^{2}\right)^{1/2}.$
Let $\|\quad\|_{\star}$ be the natural norm associated with the space
$\star\in\left\{ \mathbb{X},V,\mathbb{Q},Z,W\right\} $, where the set of
spaces will be introduced later in this section.

Given a vector field $f$ representing body
forces and a source term $g$, we consider the quasi-static Biot system of poroelasticity (\cite{biot1941general}):
\begin{align}
-\dvr\sigma(u) & =f, & \text{in \ensuremath{\Omega \times (0,T]}},\label{biot-1}\\
{K^{-1}} z+\nabla p & =0, & \text{in \ensuremath{\Omega \times (0,T]}},\label{biot-2}\\
\ddt(c_{0}p+\alpha\dvr u)+\dvr z & =g, & \text{in \ensuremath{\Omega \times (0,T]}},\label{biot-3}
\end{align}
where $u$ is the displacement, $p$ is the fluid pressure, $z$ is the
Darcy velocity, and $\sigma$ is the poroelastic stress, defined as
\begin{equation}\label{stress-comb}
\sigma = \sigma_e - \alpha p I.
\end{equation}
Here $I$ is the $d\times d$ identity matrix, $0 < \alpha \le 1$ is the
Biot-Willis constant, and $\sigma_e$ is the elastic stress satisfying the
stress-strain relationship
\begin{equation}\label{stress-strain}
A\sigma_e = \epsilon(u), \quad \epsilon(u):=\frac{1}{2}\left(\nabla u+(\nabla u)^{T}\right),
\end{equation}
where $A$ is the compliance tensor, which is a symmetric,
bounded and uniformly positive definite linear operator acting from
$\ensuremath{\mathbb{S}\to\mathbb{S}}$, extendible to $\mathbb{M}\to\mathbb{M}$. In the special case of homogeneous
and isotropic body, $A$ is given by,
\begin{equation}
A\sigma=\frac{1}{2\mu}\left(\sigma-\frac{\lambda}{2\mu+d\lambda}\operatorname{tr}(\sigma)I\right),\label{A-defn}
\end{equation}
where $\mu>0$ and $\lambda\ge0$
are the Lam\'e coefficients. In this case, $\sigma_{e}(u)=2\mu\epsilon(u)+\lambda\dvr u\,I$.
Finally, $c_{0} \ge 0$ is the mass storativity and $K$ stands for the permeability
tensor that is spatially-dependent, uniformly bounded, symmetric, and
positive definite, i.e, for constants $0 < k_{\min} \le k_{\max}
< \infty$,
\begin{equation}\label{k-spd}
\forall \text{ a.e. } x \in \O, \quad k_{\min}\zeta^T \zeta \le \zeta^T K(x) \zeta \le k_{\max}\zeta^T \zeta,
\quad \forall \zeta \in \mathbb{R}^d.
\end{equation}
To close the system, we impose the boundary conditions
\begin{align}
u & =g_{u} & \text{on \ensuremath{\Gamma_D^u \times (0,T]}} & , & \sigma n & =0
& \text{on \ensuremath{\Gamma_N^\sigma  \times (0,T]}},\label{biot-bc-1}\\
p & =g_{p} & \text{on \ensuremath{\Gamma_D^p \times (0,T]}} & , & z\cdot n
& =0 & \text{on \ensuremath{\Gamma_N^z \times (0,T]}},\label{biot-bc-2}
\end{align}
where $\ensuremath{\Gamma_D^u\cup\Gamma_N^\sigma=\Gamma_D^p\cup\Gamma_N^z=\partial\O}$ and
$n$ is the outward unit normal vector field on $\partial \O$,
along with the initial condition $p(x,0)=p_0(x)$ in $\Omega$.
Compatible initial data for the rest of the variables can be obtained from
($\ref{biot-1}$) and ($\ref{biot-2}$) at $t=0$.
Well posedness analysis for this system can be found in \cite{SHOWALTER2000310}.

We use the same mixed variational formulation for ($\ref{biot-1}$)--($\ref{biot-bc-2}$)
as discussed in \cite{dd-biot}, which was based on the approach
in \cite{Lee-Biot-five-field}. Recall that we introduced a rotation
Lagrange multiplier $\gamma:=\frac{1}{2}\left(\nabla u-\nabla u^{T}\right)\in\mathbb{N}$,
which was used to impose weakly the symmetry of the stress tensor
$\sigma$. Following the same arguments as in \cite{dd-biot}, we obtain the following mixed variational formulation for ($\ref{biot-1}$)--($\ref{biot-bc-2}$),
find $\ensuremath{(\sigma,u,\gamma,z,p):[0,T]\to\mathbb{X}\times V\times\mathbb{Q}\times Z\times W}$
such that $p(0)=p_{0}$ and for a.e. $t\in(0,T),$
\begin{align}
& \inp[A(\sigma + \alpha p I)]{\tau} + \inp[u]{\dvr{\tau}}+\inp[\gamma]{\tau}=\gnp[g_{u}]{\tau\,n}_{\Gamma_{D}^u}, & \forall\tau\in\mathbb{X},\label{eq:cts1-1}\\
& \inp[\dvr{\sigma}]{v}=-\inp[f]{v}, & \forall v\in V,\label{eq:cts2-1}\\
& \inp[\sigma]{\xi}=0, & \forall\xi\in\mathbb{Q},\label{eq:cts3-1}\\
& \inp[{K^{-1}} z]{q}-\inp[p]{\dvr{q}}=-\gnp[g_{p}]{q\cdot n}_{\Gamma_{D}^p}, & \forall q\in Z,\label{eq:cts4-1}\\
& c_{0}\inp[\partial_t{p}]{w}+\alpha\inp[\partial_t A(\sigma + \alpha pI)]{wI}
+\inp[\dvr{z}]{w}=\inp[g]{w}, & \forall w\in W,\label{eq:cts5-1}
\end{align}
where
\begin{align*}
&\mathbb{X} = \big\{ \tau\in H (\dvr;\Omega,\mathbb{M}) : \tau\,n = 0 \text{ on } \Gamma_N^\sigma  \big\},
\quad V = L^2 (\Omega, \mathbb{R}^d), \quad \mathbb{Q} = L^2 (\Omega, \mathbb{N}), \\
&Z = \big\{ q\in H (\dvr;\Omega) : q\cdot n = 0 \text{ on } \Gamma_N^z  \big\},
\quad W = L^2 (\Omega).
\end{align*}
It was shown in \cite{msfmfe-Biot} that the system
(\ref{eq:cts1-1})--(\ref{eq:cts5-1}) is well posed.

In the MFE setting, the variational formulation corresponding to (\ref{eq:cts1-1})--(\ref{eq:cts5-1})
reads as follows: find $\ensuremath{(\sigma_{h},u_{h},\gamma_{h},z_{h},p_{h}):[0,T]\to}\in\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}\times Z_{h}\times W_{h}$
such that for a.e. $t\in(0,T),$
$\ensuremath{(\sigma_{h},u_{h},\gamma_{h},z_{h},p_{h})}:[0,T]\to\mathbb{X}_{h}\times
V_{h}\times\mathbb{Q}_{h}\times Z_{h}\times W_{h}$ such that, for a.e. $t\in (0,T)$,
\begin{align}
& \inp[A(\sigma_{h} + \alpha p_h I)]{\tau}+\inp[u_{h}]{\dvr{\tau}}+\inp[\gamma_{h}]{\tau}=\gnp[g_{u}]{\tau\,n}_{\Gamma_{D}^u}, & \forall\tau\in\mathbb{X}_{h},\label{eq:fe1-1}\\
& \inp[\dvr{\sigma_{h}}]{v}=-\inp[f]{v}, & \forall v\in V_{h},\label{eq:fe2-1}\\
& \inp[\sigma_{h}]{\xi}=0, & \forall\xi\in\mathbb{Q}_{h},\label{eq:fe3-1}\\
& \inp[{K^{-1}} z_{h}]{q}-\inp[p_{h}]{\dvr{q}}=-\gnp[g_{p}]{q\cdot n}_{\Gamma_{D}^p},
& \forall q\in Z_{h},\label{eq:fe4-1}\\
& c_{0}\inp[\partial_t{p_{h}}]{w}+\alpha\inp[\partial_t A(\sigma_{h} + \alpha p_h I)]{wI}
+\inp[\dvr{z_{h}}]{w}=\inp[g]{w}, & \forall w\in W_{h},\label{eq:fe5-1}
\end{align}

where $\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}\times Z_{h}\times W_{h}\subset\mathbb{X}\times V\times\mathbb{Q}\times Z\times W$
can be chosen from any of the inf-sup stable MFE families and the
discrete initial data, $p_{h,0},$ is obtained from the elliptic projection
of the continuous initial data.

For solving the mechanics part, $\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}$
can be chosen from any of the stable triplets suitable for solving
linear elasticity problem equations with weakly imposed symmetry.
Examples of such triplets include \cite{stenberg1988family,Amara-Thomas,ArnAwaQiu,arnold2007mixed,Awanou-rect-weak,cockburn2010new,gopalakrishnan2012second,brezzi2008mixed,BBF-reduced,FarFor,msmfe-quads,lee2016towards},
where all the listed triplets satisfy the inf-sup condition,
\begin{align}
\forall v\in V_{h}, \xi\in\mathbb{Q}_{h}, \quad
\|v\|+\|\xi\| & \le C\sup_{0\ne\tau\in\mathbb{X}_{h}}\frac{\inp[v]{\dvr{\tau}}+\inp[\xi]{\tau}}{\|\tau\|_{\text{div}}}.\label{eq:inf-sup-elast-1}
\end{align}
 For the flow part, $Z_{h}\times W_{h}$ could be chosen from any
of the stable pressure-velocity pair of MFE spaces such as the Raviart-Thomas
($\mathcal{RT}$) or Brezzi-Douglas-Marini ($\mathcal{BDM}$) spaces.
These families are presented in detail in \cite{brezzi1991mixed},
where it is shown that they satisfy the following inf-sup condition,
\begin{align}
\forall w \in W_h, \quad \|w\| & \le
C\sup_{0\ne q\in Z_{h}}\frac{(\dvr q,w)}{\|q\|_{\text{div}}}\label{eq:inf-sup-darcy-1}.
\end{align}

We also note the following properties of $A$ and $K$ which will be useful in the analysis of our method, there exist constants
$0 < a_{\min} \le a_{\max} < \infty$ and $0 < k_{\min} \le k_{\max} < \infty$ such that
\begin{align}
a_{\min}\|\tau\|^{2}\le\inp[A\tau]{\tau}\le a_{\max}\|\tau\|^{2}, \quad
\forall \, \tau \in \mathbb{X},\label{eq:coercivity-elast-1}\\
k_{\min}\|q\|^{2}\le\inp[K q]{q}\le k_{\max}\|q\|^{2}, \quad \forall \, q \in Z.
\label{eq:coercivity-flow-1}
\end{align}

\subsection{Multiscale mortar domain decomposition method}

Let $\Omega=\cup_{i=1}^{N}\Omega_{i}$
be a union of non-overlapping shape regular polygonal subdomains,
where each subdomain is a union of elements of finite element partition
$\mathcal{T}_{h}$. Let $\Gamma_{i,j}={\partial\Omega}_{i}\cap{\partial\Omega}_{j},\,\Gamma=\cup_{i,j=1}^{N}\Gamma_{i,j},$
and $\Gamma_{i}={\partial\Omega}_{i}\cap\Gamma={\partial\Omega}_{i}\setminus{\partial\Omega}$ denote the interior
subdomain interfaces. The domain discretization technique we develop
in this section is the generalization of the monolithic non-overlapping
domain decomposition technique developed in \cite{dd-biot},
where the sub-domains are allowed to have multiscale non-matching
grids at their interfaces. Let $h_{i}$ be the diameter of the maximal
element in the mesh on $\Omega_{i}$ and define $h=\max_{i}h_{i}$.
For $1\le i\le N$, let $\mathbb{X}_{h,i}\times V_{h,i}\times\mathbb{\mathbb{Q}}_{h,i}\times Z_{h,i}\times W_{h,i}$
be a family of stable mixed finite elements defined on the subdomain
$\Omega_{i}$. These spaces could be chosen from any of the stable
family of spaces discussed in \cite{dd-biot}. Let the finite
element spaces $\mathbb{X}_{h,i},\,V_{h,i},\,\mathbb{Q}_{h,i},\,Z_{h,i}$,
and $W_{h,i}$ contain polynomials of degree less than or equal to
$k\ge1,\,l\ge0,\,j\ge0,\,r\ge0$, and $s\ge0$, respectively. We define
the global FE spaces, defined on $\Omega$, as follows:
\[
\mathbb{X}_{h}=\bigoplus_{1\le i\le N}\mathbb{X}_{h,i},\quad V_{h}=\bigoplus_{1\le i\le N}V_{h,i},\quad\mathbb{Q}_{h}=\bigoplus_{1\le i\le N}\mathbb{\mathbb{Q}}_{h,i},\quad Z_{h}=\bigoplus_{1\le i\le N}Z_{h,i},\quad W_{h}=\bigoplus_{1\le i\le N}W_{h,i}.
\]
The spaces $V_h$, $\mathbb{Q}_{h}$, and $W_h$ are equipped with $L^2$-norms. The spaces $\mathbb{X}_{h}$ and $Z_h$ are equipped with the norms
$$
\|\tau\|_{\mathbb{X}_h}^2 := \|\tau\|^2 + \|\dvr_h \tau\|^2 \quad \mbox{and} \quad
\|\zeta\|_{Z_h}^2 := \|\zeta\|^2 + \|\dvr_h \zeta\|^2,
$$
where for simplicity we define $\dvr_h \varphi |_{\Omega_i} := \dvr (\varphi|_{\Omega_i})$. We note that functions in $\mathbb{X}_{h}$ and $Z_{h}$ do not have continuity of the normal components across the subdomain interfaces.
This discontinuity is addressed using Lagrange multipliers defined
on suitable mortar spaces on the interface $\Gamma$. We use relatively
coarser mortar finite elements satisfying certain coarseness conditions
(which will be discussed in the later sections) to approximate the
traces of the displacement vector and the pressure at the interfaces. Let $\mathcal{T}{}_{H,i,j}$ be a shape regular quasi-uniform
finite element partition of $\Gamma_{i,j}$ constructed using a simplicial
or quadrilateral mesh in $d-1$ dimensions with maximal element diameter
$H.$ Define the global mortar fine element spaces on the union of
sub-domain interfaces, $\Gamma$, to be,
\[
\Lambda_{H}=\bigoplus_{1\le i<j\le N}\begin{pmatrix}\Lambda_{H,i,j}^{u}\\
\Lambda_{H,i,j}^{p}
\end{pmatrix},\,\Lambda_{H}^{u}=\bigoplus_{1\le i<j\le N}\Lambda_{H,i,j}^{u},\text{ and }\Lambda_{H}^{p}=\bigoplus_{1\le i<j\le N}\Lambda_{H,i,j}^{p},
\]
where $\Lambda_{H,i,j}^{u}\subset\left(L^{2}(\Gamma_{i,j})\right)^{d}$
and $\Lambda_{H,i,j}^{p}\subset L^{2}(\Gamma_{i,j})$ are mortar finite
element spaces on $\Gamma_{i,j}$ representing the displacement and
pressure Lagrange multipliers, respectively. We assume that these
mortar spaces contain either continuous or discontinuous polynomials
of degree up to $m\ge0$. Conditions on the degree and richness of
the mortar spaces in order to get a well-posed and stable method will
be discussed in the later sections.

The multiscale mortar domain decomposition formulation for the mixed
Biot problem in a semi-discrete form reads as follows: for $1\le i\le N$,
find $(\sigma_{h,i},u_{h,i},\gamma_{h,i},z_{h,i},p_{h,i},\lambda_{H}):[0,T]\to\mathbb{X}_{h,i}\times V_{h,i}\times\mathbb{Q}_{h,i}\times Z_{h,i}\times W_{h,i}\times\Lambda_{H}$
such that $p_{h,i}(0)=p_{h,0}|_{\Omega_{i}}$ and for a.e. $t\in(0,T),$
\begin{align}
& \inp[A\left(\sigma_{h,i}+\alpha p_{h,i}I\right)]{\tau}_{\Omega_{i}}+\inp[u_{h,i}]{\dvr{\tau}}_{\Omega_{i}}+\inp[\gamma_{h,i}]{\tau}_{\Omega_{i}}\nonumber \\
& \hspace{13em}=\gnp[g_{u}]{\tau\,n_{i}}_{{\partial\Omega}_{i}\cap\Gamma_D^{u}}+\gnp[\lambda_{H}^{u}]{\tau\,n_{i}}_{\Gamma_{i}}, &  & \forall\tau\in\mathbb{X}_{h,i},\label{eq:monodd-mmmfepre1}\\
& \inp[\dvr{\sigma_{h,i}}]{v}_{\Omega_{i}}=-\inp[f]{v}_{\Omega_{i}}, &  & \forall v\in V_{h,i},\label{eq:monodd-mmmfepre2}\\
& \inp[\sigma_{h,i}]{\xi}_{\Omega_{i}}=0, &  & \forall\xi\in\mathbb{Q}_{h,i},\label{eq:monodd-mmmfepre3}\\
& \inp[{K^{-1}} z_{h,i}]{\zeta}_{\Omega_{i}}-\inp[p_{h,i}]{\dvr{\zeta}}_{\Omega_{i}}=-\gnp[g_{p}]{\zeta\cdot n_{i}}_{{\partial\Omega}_{i}\cap\Gamma_D^{p}}-\gnp[\lambda_{H}^{p}]{\zeta\cdot n_{i}}_{\Gamma_{i}}, &  & \forall\zeta\in Z_{h,i},\label{eq:monodd-mmmfepre4}\\
& c_{0}\inp[\partial_t p_{h,i}]{w}_{\Omega_{i}}+\alpha\inp[\partial_t A\left(\sigma_{h,i}+\alpha p_{h,i}I\right)]{wI}_{\Omega_{i}}+\inp[\dvr{z_{h,i}}]{w}_{\Omega_{i}}=\inp[g]{w}_{\Omega_{i}}, &  & \forall w\in W_{h,i},\label{eq:monodd-mmmfepre5}\\
& \sum_{i=1}^{N}\gnp[\sigma_{h,i}\,n_{i}]{\mu^{u}}_{\Gamma_{i}}=0, &  & \forall\mu^{u}\in\Lambda_{H}^{u},\label{eq:monodd-mmmfepre6}\\
& \sum_{i=1}^{N}\gnp[z_{h,i}\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}=0, &  & \forall\mu^{p}\in\Lambda_{H}^{p},\label{eq:monodd-mmmfepre7}
\end{align}
where $n_{i}$ is the outward unit normal vector field on $\Omega_{i}$.
Note that equations (\ref{eq:monodd-mmmfepre6})$\--$(\ref{eq:monodd-mmmfepre7})
enforces a notion of weak continuity of normal components of the stress
tensor and velocity vector across the interface $\Gamma$ and that both the
flow and the elasticity problems are of Dirichlet type.

For simplicity of the analysis, we assume that $\Gamma_{D}^{u}=\Gamma_{D}^{p}={\partial\Omega}$
and $g_{u}=g_{p}=0$, to get the following reformulation of (\ref{eq:monodd-mmmfepre1})$\--$(\ref{eq:monodd-mmmfepre7}):
find $(\sigma_{h},u_{h},\gamma_{h},z_{h},p_{h},\lambda_{H}):[0,T]\to\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}\times Z_{h}\times W_{h}\times\Lambda_{H}$
such that $p_{h,i}(0)=p_{h,0}|_{\Omega_{i}}$ and for a.e. $t\in(0,T),$
\begin{align}
& \inp[A\left(\sigma_{h}+\alpha p_{h}I\right)]{\tau}+\sum_{i=1}^{N}\inp[u_{h}]{\dvr{\tau}}_{\Omega_{i}}+\inp[\gamma_{h}]{\tau}=\sum_{i=1}^{N}\gnp[\lambda_{H}^{u}]{\tau\,n_{i}}_{\Gamma_{i}}, &  & \forall\tau\in\mathbb{X}_{h},\label{eq:monodd-mmmfe1}\\
& \sum_{i=1}^{N}\inp[\dvr{\sigma_{h}}]{v}_{\Omega_{i}}=-\inp[f]{v}, &  & \forall v\in V_{h},\label{eq:monodd-mmmfe2}\\
& \inp[\sigma_{h}]{\xi}=0, &  & \forall\xi\in\mathbb{Q}_{h},\label{eq:monodd-mmmfe3}\\
& \inp[{K^{-1}} z_{h}]{\zeta}-\sum_{i=1}^{N}\inp[p_{h}]{\dvr{\zeta}}_{\Omega_{i}}=\sum_{i=1}^{N}-\gnp[\lambda_{H}^{p}]{\zeta\cdot n_{i}}_{\Gamma_{i}}, &  & \forall\zeta\in Z_{h},\label{eq:monodd-mmmfe4}\\
& c_{0}\inp[\pdrv{p_{h}}]{w}+\alpha\inp[\ddt A\left(\sigma_{h}+\alpha p_{h}I\right)]{wI}+\sum_{i=1}^{N}\inp[\dvr{z_{h}}]{w}_{\Omega_{i}}=\inp[g]{w}, &  & \forall w\in W_{h},\label{eq:monodd-mmmfe5}\\
& \sum_{i=1}^{N}\gnp[\sigma_{h}\,n_{i}]{\mu^{u}}_{\Gamma_{i}}=0, &  & \forall\mu^{u}\in\Lambda_{H}^{u},\label{eq:monodd-mmmfe6}\\
& \sum_{i=1}^{N}\gnp[z_{h}\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}=0, &  & \forall\mu^{p}\in\Lambda_{H}^{p}.\label{eq:monodd-mmmfe7}
\end{align}

\subsection{Projection and interpolation operators}

In this subsection, we discuss various interpolation and projection
operators useful in the analysis of the method.

Let $\mathcal{Q}_{h,i}^{u}:\left(L^{2}({\partial\Omega_i})\right)^{d}\to\mathbb{X}_{h,i}n_{i}$
and $\mathcal{Q}_{h,i}^{p}:L^{2}({\partial\Omega_i})\to Z_{h,i}\cdot n_{i}$
be projection operators onto the trace of the normal components of
$\mathbb{X}_{h,i}$ and $Z_{h,i}$, respectively such that for any $\phi_{u}\in\left(L^{2}({\partial\Omega_i})\right)^{d}$
and $\phi_{p}\in L^{2}({\partial\Omega_i})$,

\begin{align}
& \left<\phi_{u}-\mathcal{Q}_{h,i}^{u}\phi_{u},\,\tau n_{i}\right>_{{\partial\Omega_i}}=0, & \forall\tau\in\mathbb{X}_{h,i},\label{eq:motor-project-1}\\
& \left<\phi_{p}-\mathcal{Q}_{h,i}^{p}\phi_{p},\,\zeta\cdot n_{i}\right>_{{\partial\Omega_i}}=0, & \forall\zeta\in Z_{h,i}.\label{eq:motor-project-2}
\end{align}
We define $\mathcal{Q}_{h,i}:\left(L^{2}({\partial\Omega_i})\right)^{d}\times L^{2}({\partial\Omega_i})\to\mathbb{X}_{h,i}n_{i}\times Z_{h,i}\cdot n_{i}$
as
\begin{equation}
\mathcal{Q}_{h,i}=\begin{pmatrix}\mathcal{Q}_{h,i}^{u}\\
\mathcal{Q}_{h,i}^{p}
\end{pmatrix}.\label{eq:MrtToTraceProj}
\end{equation}
For any inf-sup stable pair of finite element spaces, $\mathbb{X}_{h,i}\times V_{h,i}$,
with $\text{div\,}\mathbb{X}_{h,i}=V_{h,i}$, there exists a mixed canonical
interpolant \cite{brezzi1991mixed}, $\Pi_{i}^{\sigma}:H^{\epsilon}(\Omega_{i},\mathbb{M})\cap\mathbb{X}_{i}\to\mathbb{X}_{h,i}$,
for any $\epsilon>0$, such that for any $\tau\in H^{\epsilon}(\Omega_{i},\mathbb{M})\cap\mathbb{X}_{h,i}$,
\begin{align}
& \left(\text{div\,}(\Pi_{i}^{\sigma}\tau-\tau),\,v\right)_{\Omega_i}=0, & \forall v\in V_{h,i},\label{eq:elast-project-1}\\
& \left<(\Pi_{i}^{\sigma}\tau-\tau)n_{i},\,\hat{\tau} n_{i}\right>_{\Gamma _{i}}=0, & \forall\hat{\tau}\in\mathbb{X}_{h,i},\label{eq:elast-project-2}\\
& \|\Pi_{i}^{\sigma}\tau\|_{\Omega_i}\le C\left(\|\tau\|_{L^{\epsilon}(\Omega_i)}+\|\dvr\tau\|_{\Omega_i}\right).
\end{align}
Similarly for any inf-sup stable pair, $Z_{h,i}\times W_{h,i}$, with
$\text{div }Z_{h,i}=W_{h,i}$, there exists a mixed canonical interpolant
$\Pi_{i}^{z}:\left(H^{\epsilon}(\Omega_{i})\right)^{d}\cap Z_{i}\to Z_{h,i}$
such that for any $\zeta\in\left(H^{\epsilon}(\Omega_{i})\right)^{d}\cap Z_{i}$,
the following holds
\begin{align}
& \left(\text{div\,}(\Pi_{i}^{z}\zeta-\zeta),\,w\right)_{\Omega_i}=0, & \forall w\in W_{h,i},\label{eq:darcy-project-1}\\
& \left<(\Pi_i^{z}\zeta-\zeta)\cdot n,\,\hat{\zeta}\cdot n\right>_{\Gamma _{i}}=0, & \forall\hat{\zeta}\in Z_{h,i},\label{eq:darcy-project-2}\\
& \|\Pi_i^{z}\zeta\|_{Z_{i}}\le C\left(\|\zeta\|_{H^{\epsilon}(\Omega_i)}+\|\dvr\zeta\|_{\Omega_i}\right).\label{eq:darcy-project-3}
\end{align}
Let ${\cal P}_{h,i}^{p}$ denote the $L^2$-orthogonal projection, ${\cal P}_{h,i}^{p}:L^{2}(\Omega_i)\to W_{h,i}$,
such that for any $w\in L^{2}(\Omega_i),$
\begin{align}
& \left({\cal P}_{h,i}^{p}w-w,\,\hat{w}\right)_{\Omega_i}=0, & \forall\hat{w}\in W_{h,i}.\label{eq:pressure-project}
\end{align}

Let ${\cal P}_{h,i}^{u}$ denote the $L^2$-orthogonal projection, ${\cal P}_{h,i}^{u}:\left(L^{2}(\Omega_i)\right)^{d}\to V_{h,i}$,
such that for any $v\in\left(L^{2}(\Omega_i)\right)^{d},$
\begin{align}
& \left({\cal P}_{h,i}^{u}v-v,\,\hat{v}\right)_{\Omega_i}=0, & \forall\hat{v}\in V_{h,i}.\label{eq:disp-project-1}
\end{align}

We also use ${\cal R}_{h,i}$ to denote the orthogonal projection,
${\cal R}_{h,i}:L^{2}(\Omega_i,\mathbb{N})\to\mathbb{Q}_{h,i}$ such that
for any $\xi\in L^{2}(\Omega_i,\mathbb{N})$ ,
\begin{align}
& \left({\cal \cal R}_{h,i}\xi-\xi,\,\hat{\xi}\right)_{\Omega_i}=0, & \forall\hat{\xi}\in\mathbb{Q}_{h,i}.\label{eq:disp-project-1-1}
\end{align}

For the analysis of the method, we will use an elliptic projection
operator, $\hat{\Pi}_{i}^{\sigma}$ onto $\mathbb{X}_{h,i}$ as defined in \cite{eldar_elastdd}. Define $\hat{\Pi}_{i}^{\sigma}:H^{\epsilon}(\Omega_{i},\mathbb{M})\cap\mathbb{X}_{i}\to\mathbb{X}_{h,i}$ as the
operator that takes $\sigma\in H^{\epsilon}(\Omega_{i},\mathbb{M})\cap\mathbb{X}_{i}$ to the finite element
approximation, $\hat{\sigma}$, of the following Neumann problem: for any $\sigma\in H^{\epsilon}(\Omega_{i},\mathbb{M})$, find $(\hat{\sigma},\hat{u},\hat{\gamma})\in\mathbb{X}_{h,i}\times V_{h,i}\times\mathbb{Q}_{h,i}$
such that
\begin{align}
  & \inp[\hat{\sigma}]{\tau}_{\Omega_{i}}+\inp[\hat{u}]{\dvr{\tau}}_{\Omega_{i}}+\inp[\hat{\gamma}]{\tau}_{\Omega_{i}}=\inp[\sigma]{\tau}_{\Omega_i}, &  & \forall\tau\in\mathbb{X}_{h,i},
 \label{ell-proj1} \\
& \inp[\dvr{\hat{\sigma}}]{v}_{\Omega_{i}}=\inp[\dvr{\sigma}]{v}_{\Omega_{i}}, &  & \forall v\in V_{h,i}, \label{ell-proj2} \\
& \inp[\hat{\sigma}]{\xi}_{\Omega_{i}}=\inp[\sigma]{\xi}_{\Omega_{i}}, &  & \forall\xi\in\mathbb{Q}_{h,i}, \label{ell-proj3} \\
& \hat{\sigma}n_{i}=(\Pi_{i}^{\sigma}\sigma)n_{i}\,\,\,\,\text{on \ensuremath{{\partial\Omega}_{i}}}. \label{ell-proj4}
\end{align}
More details on the well-posedness and properties of $\hat{\Pi}_{i}^{\sigma}$
can be found in \cite{eldar_elastdd}. In particular, the following
bounds hold
\begin{align*}
& \|\sigma-\hat{\Pi}_{i}^{\sigma}\sigma\|_{\Omega_i}\le C\|\sigma-\Pi_{i}\sigma\|_{\Omega_i}, & \sigma\in H^{1}\left(\Omega_i,\mathbb{M}\right),\\
& \|\hat{\Pi}_{i}^{\sigma}\sigma\|_{\Omega_i}\le C\left(\|\sigma\|_{H^{\epsilon}(\Omega_i)}+\|\dvr\sigma\|_{\Omega_i}\right). & \sigma\in H^{\epsilon}(\Omega_i,\mathbb{M})\cap\mathbb{X}_{i},\,0<\epsilon\le1.
\end{align*}

We also use the Scott-Zhang interpolants (see \cite{Scott-Zhang}) ${\cal I}_{H}^{u}: H^1(\Gamma) \to \Lambda_{H}^{u} \cap C(\Gamma)$ and ${\cal I}_{H}^{p}: H^1(\Gamma) \to \Lambda_{H}^{p} \cap C(\Gamma)$, defined to preserve the trace on $\partial \Gamma$ for functions that are zero on $\partial \Gamma$.

The operators defined above satisfy the following approximation bounds:

\begin{align}
  \|\psi-{\cal I}_{H}^{u}\psi\|_{t,\Gamma_{i,j}} & \le CH^{\hat{m}-t}\|\psi\|_{\hat{m},\Gamma_{i,j}}, &
  0\le t\le1,
  \,\,
t\le\hat{m}\le m+1,
  \label{eq:inter-proj-1}\\
  \|\psi-{\cal I}_{H}^{p}\psi\|_{t,\Gamma_{i,j}} & \le CH^{\hat{m}-t}\|\psi\|_{\hat{m},\Gamma_{i,j}}, &
  0\le t\le1,
  \,\,
t\le\hat{m}\le m+1,
  \label{eq:inter-proj-1a}\\
\|v-{\cal P}_{h,i}^{u}v\|_{\Omega_{i}} & \le Ch^{\hat{l}}\|v\|_{\hat{l},\Omega_i}, & 0\le\hat{l}\le l+1,\label{eq:inter-proj-3}\\
\|\zeta-{\cal P}_{h,i}^{p}\zeta\|_{\Omega_{i}} & \le Ch^{\hat{s}}\|\zeta\|_{\hat{s},\Omega_i}, & 0\le\hat{s}\le s+1,\label{eq:inter-proj-4}\\
\|\xi-{\cal R}_{h,i}\xi\|_{\Omega_i} & \le Ch^{\hat{j}}\|\xi\|_{\hat{j},\Omega_i}, & 0\le\hat{j}\le j+1,\label{eq:inter-proj-6}\\
\|\psi-{\cal Q}_{h,i}^{u}\psi\|_{\Gamma_{i,j}} & \le Ch^{\hat{k}+t}\|\psi\|_{\hat{k},\Gamma_{i,j}}, & 0\le\hat{k}\le k+1,
\label{eq:inter-proj-8}\\
\|\psi-{\cal Q}_{h,i}^{p}\psi\|_{\Gamma_{i,j}} & \le Ch^{\hat{r}+t}\|\psi\|_{\hat{r},\Gamma_{i,j}}, & 0\le\hat{r}\le r+1,
\label{eq:inter-proj-9}\\
\|\tau-\hat{\Pi}_{i}^{\sigma}\tau\|_{\Omega_{i}} & \le Ch^{\hat{k}}\|\tau\|_{\hat{k},\Omega_i}, & 0\le\hat{k}\le k+1,\label{eq:inter-proj-14}\\
\|\zeta-\Pi_{i}^{z}\zeta\|_{\Omega_{i}} & \le Ch^{\hat{r}}\|\zeta\|_{\hat{r},\Omega_i}, & 0\le\hat{r}\le r+1,\label{eq:inter-proj-15}\\
\|\text{div}(\tau-\hat{\Pi}_{i}^{\sigma}\tau)\|_{\Omega_{i}} & \le Ch^{\hat{l}}\|\text{div\,}\tau\|_{\hat{l},\Omega_i}, & 0\le\hat{l}\le l+1,\label{eq:inter-proj-5}\\
\|\text{div}(\zeta-\Pi_{i}^{z}\zeta)\|_{\Omega_{i}} & \le Ch^{\hat{s}}\|\text{div\,}\tau\|_{\hat{s},\Omega_i}, & 0\le\hat{s}\le s+1,\label{eq:inter-proj-11}
\end{align}
where the functions $\psi,\,v,\,\zeta,\,\tau,$ and $\xi$ are taken
from the domains of the operators acting on them. Bound (\ref{eq:inter-proj-1})
can be found in \cite{Scott-Zhang}, bounds (\ref{eq:inter-proj-3})$\--$(\ref{eq:inter-proj-9}) and (\ref{eq:inter-proj-5})$\--$(\ref{eq:inter-proj-11}) are standard $L^2$-projection approximation bounds \cite{ciarlet2002finite},
and bounds (\ref{eq:inter-proj-14})$\--$(\ref{eq:inter-proj-15})
can be found in \cite{brezzi1991mixed,roberts1991mixed,eldar_elastdd}.

We will also use the trace inequalities
\begin{align}
& \|\psi\|_{t,\Gamma_{i,j}} \le C\|\psi\|_{t+\frac{1}{2},\Omega_i}, \quad t>0,\label{eq:trace-1}\\
& \left\langle \psi,\tau n\right\rangle _{{\partial\Omega_i}} \le C\|\psi\|_{\frac{1}{2},{\partial\Omega_i}}\|\tau\|_{H(\text{div};\Omega_i)},
\quad
\left\langle \psi,\zeta\cdot n\right\rangle _{{\partial\Omega_i}} \le C\|\psi\|_{\frac{1}{2},{\partial\Omega_i}}\|\zeta\|_{H(\text{div};\Omega_i)},
\label{eq:trace-2}
\end{align}
which can be found in \cite{grisvard2011elliptic} and \cite{brezzi1991mixed,roberts1991mixed},
respectively.

Finally, define the projection operators $\hat{\Pi}^{\sigma},\,\Pi^{z},\,{\cal P}_{h}^{p},\,{\cal P}_{h}^{u},{\cal \,R}_{h}$
on respective spaces defined in global domain, $\Omega$, to be the
piece-wise application of $\hat{\Pi}_{i}^{\sigma},\,\Pi_{i}^{z},\,{\cal P}_{h,i}^{p},\,{\cal P}_{h,i}^{p},\,{\cal R}_{h,i}$,
respectively on subdomains $\Omega_{i}$ for $i=1,\dots,N$.

\subsection{Spaces of weakly continuous stress and velocity}

In this section, we introduce the spaces of weakly continuous stress tensors
and velocity vectors, which are defined as follows:
\[
\mathbb{X}_{h,0}=\left\{ \tau\in\mathbb{X}_{h}:\sum_{i=1}^{N}\gnp[\tau n_{i}]{\mu^{u}}_{\Gamma_{i}}=0,\,\,\,\forall\mu^{u}\in\Lambda_{H}^{u}\right\}
\]
and
\[
Z_{h,0}=\left\{ \zeta\in Z_{h}:\sum_{i=1}^{N}\gnp[\zeta\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}=0,\,\,\,\forall\mu^{p}\in\Lambda_{H}^{p}\right\} .
\]
In order to find a priori error estimates for the method (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe7})
using techniques developed for single domain system (\ref{eq:fe1-1})$\--$(\ref{eq:fe5-1})
in \cite{msfmfe-Biot}, we restate (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe7})
in terms of $\mathbb{X}_{h,0}$ and $Z_{h,0}$ as follows: find
$\ensuremath{(\sigma_{h},u_{h},\gamma_{h},z_{h},p_{h})}:[0,T]\to \left(\mathbb{X}_{h,0},V_h,\mathbb{Q}_{h},Z_{h,0},W_h \right)$
such that $p_{h}(0)=p_{h,0}$ and
\begin{align}
& \inp[A\left(\sigma_{h}+\alpha p_{h}I\right)]{\tau}+\sum_{i=1}^{N}\inp[u_{h}]{\dvr{\tau}}_{\Omega_{i}}+\inp[\gamma_{h}]{\tau}=0, & \forall\tau\in\mathbb{X}_{h,0},\label{eq:fe-mono-weak-1}\\
& \sum_{i=1}^{N}\inp[\dvr{\sigma_{h}}]{v}_{\Omega_{i}}=-\inp[f]{v}{}, & \forall v\in V_{h},\label{eq:fe-mono-weak-2}\\
& \inp[\sigma_{h}]{\xi}=0, & \forall\xi\in\mathbb{Q}_{h},\label{eq:fe-mono-weak-3}\\
& \inp[{K^{-1}} z_{h}]{\zeta}-\sum_{i=1}^{N}\inp[p_{h}]{\dvr{\zeta}}_{\Omega_{i}}=0, & \forall\zeta\in Z_{h,0},\label{eq:fe-mono-weak-4}\\
& c_{0}\inp[\partial_t p_{h}]{w}+\alpha\inp[\partial_t A\left(\sigma_{h}+\alpha p_{h}I\right)]{wI}+\sum_{i=1}^{N}\inp[\dvr{z_{h}}]{w}{}_{\Omega_{i}}=\inp[g]{w}, & \forall w\in W_{h}.\label{eq:fe-mono-weak-5}
\end{align}
Note that constructing basis functions for function spaces, $\mathbb{X}_{h,0}$
and $Z_{h,0}$, is difficult and we use the above formulation only
for the sake of error analysis. In the later sections, we will present
a reduction to interface problem approach to design the numerical
algorithm using any of the popular sub-domain spaces $\mathbb{X}_{h,i}$ and
$Z_{h,i}$ discussed earlier.

\section{Analysis of  the MMMFE Method}\label{sec:Analysis-of-MMMFE}

In this section, we present well-posedness and error analysis of the DD formulation
developed in the previous section. We start out by proving inf-sup
stability bounds for weakly continuous stress, $\mathbb{X}_{h,0}$, and velocity,
$Z_{h,0}$, spaces under appropriate conditions on the mortar space, $\Lambda_H$.
Under the same conditions, we show that the multiscale mortar DD method
is well-posed and stable. We finish the section by proving
a combined a priori error bound for all the variables in the formulation.

\subsection{Inf-sup stability for the weakly continuous spaces\label{subsec:Existence-of-inf-sup}}

In this subsection, we give inf-sup stability bounds for the weakly continuous
stress, $\mathbb{X}_{h,0}$, and velocity, $Z_{h,0}$, spaces under appropriate
conditions on the mortar space, $\Lambda_H$.

\begin{assump}
	The mortar space $\Lambda_{H}$ is chosen so that there exists a positive
	constant $C$ independent of $H$ and $h$ such that the following
	inequality holds:
	\begin{equation}
	\|\mu\|_{\Gamma_{i,j}}\le C\left(\|\mathcal{Q}_{h,i}\mu\|_{\Gamma_{i,j}}+\|\mathcal{Q}_{h,j}\mu\|_{\Gamma_{i,j}}\right),\quad\forall\mu\in\Lambda_{H},\quad1\le i<j\le n.\label{eq:mortar_assumption}
	\end{equation}

\end{assump}

\begin{remark}\label{rem:mortar_assum}

	Note that assumption \eqref{eq:mortar_assumption} implies that the
	space $\Lambda_{H}$ cannot be too rich compared to subdomain stress-velocity
	FE spaces (similar approach to \cite{APWY}) in the sense that $\Lambda_{H}^{u}$ and $\Lambda_{H}^{p}$
	are well controlled by their projections on to the normal traces of
	stress and velocity sub-domain spaces respectively. In practice, this
	condition can be easily obtained by taking a coarser mortar mesh satisfying
	$h<H\le1$(see \cite{APWY,arbogast2000mixed,pencheva2003balancing}).

\end{remark}
\begin{lem}
	Under the assumption \eqref{eq:mortar_assumption}, there exists a
	constant $\beta_{D}>0$, independent of $h$ and $H$ such that for
	any $\mu^{p}\in\Lambda_{H}^{p}$, the following holds:
	\begin{equation}
	\|\mu^{p}\|_{\Gamma}\le\beta_{D}\sup_{0\ne\zeta\in Z_{h}}\frac{\sum_{i=1}^{N}\gnp[\zeta\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}}{\|\zeta\|_{Z_{h}}}.\label{eq:mortar-p-inf-sup}
	\end{equation}
\end{lem}

\begin{proof}
	We start with any $\mu^{p}\in\Lambda_{H}^{p}$ and extend it by zero
	on ${\partial\Omega}$. Let $\phi_{i}$ be the solution to the following auxiliary
	problem
	\begin{align}
	& \dvr\nabla\phi_{i}=\overline{\Qc^{p}_{h,i}\mu^{p}}, & \text{in }\Omega_i,\label{eq:m-inf-sup-1}\\
	& \nabla\phi_{i}\cdot n_{i}=\Qc^{p}_{h,i}\mu^{p}, & \text{on }{\partial\Omega_i},\label{eq:m-inf-sup-2}
	\end{align}
	where $\overline{\Qc^{p}_{h,i}\mu^{p}}$ denotes the mean value of $\Qc^{p}_{h,i}\mu^{p}$
	on ${\partial\Omega_i}$. The above problem can be reformulated in the mixed form
	by defining $\psi_{i}=\nabla\phi_{i}$. The aforementioned elliptic
	problem is well-posed and the elliptic regularity (see \cite{grisvard2011elliptic})
	gives
	\begin{equation}
	\|\psi_{i}\|_{1/2,\Omega_i}+\|\dvr\psi\|_{\Omega_i}\le C\|\Qc^{p}_{h,i}\mu^{p}\|_{{\partial\Omega_i}}.\label{eq:m-inf-sup-3}
	\end{equation}
	Take $\zeta_{h,i}=\Pi^{z}\psi_{i}\in Z_{h,i}$ and (\ref{eq:darcy-project-2})
	combined with (\ref{eq:m-inf-sup-2}) implies that $\zeta_{h,i}\cdot n_{i}=\Qc^{p}_{h,i}\mu$
	on ${\partial\Omega_i}.$ This equality along with the definition of $\Qc^{p}_{h,i}$ imply
	\begin{align}
	& \sum_{i=1}^{N}\gnp[\zeta_{h,i}\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}=\sum_{i=1}^{N}\gnp[\Pi^{z}\psi_{i}\cdot n_{i}]{\mu^{p}}_{{\partial\Omega_i}}=\sum_{i=1}^{N}\gnp[\Pi^{z}\psi_{i}\cdot n_{i}]{\Qc^{p}_{h,i}\mu^{p}}_{{\partial\Omega_i}}\nonumber \\
	& \quad\quad\quad=\sum_{i=1}^{N}\gnp[\Qc^{p}_{h,i}\mu^{p}]{\Qc^{p}_{h,i}\mu^{p}}_{{\partial\Omega_i}}\ge C\sum_{i=1}^{N}\|\mu\|_{\Gamma_{i}},\label{eq:m-inf-sup-4}
	\end{align}
	where we have used the mortar coarseness assumption (\ref{eq:mortar_assumption}).

	Next, we note that
	\begin{equation}
	\|\zeta_{h,i}\|_{Z_i}\le C\sum_{i=1}^{N}\|\mu\|_{\Gamma_{i}},\label{eq:m-inf-sup-5}
	\end{equation}
	which follows from the stability of the canonical projection $\Pi_i^{z}$,
	(\ref{eq:darcy-project-3}), with $\epsilon=1/2$, (\ref{eq:m-inf-sup-3})
	and the stability of $\Qc^{p}_{h,i}$.

	Finally, (\ref{eq:m-inf-sup-5}) combined with (\ref{eq:m-inf-sup-4})
	and defining $\zeta:=\zeta_{h,i}$ on $\Omega_i$ completes the proof.
\end{proof}
\begin{lem}
	Under the assumption \eqref{eq:mortar_assumption}, there exists a
	constant $\beta_{E}>0$, independent of $h$ and $H$ such that for
	any $\mu^{u}\in\Lambda_{H}^{u}$, the following bound holds
	\begin{equation}
	\|\mu^{u}\|_{\Gamma}\le\beta_{E}\sup_{0\ne\tau\in\mathbb{X}_{h}}\frac{\sum_{i=1}^{N}\gnp[\tau n_{i}]{\mu^{u}}_{\Gamma_{i}}}{\|\tau\|_{\mathbb{X}_{h}}}.\label{eq:mortar-u-inf-sup}
	\end{equation}
\end{lem}

\begin{proof}
	The proof follows similar arguments as in the proof of the previous
	lemma, starting with any $\mu^{u}\in\Lambda_{H}^{u}$ and using elliptic
	regularity of corresponding elliptic problem and stability of the projections
	$\Pi_{i}^{\sigma}$ and $\Qc^{u}_{h,i}.$
\end{proof}
\begin{lem}
	\label{lem:operator_elast}Under the assumption \eqref{eq:mortar_assumption},
	there exists a linear operator $\Pi_{0}^{\sigma}:H^{\frac{1}{2}+\epsilon}(\Omega,\mathbb{M})\cap\mathbb{X}\to\mathbb{X}_{h,0}$
	for any $\epsilon>0$, such that for any $\tau\in H^{\frac{1}{2}+\epsilon}(\Omega,\mathbb{M})\cap\mathbb{X}$,
	\begin{align}
& \sum_{i=1}^{N}\left({\rm div} (\Pi_{0}^{\sigma}\tau-\tau),\,v\right)_{\Omega_{i}}=0, & \forall v\in V_h,\label{eq:pi-weak-1}\\
	& \left(\Pi_{0}^{\sigma}\tau-\tau,\,\xi\right)=0, & \forall\xi\in\mathbb{Q}_{h},\label{eq:pi-weak-2}\\
& \|\Pi_{0}^{\sigma}\tau\|\le C\left(\|\tau\|_{\frac{1}{2}+\epsilon}+\|\text{\emph{div\,}}\tau\|\right),\label{eq:pi-weak-3}\\
& \|\Pi_{0}^{\sigma}\tau-\tau\|\le C\left(\sum_{i=1}^{N} h^{\tilde{t}}\|\tau\|_{\tilde{t},\Omega_i} + h^{\tilde{k}}H^{\frac{1}{2}}\|\tau\|_{\tilde{k}+\frac{1}{2}}\right), & 0 < \tilde{t}\le k+1,\,0<\tilde{k}\le k+1,\label{eq:pi-weak-5}\\
& \|{\rm div}(\Pi_{0}^{\sigma}\tau-\tau)\|_{\Omega_{i}}\le
Ch^{\tilde{l}}\|{\rm div}\,\tau\|_{\tilde{l},\Omega_i}, & 0\le\tilde{l}\le l+1.
\label{eq:pi-weak-div}
\end{align}
\end{lem}

\begin{proof}

The proof is based on constructing $ \Pi_{0}^{\sigma}\tau|_{\partial\Omega_i} = \hat{\Pi}^\sigma_i (\tau + \delta \tau_i)$, where the correction $\delta \tau_i$ is designed to give weak continuity of the normal components.
The proof of \eqref{eq:pi-weak-1}--\eqref{eq:pi-weak-5}
is given in \cite[Lemma 4.6]{eldar_elastdd}. Bound \eqref{eq:pi-weak-div} follows from \eqref{eq:pi-weak-1} and the approximation properties of the $L^2$-projection \cite{ciarlet2002finite}.
\end{proof}
\begin{lem}
	\label{lem:operator_darcy}Under the assumption \eqref{eq:mortar_assumption},
	there exists a linear operator $\Pi_{0}^{z}:\left(H^{\frac{1}{2}+\epsilon}(\Omega)\right)^{d}\cap Z\to Z_{h,0}$
	such that for any $1\le i\le N$ and $\zeta\in\left(H^{\frac{1}{2}+\epsilon}(\Omega)\right)^{d}\cap Z,$

\begin{align}
& \sum_{i=1}^{N}\left(\text{\emph{div}}\left(\Pi_{0}^{z}\zeta-\zeta\right),\,w\right)_{\Omega_{i}}=0, & \forall w\in W_{h},\label{eq:pi-weak-d-1}\\
& \|\Pi_{0}^{z}\zeta\|_{Z_{h}}\le C\left(\|\zeta\|_{\frac{1}{2}+\epsilon} + \|\dvr\zeta\|\right),\label{eq:pi-weak-d-2}\\
& \|\Pi_{0}^{z}\zeta-\zeta\| \le C \left( \sum_{i=1}^{N} h^{\tilde{t}}\|\zeta\|_{\tilde{t},\Omega_{i}} + h^{\tilde{r}}H^{\frac{1}{2}}\|\zeta\|_{\tilde{r}+\frac{1}{2}}\right), &
  0 < \tilde{t}\le r+1,\, 0<\tilde{r}\le r+1,\label{eq:pi-weak-d-4}\\
&  \|{\rm div}(\Pi_{0}^{z}\zeta-\zeta)\|_{\Omega_{i}}\le
Ch^{\tilde{s}}\|{\rm div}\,\zeta\|_{\tilde{s},\Omega_i}, & 0\le\tilde{s}\le s+1.
\label{eq:pi-weak-d-div}
\end{align}
\end{lem}

\begin{proof}
The proof follows from the arguments given in \cite[Section 3]{arbogast2000mixed} and \cite[Section 3]{APWY}.
\end{proof}

Lemmas \ref{lem:operator_elast} and \ref{lem:operator_darcy}
along with Fortin's Lemma \cite{girault1986finite,brezzi1991mixed}
imply the following inf-sup stability
bounds with respect to the weakly continuous spaces of stress and
velocity.
\begin{lem}
	\label{thm:inf-sup-weak}Under the assumption \eqref{eq:mortar_assumption},
	there exists positive constants $C_{E}$ and $C_{D}$ independent of the discretization parameters $h$ and $H$ such that
	\begin{align}
          \forall \, v\in V_{h}, \xi\in\mathbb{Q}_{h}, \quad
	\|v\|+\|\xi\| & \le C_{E}\sup_{0\ne\tau\in\mathbb{X}_{h,0}}\frac{\sum_{i=1}^{N}\inp[v]{\dvr{\tau}}_{\Omega_{i}}+\inp[\xi]{\tau}}{\|\tau\|_{\mathbb{X}_{h}}},\label{eq:inf-sup-elast-weak}
	\end{align}
\end{lem}
\begin{align}
  \forall \, w \in W_h, \quad
  \|w\| & \le C_{D}\sup_{0\ne\zeta\in Z_{h,0}}\frac{\sum_{i=1}^{N}\left(\dvr{\zeta},\,w\right)_{\Omega_{i}}}
        {\|\zeta\|_{Z_{h}}}.\label{eq:inf-sup-darcy-weak}
\end{align}

\subsection{Well-posedness of the semi-discrete MMMFE formulation}

In this subsection, we show the existence of a unique solution
to the system of equations (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe7})
under the assumption (\ref{eq:mortar_assumption}). We follow closely
the proof for the well-posedness of the multipoint flux method for
the Biot system given in \cite{msfmfe-Biot}. We base our proof on
the theory for showing the existence of solution to a degenerate parabolic
system \cite{Showalter}. In particular, we use \cite[IV, Theorem 6.1(b)]{Showalter}
which is stated as follows:
\begin{thm}
	\label{thm:algebraic-pde}Let the linear, symmetric, and monotone
	operator $\mathcal{N}$ be given for the real vector space $E$ to
	its algebraic dual $E^{*}$, and let $E_{b}^{'}$ be the Hilbert space
	which is the dual of $E$ with the seminorm $|x|_{b}=\sqrt{\mathcal{N}x(x)}$
	for $x\in E.$ Let $\mathcal{M}\subset E\times E_{b}^{'}$ be a relation
	with the domain $D=\left\{ x\in E:\mathcal{M}(x)\ne\emptyset\right\} .$
	Assume that ${\cal M}$ is monotone and $Range({\cal N+M})=E_{b}^{'}$.
	Then for each $x_{0}\in D$ and for each ${\cal F}\in W^{1,1}\left(0,T;E_{b}^{'}\right)$,
	there is a solution $x$ of
	\begin{eqnarray*}
		\frac{d}{dt}\left(\mathcal{N}x(t)+{\cal M}\left(x(t)\right)\right)\ni{\cal F}(t), &  & \quad a.e.\,\,\,0<t<T,
	\end{eqnarray*}
	with
	\[
	{\cal N}x\in W^{1,\infty}\left(0,T;E_{b}^{'}\right),x(t)\in D,\,\text{for all }0\le t\le T,\text{ and }{\cal N}x(0)={\cal N}x_{0}.
	\]
\end{thm}

Using the above theorem, we now prove that the semi-discrete system
(\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe7}) is well-posed.
\begin{thm}
	\label{thm:well_posedness-mono-mortar}For each $\left(f,g\right)\in W^{1,\infty}\left(0,T;\left(L^{2}(\Omega)\right)^{d}\right)\times W^{1,\infty}\left(0,T;L^{2}(\Omega)\right)$
	and compatible initial data $\left(\sigma_{h,0},u_{h,0},\gamma_{h,0},z_{h,o},p_{h,0},\lambda_{H,0}\right)$,
	the system of equations \eqref{eq:monodd-mmmfe1}$\--$\eqref{eq:monodd-mmmfe7}
	has a unique solution,
	provided that the assumption \eqref{eq:mortar_assumption} holds.
\end{thm}

\begin{proof}
	We start by reformulating (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe7})
	to fit the setting of Theorem \ref{thm:algebraic-pde}. For this purpose,
	we define operators
	\begin{align*}
	& \left(A_{\sigma\s}\sigma_{h},\,\tau\right)=\left(A\sigma_{h},\,\tau\right),\quad\left(A_{\sigma p}\sigma_{h},\,w\right)=\alpha\left(A\sigma_{h},\,wI\right),\quad\left(A_{\sigma u}\sigma_{h},\,v\right)=\sum_{i=1}^{N}\left(\text{div}\,\sigma_{h,i},\,v\right),\\
	& \left(A_{\sigma\gamma}\sigma_{h},\,\xi\right)=\left(\sigma_{h},\,\xi\right),\quad\left(A_{\sigma\lambda}\sigma_{h},\,\mu^{u}\right)=\sum_{i=1}^{N}\left\langle \sigma_{h}n_{i},\,\mu^{u}\right\rangle _{\Gamma_{i}},\quad\left(A_{zz}z_{h},\,\zeta\right)=\left(K^{-1}z_{h},\,\zeta\right),\\
	& \left(A_{zp}z_{h},\,w\right)=-\sum_{i=1}^{N}\left(\text{div}\,z_{h,i},\,w\right),\quad\left(A_{z\lambda}z_{h},\,\mu^{p}\right)=\sum_{i=1}^{N}\left\langle z_{h}\cdot n_{i},\,\mu^{p}\right\rangle _{\Gamma_{i}},\\
	& \left(A_{pp}p_{h},\,w\right)=c_{0}\left(p_{h},\,w\right)+\alpha^{2}\left(Ap_{h}I,\,wI\right).
	\end{align*}

	Let us introduce the new variables $\dot{u}_h$, $\dot{\gamma}_h$, and $\dot{\lambda}^u_H$ representing $\partial_t u_h$, $\partial_t \gamma_h$, and $\partial_t \lambda^u_H$, respectively. We differentiate equation (\ref{eq:monodd-mmmfe1}) in time to get
	\begin{align}
	& \inp[\partial_t A\left(\sigma_{h}+\alpha p_{h}I\right)]{\tau}+\sum_{i=1}^{N}\inp[\dot{u}_h]{\dvr{\tau}}_{\Omega_{i}}+\inp[\dot{\gamma}_h]{\tau}=\sum_{i=1}^{N}\inp[\dot{\lambda}^u_H]{\tau n_{i}}_{\Gamma_{i}}, & \forall\tau\in\mathbb{X}_{h}.\label{eq:fe-mono-weak-diff}
	\end{align}
	Using the above definitions of operators and using (\ref{eq:fe-mono-weak-diff})
	instead of equation (\ref{eq:monodd-mmmfe1}), we can reformulate
	the problem as a system of linear equations
	\begin{eqnarray}
	\frac{d}{dt}\left(\mathcal{N}\dot{x}(t)+{\cal M}\left(\dot{x}(t)\right)\right) & ={\cal F}(t) & 0<t<T,\label{eq:eq:mod_matrix_form}
	\end{eqnarray}
	where
	\begin{align*}
	& \dot{x}=\left(\begin{array}{c}
	\sigma_{h}\\
	\dot{u}_h\\
	\dot{\gamma}_h\\
	z_{h}\\
	p_{h}\\
	\dot{\lambda}^u_H\\
	\lambda_{H}^p
	\end{array}\right),N=\left(\begin{array}{ccccccc}
	A_{\sigma\s} & 0 & 0 & 0 & A_{\sigma p}^{T} & 0 & 0\\
	0 & 0 & 0 & 0 & 0 & 0 & 0\\
	0 & 0 & 0 & 0 & 0 & 0 & 0\\
	0 & 0 & 0 & 0 & 0 & 0 & 0\\
	A_{\sigma p} & 0 & 0 & 0 & A_{pp} & 0 & 0\\
	0 & 0 & 0 & 0 & 0 & 0 & 0\\
	0 & 0 & 0 & 0 & 0 & 0 & 0
	\end{array}\right),\\
	& M=\left(\begin{array}{ccccccc}
	0 & A_{\sigma u}^{T} & A_{\sigma\gamma}^{T} & 0 & 0 & -A_{\sigma\lambda}^{T} & 0\\
	-A_{\sigma u} & 0 & 0 & 0 & 0 & 0 & 0\\
	-A_{\sigma\gamma} & 0 & 0 & 0 & 0 & 0 & 0\\
	0 & 0 & 0 & A_{zz} & A_{zp}^{T} & 0 & A_{z\lambda}^{T}\\
	0 & 0 & 0 & -A_{zp} & 0 & 0 & 0\\
	A_{\sigma\lambda} & 0 & 0 & 0 & 0 & 0 & 0\\
	0 & 0 & 0 & A_{z\lambda} & 0 & 0 & 0
	\end{array}\right),\,\mathcal{F}=\left(\begin{array}{c}
	0\\
	-f\\
	0\\
	0\\
	g\\
	0\\
	0
	\end{array}\right).
	\end{align*}
	The dual space $E_{b}^{'}$ is given by $L^{2}(\Omega,\mathbb{M})\times0\times0\times0\times L^{2}(\Omega)\times0\times0$
	and the condition ${\cal F}\in W^{1,1}\left(0,T;E_{b}^{'}\right)$
	implies that non-zero source terms can appear only in equations with
	time derivatives. This means we have to take $f=0$ in our case. We
	can fix this issue by considering an auxiliary problem that, for each
	$t\in(0,T]$, solves the system
	\begin{equation}
	\left(\begin{array}{cccc}
	A_{\sigma\s} & A_{\sigma u}^{T} & A_{\sigma\gamma}^{T} & -A_{\sigma\lambda}^{T}\\
	-A_{\sigma u} & 0 & 0 & 0\\
	-A_{\sigma\gamma} & 0 & 0 & 0\\
	A_{\sigma\lambda} & 0 & 0 & 0
	\end{array}\right)\left(\begin{array}{c}
	\sigma_{h}^{f}\\
	\partial_t u_{h}^{f}\\
	\partial_t\gamma_{h}^{f}\\
	\partial_t\lambda_{H}^{u,f}
	\end{array}\right)=\left(\begin{array}{c}
	0\\
	-f\\
	0\\
	0
	\end{array}\right).\label{eq:aux-elast-wellpos}
	\end{equation}
	Such an auxiliary system (\ref{eq:aux-elast-wellpos}) is well-posed
	and the proof can be found in \cite{eldar_elastdd}. Now we can subtract
	the solution to (\ref{eq:aux-elast-wellpos}) from the original system
	of equations (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe7})
	to obtain the modified RHS ${\cal F}=\left(A_{\sigma\s}\left(\sigma_{h}^{f}-\partial_t\sigma_{h}^{f}\right),0,0,0,q-A_{\sigma p}\partial_t\sigma_{h}^{f},0,0\right)^{T}$.

	Next we show that $Range({\cal N+M})=E_{b}^{'}$. This can be established
	by showing that the following square linear homogeneous system has
	only the trivial solution: $\ensuremath{(\hat{\sigma}_{h},\hat{u}_{h},\hat{\gamma}_{h},\hat{z}_{h},\hat{p}_{h},\hat{\lambda}_{H})}\in\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}\times Z_{h}\times W_{h}\times\Lambda_{H}$
	such that
	\begin{align}
	& \inp[A\left(\hat{\sigma}_{h}+\alpha\hat{p}_{h}I\right)]{\tau}+\sum_{i=1}^{N}\inp[\hat{u}_{h}]{\dvr{\tau}}_{\Omega_{i}}+\inp[\hat{\gamma}_{h}]{\tau}-\sum_{i=1}^{N}\gnp[\hat{\lambda}_{H}^{u}]{\tau\,n_{i}}_{\Gamma_{i}}=0, & \forall\tau\in\mathbb{X}_{h},\label{eq:well-pose-aux-1}\\
	& \sum_{i=1}^{N}\inp[\dvr{\hat{\sigma}_{h}}]{v}_{\Omega_{i}}=0, & \forall v\in V_{h},\nonumber \\
	& \inp[\hat{\sigma}_{h}]{\xi}=0, & \forall\xi\in\mathbb{Q}_{h},\nonumber \\
	& \inp[{K^{-1}}\hat{z}_{h}]{\zeta}-\sum_{i=1}^{N}\inp[\hat{p}_{h}]{\dvr{\zeta}}_{\Omega_{i}}+\sum_{i=1}^{N}\gnp[\hat{\lambda}_{H}^{p}]{\zeta\cdot n_{i}}_{\Gamma_{i}}=0, & \forall\zeta\in Z_{h},\label{eq:well-pose-aux-4}\\
	& c_{0}\inp[\partial_t{\hat{p}_{h}}]{w}+\alpha\inp[ A\left(\hat{\sigma}_{h}+\alpha\hat{p}_{h}I\right)]{wI}+\sum_{i=1}^{N}\inp[\dvr{\hat{z}_{h}}]{w}{}_{\Omega_{i}}=0, & \forall w\in W_{h},\nonumber \\
	& \sum_{i=1}^{N}\gnp[\hat{\sigma}_{h}n_{i}]{\mu^{u}}_{\Gamma_{i}}=0, & \forall\mu^{u}\in\Lambda_{H}^{u},\nonumber \\
	& \sum_{i=1}^{N}\gnp[\hat{z}_{h}\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}=0, & \forall\mu^{p}\in\Lambda_{H}^{p}.\nonumber
	\end{align}
	Taking appropriate test functions $(\tau,v,\xi,\zeta,w,\mu^{u},\mu^{p})=(\hat{\sigma}_{h},\hat{u}_{h},\hat{\gamma}_{h},\hat{z}_{h},\hat{p}_{h},\hat{\lambda}_{H}^{u},\hat{\lambda}_{H}^{p})$
	in the above system and adding the equations together gives $\|A^{\frac{1}{2}}\left(\hat{\sigma}_{h}+\alpha\hat{p}_{h}I\right)\|^{2}+c_{0}\|\hat{p}_{h}\|^{2}+\|K^{-\frac{1}{2}}\hat{z}_{h}\|^{2}=0$.
	The coercivity of $A$, (\ref{eq:coercivity-elast-1}), and $K$, (\ref{eq:coercivity-flow-1}),
	give $\hat{\sigma}_{h}+\alpha\hat{p}_{h}I=0$ and $\hat{z}_{h}=0$ respectively.
	Further, the inf-sup condition with respect to the weakly continuous
	space $Z_{h,0}$, (\ref{eq:inf-sup-darcy-weak}), implies $\hat{p}_{h}=0$
	and hence we also have $\hat{\sigma}_{h}=0$. Inf-sup condition with
	respect to the weakly continuous space $\mathbb{X}_{h,0}$, (\ref{eq:inf-sup-elast-weak}),
	implies $\hat{u}_{h}=0$ and $\hat{\gamma}_{h}=0$. Finally, (\ref{eq:mortar-p-inf-sup})
	combined with (\ref{eq:well-pose-aux-1}) implies $\hat{\lambda}_{H}^{u}=0$,
	and (\ref{eq:mortar-u-inf-sup}) combined with (\ref{eq:well-pose-aux-4})
	implies $\hat{\lambda}_{H}^{p}=0$. Similar arguments can be used
	to show that ${\cal N}$ and ${\cal M}$ are non-negative and therefore
	due to linearity, monotone.

	Now to completely satisfy the hypothesis of Theorem \ref{thm:algebraic-pde},
	we need compatible initial data $\dot{x}_{0}\in D$ which implies
	${\cal M}\dot{x}_{0}\in E_{b}^{'}$. We first construct compatible
	initial data to the continuous system (\ref{eq:cts1-1})-(\ref{eq:cts5-1}),
	$\ensuremath{(\sigma_{0},u_{0},\gamma_{0},z_{0},p_{0})}^{T}$, from
	continuous initial condition, $p_{0}$, as follows:
	\begin{enumerate}
		\item Solve equations (\ref{eq:cts1-1})$\--$(\ref{eq:cts3-1}) using $p=p_{0}$
		as given initial data to get $\sigma_{0},u_{0},\gamma_{0}$.
		\item Set $z_0=-K\nabla p_o $ and it is easy to show using integration by parts that this choice satisfies equation (\ref{eq:cts4-1}) with $p=p_{0}$.
	\end{enumerate}
	Define $\tilde{x}_{0}=\ensuremath{(\sigma_{0},u_{0},\gamma_{0},z_{0},p_{0},\lambda_{H,0}^{u},\lambda_{H,0}^{p})}^{T}$,
	where $\lambda_{H,0}^{u}=u_{0}|_{\Gamma}$ and $\lambda_{H,0}^{p}=p_{0}|_{\Gamma}$.
	Take the initial data to the system (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfepre7}),
	$x_{0}$, to be the elliptic projection of $\tilde{x}_{0}$. With
	the reduction of the problem to the case with $f=0$, we have $\left(\mathcal{N+M}\right)\tilde{x}_{0}\in E_{b}^{'}$.
	We also have
	\begin{equation}
	\left(\mathcal{N+M}\right)x_{0}=\left(\mathcal{N+M}\right)\tilde{x}_{0},\label{eq:well-pose-initial-data}
	\end{equation}
	which implies ${\cal M}x_{0}=\left(\mathcal{N+M}\right)\tilde{x}_{0}-{\cal N}x_{0}\in E_{b}^{'}$
	. Now for the modified system (\ref{eq:eq:mod_matrix_form}), we
	take the initial data, $\dot{x}_{0}$ to be $(\sigma_{h,0},0,0,z_{h,0},p_{h,0},0,0)$,
	which also satisfies ${\cal M}\dot{x}_{0}\in E_{b}^{'}$. Note that
	initial data $u_{h,0}$ , $\gamma_{h,0}$ and $\lambda_{H,0}^{u}$ are
	not needed to solve (\ref{eq:eq:mod_matrix_form}), but will be used
	later to recover solution to the original problem.

	Now we can apply Theorem \ref{thm:algebraic-pde} to prove the existence
	of a unique solution
	\begin{align*}
	&\dot{x}=\ensuremath{(\sigma_{h},\dot{u}_{h},\dot{\gamma}_{h},z_{h},p_{h},\lambda_H)},
	\end{align*}
	such that $\sigma_{h}(0)=\sigma_{h,0}$ and $p_{h}(0)=p_{h,0}.$ It is also
	easy to see that $z_{h}(0)=z_{h,0}$ by taking $t\to0$ in (\ref{eq:fe-mono-weak-4})
	and using the fact that $z_{h,0}$ and $p_{h,0}$ satisfy (\ref{eq:fe-mono-weak-4}). Finally for each $t\in[0,T]$, we define $u_{h}$ , $\gamma_{h}$,
	and $\lambda_{H}^u$ as follows:
	\[
	u_{h}(t)=u_{h,0}+\int_{0}^{t}\dot{u}_{h}(s)ds,
	\]
	\[
	\gamma_{h}(t)=\gamma_{h,0}+\int_{0}^{t}\dot{\gamma}_{h}(s)ds,
	\]
	and
	\[
	\lambda_{H}^u(t)=\lambda_{H,0}^{u}+\int_{0}^{t}\dot{\lambda}_{H}^{u}(s)ds.
	\]

	It is easy to show that $u_{h}(t)$, $\gamma_{h}(t)$ and $\lambda_{H}^u(t)$
	satisfy equation (\ref{eq:monodd-mmmfe1}). We can indeed verify this
	by integrating the differentiated version of this equation namely
	equation (\ref{eq:fe-mono-weak-diff}) with respect to time from $0$
	to any $t\in(0,T]$ and using the fact that $\sigma_{h,0},\,u_{h,0},\,\gamma_{h,0}$,
	and $\lambda_{H,0}^{u}$ satisfy equation (\ref{eq:monodd-mmmfe1}). This completes the proof that (\ref{eq:elast-project-1})$\--$(\ref{eq:monodd-mmmfepre7})
	has a solution $\ensuremath{(\sigma_{h},u_{h},\gamma_{h},z_{h},p_{h},\lambda_H)}$.
	Uniqueness of the above constructed solution follows from the stability
	bound for the solution variables which will be discussed in the next
	section.
\end{proof}

\subsection{Stability analysis for MMMFE formulation}

In this subsection, we give a stability bound for the system
(\ref{eq:monodd-mmmfe1})$\--$\eqref{eq:monodd-mmmfe7}.
\begin{thm}
	\label{thm:stability-mono-mortar} Under the assumption \eqref{eq:mortar_assumption}, there exists a constant $C>0$,
	independent of discretization parameters $h$ and $H$, and $c_{0}$
	such that the following stability bound holds for the solution of \eqref{eq:monodd-mmmfe1}$\--$\eqref{eq:monodd-mmmfe7}:
	\begin{align*}
	& \|\sigma_{h}\|_{L^{\infty}\left(0,T;\mathbb{X}_h\right)}+\|u_{h}\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|\gamma_{h}\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|z_{h}\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|p_{h}\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}\\
	& +\|\lambda_{H}^u\|_{L^{\infty}\left(0,T;L^{2}\left(\Gamma\right)\right)}+\|\lambda_{H}^p\|_{L^{\infty}\left(0,T;L^{2}\left(\Gamma\right)\right)}+\|\sigma_{h}\|_{L^{2}\left(0,T;\mathbb{X}_{h}\right)}+\|u_{h}\|_{L^{2}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|\gamma_{h}\|_{L^{2}\left(0,T;L^{2}\left(\Omega\right)\right)}
	\\
	& +\|z_{h}\|_{L^{2}\left(0,T;Z_h\right)}+\|p_{h}\|_{L^{2}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|\lambda_{H}^u\|_{L^{2}\left(0,T;L^{2}\left(\Gamma\right)\right)}+\|\lambda_{H}^p\|_{L^{2}\left(0,T;L^{2}\left(\Gamma\right)\right)}
	\\
	&  \hspace{80pt}\le C\big(\|f\|_{H^{1}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|g\|_{H^{1}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|p_{0}\|_{H^{1}(\Omega)}+\|\nabla Kp_{0}\|_{H\left(\text{\emph{div}},\Omega\right)}\big).
	\end{align*}
\end{thm}

\begin{proof}
	We start out by choosing the test functions
	\begin{equation*}
	(\tau,v,\xi,\zeta,w,\mu^{u},\mu^{p})=(\sigma_{h},\partial_t u_{h},\partial_t\gamma_{h},z_{h},p_{h},\partial_t\lambda_{H}^u,\lambda_{H}^p),
	\end{equation*}
	in equations (\ref{eq:fe-mono-weak-diff}) and (\ref{eq:monodd-mmmfe2})$\--$(\ref{eq:monodd-mmmfe7})
	and combining them to get
	\begin{align*}
	& \inp[\partial_t A\left(\sigma_{h}+\alpha p_{h}I\right)]{\sigma_{h}+\alpha p_{h}I}+c_{0}\inp[\partial_t p_{h}]{p_{h}}+\inp[{K^{-1}} z_{h}]{z_{h}}=\left(f,\partial_t u_{h}\right)+\left(q,p_{h}\right).
	\end{align*}
	The above equation
	can be rewritten as
	\begin{equation}
	\frac{1}{2}\partial_t\left(\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)\|^{2}+c_{0}\|p_{h}\|^{2}\right)+\|K^{-\frac{1}{2}}z_{h}\|^{2}=\partial_t\left(f,u_{h}\right)-\left(\partial_t f,u_{h}\right)+\left(g,p_{h}\right).\label{eq:stab-bnd-1}
	\end{equation}
	For any $t\in(0,T]$, we integrate equation (\ref{eq:stab-bnd-1})
	with respect to time from $0$ to $t$ to get
	\begin{align*}
	& \frac{1}{2}\left(\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)(t)\|^{2}+c_{0}\|p_{h}(t)\|^{2}\right)+\int_{0}^{t}\|K^{-\frac{1}{2}}z_{h}\|^{2}ds\\
	& =\frac{1}{2}\left(\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)(0)\|^{2}+c_{0}\|p_{h}(0)\|^{2}\right)+\int_{0}^{t}\left(\left(g,p_{h}\right)-\left(\partial_t f,u_{h}\right)\right)ds+\left(f,u_{h}\right)(t)-\left(f,u_{h}\right)(0).
	\end{align*}
	On applications of the Cauchy-Schwartz and Young's inequalities, we get
	\begin{align}
	& \|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)(t)\|^{2}+c_{0}\|p_{h}(t)\|^{2}+2\int_{0}^{t}\|K^{-\frac{1}{2}}z_{h}\|^{2}ds\nonumber \\
	& \le\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)(0)\|^{2}+c_{0}\|p_{h}(0)\|^{2}+\epsilon\left(\int_{0}^{t}\left(\|p_{h}\|^{2}+\|u_{h}\|^{2}\right)ds+\|u_{h}(t)\|^{2}\right)\label{eq:stab-bnd-2}\\
	& \quad+\frac{1}{\epsilon}\left(\int_{0}^{t}\left(\|g\|^{2}+\|\partial_t f\|^{2}\right)ds+\|f(t)\|^{2}\right)+\|f(0)\|^{2}+\|u_{h}(0)\|^{2}.\nonumber
	\end{align}
	Bounds for $\|u_{h}\|$ and $\|\gamma_{h}\|$ follow from the inf-sup
	condition (\ref{eq:inf-sup-elast-weak}) as follows:
	\begin{align*}
	& \|u_{h}\|+\|\gamma_{h}\|\le C_{E}\sup_{0\ne\tau\in\mathbb{X}_{h,0}}\frac{\sum_{i=1}^{N}\inp[u_{h}]{\dvr{\tau}}_{\Omega_{i}}+\inp[\gamma_{h}]{\tau}}{\|\tau\|_{\mathbb{X}_{h}}}.
	\end{align*}
	We combine the above equation along with (\ref{eq:monodd-mmmfe1})
	and the Cauchy-Schwartz inequality to obtain the bound
	\begin{align}
	& \|u_{h}\|+\|\gamma_{h}\|\le C_{E}\sup_{0\ne\tau\in\mathbb{X}_{h,0}}\frac{\inp[A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)]{A^{\frac{1}{2}}\tau}}{\|\tau\|_{\mathbb{X}_{h}}}\le C\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)\|.\label{eq:stab-bnd-3}
	\end{align}
	Equation (\ref{eq:stab-bnd-3}) also gives,
	\begin{equation}
	\int_{0}^{t}\left(\|u_{h}\|^{2}+\|\gamma_{h}\|^{2}\right)ds\le C\int_{0}^{t}\left(\|\sigma_{h}\|^{2}+\|p_{h}\|^{2}\right)ds.\label{eq:stab-bnd-4}
	\end{equation}
	Further, choose test functions $(\tau,v,\xi,\mu^u)=(\sigma_{h},u_{h},\gamma_{h},\lambda_{H}^u)$
	in (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe3}), (\ref{eq:monodd-mmmfe6})
	and combine the equations along with Cauchy-Schwartz inequality to
	get
	\[
	\|\sigma_{h}\|^{2}\le C\left(\|p_{h}^{2}\|+\epsilon\|u_{h}\|^{2}+\frac{1}{\epsilon}\|f\|^{2}\right).
	\]
	Combining the above inequality with inequality (\ref{eq:stab-bnd-4})
	yields the following bound
	\begin{equation}
	\int_{0}^{t}\left(\|u_{h}\|^{2}+\|\gamma_{h}\|^{2}\right)ds\le C\int_{0}^{t}\left(\|p_{h}^{2}\|+\epsilon\|u_{h}\|^{2}+\frac{1}{\epsilon}\|f\|^{2}\right)ds.\label{eq:stab-bnd-4-1}
	\end{equation}
	Bound for $\|p_{h}\|$ can be obtained from the inf-sup condition
	(\ref{eq:inf-sup-darcy-weak}) and equation (\ref{eq:monodd-mmmfe4})
	as follows:
	\begin{align}
	||p_{h}|| & \le C_{D}\sup_{0\ne\zeta_{h}\in Z_{h,0}}\frac{\sum_{i=1}^{N}\left(\dvr\zeta_{h},p_{h}\right)_{\Omega_{i}}}{||\zeta_{h}||_{Z_{h}}}=C_{D}\sup_{0\ne\zeta_{h}\in Z_{h,0}}\frac{\inp[{K^{-1}} z_{h}]{\zeta_h}}{||\zeta_{h}||_{Z_{h}}}\le C\|{K^{-1}} z_{h}\|,\label{eq:inf-stab-bnd-5}
	\end{align}
	where the last inequality follows from the Cauchy-Schwartz inequality.

	Further, taking test function $v = \text{div}_h\sigma_{h}$ in (\ref{eq:monodd-mmmfe2})
	and using Cauchy-Schwartz inequality yields
	\begin{equation}
	\|\text{div}_h\sigma_{h}\|^2\le\|f\|^2.\label{eq:stab-bnd-6}
	\end{equation}

	Finally, combining inequalities (\ref{eq:stab-bnd-2})$\--$(\ref{eq:stab-bnd-3})
	and (\ref{eq:stab-bnd-4-1})$\--$(\ref{eq:stab-bnd-6}) and taking
	$\epsilon$ small enough give the following
	bound
	\begin{align}
	  & \|A^{\frac{1}{2}}\left(\sigma_{h} + \alpha p_{h}I\right)(t)\|^2 +
          \|\text{div}_h\sigma_{h}(t)\|^2 + \|u_{h}(t)\|^{2}+\|\gamma_{h}(t)\|^{2}+c_{0}\|p_{h}(t)\|^{2}\nonumber \\
	& \quad + \int_{0}^{t}\left(\|\sigma_{h}\|^2 + \|\text{div}_h\sigma_{h}\|^2\right)ds + \int_{0}^{t}\left(\|u_{h}\|^{2}+\|\gamma_{h}\|^{2}+\|K^{-\frac{1}{2}}z_{h}\|^{2}+\|p_{h}\|^{2}\right)ds\nonumber \\
	& \le C\Big(\int_{0}^{t}\left(\|g(s)\|^{2}+\|\partial_t f(s)\|^{2}+\|f(s)\|^{2}\right)ds+\|f(t)\|^{2}+\|\sigma_{h}(0)\|^{2}\nonumber \\
	& \quad\quad+\|u_{h}(0)\|^{2}+\|p_{h}(0)\|^{2}+\|f(0)\|^{2}\Big).\label{eq:stab-bnd-6-0-1}
	\end{align}

Next we give bounds for $\|\text{div}_hz_{h}\|$, $\|K^{-\frac{1}{2}}z_{h}(t)\|$, and $\|p_{h}(t)\|$ for all $t\in(0,t]$, which are independent of $c_{0}$.
We start by choosing test function $w=\text{div}_hz_{h}$ in (\ref{eq:monodd-mmmfe5}) and applying the Cauchy-Schwartz inequality to obtain
\begin{equation}
  \|\text{div}_hz_{h}\| \le C\left(c_{0}\|\partial_t p_{h}\|
  + \|\partial_t A^{\frac{1}{2}}\left(\sigma_{h} + \alpha p_{h}I\right)\| + \|q\|\right).\label{eq:stab-bond-6-1}
\end{equation}
To bound the right hand side of (\ref{eq:stab-bond-6-1}), differentiate
	equations (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe4}),
	and (\ref{eq:monodd-mmmfe6})$\--$(\ref{eq:monodd-mmmfe7}) with
	respect to time and take appropriate test functions, namely
	\[
	(\tau,v,\xi,\zeta,w,\mu^{u},\mu^{p})=(\partial_t\sigma_{h},\partial_t u_{h},\partial_t\gamma_{h},z_{h},\partial_t p_{h},\partial_t\lambda_{H}^u,\lambda_{H}^p),
	\]
	in the differentiated equations and equation (\ref{eq:monodd-mmmfe5}).
	Further, combining the resulting equations and integrating in time
	from $0$ to $t\in(0,T]$ similar to what we did for equations (\ref{eq:stab-bnd-1})$\--$(\ref{eq:stab-bnd-2})
	previously, we obtain
	\begin{align}
	& 2\int_{0}^{t}\left(\|\partial_t A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)\|^{2}+c_{0}\partial_t\|p_{h}\|^{2}\right)ds+\|K^{-\frac{1}{2}}z_{h}(t)\|^{2}\nonumber \\
	& \le\epsilon\left(\int_{0}^{t}\|\partial_t u_{h}\|^{2}ds+\|p_{h}(t)\|^{2}\right)+\frac{1}{\epsilon}\left(\int_{0}^{t}\|\partial_t f\|^{2}ds+\|g(t)\|^{2}\right)\nonumber \\
	& \quad+\int_{0}^{t}\left(\|p_{h}\|^{2}+\|\partial_t g\|^{2}\right)ds+\|K^{-\frac{1}{2}}z_{h}(0)\|^{2}+\|p_{h}(0)\|^{2}+\|g(0)\|^{2}.\label{eq:stab-bnd-7}
	\end{align}
	To bound $\|\partial_t u_{h}\|$ and $\|\partial_t\gamma_{h}\|$, we use the inf-sup
	condition (\ref{eq:inf-sup-elast-weak}) and equation (\ref{eq:monodd-mmmfe1})
	as before, but now in their time differentiated forms to obtain
	\begin{equation}
	\|\partial_t u_{h}\|+\|\partial_t\gamma_{h}\|\le C\|\partial_t A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)\|.\label{eq:stab-bnd-8}
	\end{equation}
	Combining inequalities (\ref{eq:inf-stab-bnd-5}), (\ref{eq:stab-bnd-7}),
	and (\ref{eq:stab-bnd-7}), and using small enough $\epsilon$
	give
	\begin{align}
	& \int_{0}^{t}\left(\|\partial_t A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)\|^{2}+c_{0}\|\partial_t p_{h}\|^{2}+\|\partial_t u_{h}\|^{2}+\|\partial_t\gamma_{h}\|^{2}\right)ds+\|K^{-\frac{1}{2}}z_{h}(t)\|^{2}+\|p_{h}(t)\|^{2}\nonumber \\
	& \quad\le C\left(\int_{0}^{t}\left(\|\partial_t f\|^{2}+\|p_{h}\|^{2}+\|\partial_t g\|^{2}\right)ds+\|g(t)\|^{2}+\|p_{h}(0)\|^{2}+\|g(0)\|^{2}+\|z_{h}(0)\|^{2}\right).\label{eq:stab-bnd-9}
	\end{align}
	Integrating inequality (\ref{eq:stab-bond-6-1}) with respect to time
	from $0$ to $t\in(0,T]$ and combining the resulting integral inequality
	with inequalities (\ref{eq:inf-stab-bnd-5}) and (\ref{eq:stab-bnd-9})
	give
	\begin{equation} \begin{split}
	    & \|p_{h}(t)\|^{2}+\|K^{-\frac{1}{2}}z_{h}(t)\|^{2} +
            \int_{0}^{t}\|\text{div}_hz_{h}\|^{2}ds \le C\Big(\int_{0}^{t}\big(\|f(s)\|^{2}+\|\partial_t f(s)\|^{2}+\|g(s)\|^{2}\nonumber \\
	& \quad+\|\partial_t q(s)\|^{2}\big)ds+\|f(t)\|^{2}+\|g(t)\|^{2}+\|\sigma_{h}(0)\|^{2}+\|p_{h}(0)\|^{2}+\|z_{h}(0)\|^{2}+\|f(0)\|^{2}+\|g(0)\|^{2}\Big).\label{eq:stab-bnd-10}
	\end{split}
	\end{equation}
	Also, the coercivity of $A$ and $K$ given in inequalities (\ref{eq:coercivity-elast-1})
	and (\ref{eq:coercivity-flow-1}), respectively gives
	\begin{align}
	& \|z_{h}\|\le C\|K^{-\frac{1}{2}}z_{h}\|,\label{eq:stab-bnd-11}\\
	& \|\sigma_{h}\|\le C\left(\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)(t)\|+\|p_{h}\|\right).\label{eq:stab-bnd-12}
	\end{align}
	Next, we give a bound for $\|\lambda_{H}^u\|$. Combining (\ref{eq:mortar-u-inf-sup})
	and (\ref{eq:monodd-mmmfe1}) gives
	\begin{align*}
	& \|\lambda_{H}^u\|_{\Gamma}\le C\sup_{0\ne\tau\in\mathbb{X}_{h}}\frac{\sum_{i=1}^{N}\gnp[\tau n_{i}]{\lambda_{H}^u}_{\Gamma_{i}}}{\|\tau\|_{\mathbb{X}_{h}}}\\
	& \quad=C\sup_{0\ne\tau\in\mathbb{X}_{h}}\frac{1}{\|\tau\|_{\mathbb{X}_{h}}}\left(\inp[A\left(\sigma_{h}+\alpha p_{h}I\right)]{\tau}+\sum_{i=1}^{N}\inp[u_{h}]{\dvr{\tau}}_{\Omega_{i}}+\inp[\gamma_{h}]{\tau}\right)\\
	& \quad\quad\le C\left(\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)\|+\|u_{h}\|+\|\gamma_{h}\|\right).
	\end{align*}
	Squaring the above inequality and integrating from time $0$ to $t$,
	we get the following set of inequalities
	\begin{align}
	& \|\lambda_{H}^u\|_{\Gamma}^{2}\le C\left(\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)\|^{2}+\|u_{h}\|^{2}+\|\gamma_{h}\|^{2}\right),\label{eq:stab-bnd-12-1}\\
	& \int_{0}^{t}\|\lambda_{H}^u(s)\|_{\Gamma}^{2}ds\le C\int_{0}^{t}\left(\|A^{\frac{1}{2}}\left(\sigma_{h}+\alpha p_{h}I\right)\|^{2}+\|u_{h}\|^{2}+\|\gamma_{h}\|^{2}\right)ds.\label{eq:stab-bnd-12-2}
	\end{align}
	Similarly, we can bound $\|\lambda_{H}^p\|$ combining (\ref{eq:mortar-p-inf-sup})
	and (\ref{eq:monodd-mmmfe4}) to obtain
	\begin{align*}
	& \|\lambda_{H}^p\|_{\Gamma}\le C\sup_{0\ne\zeta\in Z_{h}}\frac{\sum_{i=1}^{N}\gnp[\zeta\cdot n_{i}]{\lambda_{H}^p}_{\Gamma_{i}}}{\|\zeta\|_{Z_{h}}}\\
	& \quad=C\sup_{0\ne\zeta\in Z_{h}}\frac{1}{\|\zeta\|_{Z_{h}}}\left(-\inp[{K^{-1}} z_{h}]{\zeta}+\inp[p_{h}]{\dvr{\zeta}}\right)\le C\left(\|K^{-\frac{1}{2}}z_{h}\|+\|p_{h}\|\right),
	\end{align*}
	and using the above inequality to obtain
	\begin{align}
	& \|\lambda_{H}^p\|_{\Gamma}^{2}\le C\left(\|K^{-\frac{1}{2}}z_{h}\|^{2}+\|p_{h}\|^{2}\right),\label{eq:stab-bnd-12-3}\\
	& \int_{0}^{t}\|\lambda_{H}^p\|_{\Gamma}^{2}ds\le C\int_{0}^{t}\left(\|K^{-\frac{1}{2}}z_{h}\|^{2}+\|p_{h}\|^{2}\right)ds.\label{eq:stab-bnd-12-4}
	\end{align}

	In order to bound the initial data, $\sigma_{h}(0),\,u_{h}(0),\,z_{h}(0),$
	and $p_{h}(0)$, note that we obtain the discrete initial data by
	taking elliptic projection of the continuous initial data $(\sigma_{0},u_{0},\gamma_{0},z_{0},p_{0})$
	to the continuous problem (\ref{eq:cts1-1})$\--$(\ref{eq:cts5-1}),
	see (\ref{eq:well-pose-initial-data}). Further note that the continuous
	initial data is constructed using the original pressure initial data
	$p_{0}\in H^{1}(\Omega)$ using the procedure mentioned in Section
	\ref{sec:Math-Formulation-mortar}. Following
	the arguments used in the proof so far for the steady-state version
	with $t=0$ gives
	\begin{align}
	& \|\sigma_{h}(0)\|+\|u_{h}(0)\|+\|\gamma_{h}(0)\|+\|z_{h}(0)\|+\|p_{h}(0)\|\le C\left(\|\sigma_{0}\|+\|u_{0}\|+\|\gamma_{0}\|+\|z_{0}\|+\|p_{0}\|\right)\nonumber \\
	& \quad\quad\quad\le C\left(\|p_{0}\|_{H^{1}(\Omega)}+\|K\nabla p_{0}\|_{H(\text{div};\Omega)}\right).\label{eq:stab-bnd-13}
	\end{align}

	Finally, we combine inequalities (\ref{eq:stab-bnd-6-0-1}), (\ref{eq:stab-bnd-10})$\--$(\ref{eq:stab-bnd-13})
	along with the fact that all the results derived so far hold for a
	general $t\in(0,T]$, to arrive at the stability bound in the theorem.
\end{proof}
\subsection{Error analysis}
In this subsection, we will establish a combined a priori error estimate
for all the unknowns in the formulation.
\begin{thm}
	\label{thm:Error-theorem}Let $\ensuremath{(\sigma_{h}(t),u_{h}(t),\gamma_{h}(t),z_{h}(t),p_{h}(t),\lambda_H)}\in\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}\times Z_{h}\times W_{h}\times\Lambda_{H}$
	be the solution to the system of equations \eqref{eq:monodd-mmmfe1}$\--$\eqref{eq:monodd-mmmfe7}
	under the assumption \eqref{eq:mortar_assumption} for $t\in[0,T]$,
	and suppose the solution of \eqref{eq:cts1-1}$\--$\eqref{eq:cts5-1}
	is sufficiently smooth, then there exists a positive constant $C$,
	independent of $h$, $H$ and $c_{0}$ such that the following holds:
	\begin{equation*}\begin{split}
	& \|\sigma-\sigma_{h}\|_{L^{\infty}\left(0,T;\mathbb{X}_h\right)}+\|u-u_{h}\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|\gamma-\gamma_{h}\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|z-z_{h}\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}\\
	& +\|p-p_{h}\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|u-\lambda_{H}^u\|_{L^{\infty}\left(0,T;L^{2}\left(\Gamma\right)\right)}+\|p-\lambda_{H}^p\|_{L^{\infty}\left(0,T;L^{2}\left(\Gamma\right)\right)}+\|\sigma-\sigma_{h}\|_{L^{2}\left(0,T;\mathbb{X}_h\right)}\\
	& +\|u-u_{h}\|_{L^{2}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|\gamma-\gamma_{h}\|_{L^{2}\left(0,T;L^{2}\left(\Omega\right)\right)}+\|z-z_{h}\|_{L^{2}\left(0,T;Z_h\right)}+\|p-p_{h}\|_{L^{2}\left(0,T;L^{2}\left(\Omega\right)\right)}\\
	& +\|u-\lambda_{H}^u\|_{L^{2}\left(0,T;L^{2}\left(\Gamma\right)\right)}+\|u-\lambda_{H}^p\|_{L^{2}\left(0,T;L^{2}\left(\Gamma\right)\right)} \le C\Big(h^{k_{1}}\|\sigma\|_{H^{1}\left(0,T;H^{k_{1}}\left(\Omega\right)\right)}\\
	&\quad+h^{k_{2}}H^{\frac{1}{2}}\|\sigma\|_{H^{1}\left(0,T;H^{k_{2}+\frac{1}{2}}\left(\Omega\right)\right)}+h^{l_{1}}\|\dvr\sigma\|_{L^{\infty}\left(0,T;H^{l_{1}}\left(\Omega\right)\right)}+h^{l_{2}}\|\dvr\sigma\|_{L^{2}\left(0,T;H^{l_{2}}\left(\Omega\right)\right)}\\
	& \quad +h^{l_{3}}\|u\|_{L^{2}\left(0,T;H^{l_{3}}\left(\Omega\right)\right)}+h^{l_{4}}\|u\|_{L^{\infty}\left(0,T;H^{l_{4}}\left(\Omega\right)\right)}+h^{j_{1}}\|\gamma\|_{H^{1}\left(0,T;H^{j_{1}}\left(\Omega\right)\right)}\\
	& \quad +h^{r_{1}}\|z\|_{H^{1}\left(0,T;H^{r_{1}}\left(\Omega\right)\right)}+h^{r_{2}}H^{\frac{1}{2}}\|z\|_{H^{1}\left(0,T;H^{r_{2}+\frac{1}{2}}\left(\Omega\right)\right)}+h^{s_{1}}\|\dvr z\|_{L^{2}\left(0,T;H^{s_{1}}\left(\Omega\right)\right)}\\
	& \quad +h^{s_{2}}\|p\|_{H^{1}\left(0,T;H^{s_{2}}\left(\Omega\right)\right)}+H^{m_{1}-\frac{1}{2}}\|u\|_{H^{2}\left(0,T;H^{m_{1}+\frac{1}{2}}\left(\Omega\right)\right)}+H^{m_{2}-\frac{1}{2}}\|p\|_{H^{1}\left(0,T;H^{m_{2}+\frac{1}{2}}\left(\Omega\right)\right)}\Big),\\
	& \quad\quad1\le k_{1}\le k+1,\quad\,0\le k_{2}\le k+1,\quad\,0\le l_{1},l_{2},l_{3},l_{4}\le l+1,\quad\,0\le j_{1}\le j+1,\\
	& \quad\quad\quad\quad\quad1\le r_{1},r_{2}\le r+1,\,\quad0\le s_{1},s_{2}\le s+1,\,\quad0\le m_{1},m_{2}\le m+1.
	\end{split} \end{equation*}
\end{thm}

\begin{proof} First, note that the solution to \eqref{eq:cts1-1}$\--$\eqref{eq:cts5-1} satisfies for $1\le i \le N$,
	\begin{align}
	& \inp[A(\sigma + \alpha p I)]{\tau}_{\Omega_i} + \inp[u]{\dvr{\tau}}_{\Omega_i}+\inp[\gamma]{\tau}_{\Omega_i}-\gnp[u]{\tau\,n_i}_{\Gamma_i}=0, & \forall\tau\in\mathbb{X}_i,\label{eq:cts1-sub}\\
	& \inp[\dvr{\sigma}]{v}_{\Omega_i}=-\inp[f]{v}_{\Omega_i}, & \forall v\in V_i,\label{eq:cts2-sub}\\
	& \inp[\sigma]{\xi}_{\Omega_i}=0, & \forall\xi\in\mathbb{Q}_i,\label{eq:cts3-sub}\\
	& \inp[{K^{-1}} z]{\zeta}_{\Omega_i}-\inp[p]{\dvr{\zeta}}_{\Omega_i}+\gnp[p]{\zeta\cdot n_i}_{\Gamma_i}=0, & \forall \zeta\in Z_i,\label{eq:cts4-sub}\\
	& c_{0}\inp[\partial_t{p}]{w}_{\Omega_i}+\alpha\inp[\partial_t A(\sigma + \alpha pI)]{wI}_{\Omega_i}
	+\inp[\dvr{z}]{w}_{\Omega_i}=\inp[g]{w}_{\Omega_i}, & \forall w\in W_i,\label{eq:cts5-sub}
	\end{align}
	Subtracting (\ref{eq:fe-mono-weak-1})$\--$(\ref{eq:fe-mono-weak-5})
	from (\ref{eq:cts1-sub})$\--$(\ref{eq:cts5-sub}) gives
	\begin{align}
	& \inp[A\left(\left(\sigma-\sigma_{h}\right)+\alpha\left(p-p_{h}\right)I\right)]{\tau}+\sum_{i=1}^{N}\inp[\left(u-u_{h}\right)]{\dvr{\tau}}_{\Omega_{i}} \nonumber \\
	&\qquad \hspace{3pt} +\inp[\left(\gamma-\gamma_{h}\right)]{\tau}=\sum_{i=1}^{N}\gnp[u]{\tau n_i}_{\Gamma_{i}},&  \forall\tau\in\mathbb{X}_{h,0},\label{eq:error-1}\\
	& \sum_{i=1}^{N}\inp[\dvr{\left(\sigma-\sigma_{h}\right)}]{v}_{\Omega_{i}}=0, &  \forall v\in V_{h},\label{eq:error-2}\\
	& \inp[\left(\sigma-\sigma_{h}\right)]{\xi}=0, & \forall\xi\in\mathbb{Q}_{h},\label{eq:error-3}\\
	& \inp[{K^{-1}}\left(z-z_{h}\right)]{\zeta}-\sum_{i=1}^{N}\inp[\left(p-p_{h}\right)]{\dvr{\zeta}}_{\Omega_{i}}=-\sum_{i=1}^{N}\gnp[p]{\zeta\cdot n_{i}}_{\Gamma_{i}}, &  \forall\zeta\in Z_{h,0},\label{eq:error-4}\\
	& c_{0}\inp[\partial_t{\left(p-p_{h}\right)}]{w}+\alpha\inp[\partial_t A\left(\left(\sigma-\sigma_{h}\right)+\alpha\left(p-p_{h}\right)I\right)]{wI}\nonumber \\
	&+\sum_{i=1}^{N}\inp[\dvr{\left(z-z_{h}\right)}]{w}{}_{\Omega_{i}}=0, &  \forall w\in W_{h}.\label{eq:error-5}
	\end{align}
	Next, rewrite the above error equations in terms of the approximation
	errors $\psi_{\star}$ and discretization errors $\phi_{\star}$, for $\star\in\{\sigma,u,\gamma,z,p\}$
	as follows:
\begin{align*}
	\sigma-\sigma_{h} & =(\sigma-\Pi_{0}^{\sigma}\sigma)+(\Pi_{0}^{\sigma}\sigma-\sigma_{h}):=\psi_{\sigma}+\phi_{\sigma},\\
	u-u_{h} & =(u-{\cal P}_{h}^{u}u)+({\cal P}_{h}^{u}u-u_{h}):=\psi_{u}+\phi_{u},\\
	\gamma-\gamma_{h} & =(\gamma-{\cal R}_{h}\gamma)+({\cal R}_{h}\gamma-\gamma_{h}):=\psi_{\gamma}+\phi_{\gamma},\\
	z-z_{h} & =(z-\Pi_{0}^{z}z)+(\Pi_{0}^{z}z-z_{h}):=\psi_{z}+\phi_{z},\\
	p-p_{h} & =(p-{\cal P}_{h}^{p}p)+({\cal P}_{h}^{p}p-p_{h}):=\psi_{p}+\phi_{p},\\
	u-\lambda_{H}^u & =(u-\Qc^{u}_h u)+(\Qc^{u}_h u-\lambda_{H}^u):=\psi_{\lambda^u}+\phi_{\lambda^u},\\
	p-\lambda_{H}^p & =(p-\Qc^{p}_h p)+(\Qc^{p}_h p-\lambda_{H}^p):=\psi_{\lambda^p}+\phi_{\lambda^p}.
\end{align*}
Combining (\ref{eq:error-2}) with (\ref{eq:pi-weak-1}) gives
\begin{equation}
  \text{div}_h\phi_{\sigma}=0,
  \label{eq:err-proof-1}
\end{equation}
and (\ref{eq:error-3}) combined with (\ref{eq:pi-weak-2}) gives
	\begin{equation}
	\inp[\phi_{\sigma}]{\xi}=0\,\,\,\text{for }\xi\in\mathbb{Q}_{h}.\label{eq:err-proof-2}
	\end{equation}
	We rewrite error equation (\ref{eq:error-1}) as
	\begin{align}
	  & \inp[A\left(\phi_{\sigma} + \alpha\phi_{p}I\right)]{\tau}
          + \sum_{i=1}^{N}\inp[\phi_u]{\dvr{\tau}}_{\Omega_{i}}
          + \inp[\phi_{\gamma}]{\tau} \nonumber \\
          & \quad = -\inp[A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{\tau}-\inp[\psi_{\gamma}]{\tau}
          +\sum_{i=1}^{N}\gnp[u-{\cal I}_{H}^{u}u]{\tau n_i}_{\Gamma_{i}},
          \label{eq:err-proof-2-1}
	\end{align}
	where we have used
	\begin{equation}
	\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{u}u]{\tau n_i}_{\Gamma_{i}}=0,\label{eq:cont_lagrang_u}
	\end{equation}
	for any $\tau\in\mathbb{X}_{h,0}.$ Differentiating the above equation with
	respect to time gives
	\begin{align}
	& \inp[\partial_{t}A\left(\phi_{\sigma}+\alpha\phi_{p}I\right)]{\tau}+\sum_{i=1}^{N}\inp[\partial_{t}\phi_u]{\dvr{\tau}}_{\Omega_{i}}+\inp[\phi_{\gamma}]{\tau}\nonumber \\
	& \quad \quad =-\inp[\partial_{t}A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{\tau} -\inp[\partial_{t}\psi_{\gamma}]{\tau}+\sum_{i=1}^{N}\gnp[\partial_{t}\left(u-{\cal I}_{H}^{u}u\right)]{\tau n_i}_{\Gamma_{i}}.\label{eq:err-proof-3}
	\end{align}

	Taking $\tau=\phi_{\sigma}$ in (\ref{eq:err-proof-3}) and using (\ref{eq:err-proof-1})
	and (\ref{eq:err-proof-2}) gives
\begin{align}
	& \inp[\partial_{t}A\left(\phi_{\sigma}+\alpha\phi_{p}I\right)]{\phi_{\sigma}}=-\inp[\partial_{t}A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{\phi_{\sigma}}-\inp[\partial_{t}\psi_{\gamma}]{\phi_{\sigma}}+\sum_{i=1}^{N}\inp[\partial_{t}\left(u-{\cal I}_{H}^{u}u\right)]{\phi_{\sigma}n_i}_{\Gamma_{i}}.\label{eq:err-proof-4}
\end{align}
Error equation (\ref{eq:error-5}) can be written as
\begin{align}
  c_{0}\inp[\partial_{t}\phi_{p}]{w} + \alpha\inp[\partial_{t}A\left(\phi_{\sigma} + \alpha\phi_{p}I\right)]{wI} + \sum_{i=1}^{N}\inp[\dvr{\phi_{z}}]{w}{}_{\Omega_{i}} = -\alpha\inp[\partial_{t}A\left(\psi_{\sigma} + \alpha\psi_{p}I\right)]{wI},
 \label{eq:err-proof-5}
\end{align}
where we have used \eqref{eq:pressure-project} and \eqref{eq:pi-weak-d-1}.
Taking $w=\phi_{p}$ in (\ref{eq:err-proof-5}) and combining the
resulting equation with (\ref{eq:err-proof-4}) give
\begin{align}
  & \frac{1}{2}\partial_{t}\left(\|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|^{2}+c_{0}\|\phi_{p}\|^{2}\right)+\sum_{i=1}^{N}\inp[\dvr{\phi_{z}}]{\phi_{p}}{}_{\Omega_{i}}\nonumber \\
 & \quad =-\inp[\partial_{t}A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{\phi_{\sigma}+\alpha\phi_{p}I}
-\inp[\partial_{t}\psi_{\gamma}]{\phi_{\sigma}}+\sum_{i=1}^{N}\gnp[\partial_{t}\left(u-{\cal I}_{H}^{u}u\right)]{\phi_{\sigma}n_i}_{\Gamma_{i}}.\label{eq:err-proof-6}
\end{align}
Finally error equation (\ref{eq:error-4}) can be written as
	\begin{align}
	& \inp[{K^{-1}}\phi_{z}]{\zeta}-\sum_{i=1}^{N}\inp[\phi_{p}]{\dvr{\zeta}}_{\Omega_{i}}=-\inp[{K^{-1}}\psi_{z}]{\zeta}+\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}p-p]{\zeta\cdot n_{i}}_{\Gamma_{i}},\label{eq:err-proof-8}
	\end{align}
	where we have used for $\zeta \in Z_{h,0}$,
	\begin{equation}
	\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}p]{\zeta\cdot n_{i}}_{\Gamma_{i}}=0,\label{eq:cont_lagrang_p}
	\end{equation}
	and the definition of the $L^{2}(\text{\ensuremath{\Omega})}$-projection,
	${\cal P}_{h}^{p}$ onto space $W_{h}$.

	Taking test function $\zeta=\phi_{z}$ in
	equation (\ref{eq:err-proof-8})
	and combining the resulting equation with (\ref{eq:err-proof-6})
	gives
	\begin{align}
	& \frac{1}{2}\partial_{t}\left(\|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|^{2}+c_{0}\|\phi_{p}\|^{2}\right)+\|K^{-\frac{1}{2}}\phi_{z}\|^{2}=-\inp[\partial_{t}A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{\phi_{\sigma}+\alpha\phi_{p}I}\nonumber \\
	  & \quad\quad\quad\quad-\inp[\partial_{t}\psi_{\gamma}]{\phi_{\sigma}}-\inp[{K^{-1}}\psi_{z}]{\phi_{z}}+\inp[\psi_{\sigma}]{\partial_{t}\phi_{\gamma}}
          -\sum_{i=1}^{N}\gnp[\partial_{t}({\cal I}_{H}^{u}u-u)]{\phi_{\sigma}n_i}_{\Gamma_{i}}\nonumber \\
	& \quad\quad\quad\quad+\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}p-p]{\phi_{z}\cdot n_{i}}_{\Gamma_{i}},\label{eq:err-proof-9}
	\end{align}
	where we have used (\ref{eq:err-proof-2}) to conclude that $\inp[\phi_{\sigma}]{\partial_{t}\phi_{\gamma}}=0.$

	Next, we bound the first three terms on the right hand side of (\ref{eq:err-proof-9}).
	\begin{align}
	& -\inp[\partial_{t}A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{\phi_{\sigma}+\alpha\phi_{p}I}-\inp[\partial_{t}\psi_{\gamma}]{\phi_{\sigma}}-\inp[{K^{-1}}\psi_{z}]{\phi_{z}}\nonumber \\
	& \quad\quad\quad\quad\le\|\partial_{t}A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)\|\|\phi_{\sigma}+\alpha\phi_{p}I\|+\|\partial_{t}\psi_{\gamma}\|\|\phi_{\sigma}\|+\|{K^{-1}}\psi_{z}\|\|\phi_{z}\|\nonumber \\
	& \quad\quad\quad\quad\le \frac{C}{\epsilon}\left(\|\partial_{t}\psi_{\sigma}\|^{2}+\|\partial_{t}\psi_{p}\|^{2}+\|\partial_{t}\psi_{\gamma}\|^{2}+\|\psi_{z}\|^{2}\right)+\epsilon\left(\|\phi_{\sigma}\|^{2}+\|\phi_{p}\|^{2}+\|\phi_{z}\|^{2}\right),\label{eq:err-proof-10}
	\end{align}
	where we have used the operator bounds (\ref{eq:coercivity-elast-1})
	and (\ref{eq:coercivity-flow-1}) along with Young's inequality for
	a some $\epsilon > 0$.

	Next, we give a bound on the last two boundary terms in the right
	hand side of equation (\ref{eq:err-proof-9}). For this, we note that
	the following bounds hold for any $(\tau,v)\in H(\text{div};\Omega,\mathbb M) \times H^1_0(\Omega,\mathbb R^d)$
        and $(\zeta,w)\in H(\text{div};\Omega) \times H^1_0(\Omega)$:
	\begin{align}
	\gnp[{\cal I}_{H}^{u}v-v]{\tau n_i}_{\Gamma_{i}} & = \gnp[E_{i}({\cal I}_{H}^{u}v-v)]{\tau n_{i}}_{{\partial\Omega_i}} \nonumber\\
& \le C\|E_{i}({\cal I}_{H}^{u}v-v)\|_{\frac{1}{2},{\partial\Omega_i}}\|\tau\|_{H(\text{div};\Omega_i)}
        \le C\|{\cal I}_{H}^{u}v-v\|_{\frac{1}{2},\Gamma _{i}}\|\tau\|_{H(\text{div};\Omega_i)},
        \label{eq:err-proof-10-1}
        \end{align}
\begin{align}
        \gnp[{\cal I}_{H}^{p}w-w]{\zeta \cdot n_i}_{\Gamma_{i}} &
        = \gnp[E_{i}({\cal I}_{H}^{p}\zeta-\zeta)]{\zeta \cdot n_{i}}_{{\partial\Omega_i}} \nonumber \\
	& \le C\|E_{i}({\cal I}_{H}^{p}\zeta-\zeta)\|_{\frac{1}{2},{\partial\Omega_i}}\|\zeta\|_{H(\text{div};\Omega_i)}
        \le C\|{\cal I}_{H}^{p}\zeta-\zeta\|_{\frac{1}{2},\Gamma _{i}}\|\zeta\|_{H(\text{div};\Omega_i)},\label{eq:err-proof-10-2}
\end{align}
where $E_{i}$ denotes the extension by zero from $\Gamma_i$ to
${\partial\Omega_i}$, which is continuous in the $H^{\frac12}$-norm for functions that are zero on $\partial\Gamma$,
and we have used the trace inequality (\ref{eq:trace-2}). Taking $(\tau,v)=( \phi_{\sigma},\partial_t u)$ and $(\zeta,w)=( \phi_{z},p)$ in (\ref{eq:err-proof-10-1}) and (\ref{eq:err-proof-10-2}), respectively, and using Young's inequality gives
\begin{align}
  & \gnp[\partial_{t}({\cal I}_{H}^{u}u-u)]{\phi_{\sigma}n_i}_{\Gamma_{i}}
  \le \frac{C}{\epsilon}\|\partial_{t}({\cal I}_{H}^{u}u-u)\|_{\frac{1}{2},\Gamma _{i}}^2 + \epsilon \|\phi_{\sigma}\|_{\Omega_i}^2, \label{bound-gamma-1}\\
  & \gnp[{\cal I}_{H}^{p}p-p]{\phi_{z}\cdot n_{i}}_{\Gamma_{i}}
  \le \frac{C}{\epsilon}\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma _{i}}^2
  + \epsilon(\|\phi_{z}\|_{\Omega_i}^2 + \|\dvr \phi_z\|_{\Omega_i}^2),
  \label{bound-gamma-2}
\end{align}
where we also used \eqref{eq:err-proof-1}. Combining inequalities \eqref{eq:err-proof-9}--\eqref{bound-gamma-1} respect to time from $0$ to $t\in(0,T]$ gives
	\begin{align}
	& \|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)(t)\|^{2}+c_{0}\|\phi_{p}(t)\|^{2}+\int_{0}^{t}\|{K^{-1}}\phi_{z}\|^{2}\nonumber \\
	& \quad \le C\int_{0}^{t}(\|\partial_{t}\psi_{\sigma}\|^{2}+\|\partial_{t}\psi_{p}\|^{2}+\|\partial_{t}\psi_{\gamma}\|^{2}+\|\psi_{z}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t u-\partial_t u\|_{\frac{1}{2},\Gamma _{i}}^{2}\nonumber \\
	  & \qquad +\sum_{i=1}^{N}\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma _{i}}^{2})ds+\epsilon\int_{0}^{t}\left(\|\phi_{\sigma}\|^{2}+\|\phi_{p}\|^{2}+\|\phi_{z}\|^{2}\right)ds \nonumber \\
          & \qquad + C\int_{0}^{t}\|\text{div}_h\phi_{z}\|^2 ds
	+ \int_{0}^{t}\inp[\psi_{\sigma}]{\partial_{t}\phi_{\gamma}}ds +\|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)(0)\|^{2}+c_{0}\|\phi_{p}(0)\|^{2},\label{eq:err-proof-13}
	\end{align}

	To bound the term $\displaystyle\int_{0}^{t}\inp[\psi_{\sigma}]{\partial_{t}\phi_{\gamma}}ds$,
	we use integration by parts as follows:
	\begin{align}
	& \int_{0}^{t}\inp[\psi_{\sigma}]{\partial_{t}\phi_{\gamma}}ds=\inp[\psi_{\sigma}(t)]{\phi_{\gamma}(t)}-\inp[\psi_{\sigma}(0)]{\phi_{\gamma}(0)}-\int_{0}^{t}\inp[\partial_{t}\psi_{\sigma}]{\phi_{\gamma}}ds \nonumber \\
	& \le \frac{C}{\epsilon}\left(\int_{0}^{t}\|\partial_{t}\psi_{\sigma}\|^{2}ds+\|\psi_{\sigma}(t)\|^{2}\right)+\epsilon\left(\int_{0}^{t}\|\phi_{\gamma}\|^{2}ds+\|\phi_{\gamma}(t)\|^{2}\right)+C\left(\|\psi_{\sigma}(0)\|^{2}+\|\phi_{\gamma}(0)\|^{2}\right),\label{eq-proof-13-0-1}
	\end{align}
	where we have used the Cauchy-Schwartz and Young's inequality.

	Next, we bound the errors of the form $\phi_{\star}$ for
	$\star\in\{\sigma,\gamma,u,p\}$.

	Using the inf-sup condition (\ref{eq:inf-sup-elast-weak}) and the error equation
	(\ref{eq:err-proof-1}) gives
	\begin{align}
	& \|\phi_{u}\|+\|\phi_{\gamma}\|\le C\sup_{0\ne\tau\in\mathbb{X}_{h,0}}\frac{\sum_{i=1}^{N}\inp[\phi_{u}]{\dvr{\tau}}_{\Omega_{i}}+\inp[\phi_{\gamma}]{\tau}}{\|\tau\|_{\mathbb{X}_{h}}}\nonumber \\
	& \quad\le C\sup_{0\ne\tau\in\mathbb{X}_{h,0}}\frac{1}{\|\tau\|_{\mathbb{X}_{h}}}\Big|\inp[A\left(\phi_{\sigma}+\alpha\phi_{p}I\right)]{\tau}+\inp[A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{\tau}+\inp[\psi_{\gamma}]{\tau}\nonumber \\
	& \quad-\sum_{i=1}^{N}\inp[\left({\cal I}_{H}^{u}u-u\right)]{\tau n_i}_{\Gamma_{i}}\Big|\nonumber \\
	& \quad\le C\left(\|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|+\|\psi_{\sigma}\|+\|\psi_{\gamma}\|+\|\psi_{p}\|+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}\right),\label{eq:err-proof-14}
	\end{align}
	where we have used Cauchy-Schwartz inequality, (\ref{eq:err-proof-1})
	and (\ref{eq:err-proof-10-1}) with $v=u$ to get last inequality.
	Integrating (\ref{eq:err-proof-14}) with respect to time from $0$
	to $t\in(0,T]$ gives
	\begin{align}
	& \int_{0}^{t}\left(\|\phi_{u}\|^{2}+\|\phi_{\gamma}\|^{2}\right)ds\le C\int_{0}^{t}\left(\|\phi_{\sigma}\|^{2}+\alpha\|\phi_{p}I\|^{2}\right)ds\nonumber \\
	& \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+C\int_{0}^{t}\left(\|\psi_{\sigma}\|^{2}+\|\psi_{\gamma}\|^{2}+\|\psi_{p}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}^{2}\right)ds.\label{eq:err-proof-15}
	\end{align}

	To bound $\|\phi_{p}\|$, we use the inf-sup stability condition (\ref{eq:inf-sup-darcy-weak})
	to get
	\begin{align}
	& \|\phi_{p}\|\le C\sup_{0\ne\zeta\in Z_{h,0}}\frac{\sum_{i=1}^{N}\left(\dvr\zeta,\phi_{p}\right)_{\Omega_{i}}}{\|\zeta\|_{Z_{h}}}\nonumber \\
	& \quad\quad\le C\sup_{0\ne\zeta\in Z_{h,0}}\frac{\inp[{K^{-1}}\phi_{z}]{\zeta}+\inp[{K^{-1}}\psi_{z}]{\zeta}-\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}p-p]{\zeta\cdot n_{i}}_{\Gamma_{i}}}{||\zeta||_{Z_{h}}}\nonumber \\
	& \quad\quad\le C\left(\|\psi_{z}\|+\|K^{-\frac{1}{2}}\phi_{z}\|+\sum_{i=1}^{N}\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma _{i}}\right),\label{eq:err-proof-16}
	\end{align}
	where we have used (\ref{eq:err-proof-8}) and (\ref{eq:cont_lagrang_p})
	to obtain second inequality and (\ref{eq:err-proof-10-2}) with $w=p$
	to obtain the last one. Integrating (\ref{eq:err-proof-16}) in time
	from $0$ to $t\in(0,T]$ yields
	\begin{equation}
	\int_{0}^{t}\|\phi_{p}\|^{2}ds\le C\int_{0}^{t}\left(\|\psi_{z}\|^{2}+\|K^{-\frac{1}{2}}\phi_{z}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma _{i}}^{2}\right)ds.\label{eq:err-proof-17}
	\end{equation}

	To bound the term $\displaystyle\int_{0}^{t}\| \phi_{\sigma}\|^{2}ds$, we take $\tau= \phi_{\sigma}$ in (\ref{eq:err-proof-2-1})
	and $\xi= \phi_{\gamma}$ in (\ref{eq:error-3}), and use \eqref{eq:err-proof-1}$\--$\eqref{eq:err-proof-2} and Young's inequality
to get
	\begin{align}
	& \|A^{\frac{1}{2}} \phi_{\sigma}\|^{2}=-\inp[A^{\frac{1}{2}}\alpha\phi_{p}I]{ \phi_{\sigma}}-\inp[A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{ \phi_{\sigma}}-\inp[\psi_{\gamma}]{ \phi_{\sigma}}-\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{u}u-u]{ \phi_{\sigma} n_i}_{\Gamma_{i}} \nonumber \\
	& +\inp[\psi_{\sigma}]{ \phi_{\gamma}}\le C \Big( \big(\| \phi_{p}\|+\| \psi_{\sigma}\|+\| \psi_{p}\|+\|  \psi_{\gamma} \|\big)\|  \phi_{\sigma} \| +\sum_{i=1}^{N}\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}\|  \phi_{\sigma}\|_{H(\dvr;\Omega_i)} \nonumber \\
	& \hspace{40pt}+\| \psi_{\sigma}\|\| \phi_{\gamma}\|\Big)\le \frac{C}{\epsilon}\left(\| \psi_{\sigma}\|^{2}+\| \psi_{p}\|^{2}+\| \psi_{\gamma}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}^{2}+\| \phi_{p}\|^{2}\right)\nonumber \\
	& \hspace{150pt}+\epsilon\left(\| \phi_{\sigma}\|^{2}+\| \phi_{\gamma}\|^{2}\right),\label{eq:err-proof-18}
	\end{align}
	where we have used (\ref{eq:err-proof-10-1}) with $(\tau,v)=( \phi_{\sigma},u)$
	to arrive at the first inequality and (\ref{eq:coercivity-elast-1})
	with Young's inequality to justify the last inequality. Now again
	using (\ref{eq:coercivity-elast-1}), integrating (\ref{eq:err-proof-18})
	with respect to time from $0$ to $t\in(0,T]$, and taking $\epsilon$ small enough, we get
	\begin{align}
	& \int_{0}^{t}\| \phi_{\sigma}\|^{2}ds\le C\Big(\int_{0}^{t}\big(\| \psi_{\sigma}\|^{2}+\| \psi_{p}\|^{2}+\| \psi_{\gamma}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}^{2} + \| \phi_{p}\|^{2}\big)ds\Big)\nonumber \\
	& \quad\quad\quad\quad\quad\quad+\epsilon\int_{0}^{t}\| \phi_{\gamma}\|^{2}ds.\label{eq:err-proof-19}
	\end{align}

	Combining (\ref{eq:err-proof-13})$\--$(\ref{eq:err-proof-19}),
	using (\ref{eq:err-proof-1}), and taking $\epsilon$ small enough gives
	\begin{align}
	  & \|A^{\frac{1}{2}}\left( \phi_{\sigma}+\alpha \phi_{p} I\right)\|^{2}+\| \phi_{u}\|^{2}+\| \phi_{\gamma}\|^{2}+\|c_{0}^{\frac{1}{2}} \phi_{p}\|^{2}
          +\|\dvr_h \phi_{\sigma}\|^{2}\nonumber \\
	  & +\int_{0}^{t}\left(\| \phi_{\sigma}\|^{2}+\| \phi_{u}\|^{2}+\| \phi_{\gamma}\|^{2}+\|K^{-\frac{1}{2}} \phi_{z}\|^{2}
          +\| \phi_{p}\|^{2}+\|\dvr_h \phi_{\sigma}\|^{2}\right)\nonumber \\
	& \le C\int_{0}^{t}\left(\|\partial_t \psi_{\sigma}\|^{2}+\|\partial_t \psi_{p}\|^{2}+\|\partial_t \psi_{\gamma}\|^{2}+\| \psi_{\sigma}\|^{2}+\| \psi_{p}\|^{2}+\| \psi_{\gamma}\|^{2}+\| \psi_{z}\|^{2}\right)ds+C(\| \psi_{\sigma}\|^{2}\nonumber \\
	& \quad+\| \psi_{p}\|^{2}+\| \psi_{\gamma}\|^{2})+C\int_{0}^{t}\|\text{div}_h\phi_{z}\|^{2}ds
	+C\sum_{i=1}^{N}\left({\cal I}_{H}^{u}u-u(t)\|_{\frac{1}{2},\Gamma _{i}}^{2}+\|{\cal I}_{H}^{p}p-p(t)\|_{\frac{1}{2},\Gamma _{i}}^{2}\right)\nonumber \\
	& \quad+C\sum_{i=1}^{N}\int_{0}^{t}\left(\|{\cal I}_{H}^{u}\partial_t u-\partial_t u\|_{\frac{1}{2},\Gamma _{i}}^{2}+\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}^{2}+\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma _{i}}^{2}\right)ds\nonumber \\
	& \quad+C\left(\| \phi_{\sigma}(0)\|^{2}+\| \phi_{p}(0)\|^{2}+\| \phi_{\gamma}(0)\|^{2}\right).\label{eq:err-proof-19-1}
	\end{align}
	\textbf{Bound on $\sum_{i=1}^{N}\|\dvr\phi_{z}\|_{\Omega_i}$.}

	Take $w=\phi_{z}$ in (\ref{eq:error-5}) to get, for  $i=1,2,\ldots,N,$
	\begin{align*}
	& \|\dvr\phi_{z}\|_{\Omega_i}^{2}=-\inp[c_{0}\partial_t \phi_{p}]{\dvr \phi_{z}}_{\Omega_{i}}-\inp[c_{0}\partial_t \phi_{p}]{\dvr \phi_{z}}_{\Omega_{i}}-\alpha\inp[\partial_t A\left( \phi_{\sigma}+\alpha \phi_{p} I\right)]{\dvr \phi_{z} I}_{\Omega_{i}}\\
	& -\alpha\inp[\partial_t A\left( \psi_{\sigma}+\alpha \psi_{p} I\right)]{\dvr \phi_{z} I}_{\Omega_{i}}+\inp[ \psi_{z}]{\dvr \phi_{z}}{}_{\Omega_{i}}=-\inp[c_{0}\partial_t \phi_{p}]{\dvr \phi_{z}}_{\Omega_{i}}\\
	& -\alpha\inp[\partial_t A\left( \phi_{\sigma}+\alpha \phi_{p} I\right)]{\dvr \phi_{z} I}_{\Omega_{i}}-\alpha\inp[\partial_t A\left( \psi_{\sigma}+\alpha \psi_{p} I\right)]{\dvr \phi_{z} I}_{\Omega_{i}},
	\end{align*}
	where the last equality follows from (\ref{eq:pressure-project})
	and (\ref{eq:pi-weak-d-1}). Finally, summing the above equation over
	all the subdomain indices and using (\ref{eq:coercivity-elast-1}),
	we get, for $i=1,2,\ldots,N,$
	\begin{align}
	& \|\dvr\phi_{z}\|_{\Omega_i}\le C\left(\|c_{0}^{\frac{1}{2}}\partial_t \phi_{p}\|_{\Omega_i}+\|\partial_t A^{\frac{1}{2}}\left( \phi_{\sigma}+\alpha \phi_{p} I\right)\|_{\Omega_i}+\| \psi_{p}\|_{\Omega_i}+\| \psi_{\sigma}\|_{\Omega_i}\right).\label{eq:err-proof-20}
	\end{align}
	Squaring the above equation and summing over all the subdomain indices
	give
	\begin{align}
	& \sum_{i=1}^{N}\|\dvr\phi_{z}\|_{\Omega_i}^{2}\le C\left(\|c_{0}^{\frac{1}{2}}\partial_t \phi_{p}\|^{2}+\|\partial_t A^{\frac{1}{2}}\left( \phi_{\sigma}+\alpha \phi_{p} I\right)\|^{2}+\| \psi_{p}\|^{2}+\| \psi_{\sigma}\|^{2}\right).\label{eq:err-proof-21}
	\end{align}

	In order to bound $\|c_{0}^{\frac{1}{2}}\partial_t \phi_{p}\|^{2}$ and $\|\partial_t A^{\frac{1}{2}}\left( \phi_{\sigma}+\alpha \phi_{p} I\right)\|^{2}$,
	we use arguments similar to the ones used in the stability analysis
	for the method. Differentiate (\ref{eq:err-proof-1}), (\ref{eq:err-proof-2}),
	and (\ref{eq:err-proof-8}) in time, combine them with (\ref{eq:err-proof-3}),
	and take test functions $\tau=\partial_t \phi_{\sigma}$, $\xi=\partial_t \phi_{\gamma}$, $q=\phi_{z}$,
	and $w=\partial_t \phi_{p}$ to get the following time differentiated version of
	(\ref{eq:err-proof-9}):
	\begin{align}
	& \|\partial_{t}A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|^{2}+\|c_{0}^{\frac{1}{2}}\partial_t\phi_{p}\|^{2}+\frac{1}{2}\partial_t\|K^{-\frac{1}{2}}\phi_{z}\|^{2}=-\inp[\partial_{t}A\left(\psi_{\sigma}+\alpha\psi_{p}I\right)]{\partial_t\left(\phi_{\sigma}+\alpha\phi_{p}I\right)}\nonumber \\
	& \quad-\inp[\partial_{t}\psi_{\gamma}]{\partial_t\left(\phi_{\sigma}+\alpha\phi_{p}I\right)}-\inp[\partial_t{K^{-1}}\psi_{z}]{\phi_{z}}+\inp[\partial_t\psi_{\sigma}]{\partial_{t}\phi_{\gamma}}-\sum_{i=1}^{N}\gnp[\partial_{t}\left({\cal I}_{H}^{u}u-u\right)]{\partial_t\phi_{\sigma}n_i}_{\Gamma_{i}}\nonumber \\
	& \quad+\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}\partial_t p-\partial_t p]{\phi_{z}\cdot n_{i}}_{\Gamma_{i}},\label{eq:err-proof-22}
	\end{align}
	where we have used the fact that $\inp[\partial_{t}\psi_{\gamma}]{\partial_t\alpha\phi_{p}I}=0$
	to write \begin{equation*}
	\inp[\partial_{t}\psi_{\gamma}]{\partial_t\phi_{\sigma}}=\inp[\partial_{t}\psi_{\gamma}]{\partial_t\left(\phi_{\sigma}+\alpha\phi_{p}I\right)}.
	\end{equation*}\

	To bound $\sum_{i=1}^{N}\gnp[\partial_{t}\left({\cal I}_{H}^{u}u-u\right)]{\partial_t\phi_{\sigma}n_i}_{\Gamma_{i}}$,
	we use integration by parts to rewrite
	\begin{align}
	& \sum_{i=1}^{N}\gnp[\partial_{t}\left({\cal I}_{H}^{u}u-u\right)]{\partial_t\phi_{\sigma}n_i}_{\Gamma_{i}}=\ddt\left(\sum_{i=1}^{N}\gnp[\partial_{t}\left({\cal I}_{H}^{u}u-u\right)]{\phi_{\sigma}n_i}_{\Gamma_{i}}\right) \nonumber \\
	& -\sum_{i=1}^{N}\gnp[\partial_{t}^{2}\left({\cal I}_{H}^{u}u-u\right)]{\phi_{\sigma}n_i}_{\Gamma_{i}}.\label{eq:err-proof-23}
	\end{align}
	To bound the last term on the right hand side of the above equation,
	we take $(\tau,v)=( \phi_{\sigma},\partial_t^{2}u)$ in (\ref{eq:err-proof-10-1}) to
	get
	\begin{equation}
	\Big|\sum_{i=1}^{N}\gnp[\partial_{t}^{2}\left({\cal I}_{H}^{u}u-u\right)]{\phi_{\sigma}n_i}_{\Gamma_{i}}\Big|\le C\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t^{2}u-\partial_t^{2}u\|_{\frac{1}{2},\Gamma _{i}}\|  \phi_{\sigma}\|_{H(\dvr;\Omega_i)}.\label{eq:err-proof-24}
	\end{equation}
	To bound $\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}\partial_t p-\partial_t p]{\partial_t\phi_{z}\cdot n_{i}}_{\Gamma_{i}}$,
	we take $(\zeta,w)=( \phi_{z},\partial_t p)$ in (\ref{eq:err-proof-10-2}), to
	get

	\begin{align}
	& \Big|\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}\partial_t p-\partial_t p]{ \phi_{z}\cdot n_{i}}_{\Gamma_{i}}\Big|\le C\sum_{i=1}^{N}\|{\cal I}_{H}^{p}\partial_t p-\partial_t p\|_{\frac{1}{2},\Gamma _{i}}\|  \phi_{z}\|_{H(\dvr;\Omega_i)}.\label{eq:err-proof-25}
	\end{align}

	Combining (\ref{eq:err-proof-22})$\--$(\ref{eq:err-proof-25}),
	using (\ref{eq:coercivity-elast-1})$\--$(\ref{eq:coercivity-flow-1})
	and integrating with respect to time from $0$ to $t\in(0,T]$, we
	get
	\begin{align}
	& \|K^{-\frac{1}{2}}\phi_{z}\|^{2}+\int_{0}^{t}\left(\|\partial_{t}A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|^{2}+\|c_{0}^{\frac{1}{2}}\partial_t\phi_{p}\|^{2}\right)ds\nonumber \\
	& \le C\int_{0}^{t}(\|\partial_t \psi_{\sigma}\|^{2}+\|\partial_t \psi_{p}\|^{2}+\|\partial_t \psi_{\gamma}\|^{2}+\|\partial_t \psi_{z}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t^{2}u-\partial_t^{2}u\|_{\frac{1}{2},\Gamma _{i}}^{2}\nonumber \\
	& +\sum_{i=1}^{N}\|{\cal I}_{H}^{p}\partial_t p-\partial_t p\|_{\frac{1}{2},\Gamma _{i}}^{2})ds+C\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t u-\partial_t u(t)\|_{\frac{1}{2},\Gamma _{i}}^{2}\nonumber \\
	& +\epsilon\left(\int_{0}^{t}\left(\| \phi_{\sigma}\|^{2}+\|\partial_t \phi_{\gamma}\|^{2}+\| \phi_{z}\|^{2}+\sum_{i=1}^{N}\|\dvr \phi_{z}\|_{\Omega_i}^{2}\right)ds+\| \phi_{\sigma}(t)\|^{2}\right)\nonumber \\
	& +C\left(\|K^{-\frac{1}{2}}\phi_{z}(0)\|^{2}+\| \phi_{\sigma}(0)\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t u-\partial_t u(0)\|_{\frac{1}{2},\Gamma _{i}}^{2}\right),\label{eq:err-proof-26}
	\end{align}
	where we have used Young's inequality for $\epsilon>0$ and (\ref{eq:err-proof-10-1})
	with $(\tau,v)=( \phi_{\sigma},\partial_t u)$.

	Using the inf-sup condition \eqref{eq:inf-sup-elast-weak} with $v=\partial_t \phi_{u}$, $\xi=\partial_t  \phi_{\gamma}$, the time-differentiated \eqref{eq:err-proof-2-1}, and following
	the steps similar to the ones used to arrive at (\ref{eq:err-proof-15}),
	we get
	\begin{align}

	& \int_{0}^{t}\left(\|\partial_t\phi_{u}\|^{2}+\|\partial_t\phi_{\gamma}\|^{2}\right)ds\le C\Big(\int_{0}^{t}\|\partial_{t}A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|^{2}ds \nonumber \\
	& \quad \quad  + \int_{0}^{t}\big(\|\partial_t\psi_{\sigma}\|^{2}+\|\partial_t\psi_{\gamma}\|^{2}+\|\partial_t\psi_{p}\|^{2} +\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t u-\partial_t u(t)\|_{\frac{1}{2},\Gamma _{i}}^{2}\big)ds\Big).\label{eq:err-proof-27}

	\end{align}

	Combining (\ref{eq:err-proof-26})$\--$(\ref{eq:err-proof-27}), taking $\epsilon$ small enough, and
	using (\ref{eq:err-proof-16}) implies
	\begin{align}
	& \|K^{-\frac{1}{2}}\phi_{z}(t)\|^{2}+\| \phi_{p}(t)\|^{2}+\int_{0}^{t}\left(\|\partial_{t}A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|^{2}+\|c_{0}^{\frac{1}{2}}\partial_t\phi_{p}\|^{2}\right)ds\nonumber \\
	& \le C\int_{0}^{t}(\|\partial_t \psi_{\sigma}\|^{2}+\|\partial_t \psi_{p}\|^{2}+\|\partial_t \psi_{\gamma}\|^{2}+\|\partial_t \psi_{z}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t^{2}u-\partial_t^{2}u\|_{\frac{1}{2},\Gamma _{i}}^{2}\nonumber \\
	& \quad+\sum_{i=1}^{N}\|{\cal I}_{H}^{p}\partial_t p-\partial_t p\|_{\frac{1}{2},\Omega_i}^{2})ds+C\sum_{i=1}^{N}(\|{\cal I}_{H}^{u}\partial_t u-\partial_t u(t)\|_{\frac{1}{2},\Gamma _{i}}^{2}\nonumber \\
	& \quad+\|{\cal I}_{H}^{p}p-p(t)\|_{\frac{1}{2},\Gamma _{i}}^{2})+\epsilon\left(\int_{0}^{t}\left(\| \phi_{\sigma}\|^{2}+\| \phi_{z}\|^{2}+\sum_{i=1}^{N}\|\dvr \phi_{z}\|_{\Omega_i}^{2}\right)ds\right)\nonumber \\
	& \quad+\epsilon\left(\| \phi_{\sigma}(t)\|^{2}\right)+C\left(\| \psi_{z}(t)\|^{2}+\|K^{-\frac{1}{2}}\phi_{z}(0)\|^{2}+\| \phi_{\sigma}(0)\|^{2}\right)\nonumber \\
	& \quad+C\left(\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t u-\partial_t u(0)\|_{\frac{1}{2},\Gamma _{i}}^{2}\right).\label{eq:err-proof-28}
	\end{align}

	Combining (\ref{eq:err-proof-21}) and (\ref{eq:err-proof-28}) gives
	\begin{align}
	  & \|K^{-\frac{1}{2}}\phi_{z}(t)\|^{2}+\| \phi_{p}(t)\|^{2}
          +\int_{0}^{t}\|\dvr_h\phi_{z}\|^{2}ds\nonumber \\
	& \le C\int_{0}^{t}(\|\partial_t \psi_{\sigma}\|^{2}+\|\partial_t \psi_{p}\|^{2}+\|\partial_t \psi_{\gamma}\|^{2}+\|\partial_t \psi_{z}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t^{2}u-\partial_t^{2}u\|_{\frac{1}{2},\Gamma _{i}}^{2}\nonumber \\
	& \quad+\sum_{i=1}^{N}\|{\cal I}_{H}^{p}\partial_t p-\partial_t p\|_{\frac{1}{2},\Gamma _{i}}^{2})ds+C\sum_{i=1}^{N}\Big(\|{\cal I}_{H}^{u}\partial_t u-\partial_t u\|_{\frac{1}{2},\Gamma _{i}}^{2}\nonumber \\
	& \quad+\|{\cal I}_{H}^{p}p-p(t)\|_{\frac{1}{2},\Gamma _{i}}^{2}\Big)+\epsilon\Big(\int_{0}^{t}\big(\| \phi_{\sigma}\|^{2}+\| \phi_{z}\|^{2}+\sum_{i=1}^{N}\|\dvr \phi_{z}\|_{\Omega_i}^{2}\big)ds+\| \phi_{\sigma}(t)\|^{2}\Big)\nonumber \\
	& \quad+C\Big(\| \psi_{z}(t)\|^{2}+\|K^{-\frac{1}{2}}\phi_{z}(0)\|^{2}+\| \phi_{\sigma}(0)\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}\partial_t u-\partial_t u(0)\|_{\frac{1}{2},\Gamma _{i}}^{2}\Big).\label{eq:err-proof-29}
	\end{align}

	Coercivity of $A$, (\ref{eq:coercivity-elast-1}), also implies
	\begin{equation}
	\| \phi_{\sigma}\|\le C\left(\|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|+\| \phi_{p}\|\right).\label{eq:err-proof-30}
	\end{equation}

	\textbf{Bound on $\|\lambda_{H}^u\|_{\Gamma}$ and $\|\lambda_{H}^p\|_{\Gamma}$.}

	In order to bound $\|\lambda_{H}^u\|_{\Gamma}$, we take the difference between
	equations (\ref{eq:cts1-1}) and (\ref{eq:monodd-mmmfe1}) to get
	\begin{align*}
	& \inp[A\left(\left(\sigma-\sigma_{h}\right)+\alpha\left(p-p_{h}\right)I\right)]{\tau}+\sum_{i=1}^{N}\inp[\left(u-u_{h}\right)]{\dvr{\tau}}_{\Omega_{i}}+\inp[\left(\gamma-\gamma_{h}\right)]{\tau}\\
	& \qquad-\sum_{i=1}^{N}\inp[u]{\tau n_i}_{\Gamma_{i}}=\sum_{i=1}^{N}\gnp[u-\lambda_{H}^{u}]{\tau\,n_{i}}_{\Gamma_{i}}, & \forall\tau\in\mathbb{X}_{h}.
	\end{align*}
	We can split the error terms in the above equation and use (\ref{eq:motor-project-1})
	to rewrite the above equation as
	\begin{align}
	& \sum_{i=1}^{N}\gnp[\phi_{\lHu}]{\tau\,n_{i}}_{\Gamma_{i}}=\inp[A\left( \phi_{\sigma}+\alpha \phi_{p}\right)]{\tau}+\inp[A\left( \psi_{\sigma}+\alpha \psi_{p}\right)]{\tau}+\sum_{i=1}^{N}\inp[ \phi_{u}]{\dvr{\tau}}_{\Omega_{i}}\nonumber \\
	& +\sum_{i=1}^{N}\inp[ \psi_{u}]{\dvr{\tau}}_{\Omega_{i}}+\inp[ \phi_{\gamma}]{\tau}+\inp[ \psi_{\gamma}]{\tau}+\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{u}u-u]{\tau n_i}_{\Gamma_{i}}, & \forall\tau\in\mathbb{X}_{h}.\label{eq:err-proof-31}
	\end{align}

	Inf-sup stability bound (\ref{eq:mortar-u-inf-sup}) combined with
	(\ref{eq:err-proof-31}) implies
	\begin{align*}
	& \|\phi_{\lHu}\|_{\Gamma}\le C\sup_{0\ne\tau\in\mathbb{X}_{h}}\frac{\sum_{i=1}^{N}\gnp[\tau n_{i}]{\phi_{\lHu}}_{\Gamma_{i}}}{\|\tau\|_{\mathbb{X}_{h}}}=C\sup_{0\ne\tau\in\mathbb{X}_{h}}\Big(\frac{\inp[A\left( \phi_{\sigma}+\alpha \phi_{p}\right)]{\tau}+\inp[A\left( \psi_{\sigma}+\alpha \psi_{p}\right)]{\tau}}{\|\tau\|_{\mathbb{X}_{h}}}\\
	& +\frac{\sum_{i=1}^{N}\inp[ \phi_{u}]{\dvr{\tau}}_{\Omega_{i}}+\sum_{i=1}^{N}\inp[ \psi_{u}]{\dvr{\tau}}_{\Omega_{i}}+\inp[ \phi_{\gamma}]{\tau}+\inp[ \psi_{\gamma}]{\tau}+\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{u}u-u]{\tau n_i}_{\Gamma_{i}}}{\|\tau\|_{\mathbb{X}_{h}}}\Big)\\
	& \le C\big(\|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|+\| \phi_{u}\|+\| \phi_{\gamma}\|\big)+C\Big(\| \psi_{\sigma}\|+\| \psi_{p}\|+\| \psi_{u}\|+\| \psi_{\gamma}\|\\
	& \hspace{20pt}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}\Big),
	\end{align*}
	where in the last inequality, we have used (\ref{eq:coercivity-elast-1})
	and (\ref{eq:err-proof-10-1}) with $v=u$. Squaring the above inequality
	and then integrating with respect to time from $0$ to $t$ yields
	the following two bounds
	\begin{align}
	& \|\phi_{\lHu}\|_{\Gamma}^{2}\le C\left(\|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|^{2}+\| \phi_{u}\|^{2}+\| \phi_{\gamma}\|^{2}\right)\nonumber \\
	& \quad\quad\quad\quad\quad+C\left(\| \psi_{\sigma}\|^{2}+\| \psi_{p}\|^{2}+\| \psi_{u}\|^{2}+\| \psi_{\gamma}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}^{2}\right),\label{eq:err-proof-32}\\
	& \int_{0}^{t}\|\phi_{\lHu}\|_{\Gamma}^{2}ds\le C\int_{0}^{t}\left(\|A^{\frac{1}{2}}\left(\phi_{\sigma}+\alpha\phi_{p}I\right)\|^{2}+\| \phi_{u}\|^{2}+\| \phi_{\gamma}\|^{2}\right)ds\nonumber \\
	& \quad\quad\quad\quad\quad+C\int_{0}^{t}\left(\| \psi_{\sigma}\|^{2}+\| \psi_{p}\|^{2}+\| \psi_{u}\|^{2}+\| \psi_{\gamma}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{u}u-u\|_{\frac{1}{2},\Gamma _{i}}^{2}\right)ds.\label{eq:err-proof-33}
	\end{align}

	Following similar arguments, we can bound $\|\lambda_{H}^p\|_{\Gamma}$ as well.
	To achieve this, take the difference between equations (\ref{eq:cts4-sub})
	and (\ref{eq:monodd-mmmfe4}) to get
	the following error equation
	\begin{align}
	& \inp[{K^{-1}} \phi_{z}]{\zeta}+\inp[{K^{-1}} \psi_{z}]{\zeta}-\sum_{i=1}^{N}\inp[ \phi_{p}]{\dvr{\zeta}}_{\Omega_{i}}-\sum_{i=1}^{N}\inp[ \psi_{p}]{\dvr{\zeta}}_{\Omega_{i}}\nonumber \\
	& \quad-\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}p-p]{\zeta\cdot n_{i}}_{\Gamma_{i}}=\sum_{i=1}^{N}-\gnp[\lambda_{H}^{p}]{\zeta\cdot n_{i}}_{\Gamma_{i}}=\sum_{i=1}^{N}-\gnp[\phi_{\lHp}]{\zeta\cdot n_{i}}_{\Gamma_{i}}, & \forall\zeta\in Z_{h},\label{eq:err-proof-34}
	\end{align}
	where the last equality follows from (\ref{eq:motor-project-2}).
	Inf-sup stability bound (\ref{eq:mortar-p-inf-sup}) combined with
	(\ref{eq:err-proof-34}) implies
	\begin{align}
	& \|\phi_{\lHp}\|_{\Gamma}\le C\sup_{0\ne\zeta\in Z_{h}}\frac{\sum_{i=1}^{N}\gnp[\zeta\cdot n_{i}]{\phi_{\lHp}}_{\Gamma_{i}}}{\|\zeta\|_{Z_{h}}}\nonumber \\
	& =C\sup_{0\ne\zeta\in Z_{h}}(\frac{\inp[{K^{-1}} \phi_{z}]{\zeta}+\inp[{K^{-1}} \psi_{z}]{\zeta}-\sum_{i=1}^{N}\inp[ \phi_{p}]{\dvr{\zeta}}_{\Omega_{i}}-\sum_{i=1}^{N}\inp[ \psi_{p}]{\dvr{\zeta}}_{\Omega_{i}}}{\|\zeta\|_{Z_{h}}}\nonumber \\
	&\quad  +\frac{-\sum_{i=1}^{N}\gnp[{\cal I}_{H}^{p}p-p]{\zeta\cdot n_{i}}_{\Gamma_{i}}}{\|\zeta\|_{Z_{h}}})\le C\Big(\|K^{-\frac{1}{2}} \phi_{z}\|+\| \phi_{p}\| +\| \psi_{z}\|+\| \psi_{p}\| \nonumber\\
	&\hspace{160pt}  +\sum_{i=1}^{N}\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma _{i}}\Big),\label{eq:err-proof-35}
	\end{align}
	where we have used (\ref{eq:coercivity-flow-1}) and (\ref{eq:err-proof-10-2})
	with $w=p$. Squaring the above bound and then integrating with respect
	to time from $0$ to $t$ give the following
	\begin{align}
	& \|\phi_{\lHp}\|_{\Gamma}^{2}\le C\left(\|K^{-\frac{1}{2}} \phi_{z}\|^{2}+\| \phi_{p}\|^{2}\right)+C\left(\| \psi_{z}\|^{2}+\| \psi_{p}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma _{i}}^{2}\right),\label{eq:err-proof-36}\\
	& \int_{0}^{t}\|\phi_{\lHp}\|_{\Gamma}^{2}ds\le C\int_{0}^{t}\left(\|K^{-\frac{1}{2}} \phi_{z}\|^{2}+\| \phi_{p}\|^{2}\right)ds\nonumber \\
	& \quad\quad\quad\quad\quad\quad\quad\quad+C\int_{0}^{t}\left(\| \psi_{z}\|^{2}+\| \psi_{p}\|^{2}+\sum_{i=1}^{N}\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma _{i}}^{2}\right)ds.\label{eq:err-proof-37}
\end{align}

Combining (\ref{eq:err-proof-19-1}), (\ref{eq:err-proof-29})$\--$(\ref{eq:err-proof-30}),
	(\ref{eq:err-proof-32})$\--$(\ref{eq:err-proof-33}) and (\ref{eq:err-proof-36})$\--$(\ref{eq:err-proof-37}),
	and taking $\epsilon$ small enough, we obtain
	\begin{align}
	  &
          \| \phi_{\sigma}(t)\|_{\mathbb{X}_{h}}^{2}+\| \phi_{u}(t)\|^{2}+\| \phi_{\gamma}(t)\|^{2}
          +\| \phi_{z}(t)\|^{2} + \| \phi_{p}(t)\|^{2}
          +\|\phi_{\lHu}(t)\|_{\Gamma}^{2}+\|\phi_{\lHp}(t)\|_{\Gamma}^{2}
          \nonumber \\
	  &\qquad
          +\int_{0}^{t}(\| \phi_{\sigma}\|_{\mathbb{X}_{h}}^{2}
          +\| \phi_{u}\|^{2}+\| \phi_{\gamma}\|^{2}+\| \phi_{z}\|_{Z_h}^{2}
	  +\| \phi_{p}\|^{2} +\|\phi_{\lHu}\|_{\Gamma}^{2}+\|\phi_{\lHp}\|_{\Gamma}^{2})ds
\nonumber \\
          & \quad \le C \left(\int_{0}^{t}\big(\|\partial_t \psi_{\sigma}\|^{2}+\|\partial_t \psi_{p}\|^{2}+\|\partial_t \psi_{\gamma}\|^{2}+\|\partial_t \psi_{z}\|^{2}
          +\| \psi_{\sigma}\|^{2} +\| \psi_{p}\|^{2}+\| \psi_{\gamma}\|^{2}+\| \psi_{z}\|^{2}\big)ds \right.
          \nonumber \\
	  &\qquad
          +\int_{0}^{t}\big(\|{\cal I}_{H}^{p}u-u\|_{\frac{1}{2},\Gamma}^{2}
          +\|\partial_t({\cal I}_{H}^{u} u-u)\|_{\frac{1}{2},\Gamma}^{2}
	  +\|\partial_t^{2}({\cal I}_{H}^{u}u-u)\|_{\frac{1}{2},\Gamma}^{2}\nonumber \\
          &\qquad +\|{\cal I}_{H}^{p}p-p\|_{\frac{1}{2},\Gamma}^{2}
          +\|\partial_t({\cal I}_{H}^{p}p-p)\|_{\frac{1}{2},\Gamma}^{2}\big)ds
	+ \| \psi_{\sigma}(t)\|^{2}+\| \psi_{p}(t)\|^{2}+\| \psi_{\gamma}(t)\|^{2} + \| \psi_{z}(t)\|^{2}\nonumber \\
	  & \qquad + \|({\cal I}_{H}^{u}u-u)(t)\|_{\frac{1}{2},\Gamma}^{2}
          +\|\partial_t({\cal I}_{H}^{u} u- u)(t)\|_{\frac{1}{2},\Gamma}^{2}
          +\|({\cal I}_{H}^{p}p-p)(t)\|_{\frac{1}{2},\Gamma}^{2}\nonumber \\
	  & \qquad \left. +\| \phi_{\sigma}(0)\|^{2}+\| \phi_{p}(0)\|^{2}+\|\phi_{\gamma}(0)\|^{2}+\|\phi_{z}(0)\|^{2}
          +\|\partial_t({\cal I}_{H}^{u} u-u)(0)\|_{\frac{1}{2},\Gamma}^{2}\right).
          \label{eq:err-proof-38}
	\end{align}

	To bound terms of the form $\|{\cal I}_{H}^{u}v-v\|_{\frac{1}{2},\Gamma}$
	and $\|{\cal I}_{H}^{p}w-w\|_{\frac{1}{2},\Gamma}$, we use (\ref{eq:inter-proj-1})--\eqref{eq:inter-proj-1a} and (\ref{eq:trace-1})
	to obtain
\begin{align}
& \|{\cal I}_{H}^{u}v-v\|_{\frac{1}{2},\Gamma}\le CH^{\hat{m}-\frac{1}{2}}\|v\|_{\hat{m}+\frac{1}{2},\Omega}, & \frac12\le\hat{m}\le m+1, \label{eq:err-proof-39}\\
& \|{\cal I}_{H}^{p}w-w\|_{\frac{1}{2},\Gamma}\le CH^{\hat{m}-\frac{1}{2}}\|w\|_{\hat{m}+\frac{1}{2},\Omega}, & \frac12\le\hat{m}\le m+1. \label{eq:err-proof-40}
	\end{align}

	In order to bound the initial errors, $\| \phi_{\sigma}(0)\|,\,\| \phi_{p}(0)\|,\,\|\phi_{\gamma}(0)\|$,
	and $\|\phi_{z}(0)\|$, we recall that we obtain the discrete initial
	data from the elliptic projection of the continuous initial data (see
	(\ref{eq:well-pose-initial-data})). Following the arguments similar
	to the ones used to arrive at (\ref{eq:stab-bnd-13}), we get
	\begin{equation}
	\| \phi_{\sigma}(0)\|+\| \phi_{p}(0)\|+\|\phi_{\gamma}(0)\|+\|\phi_{z}(0)\|\le C\left(\| \psi_{\sigma}(0)\|+\| \psi_{p}(0)\|+\| \psi_{\gamma}(0)\|+\| \psi_{z}(0)\|+\| \psi_{u}(0)\|\right).\label{eq:err-proof-41}
	\end{equation}

Finally, combining bounds (\ref{eq:err-proof-38})$\--$(\ref{eq:err-proof-41}) with the approximation results (\ref{eq:inter-proj-3})$\--$(\ref{eq:inter-proj-9}),
	(\ref{eq:pi-weak-5})--\eqref{eq:pi-weak-div} and (\ref{eq:pi-weak-d-4})--\eqref{eq:pi-weak-d-div} completes the proof.
\end{proof}
\begin{remark}\label{rem:err-order}

	The above theorem implies that for sufficiently smooth solution variables,
	the error in using our method is of ${\cal O}\left(h^{k+1}+h^{l+1}+h^{j+1}+h^{r+1}+h^{s+1}+H^{m+\frac{1}{2}}\right)$.
	Assuming we use inf-sup stable pairs of FE spaces containing polynomials
	of degree $l=j=s$, and $k=r$, and $l\le k$, we could choose $H={\cal O}\left(h^{\frac{l+1}{m+1/2}}\right)$
	to get a total error bound of order ${\cal O}\left(h^{l+1}\right)$.
	For example, for the choice of $l=0$ and $m=1$, we could choose
	$H={\cal O}\left(h^{\frac{2}{3}}\right)$ and for $l=0$ and $m=2$,
	we could choose $H={\cal O}\left(h^{\frac{2}{5}}\right)$ to obtain
	a total convergence rate of ${\cal O}(h).$ We will demonstrate the
	results for different choices of $H(h)$ in the numerical results
	section.

\end{remark}

\section{Implementation\label{sec:Implementation-Mortar}}

In this section, we discuss the implementation of the multiscale mortar
domain decomposition technique discussed in this paper. First, we
provide a fully discrete version of the system (\ref{eq:monodd-mmmfe1})$\--$(\ref{eq:monodd-mmmfe7})
using backward Euler time discretization. We use a related formulation
where a time differentiated elasticity equation (\ref{eq:cts1-1}) is
used. The reason for such a reformulation was discussed in the previous
paper. The fully discrete formulation of the technique is then reduced
to an interface problem which can then be solved using an iterative
solver like GMRES. Finally, we discuss the possibility of computing
and storing a multiscale basis which will help in increasing the efficiency
of the method.

\subsection{Time discretization}

For time discretization, we use the backward Euler method. Let $\{t_{n}\}_{n=0}^{N_{T}}$,
$t_{n}=n\Delta t$, $\Delta t=T/N_{T}$, be a uniform partition of
$(0,T)$. With these choices, the fully discrete problem corresponding to \eqref{eq:monodd-mmmfe1}$\--$\eqref{eq:monodd-mmmfe7}
reads as follows: for $0\le n\le N_{T}-1$ and $1\le i\le N$, find\\
$(\sigma_{h,i}^{n+1},u_{h,i}^{n+1},\gamma_{h,i}^{n+1},z_{h,i}^{n+1},p_{h,i}^{n+1},\lambda_{H}^{n+1})\in\mathbb{X}_{h,i}\times V_{h,i}\times\mathbb{Q}_{h,i}\times Z_{h,i}\times W_{h,i}\times\Lambda_{H}$
such that:
\begin{align}
& \inp[A\left(\sigma_{h,i}^{n+1}+\alpha p_{h,i}^{n+1}I\right)]{\tau}_{\Omega_i}+\inp[u_{h,i}^{n+1}]{\dvr{\tau}}_{\Omega_i}+\inp[\gamma_{h,i}^{n+1}]{\tau}_{\Omega_i}\nonumber \\
& \quad\quad=\gnp[g_{u}^{n+1}]{\tau\,n_{i}}_{{\partial\Omega_i}\cap\Gamma_D^{u}}+\gnp[\lambda_{H}^{u,n+1}]{\tau\,n_{i}}_{\Gamma_{i}}, &  & \forall\tau\in\mathbb{X}_{h,i},\label{eq:mortar-dd1-mfe1-dsc}\\
& \inp[\dvr{\sigma_{h,i}^{n+1}}]{v}_{\Omega_i}=-\inp[f^{n+1}]{v}_{\Omega_i}, &  & \forall v\in V_{h,i},\label{eq:mortar-dd1-mfe2-dsc}\\
& \inp[\sigma_{h,i}^{n+1}]{\xi}_{\Omega_i}=0, &  & \forall\xi\in\mathbb{Q}_{h,i},\label{eq:mortar-dd1-mfe3-dsc}\\
& \inp[{K^{-1}} z_{h,i}^{n+1}]{\zeta}_{\Omega_i}-\inp[p_{h,i}^{n+1}]{\dvr{\zeta}}_{\Omega_i}=-\gnp[g_{p}^{n+1}]{\zeta\cdot n_i}_{{\partial\Omega_i}\cap\Gamma_D^{p}}-\gnp[\lambda_{H}^{p,n+1}]{\zeta\cdot n_i}_{\Gamma_{i}}, &  & \forall \zeta\in Z_{h,i},\label{eq:mortar-dd1-mfe4-dsc}\\
& c_{0}\inp[\frac{p_{h,i}^{n+1}-p_{h,i}^{n}}{\Delta t}]{w}_{\Omega_i}+\alpha\inp[\frac{A\left(\sigma_{h,i}^{n+1}-\sigma_{h,i}^{n}\right)}{\Delta t}]{wI}_{\Omega_i}\nonumber \\
& \quad\quad+\alpha\inp[A\alpha\frac{p_{h,i}^{n+1}-p_{h,i}^{n}}{\Delta t}I]{wI}_{\Omega_i}+\inp[\dvr{z_{h,i}^{n+1}}]{w}_{\Omega_i}=\inp[g^{n+1}]{w}_{\Omega_i}, &  & \forall w\in W_{h,i},\label{eq:mortar-dd1-mfe5-dsc}\\
& \sum_{i=1}^{N}\gnp[\sigma_{h,i}^{n+1}\,n_{i}]{\mu^{u}}_{\Gamma_{i}}=0, &  & \forall\mu^{u}\in\Lambda_{H}^{u},\label{eq:mortar-dd1-mfe6-dsc}\\
& \sum_{i=1}^{N}\gnp[z_{h,i}^{n+1}\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}=0, &  & \forall\mu^{p}\in\Lambda_{H}^{p}.\label{eq:mortar-dd1-mfe7-dsc}
\end{align}
Note that similarly to the matching grid case in \cite{dd-biot},
the non-matching multiscale grid method also requires initial data
$p_{h,i}^{0}$ and $\sigma_{h,i}^{0}$. Such data can be obtained
by taking $p_{h,i}^{0}$ to be the $L^{2}$-projection of $p_{0}$
onto $W_{h,i}$ and solving a mixed elasticity non-matching grid non-overlapping
domain decomposition problem obtained from (\ref{eq:mortar-dd1-mfe1-dsc})$\--$(\ref{eq:mortar-dd1-mfe3-dsc})
and (\ref{eq:mortar-dd1-mfe6-dsc}) with $n=-1$ (see \cite{eldar_elastdd}).
Also, following the case of the matching grid domain decomposition
method, we will utilize a related formulation in which the elasticity
equation, (\ref{eq:cts1-1}), is differentiated in time. The reason
for this was discussed in the analysis of the resulting interface
problem in \cite{dd-biot}. We introduce new variables $\dot{u}=\partial_t u$,
$\dot{\gamma}=\partial_t\gamma$, and $\lambda_{H}^{\dot{u}}=\partial_t\lambda_{H}^u$ representing
the time derivatives of the displacement, rotation, and displacement-Lagrange
multiplier function, respectively. Using time-differentiated equation
(\ref{eq:cts1-1}) and backward Euler (see \cite{dd-biot} for details), we get

\begin{align}
& \inp[A(\sigma_{h,i}^{n+1}+\alpha p_{h,i}^{n+1}I)]{\tau}_{\Omega_i}+\Delta t\inp[\dot{u}_{h,i}^{n+1}]{\dvr{\tau}}_{\Omega_i}+\Delta t\inp[\dot{\gamma}_{h,i}^{n+1}]{\tau}_{\Omega_i}\nonumber \\
& \quad\quad=\Delta t\gnp[\partial_t g_{u}^{n+1}]{\tau\,n_{i}}_{{\partial\Omega}_{i}\cap\Gamma_D^{u}}+\Delta t\gnp[\lambda_{H}^{\dot{u},n+1}]{\tau\,n_{i}}_{\Gamma_{i}}+\inp[A(\sigma_{h,i}^{n}+\alpha p_{h,i}^{n}I)]{\tau}_{\Omega_i}, & \forall\tau\in\mathbb{X}_{h,i}.\label{eq:mortar-time-diff-elast}
\end{align}
The original variables can be recovered using\begin{equation}\label{init-recover} u_h^n = u_h^0 + \Delta t \sum_{k=1}^n \dot{u}_h^k, \quad \gamma_h^n = \gamma_h^0 + \Delta t \sum_{k=1}^n \dot{\gamma}_h^k, \quad \lambda_H^{u,n} = \lambda_H^{u,0} + \sum_{k=1}^n \dot \lambda_H^{u,k}. \end{equation}

\subsection{Reduction to an interface problem}

Similar to the non-overlapping matching grid domain decomposition
algorithm discussed in \cite{dd-biot}, we solve the system resulting
from (\ref{eq:mortar-time-diff-elast}), (\ref{eq:mortar-dd1-mfe2-dsc})$\--$(\ref{eq:mortar-dd1-mfe7-dsc})
at each time step by reducing it to an interface problem for the Lagrange
multiplier function $\lambda_{H}$. Following similar arguments as
in \cite{dd-biot}, we introduce two sets of complementary subdomain problems.

The first set of problems reads as follows: for $1\le i\le N$, find\\
$(\bar{\sigma}_{h,i}^{n+1},\bar{\dot{u}}_{h,i}^{n+1},\bar{\dot{\g}}_{h,i}^{n+1},\bar{z}_{h,i}^{n+1},\bar{p}_{h,i}^{n+1})\in\mathbb{X}_{h,i}\times V_{h,i}\times\mathbb{Q}_{h,i}\times Z_{h,i}\times W_{h,i}$
such that
\begin{align}
& \inp[A(\bar{\sigma}_{h,i}^{n+1}+\alpha\bar{p}_{h,i}^{n+1} I)]{\tau}_{\Omega_i}+\Delta t\inp[\bar{\dot{u}}_{h,i}^{n+1}]{\dvr{\tau}}_{\Omega_i}+\Delta t\inp[\bar{\dot{\g}}_{h,i}^{n+1}]{\tau}_{\Omega_i}\nonumber \\
& \quad\quad=\Delta t\gnp[g_{u}^{n+1}]{\tau\,n_i}_{{\partial\Omega_i}\cap\Gamma_D^{u}}+\inp[A(\sigma_{h,i}^{n}+\alpha p_{h,i}^{n}I)]{\tau}_{\Omega_i}, &  & \forall\tau\in\mathbb{X}_{h,i},\label{eq:mortar-bar-1-1}\\
& \inp[\dvr{\bar{\sigma}_{h,i}^{n+1}}]{v}_{\Omega_i}=-\inp[f^{n+1}]{v}_{\Omega_i}, &  & \forall v\in V_{h,i},\label{eq:mortar-bar-2-1}\\
& \inp[\bar{\sigma}_{h,i}^{n+1}]{\xi}_{\Omega_i}=0, &  & \forall\xi\in\mathbb{Q}_{h,i},\label{eq:mortar-bar-3-1}\\
& \inp[{K^{-1}}\bar{z}_{h,i}^{n+1}]{\zeta}_{\Omega_i}-\inp[\bar{p}_{h,i}^{n+1}]{\dvr{\zeta}}_{\Omega_i}=-\gnp[g_{p}^{n+1}]{\zeta\cdot n_i}_{{\partial\Omega_i}\cap\Gamma_D^{p}}, &  & \forall \zeta\in Z_{h,i},\label{eq:mortar-bar-4-1}\\
& c_{0}\inp[\bar{p}_{h,i}^{n+1}]{w}_{\Omega_i}+\alpha\inp[A(\bar{\sigma}_{h,i}^{n+1}+\alpha\bar{p}_{h,i}^{n+1} I)]{wI}_{\Omega_i}+\Delta t\inp[\dvr{\bar{z}_{h,i}^{n+1}}]{w}_{\Omega_i}\nonumber \\
& \quad\quad=\Delta t\inp[g^{n+1}]{w}_{\Omega_i}+c_{0}\inp[p_{h,i}^{n}]{w}_{\Omega_i}+\alpha\inp[A(\sigma_{h,i}^{n}+\alpha p_{h,i}^{n}I)]{wI}_{\Omega_i}, &  & \forall w\in W_{h,i}.\label{eq:mortar-bar-5-1}
\end{align}
Note that these subdomain problems have zero Dirichlet data at the
subdomain interface, the true source terms $f$ and $g$, outside
boundary conditions $g_{u}$ and $g_{p}$, and initial data $\sigma_{h,i}^{n}$
and $p_{h,i}^{n}$.

The second set of equations reads as follows: given $\lambda_{H}\in\Lambda_{H}$,
for $1\le i\le N$, find $(\sigma_{h,i}^{*,n+1}(\lambda_{H}),$ $\dot{u}_{h,i}^{*,n+1}(\lambda_{H}),\dot{\g}_{h,i}^{*,n+1}(\lambda_{H}),z_{h,i}^{*,n+1}(\lambda_{H}),p_{h,i}^{*,n+1}(\lambda_{H}))$$\in\mathbb{X}_{h,i}\times V_{h,i}\times\mathbb{Q}_{h,i}\times Z_{h,i}\times W_{h,i}$
such that:
\begin{align}
& \inp[A\big(\sigma_{h,i}^{*,n+1}(\lambda_{H})+\alphap_{h,i}^{*,n+1}(\lambda_{H})I\big)]{\tau}_{\Omega_i}+\Delta t\inp[\dot{u}_{h,i}^{*,n+1}(\lambda_{H})]{\dvr{\tau}}_{\Omega_i}\nonumber \\
& \qquad\quad+\Delta t\inp[\dot{\g}_{h,i}^{*,n+1}(\lambda_{H})]{\tau}_{\Omega_i}=\Delta t\left\langle \lambda_{H}^{\dot{u}},\tau\,n_{i}\right\rangle _{\Gamma_{i}}, & \forall\tau\in\mathbb{X}_{h,i},\label{eq:mortar-star-1-1}\\
& \left(\dvr\sigma_{h,i}^{*,n+1}(\lambda_{H}),v\right)_{\Omega_i}=0, & \forall v\in V_{h,i},\label{eq:mortar-star-2-1}\\
& \left(\sigma_{h,i}^{*,n+1}(\lambda_{H}),\xi\right)_{\Omega_i}=0, & \forall\xi\in\mathbb{Q}_{h,i},\label{eq:mortar-star-3-1}\\
& \left({K^{-1}}z_{h,i}^{*,n+1}(\lambda_{H}),\zeta\right)_{\Omega_i}-\left(p_{h,i}^{*,n+1}(\lambda_{H}),\dvr{\zeta}\right)_{\Omega_i}=-\gnp[\lambda_{H}^{p,n+1}]{\zeta\cdot n_i}_{\Gamma_{i}}, & \forall \zeta\in Z_{h,i},\label{eq:mortar-star-4-1}\\
& c_{0}\left(p_{h,i}^{*,n+1}(\lambda_{H}),w\right)+\alpha\left(A\big(\sigma_{h,i}^{*,n+1}(\lambda_{H})+\alphap_{h,i}^{*,n+1}(\lambda_{H})I\big),wI\right)_{\Omega_i}\nonumber \\
& \quad\qquad+\Delta t\left(\dvrz_{h,i}^{*,n+1}(\lambda_H)\right)=0, & \forall w\in W_{h,i}.\label{eq:mortar-star-5-1}
\end{align}
Note that this set of problems has $\lambda_{H}$ as the Dirichlet
boundary data on the interface $\Gamma$, compared to $\lambda_{h}$ in
the matching grid case. This system also has zero source terms, zero
boundary data on part of the outside boundary ${\partial\Omega}$, and zero data
from the previous time step.

Define the bilinear form $a_{H,i}^{n+1}:\lambda_H\times\lambda_H\to\mathbb{R}$, $1\le i\le N$,
$a_{H}^{n+1}:\lambda_H\times\lambda_H\to\mathbb{R}$, and the linear functional $g_{H}^{n+1}:\lambda_H\to\mathbb{R}$
for all $0\le n\le N_{T}-1$ by
\begin{align}
& a_{H,i}^{n+1}(\lambda_{H},\mu)=\gnp[\sigma_{h,i}^{*,n+1}(\lambda_{H})\,n_{i}]{\mu^{u}}_{\Gamma_{i}}-\gnp[z_{h,i}^{*,n+1}(\lambda_{H})\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}},\quad a_{H}^{n+1}(\lambda_{H},\mu)=\sum_{i=1}^{N}a_{H,i}^{n+1}(\lambda_{H},\mu),\label{eq:mortar-interface-optr-1}\\
& g_{H}^{n+1}(\mu)=\sum_{i=1}^{N}\left(-\gnp[\bar{\sigma}_{h,i}^{n+1}\,n_{i}]{\mu^{u}}_{\Gamma_{i}}+\gnp[\bar{z}_{h,i}^{n+1}\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}\right).\label{eq:mortar-intr-rhs-1}
\end{align}
It follows from (\ref{eq:mortar-dd1-mfe6-dsc})$\--$(\ref{eq:mortar-dd1-mfe7-dsc})
that for each $0\le n\le N_{T}-1$, the solution to the global problem
(\ref{eq:mortar-time-diff-elast}), (\ref{eq:mortar-dd1-mfe2-dsc})$\--$(\ref{eq:mortar-dd1-mfe7-dsc})
is equivalent to solving the interface problem for $\lambda_{H}^{n+1}\in\Lambda_{H}$:
\begin{equation}
a_{H}^{n+1}(\lambda_{H}^{n+1},\mu)=g_{H}^{n+1}(\mu),\quad\forall\mu\in\Lambda_{H},\label{eq:mortar-int-problm-1}
\end{equation}
and setting
\begin{align*}
& \sigma_{h,i}^{n+1}=\sigma_{h,i}^{*,n+1}(\lambda_{H}^{n+1})+\bar{\sigma}_{h,i}^{n+1},\quad\dot{u}_{h,i}^{n+1}=\dot{u}_{h,i}^{*,n+1}(\lambda_{H}^{n+1})+\bar{\dot{u}}_{h,i}^{n+1},\quad\dot{\gamma}_{h,i}^{n+1}=\dot{\g}_{h,i}^{*,n+1}(\lambda_{H}^{n+1})+\bar{\dot{\g}}_{h,i}^{n+1},\\
& z_{h,i}^{n+1}=z_{h,i}^{*,n+1}(\lambda_{H}^{n+1})+\bar{z}_{h,i}^{n+1},\quad p_{h,i}^{n+1}=p_{h,i}^{*,n+1}(\lambda_{H}^{n+1})+\bar{p}_{h,i}^{n+1}.
\end{align*}

\subsection{Solving the interface problem}

In order to solve the interface problem (\ref{eq:mortar-int-problm-1}),
we introduce linear operators $\mathcal{A}_{H,i}^{n+1}:\Lambda_{H,i}\to\Lambda_{H,i}$,
for $1\le i\le N$ and $\mathcal{A}_{H}^{n+1}:\Lambda_{H}\to\Lambda_{H}$ such that
for any $\lambda_H\in\Lambda_H$,
\begin{align}
& \langle\mathcal{A}_{H,i}^{n+1}\lambda_{H,i},\,\mu\rangle_{\Gamma _{i}}=\gnp[\sigma_{h,i}^{*,n+1}(\lambda_{H})\,n_{i}]{\mu^{u}}_{\Gamma_{i}}-\gnp[z_{h,i}^{*,n+1}(\lambda_{H})\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}, & \forall\mu\in\Lambda_{H,i},\label{eq:int-op-matrix-1}\\
& \mathcal{A}_{H}^{n+1}\lambda_H=\sum_{i=1}^{N}\mathcal{A}_{H,i}^{n+1}\lambda_{H,i},\label{eq:int-op-matrix-2}
\end{align}
where $\lambda_{H,i}$ and $\Lambda_{H,i}$ denote the restrictions of $\lambda_H$
and $\Lambda_{H}$ to $\Gamma _{i}$, respectively. We also define the vector
$G_{H}^{n+1}\in\Lambda_H$ as

\begin{equation}
\gnp[G_{H}^{n+1}]{\mu}_{\Gamma}=\sum_{i=1}^{N}\left(-\gnp[\bar{\sigma}_{h,i}^{n+1}\,n_{i}]{\mu^{u}}_{\Gamma_{i}}+\gnp[\bar{z}_{h,i}^{n+1}\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}\right),\qquad\forall\mu\in\Lambda_{H,i}.\label{eq:int-rhs-vector}
\end{equation}

Interface problem (\ref{eq:mortar-int-problm-1}) can now be reformulated
as finding $\lambda_H\in\Lambda_H$ such that
\begin{equation}
\mathcal{A}_{H}^{n+1}\lambda_H=G_{H}^{n+1}.\label{eq:mortar-Matrix-int-problm}
\end{equation}

Consider the $L^{2}$ orthogonal projections, $\mathcal{Q}_{h,i}^{u,T}:\mathbb{X}_{h,i}n_{i}\to\Lambda_{H}^{u}$
and $\mathcal{Q}_{h,i}^{p,T}:Z_{h,i}\cdot n_{i}:\to\Lambda_{H}^{p}$
such that for any $\tau\in\mathbb{X}_{h,i}$ and $\zeta\in Z_{h,i}$,

\begin{align*}
& \gnp[\mathcal{Q}_{h,i}^{u,T}\left(\tau\,n_{i}\right)-\tau\,n_{i}]{\mu^{u}}_{\Gamma_{i}}=0, & \forall\mu^{u}\in\Lambda_{H}^{u},\\
& \gnp[\mathcal{Q}_{h,i}^{p,T}\left(\zeta\cdot n_{i}\right)-\zeta\cdot n_{i}]{\mu^{p}}_{\Gamma_{i}}=0, & \forall\mu^{p}\in\Lambda_{H}^{p}.
\end{align*}
where $n_{i}$ is the unit outward normal to ${\partial\Omega_i}$. Define $\mathcal{Q}_{h,i}^{T}:\mathbb{X}_{h,i}n_{i}\times Z_{h,i}\cdot n_{i}\to\Lambda_{H}^{u}\times\Lambda_{H}^{p}$
such that for any $\tau\in\mathbb{X}_{h,i}$ and $\zeta\in Z_{h,i}$
\begin{equation}
\mathcal{Q}_{h,i}^{T}\begin{pmatrix}\tau\,n_{i}\\
\zeta\cdot n_{i}
\end{pmatrix}=\begin{pmatrix}\mathcal{Q}_{h,i}^{u,T}\tau\,n_{i}\\
\mathcal{Q}_{h,i}^{p,T}\zeta\cdot n_{i}
\end{pmatrix}.\label{eq:MrtToTraceProj-1}
\end{equation}

Using the above notations, we note that
\begin{align}
& \mathcal{A}_{H,i}^{n+1}\lambda_{H,i}=\mathcal{Q}_{h,i}^{T}\begin{pmatrix}\sigma_{h,i}^{*,n+1}(\lambda_{H})\,n_{i}\\
-z_{h,i}^{*,n+1}(\lambda_{H})\cdot n_{i}
\end{pmatrix},\qquadG_{H}^{n+1}=\begin{pmatrix}\sum_{i=1}^{N}-\mathcal{Q}_{h,i}^{u,T}\bar{\sigma}_{h,i}^{n+1}\,n_{i}\\
\sum_{i=1}^{N}\mathcal{Q}_{h,i}^{p,T}\bar{z}_{h,i}^{n+1}\cdot n_{i}
\end{pmatrix},\label{eq:Matrix-projection-def}
\end{align}
for $i=1,2,\ldots N.$

To solve the interface problem (\ref{eq:mortar-Matrix-int-problm}),
we use an iterative method like GMRES with an initial guess $\lambda_{H,0}\in\Lambda_H$
at each time step, $t_{n}=n\Delta t$. A detailed description is given
in Algorithm \ref{alg:Solving-interface-problem} (also see \cite{ganis2009implementation}).
\\
\\
\begin{algorithm}[h]
	\caption{\label{alg:Solving-interface-problem}Solving interface problem using
		GMRES algorithm.}

	\begin{enumerate}
		\item Solve the first set of complementary equations, (\ref{eq:mortar-bar-1-1})$\--$(\ref{eq:mortar-bar-5-1}),
		and compute $G_{H}^{n+1}$ using (\ref{eq:Matrix-projection-def}).
		\item Pick an initial guess $\lambda_{H,0}\in\Lambda_H$.
		\item Project the mortar function onto the subdomain boundaries, $\lambda_{H,0,i}\longrightarrow\mathcal{Q}_{h,i}\left(\lambda_{H,0,i}\right)$.
		\item Solve the second set of complementary equations, (\ref{eq:mortar-star-1-1})$\--$(\ref{eq:mortar-star-5-1}),
		using the projected function $\mathcal{Q}_{h,i}\left(\lambda_{H,0,i}\right)$
		as Dirichlet boundary data to obtain $\sigma_{h,i}^{*,n+1}(\lambda_{H,0})$ and $z_{h,i}^{*,n+1}(\lambda_{H,0})$.
		\item Project the solution variables to the mortar space, $\sigma_{h,i}^{*,n+1}(\lambda_{H,0})\,n_{i}\longrightarrow\mathcal{Q}_{h,i}^{u,T}\sigma_{h,i}^{*,n+1}(\lambda_{H,0})\,n_{i}$
		and $z_{h,i}^{*,n+1}(\lambda_{H,0})\cdot n_{i}\longrightarrow\mathcal{Q}_{h,i}^{p,T}z_{h,i}^{*,n+1}(\lambda_{H,0})\cdot n_{i}$.
		\item Compute the action $\mathcal{A}_{H}^{n+1}\lambda_{H,0}$ using (\ref{eq:Matrix-projection-def}).
		\item Update $\lambda_H$ using $\mathcal{A}_{H}^{n+1}\lambda_{H,0}$ in the GMRES algorithm.
		\item Repeat steps $3\--7$, with updated values of $\lambda_H$, until the residual
		for the GMRES algorithm goes below a predetermined tolerance.
	\end{enumerate}
\end{algorithm}

This method has the performance advantage over the similar method
for matching grids discussed in \cite{dd-biot} that a coarse mortar mesh could be used to obtain a smaller interface problem due to the reduction in the mortar degrees of freedom. We
implement this algorithm and study various test cases in the numerical
results section.

\subsection{Implementation with multiscale stress-flux basis (MSB)\label{subsec:Implementation-of-multiscale}}

A coarser mortar mesh can lead to a smaller interface problem, but
even in that case, the number of subdomain solves of the type (\ref{eq:mortar-star-1-1})$\--$(\ref{eq:mortar-star-5-1})
is directly proportional to both the number of mortar space degrees
of freedom and the number of time steps used. We propose the construction
and use of a multiscale stress-flux basis (MSB) which makes the number
of subdomain solves independent of the number of iterations required
for the interface problem and the number of time steps used.

Let $\left\{ \mathcal{\beta}_{H,i}^{k}\right\} _{k=0}^{N_{H}}$ be a basis for $\Lambda_{H,i}$,
where $N_{H}$ denotes the number of degrees of freedom associated
with the finite element space $\Lambda_{H,i}$. We calculate and store the
action of the interface operator of the form
\begin{align}
& \mathcal{A}_{H,i}\mathcal{\beta}_{H,i}^{k}=\mathcal{Q}_{h,i}^{T}\begin{pmatrix}\sigma_{h,i}^{*}(\mathcal{\beta}_{H,i}^{k})\,n_{i}\\
-z_{h,i}^{*}(\mathcal{\beta}_{H,i}^{k})\cdot n_{i}
\end{pmatrix},\label{eq:int-op-matrix-1-generic}
\end{align}
for $k=1,2,\ldots N_{H}$, where we obtain $\sigma_{h,i}^{*}(\mathcal{\beta}_{H,i}^{k})$
and $z_{h,i}^{*}(\mathcal{\beta}_{H,i}^{k})$ by solving (\ref{eq:mortar-star-1-1})$\--$(\ref{eq:mortar-star-5-1})
with $\mathcal{\beta}_{H,i}^{k}$ as the Dirichlet boundary data. A detailed description
of the construction of the multiscale basis elements $\left\{ \mathcal{\phi}_{H,i}^{k}\right\} _{k=0}^{N_{H}}$,
where $\mathcal{\phi}_{H,i}^{k}=\mathcal{A}_{H,i}\mathcal{\beta}_{H,i}^{k}$ is given in Algorithm \ref{alg:Construction-of-multiscale},
(also see \cite{ganis2009implementation,eldar_elastdd} for similar
constructions).

\begin{algorithm}[h]
	\caption{\label{alg:Construction-of-multiscale}Construction of multiscale
		stress-flux basis}

	\textbf{for} $k=1,\ldots,N_{H}$:
	\begin{enumerate}
		\item Project $\mathcal{\beta}_{H,i}^{k}$ onto the subdomain boundary, $\mathcal{\beta}_{H,i}^{k}\longrightarrow\mathcal{Q}_{h,i}\left(\mathcal{\beta}_{H,i}^{k}\right)$.
		\item Solve the system (\ref{eq:mortar-star-1-1})$\--$(\ref{eq:mortar-star-5-1})
		using the projected function $\mathcal{Q}_{h,i}\left(\mathcal{\beta}_{H,i}^{k}\right)$
		as the Dirichlet boundary data, to obtain $\sigma_{h,i}^{*}(\mathcal{\beta}_{H,i}^{k})$ and
		$z_{h,i}^{*}(\mathcal{\beta}_{H,i}^{k})$.
		\item Project the solution variables to the mortar space to obtain $\mathcal{\phi}_{H,i}^{k}=\begin{pmatrix}\mathcal{Q}_{h,i}^{u,T}\sigma_{h,i}^{*}(\mathcal{\beta}_{H,i}^{k})\,n_{i}\\
		-\mathcal{Q}_{h,i}^{p,T}z_{h,i}^{*}(\mathcal{\beta}_{H,i}^{k})\cdot n_{i}
		\end{pmatrix}$.
	\end{enumerate}
	\textbf{end for}
\end{algorithm}

For any $\mu\in\Lambda_{H,i}$, consider the mortar basis decomposition, $\mu=\sum_{k=0}^{N_{H}}\mu_{i}\mathcal{\beta}_{H,i}^{k}$.
We can calculate the action of the interface operator on $\mu$ as
follows:
\begin{equation}
\mathcal{A}_{H,i}\mu=\sum_{k=0}^{N_{H}}\mu_{i}\mathcal{\phi}_{H,i}^{k}.\label{eq:multi-basis-use}
\end{equation}

We can use the multiscale basis to replace steps $3\--6$ in Algorithm
\ref{alg:Solving-interface-problem}, by taking linear combinations
of the form (\ref{eq:multi-basis-use}). Note that the multiscale
stress-flux basis is computed and saved once and can be reused over
all time steps, which gains a significant performance advantage in
the case of time-dependent parabolic problems like the one we are
studying.

We further discuss and compare the efficiency of using the multiscale
stress-flux basis with other methods in Example 2 in the numerical
section of this paper.

\section{Numerical Results\label{sec:Numerical-results-mortar}}
\setcounter{figure}{3}
In this section, we report the results of various numerical tests
designed to verify and compare the well-posedness, stability, and
convergence of the multiscale mortar non-overlapping domain decomposition
method for the Biot system of poroelasticity that we have developed
in this paper. We compare the computational efficiency in different
cases, including the matching and non-matching grids on the subdomains,
and also discuss the advantages of using a multiscale basis. The numerical
schemes are implemented using deal.II finite element package \cite{dealII90,BangerthHartmannKanschat2007}.

In all the examples, we have used the FE triplet $\mathbb{X}_{h}\times V_{h}\times\mathbb{Q}_{h}=\mathcal{BDM}_{1}^{2}\times Q_{0}^{2}\times Q_{0}$
(\cite{ArnAwaQiu}) for elasticity and the FE pair $Z_{h}\times W_{h}=\mathcal{BDM}_{1}\times Q_{0}$
(\cite{brezzi1991mixed}) for Darcy on quadrilateral meshes. Here
$Q_{k}$ denotes polynomials of degree $k$ in each variable. For
the mortar spaces, $\lambda_{H}^u$ is taken to be $DQ_{m}^{2}$, and $\lambda_{H}^p$
is taken to be $DQ_{m}$ with $m=1\text{ or }2$, where $DQ_{k}$ represents the discontinuous
finite element spaces containing polynomials of degree $k$, which
lives on the subdomain interface. The degrees of FEM spaces used in this
section is given in Table \ref{tab:Degree-poly-numerics}. For solving
the interface problem, we use non-restarted unpreconditioned GMRES
with a tolerance on the relative residual $\frac{r_{k}}{r_{0}}$ as
the stopping criteria. For all the examples, this tolerance is taken
to be $10^{-6}$.

\begin{table}[h]

	\caption{\label{tab:Degree-poly-numerics}Degree of polynomials associated
		with FEM spaces used for numerical experiments.}

	\begin{centering}
		\begin{tabular}{|c|c|c|c|c|c|}
			\hline
			$\Lambda_{H}:m$ & $\mathbb{X}_{h}:k$ & $V_{h}:l$ & $\mathbb{Q}_{h}:j$ & $Z_{h}:r$ & $W_{h}:s$\tabularnewline
			\hline
			$1$ (linear) & 1 & 0 & 0 & 1 & 0\tabularnewline
			$2$ (quadratic) & 1 & 0 & 0 & 1 & 0\tabularnewline
			\hline
		\end{tabular}
		\par\end{centering}
\end{table}

In Example 1, we test and compare the convergence, stability, and
efficiency of the multiscale mortar DD method using linear ($m=1$)
and quadratic ($m=2$) mortar spaces. We do this by solving a system
of equations with a known solution on successively refined meshes.

In Example 2, we apply the multiscale DD method to solve a more practical
problem, using a highly heterogeneous medium. We compare the efficiency
of the multiscale versus fine scale methods and study the computational
advantage of constructing a multiscale stress-flux basis (MSB) discussed
in Section \ref{subsec:Implementation-of-multiscale} of this paper.

In Example 1, we solve the system of PDEs on a checkerboard
global mesh, which consists of non-matching grids on all subdomain
interfaces. In particular, the coarsest multiscale mesh in all
examples follows a subdomain mesh size ratio $\frac{1}{4}:\frac{1}{6}:\frac{1}{6}:\frac{1}{4}$
as shown in Figure \ref{fig:multi-grid-mesh}. The corresponding coarsest
mortar-interface mesh consists of two elements with mesh size $\frac{1}{2}$.

\begin{figure}[h]
	\begin{centering}
		\includegraphics[width=0.5\columnwidth]{figures/mesh_checker.jpeg}
		\par\end{centering}
	\caption{\label{fig:multi-grid-mesh}Example 1, coarsest non matching subdomain
		grid on $(0,1)^{2}.$}

\end{figure}

\subsection{Example 1: testing convergence rates}

In this example, we test the well-posedness, convergence, and stability
of the multiscale mortar DD method using linear $(m=1)$ and quadratic
$(m=2)$ mortar spaces. The global computational domain $\Omega$ is taken
to be the unit square $(0,1)^{2}$. We consider the following analytical
solution
\[
p=\exp(t)(\sin(\pi x)\cos(\pi y)+10),\quad u=\exp(t)\begin{pmatrix}x^{3}y^{4}+x^{2}+\sin((1-x)(1-y))\cos(1-y)\\
(1-x)^{4}(1-y)^{3}+(1-y)^{2}+\cos(xy)\sin(x)
\end{pmatrix}.
\]
The physical and numerical parameters are given in Table \ref{tab:Physical-Parameters_ex1-mortar}.
Using this information, we derive the right hand side and initial
conditions essential to solve the system \eqref{biot-1}$\--$\eqref{biot-bc-2}.
We partition $\Omega$ into four square subdomains with non-matching
grids as shown in Figure \ref{fig:multi-grid-mesh}.

We consider two different cases, with linear and quadratic mortar
spaces, where the former contains polynomials of degree $1$ and the
latter contains polynomials of degree $2$. To test the convergence
and verify the theoretical a priori error estimates, we successively
refine the subdomain and mortar meshes. In the linear mortar case,
we maintain a subdomain to mortar mesh ratio such that $H=Ch$, and
in the quadratic mortar case, we maintain the ratio such that $H=C\sqrt{h}$.

The convergence tables for the cases with linear and quadratic mortar
spaces with $\Delta t=10^{-4}$ and $c_{0}=1.0$ are given in Tables
\ref{tab:Example-1-lin-conv} and \ref{tab:Example-1-quad-conv},
respectively. Tables \ref{tab:Example-1-lin-conv_c0_small}$\--$\ref{tab:Example-1-quad-conv-small-c0}
present the convergence table in the case of linear mortar and quadratic
mortar spaces, respectively with $\Delta t=10^{-4}$, and $c_{0}=10^{-3}$.
We present the number of interface iterations, relative errors, and
their convergence rates in these tables. Solution plots in the case
of linear mortar with an intermediate level of refinement, $h=1/32$,
$\Delta t=10^{-3}$ and $c_{0}=1.0$ is given in Figure \ref{fig:Example-1,-computed_exact_mortar}
in order to compare with the plots obtained in \cite{dd-biot}
using monolithic domain decomposition technique using matching grids
on subdomain interfaces. Note that plots in the case of quadratic
mortar space look similar.

\begin{table}[h]
	\begin{centering}
		\caption{Example 1, physical and numerical parameters.\label{tab:Physical-Parameters_ex1-mortar}}
		\begin{tabular}{c|c} \hline Parameter & Value
			\tabularnewline \hline Permeability tensor $(K)$ & $I$
			\tabularnewline Lame coefficient $(\mu)$ & $100.0$
			\tabularnewline Lame coefficient $(\lambda)$ & $100.0$
			\tabularnewline Mass storativity $(c_{0})$ & $1.0, 10^{-3}$
			\tabularnewline Biot-Willis constant $(\alpha)$ & $1.0$
			\tabularnewline Time step $(\Delta t)$ & $10^{-3}, 10^{-4}$
			\tabularnewline Number of time steps & $100$
			\tabularnewline
			\hline
		\end{tabular}
		\par\end{centering}
\end{table}

\renewcommand{\tabcolsep}{5pt}

\begin{table}[h]
	\captionsetup{justification=centering}
	\caption{\label{tab:Example-1-lin-conv}Example 1, convergence table using
		linear mortar ($m=1$) with $H=Ch$, $\Delta t=10^{-4}$ and $c_{0}=1.0$.}

	\begin{centering}
		\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
			\hline
			$h$ & \multicolumn{2}{c|}{\# GMRES} & \multicolumn{2}{c|}{$\|\sigma-\sigma_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\dvr\,(\sigma-\sigma_{h})\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\gamma-\gamma_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|u-u_{h}\|_{L^{\infty}(L^{2})}$}\tabularnewline
			\hline
			$1/4$ & 16  & rate & 1.23e-01  & rate  & 6.09e-01  & rate  & 1.39e+00  & rate  & 5.78e-01  & rate \tabularnewline
			\cline{3-3} \cline{5-5} \cline{7-7} \cline{9-9} \cline{11-11}
			$1/8$ & 28  & -0.81  & 3.24e-02  & 1.92  & 3.11e-01  & 0.97  & 7.07e-01  & 0.97  & 2.92e-01  & 0.99 \tabularnewline
			$1/16$ & 46  & -0.72  & 8.20e-03  & 1.98  & 1.56e-01  & 0.99  & 3.55e-01  & 0.99  & 1.46e-01  & 1.00 \tabularnewline
			$1/32$ & 73  & -0.67  & 2.08e-03  & 1.98  & 7.82e-02  & 1.00  & 1.78e-01  & 1.00  & 7.31e-02  & 1.00 \tabularnewline
			$1/64$ & 122  & -0.74  & 5.39e-04  & 1.94  & 3.91e-02  & 1.00  & 8.89e-02  & 1.00  & 3.65e-02  & 1.00 \tabularnewline
			\hline
		\end{tabular}
		\par\end{centering}
	\centering{}%
	\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
		\hline
		$h$ & \multicolumn{2}{c|}{$\|z-z_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\dvr\,(z-z_{h})\|_{L^{2}(L^{2})}$} & \multicolumn{2}{c|}{$\|p-p_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|u-\lambda^{u}{}_{H}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|p-\lambda^{p}{}_{H}\|_{L^{\infty}(L^{2})}$}\tabularnewline
		\hline
		$1/4$ & 1.04e+00  & rate  & 4.15e-01  & rate  & 5.91e-02  & rate  & 7.50e-01  & rate  & 2.06e-01  & rate \tabularnewline
		\cline{3-3} \cline{5-5} \cline{7-7} \cline{9-9} \cline{11-11}
		$1/8$ & 3.72e-01  & 1.48  & 1.89e-01  & 1.14  & 2.96e-02  & 1.00  & 1.90e-01  & 1.98  & 5.30e-02  & 1.96\tabularnewline
		$1/16$ & 1.19e-01  & 1.64  & 8.50e-02  & 1.15  & 1.48e-02  & 1.00  & 4.76e-02  & 1.99  & 1.33e-02  & 2.00\tabularnewline
		$1/32$ & 3.56e-02  & 1.74  & 3.97e-02  & 1.10  & 7.39e-03  & 1.00  & 1.19e-02  & 2.00  & 3.33e-03  & 2.00\tabularnewline
		$1/64$ & 1.08e-02  & 1.72  & 1.92e-02  & 1.05  & 3.70e-03  & 1.00  & 3.04e-03  & 1.97  & 8.37e-04  & 1.99\tabularnewline
		\hline
	\end{tabular}
\end{table}

\begin{table}[h]
	\captionsetup{justification=centering}
	\caption{\label{tab:Example-1-quad-conv}Example 1, convergence table using
		quadratic mortar ($m=2$) with $H=C\sqrt{h}$, $\Delta t=10^{-4}$
		and $c_{0}=1.0$.}

	\begin{centering}
		\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
			\hline
			$h$ & \multicolumn{2}{c|}{\# GMRES} & \multicolumn{2}{c|}{$\|\sigma-\sigma_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\dvr\,(\sigma-\sigma_{h})\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\gamma-\gamma_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|u-u_{h}\|_{L^{\infty}(L^{2})}$}\tabularnewline
			\hline
			$1/4$ & 22  & rate & 1.26e-01  & rate & 6.09e-01  & rate & 1.39e+00  & rate & 5.79e-01  & rate\tabularnewline
			\cline{3-3} \cline{5-5} \cline{7-7} \cline{11-11}
			$1/16$ & 40  & -0.43  & 8.25e-03  & 1.97 & 1.56e-01  & 0.98 & 3.55e-01  & 0.99  & 1.46e-01  & 0.99 \tabularnewline
			$1/64$ & 65  & -0.35  & 5.62e-04  & 1.93 & 3.91e-02  & 1.00  & 8.89e-02  & 1.00  & 3.65e-02  & 1.00 \tabularnewline
			\hline
		\end{tabular}
		\par\end{centering}
	\centering{}%
	\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
		\hline
		$h$ & \multicolumn{2}{c|}{$\|z-z_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\dvr\,(z-z_{h})\|_{L^{2}(L^{2})}$} & \multicolumn{2}{c|}{$\|p-p_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|u-\lambda^{u}{}_{H}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|p-\lambda^{p}{}_{H}\|_{L^{\infty}(L^{2})}$}\tabularnewline
		\hline
		$1/4$ & 6.72e-01  & rate & 3.92e-01  & rate & 5.92e-02  & rate & 7.55e-01  & rate & 9.70e-02  & rate\tabularnewline
		\cline{3-3} \cline{5-5} \cline{7-7} \cline{9-9} \cline{11-11}
		$1/16$ & 8.20e-02  & 1.52  & 8.36e-02  & 1.11 & 1.48e-02  & 1.00  & 4.82e-02  & 1.99  & 6.83e-03  & 1.91\tabularnewline
		$1/64$ & 7.03e-03  & 1.77  & 1.92e-02  & 1.06 & 3.70e-03  & 1.00  & 3.31e-03  & 1.93 & 5.91e-04  & 1.77\tabularnewline
		\hline
	\end{tabular}
\end{table}

\begin{table}[h]
	\captionsetup{justification=centering}
	\caption{\label{tab:Example-1-lin-conv_c0_small}Example 1, convergence table
		for linear mortar with $H=Ch$, $\Delta t=10^{-4}$ and $c_{0}=10^{-3}$.}

	\begin{centering}
		\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
			\hline
			$h$ & \multicolumn{2}{c|}{\# GMRES} & \multicolumn{2}{c|}{$\|\sigma-\sigma_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\dvr\,(\sigma-\sigma_{h})\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\gamma-\gamma_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|u-u_{h}\|_{L^{\infty}(L^{2})}$}\tabularnewline
			\hline
			$1/4$ & 16  & rate & 1.25e-01  & rate & 6.09e-01  & rate & 1.39e+00  & rate & 5.78e-01  & rate\tabularnewline
			\cline{3-3} \cline{5-5} \cline{7-7} \cline{9-9} \cline{11-11}
			$1/8$ & 29  & -0.86  & 3.30e-02  & 1.92  & 3.11e-01  & 0.97  & 7.07e-01  & 0.97  & 2.92e-01  & 0.99 \tabularnewline
			$1/16$ & 50  & -0.79  & 8.34e-03  & 1.98  & 1.56e-01  & 0.99  & 3.55e-01  & 0.99  & 1.46e-01  & 1.00 \tabularnewline
			$1/32$ & 87  & -0.80  & 2.09e-03  & 1.99  & 7.82e-02  & 1.00  & 1.78e-01  & 1.00  & 7.31e-02  & 1.00 \tabularnewline
			$1/64$ & 157  & -0.85  & 5.38e-04  & 1.96  & 3.91e-02  & 1.00  & 8.89e-02  & 1.00  & 3.65e-02  & 1.00 \tabularnewline
			\hline
		\end{tabular}
		\par\end{centering}
	\centering{}%
	\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
		\hline
		$h$ & \multicolumn{2}{c|}{$\|z-z_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\dvr\,(z-z_{h})\|_{L^{2}(L^{2})}$} & \multicolumn{2}{c|}{$\|p-p_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|u-\lambda^{u}{}_{H}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|p-\lambda^{p}{}_{H}\|_{L^{\infty}(L^{2})}$}\tabularnewline
		\hline
		$1/4$ & 4.18e+01  & rate & 2.31e+00  & rate & 8.81e-01  & rate & 7.52e-01  & rate & 8.48e+00  & rate\tabularnewline
		\cline{3-3} \cline{5-5} \cline{7-7} \cline{9-9} \cline{11-11}
		$1/8$ & 9.68e+00  & 2.11  & 7.14e-01  & 1.69  & 2.33e-01  & 1.92  & 1.90e-01  & 1.98  & 2.11e+00  & 2.00\tabularnewline
		$1/16$ & 2.31e+00  & 2.07  & 2.00e-01  & 1.84  & 5.93e-02  & 1.98  & 4.77e-02  & 1.99  & 5.08e-01  & 2.06\tabularnewline
		$1/32$ & 5.68e-01  & 2.02  & 6.02e-02  & 1.73  & 1.62e-02  & 1.87  & 1.19e-02  & 2.00  & 1.25e-01  & 2.02\tabularnewline
		$1/64$ & 1.42e-01  & 2.00  & 2.22e-02  & 1.44  & 5.22e-03  & 1.64  & 2.98e-03  & 2.00  & 3.12e-02  & 2.00\tabularnewline
		\hline
	\end{tabular}
\end{table}

\begin{table}[h]
	\captionsetup{justification=centering}
	\caption{\label{tab:Example-1-quad-conv-small-c0}Example 1, convergence table
		for quadratic mortar with $H=C\sqrt{h}$, $\Delta t=10^{-4}$ and
		$c_{0}=10^{-3}$.}

	\begin{centering}
		\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
			\hline
			$h$ & \multicolumn{2}{c|}{\# GMRES} & \multicolumn{2}{c|}{$\|\sigma-\sigma_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\dvr\,(\sigma-\sigma_{h})\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\gamma-\gamma_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|u-u_{h}\|_{L^{\infty}(L^{2})}$}\tabularnewline
			\hline
			$1/4$ & 23  & rate & 1.28e-01  & rate  & 6.09e-01  & rate  & 1.39e+00  & rate  & 5.79e-01  & rate \tabularnewline
			\cline{3-3} \cline{5-5} \cline{7-7} \cline{11-11}
			$1/16$ & 41  & -0.41  & 8.39e-03  & 1.97  & 1.56e-01  & 0.98 & 3.55e-01  & 0.96 & 1.46e-01  & 0.99\tabularnewline
			$1/64$ & 72  & -0.41  & 5.61e-04  & 1.95  & 3.91e-02  & 1.00  & 8.89e-02  & 1.00  & 3.65e-02  & 1.00 \tabularnewline
			\hline
		\end{tabular}
		\par\end{centering}
	\centering{}%
	\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
		\hline
		$h$ & \multicolumn{2}{c|}{$\|z-z_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|\dvr\,(z-z_{h})\|_{L^{2}(L^{2})}$} & \multicolumn{2}{c|}{$\|p-p_{h}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|u-\lambda^{u}{}_{H}\|_{L^{\infty}(L^{2})}$} & \multicolumn{2}{c|}{$\|p-\lambda^{p}{}_{H}\|_{L^{\infty}(L^{2})}$}\tabularnewline
		\hline
		$1/4$ & 4.24e+01  & rate & 2.42e+00  & rate  & 9.97e-01  & rate  & 7.57e-01  & rate  & 1.07e+01  & rate\tabularnewline
		\cline{3-3} \cline{5-5} \cline{7-7} \cline{9-9} \cline{11-11}
		$1/16$ & 2.33e+00  & 2.01  & 2.01e-01  & 1.79 & 6.01e-02  & 2.06 & 4.83e-02  & 1.98 & 5.17e-01  & 2.19\tabularnewline
		$1/64$ & 1.50e-01  & 1.97  & 2.25e-02  & 1.58 & 5.40e-03  & 1.74 & 3.26e-03  & 1.95 & 3.38e-02  & 1.97\tabularnewline
		\hline
	\end{tabular}
\end{table}

\begin{figure}[h]
	\includegraphics[width=0.33\columnwidth]{\string"figures/Example 1/Stress_x\string".jpg}
	\includegraphics[width=0.33\columnwidth]{\string"figures/Example 1/Stress_Y\string".jpg}
	\includegraphics[width=0.33\columnwidth]{\string"figures/Example 1/displacement\string".jpeg}

	\includegraphics[width=0.33\columnwidth]{\string"figures/Example 1/rotation\string".jpeg}
	\includegraphics[width=0.33\columnwidth]{\string"figures/Example 1/velocity\string".jpeg}
	\includegraphics[width=0.33\columnwidth]{\string"figures/Example 1/pressure\string".jpeg}
	\captionsetup{justification=centering}
	\caption{\label{fig:Example-1,-computed_exact_mortar}Example 1, computed solution
		at final time step using a linear mortar on non-matching subdomain
		grids, top: stress x (left), stress y (middle), displacement (right), bottom: rotation (left), velocity (middle), pressure (right). Mesh size, $h=1/32$, $\Delta t=10^{-3}$ and $c_{0}=1.0$.}
\end{figure}

The numerical results that we observe are consistent with the theoretical
results from the previous sections. In particular,
we demonstrate the stability of the method over a $100$ time steps,
and Tables \ref{tab:Example-1-lin-conv} and \ref{tab:Example-1-quad-conv}
confirm convergence rates that follow from Theorem \ref{thm:Error-theorem}
and Table \ref{tab:Degree-poly-numerics}. With linear mortar $m=1$ and $H=Ch$, the interface error is ${\cal O}\left(h^{\frac{3}{2}}\right)$. With quadratic mortar $m=2$ and $H=\sqrt{h}$, the interface error is ${\cal O}\left(h^{\frac{5}{4}}\right)$. In both the cases, it is dominated by the subdomain error, which is ${\cal O}(h)$. As a result, we expect at least ${\cal O}\left(h\right)$ convergence
in both cases, which is what we observe. Comparison of the number
of interface iterations required in the case of linear and quadratic
mortars in Tables \ref{tab:Example-1-lin-conv} and \ref{tab:Example-1-quad-conv},
respectively shows that both mortar degrees result in similar accuracy
for the same level of subdomain mesh refinement. This is despite the
fact that the quadratic mortar case requires smaller number of interface
iterations compared to the linear mortar case with the same level
of subdomain mesh refinement. This is due to the choice of a coarser
mortar mesh in the case of quadratic mortar case. This points towards
a way to decrease the number of interface iterations by using a coarser
mesh and higher mortar space degree, without any loss in accuracy.
Tables \ref{tab:Example-1-lin-conv_c0_small}$\--$\ref{tab:Example-1-quad-conv-small-c0}
confirm that the stability and error bounds proved in previous sections
are not affected by smaller values of $c_{0}$. Further, Figure \ref{fig:Example-1,-computed_exact_mortar}
demonstrates the efficacy of the method in enforcing continuity of
solution variables across subdomain interfaces, weakly using coarse
mortar spaces. In fact, the solution looks almost identical to the
ones obtained using matching-subdomain grids in \cite{dd-biot},
with a smaller number of interface iterations for same level of finest
subdomain mesh refinement. This demonstrates the advantage of using
the multiscale technique we have developed over the completely matching
case that was discussed in \cite{dd-biot}.

\subsection{Example 2: heterogeneous medium}

In this example, we demonstrate the performance of our method in a
practical application with highly heterogeneous medium. First, we
compare the efficiency of our multiscale mortar method, where $H>h$,
with a fine scale method, where $H=h$. We expect the former to be
more efficient than the latter because of weaker enforcement of continuity
across subdomain interfaces using a coarser mortar space in the case
of the multiscale method. We then study the computational advantage
of using a multiscale stress-flux basis (MSB) over not using an MSB.
In the case of no-MSB, the number of subdomain solves is total \#GMRES
iterations across all time steps + $2\times$number of time steps,
where the last term comes from two extra solves required to solve the
system (\ref{eq:mortar-bar-1-1})$\--$(\ref{eq:mortar-star-5-1}) initially
and recovering the final solution once the GMRES converges. Similarly,
in the case of using MSB, total number of subdomain solves equals
the $\text{dim}(\Lambda_H)$+$2\times$number of time steps. Note that the
first term in the number of solves in the case of no-MSB is directly
proportional to the number of time steps used in time discretization,
while the same in the case of MSB method is independent of the number
of time steps used. This leads to MSB method being far more efficient
than the no-MSB method with any choice of mortar, as long as enough
number of time steps are used.

To obtain the desired level of heterogeneity in the medium, we use
the porosity and the permeability data from the Society of Petroleum Engineers
10th Comparative Solution Project (SPE10)\footnote{https://www.spe.org/web/csp/datasets/set02.htm}.
The porosity and permeability fields are given on a $60\times220$
grid and we use the rectangular region $(0,60)\times(0,220)$ as the
computational domain. We decompose the global domain into $3\times5$
subdomains consisting of identical rectangular blocks. The Young's
modulus is obtained from the porosity field data using the relation
$E=10^{2}\left(1-\frac{\phi}{c}\right)^{2.1},$ where $c=0.5$, refers
to the porosity at which the Young's modulus vanishes, see \cite{kovavcik1999correlation}
for details. These input fields are presented in Figure \ref{fig:Example-2,-input_fields_biot_mortar}.
We use parameters and boundary conditions as mentioned in Table \ref{tab:Physical-Parameters_ex2_mortar},
along with zero source terms. These conditions describe a flow from
left to right, driven by the gradient in the pressure. We use a compatible
initial condition for pressure, $p_{0}=1-x$. To obtain the essential
discrete initial data, we take the elliptic projection of the continuous
initial data, see (\ref{eq:well-pose-initial-data}). In particular,
we set $p_{h}^{0}$ to be the $L^{2}$-projection of $p_{0}$ onto
$W_{h}$ and solve a mixed elasticity domain decomposition problem
at $t=0$ to obtain $\sigma_{h}^{0}$. We also obtain $u_{h}^{0}$,
$\gamma_{h}^{0}$, and $\lambda_{H}^{u,0}$ from this solve.

We use a global $60\times220$ grid and solve the problem using both
fine scale ($H=h)$ and coarse $(H>h)$ mortar spaces. For the coarse
mortar case, we use both linear $(m=1)$ and quadratic mortars $(m=2)$
with one and two mortars per subdomain interface. The comparison of
the computed solution using different choices of mortars is given
in Figures \ref{fig:Example-2,-solution_plots-1}$\--$\ref{fig:Example-2,-solution_plots_last}.
Comparison of the number of solves required for different choices
of mortar, both in the no-MSB and MSB cases are given in Table \ref{tab:Example-2,-=000023GMRES_table}.
We report the number of subdomain solves which dominates the computational
complexity of the method.

Table \ref{tab:Example-2,-=000023GMRES_table} clearly shows that
using the multiscale mortar method requires fewer number of solves compared
to fine scale method and hence the former is computationally less
expensive than the latter. Comparison of the computed solution for
various choices of mortars in Figures \ref{fig:Example-2,-solution_plots-1}$\--$\ref{fig:Example-2,-solution_plots_last}
shows that we retain good amount of accuracy even in the case of the
coarsest mortar case with one linear mortar per interface. We also
note that using a single quadratic mortar per interface yields almost
identical results as the fine scale solution which emphasizes our
observation from the previous example that a coarse mortar can be
compensated by choosing a higher degree for mortar space $\Lambda_H$. Table
\ref{tab:Example-2,-=000023GMRES_table} also demonstrates the superiority
of using MSB for a time-dependent multiscale problem like the one
in our case. The number of solves in the case of no MSB is at least
an order of magnitude bigger than the MSB case which implies that
the construction of MSB is an excellent tool to make our multiscale
mortar method even more efficient than it already is compared to the
fine scale methods discussed in \cite{dd-biot}.

\begin{table}[h]
	\centering{}
	\caption{Example 2, parameters (top) and boundary conditions (bottom).\label{tab:Physical-Parameters_ex2_mortar}} \begin{tabular}{c|c}
		Parameter & Value\tabularnewline
		\hline
		Mass storativity $(c_{0})$ & $1.0$\tabularnewline
		Biot-Willis constant $(\alpha)$ & $1.0$\tabularnewline
		Time step $(\Delta t)$ & $10^{-3}$\tabularnewline
		Total time $(T)$ & $0.1$\tabularnewline
	\end{tabular}
	\hspace{5 mm}

	\begin{tabular}{c|c|c|c|c}
		Boundary & $\sigma$ & $u$ & $p$ & $z$\tabularnewline
		\hline
		Left & $\sigma n=-\alpha p n$ & - & $1$ & -\tabularnewline
		Bottom & $\sigma n=0$ & - & - & $z\cdot n=0$\tabularnewline
		Right & $-$ & $0$ & $0$ & -\tabularnewline
		Top & $\sigma n=0$ & - & - & $z\cdot n=0$\tabularnewline
	\end{tabular}
\end{table}

\begin{table}[h]
	\caption{Example 2, \#GMRES iterations and maximum number of subdomain solves.\label{tab:Example-2,-=000023GMRES_table}}

	\centering{}%
	\begin{tabular}{|c|c|c|c|c|}
		\hline
		mortar & Average \#GMRES & Total \#GMRES & \multicolumn{2}{c|}{Total \#Solves}\tabularnewline
		\hline
		&  &  & No MSB & MSB\tabularnewline
		\cline{4-5} \cline{5-5}
		linear fine scale & 343 & 34375 & 34575 & 968\tabularnewline
		1 linear per interface & 41 & 4149 & 4349 & 224\tabularnewline
		1 quadratic per interface & 61 & 6184 & 6384 & 236\tabularnewline
		2 linear per interface & 80 & 8010 & 8210 & 248\tabularnewline
		2 quadratic per interface & 123 & 12302 & 12502 & 272\tabularnewline
		\hline
	\end{tabular}
\end{table}

\begin{figure}[h]
	\includegraphics[width=0.33\columnwidth]{figures/Example_2/permeability.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/porosity.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/youngs_modulus.jpeg}\caption{Example 2, permeability, porosity, Young's modulus. \label{fig:Example-2,-input_fields_biot_mortar}}
\end{figure}

\begin{figure}[h]
	\includegraphics[width=0.33\columnwidth]{figures/Example_2/pressure_fine_scale.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/pressure_linear_one.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/pressure_linear_two.jpeg}\captionsetup{justification=centering}\caption{Example 2, pressure (color) and velocity (arrows): fine scale (left),
		single linear mortar per interface (middle), and two linear mortars
		per interface (right).\label{fig:Example-2,-solution_plots-1}}
\end{figure}

\begin{figure}[h]
	\centering{}\includegraphics[width=0.33\columnwidth]{figures/Example_2/pressure_quadratic_one.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/pressure_quadratic_two.jpeg}\captionsetup{justification=centering}\caption{Example 2, pressure (color) and velocity (arrows): single quadratic
		mortar per interface (left), and two quadratic mortars per interface
		(right).}
\end{figure}

\begin{figure}[h]
	\includegraphics[width=0.33\columnwidth]{figures/Example_2/velocity_fine_scale.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/velocity_linear_one.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/velocity_linear_two.jpeg}\captionsetup{justification=centering}\captionsetup{justification=centering}\caption{Example 2, velocity magnitude: fine scale (left), single linear mortar
		per interface (middle), and two linear mortars per interface (right).}
\end{figure}

\begin{figure}[h]
	\centering{}\includegraphics[width=0.33\columnwidth]{figures/Example_2/velocity_quadratic_one.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/velocity_quadratic_two.jpeg}\captionsetup{justification=centering}\caption{Example 2, velocity magnitude: single quadratic mortar per interface
		(left), and two quadratic mortars per interface (right).}
\end{figure}

\begin{figure}[h]
	\includegraphics[width=0.33\columnwidth]{figures/Example_2/displacement_fine_scale.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/displacement_linear_one.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/displacement_linear_two.jpeg}\captionsetup{justification=centering}\caption{Example 2, displacement vector (arrows) and its magnitude: fine scale
		(left), single linear mortar per interface (middle), and two linear
		mortars per interface (right).}
\end{figure}

\begin{figure}[h]
	\centering{}\includegraphics[width=0.33\columnwidth]{figures/Example_2/displacement_quadratic_one.jpeg}\includegraphics[width=0.33\columnwidth]{figures/Example_2/displacement_quadratic_two.jpeg}\captionsetup{justification=centering}\caption{Example 2, displacement vector (arrows) and its magnitude: single
		quadratic mortar per interface (left), and two quadratic mortars per
		interface (right).\label{fig:Example-2,-solution_plots_last}}
\end{figure}

\section{Conclusions}

In this paper we presented a multiscale mortar mixed finite element
technique (MMMFE) for the Biot system of poroelasticity in a five-filed
fully mixed formulation. This method is the generalization of the
monolithic domain decomposition technique discussed in the previous
paper, with the extra capability to use non-matching subdomain grids
at the interface. This capability is obtained by using composite multiscale
mortar Lagrange multiplier spaces that approximates displacement and pressure
on a coarse mortar grid at the interface. The global problem can be
reduced into a series of parallel Dirichlet type problems and an interface
problem for the composite displacement-pressure Lagrange multiplier
spaces which requires subdomain solves at each iteration. We showed
the well-posedness and stability of the method under proper assumptions.
We have also carried out an extensive error analysis of the method
to get a combined a priori error estimate for all the unknowns in
the formulation. To complete the analysis, we have done a series of
numerical experiments to put the theory into test. We observed stability
and convergence results as predicted by the theory and also demonstrated
the application of the method to a highly heterogeneous medium. We
noted that in practice, a coarser mesh with higher mortar space degree
can be used to get a smaller interface problem and hence faster convergence
without compromising the accuracy of the method. We conclude the paper
by recalling the effectiveness of the construction and use of a pre-computed
multiscale stress-flux basis (MSB), which makes the MMMFE method far
more superior than the fine scale monolithic methods, especially when
a coarse mortar mesh is used.

\bibliographystyle{abbrv}
\addcontentsline{toc}{section}{\refname}