\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}



\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]





%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}


\def\udcs{519.6} %Here the author places classificators of the paper according to Russian classification system
\def\mscs{26A33,65M60} %Here the author places  classificators according to the AMS classification list.
\setcounter{page}{144}

\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}



\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Vol. 20, No. 2, pp. 144--145 (2023) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline  %Volume, pages in Russian, filled by Editorial board
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2023.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}




\begin{document}
\renewcommand{\refname}{References}
\renewcommand{\proofname}{Proof.}
\renewcommand{\figurename}{Fig.}




\thispagestyle{empty}

\title[\sc Online GMsFEM for Poroelasticity pr.]{\large Numerical simulation of the time fractional thermoporoelasticity problem using online generalized multiscale finite element method}

\author[A.A. TYRYLGIN, U.S. KALACHIKOVA, H.BAI, A.A. ALIKHANOV]
{
{\bf \large \textsc{A.A. Tyrylgin}}\protect\orcidicon{0000-0002-2157-0921} 
{\bf \large \textsc{, U.S. Kalachikova}}\protect\orcidicon{0000-0002-1561-449X} \\  
{\bf \large \textsc{H. Bai}}\protect\orcidicon{0009-0001-5048-5339} 
{\bf \large \textsc{, A.A. Alihkanov}}\protect\orcidicon{0000-0003-0684-6667}\\ \\ 
}% 

\address{Tyrylgin Aleksei Afanasievich 
\newline\hphantom{iii} Laboratory of Computational Technologies and Artificial Intelligence, 
\newline\hphantom{iii} North-Eastern Federal University, Belinsky st., 58,
\newline\hphantom{iii} 677000, Yakutsk, Russia%
\newline\hphantom{iii} North-Caucasus Center for Mathematical Research, 
\newline\hphantom{iii} North-Caucasus Federal University, Pushkina st., 1,
\newline\hphantom{iii} 355017, Stavropol, Russia}

\email{\textcolor{blue}{aa.tyrylgin@mail.ru}}

\address{Kalachikova Uigulaana Semenovna 
\newline\hphantom{iii} Laboratory of Computational Technologies and Artificial Intelligence, 
\newline\hphantom{iii} North-Eastern Federal University, Belinsky st., 58,
\newline\hphantom{iii} 677000, Yakutsk, Russia}%
\email{\textcolor{blue}{lanasemna@mail.ru}}

\address{Huiran Bai % ?????????? ?????? ???? ??????? ? ????? ?????? ??????????? ?????? ?? ?????????? ?????
\newline\hphantom{iii} National Center for Applied Mathematics in Hunan and Hunan Key Laboratory for Computation and Simulation in Science and Engineering, 
\newline\hphantom{iii} School of Mathematics and Computational Science, Xiangtan University,
\newline\hphantom{iii} 411105, Hunan, China}%
\email{\textcolor{blue}{202331510101@smail.xtu.edu.cn}}%

\address{Anatoly Alikhanov % ?????????? ?????? ???? ??????? ? ????? ?????? ??????????? ?????? ?? ?????????? ?????
\newline\hphantom{iii} North-Caucasus Center for Mathematical Research, 
\newline\hphantom{iii} North-Caucasus Federal University, Pushkina st., 1,
\newline\hphantom{iii} 355017, Stavropol, Russia}%
\email{\textcolor{blue}{aaalikhanov@gmail.com}}%

