Convergence of the CEM-GMsFEM for compressible flow in highly heterogeneous media

Leonardo A. Poveda Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong SAR, China Shubin Fu Eastern Institute for Advanced Study, Ningbo, China Eric T. Chung Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong SAR, China Lina Zhao Department of Mathematics, City University of Hong Kong, Hong Kong SAR, China

Abstract

This paper presents and analyses a Constraint Energy Minimization Generalized Multiscale Finite Element Method (CEM-GMsFEM) for solving single-phase non-linear compressible flows in highly heterogeneous media. The construction of CEM-GMsFEM hinges on two crucial steps: First, the auxiliary space is constructed by solving local spectral problems, where the basis functions corresponding to small eigenvalues are captured. Then the basis functions are obtained by solving local energy minimization problems over the oversampling domains using the auxiliary space. The basis functions have exponential decay outside the corresponding local oversampling regions. The convergence of the proposed method is provided, and we show that this convergence only depends on the coarse grid size and is independent of the heterogeneities. An online enrichment guided by a posteriori error estimator is developed to enhance computational efficiency. Several numerical experiments on a three-dimensional case to confirm the theoretical findings are presented, illustrating the performance of the method and giving efficient and accurate numerical

Keywords: Constraint energy minimization, multiscale finite element methods, compressible flow, highly heterogeneous, local spectral problems

1 Introduction

The numerical solution of non-linear partial differential equations defined on domains with multiscale and heterogeneous properties is an active research subject in the scientific community. The subject is related to several engineering applications, such as composite materials, porous media flow, and fluid mechanics. A common feature for all these applications is that they are very computationally challenging and often impossible to solve within an acceptable tolerance using standard fine-scale approximations due to the disparity between scales that need to be represented and the inherent nonlinearities. For this reason, coarse-grid computational models are often used. These approaches are usually referred to as multiscale methods in the literature, among which we may mention: Multiscale Finite Element Method [17], the Variational Multiscale Method [18], Mixed Multiscale Finite Element Method [6], Mixed Mortar Multiscale Finite Element Method [2], the two-scale Finite Element Method [24], and the Multiscale Finite Volume method [19, 20]. The aforementioned methods share model reduction techniques, using different structures to find multiscale solutions, especially in many practical applications such as fluid flow simulations, for instance, [1, 5, 25, 21, 16, 27, 15, 26]. In particular, we consider a family of an extended version of MsFEM, known as the Generalized Multiscale Finite Element Method (GMsFEM) that was first introduced by [12, 13, 8, 11]. The main idea of GMsFEM is to construct localized basis functions by solving local spectral problems that are used to approximate the solution on a coarse grid incorporating fine-scale features. Following this, we construct an auxiliary space associated with local spectral problems in the coarse grid. The first few eigenfunctions corresponding to small eigenvalues (the convergence depends on the decay of the spectral problems [12]) are considered as the multiscale basis functions. In this paper, we extend the Constraint Energy Minimization Generalized Multiscale Finite Element Method (CEM-GMsFEM) developed in [9, 22] to solve single-phase nonlinear compressible flow. The key ideas of the method can be summarized as follows: First, we construct the auxiliary basis functions by solving local spectral problems. Then, by using oversampling techniques and localization (cf. [23]), we solve an appropriate energy subject to some constrainted oversampling regions to find the required basis functions. Finally, the resulting basis functions are shown to have exponential decay away from the target coarse element, and therefore, they are localizable. Then we rigorously analyze the convergence error estimates for the proposed scheme. Our theories indicate that the convergence rate behaves as $H/\Lambda$ , where $H$ denotes the coarse-grid size and $\Lambda$ is proportional to the minimum (taken over all coarse regions) of the eigenvalues that the corresponding eigenvectors are not included in the coarse space. Since the problems under consideration are nonlinear, some novel methodologies shall be incorporated to overcome the difficulties present in the analysis. Several numerical experiments are carried out to demonstrate the capabilities and efficiency of the proposed method.

The outline of the paper is following. In Section 2, we briefly derive a mathematical model for compressible fluid flows in porous media. Section 3 is devoted to constructing the offline multiscale space and framework of CEM-GMsFEM. Section 4 presents our convergence analysis for the proposed method. Numerical experiments are presented in Section 5. Finally, concluding remarks and future perspectives are given in Section 6.

2 Mathematical model for compressible fluid flow

In this section, we consider the governing equations of the single-phase, non-linear compressible fluid flow processes in a porous medium that are defined by

