Multicontinuum Homogenization for Coupled Flow and Transport Equations

Abstract

In this paper, we present the derivation of a multicontinuum model for the coupled flow and transport equations by applying multicontinuum homogenization. We perform the multicontinuum expansion for both flow and transport solutions and formulate novel coupled constraint cell problems to capture the multiscale property, where oversampled regions are utilized to avoid boundary effects. Assuming the smoothness of macroscopic variables, we obtain a multicontinuum system composed of macroscopic elliptic equations and convection-diffusion-reaction equations with homogenized effective properties. Finally, we present numerical results for various coefficient fields and boundary conditions to validate our proposed algorithm.

Keywords: multiscale, multicontinuum, homogenization, coupled flow and transport, miscible displacement.

1 Introduction

The coupled system of flow and transport equations is fundamental in modeling phenomena across various scientific and engineering fields, including hydrology, petroleum engineering, and environmental science. For instance, the miscible displacement of one incompressible fluid by another, which occurs in processes ranging from enhanced oil and gas recovery to the remediation of contaminated groundwater, can be governed by these equations [1]. However, a ubiquitous challenge in these scenarios is the inherent heterogeneities and high contrast in material properties. These complexities not only complicate simulation efforts but also significantly impact the accuracy of predictive models. Consequently, the adoption of upscaling techniques is necessary to accurately capture detailed multiscale properties, enabling time-saving and precise computations.

Many different methods have been developed to perform simulations on a coarse grid, one type of which is homogenization methods [2, 3, 4]. In homogenization methods, one defines upscaled coefficients for multiscale properties through local solutions of fine-grid problems and utilizes them to formulate the macroscopic flow and transport equations on the coarse grid. The homogenized solution for the transport equation can be computed after the macroscopic flow equation has been solved. Such homogenization methods have been applied in petroleum engineering [5, 6].

Another approach is multiscale methods, including the Multiscale Finite Element Method [7, 8, 9], the Multiscale Finite Volume Method [10, 11, 12], the Generalized Multiscale Finite Element Method [13, 14, 15], and others. Multiscale methods construct local basis functions via fine-grid local problems and use them to approximate the solutions. One could either apply multiscale methods to the flow equation and solve the transport equation on the fine grid, or apply these methods to both of them, depending on the media properties [11, 16]. Mixed multiscale methods have also been developed to solve the coarse transport equation [17, 18].

This paper presents the multicontinuum homogenization method following earlier works on the elliptic equation [19] using the ideas of the above two types of methods. For deriving a macroscopic equation, multicontinuum homogenization posits a multicontinuum expansion in each coarse block (or Representative Volume Element (RVE) in the case of scale separation) about macroscopic quantities and multiscale basis functions. The macroscopic quantities are introduced to represent local averages of solution in each continuum and assumed to be smooth over coarse regions. Besides, the multiscale basis functions are obtained from some oversampled locally constraint cell problems solved on a fine grid, which are inspired by [20, 21] and can be designed to promise localization.

In our study, we will apply multicontinuum homogenization to both the flow and the transport equations in the coupled system. However, challenges are that the transport equation contains heterogeneous and high-contrast convection and temporal terms, and its convection coefficient further depends on the solution of the flow equation. To address these problems, we carefully design new coupled constraint cell problems, where the one for transport reflects the multiscale property for each term. With multiscale basis functions obtained from the coupled cell problems, we perform the multicontinuum expansions in the variational forms of the flow and transport equations. We will derive a multicontinuum system composed of macroscopic elliptic equations and convection-diffusion-reaction equations, where the effective properties have been homogenized in each coarse block to capture the fine-scale details. We remark here that the properties for convection-diffusion-reaction equations will depend on the flow equation. Therefore, machine learning techniques could be used to predict effective properties, and we could, as in many model reduction methods, divide the computation into offline and online stages to save online computation. We also note that the fine-scale information can be recovered from macroscopic solutions by the multicontinuum expansions. We apply it to coupled flow and transport equations with different coefficient fields subject to different boundary conditions to validate our proposed method. We change the coarse mesh size and compare the macroscopic results with the averaged reference solutions. The relative $L^{2}$ errors are very small, and our method performs the simulations accurately.

Our main contributions consist of the following. Firstly, we formulate appropriate coupled constraint cell problems to deal with the multiscale property of the coupled equations. Secondly, we derive a novel multicontinuum coupled model about macroscopic solutions on the coarse grid. Thirdly, we present numerical results for various media fields and boundary conditions to show the validity.

The paper is organized as follows. In Section 2, we present preliminaries and review some earlier works on multiscale methods. Section 3 is contributed to the description of multicontinuum homogenization and the derivation of the multicontinuum model. We present numerical results in Section 4 and finally draw conclusions in Section 5.

2 Preliminaries

In this work, we consider the following linearized and one-way coupled system of flow and transport equations for pressure $p$ and concentration $c$ in a bounded domain $\Omega\subset\mathbb{R}^{d}$

$$\begin{split}&-\nabla\cdot(\kappa\nabla p)=g,\\ &\phi\frac{\partial c}{\partial t}+u\cdot\nabla c-\nabla\cdot(D\nabla c)=h(c),\end{split}$$ (1)

where $u=-\kappa\nabla p$, $\kappa$ is a permeability coefficient, $\phi$ is a porosity, $D$ is a diffusion coefficient, and $g$ and $h$ are source terms. We assume $\kappa$ and $D$ are heterogeneous and of high contrast. Before delving into multicontinuum homogenization, we briefly review earlier efforts on some methods for multiscale problems, including (numerical) homogenization, MsFEM, GMsFEM.

The main idea of the homogenization method is to derive macroscopic equations, with the multiscale coefficients in original equations replaced by homogenized effective properties. While the classical homogenization method relies on periodicity and posits a formal asymptotic expansion of solution, numerical homogenization defines an upscaled property for each coarse block, conversing averages in a physical sense for specific boundary conditions, without the periodicity assumption. To homogenize the flow equation in (1), for instance, we formulate local problems for each coarse block $K$

$$-\nabla\cdot(\kappa\nabla\phi_{l})=0,$$ (2)

subject to some appropriate boundary conditions, like $\phi_{l}=x_{l}$ on $\partial K$. The elements of the upscaled coefficient $\kappa^{*}$ are defined by averaging the fluxes

$$\int_{K}\kappa_{ij}^{*}\frac{\partial}{\partial x_{j}}\phi_{l}^{*}=\int_{K}\kappa_{ij}(x)\frac{\partial}{\partial x_{j}}\phi_{l},$$ (3)

and this leads to

$$\kappa_{\cdot,l}^{*}=\frac{1}{|K|}\int_{K}\kappa(x)\nabla\phi_{l}.$$ (4)

Then, the transport equation can be homogenized with a similar idea to get upscaled diffusion $D^{*}$ while the upscaled porosity $\phi^{*}$ would be its volume averaging, and consequently, we derive the following coupled macroscopic equations

$$\begin{split}&-\nabla\cdot(\kappa^{*}\nabla p^{*})=g,\\ &\phi^{*}\frac{\partial c^{*}}{\partial t}+u^{*}\cdot\nabla c^{*}-\nabla\cdot(D^{*}\nabla c^{*})=h,\end{split}$$ (5)