\thanks{\sc Tyrylgin A.A., Kalachikova U.S., Bai H., Alikhanov A.A.,
Numerical simulation of the time fractional thermoporoelasticity problem using online generalized multiscale finite element method}
\thanks{\copyright \ Tyrylgin A.A., Kalachikova U.S., Bai H., Alikhanov A.A.}
\thanks{\rm The work is supported by the grant of Russian Science Foundation no. 23-41-00037(https://rscf.ru/project/23-41-00037/) and North-Caucasus Center for Mathematical Research under agreement no. 075-02-2024-1451 with the Ministry of Science and Higher Education of the Russian Federation.}
\thanks{\it Received September, 30, 2024, Published September, 30, 2024}%


\semrtop \vspace{1cm}
\maketitle {\small

\vspace{-12pt}
\centerline{{\it Communicated by} {\sc P.P. Petrov}
%({\it filled by Editorial board})
}

\bigskip
\bigskip

\begin{quote}
\noindent{\bf Abstract:} This paper considers the time of fractional derivatives thermoporoelasticity in a fracture and heterogeneous media. The mathematical model is formulated as a related system of equations that regulate the pressure of the fluid, temperature and mechanical displacement. We use the finite element method (FEM) with a fine mesh  for spatial sampling integrated with a discrete fracture model (DFM) to capture the complexity of fractures in the heterogeneous media. Temporary sampling is achieved using an implicit scheme of final differences. To increase the effectiveness of computing technology, we use the Generalized Multiscale Finite Element Method (GMSFEM) for coarse grid approximation, effectively reducing the dimension of the problem while maintaining accuracy. A multiscale approach uses precalculated offline functions and dynamically updated online functions to process local residues, ensuring a rapid decrease in errors. Numerical experiments demonstrate the ability of the method to accurately imitate the temporal processes in a fractured porous media, reaching significant computing savings without lowering the accuracy of the solution.\medskip

\noindent{\bf Keywords:} thermoporoelasticity, fractional derivatives, multiscale methods, discrete fracture model, heterogeneous media, finite difference scheme.  
 \end{quote}
 
 \begin{quote}
\noindent{\bf Аннотация:} В данной работе рассматривается дробная производная по времени для задачи термопороупругости в трещиноватых и неоднородных средах. Математическая модель сформулирована как связанная система уравнений, которые включают давление жидкости, температуру и механическое перемещение. Мы используем метод конечных элементов (FEM) с мелкой сеткой для пространственной выборки, интегрированный с дискретной моделью трещины (DFM), чтобы охватить сложность трещин в неоднородных средах. Выборка по времени достигается с помощью неявной схемы конечных разностей. Для повышения эффективности вычислительной технологии мы используем обобщенный многомасштабный метод конечных элементов (GMSFEM) для грубой сеточной аппроксимации, эффективно уменьшая размерность задачи при сохранении точности. Многомасштабный подход использует предварительно вычисленные оффлайн функции и динамически обновляемые онлайн-функции для обработки локальной невязки, обеспечивая быстрое снижение погрешности. Численные эксперименты демонстрируют способность метода точно имитировать временные процессы в трещиноватой пористой среде, достигая значительной экономии вычислительных ресурсов без снижения точности решения.\medskip

\noindent{\bf Ключевые слова:} термопороупругость, дробные производные, многомасштабные методы, дискретная модель трещины, неоднородные среды, конечно-разностная схема. 
 \end{quote}
}




\bigskip

\section{Introduction}
Thermoporoelasticity, a field that blends thermal, hydraulic, and mechanical interactions within porous materials, plays a key role in understanding subsurface environments like geothermal systems, hydrocarbon reservoirs, and sites used for environmental purposes \cite{1,2,3,4}. It focuses on the intricate dynamics between fluid movement, heat conduction, and mechanical shifts within porous networks, which are frequently characterized by heterogeneity and the presence of fractures.

The inclusion of fractional derivatives in the mathematical models of these processes has garnered considerable attention because of their ability to represent memory effects and long-term temporal dependencies, which are frequently observed in geological formations \cite{5,6,7}. Compared to traditional integer-order derivatives, fractional derivatives provide a more refined description of history-dependent behaviors \cite{8,9,10}. This feature is especially useful for modeling phenomena like anomalous diffusion and stress relaxation, which are crucial in many engineering and geophysical contexts.

In media characterized by heterogeneity and the presence of fractures, these complexities occur across multiple scales. Traditional numerical modeling approaches often require highly detailed discretization, leading to heavy computational demands, especially in large-scale simulations. To address these issues, the Generalized Multiscale Finite Element Method (GMsFEM) has emerged as a robust technique \cite{11,12,13,14,15,16}. GMsFEM enables efficient modeling of both macroscopic behaviors and small-scale heterogeneities by utilizing coarse grid approximations enhanced by locally computed basis functions that capture the essential system features. This method is especially well-suited for fractured porous media, where it is critical to accurately model the interactions between the porous matrix and fractures \cite{17,18,19}.

A significant advancement in this study is the integration of online multiscale methods within the GMsFEM framework \cite{20,21,22}. The online multiscale approach allows for the dynamic updating of basis functions based on residual information during the simulation. This adaptive strategy enhances the accuracy and efficiency of the model by refining the multiscale basis functions as the simulation progresses, focusing computational resources on areas where they are most needed.

In this study, we develop and apply a multiscale model reduction technique for the time fractional thermoporoelasticity problem in fractured and heterogeneous media. By integrating fractional time derivatives and online multiscale methods into the GMsFEM framework, we aim to capture the memory effects and complex interactions between fluid flow, heat transfer, and mechanical deformations in a computationally efficient manner. Additionally, the use of a Discrete Fracture Model (DFM) within this framework enhances the ability to explicitly model fractures, providing a more accurate and detailed simulation of thermoporoelastic processes in fractured media.

The structure of this paper is as follows: Section 2 outlines the mathematical formulation of the thermoporoelasticity problem, incorporating fractional derivatives. Section 3 discusses the fine grid approximation using the finite element method, while Section 4 details the coarse grid approximation and the construction of multiscale basis functions using GMsFEM, including the implementation of online multiscale methods. Numerical results demonstrating the effectiveness of the proposed method are presented in Section 5. Finally, the study concludes with a summary in the concluding section.


\section{Problem formulation}
For the mechanics of the thermoporoelastic multicontinuum media, we have the following coupled system of equations 
\begin{equation}
\label{f1}
\begin{aligned}
c_{m}\frac{\partial^{\alpha_{m}}p_{m}}{\partial t^{\alpha_{m}}}+\gamma_{m}\frac{\partial^{\beta_{m}}\operatorname{div}\mathbf{u}}{\partial t^{\beta_{m}}}-\nabla\cdot(k_{m}\nabla p_{m})+\eta_{mf}(p_{m}-p_{f}) =0, \ \Omega\times(0,T),  \\
c_{f}\frac{\partial^{\alpha_{f}}p_{f}}{\partial t^{\alpha_{f}}}+\gamma_{f}\frac{\partial^{\beta_{f}}\operatorname{div}\mathbf{u}}{\partial t^{\beta_{f}}}-\nabla\cdot(k_{f}\nabla p_{f})+\eta_{mf}(p_{f}-p_{m}) =0, \ \gamma\times(0,T),  \\
s_{m}\frac{\partial^{\zeta_{m}}T_{m}}{\partial t^{\zeta_{m}}}+\delta_{m}\frac{\partial^{\beta_{m}}\operatorname{div}\mathbf{u}}{\partial t^{\beta_{m}}}+C_{w}\mathcal{V}_{m}\cdot\nabla T_{m}-\nabla\cdot(\chi_{m}\nabla T_{m})+  \\
+ \Upsilon_{mf}(T_{m}-T_{f}) =0, \ \Omega\times(0,T), \\
s_{f}\frac{\partial^{\zeta_{f}}T_{f}}{\partial t^{\zeta_{f}}}+\delta_{f}\frac{\partial^{\beta_{f}}\operatorname{div}\mathbf{u}}{\partial t^{\beta_{f}}}+C_{w}\mathcal{V}_{f}\cdot\nabla T_{f}-\nabla\cdot(\chi_{f}\nabla T_{f})+  \\
+ \Upsilon_{mf}(T_{f}-T_{m}) =0, \ \gamma\times(0,T), \\
-\operatorname{div}\boldsymbol{\sigma}(\mathbf{u})+\gamma_{m}\nabla p_{m}+\gamma_{f}\nabla p_{f}+\delta_{m}\nabla T_{m}+\delta_{f}\nabla T_{f} =0, \ \Omega\times(0,T),
\end{aligned}
\end{equation}
where $\boldsymbol{\sigma}$ denotes the stress tensor, $\mathbf{u}$ the displacement, ${\gamma _i}$ the Biot coefficient, ${M_i}$ the Biot modulus ($c_{i}={\frac{1}{M_{i}}}$) for the $i$-th component. In the case of a linear elastic stress-strain constitutive relation, we have
$$
\boldsymbol{\sigma}(\mathbf{u})=2\mu\boldsymbol{\varepsilon}(\mathbf{u})+\lambda\operatorname{div}\mathbf{u}\mathcal{I},\quad\boldsymbol{\varepsilon}(\mathbf{u})=\frac12(\nabla\mathbf{u}+\nabla\mathbf{u}^T),
$$
where $\varepsilon $ is the strain tensor, $\lambda $ and $\mu $ are the Lame’s coefficients. Here we have a volume force source that is proportional to the sum of the pressure gradients for each continuum. The first continuum describes a flow in the matrix of the porous media and the second continuum relates to the flow in low dimensional fracture networks.

\section{Fine grid approximation using FEM}

This section presents a fine grid finite element approximation for the two-dimensional time fractional thermoporoelasticity problem. For the temporal approximation, we use uniform temporal grid with nodes $t^n=n\tau(n=0,1,2,\ldots,N_{T})$, where $T$ denotes the final time for simulation and $\tau=\frac{T}{N_{T}}$ denotes temporal grid size. The values of pressures, temperature and displacement at the temporal notes $t^n=n\tau(n=0,1,2,\ldots,N_{T})$ are denoted by $(p_{i}^{n},T_{i}^{n},\mathbf{u}^{n})=(p_{i}(t^{n}),T_{i}(t^{n}),\mathbf{u}(t^{n}))$, where $p_{i}$ and $T_{i}$ denote the pressure and temperature of the $i$th continuum respectively. We begin with a variational form of the system (\ref{f1}). Then we exhibit the discrete system in matrix form. The stability and convergence of the time derivative difference scheme can be seen in this paper \cite{23}.

\textbf{Variational form.} For the spatial approximation, we use the continuous Galerkin finite element method with linear basis functions. Let us define functional spaces 
$$
W = {\left[ {{H^1}\left( \Omega  \right)} \right]^d},\quad V_m = {H^1}\left( \Omega  \right), \quad {V_f} = {H^1}\left( \gamma  \right),
$$ 
and $p_{m}\in V_m,p_{f}\in V_{f},T_{m}\in V_m,T_{f}\in V_{f},\mathbf{u}\in W.$ 

We multiply the system (\ref{f1}) by test functions
 $v_{m}\in V_m,v_{f}\in V_{f},z_{m}\in V_m,z_{f}\in V_{f},w\in W,$ respectively. The variational form of the time fractional thermoporoelasticity problem in multicontinuum media can be written as follows: given $(p_{m}^{0},p_{f}^{0},T_{m}^{0},T_{f}^{0},\mathbf{u}^{0})\in V_m\times V_{f}\times V_m\times V_{f}\times W$ iteratively find $(p_{m}^{n},p_{f}^{n},T_{m}^{n},T_{f}^{n},\mathbf{u}^{n})\in V_m\times V_{f}\times V_m\times V_{f}\times W$ such that 
$$
\xi_{\tau}^{(\alpha_{i})}m_{i}^{p}(p_{i}^{n}-p_{i}^{n-1},v_{i}) +\xi_{\tau}^{(\alpha_{i})}\sum_{j=2}^{n}\xi_{j-1}^{(\alpha_{i})}m_{i}^{p}(p_{i}^{n-j+1}-p_{i}^{n-j},v_{i}) +
$$
$$
+\xi_{\tau}^{(\beta_{i})}d_{i}^{p}(\mathbf{u}^{n}-\mathbf{u}^{n-1},v_{i})+\xi_{\tau}^{(\beta_{i})}\sum_{j=2}^{n}\xi_{j-1}^{(\beta_{i})}d_{i}^{p}(\mathbf{u}^{n-j+1}-\mathbf{u}^{n-j},v_{i}) +
$$
$$
+b_{i}^{p}(p_{i}^{n},v_{i})+\sum_{j\neq i}q_{ij}^{p}(p_{i}^{n}-p_{j}^{n},v_{i})=0,\quad\forall v_{i}\in V_{i},i=m,f, 
$$
$$
\xi_{\tau}^{(\zeta_{i})}m_{i}^{T}(T_{i}^{n}-T_{i}^{n-1},z_{i}) +\xi_{\tau}^{(\zeta_{i})}\sum_{j=2}^{n}\xi_{j-1}^{(\zeta_{i})}m_{i}^{T}(T_{i}^{n-j+1}-T_{i}^{n-j},z_{i}) +
$$
$$
+\xi_{\tau}^{(\beta_{i})}d_{i}^{T}(\mathbf{u}^{n}-\mathbf{u}^{n-1},z_{i})+\xi_{\tau}^{(\beta_{i})}\sum_{j=2}^{n}\xi_{j-1}^{(\beta_{i})}d_{i}^{T}(\mathbf{u}^{n-j+1}-\mathbf{u}^{n-j},z_{i}) +
$$
$$
+b_{i}^{T}(T_{i}^{n},z_{i})+\sum_{j\neq i}q_{ij}^{T}(T_{i}^{n}-T_{j}^{n},z_{i})=0,\quad\forall z_{i}\in V_{i},i=m,f, 
$$
$$
a(\mathbf{u}^{n},\mathbf{w})+\sum_{j}g_{i}^{p}(p_{i}^{n},\mathbf{w})+\sum_{j}g_{i}^{T}(T_{i}^{n},\mathbf{w})=0,\quad\forall\mathbf{w}\in W,
$$

For the bilinear forms, we have
$$
\begin{aligned}
&b_{i}^{p}(p_{i},v_{i})=\int_{\Omega_{i}}k_{i}\nabla p_{i}\cdot\nabla v_{i}dx,\quad m_{i}^{p}(p_{i},v_{i})=\int_{\Omega_{i}}c_{i}p_{i}v_{i}dx, \\
& q_{ij}^{p}(p_{i}-p_{j},v_{i})=\int_{\Omega_{i}}\eta_{ij}(p_{i}-p_{j})v_{i}dx,\quad  d_{i}^{p}(\mathbf{u},v_{i})=\int_{\Omega_{i}}\gamma_{i}\operatorname{div}\mathbf{u}v_{i}dx,\\
&b_{i}^{T}(T_{i},z_{i})=\int_{\Omega_{i}}\chi_{i}\nabla T_{i}\cdot\nabla z_{i}dx+\int_{\Omega_{i}}C_{w}(\mathcal{V}_{i}\cdot\nabla T_{i})z_{i}dx,\quad m_{i}^{T}(T_{i},z_{i})=\int_{\Omega_{i}}s_{i}T_{i}z_{i}dx, \\
& q_{ij}^{T}(T_{i}-T_{j},z_{i})=\int_{\Omega_{i}}\Upsilon_{ij}(T_{i}-T_{j})z_{i}dx,\quad  d_{i}^{T}(\mathbf{u},z_{i})=\int_{\Omega_{i}}\delta_{i}\operatorname{div}\mathbf{u}z_{i}dx,\\
&a(\mathbf{u},\mathbf{w})=\int_{\Omega}\boldsymbol\sigma(\mathbf{u})\cdot\boldsymbol{\varepsilon}(\mathbf{w})dx,\quad 
 g_i^{p}(p_i,\mathbf{w})=\int_{\Omega_i}\gamma_i\nabla p_i\mathbf{w}dx,\quad g_i^{T}(T_i,\mathbf{w})=\int_{\Omega_i}\delta_i\nabla T_i\mathbf{w}dx,
\end{aligned}
$$
for $i,j=m,f$, with $\Omega_{m}=\Omega,\Omega_{f}=\gamma $.

\section{Coarse grid approximation using GMsFEM}
We use the Generalized Multiscale Finite Element Method (GMsFEM) to construct the coarse grid approximation of the time fractional thermoporoelasticity problem in a fractured and heterogeneous medium. In this computational algorithm, the first four steps are offline (preprocessing) steps for a given fracture geometry and heterogeneity. 


\textit{Offline stage}
\begin{itemize}
  \item Coarse grid and local domains construction.
  \item The solution of the local problems with different boundary conditions to construct a snapshot space in each local domain.
  \item An offline space construction via the solution of the local spectral problems on the snapshot space.
  \item Generation of the projection matrix.
\end{itemize}

Next, we move on to the online stage that includes the construction of the coarse grid system using the precalculated projection matrix with offline multiscale basis functions.  

\textit{Online stage}
\begin{itemize}
  \item Construction of the coarse grid system using projection matrix.
  \item Solving the problem on the coarse grid at the current time step.
  \item The multiscale space enrichment by calculation of the online basis functions, where we solve a coupled problem in each local domain using local residual information.  We enrich the offline space and update the projection matrix using the obtained online basis functions. After that, we repeatedly solve the current time step problem on the coarse grid to update the solution. 
  \item Moving to the next time step. 
\end{itemize}

Let $\mathcal{T}^{H}$ be the partition of the domain on the coarse grid and $K_j$ be the $j$th cell of $\mathcal{T}^{H}$, we have
$$\mathcal{T}^H=\bigcup_jK_j.$$
Let $N_{v}^{H}$ denote the number of nodes of $\mathcal{T}^{H}$. For each coarse grid node, we can define local domain $\omega_{l}$($l=1,\ldots,N_{v}^{H}$), which can be regarded as a combination of the several coarse grid cells that each of them contains a $l$th node. 

We start with the construction of the offline space, where the generation of the offline basis functions for pressure, temperature and displacements are given separately. For more detailed construction of multiscale bases, see \cite{24}. The offline basis construction contains two steps: (1) snapshot space construction and (2) solution of the local spectral problem on snapshot space. 

\textbf{Coarse grid system}.
With the multiscale basis functions for pressures, temperatures, and displacements constructed as described above, we define the projection matrix:

\begin{equation}
\label{f2}
R = \begin{pmatrix}
R_p & 0 & 0 \\
0 & R_T & 0 \\
0 & 0 & R_u 
\end{pmatrix}.
\end{equation}

Then we obtain the following reduced order model:

\begin{equation}
\label{f3}
\begin{split}
\xi_\tau^{(\alpha_i)}M_i^{H,p}p_i^{H,n}+\xi_\tau^{(\beta_i)}D_i^{H,p}\mathbf{u}^{H,n}+B_i^{H,p}p_i^{H,n}+\sum_{j\neq i}Q_{ij}^{H,p}(p_i^{H,n}-p_j^{H,n})\\
=\xi_\tau^{(\alpha_i)}M_i^{H,p}p_i^{H,n-1}-\xi_\tau^{(\alpha_i)}\sum_{j=2}^n\xi_{j-1}^{(\alpha_i)}M_i^{H,p}(p_i^{H,n-j+1}-p_i^{H,n-j})\\
+\xi_\tau^{(\beta_i)}D_i^{H,p}\mathbf{u}^{H,n-1}-\xi_\tau^{(\beta_i)}\sum_{j=2}^n\xi_{j-1}^{(\beta_i)}D_i^{H,p}(\mathbf{u}^{H,n-j+1}-\mathbf{u}^{H,n-j}), \\
\xi_\tau^{(\zeta_i)}M_i^{H,T}T_i^{H,n}+\xi_\tau^{(\beta_i)}D_i^{H,T}\mathbf{u}^{H,n}+B_i^{H,T}T_i^{H,n}+\sum_{j\neq i}Q_{ij}^{H,T}(T_i^{H,n}-T_j^{H,n})\\
=\xi_\tau^{(\zeta_i)}M_i^{H,T}T_i^{H,n-1}-\xi_\tau^{(\zeta_i)}\sum_{j=2}^n\xi_{j-1}^{(\zeta_i)}M_i^{H,T}(T_i^{H,n-j+1}-T_i^{H,n-j})\\
+\xi_\tau^{(\beta_i)}D_i^{H,T}\mathbf{u}^{H,n-1}-\xi_\tau^{(\beta_i)}\sum_{j=2}^n\xi_{j-1}^{(\beta_i)}D_i^{H,T}(\mathbf{u}^{H,n-j+1}-\mathbf{u}^{H,n-j}),\\
\sum_jG_j^{H,p}p_j^{H,n}+\sum_jG_j^{H,T}T_j^{H,n}+A^H_u\mathbf{u}^{H,n}=0,
\end{split}
\end{equation}
where 
\begin{equation}
\label{f4}
\begin{aligned}
M^{H,p}_i = R M^p_i R^T, \quad 
B^{H,p}_i = R B^p_i R^T, \quad 
Q^{H,p}_{ij} = R Q^p_{ij} R^T, \quad  \\
D^{H,p}_i = R D^p_i R^T, \quad
M^{H,T}_i = R M^T_i R^T, \quad 
B^{H,T}_i = R B^T_i R^T, \quad  \\
Q^{H,T}_{ij} = R Q^T_{ij} R^T, \quad 
D^{H,T}_i = R D^T_i R^T, \quad
G^{H,p}_j = R G^p_j R^T, \quad \\
G^{H,T}_j = R G^T_j R^T, \quad 
A^H_u = R A_u R^T.
\end{aligned}
\end{equation}
Following the acquisition of solutions at the coarse scale, we proceed to reconstruct the solutions at the fine scale:
\[
p^{ms,n}_i = R^T p^{H,n}_i,  \quad
T^{ms,n}_i = R^T T^{H,n}_i, \quad
\mathbf{u}^{ms,n} = R^T \mathbf{u}^{H,n}. 
\]
In the method outlined above, it's worth noting that we exclusively store and utilize information pertaining to the coarse-grid solutions from the preceding time step.

\textbf{Online enrichment of multiscale space}
To enhance the accuracy of the proposed multiscale approximation, we introduce the development of online residual-based multiscale basis functions. These coupled online basis functions are computed after solving the coarse-scale system in the offline space, utilizing residual information during the online stage. To construct the local residual-based online multiscale basis functions, we solve the following local problems in each $\omega_l$. Find $(\Upsilon_k^{l, p},  \Upsilon_k^{l, u}, \Upsilon_k^{l, T}) \in W^h(\omega_l) \times V^h(\omega_l) \times V^h(\omega_l)$ such that
$$
\begin{aligned}
\xi_{\tau}^{(\alpha_{\omega_l})}m_{\omega_l}^{p}(\Upsilon_k^{l, p, n},v) +\xi_{\tau}^{(\alpha_{\omega_l})}\sum_{j=2}^{n}\xi_{j-1}^{(\alpha_{\omega_l})}m_{\omega_l}^{p}(\Upsilon_k^{l,p,n-j+1},v) +\xi_{\tau}^{(\beta_{\omega_l})}d_{\omega_l}^{p}(\Upsilon_k^{l, u, n},v) +  \\
+\xi_{\tau}^{(\beta_{\omega_l})}\sum_{j=2}^{n}\xi_{j-1}^{(\beta_{\omega_l})}d_{\omega_l}^{p}(\Upsilon_k^{l, u, n-j+1},v) +b_{\omega_l}^{p}(\Upsilon_k^{l, p, n},v) = r^{p,k-1}_{\omega_l},\quad\forall v\in V, \\
\xi_{\tau}^{(\zeta_{\omega_l})}m_{\omega_l}^{T}(\Upsilon_k^{l, T, n},z) +\xi_{\tau}^{(\zeta_{\omega_l})}\sum_{j=2}^{n}\xi_{j-1}^{(\zeta_{\omega_l})}m_{\omega_l}^{T}(\Upsilon_k^{l, T, n-j+1},z) +\xi_{\tau}^{(\beta_{\omega_l})}d_{\omega_l}^{T}(\Upsilon_k^{l, u, n},z) +  \\
+\xi_{\tau}^{(\beta_{\omega_l})}\sum_{j=2}^{n}\xi_{j-1}^{(\beta_{\omega_l})}d_{i}^{T}(\Upsilon_k^{l, u, n-j+1},z) + b_{\omega_l}^{T}(\Upsilon_k^{l, T, n},z)=r^{T,k-1}_{\omega_l},\quad\forall z\in V, \\
a(\Upsilon_k^{l, u, n},\mathbf{w})+\sum_{j}g_{\omega_l}^{p}(\Upsilon_k^{l, p, n},\mathbf{w})+\sum_{j}g_{\omega_l}^{T}(\Upsilon_k^{l, T, n},\mathbf{w})=r^{u,k-1}_{\omega_l},\quad\forall\mathbf{w}\in W,
\end{aligned}
$$
with
\[
{V}^h(\omega_l) = \lbrace
v \in [H^1(\omega_l)]^d: v = 0 \text{ on } \partial \omega_l
\rbrace, 
\quad 
{W}^h(\omega_l) = \lbrace
w \in H^1(\omega_l): w = 0 \text{ on }\partial \omega_l
\rbrace. 
\]


Here, we have the following bilinear forms
$$
\begin{aligned}
&b_{\omega_l}^{p}(p,v)=\int_{\omega_{l}}k_m\nabla p_{m}\cdot\nabla v dx + \int_{\gamma^{\omega_{l}}} k_f\nabla p_{f}\cdot\nabla v_f dx,\quad  \\
&m_{\omega_l}^{p}(p,v)=\int_{\omega_{l}}c_{m}p_{m}v dx + \int_{\gamma^{\omega_{l}}}c_{f}p_{f}v_f dx, \\
&d_{\omega_l}^{p}(\mathbf{u},v)=\int_{\omega_{l}}\gamma_{m}\operatorname{div}\mathbf{u}vdx + \int_{\gamma^{\omega_{l}}} \gamma_{f}\operatorname{div}\mathbf{u}v_fdx, \quad a_{
\omega_l}(\mathbf{u},\mathbf{w})=\int_{\omega_l}\boldsymbol\sigma(\mathbf{u})\cdot\boldsymbol{\varepsilon}(\mathbf{w})dx,\\
&b_{\omega_l}^{T}(T,z)=\int_{\omega_{l}}\chi_{m}\nabla T_{m}\cdot\nabla zdx+\int_{\omega_{l}}C_{w}(\mathcal{V}_{m}\cdot\nabla T_{m})zdx + \\
+ &\int_{\gamma^{\omega_{l}}}\chi_{f}\nabla T_{f}\cdot\nabla z_fdx +\int_{\gamma^{\omega_{l}}}C_{w}(\mathcal{V}_{f}\cdot\nabla T_{f})z_fdx, \\
&d_{\omega_l}^{T}(\mathbf{u},z)=\int_{\omega_{l}}\delta_{m}\operatorname{div}\mathbf{u}z dx + \int_{\gamma^{\omega_{l}}}\delta_{f}\operatorname{div}\mathbf{u}z_f dx, \quad \\
&g_{\omega_l}^{p}(p,\mathbf{w})=\int_{\Omega_i}\gamma_m\nabla p_m\mathbf{w}dx + \int_{\gamma^{\omega_{l}}}\gamma_f\nabla p_f\mathbf{w_f}dx, \\
&g_{\omega_l}^{T}(T,\mathbf{w})=\int_{\omega_l}\delta_m\nabla T_m\mathbf{w}dx + \int_{\gamma^{\omega_{l}}}\delta_f\nabla T_f\mathbf{w_f}dx, \quad m_{\omega_l}^{T}(T,z)=\int_{\omega_{l}}s_{m}T_{m}zdx + \\
\end{aligned}
$$
for $i,j=m,f$, with $\Omega_{m}=\Omega,\Omega_{f}=\gamma $.
and the right-hand side is based on the local residual information
$$
\begin{aligned}
r^{p, k}_{\omega_l} =  - \xi_{\tau}^{(\alpha_{\omega_l})}m_{\omega_l}^{p}(p^{k,n-1}_{ms}-\hat{p}_{ms}^{n-1},v) - \xi_{\tau}^{(\alpha_{\omega_l})}\sum_{j=2}^{n}\xi_{j-1}^{(\alpha_{\omega_l})}m_{\omega_l}^{p}(p^{k,n-j}_{ms}-\hat{p}_{ms}^{n-j},v)\\
- \xi_{\tau}^{(\beta_{\omega_l})}d_{\omega_l}^{p}(u^{k,n-1}_{ms}-\hat{u}_{ms}^{n-1},v) - \xi_{\tau}^{(\beta_{\omega_l})}\sum_{j=2}^{n}\xi_{j-1}^{(\beta_{\omega_l})}d_{\omega_l}^{p}(u^{k,n-j}_{ms}-\hat{u}_{ms}^{n-j},v) - b_{\omega_l}^{p}(p_{ms}^{k,n},v), \\ 
\end{aligned}
$$
$$
\begin{aligned}
r^{T, k}_{\omega_l} =  - \xi_{\tau}^{(\alpha_{\omega_l})}m_{\omega_l}^{T}(T^{k,n-1}_{ms}-\hat{T}_{ms}^{n-1},z) - \xi_{\tau}^{(\alpha_{\omega_l})}\sum_{j=2}^{n}\xi_{j-1}^{(\alpha_{\omega_l})}m_{\omega_l}^{T}(T^{k,n-j}_{ms}-\hat{T}_{ms}^{n-j},z) \\
 - \xi_{\tau}^{(\beta_{\omega_l})}d_{\omega_l}^{T}(u^{k,n-1}_{ms}-\hat{u}_{ms}^{n-1},z) - \xi_{\tau}^{(\beta_{\omega_l})}\sum_{j=2}^{n}\xi_{j-1}^{(\beta_{\omega_l})}d_{\omega_l}^{T}(u^{k,n-j}_{ms}-\hat{u}_{ms}^{n-j},z) - b_{\omega_l}^{T}(T_{ms}^{k,n},z), \\
\end{aligned}
$$
$$
\begin{aligned}
r^{u,k} = - a(u_{ms}^{k, n},\mathbf{w}) - \sum_{j}g_{\omega_l}^{p}(p_{ms}^{k,n},\mathbf{w}) - \sum_{j}g_{\omega_l}^{T}(T_{ms}^{k, n},\mathbf{w}),
\end{aligned}
$$

Using the constructed coupled online basis functions, we enrich the offline spaces $V_{ms}$, $V_{ms}$ and $W_{ms}$ by adding $\Upsilon_k^{l, p}$, $\Upsilon_k^{l, T}$ and $\Upsilon_k^{l, u}$
\begin{equation}
\label{f5}
\begin{split}
V_{ms} = \text{span} \lbrace 
\phi^{l,j},  \, 
\Upsilon_k^{l, p},
\, l = 1,...,N^H_v,  \, j = 1,...,M^{l,p},  \, k = 1,2, ...,
\rbrace, \\
V_{ms} = \text{span} \lbrace 
\phi^{l,j},  \, 
\Upsilon_k^{l, T},
\, l = 1,...,N^H_v,  \, j = 1,...,M^{l,T},  \, k = 1,2, ...,
\rbrace, \\
W_{ms} = \text{span} \lbrace 
\Phi^{l,j},  \,
\Upsilon_k^{l, u},  
\,  l = 1,...,N^H_v,  \, j = 1,...,M^{l,u}, \, k = 1,2, ...
\rbrace, 
\end{split}
\end{equation}
where $k$ is the number of online iteration for the current time step. We will update online basis functions for some time steps.    


Next, we present an algorithm for the multiscale method with the online residual-based multiscale basis functions. 
Let $R^{off}$ be the projection matrix constructed using offline multiscale basis functions
 \[
R^{off} = \begin{pmatrix}
R^{off}_p & 0 & 0 \\ 
0 & R^{off}_T & 0 \\
0 & 0 & R^{off}_u
\end{pmatrix}, 
\]\[
R^{off}_p = (\phi^{1,1}, ... , \phi^{N^H_v, M^{N^H_v,p}})^T.
\]\[
R^{off}_T = (\phi^{1,1}, ... , \phi^{N^H_v, M^{N^H_v,T}})^T.
\]\[
R^{off}_u = (\Phi^{1,1}, ... ,\Phi^{N^H_v, M^{N^H_v,u}})^T,
\]

\textit{Multiscale algorithm with online enrichment} 
\begin{itemize}
\item Define projection matrix $R = \check{R}$ for the current time step $n$ with $R = R^{off}$ for $n=0$ ($\check{R}$ is the projection matrix from the previous time step). 
\item Construct and solve the coarse-scale problem at the current time step.
\begin{itemize}
\item If we want to add/update online basis functions for the current time step, we solve
$$
\begin{aligned}
&\xi_\tau^{(\alpha_i)}M_i^{H,k-1,p}p_i^{H,k,n}+\xi_\tau^{(\beta_i)}D_i^{H,k-1,p}\mathbf{u}^{H,k,n}+B_i^{H,k-1,p}p_i^{H,k,n}\\
&+\sum_{j\neq i}Q_{ij}^{H,k-1,p}(p_i^{H,k,n}-p_j^{H,k-1,n}) \\
&=\xi_\tau^{(\alpha_i)}M_i^{H,k-1,p}p_i^{H,k,n-1}-\xi_\tau^{(\alpha_i)}\sum_{j=2}^n\xi_{j-1}^{(\alpha_i)}M_i^{H,k-1,p}(p_i^{H,k,n-j+1}-p_i^{H,k-1,n-j})\\
&+\xi_\tau^{(\beta_i)}D_i^{H,k-1,p}\mathbf{u}^{H,k,n-1}-\xi_\tau^{(\beta_i)}\sum_{j=2}^n\xi_{j-1}^{(\beta_i)}D_i^{H,k-1,p}(\mathbf{u}^{H,k,n-j+1}-\mathbf{u}^{H,k-1,n-j}),\\
&\xi_\tau^{(\zeta_i)}M_i^{H,k-1,T}T_i^{H,k,n}+\xi_\tau^{(\beta_i)}D_i^{H,k-1,T}\mathbf{u}^{H,k,n}+B_i^{H,k-1,T}T_i^{H,k,n}+\\
& + \sum_{j\neq i}Q_{ij}^{H,k-1,T}(T_i^{H,k,n}-T_j^{H,k-1,n}) \\
&=\xi_\tau^{(\zeta_i)}M_i^{H,k-1,T}T_i^{H,k,n-1}-\xi_\tau^{(\zeta_i)}\sum_{j=2}^n\xi_{j-1}^{(\zeta_i)}M_i^{H,k-1,T}(T_i^{H,k,n-j+1}-T_i^{H,k-1,n-j})\\
&+\xi_\tau^{(\beta_i)}D_i^{H,k-1,T}\mathbf{u}^{H,k,n-1}-\xi_\tau^{(\beta_i)}\sum_{j=2}^n\xi_{j-1}^{(\beta_i)}D_i^{H,k-1,T}(\mathbf{u}^{H,k,n-j+1}-\mathbf{u}^{H,k-1,n-j}),\\
&\sum_jG_j^{H,k-1,p}p_j^{H,k,n}+\sum_jG_j^{H,k-1,T}T_j^{H,k,n}+A^{H,k-1}_u\mathbf{u}^{H,k,n}=0,
\end{aligned}
$$
for $k = 1,2,...$ with  
\[
M^{H,p}_i = R M^p_i R^T, \quad 
B^{H,p}_i = R B^p_i R^T, \quad 
Q^{H,p}_{ij} = R Q^p_{ij} R^T, \quad 
D^{H,p}_i = R D^p_i R^T,
\]
\[
M^{H,T}_i = R M^T_i R^T, \quad 
B^{H,T}_i = R B^T_i R^T, \quad 
Q^{H,T}_{ij} = R Q^T_{ij} R^T, \quad 
D^{H,T}_i = R D^T_i R^T, 
\]
\[
G^{H,p}_j = R G^p_j R^T, \quad 
G^{H,T}_j = R G^T_j R^T, \quad 
A^H_u = R A_u R^T,
\]
and \[
p^{ms,n}_i = R^T p^{H,n}_i,  \quad
T^{ms,n}_i = R^T T^{H,n}_i, \quad
\mathbf{u}^{ms,n} = R^T \mathbf{u}^{H,n}. 
\]

For the projection matrix, we have $R^k = R^{off}$ for $k = 0$ and 
 \[
R^{k} = \begin{pmatrix}
R^{k}_p & 0 & 0 \\
0 & R^{k}_T & 0 \\
0 & 0 & R^{k}_u
\end{pmatrix}, 
\]\[
R^{k}_p = (\phi^{1,1}, ... , \phi^{N^H_v, M^{N^H_v,p}},   
\Upsilon_1^{1,  p}, ..., \Upsilon_1^{N^H_v, p},...,
\Upsilon_k^{1,  p}, ..., \Upsilon_k^{N^H_v, p} )^T.
\]\[
R^{k}_T = (\phi^{1,1}, ... , \phi^{N^H_v, M^{N^H_v,T}},   
\Upsilon_1^{1,  T}, ..., \Upsilon_1^{N^H_v, T},...,
\Upsilon_k^{1,  T}, ..., \Upsilon_k^{N^H_v, T} )^T.
\]\[
R^{k}_u = (\Phi^{1,1}, ... ,\Phi^{N^H_v, M^{N^H_v,u}},   
\Upsilon_1^{1, u}, ..., \Upsilon_1^{N^H_v, u}, ...,
\Upsilon_k^{1, u}, ..., \Upsilon_k^{N^H_v, u} )^T, 
\]
for $k=1,2,...$.   

Here, online basis functions $(\Upsilon_k^{l, p},\Upsilon_k^{l, T}, \Upsilon_k^{l, u})$ are calculated using the solution from previous iteration $(p_{ms}^{k-1}, T_{ms}^{k-1}, u_{ms}^{k-1})$.
\item Else, we solve equation \ref{f3} with \ref{f4}.
\end{itemize}
\item Move to the next time step. % and repeating the online stage steps.
\end{itemize}
In general, online basis functions can be adaptively added only in certain local domains with a significant residual \cite{20, 21}. Next, we provide numerical results for heterogeneous and fractured poroelastic media.

\section{Numerical results}

\begin{figure}[t!]
\centering
\includegraphics[width=0.30\textwidth]{pics/meshfrac1.png}
\includegraphics[width=0.30\textwidth]{pics/meshfrac2.png}
\caption{Computation domain and grids.  Coarse grid (blue color), fine grid (green), and fractures (red).}
\label{fig1}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=0.32\linewidth]{pics/E.png}
\includegraphics[width=0.32\linewidth]{pics/kdmu.png}
\includegraphics[width=0.32\linewidth]{pics/Dif.png}
\caption{Elasticity coefficient $E$(left), heterogeneous permeability $k_m$(center) and thermal conductivity $\chi_m$(right).}
\label{fig2}
\end{figure}

In this section, we present the numerical results related to thermoporoelasticity problems in heterogeneous and fractured media that incorporate fractional derivatives. The coarse grid is characterized by uniform rectangular cells. Figure \ref{fig1} illustrates both the coarse and fine computational grids. The fine grid consists of $25,846$ cells and $12,944$ vertices, whereas the coarse grid has $121$ vertices and $100$ cells. Our analysis centers on the time-fractional diffusion equation relevant to the thermoporoelasticity issue within $\Omega = (0, 50)^2$, specifically focusing on thermoporoelasticity in fractured media. For the coefficients that represent the properties of the matrix and fractures, we set $\gamma_m = 0.1, \gamma_f = 0, \delta_m = 0.1, \delta_f = 0.1, s_m = 10^2, s_f = 10^2, k_f = 1.0, M_m = 10, M_f = 10^3, \nu = 0.3$, $\chi_m = 400$, and $\chi_f = 1.0$. The computations are performed with $T_{max} = 86,400$, using a time step of $\tau = 8,640$, $\eta_{mf} = k_f$, and $\Upsilon_{mf} = \chi_f$. The heterogeneous coefficients for Young's modulus, permeability, and thermal conductivity are depicted in Fig. \ref{fig2}. A numerical solution is provided under the boundary conditions $\mathbf{u}_x = 0, \boldsymbol\sigma_y = 0, x \in \Gamma_L \cup \Gamma_R, \mathbf{u}_y = 0, \boldsymbol\sigma_x = 0, x \in \Gamma_T \cup \Gamma_B$ for displacement, along with the initial conditions $p^0 = 1$ and $T^0 = 40$ for pressure and temperature, respectively.

\begin{table}[h!]
\centering
\begin{tabular}{ | c | c | c | c | c | c | c | c |}
\hline
$M$ & $DOF_H$ & $e_{L_2}^u$ (\%) & $e_{H_1}^u$  (\%) & $e_{L_2}^{p}$  (\%) & $e_{H_1}^{p}$  (\%) & $e_{L_2}^{T}$  (\%) & $e_{H_1}^{T}$  (\%) \\ \hline
\multicolumn{7}{|c|}{Offline Basis} \\
\hline
1 & 484 & 9.191 & 32.819 & 5.121 & 79.421 & 0.258 & 57.072\\
2 & 968 & 5.078 & 18.778 & 3.880 & 66.501 & 0.203 & 47.761\\
4 & 1936 & 2.109 & 11.579 & 2.142 & 46.316 & 0.098 & 30.354\\
8 & 3872 & 0.978 & 7.360 & 1.185 & 33.598 & 0.046 & 19.767\\
12 & 5808 & 0.716 & 5.877 & 0.824 & 27.569 & 0.032 & 15.951 \\
16 & 7744 & 0.614 & 5.221 & 0.695 & 25.133 & 0.028 & 14.716\\
\hline
\multicolumn{7}{|c|}{Online basis} \\
\hline
1+1 & 968 & 6.328 & 26.291 & 2.525 & 46.063 & 0.144 & 30.067\\
2+1 & 1452 & 3.101 & 13.439 & 2.166 & 40.867 & 0.138 & 30.273\\
4+1 & 2420 & 1.114 & 7.071 & 1.007 & 25.237 & 0.053 & 16.027\\
8+1 & 4356 & 0.343 & 3.063 & 0.451 & 15.412 & 0.017 & 7.620\\
12+1 & 6292 & 0.190 & 2.027 & 0.281 & 11.627 & 0.007 & 4.134 \\
16+1 & 8228 & 0.149 & 1.518 & 0.238 & 10.346 & 0.007 & 3.978\\
\hline
\end{tabular}
\caption{Relative errors for displacement, pressure and temperature with fractional order derivative $\alpha = 0.8$.}
\label{t1}
\end{table}

\begin{table}[h!]
\centering
\begin{tabular}{ | c | c | c | c | c | c | c | c |}
\hline
$M$ & $DOF_H$ & $e_{L_2}^u$ (\%) & $e_{H_1}^u$  (\%) & $e_{L_2}^{p}$  (\%) & $e_{H_1}^{p}$  (\%) & $e_{L_2}^{T}$  (\%) & $e_{H_1}^{T}$  (\%) \\ \hline
\multicolumn{7}{|c|}{Offline Basis} \\
\hline
1 & 484 & 8.390 & 27.900 & 7.463 & 69.941 & 0.570 & 51.236\\
2 & 968 & 4.376 & 15.820 & 3.880 & 56.429 & 0.388 & 40.967\\
4 & 1936 & 1.854 & 9.029 & 2.274 & 36.302 & 0.170 & 25.355\\
8 & 3872 & 0.813 & 5.387 & 1.054 & 24.103 & 0.070 & 15.639\\
12 & 5808 & 0.572 & 4.304 & 0.700 & 19.316 & 0.047 & 12.614 \\
16 & 7744 & 0.484 & 3.809 & 0.604 & 17.792 & 0.041 & 11.699\\
\hline
\multicolumn{7}{|c|}{Online basis} \\
\hline
1+1 & 968 & 4.858 & 21.841 & 3.762 & 44.260 & 0.290 & 29.900\\
2+1 & 1452 & 2.674 & 11.409 & 2.818 & 38.431 & 0.250 & 27.642\\
4+1 & 2420 & 0.807 & 5.493 & 1.286 & 22.869 & 0.085 & 14.166\\
8+1 & 4356 & 0.228 & 2.193 & 0.446 & 12.219 & 0.025 & 6.334\\
12+1 & 6292 & 0.117 & 1.492 & 0.245 & 8.382 & 0.009 & 3.341 \\
16+1 & 8228 & 0.084 & 1.068 & 0.221 & 7.701 & 0.009 & 3.210\\
\hline
\end{tabular}
\caption{Relative errors for displacement, pressure and temperature with fractional order derivative $\alpha = 0.9$.}
\label{t2}
\end{table}

\begin{table}[h!]
\centering
\begin{tabular}{ | c | c | c | c | c | c | c | c |}
\hline
$M$ & $DOF_H$ & $e_{L_2}^u$ (\%) & $e_{H_1}^u$  (\%) & $e_{L_2}^{p}$  (\%) & $e_{H_1}^{p}$  (\%) & $e_{L_2}^{T}$  (\%) & $e_{H_1}^{T}$  (\%) \\ \hline
\multicolumn{7}{|c|}{Offline Basis} \\
\hline
1 & 484 & 9.060 & 25.4 & 8.677 & 63.328 & 1.300 & 48.331\\
2 & 968 & 4.502 & 14.211 & 4.900 & 48.362 & 0.772 & 37.176\\
4 & 1936 & 1.616 & 7.248 & 1.931 & 28.852 & 0.284 & 21.761\\
8 & 3872 & 0.773 & 4.037 & 0.792 & 17.814 & 0.116 & 13.437\\
12 & 5808 & 0.596 & 3.255 & 0.517 & 14.176 & 0.082 & 11.060 \\
16 & 7744 & 0.514 & 2.901 & 0.457 & 13.169 & 0.073 & 10.310\\
\hline
\multicolumn{7}{|c|}{Online basis} \\
\hline
1+1 & 968 & 4.007 & 18.806 & 4.550 & 43.186 & 0.613 & 30.385\\
2+1 & 1452 & 2.730 & 10.299 & 3.001 & 34.879 & 0.512 & 27.055\\
4+1 & 2420 & 0.584 & 4.341 & 1.154 & 19.103 & 0.131 & 12.230\\
8+1 & 4356 & 0.147 & 1.614 & 0.321 & 8.991 & 0.036 & 5.465\\
12+1 & 6292 & 0.075 & 1.125 & 0.162 & 5.832 & 0.014 & 2.937 \\
16+1 & 8228 & 0.050 & 0.826 & 0.148 & 5.388 & 0.012 & 2.586\\
\hline
\end{tabular}
\caption{Relative errors for displacement, pressure and temperature with fractional order derivative $\alpha = 1.0$.}
\label{t3}
\end{table}

\begin{figure}[t!]
\centering
\includegraphics[width=0.24\linewidth]{pics/08exactp}
\includegraphics[width=0.24\linewidth]{pics/08exactt}  
\includegraphics[width=0.24\linewidth]{pics/08exactux}
\includegraphics[width=0.24\linewidth]{pics/08exactuy} \\
\includegraphics[width=0.24\linewidth]{pics/08multip}
\includegraphics[width=0.24\linewidth]{pics/08multit}
\includegraphics[width=0.24\linewidth]{pics/08multiux}
\includegraphics[width=0.24\linewidth]{pics/08multiuy} \\
\includegraphics[width=0.24\linewidth]{pics/08onlinep}
\includegraphics[width=0.24\linewidth]{pics/08onlinet}  
\includegraphics[width=0.24\linewidth]{pics/08onlineux}
\includegraphics[width=0.24\linewidth]{pics/08onlineuy} \\
\caption{Thermoporoelasticity problem in fractured and heterogeneous media: Distribution of pressure, temperature and displacement along $X$ and $Y$ at final time for fractional order derivative $\alpha = 0.8$ (from left to right).}
\label{fig3}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=0.24\linewidth]{pics/09exactp}
\includegraphics[width=0.24\linewidth]{pics/09exactt}  
\includegraphics[width=0.24\linewidth]{pics/09exactux}
\includegraphics[width=0.24\linewidth]{pics/09exactuy} \\
\includegraphics[width=0.24\linewidth]{pics/09multip}
\includegraphics[width=0.24\linewidth]{pics/09multit}
\includegraphics[width=0.24\linewidth]{pics/09multiux}
\includegraphics[width=0.24\linewidth]{pics/09multiuy} \\
\includegraphics[width=0.24\linewidth]{pics/09onlinep}
\includegraphics[width=0.24\linewidth]{pics/09onlinet}  
\includegraphics[width=0.24\linewidth]{pics/09onlineux}
\includegraphics[width=0.24\linewidth]{pics/09onlineuy} \\
\caption{Thermoporoelasticity problem in fractured and heterogeneous media: Distribution of pressure, temperature and displacement along $X$ and $Y$ at final time for fractional order derivative $\alpha = 0.9$ (from left to right).}
\label{fig4}
\end{figure}

