1 Introduction↩︎

Numerical simulation of incompressible and immiscible two-phase flow in porous media is an important problem in petroleum reservoir engineering, hydrology, and subsurface energy applications. The governing model is derived from mass conservation, Darcy’s law, the saturation constraint, and the capillary pressure relation. In realistic geological media, the absolute permeability usually contains multiple spatial scales, high-conductivity channels, inclusions, and very high contrast. These heterogeneous features strongly affect the pressure and velocity fields, while the saturation equation may develop sharp fronts or discontinuous profiles. Direct fine-grid simulation can resolve these effects, but it is computationally expensive when the pressure-velocity system has to be solved repeatedly during a long-time two-phase flow simulation [1]–[4].

For the temporal discretization of two-phase flow, fully implicit methods are robust but require solving large nonlinear systems at each time step. The implicit-pressure explicit-saturation (IMPES) method is more efficient since it separates the pressure solve from the saturation update [5]–[7]. However, standard IMPES schemes may suffer from restrictive stability conditions, loss of mass conservation for one phase, phase bias, and violation of saturation bounds. In heterogeneous media with capillary pressure, saturation discontinuities across material interfaces introduce further difficulties. Several improved IMPES, discontinuous Galerkin, and capillarity-treatment strategies have been proposed to address these issues [8]–[11].

Recently, Chen and Sun proposed a physics-preserving IMPES scheme, hereafter referred to as P-IMPES, for incompressible and immiscible two-phase flow in heterogeneous porous media [12]. The main idea is to rewrite the phase Darcy velocities in terms of the total velocity and an auxiliary velocity related to the capillary potential gradient. The resulting scheme is locally mass conservative for both phases, preserves the normal continuity of the total velocity, treats the wetting and non-wetting phases in an unbiased way, and keeps the saturation within physical bounds under a suitable CFL condition. These properties are especially important because the computed velocity field is directly used in the saturation transport equation.

Besides time discretization, the spatial approximation of the pressure–velocity system is another major computational challenge. Classical upscaling methods reduce fine-scale heterogeneity to effective coarse-scale coefficients [13]–[16]. Multiscale finite element and multiscale finite volume methods provide an alternative reduced-order strategy by constructing coarse-scale basis functions from local fine-scale problems [17]–[20]. Since the saturation equation is very sensitive to the quality of the velocity approximation, local mass conservation is essential for multiscale simulations of porous-media flow. This motivates mixed multiscale finite element methods, multiscale mortar methods, and conservative flux postprocessing techniques [21]–[24].

The Generalized Multiscale Finite Element Method (GMsFEM) improves the flexibility of multiscale methods by constructing multiple basis functions through local spectral decompositions [25]–[28]. The GMsFEM framework has been further developed in adaptive enrichment, oversampling, randomized and sparse model reduction, and wave and space-time applications [29]–[35]. For high-contrast media, however, important features such as high-conductivity channels may be nonlocal and cannot be captured accurately by standard localized basis functions. Related localized and operator-adapted multiscale ideas for rough-coefficient elliptic problems have also been developed in [36], [37]. To overcome the difficulty caused by high contrast and nonlocal features, the Constraint Energy Minimizing Generalized Multiscale Finite Element Method (CEM-GMsFEM) was introduced in [38]. The method constructs an auxiliary space from local spectral problems and then computes constraint energy minimizing basis functions in oversampling regions. By including the eigenfunctions associated with small contrast-dependent eigenvalues, CEM-GMsFEM captures channelized features and achieves contrast-independent convergence with respect to the coarse mesh size. Its mixed formulation was developed in [39], where auxiliary pressure basis functions are used to construct localized velocity basis functions. The resulting method provides fine-grid mass-conservative velocity approximations and first-order convergence independent of the contrast, provided that enough oversampling layers are used.

Although P-IMPES provides a physics-preserving time discretization and mixed CEM-GMsFEM provides a contrast-robust spatial discretization, their combination for two-phase flow still faces a time-dependent coefficient issue. At time level \(t_n\), the coefficient in the pressure–velocity system is \(\kappa_n=\lambda_t(S_w^n)K,\) where \(K\) is the absolute permeability and \(\lambda_t(S_w^n)\) is the total mobility depending on the current wetting-phase saturation. Hence the coefficient changes as the saturation evolves. If the multiscale spaces are constructed only from the initial coefficient and kept fixed, the resulting basis functions may no longer be well adapted to the current coefficient field, which can deteriorate the velocity approximation and the saturation front.

Compared with the P-IMPES-MsFEM method in [40], which used residual-driven enrichment and conservative postprocessing and established phasewise conservation, unbiasedness, conditional bounds preservation, and velocity–saturation error relations, the present work addresses mobility-induced coefficient obsolescence with mixed CEM-GMsFEM, whose oversampled energy-minimizing bases are designed for contrast-robust approximation. Its distinct ingredients are coefficient-variation-triggered global space reconstruction, fine-cell-residual marking for postprocessing instead of source-based marking, and an error estimate separating contributions from the update tolerance, coarse mesh, spectral truncation, and front layer. The P-IMPES properties are therefore inherited and reverified, rather than claimed as new.

In this paper, we propose an adaptive mixed CEM-GMsFEM for the physics-preserving IMPES simulation of incompressible and immiscible two-phase flow in high-contrast porous media. The method couples the P-IMPES formulation with mixed CEM-GMsFEM spaces for the total velocity, auxiliary capillary velocity, and phase pressure. To reduce the error caused by coefficient variation, we propose an error indicator \(\eta_n\) to monitor the saturation-induced variation of the mobility-weighted coefficient field. When \(\eta_n\) exceeds a prescribed tolerance, the multiscale spaces are regenerated using the current coefficient; otherwise, the existing spaces are reused. A local postprocessing step is also applied to recover fine-grid mass conservation of the total velocity.

The rest of the paper is organized as follows. Section 2 introduces the two-phase flow model and the P-IMPES discretization. Section 3 presents the mixed CEM-GMsFEM construction and the adaptive update strategy. Section 4 gives the analysis of conservation, unbiased property, bounds preservation, and error estimates. Numerical results are reported in Section 5, and conclusions are given in Section 6.

2 Preliminaries↩︎

In this section, we introduce the foundational mathematical model for incompressible and immiscible two-phase flow within a porous medium. We denote the wetting and non-wetting phases by subscripts \(w\) and \(n\), respectively. The mathematical model is derived from the conservation of mass, Darcy’s law, the capillary pressure relationship, and the saturation constraint. Following the model formulation, we introduce the P-IMPES scheme, which is unbiased with regard to the two phases and satisfies the mass conservation for both phases.

2.1 Mathematical model↩︎

We consider a model in porous media \(\Omega \subset \mathbb{R}^2\) given as follows, \[\begin{align} \phi\frac{\partial S_\alpha}{\partial t} + \nabla \cdot \mathbf{u}_\alpha = F_\alpha, &\quad \text{in } \Omega, \quad \alpha = w, n, \nonumber\label{2461}\\ \mathbf{u}_\alpha = -\frac{k_{r\alpha}}{\mu_\alpha}\mathbf{K}\nabla p_\alpha , &\quad \text{in } \Omega, \quad \alpha = w, n,\\ S_n + S_w = 1, &\quad \text{in } \Omega, \nonumber \\ p_c(S_w) = p_n - p_w, &\quad \text{in } \Omega. \nonumber \end{align}\tag{1}\]