\left\{\begin{split}\partial_{t}(\phi\rho)-\nabla\cdot\left(\frac{\kappa}{\mu}\rho\nabla p\right)&=q,\quad\mbox{in }\Omega_{T}:=\Omega\times(0,T],\quad T>0,\\ \frac{\kappa}{\mu}\rho\nabla p\cdot n&=0,\quad\mbox{on }\Gamma_{N}\times(0,T],\\ p&=p^{D},\quad\mbox{on }\Gamma_{D}\times(0,T],\\ p&=p_{0},\quad\mbox{on }\Omega\times\{t=0\}.\end{split}\right.

(2.1)

For simplicity of presentation, let $\Omega\subset\mathbb{R}^{d}$ be the computational domain with a boundary defined by $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$ . We will henceforth neglect the gravity effects and capillary forces and assume that $\phi$ , the porosity of the medium, is assumed to be a constant. We aim to seek the fluid pressure $p$ , $\kappa$ denotes the permeability field that may be highly heterogeneous, such that $\kappa_{0}\leq\kappa\leq\kappa_{1}$ , where $0<\kappa_{0}<\kappa_{1}<\infty$ ; and $\mu$ is the constant fluid viscosity. The fluid density $\rho$ is a function of the fluid pressure $p$ defined as

\rho(p)=\rho_{\mathrm{ref}}e^{c(p-p_{\mathrm{ref}})},

(2.2)

where $\rho_{\mathrm{ref}}$ is the given reference density and $p_{\mathrm{ref}}$ is the reference pressure. Finally, $n$ denotes the outward unit-normal vector on $\partial\Omega$ .

For a sub-domain $D\subset\Omega$ , let $\mathrm{V}:=\mathrm{H}^{1}(D)$ be the standard Sobolev spaces endowed with the norm $\|\cdot\|_{1,D}$ . We further denote by $(\cdot,\cdot)$ and $\|\cdot\|_{0,D}$ the inner product and norm, respectively in $\mathrm{L}^{2}(D)$ . The subscript $D$ will be omitted whenever $D=\Omega$ . In addition, we use the space $\mathrm{V}_{0}:=\mathrm{H}^{1}_{0}(D)$ , which is a subspace of $\mathrm{V}$ made of functions that vanish at $\Gamma_{D}$ . Finally, let $\mathrm{L}^{2}(0,T;\mathrm{L}^{2}(D))$ denote the set of functions with norm

\|v\|_{\mathrm{L}^{2}(0,T;\mathrm{L}^{2}(D))}=\left(\int_{0}^{T}\|v(\cdot,t)\|^{2}_{0,D}dt\right)^{1/2}.

Throughout the paper, $a\preceq b$ means there exists a positive constant $C$ independent of the mesh size such that $a\leq Cb$ .

2.1 A finite element approximation

In this subsection, we introduce the notions of fine and coarse grids to discretize the problem (2.1). Let ${\cal T}^{H}$ be a usual conforming partition of the computational domain $\Omega$ into coarse block $K\in{\cal T}^{H}$ with diameter $H$ . Then, we denote this partition as the coarse grid and assume that each coarse element is partitioned into a connected union of fine-grid blocks. In this case, the fine grid partition will be denoted by ${\cal T}^{h}$ , and is, by definition, a refinement of the coarse grid ${\cal T}^{H}$ , such that $h\ll H$ . We shall denote $\{x_{i}\}_{i=1}^{N_{\tiny{\rm{coarse}}}}$ as the vertices of the coarse grid ${\cal T}^{H}$ , where $N_{\tiny{\rm{coarse}}}$ denotes the number of coarse nodes. We define the neighborhood of the node $x_{i}$ by

\omega_{i}=\bigcup\{K_{j}\in{\cal T}^{H}:x_{i}\in\overline{K}_{j}\}.

In addition, for CEM-GMsFEM considered in this paper, we have that given a coarse block $K_{i}$ , we represent the oversampling region $K_{i,m}\subset\Omega$ obtained by enlarging $K_{i}$ with $m\geq 1$ coarse grid layers, see Fig. 1.

Figure 1: Illustration of the

2

D multiscale grid with a typical coarse element

K_{i}

and oversampling domain

K_{i,2}

, the fine grid element and neighborhood

\omega_{i}

of the node

x_{i}

We consider the linear finite element space $\mathrm{V}^{h}$ associated with the grid ${\cal T}^{h}$ , where the basis functions in this space are the standard Lagrange basis functions defined as $\{\eta^{i}\}_{i=1}^{N_{\tiny{\rm{fine}}}}$ . Then, the semi-discrete finite element approximation to (2.1) on the fine grid is to find $p^{h}\in\mathrm{V}^{h}$ such that

(\phi\partial_{t}\rho(p^{h}),v)+\left(\tfrac{\kappa}{\mu}\rho(p^{h})\nabla p^{h},\nabla v\right)=(q,v),\quad\mbox{for each }v\in\mathrm{V}^{h}.

(2.3)

We can now define the fully-discrete scheme for the discrete formulation (2.3). Let $0=t_{0}<t_{1}<\cdots<t_{N_{\tiny{\rm{time}}}-1}<t_{N_{\tiny{\rm{time}}}}=T$ be a partition of the interval $[0,T]$ , with time-step size given by $\Delta_{n}=t_{n}-t_{n-1}$ , for $n=1,\dots,N_{\tiny{\rm{time}}}$ , where $N_{\tiny{\rm{time}}}$ is an integer. The backward Euler time integration leads to finding $p^{h}_{n}$ such that

(\phi\rho(p^{h}_{n}),v)-(\phi\rho(p^{h}_{n-1}),v)+\Delta_{n}\left(\tfrac{\kappa}{\mu}\rho(p^{h}_{n})\nabla p^{h}_{n},\nabla v\right)=\Delta_{n}(q,v),\quad\mbox{for each }v\in\mathrm{V}^{h}.

(2.4)

Linearization of (2.4) via Newton-Raphson iteration yields a iterative linear matrix problem,

{\mathbf{J}}^{n,k}\boldsymbol{\delta}_{p^{n,k}}=-{\mathbf{F}}^{n,k},

where ${\mathbf{J}}^{n,k}:=[J_{ji}^{n,k}]_{i,j=1}^{N_{\tiny{\rm{fine}}}}$ denotes the Jacobi matrix, with entries

J_{ji}^{n,k}:=(\phi\rho(p^{h,k}_{n})\eta_{i},\eta_{j})+\Delta_{n}\left(\frac{\kappa}{\mu}\rho(p^{h,k}_{n})\nabla\eta_{i},\nabla\eta_{j}\right)+\Delta_{n}\left(c\frac{\kappa}{\mu}\eta_{i}\rho(p^{h,k}_{n})\sum_{i}p^{i,k}_{n}\nabla\eta_{i},\nabla\eta_{j}\right),

${\mathbf{F}}^{n,k}:=[F_{j}^{n,k}]_{j=1}^{N_{\tiny{\rm{fine}}}}$ is the residual with entries

	$\displaystyle F_{j}^{n,k}$	$\displaystyle=\left(\phi\rho\left(\sum_{i=1}^{N_{\tiny{\rm{fine}}}}p^{i,k}_{n}\eta_{i}\right),\eta_{j}\right)-\left(\phi\rho\left(\sum_{i=1}^{N_{\tiny{\rm{fine}}}}p^{i}_{n-1}\eta_{i}\right),\eta_{j}\right)+\Delta_{n}\left(\frac{\kappa}{\mu}\rho\left(\sum_{i=1}^{N_{\tiny{\rm{fine}}}}p^{i,k}_{n}\eta_{i}\right)\sum_{i=1}^{N_{\tiny{\rm{fine}}}}p^{i,k}_{n}\nabla\eta_{i},\nabla\eta_{j}\right)$
		$\displaystyle\quad-\Delta_{n}(q,\eta_{j})=0,$

and $p^{k+1}_{n}=p^{k}_{n}+\boldsymbol{\delta}_{p^{n,k}}$ , where $p^{h,k}_{n}=\sum_{i}p^{i,k}_{n}\eta_{i}$ and $p^{h}_{n-1}=\sum_{i}p^{i}_{n-1}\eta_{i}$ , with $k$ and $k-1$ the new and old Newton iteration. Here $\{\eta^{i}\}_{i=1}^{N_{\tiny{\rm{fine}}}}$ represents the finite element basis functions for $\mathrm{V}^{h}$ .

3 Construction of CEM-GMsFEM basis functions

This section will describe the construction of CEM-GMsFEM basis functions using the framework of [9] and [22]. This procedure can be divided into two stages. The first stage involves constructing the auxiliary spaces by solving a local spectral problem in each coarse element $K$ , see [12]. The second stage is to provide the multiscale basis functions by solving some local constraint energy minimization problems in oversampling regions.

3.1 Auxiliary basis function

In this subsection, we present the construction of the auxiliary multiscale basis functions by solving the local eigenvalue problem for each coarse element $K_{i}$ . We consider $\mathrm{V}(K_{i}):=\mathrm{V}\scriptstyle\big{|}_{K_{i}}$ the restriction of the space $\mathrm{V}$ to the coarse element $K_{i}$ . We solve the following local eigenvalue problem: find $\{\lambda^{(i)}_{j},\varphi^{(i)}_{j}\}$ such that

a_{i}(\varphi^{(i)}_{j},w)=\lambda^{(i)}_{j}s_{i}(\varphi^{(i)}_{j},w),\quad\mbox{for each }w\in\mathrm{V}(K_{i}),

(3.1)

where

a_{i}(v,w):=\int_{K_{i}}\kappa\rho(p_{0})\nabla v\cdot\nabla wd\mathrm{x},\quad s_{i}(v,w):=\int_{K_{i}}\widetilde{\kappa}vwd\mathrm{x}.

Here $\widetilde{\kappa}=\rho(p_{0})\kappa\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}|\nabla\chi_{i}|^{2}$ , where $N_{\tiny{\rm{coarse}}}$ is the total number of neighborhoods, $p_{0}$ is the pressure $p$ at the initial time and $\{\chi_{i}\}$ is a set of partition of unity functions for the coarse grid ${\cal T}^{H}$ , see [3]. The problem defined above is solved on the fine grid in the actual computation. We assume that the eigenfunctions satisfy the normalized condition $s_{i}(\varphi^{(i)}_{j},\varphi^{(i)}_{j})=1$ . We let $\lambda^{(i)}_{j}$ be the eigenvalues of (3.1) arranged in ascending order. We shall use the first $L_{i}$ eigenfunctions to construct the local auxiliary multiscale space $\mathrm{V}_{\mathrm{aux}}^{(i)}:=\{\varphi^{(i)}_{j}:1\leq j\leq L_{i}\}$ . We can define the global auxiliary multiscale space as $\mathrm{V}_{\mathrm{aux}}:=\bigoplus_{i=1}^{N_{\tiny{\rm{coarse}}}}\mathrm{V}_{\mathrm{aux}}^{(i)}$ .

For the local auxiliary space $\mathrm{V}_{\mathrm{aux}}^{(i)}$ , the bilinear form $s_{i}$ given above defines an inner product with norm $\|v\|_{s(K_{i})}=s_{i}(v,v)^{1/2}$ . Then, we can define the inner product and norm for the global auxiliary multiscale space $\mathrm{V}_{\mathrm{aux}}$ , which are defined by

s(v,w)=\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}s_{i}(v,w),\quad\|v\|_{s}:=s(v,v)^{1/2},\quad\mbox{for each }v,w\in\mathrm{V}_{\mathrm{aux}}.

To construct the CEM-GMsFEM basis functions, we use the following definition.

Definition 3.1 ([9]).

Given a function $\varphi^{(i)}_{j}\in\mathrm{V}_{\mathrm{aux}}$ , if a function $\psi\in\mathrm{V}$ satisfies

s(\psi,\varphi^{(i)}_{j}):=1,\quad s(\psi,\varphi_{j^{\prime}}^{(i^{\prime})})=0,\quad\mbox{if }j^{\prime}\neq j\mbox{ or }i^{\prime}\neq i,

then, we say that is $\varphi^{(i)}_{j}$ -orthogonal where $s(v,w)=\sum_{i=1}^{N}s_{i}(v,w)$ .

Now, we define $\pi:\mathrm{V}\to\mathrm{V}_{\mathrm{aux}}$ to be the projection with respect to the inner product $s(v,w)$ . So, $\pi$ is defined by

\pi(v):=\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\pi_{i}(v)=\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\sum_{j=1}^{L_{i}}s_{i}(v,\varphi^{(i)}_{j})\varphi^{(i)}_{j},\quad\mbox{for each }v\in\mathrm{V},

where $\pi_{i}:L^{2}(K_{i})\to\mathrm{V}_{\mathrm{aux}}^{(i)}$ denotes the projection with respect to inner product $s_{i}(\cdot,\cdot)$ . The null space of the operator $\pi$ is defined by $\widetilde{\mathrm{V}}=\{v\in\mathrm{V}:\pi(v)=0\}$ . Now, we will construct the multiscale basis functions. Given a coarse block $K_{i}$ , we denote the oversampling region $K_{i,m}\subset\Omega$ obtained by enlarging $K_{i}$ with an arbitrary number of coarse grid layers $m\geq 1$ , see Figure 1. Let $\mathrm{V}_{0}(K_{i,m}):=\mathrm{H}^{1}_{0}(K_{i,m})$ . Then, we define the multiscale basis function

\psi_{j,\mathrm{ms}}^{(i)}=\mathrm{argmin}\{a(\psi,\psi):\psi\in\mathrm{V}_{0}(K_{i,m}),\,\psi\mbox{ is $\varphi^{(i)}_{j}$-orthogonal}\},

(3.2)

where $\mathrm{V}(K_{i,m})$ is the restriction of $\mathrm{V}_{0}$ in $K_{i,m}$ and $\mathrm{V}_{0}(K_{i,m})$ is the subspace of $\mathrm{V}(K_{i,m})$ with zero trace on $\partial K_{i,m}$ . The multiscale finite element space $\mathrm{V}^{\mathrm{ms}}$ is defined by

\mathrm{V}^{\mathrm{ms}}=\mathrm{span}\{\psi_{j,\mathrm{ms}}^{(i)}:1\leq j\leq L_{i},1\leq i\leq N_{\tiny{\rm{coarse}}}\}.