Tables \ref{t1}-\ref{t3} enhance our analysis by displaying the relative $??_2$ and energy $H_1$ errors associated with the use of offline and online basis functions. The comparative evaluation of errors between fine-scale and multiscale solutions, considering the different types of basis functions, reveals a significant trend: a notable reduction in error as the number of basis functions increases for each fractional order derivative. This demonstrates the effectiveness of both offline and online approaches in improving solution accuracy.

Figures \ref{fig3}–\ref{fig5} provide a comprehensive view of the spatial distribution of pressure, temperature, and displacement at the final time step. The first row presents the fine-scale solution obtained through the finite element method, serving as a benchmark for accuracy. The second row displays multiscale solutions using eight basis functions within the GMsFEM framework, showcasing the method's effectiveness in capturing complex patterns along the $X$ and $Y$ directions. Additionally, the third row introduces the multiscale solution incorporating eight multiscale basis functions one online residual-based multiscale basis function. This enhancement further refines the solution, illustrating the method's adaptability and improved precision. This visual comparison not only validates the accuracy and robustness of the proposed approach but also highlights its applicability in real-world scenarios where understanding the detailed spatial distribution of pressure, temperature, and displacement is crucial.

The analysis focused on evaluating how the number of multiscale basis functions affects the accuracy of the solution. The results indicate that as the number of basis functions increases, there is a significant reduction in relative $L_2$ and energy $H_1$ errors for displacement, pressure, and temperature. For example, with the fractional order derivative set at $\alpha = 0.8$, increasing the number of offline basis functions led to a marked decrease in errors, with the relative error for displacement reducing from approximately $9.1\%$ to less than $0.6\%$, and similar trends observed for pressure and temperature.

