\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{multirow,hhline}

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

\begin{document}
\title[\hfilneg EJDE-2017/137\hfil Fully coupled thermo-chemo-poroelastic system]
{Well-Posedness of a fully coupled thermo-chemo-poroelastic system with
applications to petroleum rock mechanics}

\author[T. Malysheva, L. W. White \hfil EJDE-2017/137\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}

\dedicatory{Communicated by Ralph Showalter}

\thanks{Submitted February 27, 2017. Published May 24, 2017.}
\subjclass[2010]{35D30, 35E99, 35G16, 35Q74, 35Q86}
\keywords{Parabolic-elliptic system; poroelasticity; thermo-poroelasticity; 
\hfill\break\indent  thermo-chemo-poroelasticity; existence; uniqueness; 
well-posedness}

\begin{abstract}
 We consider a system of fully coupled parabolic and elliptic equations
 constituting the general model of chemical thermo-poroelasticity for a
 fluid-saturated porous media. The main result of this paper is the
 developed well-posedness theory for the corresponding initial-boundary
 problem arising from petroleum rock mechanics applications.
 Using the proposed pseudo-decoupling method, we establish, subject to
 some natural assumptions imposed on matrices of diffusion coefficients,
 the existence, uniqueness, and continuous dependence on initial and
 boundary data of a weak solution to the problem. Numerical experiments
 confirm the applicability of the obtained well-posedness results for
  thermo-chemo-poroelastic models with real-data parameters.
\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}
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks

 \section{Introduction}

In this article we investigate a model describing fully coupled processes
of quasi-static elastic deformation and thermal, solute, and fluid diffusions
in porous media. This work is motivated by petroleum rock mechanics applications
dealing with drilling and borehole stability in high-temperature,
high-pressure chemically active rock formations. The system of coupled
diffusion and deformation equations is taken from Diek, White, and
Roegiers \cite{Diek1} and constitutes the general thermo-chemo-poroelasticity
theory of porous media saturated by a compressible and thermally expansible fluid.

The underlying equations are formulated in terms of 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)$ and will henceforth be referred to as
the TCpu system. In the case of a homogeneous and isotropic medium the
equations of the TCpu system take the  form:\\
thermal diffusion
\begin{equation}
\Lambda\dot{T}+\Sigma\dot{C}+\Phi\dot{p}-\frac{k^\text{T}}{T_F}\nabla^2T
-C_FD^\text{T}\rho_f\tilde{\Omega}\nabla^2C+K^\text{T}\nabla^2p
=-\zeta(\nabla\cdot\dot{\mathbf{u}}), \label{eq1.1}
\end{equation}
solute diffusion
\begin{equation}
\phi\dot{C}-C_FD^\text{T}\nabla^2T-D\nabla^2C+\frac{k\mathfrak{R}}{\eta}\nabla^2p=0,
 \label{eq1.2}
\end{equation}
fluid diffusion
\begin{equation}
\Gamma\dot{T}+\chi\dot{C}+\Psi\dot{p}+K^\text{T}\nabla^2T
+\frac{k \mathfrak{R}}{\eta} \rho_f\tilde{\Omega}\nabla^2C
-\frac{k}{\eta}\nabla^2p=-\alpha(\nabla\cdot\dot{\mathbf{u}}), \label{eq1.3}
\end{equation}
and the Navier-type elastic equation
\begin{equation}
\Big(K+\frac{G}{3}\Big)\nabla(\nabla\cdot\dot{\mathbf{u}})
+G\nabla^2\dot{\mathbf{u}}
=\tilde{\zeta}\nabla\dot{T}-\xi\nabla\dot{C}+\tilde{\alpha}\nabla\dot{p},
 \label{eq1.4}
\end{equation}
where the superscript dot $(\dot{~})$ denotes a time derivative. 
The description and values of physical constants are taken from \cite{Diek1} 
and presented in the appendix of the paper. These equations supplemented 
by appropriate initial and boundary conditions constitute an initial-boundary 
value problem (IBVP) defined in an open region $\Omega \subset \mathbb{R}^n, ~n=2,3$, 
exterior to the borehole with a sufficiently smooth boundary $\Gamma$. 
We assume that the outer (far-field) boundary $\Gamma_F \subset \Gamma$ 
of the region has a nonempty interior relative to $\Gamma$ and is specified by 
the following conditions: (i) the absolute temperature, solute mass fraction, 
pore pressure, and displacements are time-independent, and (ii) 
displacements and their velocities are negligibly small.

For the sake of convenience, we combine thermal, solute, and fluid diffusion
 equations into a single vector diffusion equation. To that end, we 
introduce the vector $\bar{\mathbf{V}}= [T ~C ~p]^T$, with the superscript $T$ 
meaning transpose, the matrices of diffusion coefficients
\begin{equation}
M=\begin{bmatrix} 
\Lambda & \Sigma & \Phi \\
0 & \phi & 0 \\
\Gamma & \chi & \Psi
\end{bmatrix}, \quad
A=\begin{bmatrix} 
\frac{k^\text{T}}{T_F} & C_FD^\text{T}\rho_f\tilde{\Omega} & -K^\text{T} \\
C_FD^\text{T} & D & -\frac{k\mathfrak{R}}{\eta} \\
-K^\text{T} & -\frac{k \mathfrak{R}}{\eta} \rho_f\tilde{\Omega} & \frac{k}{\eta}
\end{bmatrix} \label{eq1.5}
\end{equation}
and the coupling vectors
\begin{equation}
\mathbf{{b}}_0=\begin{bmatrix} 
\zeta \\
0 \\
\alpha
\end{bmatrix}, \quad
\mathbf{b}_1=\begin{bmatrix} 
\tilde{\zeta} \\
-\xi \\
\tilde{\alpha}
\end{bmatrix} = \begin{bmatrix} 
\zeta \\
0\\
\alpha
\end{bmatrix} + \begin{bmatrix} 
s_F\tilde{\omega} \\
-\xi \\
-\frac{\tilde{\omega}}{\rho_f}
\end{bmatrix} =:\mathbf{b}_0 +\mathbf{b}_d \,. \label{eq1.6}
\end{equation}
With the above notation, the IBVP for the fully coupled TCpu 
system \eqref{eq1.1}-\eqref{eq1.4} has the  form
\begin{gather}
M \dot{\bar{\mathbf{V}}} - A \nabla^{2} \bar{\mathbf{V}} 
=-\mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}), \quad \text{in } 
\Omega \times (0,t_f), \label{eq1.7}\\
\Big(K+\frac{G}{3}\Big)\nabla(\nabla\cdot\dot{\mathbf{u}})
+G\nabla^2\dot{\mathbf{u}}=\nabla(\mathbf{b}_1 \cdot \dot{\bar{\mathbf{V}}}), 
\quad \text{in } \Omega \times (0,t_f), \label{eq1.8}\\
\bar{\mathbf{V}}(\mathbf{x},t) = \mathbf{V}_B(\mathbf{x},t), \quad 
\text{on } \Gamma \times [0,t_f), \label{eq1.9}  \\
\mathbf{u}(\mathbf{x},t) = \mathbf{0}, \quad \text{on } \Gamma_F 
\times (0, t_f), \label{eq1.10}\\
\dot{\mathbf{u}}(\mathbf{x},t) = \mathbf{0}, \quad \text{on } 
\Gamma_F \times (0, t_f),\label{eq1.11} \\
\dot{\boldsymbol{\tau}} \mathbf{n} 
=\big((\mathbf{b}_1 \cdot \dot{\mathbf{V}})I
+\dot{\hat{\boldsymbol{\sigma}}}\big)
\mathbf{n}, \quad \text{on } \Gamma \setminus \Gamma_F \times (0,t_f),
 \label{eq1.12}\\
\bar{\mathbf{V}}(\mathbf{x},0) = \mathbf{V}_I(\mathbf{x}), \quad \text{in } 
\Omega, \label{eq1.13}
\end{gather}
where $t_f \in (0, \infty)$ stands for a final time, $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  time-independent 
on the far-field boundary $\Gamma_F$.

The aim of this paper is to establish existence, uniqueness, and continuous 
dependence on initial and boundary data of a weak solution to the fully 
coupled IBVP \eqref{eq1.7}-\eqref{eq1.13}.

One can find an enormous amount of literature concerning numerical techniques 
and simulations for models describing coupled thermal, chemical, hydraulic, 
and quasi-static elastic deformation processes in porous media relevant 
to petroleum rock mechanics applications. However, there are very 
few papers dealing with analytical solutions and well-posedness of such 
coupled problems. Typically, analytical solutions are derived under 
assumptions that some of the diffusion processes or couplings can be 
neglected \cite{Abousleiman, Belotserkovets, Coussy, Ekbote, Ghassemi1, 
Kodashima1, Kodashima2, McTigue}, and most well-posedness results are 
obtained for poroelastic and thermoelastic  systems only. The existence, 
 uniqueness, and regularity theory for linear Biot's consolidation models 
in poroelasticity and a coupled quasi-static problem in thermoelasticity 
was developed by Showalter \cite{Showalter1, Showalter2} from the theory of 
linear degenerate evolution equations in a Hilbert space. 
It should be noted that, with slight modifications, these results can be 
extended to a fully coupled thermo-poroelastic model. In the papers 
of Barucq, Madaune-Tort, and Saint-Macary \cite{Barucq1, Barucq2}, 
the existence and uniqueness of weak solutions to non-linear fully dynamic 
and quasi-static Biot's consolidation models of poroelasticity for either 
Newtonian or non-Newtonian fluid were established using Galerkin approximants.

The main difficulty of the problem under consideration is related to the 
complexity of cross-coupling mechanisms involved in the fully coupled 
model of thermo-chemo-poroelasticity. In contrast with the fully 
coupled thermo-poroelastic system, the matrices of diffusion coefficients
 $M$ and $A$ are non-symmetric and the diffusion equation \eqref{eq1.7}
cannot be rescaled to make the coupling vectors $\mathbf{b}_0$ and 
$\mathbf{b}_1$ equal; concsequently, the techniques presented in 
\cite{Barucq1, Barucq2, Showalter1, Showalter2} are not applicable to 
the problem \eqref{eq1.7}-\eqref{eq1.13}. In our preceding
paper \cite{Malysheva} we obtained sufficient conditions for 
Hadamard well-posedness of the system \eqref{eq1.7}-\eqref{eq1.13}.
 However, the question of finding necessary conditions for existence and 
uniqueness of a solution and its continuous dependence on data remains 
open and therefore, the well-posedness of the fully coupled TCpu 
problem \eqref{eq1.7}-\eqref{eq1.13} requires further investigation.

The novel contribution of the present paper consists of the following:

(i)
The well-posedness theory is developed for the IBVP 
\eqref{eq1.7}-\eqref{eq1.13} for the system of fully coupled partial
differential equations constituting the general thermo-chemo-poroelasticity 
theory of porous media saturated by a compressible and thermally expansible 
fluid. It is shown that, subject to some natural assumptions imposed on
 matrices of diffusion coefficients, the system \eqref{eq1.7}-\eqref{eq1.13}
admits a unique weak solution and this solution depends continuously on 
initial and boundary data.

(ii) A novel method that allows one to transform the coupled parabolic-elliptic
 IBVP \eqref{eq1.7}-\eqref{eq1.13} to an IBVP for a single implicit
 equation is presented. This method will be called the pseudo-decoupling 
method because it is based on the construction of operators that fully 
eliminate coupling terms by implicitly incorporating the elliptic IBVP 
for the elasticity system into the parabolic equation representing diffusion 
processes. This approach differs from those used before in the 
literature \cite{Showalter1, Showalter2} in the construction of operators 
that transform a coupled parabolic-elliptic system into an implicit system.
 In contrast with \cite{Showalter1, Showalter2} dealing with anti-symmetric
 coupling terms, the proposed method makes no assumption on the relation 
between the coupling vectors $\mathbf{b}_0$ and $\mathbf{b}_1$ and therefore, 
it is applicable to a wider variety of coupling mechanisms.


The article is organized as follows. In Section 2 we introduce the 
pseudo-decoupling method to transform the fully coupled
 parabolic-elliptic IBVP \eqref{eq1.7}-\eqref{eq1.13} to an IBVP
for a single implicit equation. The obtained IBVP for the implicit equation 
is then decomposed into a parabolic IBVP and an implicit IBVP with 
homogeneous boundary and initial conditions. The well-posedness of the 
parabolic IBVP is discussed in Section 2. In Section 3 we prove the 
well-posedness in a weak sense of the implicit IBVP with homogeneous 
boundary and initial conditions and a generalized source term. 
Section 4 contains the main results of this paper summarized in 
Theorem  \ref{thrm_main} that establishes the existence, uniqueness, 
and continuous dependence on initial and boundary data of a weak 
solution to the fully coupled IBVP \eqref{eq1.7}-\eqref{eq1.13}
for the TCpu system. In Section 5 we provide a numerical example illustrating 
the applicability of the obtained well-posedness results for the TCpu model 
with real data.


\subsection{Notation}

Let $\Omega$ be a bounded open domain in $\mathbb{R}^n$, $n=2,3$, with a 
sufficiently smooth boundary $\Gamma$, $\bar{\Omega}=\Omega \cup \Gamma$, 
and $\Gamma_F \subset \Gamma$ have a nonempty interior relative to $\Gamma$. 
Throughout this article, $\mathbf{x}=(x_1, x_2, \dots, x_n)\in\bar{\Omega}$, 
$\mathbf{u}=[u_1 ~u_2 \ldots u_n]^T$ is the general notation for an 
$\mathbb{R}^n$-valued function, the superscript $T$ means transpose, and 
$\partial_i$ stands for the partial derivative with respect to $x_i$. 
We introduce the following spaces of vector-valued functions:
$\mathbb{H}^n=L^2(\Omega)^n$ and $\mathbb{V}^n=H^1(\Omega)^n$ 
endowed with their standard norms denoted by $\|\cdot\|_{\mathbb{H}^n}$ and 
$\|\cdot\|_{\mathbb{V}^n}$, respectively; the space 
$\mathbb{V}_0^n=H_0^1(\Omega)^n$ with the norm
\begin{equation*}
\|\mathbf{u}\|_{0,n}=\Big[\sum_{k=1}^n\int_{\Omega} |\nabla u_k|^2
\mathrm{d}\Omega\Big]^{1/2}
\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}_0^n, ~n \in \mathbb{N}$.
 Let $-\infty \leq a < b \leq \infty$ and $\mathrm{X}$ be a Hilbert space. 