By introducing the Lagrange multiplier, the problem (3.2) is equivalent to the explicit form: find $\psi_{j,\mathrm{ms}}^{(i)}\in\mathrm{V}_{0}(K_{i,m})$ , $\lambda\in\mathrm{V}_{\mathrm{aux}}^{(i)}(K_{i})$ such that

\begin{cases}a(\psi_{j,\mathrm{ms}}^{(i)},\eta)+s(\eta,\lambda)&=0,\quad\mbox{for all }\eta\in\mathrm{V}(K_{i,m}),\\ s(\psi_{j,\mathrm{ms}}^{(i)}-\varphi^{(i)}_{j},\nu)&=0,\quad\mbox{for all }\nu\in\mathrm{V}_{\mathrm{aux}}^{(i)}(K_{i,m}),\end{cases}

where $\mathrm{V}_{\mathrm{aux}}^{(i)}(K_{i,m})$ is the union of all local auxiliary spaces for $K_{i}\subset K_{i,m}$ . Thus, the semi-discrete multiscale approximation reads as follows: find $p^{\mathrm{ms}}_{n}\in\mathrm{V}^{\mathrm{ms}}$ such that

(\phi\partial_{t}\rho(p^{\mathrm{ms}}_{n}),v)+\left(\tfrac{\kappa}{\mu}\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n},\nabla v\right)=(q,v),\quad\mbox{for each }v\in\mathrm{V}^{\mathrm{ms}}.

(3.3)

Using the backward Euler time-stepping scheme, we have a full-discrete formulation: find $p^{\mathrm{ms}}_{n}\in\mathrm{V}^{\mathrm{ms}}$ such that

(\phi\rho(p^{\mathrm{ms}}_{n}),v)-(\phi\rho(p^{\mathrm{ms}}_{n-1}),v)+\Delta_{n}\left(\tfrac{\kappa}{\mu}\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n},\nabla v\right)=\Delta_{n}(q,v),\quad\mbox{for each }v\in\mathrm{V}^{\mathrm{ms}}.

(3.4)

4 Convergence analysis

In this section, we establish the estimates of the convergence order of the proposed method.

4.1 Error estimates

In this subsection, we will present the convergence error estimates for the semi-discrete scheme (3.3). The analysis consists of two main steps. First, we derive the error estimate for the difference between the exact solution and its corresponding elliptic projection. Second, we estimate the difference between the solution of (2.1) and solution of (3.3) by the difference between exact solution and the elliptic projection solution of problem (2.1).

To begin, we let $\hat{p}\in\mathrm{V}^{\mathrm{ms}}$ be the elliptic projection of the function $p\in\mathrm{V}$ that is defined by

(\tfrac{\kappa}{\mu}\rho(p_{0})\nabla(p-\hat{p}),\nabla w)=0,\quad\mbox{for each }w\in\mathrm{V}^{\mathrm{ms}}.

(4.1)

The following lemma gives us the error bound of $\hat{p}$ for the nonlinear parabolic problem.

Lemma 4.1.

Let $p$ be the solution of (2.3). For each $t>0$ , we define the elliptic projection $\hat{p}\in\mathrm{V}^{\mathrm{ms}}$ that satisfies (4.1). Then,

\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})(t)\|_{0}\preceq\Lambda^{-1/2}H\|(\tfrac{\kappa}{\mu})^{1/2}(q-\partial_{t}\rho(p))(t)\|_{0},

where $\Lambda=\min_{1\leq i\leq N}\lambda^{(i)}_{L_{i}+1}$ .

Proof.

Let $\hat{p}\in\mathrm{V}^{\mathrm{ms}}$ be the projection of $p$ . By boundedness of $\rho$ and orthogonality property, we can write

	$\displaystyle\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\\|^{2}_{0}\preceq\int_{\Omega}(\tfrac{\kappa}{\mu})\rho(p_{0})\|\nabla(p-\hat{p})\|^{2}d\mathrm{x}$	$\displaystyle=(\tfrac{\kappa}{\mu}\rho(p_{0})\nabla(p-\hat{p}),\nabla(p-\hat{p}))$
		$\displaystyle=(\tfrac{\kappa}{\mu}\rho(p_{0})\nabla p,\nabla(p-\hat{p})).$

Invoking again the boundedness of $\rho$ , we have

|(\tfrac{\kappa}{\mu}\rho(p_{0})\nabla p,\nabla(p-\hat{p}))|\preceq|(\tfrac{\kappa}{\mu}\rho(p)\nabla p,\nabla(p-\hat{p}))|.

Now, from problem (2.3), we get that

(\tfrac{\kappa}{\mu}\rho(p)\nabla p,\nabla(p-\hat{p}))=(q-\phi\partial_{t}\rho(p),\nabla(p-\hat{p})),\quad\mbox{for all }t>0.

Therefore, we arrive at

\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\|^{2}_{0}\preceq(q-\phi\partial_{t}\rho(p),p-\hat{p})\leq\|\widetilde{\kappa}^{-1/2}(q-\phi\partial_{t}\rho(p))\|_{0}\|p-\hat{p}\|_{s},\quad\mbox{for all }t>0.

(4.2)

Since $p-\hat{p}\in\widetilde{\mathrm{V}}$ , implies that $\pi(p-\hat{p})=0$ . According to [9], the coarse blocks $K_{i}$ with $i=1,\dots,N_{\tiny{\rm{coarse}}}$ are disjoint, so we obtain that $\pi_{i}(p-\hat{p})=0$ , for all $i=1,2,\dots,N_{\tiny{\rm{coarse}}}$ . Thus, the local spectral problem (3.1) yields that

\|p-\hat{p}\|^{2}_{s}=\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\|p-\hat{p}\|^{2}_{s_{i}}=\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\|(I-\pi_{i})(p-\hat{p})\|^{2}_{s_{i}}\preceq\frac{1}{\Lambda}\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\|^{2}_{0,K_{i}},

(4.3)

where $\Lambda=\min_{1\leq i\leq N}\lambda^{(i)}_{L_{i}+1}$ . Therefore, by combining (4.2) and (4.3), and using the fact $|\nabla\chi_{i}|={\cal O}(H^{-1})$ , we obtain

\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\|_{0}\preceq\Lambda^{-1/2}H\|\kappa^{-1/2}(q-\phi\partial_{t}\rho(p))\|_{0}.

This completes the proof. ∎

The above estimate is the essence of the following result.

Lemma 4.2.

Under Assumptions of Lemma 4.1, the following estimates hold

	$\displaystyle\\|p-\hat{p}\\|_{0}$	$\displaystyle\preceq\Lambda^{-1}H^{2}\\|\kappa^{-1/2}(q-\phi\partial_{t}\rho(p))\\|_{0},$
	$\displaystyle\\|\partial_{t}(p-\hat{p})\\|_{0}$	$\displaystyle\preceq\Lambda^{-1}H^{2}\\|\kappa^{-1/2}\partial_{t}(q-\phi\partial_{t}\rho(p))\\|_{0}.$

Proof.

First, we will invoke the duality argument. For each $t>0$ , we define $w\in\mathrm{V}_{0}$ by

\int_{\Omega}\kappa\rho(p_{0})\nabla w\cdot\nabla vd\mathrm{x}=\int_{\Omega}(p-\hat{p})vd\mathrm{x},\quad\mbox{for each }v\in\mathrm{V}_{0},

and consider $\hat{w}$ as elliptic projection of $w$ in $\mathrm{V}^{\mathrm{ms}}$ . By Lemma 4.1, for $v=p-\hat{p}$ , we have

	$\displaystyle\\|p-\hat{p}\\|^{2}_{0}=\int_{\Omega}\kappa\rho(p_{0})\nabla w\cdot\nabla(p-\hat{p})d\mathrm{x}$	$\displaystyle=\int_{\Omega}\kappa\rho(p_{0})\nabla(w-\hat{w})\cdot\nabla(p-\hat{p})d\mathrm{x}$
		$\displaystyle\preceq\int_{\Omega}\tfrac{\kappa}{\mu}\rho(p_{0})\nabla(w-\hat{w})\cdot\nabla(p-\hat{p})d\mathrm{x}$
		$\displaystyle\leq\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(w-\hat{w})\\|_{0}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\\|_{0}$
		$\displaystyle\preceq\left(H\Lambda^{-1/2}\max\{\kappa^{-1/2}\}\\|p-\hat{p}\\|_{0}\right)$
		$\displaystyle\quad\times\left(H\Lambda^{-1/2}\\|\kappa^{-1/2}(q-\phi\partial_{t}\rho(p))\\|_{0}\right).$

Hence, we have

\|p-\hat{p}\|_{0}\preceq\Lambda^{-1}H^{2}\|\kappa^{-1/2}(q-\phi\partial_{t}\rho(p))\|_{0}.

By a similar computation, we can obtain the second estimate. This completes the proof. ∎

We will derive an error estimate for the difference between the solution of (2.1) and the CEM-GMsFEM solution of (3.3) using the framework of [15].

Theorem 4.3.

