\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 15, pp. 1--17.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2016/15\hfil Well-posedness of a coupled TCP system]
{Sufficient conditions for Hadamard well-posedness of a
coupled thermo-chemo-poroelastic system}

\author[T. Malysheva, L. W. White \hfil EJDE-2016/15\hfilneg]
{Tetyana Malysheva, Luther W. White}

\address{Tetyana Malysheva \newline
Department of Natural \& Applied Sciences \\
University of Wisconsin-Green Bay, Green Bay, WI 54311-7001, USA}
\email{malyshet@uwgb.edu}

\address{Luther W. White \newline
Department of Mathematics \\
University of Oklahoma, Norman, OK 73019-3103, USA}
\email{lwhite@ou.edu}

\thanks{Submitted December 15, 2015. Published January 8, 2016.}
\subjclass[2010]{35D30, 35E99, 35G16, 35Q74, 35Q86}
\keywords{Parabolic-elliptic system; poroelasticity; thermo-poroelasticity,
\hfill\break\indent thermo-chemo-poroelasticity; Hadamard well-posedness}

\begin{abstract}
 This article addresses the well-posedness of a coupled parabolic-elliptic
 system modeling fully coupled thermal, chemical, hydraulic, and mechanical
 processes in porous formations that impact drilling and borehole stability.
 The underlying thermo-chemo-poroelastic model is a system of time-dependent
 parabolic equations describing thermal, solute, and fluid diffusions coupled
 with Navier-type elliptic equations that attempt to capture the elastic
 behavior of rock around a borehole. An existence and uniqueness theory for
 a corresponding initial-boundary value problem is an open problem in the field.
 We give sufficient conditions for the well-posedness in the sense of Hadamard
 of a weak solution to a fully coupled parabolic-elliptic initial-boundary value
 problem describing homogeneous and isotropic media.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{assumption}[theorem]{Assumption}
\allowdisplaybreaks

 \section{Introduction}

In this article, we are concerned with the well-posedness of a coupled
parabolic-elliptic system arising in petroleum- and geothermal-related
applications of rock mechanics. This work is motivated by the problems
of drilling and borehole stability in porous formations that involve the
modeling of fully coupled thermal, chemical, hydraulic, and mechanical
(elastic deformation) processes.
		
	\subsection{Literature review}
The poroelasticity theory describing the coupled processes of elastic deformation
and pore fluid diffusion in fluid-saturated isothermal porous media can
 be traced back to the pioneering works of von Terzaghi
\cite{Terzaghi1, Terzaghi2} and Biot \cite{Biot1}.
Biot \cite{Biot2} pointed out a complete mathematical analogy between
poroelasticity and thermoelasticity with the temperature playing the
same role as the fluid pressure and heat conduction corresponding to fluid flow.
A complete well-posedness analysis of a general initial-boundary value problem
for a system of coupled partial differential equations that describes the Biot
consolidation model \cite{Biot1} in poroelasticity, as well as a coupled
quasi-static problem in thermoelasticity, has been carried out by
Showalter \cite{Showalter2, Showalter3}. Based on the theory of linear
degenerate evolution equations in Hilbert spaces, the existence and uniqueness
of strong and weak solutions, as well as regularity theory for the
initial-boundary value problem, were developed.
	
The Biot theory for isothermal systems was first extended by
 Schiffman \cite{Schiffman} to account for the effects of thermal expansion
of both the pore fluid and the elastic matrix. Since then, a substantial
literature on thermo-poroelasticity theory and the modeling of the coupled
hydro-thermo-mechanical behavior of a fluid-saturated porous media has
been developed for rock mechanics, including petroleum and geothermal
borehole stability
\cite{Abousleiman, Belotserkovets, Coussy, Diek2, Diek3, Diek4, Kodashima1,
Kodashima2, Kurashige, McTigue, Pao, Wang}.
Considerable attention has been recently placed on the impact of chemical
processes in porous media on drilling and borehole stability
\cite{Diek1, Ghassemi3, Rafieepour, Zhou}. Due to the complexity of
cross-coupling mechanisms involved in thermo-poroelastic and
chemo-thermo-poroelastic models, the general question of existence and
uniqueness of solutions to the corresponding initial-boundary value problems
remains open and very few analytical solutions are currently available.
Typically, the solutions are derived under the assumptions that some of
the couplings can be neglected
\cite{Coussy, Ghassemi3, Kodashima1, Kodashima2, McTigue, Rafieepour}.
An exact unique analytical solution for a specific case of the fully
coupled thermo-hydro-mechanical response of a fluid-saturated porous
sphere under mechanical pulse load was developed by Belotserkovets and
Prevost \cite{Belotserkovets} using the Laplace transformation and the
residue theorem.
	
The model that we study is based on equations derived by Diek \cite{Diek1}
and constitutes the general theory of fully coupled chemical
thermo-poroelasticity for porous media saturated by a compressible fluid.
The theory satisfies the first and second laws of thermodynamics and is
based on concepts of irreversible thermodynamics, a novel rock constitutive
relation, and Onsager's transport phenomenology \cite{Diek1}.
Our objective here is to establish the well-posedness theory for the
fully coupled thermo-chemo-poroelastic (TCP) system describing homogeneous
and isotropic fluid-sa\-turated porous media. The main result of this work is
Theorem \ref{main} which gives a sufficient condition for the well-posedness
in the sense of Hadamard of a weak solution to the coupled linear parabolic-elliptic
initial-boundary value problem describing the fully coupled TCP model.
		
		\subsection{Notation}
Let $\Omega \subset \mathbb{R}^n, ~n=2,3$, be a bounded open domain with a
 sufficiently smooth boundary $\Gamma$ and $\bar{\Omega}=\Omega \cup \Gamma$.
 We write $\mathbf{x}$ for $(x,y,z) \in \bar{\Omega} \subset \mathbb{R}^3$
or $(x,y) \in \bar{\Omega} \subset \mathbb{R}^2$. Let $\Gamma_F \subset \Gamma$
 and $\Gamma_F$ have a nonempty interior relative to $\Gamma$.
The following notation will be used for the spaces of vector-valued functions:
$\mathbb{H}^n=L^2(\Omega)^n$, $\mathbb{V}^n=H^1(\Omega)^n$,
$\mathbb{V}_0^n=H_0^1(\Omega)^n$ with the norm
\begin{equation*}
\|\mathbf{u}\|_{\mathbb{V}_0^n}=\Big[\sum_{k=1}^n\int_{\Omega} |\nabla u_k|^2
 \mathrm{d}\Omega\Big]^{1/2}, \quad \mathbf{u}= [u_1 \ldots u_n]^T