We denote by $L^2(a,b;\mathrm{X})$ the space of $L^2$-integrable 
functions from $[a,b]$ into $\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*}
The space $L^{\infty}(a,b;\mathrm{X})$ is the space of essentially 
bounded functions from $[a,b]$ into $\mathrm{X}$ equipped with the norm
\begin{equation*}
\|u\|_{L^\infty(a,b;\mathrm{X})}= \operatorname{ess\,sup}_{[a,b]}
\|u(t)\|_{X}\,.
\end{equation*}
We will occasionally use the Einstein summation convention: whenever
 an index is repeated once in the same term, it implies summation over 
the specified range of the index. For example, 
$\sum_{k=1}^n \sum_{l=1}^n a_{ijkl}\varepsilon_{kl}\equiv a_{ijkl}\varepsilon_{kl}$.

\section{Pseudo-decoupling method}

In this section we present the pseudo-decoupling method to transform the 
coupled parabolic-elliptic IBVP \eqref{eq1.7}-\eqref{eq1.13}
to an IBVP for a single implicit equation and discuss further decomposition 
of the implicit system.

\subsection{Elastic system}

Elastic deformation of the region $\Omega$ is characterized using 
the linearized strain tensor
\begin{equation}
\boldsymbol{\varepsilon}(\mathbf{u})
= \frac{1}{2}\big(\nabla\mathbf{u}+\nabla\mathbf{u}^T\big), \label{eq2.1}
\end{equation}
where $\mathbf{u}$ is the displacement vector. The material properties 
of rock are described by the relation between the stress tensor 
$\boldsymbol{\tau}=[\tau_{ij}]$ and the strain tensor 
$\boldsymbol{\varepsilon}=[\varepsilon_{ij}]$, the generalized Hooke's law,
\begin{equation}
\tau_{ij}=a_{ijkl}\varepsilon_{kl}(\mathbf{u}), \label{eq2.2}
\end{equation}
where $a_{ijkl}$ are the coefficients of elasticity independent of the 
strain tensor, with the properties of symmetry
\begin{equation}
a_{ijkl}=a_{jilk}=a_{klij} \label{eq2.3}
\end{equation}
and of ellipticity: there exists a constant $\alpha_0>0$ such that
\begin{equation*}
a_{ijkl}\epsilon_{ij}\epsilon_{kl} \geq \alpha_0 \epsilon_{ij}\epsilon_{ij}
\end{equation*}
for all symmetric $\mathbf{\epsilon}\in\mathbb{R}^{n\times n}$. 
In the case of a homogeneous and isotropic medium, the stress-strain 
relation \eqref{eq2.2} in terms of the bulk modulus $K$ and the
shear modulus $G$ takes the form
\begin{equation}
\tau_{ij} = 2G\varepsilon_{ij} + \Big(K-\frac{2G}{3}\Big)
\varepsilon_{kk}\delta_{ij}\,. \label{eq2.4}
\end{equation}

According to the principle of minimum total potential energy, 
among all admissible displacements satisfying the boundary conditions
 \eqref{eq1.10} and \eqref{eq1.11}, the actual displacement
that the region $\Omega$ undergoes is the one that minimizes the total 
potential energy $\mathcal{V}$ of the elastic system 
\eqref{eq1.8}, \eqref{eq1.10}-\eqref{eq1.12}:
\begin{equation*}
\mathcal{V}(\mathbf{u})=\mathcal{V}_S(\mathbf{u})-W_b(\mathbf{u})-W_S(\mathbf{u}),
\end{equation*}
where
\begin{equation*}
\mathcal{V}_S(\mathbf{u})=\frac{1}{2}\int_\Omega \tau_{ij}(\mathbf{u})
\varepsilon_{ij}(\mathbf{u}) \mathrm{d} \Omega
\end{equation*}
is the elastic energy of the system;
\begin{equation*}
W_b(\mathbf{u})=\int_\Omega (\mathbf{b}_1 \cdot \bar{\mathbf{V}}) 
(\nabla \cdot \mathbf{u}) \mathrm{d} \Omega
\end{equation*}
is the work done by body forces due to the absolute temperature, 
solute mass fraction, and pore pressure; and
\begin{equation*}
W_S(\mathbf{u})=\int_{\Gamma}(\boldsymbol{\hat{\sigma}}\mathbf{n})
\cdot\mathbf{u} \mathrm{d} \Gamma
\end{equation*}
is the  work done by applied boundary stress.

Let us define a bilinear form 
$a_E:\mathbb{V}^n \times \mathbb{V}^n \to \mathbb{R}$ by
\begin{equation}
a_E(\mathbf{u},\boldsymbol{\Phi}) 
= \int_\Omega \tau_{ij}(\mathbf{u})\varepsilon_{ij}(\boldsymbol{\Phi}) 
\mathrm{d}\Omega\,. \label{eq2.5}
\end{equation}
Then the total potential energy takes the form
\begin{equation*}
\mathcal{V}(\mathbf{u})=\frac{1}{2}a_E(\mathbf{u},\mathbf{u})
-\int_\Omega (\mathbf{b}_1 \cdot \bar{\mathbf{V}}) (\nabla \cdot \mathbf{u}) 
\mathrm{d} \Omega - \int_{\Gamma}(\boldsymbol{\hat{\sigma}}\mathbf{n})
\cdot\mathbf{u} \mathrm{d} \Gamma\,.
\end{equation*}

Using the principle of minimum total potential energy, it was shown 
in \cite{Malysheva} that the elastic system \eqref{eq1.8},
\eqref{eq1.10}-\eqref{eq1.12} is equivalent to
\begin{gather*}
\Big(K+\frac{G}{3}\Big)\nabla(\nabla \cdot \mathbf{u})+G\nabla^2\mathbf{u} 
= \nabla(\mathbf{b}_1 \cdot \bar{\mathbf{V}}), \quad 
\text{in } \Omega \times (0,t_f), \\
\mathbf{u} =\mathbf{0}, \quad \text{on } \Gamma_F \times (0,t_f) \\
\boldsymbol{\tau} \mathbf{n}=\big((\mathbf{b}_1 \cdot \bar{\mathbf{V}})I
+\boldsymbol{\hat{\sigma}}\big)\mathbf{n}, \quad \text{on } 
\Gamma \setminus \Gamma_F \times (0,t_f),
\end{gather*}
and the variational form of the above system is 
\begin{equation}
a_E(\mathbf{u},\boldsymbol{\Phi})-\int_\Omega (\mathbf{b}_1 
\cdot \bar{\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{eq2.6}
\end{equation}

We will be working under the following assumption on the applied boundary 
stress tensor $\boldsymbol{\hat{\sigma}}$.

\begin{assumption}  \label{assm2.1} \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}

\subsection{Diffusion system}

We observe from \eqref{eq1.5} that the matrix $A$ can be written
 as a product of a symmetric matrix $A_0$ and a diagonal matrix $R$ in the form
\begin{equation*}
A=A_0R,
\end{equation*}
where
\begin{equation}
A_0=\begin{bmatrix}
\frac{k^\text{T}}{T_F} & C_FD^\text{T} & -K^\text{T} \\[4pt]
C_FD^\text{T} & \frac{D}{\rho_f\tilde{\Omega}} & -\frac{k\mathfrak{R}}{\eta} \\[4pt]
-K^\text{T} & -\frac{k \mathfrak{R}}{\eta}  & \frac{k}{\eta}
\end{bmatrix}, 
\quad
R= \begin{bmatrix} 
1 & 0 & 0 \\
0 & \rho_f\tilde{\Omega} & 0 \\
0 & 0  & 1
\end{bmatrix}, \quad \rho_f\tilde{\Omega}>0\,. \label{eq2.7}
\end{equation}
Define the vector
\begin{equation}
\bar{\mathbf{V}}_R=R\bar{\mathbf{V}}\,. \label{eq2.8}
\end{equation}
This transformation leads to the IBVP equivalent to the diffusion 
system \eqref{eq1.7}, \eqref{eq1.9}, and \eqref{eq1.13}:
\begin{gather}
M_R \dot{\bar{\mathbf{V}}}_R - A_0 \nabla^{2} \bar{\mathbf{V}}_R 
+ \mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}) = \mathbf{0}, \quad 
\text{in } \Omega \times (0,t_f), \label{eq2.9}\\
\bar{\mathbf{V}}_R(\mathbf{x},t) = R\mathbf{V}_B(\mathbf{x},t), \quad 
\text{on } \Gamma \times [0,t_f), \label{eq2.10}  \\
\bar{\mathbf{V}}_R(\mathbf{x},0) = R\mathbf{V}_I(\mathbf{x}), \quad \text{in } 
\Omega, \label{eq2.11}
\end{gather}
where $M_R=MR^{-1}$ is a non-symmetric matrix.

We shall be using the following assumption imposed on the matrices of 
diffusion coefficients $M_R$ and $A_0$.

\begin{assumption}  \label{assm2.2} \rm
The matrices $M_R=[m_{ij}]_{i,j=1}^3$ and $A_0=[a_{ij}]_{i,j=1}^3$ 
satisfy the following conditions:
\begin{itemize}
\item[(i)] $0 < m_{ii} < b_{Ri}^2, ~~i=1,2,3$, where $\mathbf{b}_R^T = [b_{R1} ~b_{R2} ~b_{R3}]^T = \mathbf{b}_d^T R^{-1}$.
\item[(ii)] $0 < (m_{ij}+m_{ji})^2 < m_{ii}m_{jj}, ~~i,j=1,2,3.$
\item[(iii)] The matrices $\frac{1}{2}\big(M_R+M_R^T\big)$ and $A_0$ are 
positive definite.
\end{itemize}
\end{assumption}

\subsection{Pseudo-decoupling method}

One may observe from \eqref{eq2.1}-\eqref{eq2.3} and \eqref{eq2.5} that
\begin{equation*}
a_E(\mathbf{u},\mathbf{v})=a_E(\mathbf{v},\mathbf{u}), \quad \forall 
\mathbf{u}, \mathbf{v} \in \mathbb{V}^n
\end{equation*}
and, for every $\mathbf{u}, \mathbf{v} \in \tilde{\mathbb{V}}_0^n$,
\begin{equation*}
|a_E(\mathbf{u},\mathbf{v})|\leq n \max_{i,j,k,l} 
\{a_{ijkl}\} \|\mathbf{u}\|_{0,n}\|\mathbf{v}\|_{0,n}\,.
\end{equation*}
Also, Korn's inequality and its consequence 
\cite[Chapter III, Theorems 3.1 and 3.3]{Duvaut_Lions} 
imply that there exists a constant $\gamma=\gamma(\Omega)>0$ such that
\begin{equation}
a_E(\mathbf{u},\mathbf{u}) \geq \gamma \|\mathbf{u}\|_{0,n}^2, \quad 
\forall \mathbf{u}\in \tilde{\mathbb{V}}_0^n \,. \label{eq2.22}
\end{equation}
Thus, the bilinear form $a_E(\cdot, \cdot)$ is symmetric, continuous, 
and coercive on $\tilde{\mathbb{V}}_0^n$. 

The above observation allows us to define the linear operator 
$E_0:\mathbb{H}\to \tilde{\mathbb{V}}_0^n$ by
\begin{equation}
a_E(E_0\varphi,\boldsymbol{\Phi})
=\int_\Omega \varphi(\nabla\cdot\boldsymbol{\Phi}) \mathrm{d} \Omega, 
\quad \forall \boldsymbol{\Phi}\in\tilde{\mathbb{V}}_0^n \label{eq2.23}
\end{equation}
and the function $\mathbf{u_B}=\mathbf{u_B}(\boldsymbol{\hat{\sigma}})
\in\tilde{\mathbb{V}}_0^n$ by
\begin{equation}
a_E(\mathbf{u}_B,\boldsymbol{\Phi}) 
= \int_{\Gamma} (\boldsymbol{\hat{\sigma}}\mathbf{n})\cdot\boldsymbol{\Phi}
 \mathrm{d} \Gamma, \quad \forall \boldsymbol{\Phi}\in\tilde{\mathbb{V}}_0^n.
 \label{eq2.24}
\end{equation}
The operator $E_0$ is continuous. Indeed, from \eqref{eq2.22}
and \eqref{eq2.23}, for every $\varphi\in\mathbb{H}$,
\begin{equation*}
\gamma \|E_0\varphi\|_{0,n}^2 \leq a_E(E_0\varphi,E_0\varphi) 
\leq \|\varphi\|_{\mathbb{H}} \|\nabla\cdot E_0\varphi\|_{\mathbb{H}} 
\leq \sqrt{n}\|\varphi\|_{\mathbb{H}}\|E_0\varphi\|_{0,n}\,.
\end{equation*}
Thus,
\begin{equation}
\|E_0\varphi\|_{0,n} \leq \frac{\sqrt{n}}{\gamma}\|\varphi\|_{\mathbb{H}}, 
\quad \forall \varphi\in\mathbb{H}\,. \label{eq2.25}
\end{equation}