Let $p$ be the solution obtained from (2.3), $p^{\mathrm{ms}}\in\mathrm{V}^{\mathrm{ms}}$ be the multiscale solution of (3.3) using CEM-GMsFEM and $\hat{p}$ be an elliptic projection of $p$ in $\mathrm{V}^{\mathrm{ms}}$ . Then, the following error estimate holds

	$\displaystyle\\|(p-p^{\mathrm{ms}})(t)\\|_{0}^{2}+\int_{0}^{T}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-p^{\mathrm{ms}})\\|^{2}_{0}dt$	$\displaystyle\preceq\Lambda^{-1}H^{2}\left(\\|\kappa^{-1/2}\partial_{t}(q-\phi\partial_{t}\rho(p))(t)\\|_{0}^{2}\right.$
		$\displaystyle\quad+\left.\int_{0}^{T}\\|\kappa^{-1/2}\partial_{t}(q-\phi\partial_{t}\rho(p))\\|^{2}_{0}dt\right)$
		$\displaystyle\quad+\\|(\hat{p}-p^{\mathrm{ms}})(0)\\|^{2}_{0}.$

Proof.

Subtracting (3.3) from (2.3), and using (4.7), we have that

(\phi\partial_{t}\rho(p),v)+\left(\tfrac{\kappa}{\mu}\rho(p)\nabla p,\nabla v\right)-(\phi\partial_{t}\rho(p^{\mathrm{ms}}),v)-\left(\tfrac{\kappa}{\mu}\rho(p^{\mathrm{ms}})\nabla p^{\mathrm{ms}},\nabla v\right)=0,\quad\mbox{for each }v\in\mathrm{V}^{\mathrm{ms}}.

Since $\hat{p}\in\mathrm{V}^{\mathrm{ms}}$ , we put $v=\hat{p}-p^{\mathrm{ms}}$ , then follows that

\underset{I_{1}}{\underbrace{(\phi\partial_{t}(\rho(p)-\rho(p^{\mathrm{ms}})),\hat{p}-p^{\mathrm{ms}})}}+\underset{I_{2}}{\underbrace{\left(\tfrac{\kappa}{\mu}(\rho(p)\nabla p-\rho(p^{\mathrm{ms}})\nabla p^{\mathrm{ms}}),\nabla(\hat{p}-p^{\mathrm{ms}})\right)}}=0.

(4.4)

About $I_{1}$ , we can rewrite

I_{1}=\underset{I_{3}}{\underbrace{(\phi\partial_{t}(\rho(\hat{p})-\rho(p^{\mathrm{ms}})),\hat{p}-p^{\mathrm{ms}})}}+\underset{I_{4}}{\underbrace{(\phi\partial_{t}(\rho(p)-\rho(\hat{p})),\hat{p}-p^{\mathrm{ms}})}}.

(4.5)

For $I_{3}$ , we obtain

	$\displaystyle I_{3}$	$\displaystyle=\frac{d}{dt}\int_{\Omega}\phi\int_{0}^{\hat{p}-p^{\mathrm{ms}}}\rho^{\prime}(\hat{p}+\xi)\xi d\xi d\mathrm{x}-\underset{I_{5}}{\underbrace{\int_{\Omega}\phi\int_{0}^{\hat{p}-p^{\mathrm{ms}}}\rho^{\prime\prime}(\hat{p}+\xi)\partial_{t}\hat{p}\xi d\xi d\mathrm{x}}}$
		$\displaystyle\quad+\underset{I_{6}}{\underbrace{\int_{\Omega}\phi\int_{\Omega}(\rho^{\prime}(p^{\mathrm{ms}})-\rho^{\prime}(\hat{p}))\partial_{t}\hat{p}(\hat{p}-p^{\mathrm{ms}})d\mathrm{x}}}.$

Following [25, 21], we have that the terms $I_{5}$ and $I_{6}$ are bounded by $\|\hat{p}-p^{\mathrm{ms}}\|^{2}_{0}$ . We deduce that

I_{3}\geq\frac{d}{dt}\int_{\Omega}\phi\int_{0}^{\hat{p}-p^{\mathrm{ms}}}\rho^{\prime}(\hat{p}+\xi)\xi d\xi d\mathrm{x}-C_{1}\|\hat{p}-p^{\mathrm{ms}}\|^{2}_{0},

where $C_{1}$ is a positive constant independent of the mesh size. In virtue of $\rho^{\prime}$ being bounded below positively, we have

\int_{\Omega}\phi\int_{0}^{\hat{p}-p^{\mathrm{ms}}}\rho^{\prime}(\hat{p}+\xi)\xi d\xi d\mathrm{x}\geq C_{2}\|\hat{p}-p^{\mathrm{ms}}\|^{2}_{0}.

Then, for $I_{3}$ , we obtain

I_{3}=(\phi\partial_{t}(\rho(\hat{p})-\rho(p^{\mathrm{ms}})),\hat{p}-p^{\mathrm{ms}})\geq\frac{d}{dt}\|\hat{p}-p^{\mathrm{ms}}\|^{2}_{0}-C_{3}\|\hat{p}-p^{\mathrm{ms}}\|^{2}_{0}.

(4.6)

For $I_{4}$ , by using the chain rule and Young’s inequality, one can get

	$\displaystyle I_{4}$	$\displaystyle=(\phi(\rho^{\prime}(p)-\rho^{\prime}(\hat{p}))\partial_{t}\hat{p},\hat{p}-p^{\mathrm{ms}})+(\phi\rho^{\prime}(p)(\partial_{t}p-\partial_{t}\hat{p}),\hat{p}-p^{\mathrm{ms}})$
		$\displaystyle\preceq\\|p-\hat{p}\\|^{2}_{0}+\\|\partial_{t}(p-\hat{p})\\|^{2}_{0}+\\|\hat{p}-p^{\mathrm{ms}}\\|^{2}_{0}.$

Now, for $I_{2}$ , we get

I_{2}=\underset{I_{7}}{\underbrace{\left(\tfrac{\kappa}{\mu}(\rho(p^{\mathrm{ms}})\nabla(\hat{p}-p^{\mathrm{ms}}),\nabla(\hat{p}-p^{\mathrm{ms}})\right)}}+\underset{I_{8}}{\underbrace{\left(\tfrac{\kappa}{\mu}\left(\rho(p)\nabla p-\rho(p^{\mathrm{ms}})\nabla\hat{p}\right),\nabla(\hat{p}-p^{\mathrm{ms}})\right)}}

Then, for $I_{7}$ we have

I_{7}\geq C\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(\hat{p}-p^{\mathrm{ms}})\|^{2}_{0},

and for $I_{8}$ by invoking Young’s inequality, one obtains

	$\displaystyle I_{8}$	$\displaystyle=\left(\tfrac{\kappa}{\mu}(\rho(p)-\rho(p^{\mathrm{ms}}))\nabla p,\nabla(\hat{p}-p^{\mathrm{ms}})\right)+\left(\tfrac{\kappa}{\mu}(\rho(p^{\mathrm{ms}})\nabla p-\rho(p^{\mathrm{ms}})\nabla\hat{p}),\nabla(\hat{p}-p^{\mathrm{ms}})\right)$
		$\displaystyle\leq C\left(\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\\|^{2}_{0}+\\|p-p^{\mathrm{ms}}\\|^{2}_{0}\right)+\epsilon\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(\hat{p}-p^{\mathrm{ms}})\\|^{2}_{0},$

where in the last inequality we use the boundedness of $\rho$ and $\rho^{\prime}$ , and $p\in W^{1,\infty}(\Omega)$ . Combining the above estimates and taking $\epsilon$ small enough, we can obtain

	$\displaystyle\frac{d}{dt}\\|\hat{p}-p^{\mathrm{ms}}\\|^{2}_{0}+\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p^{\mathrm{ms}}-\hat{p})\\|^{2}_{0}$	$\displaystyle\preceq\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\\|^{2}_{0}+\\|p-\hat{p}\\|^{2}_{0}+\\|\partial_{t}(p-\hat{p})\\|^{2}_{0}$
		$\displaystyle\quad+\\|\hat{p}-p^{\mathrm{ms}}\\|^{2}_{0}.$

Integrating with respect to time $t$ and invoking the continuous Gronwall’s inequality [4], we can infer that

	$\displaystyle\\|(\hat{p}-p^{\mathrm{ms}})(t)\\|^{2}_{0}+\int_{0}^{T}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p^{\mathrm{ms}}-\hat{p})\\|^{2}_{0}dt$	$\displaystyle\preceq\\|(\hat{p}-p^{\mathrm{ms}})(0)\\|^{2}_{0}+\int_{0}^{T}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\\|^{2}_{0}dt$
		$\displaystyle\quad+\int_{0}^{T}\left(\\|p-\hat{p}\\|^{2}_{0}+\\|\partial_{t}(p-\hat{p})\\|^{2}_{0}\right)dt.$

Thus, we use the triangle inequality to get

	$\displaystyle\\|(p-p^{\mathrm{ms}})(t)\\|^{2}_{0}+\int_{0}^{T}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-p^{\mathrm{ms}})\\|^{2}_{0}dt$	$\displaystyle\preceq\\|(\hat{p}-p^{\mathrm{ms}})(0)\\|^{2}_{0}+\int_{0}^{T}(\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\\|^{2}_{0}dt$
		$\displaystyle\quad+\int_{0}^{T}\left(\\|p-\hat{p}\\|^{2}_{0}+\\|\partial_{t}(p-\hat{p})\\|^{2}_{0}\right)dt+\\|(p-\hat{p})(t)\\|^{2}_{0}.$

Finally, the proof is completed by using Lemmas 4.1 and 4.2. ∎

4.2 A posteriori error estimate

We shall give an a posteriori error estimate, which provides a computable error bound to access the quality of the numerical solution. To begin, notice that since $\mathrm{V}^{\mathrm{ms}}\subset\mathrm{V}^{h}$ , we can derive from the fully-discrete approximation a residual expression defined by

r^{\mathrm{ms}}_{n}(v):=(\phi\rho(p^{\mathrm{ms}}_{n}),v)-(\phi\rho(p^{\mathrm{ms}}_{n-1}),v)+\Delta_{n}\left(\tfrac{\kappa}{\mu}\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n},\nabla v\right)-\Delta_{n}(q,v),\quad\mbox{for each }v\in\mathrm{V}^{\mathrm{ms}}.