The use of online basis functions, which were updated dynamically based on residual information, further improved accuracy. The introduction of even a single online basis function per local domain substantially reduced errors compared to the offline basis function approach alone. For instance, the relative error for displacement was reduced to approximately $0.3\%$ when four online basis functions were used in conjunction with the offline basis set. This highlights the efficiency of the online enrichment approach in capturing local solution features more effectively. Finally, the computational efficiency of the method was assessed by comparing the offline and online approaches. The results showed that while the offline method required less computational effort, the online enrichment approach provided superior accuracy with only a moderate increase in computational cost. This balance between accuracy and efficiency makes the proposed method well-suited for large-scale simulations of thermoporoelasticity in fractured media.

\begin{figure}[h!]
\centering
\includegraphics[width=0.24\linewidth]{pics/10exactp}
\includegraphics[width=0.24\linewidth]{pics/10exactt}  
\includegraphics[width=0.24\linewidth]{pics/10exactux}
\includegraphics[width=0.24\linewidth]{pics/10exactuy} \\
\includegraphics[width=0.24\linewidth]{pics/10multip}
\includegraphics[width=0.24\linewidth]{pics/10multit}
\includegraphics[width=0.24\linewidth]{pics/10multiux}
\includegraphics[width=0.24\linewidth]{pics/10multiuy} \\
\includegraphics[width=0.24\linewidth]{pics/10onlinep}
\includegraphics[width=0.24\linewidth]{pics/10onlinet}  
\includegraphics[width=0.24\linewidth]{pics/10onlineux}
\includegraphics[width=0.24\linewidth]{pics/10onlineuy} \\
\caption{Thermoporoelasticity problem in fractured and heterogeneous media: Distribution of pressure, temperature and displacement along $X$ and $Y$ at final time for fractional order derivative $\alpha = 1.0$ (from left to right).}
\label{fig5}
\end{figure}