Applying \eqref{eq2.8}, \eqref{eq2.23}, and \eqref{eq2.24} to
\eqref{eq2.6} gives, for every $\boldsymbol{\Phi}\in\tilde{\mathbb{V}}_0^n$,
\begin{align*}
a_E(\mathbf{u},\boldsymbol{\Phi}) 
= & \int_\Omega (\mathbf{b}_1^T R^{-1}\bar{\mathbf{V}}_R)(\nabla 
\cdot \boldsymbol{\Phi})\mathrm{d} \Omega 
+ \int_{\Gamma \setminus \Gamma_F}(\boldsymbol{\hat{\sigma}}
\mathbf{n})\cdot\boldsymbol{\Phi}\mathrm{d} \Gamma \\
= &a_E\big(E_0(\mathbf{b}_1^T R^{-1}\bar{\mathbf{V}}_R),
\boldsymbol{\Phi}\big) + a_E(\mathbf{u}_B,\boldsymbol{\Phi})\,.
\end{align*}
It follows that
\begin{equation}
\mathbf{u}=E_0(\mathbf{b}_1^T R^{-1}\bar{\mathbf{V}}_R) 
+ \mathbf{u}_B, \label{eq2.26}
\end{equation}
and hence,
\begin{equation}
\nabla\cdot\dot{\mathbf{u}}=\nabla\cdot E_0(\mathbf{b}_1^T R^{-1}
\dot{\bar{\mathbf{V}}}_R) +\nabla\cdot\dot{\mathbf{u}}_B \label{eq2.27}
\end{equation}
We now introduce the linear operator $\hat{E}:\mathbb{H}^3\to\mathbb{H}^3$ by
\begin{equation}
\hat{E}\boldsymbol{\psi}=\begin{bmatrix}
\nabla\cdot (E_0\psi_1) \\
\nabla\cdot (E_0\psi_2) \\
\nabla\cdot (E_0\psi_3)
\end{bmatrix}\,. \label{eq2.28}
\end{equation}
The continuity of $E_0$ implies the continuity of $\hat{E}$, 
and from \eqref{eq2.25} it is straightforward to show that
\begin{equation}
\|\hat{E}\boldsymbol{\psi}\|_{\mathbb{H}^3} 
\leq \frac{n}{\gamma}\|\boldsymbol{\psi}\|_{\mathbb{H}^3}, \quad 
\forall \boldsymbol{\psi}\in\mathbb{H}^3\,. \label{eq2.29}
\end{equation}
With the above notation, \eqref{eq2.27} can be rewritten as
\begin{equation}
\nabla\cdot\dot{\mathbf{u}}=\mathbf{b}_1^T R^{-1}\hat{E}\dot{\bar{\mathbf{V}}}_R
 + \nabla\cdot\dot{\mathbf{u}}_B\,. \label{eq2.30}
\end{equation}
Substituting \eqref{eq2.30} into \eqref{eq2.9}, we obtain the  implicit
equation
\begin{equation}
\big[M_R+\mathbf{b}_0 \mathbf{b}_1^TR^{-1}\hat{E}\big]\dot{\bar{\mathbf{V}}}_R
-A_0\nabla^2 \bar{\mathbf{V}}_R + \mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}_B)
 = \mathbf{0}, \quad \text{in } \Omega \times (0,t_f)\,. \label{eq2.31A}
\end{equation}
The above equation supplemented with the boundary and initial conditions
 \eqref{eq2.10} and \eqref{eq2.11} yields the IBVP for a single 
implicit equation,
\begin{gather*}
\big[M_R+\mathbf{b}_0 \mathbf{b}_1^TR^{-1}\hat{E}\big]
\dot{\bar{\mathbf{V}}}_R-A_0\nabla^2 \bar{\mathbf{V}}_R 
+ \mathbf{b}_0 (\nabla \cdot \dot{\mathbf{u}}_B) 
= \mathbf{0}, \quad \text{in }\Omega \times (0,t_f)\,, \\ %\quad{\eqref{eq2.31A}}\\
\bar{\mathbf{V}}_R(\mathbf{x},t) = R\mathbf{V}_B(\mathbf{x},t), \quad 
\text{on } \Gamma \times [0,t_f)\,, \\ %\quad {\eqref{eq2.10}} \\
\bar{\mathbf{V}}_R(\mathbf{x},0) = R\mathbf{V}_I(\mathbf{x}), 
\quad \text{in } \Omega \\ %\quad{\eqref{eq2.11}}
\end{gather*}
that is equivalent to the fully coupled parabolic-elliptic IBVP 
\eqref{eq1.7}-\eqref{eq1.13}. We will denote this IBVP by $\mathcal{P}0$.

By the superposition principle, the problem $\mathcal{P}0$ can be 
decomposed into the following two subproblems: the (autonomous)
 parabolic problem, denoted by $\mathcal{P}1$,
\begin{equation} \label{eq2.16}
\begin{gathered}
M_R\dot{\bar{\mathbf{W}}}-A_0\nabla^2\bar{\mathbf{W}}
= \mathbf{0}, \quad \text{in } \Omega \times (0,t_f),  \\
\bar{\mathbf{W}}(\mathbf{x},t)=R\mathbf{V}_B(\mathbf{x},t), \quad
\text{on } \Gamma \times [0,t_f),  \\
\bar{\mathbf{W}}(\mathbf{x},0)=R\mathbf{V}_I(\mathbf{x}), \quad \text{in }
\Omega,
\end{gathered}
\end{equation}
and the implicit problem, denoted by $\mathcal{P}2$,
\begin{gather}
\Big[M_R+\mathbf{b}_0 \mathbf{b}_1^TR^{-1}\hat{E}\Big] \dot{\mathbf{V}}
 -A_0\nabla^2 \mathbf{V} = \mathbf{F}, \quad \text{in } \Omega \times (0,t_f),
 \label{eq2.31}\\
\mathbf{V}(\mathbf{x},t) = \mathbf{0}, \quad \text{on } \Gamma \times [0,t_f),
\nonumber \\
\mathbf{V}(\mathbf{x},0) = \mathbf{0}, \quad \text{in } \Omega ,\label{eq2.14}
\end{gather}
where
\begin{equation}
\mathbf{F}=-\mathbf{b}_0 \mathbf{b}_1^T R^{-1}\hat{E}\dot{\bar{\mathbf{W}}}
-\mathbf{b}_0 (\nabla\cdot\dot{\mathbf{u}}_B), \label{eq2.32}
\end{equation}
$\dot{\bar{\mathbf{W}}}$ is the time derivative of a solution to the problem
$\mathcal{P}1$ and $\mathbf{u}_B$ is defined by \eqref{eq2.24}.
Then the solution to the problem $\mathcal{P}0$ takes the form
\begin{equation}
\bar{\mathbf{V}}_R=\mathbf{V}+\bar{\mathbf{W}} \label{eq2.18}
\end{equation}
and from \eqref{eq2.8}, \eqref{eq2.26}, and \eqref{eq2.18}, the solution
$(\bar{\mathbf{V}}, \mathbf{u})$ to the IBVP \eqref{eq1.7}-\eqref{eq1.13} is
\begin{gather}
\bar{\mathbf{V}} = R^{-1}(\mathbf{V}+\bar{\mathbf{W}}), \label{eq2.33}\\
\mathbf{u}=E_0(\mathbf{b}_1^T R^{-1}\mathbf{V})+ E_0(\mathbf{b}_1^T
 R^{-1}\bar{\mathbf{W}})+\mathbf{u}_B\,. \label{eq2.26A}
\end{gather}

\begin{remark} \rm \label{auxiliaryBC}
Equations \eqref{eq2.10}, \eqref{eq2.16}, and \eqref{eq2.18} yield the 
auxiliary boundary condition for the problem $\mathcal{P}2$:
\begin{equation}
\dot{\mathbf{V}}(\mathbf{x},t) = \mathbf{0}, \quad \text{on } 
\Gamma \times [0,t_f)\,.  \label{eq2.19}
\end{equation}
\end{remark}

\begin{remark} \rm
The autonomous IBVP $\mathcal{P}1$ can further be decomposed into the 
elliptic boundary value problem, denoted by  $\mathcal{P}1.1$, with time
 $t\in[0,t_f]$ as a parameter
\begin{gather*}
-A_0\nabla^2\mathbf{W}_0 = \mathbf{0}, \quad \text{in } \Omega, \\
\mathbf{W}_0(\mathbf{x},t) = R\mathbf{V}_B(\mathbf{x},t), \quad \text{on } \Gamma,
\end{gather*}
where $\dot{\mathbf{W}}_0$ satisfies
\begin{gather}
-A_0\nabla^2\dot{\mathbf{W}}_0 = \mathbf{0}, \quad \text{in } \Omega, \label{eq2.36}\\
\dot{\mathbf{W}}_0(\mathbf{x},t) = R\dot{\mathbf{V}}_B(\mathbf{x},t), \quad
 \text{on } \Gamma , \label{eq2.37}
\end{gather}
and the parabolic IBVP, denoted by $\mathcal{P}1.2$:
\begin{gather*}
M_R\dot{\mathbf{W}}-A_0\nabla^2\mathbf{W} 
= -M_R\dot{\mathbf{W}}_0, \quad \text{in }\Omega \times (0,t_f), \\
\mathbf{W}(\mathbf{x},t)=\mathbf{0}, \quad \text{on }\Gamma \times [0,t_f),  \\
\mathbf{W}(\mathbf{x},0)=R\mathbf{V}_I(\mathbf{x})-\mathbf{W}_0(\mathbf{x},0), 
\quad \text{in }\Omega,
\end{gather*}
where $\dot{\mathbf{W}}_0$ is the solution to \eqref{eq2.36}-\eqref{eq2.37}. 
Then the solution to problem $\mathcal{P}1$ is 
\begin{equation}
\bar{\mathbf{W}}=\mathbf{W}_0+\mathbf{W}\,. \label{eq2.41}
\end{equation}
The existence and uniqueness of weak solutions to the subproblems
 $\mathcal{P}1.1$ and $\mathcal{P}1.2$ and continuous dependence of the
 solutions on the initial and boundary data follow from the standard 
elliptic and parabolic theory, yielding the well-posedness in a weak 
sense of the problem $\mathcal{P}1$. Thus, the question of the
 well-posedness of the fully coupled IBVP \eqref{eq1.7}-\eqref{eq1.13} 
amounts to establishing the well-posedness of the problem $\mathcal{P}2$.
\end{remark}

\subsection{Problem $\mathcal{P}1$}
Let us now focus our attention on the properties of a weak solution to 
the problem $\mathcal{P}1$. Henceforth, $\hat{C}>0$ denotes a generic constant
 independent of functions to be estimated. We shall start with two 
lemmas concerning weak solutions to the subproblems $\mathcal{P}1.1$ and
 $\mathcal{P}1.2$. We omit their proofs because they follow the 
standard Galerkin methods for elliptic and parabolic problems, respectively.

\begin{lemma} \label{lemma2.5}
Given $\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)$, under
Assumption \ref{assm2.2} (iii) on the matrix $A_0$, there exists a unique 
weak solution $\mathbf{W}_0\in L^2\big(0,t_f;\mathbb{V}^3\big)$ to 
the problem $\mathcal{P}1.1$ such that 
$\dot{\mathbf{W}}_0\in L^2\big(0,t_f;\mathbb{V}^3\big)$ and the following 
estimates hold
\begin{gather*}
\|\mathbf{W}_0\|_{L^2(0,t_f;\mathbb{V}^3)} 
\leq \hat{C} \|\mathbf{V}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}, \\
\|\dot{\mathbf{W}}_0\|_{L^2(0,t_f;\mathbb{V}^3)} 
\leq \hat{C} \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}\,.
\end{gather*}
Moreover, for each $t\in[0,t_f]$,
\begin{gather*}
\|\mathbf{W}_0(t)\|_{\mathbb{V}^3} 
\leq \hat{C} \|\mathbf{V}_B(t)\|_{H^{1/2}(\Gamma)^3},  \\
\|\dot{\mathbf{W}}_0(t)\|_{\mathbb{V}^3} 
\leq \hat{C} \|\dot{\mathbf{V}}_B(t)\|_{H^{1/2}(\Gamma)^3}\,.
\end{gather*}
\end{lemma}

\begin{lemma} \label{lemma2.6}
Given $\dot{\mathbf{W}}_0\in L^2(0,t_f;\mathbb{V}^3)$, 
$\mathbf{V}_I \in \mathbb{V}^3$, and $\mathbf{W}_0(0)\in \mathbb{V}^3$, 
under Assumption \ref{assm2.2} (iii),  problem $\mathcal{P}1.2$ admits a 
unique weak solution $\mathbf{W}\in L^\infty\big(0,t_f;\mathbb{V}_0^3\big)$ 
with $\dot{\mathbf{W}}\in L^2\big(0,t_f;\mathbb{H}^3\big)$ and the solution 
depends continuously on the data; that is,
\begin{gather*}
\|\mathbf{W}\|_{L^\infty(0,t_f;\mathbb{V}_0^3)} 
\leq \hat{C}\big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
+\|\mathbf{W}_0(0)\|_{\mathbb{V}^3}+\|\dot{\mathbf{W}}_0\|_{L^2(0,t_f;\mathbb{V}^3)}
\big),  \\
\|\dot{\mathbf{W}}\|_{L^2(0,t_f;\mathbb{H}^3)} 
\leq \hat{C}\big( \|\mathbf{V}_I\|_{\mathbb{V}^3}
+\|\mathbf{W}_0(0)\|_{\mathbb{V}^3}+\|\dot{\mathbf{W}}_0\|_{L^2(0,t_f;\mathbb{V}^3)}
\big)\,.
\end{gather*}
\end{lemma}

As a direct consequence of Lemma \ref{lemma2.5}, Lemma \ref{lemma2.6}, 
and \eqref{eq2.41} we have the following results regarding the well-posedness 
and properties of a weak solution to the  problem $\mathcal{P}1$.

\begin{corollary} \label{cor1}
Given the boundary and initial 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)$ and 
$\mathbf{V}_I \in \mathbb{V}^3$, under Assumption \ref{assm2.2} (iii), 
problem $\mathcal{P}1$ admits a unique weak solution
\begin{equation}
\bar{\mathbf{W}}\in L^2(0,t_f;\mathbb{V}^3) \quad \text{with} \quad 
\dot{\bar{\mathbf{W}}}\in L^2(0,t_f;\mathbb{H}^3) \label{eq2.48}
\end{equation}
and the solution depends continuously on the data in the sense that the 
following estimates hold
\begin{gather}
\begin{aligned}
\|\bar{\mathbf{W}}\|_{L^2(0,t_f; \mathbb{V}^3)} 
&\leq \hat{C}\big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
 +\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3}  \\