(4.7)

We also consider, the local residuals. For each coarse node $x_{i}$ , we define $\omega_{i}$ be the set of coarse blocks having the vertex $x_{i}$ . For each coarse neighborhood $\omega_{i}$ , we define the local residual functional $r_{i}:\mathrm{V}\to\mathbb{R}$ by

r_{n}^{(i)}(v)=r(\chi_{i}v)=(\phi\rho(p^{\mathrm{ms}}_{n}),\chi_{i}v)-(\phi\rho(p^{\mathrm{ms}}_{n-1}),\chi_{i}v)+\Delta_{n}\left(\tfrac{\kappa}{\mu}\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n},\nabla\chi_{i}v\right)-\Delta_{n}(q,\chi_{i}v),

for all $v\in\mathrm{V}$ . The local residual $r_{i}$ gives a measure of the error $p-p^{\mathrm{ms}}_{n}$ in the coarse neighborhood $\omega_{i}$ .

Theorem 4.4.

Let $p_{n}$ be the solution obtained from (2.1) at $t_{n}$ and $p^{\mathrm{ms}}_{n}\in\mathrm{V}^{\mathrm{ms}}$ denote the CEM-GMsFEM solution of the fully discrete scheme of (3.4) at $t_{n}$ . Then, there exists a positive constant $C$ independent of the mesh size such that

\|p_{N_{\tiny{\rm{time}}}}-p^{\mathrm{ms}}_{N_{\tiny{\rm{time}}}}\|^{2}_{0}+\sum_{n=1}^{N_{\tiny{\rm{time}}}}\Delta_{n}\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p_{n}-p^{\mathrm{ms}}_{n})\|^{2}_{0}\preceq(1+\Lambda^{-1})\sum_{n=1}^{N_{\tiny{\rm{time}}}}\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\|\widetilde{r}^{n}_{i}\|^{2}_{\mathrm{V}^{*}_{i}}+\|p_{0}-p^{\mathrm{ms}}_{0}\|^{2}_{0},

where

\widetilde{r}_{n}^{(i)}=\Delta_{n}\int_{\omega_{i}}q_{n}vd\mathrm{x}-\int_{\omega_{i}}\phi(\rho(p^{\mathrm{ms}}_{n})-\rho(p^{\mathrm{ms}}_{n-1}))vd\mathrm{x}-\Delta_{n}\int_{\omega_{i}}\tfrac{\kappa}{\mu}\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n}\cdot\nabla vd\mathrm{x},

and the residual norm is defined by

\|\widetilde{r}_{n}^{(i)}\|_{\mathrm{V}^{*}_{i}}=\sup_{v\in\mathrm{L}^{2}(t_{n},t_{n+1};\mathrm{H}^{1}_{0}(\omega_{i}))}\frac{\widetilde{r}_{n}^{(i)}}{\|v\|_{\mathrm{V}_{i}}}.

Proof.

Subtracting (3.4) from (2.4), we get for $p\in\mathrm{V}$ at $t_{n}$

(\phi(\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})),v)-(\phi(\rho(p_{n-1})-\rho(p^{\mathrm{ms}}_{n-1}),v)+\Delta_{n}\left(\tfrac{\kappa}{\mu}(\rho(p_{n})\nabla p_{n}-\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n},\nabla v\right)=0,

(4.8)

for each $v\in\mathrm{V}^{\mathrm{ms}}$ Putting $v=p_{n}-p^{\mathrm{ms}}_{n}$ and using the fact that $\rho$ is bounded below positively, we easily obtain that

(\phi(\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})),p_{n}-p^{\mathrm{ms}}_{n})\geq C\|p_{n}-p^{\mathrm{ms}}_{n}\|^{2}_{0}.

Similarly, for the second term of (4.8), we can use the boundedness of $\rho$ and the Young’s inequality to yield

(\phi(\rho(p_{n-1})-\rho(p^{\mathrm{ms}}_{n-1}),p_{n}-p^{\mathrm{ms}}_{n})\leq C\|p_{n-1}-p^{\mathrm{ms}}_{n-1}\|^{2}_{0}+\epsilon\|p_{n}-p^{\mathrm{ms}}_{n}\|^{2}_{0}.

Gathering the above inequalities and for $\epsilon$ small enough, we arrive to

\begin{split}(\phi(\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})),p_{n}-p^{\mathrm{ms}}_{n})-(\phi(\rho(p_{n-1})-\rho(p^{\mathrm{ms}}_{n-1}),p_{n}-p^{\mathrm{ms}}_{n})&\geq C\left(\|p_{n}-p^{\mathrm{ms}}_{n}\|^{2}_{0}\right.\\ &\quad-\left.\|p_{n-1}-p^{\mathrm{ms}}_{n-1}\|^{2}_{0}\right).\end{split}

For third term of (4.8), we have that

\left(\tfrac{\kappa}{\mu}\left(\rho(p_{n})\nabla p_{n}-\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n}\right),\nabla(p_{n}-p^{\mathrm{ms}}_{n})\right)\geq C\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p_{n}-p^{\mathrm{ms}}_{n})\|_{0}^{2}.

Thus, these above inequalities drive us the expression

	$\displaystyle\\|p_{n}-p^{\mathrm{ms}}_{n}\\|^{2}_{0}$	$\displaystyle-\\|p_{n-1}-p^{\mathrm{ms}}_{n-1}\\|^{2}_{0}+\Delta_{n}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p_{n}-p^{\mathrm{ms}}_{n})\\|^{2}_{0}\preceq(\phi(\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})),p_{n}-p^{\mathrm{ms}}_{n})$
		$\displaystyle\quad-(\phi(\rho(p_{n-1})-\rho(p^{\mathrm{ms}}_{n-1})),p_{n}-p^{\mathrm{ms}}_{n})$
		$\displaystyle\quad+\Delta_{n}\left(\tfrac{\kappa}{\mu}\left(\rho(p_{n})\nabla p_{n}-\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n}\right),\nabla(p_{n}-p^{\mathrm{ms}}_{n})\right)$
		$\displaystyle=(\phi(\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})),p_{n}-p^{\mathrm{ms}}_{n})-(\phi(\rho(p_{n-1})-\rho(p^{\mathrm{ms}}_{n-1})),p_{n}-p^{\mathrm{ms}}_{n})$
		$\displaystyle\quad+\Delta_{n}\left(\tfrac{\kappa}{\mu}\left(\rho(p_{n})\nabla p_{n}-\rho(p_{n})\nabla p^{\mathrm{ms}}_{n}\right),\nabla(p_{n}-p^{\mathrm{ms}}_{n})\right)$
		$\displaystyle\quad+\Delta_{n}\left(\tfrac{\kappa}{\mu}\left(\rho(p_{n})\nabla p^{\mathrm{ms}}_{n}-\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n}\right),\nabla(p_{n}-p^{\mathrm{ms}}_{n})\right).$

Reorganizing the terms and using the definition of the weak formulation (2.1), we get that

\begin{split}&\|p_{n}-p^{\mathrm{ms}}_{n}\|^{2}_{0}-\|p_{n-1}-p^{\mathrm{ms}}_{n-1}\|^{2}_{0}+\Delta_{n}\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p_{n}-p^{\mathrm{ms}}_{n})\|^{2}_{0}\preceq\Delta_{n}(q^{n},p_{n}-p^{\mathrm{ms}}_{n})\\ &\quad-(\phi(\rho(p^{\mathrm{ms}}_{n})-\rho(p^{\mathrm{ms}}_{n-1})),p_{n}-p^{\mathrm{ms}}_{n})-\Delta_{n}\left(\tfrac{\kappa}{\mu}\rho(p_{n})\nabla p^{\mathrm{ms}}_{n},\nabla(p_{n}-p^{\mathrm{ms}}_{n})\right)\\ &\quad+\Delta_{n}\left(\tfrac{\kappa}{\mu}\left(\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})\right)\nabla p^{\mathrm{ms}}_{n},\nabla(p_{n}-p^{\mathrm{ms}}_{n})\right).\end{split}

(4.9)

We will limit the right-hand side of (4.9). In light of residual expression (4.7), we have that

r^{\mathrm{ms}}(v)=0,\quad\mbox{for each }v\in\mathrm{V}^{\mathrm{ms}}.

Denote $w=p^{\mathrm{ms}}_{n}-p_{n}$ and use $\hat{w}\in\mathrm{V}^{\mathrm{ms}}$ , the elliptic projection of $w$ . Thus,

	$\displaystyle r^{\mathrm{ms}}_{n}(w)=r^{\mathrm{ms}}_{n}(w-\hat{w})$	$\displaystyle=\Delta_{n}(q_{n},w-\hat{w})-(\phi(\rho(p^{\mathrm{ms}}_{n})-\rho(p^{\mathrm{ms}}_{n-1})),w-\hat{w})$
		$\displaystyle\quad-\Delta_{n}\left(\tfrac{\kappa}{\mu}\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n},\nabla(w-\hat{w})\right).$