Here \(\phi\) is the porosity of the medium, which is a constant. Let \(\mathbf{K}\), \(S_\alpha\), \(\mathbf{u}_\alpha\) and \(p_\alpha\) be the absolute permeability tensor, saturation, Darcy’s velocity and pressure corresponding to phase \(\alpha\). Total velocity is denoted by \(\mathbf{u}_t = \mathbf{u}_w + \mathbf{u}_n\). Define \(k_{r\alpha}\), \(\mu_\alpha\) and \(F_\alpha\) as relative permeability, viscosity and sink/source term of phase \(\alpha\). Besides, \(F_t = F_w + F_n\). Phase mobility is denoted by \(\lambda_\alpha = \frac{k_{r\alpha}}{\mu_\alpha}\). Total mobility is further defined by \(\lambda_t = \lambda_w + \lambda_n\). Then we define fractional flow functions as \(f_w = \lambda_w/\lambda_t\), \(f_n = \lambda_n/\lambda_t\).

We let the boundary be \(\Gamma := \partial \Omega\) and it can be partitioned as \(\Gamma = \Gamma_D \cup \Gamma_N\), where \(\Gamma_D\) and \(\Gamma_N\) are Dirichlet and Neumann boundaries to solve \(\mathbf{u}_\alpha\). Moreover, to solve \(S_\alpha\), we let \(\Gamma = \Gamma_{\mathrm{in}} \cup \Gamma_{\mathrm{out}}\), where \(\Gamma_{\mathrm{in}} = \{x \in \Gamma : \mathbf{u}_t(x) \cdot \mathbf{n}(x) < 0\}\) is the inflow boundary, \(\Gamma_{\mathrm{out}} = \{x \in \Gamma : \mathbf{u}_t(x) \cdot \mathbf{n}(x) \ge 0\}\) is the outflow boundary, and \(\mathbf{n}\) is the unit outer normal vector to \(\Gamma\). The initial and boundary conditions are \(S_\alpha = S_\alpha^0\), for \(t = 0\), \(p_\alpha = p_\alpha^B\) on \(\Gamma_D\) and \(\mathbf{u}_\alpha \cdot \mathbf{n} = g_\alpha^N\) on \(\Gamma_N\), where \(\alpha = w, n\).

2.2 Physics-preserving IMPES scheme (P-IMPES)↩︎

To compute the reference solutions, we solve Darcy flow on a fine-scale mesh \(\mathcal{T}_h\). We use the lowest Raviart-Thomas vector field (RT0) denoted by \(V_h\) for fine scale velocity bases. Besides, we define \(Q_h\) as the pressure space spanned by piecewise constant functions and each basis corresponds to a particular fine-scale element. And we use \(V_h(D)\) and \(Q_h(D)\) to represent \(V_h\) and \(Q_h\) restricted in \(D\). The P-IMPES scheme results in the following equation system.
For any test functions \(\mathbf{v} \in V_h\) and \(q \in Q_h\), the scheme satisfies: \[\begin{align} & (\kappa_n^{-1}\mathbf{u}_t^{h,n+1}, \mathbf{v}) - (p_w^{h,n+1}, \nabla \cdot \mathbf{v}) = (\kappa_n^{-1} f_n(S_w^{h,n})\xi_c^{h,n+1}, \mathbf{v}) - \int_{\Gamma_D} p_w^B \mathbf{v} \cdot \mathbf{n},\tag{2}\\ & \sum_\alpha \beta_\alpha(\mathbf{u}_t^{h,n+1}, q; S_w^{h,n}) = (F_t, q),\tag{3}\\ & \left(\phi \frac{S_\alpha^{h,n+1} - S_\alpha^{h,n}}{t_{n+1} - t_n}, q\right) + \beta_\alpha(\mathbf{u}_t^{h,n+1}, q; S_w^{h,n}) = (F_\alpha, q) + \sigma_\alpha \beta_c(\xi_c^{h,n+1}, q; S_w^{h,n}),\tag{4}\\ & (\kappa_n^{-1} \xi_c^{h,n+1}, \mathbf{v}) = (p_c(S_w^{h,n}), \nabla \cdot \mathbf{v}) - \int_{\Gamma_D} (p_n^B - p_w^B)\mathbf{v} \cdot \mathbf{n},\\ & (S_n^{h,n+1} + S_w^{h,n+1}, q) = (1, q),\\ & (p_n^{h,n+1} - p_w^{h,n+1}, q) = (p_c(S_w^{h,n}), q). \end{align}\]

Define \(t_i\) as the \(i\)-th time step in a uniform partition of \([0, T]\), where \(T\) is the final time \(\Delta t = t_{n+1} - t_n\) for each \(n\). Set \(\sigma_w = 1\) and \(\sigma_n = -1\). Let \(\mathbf{u}_t^{h, n} \in V_h(\Omega)\) and \(\xi_c^{h, n} \in V_h(\Omega)\) be reference velocity solutions to (2.3) and (2.6) at time \(t_n\). \(p_w^{h, n}\) is the reference pressure of wetting phase at time \(t_n\) which is the solution to (2.3). And \(p_n^{h, n}\) is defined to be the pressure of non-wetting phase at time \(t_n\) and it can be solved by (2.8).

Define \(\kappa_n = \lambda_t(S_w^{h, n})\mathbf{K}\). After we get the solutions for \(\mathbf{u}_t^{h,n+1}\) and \(\mathbf{\xi}_c^{h,n+1}\), we can update the wetting and non-wetting velocities \(\mathbf{u}_w^{h,n+1}\) and \(\mathbf{u}_n^{h,n+1}\) on each element as follows: \[\begin{align} \mathbf{u}_w^{h,n+1} &= f_w(S_w^{h,n})\mathbf{u}_t^{h,n+1} - f_w(S_w^{h,n})f_n(S_w^{h,n}){\xi}_c^{h,n+1}, \\ \mathbf{u}_n^{h,n+1} &= f_n(S_w^{h,n})\mathbf{u}_t^{h,n+1} + f_w(S_w^{h,n})f_n(S_w^{h,n}){\xi}_c^{h,n+1}. \end{align}\] For \(q \in Q_h\) is piecewise constant, \[\begin{align} \beta_\alpha(\mathbf{v}, q; S_w^h) &= \sum_{K^h \in \mathcal{T}_h} \int_{\partial K^h} f_\alpha(S_{w,\alpha}^{*,h})\mathbf{v} \cdot \mathbf{n} q, \quad \alpha = w, n,\\ \beta_c(\mathbf{v}, q; S_w^h) &= \sum_{K^h \in \mathcal{T}_h} \int_{\partial K^h} f_n(S_{w,n}^{*,h}) f_w(S_{w,w}^{*,h})\mathbf{v} \cdot \mathbf{n} q, \end{align}\] where the upwind value \(S_{w,\alpha}^{*,h}\) on \(e \subset \partial K^h\) in the function \(f_\alpha(S_{w,\alpha}^{*,h})\) is defined as follows: \[S_{\alpha}^{*,h}|_e = \begin{cases} S_{\alpha}^h|_{K^h}, & \text{if } \{\mathbf{u}_{\alpha}^h \cdot \mathbf{n}\}_e \ge 0, \\ S_{\alpha}^h|_{K^{h,1}}, & \text{if } \{\mathbf{u}_{\alpha}^h \cdot \mathbf{n}\}_e < 0, \end{cases} \qquad S_{w,\alpha}^{*,h} = \begin{cases} S_w^{*,h}, & \alpha = w, \\ 1 - S_n^{*,h}, & \alpha = n. \end{cases}\] Here \(K^h \cap K^{h,1} = e\) and \(\mathbf{n}\) is an outward normal vector to \(K^h\). For \(e \subset \Gamma_{\mathrm{in}}\), \(S_{w,\alpha}^{*,h}|_e = S_w^h|_{K^h}\). We further define \(\beta_t = \beta_w + \beta_n\). At the initial time \(t=0\), using the initial saturation \(S_w^0\), the numerical flux terms reduce to \[\begin{align} \beta_{\alpha}(\mathbf{v},q;S_w^{0})&=\sum_{K^h \in \mathcal{T}_h} \int_{K^h} f_{\alpha}(S_w^0)q \nabla \cdot \mathbf{v} \, dx , \qquad\alpha=w,n, \\ \beta_c(\mathbf{v},q;S_w^{0}) &=\sum_{K^h \in \mathcal{T}_h} \int_{K^h} f_w(S_w^0)f_n(S_w^0)q \nabla \cdot \mathbf{v} \, dx \end{align}\]

