Efficient Iterative Decoupling Methods for Thermo-Poroelasticity Based on a Four-Field Formulation

Cai Mingchao Morgan State University Li Jingzhi South University of Science and Technology Li Ziliang ²²footnotemark: 2 Liu Qiang³³footnotemark: 3 Shenzhen University

Abstract

This paper studies the thermo-poroelasticity model. By introducing an intermediate variable, we transform the original three-field model into a four-field model. Building upon this four-field model, we present both a coupled finite element method and a decoupled iterative finite element method. We prove the stability and optimal convergence of the coupled finite element method. Furthermore, we establish the convergence of the decoupled iterative method. This paper focuses primarily on analyzing the iterative decoupled algorithm. It demonstrates that the algorithm’s convergence does not require any additional assumptions about physical parameters or stabilization parameters. Numerical results are provided to demonstrate the effectiveness and theoretical validity of these new methods.

Keywords thermo-poroelasticity, decoupled algorithms, finite element method.

1 Introduction

The thermo-poroelasticity model provides a comprehensive framework for describing the interplay between mechanical deformation, fluid flow, and heat transfer in porous media. By integrating Biot’s poroelasticity theory with the heat equation, this model captures the effects of temperature variations on stress and displacement, accounting for the compressibility and thermal expansion of both fluid and solid components. The foundational work in poroelasticity was pioneered by Terzaghi and later expanded by Biot [22, 3, 4], who developed a coupled model to describe solid-fluid and fluid-solid interactions. Extending these concepts to non-isothermal scenarios, the thermo-poroelasticity model comprises a force balance equation, a mass conservation equation, and a heat conduction equation, effectively representing the complex interactions among these physical processes. This model’s broad applicability across various engineering and environmental disciplines has driven significant research interest. Applications span various areas such as soil mechanics, pharmaceutical tablet analysis, nuclear waste management [24], geothermal energy systems [23], biomedical engineering, and carbon dioxide sequestration [21].

In recent years, there has been significant progress in computational research on the thermo-poroelasticity model. Studies such as [9] have demonstrated the existence and uniqueness of solutions for nonlinear thermo-poroelasticity, along with deriving a priori energy estimates and regularity properties. To overcome computational challenges, the classical three-field formulation has been expanded to a four-field model by introducing the divergence of displacement as an auxiliary variable [10], which has also been extended to nonlinear thermo-poroelasticity incorporating convective transport [13]. Various numerical methods have been proposed for these models, including Galerkin methods and combined Galerkin-mixed finite element methods for nonlinear formulations [27, 25]. Additional approaches include a hybrid finite element method utilizing mixed finite elements for pressure, characteristic finite elements for temperature, and Galerkin finite elements for displacement [26], as well as discontinuous Galerkin methods and related algorithms for fully coupled nonlinear problems [5, 1].

Among recent numerical algorithms, iteratively decoupled algorithms have been extensively explored. For instance, unconditionally stable sequential methods for all-ways coupled thermo-poromechanics were proposed in [17]. Splitting schemes for decoupling poroelasticity and thermoelasticity [18], and monolithic, partially decoupled, and fully decoupled schemes based on a five-field formulation incorporating heat and Darcy fluxes [8] have also been investigated. Furthermore, a decoupled iterative solution approach combined with reduced-order modeling has been developed to improve computational efficiency [2]. While convergence proofs for coupled and decoupled iterative schemes exist [18, 8], analysis of the convergence order of these methods remains limited.

This work focuses on developing an unconditionally stable and optimally convergent algorithm rooted in a novel four-field formulation of the thermo-poroelasticity model. To mitigate the issue of Poisson locking, we introduce an auxiliary variable $\xi=-\lambda\nabla\cdot\bm{u}+\alpha p+\beta T$ , transforming the original three-field model into a robust four-field framework. Building upon this reformulation, we design and rigorously analyze an iteratively decoupled method tailored for linear thermo-poroelastic problems. Theoretical analysis establishes the method’s unconditional stability and optimal convergence. Unlike earlier approaches, such as time-extrapolation-based decoupling algorithms for Biot’s model [16], which often face stability limitations and reduced accuracy due to explicit extrapolation, our method overcomes these drawbacks. By iteratively solving each subproblem while incorporating time extrapolation, we ensure not only stability but also maintain optimal convergence rates. This approach provides a more reliable and accurate computational framework for addressing coupled thermo-poroelastic processes.

The structure of the paper is as follows: In Section 2, we provide an overview of the thermo-poroelastic problem. We introduce a four-field model and give its a priori estimates. In Section 3, we present two algorithms: the coupled algorithm and the decoupled algorithm. In Section 4, we analyze the convergence of the decoupled algorithm. Furthermore, we derive the convergence rate for the discrete coupled problem relative to the continuous coupled problem. Finally, Section 5 showcases the numerical results.

2 Mathematical model and its reformulation

We introduce some notations that will be utilized throughout the subsequent text. Let $\Omega\in\mathbb{R}^{d},~{}d=2,3$ be a bounded domain and $J=(0,\tau)$ with $\tau$ being the final time. For $1\leq s<\infty$ , let $L^{s}(\Omega)=\{v:\int_{\Omega}|v|^{s}<\infty\}$ represent the standard Banach space, with the associated norm $\|v\|_{L^{s}(\Omega)}=(\int_{\Omega}|v|^{s})^{\frac{1}{s}}$ . Usually, without special emphasis, $\|\cdot\|_{L^{2}(\Omega)}$ can alternatively be substituted with $\|\cdot\|_{L^{2}(\Omega)}$ . We use $(\cdot,\cdot)$ to denote the $L^{2}$ -inner product. For $s=\infty$ , let $L^{\infty}(\Omega)=\{v:\Omega\to\mathbb{R}:\text{ess sup}_{x\in\Omega}|v(x)|<\infty\}$ be the space of uniformly bounded measurable functions on $\Omega$ , with norm $\|v\|_{L^{\infty}(\Omega)}=\text{ess sup}_{x\in\Omega}|v|$ . Let $W_{s}^{m}(\Omega)=\{v:D^{\alpha}v\in L^{s}(\Omega),0\leq|\alpha|\leq m,\|v\|\leq\infty\}$ denote the standard Sobolev space with the associated norm $\|v\|_{W_{s}^{k}(\Omega)}=(\sum_{|\alpha|\leq k}\|D^{\alpha}v\|^{s}_{L^{s}(\Omega)})^{1/s}$ . In particular, define $H^{m}(\Omega):=W_{2}^{m}(\Omega)$ as the special Sobolev space with $p=2$ and define $H_{0}^{m}(\Omega)=\{v\in H^{1}(\Omega):v|_{\partial\Omega}=0\}$ . For a Banach space $X$ , we give $L^{s}(J;X)=\{v:(\int_{J}\|v\|_{X}^{s})^{\frac{1}{s}}<\infty\}$ . For $v\in L^{s}(J;X)$ , define the associated norm $\|v\|_{L^{s}(J;X)}=(\int_{J}\|v\|_{X}^{s})^{\frac{1}{s}}$ . For $v\in L^{2}(J,L^{2}(\Omega))$ , it is a function in space $L^{2}(\Omega)$ square-integrable in time. Moreover, we define the space $L^{\infty}(J;L^{2}(\Omega))=\{v(t):\text{ess sup}_{t\in J}\|v(t)\|_{L^{2}(\Omega)}<\infty\}$ . We also use $\bm{v}_{t}$ and $v_{t}$ to denote the partial derivative of the function $\bm{v}$ in $[H^{1}(\Omega)]^{d}$ and $\bm{v}$ in $H^{1}(\Omega)$ with respect to $t$ , respectively.

2.1 A linear thermo-poroelasticity model

We consider the following linear thermo-poroelasticity problem [7], where the goal is to determine the displacement $\bm{u}$ , fluid pressure $p$ , and temperature $T$ . The governing equations are:

$\displaystyle-\nabla\cdot(2\mu\bm{\varepsilon}(\bm{u})+\lambda\nabla\cdot\bm{u}\bm{I})+\alpha\nabla p+\beta\nabla T$	$\displaystyle=\bm{f},\quad\text{in }\Omega\times J,$	(2.1)
$\displaystyle\frac{\partial}{\partial t}(c_{0}p-b_{0}T+\alpha\nabla\cdot\bm{u})-\nabla\cdot(\bm{K}\nabla p)$	$\displaystyle=g,\quad\text{in }\Omega\times J,$
$\displaystyle\frac{\partial}{\partial t}(a_{0}T-b_{0}p+\beta\nabla\cdot\bm{u})-\nabla\cdot(\bm{\Theta}\nabla T)$	$\displaystyle=H_{s},\quad\text{in }\Omega\times J,$

subject to the initial conditions:

\bm{u}(\cdot,0)=\bm{u}^{0},\quad p(\cdot,0)=p^{0},\quad T(\cdot,0)=T^{0},

(2.2)

and the homogeneous Dirichlet boundary conditions:

\bm{u}=0,\quad p=0,\quad T=0,\quad\text{on }\partial\Omega\times J.

(2.3)

The operator $\bm{\varepsilon}$ is defined as $\bm{\varepsilon}(\bm{u})=\frac{1}{2}(\nabla\bm{u}+\nabla\bm{u}^{T})$ and the matrix $\bm{I}$ is the identity tensor. The right-hand side functions $\bm{f},g$ , and $H_{s}$ denote the body force, the mass source, and the heat source respectively. $\alpha$ is the Biot-Willis constant and $\beta$ is the thermal stress coefficient. $a_{0},b_{0}$ , and $c_{0}$ represent the effective thermal capacity, the thermal dilation coefficient, and the specific storage coefficient, respectively. The matrix parameter $\bm{K}=(K_{ij})^{d}_{ij=1}$ denotes the permeability divided by fluid viscosity and matrix parameter $\bm{\Theta}=(\Theta_{ij})^{d}_{ij=1}$ denotes the effective thermal conductivity. $\mu$ and $\lambda$ represent the Lamé parameters, which can be formulated in terms of Young’s modulus $E$ and Poisson’s ratio $\nu$ :

\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)},\qquad\mu=\frac{E}{2(1+\nu)}.

Note that in this paper considered coefficients and variables in the model are without units which is also known as dimensionless form.

We now introduce the four-field formulation for the above model (2.1). More clearly, following the methodology for handling Biot’s problem in [19, 20], we define an intermediate auxiliary variable to represent the volumetric contribution to the total stress:

\xi=-\lambda\nabla\cdot\bm{u}+\alpha p+\beta T,

which is commonly referred to as the pseudo-total pressure [1]. Substituting this variable into equation (2.1) and using parameter transformations above, the four-field thermo-poroelasticity problem can be expressed as:

$\displaystyle-2\mu\nabla\cdot\bm{\varepsilon}(\bm{u})+\nabla\xi$	$\displaystyle=\bm{f},$	(2.4)
$\displaystyle-\lambda\nabla\cdot\bm{u}-\xi+\alpha p+\beta T$	$\displaystyle=0,$
$\displaystyle-\frac{\alpha}{\lambda}\frac{\partial}{\partial t}\xi+(c_{0}+\frac{\alpha^{2}}{\lambda})\frac{\partial}{\partial t}p-\nabla\cdot(\bm{K}\nabla p)+(\frac{\alpha\beta}{\lambda}-b_{0})\frac{\partial}{\partial t}T$	$\displaystyle=g,$
$\displaystyle-\frac{\beta}{\lambda}\frac{\partial}{\partial t}\xi+(\frac{\alpha\beta}{\lambda}-b_{0})\frac{\partial}{\partial t}p+(a_{0}+\frac{\beta^{2}}{\lambda})\frac{\partial}{\partial t}T-\nabla\cdot(\bm{\Theta}\nabla T)$	$\displaystyle=H_{s}.$

We still use the initial condition (2.2) and the Dirichlet boundary conditions (2.3) for (2.4). In practical situations, nonhomogeneous Dirichlet and Neumann boundary conditions are commonly encountered. The analysis performed for homogeneous boundary conditions can be straightforwardly extended to accommodate these cases.

As outlined in [8], we make the following Assumptions for the model parameters throughout this paper (similar assumptions are also discussed in [10, 27]):

(A1) $\bm{K}$ and $\bm{\Theta}$ are constant in time, symmetric, definite, and positive. There exist $k_{m}>0$ and $k_{M}>0$ such that

k_{m}|\bm{v}|^{2}\leq\bm{v}^{T}\bm{K}(x)\bm{v},\quad|\bm{K}(x)\bm{v}|\leq k_{M}|\bm{v}|,\quad\forall\bm{v}\in\mathbb{R}^{d}\backslash\{\bm{0}\},

where $|\cdot|$ denotes the $l^{2}$ -norm of a vector in $\mathbb{R}^{d}$ . There exist $\theta_{m}>0$ and $\theta_{M}>0$ such that

\theta_{m}|\bm{v}|^{2}\leq\bm{v}^{T}\bm{\Theta}(x)\bm{v},\quad|\bm{\Theta}(x)\bm{v}|\leq\theta_{M}|\bm{v}|,\quad\forall\bm{v}\in\mathbb{R}^{d}\backslash\{\bm{0}\}.

(A2) The constants $\lambda$ , $\mu$ , $\alpha$ and $\beta$ are strictly positive constants.

(A3) The constants $a_{0},b_{0}$ and $c_{0}$ are such that $a_{0},~{}c_{0}\geq b_{0}\geq 0$ .

In addition to Assumptions (A1)-(A3), the following Assumptions will also be maintained throughout the remainder of this paper.

We assume the initial conditions satisfy $\bm{u}^{0}\in[H^{1}(\Omega)]^{d}$ and $p^{0},T^{0}\in L^{2}(\Omega)$ , while the source terms satisfy $\bm{f}\in[L^{2}(\Omega)]^{d}$ , $g\in L^{2}(\Omega)$ , and $H_{s}\in L^{2}(\Omega)$ . For simplicity, we further assume that $\bm{f}$ , $g$ , and $H_{s}$ are time-independent.

Define the following spaces

\bm{V}=\bm{H}^{1}_{0}(\Omega),~{}Q=L^{2}(\Omega),~{}W=H_{0}^{1}(\Omega).

and their dual spaces are denoted as $\bm{V}^{\prime}$ , $Q^{\prime}$ and $W^{\prime}$ . Given $\tau>0$ , a 4-tuple $(\bm{u},\xi,p,T)$ with

\bm{u}\in L^{\infty}(J;\bm{V}),~{}\xi\in L^{\infty}(J;Q),~{}p,~{}T\in L^{\infty}(J;W),

\nabla u_{t},~{}\xi_{t},~{}p_{t},~{}T_{t}\in L^{2}(J;W^{\prime}).

is called a weak solution of problem (2.4) if for almost every $t\in(0,\tau]$ there holds,

$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}),\bm{\varepsilon}(\bm{v}))-(\nabla\cdot\bm{v},\xi)$	$\displaystyle=(\bm{f},\bm{v}),~{}\forall\bm{v}\in\bm{V},$	(2.5)
$\displaystyle-(\nabla\cdot\bm{u},\phi)-\frac{1}{\lambda}(\xi,\phi)+\frac{\alpha}{\lambda}(p,\phi)+\frac{\beta}{\lambda}(T,\phi)$	$\displaystyle=0,~{}\forall\phi\in Q,$
$\displaystyle-\frac{\alpha}{\lambda}(\xi_{t},q)+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{t},q)+\bm{(}\bm{K}\nabla p,\nabla q)+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{t},q)$	$\displaystyle=(g,q),~{}\forall q\in W,$
$\displaystyle-\frac{\beta}{\lambda}(\xi_{t},S)+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{t},S)+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{t},S)+(\bm{\Theta}\nabla T,\nabla S)$	$\displaystyle=(H_{s},S),~{}\forall S\in W.$

2.2 A priori estimates of the weak solution

In this subsection, we provide an energy estimate and a priori estimates. Using these estimates, we establish the existence and uniqueness theorem in Theorem 2.2. Firstly, we present Korn’s inequality([6]), which elucidates the relationship between $\|\bm{u}\|_{H^{1}(\Omega)}$ and $\|\bm{\varepsilon}(\bm{u})\|_{L^{2}(\Omega)}$ .

Lemma 2.1.

There exists a constant $C$ such that the Korn’s inequality

C\|\bm{u}\|_{H^{1}(\Omega)}^{2}\leq\|\bm{\varepsilon}(\bm{u})\|_{L^{2}(\Omega)}^{2},\quad\forall\bm{u}\in[H^{1}_{0}(\Omega)]^{d},~{}d=2,3

holds true.

Then, we give an energy estimate.

Lemma 2.2.

Assume that the Assumptions of (A1)-(A3) hold. Every weak solution $(\bm{u},\xi,p,T)$ of problem (2.5) satisfies the following energy law:

\displaystyle E(t)+\int_{0}^{t}\bm{K}(\nabla p,\nabla p)ds+\int_{0}^{t}\bm{\Theta}(\nabla T,\nabla T)ds-\int_{0}^{t}(g,p)ds-\int_{0}^{t}(H_{s},T)ds=E(0),

(2.6)

for all $t\in J$ , where

	$\displaystyle E(t):=$	$\displaystyle\mu\\|\bm{\varepsilon}(\bm{u})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p+\beta T-\xi\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p\\|_{L^{2}(\Omega)}^{2}$		(2.7)
		$\displaystyle+\frac{a_{0}-b_{0}}{2}\\|T\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\\|p-T\\|_{L^{2}(\Omega)}^{2}-(\bm{f},\bm{u}).$		(2.7)

Furthermore,

\displaystyle\|\xi(t)\|_{L^{2}(\Omega)}\leq C\big{(}2\mu\|\bm{\varepsilon}(u)\|_{L^{2}(\Omega)}+\|\bm{f}\|_{L^{2}(\Omega)}\big{)},

(2.8)

where $C$ is a constant depending only on $\Omega$ .

Proof.

Differentiating the second equation of system (2.5) with respect to $t$ , taking $\bm{v}=\bm{u}_{t},\phi=-\xi,q=p$ and $S=T$ in (2.5), we have

$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}),\bm{\varepsilon}(\bm{u}_{t}))-(\nabla\cdot\bm{u}_{t},\xi)$	$\displaystyle=(\bm{f},\bm{u}_{t}),$	(2.9)
$\displaystyle(\nabla\cdot\bm{u}_{t},\xi)+\frac{1}{\lambda}(\xi_{t},\xi)-\frac{\alpha}{\lambda}(p_{t},\xi)-\frac{\beta}{\lambda}(T_{t},\xi)$	$\displaystyle=0,$
$\displaystyle-\frac{\alpha}{\lambda}(\xi_{t},p)+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{t},p)+(\bm{K}\nabla p,\nabla p)+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{t},p)$	$\displaystyle=(g,p),$
$\displaystyle-\frac{\beta}{\lambda}(\xi_{t},T)+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{t},T)+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{t},T)+(\bm{\Theta}\nabla T,\nabla T)$	$\displaystyle=(H_{s},T).$

Adding the above four equations together, we have

		$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}),\bm{\varepsilon}(\bm{u}_{t}))+\frac{1}{\lambda}(\xi_{t},\xi)-\frac{\alpha}{\lambda}(p_{t},\xi)-\frac{\beta}{\lambda}(T_{t},\xi)$		(2.10)
		$\displaystyle-\frac{\alpha}{\lambda}(\xi_{t},p)+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{t},p)+(\bm{K}\nabla p,\nabla p)+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{t},p)$
		$\displaystyle-\frac{\beta}{\lambda}(\xi_{t},T)+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{t},T)+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{t},T)+(\bm{\Theta}\nabla T,\nabla T)$
		$\displaystyle-(\bm{f},\bm{u}_{t})-(g,p)-(H_{s},T)=0,$

From the expression of (2.7), we note that

	$\displaystyle\frac{d}{dt}\Big{(}\frac{1}{2\lambda}\\|\alpha p+\beta T-\xi\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\\|p-T\\|_{L^{2}(\Omega)}^{2}\Big{)}$	(2.11)
$\displaystyle=$	$\displaystyle\frac{1}{\lambda}(\alpha p+\beta T-\xi,\alpha p_{t}+\beta T_{t}-\xi_{t})+(c_{0}-b_{0})(p,p_{t})+(a_{0}-b_{0})(T,T_{t})$
	$\displaystyle+b_{0}(p-T,p_{t}-T_{t})$
$\displaystyle=$	$\displaystyle\frac{1}{\lambda}(\xi_{t},\xi)-\frac{\alpha}{\lambda}(p_{t},\xi)-\frac{\beta}{\lambda}(T_{t},\xi)-\frac{\alpha}{\lambda}(\xi_{t},p)+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{t},p)$
	$\displaystyle+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{t},p)-\frac{\beta}{\lambda}(\xi_{t},T)+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{t},T)$
	$\displaystyle+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{t},T).$

Integrating (2.10) from $0$ to $t$ and using (2.11), we derive (2.6). The bound for $\xi$ follows from the inf-sup condition between $\bm{V}$ and $Q$ and the Korn’s inequality. We have the following inequality holds

$\displaystyle\zeta_{0}\\|\xi\\|_{L^{2}(\Omega)}=$	$\displaystyle\sup_{\bm{v}\in\bm{V}}\frac{\|(\text{div }\bm{v},\xi)\|}{\\|\bm{v}\\|_{H^{1}(\Omega)}}$	(2.12)
$\displaystyle\leq$	$\displaystyle\sup_{\bm{v}\in\bm{V}}\frac{\|2\mu(\varepsilon(\bm{u}),\varepsilon(\bm{v}))+(\bm{f},\bm{v})\|}{\\|\bm{v}\\|_{H^{1}(\Omega)}}$
$\displaystyle\leq$	$\displaystyle C_{K}\big{(}2\mu\\|\bm{\varepsilon}(u)\\|_{L^{2}(\Omega)}+\\|\bm{f}\\|_{L^{2}(\Omega)}\big{)},$

where constant $\zeta_{0}$ is from the inf-sup condition and $C_{k}$ is from the Korn’s inequality. This deduces (2.8). ∎

Then, the energy law 2.2 implies the following priori estimates.

Theorem 2.1.

Assume that the Assumptions of (A1)-(A3) hold. There exists a positive constant $C=C(\|\bm{\varepsilon}(\bm{u}^{0})\|_{L^{2}(\Omega)}$ , $\|\xi^{0}\|_{L^{2}(\Omega)}$ , $\|p^{0}\|_{L^{2}(\Omega)}$ , $\|T^{0}\|_{L^{2}(\Omega)}$ , $\|\bm{f}\|_{L^{2}(\Omega)}$ , $\|g\|_{L^{2}(\Omega)}$ , $\|H_{s}\|_{L^{2}(\Omega)})$ such that

		$\displaystyle\sqrt{2\mu}\\|\bm{\varepsilon}(\bm{u})\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{\frac{1}{\lambda}}\\|\alpha p+\beta T-\xi\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{c_{0}-b_{0}}\\|p\\|_{L^{\infty}(J;L^{2}(\Omega))}$		(2.13)
		$\displaystyle+\sqrt{a_{0}-b_{0}}\\|T\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{b_{0}}\\|p-T\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{k_{m}}\\|\nabla p\\|_{L^{2}(J;L^{2}(\Omega))}$
		$\displaystyle+\sqrt{\theta_{m}}\\|\nabla T\\|_{L^{2}(J;L^{2}(\Omega))}\leq C.$

Proof.

From (2.10), by using the Assumption (A1), the Cauchy-Schwarz inequality, and Young’s inequality, we have

		$\displaystyle\mu\frac{d}{dt}\\|\bm{\varepsilon}(\bm{u})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\frac{d}{dt}\\|\alpha p+\beta T-\xi\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\frac{d}{dt}\\|p\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\frac{d}{dt}\\|T\\|_{L^{2}(\Omega)}^{2}$		(2.14)
		$\displaystyle+\frac{b_{0}}{2}\frac{d}{dt}\\|p-T\\|_{L^{2}(\Omega)}^{2}+k_{m}\\|\nabla p\\|_{L^{2}(\Omega)}^{2}+\theta_{m}\\|\nabla T\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle\leq(\bm{f},\bm{u}_{t})+(g,p)+(H_{s},T).$

Integrating (2.14) from $0$ to $t$ , applying the Cauchy-Schwarz inequality, Young’s inequality, we can find a constant $\mu^{*}$ such that

		$\displaystyle\mu\\|\bm{\varepsilon}(\bm{u})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p+\beta T-\xi\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T\\|_{L^{2}(\Omega)}^{2}$		(2.15)
		$\displaystyle+\frac{b_{0}}{2}\\|p-T\\|_{L^{2}(\Omega)}^{2}+\int_{0}^{t}k_{m}\\|\nabla p\\|_{L^{2}(\Omega)}^{2}ds+\int_{0}^{t}\theta_{m}\\|\nabla T\\|_{L^{2}(\Omega)}^{2}ds$
		$\displaystyle\leq\mu\\|\bm{\varepsilon}(\bm{u}^{0})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p^{0}+\beta T^{0}-\xi^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T^{0}\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\frac{b_{0}}{2}\\|p^{0}-T^{0}\\|_{L^{2}(\Omega)}^{2}+2\mu^{}\\|\bm{f}\\|_{L^{2}(\Omega)}^{2}+\frac{1}{4\mu^{}}\\|\bm{u}\\|_{L^{2}(\Omega)}^{2}+\frac{1}{4\mu^{*}}\\|\bm{u}^{0}\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\frac{1}{2}\int_{0}^{t}(\\|g\\|_{L^{2}(\Omega)}^{2}+\\|H_{s}\\|_{L^{2}(\Omega)}^{2}+\\|p\\|_{L^{2}(\Omega)}^{2}+\\|T\\|_{L^{2}(\Omega)}^{2}).$

Then, by the Poincaré inequality, Gronwall’s inequality, and Korn’s inequality, choosing $\mu^{*}=\frac{1}{2\mu C_{K}}$ , we get

		$\displaystyle\mu\\|\bm{\varepsilon}(\bm{u})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p+\beta T-\xi\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T\\|_{L^{2}(\Omega)}^{2}$		(2.16)
		$\displaystyle+\frac{b_{0}}{2}\\|p-T\\|_{L^{2}(\Omega)}^{2}+\int_{0}^{t}k_{m}\\|\nabla p\\|_{L^{2}(\Omega)}^{2}ds+\int_{0}^{t}\theta_{m}\\|\nabla T\\|_{L^{2}(\Omega)}^{2}ds$
		$\displaystyle\lesssim\frac{3\mu}{2}\mu\\|\bm{\varepsilon}(\bm{u}^{0})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p^{0}+\beta T^{0}-\xi^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T^{0}\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\frac{b_{0}}{2}\\|p^{0}-T^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{1}{\mu}\\|\bm{f}\\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{2}\\|\bm{\varepsilon}(\bm{u})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\int_{0}^{t}(\\|g\\|_{L^{2}(\Omega)}^{2}+\\|H_{s}\\|_{L^{2}(\Omega)}^{2}),$

which implies that (2.13) holds. ∎

Lemma 2.3.

Assume that the Assumptions of (A1)-(A3) hold. Then, there exist positive constants $C_{1}=C_{1}(\|\bm{\varepsilon}(\bm{u}^{0})\|_{L^{2}(\Omega)},~{}\|\xi^{0}\|_{L^{2}(\Omega)},~{}\|p^{0}\|_{L^{2}(\Omega)},~{}\|T^{0}\|_{L^{2}(\Omega)}$ , $\|\bm{f}\|_{L^{2}(\Omega)}$ , $\|g\|_{L^{2}(\Omega)},\|H_{s}\|_{L^{2}(\Omega)}$ , $\|\nabla p(0)\|_{L^{2}(\Omega)}$ , $\|\nabla T(0)\|_{L^{2}(\Omega)})$ and $C_{2}=C_{2}(C_{1},\|\bm{u}^{0}\|_{H^{2}(\Omega)},~{}\|p^{0}\|_{H^{2}(\Omega)},~{}\|T^{0}\|_{H^{2}(\Omega)})$ such that

		$\displaystyle\sqrt{2\mu}\\|\bm{\varepsilon}(\bm{u}_{t})\\|_{L^{2}(J;L^{2}(\Omega))}+\sqrt{\frac{1}{\lambda}}\\|\alpha p_{t}+\beta T_{t}-\xi_{t}\\|_{L^{2}(J;L^{2}(\Omega))}+\sqrt{c_{0}-b_{0}}\\|p_{t}\\|_{L^{2}(J;L^{2}(\Omega))}$		(2.17)
		$\displaystyle+\sqrt{a_{0}-b_{0}}\\|T_{t}\\|_{L^{2}(J;L^{2}(\Omega))}+\sqrt{b_{0}}\\|p_{t}-T_{t}\\|_{L^{2}(J;L^{2}(\Omega))}+\sqrt{k_{m}}\\|\nabla p_{t}\\|_{L^{\infty}(J;L^{2}(\Omega))}$
		$\displaystyle\sqrt{\theta_{m}}\\|\nabla T_{t}\\|_{L^{\infty}(J;L^{2}(\Omega))}\leq C_{1},$

and

		$\displaystyle\sqrt{2\mu}\\|\bm{\varepsilon}(\bm{u}_{t})\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{\frac{1}{\lambda}}\\|\alpha p_{t}+\beta T_{t}-\xi_{t}\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{c_{0}-b_{0}}\\|p_{t}\\|_{L^{\infty}(J;L^{2}(\Omega))}$		(2.18)
		$\displaystyle+\sqrt{a_{0}-b_{0}}\\|T_{t}\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{b_{0}}\\|p_{t}-T_{t}\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{2k_{m}}\\|\nabla p_{t}\\|_{L^{2}(J;L^{2}(\Omega))}$
		$\displaystyle\sqrt{2\theta_{m}}\\|\nabla T_{t}\\|_{L^{2}(J;L^{2}(\Omega))}\leq C_{2}.$

Proof.

Differentiating the first and the second equation of system (2.5) with respect to $t$ , taking $\bm{v}=\bm{u}_{t},~{}\phi=-\xi_{t},~{}q=p_{t}$ , and $S=T_{t}$ respectively, and adding the resulting equations together we have

		$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}_{t}),\bm{\varepsilon}(\bm{u}_{t}))+\frac{1}{\lambda}(\xi_{t},\xi_{t})-\frac{\alpha}{\lambda}(p_{t},\xi_{t})-\frac{\beta}{\lambda}(T_{t},\xi_{t})$		(2.19)
		$\displaystyle-\frac{\alpha}{\lambda}(\xi_{t},p_{t})+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{t},p_{t})+(\bm{K}\nabla p,\nabla p_{t})+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{t},p_{t})$
		$\displaystyle-\frac{\beta}{\lambda}(\xi_{t},T_{t})+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{t},T_{t})+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{t},T_{t})+(\bm{\Theta}\nabla T,\nabla T_{t})$
		$\displaystyle=(g,p_{t})+(H_{s},T_{t}).$

Based on some algebraic operations, we have

		$\displaystyle 2\mu\\|\bm{\varepsilon}(\bm{u}_{t})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{\lambda}\\|\alpha p_{t}+\beta T_{t}-\xi_{t}\\|_{L^{2}(\Omega)}^{2}+(c_{0}-b_{0})\\|p_{t}\\|_{L^{2}(\Omega)}^{2}+(a_{0}-b_{0})\\|T_{t}\\|_{L^{2}(\Omega)}^{2}$		(2.20)
		$\displaystyle+b_{0}\\|p_{t}-T_{t}\\|_{L^{2}(\Omega)}^{2}+(\bm{K}\nabla p,\nabla p_{t})+(\bm{\Theta}\nabla T,\nabla T_{t})$
		$\displaystyle=(g,p_{t})+(H_{s},T_{t}).$

Integrating equation (2.20) from $0$ to $t$ , by Assumption (A1) we have

		$\displaystyle\int_{0}^{t}\Big{(}2\mu\\|\bm{\varepsilon}(\bm{u}_{t})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{\lambda}\\|\alpha p_{t}+\beta T_{t}-\xi_{t}\\|_{L^{2}(\Omega)}^{2}+(c_{0}-b_{0})\\|p_{t}\\|_{L^{2}(\Omega)}^{2}+(a_{0}-b_{0})\\|T_{t}\\|_{L^{2}(\Omega)}^{2}$		(2.21)
		$\displaystyle+b_{0}\\|p_{t}-T_{t}\\|_{L^{2}(\Omega)}^{2}\Big{)}ds+k_{m}\\|\nabla p\\|_{L^{2}(\Omega)}^{2}+\theta_{m}\\|\nabla T\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle\leq(g,p(t)-p(0))+(H_{s},T(t)-T(0))+k_{M}\\|\nabla p(0)\\|_{L^{2}(\Omega)}^{2}+\theta_{M}\\|\nabla T(0)\\|_{L^{2}(\Omega)}^{2},$

which, when combined with equation (2.13) in lemma 2.2 and Gronwall’s inequality, implies that (2.17) holds.