&\quad + \|\mathbf{V}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)} 
 + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}\big)
\end{aligned}, \label{eq2.49} \\
\|\dot{\bar{\mathbf{W}}}\|_{L^2(0,t_f; \mathbb{H}^3)} 
\leq \hat{C}\big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
+\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3} 
+ \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}\big)\,. \label{eq2.50}
\end{gather}
\end{corollary}

\section{Well-posedness of problem $\mathcal{P}2$ with a generalized source term} 

In this section we address the questions of existence, uniqueness, 
and continuous dependence on data of a weak solution to the implicit problem 
$\mathcal{P}2$ with a generalized source term $\mathbf{F}$. 
Throughout the section we do not impose the relationship \eqref{eq2.32} 
between $\mathbf{F}$, the solution $\bar{\mathbf{W}}$ to  problem $\mathcal{P}1$,
 and the function $\mathbf{u_B}(\boldsymbol{\hat{\sigma}})$.

\subsection{A priori estimates}
We begin by obtaining formal a priori energy estimates for  system 
$\mathcal{P}2$. An important feature of this system is that the elements 
of the matrices $M_R$ and $A_0$ and the coupling vectors $\mathbf{b}_0$ 
and $\mathbf{b}_1$ differ by more than 20 orders of magnitude. 
Therefore, a more refined element-wise approach that allows one to determine 
precise relationships between components of similar magnitude needs to be applied.

We first formally multiply \eqref{eq2.31} by $\dot{\mathbf{V}}^T$, integrate 
over $\Omega$, and split the first integral on the left-hand side into the 
sum of two integrals:
\begin{equation}
\int_\Omega \dot{\mathbf{V}}^T M_R\dot{\mathbf{V}}\mathrm{d}\Omega 
+ \int_\Omega \dot{\mathbf{V}}^T \mathbf{b}_0 \mathbf{b}_1^T R^{-1}\hat{E}
\dot{\mathbf{V}}\mathrm{d}\Omega 
- \int_\Omega \dot{\mathbf{V}}^T A_0\nabla^2 \mathbf{V}\mathrm{d}\Omega 
= \int_\Omega \dot{\mathbf{V}}^T\mathbf{F} \mathrm{d}\Omega\,. \label{eq3.1}
\end{equation}
We shall estimate each term of \eqref{eq3.1} separately. For the first term 
on the left-hand side we have
\begin{align*}
&\Big|\int_\Omega \Big[\dot{\mathbf{V}}^T M_R \dot{\mathbf{V}}
- \big(m_{11}\dot{V}_1^2+m_{22}\dot{V}_2^2+m_{33}\dot{V}_3^2\big)\Big]
 \mathrm{d}\Omega\Big| \\
& \leq \frac{|m_{12}+m_{21}|}{2} 2\Big|\int_\Omega \dot{V}_1\dot{V}_2
\mathrm{d}\Omega \Big| + \frac{|m_{23}+m_{32}|}{2} 2\Big|
\int_\Omega \dot{V}_2\dot{V}_3\mathrm{d}\Omega\Big|\\
&\quad + \frac{|m_{13}+m_{31}|}{2} 2\Big|\int_\Omega \dot{V}_1\dot{V}_3\mathrm{d}
 \Omega \Big| \\
& \leq \frac{|m_{12}+m_{21}|}{2} 
 \Big[\frac{1}{\varepsilon_1}\|\dot{V}_1\|_{\mathbb{H}}^2 
 + \varepsilon_1 \|\dot{V}_2\|_{\mathbb{H}}^2 \Big] 
 + \frac{|m_{23}+m_{32}|}{2} \Big[\frac{1}{\varepsilon_2}\|\dot{V}_2\|_{\mathbb{H}}^2 
 + \varepsilon_2 \|\dot{V}_3\|_{\mathbb{H}}^2 \Big] \\
&\quad + \frac{|m_{13}+m_{31}|}{2} \Big[\frac{1}{\varepsilon_3}
 \|\dot{V}_3\|_{\mathbb{H}}^2 + \varepsilon_3 \|\dot{V}_1\|_{\mathbb{H}}^2 \Big]\,,
\end{align*}
where $\varepsilon_i>0$, $i=1,2,3$, and therefore,
\begin{equation} \label{eq3.2}
\begin{aligned}
\int_\Omega \dot{\mathbf{V}}^T M_R \dot{\mathbf{V}}\mathrm{d}\Omega
& \geq \Big[m_{11} - \frac{|m_{12}+m_{21}|}{2}\frac{1}{\varepsilon_1}
 - \frac{|m_{13}+m_{31}|}{2}\varepsilon_3\Big]\|\dot{V}_1\|_{\mathbb{H}}^2  \\
&\quad + \Big[m_{22} - \frac{|m_{23}+m_{32}|}{2} \frac{1}{\varepsilon_2}
 - \frac{|m_{21}+m_{12}|}{2}\varepsilon_1\Big]\|\dot{V}_2\|_{\mathbb{H}}^2  \\
&\quad + \Big[m_{33} - \frac{|m_{31}+m_{13}|}{2}\frac{1}{\varepsilon_3}
  - \frac{|m_{32}+m_{23}|}{2}\varepsilon_2\Big]\|\dot{V}_3\|_{\mathbb{H}}^2\,.
\end{aligned}
\end{equation}
Using \eqref{eq2.23} and \eqref{eq2.28} together with \eqref{eq1.6}
and the definition \eqref{eq2.7} of the matrix $R$, the second term on
 the left-hand side of \eqref{eq3.1} can be written as
\begin{equation} \label{eq3.3}
\begin{aligned}
&\int_\Omega \dot{\mathbf{V}}^T \mathbf{b}_0 \mathbf{b}_1^T R^{-1}\hat{E}
 \dot{\mathbf{V}}\mathrm{d}\Omega  \\
&= \int_\Omega (\mathbf{b}_0^T \dot{\mathbf{V}}) \nabla
 \cdot E_0(\mathbf{b}_1^T R^{-1}\dot{\mathbf{V}}) \mathrm{d}\Omega  \\
& = a_E\big(E_0(\mathbf{b}_0^T \dot{\mathbf{V}}), E_0(\mathbf{b}_1^T R^{-1}
 \dot{\mathbf{V}})\big) \\
& = a_E\big(E_0(\mathbf{b}_0^T \dot{\mathbf{V}}), E_0((\mathbf{b}_0
 +\mathbf{b}_d)^T R^{-1}\dot{\mathbf{V}})\big)  \\
& = a_E\big(E_0(\mathbf{b}_0^T \dot{\mathbf{V}}),
 E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\big)
 + a_E\big(E_0(\mathbf{b}_0^T \dot{\mathbf{V}}),
 E_0(\mathbf{b}_R^T \dot{\mathbf{V}})\big),
\end{aligned}
\end{equation}
where $\mathbf{b}_R^T = [b_{R1} ~b_{R2} ~b_{R3}]^T = \mathbf{b}_d^T R^{-1}$.
 From \eqref{eq2.22},
\begin{equation}
a_E\big(E_0(\mathbf{b}_0^T \dot{\mathbf{V}}), E_0(\mathbf{b}_0^T
\dot{\mathbf{V}})\big)
\geq \gamma \|E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\|_{0,n}^2  \label{eq3.4}
\end{equation}
and the last term in \eqref{eq3.3} is estimated as follows:
\begin{align*}
&\big|a_E\big(E_0 (\mathbf{b}_0^T \dot{\mathbf{V}}),
 E_0(\mathbf{b}_R^T \dot{\mathbf{V}})\big)\big| \\
&=\big|a_E\big(E_0(\mathbf{b}_R^T \dot{\mathbf{V}}),
  E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\big)\big| \\
& = \big|\int_\Omega (\mathbf{b}_R^T \dot{\mathbf{V}})
 \nabla\cdot E_0(\mathbf{b}_0^T \dot{\mathbf{V}}) \mathrm{d}\Omega\big| \\
&\leq \|\mathbf{b}_R^T \dot{\mathbf{V}}\|_{\mathbb{H}}
 \|\nabla\cdot E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\|_{\mathbb{H}}\\
& \leq\frac{\varepsilon}{2}\|\mathbf{b}_R^T \dot{\mathbf{V}}\|_{\mathbb{H}}^2
 + \frac{n}{2\varepsilon}\|E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\|_{0,n}^2\\
& \leq \frac{3\varepsilon}{2}\big(b_{R1}^2\|\dot{V}_1\|_{\mathbb{H}}^2
 + b_{R2}^2\|\dot{V}_2\|_{\mathbb{H}}^2 + b_{R3}^2\|\dot{V}_3\|_{\mathbb{H}}^2 \big)
 + \frac{n}{2\varepsilon}\|E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\|_{0,n}^2\,.
\end{align*}
Thus,
\begin{equation}
\begin{aligned}
a_E\big(E_0(\mathbf{b}_0^T \dot{\mathbf{V}}),
 E_0(\mathbf{b}_R^T \dot{\mathbf{V}})\big) 
&\geq  - \frac{3\varepsilon}{2}\big(b_{R1}^2\|\dot{V}_1\|_{\mathbb{H}}^2
  + b_{R2}^2\|\dot{V}_2\|_{\mathbb{H}}^2 + b_{R3}^2\|\dot{V}_3\|_{\mathbb{H}}^2 \big)
 \\
 &\quad  - \frac{n}{2\varepsilon}\|E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\|_{0,n}^2\,.
\end{aligned} \label{eq3.5}
\end{equation}
To handle the the last term on the left-hand side of \eqref{eq3.1}, we use
the divergence theorem and the auxiliary boundary condition \eqref{eq2.19} to obtain
\begin{equation}
\begin{aligned}
&- \int_\Omega \dot{\mathbf{V}}^T A_0\nabla^2 \mathbf{V}\mathrm{d}\Omega \\
& = \int_\Omega \Big(\nabla \dot{V}_1 \cdot (a_{11}\nabla V_1 + a_{12}\nabla V_2
  + a_{13}\nabla V_3)  + \nabla \dot{V}_2 \cdot \big(a_{21}\nabla V_1
 + a_{22}\nabla V_2 \\
&\quad + a_{23}\nabla V_3\big) 
 + \nabla \dot{V}_3 \cdot (a_{31}\nabla V_1 + a_{32}\nabla V_2
 + a_{33}\nabla V_3)\Big) \mathrm{d}\Omega   \\
&= \int_\Omega \sum_{i=1}^n\partial_i \dot{\mathbf{V}}^T A_0 \partial_i
 \mathbf{V} \mathrm{d} \Omega\,.
\end{aligned}\label{eq3.6}
\end{equation}
Let us define the bilinear form
$\bar{a}:\mathbb{V}_0^3\times\mathbb{V}_0^3\to \mathbb{R}$ by
\begin{equation}
\bar{a}(\boldsymbol{\psi}, \boldsymbol{\Phi})
= \int_\Omega \sum_{i=1}^n\partial_i \mathbf{\boldsymbol{\psi}}^T A_0
\partial_i \boldsymbol{\Phi} \mathrm{d} \Omega \,. \label{eq3.7}
\end{equation}
Taking into account the symmetry of $A_0$, \eqref{eq3.6} and \eqref{eq3.7} yield
\begin{equation}
- \int_\Omega \dot{\mathbf{V}}^T A_0\nabla^2 \mathbf{V}\mathrm{d}\Omega
= \frac{1}{2} \frac{d}{dt} \bar{a}(\mathbf{V}, \mathbf{V})\,. \label{eq3.8}
\end{equation}

\begin{remark}  \label{a_properties}\rm
The symmetric bilinear form 
$\bar{a}:\mathbb{V}_0^3\times\mathbb{V}_0^3\to \mathbb{R}$ 
is coercive and continuous. Indeed, Assumption \ref{assm2.2} (iii) 
on the matrix $A_0$ guarantees that there exists $\hat{\alpha}>0$ such that
\begin{equation}
\bar{a}(\boldsymbol{\Phi}, \boldsymbol{\Phi}) 
\geq \hat{\alpha}\|\boldsymbol{\Phi}\|_{0,3}^2, \quad \forall 
\boldsymbol{\Phi}\in\mathbb{V}_0^3 \,.\label{eq3.9}
\end{equation}
On the other hand, for every 
$\boldsymbol{\psi}, \boldsymbol{\Phi} \in \mathbb{V}_0^3$,
\begin{align*}
|\bar{a}(\boldsymbol{\psi}, \boldsymbol{\Phi})|
 & = \big| \int_\Omega \sum_{i=1}^n\partial_i \mathbf{\boldsymbol{\psi}}^T A_0 
 \partial_i \boldsymbol{\Phi} \mathrm{d} \Omega \big| \\
&\leq \underset{1\leq j,k\leq 3}{\max} \{|a_{jk}|\} \sum_{i=1}^n 
 \sum_{j,k=1}^3 \|\partial_i \psi_j\|_{\mathbb{H}}
 \|\partial_i \Phi_k\|_{\mathbb{H}} \\
& \leq \underset{1\leq j,k\leq 3}{\max} \{|a_{jk}|\} \sum_{j=1}^3 
\sum_{i=1}^n \|\partial_i \psi_j\|_{\mathbb{H}} \sum_{k=1}^3 
 \sum_{i=1}^n \|\partial_i \Phi_k\|_{\mathbb{H}} \\
& \leq 3n \underset{1\leq j,k\leq 3}{\max} \{|a_{jk}|\} 
 \|\boldsymbol{\psi}\|_{0,3} \|\boldsymbol{\Phi}\|_{0,3}\,,
\end{align*}
and the result follows.
\end{remark}
Applying \eqref{eq3.2}, \eqref{eq3.3}-\eqref{eq3.5}, and \eqref{eq3.8} 
to \eqref{eq3.1} we have
\begin{equation} \label{eq3.10}
\begin{aligned}
& 2\Big[m_{11} - \frac{|m_{12}+m_{21}|}{2}\frac{1}{\varepsilon_1}
- \frac{|m_{13}+m_{31}|}{2} \varepsilon_3
  - \frac{3\varepsilon}{2}b_{R1}^2\Big]\|\dot{V}_1\|_{\mathbb{H}}^2   \\