\end{equation*}
and $\tilde{\mathbb{V}}_0^n=\big\{ \boldsymbol{\varphi}\in \mathbb{V}^n:
\boldsymbol{\varphi}\big|_{\Gamma_F} =\mathbf{0}\big\}$
with the norm inherited from $\mathbb{V}^n$ or from
$\mathbb{V}_0^n$, $n \in \mathbb{N}$, as above. We denote
by $L^2(a,b;\mathrm{X})$ the space of $L^2$-integrable functions from
$[a,b]\subset \mathbb{R}$ into a Hilbert space $\mathrm{X}$ with the norm
\begin{equation*}
\|u\|_{L^2(a,b;\mathrm{X})}=\Big[\int_a^b\|u(t)\|_{\mathrm{X}}^2
\mathrm{d} t\Big]^{1/2}
\end{equation*}
We also introduce the Hilbert space
\begin{equation*}
\mathcal{W}(a,b;\mathrm{X})=\big\{u:u\in L^2(a,b;\mathrm{X}),\,
\dot{u}\in L^2(a,b;\mathrm{X}')\big\}
\end{equation*}
with the norm
\begin{equation*}
\|u\|_{\mathcal{W}}=\Big[\|u\|_{L^2(a,b;\mathrm{X})}^2
+ \|\dot{u}\|_{L^2(a,b;\mathrm{X}')}^2\Big]^{1/2}
\end{equation*}
where the superscript dot $(~\dot{}~)$ denotes a time derivative and
$\mathrm{X}'$ is the dual of $\mathrm{X}$.

Let $\mathbf{f}= [f_1 \ldots f_n]^T \in \mathbb{V}^n$. We will use the
notation $\nabla \mathbf{f}=[\nabla f_1 \ldots \nabla f_n]^T$ and
\begin{equation*}
(\nabla\mathbf{f}, \nabla\mathbf{g})_{\mathbb{H}^n}
=\sum_{k=1}^n(\nabla f_k, \nabla g_k)_{\mathbb{H}^n}, \quad
\mathbf{f}, \mathbf{g} \in\mathbb{V}^n
\end{equation*}
Throughout this article, the symbol $c$ will be used to denote the
Poincar\'e constant: $c=c(\Omega)>0$,
$\|\mathbf{u}\|_{\mathbb{H}^n}^2 \leq c\|\mathbf{u}\|_{\mathbb{V}_0^n}^2$,
$\mathbf{u}\in \mathbb{V}_0^n$.

\section{The coupled TCP model: underlying equations}
   \label{model}

	The underlying TCP model is a system of time-dependent parabolic partial
differential equations coupled with Navier-type elliptic equations with time,
 $t \in (0, t_f)$, as a parameter. The parabolic equations represent heat,
 solute, and fluid diffusions, and the Navier-type elliptic equations attempt
to capture the elastic behavior of rock, while incorporating thermal,
chemical, and porous media effects. The equations are developed in terms of
the following variables: the absolute temperature $T(\mathbf{x},t)$, the
solute mass fraction $C(\mathbf{x},t)$, the pore pressure $p(\mathbf{x},t)$,
 and the vector of solid displacements $\mathbf{u}(\mathbf{x},t)$.
	
The coupled partial differential equations supplemented by the appropriate
initial and boundary conditions constitute an initial-boundary-value problem
with constant coefficients defined in an open region $\Omega$ exterior
to the borehole. The region is specified with a sufficiently smooth boundary
$\Gamma$ such that, without loss of generality, we may assume that, on the
outer (far-field) boundary $\Gamma_F$, the absolute temperature, solute mass
fraction, pore pressure, and displacements are time-independent and displacements
 and their velocities are negligibly small. Specifically, the region is
an elliptical annulus in a two-dimensional case and a vertical or inclined
finite cylinder in a three-dimensional case.
	
	It is convenient to consider the thermal diffusion, solute diffusion,
and fluid diffusion as a system (diffusion system). Hence, we introduce
the vector $\mathbf{V}= [T ~C ~p]^T$ and, with this notation, the
initial-boundary value problem describing the underlying TCP model has the form
\begin{gather}
M \dot{\mathbf{V}} - A \nabla^{2} \mathbf{V}
=-\mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}), \quad
\text{in } \Omega \times (0,t_f), \label{eq1}\\
%
\Big(K+\frac{G}{3}\Big)\nabla(\nabla\cdot\dot{\mathbf{u}})+G\nabla^2
\dot{\mathbf{u}}=\nabla(\mathbf{b}_1 \cdot \dot{\mathbf{V}}), \quad \text{in }
\Omega \times (0,t_f) \label{eq2}
\end{gather}
with boundary conditions
\begin{gather}
\mathbf{V}(\mathbf{x},t)= \mathbf{V}_B(\mathbf{x},t), \quad \text{on }
\Gamma \times (0,t_f) \label{eq3} \\
%
\mathbf{u}(\mathbf{x},t) \approx \mathbf{0}, \quad \text{on }
\Gamma_F \times (0, t_f) \label{eq4}\\
%
\dot{\mathbf{u}}(\mathbf{x},t) \approx \mathbf{0}, \quad \text{on }
\Gamma_F \times (0, t_f) \label{eq5}\\
%
\boldsymbol{\dot{\tau}} \mathbf{n}=\big((\mathbf{b}_1 \cdot
\dot{\mathbf{V}})I+\boldsymbol{\dot{\hat{\sigma}}}\big)\mathbf{n},
\quad \text{on } \Gamma \setminus \Gamma_F \times (0,t_f) \label{eq6}
\end{gather}
and initial conditions
\begin{equation}
\mathbf{V}(\mathbf{x},0) =
\begin{cases} 
 \mathbf{V}_I(\mathbf{x}), &\quad \text{in }\Omega  \\
 \mathbf{V}_B(\mathbf{x},0), &\quad \text{on }\Gamma
\end{cases} \label{eq7}
\end{equation}
%
where $M$ and $A$ are $3\times3$ matrices of diffusion coefficients,
$\mathbf{b}_0$ and $\mathbf{b}_1$ are $3\times1$ coupling vectors
determined by physical properties and input parameters of the rock/fluid system,
$K$ and $G$ are the bulk and shear moduli, respectively,
 $\boldsymbol{\tau}$ is the stress tensor, $\mathbf{n}$ is the outward unit
normal vector on the boundary, $I$ is the $n \times n$ identity matrix,
 $n=2$ or $3$ is the dimension of the problem, and $\boldsymbol{\hat{\sigma}}$
is the applied boundary stress tensor. The boundary function $\mathbf{V}_B$
is spatially independent on the inner (borehole) boundary $\Gamma_B$ and
time-independent on the far-field boundary $\Gamma_F$.

\begin{remark} \label{rmk2.1} \rm
The matrices $M=[m_{ij}]$ and $A=[a_{ij}], ~i,j=1,2,3$, of diffusion
coefficients are non-symmetric with the entries specified as follows:
$m_{21}=m_{23}=0$ and no other entries of $ M$ are zero, and
$ a_{ij}\neq 0$, $i,j = 1, 2, 3$. The second entry in the coupling vector
$\mathbf{b}_0$ is zero and all the entries of the coupling vector $\mathbf{b}_1$
are non-zero. This implies that the TCP model is non-symmetric, and the
diffusion equation \eqref{eq1} cannot be rescaled to make the coupling
vectors equal. Therefore, the methods presented in \cite{Showalter3}
are not applicable to the fully coupled TCP problem \eqref{eq1}-\eqref{eq7}.
\end{remark}

We begin with a lemma that provides an alternative formulation of the
Navier-type elastic system \eqref{eq2}, \eqref{eq4}-\eqref{eq6}.
Its proof, based on the principle of minimum total potential energy,
will play an important role in well-posedness analysis presented
in the next section.

\begin{lemma} \label{lemma1}
The Navier-type elastic system \eqref{eq2}, \eqref{eq4}-\eqref{eq6}
is equivalent to the system
\begin{gather}
\big(K+\frac{G}{3}\big)\nabla(\nabla \cdot \mathbf{u})+G\nabla^2\mathbf{u}
= \nabla(\mathbf{b}_1 \cdot \mathbf{V}), \quad
\text{in } \Omega \times (0,t_f) \label{eq8} \\
\mathbf{u} \approx \mathbf{0}, \quad \quad \quad
\text{on } \Gamma_F \times (0,t_f) \tag{\ref{eq4}}\\
\boldsymbol{\tau} \mathbf{n}=\big((\mathbf{b}_1 \cdot
\mathbf{V})I+\boldsymbol{\hat{\sigma}}\big)\mathbf{n}, \quad
\text{on } \Gamma \setminus \Gamma_F \times (0,t_f) \label{eq9}
\end{gather}
\end{lemma}

\begin{proof}
For the sake of simplicity, we restrict the proof to the two-dimensional
case $\Omega \subset \mathbb{R}^2.$ Directly analogous arguments can be
developed for the three-dimensional case. By the principle of minimum
total potential energy, the region $\Omega$ shall displace to a position
 that minimizes the total potential energy of the elastic system,
\begin{equation}
\mathcal{V}(\mathbf{u})=\mathcal{V}_S(\mathbf{u})-W_b(\mathbf{u})-W_S(\mathbf{u})
 \label{eq10}
\end{equation}
where $\mathcal{V}_S(\mathbf{u})$ is the elastic energy of the system;
$W_b(\mathbf{u})$ is work done by body forces due to the absolute temperature,
solute mass fraction, and pore pressure; and $W_S(\mathbf{u})$ is work done
by applied boundary stress.

Now we will specify $\mathcal{V}_S(\mathbf{u})$, $W_b(\mathbf{u})$, and
$W_S(\mathbf{u})$. Here and in the following we will suppress the time
 dependence of the displacement vector $\mathbf{u}$ for convenience.
Given the displacement
$\mathbf{u}(\mathbf{x})=u(\mathbf{x})\mathbf{i}+v(\mathbf{x})\mathbf{j}$,
the linearized strain is the second order symmetric tensor
\begin{equation*}
\boldsymbol{\varepsilon}(\mathbf{u})
=\frac{1}{2}\big(\nabla\mathbf{u}+\nabla\mathbf{u}^T\big)
\end{equation*}
where
\begin{equation*}
\varepsilon_{11}=u_x, \quad \varepsilon_{12}=\varepsilon_{21}
=\frac{1}{2}\big(u_y+v_x\big), \quad \varepsilon_{22}=v_y
\end{equation*}
In linear elasticity, a relation between the stress tensor
$\boldsymbol{\tau}$ and the linearized strain tensor
$\boldsymbol{\varepsilon}(\mathbf{u})$ is given by the generalized Hooke's law
\begin{equation}
\tau_{ij}(\mathbf{u})=a_{ijkl}\varepsilon_{kl}(\mathbf{u}) \label{eq11}
\end{equation}
where Einstein summation convention is used, indices run from $1$ to $n$,
and $a_{ijkh}$ are the coefficients of elasticity independent of the strain
tensor with the properties of symmetry
\begin{equation}
a_{ijkh}=a_{jihk}=a_{khij} \label{eq12}
\end{equation}
and of ellipticity: there exists a constant $ \alpha>0$ such that
\begin{equation*}
a_{ijkh}\varepsilon_{ij}\varepsilon_{kh}
\geq \alpha \varepsilon_{ij}\varepsilon_{ij}, \quad \forall \varepsilon_{ij}.
\end{equation*}
Specifically, for a homogeneous elastic isotropic medium, the stress-strain
relation in terms of the bulk modulus $K$ and the shear modulus $G$ has the form
\begin{equation}
\boldsymbol{\tau} = 2G\boldsymbol{\varepsilon}
+ \big(K-\frac{2G}{3}\big)(\text{tr}\boldsymbol{\varepsilon})I \label{eq13}
\end{equation}
With this notation, the elastic strain energy of the system is
\begin{equation}
\mathcal{V}_S(\mathbf{u}) = \frac{1}{2}\int_\Omega \tau_{11}(\mathbf{u})
\varepsilon_{11}(\mathbf{u}) + 2\tau_{12}(\mathbf{u})\varepsilon_{12}
(\mathbf{u})+\tau_{22}(\mathbf{u})\varepsilon_{22}(\mathbf{u})\mathrm{d} \Omega;
 \label{eq14}
\end{equation}
the work done by body forces due to the absolute temperature, solute mass
fraction, and pore pressure is
\begin{equation}
W_b(\mathbf{u})=\int_\Omega (\mathbf{b}_1 \cdot
\mathbf{V})\big(\text{tr}\boldsymbol{\varepsilon}(\mathbf{u})\big) \mathrm{d}
 \Omega; \label{eq15}
\end{equation}
and the work done by applied boundary stress is
\begin{equation}
W_S(\mathbf{u})=\int_{\Gamma}(\boldsymbol{\hat{\sigma}}\mathbf{n})
\cdot\mathbf{u} \mathrm{d} \Gamma \label{eq16}
\end{equation}
At this point, we wish to express the elastic strain energy
$\mathcal{V}_S(\mathbf{u})$ as a functional on $\mathbb{V}^2$.
Let $\boldsymbol{\Phi}(\mathbf{x})=\phi(\mathbf{x})\mathbf{i}
+\psi(\mathbf{x})\mathbf{j}$. Define the bilinear form
$a_E:\mathbb{V}^2 \times \mathbb{V}^2 \to \mathbb{R}$ by
\begin{equation}
a_E(\mathbf{u},\boldsymbol{\Phi}) = \int_\Omega \tau_{11}(\mathbf{u})
\varepsilon_{11}(\boldsymbol{\Phi})+2\tau_{12}(\mathbf{u})
\varepsilon_{12}(\boldsymbol{\Phi}) +\tau_{22}(\mathbf{u})
\varepsilon_{22}(\boldsymbol{\Phi}) \mathrm{d} \Omega,  \label{eq17}
\end{equation}
for $\mathbf{u}, \boldsymbol{\Phi} \in \mathbb{V}^2$.
From \eqref{eq10} and \eqref{eq14}-\eqref{eq17}, the total potential
 energy of the system has the form
\begin{equation}
\mathcal{V}(\mathbf{u})=\frac{1}{2}a_E(\mathbf{u},\mathbf{u})
-\int_\Omega (\mathbf{b}_1 \cdot \mathbf{V}) \nabla
\cdot \mathbf{u} \mathrm{d} \Omega
- \int_{\Gamma}(\boldsymbol{\hat{\sigma}}\mathbf{n})\cdot\mathbf{u}
\mathrm{d} \Gamma \label{eq18}
\end{equation}
Define the following vectors:
\begin{equation*}
\boldsymbol{\tau}_1=[\tau_{11} \quad \tau_{12}]^T, \quad
\boldsymbol{\tau}_2=[\tau_{21} \quad \tau_{22}]^T
\end{equation*}
Then
\begin{align*}
a_E(\mathbf{u},\boldsymbol{\Phi})
=& \int_\Omega \tau_{11}(\mathbf{u})\phi_x+\tau_{12}(\mathbf{u})\phi_y
 +\tau_{21}(\mathbf{u})\psi_x+\tau_{22}(\mathbf{u})\psi_y \mathrm{d} \Omega\\
=& \int_\Omega \boldsymbol{\tau}_1(\mathbf{u}) \cdot \nabla \phi
 + \boldsymbol{\tau}_2(\mathbf{u}) \cdot \nabla \psi \mathrm{d} \Omega \\
=&\int_{\Gamma} \big(\boldsymbol{\tau}(\mathbf{u})\mathbf{n}\big)
 \cdot \boldsymbol{\Phi}\mathrm{d} \Gamma
 - \int_\Omega [\nabla \cdot \boldsymbol{\tau}_1
\quad \nabla \cdot \boldsymbol{\tau}_2]^T(\mathbf{u})
 \cdot \boldsymbol{\Phi} \mathrm{d} \Omega
\end{align*}
and we obtain the  Green's formula:
\begin{equation}
a_E(\mathbf{u},\boldsymbol{\Phi})=\int_{\Gamma}
(\boldsymbol{\tau}(\mathbf{u})\mathbf{n}) \cdot \boldsymbol{\Phi}\mathrm{d} \Gamma
- \int_\Omega [\nabla \cdot \boldsymbol{\tau}_1 \quad
\nabla \cdot \boldsymbol{\tau}_2]^T(\mathbf{u}) \cdot \boldsymbol{\Phi}
\mathrm{d} \Omega, \quad \forall \mathbf{u}, \boldsymbol{\Phi}
\in \mathbb{V}^2  \label{eq19}
\end{equation}
%
Referring to \eqref{eq4} and \eqref{eq5}, the set of admissible displacements,
in general, is
\begin{equation}
\mathcal{U}_{ad}=\big\{\mathbf{u} \in \tilde{\mathbb{V}}_0^n:
\dot{\mathbf{u}} \in \tilde{\mathbb{V}}_0^n\big\}, \quad n=2 \text{ or }
 3 \label{eq20}
\end{equation}
Using the principle of minimum total potential energy, the displacement
$\mathbf{u} \in \mathcal{U}_{ad}$ that the region $\Omega \subset \mathbb{R}^2$
undergoes is given by
\begin{equation}
D\mathcal{V}(\mathbf{u})\boldsymbol{\Phi}=0, \quad \forall
\boldsymbol{\Phi} \in \tilde{\mathbb{V}}_0^2 \label{eq21}
\end{equation}
where
\begin{equation*}
D\mathcal{V}(\mathbf{u})\boldsymbol{\Phi}=\frac{d}{d\delta}
 \mathcal{V}(\mathbf{u}+\delta\boldsymbol{\Phi})\Big\vert_{\delta=0}
\end{equation*}
is the G\^{a}teaux differential of $\mathcal{V}$ with increment
$\boldsymbol{\Phi}$.
From \eqref{eq18} and \eqref{eq21},
\begin{equation}
a_E(\mathbf{u},\boldsymbol{\Phi})-\int_\Omega (\mathbf{b}_1
\cdot \mathbf{V})\nabla \cdot \boldsymbol{\Phi}\mathrm{d} \Omega
- \int_{\Gamma \setminus \Gamma_F}(\boldsymbol{\hat{\sigma}}
\mathbf{n})\cdot\boldsymbol{\Phi}\mathrm{d} \Gamma = 0,
\quad \forall \boldsymbol{\Phi} \in \tilde{\mathbb{V}}_0^2 \label{eq22}
\end{equation}
and applying the divergence theorem we have
\begin{equation}
\begin{aligned}
&a_E(\mathbf{u},\boldsymbol{\Phi})-\int_{\Gamma \setminus \Gamma_F}
\big((\mathbf{b}_1 \cdot \mathbf{V})I\mathbf{n}) \cdot \boldsymbol{\Phi} \mathrm{d}
\Gamma \\
&+ \int_\Omega \nabla (\mathbf{b}_1 \cdot \mathbf{V}) \cdot
\boldsymbol{\Phi} \mathrm{d} \Omega
 - \int_{\Gamma \setminus \Gamma_F}(\boldsymbol{\hat{\sigma}}
\mathbf{n})\cdot\boldsymbol{\Phi}\mathrm{d} \Gamma = 0,
\quad \forall \boldsymbol{\Phi} \in \tilde{\mathbb{V}}_0^2
\end{aligned}\label{eq23}
\end{equation}
Green's formula \eqref{eq19} and \eqref{eq23} yield
\begin{gather}
[\nabla \cdot \boldsymbol{\tau}_1 \quad \nabla \cdot \boldsymbol{\tau}_2]^T
= \nabla (\mathbf{b}_1 \cdot \mathbf{V}), \quad \text{in } \Omega \label{eq24} \\
\boldsymbol{\tau} \mathbf{n}=\big((\mathbf{b}_1 \cdot
\mathbf{V})I+\boldsymbol{\hat{\sigma}}\big)\mathbf{n}, \quad \text{on }
\Gamma \setminus \Gamma_F, \label{eq25}
\end{gather}
and the stress-strain relation \eqref{eq13} gives
\begin{equation}
[\nabla \cdot \boldsymbol{\tau}_1 \quad \nabla
\cdot \boldsymbol{\tau}_2]^T = \big(K+\frac{G}{3}\big)
\nabla(\nabla \cdot \mathbf{u})+G\nabla^2\mathbf{u} \label{eq26}
\end{equation}
From \eqref{eq24}-\eqref{eq26} and \eqref{eq4}, we obtain the system
 \eqref{eq8}, \eqref{eq9}, and \eqref{eq4}.

On the other hand,  since $\dot{\mathbf{u}} \in \tilde{\mathbb{V}}_0^2$,
differentiating \eqref{eq22} with respect to time, we have
\begin{equation*}
a_E(\dot{\mathbf{u}},\boldsymbol{\Phi})-\int_\Omega (\mathbf{b}_1 \cdot
\dot{\mathbf{V}})\nabla \cdot \boldsymbol{\Phi}\mathrm{d} \Omega
 - \int_{\Gamma \setminus \Gamma_F}(\dot{\hat{\boldsymbol{\sigma}}}
\mathbf{n})\cdot\boldsymbol{\Phi}\mathrm{d}\Gamma = 0, \quad
\forall \boldsymbol{\Phi} \in \tilde{\mathbb{V}}_0^2
\end{equation*}
Using the same argument as above leads to the equivalence of the systems
\eqref{eq8}, \eqref{eq9}, \eqref{eq4} and \eqref{eq2}, \eqref{eq4}-\eqref{eq6}.
\end{proof}

\begin{remark} \label{rmk2.3} \rm
From the condition of mechanical equilibrium
\begin{equation*}
\nabla \cdot \boldsymbol{\sigma}=0
\end{equation*}
where $ \boldsymbol{\sigma}=\boldsymbol{\tau}-(\mathbf{b}_1 \cdot \mathbf{V})I$
 is the poroelastic stress tensor, and from the stress-strain relation \eqref{eq13},
 it follows that \eqref{eq8} constitutes the equation of equilibrium.
\end{remark}

		  \section{Well-posedness of the diffusion and elastic systems}
   \label{analysis}

In this section, we discuss the well-posedness of the parabolic initial-boundary
value problem for the diffusion system and the elliptic initial-boundary
value problem for the elastic system considering coupling terms as data.

\subsection{Diffusion system}
The following assumptions will be made on the matrices of diffusion
coefficients and the boundary and initial functions.

\begin{assumption} \label{assm1} \rm
The matrix $M$ is non-singular and all the eigenvalues of the matrix
$M^{-1}A$ are positive: $\lambda_i>0$, $i=1,2,3$.
\end{assumption}

Under Assumption \ref{assm1}, the matrix $M^{-1}A$ admits the eigendecomposition
\begin{equation}
M^{-1}A = PDP^{-1} \label{eq27}
\end{equation}
where $P=[\mathbf{e}_1 \quad \mathbf{e}_2 \quad \mathbf{e}_3], ~\mathbf{e}_i$
is the eigenvector of $M^{-1}A$ corresponding to the eigenvalue $\lambda_i$,
$\|\mathbf{e}_i\|=1, ~i=1,2,3$, and
$D=\operatorname{diag}(\lambda_1, \lambda_2, \lambda_3)$.
	
\begin{assumption} \label{assm2} \rm
	$\mathbf{V}_B \in L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big)$,
$\dot{\mathbf{V}}_B \in L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big)$,
and $\mathbf{V}_I \in \mathbb{V}^3$.
\end{assumption}
	
We transform the parabolic diffusion system \eqref{eq1}, \eqref{eq3},
and \eqref{eq7} to the equivalent diagonalized system with homogeneous
boundary conditions. The inverse trace theorem \cite{Steinbach}
yields that there exists a continuous mapping
\begin{equation*}
\gamma_0^{-}:L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big) \to
 L^2\big(0,t_f; \mathbb{V}^3\big)