3 Construction of Mixed CEM-GMsFEM↩︎

In this section, we will introduce the framework of the mixed CEM-GMsFEM. The multiscale finite element method consists of two steps. First, we construct a multiscale space \(Q_{ms}\) for approximating the pressure. Based on the space \(Q_{ms}\), we construct another multiscale space \(V_{ms}\) for the velocity.

First, we introduce some notations that will be used later. Given a subset \(S \subset \Omega\), we define \(V_{h,0}(S) := \{v \in V_h \cap H(\operatorname{div}; S) : v \cdot \mathbf{n}_S = 0 \text{ on } \partial S\}\) and \(Q_h(S) := Q_h \cap L^2(S)\). Denote \(\kappa_n^H=\lambda_t(S_w^{H,n}) \mathbf{K},\;\;p_c^n=p_c(S_w^{H,n})\)

Define the bilinear forms as follows: \[\begin{align} &a_n(v,w)=\int_\Omega \kappa_n^{-1} v \cdot w \, dx, \quad a_n^H(v, w) = \int_\Omega (\kappa_n^H)^{-1} v \cdot w \, dx,\\ & b(v, q) = \int_\Omega q \nabla \cdot v \, dx, \quad \quad b_D(v,p^B)=\int_{\Gamma_D} p^B v\cdot n \,ds, \end{align}\]

Note that the pressure basis functions are locally supported on individual coarse elements. In contrast, the associated velocity basis functions are supported on oversampled regions formed by adding several coarse grid layers to the original element (see Figure 1). The localized feature of these velocity bases is the key of our approach.

3.1 Construction of pressure multiscale space↩︎

In this section, we outline the procedure for generating the pressure multiscale approximation space, denoted as \(Q_{\mathrm{ms}}\). For each coarse grid block \(K_i \in \mathcal{T}^H\), we establish a localized eigenvalue problem. Specifically, we seek the eigenpairs \((\phi_j^i, p_j^i) \in V_{h,0}(K_i) \times Q_h(K_i)\) and the corresponding eigenvalues \(\lambda_j^i \in \mathbb{R}\) satisfying: \[\begin{align} a_0^H(\phi_j^i, v) - b(v, p_j^i) &= 0 && \forall v \in V_{h,0}(K_i), \\ b(\phi_j^i, q) &= \lambda_j^i s_i^H(p_j^i, q) && \forall q \in Q_h(K_i), \end{align}\] where \(j = 1, \dots, L_i\), with the integer \(L_i \in \mathbb{N}^+\) representing the total number of localized degrees of freedom determined by the fine and coarse meshes. The localized bilinear operator \(s_i: Q_h(K_i) \times Q_h(K_i) \to \mathbb{R}\) is formulated as follows: \[s_i^H(p, q) = \int_{K_i} \tilde{\kappa}^H p q \, dx, \quad \text{where } \tilde{\kappa}^H = \kappa_0^H \sum_{j=1}^{N_c} |\nabla \chi_j|^2.\]

Here, the basis functions \(\{\chi_j\}_{j=1}^{N_c}\) constitute a specialized multiscale partition of unity (POU). To be precise, for any internal coarse node \(x_j\) and its associated overlapping neighborhood \(\omega_j = \bigcup \{K \in \mathcal{T}^H : x_j \in \partial K\}\), the POU function \(\chi_j\) is computed by solving the following local elliptic PDEs within each \(K \subset \omega_j\): \[\begin{align} -\nabla \cdot (\kappa_0^H \nabla \chi_j) &= 0 && \text{in the interior of } K \subset \omega_j, \\ \chi_j &= g_j && \text{on } \partial K \setminus \partial \omega_j \;(\text{for all } K \subset \omega_j), \\ \chi_j &= 0 && \text{on the boundary } \partial \omega_j. \end{align}\] In the boundary conditions above, \(g_j\) is a piecewise linear and continuous function along the edges of the coarse element.

Based on the normalization, we set the constraint \(s_i^H(p_j^i, p_j^i) = 1\). Then ordering the eigenvalues such that \(0 \le \lambda_1^i \le \lambda_2^i \le \dots \le \lambda_{L_i}^i\) and truncate the spectrum, retain only the first \(J_i\) (where \(1 \le J_i \le L_i\)) eigenfunctions that correspond to the smallest eigenvalues. The global multiscale space for the pressure is defined as the linear span of these selected basis functions: \[Q_{\mathrm{ms}} = \operatorname{span}\{p_j^i : i = 1, \dots, N, \;j = 1, \dots, J_i\}.\]

3.2 Construction of velocity multiscale functions↩︎

In this subsection, we detail the derivation of the multiscale velocity space, denoted as \(V_{\mathrm{ms}}\). First, we define a global projection operator \(\pi: Q_h \to Q_{\mathrm{ms}}\) , which can be explicitly formulated as: \[\pi_s q = \sum_{i=1}^N \sum_{j=1}^{J_i} s_i^H(p_j^i, q) p_j^i, \quad \forall q \in Q_h.\]

First of all, we define a global bilinear form \(s: Q_h \times Q_h \to \mathbb{R}\) by the local bilinear form \(s_i\) as \(s(p, q) = \sum_{i=1}^N s_i^H(p, q)\). From this perspective, the operator \(\pi_s\) is an orthogonal projection mapping the fine grid space \(Q_h\) onto \(Q_{\mathrm{ms}}\) with respect to the \(s(\cdot, \cdot)\) inner product.

Global energy-minimizing solutions have decay property and are highly localized. This feature allows us to approximate the global velocity basis within a truncated oversampled region. For any given coarse grid block \(K_i \in \mathcal{T}^H\), we expand its boundary by \(\ell\) (\(\ell \in \mathbb{N}^+\)) coarse grid layers, see figure 1 below. The resulting oversampled domains are defined as: \[K_{i,0} = K_i, \quad K_{i,\ell} = \bigcup \left\{ K \in \mathcal{T}^H : K \cap \overline{K_{i,\ell-1}} \neq \emptyset \right\}, \quad \text{for } \ell = 1, 2, \dots\] For notational simplicity, we let \(K_i^+\) denote the target localized region \(K_{i,\ell}\) with a sufficiently large layer parameter \(\ell\).