where $u^{*}=-\kappa^{*}\nabla p^{*}$, $\kappa^{*}$ and $D^{*}$ are anisotropic upscaled coefficients. The macroscopic solutions can be obtained by solving the coupled system above on the coarse grid.

The Multiscale Finite Element Method (MsFEM), rather than deriving a macroscopic equation, establishes local multiscale basis functions via local problems formulated on coarse blocks to approximate the solution. Specifically, with multiscale basis functions $\phi^{\omega_{i},p}$ for pressure and $\phi^{\omega_{i},c}$ for concentration at each coarse node $i$ supported in coarse neighborhood $\omega_{i}$, we seek approximate solutions $p_{H}$ and $c_{H}$ that lie in $V_{\text{ms}}^{p}=\text{span}_{i}\{\phi^{\omega_{i},p}\}$ and $V_{\text{ms}}^{c}=\text{span}_{i}\{\phi^{\omega_{i},c}\}$ respectively, and

$$p_{H}=\sum_{i}p_{i}\phi^{\omega_{i},p},\quad c_{H}=\sum_{i}c_{i}\phi^{\omega_{i},c}.$$ (6)

Now, we illustrate how to construct local problems. Since the transport equation in (1) lacks temporal heterogeneity, $\phi^{\omega_{i},p}$ and $\phi^{\omega_{i},c}$ can be chosen to satisfy

$$\begin{split}&-\nabla\cdot(\kappa\nabla\phi^{\omega_{i},p})=0,\\ &(-\kappa\sum_{i}p_{i}\nabla\phi^{\omega_{i},p})\cdot\nabla\phi^{\omega_{i},c}-\nabla\cdot(D\nabla\phi^{\omega_{i},c})=0\end{split}$$ (7)

in each coarse block K such that $\phi^{\omega_{i},p}=\phi^{\omega_{i},c}=\phi_{0}^{\omega_{i}}$ on $\partial K$, where $\phi_{0}^{\omega_{i}}$ represents the nodal basis at node $i$ of the coarse finite element space. Indeed, the mesh grids for pressure and concentration are not necessarily the same, and the boundary conditions for local problems can be chosen differently. Note also that local problems can be defined in representative volume element (RVE) in the case of scale separation, and fine-scale information could be recovered, which is unachievable by most homogenization methods.

The Generalized Multiscale Finite Element Method (GMsFEM), evolving from MsFEM, adds additional local degrees of freedom as necessary and thus can systematically enrich coarse spaces and capture multiscale features more precisely. The first step is to construct the snapshot space $V_{\text{snap }}^{(i),p}$ and $V_{\text{snap }}^{(i),c}$ for each generic coarse region $\omega_{i}$, which can be composed of, for example, harmonic extensions of fine-grid functions on $\partial\omega^{i}$, that is,

$$\begin{split}&-\nabla\cdot(\kappa\nabla\psi_{l}^{\omega_{i},p})=0,\\ &(-\kappa\nabla p_{H})\cdot\nabla\psi_{l}^{\omega_{i},c}-\nabla\cdot(D\nabla\psi_{l}^{\omega_{i},c})=0\end{split}$$ (8)

in $\omega_{i}$ with $\psi_{l}^{\omega_{i},p}=\psi_{l}^{\omega_{i},c}=\delta_{l}^{h}$ on $\partial\omega^{i}$. Here $\delta_{l}^{h}$ is a piecewise linear function on $\partial\omega^{i}$ satisfying $\delta_{l}^{h}(x_{k})=\delta_{l,k}$ for any fine-grid boundary node $x_{k}$. Subsequently, we only keep the dominant component subspace of each local snapshot space and form a multiscale space. A local spectral decomposition for each snapshot space $V_{\text{snap }}^{(i)}$ is constructed for model reduction purposes. This involves finding $\lambda_{j}^{\omega_{i}}\in\mathbb{R}$ and $\phi_{j}^{\omega_{i}}\in V_{\text{snap }}^{(i)}$ such that

$$a_{\omega_{i}}(\phi_{j}^{\omega_{i}},w)=\lambda_{j}^{\omega_{i}}s_{\omega_{i}}(\phi_{j}^{\omega_{i}},w),\quad\forall w\in V_{\text{snap }}^{(i)}.$$ (9)

We select the eigenvectors corresponding to the smallest eigenvalues, denoted as $\phi_{l}^{\omega_{i},p}$ and $\phi_{m}^{\omega_{i},c}$, and multiply them by multiscale partition of unity $\chi_{i}$ to form the multiscale basis, i.e., $V_{\text{ms}}^{p}=\text{span}_{i,l}\{\chi_{i}\phi_{l}^{\omega_{i},p}\}$ and $V_{\text{ms}}^{c}=\text{span}_{i,m}\{\chi_{i}\phi_{m}^{\omega_{i},c}\}$. We remark that the degrees of freedom can be added adaptively according to a-posteriori error estimates. For more details on the methods above, one may refer to [4, 7, 13, 22].

3 Multicontinuum homogenization and multicontinuum model

In this section, we will present the main ideas and concepts of multicontinuum homogenization and apply it to the coupled flow and transport equations to derive a multicontinuum model.

3.1 Main ideas

The multicontinuum homogenization method, similar to classical homogenization, also aims for macroscopic equations. It introduces macroscopic quantities to represent the local averages of solution and posits the general multicontinuum expansion in each local region, like RVE,

$$u=\phi_{i}U_{i}+\phi_{i}^{m}\nabla_{m}U_{i}+\phi_{i}^{mn}\nabla_{mn}^{2}U_{i}+\cdots,$$ (10)

where $U_{i}$ is the macroscopic quantity for continuum $i$ with physical meaning, $\phi_{i}$, $\phi_{i}^{m}$, $\phi_{i}^{mn}$, $\ldots$ are auxiliary basis functions obtained via some carefully constructed cell problems, which minimize local energy under certain constraints. The summation over repeated indices is taken. Here, the choice of cell problems is crucial for ensuring localization of the local basis. In the following, we will ignore the terms after the second term and only use the two-term expansion for simplicity. We apply the above expansion to both the solution and test function and substitute them into the variational form, where RVE concepts could be utilized. Using the smoothness of macroscopic quantities, we can take them out of the local integrals and, after some calculation, arrive at a macroscopic equation about $U_{i}$, which can recover the fine-scale information by expansion (10). Note that the smoothness of macroscopic quantities is reasonable and crucial for multicontinuum homogenization. We also remark that the continua are not necessarily defined apriori but could be achieved by local spectral problems.

Next, we will describe the multicontinuum homogenization method for (1) in detail. Since the pressure field’s equation does not depend on the concentration, we can perform homogenization in a split way.

3.2 Flow homogenization

Let us consider the following variational formulation of the flow equation

$$\int_{\Omega}\kappa\nabla p\cdot\nabla q=\int_{\Omega}gq,\quad\forall q\in H_{0}^{1}(\Omega).$$ (11)

Cell problems.

We assume that our computational domain $\Omega$ is partitioned into coarse-grid elements $\omega$’s, whose sizes are larger than the scale of heterogeneities. Also, we assume there is a Representative Volume Element (RVE) $R_{\omega}$ inside each coarse element $\omega$. These RVEs have the following properties:

• $R_{\omega}$ can represent the entire $\omega$ in terms of its heterogeneities;

• $R_{\omega}$ contains N continua (or components). For each continuum of $R_{\omega}$, we introduce the following characteristic function $\psi_{j}$ ($j=1,...,N$)

$$\psi_{j}=\begin{cases}1,&\text{in continuum}\leavevmode\nobreak\ j,\\ 0,&\text{otherwise}.\end{cases}$$

Figure 1: Illustration of computational domain $\Omega$, coarse block $\omega$, RVE $R_{\omega}$ and oversampled RVE $R_{\omega}^{+}$

For performing homogenization, we introduce two cell problems in the oversampled RVE $R_{\omega}^{+}$ constructed around $R_{\omega}$ and consisted of several RVEs $R_{\omega}^{k}$, where $R_{\omega}^{k_{0}}=R_{\omega}$. The relationship between the computational domain $\Omega$, coarse block $\omega$, RVE $R_{\omega}$ and oversampled RVE $R_{\omega}^{+}$ has been illustrated in Figure 1. We would like to minimize the local energy in ${R_{\omega}^{+}}$ under different constraints, where Lagrange multipliers are applied. The first cell problem (12) considers gradient effects and imposes constraints to represent the linear functions in the average behavior of each continuum

$$\begin{split}&\int_{R_{\omega}^{+}}\kappa\nabla\varphi^{m,p}_{i}\cdot\nabla q-\sum_{j,k}{\beta_{ij}^{mk}\over\int_{R_{\omega}^{k}}\psi_{j}^{k}}\int_{R_{\omega}^{k}}\psi_{j}^{k}q=0,\\ &\int_{R_{\omega}^{k}}\varphi^{m,p}_{i}\psi_{j}^{k}=\delta_{ij}\int_{R_{\omega}^{k}}(x_{m}-\tilde{x}_{m})\psi_{j}^{k},\\ &\int_{R_{\omega}^{k_{0}}}(x_{m}-\tilde{x}_{m})\psi_{j}^{k_{0}}=0,\\ \end{split}$$ (12)

where $\tilde{x}_{m}$ is a constant. The second cell problem considers different averages in each continuum and represents the constants in the average behavior

$$\begin{split}&\int_{R_{\omega}^{+}}\kappa\nabla\varphi_{i}^{p}\cdot\nabla q-\sum_{j,k}{\beta_{ij}^{k}\over\int_{R_{\omega}^{k}}\psi_{j}^{k}}\int_{R_{\omega}^{k}}\psi_{j}^{k}q=0,\\ &\int_{R_{\omega}^{k}}\varphi_{i}^{p}\psi_{j}^{k}=\delta_{ij}\int_{R_{\omega}^{k}}\psi_{j}^{k}.\end{split}$$ (13)

We remark that the solution in $R_{\omega}$ should be independent of the oversampling size as a result of the construction of our cell problems.

Coarse-grid model.

Next, we introduce the following multicontinuum expansion of $p$ in each $R_{\omega}$

$$\begin{split}p\approx\varphi_{i}^{p}P_{i}+\varphi_{i}^{m,p}\nabla_{m}P_{i},\end{split}$$ (14)

where $P_{i}$ is a smooth function in a macroscopic sense representing the homogenized solution for the $i$-th continuum. In our case, we have $P_{i}(x_{\omega}^{*})=\dfrac{\int_{R_{\omega}}p\psi_{i}}{\int_{R_{\omega}}\psi_{i}}$ for some $x_{\omega}^{*}\in R_{\omega}$.

Using the definition of RVE, we obtain the following approximation of the variational formulation

$$\begin{split}&\int_{\Omega}gq=\int_{\Omega}\kappa\nabla p\cdot\nabla q=\sum_{\omega}\int_{\omega}\kappa\nabla p\cdot\nabla q\approx\sum_{\omega}{|\omega|\over|R_{\omega}|}\int_{R_{\omega}}\kappa\nabla p\cdot\nabla q.\end{split}$$ (15)

Next, if we replace $p$ with its multicontinuum expansion (14), we get

$$\begin{split}\int_{R_{\omega}}\kappa\nabla p\cdot\nabla q\approx\int_{R_{\omega}}\kappa\nabla(\varphi_{i}^{p}P_{i})\cdot\nabla q+\int_{R_{\omega}}\kappa\nabla(\varphi_{i}^{m,p}\nabla_{m}P_{i})\cdot\nabla q.\end{split}$$ (16)

Additionally, we define the following multicontinuum expansion of $q$

$$q\approx\varphi_{s}^{p}Q_{s}+\varphi_{s}^{k,p}\nabla_{k}Q_{s}.$$ (17)

According to our assumption on the smoothness of the macroscopic variables, the variations of $P_{i}$ and $\nabla_{m}P_{i}$ are relatively small when compared to the variations of $\varphi_{i}^{p}$ and $\varphi_{i}^{m,p}$. Similarly, this applies to $Q_{s}$ and $\nabla_{k}Q_{s}$.

Then, we can approximate the first term of (16) as follows

$$\begin{split}&\int_{R_{\omega}}\kappa\nabla(\varphi_{i}^{p}P_{i})\cdot\nabla q\approx P_{i}(x_{\omega})\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{p}\cdot\nabla q\approx\\ &P_{i}(x_{\omega})Q_{s}(x_{\omega})\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{p}\cdot\nabla\varphi_{s}^{p}+P_{i}(x_{\omega})\nabla_{m}Q_{s}(x_{\omega})\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{p}\cdot\nabla\varphi_{s}^{m,p}=\\ &P_{i}(x_{\omega})\beta_{is}^{*}Q_{s}(x_{\omega})+\beta_{is}^{m*}P_{i}(x_{\omega})\nabla Q_{s}(x_{\omega}),\end{split}$$ (18)

where

$$\beta_{ik}^{m*}=\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{p}\cdot\nabla\varphi_{k}^{m,p},\ \ \beta_{ik}^{*}=\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{p}\cdot\nabla\varphi_{k}^{p}.$$ (19)

The second term of (16) is approximated in a similar way

$$\begin{split}&\int_{R_{\omega}}\kappa\nabla(\varphi_{i}^{m,p}\nabla_{m}P_{i})\cdot\nabla q\approx\nabla_{m}P_{i}(x_{\omega})\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{m,p}\cdot\nabla q\approx\\ &\nabla_{m}P_{i}(x_{\omega})\nabla_{k}Q_{s}(x_{\omega})\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{m,p}\cdot\nabla\varphi_{s}^{k,p}+\nabla_{m}P_{i}(x_{\omega})Q_{s}(x_{\omega})\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{m,p}\cdot\nabla\varphi_{s}^{p}=\\ &\nabla_{m}P_{i}(x_{\omega})\nabla_{k}Q_{s}(x_{\omega})\alpha_{is}^{km}+\nabla_{m}P_{i}(x_{\omega})Q_{s}(x_{\omega})\beta_{is}^{m*},\end{split}$$ (20)

where

$$\alpha_{is}^{km}=\int_{R_{\omega}}\kappa\nabla\varphi_{i}^{m,p}\cdot\nabla\varphi_{s}^{k,p}.$$ (21)

Therefore, by utilizing continuous approximations for $P_{i}$ and $Q_{j}$, we have