\end{equation*}
such that $\gamma_0\gamma_0^{-}w=w$,
$w\in L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big)$, where
\begin{gather*}
\gamma_0:L^2\big(0,t_f; \mathbb{V}^3\big) \to
 L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big)\\
\quad  \gamma_0(v)=v\big|_{\Gamma}, \quad v\in L^2\big(0,t_f; \mathbb{V}^3\big)
\end{gather*}
is a linear and continuous trace mapping. We define the vector
$\mathbf{W} \in L^2\big(0,t_f; \mathbb{V}^3\big)$ by
\begin{equation}
\mathbf{W}(\mathbf{x},t)=\gamma_0^{-}\mathbf{V}_B(\mathbf{x},t)  \label{eq28}
\end{equation}
and the vector
	\begin{equation}
	\mathbf{U}(\mathbf{x},t) = P^{-1}\big(\mathbf{V}(\mathbf{x},t)
-\mathbf{W}(\mathbf{x},t)\big) \label{eq29}
	\end{equation}
The transformation given by \eqref{eq27}-\eqref{eq29} leads to the
following initial-boundary value problem equivalent to the parabolic
diffusion system \eqref{eq1}, \eqref{eq3}, and \eqref{eq7}:
\begin{gather}
\dot{\mathbf{U}} - D \nabla^{2} \mathbf{U}=-P^{-1}\dot{\mathbf{W}}
+DP^{-1}\nabla^{2} \mathbf{W}-P^{-1}M^{-1}\mathbf{b}_0
(\nabla \cdot \dot{\mathbf{u}}), \quad \text{in } \Omega \times (0, t_f) \label{eq30}
 \\
\mathbf{U}(\mathbf{x},t)= \mathbf{0},  \quad \text{on }
\Gamma \times (0,t_f) \label{eq31} \\
\mathbf{U}(\mathbf{x},0)=\mathbf{U}_0(\mathbf{x}), \quad \text{in }
\bar{\Omega} \label{eq32}
\end{gather}
where
\begin{equation}
\mathbf{U}_0(\mathbf{x})=
\begin{cases}
P^{-1}\big(\mathbf{V}_I(\mathbf{x})-\mathbf{W}(\mathbf{x},0)\big), & \text{in } \Omega\\
\mathbf{0}, & \text{on } \Gamma
\end{cases} \label{eq33}
\end{equation}

\begin{remark} \label{rmk3.3} \rm
From \eqref{eq20}, $\mathbf{u} \in \mathcal{U}_{ad}$ satisfies
$\nabla \cdot \dot{\mathbf{u}} \in L^2(0,t_f; L^2(\Omega))$.
Also, Assumption \ref{assm2} and \eqref{eq28} yield
$\dot{\mathbf{W}} \in L^2\big(0,t_f; \mathbb{V}^3\big)$ and, with the use
of \eqref{eq33}, it follows that $\mathbf{U}_0 \in \mathbb{V}_0^3$.
\end{remark}

Define a bilinear form $a: \mathbb{V}_0^3 \times \mathbb{V}_0^3 \to \mathbb{R}$ by
\begin{equation*}
a(\mathbf{v},\mathbf{w})=\sum_{k=1}^3 \lambda_k \int_{\Omega}
\nabla v_k \cdot\nabla w_k \mathrm{d}\Omega,
\end{equation*}
where $\mathbf{v}=[v_1 ~v_2 ~v_3]^T$ and ${\mathbf{w}}=[w_1 ~w_2 ~w_3]^T$.
The form $a(\cdot, \cdot)$ defines a scalar product on $\mathbb{V}_0^3$
with the norm equivalent to the standard norm on $\mathbb{V}_0^3$ and the
following inequalities hold:
\begin{gather}
|a(\mathbf{v}, \mathbf{w})| \leq 3\lambda_{\rm max}\|\mathbf{v}\|_{\mathbb{V}_0^3}\|\mathbf{w}\|_{\mathbb{V}_0^3}, \quad \mathbf{v}, \mathbf{w} \in \mathbb{V}_0^3 \label{eq34} \\
\lambda_{\rm min}\|\mathbf{v}\|_{\mathbb{V}_0^3}^2
\leq a(\mathbf{v},\mathbf{v})\leq \lambda_{\rm max}
\|\mathbf{v}\|_{\mathbb{V}_0^3}^2, \quad \mathbf{v} \in \mathbb{V}_0^3 \label{eq35}
\end{gather}
We denote this scalar product by $(\cdot, \cdot)_a$; that is,
\begin{equation}
a(\mathbf{v},\mathbf{w})=(\mathbf{v},\mathbf{w})_a \label{eq36}
\end{equation}
and the corresponding norm is $\|\cdot\|_a=(\cdot, \cdot)_a^{1/2}$.

We consider the following weak formulation of the problem \eqref{eq30}-\eqref{eq33}.