Differentiating the first, the third, and the fourth equations of (2.5) with respect to $t$ , taking $\bm{v}=\bm{u}_{tt},\phi=-\xi_{t},q=p_{t}$ , and $S=T_{t}$ respectively, and adding the resulting equations we have

		$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}_{tt}),\bm{\varepsilon}(\bm{u}_{t}))+\frac{1}{\lambda}(\xi_{tt},\phi_{t})-\frac{\alpha}{\lambda}(p_{tt},\phi_{t})-\frac{\beta}{\lambda}(T_{tt},\phi_{t})$		(2.22)
		$\displaystyle-\frac{\alpha}{\lambda}(\xi_{tt},p_{t})+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{tt},p_{t})+(\bm{K}\nabla p_{t},\nabla p_{t})+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{tt},p_{t})$
		$\displaystyle-\frac{\beta}{\lambda}(\xi_{tt},T_{t})+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{tt},T_{t})+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{tt},T_{t})+(\bm{\Theta}\nabla T_{t},\nabla T_{t})$
		$\displaystyle=0.$

By the Assumption (A1), the Cauchy-Schwarz inequality and Young’s inequality, we have

		$\displaystyle\mu\frac{d}{dt}\\|\bm{\varepsilon}(\bm{u}_{t})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\frac{d}{dt}\\|\alpha p_{t}+\beta T_{t}-\xi_{t}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\frac{d}{dt}\\|p_{t}\\|_{L^{2}(\Omega)}^{2}$		(2.23)
		$\displaystyle+\frac{a_{0}-b_{0}}{2}\frac{d}{dt}\\|T_{t}\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\frac{d}{dt}\\|p_{t}-T_{t}\\|_{L^{2}(\Omega)}^{2}+k_{m}\\|\nabla p_{t}\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\theta_{m}\\|\nabla T_{t}\\|_{L^{2}(\Omega)}^{2}\leq 0.$

Integrating equation (2.23) from 0 to $t$ , we have

		$\displaystyle\mu\\|\bm{\varepsilon}(\bm{u}_{t})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p_{t}+\beta T_{t}-\xi_{t}\\|_{L^{2}(\Omega)}^{2}+\frac{(c_{0}-b_{0})}{2}\\|p_{t}\\|_{L^{2}(\Omega)}^{2}+\frac{(a_{0}-b_{0})}{2}\\|T_{t}\\|_{L^{2}(\Omega)}^{2}$		(2.24)
		$\displaystyle+\frac{b_{0}}{2}\\|p_{t}-T_{t}\\|_{L^{2}(\Omega)}^{2}+\int_{0}^{t}k_{m}\\|\nabla p_{t}\\|_{L^{2}(\Omega)}^{2}+\int_{0}^{t}\theta_{m}\\|\nabla T_{t}\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle\leq\mu\\|\bm{\varepsilon}(\bm{u}_{t}(0))\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p_{t}(0)+\beta T_{t}(0)-\xi_{t}(0)\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p_{t}(0)\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\frac{a_{0}-b_{0}}{2}\\|T_{t}(0)\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\\|p_{t}(0)-T_{t}(0)\\|_{L^{2}(\Omega)}^{2},$

which implies that (2.18) holds. ∎

With the assistance of Theorem 2.1, Lemma 2.2, and Lemma 2.3, we can establish the existence and uniqueness of the weak solution to the problem (2.4). Since the system is linear, the proof can be easily obtained by employing the Galerkin method and the compactness argument, along with the insights presented in [9].

Theorem 2.2.

Assume that the Assumptions of (A1)-(A3) hold. Let $\bm{u}^{0}\in H^{1}(\Omega)$ , $p^{0}\in L^{2}(\Omega)$ , $T^{0}\in L^{2}(\Omega)$ , $\bm{f}\in L^{2}(\Omega)$ , $g\in L^{2}(\Omega)$ , and $H_{s}\in L^{2}(\Omega)$ . Then a unique weak solution exists to the problem (2.4). Furthermore, if $\bm{u}^{0}\in H^{2}(\Omega)$ , $p^{0}\in H^{2}(\Omega)$ , $T^{0}\in H^{2}(\Omega)$ , $\bm{f}\in L^{2}(\Omega)$ , $g\in L^{2}(\Omega)$ , and $H_{s}\in L^{2}(\Omega)$ , then the energy estimations (2.17) and (2.18) for the solution hold true.

3 Numerical algorithms

This section introduces numerical algorithms derived from the four-field formulation previously discussed. Let $h$ denote the maximum diameter of all elements in the mesh, and for $h>0$ , let $\mathcal{T}_{h}$ represent a family of triangulations of the domain $\Omega$ composed of triangular elements. The triangulations are assumed to be shape-regular and quasi-uniform. The discrete spaces for $\bm{u}$ and $\xi$ are denoted by $\bm{V}_{h}\times Q_{h}$ , which are required to satisfy the following fundamental stability condition.

\displaystyle\inf_{\xi\in Q_{h}}\sup_{\bm{v}\in\bm{V}_{h}}\frac{(\nabla\cdot\bm{v},\xi)}{\|\bm{v}\|_{H^{1}(\Omega)}\|\xi\|_{L^{2}(\Omega)}}\geq\gamma>0,

(3.1)

where $\gamma$ is independent of $h$ , which means $V_{h}\times Q_{h}$ is a stable Stokes pair. As an example, we utilize the Taylor-Hood elements for the pair $(\bm{u},\xi)$ , specifically $(\bm{P}_{2},P_{1})$ Lagrange finite elements, and $P_{1}$ Lagrange finite elements for both the fluid pressure $p$ and temperature $T$ . The finite element spaces defined on $\mathcal{T}_{h}$ are as follows:

		$\displaystyle\bm{V}_{h}=\left\{\bm{v}_{h}\in[H_{0}^{1}(\Omega)]^{2}\cap\bm{C}^{0}(\bar{\Omega}):\bm{v}_{h}\|_{K}\in\mathbb{P}_{2}(K),\,\forall K\in\mathcal{T}_{h}\right\},$
		$\displaystyle Q_{h}=\left\{\phi_{h}\in L^{2}(\Omega)\cap C^{0}(\bar{\Omega}):\phi_{h}\|_{K}\in\mathbb{P}_{1}(K),\,\forall K\in\mathcal{T}_{h}\right\},$
		$\displaystyle W_{h}=\left\{q_{h}\in H_{0}^{1}(\Omega)\cap C^{0}(\bar{\Omega}):q_{h}\|_{K}\in\mathbb{P}_{1}(K),\,\forall K\in\mathcal{T}_{h}\right\}.$

It is worth noting that other stable Stokes element pairs, such as the MINI element, could also be used.

For time discretization, we consider an equidistant partition consisting of a series of points from $0$ to $T$ with a constant step size of $\Delta t$ . Specifically, we define $\bm{u}^{n}$ as $\bm{u}(t_{n})$ , where $t_{n}$ represents the $n$ -th point in the partition. Similarly, we define $\xi^{n}=\xi(t_{n})$ , $p^{n}=p(t_{n})$ and $T^{n}=T(t_{n})$ .

3.1 A coupled algorithm

Suppose that an initial value $(\bm{u}_{h}^{0},~{}\xi_{h}^{0},~{}p_{h}^{0},~{}T_{h}^{0})\in\bm{V}_{h}\times Q_{h}\times W_{h}\times W_{h}$ is provided. Using the back Euler method to (2.5), for the given finite element spaces $\bm{V}_{h}\subset\bm{V}$ , $Q_{h}\subset Q$ , $W_{h}\subset W$ , the coupled algorithm (Alg. 1) reads as, find $(\bm{u}_{h}^{n},~{}\xi_{h}^{n},~{}p_{h}^{n},~{}T_{h}^{n})\in\bm{V}_{h}\times Q_{h}\times W_{h}\times W_{h}$ such that for all $(\bm{v}_{h},~{}\phi_{h},~{}q_{h},~{}S_{h})\in\bm{V}_{h}\times Q_{h}\times W_{h}\times W_{h}$ ,


$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}^{n}_{h}),\bm{\varepsilon}(\bm{v}_{h}))-(\nabla\cdot\bm{v}_{h},\xi^{n}_{h})$	$\displaystyle=(\bm{f},\bm{v}_{h})$	(3.2a)

$\displaystyle-(\nabla\cdot\bm{u}^{n}_{h},\phi_{h})-\frac{1}{\lambda}(\xi^{n}_{h},\phi_{h})+\frac{\alpha}{\lambda}(p^{n}_{h},\phi_{h})+\frac{\beta}{\lambda}(T^{n}_{h},\phi_{h})$	$\displaystyle=0,$	(3.2b)

$\displaystyle-\frac{\alpha}{\lambda}(\xi_{h}^{n},q_{h})+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{h}^{n},~{}q_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{h}^{n},~{}q_{h})+\Delta t(\bm{K}\nabla p^{n}_{h},~{}\nabla q_{h})$	$\displaystyle=$	(3.2c)
$\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{h}^{n-1},~{}q_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{h}^{n-1},~{}q_{h})-\frac{\alpha}{\lambda}(\xi_{h}^{n-1},~{}q_{h})$	$\displaystyle+\Delta t(g,~{}q_{h}),$	(3.2c)

$\displaystyle-\frac{\beta}{\lambda}(\xi_{h}^{n},S_{h})+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{h}^{n},S_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{h}^{n},~{}S_{h})+\Delta t(\bm{\Theta}\nabla T^{n}_{h},~{}\nabla S_{h})$	$\displaystyle=$	(3.2d)
$\displaystyle(a_{0}+\frac{\beta^{2}}{\lambda})(T_{h}^{n-1},~{}S_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{h}^{n-1},~{}S_{h})-\frac{\beta}{\lambda}(\xi_{h}^{n-1},~{}S_{h})$	$\displaystyle+\Delta t(H_{s},~{}S_{h}).$	(3.2d)

To avoid solving the above problem in a fully coupled manner, the system (3.2a)-(3.2d) can be reformulated into two decoupled subproblems. If the term $\frac{\alpha}{\lambda}(p^{n}_{h},\phi_{h})+\frac{\beta}{\lambda}(T^{n}_{h},\phi_{h})$ in (3.2b) were moved to the right-hand side, equations (3.2a) and (3.2b) become equivalent to the variational form of the generalized Stokes problem for the variables $\bm{u}_{h}^{n}$ and $\xi_{h}^{n}$ . Similarly, moving the terms $\frac{\alpha}{\lambda}(\xi_{h}^{n},q_{h})$ in (3.2c) and $\frac{\beta}{\lambda}(\xi_{h}^{n},S_{h})$ in (3.2d) to the right-hand sides yields variational forms corresponding to reaction-diffusion equations for $p_{h}^{n}$ and $T_{h}^{n}$ respectively. These observations facilitate the design of two time-extrapolation-based decoupled algorithms: The first algorithm solves for $(\bm{u}_{h}^{n},\xi_{h}^{n})$ first, followed by solving the reaction-diffusion equations for $(p_{h}^{n},T_{h}^{n})$ . Conversely, the second algorithm begins with $(p_{h}^{n},T_{h}^{n})$ and subsequently computes $(\bm{u}_{h}^{n},\xi_{h}^{n})$ .

Based on prior experience with Biot’s model [12, 16], these decoupling strategies typically impose stability constraints that necessitate small time step sizes. Specifically, the time step $\Delta t$ must generally be of order $O(h^{2})$ . Consequently, time-extrapolation-based decoupled algorithms may fail to ensure stability or accuracy if the time step size is too large.

3.2 An iteratively decoupled algorithm

To overcome the stability constraints, we propose an iterative decoupling algorithm for solving $(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})$ . We have the following Iteratively Decoupled Algorithm (Algorithm 2):

For each time step $t_{n}$ ( $n\geq 1$ ), let the previous step states $(p_{h}^{n,i-1},T_{h}^{n,i-1},\xi_{h}^{n,i-1})$ be given, with the initial values set as $p_{h}^{n,0}=p_{h}^{n-1}$ , $T_{h}^{n,0}=T_{h}^{n-1}$ , and $\xi_{h}^{n,0}=\xi_{h}^{n-1}$ . Let $i$ be the iteration index. For a given $n$ , the $i$ -th iteration is formulated as:

step 1: Given $\xi_{h}^{n,i-1}$ , find $(p_{h}^{n,i},T_{h}^{n,i})$ such that

$\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{h}^{n,i},q_{h})$	$\displaystyle+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{h}^{n,i},q_{h})+\Delta t(\bm{K}\nabla p^{n,i}_{h},\nabla q_{h})$	(3.3)
$\displaystyle=$	$\displaystyle\Delta t(g,q_{h})+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{h}^{n-1},q_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{h}^{n-1},q_{h})$
	$\displaystyle+\frac{\alpha}{\lambda}(\xi_{h}^{n,i-1},q_{h})-\frac{\alpha}{\lambda}(\xi_{h}^{n-1},q_{h}),$
$\displaystyle(a_{0}+\frac{\beta^{2}}{\lambda})(T_{h}^{n,i},S_{h})$	$\displaystyle+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{h}^{n,i},S_{h})+\Delta t(\bm{\Theta}\nabla T^{n,i}_{h},\nabla S_{h})$
$\displaystyle=$	$\displaystyle\Delta t(H_{s},S_{h})+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{h}^{n-1},S_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{h}^{n-1},S_{h})$
	$\displaystyle+\frac{\beta}{\lambda}(\xi_{h}^{n,i-1},S_{h})-\frac{\beta}{\lambda}(\xi_{h}^{n-1},S_{h}),$

step 2: Given $p_{h}^{n,i},T_{h}^{n,i}$ ,find $(\bm{u}_{h}^{n,i},\xi_{h}^{n,i})$ such that

	$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}^{n,i}_{h}),\bm{\varepsilon}(\bm{v}_{h}))-(\nabla\cdot\bm{v}_{h},\xi^{n,i}_{h})$	$\displaystyle=(\bm{f},\bm{v}_{h}),$		(3.4)
	$\displaystyle(\nabla\cdot\bm{u}^{n,i}_{h},\phi_{h})+\frac{1}{\lambda}(\xi^{n,i}_{h},\phi_{h})$	$\displaystyle=\frac{\alpha}{\lambda}(p^{n,i}_{h},\phi_{h})+\frac{\beta}{\lambda}(T^{n,i}_{h},\phi_{h}).$		(3.4)

The algorithm proceeds as follows: 1. Compute iteration Step 1 and Step 2 in succession. 2. Repeat the iteration process until the sequence $(\bm{u}_{h}^{n,i},\xi_{h}^{n,i},p_{h}^{n,i},T_{h}^{n,i})$ converges to within a prescribed tolerance.

For simplicity, the backward Euler scheme has been employed for time discretization in both Alg. 1 and Alg. 2. However, it is important to note that alternative higher-order time-stepping methods could also be applied effectively in this context.

4 Convergence analysis

This section establishes the convergence of Alg. 2, demonstrating that the iterative sequences $(\bm{u}_{h}^{n,i},\xi_{h}^{n,i},p_{h}^{n,i},T_{h}^{n,i})$ converge to the solution $(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})$ of the coupled algorithm Alg. 1 as $i\to\infty$ . Additionally, we present the error analysis for Alg. 1, confirming its unconditional stability and convergence. Assuming the discrete spaces $\bm{V}_{h}$ , $Q_{h}$ , and $W_{h}$ , we establish that the time discretization error is of order $O(\Delta t)$ , while the energy-norm errors for $\bm{u}$ and $\xi$ are of order $O(h^{2})$ . To begin, we review an essential lemma that aids us in proving the theorem.

Lemma 4.1.

There exists a constant $C_{inv}>0$ such that the following inverse inequality is satisfied:

\|\bm{v}\|_{H^{1}(\Omega)}\leq C_{inv}h^{-1}\|\bm{v}\|_{L^{2}(\Omega)},\quad\forall\bm{v}\in\bm{V}_{h},

where $h=\max_{K\in\mathcal{T}_{h}}h_{K}$ , as established in [11].

To proceed the analysis, in the remaining part of the article, we let $e_{\bm{u}}^{i},e_{\xi}^{i},e_{p}^{i},e_{T}^{i}$ be differences between the solutions at the iteration $i$ of problem (3.3)-(3.4) and the solutions to problem (3.2a)-(3.2d), i.e.

\displaystyle(e_{\bm{u}}^{i},e_{\xi}^{i},e_{p}^{i},e_{T}^{i}):=(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})-(\bm{u}_{h}^{n,i},\xi_{h}^{n,i},p_{h}^{n,i},T_{h}^{n,i}).

(4.1)

Now, we present the following key result, establishing the convergence of Alg.2.

Theorem 4.1.

Let $(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})$ be the solution of the coupled problem (3.2a)-(3.2d). Let $(\bm{u}_{h}^{n,i},\xi_{h}^{n,i},p_{h}^{n,i},T_{h}^{n,i})$ be the solution of the iteratively decoupled problem (3.3)-(3.4). Assuming the Assumptions (A1)-(A3) hold, then, the solution sequence $\{(\bm{u}_{h}^{n,i},\xi_{h}^{n,i},p_{h}^{n,i},T_{h}^{n,i})\}_{i=1}^{\infty}$ generated by the iteratively decoupled algorithm converges to the solution $(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})$ of the coupled problem.