For each selected pressure eigenfunction \(p_j^i \in Q_{ms}\), the corresponding multiscale velocity basis function \(\psi_{j,{ms}}^i \in V_{h,0}(K_i^+)\) along with \(q_{j,{ms}}^i \in Q_h(K_i^+)\) are determined by solving the following localized system: \[\begin{align} a_0^H(\psi_{j,{ms}}^i, v) - b(v, q_{j,{ms}}^i) &= 0 && \forall v \in V_{h,0}(K_i^+), \\ s(\pi_s q_{j,{ms}}^i, \pi q) + b(\psi_{j,{ms}}^i, q) &= s(p_j^i, q) && \forall q \in Q_h(K_i^+). \end{align}\]

Then, the global multiscale velocity space is defined by taking the linear span of all these localized velocity basis functions: \[V_{{ms}} = \operatorname{span} \left\{ \psi_{j,{ms}}^i : i = 1, \dots, N, \;j = 1, \dots, J_i \right\}.\]

3.3 The adaptive method and discretizations↩︎

From the previous discussion, the multiscale spaces for pressure and velocity have been constructed as \(Q_{ ms}=\operatorname{span}\{p_1,\ldots,p_{N_{\rm ms}}\}\), \(V_{ ms}=\operatorname{span}\{\psi_1,\ldots,\psi_{N_{\rm ms}}\}\), where \(N_{ ms}\) denotes the number of multiscale basis functions in the final space. In the standard implementation, these multiscale spaces are constructed at the initial time and then kept fixed throughout the whole simulation. However, since the effective permeability \(\kappa_n=\lambda_t(S_w^n)K\) depends on the saturation, the coefficient entering the pressure and velocity equations changes as the saturation updates. Therefore, the multiscale basis functions constructed from the initial coefficient may become outdated at later time levels.

To reduce this coefficient-outdated error, we introduce an adaptive update strategy for the multiscale spaces. The main idea is to monitor the variation of the coefficient field by an error indicator and update the CEM multiscale spaces once this indicator exceeds a prescribed tolerance. More precisely, when the update criterion is satisfied, the multiscale spaces are reconstructed by replacing the initial coefficient used in the bilinear forms in Subsections 3.1 and 3.2 with the current coefficient \(\kappa_n^H\).

We first introduce an error indicator to measure the variation of the coefficient field. Let \(i\) denote the last time level at which the multiscale space was updated, and let \(n\) be the current time level. We define \[\eta_n = \left\| \left(\kappa_i^H\right)^{-\frac{1}{2}} - \left(\kappa_n^H\right)^{-\frac{1}{2}} \right\|_{L^2(\Omega)}\] Then the adaptive multiscale solution can be obtained by the following steps.
For given \(S_w^{H,n}\) at the time step \(t_n\) and most recent updated space \(Q_{ms}^i\) and \(V_{ms}^i\), we seek the multiscale solutions of the P-IMPES scheme at the time step \(t_{n+1}\) as follows:

Step 1. First, find the capillary velocity \(\xi_c^{H,n+1}\) by solving \[\begin{align} a_n^H( \xi_c^{H,n+1},w)=b(w,p_c(S_w^{H,n}))-b_D(w,p_n^B-p_w^B) && \forall w\in V_{ms}^i\label{3461} \end{align}\tag{5}\]

Step 2. Then seek \(u_t^{H,n+1} \in V_{ms}^i\) and \(p_w^{H,n+1} \in Q_{ms}^i\) satisfying \[\begin{align} a_n^H(u_t^{H,n+1},w) - b(w,p_w^{H,n+1}) &= a_n^H(f_n(S_w^{H,n})\xi_c^{H,n+1} ,w) - b_D(w,p_w^B) \quad \forall w \in V_{ms}^i\tag{6} \\ \sum_\alpha \beta_\alpha({u}_t^{H,n+1}, q; S_w^{H,n})& = (F_t,q) \qquad \qquad \forall q \in Q_{ms}^i,\tag{7} \end{align}\] Moreover, to guarantee the uniqueness of the solution, we add an additional constraint as \(\int_{\Omega} p_w^{H,n+1} = 0\) in step 2.

Step 3. Update the two phase saturation \(S_w^{H,n+1}\) and \(S_n^{H,n+1}\) in every fine element \(K^h \in \mathcal{T}_h\) by \[\begin{align} (\phi \frac{S_w^{H,n+1} - S_w^{H,n}}{\Delta t},\mathbf{1}_{K^h})& = -\beta_w({u}_t^{H,n+1}, \mathbf{1}_{K^h}; S_w^{H,n}) +\beta_c(\xi_c^{H,n+1},\mathbf{1}_{K^h}; S_w^{H,n}) +( F_w , \mathbf{1}_{K^h})\label{3464} \\ S_n^{H,n+1}&=1-S_w^{H,n+1} \end{align}\tag{8}\]

Step 4. Then update the effective permeability by \(\kappa_{n+1}^H=\lambda_t(S_w^{H,n+1})K\), and we compute the error indicator \[\eta_{n+1} = \left\| \left(\kappa_i^H\right)^{-\frac{1}{2}} - \left(\kappa_{n+1}^H\right)^{-\frac{1}{2}} \right\|_{L^2(\Omega)} ,\] For a given prescribed tolerance \(\varepsilon_{\rm }>0\).

If \(\eta_{n+1}>\varepsilon\), we regenerate the multiscale spaces in Subsections 3.1 and 3.2 using the coefficient \(\kappa_{n+1}^H\). The newly constructed spaces are denoted by \(Q_{ ms}^{n+1}\), \(V_{ ms}^{n+1}\). Then we set \(i=n+1\) .

Otherwise, if \(\eta_{n+1}\leq \varepsilon\), the current multiscale spaces are reused, namely \(Q_{ ms}^{i}\), \(V_{ ms}^{i}\), which are kept unchanged for the next time step.
0◻

Although our method is formulated by focusing on the wetting phase (\(\alpha = w\)), in fact, this P-IMPES scheme satisfies an unbiased property. It means we can get exactly the same solution using \(\alpha = n\). We will theoretically prove this property in the next section and numerically check it in the numerical experiment.