$$\begin{split}\int_{R_{\omega}}\kappa\nabla p\cdot\nabla q\approx P_{i}\beta_{ij}^{n*}\nabla_{n}Q_{j}+P_{i}\beta_{ij}^{*}Q_{j}+\nabla_{m}P_{i}\alpha_{ij}^{mn}\nabla_{n}Q_{j}+\nabla_{m}P_{i}\beta_{ij}^{m*}Q_{j}.\end{split}$$ (22)

According to [19], we have the following estimates for the solutions of the local cell problems and their gradients

$$\begin{split}&\|\varphi_{i}^{p}\|=O(1),\ \|\nabla\varphi_{i}^{p}\|=O({1\over\epsilon}),\\ &\|\varphi_{i}^{m,p}\|=O(\epsilon),\ \|\nabla\varphi_{i}^{m,p}\|=O(1),\\ \end{split}$$ (23)

where $\epsilon$ is the size of RVE. Then, we can obtain the following scalings

$$\beta_{ij}^{*}=O(\cfrac{|R_{\omega}|}{\epsilon^{2}}),\ \beta_{ij}^{m*}=O(\cfrac{|R_{\omega}|}{\epsilon}),\ \alpha_{ij}^{mn}=O(|R_{\omega}|).\$$

We introduce $\widehat{\beta_{ij}}$, $\widehat{\beta_{ij}^{m*}}$, and $\widehat{\alpha_{ij}^{mn}}$ as follows

$$\begin{split}\widehat{\beta_{ij}}=\cfrac{\epsilon^{2}}{|R_{\omega}|}\beta_{ij}^{*},\ \widehat{\beta_{ij}^{m*}}=\cfrac{\epsilon}{|R_{\omega}|}\beta_{ij}^{m*},\ \widehat{\alpha_{ij}^{mn}}=\frac{1}{|R_{\omega}|}\alpha_{ij}^{mn}.\ \end{split}$$ (24)

Using these scalings, we obtain

$$\begin{split}\int_{\Omega}\kappa\nabla p\cdot\nabla q\approx\int_{\Omega}\widehat{\alpha_{ij}^{mn}}\nabla_{m}P_{i}\nabla_{n}Q_{j}+{1\over\epsilon}\int_{\Omega}\widehat{\beta_{ij}^{m}}\nabla_{m}P_{i}Q_{j}+{1\over\epsilon}\int_{\Omega}\widehat{\beta_{ij}^{m}}P_{i}\nabla_{m}Q_{j}+{1\over\epsilon^{2}}\int_{\Omega}\widehat{\beta_{ij}}P_{i}Q_{j}.\end{split}$$ (25)

By applying the integration by parts formula, one can find out that that the sum of the second and third terms is negligible. Finally, we obtain the following multicontinuum flow equations

$$\begin{split}-\nabla_{n}(\widehat{\alpha_{ij}^{mn}}\nabla_{m}P_{j})+{1\over\epsilon^{2}}\widehat{\beta_{ij}}P_{j}=g_{i}.\end{split}$$ (26)

3.3 Transport homogenization

Let us consider the following variational formulation of the transport equation

$$\int_{\Omega}\phi\frac{\partial c}{\partial t}v+\int_{\Omega}(u\cdot\nabla c)v+\int_{\Omega}D\nabla c\cdot\nabla v=\int_{\Omega}h(c)v,\quad\forall v\in H_{0}^{1}(\Omega).$$ (27)

Cell problems.

To formulate and solve the cell problems, we need to linearize the convective term using known $P_{i}$. One can obtain these functions by solving the flow problem prior to that. Another approach is to consider many possible variants of $P_{i}$ to approximate the effective properties using machine learning techniques.

Again, the first cell problem (28) takes into account gradient effects

$$\begin{split}&\int_{R_{\omega}^{+}}\left[-\kappa\nabla(\varphi_{i}^{p}P_{i}+\varphi_{i}^{m,p}\nabla_{m}P_{i})\right]\cdot\nabla\varphi^{m,c}_{i}v+\int_{R_{\omega}^{+}}D\nabla\varphi^{m,c}_{i}\cdot\nabla v-\sum_{j,k}{\theta_{ij}^{mk}\over\int_{R_{\omega}^{k}}\psi_{j}^{k}}\int_{R_{\omega}^{k}}\psi_{j}^{k}v=0,\\ &\int_{R_{\omega}^{k}}\varphi^{m,c}_{i}\psi_{j}^{k}=\delta_{ij}\int_{R_{\omega}^{k}}(x_{m}-\tilde{x}_{m})\psi_{j}^{k},\\ &\int_{R_{\omega}^{k_{0}}}(x_{m}-\tilde{x}_{m})\psi_{j}^{k_{0}}=0.\\ \end{split}$$ (28)

The second cell problem considers different averages in each continuum

$$\begin{split}&\int_{R_{\omega}^{+}}\left[-\kappa\nabla(\varphi_{i}^{p}P_{i}+\varphi_{i}^{m,p}\nabla_{m}P_{i})\right]\cdot\nabla\varphi^{c}_{i}v+\int_{R_{\omega}^{+}}D\nabla\varphi_{i}^{c}\cdot\nabla v-\sum_{j,k}{\theta_{ij}^{k}\over\int_{R_{\omega}^{k}}\psi_{j}^{k}}\int_{R_{\omega}^{k}}\psi_{j}^{k}v=0,\\ &\int_{R_{\omega}^{k}}\varphi_{i}^{c}\psi_{j}^{k}=\delta_{ij}\int_{R_{\omega}^{k}}\psi_{j}^{k}.\end{split}$$ (29)

Coarse-grid model.

We apply the following multicontinuum expansion of $c$ in $R_{\omega}$ using $\varphi_{i}^{c}$ and $\varphi_{i}^{m,c}$ as follows

$$\begin{split}c\approx\varphi_{i}^{c}C_{i}+\varphi_{i}^{m,c}\nabla_{m}C_{i},\end{split}$$ (30)

where $C_{i}$ is a smooth function representing the homogenized solution for the $i$-th continuum.

We obtain the following approximation of the variational formulation

$$\begin{split}&\int_{\Omega}\phi{\partial c\over\partial t}v+\int_{\Omega}(u\cdot\nabla c)v+\int_{\Omega}D\nabla c\cdot\nabla v=\sum_{\omega}\int_{\omega}\phi{\partial c\over\partial t}v+\sum_{\omega}\int_{\omega}(u\cdot\nabla c)v+\sum_{\omega}\int_{\omega}D\nabla c\cdot\nabla v\approx\\ &\sum_{\omega}{|\omega|\over|R_{\omega}|}\int_{R_{\omega}}\phi{\partial c\over\partial t}v+\sum_{\omega}{|\omega|\over|R_{\omega}|}\int_{R_{\omega}}(u\cdot\nabla c)v+\sum_{\omega}{|\omega|\over|R_{\omega}|}\int_{R_{\omega}}D\nabla c\cdot\nabla v=\sum_{\omega}{|\omega|\over|R_{\omega}|}\int_{R_{\omega}}h(c)v\end{split}$$ (31)

The homogenization of the diffusion term is similar to the flow case. We can represent it as follows