Given $\mathbf{W} \in L^2\big(0,t_f; \mathbb{V}^3\big)$,
$\dot{\mathbf{W}} \in L^2\big(0,t_f; \mathbb{V}^3\big)$,
$\nabla \cdot \dot{\mathbf{u}} \in L^2(0,t_f; L^2(\Omega))$, and
$\mathbf{U}_0\in \mathbb{H}^3$, find $\mathbf{U}\in L^2(0,t_f;\mathbb{V}_0^3)$
such that, for all $\boldsymbol{\varphi}\in\mathbb{V}_0^3$,
\begin{equation}
\begin{aligned}
\frac{d}{dt}(\mathbf{U},\boldsymbol{\varphi})_{\mathbb{H}^3}
+(\mathbf{U},\boldsymbol{\varphi})_a
&= (-P^{-1}\dot{\mathbf{W}},\boldsymbol{\varphi})_{\mathbb{H}^3}
+ (DP^{-1}\nabla^{2} \mathbf{W}, \boldsymbol{\varphi})_{\mathbb{H}^3}  \\
&\quad+ (-P^{-1}M^{-1}\mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}),
\boldsymbol{\varphi})_{\mathbb{H}^3}
\end{aligned}  \label{eq37}
\end{equation}
and
\begin{equation}
\mathbf{U}(0)=\mathbf{U}_0 \label{eq38}
\end{equation}
The next lemma establishes the Hadamard well-posedness of the problem \eqref{eq37}
 and \eqref{eq38}.

\begin{lemma} \label{lem1}
Problem \eqref{eq37} and \eqref{eq38} admits a unique solution
$\mathbf{U}\in \mathcal{W}(0,t_f;\mathbb{V}_0^3)$ and this solution depends
continuously on the data, that is the mapping
\begin{equation*}
\mathbf{W}, \dot{\mathbf{W}}, \nabla \cdot \dot{\mathbf{u}},
\mathbf{U}_0 \mapsto \mathbf{U}
\end{equation*}
from $L^2\big(0,t_f; \mathbb{V}^3\big) \times L^2\big(0,t_f; \mathbb{V}^3\big)
\times L^2(0,t_f; L^2(\Omega)) \times \mathbb{H}^3$ to
$\mathcal{W}(0,t_f;\mathbb{V}_0^3)$ is continuous.
\end{lemma}

\begin{proof}
The result is based on the standard Faedo-Galerkin approximation
technique and we omit the proof of the existence and uniqueness of a solution,
 as it essentially follows the proofs of \cite[Theorem 1.2, Chapter III]{Lions1}
and \cite[Theorem 1.1, Chapter III]{Temam}.

We will show that the solution $\mathbf{U}\in \mathcal{W}(0,t_f;\mathbb{V}_0^3)$
depends continuously on the data
$\mathbf{W}, \dot{\mathbf{W}}, \nabla \cdot \dot{\mathbf{u}}$ and
$\mathbf{U}_0$. For each $m \in \mathbb{Z}_+$, define an approximate solution
$\mathbf{U}_m$ as
\begin{equation*}
\mathbf{U}_m(\mathbf{x},t)=\sum_{i=1}^m g_{im}(t)\mathbf{w}_i(\mathbf{x})
\end{equation*}
where $g_{im}(t), 1\leq i\leq m$, are scalar functions defined on
$[0,t_f]$, $ \mathbf{w}_1, \ldots, \mathbf{w}_m, \ldots$ is a countable
set of functions which is dense in $\mathbb{V}_0^3$,
$\mathbf{U}_m(\mathbf{x},t)$ satisfies
\begin{equation}
\begin{aligned}
&(\dot{\mathbf{U}}_m(t),\mathbf{w}_i)_{\mathbb{H}^3}
 +(\mathbf{U}_m(t),\mathbf{w}_i)_a \\
& = (-P^{-1}\dot{\mathbf{W}}(t),\mathbf{w}_i)_{\mathbb{H}^3} + \big(\nabla(-DP^{-1}\mathbf{W}(t)), \nabla \mathbf{w}_i\big)_{\mathbb{H}^3}  \\
&\quad + \big(-P^{-1}M^{-1}\mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}(t)),
\mathbf{w}_i\big)_{\mathbb{H}^3}, \quad 1\leq i \leq m
\end{aligned} \label{eq39}
\end{equation}
and
\begin{equation}
\mathbf{U}_m(\mathbf{x},0)= \mathbf{U}_{0m}\underset{m\to \infty}
{\longrightarrow}\mathbf{U}_0 \quad \text{in }\mathbb{H}^3 \label{eq40}
\end{equation}
Multiplying \eqref{eq39} by $g_{jm}(t)$, adding for $j$, and using
\begin{equation*}
(\dot{\mathbf{U}}_m(t),\mathbf{U}_m(t))_{\mathbb{H}^3}
=\frac{1}{2}\|\dot{\mathbf{U}}_m(t)\|_{\mathbb{H}^3}^2
\end{equation*}
we have
\begin{align*}
&\frac{d}{dt}\|\mathbf{U}_m(t)\|_{\mathbb{H}^3}^2+2\|\mathbf{U}_m(t)\|_a^2 \\
& \leq 2|(P^{-1}\dot{\mathbf{W}}(t),\mathbf{U}_m(t))_{\mathbb{H}^3}|
 + 2|\big(\nabla(DP^{-1}\mathbf{W}(t)), \nabla \mathbf{U}_m(t)\big)_{\mathbb{H}^3}| \\
& \quad+ 2|\big(P^{-1}M^{-1}\mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}(t)),
\mathbf{U}_m(t)\big)_{\mathbb{H}^3}|
\end{align*}
 Applying \eqref{eq35}, \eqref{eq36}, and the Poincar\'e inequality gives
\begin{align*}
&\frac{d}{dt}\|\mathbf{U}_m(t)\|_{\mathbb{H}^3}^2
 +2\lambda_{\rm min}\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}^2\\
& \leq \frac{3c}{\lambda_{\rm min}}\|P^{-1}\|^2
 \|\dot{\mathbf{W}}\|_{\mathbb{V}^3}^2
 +\frac{\lambda_{\rm min}}{3c}c\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}^2 \\
& \quad +\frac{3}{\lambda_{\rm min}}\|DP^{-1}\|^2 \|\mathbf{W}\|_{\mathbb{V}^3}^2
 +\frac{\lambda_{\rm min}}{3}\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}^2  \\
& \quad +\frac{3c}{\lambda_{\rm min}}|P^{-1}M^{-1}\mathbf{b}_0|^2\|\nabla
 \cdot \dot{\mathbf{u}}(t)\|_{L^2(\Omega)}^2
 + \frac{\lambda_{\rm min}}{3c}c\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}^2
\end{align*}
and therefore
\begin{align*}
\frac{d}{dt}\|\mathbf{U}_m(t)\|_{\mathbb{H}^3}^2
 +\lambda_{\rm min}\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}^2
&\leq \frac{3c}{\lambda_{\rm min}}\|P^{-1}\|^2\|
 \dot{\mathbf{W}}\|_{\mathbb{V}^3}^2
 + \frac{3\lambda_{\rm max}^2}{\lambda_{\rm min}}\|P^{-1}\|^2
  \|\mathbf{W}\|_{\mathbb{V}^3}^2 \\
& \quad +\frac{3c}{\lambda_{\rm min}}|P^{-1}M^{-1}\mathbf{b}_0|^2
 \|\nabla \cdot \dot{\mathbf{u}}(t)\|_{L^2(\Omega)}^2
\end{align*}
Integrating the last inequality over $(0,t_f)$ and applying \eqref{eq40},
we obtain
\begin{equation}
\begin{aligned}
\|\mathbf{U}_m\|_{L^2(0,t_f;\mathbb{V}_0^3)}^2
&\leq \frac{1}{\lambda_{\rm min}}\|\mathbf{U}_0\|_{\mathbb{H}^3}^2
 + \frac{3c}{\lambda_{\rm min}^2}\|P^{-1}\|^2\|
 \dot{\mathbf{W}}\|_{L^2(0,t_f;\mathbb{V}^3)}^2  \\
&\quad +\frac{3\lambda_{\rm max}^2}{\lambda_{\rm min}^2}\|P^{-1}\|^2
 \|\mathbf{W}\|_{L^2(0,t_f;\mathbb{V}^3)}^2  \\
&\quad + \frac{3c}{\lambda_{\rm min}^2}|P^{-1}M^{-1}\mathbf{b}_0|^2
\|\nabla \cdot \dot{\mathbf{u}}\|_{L^2(0,t_f;L^2(\Omega))}^2
\end{aligned} \label{eq41}
\end{equation}
The inequality \eqref{eq41} holds for each $m\in \mathbb{Z}_+$ and therefore,
\begin{equation}
\begin{aligned}
\|\mathbf{U}\|_{L^2(0,t_f;\mathbb{V}_0^3)}^2
 &\leq \frac{1}{\lambda_{\rm min}}\|\mathbf{U}_0\|_{\mathbb{H}^3}^2
 + \frac{3c}{\lambda_{\rm min}^2}\|P^{-1}\|^2
 \|\dot{\mathbf{W}}\|_{L^2(0,t_f;\mathbb{V}^3)}^2  \\
&\quad +\frac{3\lambda_{\rm max}^2}{\lambda_{\rm min}^2}\|P^{-1}\|^2
 \|\mathbf{W}\|_{L^2(0,t_f;\mathbb{V}^3)}^2  \\
&\quad + \frac{3c}{\lambda_{\rm min}^2}|P^{-1}M^{-1}\mathbf{b}_0|^2
 \|\nabla \cdot \dot{\mathbf{u}}\|_{L^2(0,t_f;L^2(\Omega))}^2 \label{eq42}
\end{aligned}
\end{equation}
To estimate $\dot{\mathbf{U}}_m(t)$, note that
$\dot{\mathbf{U}}_m(t) = {\sum_{i=1}^m} \dot{g}_{im}(t)\mathbf{w}_i(\mathbf{x})
\in \text{span}\{\mathbf{w}_1, \ldots,\mathbf{w}_m\}=E_m \subset \mathbb{V}_0^3$
and
\begin{equation}
\|\dot{\mathbf{U}}_m(t)\|_{H^{-1}(\Omega)^3}
=\sup_{\mathbf{v}\in E_m\setminus \{0\}} \frac{|(\dot{\mathbf{U}}_m(t),
\mathbf{v})_{\mathbb{H}^3}|}{\|\mathbf{v}\|_{\mathbb{V}_0^3}} \label{eq43}
\end{equation}
Using \eqref{eq34}, \eqref{eq36}, \eqref{eq39}, the Poincar\'e inequality,
and linearity argument yields
\begin{equation}
\begin{aligned}
(\dot{\mathbf{U}}_m(t), \boldsymbol{\varphi})_{\mathbb{H}^3}
&\leq |a(\mathbf{U}_m(t),\boldsymbol{\varphi})|
+ |(P^{-1}\dot{\mathbf{W}},\boldsymbol{\varphi})_{\mathbb{H}^3}|
 + |(\nabla(DP^{-1}\mathbf{W}), \nabla\boldsymbol{\varphi})_{\mathbb{H}^3}|  \\
&\quad+ |(P^{-1}M^{-1}\mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}),
\boldsymbol{\varphi})_{\mathbb{H}^3}|  \\
& \leq \Big(3\lambda_{\rm max}\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}
 +c\|P^{-1}\| \|\dot{\mathbf{W}}\|_{\mathbb{V}^3}+\lambda_{\rm max}\|P^{-1}\|
\|\mathbf{W}\|_{\mathbb{V}^3}  \\
&\quad + c |P^{-1}M^{-1}\mathbf{b}_0| \|\nabla \cdot
\dot{\mathbf{u}}\|_{L^2(\Omega)} \Big)\|\boldsymbol{\varphi}\|_{\mathbb{V}_0^3},
\quad \forall \boldsymbol{\varphi}\in \mathbb{V}_0^3
\end{aligned}\label{eq44}
\end{equation}
From \eqref{eq43} and \eqref{eq44}, we have
\begin{align*}
\|\dot{\mathbf{U}}_m(t)\|_{H^{-1}(\Omega)^3}
&\leq 3\lambda_{\rm max}\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}
+c\|P^{-1}\| \|\dot{\mathbf{W}}\|_{\mathbb{V}^3}
 +\lambda_{\rm max}\|P^{-1}\| \|\mathbf{W}\|_{\mathbb{V}^3} \\