Proof.

Subtracting equation (3.3)-(3.4) from equation (3.2a)-(3.2d) respectively, we obtain


$\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})(e_{p}^{i},q_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(e_{T}^{i},q_{h})$	$\displaystyle+\Delta t(\bm{K}\nabla e_{p}^{i},\nabla q_{h})=\frac{\alpha}{\lambda}(e_{\xi}^{i-1},q_{h}),$	(4.2a)

$\displaystyle(a_{0}+\frac{\beta^{2}}{\lambda})(e_{T}^{i},S_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(e_{p}^{i},S_{h})+\Delta t(\bm{\Theta}\nabla e_{T}^{i},\nabla S_{h})=$	$\displaystyle\frac{\beta}{\lambda}(e_{\xi}^{i-1},S_{h}),$	(4.2b)

$\displaystyle 2\mu(\bm{\varepsilon}(e_{\bm{u}}^{i}),\bm{\varepsilon}(\bm{v}_{h}))-(\nabla\cdot\bm{v}_{h},e_{\xi}^{i})=$	$\displaystyle 0,$	(4.2c)

$\displaystyle(\nabla\cdot e_{\bm{u}}^{i},\phi_{h})+\frac{1}{\lambda}(e_{\xi}^{i},\phi_{h})=$	$\displaystyle\frac{\alpha}{\lambda}(e_{p}^{i},\phi_{h})+\frac{\beta}{\lambda}(e_{T}^{i},\phi_{h}).$	(4.2d)

Taking $q_{h}=e_{p}^{i}$ , $S_{h}=e_{T}^{i}$ , $\bm{v}_{h}=e_{\bm{u}}^{i}$ and $\phi_{h}=e_{\xi}^{i}$ in equation (4.2a), (4.2b), (4.2c) and (4.2d), we have

\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})\|e_{p}^{i}\|_{L^{2}(\Omega)}^{2}+(\frac{\alpha\beta}{\lambda}-b_{0})(e_{T}^{i},e_{p}^{i})

\displaystyle+\Delta t(\bm{K}\nabla e_{p}^{i},\nabla e_{p}^{i})=\frac{\alpha}{\lambda}(e_{\xi}^{i-1},e_{p}^{i}),

(4.3)

\displaystyle(a_{0}+\frac{\beta^{2}}{\lambda})\|e_{T}^{i}\|_{L^{2}(\Omega)}^{2}+(\frac{\alpha\beta}{\lambda}-b_{0})(e_{T}^{i},e_{p}^{i})+\Delta t(\bm{\Theta}\nabla e_{T}^{i},\nabla e_{T}^{i})=

\displaystyle\frac{\beta}{\lambda}(e_{\xi}^{i-1},e_{T}^{i}).

(4.4)

\displaystyle 2\mu(\bm{\varepsilon}(e_{\bm{u}}^{i}),\bm{\varepsilon}(e_{\bm{u}}^{i}))-(\nabla\cdot e_{\bm{u}}^{i},e_{\xi}^{i})=

\displaystyle 0.

(4.5)

\displaystyle(\nabla\cdot e_{\bm{u}}^{i},e_{\xi}^{i})+\frac{1}{\lambda}(e_{\xi}^{i},e_{\xi}^{i})=

\displaystyle\frac{\alpha}{\lambda}(e_{p}^{i},e_{\xi}^{i})+\frac{\beta}{\lambda}(e_{T}^{i},e_{\xi}^{i}).

(4.6)

Summing up equation (4.3) and equation (4.4), we have

	$\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})\\|e_{p}^{i}\\|_{L^{2}(\Omega)}^{2}+(a_{0}+\frac{\beta^{2}}{\lambda})\\|e_{T}^{i}\\|_{L^{2}(\Omega)}^{2}+2(\frac{\alpha\beta}{\lambda}-b_{0})(e_{p}^{i},e_{T}^{i})$	(4.7)
	$\displaystyle+\Delta t(\bm{K}\nabla e_{p}^{i},\nabla e_{p}^{i})+\Delta t(\bm{\Theta}\nabla e_{T}^{i},\nabla e_{T}^{i})$
$\displaystyle=$	$\displaystyle\frac{\alpha}{\lambda}(e_{p}^{i},e_{\xi}^{i-1})+\frac{\beta}{\lambda}(e_{T}^{i},e_{\xi}^{i-1})=(\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i},e_{\xi}^{i-1}).$

By the Cauchy-Schwarz inequality and (4.7), we have

	$\displaystyle(c_{0}-b_{0})\\|e_{p}^{i}\\|_{L^{2}(\Omega)}^{2}+(a_{0}-b_{0})\\|e_{T}^{i}\\|_{L^{2}(\Omega)}^{2}+\lambda\\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\\|_{L^{2}(\Omega)}^{2}$	(4.8)
$\displaystyle\leq$	$\displaystyle c_{0}\\|e_{p}^{i}\\|_{L^{2}(\Omega)}^{2}-2b_{0}(e_{p}^{i},e_{T}^{i})+a_{0}\\|e_{T}^{i}\\|_{L^{2}(\Omega)}^{2}+\frac{\alpha^{2}}{\lambda}\\|e_{p}^{i}\\|_{L^{2}(\Omega)}^{2}+2\frac{\alpha\beta}{\lambda}(e_{p}^{i},e_{T}^{i})+\frac{\beta^{2}}{\lambda}\\|e_{T}^{i}\\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq$	$\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})\\|e_{p}^{i}\\|_{L^{2}(\Omega)}^{2}+(a_{0}+\frac{\beta^{2}}{\lambda})\\|e_{T}^{i}\\|_{L^{2}(\Omega)}^{2}+2(\frac{\alpha\beta}{\lambda}-b_{0})(e_{p}^{i},e_{T}^{i})$
	$\displaystyle+\Delta t(\bm{K}\nabla e_{p}^{i},\nabla e_{p}^{i})+\Delta t(\bm{\Theta}\nabla e_{T}^{i},\nabla e_{T}^{i})$
$\displaystyle=$	$\displaystyle(\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i},e_{\xi}^{i-1})$
$\displaystyle\leq$	$\displaystyle\\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\\|_{L^{2}(\Omega)}\\|e_{\xi}^{i-1}\\|_{L^{2}(\Omega)}.$

By simultaneously dividing both sides of the equation (4.8) by $\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\|_{L^{2}(\Omega)}$ which is positive, otherwise we yields $\|e_{p}^{i}\|_{L^{2}(\Omega)}=\|e_{T}^{i}\|_{L^{2}(\Omega)}$ , we obtain

\displaystyle\frac{(c_{0}-b_{0})\|e_{p}^{i}\|_{L^{2}(\Omega)}^{2}+(a_{0}-b_{0})\|e_{T}^{i}\|_{L^{2}(\Omega)}^{2}}{\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\|_{L^{2}(\Omega)}^{2}}\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\|_{L^{2}(\Omega)}+\lambda\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\|_{L^{2}(\Omega)}

\displaystyle\leq\|e_{\xi}^{i-1}\|_{L^{2}(\Omega)}.

Let

\mathbb{D}:=\frac{1}{\lambda^{2}}\frac{(c_{0}-b_{0})\|e_{p}^{i}\|_{L^{2}(\Omega)}^{2}+(a_{0}-b_{0})\|e_{T}^{i}\|_{L^{2}(\Omega)}^{2}}{\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\|_{L^{2}(\Omega)}^{2}}=\frac{(c_{0}-b_{0})\|e_{p}^{i}\|_{L^{2}(\Omega)}^{2}+(a_{0}-b_{0})\|e_{T}^{i}\|_{L^{2}(\Omega)}^{2}}{\|\alpha e_{p}^{i}+\beta e_{T}^{i}\|_{L^{2}(\Omega)}^{2}},

which is positive and independent of $\lambda$ , we can next derive

	$\displaystyle\lambda^{2}\mathbb{D}\\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\\|_{L^{2}(\Omega)}+\lambda\\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\\|_{L^{2}(\Omega)}$	$\displaystyle\leq\\|e_{\xi}^{i-1}\\|_{L^{2}(\Omega)}$		(4.9)
	$\displaystyle\\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\\|_{L^{2}(\Omega)}$	$\displaystyle\leq\frac{1}{\lambda^{2}\mathbb{D}+\lambda}\\|e_{\xi}^{i-1}\\|_{L^{2}(\Omega)}.$		(4.9)

Summing up equation (4.5) and equation (4.6), we obtain

\displaystyle 2\mu(\bm{\varepsilon}(e_{\bm{u}}^{i}),\bm{\varepsilon}(e_{\bm{u}}^{i}))+\frac{1}{\lambda}\|e_{\xi}^{i}\|_{L^{2}(\Omega)}^{2}=

\displaystyle(\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i},e_{\xi}^{i}).

(4.10)

Dropping the first positive term, and using the conclusion of (4.9), we have

\displaystyle\|e_{\xi}^{i}\|_{L^{2}(\Omega)}^{2}\leq

\displaystyle\lambda(\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i},e_{\xi}^{i})\leq\lambda\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\|_{L^{2}(\Omega)}\|e_{\xi}^{i}\|_{L^{2}(\Omega)},

(4.11)

and then

\displaystyle\|e_{\xi}^{i}\|_{L^{2}(\Omega)}\leq

\displaystyle\lambda\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\|_{L^{2}(\Omega)}\leq\frac{\lambda}{\lambda^{2}\mathbb{D}+\lambda}\|e_{\xi}^{i-1}\|_{L^{2}(\Omega)}.

(4.12)

Therefore, the convergence of $\xi_{h}^{n,i}$ is proved. Applying Lemma 2.1 to equation (4.5), we see that

	$\displaystyle 2\mu\\|\bm{\varepsilon}(e_{\bm{u}}^{i})\\|_{L^{2}(\Omega)}^{2}$	$\displaystyle=2\mu(\bm{\varepsilon}(e_{\bm{u}}^{i}),\bm{\varepsilon}(e_{\bm{u}}^{i}))=(\nabla\cdot e_{\bm{u}}^{i},e_{\xi}^{i})$		(4.13)
		$\displaystyle\leq\\|\nabla\cdot e_{\bm{u}}^{i}\\|_{L^{2}(\Omega)}\\|e_{\xi}^{i}\\|_{L^{2}(\Omega)}\leq C_{K}\\|\bm{\varepsilon}(e_{\bm{u}}^{i})\\|_{L^{2}(\Omega)}\\|e_{\xi}^{i}\\|_{L^{2}(\Omega)},$		(4.13)

which yields the convergence of $\bm{u}_{h}^{n,i}$ . By neglecting the third term on the left-hand side of (4.8) and utilizing the result of (4.9) on the right-hand side of (4.8), we obtain

	$\displaystyle(c_{0}-b_{0})\\|e_{p}^{i}\\|_{L^{2}(\Omega)}^{2}+(a_{0}-b_{0})\\|e_{T}^{i}\\|_{L^{2}(\Omega)}^{2}$	$\displaystyle\leq\\|\frac{\alpha}{\lambda}e_{p}^{i}+\frac{\beta}{\lambda}e_{T}^{i}\\|_{L^{2}(\Omega)}\\|e_{\xi}^{i-1}\\|_{L^{2}(\Omega)}$		(4.14)
		$\displaystyle\leq\frac{1}{\lambda^{2}\mathbb{D}+\lambda}\\|e_{\xi}^{i-1}\\|_{L^{2}(\Omega)}^{2},$		(4.14)

which implies the convergence of $p_{h}^{n,i}$ and $T_{h}^{n,i}$ . ∎

Remark 1.

The theorem remains valid in the special cases where either $c_{0}-b_{0}=0$ or $a_{0}-b_{0}=0$ , provided the other term is nonzero, as the proof remains unaffected under these conditions.

If both $c_{0}-b_{0}=0$ and $a_{0}-b_{0}=0$ , we can still demonstrate the convergence of the iterative decoupling algorithm. Assuming $c_{0}-b_{0}=a_{0}-b_{0}=0$ , it follows that $\mathbb{D}=0$ . From inequality (4.12), it is evident that $\|e_{\xi}^{i}\|_{L^{2}(\Omega)}$ forms a monotonically non-increasing sequence bounded below by zero, ensuring its convergence. Let

\lim_{i\to\infty}\|e_{\xi}^{i}\|_{L^{2}(\Omega)}=\rho>0.

From (4.10) and (4.12), we derive:

2\mu(\bm{\varepsilon}(e_{\bm{u}}^{i}),\bm{\varepsilon}(e_{\bm{u}}^{i}))+\frac{1}{\lambda}\|e_{\xi}^{i}\|_{L^{2}(\Omega)}^{2}\leq\|e_{\xi}^{i}\|_{L^{2}(\Omega)}\|e_{\xi}^{i-1}\|_{L^{2}(\Omega)}.

Given $\rho>0$ , this implies $\lim_{i\to\infty}\|\bm{\varepsilon}(e_{\bm{u}}^{i})\|_{L^{2}(\Omega)}=0$ as $i\to\infty$ . Utilizing the discrete inf-sup condition and referencing (4.2c), we obtain:

\gamma\|e_{\xi}^{i}\|_{L^{2}(\Omega)}\leq\sup_{\bm{v}_{h}\in\bm{V}_{h}}\frac{|(\nabla\cdot\bm{v}_{h},e_{\xi}^{i})|}{\|\bm{v}_{h}\|_{H^{1}(\Omega)}}=\frac{|2\mu(\bm{\varepsilon}(\bm{v}_{h}),\bm{\varepsilon}(e_{\bm{u}}^{i}))|}{\|\bm{v}_{h}\|_{H^{1}(\Omega)}}\lesssim\|\bm{\varepsilon}(e_{\bm{u}}^{i})\|_{L^{2}(\Omega)}.

This leads to $\rho=0$ , which contradicts the assumption $\rho>0$ . Thus, $\|e_{\xi}^{i}\|_{L^{2}(\Omega)}\to 0$ as $i\to\infty$ . Consequently, it follows that $\|e_{\bm{u}}^{i}\|_{L^{2}(\Omega)}\to 0$ , $\|e_{p}^{i}\|_{L^{2}(\Omega)}\to 0$ , and $\|e_{T}^{i}\|_{L^{2}(\Omega)}\to 0$ as $i\to\infty$ .

Here, we denote $d_{t}\eta_{n}:=(\eta^{n}-\eta^{n-1})/\Delta t$ , where $\eta$ can be a vector or a scalar.

Lemma 4.2.

Let $\{(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})\}$ be solution of the coupled algorithm, then the following identity holds:

\displaystyle J_{h}^{l}+S_{h}^{l}=J_{h}^{0},\text{ for }l\geq 1,

(4.15)

where for $l\geq 0$ ,

	$\displaystyle J_{h}^{l}:=$	$\displaystyle\mu\\|\bm{\varepsilon}(\bm{u}_{h}^{l})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p_{h}^{l}+\beta T_{h}^{l}-\xi_{h}^{l}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p_{h}^{l}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T_{h}^{l}\\|_{L^{2}(\Omega)}^{2}$		(4.16)
		$\displaystyle+\frac{b_{0}}{2}\\|p_{h}^{l}-T_{h}^{l}\\|_{L^{2}(\Omega)}^{2}-(\bm{f},\bm{u}_{h}^{l}),$		(4.16)

and

$\displaystyle S_{h}^{l}:=$	$\displaystyle\Delta t\sum_{n=1}^{l}\Big{[}\Delta t\Big{(}\mu\\|\bm{\varepsilon}(d_{t}\bm{u}_{h}^{l})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|d_{t}(\alpha p_{h}^{l}+\beta T_{h}^{l}-\xi_{h}^{l})\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|d_{t}p_{h}^{l}\\|_{L^{2}(\Omega)}^{2}$	(4.17)
	$\displaystyle+\frac{a_{0}-b_{0}}{2}\\|d_{t}T_{h}^{l}\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\\|p_{h}^{l}-T_{h}^{l}\\|_{L^{2}(\Omega)}^{2}\Big{)}$
	$\displaystyle+\bm{K}(\nabla p_{h}^{n},\nabla p_{h}^{n})+\bm{\Theta}(\nabla T_{h}^{n},\nabla T_{h}^{n})-(g,\nabla p_{h}^{n})-(H_{s},\nabla T_{h}^{n})\Big{]}.$

Proof.

Using operator $d_{t}$ for the equation (3.2b), taking $\bm{v}_{h}=d_{t}\bm{u}_{h}^{n},\phi_{h}=-\xi_{h}^{n},q_{h}=p_{h}^{n}$ and $S_{h}=T_{h}^{n}$ , we have