$$\begin{split}\int_{R_{\omega}}D\nabla c\cdot\nabla v\approx C_{i}\theta_{ij}^{n*}\nabla_{n}V_{j}+C_{i}\theta_{ij}^{*}V_{j}+\nabla_{m}C_{i}\eta_{ij}^{mn}\nabla_{n}V_{j}+\nabla_{m}C_{i}\theta_{ij}^{m*}V_{j},\end{split}$$ (32)

where

$$\eta_{is}^{km}=\int_{R_{\omega}}D\nabla\varphi_{i}^{m,c}\cdot\nabla\varphi_{s}^{k,c},\ \ \theta_{ik}^{m*}=\int_{R_{\omega}}D\nabla\varphi_{i}^{c}\cdot\nabla\varphi_{k}^{m,c},\ \ \theta_{ik}^{*}=\int_{R_{\omega}}D\nabla\varphi_{i}^{c}\cdot\nabla\varphi_{k}^{c}.$$ (33)

We can also approximate the term with time derivative

$$\begin{split}&\int_{R_{\omega}}\phi{\partial c\over\partial t}v\approx\int_{R_{\omega}}\phi{\partial(\varphi_{i}^{c}C_{i}+\varphi_{i}^{m,c}\nabla_{m}C_{i})\over\partial t}(\varphi_{j}^{c}V_{j}+\varphi_{j}^{n,c}\nabla_{n}V_{j})\approx\int_{R_{\omega}}\phi{\partial(\varphi_{i}^{c}C_{i})\over\partial t}\varphi_{j}^{c}V_{j}\approx\\ &{\partial C_{i}(x_{\omega})\over\partial t}V_{j}(x_{\omega})\int_{R_{\omega}}\phi\varphi_{i}^{c}\varphi_{j}^{c}={\partial C_{i}(x_{\omega})\over\partial t}V_{j}(x_{\omega})\gamma_{ij},\end{split}$$ (34)

where

$$\gamma_{ij}=\int_{R_{\omega}}\phi\varphi_{i}^{c}\varphi_{j}^{c}.$$

Here for the second approximate equality, we use the fact that $\varphi_{i}^{m,c}$ is of the order $\epsilon$ and $\varphi_{i}^{c}$ is of the order $1$ as will be stated in (39). Again we make use of continuous approximations for $C_{i}$ and $V_{j}$, allowing us to derive the following approximation

$$\int_{R_{\omega}}\phi{\partial c\over\partial t}v\approx{\partial c_{i}\over\partial t}\gamma_{ij}v_{j}.$$ (35)

Next, we need to perform homogenization for the convection term. With our previous assumptions and estimates for the solutions of cell problems, we write

$$\begin{split}&\int_{R_{\omega}}(u\cdot\nabla c)v\approx\int_{R_{\omega}}\left[-\kappa\nabla(\varphi_{s}^{p}P_{s}+\varphi_{s}^{l,p}\nabla_{l}P_{s})\cdot\nabla(\varphi_{i}^{c}C_{i}+\varphi_{i}^{m,c}\nabla_{m}C_{i})\right]\left[\varphi_{j}^{c}V_{j}+\varphi_{j}^{n,c}\nabla_{n}V_{j}\right]\approx\\ &\int_{R_{\omega}}\left[-\kappa\nabla(\varphi_{s}^{p}P_{s}+\varphi_{s}^{l,p}\nabla_{l}P_{s})\cdot\nabla(\varphi_{i}^{c}C_{i}+\varphi_{i}^{m,c}\nabla_{m}C_{i})\right]\varphi_{j}^{c}V_{j}=\\ &\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{p}P_{s}\cdot\nabla\varphi_{i}^{c}C_{i})\varphi_{j}^{c}V_{j}+\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{p}P_{s}\cdot\nabla\varphi_{i}^{m,c}\nabla_{m}C_{i})\varphi_{j}^{c}V_{j}+\\ &\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{l,p}\nabla_{l}P_{s}\cdot\nabla\varphi_{i}^{c}c_{i})\varphi_{j}^{c}V_{j}+\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{l,p}\nabla_{l}P_{s}\cdot\nabla\varphi_{i}^{m,c}\nabla_{m}C_{i})\varphi_{j}^{c}V_{j}\approx\\ &P_{s}(x_{\omega})C_{i}(x_{\omega})V_{j}(x_{\omega})\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{p}\cdot\nabla\varphi_{i}^{c})\varphi_{j}^{c}+P_{s}(x_{\omega})\nabla_{m}C_{i}(x_{\omega})V_{j}(x_{\omega})\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{p}\cdot\nabla\varphi_{i}^{m,c})\varphi_{j}^{c}+\\ &\nabla_{l}P_{s}(x_{\omega})C_{i}(x_{\omega})V_{j}(x_{\omega})\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{l,p}\cdot\nabla\varphi_{i}^{c})\varphi_{j}^{c}+\nabla_{l}P_{s}(x_{\omega})\nabla_{m}C_{i}(x_{\omega})V_{j}(x_{\omega})\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{l,p}\cdot\nabla\varphi_{i}^{m,c})\varphi_{j}^{c}=\\ &P_{s}(x_{\omega})C_{i}(x_{\omega})V_{j}(x_{\omega})\zeta_{sij}+P_{s}(x_{\omega})\nabla_{m}C_{i}(x_{\omega})V_{j}(x_{\omega})\chi_{sij}^{m}+\\ &\nabla_{l}P_{s}(x_{\omega})C_{i}(x_{\omega})V_{j}(x_{\omega})\Upsilon_{sij}^{l}+\nabla_{l}P_{s}(x_{\omega})\nabla_{m}C_{i}(x_{\omega})V_{j}(x_{\omega})\iota_{sij}^{lm},\end{split}$$ (36)

where

$$\begin{split}&\zeta_{sij}=\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{p}\cdot\nabla\varphi_{i}^{c})\varphi_{j}^{c},\quad\chi_{sij}^{m}=\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{p}\cdot\nabla\varphi_{i}^{m,c})\varphi_{j}^{c},\\ &\Upsilon_{sij}^{l}=\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{l,p}\cdot\nabla\varphi_{i}^{c})\varphi_{j}^{c},\quad\iota_{sij}^{lm}=\int_{R_{\omega}}(-\kappa)(\nabla\varphi_{s}^{l,p}\cdot\nabla\varphi_{i}^{m,c})\varphi_{j}^{c}.\end{split}$$ (37)

Applying the continuous approximations for $C_{i}$ and $V_{j}$ gives us

$$\begin{split}&\int_{R_{\omega}}(u\cdot\nabla c)v\approx P_{s}C_{i}V_{j}\zeta_{sij}+P_{s}\nabla_{m}C_{i}V_{j}\chi_{sij}^{m}+\nabla_{l}P_{s}C_{i}V_{j}\Upsilon_{sij}^{l}+\nabla_{l}P_{s}\nabla_{m}C_{i}V_{j}\iota_{sij}^{lm}.\end{split}$$ (38)

We can establish the following estimates

$$\begin{split}&\|\varphi_{i}^{c}\|=O(1),\ \|\nabla\varphi_{i}^{c}\|=O({1\over\epsilon}),\\ &\|\varphi_{i}^{m,c}\|=O(\epsilon),\ \|\nabla\varphi_{i}^{m,c}\|=O(1).\\ \end{split}$$ (39)