\section{Conclusions}

In this study, we developed and applied a multiscale model reduction technique for the time fractional thermoporoelasticity problem in heterogeneous and fractured media. By integrating fractional derivatives with the Generalized Multiscale Finite Element Method (GMsFEM) and implementing online residual-based basis function enrichment, we were able to achieve a highly accurate and computationally efficient simulation approach.

Our numerical experiments demonstrated that the use of multiscale basis functions significantly improves the accuracy of the solution, particularly when enhanced by online updates that adaptively target areas with high residuals. The results showed that increasing the number of basis functions effectively reduces errors, and that the incorporation of online basis functions yields even greater improvements in accuracy, especially in capturing the complex interactions within fractured porous media.

Overall, the proposed method offers a robust and flexible framework for modeling complex thermoporoelastic processes in fractured and heterogeneous media. Its ability to balance computational efficiency with high accuracy makes it an attractive tool for large-scale simulations in geomechanics, geothermal energy, and related fields. 

\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).

\bibitem{1} J. Bear, M. Corapcioglu, 
\href{https://doi.org/10.1029/WR017i003p00723}{\it A mathematical model for consolidation in a thermoelastic aquifer due to hot water injection or pumping}, Water Resources Research, \textbf{17}:3 (1981), 723--736.

\bibitem{2} D. A. Ammosov, M. V. Vasilyeva, E.T. Chung, 
\href{https://doi.org/10.1016/j.cam.2021.113995}{\it Generalized multiscale finite element method for thermoporoelasticity problems in heterogeneous and fractured media}, Water Resources Research, {\bf 1} (2022), 113995.

\bibitem{3} D. W. Smith, J. R. Booker,  
\href{https://10.1002/nag.1610170302}{\it Green’s functions for a fully coupled thermoporoelastic material}, International Journal for Numerical and Analytical Methods in
Geomechanics, \textbf{17}:3 (1993), 139--163.

\bibitem{4} M. Kurashige,  
\href{https://10.1016/0020-7683(89)90020-6}{\it A thermoelastic theory of fluid-filled porous materials}, International Journal of Solids and Structures, \textbf{25}:9 (1989), 1039--1052.

\bibitem{5} J. Hu, A. A. Alikhanov, Y. Efendiev, W.T. Leung,  
\href{https://10.1007/s13540-022-00092-6}{\it Partially explicit time discretization for time fractional diffusion equation}, Fractional Calculus and Applied Analysis, \textbf{25}:5 (2022), 1908--1924.

\bibitem{6} G. Alaimo, V. Piccolo, A. Cutolo, L. Deseri, M. Fraldi, M. Zingales,  
\href{https://10.1016/j.mechrescom.2019.103395}{\it A fractional order theory of poroelasticity}, Mechanics Research Communications,  {\bf 100} (2019), 103395.

\bibitem{7} A. A. Tyrylgin, M. V. Vasilyeva, A. A. Alikhanov, D. Sheen, 
\href{https://10.1016/j.cam.2022.114670}{\it A computational
macroscale model for the time fractional poroelasticity problem in fractured and heterogeneous media}, Journal of Computational and Applied Mathematics,  {\bf 418} (2023), 114670.

\bibitem{8} M. Caputo, F. Mainardi, 
\href{https://10.1007/BF02820620}{\it Linear models of dissipation in anelastic solids}, La Rivista del Nuovo Cimento, \textbf{1}:2 (1971), 161--198.

\bibitem{9} M. Caputo,
\href{https://10.1121/1.1903344}{\it Vibrations of an infinite viscoelastic layer with a dissipative memory}, The Journal of the Acoustical Society of America, \textbf{56}:3 (1974), 897--904.

\bibitem{10} M. Caputo,
\href{https://10.1029/1999WR900299}{\it Models of flux in porous media with memory}, Water Resources Research, \textbf{36}:3 (2000), 693--705.

\bibitem{11} Y. Efendiev, J. Galvis, T.Y. Hou,
\href{https://10.1016/j.jcp.2013.04.045}{\it Generalized multiscale finite element methods}, Journal of computational physics, {\bf 251} (2013), 116--135.

\bibitem{12} D. L. Brown, M. V. Vasilyeva,
\href{https://10.1016/j.cam.2015.08.007}{\it A generalized multiscale finite element method for poroelasticity problems i: linear problems}, Journal of Computational and Applied
Mathematics, {\bf 294} (2016), 372--388.

\bibitem{13} E. T. Chung, Y. Efendiev, W.T. Leung,
\href{https://10.1016/j.cma.2018.04.010}{\it Constraint energy minimizing generalized multiscale finite element method}, Computer Methods in Applied Mechanics
and Engineering, {\bf 339} (2018), 298--319.

\bibitem{14} Q. Li, Y. Wang, M.V. Vasilyeva,
\href{https://10.1016/j.cam.2017.12.040}{\it Multiscale model reduction for fluid infiltration simulation through dual-continuum porous media with localized uncertainties}, Journal of Computational and Applied Mathematics, {\bf 336} (2018), 127--146.

\bibitem{15} V.N. Alekseev, U.S. Kalachikova, Y. Yang,
\href{https://10.1134/S1995080223100037}{\it Partial Learning Using Partially Explicit Discretization for Heterogeneous Transport Problem Simulation}, Lobachevskii Journal of Mathematics, \textbf{44}:10 (2023), 4103--4115.

\bibitem{16} E.T. Chung, U.S. Kalachikova, M. V. Vasilyeva, V.N. Alekseev,
\href{https://10.1016/j.matcom.2021.11.001}{\it Generalized multiscale discontinuous Galerkin method for convection–diffusion equation in perforated media}, Mathematics and Computers in Simulation, {\bf 193} (2023), 666--688.

\bibitem{17} Y. Efendiev, S. Li, G. Yao, N. Zhang,
\href{https://10.1007/s13137-015-0075-7}{\it Hierarchical multiscale modeling for flows in fractured media using generalized multiscale finite element method}, GEM-International Journal on Geomathematics, {\bf 6} (2015), 141--162.

\bibitem{18} I. Y. Akkutlu, Y. Efendiev, M. V. Vasilyeva, Y. Wang,
\href{https://10.1016/j.jngse.2017.02.040}{\it Multiscale model reduction for shale gas transport in a coupled discrete fracture and dual-continuum porous media}, Journal of Natural Gas Science and Engineering, {\bf 48} (2017), 65--76.

\bibitem{19} I. Y. Akkutlu, Y. Efendiev, M. V. Vasilyeva,
\href{https://10.1007/s10596-016-9571-6}{\it Multiscale model reduction for shale gas transport in fractured media}, Computational Geosciences, {\bf 20} (2016), 953--973.

\bibitem{20} E. T. Chung, Y. Efendiev, W.T. Leung, M.V. Vasilyeva, Y. Wang,
\href{https://10.1080/00036811.2016.1199799}{\it Online adaptive local multiscale model reduction for heterogeneous problems in perforated domains}, Applicable Analysis, \textbf{96}:12 (2017), 2002--2031.

\bibitem{21} E. T. Chung, Y. Efendiev, W.T. Leung,
\href{https://10.1016/j.jcp.2015.07.068}{\it Residual-driven online generalized multiscale finite element methods}, Journal of Computational Physics, {\bf 302} (2015), 176--190.

\bibitem{22} D. A. Spiridonov, M. V. Vasilyeva, A.A. Tyrylgin, E.T. Chung,
\href{https://10.3390/math9121382}{\it An online generalized multiscale finite element method for unsaturated filtration problem in fractured media}, Mathematics, \textbf{9}:12 (2021), 1382.

\bibitem{23} A. A. Alikhanov, C. Huang,
\href{https://10.1016/j.amc.2021.126545}{\it A high-order l2 type difference scheme for the timefractional diffusion equation}, Applied Mathematics and Computation, {\bf 411} (2021), 126545.

\bibitem{24} A. A. Alikhanov, H. Bai, J. Huang, A. A. Tyrylgin, Y. Yang,
\href{https://10.1016/j.cam.2024.116157}{\it Multiscale model reduction for the time fractional thermoporoelasticity problem in fractured and heterogeneous media}, Journal of Computational and Applied Mathematics, {\bf 455} (2025), 116157.


\end{thebibliography}
\end{document}