&\quad + c |P^{-1}M^{-1}\mathbf{b}_0| \|\nabla \cdot
\dot{\mathbf{u}}\|_{L^2(\Omega)}
\end{align*}
Squaring the last inequality and integrating it over $(0,t_f)$ gives
\begin{align*}
\|\dot{\mathbf{U}}_m\|_{L^2(0,t_f;H^{-1}(\Omega)^3)}^2
&\leq 36\lambda_{\rm max}^2\|\mathbf{U}_m\|_{L^2(0,t_f;
 \mathbb{V}_0^3)}^2\\
&\quad +4c^2\|P^{-1}\|^2 \|\dot{\mathbf{W}}
 \|_{L^2(0,t_f;\mathbb{V}^3)}^2
 + 4\lambda_{\rm max}^2\|P^{-1}\|^2 \|\mathbf{W}\|_{L^2(0,t_f;\mathbb{V}^3)}^2 \\
&\quad + 4c^2 |P^{-1}M^{-1}\mathbf{b}_0|^2 \|\nabla \cdot
\dot{\mathbf{u}}\|_{L^2(0,t_f;L^2(\Omega))}^2
\end{align*}
Applying \eqref{eq41} to the first term on the right-hand side of this
inequality, we obtain
\begin{align*}
\|\dot{\mathbf{U}}_m\|_{L^2(0,t_f;H^{-1}(\Omega)^3)}^2
&\leq \frac{36\lambda_{\rm max}^2}{\lambda_{\rm min}}
\|\mathbf{U}_0\|_{\mathbb{H}^3}^2\\
&\quad+\Big(\frac{108c\lambda_{\rm max}^2}{\lambda_{\rm min}^2}
+4c^2\Big)\|P^{-1}\|^2 \|\dot{\mathbf{W}}\|_{L^2(0,t_f;\mathbb{V}^3)}^2 \\
&\quad + \Big(\frac{108\lambda_{\rm max}^4}{\lambda_{\rm min}^2}
+ 4\lambda_{\rm max}^2\Big)\|P^{-1}\|^2 \|\mathbf{W}\|_{L^2(0,t_f;\mathbb{V}^3)}^2 \\
&\quad + \Big(\frac{108c\lambda_{\rm max}^2}{\lambda_{\rm min}^2}
+ 4c^2\Big) |P^{-1}M^{-1}\mathbf{b}_0|^2 \|\nabla \cdot
\dot{\mathbf{u}}\|_{L^2(0,t_f;L^2(\Omega))}^2
\end{align*}
The last inequality holds for each $m\in \mathbb{Z}_+$ and therefore,
\begin{equation}
\begin{aligned}
&\|\dot{\mathbf{U}}\|_{L^2(0,t_f;H^{-1}(\Omega)^3)}^2 \\
&\leq \frac{36\lambda_{\rm max}^2}{\lambda_{\rm min}}
 \|\mathbf{U}_0\|_{\mathbb{H}^3}^2
 +\Big(\frac{108c\lambda_{\rm max}^2}{\lambda_{\rm min}^2}
 +4c^2\Big)\|P^{-1}\|^2 \|\dot{\mathbf{W}}\|_{L^2(0,t_f;\mathbb{V}^3)}^2  \\
&\quad + \Big(\frac{108\lambda_{\rm max}^4}{\lambda_{\rm min}^2}
 + 4\lambda_{\rm max}^2\Big)\|P^{-1}\|^2 \|\mathbf{W}\|_{L^2(0,t_f;\mathbb{V}^3)}^2  \\
&\quad + \Big(\frac{108c\lambda_{\rm max}^2}{\lambda_{\rm min}^2}
 + 4c^2\Big) |P^{-1}M^{-1}\mathbf{b}_0|^2 \|\nabla \cdot
 \dot{\mathbf{u}}\|_{L^2(0,t_f;L^2(\Omega))}^2
\end{aligned} \label{eq45}
\end{equation}
The result follows from \eqref{eq42} and \eqref{eq45}.
\end{proof}

Lemma \ref{lem1} yields the Hadamard well-posedness in the weak sense of the
diffusion system \eqref{eq1}, \eqref{eq3}, and \eqref{eq7}, namely,
the existence and uniqueness of a weak solution to the model and its
 continuous dependence on the boundary data, initial data, and the
divergence of the rock deformation velocity field.

\begin{theorem}[Well-posedness of the diffusion system] \label{theorem1}
	Given boundary data \\
 $\mathbf{V}_B \in L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big)$ with
$\dot{\mathbf{V}}_B \in L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big)$,
initial data $\mathbf{V}_I \in \mathbb{H}^3$, and
$\nabla \cdot \dot{\mathbf{u}} \in L^2\big(0,t_f;L^2(\Omega)\big)$,
the diffusion system \eqref{eq1}, \eqref{eq3}, and \eqref{eq7} admits a
unique weak solution $\mathbf{V}\in L^2\big(0,t_f;\mathbb{V}^3\big)$ with
$\dot{\mathbf{V}}\in L^2\big(0,t_f;H^{-1}(\Omega)^3\big)$ and this solution
depends continuously on the data
$\mathbf{V}_B, \dot{\mathbf{V}}_B, \mathbf{V}_B(0), \mathbf{V}_I$, and
$\nabla \cdot \dot{\mathbf{u}}$. That is, the mapping
\begin{equation*}
\mathbf{V}_B, \dot{\mathbf{V}}_B, \mathbf{V}_B(0), \mathbf{V}_I,
\nabla \cdot \dot{\mathbf{u}} \mapsto \mathbf{V}, \dot{\mathbf{V}}
\end{equation*}
from $L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big) \times
L^2\big(0,t_f; H^{1/2}(\Gamma)^3\big) \times H^{1/2}(\Gamma)^3
\times \mathbb{H}^3 \times L^2\big(0,t_f;L^2(\Omega)\big)$ to
$L^2\big(0,t_f;\mathbb{V}^3\big) \times L^2\big(0,t_f;H^{-1}(\Omega)^3\big)$
is continuous.
\end{theorem}

\begin{proof}
The existence and uniqueness of a weak solution follow immediately from the
transformation \eqref{eq29}. It remains to prove the continuous dependence
on data. From \eqref{eq28}, \eqref{eq29} and the Poincar\'e inequality, we have
\begin{gather}
\|\mathbf{V}\|_{L^2(0,t_f;\mathbb{V}^3)}^2
\leq 2(1+c)\|P\|^2\|\mathbf{U}\|_{L^2(0,t_f;\mathbb{V}_0^3)}^2
+ 2\hat{C}\|\mathbf{V}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}^2 \label{eq46}
\\
\|\dot{\mathbf{V}}\|_{L^2(0,t_f;H^{-1}(\Omega)^3)}^2
\leq 2\|P\|^2\|\dot{\mathbf{U}}\|_{L^2(0,t_f;H^{-1}(\Omega)^3)}^2
+ 2\hat{C}\|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}^2 \label{eq47}
\end{gather}
where $\hat{C}>0$ is a constant. Also, the transformation \eqref{eq33} gives
\begin{equation}
\|\mathbf{U}_0\|_{\mathbb{H}^3}^2 \leq 2\|P^{-1}\|^2\|\mathbf{V}_I\|_{\mathbb{H}^3}^2
+ 2\|P^{-1}\|^2\hat{C}\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3}^2 \label{eq48}
\end{equation}
Applying \eqref{eq48} to \eqref{eq42} and \eqref{eq45} and substituting the
results together with
\begin{gather*}
\|\mathbf{W}\|_{L^2(0,t_f;\mathbb{V}^3)}^2
\leq \hat{C}\|\mathbf{V}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}^2 \\
\|\dot{\mathbf{W}}\|_{L^2(0,t_f;\mathbb{V}^3)}^2
\leq \hat{C}\|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}^2
\end{gather*}
into \eqref{eq46} and \eqref{eq47}, respectively, complete the proof.
\end{proof}

The next lemma further characterizes the solution $\mathbf{V}$ to the diffusion
system \eqref{eq1}, \eqref{eq3}, and \eqref{eq7} dealing with the initial
data $\mathbf{V}_I \in \mathbb{V}^3$. This result will be needed in the
later discussion of the coupled TCP model.

\begin{lemma} \label{lemma2}
Under assumption $\mathbf{V}_I \in \mathbb{V}^3$, the weak solution
$\mathbf{V}$ of the diffusion system \eqref{eq1}, \eqref{eq3}, and
\eqref{eq7} satisfies  $\dot{\mathbf{V}}\in L^2(0,t_f;\mathbb{H}^3)$ and
the following a priori estimate holds
\begin{equation}
\|\dot{\mathbf{V}}\|_{L^2(0,t_f;\mathbb{H}^3)}^2
\leq \Big(4+\frac{6}{\lambda_{\rm min}^2}\Big)\|P\|^2 |P^{-1}M^{-1}
\mathbf{b}_0|^2 \|\nabla \cdot \dot{\mathbf{u}}\|_{L^2(0,t_f;L^2(\Omega))}^2
+ \theta \label{eq49}
\end{equation}
where $\theta>0$ is a constant that does not depend on
$\nabla \cdot \dot{\mathbf{u}}$.
\end{lemma}

\begin{proof}
As before, we use the standard Faedo-Galerkin approximation method applied
to the problem \eqref{eq37} and \eqref{eq38}. Note that the assumption
 $\mathbf{V}_I \in \mathbb{V}^3$ and \eqref{eq33} imply
$\mathbf{U}_0 \in \mathbb{V}_0^3$.  Multiplying \eqref{eq39} by
$\dot{g}_{jm}(t)$, adding for $j$, and using
$(\mathbf{U}_m(t),\dot{\mathbf{U}}_m(t))_a
=\frac{1}{2}\|\dot{\mathbf{U}}_m(t)\|_a$
yield the  inequality
\begin{align*}
&2\|\dot{\mathbf{U}}_m(t)\|_{\mathbb{H}^3}^2 
 + \frac{d}{dt}\|\mathbf{U}_m(t)\|_a^2\\
& \leq 2|(P^{-1}\dot{\mathbf{W}},\dot{\mathbf{U}}_m(t))_{\mathbb{H}^3}| 
 +2\frac{d}{dt}\big(\nabla(-DP^{-1}\mathbf{W}),
 \nabla\mathbf{U}_m(t)\big)_{\mathbb{H}^3} \\
&\quad  +2|\big(\nabla(DP^{-1}\dot{\mathbf{W}}), 
 \nabla\mathbf{U}_m(t)\big)_{\mathbb{H}^3}| 
 + 2|\big(P^{-1}M^{-1}\mathbf{b}_0 (\nabla \cdot
\dot{\mathbf{u}}), \dot{\mathbf{U}}_m(t)\big)_{\mathbb{H}^3}|
\end{align*}
which, in turn, gives
\begin{align*}
&2\|\dot{\mathbf{U}}_m(t)\|_{\mathbb{H}^3}^2
+ \frac{d}{dt}\|\mathbf{U}_m(t)\|_a^2 \\
&\leq 2\|P^{-1}\|^2 \|\dot{\mathbf{W}}\|_{\mathbb{V}^3}^2
 + \frac{1}{2}\|\dot{\mathbf{U}}_m(t)\|_{\mathbb{H}^3}^2 \\
& + 2\frac{d}{dt}\big(\nabla(-DP^{-1}\mathbf{W}),
 \nabla\mathbf{U}_m(t)\big)_{\mathbb{H}^3} +c\|DP^{-1}\|^2
  \|\dot{\mathbf{W}}\|_{\mathbb{V}^3}^2  \\
&+\frac{1}{c}\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}^2
  + 2|P^{-1}M^{-1}\mathbf{b}_0|^2 \|\nabla \cdot \dot{\mathbf{u}}\|_{L^2(\Omega)}^2
+ \frac{1}{2}\|\dot{\mathbf{U}}_m(t)\|_{\mathbb{H}^3}^2
\end{align*}
and then
\begin{align*}
\|\dot{\mathbf{U}}_m(t)\|_{\mathbb{H}^3}^2
+ \frac{d}{dt}\|\mathbf{U}_m(t)\|_a^2
 &\leq (2+c\lambda_{\rm max}^2)\|P^{-1}\|^2 \|\dot{\mathbf{W}}\|_{\mathbb{V}^3}^2
 + \frac{1}{c}\|\mathbf{U}_m(t)\|_{\mathbb{V}_0^3}^2  \\