We note that the velocity obtained from 6 –7 is conservative only with respect to the multiscale pressure space. To avoid unnecessary local solves, we postprocess only those coarse elements containing a nonzero fine-grid mass residual. Let \(\widehat{\mathbf{u}}_t^{H,n+1}\) denote the multiscale velocity before postprocessing. For each fine element \(K^h\in\mathcal{T}_h\), define \(R_{K^h}^{n+1} := (F_t,\mathbf{1}_{K^h}) - \sum_{\alpha} \beta_\alpha \left( \widehat{\mathbf{u}}_t^{H,n+1}, \mathbf{1}_{K^h};S_w^{H,n} \right),\) and introduce the marked set \(\mathcal{T}_{H,\mathrm{pt}}^{n+1} = \left\{ K\in\mathcal{T}_H: \exists K^h\in K,R_{K^h}^{n+1}\neq0 \right\}.\) For each \(K\in\mathcal{T}_{H,\mathrm{pt}}^{n+1}\), we seek \(\mathbf{u}_{t,p,K}^{H,n+1}\in V_h(K)\) and \(p_{w,p,K}^{H,n+1}\in Q_h(K)/\mathbb{R}\) satisfying \[\begin{align} a_{n,K}^H(\mathbf{u}_{t,p,K}^{H,n+1},\mathbf{v}) -b_K(\mathbf{v},p_{w,p,K}^{H,n+1}) &= a_{n,K}^H \left( f_n(S_w^{H,n})\boldsymbol{\xi}_c^{H,n+1}, \mathbf{v} \right), &&\forall\mathbf{v}\in V_{h,0}(K),\\ \sum_\alpha \beta_{\alpha,K} \left( \mathbf{u}_{t,p,K}^{H,n+1}, q;S_w^{H,n} \right) &=(F_t,q)_K, &&\forall q\in Q_h(K), \end{align}\] subject to \[\mathbf{u}_{t,p,K}^{H,n+1}\cdot\mathbf{n}_K = \widehat{\mathbf{u}}_t^{H,n+1}\cdot\mathbf{n}_K \quad\text{on }\partial K, \qquad (p_{w,p,K}^{H,n+1},\mathbf{1}_K)=0.\] The upwind directions are kept fixed during each local solve. On the unmarked coarse elements, no local problem is solved, and we set \(\mathbf{u}_{t,p}^{H,n+1}|_K = \widehat{\mathbf{u}}_t^{H,n+1}|_K,\text{ for }K\notin\mathcal{T}_{H,\mathrm{pt}}^{n+1}.\) Since the coarse-element constants belong to \(Q_{\mathrm{ms}}^i\), equation 7 implies \(\sum_{K^h\subset K}R_{K^h}^{n+1}=0,\) which gives the compatibility condition for every marked local problem. Moreover, the normal trace is unchanged on each coarse-element boundary, and hence \(\mathbf{u}_{t,p}^{H,n+1}\in V_h\). We introduce the cumulative marked set \(\mathcal{A}_{H,\mathrm{pt}}^{i,n+1} = \bigcup_{k=i+1}^{n+1} \mathcal{T}_{H,\mathrm{pt}}^{k},\) where \(i\) denotes the most recent time level at which the CEM space was updated. We then define the selectively enriched velocity space by \(\widetilde{V}_{\mathrm{ms}}^{i,n+1} = V_{\mathrm{ms}}^i + \bigoplus_{K\in\mathcal{A}_{H,\mathrm{pt}}^{i,n+1}} V_{h,0}(K) \subset V_h.\) The postprocessed pair \((\mathbf{u}_{t,p}^{H,n+1},p_{w,p}^{H,n+1}) \in \widetilde{V}_{\mathrm{ms}}^{i,n+1} \times Q_h/\mathbb{R}\) is characterized by \[\begin{align} a_n^H \left( \mathbf{u}_{t,p}^{H,n+1},\mathbf{w} \right) - b \left( \mathbf{w},p_{w,p}^{H,n+1} \right) &= a_n^H \left( f_n(S_w^{H,n})\boldsymbol{\xi}_c^{H,n+1}, \mathbf{w} \right) - b_D(\mathbf{w},p_w^B), && \forall\mathbf{w}\in \widetilde{V}_{\mathrm{ms}}^{i,n+1}, \tag{9} \\ \sum_\alpha \beta_\alpha \left( \mathbf{u}_{t,p}^{H,n+1}, q;S_w^{H,n} \right) &= (F_t,q), && \forall q\in Q_h. \tag{10} \end{align}\] We therefore set \(\mathbf{u}_t^{H,n+1} = \mathbf{u}_{t,p}^{H,n+1}, p_w^{H,n+1} = p_{w,p}^{H,n+1}\) in the subsequent analysis and saturation update.

Next, we will discuss the full discretization of the problem. Let the fine-scale spaces for pressure and velocity be spanned by \(Q_h = \text{span}\{q_1, \ldots, q_{N_f}\}\) and \(V_h = \text{span}\{\phi_1, \ldots, \phi_{N_e}\}\), and the corresponding multiscale spaces be \(Q_{ms}^i = \text{span}\{p_1, \ldots, p_{N_{ms}}\}\) and \(V_{ms}^i = \text{span}\{\psi_1, \ldots, \psi_{N_{ms}}\}\). To relate these spaces algebraically, we construct the mapping matrices. Since for every \(\psi_i \in V_{ms} \subset V_h\), we can express \(\psi_i\) by the linear span of fine scale basis as \(\psi_i=\sum\limits_{j=1}^{N_e} {v}_j \phi_j\). Collect all these representative vectors \(\mathbf{v}_i=(v_j)_{N_e \times 1}\) and assemble them into a mapping matrix denoted by \(V_i=[\mathbf{v}_1, \ldots, \mathbf{v}_{N_{ms}}] \in \mathbf{R}^{N_e \times N_{ms}}\). Similarly, we have the mapping matrix for the multiscale pressure space \(P_i=[\mathbf{p}_1, \ldots,\mathbf{p}_{N_{ms}}]\). Recall that \(i\) is the last time level when the most recent multiscale space was updated.

First, we give the following definitions of the matrices: \[\begin{align} &A^n=\left(a_n^H(\phi_i,\phi_j) \right)_{N_{e} \times N_{e}}, \quad A^{H,n}=V_i^TA^nV_i\\ &B^n=\left( b(\phi_i,q_j) \right)_{N_{e} \times N_{f}}, \qquad B^{H,n}=V^T_iB^nP_i \\ & B_{\alpha}^n=\left( \beta_\alpha(\phi_i, q_j; S_w^{H,n})\right)_{N_{e} \times N_f}, \quad B_{\alpha}^{H,n}=V^T_iB_{\alpha}^{n}P_i \\ &B_c^n=\left(\beta_c(\phi_i, q_j; S_w^{H,n})\right)_{N_{e} \times N_{f}}, \quad B_{c}^{H,n}=V^T_iB_{c}^{n}P_i \\ &A_c^n=\left(a_n^H(f_n(S_w^{H,n})\phi_i,\phi_j) \right)_{N_{e} \times N_{e}}, \quad A_c^{H,n}=V^T_iA_c^nV_i \\ &D =\left( b_D(\phi_i,1) \right)_{N_{e} \times1} ,\qquad \qquad D^H=V^T_iD \\ & \mathbf{F}_{\alpha}=(F_{\alpha},q_i)_{N_{f} \times1} , \quad \mathbf{F}_t=\mathbf{F}_w +\mathbf{F}_n , \quad \Pi_i\mathbf{F}_{\alpha}=P_i^T \mathbf{F}_{\alpha} \end{align}\]