&+ 2\Big[m_{22} - \frac{|m_{23}+m_{32}|}{2} \frac{1}{\varepsilon_2}
  - \frac{|m_{21}+m_{12}|}{2}\varepsilon_1- \frac{3\varepsilon}{2}b_{R2}^2\Big]
 \|\dot{V}_2\|_{\mathbb{H}}^2   \\
&+  2\Big[m_{33} - \frac{|m_{31}+m_{13}|}{2}\frac{1}{\varepsilon_3}
 - \frac{|m_{32}+m_{23}|}{2}\varepsilon_2
 - \frac{3\varepsilon}{2}b_{R3}^2\Big]\|\dot{V}_3\|_{\mathbb{H}}^2   \\
&+ 2\Big(\gamma-\frac{n}{2\varepsilon}\Big)\|E_0(\mathbf{b}_0^T \dot{\mathbf{V}})
 \|_{0,n}^2 + \frac{d}{dt} \bar{a}(\mathbf{V}, \mathbf{V}) \\
&\leq 2 \int_\Omega \dot{\mathbf{V}}^T\mathbf{F} \mathrm{d}\Omega\,.
\end{aligned}
\end{equation}
The right-hand side of \eqref{eq3.10} is majorized by
\begin{equation}
\begin{aligned}
2 \int_\Omega \dot{\mathbf{V}}^T\mathbf{F} \mathrm{d}\Omega
& \leq m_{11}\|\dot{V}_1\|_{\mathbb{H}}^2 +\frac{1}{m_{11}}\|F_1\|_{\mathbb{H}}^2
 + m_{22}\|\dot{V}_2\|_{\mathbb{H}}^2 +\frac{1}{m_{22}}\|F_2\|_{\mathbb{H}}^2  \\
&\quad + m_{33}\|\dot{V}_3\|_{\mathbb{H}}^2 +\frac{1}{m_{33}}
\|F_3\|_{\mathbb{H}}^2\,.
\end{aligned}  \label{eq3.11}
\end{equation}
From \eqref{eq3.10} and \eqref{eq3.11} we obtain
\begin{align*}
& \Big[m_{11} - |m_{12}+m_{21}|\frac{1}{\varepsilon_1}
 - |m_{13}+m_{31}| \varepsilon_3
 - 3\varepsilon b_{R1}^2\Big]\|\dot{V}_1\|_{\mathbb{H}}^2 \\
&+ \Big[m_{22} - |m_{23}+m_{32}| \frac{1}{\varepsilon_2}
 - |m_{21}+m_{12}| \varepsilon_1- 3\varepsilon b_{R2}^2\Big]
 \|\dot{V}_2\|_{\mathbb{H}}^2 \\
&+ \Big[m_{33} - |m_{31}+m_{13}|\frac{1}{\varepsilon_3}
 - |m_{32}+m_{23}|\varepsilon_2 - 3\varepsilon b_{R3}^2\Big]
 \|\dot{V}_3\|_{\mathbb{H}}^2 \\
&+ 2\Big(\gamma-\frac{n}{2\varepsilon}\Big)\|E_0(\mathbf{b}_0^T
 \dot{\mathbf{V}})\|_{0,n}^2 +  \frac{d}{dt} \bar{a}(\mathbf{V}, \mathbf{V}) \\
&\leq \frac{1}{m_{11}}\|F_1\|_{\mathbb{H}}^2
 + \frac{1}{m_{22}}\|F_2\|_{\mathbb{H}}^2 
+ \frac{1}{m_{33}}\|F_3\|_{\mathbb{H}}^2\,.
\end{align*}
Integrating the last inequality with respect to time and applying
\eqref{eq2.14} and \eqref{eq3.9} yield
\begin{equation}
\hat{\beta}\int_0^t\|\dot{\mathbf{V}}\|_{\mathbb{H}^3}^2 \mathrm{d}s
 + \hat{\gamma}\int_0^t\|E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\|_{0,n}^2 \mathrm{d}s
+ \hat{\alpha}\|\mathbf{V}(t)\|_{0,3}^2
\leq \frac{1}{\underline {m}}\int_0^t\|\mathbf{F}\|_{\mathbb{H}^3}^2
\mathrm{d}s \label{eq3.12}
\end{equation}
for every $t\in[0,t_f]$, where
$\hat{\beta}=\underset{1\leq i \leq 3}{\min}~\beta_i$ with
\begin{gather}
\beta_i=m_{ii} - |m_{ij}+m_{ji}|\frac{1}{\varepsilon_i}
 - |m_{ik}+m_{ki}|\varepsilon_k - 3\varepsilon b_{Ri}^2\,, \label{eq3.13} \\
i=1,2,3, \; j=(i\bmod{3})+1, ~k=(j\bmod{3})+1 \,,\nonumber \\
\hat{\gamma}=2\Big(\gamma-\frac{n}{2\varepsilon}\Big) \label{eq3.14}
\end{gather}
and $\underline{m}=\min_{1\leq i \leq 3}  m_{ii} >0$. The values
$\varepsilon_i>0$, $i=1,2,3$, and $\varepsilon >0$ must be chosen to satisfy
$\hat{\beta}>0$ and $\hat{\gamma} \geq 0$ and we have the following result.

\begin{lemma} \label{lemma3.2}
Under Assumption \ref{assm2.2} (i)-(ii), the a priori energy estimate 
\eqref{eq3.12} holds with $\hat{\beta}>0$ and $\hat{\gamma} \geq 0$
 if $\varepsilon_i$, $i=1,2,3$, and $\varepsilon$ satisfy the following conditions.
\begin{gather}
\frac{p_i|m_{ij}+m_{ji}|}{m_{ii}} \leq \varepsilon_{i} 
\leq \frac{m_{jj}}{q_j|m_{ij}+m_{ji}|}, \quad i=1,2,3, \; 
j=(i\bmod{3})+1\,, \label{eq3.15}\\
\frac{n}{2\gamma} \leq \varepsilon \leq \underset{1\leq i \leq 3}{\mathrm{min}}
\frac{m_{ii}}{3r_ib_{Ri}^2}\,, \label{eq3.16}
\end{gather}
where $p_i>0$, $q_i>0$, $r_i>0$, $\frac{1}{p_i}+\frac{1}{q_i}+\frac{1}{r_i}<1$, 
$i=1,2,3$, $\gamma=\gamma(\Omega)$ is given by \eqref{eq2.22}, and 
$n$ is the dimension of the problem.
\end{lemma}

\begin{proof}
The condition $\hat{\beta}>0$ is equivalent to $\beta_i>0, ~i=1,2,3$. 
From \eqref{eq3.13} we observe that the latter condition holds true if 
we assume that there exist $p_i>0$, $q_i>0$, and $r_i>0$, 
$\frac{1}{p_i}+\frac{1}{q_i}+\frac{1}{r_i}<1$, $i=1,2,3$, such that
\begin{equation}
\begin{gathered}
 |m_{ij}+m_{ji}|\frac{1}{\varepsilon_i} \leq \frac{1}{p_i}m_{ii}\,,\\
 |m_{ik}+m_{ki}|\varepsilon_k \leq \frac{1}{q_i}m_{ii}\,, \\
3\varepsilon b_{Ri}^2 \leq \frac{1}{r_i}m_{ii},
\end{gathered} \label{eq3.17}
\end{equation}
where $i=1,2,3$, $j=(i\bmod{3})+1$, $k=(j\bmod{3})+1$.
The first two inequalities in \eqref{eq3.17} give the following conditions
 on $\varepsilon_1$:
\begin{equation*}
|m_{12}+m_{21}|\frac{1}{\varepsilon_1} 
\leq \frac{1}{p_1}m_{11}, \quad |m_{12}+m_{21}|\varepsilon_1 
\leq \frac{1}{q_2}m_{22}
\end{equation*}
which implies
\begin{equation}
\frac{p_1|m_{12}+m_{21}|}{m_{11}} 
\leq \varepsilon_1 \leq \frac{m_{22}}{q_2|m_{12}+m_{21}|}. \label{eq3.18}
\end{equation}
Similarly, we obtain
\begin{equation}
\frac{p_2|m_{23}+m_{32}|}{m_{22}} 
\leq \varepsilon_2 \leq \frac{m_{33}}{q_3|m_{23}+m_{32}|} \label{eq3.19}
\end{equation}
and
\begin{equation}
\frac{p_3|m_{13}+m_{31}|}{m_{33}} \leq \varepsilon_3 
\leq \frac{m_{11}}{q_1|m_{13}+m_{31}|}. \label{eq3.20}
\end{equation}
The existence of positive intervals 
\[
\big[ \frac{p_i|m_{ij}+m_{ji}|}{m_{ii}}, 
\frac{m_{jj}}{q_j|m_{ij}+m_{ji}|}\big], \quad i=1,2,3,\;
j=(i\bmod{3})+1,
\]
 is guaranteed by Assumption \ref{assm2.2} (ii); hence, 
 \eqref{eq3.18}-\eqref{eq3.20} yield \eqref{eq3.15}.

The condition \eqref{eq3.16} on $\varepsilon$ follows immediately 
from \eqref{eq3.14} and the last inequality in \eqref{eq3.17}.
 Assumption \ref{assm2.2} (i) together with \eqref{eq2.4}, \eqref{eq2.5}, 
\eqref{eq2.22}, and the results in Horgan \cite{Horgan} concerning 
the Korn's constant guarantee the existence of a positive interval 
for $\varepsilon$.
\end{proof}

\begin{remark}\rm \label{solution_spaces}
From formal a priori energy estimate \eqref{eq3.12} we conclude that a weak solution to problem $\mathcal{P}1$ is expected to be
\begin{equation*}
\mathbf{V} \in L^{\infty}(0, t_f; \mathbb{V}_0^3), \quad 
\dot{\mathbf{V}}\in L^2(0,t_f;\mathbb{H}^3)
\end{equation*}
provided $\mathbf{F}\in L^2(0,t_f;\mathbb{H}^3)$.
\end{remark}

\subsection{Abstract formulation}
The preceding remark suggests the weak formulation of the problem $\mathcal{P}2$ 
as follows:
Given $\mathbf{F}\in L^2(0,t_f;\mathbb{H}^3)$, find
 $\mathbf{V} \in L^{\infty}(0, t_f; \mathbb{V}_0^3)$ with 
$\dot{\mathbf{V}}\in L^2(0,t_f;\mathbb{H}^3)$ such that, for all 
$\boldsymbol{\psi}\in\mathbb{V}_0^3$,
\begin{gather}
\begin{aligned}
&\int_\Omega \boldsymbol{\psi}^T M_R\dot{\mathbf{V}}\mathrm{d}\Omega 
+ \int_\Omega \boldsymbol{\psi}^T \mathbf{b}_0 \mathbf{b}_1^T R^{-1}\hat{E}
\dot{\mathbf{V}}\mathrm{d}\Omega 
- \int_\Omega \boldsymbol{\psi}^T A_0\nabla^2 \mathbf{V}\mathrm{d}\Omega \\
&= \int_\Omega \boldsymbol{\psi}^T\mathbf{F} \mathrm{d}\Omega,
\end{aligned} \label{eq3.21} \\
\mathbf{V}(\mathbf{x},0)=\mathbf{0}\,. \label{eq3.22}
\end{gather}
Let us define the continuous bilinear form 
$m:\mathbb{H}^3\times\mathbb{H}^3\to\mathbb{R}$ by
\begin{equation*} %\label{eq3.23}
m(\boldsymbol{\psi}, \boldsymbol{\Phi}) 
= \int_\Omega \boldsymbol{\psi}^T M_R \boldsymbol{\Phi}\mathrm{d}\Omega\,.
\end{equation*}
Then the first integral on the left-had side of \eqref{eq3.21} takes the form
\begin{equation}
\int_\Omega \boldsymbol{\psi}^T M_R\dot{\mathbf{V}}\mathrm{d}\Omega 
= m(\boldsymbol{\psi}, \dot{\mathbf{V}})\,. \label{eq3.24}
\end{equation}
Applying \eqref{eq2.23} and \eqref{eq2.28} to the second term on the 
left-hand side of \eqref{eq3.21} gives
\begin{equation}
\int_\Omega \boldsymbol{\psi}^T \mathbf{b}_0 \mathbf{b}_1^T R^{-1}
\hat{E}\dot{\mathbf{V}}\mathrm{d}\Omega 
= a_E(E_0(\mathbf{b}_0^T\boldsymbol{\psi}), E_0(\mathbf{b}_1^TR^{-1}
\dot{\mathbf{V}}))\,. \label{eq3.25}
\end{equation}
Furthermore, from the definition of the linear operator $E_0$ \eqref{eq2.23} 
and the continuity of $a_E(\cdot,\cdot)$ and $E_0$, we observe that
\begin{equation*}
\boldsymbol{\psi}, \boldsymbol{\Phi} 
\mapsto a_E(E_0(\mathbf{b}_0^T\boldsymbol{\psi}), 
E_0(\mathbf{b}_1^TR^{-1}\boldsymbol{\Phi}))
\end{equation*}
is a bilinear continuous map from 
$\mathbb{H}^3\times\mathbb{H}^3$ to $\mathbb{R}$.
Thus, we can define a continuous bilinear form 
$\bar{l}:\mathbb{H}^3\times\mathbb{H}^3 \to \mathbb{R}$ by
\begin{equation}
\bar{l}(\boldsymbol{\psi}, \boldsymbol{\Phi}) 
= m(\boldsymbol{\psi}, \boldsymbol{\Phi}) 
+ a_E(E_0(\mathbf{b}_0^T\boldsymbol{\psi}), 
E_0(\mathbf{b}_1^TR^{-1}\boldsymbol{\Phi}))\,. \label{eq3.26}
\end{equation}
Combining \eqref{eq3.24}-\eqref{eq3.26} leads to
\begin{equation}
\int_\Omega \boldsymbol{\psi}^T M_R\dot{\mathbf{V}}\mathrm{d}\Omega 
+ \int_\Omega \boldsymbol{\psi}^T \mathbf{b}_0 \mathbf{b}_1^T R^{-1}\hat{E}
\dot{\mathbf{V}}\mathrm{d}\Omega = \bar{l}(\boldsymbol{\psi},
\dot{\mathbf{V}})\,. \label{eq3.27}
\end{equation}