& \quad + 2\frac{d}{dt}\big(\nabla(-DP^{-1}\mathbf{W}),
 \nabla\mathbf{U}_m(t)\big)_{\mathbb{H}^3} \\
&\quad + 2|P^{-1}M^{-1}\mathbf{b}_0|^2 \|\nabla \cdot
\dot{\mathbf{u}}\|_{L^2(\Omega)}^2
\end{align*}
Integrating the above inequality over $(0,t_f)$ and using \eqref{eq35}
 and \eqref{eq36}, we have
\begin{align*}
&\|\dot{\mathbf{U}}_m\|_{L^2(0,t_f;\mathbb{H}^3)}^2
 + \lambda_{\rm min}\|\mathbf{U}_m(t_f)\|_{\mathbb{V}_0^3}^2 \\
&\leq \lambda_{\rm max}\|\mathbf{U}_m(0)\|_{\mathbb{V}_0^3}^2
 + (2+c\lambda_{\rm max}^2)\|P^{-1}\|^2 \|\dot{\mathbf{W}}\|_{L^2(0,t_f
 ;\mathbb{V}^3)}^2 \\
& \quad +\lambda_{\rm max}\|P^{-1}\|^2\|\mathbf{W}(0)\|_{\mathbb{V}^3}^2
 + \lambda_{\rm max} \|\mathbf{U}_m(0)\|_{\mathbb{V}_0^3}^2\\
& \quad +\frac{\lambda_{\rm max}^2}{\lambda_{\rm min}}
 \|P^{-1}\|^2\|\mathbf{W}(t_f)\|_{\mathbb{V}^3}^2
 +\lambda_{\rm min}\|\mathbf{U}_m(t_f)\|_{\mathbb{V}_0^3}^2\\
& \quad +\frac{1}{c}\|\mathbf{U}_m\|_{L^2(0,t_f;\mathbb{V}_0^3)}^2
 + 2|P^{-1}M^{-1}\mathbf{b}_0|^2 \|\nabla \cdot
 \dot{\mathbf{u}}\|_{L^2(0,t_f;L^2(\Omega))}^2
\end{align*}
Applying \eqref{eq41} and the Poincar\'e inequality to the last inequality, we obtain
\begin{align*}
\|\dot{\mathbf{U}}_m\|_{L^2(0,t_f;\mathbb{H}^3)}^2
&\leq \Big(\frac{1}{\lambda_{\rm min}}+2\lambda_{\rm max}\Big)
 \|\mathbf{U}_0\|_{\mathbb{V}_0^3}^2
 + \lambda_{\rm max}\|P^{-1}\|^2\|\mathbf{W}(0)\|_{\mathbb{V}^3}^2 \\
& \quad + \frac{\lambda_{\rm max}^2}{\lambda_{\rm min}}
 \|P^{-1}\|^2\|\mathbf{W}(t_f)\|_{\mathbb{V}^3}^2
 + \frac{3\lambda_{\rm max}^2}{c\lambda_{\rm min}^2} \|P^{-1}\|^2
  \|\mathbf{W}\|_{L^2(0,t_f;\mathbb{V}^3)}^2 \\
& \quad + \Big(2+c\lambda_{\rm max}^2
 +\frac{3}{\lambda_{\rm min}^2}\Big)\|P^{-1}\|^2
 \|\dot{\mathbf{W}}\|_{L^2(0,t_f;\mathbb{V}^3)}^2 \\
& \quad + \Big(2+\frac{3}{\lambda_{\rm min}^2}\Big)
  |P^{-1}M^{-1}\mathbf{b}_0|^2 \|\nabla \cdot \dot{\mathbf{u}}
 \|_{L^2(0,t_f;L^2(\Omega))}^2, \quad \forall m\in\mathbb{Z}_+
\end{align*}
This shows that the sequence $\dot{\mathbf{U}}_m$ ranges in a bounded set
of $L^2(0,t_f;\mathbb{H}^3)$ and
\begin{equation}
\|\dot{\mathbf{U}}\|_{L^2(0,t_f;\mathbb{H}^3)}^2
\leq \Big(2+\frac{3}{\lambda_{\rm min}^2}\Big)
|P^{-1}M^{-1}\mathbf{b}_0|^2 \|\nabla \cdot
\dot{\mathbf{u}}\|_{L^2(0,t_f;L^2(\Omega))}^2 + \hat{\theta}\label{eq50}
\end{equation}
where $\hat{\theta}>0$ is a constant independent of
 $\nabla \cdot \dot{\mathbf{u}}$. From \eqref{eq28} and \eqref{eq29},
\begin{equation*}
\|\dot{\mathbf{V}}\|_{L^2(0,t_f;\mathbb{H}^3)}^2
\leq 2\|P\|^2\|\dot{\mathbf{U}}\|_{L^2(0,t_f;\mathbb{H}^3)}^2
+ 2\hat{C}\|\dot{\mathbf{V}}_B\|_{L^2(0,t_f;H^{1/2}(\Gamma)^3)}
\end{equation*}
where $\hat{C}>0$ is a constant, and applying \eqref{eq50} the result follows.
\end{proof}


	\subsection{Elastic system}
From now on, we will consider the Navier-type elastic system
\eqref{eq8}, \eqref{eq9}, and \eqref{eq4}. The following assumption
 will be made on the applied boundary stress tensor $\boldsymbol{\hat{\sigma}}$.

\begin{assumption} \label{assm3} \rm
\begin{equation*}
\hat{\boldsymbol{\sigma}} \in L^2\big(0,t_f;L^2(\Gamma)^{n\times n}\big) \quad
\text{and} \quad \dot{\hat{\boldsymbol{\sigma}}} \in L^2\big(0,t_f;L^2(\Gamma)^{n\times n}\big)
\end{equation*}
where $n=2$ or $3$ is the dimension of the problem.
\end{assumption}

In solid mechanics the weak or variational formulation of a boundary value
problem is equivalent to the principle of minimum total potential energy.
 Utilizing the notation introduced in Section \ref{model}, let
 $a_E:\mathbb{V}^n \times \mathbb{V}^n \to \mathbb{R}$, $n=2$ or $3$,
 be a bilinear form defined as
\begin{equation}
a_E(\mathbf{u},\boldsymbol{\Phi})
= \int_\Omega \sum_{i,j=1}^{n} \tau_{ij}(\mathbf{u})
\varepsilon_{ij}(\boldsymbol{\Phi})\mathrm{d} \Omega, \quad
\mathbf{u}, \boldsymbol{\Phi} \in \mathbb{V}^n \label{eq51}
\end{equation}
where $\boldsymbol{\tau} = [\tau_{ij}]$ and
$\boldsymbol{\varepsilon} = [\varepsilon_{ij}]$, $i,j =1, \ldots, n$,
are the stress and strain tensors, respectively.
From \eqref{eq18}, \eqref{eq20}, and \eqref{eq21}, under
Assumption \ref{assm3}, we obtain the following weak formulation of
problem \eqref{eq8}, \eqref{eq9}, and \eqref{eq4}:

Given $\mathbf{V}\in L^2\big(0,t_f;\mathbb{V}^3\big)$ with
$\dot{\mathbf{V}}\in L^2\big(0,t_f;\mathbb{H}^3\big)$ and
$\hat{\boldsymbol{\sigma}}\in L^2\big(0,t_f;L^2(\Gamma)^{n\times n}\big)$
with $\dot{\hat{\boldsymbol{\sigma}}}\in L^2\big(0,t_f;L^2(\Gamma)^{n\times n}\big)$,
 find $\mathbf{u}\in L^2(0,t_f;\tilde{\mathbb{V}}_0^n)$ such that
$ \dot{\mathbf{u}}\in L^2(0,t_f;\tilde{\mathbb{V}}_0^n)$ and,
for each $t \in [0,t_f)$,
\begin{equation}
a_E(\mathbf{u},\boldsymbol{\Phi})-\int_\Omega (\mathbf{b}_1 \cdot
\mathbf{V})\nabla \cdot \boldsymbol{\Phi}\mathrm{d} \Omega
 - \int_{\Gamma \setminus \Gamma_F}(\boldsymbol{\hat{\sigma}}
\mathbf{n})\cdot\boldsymbol{\Phi}\mathrm{d} \Gamma = 0, \quad
\forall \boldsymbol{\Phi} \in \tilde{\mathbb{V}}_0^n \label{eq52}
\end{equation}
where $\mathbf{n}$ is the outward unit normal vector on the boundary.

Next we will show that problem \eqref{eq52} is well-posed in the sense
 of Hadamard. Define continuous linear functionals $F$ and $\dot{F}$ on
$\mathbb{V}^n$ by the pairings
\begin{equation*}
\boldsymbol{\Phi}\mapsto \langle F, \boldsymbol{\Phi}\rangle \quad \text{and} \quad
\boldsymbol{\Phi}\mapsto \langle \dot{F}, \boldsymbol{\Phi}\rangle
\end{equation*}
where
\begin{gather}
\langle F,\boldsymbol{\Phi}\rangle
= \int_\Omega (\mathbf{b}_1 \cdot \mathbf{V})\nabla \cdot \boldsymbol{\Phi}
\mathrm{d} \Omega
+ \int_{\Gamma \setminus \Gamma_F}(\hat{\boldsymbol{\sigma}}\mathbf{n})
\cdot\boldsymbol{\Phi}\mathrm{d} \Gamma,  \label{eq53}
\\
\langle \dot{F},\boldsymbol{\Phi}\rangle
= \int_\Omega (\mathbf{b}_1 \cdot \dot{\mathbf{V}})\nabla \cdot
\boldsymbol{\Phi}\mathrm{d} \Omega
+ \int_{\Gamma \setminus \Gamma_F}(\dot{\hat{\boldsymbol{\sigma}}}
\mathbf{n})\cdot\boldsymbol{\Phi}\mathrm{d} \Gamma, \quad
\forall \boldsymbol{\Phi} \in \mathbb{V}^n \label{eq54}
\end{gather}
From \eqref{eq52}-\eqref{eq54}, it follows that, for each $t \in [0,t_f)$,
\[
a_E(\mathbf{u},\boldsymbol{\Phi})=\langle F,\boldsymbol{\Phi}\rangle ,  \quad
a_E(\dot{\mathbf{u}},\boldsymbol{\Phi})=\langle \dot{F},\boldsymbol{\Phi}\rangle,
\quad \forall \boldsymbol{\Phi} \in \tilde{\mathbb{V}}_0^n
\]
and, using the trace theorem, for each
$\boldsymbol{\Phi} \in \tilde{\mathbb{V}}_0^n$ and $ t \in [0,t_f)$,
\begin{gather}
\begin{aligned}
|a_E(\mathbf{u},\boldsymbol{\Phi})|
& \leq \|\mathbf{b}_1 \cdot \mathbf{V}\|_{L^2(\Omega)}
 \|\nabla \cdot \boldsymbol{\Phi}\|_{L^2(\Omega)}
 + \|\boldsymbol{\hat{\sigma}}\mathbf{n}\|_{L^2(\Gamma)^n}
 \|\boldsymbol{\Phi}\|_{L^2(\Gamma)^n}  \\
& \leq \big(\sqrt{n}|\mathbf{b}_1| \|\mathbf{V}\|_{\mathbb{V}^n}
+\hat{C} \|\boldsymbol{\hat{\sigma}}\|_{L^2(\Gamma)^{n\times n}}\big)
\|\boldsymbol{\Phi}\|_{\mathbb{V}^n}
\end{aligned}\label{eq55}
\\
\begin{aligned}
|a_E(\dot{\mathbf{u}},\boldsymbol{\Phi})|
& \leq \|\mathbf{b}_1 \cdot \dot{\mathbf{V}}\|_{L^2(\Omega)}
\|\nabla \cdot \boldsymbol{\Phi}\|_{L^2(\Omega)}
+ \|\dot{\hat{\boldsymbol{\sigma}}}\mathbf{n}
 \|_{L^2(\Gamma)^n}\|\boldsymbol{\Phi}\|_{L^2(\Gamma)^n}  \\