$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}^{n}_{h}),\bm{\varepsilon}(d_{t}\bm{u}_{h}^{n}))-(\nabla\cdot d_{t}\bm{u}_{h}^{n},\xi^{n}_{h})$	$\displaystyle=(\bm{f},d_{t}\bm{u}_{h}^{n}),$	(4.18)
$\displaystyle(\nabla\cdot d_{t}\bm{u}^{n}_{h},\xi_{h}^{n})+\frac{1}{\lambda}(d_{t}\xi^{n}_{h},\xi_{h}^{n})-\frac{\alpha}{\lambda}(d_{t}p^{n}_{h},\xi_{h}^{n})-\frac{\beta}{\lambda}(d_{t}T^{n}_{h},\xi_{h}^{n})$	$\displaystyle=0,$
$\displaystyle-\frac{\alpha}{\lambda}(\xi_{h}^{n},p_{h}^{n})+(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{h}^{n},p_{h}^{n})+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{h}^{n},p_{h}^{n})+\Delta t(\bm{K}\nabla p^{n}_{h},\nabla p_{h}^{n})$	$\displaystyle=$
$\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})(p_{h}^{n-1},p_{h}^{n})+(\frac{\alpha\beta}{\lambda}-b_{0})(T_{h}^{n-1},p_{h}^{n})-\frac{\alpha}{\lambda}(\xi_{h}^{n-1},p_{h}^{n})$	$\displaystyle+\Delta t(g,p_{h}^{n}),$
$\displaystyle-\frac{\beta}{\lambda}(\xi_{h}^{n},T_{h}^{n})+(a_{0}+\frac{\beta^{2}}{\lambda})(T_{h}^{n},T_{h}^{n})+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{h}^{n},T_{h}^{n})+\Delta t(\bm{\Theta}\nabla T^{n}_{h},\nabla T_{h}^{n})$	$\displaystyle=$
$\displaystyle(a_{0}+\frac{\beta^{2}}{\lambda})(T_{h}^{n-1},T_{h}^{n})+(\frac{\alpha\beta}{\lambda}-b_{0})(p_{h}^{n-1},T_{h}^{n})-\frac{\beta}{\lambda}(\xi_{h}^{n-1},T_{h}^{n})$	$\displaystyle+\Delta t(H_{s},T_{h}^{n}).$

Summing up the four equations of system (4.18) and by the definition of $d_{t}$ , we have

		$\displaystyle 2\mu(\bm{\varepsilon}(\bm{u}_{h}^{n}),\bm{\varepsilon}(d_{t}\bm{u}_{h}^{n}))+\frac{1}{\lambda}(d_{t}\xi_{h}^{n},\xi_{h}^{n})-\frac{\alpha}{\lambda}(d_{t}p_{h}^{n},\xi_{h}^{n})-\frac{\beta}{\lambda}(d_{t}T_{h}^{n},\xi_{h}^{n})$		(4.19)
		$\displaystyle-\frac{\alpha}{\lambda}(d_{t}\xi_{h}^{n},p_{h}^{n})+(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}p_{h}^{n},p_{h}^{n})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}T_{h}^{n},p_{h}^{n})+(\bm{K}\nabla p_{h}^{n},\nabla p_{h}^{n})$
		$\displaystyle-\frac{\beta}{\lambda}(d_{t}\xi_{h}^{n},T_{h}^{n})+(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}T_{h}^{n},T_{h}^{n})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}p_{h}^{n},T_{h}^{n})+(\bm{\Theta}\nabla T_{h}^{n},\nabla T_{h}^{n})$
		$\displaystyle=(\bm{f},d_{t}\bm{u}_{h}^{n})+(g,p_{h}^{n})+(H_{s},T_{h}^{n}).$

Then, using the identity

2(\eta^{n},d_{t}\eta^{n})=d_{t}\|\eta^{n}\|_{L^{2}(\Omega)}^{2}+\Delta t\|d_{t}\eta^{n}\|_{L^{2}(\Omega)}^{2},

(4.20)

we have

		$\displaystyle\mu d_{t}\\|\bm{\varepsilon}(\bm{u}_{h}^{n})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}d_{t}\\|\alpha p_{h}^{n}+\beta T_{h}^{n}-\xi_{h}^{n}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}d_{t}\\|p_{h}^{n}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}d_{t}\\|T_{h}^{n}\\|_{L^{2}(\Omega)}^{2}$		(4.21)
		$\displaystyle+\frac{b_{0}}{2}d_{t}\\|p_{h}^{n}-T_{h}^{n}\\|_{L^{2}(\Omega)}^{2}+\Delta t\Big{(}\mu\\|d_{t}\bm{\varepsilon}(\bm{u}_{h}^{n})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|d_{t}(\alpha p_{h}^{n}+\beta T_{h}^{n}-\xi_{h}^{n})\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\frac{c_{0}-b_{0}}{2}\\|d_{t}p_{h}^{n}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|d_{t}T_{h}^{n}\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\\|d_{t}(p_{h}^{n}-T_{h}^{n})\\|_{L^{2}(\Omega)}^{2}\Big{)}+(\bm{K}\nabla p_{h}^{n},\nabla p_{h}^{n})$
		$\displaystyle+(\bm{\Theta}\nabla T_{h}^{n},\nabla T_{h}^{n})=(\bm{f},d_{t}\bm{u}_{h}^{n})+(g,p_{h}^{n})+(H_{s},T_{h}^{n}).$

Applying the summation operator $\Delta t\sum_{n=1}^{l}$ to both sides of the above equation, we obtain (4.15). ∎

For the discrete finite element spaces $\bm{V}_{h}$ , $Q_{h}$ , and $W_{h}$ defined in Section 3.1, we introduce the projection operators $\Pi^{\bm{V}_{h}}:\bm{V}\to\bm{V}_{h}$ , $\Pi^{Q_{h}}:Q\to Q_{h}$ , $\Pi^{W_{p,h}}:W\to W_{h}$ , and $\Pi^{W_{T,h}}:W\to W_{h}$ , which satisfy the following equations: for all $(\bm{v}_{h},\phi_{h},p_{h},T_{h})\in\bm{V}_{h}\times Q_{h}\times W_{h}\times W_{h}$ ,

2\mu(\bm{\varepsilon}(\Pi^{\bm{V}_{h}}\bm{u}),\bm{\varepsilon}(\bm{v}_{h}))-(\nabla\cdot\bm{v}_{h},\Pi^{Q_{h}}\xi)=2\mu(\bm{\varepsilon}(\bm{u}),\bm{\varepsilon}(\bm{v}_{h}))-(\nabla\cdot\bm{v}_{h},\xi),

(4.22)

(\nabla\cdot\Pi^{\bm{V}_{h}}\bm{u},\phi_{h})+\frac{1}{\lambda}(\Pi^{Q_{h}}\xi,\phi_{h})=(\nabla\cdot\bm{u},\phi_{h})+\frac{1}{\lambda}(\xi,\phi_{h})

(4.23)

\bm{K}(\nabla\Pi^{W_{p,h,}}p,\nabla q_{h})=\bm{K}(\nabla p,\nabla q_{h}),

(4.24)

\bm{\Theta}(\nabla\Pi^{W_{T,h}}T,\nabla S_{h})=\bm{\Theta}(\nabla T,\nabla S_{h}).

(4.25)

Furthermore, for the operators $\Pi^{\bm{V}_{h}}$ , $\Pi^{Q_{h}}$ , $\Pi^{W_{p,h,}}$ and $\Pi^{W_{T,h}}$ we have following properties.

\|\nabla(\Pi^{\bm{V}_{h}}\bm{u}-\bm{u})\|_{L^{2}(\Omega)}+\|\Pi^{Q_{h}}\xi-\xi\|_{L^{2}(\Omega)}\lesssim h^{2}\|\bm{u}\|_{H^{3}(\Omega)}+h^{2}\|\xi\|_{H^{2}(\Omega)},

(4.26)

\|\Pi^{W_{p,h,}}p-p\|_{L^{2}(\Omega)}+h\|\nabla(\Pi^{W_{p,h,}}p-p)\|_{L^{2}(\Omega)}\lesssim h^{2}\|p\|_{H^{2}(\Omega)},

(4.27)

and

\|\Pi^{W_{T,h}}T-T\|_{L^{2}(\Omega)}+h\|\nabla(\Pi^{W_{T,h}}T-T)\|_{L^{2}(\Omega)}\lesssim h^{2}\|T\|_{H^{2}(\Omega)}.

(4.28)

For convenience, we introduce the following notations:

		$\displaystyle e_{\bm{u}}^{n}=\bm{u}^{n}-\bm{u}_{h}^{n}=(\bm{u}^{n}-\Pi^{\bm{V}_{h}}\bm{u}^{n})+(\Pi^{\bm{V}_{h}}\bm{u}^{n}-\bm{u}_{h}^{n}):=e_{\bm{u}}^{I,n}+e_{\bm{u}}^{h,n},$		(4.29)
		$\displaystyle e_{\xi}^{n}=\xi^{n}-\xi_{h}^{n}=(\xi^{n}-\Pi^{Q_{h}}\xi^{n})+(\Pi^{Q_{h}}\xi^{n}-\xi_{h}^{n}):=e_{\xi}^{I,n}+e_{\xi}^{h,n},$
		$\displaystyle e_{p}^{n}=p^{n}-p_{h}^{n}=(p^{n}-\Pi^{W_{p,h}}p^{n})+(\Pi^{W_{p,h}}p^{n}-p_{h}^{n}):=e_{p}^{I,n}+e_{p}^{h,n},$
		$\displaystyle e_{T}^{n}=T^{n}-T_{h}^{n}=(T^{n}-\Pi^{W_{p,h}}T^{n})+(\Pi^{W_{p,h}}T^{n}-T_{h}^{n}):=e_{T}^{I,n}+e_{T}^{h,n}.$

Next, we formulate a discrete energy law, which serves as a discretized counterpart to the continuous statement presented in Lemma 2.2.

Lemma 4.3.

Let $\{(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})\}$ be defined by the coupled algorithm, then the following identity holds:

	$\displaystyle E_{h}^{l}+\Delta t\sum_{n=1}^{l}\big{(}(\bm{K}\nabla e_{p}^{h,n},\nabla e_{p}^{h,n})+(\bm{\Theta}\nabla e_{T}^{h,n},\nabla e_{T}^{h,n})\big{)}$	(4.30)
	$\displaystyle+\Delta t^{2}\sum_{n=1}^{l}\Big{(}\mu\\|d_{t}\bm{\varepsilon}(e_{\bm{u}}^{h,n})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|d_{t}(\alpha e_{p}^{h,n}+\beta e_{T}^{h,n}-e_{\xi}^{h,n})\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|d_{t}e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2}$
	$\displaystyle+\frac{a_{0}-b_{0}}{2}\\|d_{t}e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\\|e_{p}^{h,n}-e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}\Big{)}$
$\displaystyle=$	$\displaystyle E_{h}^{0}+\Delta t\sum_{n=1}^{l}\Big{(}(\nabla\cdot(d_{t}\bm{u}^{n}-\bm{u}_{t}^{n}),e_{\xi}^{h,n})+\frac{1}{\lambda}(d_{t}\xi^{n}-\xi_{t}^{n},e_{\xi}^{h,n})-\frac{\alpha}{\lambda}(d_{t}p^{n}-p_{t}^{n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}T^{n}-T_{t}^{n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\alpha}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{p}^{h,n})+(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{p}^{h,n})$
	$\displaystyle+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{p}^{h,n}),$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{T}^{h,n})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{T}^{h,n})$
	$\displaystyle+(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{T}^{h,n})\Big{)},$

where

	$\displaystyle E_{h}^{n}$	$\displaystyle=\mu\\|\bm{\varepsilon}(e_{\bm{u}}^{h,n})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha e_{p}^{h,n}+\beta e_{T}^{h,n}-e_{\xi}^{h,n}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}$		(4.31)
		$\displaystyle+\frac{b_{0}}{2}\\|e_{p}^{h,n}-e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}.$		(4.31)

Proof.

Utilizing the first equation of the continuous weak formulation (2.5), the finite element formulations (3.2a) and (4.22), along with the definitions of $e_{\bm{u}}^{h,n}$ and $e_{\xi}^{h,n}$ as given in equation (4.29), we derive the following relation:

2\mu(\bm{\varepsilon}(e_{\bm{u}}^{h,n}),\bm{\varepsilon}(\bm{v}_{h}))-(\nabla\cdot\bm{v}_{h},e_{\xi}^{h,n})=0.

(4.32)

Furthermore, by combining the second equation of (2.5), the finite element formulation (3.2b), and the projection property (4.23), we obtain

	$\displaystyle(\nabla\cdot d_{t}e_{\bm{u}}^{h,n},\phi_{h})+\frac{1}{\lambda}(d_{t}e_{\xi}^{h,n},\phi_{h})-\frac{\alpha}{\lambda}(d_{t}e_{p}^{h,n},\phi_{h})-\frac{\beta}{\lambda}(d_{t}e_{T}^{h,n},\phi_{h})$	(4.33)
$\displaystyle=$	$\displaystyle(\nabla\cdot(d_{t}\bm{u}^{n}-\bm{u}_{t}^{n}),\phi_{h})+\frac{1}{\lambda}(d_{t}\xi^{n}-\xi_{t}^{n},\phi_{h})-\frac{\alpha}{\lambda}(d_{t}\Pi^{W_{p,h}}p^{n}-p_{t}^{n},\phi_{h})$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},\phi_{h}).$

Using the last two equations of system (2.5) and of system (3.2a)-(3.2d), equation (4.24) and (4.25), we have

	$\displaystyle-\frac{\alpha}{\lambda}(d_{t}e_{\xi}^{h,n},q_{h})+(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}e_{p}^{h,n},q_{h})+\bm{(}\bm{K}\nabla e_{p}^{h,n},\nabla q_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}e_{T}^{h,n},q_{h})$	(4.34)
$\displaystyle=$	$\displaystyle-\frac{\alpha}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},q_{h})+(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},q_{h})$
	$\displaystyle+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},q_{h}),$

and

	$\displaystyle-\frac{\beta}{\lambda}(d_{t}e_{\xi}^{h,n},S_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}e_{p}^{h,n},S_{h})+(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}e_{T}^{h,n},S_{h})+(\bm{\Theta}\nabla e_{T}^{h,n},\nabla S_{h})$	(4.35)
$\displaystyle=$	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},S_{h})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},S_{h})$
	$\displaystyle+(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},S_{h}).$

Setting $\bm{v}_{h}=d_{t}e_{\bm{u}}^{h,n}$ in (4.32), $\phi_{h}=e_{\xi}^{h,n}$ in (4.33) and $q_{h}=e_{p}^{h,n}$ , $S_{h}=e_{T}^{h,n}$ in (4.34)-(4.35) and adding the resulting equations together, we derive

	$\displaystyle 2\mu(\bm{\varepsilon}(e_{\bm{u}}^{h,n}),\bm{\varepsilon}(d_{t}e_{\bm{u}}^{h,n}))-(\nabla\cdot d_{t}e_{\bm{u}}^{h,n},e_{\xi}^{h,n})$	(4.36)
	$\displaystyle+(\nabla\cdot d_{t}e_{\bm{u}}^{h,n},e_{\xi}^{h,n})+\frac{1}{\lambda}(d_{t}e_{\xi}^{h,n},e_{\xi}^{h,n})-\frac{\alpha}{\lambda}(d_{t}e_{p}^{h,n},e_{\xi}^{h,n})-\frac{\beta}{\lambda}(d_{t}e_{T}^{h,n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\alpha}{\lambda}(d_{t}e_{\xi}^{h,n},e_{p}^{h,n})+(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}e_{p}^{h,n},e_{p}^{h,n})+(\bm{K}\nabla e_{p}^{h,n},\nabla e_{p}^{h,n})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}e_{T}^{h,n},e_{p}^{h,n})$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}e_{\xi}^{h,n},e_{T}^{h,n})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}e_{p}^{h,n},e_{T}^{h,n})+(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}e_{T}^{h,n},e_{T}^{h,n})+(\bm{\Theta}\nabla e_{T}^{h,n},\nabla e_{T}^{h,n})$
$\displaystyle=$	$\displaystyle(\nabla\cdot(d_{t}\bm{u}^{n}-\bm{u}_{t}^{n}),e_{\xi}^{h,n})+\frac{1}{\lambda}(d_{t}\xi^{n}-\xi_{t}^{n},e_{\xi}^{h,n})-\frac{\alpha}{\lambda}(d_{t}\Pi^{W_{p,h}}p^{n}-p_{t}^{n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\alpha}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{p}^{h,n})+(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{p}^{h,n})$
	$\displaystyle+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{p}^{h,n})$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{T}^{h,n})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{T}^{h,n})$
	$\displaystyle+(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{T}^{h,n}).$