Subsequently, we derive the following scalings

$$\begin{split}&\theta_{ij}^{*}=O(\cfrac{|R_{\omega}|}{\epsilon^{2}}),\ \theta_{ij}^{m*}=O(\cfrac{|R_{\omega}|}{\epsilon}),\ \eta_{ij}^{mn}=O(|R_{\omega}|),\ \gamma_{ij}=O(|R_{\omega}|),\\ &\zeta_{sij}=O(\cfrac{|R_{\omega}|}{\epsilon^{2}}),\ \chi_{sij}^{m}=O(\cfrac{|R_{\omega}|}{\epsilon}),\ \Upsilon_{sij}^{l}=O(\cfrac{|R_{\omega}|}{\epsilon}),\ \iota_{sij}^{lm}=O(|R_{\omega}|).\end{split}$$

We introduce the scaled effective properties

$$\begin{split}&\widehat{\theta_{ij}}=\cfrac{\epsilon^{2}}{|R_{\omega}|}\beta_{ij}^{*},\ \widehat{\theta_{ij}^{m*}}=\cfrac{\epsilon}{|R_{\omega}|}\beta_{ij}^{m*},\ \widehat{\eta_{ij}^{mn}}=\frac{1}{|R_{\omega}|}\alpha_{ij}^{mn},\ \widehat{\gamma_{ij}}=\frac{1}{|R_{\omega}|}\gamma_{ij},\\ &\widehat{\zeta_{sij}}=\cfrac{\epsilon^{2}}{|R_{\omega}|}\zeta_{sij},\ \widehat{\chi_{sij}^{m}}=\cfrac{\epsilon}{|R_{\omega}|}\chi_{sij}^{m},\ \widehat{\Upsilon_{sij}^{l}}=\cfrac{\epsilon}{|R_{\omega}|}\Upsilon_{sij}^{l},\ \widehat{\iota_{sij}^{lm}}=\frac{1}{|R_{\omega}|}\iota_{sij}^{lm}.\end{split}$$ (40)

With these scalings, we have

$$\begin{split}&\int_{\Omega}\phi{\partial c\over\partial t}v+\int_{\Omega}(u\cdot\nabla c)v+\int_{\Omega}D\nabla c\cdot\nabla v\approx\\ &\int_{\Omega}\widehat{\gamma_{ij}}{\partial C_{i}\over\partial t}V_{j}+\int_{\Omega}\widehat{\eta_{ij}^{mn}}\nabla_{m}C_{i}\nabla_{n}V_{j}+{1\over\epsilon}\int_{\Omega}\widehat{\theta_{ij}^{m}}\nabla_{m}C_{i}V_{j}+{1\over\epsilon}\int_{\Omega}\widehat{\theta_{ij}^{m}}C_{i}\nabla_{m}V_{j}+{1\over\epsilon^{2}}\int_{\Omega}\widehat{\theta_{ij}}C_{i}V_{j}+\\ &{1\over\epsilon^{2}}\int_{\Omega}P_{s}C_{i}V_{j}\widehat{\zeta_{sij}}+{1\over\epsilon}\int_{\Omega}P_{s}\nabla_{m}C_{i}V_{j}\widehat{\chi_{sij}^{m}}+{1\over\epsilon}\int_{\Omega}\nabla_{l}P_{s}C_{i}V_{j}\widehat{\Upsilon_{sij}^{l}}+\int_{\Omega}\nabla_{l}P_{s}\nabla_{m}C_{i}V_{j}\widehat{\iota_{sij}^{lm}}.\end{split}$$ (41)

One can see that the sum of the third and fourth terms in (41) is negligible. It can be shown by using integration by parts. Consequently, we have

$$\begin{split}&\int_{\Omega}\phi{\partial c\over\partial t}v+\int_{\Omega}(u\cdot\nabla c)v+\int_{\Omega}D\nabla c\cdot\nabla v\approx\\ &\int_{\Omega}\widehat{\gamma_{ij}}{\partial C_{i}\over\partial t}V_{j}+\int_{\Omega}\widehat{\eta_{ij}^{mn}}\nabla_{m}C_{i}\nabla_{n}V_{j}+{1\over\epsilon^{2}}\int_{\Omega}(\widehat{\theta_{ij}}+P_{s}\widehat{\zeta_{sij}}+\epsilon\nabla_{l}P_{s}\widehat{\Upsilon_{sij}^{l}})C_{i}V_{j}+\\ &{1\over\epsilon}\int_{\Omega}(P_{s}\widehat{\chi_{sij}^{m}}+\epsilon\nabla_{l}P_{s}\widehat{\iota_{sij}^{lm}})\nabla_{m}C_{i}V_{j}=\\ &\int_{\Omega}\widehat{\gamma_{ij}}{\partial C_{i}\over\partial t}V_{j}+\int_{\Omega}\widehat{\eta_{ij}^{mn}}\nabla_{m}C_{i}\nabla_{n}V_{j}+{1\over\epsilon}\int_{\Omega}\widehat{\xi_{ij}^{m}}\nabla_{m}C_{i}V_{j}+{1\over\epsilon^{2}}\int_{\Omega}\widehat{\Theta_{ij}}C_{i}V_{j},\end{split}$$ (42)

where

$$\widehat{\xi_{ij}^{m}}=P_{s}\widehat{\chi_{sij}^{m}}+\epsilon\nabla_{l}P_{s}\widehat{\iota_{sij}^{lm}},\quad\widehat{\Theta_{ij}}=\widehat{\theta_{ij}}+P_{s}\widehat{\zeta_{sij}}+\epsilon\nabla_{l}P_{s}\widehat{\Upsilon_{sij}^{l}}.$$

Therefore, we obtain the following multicontinuum equations

$$\begin{split}\widehat{\gamma_{ij}}(x,P){\partial C_{j}\over\partial t}-\nabla_{n}(\widehat{\eta_{ij}^{mn}}(x,P)\nabla_{m}C_{j})+{1\over\epsilon}\widehat{\xi_{ij}^{m}}(x,P)\nabla_{m}C_{j}+{1\over\epsilon^{2}}\widehat{\Theta_{ij}}(x,P)C_{j}=h_{i}(x,C,P).\end{split}$$ (43)

The effective properties depend on $x$ because they can be different in each coarse block. Moreover, the coefficients of the time term and the source term depend on $P=(P_{1},P_{2},...,P_{N})$ because the cell problems’ solutions $\phi_{i}^{c}$ and $\phi_{i}^{m,c}$ depend on $P$.

3.4 Coupled multicontinuum model

Finally, we have the following coupled multicontinuum model described by macroscopic elliptic equations and convection-diffusion-reaction equations

$$\begin{split}&-\nabla_{n}(\widehat{\alpha_{ij}^{mn}}\nabla_{m}P_{j})+{1\over\epsilon^{2}}\widehat{\beta_{ij}}P_{j}=g_{i},\\ &\widehat{\gamma_{ij}}(x,P){\partial C_{j}\over\partial t}-\nabla_{n}(\widehat{\eta_{ij}^{mn}}(x,P)\nabla_{m}C_{j})+{1\over\epsilon}\widehat{\xi_{ij}^{m}}(x,P)\nabla_{m}C_{j}+{1\over\epsilon^{2}}\widehat{\Theta_{ij}}(x,P)C_{j}=h_{i}(x,C,P).\end{split}$$ (44)

