\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{\ }]


\newcommand{\grad}{\mathop{\rm grad}\nolimits}
\renewcommand{\div}{\mathop{\rm div}\nolimits}
\newcommand{\const}{\mathop{\rm const}\nolimits}

\DeclareMathOperator*{\argmin}{arg\,min}


%%%%%---------- 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,35B27,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{aaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}



\begin{document}
\renewcommand{\refname}{References}
\renewcommand{\proofname}{Proof.}
\renewcommand{\figurename}{Fig.}
\renewcommand{\tablename}{Table}
\newcommand{\old}{\textcolor{red}}


\thispagestyle{empty}

\title[\sc Multicontinuum homogenization]{\large Multicontinuum homogenization for time-fractional diffusion equation}
\author[U.S. KALACHIKOVA, D.A. AMMOSOV et al.]
{{\bf \large \textsc{U.S. Kalachikova}}\protect\orcidicon{0000-0002-1561-449X} 
{\bf \large \textsc{, D.A. Ammosov}}\protect\orcidicon{0000-0002-8616-9397}\\ 
{\bf \large \textsc{A.A. Tyrylgin}}\protect\orcidicon{0000-0002-2157-0921} 
{\bf \large \textsc{, H. Bai}}\protect\orcidicon{0009-0001-5048-5339} 
{\bf \large \textsc{, A.A. Alikhanov}}\protect\orcidicon{0000-0003-0684-6667}\\ \\ 
}% 


\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{Ammosov Dmitry Andreevich 
\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}{dmitryammosov@gmail.com}}%