\begin{lemma} \label{lemma_3.4}
The bilinear form $\bar{l}:\mathbb{H}^3\times\mathbb{H}^3 \to \mathbb{R}$ 
is coercive and hence non-degenerate.
\end{lemma}

\begin{proof}
From \eqref{eq3.2}-\eqref{eq3.5}, and \eqref{eq3.27}, for every 
$\dot{\mathbf{V}} \in \mathbb{H}^3$, we have
\begin{align*}
\bar{l}(\dot{\mathbf{V}}, \dot{\mathbf{V}}) 
& \geq \Big[m_{11} - \frac{|m_{12}+m_{21}|}{2}\frac{1}{\varepsilon_1} 
- \frac{|m_{13}+m_{31}|}{2} \varepsilon_3 
- \frac{3\varepsilon}{2}b_{R1}^2\Big]\|\dot{V}_1\|_{\mathbb{H}}^2 \\
&\quad + \Big[m_{22} - \frac{|m_{23}+m_{32}|}{2} \frac{1}{\varepsilon_2} 
 - \frac{|m_{21}+m_{12}|}{2}\varepsilon_1
 - \frac{3\varepsilon}{2}b_{R2}^2\Big]\|\dot{V}_2\|_{\mathbb{H}}^2 \\
&\quad  + \Big[m_{33} - \frac{|m_{31}+m_{13}|}{2}\frac{1}{\varepsilon_3}
 - \frac{|m_{32}+m_{23}|}{2}\varepsilon_2
 - \frac{3\varepsilon}{2}b_{R3}^2\Big]\|\dot{V}_3\|_{\mathbb{H}}^2 \\
&\quad + \Big(\gamma-\frac{n}{2\varepsilon}\Big)\|E_0(\mathbf{b}_0^T 
 \dot{\mathbf{V}})\|_{0,n}^2 \\
& \geq \frac{1}{2} \sum_{i=1}^3 (m_{ii}+\beta_i)\|\dot{V}_i\|_{\mathbb{H}}^2
 + \frac{1}{2}\hat{\gamma}\|E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\|_{0,n}^2\,,
\end{align*}
where $\beta_i >0$, $i=1,2,3$, and $\hat{\gamma} \geq 0$ are given by 
\eqref{eq3.13}-\eqref{eq3.16}. It follows that
\begin{equation*}
\bar{l}(\dot{\mathbf{V}}, \dot{\mathbf{V}}) 
\geq \frac{1}{2}\underline{m}\|\dot{\mathbf{V}}\|_{\mathbb{H}^3}^2
\end{equation*}
where $\underline{m}=\min_{1\leq i \leq 3} m_{ii}>0$, which completes the 
proof.
\end{proof}

To handle the last term on the left-hand side of \eqref{eq3.21}, we use 
the same argument as in deriving formal a priori estimates
 (see  \eqref{eq3.6} and \eqref{eq3.7} above) and obtain
\begin{equation}
- \int_\Omega \boldsymbol{\psi}^T A_0\nabla^2 \mathbf{V}\mathrm{d}\Omega 
= \bar{a}(\boldsymbol{\psi},\mathbf{V})\,. \label{eq3.28}
\end{equation}
Substituting \eqref{eq3.27} and \eqref{eq3.28} into \eqref{eq3.21} yields 
the following abstract formulation of  problem $\mathcal{P}2$ equivalent 
to \eqref{eq3.21}-\eqref{eq3.22}:
Given $\mathbf{F}\in L^2(0,t_f;\mathbb{H}^3)$, find 
$\mathbf{V} \in L^{\infty}(0, t_f; \mathbb{V}_0^3)$ such that 
$\dot{\mathbf{V}}\in L^2(0,t_f;\mathbb{H}^3)$ and, for all 
$\boldsymbol{\psi}\in\mathbb{V}_0^3$,
\begin{gather}
\bar{l}(\boldsymbol{\psi},\dot{\mathbf{V}}) + \bar{a}(\boldsymbol{\psi},\mathbf{V}) 
= (\boldsymbol{\psi}, \mathbf{F})_{\mathbb{H}^3}\,, \label{eq3.29} \\
\mathbf{V}(\mathbf{x},0)=\mathbf{0}\,. \label{eq3.30}
\end{gather}

\begin{remark}  \label{alternative_form} \rm
The bilinear forms $\bar{l}:\mathbb{H}^3\times\mathbb{H}^3 \to \mathbb{R}$ and 
$\bar{a}:\mathbb{V}_0^3\times\mathbb{V}_0^3\to \mathbb{R}$ are continuous and
 coercive and besides, $\bar{a}(\cdot, \cdot)$ is symmetric. 
Therefore, we can associate with $\bar{l}(\cdot, \cdot)$ and $\bar{a}(\cdot, \cdot)$ 
the linear bijective operators $\mathcal{L}:\mathbb{H}^3 \to \mathbb{H}^3$ and
 $\mathcal{A}:\mathbb{V}_0^3\to H^{-1}(\Omega)^3$, respectively, by setting
\begin{gather*}
\bar{l}(\boldsymbol{\psi},\boldsymbol{\Phi})
=(\boldsymbol{\psi}, \mathcal{L}\boldsymbol{\Phi})_{\mathbb{H}^3}, \quad 
\boldsymbol{\psi},\boldsymbol{\Phi} \in \mathbb{H}^3, \\
\bar{a}(\boldsymbol{\psi},\boldsymbol{\Phi})
=\langle \boldsymbol{\psi},\mathcal{A}\boldsymbol{\Phi}\rangle, \quad 
\boldsymbol{\psi},\boldsymbol{\Phi} \in \mathbb{V}_0^3\,,
\end{gather*}
where $\langle \cdot, \cdot \rangle$ is the duality pairing between 
$\mathbb{V}_0^3$ and $H^{-1}(\Omega)^3$ and 
$\langle \boldsymbol{\psi}, \boldsymbol{\phi}\rangle
 = (\boldsymbol{\psi},\boldsymbol{\phi})_{\mathbb{H}^3}$ for every 
$\boldsymbol{\psi}\in \mathbb{V}_0^3$ and $\boldsymbol{\phi}\in\mathbb{H}^3$. 
Then \eqref{eq3.29} can be written as
\begin{equation*}
\langle \boldsymbol{\psi},\mathcal{L}\dot{\mathbf{V}} \rangle 
+ \langle \boldsymbol{\psi},\mathcal{A}\mathbf{V} \rangle 
= \langle \boldsymbol{\psi},\mathbf{F} \rangle, \quad 
\forall \boldsymbol{\psi}\in \mathbb{V}_0^3
\end{equation*}
and we obtain an alternative formulation of the problem \eqref{eq3.29}-\eqref{eq3.30}:
Given $\mathbf{F}\in L^2(0,t_f;\mathbb{H}^3)$, find 
$\mathbf{V} \in L^{\infty}(0,t_f; \mathbb{V}_0^3)$ such that 
$\dot{\mathbf{V}}\in L^2(0,t_f;\mathbb{H}^3)$ and
\begin{gather*}
\mathcal{L}\dot{\mathbf{V}} + \mathcal{A}\mathbf{V} = \mathbf{F}\,, \\
\mathbf{V}(\mathbf{x},0)=\mathbf{0}\,.
\end{gather*}
\end{remark}

\subsection{Well-posedness in a weak sense}
The main result of this section is formulated in the next theorem.

\begin{theorem} \label{thrm3.6}
Given $\mathbf{F}\in L^2(0,t_f;\mathbb{H}^3)$, under Assumption \ref{assm2.2}, 
there exists a unique weak solution
\begin{equation}
\mathbf{V} \in L^{\infty}(0,t_f; \mathbb{V}_0^3) \quad \text{with} \quad 
\dot{\mathbf{V}}\in L^2(0,t_f;\mathbb{H}^3) \label{eq3.31}
\end{equation}
in the sense of \eqref{eq3.29}-\eqref{eq3.30} to the problem $\mathcal{P}2$ and 
 the solution depends continuously on the data $\mathbf{F}$; that is, the mapping
\begin{equation*}
\mathbf{F} \mapsto \mathbf{V}, \dot{\mathbf{V}}
\end{equation*}
from $L^2(0,t_f;\mathbb{H}^3)$ to 
$L^{\infty}(0,t_f; \mathbb{V}_0^3) \times L^2(0,t_f;\mathbb{H}^3)$ is continuous.
\end{theorem}

\begin{remark} \rm
The proof of the existence and uniqueness of a weak solution to  problem 
$\mathcal{P}2$ can be completed by the standard Galerkin method. We omit it
 and focus on continuous dependence results. The key to the proof of 
existence and uniqueness are energy estimates arising from \eqref{eq3.12}, 
the continuity and coercivity of the bilinear forms $\bar{l}$ and $\bar{a}$, 
as well as the symmetry of $\bar{a}$.
\end{remark}

\begin{proof}
The continuous dependence on the data $\mathbf{F}$ follows from the a priori 
energy estimate \eqref{eq3.12}:  for every $t\in [0,t_f]$,
\begin{equation*}
\hat{\beta}\int_0^t\|\dot{\mathbf{V}}\|_{\mathbb{H}^3}^2 \mathrm{d}s 
+ \hat{\gamma}\int_0^t\|E_0(\mathbf{b}_0^T \dot{\mathbf{V}})\|_{0,n}^2 \mathrm{d}s 
+ \hat{\alpha}\|\mathbf{V}(t)\|_{0,3}^2 
\leq \frac{1}{\underline {m}}\int_0^t\|\mathbf{F}\|_{\mathbb{H}^3}^2 
\mathrm{d}s \quad{\eqref{eq3.12}},
\end{equation*}
where $\hat{\beta}$, $\hat{\alpha}$ and $\underline{m}$ are positive constants 
and $\hat{\gamma}\geq 0$. The above inequality implies
\begin{equation*}
\|\mathbf{V}(t)\|_{0,3}^2 \leq \frac{1}{\hat{\alpha}\underline{m}} 
\int_0^{t_f} \|\mathbf{F}\|_{\mathbb{H}^3}^2  \mathrm{d}s
\end{equation*}
for every $t\in [0,t_f]$, and thus,
\begin{equation}
\|\mathbf{V}\|_{L^\infty(0,t_f;\mathbb{V}_0^3)} 
\leq \hat{C} \|\mathbf{F}\|_{L^2(0,t_f;\mathbb{H}^3)}\,. \label{eq3.32}
\end{equation}
Taking $t=t_f$ in \eqref{eq3.12} also yields
\begin{equation}
\|\dot{\mathbf{V}}\|_{L^2(0,t_f;\mathbb{H}^3)} 
\leq \hat{C} \|\mathbf{F}\|_{L^2(0,t_f;\mathbb{H}^3)} \label{eq3.33}
\end{equation}
which completes the proof.
\end{proof}

\section{Main results}

The main results of the paper are summarized in the following theorem 
that establishes the well-posedness in a weak sense of the fully coupled 
parabolic-elliptic IBVP \eqref{eq1.7}-\eqref{eq1.13}.

\begin{theorem} \label{thrm_main}
Given the initial data $\mathbf{V}_I \in \mathbb{V}^3$ and the 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)$ 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)$,
 under Assumption \ref{assm2.2} on the matrices of diffusion coefficients, 
the IBVP \eqref{eq1.7}-\eqref{eq1.13} for the fully coupled TCpu system 
\eqref{eq1.1}-\eqref{eq1.4} admits a unique weak solution
\begin{equation}
(\bar{\mathbf{V}}, \mathbf{u}) \in L^2(0,t_f;\mathbb{V}^3) 
\times L^2(0,t_f;\tilde{\mathbb{V}}_0^n) \label{eq4.1}
\end{equation}
with
\begin{equation}
(\dot{\bar{\mathbf{V}}}, \dot{\mathbf{u}}) \in L^2(0,t_f;\mathbb{H}^3) 
\times L^2(0,t_f;\tilde{\mathbb{V}}_0^n) \label{eq4.2}
\end{equation}
and this solution depends continuously on the data $\mathbf{V}_I$, 
$\mathbf{V}_B(0)$, $\mathbf{V}_B$, $\dot{\mathbf{V}_B}$, 
$\hat{\boldsymbol{\sigma}}$, and $\dot{\hat{\boldsymbol{\sigma}}}$.
\end{theorem}

\begin{proof}
(i) \emph{Existence and uniqueness.} From the results of Section 2, we know 
that the solution $(\bar{\mathbf{V}}, \mathbf{u})$ to 
 problem \eqref{eq1.7}-\eqref{eq1.13} is given by 
\eqref{eq2.33} and \eqref{eq2.26A}:
\begin{gather*}
\bar{\mathbf{V}} = R^{-1}(\mathbf{V}+\bar{\mathbf{W}}) ,\\
\mathbf{u}=E_0(\mathbf{b}_1^T R^{-1}\mathbf{V})+ E_0(\mathbf{b}_1^T
 R^{-1}\bar{\mathbf{W}})+\mathbf{u}_B ,
\end{gather*}
where $\mathbf{V}$ and $\bar{\mathbf{W}}$ are the solutions to the problems
 $\mathcal{P}2$ and $\mathcal{P}1$, respectively, and $\mathbf{u}_B$ is 
defined by \eqref{eq2.24}. The well-posedness of a weak solution to the 
problem $\mathcal{P}1$ is stated in Corollary \ref{cor1}, and 
Theorem \ref{thrm3.6} provides the existence, uniqueness and continuous 
dependence on data of a weak solution to the problem $\mathcal{P}2$ with a 
generalized source term $\mathbf{F}$.