& \leq \big(\sqrt{n}|\mathbf{b}_1| \|\dot{\mathbf{V}}\|_{\mathbb{H}^n}
+\hat{C} \|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(\Gamma)^{n\times n}}\big)
\|\boldsymbol{\Phi}\|_{\mathbb{V}^n}
\end{aligned} \label{eq56}
\end{gather}
where $\hat{C}>0$ is a constant.

To prove the existence and uniqueness of a solution to the problem \eqref{eq52},
we will  use of the Korn inequality \cite[Theorem 3.1, Chapter III]{Duvaut_Lions}
and its consequence \cite[Theorem 3.3, Chapter III]{Duvaut_Lions}
which imply that there exists a constant $c_E>0$ such that
\begin{equation}
a_E(\mathbf{v},\mathbf{v}) \geq c_E\|\mathbf{v}\|_{\mathbb{V}^n}^2, \quad
\forall \mathbf{v} \in \tilde{\mathbb{V}}_0^n, ~t \in [0,t_f) \label{eq57}
\end{equation}
Also, from \eqref{eq11}, \eqref{eq12}, and \eqref{eq51}, we have
\begin{equation}
a_E(\mathbf{u},\mathbf{v}) = a_E(\mathbf{v},\mathbf{u}), \quad
\forall \mathbf{u}, \mathbf{v} \in \mathbb{V}^n, \quad
\forall \mathbf{v} \in \tilde{\mathbb{V}}_0^n\label{eq58}
\end{equation}
and, for each $t\in [0,t_f)$,
\begin{equation}
\begin{aligned}
|a_E(\mathbf{u},\mathbf{v})|
&\leq \max_{i,j,k,l} \{a_{ijkl}\}\sum_{i,j,k,l=1}^n
\big| \int_\Omega \varepsilon_{ij}(\mathbf{u})\varepsilon_{kl}
 (\mathbf{v})\mathrm{d} \Omega\big|  \\
&\leq n\cdot \max_{i,j,k,l} \{a_{ijkl}\} \|\mathbf{u}\|_{\mathbb{V}^n}
\|\mathbf{v}\|_{\mathbb{V}^n}
\end{aligned} \label{eq59}
\end{equation}
Equations \eqref{eq57}-\eqref{eq59} yield that $a_E(\cdot, \cdot)$
is symmetric and continuous on $\tilde{\mathbb{V}}_0^n$ and there exist
constants $c_E>0$ and $\beta>0$ such that
\begin{equation*}
c_E\|\mathbf{u}\|_{\mathbb{V}^n}^2 \leq a_E(\mathbf{u},\mathbf{u})
\leq \beta \|\mathbf{u}\|_{\mathbb{V}^n}^2, \quad \forall \mathbf{u}
\in \tilde{\mathbb{V}}_0^n, \; t \in [0,t_f)
\end{equation*}
Thus, for each $t \in [0,t_f)$, $a_E(\cdot ,\cdot)$ is an inner product
 on $\tilde{\mathbb{V}}_0^n$ that is associated with a topology equivalent
to the standard topology on $\tilde{\mathbb{V}}_0^n$. From the Riesz
representation theorem and  \eqref{eq53}-\eqref{eq57} we have the following result.

\begin{lemma} \label{lem3.8}
For each $t \in [0,t_f)$, there exists a unique function
$\mathbf{u} \in \tilde{\mathbb{V}}_0^n$ with
$\dot{\mathbf{u}} \in \tilde{\mathbb{V}}_0^n$ such that \eqref{eq52} holds
for all $\boldsymbol{\Phi} \in \tilde{\mathbb{V}}_0^n$.
Furthermore,
\begin{gather*}
\|\mathbf{u}\|_{\mathbb{V}^n}
\leq \frac{\sqrt{n}}{c_E}|\mathbf{b}_1| \|\mathbf{V}\|_{\mathbb{V}^n}
+\frac{\hat{C}}{c_E}\|\hat{\boldsymbol{\sigma}}\|_{L^2(\Gamma)^{n\times n}}
\\
\|\dot{\mathbf{u}}\|_{\mathbb{V}^n}
\leq \frac{\sqrt{n}}{c_E}|\mathbf{b}_1| \|\dot{\mathbf{V}}\|_{\mathbb{H}^n}
+\frac{\hat{C}}{c_E}\|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(\Gamma)^{n\times n}}
\end{gather*}
\end{lemma}

As a consequence, we obtain the Hadamard well-posedness in the weak sense of
the Navier-type elastic system \eqref{eq8}, \eqref{eq9}, and \eqref{eq4}:

\begin{corollary}[Well-posedness of the elastic system] \label{corollary1}
Assume $\mathbf{V}\in L^2\big(0,t_f;\mathbb{V}^3\big)$ \\
with $\dot{\mathbf{V}}\in L^2\big(0,t_f;\mathbb{H}^3\big)$ and
$\hat{\boldsymbol{\sigma}}\in L^2\big(0,t_f;L^2(\Gamma)^{n\times n}\big)$
 with $\dot{\hat{\boldsymbol{\sigma}}}\in L^2\big(0,t_f;L^2(\Gamma)^{n\times n}\big)$.
Then  \eqref{eq8}, \eqref{eq9}, and \eqref{eq4} admits a unique
weak solution
 $\mathbf{u} \in L^2\big(0,t_f;\tilde{\mathbb{V}}_0^n\big)$ with
$\dot{\mathbf{u}} \in L^2\big(0,t_f;\tilde{\mathbb{V}}_0^n\big)$ and this solution
depends continuously on the data. That is, the mapping
\begin{equation*}
\mathbf{V}, \dot{\mathbf{V}}, \hat{\boldsymbol{\sigma}}, ~\dot{\hat{\boldsymbol{\sigma}}} \mapsto \mathbf{u}, \dot{\mathbf{u}}
\end{equation*}
from $L^2\big(0,t_f;\mathbb{V}^3\big) \times L^2\big(0,t_f;\mathbb{H}^3\big)
\times L^2\big(0,t_f;L^2(\Gamma)^{n\times n}\big)
\times L^2\big(0,t_f;L^2(\Gamma)^{n\times n}\big)$ to
$L^2\big(0,t_f;\tilde{\mathbb{V}}_0^n\big)
\times L^2\big(0,t_f;\tilde{\mathbb{V}}_0^n\big)$ is continuous.
\end{corollary}

\section{Main results}
In this section, we present the main result of this paper, Theorem \ref{main},
and an example of its application. Theorem \ref{main} gives a sufficient
condition for the well-posedness in the weak sense of the coupled
parabolic-elliptic initial-boundary value problem \eqref{eq1}-\eqref{eq7}
describing the fully coupled TCP model. This condition depends on physical
 parameters of the system, coupling vectors, and the Korn constant which in
turn depends only on the shape of a domain. Evaluation of the Korn constant
becomes necessary for each specific domain and is essential for the proposed
well-posedness theory. For instance, in the case of a two-dimensional
annular region that does not experience pure rotation, the Korn constant
can be expressed explicitly in terms of the ratio of the inner and outer radii
of the domain. As a consequence, for such regions, the sufficient condition
for the well-posedness of the problem can be formulated in terms of physical
parameters, coupling vectors, and the ratio of the inner and outer radii.
This result is presented in Corollary \ref{corollary2}. The sufficient conditions
presented in Theorem \ref{main} and Corollary \ref{corollary2} also provide
 compatibility conditions on the coupling vectors.

\begin{theorem} \label{main}
Under Assumptions \ref{assm1}, \ref{assm2}, and \ref{assm3},
 the coupled TCP system \eqref{eq1}-\eqref{eq7} admits a unique weak solution
\begin{equation}
\begin{gathered}
(\mathbf{V}, \mathbf{u}) \in L^2\big(0,t_f; \mathbb{V}^3\big)
\times L^2\big(0,t_f; \tilde{\mathbb{V}}_0^n\big)  \\
\text{with } (\dot{\mathbf{V}}, \dot{\mathbf{u}})
\in L^2\big(0,t_f; \mathbb{H}^3\big)
\times L^2\big(0,t_f; \tilde{\mathbb{V}}_0^n\big)
\end{gathered}\label{eq60}
\end{equation}
if
\begin{equation}
c_E-2\Big(2n+\frac{3n}{\lambda_{\rm min}^2}\Big)^{1/2}\|P\|
\, |P^{-1}M^{-1}\mathbf{b}_0| \, |\mathbf{b}_1| >0 \label{eq61}
\end{equation}
where $c_E=c_E(\Omega)>0$ is the Korn constant, $n$ is the dimension of
the problem, $M^{-1}A=PDP^{-1}$,
$D=\operatorname{diag}(\lambda_1, \lambda_2, \lambda_3)$, and
$\mathbf{b}_0$ and $\mathbf{b}_1$ are coupling vectors.
Moreover, this solution depends continuously on the boundary data
$\mathbf{V}_B$ and $\dot{\mathbf{V}}_B$, the initial data
$\mathbf{V}_I$ and $\mathbf{V}_B(0)$, the applied boundary stress
 $\hat{\boldsymbol{\sigma}}$, and its time derivative
$\dot{\hat{\boldsymbol{\sigma}}}$.
\end{theorem}

\begin{proof}
Property \eqref{eq60} of the solution $(\mathbf{V}, \mathbf{u})$ and the
continuous dependence of the solution on the data follow immediately
from Theorem \ref{theorem1}, Lemma \ref{lemma2}, and Corollary \ref{corollary1}.
 We only need to prove that  condition \eqref{eq61} is sufficient for
 well-posedness of the system.

Since, for each $t \in [0,t_f)$, $\dot{\mathbf{u}}\in \tilde{\mathbb{V}}_0^n$,
from \eqref{eq52}, the trace theorem, and the Poincar\'e inequality we have
\begin{align*}
a_E(\dot{\mathbf{u}},\dot{\mathbf{u}})
&=\int_{\Omega} (\mathbf{b}_1 \cdot \dot{\mathbf{V}})\nabla \cdot\dot{\mathbf{u}}\mathrm{d}\Omega + \int_{\Gamma \setminus \Gamma_F}(\dot{\hat{\boldsymbol{\sigma}}}\mathbf{n})\cdot \dot{\mathbf{u}}\mathrm{d} \Gamma  \\
%
&\leq \|\mathbf{b}_1 \cdot \dot{\mathbf{V}}\|_{L^2(\Omega)}\|\nabla
\cdot\dot{\mathbf{u}}\|_{L^2(\Omega)}+\|\dot{\hat{\boldsymbol{\sigma}}}
\mathbf{n}\|_{L^2(\Gamma)^n}\|\dot{\mathbf{u}}\|_{L^2(\Gamma)^n} \\
& \leq \sqrt{n}|\mathbf{b}_1\||\dot{\mathbf{V}}\|_{\mathbb{H}^3}
 \|\nabla \dot{\mathbf{u}}\|_{\mathbb{H}^n} 
 + \hat{C}\|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(\Gamma)^{n\times n}}
 \|\nabla \dot{\mathbf{u}}\|_{\mathbb{H}^n}
\end{align*}
where $\hat{C}>0$ is a constant. Utilizing \eqref{eq57}, the last inequality yields
\begin{equation*}
c_E\|\nabla \dot{\mathbf{u}}\|_{\mathbb{H}^n} \leq \sqrt{n}|\mathbf{b}_1|
\, \|\dot{\mathbf{V}}\|_{\mathbb{H}^3} 
+ \hat{C}\|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(\Gamma)^{n\times n}}, \quad \forall t\in [0,t_f)
\end{equation*}
and therefore,
\begin{equation}
c_E\|\nabla \dot{\mathbf{u}}\|_{L^2(0,t_f;\mathbb{H}^n)} 
\leq \sqrt{2n}|\mathbf{b}_1\||\dot{\mathbf{V}}\|_{L^2(0,t_f;\mathbb{H}^3)} 
+ \sqrt{2}\hat{C}\|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(0,t_f;L^2(\Gamma
 )^{n\times n})} \label{eq62}