Let us rewrite $r^{\mathrm{ms}}_{n}(w-\hat{w})=\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\widetilde{r}_{n}^{(i)}(\chi_{i}(w-\hat{w}))$ [10, 22], then

\begin{split}\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\widetilde{r}_{n}^{(i)}(\chi_{i}(w-\hat{w}))&=\Delta_{n}\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\left(\int_{\omega_{i}}q_{n}\chi_{i}(w-\hat{w})d\mathrm{x}-\int_{\omega_{i}}\phi\frac{\rho(p^{\mathrm{ms}}_{n})-\rho(p^{\mathrm{ms}}_{n-1})}{\Delta_{n}}\chi_{i}(w-\hat{w})d\mathrm{x}\right.\\ &\quad-\left.\int_{\omega_{i}}\tfrac{\kappa}{\mu}\rho(p^{\mathrm{ms}}_{n})\nabla p^{\mathrm{ms}}_{n}\cdot\nabla\chi_{i}(w-\hat{w})d\mathrm{x}\right).\end{split}

(4.10)

Note that

\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\widetilde{r}^{n}_{i}(\chi_{i}(w-\hat{w}))\leq\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\|\widetilde{r}^{n}_{i}\|_{\mathrm{V}^{*}_{i}}\|\chi_{i}(w-\hat{w})\|_{\mathrm{V}_{i}},

(4.11)

where $\|v\|_{\mathrm{V}_{i}}=\|v\|^{2}_{0,\omega_{i}}+\Delta_{n}\|(\tfrac{\kappa}{\mu})^{1/2}\nabla v\|^{2}_{0,\omega_{i}}$ , where $\|\cdot\|_{0,\omega_{i}}$ denotes the $\mathrm{L}^{2}$ -norm restricted to $\omega_{i}$ . Notice also that,

\|(\tfrac{\kappa}{\mu})^{1/2}\nabla\chi_{i}(w-\hat{w})\|_{0,\omega_{i}}\preceq\left(\|w-\hat{w}\|^{2}_{s(\omega_{i})}+\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(w-\hat{w})\|^{2}_{0,\omega_{i}}\right)^{1/2},

(4.12)

where $\|\cdot\|_{s(\omega_{i})}$ represents the $s$ -norm restricted to $\omega_{i}$ . For the second term on the right-hand side, by using the orthogonality property, i.e., $((\tfrac{\kappa}{\mu})^{1/2}\rho(p_{0})\nabla(w-\hat{w}),\nabla v)=0$ , for all $v\in\mathrm{V}^{\mathrm{ms}}$ , we get

\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(w-\hat{w})\|^{2}_{0,\omega_{i}}\preceq\|(\tfrac{\kappa}{\mu})^{1/2}\nabla w\|^{2}_{0,\omega_{i}}.

(4.13)

Now, for the first term on the right-hand side of (4.12), we shall use the duality argument. Let $g=\widetilde{\kappa}(w-\hat{w})$ and $z\in\mathrm{V}_{0}$ be the solution of problem below

\int_{\omega_{i}}\kappa\rho(p_{0})\nabla z\cdot\nabla vd\mathrm{x}=\int_{\omega_{i}}gvd\mathrm{x},\quad\mbox{for each }v\in\mathrm{V}_{0}.

Putting $v=w-\hat{w}$ , using the Cauchy-Schwarz inequality and equation (4.13), we arrives to

	$\displaystyle\\|w-\hat{w}\\|^{2}_{s(\omega_{i})}$	$\displaystyle=\int_{\omega_{i}}g(w-\hat{w})d\mathrm{x}=\int_{\omega_{i}}\kappa\rho(p_{0})\nabla z\cdot\nabla(w-\hat{w})d\mathrm{x}$
		$\displaystyle=\int_{\omega_{i}}\kappa\rho(p_{0})\nabla(z-\hat{z})\cdot\nabla(w-\hat{w})d\mathrm{x}$
		$\displaystyle\leq\\|(\kappa\rho(p_{0}))^{1/2}\nabla(z-\hat{z})\\|_{0,\omega_{i}}\\|(\kappa\rho(p_{0}))^{1/2}\nabla(w-\hat{w})\\|_{0,\omega_{i}}$
		$\displaystyle\preceq\Lambda^{-1/2}\\|\widetilde{\kappa}^{-1/2}g\\|_{\mathrm{L}^{2}(\omega_{i})}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(w-\hat{w})\\|_{0,\omega_{i}}$
		$\displaystyle\preceq\Lambda^{-1/2}\\|w-\hat{w}\\|_{s(\omega_{i})}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla w\\|_{0,\omega_{i}}.$

So, we have

\|w-\hat{w}\|_{s(\omega_{i})}\preceq\Lambda^{-1/2}\|(\tfrac{\kappa}{\mu})^{1/2}\nabla w\|_{0,\omega_{i}}.

(4.14)

Gathering (4.12)-(4.14), we arrive to

\|(\tfrac{\kappa}{\mu})^{1/2}\nabla\chi_{i}(w-\hat{w})\|_{0,\omega_{i}}\preceq(1+\Lambda^{-1})^{1/2}\|(\tfrac{\kappa}{\mu})^{1/2}\nabla w\|_{0,\omega_{i}}.

Analogously, we estimate $\|\chi_{i}(w-\hat{w})\|_{0,\omega_{i}}\preceq\|(\tfrac{\kappa}{\mu})^{1/2}\nabla w\|_{0,\omega_{i}}$ . Therefore,

\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\|\widetilde{r}_{n}^{(i)}\|_{\mathrm{V}^{*}_{i}}\|\chi_{i}(w-\hat{w})\|_{\mathrm{V}_{i}}\preceq(1+\Lambda^{-1})^{1/2}\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\|\widetilde{r}_{n}^{(i)}\|_{\mathrm{V}^{*}_{i}}\|(\tfrac{\kappa}{\mu})^{1/2}\nabla w\|_{0,\omega_{i}}.

(4.15)

For the last term of the right-hand side of (4.9), we have

\Delta_{n}\left(\tfrac{\kappa}{\mu}\left(\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})\right)\nabla p^{\mathrm{ms}}_{n},\nabla(p_{n}-p^{\mathrm{ms}}_{n})\right)\preceq\Delta_{n}\|\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})\|_{0}\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p_{n}-p^{\mathrm{ms}}_{n})\|_{0}.

(4.16)

Combining (4.15) and (4.16), and using Young’s inequality, one can express for (4.9) by summing over all $n$

	$\displaystyle\\|p_{N_{\tiny{\rm{time}}}}-p^{\mathrm{ms}}_{N_{\tiny{\rm{time}}}}\\|^{2}_{0}+\sum_{n=1}^{N_{\tiny{\rm{time}}}}\Delta_{n}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p_{n}-p^{\mathrm{ms}}_{n})\\|^{2}_{0}$	$\displaystyle\preceq(1+\Lambda^{-1})\sum_{n=1}^{N_{\tiny{\rm{time}}}}\sum_{i=1}^{N_{\tiny{\rm{coarse}}}}\\|\widetilde{r}_{n}^{(i)}\\|^{2}_{\mathrm{V}^{*}_{i}}$
		$\displaystyle\quad+\sum_{n=1}^{N_{\tiny{\rm{time}}}}\Delta_{n}\\|\rho(p_{n})-\rho(p^{\mathrm{ms}}_{n})\\|_{0}+\\|p_{0}-p^{\mathrm{ms}}_{0}\\|^{2}_{0}.$

The proof is completed by using the discrete Gronwall inequality. ∎

5 Numerical results

We now present numerical results of the nonlinear single-phase compressible flow in highly heterogeneous porous media with the performance of the CEM-GMsFEM that are summarized in the following three separate experiments. All parameters on the flow model, boundary, and initial conditions in each numerical experiment are described in detail to allow their proper reproduction of them. The main aim of the simulation is to demonstrate the viability of the proposed numerical approximation and improve the convergence rate shown in Section 4. We implement the CEM-GMsFEM in Matlab language and use the numerical experiments presented in [27, 15] as a reference guide to our three-dimensional experiments. We will use the Euler-backward for the time discretization and a Newton-Raphson method with a tolerance of $10^{-6}$ for the non-linear problem. Only $2-4$ Newton’s iterations are needed in the computations presented below.

We consider three high-contrast permeability fields that are the disjoint union of a background region with $10^{5}$ millidarcys and other regions of $10^{9}$ millidarcys (see Figure 2). We also consider a fractured porous medium. In this case, the permeability value in fractures is much larger than in the surrounding medium. Finally, we employ the first $30$ layers of SPE10 $3$ D dataset from [7], which is widely used in the reservoir simulation community to test multiscale approaches. All experiments employ parameters of viscosity $\mu=5$ cP, porosity $\phi=500$ , fluid compressibility $c=1.0\times 10^{-8}\text{Pa}^{-1}$ , the reference pressure $p_{\mathrm{ref}}=2.00\times 10^{7}$ Pa, and the reference density $\rho_{\mathrm{ref}}=850$ kg/ $\mbox{m}^{3}$ .

Example 5.1.