The next step toward the proof of Theorem \ref{thrm_main} is to show that 
the source term, \eqref{eq2.32}, 
\begin{equation*}
\mathbf{F}=-\mathbf{b}_0 \mathbf{b}_1^T R^{-1}\hat{E}\dot{\bar{\mathbf{W}}}
-\mathbf{b}_0 (\nabla\cdot\dot{\mathbf{u}}_B) 
\end{equation*}
corresponding to the IBVP \eqref{eq1.7}-\eqref{eq1.13} satisfies the assumption 
$\mathbf{F}\in L^2(0,t_f;\mathbb{H}^3)$ of Theorem \ref{thrm3.6}. 
The definition \eqref{eq2.28} of the operator $\hat{E}$ and Corollary 
\ref{cor1} yield
\begin{equation}
\mathbf{b}_0 \mathbf{b}_1^T R^{-1}\hat{E}
\dot{\bar{\mathbf{W}}}\in L^2(0,t_f;\mathbb{H}^3)\,. \label{eq4.3}
\end{equation}
Assumption \ref{assm2.1} and  \eqref{eq2.24} imply
\begin{gather}
\mathbf{u}_B\in L^2(0,t_f;\tilde{\mathbb{V}}_0^n)\,, \label{eq4.4} \\
\dot{\mathbf{u}}_B\in L^2(0,t_f;\tilde{\mathbb{V}}_0^n)\,, \label{eq4.5}
\end{gather}
and thus,
\begin{equation}
\mathbf{b}_0 (\nabla\cdot\dot{\mathbf{u}}_B) \in L^2(0,t_f;\mathbb{H}^3)\,. 
\label{eq4.6}
\end{equation}
From \eqref{eq2.32}, \eqref{eq4.3}, and \eqref{eq4.6} we conclude that 
$\mathbf{F}\in L^2(0,t_f;\mathbb{H}^3)$; hence, Theorem \ref{thrm3.6} 
holds for the problem $\mathcal{P}2$ associated with the 
IBVP \eqref{eq1.7}-\eqref{eq1.13}. The existence and uniqueness of a 
weak solution $(\bar{\mathbf{V}}, \mathbf{u})$ satisfying \eqref{eq4.1} 
and \eqref{eq4.2} now follow immediately from \eqref{eq2.33}, 
\eqref{eq2.48}, \eqref{eq3.31}, \eqref{eq4.4}, \eqref{eq4.5}, 
and the definition \eqref{eq2.23} of the operator $E_0$ applied to 
 \eqref{eq2.26A} and its time derivative.

(ii) \emph{Continuous dependence on data.} From \eqref{eq2.33}, \eqref{eq2.49},
 \eqref{eq3.32}, and \eqref{eq4.1} we have
\begin{equation} \label{eq4.7}
\begin{aligned}
\|\bar{\mathbf{V}}\|_{L^2(0,t_f;\mathbb{V}^3)}
&\leq \hat{C} \big(\|\mathbf{V}\|_{L^\infty (0,t_f;\mathbb{V}_0^3)}
+ \|\bar{\mathbf{W}}\|_{L^2 (0,t_f;\mathbb{V}^3)}\big) \\
& \leq \hat{C} \big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
+\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3}
+ \|\mathbf{V}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}  \\
& \quad + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}
+\|\mathbf{F}\|_{L^2(0,t_f;\mathbb{H}^3)}\big)\,.
\end{aligned}
\end{equation}
Similarly, the time derivative of \eqref{eq2.33} together with
\eqref{eq2.50}, \eqref{eq3.33}, and \eqref{eq4.2} give
\begin{equation} \label{eq4.8}
\begin{aligned}
\|\dot{\bar{\mathbf{V}}}\|_{L^2(0,t_f;\mathbb{H}^3)}
&\leq \hat{C} \big(\|\dot{\mathbf{V}}\|_{L^\infty (0,t_f;\mathbb{H}^3)}
+ \|\dot{\bar{\mathbf{W}}}\|_{L^2 (0,t_f;\mathbb{H}^3)}\big) \\
& \leq \hat{C} \big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
 +\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3}  \\
& \quad + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}
+\|\mathbf{F}\|_{L^2(0,t_f;\mathbb{H}^3)}\big)\,.
\end{aligned}
\end{equation}
Next we estimate the source term on the right-hand side of \eqref{eq4.7}
and \eqref{eq4.8}. Applying the continuity property \eqref{eq2.29} of
$\hat{E}$ and \eqref{eq2.50} to \eqref{eq2.32} gives
\begin{equation} \label{eq4.9}
\begin{aligned}
\|\mathbf{F}\|_{L^2(0,t_f;\mathbb{H}^3)}
& \leq \hat{C} \big(\|\dot{\bar{\mathbf{W}}}\|_{L^2(0,t_f;\mathbb{H}^3)}
+ \|\nabla \cdot \dot{\mathbf{u}}_B\|_{L^2(0,t_f;L^2(\Omega))}\big)   \\
& \leq \hat{C} \big(\|\mathbf{V}_I\|_{\mathbb{V}^3}+\|\mathbf{V}_B(0)\|_{H^{1/2}
 (\Gamma)^3}   \\
& \quad + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}
+ \|\nabla \cdot \dot{\mathbf{u}}_B\|_{L^2(0,t_f;\mathbb{H})}\big)\,.
\end{aligned}
\end{equation}
Now we evaluate the function $\mathbf{u}_B$ and its time derivative.
Taking $\boldsymbol{\Phi}=\mathbf{u}_B$ in \eqref{eq2.24} and using
\eqref{eq2.22} and Assumption \ref{assm2.1}, we have, for every $t\in[0,t_f]$,
\begin{equation} \label{eq4.10}
\begin{aligned}
\gamma \|\mathbf{u}_B\|_{0,n}^2
& \leq \big|a_E(\mathbf{u}_B,\mathbf{u}_B)\big|
= \Big|\int_{\Gamma} (\boldsymbol{\hat{\sigma}}\mathbf{n})
 \cdot\mathbf{u}_B \mathrm{d} \Gamma\Big|  \\
& \leq \|\boldsymbol{\hat{\sigma}}\mathbf{n}\|_{L^2(\Gamma)^n}
 \|\mathbf{u}_B\|_{L^2(\Gamma)^n}
 \leq \hat{C} \|\boldsymbol{\hat{\sigma}}\|_{L^2(\Gamma)^{n\times n}}
\|\mathbf{u}_B\|_{0,n}\,.
\end{aligned}
\end{equation}
The above inequality follows from the trace theorem and Poincar\'e-Friedrichs'
inequality. Equation \eqref{eq4.10} implies, for every $t\in[0,t_f]$,
\begin{equation}
\|\mathbf{u}_B(t)\|_{0,n}
\leq \hat{C} \|\boldsymbol{\hat{\sigma}}(t)\|_{L^2(\Gamma)^{n\times n}}\,.
\label{eq4.11}
\end{equation}
The time derivative of \eqref{eq2.24} and Assumption \ref{assm2.1}
 yield $\dot{\mathbf{u}}_B\in\tilde{\mathbb{V}}_0^n$.
Therefore, differentiating \eqref{eq2.24} with respect to time and taking
$\boldsymbol{\Phi}=\dot{\mathbf{u}}_B$, the above argument also gives
\begin{equation}
\|\dot{\mathbf{u}}_B(t)\|_{0,n}
\leq \hat{C} \|\dot{\hat{\boldsymbol{\sigma}}}(t)\|_{L^2(\Gamma)^{n\times n}}
\label{eq4.11A}
\end{equation}
for every $t\in[0,t_f]$, and consequently,
\begin{equation}
\|\nabla \cdot \dot{\mathbf{u}}_B(t)\|_{\mathbb{H}}
\leq \hat{C} \|\dot{\hat{\boldsymbol{\sigma}}}(t)\|_{L^2(\Gamma)^{n\times n}}\,.
\label{eq4.12}
\end{equation}
Integrating \eqref{eq4.11}-\eqref{eq4.12} over $[0,t_f]$, we get
\begin{gather}
\|\mathbf{u}_B\|_{L^2(0,t_f;\tilde{\mathbb{V}}_0^n)}
\leq \hat{C} \|\boldsymbol{\hat{\sigma}}\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}\,,
\label{eq4.13} \\
\|\dot{\mathbf{u}}_B\|_{L^2(0,t_f;\tilde{\mathbb{V}}_0^n)}
\leq \hat{C} \|\dot{\hat{\boldsymbol{\sigma}}}
\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}\,,
\label{eq4.15} \\
\|\nabla \cdot \dot{\mathbf{u}}_B\|_{L^2(0,t_f;\mathbb{H})}
\leq \hat{C} \|\dot{\hat{\boldsymbol{\sigma}}}
\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}\,. \label{eq4.16}
\end{gather}
Substituting \eqref{eq4.16} into \eqref{eq4.9} gives
\begin{equation} \label{eq4.17}
\begin{aligned}
\|\mathbf{F}\|_{L^2(0,t_f;\mathbb{H}^3)}
& \leq \hat{C} \big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
 +\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3}   \\
& \quad + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}
 + \|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}\big)\,.
\end{aligned}
\end{equation}
Applying \eqref{eq4.17} to \eqref{eq4.7} and \eqref{eq4.8}, we obtain respectively
\begin{equation}
\begin{aligned}
\|\bar{\mathbf{V}}\|_{L^2(0,t_f;\mathbb{V}^3)}
& \leq \hat{C} \big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
 +\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3}
 + \|\mathbf{V}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)} \\
& \quad + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}
 + \|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}\big)
\end{aligned} \label{eq4.18}
\end{equation}
and
\begin{equation}
\begin{aligned}
\|\dot{\bar{\mathbf{V}}}\|_{L^2(0,t_f;\mathbb{H}^3)}
& \leq \hat{C} \big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
 +\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3}  \\
& \quad + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}
 + \|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}\big)\,.
\end{aligned} \label{eq4.19}
\end{equation}
Next we show that the $\mathbf{u}$-component of the solution depends
 continuously on data. Applying \eqref{eq2.8} and the continuity of
$E_0$ \eqref{eq2.25} to \eqref{eq2.26} and combining the result with \eqref{eq4.13},
 and \eqref{eq4.18} yields
\begin{align*}
\|\mathbf{u}\|_{L^2(0,t_f;\tilde{\mathbb{V}}_0^n)}
& \leq \hat{C}\big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
 +\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3}  \\
& \quad + \|\mathbf{V}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}
 + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}   \\
& \quad + \|\hat{\boldsymbol{\sigma}}\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}
 + \|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}\big)\,.
\end{align*}
Finally, differentiating \eqref{eq2.26} with respect to time and using
\eqref{eq2.25}, \eqref{eq4.15}, and \eqref{eq4.19}, we obtain
\begin{align*}
\|\dot{\mathbf{u}}\|_{L^2(0,t_f;\tilde{\mathbb{V}}_0^n)}
& \leq \hat{C}\big(\|\mathbf{V}_I\|_{\mathbb{V}^3}
 +\|\mathbf{V}_B(0)\|_{H^{1/2}(\Gamma)^3} \nonumber \\
& \quad + \|\dot{\mathbf{V}}_B\|_{L^2(0,t_f; H^{1/2}(\Gamma)^3)}
 + \|\dot{\hat{\boldsymbol{\sigma}}}\|_{L^2(0,t_f;L^2(\Gamma)^{n\times n})}\big)
\end{align*}
which completes the proof.
\end{proof}

\section{Numerical Example}

In this section we provide the results of numerical experiments for the 
TCpu system with real data illustrating the applicability of the 
proposed well-posedness theory. The calculations were performed in MATLAB R2015b.

We focus on the two-dimensional case and consider the IBVP 
\eqref{eq1.7}-\eqref{eq1.13} for the TCpu system defined in an annular 
region $\Omega \subset \mathbb{R}^2$ with the inner (borehole) radius $R_B$ 
and the outer (far-field) radius $R_F$. The following assumptions are made 
on the region $\Omega$.

\begin{assumption} \label{assm5.1} \rm
\begin{itemize}
\item[(i)] $\delta=\frac{R_B}{R_F} \leq e^{-1}$
\item[(ii)] The region $\Omega$ does not experience pure rotation.
\end{itemize}
\end{assumption}

To eliminate pure rotation, we impose the following 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{eq5.1}
\end{equation}
Using the results obtained by Horgan \cite{Horgan}, it was shown in
 \cite{Malysheva} that under Assumption \ref{assm5.1} and the side 
condition \eqref{eq5.1} the Korn constant $\gamma$ appearing in 
\eqref{eq2.22} is given by
\begin{equation}
\gamma=\frac{G}{2}\Big[1-\Big(\frac{3}{\delta^{-2}+1+\delta^2}\Big)^{1/2}\Big]\,, 
\label{eq5.2}
\end{equation}
where $G$ is the shear modulus; hence, for $0< \delta \leq e^{-1}$,
\begin{equation}
\frac{G}{2}\Big[1-\Big(\frac{3}{e^2+1+e^{-2}}\Big)^{1/2}\Big] 
\leq \gamma < \frac{G}{2}\,.  \label{eq5.3}
\end{equation}


\begin{example} \rm
In this example, the values of the TCpu model parameters describing the
 physical properties of a rock/fluid system with various values of the
 permeability/fluid viscosity ratio, $k/\eta$, and the solute reflected 