\end{equation}
Substituting the a priori estimate \eqref{eq49} into \eqref{eq62}, we have
\begin{equation}
\begin{aligned}
&c_E\|\nabla \dot{\mathbf{u}}\|_{L^2(0,t_f;\mathbb{H}^n)} \\
&\leq 2\Big(2n+\frac{3n}{\lambda_{\rm min}^2}\Big)^{1/2}
\|P\|\,|P^{-1}M^{-1}\mathbf{b}_0|\,|\mathbf{b}_1|
\, \|\nabla \dot{\mathbf{u}}\|_{L^2(0,t_f;\mathbb{H}^n)} + \Theta 
\end{aligned} \label{eq63}
\end{equation}
where $\Theta>0$ is a constant independent of $\mathbf{u}$ and $\dot{\mathbf{u}}$.
Applying the Hadamard well-posedness condition for operator equations 
\cite{Ivanov} to \eqref{eq63} yields the sufficient condition \eqref{eq61}.
\end{proof}

\begin{corollary} \label{corollary2}
For an annular domain $ \Omega \subset \mathbb{R}^2$ of inner radius $R_B$ 
and outer radius $R_F$ with $ \delta=\frac{R_B}{R_F} \leq e^{-1}$, 
under the side condition on the displacement $\mathbf{u}=u\mathbf{i}+v\mathbf{j}$,
\begin{equation}
\int_{\Omega} \Big( \frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}\Big) 
\mathrm{d}\Omega =0 \label{eq64}
\end{equation}
the coupled TCP system \eqref{eq1}-\eqref{eq7} has a unique weak solution if
\begin{equation}
\frac{G}{2}\Big[1-\Big(\frac{3}{\delta^{-2}+1+\delta^2}\Big)^{1/2}\Big] 
- 2\Big(4+\frac{6}{\lambda_{\rm min}^2}\Big)^{1/2} 
\|P\|\, |P^{-1}M^{-1}\mathbf{b}_0| \, |\mathbf{b}_1| >0 \label{eq65}
\end{equation}
where $G$ is the shear modulus, $M^{-1}A=PDP^{-1}$,
 $D=\operatorname{diag}(\lambda_1, \lambda_2, \lambda_3)$, and $\mathbf{b}_0$ 
and $\mathbf{b}_1$ are coupling vectors.
\end{corollary}

\begin{proof}
Under the given conditions, the following inequality holds \cite{Horgan}:
\begin{equation*}
\frac{D(\dot{\mathbf{u}})}{S(\dot{\mathbf{u}})} 
\leq 4\Big[1-\Big(\frac{3}{\delta^{-2}+1+\delta^2}\Big)^{1/2}\Big]^{-1}
\end{equation*}
where $D(\dot{\mathbf{u}}) = \|\nabla\dot{\mathbf{u}}\|_{\mathbb{H}^2}^2$ and 
$S(\dot{\mathbf{u}})= {\int_{\Omega}} \sum_{i,j=1}^2\varepsilon_{ij}
(\dot{\mathbf{u}})^2 \mathrm{d} \Omega$. Therefore,
\begin{equation}
S(\dot{\mathbf{u}}) \geq \frac{1}{4}
\Big[1-\Big(\frac{3}{\delta^{-2}+1+\delta^2}\Big)^{1/2}\Big] 
\|\nabla\dot{\mathbf{u}}\|_{\mathbb{H}^2}^2 \label{eq66}
\end{equation}
For homogeneous and isotropic medium, the coefficients of elasticity 
$ a_{ijkl}$, $i,j,k,l=1,2$, in terms of the Lam\'e parameters 
$\lambda$ and $\mu$ are $a_{1111}=a_{2222}= 2\mu + \lambda$, 
$a_{1122}=a_{2211}=\lambda$, $a_{1212}=a_{2121}=a_{1221}=a_{2112}=\mu$ and, 
from \eqref{eq11} and \eqref{eq12}, we have
\begin{equation}
\begin{aligned}
a_E(\dot{\mathbf{u}},\dot{\mathbf{u}})
&= \int_{\Omega} 2\mu(\varepsilon_{11}^2+\varepsilon_{12}^2
 +\varepsilon_{21}^2+\varepsilon_{22}^2)
 + \lambda (\varepsilon_{11}+\varepsilon_{22})^2 \mathrm{d} \Omega  \\
& \geq 2\mu\int_{\Omega} \varepsilon_{11}^2+\varepsilon_{12}^2
 +\varepsilon_{21}^2+\varepsilon_{22}^2 \mathrm{d}\Omega  \\
& = 2\mu S(\dot{\mathbf{u}})
\end{aligned} \label{eq67}
\end{equation}
where $\mu=G$ is the shear modulus. Equations \eqref{eq57}, \eqref{eq66},
and \eqref{eq67} give the Korn constant
\begin{equation*}
c_E=\frac{G}{2}\Big[1-\Big(\frac{3}{\delta^{-2}+1+\delta^2}\Big)^{1/2}\Big]
\end{equation*}
Applying this constant to the sufficient condition \eqref{eq61} yields
 \eqref{eq65}.
\end{proof}

\begin{remark} \rm
The purpose of the condition \eqref{eq64} is to eliminate a pure rotation.
\end{remark}

\begin{thebibliography}{00}

\bibitem{Abousleiman} Y. Abousleiman, S. Ekbote;
\emph{Solutions for the inclined borehole in a porothermoelastic transversely 
isotropic medium}, ASME J. Appl. Mech. 72 (2005), 102-114.

\bibitem{Belotserkovets} A. Belotserkovets, J. Prevost;
\emph{Thermoporoelastic response of a fluid-saturated porous sphere: 
An analytical solution}, Int. J. Eng. Sci. 49 (2011), 1415-1423.

\bibitem{Biot1} M. A. Biot;
\emph{General theory of three-dimensional consolidation}, 
J. Appl. Phys. 12 (1941), 155-169.

\bibitem{Biot2} M. A. Biot;
\emph{Thermoelasticity and irreversible thermodynamics}, 
J. Appl. Phys. 28 (1956) 240-253.

\bibitem{Coussy} J. Coussy;
\emph{Poromechanics}, John Wiley \&  Sons Ltd., Chichester, 2004.

\bibitem{Diek1} A. Diek, L. White, J. C. Roegiers, K. Bartko, F. Chang;
\emph{A fully coupled thermoporoelastic model for drilling in chemically 
active formations}, 45th US Rock Mechanics, 
Geomechanics Symposium 2011, 2 (2011), 920-931.

\bibitem{Diek2} A. Diek, L. White, J. C. Roegiers;
\emph{A fully Coupled thermoporoelastic model for drilling in geological 
formations}, Geothermal Resources Council Transactions 35 (2011) 159-164.

\bibitem{Diek3} A. Diek, L. White, J. C. Roegiers, D. Blankenship;
\emph{A fully coupled thermoporoelastic model for drilling in HPHT formations,
 Harmonising Rock Engineering and the Environment}, Proceedings of the 
12th ISRM International Congress on Rock Mechanics< (2012) 1303-1308.

\bibitem{Diek4} A. Diek, L. White, J. C. Roegiers, J. Moore, J. D. McLennan;
\emph{Borehole preconditioning of geothermal wells for enhanced geothermal
 system reservoir development}, Proceedings, 37 Workshop on Geothermal Reservoir 
Eng. (2012), 1-11.

\bibitem{Duvaut_Lions} G. Duvaut, J. L. Lions;
\emph{Inequalities in Mechanics and Physics}, Springer, Berlin, 1976.

\bibitem{Ghassemi3} A. Ghassemi, A. Diek;
\emph{Linear chemo-poroelasticity for swelling shales: Theory and application}, 
J. Petrol. Sci. Eng., 38 (2003), 199-212.

\bibitem{Horgan} C. O. Horgan;
\emph{Korn's inequalities and their applications in continuum mechanics}, 
SIAM Review, 37 (1995), 491-511.

\bibitem{Ivanov} V. K. Ivanov, V. V. Vasin, V. P. Tanana;
\emph{Theory of Linear Ill-Posed Problems and Its Applications}, 
Inverse and Ill-Posed Problems Series, VSP, Utrecht, 2002.

\bibitem{Kodashima1} T. Kodashima, M. Kurashige;
\emph{Thermal stresses in a fluid-saturated poroelastic hollow sphere}, 
J. Therm. Stresses 19 (1996) 139-151.

\bibitem{Kodashima2} T. Kodashima, M. Kurashige;
\emph{Active cooling and thermal stresses reduction in a poroelastic 
hollow sphere}, J. Therm. Stresses, 20 (1997), 389-405.

\bibitem{Kurashige} M. Kurashige;
\emph{A thermoelastic theory of fluid-filled porous materials}, 
Int. J. Solids Struct. 25 (1989), 1039-1052.

\bibitem{Lions1} J.-L. Lions;
\emph{Optimal Control of Systems Governed by Partial Differential Equations}, 
Springer-Verlag, New York - Heidelberg - Berlin, 1971.

\bibitem{McTigue} D. McTigue;
\emph{Thermoelastic response of fluid-saturated porous rock}, 
J. Geophys. Res. 91 (1986) 9533-9542.

\bibitem{Pao} W. K. S. Pao, R. W. Lewis, I. Masters;
\emph{A fully coupled hydro-thermo-poro-mechanical model for black oil reservoir 
simulation}, Int. J. Numer. Anal. Meth. Geomech., 25 (2001), 1229-1256.

\bibitem{Rafieepour} S. Rafieepour, H. Jalayeri, C. Ghotbi, M. R. Pishvaie;
\emph{Simulation of wellbore stability with thermo-hydro-chemo-mechanical
 coupling in troublesome formations: An example from Ahwaz oil field},
 SW Iran, Arab. J. Geosci., 8 (2015), 379-396.

\bibitem{Schiffman} R. L. Schiffman;
\emph{A thermoelastic theory of consolidation}, Environmental and Geophysical 
Heat Transfer, 4 (1971), 78-84.

\bibitem{Showalter2} R. E. Showalter;
\emph{Diffusion in deformable media}, IMA Volumes in Mathematics and 
its Applications, 131 (2000) 115-130.

\bibitem{Showalter3} R. E. Showalter;
\emph{Diffusion in poro-elastic media}, J. Math. Anal. Appl., 251 (2000), 310-340.

\bibitem{Steinbach} O. Steinbach;
\emph{Numerical Approximation Methods for Elliptic Boundary Value Problems: 
Finite and Boundary Elements}, Springer, New York, 2008.

\bibitem{Temam} R. Temam;
\emph{Navier-Stokes Equations: Theory and Numerical Analysis}, 
North-Holland, Amsterdam - New York - Oxford, 1977.

\bibitem{Terzaghi1} K. Terzaghi;
\emph{Die berechnung der durchlassigkeitzifer des tones aus
 dem verlauf der hydrodynamischen spannungserscheinungen}, 
Akademie der Wissenschaften, Mathematisch-naturwissenschaftliche Klasse, 
Part IIa 132 (1923), 125-138.

\bibitem{Terzaghi2} K. von Terzaghi;
\emph{Theoretical Soil Mechanics}, Wiley, New York, 1943.

\bibitem{Wang} X. Wang, A. Ghassemi;
\emph{A 3D thermal-poroelastic model for geothermal reservoir stimulation}, 
Proceedings of the 37th Workshop on Geothermal Reservoir Engineering, 2012 (2012),
 SGP-TR-194.

\bibitem{Zhou} X. Zhou, A. Ghassemi;
\emph{Finite element analysis of coupled chemo-poro-thermo-mechanical 
effects around a wellbore in swelling shale}, Int. J. Rock Mech. Min. Sci.,
 46 (2009), 769-778.

\end{thebibliography}

\end{document}