(Matrix–Vector Formulation) For given \(\mathbf{S}_w^{H,n}\in \mathbb{R}^{N_f}\) at the time step \(t_n\). We seek the multiscale solutions of the P-IMPES scheme at the time step \(t_{n+1}\) as follows:
Step 1. First, discretize the capillary velocity \(\xi_c^{H,n+1}=\sum\limits_{i=1}^{N_{ms}} x_i^{n+1} \psi_i\) and obtain \(x_c^{n+1}=(x_i^{n+1})_{N_{ms}\times 1}\) by solving \[\begin{align} A^{H,n}x_c^{n+1}=V^T_iB^np_c(\mathbf{S}_w^{H,n})-(p_n^B-p_w^B)D^H \end{align}\] Step 2. We further seek \(\mathbf{u}_t^{H,n+1} \in \mathbb{R}^{N_{ms}}\) , \(\mathbf{p}_w^{H,n+1}\in \mathbb{R}^{N_{ms}}\) by \[\begin{bmatrix} {A}^{H,n} & -B^{H,n} \\ -(B_{w}^{H,n}+B_n^{H,n})^T & 0 \end{bmatrix} \begin{bmatrix} \mathbf{u}_t^{H,n+1} \\ \mathbf{p}_w^{H,n+1} \end{bmatrix} = \begin{bmatrix} {A}_c^{H,n} x_c^{n+1} - p_w^B{D^H} \\ -\Pi_i\mathbf{F}_t \end{bmatrix}.\] Then apply the local postprocessing procedure in every coarse element \(K\in \mathcal{A}_{H,\mathrm{pt}}^{i,n+1}\) and assemble the postprocessed fine-grid velocity \(\mathbf{u}_{t,p}^{H,n+1}\).
Step 3. Finally, update the two phase saturation \(\mathbf{S}_w^{H,n+1}\) and \(\mathbf{S}_n^{H,n+1}\) by \[\begin{align} \phi \frac{\mathbf{S}_w^{H,n+1} - \mathbf{S}_w^{H,n}}{\Delta t}& = -(B_w^{n})^T\mathbf{u}_{t,p}^{H,n+1} +(B_c^{n})^TV_ix_c^{n+1}+\mathbf{F}_w \\ \mathbf{S}_n^{H,n+1}&=1-\mathbf{S}_w^{H,n+1} \end{align}\] Step 4. Adaptive update of the multiscale spaces. After obtaining \(\mathbf{S}_w^{H,n+1}\), set \(\kappa_{n+1}^H=\lambda_t(\mathbf{S}_w^{H,n+1})K\) and compute \(\eta_{n+1} = \| \left(\kappa_i^H\right)^{-\frac{1}{2}} - \left(\kappa_{n+1}^H\right)^{-\frac{1}{2}} \|_2 .\)

If \(\eta_{n+1}>\varepsilon\), then regenerate the CEM multiscale spaces with the current coefficient \(\kappa_{n+1}^H\), and update \[Q_{ms}^{i}=Q_{ ms}^{n+1},\qquad V_{ms}^{i}=V_{ ms}^{n+1},\qquad P_i=P_{n+1},\qquad V_i=V_{n+1},\qquad i=n+1.\] Otherwise, the current spaces and mapping matrices are kept unchanged.

We remark that Step 2 is solved together with the condition \(\int_{\Omega} p = 0\), so we will add Lagrange multiplier in (3.9) to satisfy the constraint.

4 Analysis↩︎

In this section, we present the theoretical analysis of the proposed method in two main parts. First, following the analytical framework in [12], we establish several key properties of the scheme. Subsequently, we derive an essential relationship between the saturation error and the velocity error.

4.1 Local conservation of mass and unbiased property↩︎

In this subsection, we prove two physical properties of the proposed multiscale P-IMPES scheme: local mass conservation and the unbiased property. Although only the wetting-phase saturation is updated explicitly, the non-wetting-phase saturation is determined by the saturation constraint and also satisfies a conservative discrete equation. This follows from the fine-grid conservation of the total velocity and the wetting-phase conservation equation. We also show that the same discrete solution is obtained whether the saturation equation is written for the wetting phase or for the non-wetting phase.

Next We will prove the result in two parts.

(i) Local conservation of mass. The scheme updates only the wetting phase via \[\begin{align} \left(\phi\frac{S_w^{H,n+1}-S_w^{H,n}}{t_{n+1}-t_n},q\right) +\beta_w(\mathbf{u}_t^{H,n+1},q;S_w^{H,n}) &=(F_w,q)+\beta_c(\xi_c^{H,n+1},q;S_w^{H,n}), \\ S_n^{H,n+1}&=1-S_w^{H,n+1} \end{align}\] The two phase mass conservation means that both phases satisfy: \[\begin{align} \left(\phi\frac{S_w^{H,n+1}-S_w^{H,n}}{t_{n+1}-t_n},q\right) +\beta_w(\mathbf{u}_t^{H,n+1},q;S_w^{H,n}) &=(F_w,q)+\beta_c(\xi_c^{H,n+1},q;S_w^{H,n}), \tag{11}\\ \left(\phi\frac{S_n^{H,n+1}-S_n^{H,n}}{t_{n+1}-t_n},q\right) +\beta_n(\mathbf{u}_t^{H,n+1},q;S_w^{H,n}) &=(F_n,q)-\beta_c(\xi_c^{H,n+1},q;S_w^{H,n}). \tag{12} \end{align}\] Since (4.1) and (4.2) coincide, it suffices to show that \(S_n^{H,n+1}\) from \(S_n^{H,n+1}=1-S_w^{H,n+1}\) agrees with the solution of (4.3). On the one hand, adding (4.2) and (4.3), using 10 we have \[\sum_\alpha\left(\phi\frac{S_\alpha^{H,n+1}-S_\alpha^{H,n}}{t_{n+1}-t_n},q\right)=0.\] Choosing \(q=\mathbf{1}_{K^h}\) gives \(\sum_\alpha(S_\alpha^{H,n+1}-S_\alpha^{H,n})=0\) on each \(K^h\). Combined with the initial condition \(S_w^{H,0}+S_n^{H,0}=1\), this derives \(S_w^{H,n+1}+S_n^{H,n+1}=1\) for all \(n\). On the other hand, subtracting (4.1) by (3.7) and replacing \(S_w^{H}=1-S_n^{H}\) in the first term left, we can obtain (4.3) directly.

(ii) Unbiased property for velocity. From part (i), \(S_w^{H,n+1}\) and \(S_n^{H,n+1}\) are unbiased. The total velocity \(\mathbf{u}_t^{H,n+1}\) and the wetting pressure \(p_w^{H,n+1}\) are uniquely determined by 6 and 7 ; \(\xi_c^{H,n+1}\) is uniquely solved by 5 . The phase velocities \(\mathbf{u}_\alpha^{H,n+1}\) are unbiased via \[\begin{align} \mathbf{u}_w^{H,n+1}&=f_w(S_w^{H,n})\mathbf{u}_t^{H,n+1}-f_w(S_w^{H,n})f_n(S_w^{H,n})\xi_c^{H,n+1},\\ \mathbf{u}_n^{H,n+1}&=f_n(S_w^{H,n})\mathbf{u}_t^{H,n+1}+f_w(S_w^{H,n})f_n(S_w^{H,n})\xi_c^{H,n+1}, \end{align}\]

4.2 Bounds-preserving property for both phases↩︎

This part is motivated by the arguments in [12]. We first give two auxiliary lemmas that will be used to establish the bounds-preserving property.

This lemma shows that each fractional flow function can be bounded above by a linear function of the corresponding saturation.

The above two lemmas have been proved in [12]. Based on these results, we now state the theorem that guarantees the bounds-preserving property.