By the identity (4.20), we have

	$\displaystyle d_{t}\Big{(}\mu\\|\bm{\varepsilon}(e_{\bm{u}}^{h,n})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha e_{p}^{h,n}+\beta e_{T}^{h,n}-e_{\xi}^{h,n}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}$	(4.37)
	$\displaystyle+\frac{b_{0}}{2}\\|e_{p}^{h,n}-e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}\Big{)}+(\bm{K}\nabla e_{p}^{h,n},\nabla e_{p}^{h,n})+(\bm{\Theta}\nabla e_{T}^{h,n},\nabla e_{T}^{h,n})$
	$\displaystyle+\Delta t\Big{(}\mu\\|d_{t}\bm{\varepsilon}(e_{\bm{u}}^{h,n})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|d_{t}(\alpha e_{p}^{h,n}+\beta e_{T}^{h,n}-e_{\xi}^{h,n})\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|d_{t}e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2}$
	$\displaystyle+\frac{a_{0}-b_{0}}{2}\\|d_{t}e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\\|e_{p}^{h,n}-e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}\Big{)}$
$\displaystyle=$	$\displaystyle(\nabla\cdot(d_{t}\bm{u}^{n}-\bm{u}_{t}^{n}),e_{\xi}^{h,n})+\frac{1}{\lambda}(d_{t}\xi^{n}-\xi_{t}^{n},e_{\xi}^{h,n})-\frac{\alpha}{\lambda}(d_{t}\Pi^{W_{p,h}}p^{n}-p_{t}^{n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\alpha}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{p}^{h,n})+(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{p}^{h,n})$
	$\displaystyle+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{p}^{h,n})$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{T}^{h,n})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{T}^{h,n})$
	$\displaystyle+(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{T}^{h,n}).$

Then, applying the summation operator $\Delta t\sum_{n=1}^{l}$ to both sides of the above equation, we get (4.30). ∎

Divide the interval $[0,\tau]$ into $l$ equal parts and let $t_{n}$ represent the time value at the n-th time step, $t_{0}=0$ and $t_{l}=\tau$ . For simplicity, $X\lesssim Y$ is used to denote an inequality $X\leq CY$ , where $C$ is a positive constant independent of mesh sizes $h$ .

Theorem 4.2.

Let $\{(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})\}$ be the solution of the coupled algorithm, then the following estimates hold:

	$\displaystyle\max_{0\leq n\leq l}\Big{[}\mu\\|\bm{\varepsilon}(e_{\bm{u}}^{h,n})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha e_{p}^{h,n}+\beta e_{T}^{h,n}-e_{\xi}^{h,n}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2}$	(4.38)
	$\displaystyle+\frac{a_{0}-b_{0}}{2}\\|e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}+\frac{b_{0}}{2}\\|e_{p}^{h,n}-e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}\Big{]}$
	$\displaystyle+\Delta t\sum_{n=1}^{l}\big{(}(\bm{K}\nabla e_{p}^{h,n},\nabla e_{p}^{h,n})+(\bm{\Theta}\nabla e_{T}^{h,n},\nabla e_{T}^{h,n})\big{)}$
$\displaystyle\leq$	$\displaystyle\tilde{C}_{1}\Delta t^{2}+\tilde{C}_{2}h^{4},$

where

		$\displaystyle\tilde{C}_{1}=\tilde{C}_{1}(\int_{t_{0}}^{\tau}\\|\bm{u}_{tt}\\|_{H^{1}(\Omega)}^{2},\int_{t_{0}}^{\tau}\\|\xi_{tt}\\|_{L^{2}(\Omega)}^{2},\int_{t_{0}}^{\tau}\\|p_{tt}\\|_{L^{2}(\Omega)}^{2},\int_{t_{0}}^{\tau}\\|T_{tt}\\|_{L^{2}(\Omega)}^{2})$		(4.39)
		$\displaystyle\tilde{C}_{2}=\tilde{C}_{2}(\int_{t_{0}}^{\tau}\\|\xi_{t}\\|_{H^{2}(\Omega)}^{2},\int_{t_{0}}^{\tau}\\|p_{t}\\|_{H^{2}(\Omega)}^{2},\int_{t_{0}}^{\tau}\\|T_{t}\\|_{H^{2}(\Omega)}^{2}),$		(4.39)

Proof.

For the initial conditions, we set $\bm{u}_{h}^{0}=\Pi^{\bm{V}_{h}}\bm{u}^{0}$ , $\xi_{h}^{0}=\Pi^{Q_{h}}\xi^{0}$ , $p_{h}^{0}=\Pi^{W_{p,h}}p^{0}$ , and $T_{h}^{0}=\Pi^{W_{T,h}}T^{0}$ . Using these definitions, the following inequality is obtained from (4.30).

	$\displaystyle E_{h}^{l}+\Delta t\sum_{n=1}^{l}\big{(}(\bm{K}\nabla e_{p}^{h,n},\nabla e_{p}^{h,n})+(\bm{\Theta}\nabla e_{T}^{h,n},\nabla e_{T}^{h,n})\big{)}$	(4.40)
$\displaystyle\leq$	$\displaystyle\Delta t\sum_{n=1}^{l}\Big{(}(\nabla\cdot(d_{t}\bm{u}^{n}-\bm{u}_{t}^{n}),e_{\xi}^{h,n})+\frac{1}{\lambda}(d_{t}\xi^{n}-\xi_{t}^{n},e_{\xi}^{h,n})-\frac{\alpha}{\lambda}(d_{t}p^{n}-p_{t}^{n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}T^{n}-T_{t}^{n},e_{\xi}^{h,n})$
	$\displaystyle-\frac{\alpha}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{p}^{h,n})+(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{p}^{h,n})$
	$\displaystyle+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{p}^{h,n}),$
	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{T}^{h,n})+(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{T}^{h,n})$
	$\displaystyle+(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{T}^{h,n})\Big{)}$
$\displaystyle:=$	$\displaystyle D_{1}+D_{2}+D_{3}+D_{4}+D_{5}+D_{6}+D_{7}+D_{8}+D_{9}+D_{10}.$