We can obtain the macroscopic solutions by solving (44) on the coarse grid. In the next section, we present numerical results to verify the obtained model.

4 Numerical examples

We will apply our proposed multicontinuum homogenization approach to solve the system of equations  (1) subject to certain boundary conditions, where the spatial domain $\Omega$ is chosen to be a unit square $[0,1]\times[0,1]$. We assume that the permeability coefficient $\kappa$ and the diffusion coefficient $D$ are of high contrast, and we will consider different coefficient fields in the following examples.

We will divide the computational domain $\Omega$ into $M\times M$ square coarse blocks of the same size and define the coarse mesh size to be $H=1/M$. Each whole coarse block will be taken as an RVE for itself. The oversampling RVE is defined as an extension of itself by $l$ layers of coarse blocks. We denote the regions of low and high values by $\Omega_{1}$ and $\Omega_{2}$, respectively, and define the relative $L^{2}$ errors of solution in $\Omega_{1}$ and $\Omega_{2}$ at a specific time $t$ by

$$e^{(i)}(t)=\sqrt{\frac{\sum_{K}|\frac{1}{|K|}\int_{K}C_{i}(x,t)dx-\frac{1}{|K\cap\Omega_{i}|}\int_{K\cap\Omega_{i}}c(x,t)dx|^{2}}{\sum_{K}|\frac{1}{|K\cap\Omega_{i}|}\int_{K\cap\Omega_{i}}c(x,t)dx|^{2}}},$$ (45)

where $i=1,2$, and $K$ denotes the RVE, which is taken to be the coarse block.

4.1 Example 1: Layered field

4.1.1 Case 1

First, we consider the case when the permeability field $\kappa$ and the diffusion field $D$ consist of horizontal layers, as shown in Figure 2, and

$$\kappa(x)=D(x)=\begin{cases}10^{-4},\quad&x\in\Omega_{1},\\ 1,\quad&x\in\Omega_{2}.\end{cases}$$ (46)

We assume $\phi=1$, $g(x)=e^{-40\left((x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}\right)}$ and $h(x)=0.1e^{-40\left((x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}\right)}$ for any $x=(x_{1},x_{2})\in\Omega$. We set the following homogeneous Dirichlet boundary conditions

$$\begin{split}p(x,t)=0,\quad&x\in\partial\Omega,\quad\forall t\in[0,\infty),\\ c(x,t)=0,\quad&x\in\partial\Omega,\quad\forall t\in[0,\infty).\end{split}$$ (47)

Figure 2: Layered field for Example 1 ($\Omega_{1}$: blue regions; $\Omega_{2}$: yellow regions)

To discretize the spatial domain, we fix the fine mesh size to be $1/400$ and take the coarse mesh size to be $1/20$ and $1/40$, respectively. We use an implicit Euler scheme for the time discretization and choose the time step $\tau=0.001$. The number of oversampling layers $l$ is taken to be $\lceil-2\log(H)\rceil$, that is, $l=6$ for $H=1/20$ and $l=8$ for $H=1/40$, which is inspired by the analysis in [20, 19].

Figure 3: Initial concentration when $t=0$

We depict the concentration at the initial time $t=0$ in Figure 3. The relative $L^{2}$ errors of our proposed algorithm are shown in Table 1. We also depict the results when $H=1/40$ in Figure 4 and Figure 5. Figure 4 demonstrates the multiscale solution $C_{1}$ and the reference averaged solution in $\Omega_{1}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$; Figure 5 demonstrates the multiscale solution $C_{2}$ and its corresponding reference averaged solution.

In this paragraph, we examine the obtained results. First, we can clearly see that the errors at different time $t$ between the multiscale and reference solutions are very small and our method can approximate the results accurately. Furthermore, decreasing the coarse mesh size $H$ makes the results more precise. Also, one can observe that the averaged concentration in $\Omega_{1}$ transports more slowly than that in $\Omega_{2}$. This is reasonable because the diffusion coefficient is much higher in $\Omega_{2}$ than in $\Omega_{1}$, as defined in (46).

Figure 4: Solutions of concentration when $H=1/40$ for Case 1 in Example 1. First column: multiscale solution $C_{1}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$. Second column: reference averaged solution in $\Omega_{1}$ at the corresponding time instants.

Figure 5: Solutions of concentration when $H=1/40$ for Case 1 in Example 1. First column: multiscale solution $C_{2}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$. Second column: reference averaged solution in $\Omega_{2}$ at the corresponding time instants.

Table 1: Relative $L^{2}$ errors at $t=0.02$, $0.1$, $0.5$, $1$, $2$ for Case 1 in Example 1. Left: $H=1/20$ and $l=6$. Right: $H=1/40$ and $l=8$.

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$4.63\%$	$2.13\%$
$0.1$	$3.03\%$	$2.89\%$
$0.5$	$2.54\%$	$1.90\%$
$1.0$	$2.39\%$	$1.53\%$
$2.0$	$1.77\%$	$1.20\%$

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$2.43\%$	$2.20\%$
$0.1$	$0.76\%$	$1.32\%$
$0.5$	$0.39\%$	$0.22\%$
$1.0$	$0.55\%$	$0.51\%$
$2.0$	$0.60\%$	$0.82\%$

4.1.2 Case 2

For Case 2, we change the boundary condition for the pressure equation to be an inhomogeneous Dirichlet boundary condition while keeping the homogeneous Dirichlet condition for the concentration equation, that is,

$$\begin{split}p(x,t)=x,\quad&x\in\partial\Omega,\quad\forall t\in[0,\infty),\\ c(x,t)=0,\quad&x\in\partial\Omega,\quad\forall t\in[0,\infty).\end{split}$$ (48)

All the other assumptions are the same as in Case 1.

We list the relative $L^{2}$ errors when $H=1/20$ and $H=1/40$ in Table 2. The numerical solutions when $H=1/40$ are presented in Figures 6, 7, and 8. One can notice that our method yields results that closely align with the reference solutions in this case. We can also observe from Figure 6 that the concentration exhibits faster transport in the $-x$-direction compared to the $x$-direction. Indeed, due to the imposed boundary condition (48), the pressure generally increases along the $x$-direction on average. According to Darcy’s law, the fluid predominantly flows toward regions of decreasing pressure, thereby carrying substances in the $-x$-direction and offsetting the effects of diffusion.

Figure 6: Snapshots of $C_{1}$ and $C_{2}$ when $t=1$ for Case 2 in Example

Figure 7: Solutions of concentration when $H=1/40$ for Case 2 in Example 1. First column: multiscale solution $C_{1}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$. Second column: reference averaged solution in $\Omega_{1}$ at the corresponding time instants.

Figure 8: Solutions of concentration when $H=1/40$ for Case 2 in Example 1. First column: multiscale solution $C_{2}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$. Second column: reference averaged solution in $\Omega_{2}$ at the corresponding time instants.