Remark. By the proof of Theorem 4.1, the stability restriction on the time step can be written in a phase-flux form as \[\max_{\alpha=w,n} E_\alpha^{n+1}\frac{\Delta t}{h}+C\Delta t<\phi,\] where \(0<\phi<1\), \(C\) depends on the sink term, and \[E_\alpha^{n+1} = \max_{K^h\in\mathcal{T}_h} \sum_{F\subset\partial K^h} (f_\alpha(S_\alpha^{*,H,n}) \mathbf{v}_\alpha^{H,n+1}\cdot\mathbf{n} )^+ , \qquad \alpha=w,n .\] Here \(a^+=\max\{a,0\}\), and \(\mathbf{v}_\alpha^{H,n+1}\) denotes the effective phase velocity appearing in the saturation equation. In the numerical experiments, we use the following CFL condition, which includes both the convective flux induced by the total velocity and the capillary diffusion flux induced by \(\boldsymbol{\xi}_c^{H,n+1}\): \[\begin{align} \theta_{cfl} &= \max_{n} \frac{\Delta t}{\phi h} \Bigg( \|\mathbf{u}_t^{H,n+1}\|_{L^\infty(\Omega)} \|f_w'\|_{L^\infty(0,1)} + \|{\xi}_c^{H,n+1}\|_{L^\infty(\Omega)} \|(f_w f_n)'\|_{L^\infty(0,1)} \Bigg) < 1 . \end{align}\] The time step size \(\Delta t\) is chosen such that \(\theta_{cfl}<1\).

4.3 Error estimate↩︎

In this subsection, we establish the error estimates. To obtain the numerical errors induced by the multiscale velocity approximation, we restrict our analysis to the zero-capillary form, namely the advection-dominated Buckley-Leverett flow. Incorporating capillary pressure leads to a degenerate parabolic–hyperbolic problem and requires additional estimates for the capillary velocity, which are left for future work. Under this condition, we set zero Dirichlet and zero Neumann boundary condition, then 2 , 3 can be reduced to: \[\begin{align} a_n(\mathbf{u}_t^{h, n+1},v) - b(v,p_w^{h, n+1} ) &= 0, && \forall v\in V_h \tag{13} \\ b(\mathbf{u}_t^{h, n+1}, q) &= (F_t, q) && \forall q\in Q_h \tag{14} \end{align}\]

Similarly for 9 , 10 , the postprocessed multiscale solution satisfies: \[\begin{align} a_n^H(\mathbf{u}_t^{H, n+1},w) - b(w,p_w^{H, n+1} ) &= 0 && \forall w\in \widetilde{V}_{ms}^{i,n+1}\tag{15} \\ b(\mathbf{u}_t^{H, n+1}, q) &= (F_t, q) && \forall q\in Q_{h}\tag{16} \end{align}\]

Besides, define a local version of the norms \(\|\cdot\|_{a_i(D)}\), \(\|\cdot\|_{V_i(D)}\) and \(\|\cdot\|_{s_i(D)}\) at corresponding time \(t=t_i\). \[\begin{align} \|v\|_{V_i(D)}^2 &= \int_D \tilde{\kappa}^{-1}_i|\nabla \cdot v|^2 + \int_D \kappa^{-1}_i|v|^2, \quad \|q\|_{s_i(D)}^2 = \int_D \tilde{\kappa}_i|q|^2\\ \|v\|_{a_i(D)}^2 &= \int_D \kappa_i^{-1}|v|^2 , \quad \|v\|_{a_i^H(D)}^2 = \int_D (\kappa_i^H)^{-1}|v|^2. \end{align}\]

Next, we will introduce some assumptions used in the error analysis.

4.3.0.1 Same-direction assumption.

For every fine edge \(F\subset \partial K_h\), \[\left(u_t^{h,n+1}\cdot \mathbf{n}\right) \left(u_t^{H,n+1}\cdot \mathbf{n}\right)\big|_F \ge 0.\] Thus the fine-scale and multiscale velocities determine the same upwind direction on each fine edge.

4.3.0.2 Front-layer decomposition.

For each time level \(t_n\), we decompose the fine mesh into a smooth region and a front layer, \[\mathcal{T}_h=\mathcal{G}_h^n\cup\mathcal{B}_h^n, \qquad \mathcal{G}_h^n\cap\mathcal{B}_h^n=\emptyset.\] The set \(\mathcal{B}_h^n\) contains the cells close to the saturation fronts, shocks, and possible mismatch strips between the fine-scale and multiscale fronts. We assume that \[|\Omega_{\mathcal{B}}^n| \le \delta_h , \qquad 0\le n\le N.\] On the smooth region, the two saturation profiles are locally smooth in the following discrete sense: for neighboring fine cells \(K^h,K^{h,1}\in\mathcal{G}_h^n\) which share an edge \(F\), we define \[\begin{align} \zeta_h = \max_{0\le n\le N} \max_{\substack{ K^h,K^{h,1}\in \mathcal{G}_h^n\\ F=K^h\cap K^{h,1} }} \Big[ \omega_{n,F}^{\frac{1}{2}} |S_{w,K^h}^{h,n}-S_{w,K^{h,1}}^{h,n}| + (\omega_{n,F}^H)^{\frac{1}{2}} |S_{w,K^h}^{H,n}-S_{w,K^{h,1}}^{H,n}| \Big], \end{align} \label{44611}\tag{17}\] where \(\omega_{n,F}=\max(\kappa_n|_{K^h},\kappa_n|_{K^{h,1}}), \quad \omega_{n,F}^H=\max(\kappa_n^H|_{K^h},\kappa_n^H|_{K^{h,1}}).\) On the smooth region, we assume that \[\frac{\zeta_h}{h}=O(1) \qquad \text{as } h\to0.\] Moreover, the front layer is enlarged by a fixed number of fine-cell layers so that the following stencil closure property holds: if \(K^h\in\mathcal{G}_h^{n+1}\), then every fine cell involved in the upwind update of \(K^h\) at time \(t_n\) belongs to \(\mathcal{G}_h^n\).

For a velocity \(u\), define the upwind update operator \(\mathcal{L}_u\) by \[\begin{align} (\mathcal{L}_u S)_{K^h} = S_{K^h} - \frac{\Delta t}{\phi |K^h|} \sum_{F\subset\partial K^h} |F|\left(u_{K^{h,F}}^{+} f_w(S_{K^{h}})+ u_{K^{h,F}}^{-} f_w(S_{K^{h,i}}) \right) +\frac{\Delta t}{\phi |K^h|}(F_w,1_{K^h}), \label{44612} \end{align}\tag{18}\] Here \(u^{+}=\max\{u,0\}\) denotes the outgoing flux, \(u^{-}=\min\{u,0\}\) denotes the incoming flux.

Remark that to obtain the \(O(H)\) convergence in \(\|\mathbf{u}_t^{h,i+1} - \mathbf{u}_t^{H,i+1} \|_{a_i^H}\), we need to choose the size of the oversampling domain \(l=O(log(\mathbf{K} /{H^2}))\). For a detailed derivation of these contents, refer to [39].

5 Numerical experiments↩︎