For our first example, we set a fine grid resolution of $64^{3}$ , with a size of $h=20$ m, and a coarse grid resolution of $8^{3}$ of size $H=8h$ . The coarsening factor is chosen due this coarse grid provides the most computationally efficient performance for the method. For the CEM-GMsFEM, we use $4$ basis functions and $4$ oversampling layers. We know well that the number of bases is sufficient to improve the accuracy of CEM-GMsFEM; see [9]. Then, we have a coarse system with dimension $4916\,(=729\times\mbox{number of basis functions})$ , and the fine-scale system has a dimension of $274625$ . The permeability field $\kappa_{1}$ used in this experiment is depicted in Figure 2a. We define a model configuration as follows: four vertical injectors in each corner and a unit sink in the center of the domain to drive the flow, and employ the full zero Neumann boundary condition and an initial pressure field $p_{0}=2.16\times 10^{7}$ Pa. We consider a fine grid resolution of $64^{3}$ , whose fine grid size is given by $h=20$ m, meanwhile the coarse grid resolution of $8^{3}$ , whose coarse size is given by $H=8h$ . The time step $\Delta_{n}$ is $7$ days, and the total time simulation will be $T=25\Delta_{n}(=175\mbox{ days})$ . Figure 3 shows the pressure profiles with the sink term and zero Neumann boundary conditions in Fig. 3 for different instants at day $t=77$ and $t=147$ . In this case, we obtain a relative $L^{2}$ - and $H^{1}$ -error of 2.1138E-03 and 3.8058E-02 respectively.

Example 5.2.

We consider a combination of zero Neumann and nonzero Dirichlet boundary conditions as in [27, 15]. We set a fine grid resolution of $64^{3}$ , with a size of $h=20$ m, and different coarse grid resolutions of $4^{3},8^{3}$ and $16^{3}$ . The time step $\Delta_{n}$ and total time simulation are the same as Example 5.1. We impose zero Neumann condition on boundaries of planes $xy$ and $xz$ and let $p=2.16\times 10^{7}$ Pa in the first $yz$ plane and $p=2.00\times 10^{7}$ Pa in the last $yz$ plane for all time instants, no additional source is imposed. The permeability field used is $\kappa_{2}$ (Figure 2b). The pressure difference will drive the flow, and the initial field $p_{0}$ linearly decreases along the $x$ axis and is fixed in the $yz$ plane. Table 1 shows that numerical results use $4$ basis functions on each coarse block with different coarse grid sizes ( $H=4h,8h$ and $16h$ ), where $\varepsilon_{0}$ and $\varepsilon_{1}$ denote the relative $L^{2}$ and energy error estimate between the reference solution and CEM-GMsFEM solution defined by

\varepsilon_{0}=\left(\frac{\sum_{i=1}^{N_{\tiny{\rm{time}}}}(p^{h}_{i}-p^{\mathrm{ms}}_{i})^{2}}{\sum_{i=1}^{N_{\tiny{\rm{time}}}}(p^{h}_{i})^{2}}\right)^{1/2},\quad\varepsilon_{1}=\left(\frac{\sum_{i=1}^{N_{\tiny{\rm{time}}}}(\frac{\kappa}{\mu})^{1/2}\rho(p_{0})\nabla(p^{h}_{i}-p^{\mathrm{ms}}_{i})^{2}}{\sum_{i=1}^{N_{\tiny{\rm{time}}}}(p^{h}_{i})^{2}}\right)^{1/2},

where $p^{h}_{i}$ denotes the references solution and $p^{\mathrm{ms}}_{i}$ is the CEM-GMsFEM approximation for $i=1,\dots,N_{\tiny{\rm{time}}}$ . For instance, for a coarse grid size of $H=8h$ , we obtain the relative errors $\varepsilon_{0}=$ 4.5581E-04 and $\varepsilon_{1}=$ 2.5249E-01. In Figure 4, we depict the numerical solution profiles with a fine grid resolution of $64^{3}$ and coarse grid resolution of $8^{3}$ at day $t=105$ and $t=140$ , which is hard to find any difference between the reference solution and CEM-GMsFEM solution. Therefore, we have a good agreement.

Number of basis	$H$	Number of oversampling layers $m$	$\varepsilon_{0}$	$\varepsilon_{1}$
4	$4h$	$3$	2.4004E-03	4.4227E-01
4	$8h$	$4$	4.5581E-04	2.5249E-01
4	$16h$	$5$	1.4257E-04	1.258E-01

Table 1: Convergence rate for Example 5.2 with different numbers of oversampling layers (

m

) with a combination of zero Neumann and nonzero Dirichlet boundary conditions.

Example 5.3.

For the third experiment, we consider the combination of zero Neumann and nonzero Dirichlet boundary conditions as in Example 5.2. We set a fine grid resolution of $32^{3}$ (fine-scale system with dimension $35937$ ), with a size of $h=20$ m, and different coarse grid resolutions of $4^{3},8^{3}$ and $16^{3}$ (coarse-scale system with dimension $500$ , $2916$ and $19652$ respectively). The time step $\Delta_{n}$ and total time simulation are the same as Example 5.1. The fractured medium $\kappa_{3}$ used is depicted in Figure 2c. For this experiment, we employ the framework from [14] and apply the CEM-GMsFEM to the $3$ D model. The domain $\Omega$ can be represented by

\Omega=\Omega_{0}\cup(\cup_{i}\Omega_{\mathrm{frac},i}),

where $\Omega_{0}$ represents the matrix and subscript $\mathrm{frac}$ denotes the fracture regions. Then, we can write the finite element discretization corresponding to equation (2.3)

	$\displaystyle(\phi\partial_{t}\rho(p^{h}),v)+\left(\tfrac{\kappa}{\mu}\rho(p^{h})\nabla p^{h},\nabla v\right)$	$\displaystyle=(\phi\partial_{t}\rho(p^{h}),v)_{\Omega_{0}}+\sum_{i}(\phi\partial_{t}\rho(p^{h}),v)_{\Omega_{\mathrm{frac},i}}$
		$\displaystyle\quad+\left(\tfrac{\kappa}{\mu}\rho(p^{h})\nabla p^{h},\nabla v\right)_{\Omega_{0}}+\sum_{i}\left(\tfrac{\kappa}{\mu}\rho(p^{h})\nabla p^{h},\nabla v\right)_{\Omega_{\mathrm{frac},i}}$
		$\displaystyle=(q,v),\quad\mbox{for each }v\in\mathrm{V}^{h},$

In Table 2, we give the convergence rate $T=25\Delta_{n}(=175\text{ days})$ with different coarse-grid sizes $H$ . We notice that the error significantly decreases as the size of the coarse grid is finer. Then, we have that the CEM-GMsFEM gives a good approximation of the solution for the case of the fractured medium. Figure 5 shows the numerical solutions at $t=70$ and $t=140$ .

Number of basis	$H$	Number of oversampling layers $m$	$\varepsilon_{0}$	$\varepsilon_{1}$
4	$4h$	$3$	1.6111E-03	2.9441E-01
4	$8h$	$4$	1.9532E-04	1.1653E-01
4	$16h$	$5$	1.0160E-04	5.1492E-02

Table 2: Convergence rate of Example 5.3 with different numbers of oversampling layers (

m

) with a combination of zero Neumann and nonzero Dirichlet boundary conditions.

Example 5.4.

For the last experiment, we consider the same boundary conditions as Example 5.2. The permeability field used is $\kappa_{4}$ (see Figure LABEL:fig:perm4) and the time step $\Delta_{n}$ is $7$ days, with total time simulation will be $T=26\Delta_{n}(=186\mbox{ days})$ . We set a fine grid resolution of $220\times 60\times 30$ (fine-scale system with a dimension of $417911$ ), with a size of $h=20$ m, and coarse grid resolution of $10^{3}$ (coarse-scale system with a dimension $5324$ ). We show the pressure profiles comparison in Figure 6. This experiment obtains an error estimate of $\varepsilon_{0}=$ 1.8377E-03 and $\varepsilon_{1}=$ 3.4547E-01, respectively.

6 Concluding remarks

This paper studies the convergence of the numerical approximations to the highly heterogeneous nonlinear single-phase compressible flow by CEM-GMsFEM. We first build an auxiliary for the proposed method by solving spectral problems. Then, we construct a multiscale basis function by solving some constraint energy minimization problems in the oversampling local regions. So, we obtain multiscale basis functions for the pressure. This work defines the elliptic projection in the multiscale space spanned by CEM-GMsFEM basis functions for convergence analysis. Thus, we present the convergence of the semi-discrete formulation. The convergence depends on the coarse mesh size and the decay of eigenvalues of local spectral problems; an a posteriori error estimate is derived underlining discretization. Some numerical examples have been presented to verify the feasibility of the proposed method concerning convergence and stability. We observe that the CEM-GMsFEM is shown to have a second-order convergence rate in the $L^{2}$ -norm and a first-order convergence rate in the energy norm concerning the coarse grid size.

A foreseeable result in ongoing research is to boost the performance of the coarse-grid simulation, mainly where the source term is singular; one may need to further improve the accuracy of the approximation without additional refinement in the grid. We can enrich the multiscale space for such a goal by adding additional basis functions in the online stage [10]. These new multiscale basis functions are constructed on the oversampling technique and the information on local residuals. Consequently, we could present an adaptive enrichment algorithm to reduce error in some regions with large residuals.

Acknowledgement

The research of Eric Chung is partially supported by the Hong Kong RGC General Research Fund (Projects: 14305222 and 14304021).