Table 2: Relative $L^{2}$ errors at $t=0.02$, $0.1$, $0.5$, $1$, $2$ for Case 2 in Example 1. Left: $H=1/20$ and $l=6$. Right: $H=1/40$ and $l=8$.

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$4.64\%$	$2.17\%$
$0.1$	$3.03\%$	$2.89\%$
$0.5$	$2.54\%$	$1.87\%$
$1.0$	$2.39\%$	$1.50\%$
$2.0$	$1.75\%$	$1.16\%$

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$2.44\%$	$2.25\%$
$0.1$	$0.76\%$	$1.34\%$
$0.5$	$0.39\%$	$0.23\%$
$1.0$	$0.55\%$	$0.52\%$
$2.0$	$0.60\%$	$0.82\%$

4.1.3 Case 3

For Case 3, we consider a homogeneous Neumann condition for the concentration equation and the boundary conditions will be

$$\begin{split}p(x,t)=x,\quad&x\in\partial\Omega,\quad\forall t\in[0,\infty),\\ \nabla c(x,t)\cdot\nu=0,\quad&x\in\partial\Omega,\quad\forall t\in[0,\infty),\end{split}$$ (49)

where $\nu$ is the outward normal vector defined on $\partial\Omega$. We keep all the other assumptions the same as in Case 1.

The relative $L^{2}$ errors when $H=1/20$ and $H=1/40$ are presented in Table 3. The numerical results when $H=1/40$ are shown in Figures 9 and 10. Our method ensures negligible discrepancy from the reference solution, and the convection can still be observed.

Figure 9: Solutions of concentration when $H=1/40$ for Case 3 in Example 1. First column: multiscale solution $C_{1}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$. Second column: reference averaged solution in $\Omega_{1}$ at the corresponding time instants.

Figure 10: Solutions of concentration when $H=1/40$ for Case 3 in Example 1. First column: multiscale solution $C_{2}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$. Second column: reference averaged solution in $\Omega_{2}$ at the corresponding time instants.

Table 3: Relative $L^{2}$ errors at $t=0.02$, $0.1$, $0.5$, $1$, $2$ for Case 3 in Example 1. Left: $H=1/20$ and $l=6$. Right: $H=1/40$ and $l=8$.

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$4.64\%$	$2.90\%$
$0.1$	$3.04\%$	$3.23\%$
$0.5$	$2.40\%$	$1.99\%$
$1.0$	$2.11\%$	$1.92\%$
$2.0$	$1.77\%$	$1.72\%$

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$2.44\%$	$3.00\%$
$0.1$	$0.83\%$	$1.82\%$
$0.5$	$0.46\%$	$0.24\%$
$1.0$	$0.41\%$	$0.20\%$
$2.0$	$0.24\%$	$0.20\%$

4.2 Example 2: Circular field

In this example, we consider a different structure of permeability field $\kappa$ and the diffusion field $D$ as illustrated in Figure 11 and

$$\kappa(x)=D(x)=\begin{cases}10^{-4},\quad&x\in\Omega_{1},\\ 1,\quad&x\in\Omega_{2}.\end{cases}$$ (50)

All other assumptions remain unchanged.

Again, we apply our algorithm to different boundary conditions corresponding to Cases 1, 2, and 3 in Example 1 and present the errors in Tables 4, 5, and 6, respectively. Numerical solutions for Case 3 when $H=1/40$ are depicted in Figures 12 and 13, where the color bars have been adjusted to display the results more clearly. One can observe that the transport of the concentration for Case 3 in Example 2 is isotropic and more rapid than in Example 1, which is a consequence of the symmetry and porosity of the circular field. We also note that the numerical results for all three cases are in tight agreement with reference solutions.

Figure 11: Circular field for Example 2 ($\Omega_{1}$: blue regions; $\Omega_{2}$: yellow regions)

Table 4: Relative $L^{2}$ errors at $t=0.02$, $0.1$, $0.5$, $1$, $2$ for Case 1 in Example 2. Left: $H=1/20$ and $l=6$. Right: $H=1/40$ and $l=8$.

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$4.04\%$	$1.90\%$
$0.1$	$2.07\%$	$1.75\%$
$0.5$	$1.21\%$	$0.91\%$
$1.0$	$0.87\%$	$0.79\%$
$2.0$	$0.60\%$	$0.45\%$

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$1.74\%$	$0.44\%$
$0.1$	$0.19\%$	$0.73\%$
$0.5$	$0.41\%$	$0.21\%$
$1.0$	$0.28\%$	$0.22\%$
$2.0$	$0.20\%$	$0.13\%$

Table 5: Relative $L^{2}$ errors at $t=0.02$, $0.1$, $0.5$, $1$, $2$ for Case 2 in Example 2. Left: $H=1/20$ and $l=6$. Right: $H=1/40$ and $l=8$.

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$4.04\%$	$1.91\%$
$0.1$	$2.07\%$	$1.76\%$
$0.5$	$1.20\%$	$0.90\%$
$1.0$	$0.86\%$	$0.78\%$
$2.0$	$0.59\%$	$0.45\%$

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$1.74\%$	$0.45\%$
$0.1$	$0.20\%$	$0.73\%$
$0.5$	$0.41\%$	$0.21\%$
$1.0$	$0.28\%$	$0.22\%$
$2.0$	$0.21\%$	$0.13\%$

Table 6: Relative $L^{2}$ errors at $t=0.02$, $0.1$, $0.5$, $1$, $2$ for Case 3 in Example 2. Left: $H=1/20$ and $l=6$. Right: $H=1/40$ and $l=8$.

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$4.04\%$	$1.90\%$
$0.1$	$2.06\%$	$1.80\%$
$0.5$	$0.65\%$	$0.38\%$
$1.0$	$0.22\%$	$0.20\%$
$2.0$	$0.17\%$	$0.16\%$

$t$	$e^{(1)}(t)$	$e^{(2)}(t)$
$0.02$	$1.74\%$	$0.48\%$
$0.1$	$0.21\%$	$0.78\%$
$0.5$	$0.22\%$	$0.08\%$
$1.0$	$0.07\%$	$0.05\%$
$2.0$	$0.04\%$	$0.04\%$

Figure 12: Solutions of concentration when $H=1/40$ for Case 3 in Example 2. First column: multiscale solution $C_{1}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$. Second column: reference averaged solution in $\Omega_{1}$ at the corresponding time instants.

Figure 13: Solutions of concentration when $H=1/40$ for Case 3 in Example 2. First column: multiscale solution $C_{2}$ at $t=0.02$, $0.1$, $0.5$, $1$, $2$. Second column: reference averaged solution in $\Omega_{2}$ at the corresponding time instants.

5 Conclusion

In this work, we have applied the multicontinuum homogenization method to develop a new multicontinuum model for the coupled flow and transport equations. Coupled constraint cell problems were carefully formulated to obtain the localizable multiscale basis functions, which can capture the heterogeneity and high contrast properties. Macroscopic variables were introduced to represent the local averages of solutions in each continuum. By performing the multicontinuum expansions to both the flow and transport equations, we have achieved a macroscopic system consisting of homogenized elliptic equations and convection-diffusion-reaction equations.

We have conducted a number of numerical experiments to verify the obtained multicontinuum model. We have considered high-contrast layered and circular heterogeneous media with homogeneous and inhomogeneous Dirichlet boundary conditions and mixed Dirichlet-Neumann boundary conditions. Numerical results have shown that the proposed model allows us to approximate the reference averaged solutions for all considered cases accurately. The computed errors demonstrate the convergence of the proposed method on the coarse grid size.