We now estimate each term on the right-hand side of (4.40). Using the Cauchy-Schwarz inequality, Taylor series expansion, and Proposition 3.1. in [14], we can bound $D_{1}$ by

	$\displaystyle D_{1}=$	$\displaystyle\Delta t\sum_{n=1}^{l}(\nabla\cdot(d_{t}\bm{u}^{n}-\bm{u}_{t}^{n}),e_{\xi}^{h,n})$
	$\displaystyle=$	$\displaystyle\sum_{n=1}^{l}(\nabla\cdot(\bm{u}^{n}-\bm{u}^{n-1}-\Delta t\bm{u}_{t}^{n}),e_{\xi}^{h,n})$
	$\displaystyle\lesssim$	$\displaystyle\frac{1}{\Delta t}\sum_{n=1}^{l}\\|\bm{u}^{n}-\bm{u}^{n-1}-\Delta t\bm{u}_{t}^{n}\\|_{H^{1}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{\xi}^{h,n}\\|_{L^{2}(\Omega)}^{2}$
	$\displaystyle=$	$\displaystyle\frac{1}{\Delta t}\sum_{n=1}^{l}\\|\int_{t_{n-1}}^{t_{n}}\bm{u}_{tt}(s)(t_{n-1}-s)ds\\|_{H^{1}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{\xi}^{h,n}\\|_{L^{2}(\Omega)}^{2}$
	$\displaystyle\lesssim$	$\displaystyle\Delta t^{2}\int_{t_{0}}^{\tau}\\|\bm{u}_{tt}\\|_{H^{1}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{\xi}^{h,n}\\|_{L^{2}(\Omega)}^{2}$

Similarly, $D_{2}$ , $D_{3}$ and $D_{4}$ can be bounded as follows:

D_{2}=\Delta t\sum_{n=1}^{l}\frac{1}{\lambda}(d_{t}\xi^{n}-\xi_{t}^{n},e_{\xi}^{h,n})\lesssim\frac{1}{\lambda}(\Delta t^{2}\int_{t_{0}}^{\tau}\|\xi_{tt}\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\|e_{\xi}^{h,n}\|_{L^{2}(\Omega)}^{2}),\\

D_{3}=\Delta t\sum_{n=1}^{l}\frac{\alpha}{\lambda}(d_{t}p^{n}-p_{t}^{n},e_{\xi}^{h,n})\lesssim\frac{\alpha}{\lambda}(\Delta t^{2}\int_{t_{0}}^{\tau}\|p_{tt}\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\|e_{\xi}^{h,n}\|_{L^{2}(\Omega)}^{2}),\\

and

D_{4}=\Delta t\sum_{n=1}^{l}\frac{\beta}{\lambda}(d_{t}T^{n}-T_{t}^{n},e_{\xi}^{h,n})\lesssim\frac{\beta}{\lambda}(\Delta t^{2}\int_{t_{0}}^{\tau}\|T_{tt}\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\|e_{\xi}^{h,n}\|_{L^{2}(\Omega)}^{2}).\\

By use of estimate (4.26), the Cauchy-Schwarz inequality, Taylor series expansion, and Proposition 3.1. in [14], we see that $D_{5}$ satisfies

	$\displaystyle D_{5}=$	$\displaystyle-\Delta t\sum_{n=1}^{l}\frac{\alpha}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{p}^{h,n})$
	$\displaystyle=$	$\displaystyle-\sum_{n=1}^{l}\frac{\alpha}{\lambda}(\Pi^{Q_{h}}\xi^{n}-\Pi^{Q_{h}}\xi^{n-1}-\Delta t\xi_{t}^{n},e_{p}^{h,n})$
	$\displaystyle=$	$\displaystyle-\sum_{n=1}^{l}\frac{\alpha}{\lambda}(\Pi^{Q_{h}}\xi^{n}-\Pi^{Q_{h}}\xi^{n-1}-(\xi^{n}-\xi^{n-1})+(\xi^{n}-\xi^{n-1})-\Delta t\xi_{t}^{n},e_{p}^{h,n})$
	$\displaystyle\lesssim$	$\displaystyle\frac{\alpha}{\lambda}(h^{4}\frac{1}{\Delta t}\sum_{n=1}^{l}\\|\xi^{n}-\xi^{n-1}\\|_{H^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|\xi_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2})$
	$\displaystyle\lesssim$	$\displaystyle\frac{\alpha}{\lambda}(h^{4}\int_{t_{0}}^{\tau}\\|\xi_{t}\\|_{H^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|\xi_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2}),$

and similarly,

	$\displaystyle D_{6}=$	$\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{p}^{h,n})$
	$\displaystyle\lesssim$	$\displaystyle(c_{0}+\frac{\alpha^{2}}{\lambda})(h^{4}\int_{t_{0}}^{\tau}\\|p_{t}\\|_{H^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|p_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2}),$

	$\displaystyle D_{7}=$	$\displaystyle(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{p}^{h,n})$
	$\displaystyle\lesssim$	$\displaystyle(\frac{\alpha\beta}{\lambda}-b_{0})(h^{4}\int_{t_{0}}^{\tau}\\|T_{t}\\|_{H^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|T_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{p}^{h,n}\\|_{L^{2}(\Omega)}^{2}),$

	$\displaystyle D_{8}=$	$\displaystyle-\frac{\beta}{\lambda}(d_{t}\Pi^{Q_{h}}\xi^{n}-\xi_{t}^{n},e_{T}^{h,n})$
	$\displaystyle\lesssim$	$\displaystyle\frac{\beta}{\lambda}(h^{4}\int_{t_{0}}^{\tau}\\|\xi_{t}\\|_{H^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|\xi_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}),$

	$\displaystyle D_{9}=$	$\displaystyle(\frac{\alpha\beta}{\lambda}-b_{0})(d_{t}\Pi^{W_{p,h,}}p^{n}-p_{t}^{n},e_{T}^{h,n})$
	$\displaystyle\lesssim$	$\displaystyle\|\frac{\alpha\beta}{\lambda}-b_{0}\|(h^{4}\int_{t_{0}}^{\tau}\\|p_{t}\\|_{H^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|p_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2})$

	$\displaystyle D_{10}=$	$\displaystyle(a_{0}+\frac{\beta^{2}}{\lambda})(d_{t}\Pi^{W_{T,h}}T^{n}-T_{t}^{n},e_{T}^{h,n})$
	$\displaystyle\lesssim$	$\displaystyle(a_{0}+\frac{\beta^{2}}{\lambda})(h^{4}\int_{t_{0}}^{\tau}\\|T_{t}\\|_{H^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|T_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t\sum_{n=1}^{l}\\|e_{T}^{h,n}\\|_{L^{2}(\Omega)}^{2}).$

The above bounds lead to

	$\displaystyle E_{h}^{l}+\Delta t\sum_{n=1}^{l}\big{(}(\bm{K}\nabla e_{p}^{h,n},\nabla e_{p}^{h,n})+(\bm{\Theta}\nabla e_{T}^{h,n},\nabla e_{T}^{h,n})\big{)}$	(4.41)
$\displaystyle\lesssim$	$\displaystyle\Delta t^{2}\int_{t_{0}}^{\tau}\\|\bm{u}_{tt}\\|_{H^{1}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|\xi_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|p_{tt}\\|_{L^{2}(\Omega)}^{2}+\Delta t^{2}\int_{t_{0}}^{\tau}\\|T_{tt}\\|_{L^{2}(\Omega)}^{2}$
	$\displaystyle+h^{4}\int_{t_{0}}^{\tau}\\|\xi_{t}\\|_{H^{2}(\Omega)}^{2}+h^{4}\int_{t_{0}}^{\tau}\\|p_{t}\\|_{H^{2}(\Omega)}^{2}+h^{4}\int_{t_{0}}^{\tau}\\|T_{t}\\|_{H^{2}(\Omega)}^{2},$

which by the discrete Gronwall’s inequality implies that (4.38) holds. ∎

Theorem 4.3.

Let $\{(\bm{u}_{h}^{n},\xi_{h}^{n},p_{h}^{n},T_{h}^{n})\}$ be defined by the coupled algorithm,then the following estimates hold:

		$\displaystyle\max_{0\leq n\leq l}\Big{[}\mu\\|\bm{\varepsilon}(e_{\bm{u}}^{n})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|e_{\xi}^{n}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|e_{p}^{n}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|e_{T}^{n}\\|_{L^{2}(\Omega)}^{2}\Big{]}$		(4.42)
	$\displaystyle\leq$	$\displaystyle\tilde{C}_{1}\Delta t^{2}+\tilde{C}_{2}h^{4},$		(4.42)

where $\tilde{C}_{1}$ and $\tilde{C}_{2}$ are defined the same as in Theorem 4.2.

Proof.

By use of the (4.29) and Theorem 4.2, we obtain the estimate (4.42).

∎

5 Numerical experiments

This section presents numerical experiments to assess the computational accuracy and efficiency of the proposed algorithms in Section 3. The objective is to evaluate the performance of the two algorithms under different physical parameter settings. A uniform triangular partition, $\mathcal{T}_{h}$ , is employed to generate the finite element mesh over the domain $\Omega=[0,1]\times[0,1]$ , with the terminal time set to $\tau=0.01$ . The time interval is defined as $J=[0,\tau]$ . The experiments are conducted using the following benchmark problem. All computations are carried out using the FreeFEM++ software [15].

Example 1: We appropriately select the body force $\bm{f}$ , mass source $g$ , and heat source $h$ to ensure that the exact solution is

		$\displaystyle\bm{u}(x,y,t)=e^{-t}\left(\begin{array}[]{c}\sin(2\pi y)(\cos(2\pi x)-1)+\frac{1}{\mu+\lambda}\sin(\pi x)\sin(\pi y)\\ \sin(2\pi x)(1-\cos(2\pi y))+\frac{1}{\mu+\lambda}\sin(\pi x)\sin(\pi y)\\ \end{array}\right),$
		$\displaystyle p(x,y,t)=e^{-t}\sin(\pi x)\sin(\pi y),$
		$\displaystyle T(x,y,t)=e^{-t}\sin(\pi x)\sin(\pi y).$

The corresponding boundary conditions and initial conditions are given as:

\bm{u}(x,y,t)=\bm{0},~{}p(x,y,t)=0,~{}T(x,y,t)=0,\quad\mbox{on}~{}\partial\Omega,

and

		$\displaystyle\bm{u}(x,y,0)=\left(\begin{array}[]{c}\sin(2\pi y)(\cos(2\pi x)-1)+\frac{1}{\mu+\lambda}\sin(\pi x)\sin(\pi y)\\ \sin(2\pi x)(1-\cos(2\pi y))+\frac{1}{\mu+\lambda}\sin(\pi x)\sin(\pi y)\\ \end{array}\right),$
		$\displaystyle p(x,y,0)=\sin(\pi x)\sin(\pi y),$
		$\displaystyle T(x,y,0)=\sin(\pi x)\sin(\pi y).$

In the experiments, we use uniform triangular meshes with sizes $h=1/16$ , $h=1/32$ , $h=1/64$ , and $h=1/128$ . These meshes are refined by bisecting each triangle through its midpoints. At the terminal time $\tau$ , we report the computed $L^{2}$ -norm errors for $\xi$ and $H^{1}$ -norm errors for $\bm{u}$ , $p$ , and $T$ , along with their corresponding convergence rates. The term ”iters” denotes the number of iterations used in the iterative decoupling algorithm. Larger time step sizes are intentionally selected in the tests to emphasize the decoupled algorithm’s efficiency and robustness while demonstrating its stability in the absence of restrictive constraints.

For the coupled algorithm (Alg. 1), the time step size is set to $\Delta t=10^{-3}$ , requiring 10 steps to reach the terminal time $\tau$ . This means that the coupled system (3.2a)–(3.2d) is solved 10 times. For the decoupled algorithm (Alg. 2), two scenarios are considered: $\Delta t=5\times 10^{-3}$ and $\Delta t=10^{-2}$ , corresponding to 2 and 1 time steps, respectively, to reach $\tau$ . In these cases, the number of iterations is set to 5 and 10, respectively. Notably, the total number of system solves for (3.3)–(3.4) matches that of Alg. 1, ensuring an equivalent computational cost. This balance is achieved by appropriately selecting $\Delta t$ values and iteration counts. Despite this equivalence in system solves, the computational time for the decoupled algorithm is significantly reduced because solving smaller, independent subproblems is more efficient than solving the coupled system.

Three key observations emerge from the results: 1. Both the coupled and decoupled algorithms achieve optimal convergence rates in energy norm errors, as shown in Tables LABEL:tab:K_Theta_0.1_nu_0.3 through LABEL:tab:a0=b0=c0=0. The decoupled algorithm matches the accuracy of the coupled algorithm while requiring less computational time. 2. The last entries in Tables LABEL:tab:K_Theta_0.1_nu_0.3 through LABEL:tab:a0=b0=c0=0 reveal that increasing the number of iterations from 5 to 10 in the decoupled algorithm leads to slightly lower errors and higher convergence rates. This demonstrates that the algorithm’s accuracy improves with additional iterations. 3. Both algorithms maintain their effectiveness under extreme parameter values, as demonstrated in the subsequent subsections.

These results validate the computational efficiency and stability of the proposed algorithms across varying scenarios.

5.1 Tests for the parameter $\nu$

Table 1: Example 1: The errors and convergence orders with decreasing mesh sizes of coupled and iterative decoupling algorithms.

\nu=0.3

$h$	$\\|\bm{u}\\|_{H^{1}}$	rate	$\\|\xi\\|_{L^{2}}$	rate	$\\|p\\|_{H^{1}}$	rate	$\\|T\\|_{H^{1}}$	rate
Alg.1, $\Delta t=10^{-3}$
$1/16$	1.01028e-01	0.00	5.96583e-03	0.00	2.29827e-01	0.00	2.29827e-01	0.00
$1/32$	2.54166e-02	1.99	1.47953e-03	2.01	1.09760e-01	1.07	1.09760e-01	1.07
$1/64$	6.36417e-03	2.00	3.69057e-04	2.00	5.42034e-02	1.02	5.42034e-02	1.02
$1/128$	1.59166e-03	2.00	9.20381e-05	2.00	2.70161e-02	1.00	2.70161e-02	1.00
Alg.2, $\Delta t=5\times 10^{-3}$ ,iters=5
$1/16$	1.01030e-01	0.00	5.98023e-03	0.00	2.31732e-01	0.00	2.31732e-01	0.00
$1/32$	2.54169e-02	1.99	1.48305e-03	2.01	1.10015e-01	1.07	1.10015e-01	1.07
$1/64$	6.36426e-03	2.00	3.69788e-04	2.00	5.42370e-02	1.02	5.42370e-02	1.02
$1/128$	1.59168e-03	2.00	9.21223e-05	2.01	2.70211e-02	1.01	2.70211e-02	1.01
Alg.2, $\Delta t=10^{-2}$ ,iters=10
$1/16$	1.01031e-01	0.00	5.99684e-03	0.00	2.34369e-01	0.00	2.34369e-01	0.00
$1/32$	2.54172e-02	1.99	1.48591e-03	2.01	1.10312e-01	1.09	1.10312e-01	1.09
$1/64$	6.36419e-03	2.00	3.69110e-04	2.01	5.42412e-02	1.02	5.42412e-02	1.02
$1/128$	1.59156e-03	2.00	9.07679e-05	2.02	2.70075e-02	1.01	2.70075e-02	1.01

Table 2: Example 1: The errors and convergence orders with decreasing mesh sizes of coupled and iterative decoupling algorithms.

\nu=0.49999

$h$	$\\|\bm{u}\\|_{H^{1}}$	rate	$\\|\xi\\|_{L^{2}}$	rate	$\\|p\\|_{H^{1}}$	rate	$\\|T\\|_{H^{1}}$	rate
Alg.1, $\Delta t=10^{-3}$
$1/16$	1.00295e-01	0.00	9.63830e-03	0.00	2.15490e-01	0.00	2.15490e-01	0.00
$1/32$	2.52279e-02	1.99	2.38247e-03	2.02	1.07906e-01	1.00	1.07906e-01	1.00
$1/64$	6.31671e-03	2.00	5.94250e-04	2.00	5.39732e-02	1.00	5.39732e-02	1.00
$1/128$	1.57979e-03	2.00	1.48471e-04	2.00	2.69891e-02	1.00	2.69891e-02	1.00
Alg.2, $\Delta t=5\times 10^{-3}$ , iters=5
$1/16$	1.00295e-01	0.00	9.63812e-03	0.00	2.15492e-01	0.00	2.15492e-01	0.00
$1/32$	2.52279e-02	1.99	2.38236e-03	2.02	1.07907e-01	1.00	1.07907e-01	1.00
$1/64$	6.31671e-03	2.00	5.94164e-04	2.00	5.39735e-02	1.00	5.39735e-02	1.00
$1/128$	1.57979e-03	2.00	1.48397e-04	2.00	2.69893e-02	1.00	2.69893e-02	1.00
Alg.2, $\Delta t=10^{-2}$ , iters=10
$1/16$	1.00295e-01	0.00	9.63791e-03	0.00	2.15495e-01	0.00	2.15495e-01	0.00
$1/32$	2.52279e-02	1.99	2.38224e-03	2.02	1.07907e-01	1.00	1.07907e-01	1.00
$1/64$	6.31671e-03	2.00	5.94068e-04	2.00	5.39738e-02	1.00	5.39738e-02	1.00
$1/128$	1.57979e-03	2.00	1.48328e-04	2.00	2.69895e-02	1.00	2.69895e-02	1.00

In this subsection, we evaluate the performance of the algorithms under various settings of the Poisson ratio $\nu$ . For the other parameters, we fix $E=1$ , $c_{0}=a_{0}=0.2$ , $b_{0}=0.1$ , $\alpha=\beta=0.1$ , and $\bm{K}=\bm{\Theta}=0.1\bm{I}$ , where $\bm{I}$ denotes the $2\times 2$ identity matrix. The data shown in Table LABEL:tab:K_Theta_0.1_nu_0.3, along with the subsequent tables LABEL:tab:_K_1e-6-LABEL:tab:a0=b0=c0=0, are based on a setting where $\nu=0.3$ . These tables represent the situation where the thermo-poroelastic material is compressible.

In Table LABEL:tab:_nu_0.49999, we set a constant Poisson ratio of $\nu=0.49999$ , while keeping all other physical parameters fixed. Since the Poisson ratio $\nu$ is close to $0.5$ which in other words means $\lambda\to\infty$ , the thermo-poroelastic material behaves nearly incompressibly, leading the mixed linear elasticity model to approach the incompressible Stokes model. The results presented in Table LABEL:tab:_nu_0.49999 demonstrate that the energy-norm errors for all algorithms exhibit optimal orders as the material approaches incompressibility, thereby emphasizing the robustness of these algorithms with respect to the parameter $\nu$ .

5.2 Tests for the parameters $K$ and $\Theta$

Table 3: Example 1: The errors and convergence orders with decreasing mesh sizes of coupled and iterative decoupling algorithms.

\bm{K}=10^{-6}\bm{I}

$h$	$\\|\bm{u}\\|_{H^{1}}$	rate	$\\|\xi\\|_{L^{2}}$	rate	$\\|p\\|_{H^{1}}$	rate	$\\|T\\|_{H^{1}}$	rate
Alg.1, $\Delta t=10^{-3}$
$1/16$	1.01039e-01	0.00	6.07933e-03	0.00	2.66261e-01	0.00	2.35294e-01	0.00
$1/32$	2.54193e-02	1.99	1.50919e-03	2.01	1.18339e-01	1.17	1.10485e-01	1.09
$1/64$	6.36485e-03	2.00	3.76523e-04	2.00	5.62194e-02	1.07	5.42944e-02	1.02
$1/128$	1.59182e-03	2.00	9.38336e-05	2.00	2.74790e-02	1.03	2.70270e-02	1.01
Alg.2, $\Delta t=5\times 10^{-3}$ , iters=5
$1/16$	1.01040e-01	0.00	6.08908e-03	0.00	2.68502e-01	0.00	2.37202e-01	0.00
$1/32$	2.54196e-02	1.99	1.51189e-03	2.01	1.18668e-01	1.18	1.10762e-01	1.10
$1/64$	6.36496e-03	2.00	3.77393e-04	2.00	5.62122e-02	1.08	5.43382e-02	1.03
$1/128$	1.59189e-03	2.00	9.43212e-05	2.00	2.74213e-02	1.04	2.70369e-02	1.01
Alg.2, $\Delta t=10^{-2}$ , iters=10
$1/16$	1.01041e-01	0.00	6.10036e-03	0.00	2.71376e-01	0.00	2.39870e-01	0.00
$1/32$	2.54196e-02	1.99	1.51238e-03	2.01	1.18975e-01	1.19	1.11050e-01	1.11
$1/64$	6.36474e-03	2.00	3.75044e-04	2.01	5.60370e-02	1.09	5.43245e-02	1.03
$1/128$	1.59163e-03	2.00	9.15599e-05	2.03	2.72366e-02	1.04	2.70131e-02	1.01

Table 4: Example 1: The errors and convergence orders with decreasing mesh sizes of coupled and iterative decoupling algorithms.

\bm{\Theta}=10^{-6}\bm{I}

$h$	$\\|\bm{u}\\|_{H^{1}}$	rate	$\\|\xi\\|_{L^{2}}$	rate	$\\|p\\|_{H^{1}}$	rate	$\\|T\\|_{H^{1}}$	rate
Alg.1, $\Delta t=10^{-3}$
$1/16$	1.01039e-01	0.00	6.07933e-03	0.00	2.35294e-01	0.00	2.66261e-01	0.00
$1/32$	2.54193e-02	1.99	1.50919e-03	2.01	1.10485e-01	1.09	1.18339e-01	1.17
$1/64$	6.36485e-03	2.00	3.76523e-04	2.00	5.42944e-02	1.02	5.62194e-02	1.07
$1/128$	1.59182e-03	2.00	9.38336e-05	2.00	2.70270e-02	1.01	2.74790e-02	1.03
Alg.2, $\Delta t=5\times 10^{-3}$ , iters=5
$1/16$	1.01040e-01	0.00	6.08908e-03	0.00	2.37202e-01	0.00	2.68502e-01	0.00
$1/32$	2.54196e-02	1.99	1.51189e-03	2.01	1.10762e-01	1.10	1.18668e-01	1.18
$1/64$	6.36496e-03	2.00	3.77393e-04	2.00	5.43382e-02	1.03	5.62122e-02	1.08
$1/128$	1.59189e-03	2.00	9.43212e-05	2.00	2.70369e-02	1.01	2.74213e-02	1.04
Alg.2, $\Delta t=10^{-2}$ , iters=10
$1/16$	1.01041e-01	0.00	6.10036e-03	0.00	2.39870e-01	0.00	2.71376e-01	0.00
$1/32$	2.54196e-02	1.99	1.51238e-03	2.01	1.11050e-01	1.11	1.18975e-01	1.19
$1/64$	6.36474e-03	2.00	3.75044e-04	2.01	5.43245e-02	1.03	5.60370e-02	1.09
$1/128$	1.59163e-03	2.00	9.15599e-05	2.03	2.70131e-02	1.01	2.72366e-02	1.04

Table 5: Example 1: The errors and convergence orders with decreasing mesh sizes of coupled and iterative decoupling algorithms.

\bm{K}=10^{-6}\bm{I},\bm{\Theta}=10^{-6}\bm{I}

$h$	$\\|\bm{u}\\|_{H^{1}}$	rate	$\\|\xi\\|_{L^{2}}$	rate	$\\|p\\|_{H^{1}}$	rate	$\\|T\\|_{H^{1}}$	rate
Alg.1, $\Delta t=10^{-3}$
$1/16$	1.01060e-01	0.00	6.30251e-03	0.00	3.13011e-01	0.00	3.13011e-01	0.00
$1/32$	2.54251e-02	1.99	1.57062e-03	2.00	1.32705e-01	1.24	1.32705e-01	1.24
$1/64$	6.36635e-03	2.00	3.92445e-04	2.00	6.01621e-02	1.14	6.01621e-02	1.14
$1/128$	1.59219e-03	2.00	9.77413e-05	2.01	2.84617e-02	1.08	2.84617e-02	1.08
Alg.2, $\Delta t=5\times 10^{-3}$ , iters=5
$1/16$	1.01060e-01	0.00	6.29942e-03	0.00	3.12563e-01	0.00	3.12563e-01	0.00
$1/32$	2.54252e-02	1.99	1.57117e-03	2.00	1.32870e-01	1.23	1.32870e-01	1.23
$1/64$	6.36653e-03	2.00	3.94108e-04	2.00	6.05260e-02	1.13	6.05260e-02	1.13
$1/128$	1.59241e-03	2.00	9.99148e-05	1.98	2.89124e-02	1.07	2.89124e-02	1.07
Alg.2, $\Delta t=10^{-2}$ , iters=10
$1/16$	1.01060e-01	0.00	6.29765e-03	0.00	3.12246e-01	0.00	3.12246e-01	0.00
$1/32$	2.54246e-02	1.99	1.56551e-03	2.01	1.31875e-01	1.24	1.31875e-01	1.24
$1/64$	6.36589e-03	2.00	3.87375e-04	2.01	5.93113e-02	1.15	5.93113e-02	1.15
$1/128$	1.59179e-03	2.00	9.32757e-05	2.05	2.77070e-02	1.10	2.77070e-02	1.10

In this subsection, we investigate the accuracy of the algorithms under varying settings of hydraulic conductivity $\bm{K}$ and effective thermal conductivity $\bm{\Theta}$ . For reference, we have already tested the case $\bm{K}=\bm{\Theta}=0.1\bm{I}$ in Table LABEL:tab:K_Theta_0.1_nu_0.3 in previous experiments. In the current tests (Table LABEL:tab:_K_1e-6-Table LABEL:tab:K_Theta_1e-6), we set $(\bm{K},\bm{\Theta})=(10^{-6}\bm{I},0.1\bm{I})$ , $(\bm{K},\bm{\Theta})=(0.1\bm{I},10^{-6}\bm{I})$ , and $(\bm{K},\bm{\Theta})=(10^{-6}\bm{I},10^{-6}\bm{I})$ , respectively, to observe the impact of significantly smaller conductivity values on the solution accuracy. The remaining key parameters are fixed as $E=1$ , $c_{0}=a_{0}=0.2$ , $b_{0}=0.1$ , $\alpha=\beta=0.1$ , and $\nu=0.3$ . We present the numerical results obtained using Alg. 1 (the coupled algorithm) and Alg. 2 (the iterative decoupling algorithm), with varying numbers of iterations, in Tables LABEL:tab:_K_1e-6 through LABEL:tab:K_Theta_1e-6.

By comparing the results from Tables LABEL:tab:_K_1e-6 and LABEL:tab:K_Theta_1e-6 with those in Table LABEL:tab:K_Theta_0.1_nu_0.3, we observe that all algorithms continue to produce solutions with optimal convergence rates. This indicates that the algorithms remain effective and robust, even when hydraulic and thermal conductivities are significantly reduced, further validating their performance under a wide range of physical conditions.

5.3 Tests for the parameter $a_{0},b_{0},c_{0}$

In this subsection, we investigate the influence of the effective thermal parameter $a_{0}$ , the thermal dilation coefficient $b_{0}$ , and the specific storage coefficient $c_{0}$ on the accuracy of the algorithms. As discussed in Remark 1, the case $c_{0}-b_{0}=a_{0}-b_{0}=0$ may affect the convergence rate of the iterative decoupling algorithms. To evaluate this, Table LABEL:tab:a0=b0=c0=0 presents the limiting scenario where the coefficients $a_{0}$ , $b_{0}$ , and $c_{0}$ are all set to zero. For the remaining parameters, we fix $E=1$ , $\alpha=\beta=0.1$ , and $\bm{K}=\bm{\Theta}=0.1\bm{I}$ .

Table 6: Example 1: The errors and convergence orders with decreasing mesh sizes of coupled and iterative decoupling algorithms.

a_{0}=b_{0}=c_{0}=0

$h$	$\\|\bm{u}\\|_{H^{1}}$	rate	$\\|\xi\\|_{L^{2}}$	rate	$\\|p\\|_{H^{1}}$	rate	$\\|T\\|_{H^{1}}$	rate
Alg.1, $\Delta t=10^{-3}$
$1/16$	1.01017e-01	0.00	5.86499e-03	0.00	2.18610e-01	0.00	2.18610e-01	0.00
$1/32$	2.54138e-02	1.99	1.45400e-03	2.01	1.08332e-01	1.01	1.08332e-01	1.01
$1/64$	6.36347e-03	2.00	3.62641e-04	2.00	5.40235e-02	1.00	5.40235e-02	1.00
$1/128$	1.59148e-03	2.00	9.04325e-05	2.00	2.69935e-02	1.00	2.69935e-02	1.00
Alg.2, $\Delta t=5\times 10^{-3}$ , iters=5
$1/16$	1.01050e-01	0.00	6.19651e-03	0.00	2.48520e-01	0.00	2.48520e-01	0.00
$1/32$	2.54262e-02	1.99	1.57535e-03	1.98	1.14200e-01	1.12	1.14200e-01	1.12
$1/64$	6.37219e-03	2.00	4.59498e-04	1.78	5.61495e-02	1.02	5.61495e-02	1.02
$1/128$	1.60655e-03	1.99	2.39721e-04	0.94	2.88472e-02	0.96	2.88472e-02	0.96
Alg.2, $\Delta t=10^{-2}$ , iters=10
$1/16$	1.01126e-01	0.00	6.95719e-03	0.00	3.08253e-01	0.00	3.08253e-01	0.00
$1/32$	2.54412e-02	1.99	1.73121e-03	2.01	1.21121e-01	1.35	1.21121e-01	1.35
$1/64$	6.36954e-03	2.00	4.24246e-04	2.03	5.55013e-02	1.13	5.55013e-02	1.13
$1/128$	1.59228e-03	2.00	9.84554e-05	2.11	2.70932e-02	1.03	2.70932e-02	1.03

A comparison of the second entry in Table LABEL:tab:a0=b0=c0=0 with the corresponding entry in Table LABEL:tab:K_Theta_0.1_nu_0.3 reveals that setting $a_{0}=b_{0}=c_{0}=0$ leads to a slight degradation in the errors and convergence rates for the iterative decoupling algorithm (Alg. 2). However, when the number of iterations increases to $10$ , errors and rates reach optimal levels, as shown in the third entry of Table LABEL:tab:a0=b0=c0=0. These results confirm that the energy norm errors for all algorithms maintain optimal convergence orders, underscoring their robustness with respect to the parameters $a_{0}$ , $b_{0}$ , and $c_{0}$ .

6 Conclusions

In this paper, we propose a novel algorithm designed to decouple the computations of the linear thermo-poroelastic model, based on a newly developed four-field formulation. The algorithm effectively separates the numerical computations into two sub-models, streamlining the overall computational process. We perform a comprehensive error analysis to establish the algorithm’s unconditional stability and optimal convergence. Furthermore, numerical experiments are performed to validate the theoretical error estimates. The results demonstrate that the proposed algorithm is unconditionally stable, computationally efficient, and achieves optimal convergence rates even under widely varying parameter settings or in the limiting cases of parameters, highlighting its robustness and reliability for addressing thermo-poroelastic problems.

Acknowledgments

The work of M. Cai is partially supported by the NIH-RCMI award (Grant No. 347U54MD013376) and the affiliated project award from the Center for Equitable Artificial Intelligence and Machine Learning Systems (CEAMLS) at Morgan State University (project ID 02232301). The work of J. Li and Z. Li are partially supported by the Shenzhen Sci-Tech Fund No. RCJC20200714114556020, Guangdong Basic and Applied Research Fund No. 2023B1515250005. The work of Q. Liu is partially supported by the NSFC (12271186, 12271178), the Guangdong Basic and Applied Basic Research Foundation (2022B1515120009), and the Science and Technology Program of Shenzhen, China(20231121110406001).

References

[1] Paola F Antonietti, Stefano Bonetti, and Michele Botti. Discontinuous Galerkin approximation of the fully coupled thermo-poroelastic problem. SIAM Journal on Scientific Computing, 45(2):A621–A645, 2023.
[2] Francesco Ballarin, Sanghyun Lee, and Son-Young Yi. Projection-based reduced order modeling of an iterative scheme for linear thermo-poroelasticity. Results in Applied Mathematics, 21:100430, 2024.
[3] Maurice A Biot. General theory of three-dimensional consolidation. Journal of applied physics, 12(2):155–164, 1941.
[4] Maurice A Biot. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of applied physics, 26(2):182–185, 1955.
[5] Stefano Bonetti, Michele Botti, and Paola F Antonietti. Robust discontinuous Galerkin-based scheme for the fully-coupled nonlinear thermo-hydro-mechanical problem. IMA Journal of Numerical Analysis, page drae045, 2024.
[6] Susanne C Brenner. Korn’s inequalities for piecewise $H^{1}$ vector fields. Mathematics of Computation, pages 1067–1087, 2004.
[7] Mats K Brun, Inga Berre, Jan M Nordbotten, and Florin A Radu. Upscaling of the coupling of hydromechanical and thermal processes in a quasi-static poroelastic medium. Transport in Porous Media, 124:137–158, 2018.
[8] Mats Kirkesæther Brun, Elyes Ahmed, Inga Berre, Jan Martin Nordbotten, and Florin Adrian Radu. Monolithic and splitting solution schemes for fully coupled quasi-static thermo-poroelasticity with nonlinear convective transport. Computers & Mathematics with Applications, 80(8):1964–1984, 2020.
[9] Mats Kirkesæther Brun, Elyes Ahmed, Jan Martin Nordbotten, and Florin Adrian Radu. Well-posedness of the fully coupled quasi-static thermo-poroelastic equations with nonlinear convective transport. Journal of Mathematical Analysis and Applications, 471(1-2):239–266, 2019.
[10] Yuxiang Chen and Zhihao Ge. Multiphysics finite element method for quasi-static thermo-poroelasticity. Journal of Scientific Computing, 92(2):43, 2022.
[11] Philippe G Ciarlet. The finite element method for elliptic problems. SIAM, 2002.
[12] Xiaobing Feng, Zhihao Ge, and Yukun Li. Analysis of a multiphysics finite element method for a poroelasticity model. IMA Journal of Numerical Analysis, 38(1):330–359, 03 2017.
[13] Zhihao Ge and Dandan Xu. Analysis of multiphysics finite element method for quasi-static thermo-poroelasticity with a nonlinear convective transport term. arXiv preprint arXiv:2310.05084, 2023.
[14] H. Gu, M. Cai, Jingzhi Li, and Guoliang Ju. A priori error estimates of two monolithic schemes for biot’s consolidation model. Numerical Methods for Partial Differential Equations, 2023.
[15] Frédéric Hecht. New development in freefem++. In J. Num. Math., 2012.
[16] Guoliang Ju, Mingchao Cai, Jingzhi Li, and Jing Tian. Parameter-robust multiphysics algorithms for biot model with application in brain edema simulation. Mathematics and Computers in Simulation (MATCOM), 177, 2020.
[17] Jihoon Kim. Unconditionally stable sequential schemes for all-way coupled thermoporomechanics: Undrained-adiabatic and extended fixed-stress splits. Computer Methods in Applied Mechanics and Engineering, 341:93–112, 2018.
[18] Alexandr E Kolesov, Petr N Vabishchevich, and Maria V Vasilyeva. Splitting schemes for poroelasticity and thermoelasticity problems. Computers & Mathematics with Applications, 67(12):2185–2198, 2014.
[19] Jeonghun J Lee, Kent-Andre Mardal, and Ragnar Winther. Parameter-robust discretization and preconditioning of Biot’s consolidation model. SIAM Journal on Scientific Computing, 39(1):A1–A24, 2017.
[20] Ricardo Oyarzúa and Ricardo Ruiz-Baier. Locking-free finite element methods for poroelasticity. SIAM Journal on Numerical Analysis, 54(5):2951–2973, 2016.
[21] J. Rutqvist, Y. S. Wu, C. F. Tsang, and G. Bodvarsson. A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock. International Journal of Rock Mechanics and Mining Sciences, 39(4):429–442, 2002.
[22] Karl Terzaghi. Theoretical soil mechanics. 1943.
[23] N. Watanabe, B. Zehner, W. Wang, C. I. Mcdermott, and O. Kolditz. Numerical analysis and visualization of uncertainty effects in thermo-hydro-mechanical coupled processes in a hot-dry-rock geothermal reservoir. IAHS-AISH publication, pages 57–62, 2011.
[24] J. L Yow and J. R Hunt. Coupled processes in rock mass performance with emphasis on nuclear waste isolation. International Journal of Rock Mechanics and Mining Sciences, 39(2):143–150, 2002.
[25] Jing Zhang and Hongxing Rui. Galerkin method for the fully coupled quasi-static thermo-poroelastic problem. Computers & Mathematics with Applications, 118:95–109, 2022.
[26] Jing Zhang and Hongxing Rui. The mfe-cfe-gfe method for the fully coupled quasi-static thermo-poroelastic problem. Numer. Math. Theor. Meth. Appl., 16:792–819, 2023.
[27] Jing Zhang and Hongxing Rui. A coupling of Galerkin and mixed finite element methods for the quasi-static thermo-poroelasticity with nonlinear convective transport. Journal of Computational and Applied Mathematics, 441:115672, 2024.

$\displaystyle\zeta_{0}\\|\xi\\|_{L^{2}(\Omega)}=$	$\displaystyle\sup_{\bm{v}\in\bm{V}}\frac{\|(\text{div }\bm{v},\xi)\|}{\\|\bm{v}\\|_{H^{1}(\Omega)}}$	(2.12)
$\displaystyle\leq$	$\displaystyle\sup_{\bm{v}\in\bm{V}}\frac{\|2\mu(\varepsilon(\bm{u}),\varepsilon(\bm{v}))+(\bm{f},\bm{v})\|}{\\|\bm{v}\\|_{H^{1}(\Omega)}}$
$\displaystyle\leq$	$\displaystyle C_{K}\big{(}2\mu\\|\bm{\varepsilon}(u)\\|_{L^{2}(\Omega)}+\\|\bm{f}\\|_{L^{2}(\Omega)}\big{)},$

		$\displaystyle\sqrt{2\mu}\\|\bm{\varepsilon}(\bm{u})\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{\frac{1}{\lambda}}\\|\alpha p+\beta T-\xi\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{c_{0}-b_{0}}\\|p\\|_{L^{\infty}(J;L^{2}(\Omega))}$		(2.13)
		$\displaystyle+\sqrt{a_{0}-b_{0}}\\|T\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{b_{0}}\\|p-T\\|_{L^{\infty}(J;L^{2}(\Omega))}+\sqrt{k_{m}}\\|\nabla p\\|_{L^{2}(J;L^{2}(\Omega))}$
		$\displaystyle+\sqrt{\theta_{m}}\\|\nabla T\\|_{L^{2}(J;L^{2}(\Omega))}\leq C.$

		$\displaystyle\mu\\|\bm{\varepsilon}(\bm{u})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p+\beta T-\xi\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T\\|_{L^{2}(\Omega)}^{2}$		(2.15)
		$\displaystyle+\frac{b_{0}}{2}\\|p-T\\|_{L^{2}(\Omega)}^{2}+\int_{0}^{t}k_{m}\\|\nabla p\\|_{L^{2}(\Omega)}^{2}ds+\int_{0}^{t}\theta_{m}\\|\nabla T\\|_{L^{2}(\Omega)}^{2}ds$
		$\displaystyle\leq\mu\\|\bm{\varepsilon}(\bm{u}^{0})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p^{0}+\beta T^{0}-\xi^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T^{0}\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\frac{b_{0}}{2}\\|p^{0}-T^{0}\\|_{L^{2}(\Omega)}^{2}+2\mu^{}\\|\bm{f}\\|_{L^{2}(\Omega)}^{2}+\frac{1}{4\mu^{}}\\|\bm{u}\\|_{L^{2}(\Omega)}^{2}+\frac{1}{4\mu^{*}}\\|\bm{u}^{0}\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\frac{1}{2}\int_{0}^{t}(\\|g\\|_{L^{2}(\Omega)}^{2}+\\|H_{s}\\|_{L^{2}(\Omega)}^{2}+\\|p\\|_{L^{2}(\Omega)}^{2}+\\|T\\|_{L^{2}(\Omega)}^{2}).$

		$\displaystyle\mu\\|\bm{\varepsilon}(\bm{u})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p+\beta T-\xi\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T\\|_{L^{2}(\Omega)}^{2}$		(2.16)
		$\displaystyle+\frac{b_{0}}{2}\\|p-T\\|_{L^{2}(\Omega)}^{2}+\int_{0}^{t}k_{m}\\|\nabla p\\|_{L^{2}(\Omega)}^{2}ds+\int_{0}^{t}\theta_{m}\\|\nabla T\\|_{L^{2}(\Omega)}^{2}ds$
		$\displaystyle\lesssim\frac{3\mu}{2}\mu\\|\bm{\varepsilon}(\bm{u}^{0})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\lambda}\\|\alpha p^{0}+\beta T^{0}-\xi^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{c_{0}-b_{0}}{2}\\|p^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{a_{0}-b_{0}}{2}\\|T^{0}\\|_{L^{2}(\Omega)}^{2}$
		$\displaystyle+\frac{b_{0}}{2}\\|p^{0}-T^{0}\\|_{L^{2}(\Omega)}^{2}+\frac{1}{\mu}\\|\bm{f}\\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{2}\\|\bm{\varepsilon}(\bm{u})\\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\int_{0}^{t}(\\|g\\|_{L^{2}(\Omega)}^{2}+\\|H_{s}\\|_{L^{2}(\Omega)}^{2}),$

		$\displaystyle\sqrt{2\mu}\\|\bm{\varepsilon}(\bm{u}_{t})\\|_{L^{2}(J;L^{2}(\Omega))}+\sqrt{\frac{1}{\lambda}}\\|\alpha p_{t}+\beta T_{t}-\xi_{t}\\|_{L^{2}(J;L^{2}(\Omega))}+\sqrt{c_{0}-b_{0}}\\|p_{t}\\|_{L^{2}(J;L^{2}(\Omega))}$		(2.17)
		$\displaystyle+\sqrt{a_{0}-b_{0}}\\|T_{t}\\|_{L^{2}(J;L^{2}(\Omega))}+\sqrt{b_{0}}\\|p_{t}-T_{t}\\|_{L^{2}(J;L^{2}(\Omega))}+\sqrt{k_{m}}\\|\nabla p_{t}\\|_{L^{\infty}(J;L^{2}(\Omega))}$
		$\displaystyle\sqrt{\theta_{m}}\\|\nabla T_{t}\\|_{L^{\infty}(J;L^{2}(\Omega))}\leq C_{1},$

Efficient Iterative Decoupling Methods for Thermo-Poroelasticity Based on a Four-Field Formulation

Abstract

1 Introduction

2 Mathematical model and its reformulation

2.1 A linear thermo-poroelasticity model

2.2 A priori estimates of the weak solution

Lemma 2.1.

Lemma 2.2.

Proof.

Theorem 2.1.

Proof.

Lemma 2.3.

Proof.

Theorem 2.2.

3 Numerical algorithms

3.1 A coupled algorithm

3.2 An iteratively decoupled algorithm

4 Convergence analysis

Lemma 4.1.

Theorem 4.1.

Proof.

Remark 1.

Lemma 4.2.

Proof.

Lemma 4.3.

Proof.

Theorem 4.2.

Proof.

Theorem 4.3.

Proof.

5 Numerical experiments

5.1 Tests for the parameter ν\nu

5.2 Tests for the parameters KK and Θ\Theta

5.3 Tests for the parameter a0,b0,c0a_{0},b_{0},c_{0}

6 Conclusions

Acknowledgments

References

5.1 Tests for the parameter $\nu$

5.2 Tests for the parameters $K$ and $\Theta$

5.3 Tests for the parameter $a_{0},b_{0},c_{0}$