In this section, we present several numerical experiments to evaluate the performance of the proposed multiscale method. In particular, we focus on its mass conservation property, saturation approximation, and error convergence behavior in heterogeneous and high-contrast porous media. In the following experiments, the computational domain is \(\Omega = [0, 1]^2\), the capillary pressure is modeled by the Hoteit–Firoozabadi formula \(p_c(S_w)=-\frac{B_c}{\sqrt{\mathbf{K}}}\log \bar S_w\), where \(B_c \geq 0\) is a capillary parameter, and \(\bar S_w\) is the effective saturation defined as \(\bar S_w=\frac{S_w-S_{rw}}{1-S_{rw}-S_{rn}}\), here we choose \(S_{rw} = S_{rn} = 10^{-6}\). The porosity is fixed as \(\phi=0.2\) throughout all numerical tests. Moreover, the relative permeabilities of the wetting and non-wetting phases are given by \(k_{rw}=\bar S_w^2,\quad k_{rn}=(1-\bar S_w)^2\), the viscosity for two phases is set as \(\mu_w=1\) and \(\mu_n=5\). We use the following relative errors to evaluate the performance of the proposed method: \[e_{u} = \frac{\|\mathbf{u}_h-\mathbf{u}_{\mathrm{ms}}\|_{\kappa^{-1}}}{\|\mathbf{u}_h\|_{\kappa^{-1}}}, \quad e_{S} = \frac{\|S_h-S_{ms}\|_{L^2(\Omega)}}{\|S_h\|_{L^2(\Omega)}},\] where \(\|u\|_{\kappa^{-1}}^2 = \int_{\Omega} \mathbf{K}^{-1}|u|^2\), \((\mathbf{u}_h,S_h)\) is the fine-scale solution as reference solution, and \((\mathbf{u}_{ms},S_{ms})\) is the multiscale solution from our scheme.
Example 5.1. (Two phase conservation) In this example we will check the two phase conservation and unbiased property using the multiscale basis obtained in a local region of the high-contrast medium \(\kappa_1\). As mentioned in the subsection 4.1, there are two ways 4.1 and 11 12 to obtain the non-wetting phase saturation, denoted by \(S_{n,1}\) and \(S_{n,2}\) respectively. If \(S_{n,1}=S_{n,2}\), the local mass conservation is attained for both two phases. We use a high contrast permeability field \(\kappa_1\) to check the property. See Figure 2 for the illustration of the source term \(F_t\) and the permeability field \(\kappa_1\). We choose only \(1\) basis in every coarse element, it is sufficient to conclude that multiscale solutions can well satisfy this conservation property since solutions computed with more multiscale bases will obtain higher accuracy. We set the fine mesh size to \(100 \times 100\), the tolerance \(\varepsilon=0.1\) in adaptive algorithm and the capillary parameter \(B_c=0\). The source/sink term function \(F_t\) is set to be \[F_t(\boldsymbol{x})= \begin{cases} 1, & \boldsymbol{x}\in [0,0.01]^2,\\ -1, & \boldsymbol{x}\in [0.99,1]^2,\\ 0, & \text{otherwise}, \end{cases}\]

Figure 3.1-3.3 shows the non-wetting phase saturation computed by the two different formulations at \(T=1,2,3\), with their absolute differences. It can be seen that \(S_{n,1}\) and \(S_{n,2}\) produce nearly the same saturation distributions during the whole simulation. The difference \(|S_{n,1}-S_{n,2}|\) remains at the order of about \(10^{-13}\) throughout the simulation, indicating that the two formulations produce almost the same outcomes.
Example 5.2. (Capillary pressure effect) In this example, we investigate the influence of the capillary pressure on the two-phase flow simulation. Recall that the capillary pressure is modeled by the Hoteit-Firoozabadi formula \(p_c(S_w) = -\frac{B_c}{\sqrt{\mathbf{K}}}\log \bar{S}_w\), the parameter \(B_c \ge 0\) controls the strength of the capillary effect. We test three representative cases: \(B_c = 0\), \(B_c = 10^{-3}\) and \(B_c = 10^{-2}\) at different time \(T=1,2,3\). For other parameters, we keep the same as in Example 5.1.

Figure 4.1-4.3 illustrate the influence of capillary pressure on the saturation evolution. When \(B_c=0\), the displacement is mainly dominated by the advection and the heterogeneous permeability field, so the wetting phase mainly moves through the high-permeability channels and the fronts remain relatively sharp. As \(B_c\) increases to \(10^{-3}\) and \(10^{-2}\), the capillary pressure introduces an additional diffusion-like effect, which smooths the saturation profile and broadens the transition region near the front. In particular, for \(B_c=10^{-2}\), the channelized displacement becomes less pronounced and the wetting phase spreads more uniformly into the surrounding region. These results show that stronger capillary pressure regularizes the saturation distribution and weakens sharp saturation discontinuities.

Example 5.3. In this example, we study the influence of the adaptive update tolerance \(\varepsilon\), coarse element size \(H\) and the number of multiscale basis functions on the velocity and saturation errors. We consider two representative high-contrast permeability fields. The first one is a high-permeability channel field, denoted by \(\kappa_2\), the second one is taken from the top layer of the SPE10 permeability field and is denoted by \(\kappa_3\). These two permeability fields are shown in the figure below. The capillary pressure parameter is set to be \(B_c=0\) and fine-mesh size is still \(h=0.01\). For other parameters, we keep the same as before.

First, we investigate the influence of the coarse mesh size \(H\) and the number of multiscale basis functions on the velocity and pressure errors at the initial time. The results are shown in Table 1 and Table 2, respectively.

4.5pt

4.5pt

Next, we investigate the evolution of the saturation error for the two permeability fields \(\kappa_2\) and \(\kappa_3\). In particular, we study how the saturation error \(e_S\) depends on the adaptive update tolerance \(\varepsilon\). To balance computational time and accuracy, we select the parameters \(H\) and \(n\) according to the Table 1 and Table 2, where the corresponding velocity error \(e_u\) is sufficiently small. These results are performed at different final times and for different permeability fields.

Next, we use the SPE10 permeability field \(\kappa_3\) to test the influence of different adaptive tolerances \(\varepsilon\) at different final times.

The numerical results show that the proposed adaptive multiscale method can provide accurate approximations for both permeability fields \(\kappa_2\) and \(\kappa_3\). In general, smaller coarse mesh size \(H\) and more multiscale basis functions lead to smaller velocity and pressure errors. For the saturation equation, the numerical results show that the saturation error \(e_S\) decreases as the adaptive tolerance \(\varepsilon\) becomes smaller. This indicates that updating of the multiscale space appropriately can improve the accuracy of the saturation approximation.

6 Conclusion↩︎

In this paper, we developed an adaptive mixed CEM-GMsFEM for incompressible and immiscible two-phase flow in heterogeneous porous media. The method combines with a physics-preserving IMPES scheme, and the multiscale spaces are adaptively regenerated according to the variation of the effective permeability induced by the saturation evolution. We proved the local mass conservation for both phases, the unbiased property of the formulation, and the bounds-preserving property under a suitable CFL condition. For the zero-capillary case, velocity and saturation error estimates were derived, showing the influence of the adaptive tolerance, coarse mesh size, spectral approximation error, and front layer. Numerical results confirm the theoretical properties and show that the capillary pressure enhances the diffusive nature of the advection-dominated motion. Moreover, smaller adaptive tolerance generally leads to smaller saturation errors, while suitable coarse mesh size and number of multiscale basis functions provide accurate velocity and pressure approximations with reduced computational cost.

Acknowledgement↩︎

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

References↩︎