\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{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 Kalachikova U.S., Ammosov D.A., Tyrylgin A.A., Bai H., Alikhanov A.A.,
Multicontinuum homogenization for time-fractional diffusion equation}
\thanks{\copyright \ 2024 Kalachikova U.S., Ammosov D.A., Tyrylgin A.A., Bai H., Alikhanov A.A.}
\thanks{\rm The work is supported by the Russian government project Science and Universities (project FSRG-2024-0003) aimed at supporting junior laboratories and the Russian Science Foundation grant 23-71-30013 (https://rscf.ru/project/23-71-30013/)}
\thanks{\it Received April, 1, 2024, Published August, 23, 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:} 
In this paper, we derive a multicontinuum time-fractional diffusion equations based on Caputo fractional derivative using a multicontinuum homogenization approach. For this purpose, we formulate cell problems with constraints considering various effects. As a result, we obtain a decomposition of the solution into macroscopic variables (continua). Assuming the smoothness of these macroscopic variables, we derive multicontinuum equations for the general case. Then, we consider a particular case of a dual-continuum model in an isotropic medium. We present numerical experiments for two-dimensional model problems with different fractional derivative orders, demonstrating the high efficiency of the proposed approach.
\medskip

\noindent{\bf Keywords:} multicontinuum homogenization, time-fractional diffusion equation, heterogeneous media, fractional derivatives, multiscale modeling.
 \end{quote}

 \begin{quote}
  \noindent{\bf Аннотация:} 
  В данной статье выводятся многоконтинуальные уравнения диффузии с дробной производной по времени на основе дробной производной Капуто с использованием метода многоконтинуального усреднения. Для этого формулируются задачи на ячейках с ограничениями, учитывающие различные эффекты. В результате строится разложение решения на макроскопические переменные (континуумы). Предполагая гладкость данных макроскопических переменных, выводятся многоконтинуальные уравнения для общего случая. Затем рассматривается частный случай модели с двойным континуумом в изотропной среде. Приводятся численные эксперименты для двумерных модельных задач с различными порядками дробной производной, которые демонстрируют высокую эффективность предлагаемого подхода.
  \medskip
  
  \noindent{\bf Ключевые слова:} многоконтинуальное усреднение, уравнение диффузии с дробной производной по времени, неоднородные среды, дробные производные, многомасштабное моделирование.
   \end{quote}
}

\bigskip

\section{Introduction}

% about fractional (U)
Time-fractional diffusion equations have developed as an effective mathematical tool for describing anomalous diffusion processes in which the rate of diffusion differs from the traditional Fickian behavior. In contrast to conventional diffusion equations that utilize integer-order derivatives, time-fractional diffusion equations integrate derivatives of non-integer order, thus encapsulating the memory and heredity characteristics of the system. In this context, we employ the Caputo derivative, which is especially appropriate for starting value issues and facilitates the integration of initial conditions in a physically significant manner \cite{caputo1967linear}. These equations are widely used to describe transport phenomena in various complex media, including porous materials and fractured reservoirs, where diffusive behavior exhibits long-term dependence and non-local effects \cite{metzler2000random, podlubny1998fractional}. The fractional derivative, generally of order $0<\alpha<1$, incorporates a time-dependent kernel that considers the history of the diffusive process, resulting in enhanced explanations of anomalous transport \cite{mainardi2022fractional}.

% about homogenization
It should be noted that many real applications are complicated by heterogeneous multiscale media. For accurate modeling in such cases, it is necessary to use detailed computational grids, which significantly increase the computational cost of the problem. Moreover, the resulting huge matrices may be ill-conditioned and require special treatment. The presence of fractional derivatives in time additionally complicates the calculations. An alternative approach is to apply various homogenization methods that allow computations on coarse computational grids, reducing the computational costs \cite{hu2020homogenization, kawamoto2023homogenization, efendiev2004numerical}. The main idea of the homogenization methods is to compute effective properties at macroscopic points by solving cell problems.

However, standard homogenization approaches can lead to inaccurate solutions in cases with high contrast. This is because having only one homogenized coefficient per macroscopic point may be not sufficient. Multicontinuum modeling approaches identify several continua in the medium and introduce separate effective properties \cite{rubinvstein1948question, barenblatt1960basic, arbogast1990derivation, panasenko2018multicontinuum, alekseev2023partial}. Thus, both low-conductivity and high-conductivity regions are taken into account. One of the new methods in this approach is the multicontinuum homogenization method \cite{efendiev2023multicontinuum, chung2024multicontinuum, leung2024some}. This method provides a rigorous and, at the same time, flexible methodology for deriving multicontinuum models. The main idea of this method is to construct special cell problems to account for different solution effects. As a result, we obtain a function decomposition over continua. Using this decomposition, we can rigorously derive the multicontinuum equations. This method has already been applied to the derivation of various multicontinuum models \cite{xie2024multicontinuum, ammosov2024multicontinuum, ammosov2023multicontinuum, kalachikova2024advancing}.

Note that one can also use multiscale finite element approaches to solve the problems with high-contrast coefficients and fractional derivatives in time on a coarse grid \cite{tyrylgin2023computational, alikhanov2025multiscale, li2022partially, hu2022partially}. These methods provide accurate solutions and can handle complex high-contrast heterogeneous media. However, they result in discrete coarse-grid models and not in the form of macroscopic laws.

% about this paper
In this paper, we apply the multicontinuum homogenization method to derive a multicontinuum time-fractional diffusion model based on Caputo fractional derivative. We construct cell problems to account for gradient effects and averages. After solving these problems, we obtain the expansion of the solution over continua. We then derive a general multicontinuum model for an arbitrary number of continua. As a particular case, we present a dual-continuum model with isotropic effective properties. To verify the proposed approach, we consider two-dimensional model problems with a high-contrast diffusion coefficient. We consider different orders of fractional derivatives. We compute the relative $L_2$ errors to evaluate the accuracy of the multicontinuum model. Numerical results show that the proposed multicontinual model provides high accuracy for different orders of the fractional derivative, significantly reducing the computational cost.

% paper structure
The paper has the following structure. In Section \ref{sec:derivation}, we present the derivation of a multicontinuum time-fractional diffusion model in general form using a multicontinuum homogenization approach.
Section \ref{sec:dual_continuum} presents a dual-continuum model with isotropic coefficients. In Section \ref{sec:numerical_results}, we conduct numerical experiments to check the proposed multicontinuum model.
Finally, we present conclusions in Section \ref{sec:conclusion}.

\section{Multicontinuum time-fractional diffusion model derivation}\label{sec:derivation}

This section presents the derivation of a multicontinuum time-fractional diffusion model in general form using a multicontinuum homogenization approach. First, let us consider the following time-fractional diffusion equation in the domain $\Omega \subset \mathbb{R}^d$ ($d = 2, 3$)
\begin{equation}\label{eq:original0}
\frac{\partial^\alpha u}{\partial t^\alpha} - \nabla \cdot \left( \kappa \nabla u  \right) = f, \quad x \in \Omega, \quad t \in (0, t_{max}].
\end{equation}
Here, $u := u(x, t)$ is a sought field, $\kappa := \kappa(x)$ is a heterogeneous coefficient, and $f := f(x)$ is a source/sink term. Note that $\kappa$ possesses high contrast such that $\text{max}_{x \in \Omega} \kappa(x) / \text{min}_{x \in \Omega} \kappa(x) \gg 1$. In fluid filtration applications, $u$ can be a pressure field, and $\kappa$ can be a permeability coefficient of the porous medium with highly conductive channels.

For the time derivative term, ${\partial^\alpha \over \partial t^\alpha}$ denotes the Caputo derivative that can be expressed as
\begin{equation}
  {\partial^\alpha u(t) \over \partial t^\alpha} = \frac{1}{\Gamma (1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha} \frac{\partial u}{\partial s}(s)ds, \quad 0<\alpha\leq 1,
\end{equation}
where $\Gamma$ is the gamma function.

The multicontinuum homogenization procedure starts with the variational formulation of the original fine-scale model. Let us complement \eqref{eq:original0} with zero Dirichlet boundary conditions and set some appropriate initial conditions. Then, we can obtain the following variational formulation
\begin{equation}\label{eq:original1}
\int_\Omega \frac{\partial^\alpha u}{\partial t^\alpha} v + \int_\Omega \kappa \nabla u \cdot \nabla v = \int_\Omega f v, \quad \forall v \in H_0^1(\Omega).
\end{equation}

Next, let us make the following assumptions. First, let us assume that the computational domain $\Omega$ is partitioned into coarse blocks $\omega$. We suppose that we can define Representative Volume Elements $R_{\omega}^p$ (RVEs) inside each coarse block $\omega$. Each RVE can approximate the whole coarse block in terms of heterogeneity. Moreover, we can define an oversampled RVE $R_\omega^+$ by extending a target (central) RVE $R_\omega = R_\omega^{p_0}$ with other RVEs $R_\omega^{p}$ around it (see Figure \ref{fig:ill}). We assume that, within each RVE, there are $N$ distinct average states, which we call continua. For each continuum $j$, we define a characteristic function $\psi_j$ such that
\[
\psi_j = \begin{cases}
  1, & \text{in continuum} ~j,\\
  0, & \text{otherwise}.
\end{cases}
\]

\begin{figure}
\centering
\includegraphics[width=0.6\textwidth]{figs/scheme_homogenization}
\caption{Illustration of $\Omega$, $\omega$, $R_\omega^+$, and $R_\omega$}
\label{fig:ill}
\end{figure}

\subsection*{Multicontinuum expansion}

Inside each RVE, we define $N$ macrosopic variables, representing the averages of the solution in the corresponding subregions (continua) as follows
\begin{equation*}
U_i (x_\omega, t) = \frac{\int_{R_\omega} u(x, t) \psi_i}{\int_{R_\omega} \psi_i},
\end{equation*}
where $x_\omega$ is a point in $R_\omega$. We assume that $U_i$ are smooth functions, if we take them over all RVEs.

Then, we can define the following expansion of the solution $u$ into these macroscopic variables in each $R_\omega$ \cite{chung2024multicontinuum}
\begin{equation*}
u = \phi_i U_i + \phi_i^m \nabla_m U_i + \phi_i^{mn} \nabla_{mn}^2 U_i + ...,
\end{equation*}
where $\phi$ function are solutions of cell problems. Note that we suppose Einstein summation convention over repeated indices.

However, in our work, we consider the following expansion
\begin{equation}\label{eq:approx1}
u \approx \phi_i U_i + \phi_i^m \nabla_m U_i.
\end{equation}
Therefore, we consider only the averages and gradients.

\subsection*{Cell problems}

We formulate the cell problems in the oversampled RVEs $R_\omega^+$ with contstraints to consider various homogenized effects of the solution. The first cell problem considers the averages of the solution
\begin{equation}
\label{eq:cell_aver_aver}
\begin{split}
&\int_{R_\omega^+}\kappa \nabla \phi_{i}\cdot \nabla v -
\sum_{j,p} {B_{ij}^p \over \int_{R_\omega^p} \psi_j^p} \int_{R_\omega^p} \psi_j^p v = 0,\\
&\int_{R_\omega^p}  \phi_{i} \psi_j^p = \delta_{ij} \int_{R_\omega^p} \psi_j^p.
\end{split}
\end{equation}

The second cell problem imposes constraints to consider gradient effects
\begin{equation}
\label{eq:cell_grad_aver}
\begin{split}
&\int_{R_\omega^+} \kappa \nabla \phi^m_{i} \cdot \nabla v -  \sum_{j,p} {B_{ij}^{mp} \over \int_{R_\omega^p}\psi_j^p} \int_{R_\omega^p}\psi_j^p v =0,\\
&\int_{R_\omega^p}  \phi^m_{i} \psi_j^p = \delta_{ij} \int_{R_\omega^p} (x_m-c_{mj})\psi_j^p ,\\
&\int_{R_\omega^{p_0}} (x_m-c_{mj})\psi_j^{p_0} =0\ \text{condition for $c$}. \\
\end{split}
\end{equation}

Note that we can obtain the following estimates of the cell problems' solutions \cite{efendiev2023multicontinuum}
\begin{equation}
\label{eq:phi_estimates}
\begin{split}
&\|\phi_i\|=O(1),\ \|\nabla \phi_i\|=O({1\over\epsilon})\\
&\|\phi_i^m\|=O(\epsilon),\ \|\nabla \phi_i^m\|=O(1),\\
\end{split}
\end{equation}
where $\epsilon$ is a diameter of RVE.

\subsection*{Derivation}

By definition, $R_\omega$ can represent the whole coarse block $\omega$ in terms of heterogeneity. In this way, one can obtain the following approximation of the variational formulation \eqref{eq:original1}
\begin{equation}
\label{eq:assump1}
\begin{split}
&\int_\Omega f v = \int_\Omega {\partial^\alpha u \over \partial t^\alpha} v + \int_\Omega \kappa \nabla u \cdot \nabla v = \sum_\omega \int_\omega {\partial^\alpha u \over \partial t^\alpha} v + \sum_\omega\int_\omega \kappa \nabla u \cdot \nabla v \approx\\
&\sum_\omega {|\omega|\over |R_\omega|} \int_{R_\omega}
 {\partial^\alpha u \over \partial t^\alpha} v +
\sum_\omega {|\omega|\over |R_\omega|} \int_{R_\omega} \kappa \nabla u \cdot \nabla v.
\end{split}
\end{equation}

Let us substitute the multicontinuum expansion \eqref{eq:approx1} into the diffusion term
\begin{equation}\label{eq:approx_mc_elliptic}
\begin{split}
\int_{R_\omega} \kappa \nabla u \cdot \nabla v \approx
\int_{R_\omega} \kappa \nabla (\phi_i U_i) \cdot\nabla v  +
\int_{R_\omega} \kappa \nabla (\phi_i^m \nabla_m U_i) \cdot\nabla v. 
\end{split}
\end{equation}

According to our assumptions, macroscopic variables $U_i$ are smooth functions. Hence, the variations of $U_i$ and $\nabla_m U_i$ are minor compared to the variations of $\phi_i$ and $\phi_i^m$. That is why we assume that the following approximations hold $\int_{R_\omega} \kappa \nabla (\phi_i U_i) \cdot \nabla v \approx \int_{R_\omega} \kappa U_i \nabla \phi_i \cdot\nabla v$ and $\int_{R_\omega} \kappa \nabla (\phi_i^m \nabla_m U_i)  \cdot\nabla v \approx \int_{R_\omega} \kappa \nabla_m U_i \nabla \phi_i^m \cdot \nabla v$.

Moreover, let us introduce the expansion of $v$
\begin{equation}\label{eq:approx_v}
v \approx \phi_s V_s + \phi_s^k \nabla_k V_s.
\end{equation}

Next, let us consider each term of the right-hand side of \eqref{eq:approx_mc_elliptic} separately. If we consider the first term and apply our assumptions above and the expansion \eqref{eq:approx_v}, we obtain the following approximation
\begin{equation}\label{eq:mc_elliptic_appr_1}
\begin{split}
&\int_{R_\omega} \kappa \nabla (\phi_i U_i) \cdot \nabla v \approx
U_i(x_{\omega}) \int_{R_\omega} \kappa \nabla \phi_i \cdot \nabla v \approx\\
&U_i (x_{\omega}) V_s(x_\omega) \int_{R_\omega} \kappa \nabla \phi_i \cdot \nabla \phi_s +
U_i(x_{\omega}) \nabla_m V_s(x_\omega) \int_{R_\omega} \kappa \nabla \phi_i \cdot \nabla \phi_s^m =\\
&U_i(x_{\omega})  B_{is} V_s(x_\omega) + B_{is}^{m} U_i(x_{\omega}) \nabla  V_s(x_\omega),
\end{split}
\end{equation}
where
\begin{equation}\label{eq:B_defintions}
B_{ik}^{m}=\int_{R_\omega} \kappa \nabla \phi_i^m \cdot \nabla \phi_k, \ \ 
B_{ik}=\int_{R_\omega} \kappa \nabla \phi_i \cdot \nabla \phi_k.
\end{equation}
Note that we use midpoint approximations of macroscopic variables in \eqref{eq:mc_elliptic_appr_1}, according to our assumptions about their smoothness.

We can appropriate the second term of \eqref{eq:approx_mc_elliptic} in a similar way.
\begin{equation}
\begin{split}
&\int_{R_\omega} \kappa \nabla (\phi_i^m \nabla_m U_i)  \cdot \nabla v \approx
\nabla_m U_i (x_\omega) \int_{R_\omega^+} \kappa \nabla \phi_i^m\cdot \nabla v \approx\\
&\nabla_m U_i (x_\omega)  \nabla_k V_s (x_\omega) \int_{R_\omega} \kappa \nabla \phi_i^m \cdot \nabla \phi_s^k +
\nabla_m U_i (x_\omega)  V_s (x_\omega) \int_{R_\omega} \kappa \nabla \phi_i^m \cdot \nabla \phi_s =\\
&\nabla_m U_i (x_\omega)  \nabla_k V_s (x_\omega) A_{is}^{km} + \nabla_m U_i (x_\omega)  V_s (x_\omega) B_{is}^{m},
\end{split}
\end{equation}
where
\begin{equation}\label{eq:A_definition}
A_{is}^{km}= \int_{R_\omega} \kappa \nabla \phi_i^m \cdot \nabla \phi_s^k.
\end{equation}

Then, we apply continuous approximations for $U_i$ and $V_j$ and obtain the following approximation of the diffusion term
\begin{equation}\label{eq:mc_elliptic_appr_final}
\begin{split}
\int_{R_\omega} \kappa \nabla u  \cdot \nabla v \approx
U_i  B_{ij}^{n} \nabla_n V_j + U_i  B_{ij} V_j +
\nabla_m U_i A_{ij}^{mn} \nabla_n V_j + \nabla_m U_i B_{ij}^{m} V_j.
\end{split}
\end{equation}

Before we proceed to the time-fractional derivative term, let us make some remarks. First, as one can see from the cell problems \eqref{eq:cell_aver_aver} and \eqref{eq:cell_grad_aver}, $\phi_j$ and $\phi_j^m$ do not depend on time, since we suppose that $\psi_j$ does not change in time. Next, the following approximations hold $u \approx \phi_i U_i + \phi_i^m \nabla_m U_i \approx \phi_i U_i$ and $v \approx \phi_j V_j + \phi_j^n \nabla_n V_j \approx \phi_j V_j$ due to the estimates \eqref{eq:phi_estimates}. Then, we can approximate the time-fractional derivative term as follows
\begin{equation}
\begin{split}
&\int_{R_\omega} {\partial^\alpha u \over \partial t^\alpha} v \approx
\int_{R_\omega} 
{\partial^\alpha ( \phi_i U_i  + \phi_i^m \nabla_m U_i) \over \partial t^\alpha}
( \phi_j V_j  + \phi_j^n \nabla_n V_j) \approx\\
&\int_{R_\omega} {\partial^\alpha ( \phi_i U_i)   \over \partial t^\alpha} \phi_j V_j \approx
{\partial^\alpha U_i(x_\omega) \over \partial t^\alpha} V_j(x_\omega) \int_{R_\omega} \phi_i \phi_j = 
{\partial^\alpha U_i(x_\omega) \over \partial t^\alpha} V_j(x_\omega) C_{ij},
\end{split}
\end{equation}
where
\begin{equation}\label{eq:C_definition}
C_{ij} = \int_{R_\omega} \phi_i \phi_j.
\end{equation}

Again, we apply continuous approximations for $U_i$ and $V_j$ and obtain the following approximation of the time-fractional derivative term
\begin{equation}\label{eq:mc_time_appr_final}
\int_{R_\omega} {\partial^\alpha u \over \partial t^\alpha} v \approx
{\partial^\alpha U_i \over \partial t^\alpha} C_{ij} V_j
\end{equation}

Next, let us estimate $B_{ij}$, $B_{ij}^{m}$, $A_{ij}^{mn}$, and $C_{ij}$ using their definitions \eqref{eq:B_defintions}, \eqref{eq:A_definition}, and \eqref{eq:C_definition} and the estimates of $\phi_i$ and $\phi_i^m$ \eqref{eq:phi_estimates}
\begin{equation*}
B_{ij} =O(\cfrac{|R_\omega|}{\epsilon^{2}}),\
B_{ij}^{m}=O(\cfrac{|R_\omega|}{\epsilon}),\ 
A_{ij}^{mn}=O(|R_\omega|),\
C_{ij}=O(|R_\omega|).
\end{equation*}

Using these estimates, we introduce the scaled effective properties $\widehat{B_{ij}}$, $\widehat{B_{ij}^{m}}$, $\widehat{A_{ij}^{mn}}$, and $\widehat{C_{ij}}$ in the following way
\begin{equation}\label{eq:rescaled_coefficients}
\begin{split}
\widehat{B_{ij}}  =\cfrac{\epsilon^{2}}{|R_\omega|}B_{ij},\ 
\widehat{B_{ij}^{m}}  =\cfrac{\epsilon}{|R_\omega|}B_{ij}^{m},\ 
\widehat{A_{ij}^{mn}} = \frac{1}{|R_\omega|}A_{ij}^{mn},\
\widehat{C_{ij}} = \frac{1}{|R_\omega|}C_{ij}.
\end{split}
\end{equation}

Then, with these effective properties, we obtain the following approximation of the variational formulation \eqref{eq:original1}
\begin{equation}
\label{eq:mc1}
\begin{split}
&\int_\Omega {\partial^\alpha u \over \partial t^\alpha} v +
\int_\Omega \kappa \nabla u \cdot \nabla v \approx
\int_\Omega \widehat{C_{ij}} {\partial^\alpha U_i \over \partial t^\alpha} V_j+
 \int_\Omega \widehat{A_{ij}^{mn}} \nabla_m U_i \nabla_n V_j +\\
&{1\over\epsilon}\int_\Omega \widehat{B_{ij}^m} \nabla_m U_i V_j +
{1\over\epsilon}\int_\Omega  \widehat{B_{ij}^m} U_i \nabla_m V_j +
{1\over\epsilon^2}\int_\Omega \widehat{B_{ij}} U_i V_j.
\end{split}
\end{equation}

Using integration by parts, one can show that the sum of the third and fourth terms is negligible \cite{efendiev2023multicontinuum}. Thus, we obtain the following multicontinuum time-fractional model in a strong form
\begin{equation}\label{eq:mc_index}
\begin{split}
\widehat{C_{ij}} {\partial^\alpha U_j \over \partial t^\alpha} -\nabla_n (\widehat{A_{ij}^{mn}} \nabla_m U_j) +
{1\over \epsilon^2}\widehat{B_{ij}} U_j = f_i.
\end{split}
\end{equation}

As one can see, the reaction terms will dominant unless we have high effective diffusivity.


\section{Dual-continuum time-fractional diffusion model}\label{sec:dual_continuum}

Suppose that we have two continua in our medium. Let the first continuum be low conductive, and the second continuum be highly conductive. Next, we expand the general model \eqref{eq:mc_index}
\begin{equation*}
\begin{split}
&\widehat{C_{11}} {\partial^\alpha U_1 \over \partial t^\alpha} +
\widehat{C_{12}} {\partial^\alpha U_2 \over \partial t^\alpha} -
\nabla_n (\widehat{A_{11}^{mn}} \nabla_m U_1) - 
\nabla_n (\widehat{A_{12}^{mn}} \nabla_m U_2) +\\
&{1\over \epsilon^2} \widehat{B_{11}} U_1 +
{1\over \epsilon^2} \widehat{B_{12}} U_2  =
f_1, \\
&\widehat{C_{21}} {\partial^\alpha U_1 \over \partial t^\alpha} +
\widehat{C_{22}} {\partial^\alpha U_2 \over \partial t^\alpha} -
\nabla_n (\widehat{A_{21}^{mn}} \nabla_m U_1) - 
\nabla_n (\widehat{A_{22}^{mn}} \nabla_m U_2) +\\
&{1\over \epsilon^2} \widehat{B_{21}} U_1 +
{1\over \epsilon^2} \widehat{B_{22}} U_2  =
f_2.
\end{split}
\end{equation*}

Let us assume that the microstructure of the medium results in isotropic effective properties. Therefore, we can approximate the diffusion terms as $\widehat{A_{ij}^{mn}} \approx D_{ij} \delta_{mn}$. Next, note that we have the following properties of $\widehat{B_{ij}}$ \cite{ammosov2023multicontinuum}
\begin{equation}\label{eq:B_properties_2}
\widehat{B_{12}} \approx -\widehat{B_{11}}, \quad 
\widehat{B_{21}} \approx -\widehat{B_{22}}, \quad 
\widehat{B_{12}} = \widehat{B_{21}}, \quad
\widehat{B_{11}} \approx \widehat{B_{22}} \approx \widehat{B} > 0.
\end{equation}

Therefore, we obtain the following dual-continuum time-fractional diffusion model in the isotropic case
\begin{equation}\label{eq:dual_continuum_model}
\begin{split}
&\widehat{C_{11}} {\partial^\alpha U_1 \over \partial t^\alpha} +
\widehat{C_{12}} {\partial^\alpha U_2 \over \partial t^\alpha} -
\nabla \cdot (\widehat{D_{11}} \nabla U_1) - 
\nabla \cdot (\widehat{D_{12}} \nabla U_2) +
{1\over \epsilon^2} \widehat{B} (U_1 - U_2)  =
f_1, \\
&\widehat{C_{21}} {\partial^\alpha U_1 \over \partial t^\alpha} +
\widehat{C_{22}} {\partial^\alpha U_2 \over \partial t^\alpha} -
\nabla \cdot (\widehat{D_{21}} \nabla U_1) - 
\nabla \cdot (\widehat{D_{22}} \nabla U_2) +
{1\over \epsilon^2} \widehat{B} (U_2 - U_1)  =
f_2.
\end{split}
\end{equation}

\subsection*{Approximation of the time-fractional derivatives} We approximate the time-fractional derivatives of the macroscopic functions $U_i$ by a numerical scheme based on a finite-difference approach, which captures the non-local and memory effects inherent in fractional derivatives. One can find the stability and convergence of the time-fractional derivative difference scheme in \cite{alikhanov2015new, alikhanov2021high}.

Let $N_t$ be a count of time steps, and $\tau = t_{max} / N_t$ be a time step size. The time-fractional derivative difference scheme for has the following form
\begin{equation*}
  \frac{\partial^{\alpha} U_i^k}{\partial t^{\alpha}} \approx \xi^{(\alpha)}_{\tau} \left( U_i^k - U_i^{k-1} + \sum_{q=2}^{k}\xi^{(\alpha)}_{q-1} (U_i^{k-q+1}-U_i^{k-q})\right),
\end{equation*}
where $U_i^k \approx U_i(t_k)$, $t_k = k \tau$, and the coefficients $\xi^{(\alpha)}_{\tau}$ and $\xi^{(\alpha)}_{q-1}$ are defined as follows
\begin{equation*}
\xi^{(\alpha)}_{\tau} = \frac{1}{\tau^{\alpha} \Gamma (2-\alpha)}, \quad \xi^{(\alpha)}_{q-1} = q^{1-\alpha} - (q - 1)^{1-\alpha}.
\end{equation*}
In the next section, we present numerical results for the obtained dual-continuum model.


\section{Numerical results}\label{sec:numerical_results}

Let us consider a computational domain $\Omega = \Omega_1 \cup \Omega_2 = [0, 1] \times [0, 1]$ (see Figure \ref{fig:domain}). We construct a fine grid by partitioning $\Omega$ into $400 \times 400$ square and then dividing each square into two triangles. We consider two coarse grids constructed similarly but using $20 \times 20$ and $40 \times 40$ squares (coarse blocks). 

We take each coarse block as an RVE, while the oversampled RVE is an extension of the coarse block by $l$ layers of other coarse blocks. For choosing the number of layes, we use the formula $l = \lceil -2 \text{ log}(H) \rceil$, where $H$ is a coarse grid size \cite{efendiev2023multicontinuum}. Therefore, we have $l = 6$ for the coarse grid $20 \times 20$ and $l = 8$ for the coarse grid $40 \times 40$.

\begin{figure}[!hbt]
\centering
\includegraphics[width=0.35\textwidth]{figs/microstructure.png}
\caption{Computational domain ($\Omega_1$ is blue, $\Omega_2$ is red)}
\label{fig:domain}
\end{figure}

We take the following high-contrast heterogeneous coefficient $\kappa$
\begin{equation*}
\kappa = \begin{cases}
  10^{-4}, & x \in \Omega_1,\\
  1, & x \in \Omega_2.
\end{cases}
\end{equation*}

For boundary conditions, we set zero Dirichlet boundary conditions on all the boundaries. We set the duration of the simulated process as $t_{max} = 1$ with $N_t = 50$ time steps. For the initial condition, we set $u_0 = 0$. We set the following source term
\begin{equation*}
f = \text{exp}( -40 [ (x_1 - 0.5)^2 + (x_2 - 0.5)^2 ] ).
\end{equation*}

To check the effectivenes of the proposed multicontinuum approach, we define the following relative $L_2$ errors
\[
e_{2}^{(i)} = 
\sqrt{
\frac{\sum_K |\frac{1}{|K|} \int_K U_i dx - \frac{1}{|K \cap \Omega_i|} \int_{K \cap \Omega_i} u dx|^2}
{\sum_K |\frac{1}{|K \cap \Omega_i|} \int_{K \cap \Omega_i} u dx|^2}  
} \cdot 100\%, 
\]
where $i = 1, 2$, and $K$ is the RVE (taken to be $\omega$).

For spatial approximation, we use a finite element method with standard piecewise linear basis functions. The numerical implementation is based on the FEniCS computational package \cite{logg2012automated}. Visualization of the obtained results are performed using the ParaView software \cite{ahrens2005paraview}.

Figure \ref{fig:fine_res} depicts distributions of the fine-grid solution at different time steps for the fractional order derivative $\alpha = 0.8$. We observe an essential influence of the heterogeneous coefficient $\kappa$ on the solution. The solution field diffuses significantly faster in the high-conductivity regions ($\Omega_2$) than in the low-conductivity regions ($\Omega_1$). In general, we observe that the solution field gradually increases due to the source and diffuses over the domain.

\begin{figure}[!hbt]
  \centering
  \includegraphics[width=1.0\textwidth]{figs/u_fine.png}
  \caption{Distributions of the fine-grid solution at different time steps $t_k, k=1, 25, 50$ (from left to right) for the fractional order derivative $\alpha = 0.8$.}
  \label{fig:fine_res}
  \end{figure}

We present the effective properties computed for different coarse grids in Table \ref{tab:effective_properties}. One can see that the effective diffusion properties are isotropic. Hence, in our dual-continuum model \eqref{eq:dual_continuum_model}, we can take $\widehat{D_{ij}} = \widehat{A_{ij}^{11}} = \widehat{A_{ij}^{22}}$. Note that the reaction terms are consistent with the properties \eqref{eq:B_properties_2}. Thus, we can take $\widehat{B} = \widehat{B_{11}} = \widehat{B_{22}}$.

\begin{table}[h!]
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\multicolumn{3}{|c|}{Coarse grid $20 \times 20$} \\ 
\hline
$A_{11}^{11} / |R_\omega|$ & $A_{11}^{12} / |R_\omega|$ & $A_{11}^{22} / |R_\omega|$\\
\hline
0.000925  & $\approx 0$ & 0.000925 \\
\hline
$A_{12}^{11} / |R_\omega|$ & $A_{12}^{12} / |R_\omega|$ & $A_{12}^{22} / |R_\omega|$\\
\hline
-0.000759  & $\approx 0$ & -0.000759 \\
\hline
$A_{22}^{11} / |R_\omega|$ & $A_{22}^{12} / |R_\omega|$ & $A_{22}^{22} / |R_\omega|$\\
\hline
0.21946  & $\approx 0$ & 0.21946 \\
\hline
$B_{11} / |R_\omega|$ & $B_{12} / |R_\omega|$ & $B_{22} / |R_\omega|$\\
\hline
4.78677  & -4.78677 & 4.78677 \\
\hline
$C_{11} / |R_\omega|$ & $C_{12} / |R_\omega|$ & $C_{22} / |R_\omega|$\\
\hline
0.900799  & -0.260799 & 0.620799 \\
\hline
\end{tabular}
\quad 
\begin{tabular}{|c|c|c|}
\hline
\multicolumn{3}{|c|}{Coarse grid $40 \times 40$} \\ 
\hline
$A_{11}^{11} / |R_\omega|$ & $A_{11}^{12} / |R_\omega|$ & $A_{11}^{22} / |R_\omega|$\\
\hline
0.000463  & $\approx 0$ & 0.000463 \\
\hline
$A_{12}^{11} / |R_\omega|$ & $A_{12}^{12} / |R_\omega|$ & $A_{12}^{22} / |R_\omega|$\\
\hline
-0.000297  & $\approx 0$ & -0.000297 \\
\hline
$A_{22}^{11} / |R_\omega|$ & $A_{22}^{12} / |R_\omega|$ & $A_{22}^{22} / |R_\omega|$\\
\hline
0.219039  & $\approx 0$ & 0.219039 \\
\hline
$B_{11} / |R_\omega|$ & $B_{12} / |R_\omega|$ & $B_{22} / |R_\omega|$\\
\hline
4.78677  & -4.78677 & 4.78677 \\
\hline
$C_{11} / |R_\omega|$ & $C_{12} / |R_\omega|$ & $C_{22} / |R_\omega|$\\
\hline
0.900799  & -0.260799 & 0.620799 \\
\hline
\end{tabular}
\end{center}
\caption{Effective properties computed for different coarse grids.}
\label{tab:effective_properties}
\end{table}

Figures \ref{fig:res1}-\ref{fig:res2} present distributions of the numerical solution with the fractional order derivative $\alpha = 0.8$ at different time steps for the first and second continua, respectively. From top to bottom, we depict the reference averaged and multiscale averaged solutions, respectively. One can see that the results are very similar, indicating that our multicontinuum approach can provide high accuracy.

\begin{figure}[!hbt]
\centering
\includegraphics[width=1.0\textwidth]{figs/Uc_1}
\caption{Distributions of the solution for the first continuum at different time steps $t_k, k=1, 25, 50$ (from left to right) for the fractional order derivative $\alpha = 0.8$. First row: reference averaged solution. Second row: multiscale averaged solution. Coarse grid $40 \times 40$, relative $L_2$ errors: 2.1231\%, 1.8163\% and 1.8509\%.}
\label{fig:res1}
\end{figure}

In terms of the simulated process, it generally coincides with the dynamics of the fine-grid solution. We observe how the solution fields gradually increase due to the source in the middle and diffuse over the domain. As expected, the solution of the second continuum spreads faster than the first one. In terms of physical meaning, the solution of the first continuum represents the average in the low conductive region ($\Omega_1$), and the solution of the second continuum represents the average in the highly conductive region ($\Omega_2$).

\begin{figure}[!hbt]
\centering
\includegraphics[width=1.0\textwidth]{figs/Uc_2}
\caption{Distributions of the solution for the second continuum at different time steps $t_k, k=1, 25, 50$ (from left to right) for the fractional order derivative $\alpha = 0.8$. First row: reference averaged solution. Second row: multiscale averaged solution. Coarse grid $40 \times 40$, relative $L_2$ errors: 1.2931\%, 1.1787\% and 1.1654\%.}
\label{fig:res2}
\end{figure}

Tables \ref{tab:t20}-\ref{tab:t40} present relative $L_2$ errors (\%) for coarse grids $20 \times 20$ and $40 \times 40$, respectively. In both tables, we consider three cases of $\alpha$ (0.8, 0.9, and 1.0). Let us consider the errors of the coarse grid $20 \times 20$. We can see that the errors gradually decrease with time for all the fractional order derivatives. At the initial time steps, the higher order derivatives have slightly larger errors. However, after some time, the difference between different $\alpha$ cases becomes very small. At the final time, in all $\alpha$ cases, the relative $L_2$ errors for the first continuum are less than $1.86\%$ and the relative errors for the second continuum are less than $1.17\%$.

\begin{table}[h!]
\begin{center}
\begin{tabular}{c|c|c|c|c|c|c}
\hline
& \multicolumn{2}{c|}{$\alpha = 0.8$} & \multicolumn{2}{c|}{$\alpha = 0.9$} & \multicolumn{2}{c}{$\alpha = 1.0$}\\ 
\hline
k & $e_2^{(1)}$ & $e_2^{(2)}$ & $e_2^{(1)}$ & $e_2^{(2)}$ & $e_2^{(1)}$ & $e_2^{(2)}$\\
\hline
1  & 2.1231 & 1.2931 & 2.3511 & 1.3342 & 2.6531 & 1.3651 \\
10 & 1.8385 & 1.2524 & 1.8694 & 1.2995 & 1.9599 & 1.3611 \\
15 & 1.7997 & 1.2157 & 1.772  & 1.2504 & 1.7564 & 1.3001 \\
20 & 1.8046 & 1.1923 & 1.769  & 1.2127 & 1.7241 & 1.248  \\
25 & 1.8163 & 1.1787 & 1.787  & 1.1877 & 1.7428 & 1.2073 \\
30 & 1.8271 & 1.1713 & 1.8064 & 1.173  & 1.7727 & 1.1797 \\
35 & 1.8355 & 1.1674 & 1.8225 & 1.165  & 1.8006 & 1.1633 \\
40 & 1.8421 & 1.1657 & 1.8347 & 1.1613 & 1.8228 & 1.155  \\
45 & 1.8471 & 1.1652 & 1.8437 & 1.1601 & 1.8392 & 1.1518 \\
50 & 1.8509 & 1.1654 & 1.8502 & 1.1604 & 1.8508 & 1.1517 \\
\hline
\end{tabular}
\end{center}
\caption{Relative $L_2$ errors (\%) for the coarse grid $20 \times 20$. The index $k$ denotes the time step's number.}
\label{tab:t20}
\end{table}

Let us consider the relative $L_2$ errors (\%) for the coarse grid $40 \times 40$. One can see that the errors are significantly smaller than in the coarse grid $20 \times 20$. Again, we see that at the initial time steps, the higher order derivatives possesses larger errors. However, with more time steps, the difference becomes less noticeable. We observe the gradual decrease of the errors with time for all the cases of $\alpha$. At the final time, we have the relative $L_2$ errors for the first continuum less than $0.95\%$ and the relative $L_2$ errors for the second continuum less than $0.55\%$.

\begin{table}[h!]
  \begin{center}
  \begin{tabular}{|c|c|c|c|c|c|c|}
  \hline
  & \multicolumn{2}{c|}{$\alpha = 0.8$} & \multicolumn{2}{c|}{$\alpha = 0.9$} & \multicolumn{2}{c|}{$\alpha = 1.0$}\\ 
  \hline
k & $e_2^{(1)}$ & $e_2^{(2)}$ & $e_2^{(1)}$ & $e_2^{(2)}$ & $e_2^{(1)}$ & $e_2^{(2)}$\\
  \hline
 1  & 1.337  & 0.7067 & 1.5689 & 0.7838 & 1.876  & 0.9016 \\
 10 & 1.0533 & 0.6466 & 1.1047 & 0.6844 & 1.203  & 0.7377 \\
 15 & 0.9882 & 0.6131 & 0.9944 & 0.6435 & 1.0132 & 0.6851 \\
 20 & 0.9651 & 0.5888 & 0.9589 & 0.6091 & 0.9556 & 0.6425 \\
 25 & 0.9547 & 0.5724 & 0.9451 & 0.5834 & 0.9344 & 0.6042 \\
 30 & 0.9492 & 0.5614 & 0.9396 & 0.5658 & 0.9272 & 0.5755 \\
 35 & 0.9461 & 0.5538 & 0.9374 & 0.5540 & 0.9261 & 0.5562 \\
 40 & 0.9439 & 0.5485 & 0.9365 & 0.5462 & 0.9272 & 0.5436 \\
 45 & 0.9423 & 0.5446 & 0.9361 & 0.5408 & 0.9289 & 0.5357 \\
 50 & 0.9411 & 0.5418 & 0.9357 & 0.5372 & 0.9304 & 0.5307 \\
  \hline
  \end{tabular}
  \end{center}
  \caption{Relative $L_2$ errors (\%) for the coarse grid $40 \times 40$. The index $k$ denotes the time step's number.}
  \label{tab:t40}
  \end{table}

The obtained numerical results demonstrate that our proposed multicontinuum approach provides high accuracy with the significant reduction of the computational cost for different orders of time-fractional derivatives.

\section{Conclusion}\label{sec:conclusion}
In this paper, we have considered a time-fractional diffusion model based on Caputo fractional derivative. Using the multicontinuum homogenization approach, we have derived a multicontinuum time-fractional diffusion model. We have formulated special cell problems that account for the different homogenized effects of the solution. As a result, we have obtained a general multicontinuum model. Based on it, we have derived a particular case of a dual-continuum model in an isotropic medium.

To verify the proposed multicontinuum approach, we have performed numerical experiments. For this purpose, we have solved two-dimensional model problems for different fractional derivative order cases. The results demonstrate the high accuracy of the multicontinuum model with a significant reduction in computational costs. The proposed multicontinuum approach allows us to consider different cases of fractional derivative orders.

\begin{thebibliography}{10}

  \bibitem{caputo1967linear}
  Caputo M.,
  \href{https://doi.org/10.1111/j.1365-246X.1967.tb02303.x}{\it Linear models of dissipation whose q is almost frequency independent-ii}, Geophysical journal international, \textbf{13}:5, (1967), 529--539.
  
  \bibitem{metzler2000random}
  Metzler R., Klafter J., 
  \href{https://doi.org/10.1016/S0370-1573(00)00070-3}{\it The random walk's guide to anomalous diffusion: a fractional dynamics approach}, Physics reports, \textbf{339}:1, (2000), 1--77.
  
  \bibitem{podlubny1998fractional}
  Podlubny I., 
  \href{https://www.google.com/books?hl=ru&lr=&id=K5FdXohLto0C&oi=fnd&pg=PP1&ots=lKYGXKXmbR&sig=V7hml1qlr7I6a1eP4nafRkjqZ40}{\it Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications},
  Elsevier, 1998.
  
  \bibitem{mainardi2022fractional}
  Mainardi F., 
  \href{http://www.fracalmo.org/mainardi/MAINARDI_BOOK-COVER-PREFACE-CONTENTS.pdf}{\it Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models},
  World Scientific, 2022.
  
  \bibitem{hu2020homogenization}
  Hu J., Li G., 
  \href{https://doi.org/10.1016/j.jcp.2020.109231}{\it Homogenization of time-fractional diffusion equations with periodic coefficients},
  Journal of Computational Physics, \textbf{408}, 2020, 109231.
  
  \bibitem{kawamoto2023homogenization}
  Kawamoto A., Machida M., Yamamoto M., 
  \href{https://doi.org/10.1007/s13540-023-00195-8}{\it Homogenization and inverse problems for fractional diffusion equations},
  Fractional Calculus and Applied Analysis, \textbf{26}:5, (2023), 2118--2165.
  
  \bibitem{efendiev2004numerical}
  Efendiev Y., Pankov A., 
  \href{https://doi.org/10.1137/030600266}{\it Numerical homogenization of nonlinear random parabolic operators},
  Multiscale Modeling \& Simulation, \textbf{2}:2, (2004), 237--268.
  
  \bibitem{rubinvstein1948question}
  Rubin{\v{s}}te{\i}n L.I.,
  {\it On a question about the propagation of heat in heterogeneous media}.(russian) izvestiya akad.,
  Nauk SSSR. Ser. Geograf. Geofiz, \textbf{12}:1, (1948), 27--45.
  
  \bibitem{barenblatt1960basic}
  Barenblatt G. I., Zheltov I. P., Kochina I. N., 
  \href{https://blasingame.engr.tamu.edu/z_zCourse_Archive/P620_14C/P620_14C_zReference/(Barenblatt_Zheltov_Kochina)_Basic_Concepts_in_the_Theory_of_Seepage_of_Homogeneous_Liquids_in_Fissured_Rocks.pdf}{\it Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]},
  Journal of applied mathematics and mechanics, \textbf{24}:5, (1960), 1286--1303.
  
  \bibitem{arbogast1990derivation}
  Arbogast T., Douglas, Jr J., Hornung U., 
  \href{https://doi.org/10.1137/0521046}{\it
  Derivation of the double porosity model of single phase flow via homogenization theory},
  SIAM Journal on Mathematical Analysis, \textbf{21}:4, (1990), 823--836.
  
  \bibitem{panasenko2018multicontinuum}
  Panasenko G. P.,
  \href{https://doi.org/10.1007/s10958-018-3889-7}{\it
  Multicontinuum wave propagation in a laminated beam with contrasting stiffness and density of layers},
  Journal of Mathematical Sciences, \textbf{232}, (2018), 503--515.
  
  \bibitem{alekseev2023partial}
  Alekseev V. N., Kalachikova U. S., Yang Y.,
  \href{https://doi.org/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{efendiev2023multicontinuum}
  Efendiev Y., Leung W. T., 
  \href{https://doi.org/10.1016/j.jcp.2022.111761}{\it 
  Multicontinuum homogenization and its relation to nonlocal multicontinuum theories},
  Journal of Computational Physics, \textbf{474}, (2023), 111761.
  
  \bibitem{chung2024multicontinuum}
  Chung E., Efendiev Y., Galvis J., Leung W. T.,
  \href{https://doi.org/10.1016/j.jcp.2024.112980}{\it
  Multicontinuum homogenization. general theory and applications},
  Journal of Computational Physics, \textbf{510}, (2024), 112980.
  
  \bibitem{leung2024some}
  Leung W. T.,
  \href{https://doi.org/10.48550/arXiv.2401.12799}{\it
  Some convergence analysis for multicontinuum homogenization},
  arXiv preprint arXiv:2401.12799, (2024).
  
  \bibitem{xie2024multicontinuum}
  Xie W., Efendiev Y., Huang Y., Leung W. T., Yang Y., 
  \href{https://doi.org/10.48550/arXiv.2404.17471}{\it
  Multicontinuum homogenization in perforated domains},
  arXiv preprint arXiv:2404.17471, (2024).
  
  \bibitem{ammosov2024multicontinuum}
  Ammosov D., Leung W. T., Shan B., Huang J.,
  \href{https://doi.org/10.48550/arXiv.2405.14572}{\it
  Multicontinuum homogenization for coupled flow and transport equations},
  arXiv preprint arXiv:2405.14572, (2024).
  
  \bibitem{ammosov2023multicontinuum}
  Ammosov D., Stepanov S., Spiridonov D., Li W.,
  \href{https://doi.org/10.1515/rnam-2023-0016}{\it
  Multicontinuum homogenization for richards’ equation: The derivation and numerical experiments},
  Russian Journal of Numerical Analysis and Mathematical Modelling, \textbf{38}:4, (2023), 207--218.
  
  \bibitem{kalachikova2024advancing}
  Kalachikova U., Ammosov D.,
  \href{https://doi.org/10.1016/j.cam.2024.115755}{\it
  Advancing wave equation analysis in dual-continuum systems: A partial learning approach with discrete empirical interpolation and deep neural networks},
  Journal of Computational and Applied Mathematics, \textbf{443}, (2024), 115755.
  
  \bibitem{tyrylgin2023computational}
  Tyrylgin A., Vasilyeva M., Alikhanov A., Sheen D.,
  \href{https://doi.org/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, \textbf{418}, (2023), 114670.
  
  \bibitem{alikhanov2025multiscale}
  Alikhanov A., Bai H., Huang J., Tyrylgin A., Yang Y.,
  \href{https://doi.org/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, \textbf{455}, (2025), 116157.
  
  \bibitem{li2022partially}
  Li W., Alikhanov A., Efendiev Y., Leung W. T.,
  \href{https://doi.org/10.1016/j.cnsns.2022.106440}{\it
  Partially explicit time discretization for nonlinear time fractional diffusion equations},
  Communications in Nonlinear Science and Numerical Simulation, \textbf{113}, (2022), 106440.
  
  \bibitem{hu2022partially}
  Hu J., Alikhanov A., Efendiev Y., Leung W. T.,
  \href{https://doi.org/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{alikhanov2015new}
  Alikhanov A.A.,
  \href{https://doi.org/10.1016/j.jcp.2014.09.031}{\it
  A new difference scheme for the time fractional diffusion equation},
  Journal of Computational Physics, \textbf{280}, (2015), 424--438.
  
  \bibitem{alikhanov2021high}
  Alikhanov A. A., Huang C.,
  \href{https://doi.org/10.1016/j.amc.2021.126545}{\it
  A high-order l2 type difference scheme for the time-fractional diffusion equation},
  Applied Mathematics and Computation, \textbf{411}, (2021), 126545.
  
  \bibitem{logg2012automated}
  Logg A., Mardal K. A., Wells G.,
  \href{https://books.google.ru/books?hl=ru&lr=&id=ASWN_VRr1NQC&oi=fnd&pg=PR3&dq=Anders+Logg,+Kent-Andre+Mardal,+and+Garth+Wells.+++Automated+solution+of+differential+equations+by+the+finite+element+method&ots=F_7XrMApW8&sig=Ed87-aWBs__tqlbaI6o38_hv_d0&redir_esc=y#v=onepage&q=Anders%20Logg%2C%20Kent-Andre%20Mardal%2C%20and%20Garth%20Wells.%20%20%20Automated%20solution%20of%20differential%20equations%20by%20the%20finite%20element%20method&f=false}{\it
  Automated solution of differential equations by the finite element method: The FEniCS book},
  Springer Science \& Business Media, \textbf{84}, (2012).
  
  \bibitem{ahrens2005paraview}
  Ahrens J., Geveci B., Law C., Hansen C., Johnson C.,
  \href{https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=37289ca24fde4c63281c061815d625c7acee1e6a#page=736}{\it
  Paraview: An end-user tool for large data visualization.
  The visualization handbook}, \textbf{717}:8, (2005).
  
  \end{thebibliography}
  
  
\end{document}