fraction, $\mathfrak{R}$, are taken from \cite{Diek1} and represented 
in Table \ref{table:data1} in the appendix of the paper. With these data, 
direct calculations show that parts (i) and (ii) of Assumption \ref{assm2.2} 
are satisfied as follows regardless of the values of $k/\eta$ 
and $\mathfrak{R}$ because the matrix $M_R$ and the vector $\mathbf{b}_R$ 
are independent of those parameters:
\begin{gather*}
0 < m_{11} = 7.1188\mathrm{e}{+03} < 2.7196\mathrm{e}{+11} = b_{R1}^2 \\
0< m_{22} = 7.0722\mathrm{e}{-10} < 2.9166\mathrm{e}{-03} = b_{R2}^2 \\
0 < m_{33} = 4.9082\mathrm{e}{-11} <  1.2659\mathrm{e}{-02} = b_{R3}^2
\end{gather*}
and
\begin{gather*}
0 < (m_{12}+m_{21})^2 = 1.6800\mathrm{e}{-12} 
< 5.0346\mathrm{e}{-06}=m_{11}m_{22} \,,\\
0 < (m_{23}+m_{32})^2 = 6.9032\mathrm{e}{-25} 
< 3.4712\mathrm{e}{-20}=m_{22}m_{33} \,,\\
0 < (m_{13}+m_{31})^2 = 6.1784\mathrm{e}{-09} 
< 3.4941\mathrm{e}{-07}=m_{11}m_{33}\,.
\end{gather*}

From Table \ref{example1_eigenvalues} we observe that the eigenvalues 
$\{\mu_i\}$ of the matrix $\frac{1}{2}\big(M_R+M_R^T\big)$ and the 
eigenvalues $\{\lambda_i\}$ of the matrix $A_0$ are distinct and positive,
 which immediately implies the fulfillment of Assumption \ref{assm2.2} (iii).


Next we determine the values $\varepsilon_i$, $i=1,2,3$, and
 $\varepsilon$ that satisfy conditions \eqref{eq3.15} and \eqref{eq3.16} 
and consequently, guarantee the energy estimates \eqref{eq3.12} and 
the well-posedness of the IBVP \eqref{eq1.7}-\eqref{eq1.13}.
Calculation results for different values of the ratio of radii $\delta$ 
and numbers $p_i$, $q_i$, and $r_i$, $i=1,2,3$, are summarized in 
Table \ref{example1_epsilon} and are the same for all given values 
of $k/\eta$ and $\mathfrak{R}$.

The results of numerical experiments confirm the applicability of 
the proposed well-posdness theory to the TCpu model with real data parameters.
\end{example}

\begin{table}[ht]
\caption{Eigenvalues of the matrices $\frac{1}{2}\big(M_R+M_R^T\big)$ and $A_0$}
\renewcommand{\arraystretch}{1.3}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
\multirow{2}{*}{Eigenvalue} &
	$k/\eta = 10^{-16}$ & $k/\eta = 10^{-16}$ & $k/\eta = 10^{-17}$ 
&$k/\eta = 10^{-17}$\\
& $\mathfrak{R} = 0.01$ & $\mathfrak{R} = 0.15$ & $\mathfrak{R} = 0.01$ 
& $\mathfrak{R} = 0.15$\\
\hline
$\mu_1$ & $7.1188\mathrm{e}{+03}$ & $7.1188\mathrm{e}{+03}$ 
& $7.1188\mathrm{e}{+03}$ & $7.1188\mathrm{e}{+03}$ \\
\hline
$\mu_2$ & $7.0681\mathrm{e}{-10}$ & $7.0681\mathrm{e}{-10}$ 
& $7.0681\mathrm{e}{-10}$ & $7.0681\mathrm{e}{-10}$ \\
\hline
$\mu_3$ & $4.8865\mathrm{e}{-11}$ & $4.8865\mathrm{e}{-11}$ 
& $4.8865\mathrm{e}{-11}$ & $4.8865\mathrm{e}{-11}$ \\
\hline
$\lambda_1$ & $9.8003\mathrm{e}{-03}$ & $9.8003\mathrm{e}{-03}$ 
& $9.8003\mathrm{e}{-03}$ & $9.8003\mathrm{e}{-03}$ \\
\hline
$\lambda_2$ & $1.0000\mathrm{e}{-16}$ & $1.0231\mathrm{e}{-16}$ 
& $9.9910\mathrm{e}{-18}$ & $1.0374\mathrm{e}{-17}$ \\
\hline
$\lambda_3$ & $4.8252\mathrm{e}{-18}$ & $3.2175\mathrm{e}{-18}$ 
& $4.6728\mathrm{e}{-18}$ & $4.2899\mathrm{e}{-18}$ \\
\hline
\end{tabular}
\end{center}
\label{example1_eigenvalues}
\end{table}

\begin{table}[b]
\caption{Values of $\varepsilon_i$, $i=1,2,3$, and $\varepsilon$ }
\renewcommand{\arraystretch}{1.3}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\multirow{6}{*}
{$\begin{array}{c}
	p_i=\frac{1}{4}\\
	q_i=\frac{1}{4}\\
	r_i=\frac{1}{4}
	\end{array}
	$} & & $4.5518\mathrm{e}{-11} \leq \varepsilon_1 \leq 4.6292\mathrm{e}{+05}$ \\
	& & $1.3848\mathrm{e}{-12} \leq \varepsilon_2 \leq 2.3630\mathrm{e}{+02}$ \\
	& & $4.0036\mathrm{e}{+05} \leq \varepsilon_3 \leq 3.6227\mathrm{e}{+08}$ \\
	\cline{2-3}
	& $\delta=0.01$ & $1.1488\mathrm{e}{-10} \leq \varepsilon \leq 5.1698\mathrm{e}{-09}$\\
	\cline{2-3}
	& $\delta=0.1$ & $1.3639\mathrm{e}{-10} \leq \varepsilon \leq 5.1698\mathrm{e}{-09}$\\
	\cline{2-3}
	& $\delta=e^{-1}$ & $2.7753\mathrm{e}{-10} \leq \varepsilon \leq 5.1698\mathrm{e}{-09}$\\
\hhline{|=|=|=|}
%
\multirow{5}{*}
{$\begin{array}{c}
	p_i=\frac{1}{8}\\
	q_i=\frac{1}{3}\\
	r_i=\frac{1}{2}
	\end{array}
	$} & & $2.2759\mathrm{e}{-11} \leq \varepsilon_1 \leq 3.4719\mathrm{e}{+05}$ \\
	& & $6.9238\mathrm{e}{-13} \leq \varepsilon_2 \leq 1.7722\mathrm{e}{+02}$ \\
	& & $2.0018\mathrm{e}{+05} \leq \varepsilon_3 \leq 2.7170\mathrm{e}{+08}$ \\
	\cline{2-3}
	& $\delta=0.01$ & $1.1488\mathrm{e}{-10} \leq \varepsilon \leq 2.5849\mathrm{e}{-09}$\\
	\cline{2-3}
	& $\delta=0.1$ & $1.3639\mathrm{e}{-10} \leq \varepsilon \leq 2.5849\mathrm{e}{-09}$\\
	\cline{2-3}
	& $\delta=e^{-1}$ & $2.7753\mathrm{e}{-10} \leq \varepsilon \leq 2.5849\mathrm{e}{-09}$\\
\hhline{|=|=|=|}
%
\multirow{5}{*}
{$\begin{array}{c}
	p_i=\frac{1}{10}\\
	q_i=\frac{1}{100}\\
	r_i=\frac{1}{1000}
	\end{array}
	$} & & $1.8207\mathrm{e}{-11} \leq \varepsilon_1 \leq 1.1573\mathrm{e}{+07}$ \\
	& & $5.5390\mathrm{e}{-13} \leq \varepsilon_2 \leq 5.9075\mathrm{e}{+03}$ \\
	& & $1.6014\mathrm{e}{+05} \leq \varepsilon_3 \leq 9.0567\mathrm{e}{+09}$ \\
	\cline{2-3}
	& $\delta=0.01$ & $1.1488\mathrm{e}{-10} \leq \varepsilon \leq 1.2924\mathrm{e}{-06}$\\
	\cline{2-3}
	& $\delta=0.1$ & $1.3639\mathrm{e}{-10} \leq \varepsilon \leq 1.2924\mathrm{e}{-06}$\\
	\cline{2-3}
	& $\delta=e^{-1}$ & $2.7753\mathrm{e}{-10} \leq \varepsilon \leq 1.2924\mathrm{e}{-06}$\\
	\hline
\end{tabular}
\end{center}
\label{example1_epsilon}
\end{table}

\section{Appendix: Physical parameters of the coupled TCpu system}

The elements of the matrices of diffusion coefficients $M$ and $A$ \eqref{eq1.5}, 
the coupling vectors $\mathbf{b}_0$ and $\mathbf{b}_1$ \eqref{eq1.6}, 
the bulk modulus $K$, and the shear modulus $G$ are defined as follows:

\begin{gather*}
\Lambda = \frac{c}{T_F}-s_F\beta_s\tilde{\omega}, \quad
c=(1-\phi)\rho_sc_s+\phi\rho_fc_f, \quad
s_F=c_f-\frac{\beta_fp_F}{\rho_f}, \quad
\tilde{\omega}=\frac{\omega}{1-C_F}, \\
\Sigma=\xi\beta_s, \quad
\xi=\frac{RT_F\tilde{\omega}(1-2C_F)}{M^sC_F}, \quad
\Phi=\Xi-\phi \beta_f+\frac{\beta_s\tilde{\omega}}{\rho_f}, \quad
\Xi=(\alpha-\phi)\beta_s, \\
\alpha=1-\frac{K}{K_s}, \quad K=\frac{E}{3(1-2\nu)}, \quad
\Gamma = \Xi-\phi\beta_f-\frac{s_F\tilde{\omega}}{K_s}, \quad
\chi=-\frac{\xi}{K_s}, \quad \psi=\frac{\alpha-\phi}{K_s},\\
\Psi= \psi +\frac{\phi}{K_f}-\frac{\tilde{\omega}}{\rho_fK_s},  \quad
\tilde{\Omega} = \frac{RT_F}{M^s(1-C_F)C_F}, \quad
\zeta=K\beta_s, \quad G= \frac{E}{2(1+\nu)}\,.
\end{gather*}
Other physical properties of the rock/fluid system are given in 
Table \ref{table:data1}.

\begin{table}[ht]
\caption{Physical properties of the rock/fluid system (data taken from \cite{Diek1})}
\renewcommand{\arraystretch}{1.3}
\begin{center}
\begin{tabular}{|l|c|}
\hline
\textbf{Property} & \textbf{Value} \\
\hline
Drained elastic modulus, $~E$ & $45\times 10^9$ Pa \\
\hline
Solid bulk modulus, $~K_s$ & $65\times 10^9$ Pa \\
\hline
Drained Poisson's ratio, $~\nu$ & 0.27 \\
\hline
Permeability/Fluid viscosity ratio, $k/\eta$ & $10^{-16}$, 
$10^{-17}$ m$^2$/(Pa$\cdot$s) \\
\hline
Reference porosity, $\phi$ & 0.15 \\
\hline
Fluid density, $\rho_f$ & 1111 kg/m$^3$ \\
\hline
Fluid specific heat capacity, $~c_f$ & 4186 J/(Kg$^{\circ}$K) \\
\hline
Fluid bulk modulus, $~K_f$ & $3.3 \times 10^9$ Pa \\
\hline
Fluid volumetric thermal expansion coefficient, $~\beta_f$ 
& $3\times 10^{-4 ~\circ}$K$^{-1}$ \\
\hline
Solid density, $~\rho_s$ & $2.83 \times 10^3$ kg/m$^3$ \\
\hline
Solid specific heat capacity, $~c_s$ & 920 J/(Kg$^{\circ}$K) \\
\hline
Solid volumetric thermal expansion coefficient, $~\beta_s$
 & $2.4 \times 10^{-5 ~\circ}$K$^{-1}$ \\
\hline
Rock thermal conductivity coefficient, $~k^T$ & 4 W/(m$^{\circ}$K) \\
\hline
Thermal osmosis coefficient, $~K^T$ & $10^{-11}$ m$^2$/(s$^{\circ}$K) \\
\hline
Solute molar mass, $~M^s$ & 0.1111 kg/m$^3$ \\
\hline
Solute reflected fraction, $~\mathfrak{R}$ & 0.01, 0.15 \\
\hline
Solute chemical diffusion coefficient, $~D$ & $10^{-9}$ m$^2$/s \\
\hline
Solute thermal diffusion coefficient, $~D^T$ & $10^{-10}$ m$^2$/(s$^{\circ}$K) \\
\hline
Chemical stress coupling parameter, $~\omega$ & 100 kg/m$^3$ \\
\hline
Universal gas constant, $~R$ & 8.3144598 JK$^{-1}$/mol$^{-1}$ \\
\hline
In situ formation pore pressure, $~p_F$ & 52 MPa \\
\hline
Reference formation solute mass fraction, $~C_F$ & 0.1 \\
\hline
In situ formation temperature, $~T_F$ & 135$^\circ$C \\
\hline
\end{tabular}
\end{center}
\label{table:data1}
\end{table}


\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{Barucq1} H. Barucq, M. Madaune-Tort, P. Saint-Macary;
\emph{Some existence-uniqueness results for a class of one-dimensional nonlinear 
Biot models}, Nonlinear Analysis, 61 (2005), 591-612.

\bibitem{Barucq2} H. Barucq, M. Madaune-Tort, P. Saint-Macary;
\emph{On nonlinear Biot’s consolidation models},
 Nonlinear Analysis, 63 (2005), 985-995.

\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{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{Duvaut_Lions} G. Duvaut, J. L. Lions;
\emph{Inequalities in Mechanics and Physics}, Springer, Berlin, 1976.

\bibitem{Ekbote} S. Ekbote, Y. Abousleiman;
\emph{Porochemothermoelastic solution for an inclined borehole in transversely
 isotropic formation}, J. Eng. Mech., 131 (2005), 522-533.

\bibitem{Ghassemi1} 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{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{Malysheva} T.	Malysheva, L. W. White;
\emph{Sufficient condition for Hadamard well-posedness of a coupled 
thermo-chemo-poroelastic system}, Electron. J. Differential Equations, Vol. 2016, 
No. 15 (2016), 1-17.

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

\bibitem{Showalter1} R. E. Showalter;
\emph{Diffusion in deformable media, IMA Volumes in Mathematics and its 
Applications}, 131 (2000), 115-130.
%
\bibitem{Showalter2} R. E. Showalter;
\emph{Diffusion in poro-elastic media}, J. Math. Anal. Appl., 251 (2000), 310-340.

\end{thebibliography}

\end{document}