References

[1] Manal Alotaibi, Huangxin Chen, and Shuyu Sun. Generalized multiscale finite element methods for the reduced model of darcy flow in fractured porous media. Journal of Computational and Applied Mathematics, 413:114305, 2022.
[2] T. Arbogast, G. Pencheva, M. F. Wheeler, and I. Yotov. A multiscale mortar mixed finite element method. Multiscale Model. Simul., 6(1):319–346, 2007.
[3] I. Babuška and J. M. Melenk. The partition of unity method. Internat. J. Numer. Methods Engrg., 40(4):727–758, 1997.
[4] John R. Cannon, Richard E. Ewing, Yinnian He, and Yanping Lin. A modified nonlinear Galerkin method for the viscoelastic fluid motion equations. Internat. J. Engrg. Sci., 37(13):1643–1662, 1999.
[5] Jie Chen, Eric T Chung, Zhengkang He, and Shuyu Sun. Generalized multiscale approximation of mixed finite elements with velocity elimination for subsurface flow. Journal of Computational Physics, 404:109133, 2020.
[6] Zhiming Chen and Thomas Y. Hou. A mixed multiscale finite element method for elliptic problems with oscillating coefficients. Math. Comp., 72(242):541–576, 2003.
[7] M. A. Christie and M. J. Blunt. Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques. SPE Reservoir Evaluation & Engineering, 4(04):308–317, August 2001. _eprint: https://onepetro.org/REE/article-pdf/4/04/308/2586053/spe-72469-pa.pdf.
[8] Eric T. Chung, Yalchin Efendiev, and Shubin Fu. Generalized multiscale finite element method for elasticity equations. GEM Int. J. Geomath., 5(2):225–254, 2014.
[9] Eric T. Chung, Yalchin Efendiev, and Wing Tat Leung. Constraint energy minimizing generalized multiscale finite element method. Comput. Methods Appl. Mech. Engrg., 339:298–319, 2018.
[10] Eric T. Chung, Yalchin Efendiev, and Wing Tat Leung. Fast online generalized multiscale finite element method using constraint energy minimization. J. Comput. Phys., 355:450–463, 2018.
[11] Eric T. Chung, Yalchin Efendiev, and Guanglian Li. An adaptive GMsFEM for high-contrast flow problems. J. Comput. Phys., 273:54–76, 2014.
[12] Yalchin Efendiev, Juan Galvis, and Thomas Y. Hou. Generalized multiscale finite element methods (GMsFEM). J. Comput. Phys., 251:116–135, 2013.
[13] Yalchin Efendiev, Juan Galvis, Guanglian Li, and Michael Presho. Generalized multiscale finite element methods: Oversampling strategies. International Journal for Multiscale Computational Engineering, 12(6), 2014.
[14] Yalchin Efendiev, Seong Lee, Guanglian Li, Jun Yao, and Na Zhang. Hierarchical multiscale modeling for flows in fractured media using generalized multiscale finite element method. GEM Int. J. Geomath., 6(2):141–162, 2015.
[15] Shubin Fu, Eric Chung, and Lina Zhao. Generalized Multiscale Finite Element Method for Highly Heterogeneous Compressible Flow. Multiscale Model. Simul., 20(4):1437–1467, 2022.
[16] Hadi Hajibeygi and Patrick Jenny. Multiscale finite-volume method for parabolic problems arising from compressible multiphase flow in porous media. J. Comput. Phys., 228(14):5129–5147, 2009.
[17] Thomas Y. Hou and Xiao-Hui Wu. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys., 134(1):169–189, 1997.
[18] Thomas J. R. Hughes, Gonzalo R. Feijóo, Luca Mazzei, and Jean-Baptiste Quincy. The variational multiscale method—a paradigm for computational mechanics. Comput. Methods Appl. Mech. Engrg., 166(1-2):3–24, 1998.
[19] Patrick Jenny, SH Lee, and Hamdi A Tchelepi. Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. Journal of computational physics, 187(1):47–67, 2003.
[20] Patrick Jenny, SH Lee, and Hamdi A Tchelepi. Adaptive multiscale finite-volume method for multiphase flow and transport in porous media. Multiscale Model. Simul., 3(1):50–64, 2004.
[21] Mi-Young Kim, Eun-Jae Park, Sunil G. Thomas, and Mary F. Wheeler. A multiscale mortar mixed finite element method for slightly compressible flows in porous media. J. Korean Math. Soc., 44(5):1103–1119, 2007.
[22] Mengnan Li, Eric Chung, and Lijian Jiang. A constraint energy minimizing generalized multiscale finite element method for parabolic equations. Multiscale Model. Simul., 17(3):996–1018, 2019.
[23] Axel Må lqvist and Daniel Peterseim. Localization of elliptic multiscale problems. Math. Comp., 83(290):2583–2603, 2014.
[24] Ana M. Matache and Christoph Schwab. Homogenization via $p$ -FEM for problems with microstructure. In Proceedings of the Fourth International Conference on Spectral and High Order Methods (ICOSAHOM 1998) (Herzliya), volume 33, pages 43–59, 2000.
[25] Eun-Jae Park. Mixed finite element methods for generalized Forchheimer flow in porous media. Numer. Methods Partial Differential Equations, 21(2):213–228, 2005.
[26] Maria Vasilyeva, Eric T Chung, Siu Wun Cheung, Yating Wang, and Georgy Prokopev. Nonlocal multicontinua upscaling for multicontinua flow problems in fractured porous media. Journal of Computational and Applied Mathematics, 355:258–267, 2019.
[27] Yixuan Wang, Hadi Hajibeygi, and Hamdi A Tchelepi. Algebraic multiscale solver for flow in heterogeneous porous media. Journal of Computational Physics, 259:284–303, 2014.

	$\displaystyle\\|p-\hat{p}\\|_{0}$	$\displaystyle\preceq\Lambda^{-1}H^{2}\\|\kappa^{-1/2}(q-\phi\partial_{t}\rho(p))\\|_{0},$
	$\displaystyle\\|\partial_{t}(p-\hat{p})\\|_{0}$	$\displaystyle\preceq\Lambda^{-1}H^{2}\\|\kappa^{-1/2}\partial_{t}(q-\phi\partial_{t}\rho(p))\\|_{0}.$

	$\displaystyle\\|p-\hat{p}\\|^{2}_{0}=\int_{\Omega}\kappa\rho(p_{0})\nabla w\cdot\nabla(p-\hat{p})d\mathrm{x}$	$\displaystyle=\int_{\Omega}\kappa\rho(p_{0})\nabla(w-\hat{w})\cdot\nabla(p-\hat{p})d\mathrm{x}$
		$\displaystyle\preceq\int_{\Omega}\tfrac{\kappa}{\mu}\rho(p_{0})\nabla(w-\hat{w})\cdot\nabla(p-\hat{p})d\mathrm{x}$
		$\displaystyle\leq\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(w-\hat{w})\\|_{0}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-\hat{p})\\|_{0}$
		$\displaystyle\preceq\left(H\Lambda^{-1/2}\max\{\kappa^{-1/2}\}\\|p-\hat{p}\\|_{0}\right)$
		$\displaystyle\quad\times\left(H\Lambda^{-1/2}\\|\kappa^{-1/2}(q-\phi\partial_{t}\rho(p))\\|_{0}\right).$

	$\displaystyle\\|(p-p^{\mathrm{ms}})(t)\\|_{0}^{2}+\int_{0}^{T}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(p-p^{\mathrm{ms}})\\|^{2}_{0}dt$	$\displaystyle\preceq\Lambda^{-1}H^{2}\left(\\|\kappa^{-1/2}\partial_{t}(q-\phi\partial_{t}\rho(p))(t)\\|_{0}^{2}\right.$
		$\displaystyle\quad+\left.\int_{0}^{T}\\|\kappa^{-1/2}\partial_{t}(q-\phi\partial_{t}\rho(p))\\|^{2}_{0}dt\right)$
		$\displaystyle\quad+\\|(\hat{p}-p^{\mathrm{ms}})(0)\\|^{2}_{0}.$

	$\displaystyle\\|w-\hat{w}\\|^{2}_{s(\omega_{i})}$	$\displaystyle=\int_{\omega_{i}}g(w-\hat{w})d\mathrm{x}=\int_{\omega_{i}}\kappa\rho(p_{0})\nabla z\cdot\nabla(w-\hat{w})d\mathrm{x}$
		$\displaystyle=\int_{\omega_{i}}\kappa\rho(p_{0})\nabla(z-\hat{z})\cdot\nabla(w-\hat{w})d\mathrm{x}$
		$\displaystyle\leq\\|(\kappa\rho(p_{0}))^{1/2}\nabla(z-\hat{z})\\|_{0,\omega_{i}}\\|(\kappa\rho(p_{0}))^{1/2}\nabla(w-\hat{w})\\|_{0,\omega_{i}}$
		$\displaystyle\preceq\Lambda^{-1/2}\\|\widetilde{\kappa}^{-1/2}g\\|_{\mathrm{L}^{2}(\omega_{i})}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla(w-\hat{w})\\|_{0,\omega_{i}}$
		$\displaystyle\preceq\Lambda^{-1/2}\\|w-\hat{w}\\|_{s(\omega_{i})}\\|(\tfrac{\kappa}{\mu})^{1/2}\nabla w\\|_{0,\omega_{i}}.$