An HMM--ELLAM scheme on generic polygonal meshes for miscible incompressible flows in porous media

TL;DR

This paper presents a novel numerical scheme combining characteristic and finite volume methods for simulating miscible incompressible flows in porous media on complex polygonal meshes, improving accuracy and conservation.

Contribution

It introduces a Darcy velocity reconstruction with local consistency, addresses implementation challenges with distorted cells, and enhances treatment near injection wells for better fluid conservation.

Findings

Applicable to generic polygonal meshes, including non-conforming ones.

Ensures local consistency in Darcy velocity reconstruction.

Improves fluid conservation near injection wells.

Abstract

We design a numerical approximation of a system of partial differential equations modelling the miscible displacement of a fluid by another in a porous medium. The advective part of the system is discretised using a characteristic method, and the diffusive parts by a finite volume method. The scheme is applicable on generic (possibly non-conforming) meshes as encountered in applications. The main features of our work are the reconstruction of a Darcy velocity, from the discrete pressure fluxes, that enjoys a local consistency property, an analysis of implementation issues faced when tracking, via the characteristic method, distorted cells, and a new treatment of cells near the injection well that accounts better for the conservativity of the injected fluid.

Tables

Table 1. Table 1. Relative errors in velocity reconstruction.

Mesh	$\frac{\Vert V-V_{\mathrm{KR}}\Vert }{\Vert V\Vert }$	$\frac{\Vert V-V_C\Vert }{\Vert V\Vert }$
Cartesian	2.9038e-15	5.9529e-16
Hexahedral	0.0371	5.8098e-15
Non-conforming	0.0348	6.7103e-15
Kershaw	0.3151	4.6738e-15

Table 2. Table 2. Regularity parameter of the meshes and nb of points to approximate the trace-back cells.

Mesh	$m_{ℳ\text{reg}}$	${\log }_2\left(m_{ℳ\text{reg}}\right)$	points per edge
Cartesian	2	1	1
Hexahedral	5.4772	2.4534	3
Non-conforming	2.7619	1.4657	2
Kershaw	32.0274	5.0012	6

Equations

$$\nabla \cdot \mathbf{u}$$

$$\mathbf{u}$$

$$\varphi \frac{\partial c}{\partial t}+\nabla \cdot (\mathbf{u}c-\mathbf{D}(\mathbf{x},\mathbf{u})\nabla c)=q^{+}-c{q}^{-}:=q_c\text{mbox}on\Omega \times [0,T]$$

$$\mathbf{D}(\mathbf{x},\mathbf{u})=\varphi (\mathbf{x})\left[d_m\mathbf{I}+d_l|\mathbf{u}|\mathcal{P}(\mathbf{u})+d_t|\mathbf{u}|\left(\mathbf{I}-\mathcal{P}(\mathbf{u})\right)\right]\text{mbox}with\mathcal{P}(\mathbf{u})={\left(\frac{u_i{u}_j}{|\mathbf{u}{|}^2}\right)}_{i,j}.$$

$$\mathbf{u}\cdot\mathbf{n}=(\mathbf{D}\nabla c)\cdot\mathbf{n}=0\mbox{ on }\partial\Omega\times[0,T]\,,\qquad c(\cdot,0)=0\mbox{ in $\Omega$}.$$

$${\int }_{\Omega }p(\mathbf{x},t)\mathrm{d}\mathbf{x}=0\forall t\in [0,T]$$

$${\int }_{\Omega }\frac{\mathbf{K}}{\mu (c)}\nabla p\cdot \nabla \xi ={\int }_{\Omega }q(t^{(n+1)})\xi ,\forall \xi \in H^1(\Omega ).$$

$$X_{\mathcal{T}}:=\{w=((w_K{)}_{K\in \mathcal{M}},(w_{\sigma }{)}_{\sigma \in {\mathcal{E}}_K})\}.$$

$$\forall w\in X_{\mathcal{T}},\forall \mathbf{x}\in T_{K,\sigma },$$

$${\nabla }_{H}w(\mathbf{x})={\overline{\nabla }}_{K}w+\frac{\sqrt{2}}{d_{K,\sigma }}[w_{\sigma }-w_K-{\overline{\nabla }}_{K}w\cdot ({\overline{\mathbf{x}}}_{\sigma }-{\mathbf{x}}_K)]{\mathbf{n}}_{K,\sigma },$$

$$\text{mbox}where{\overline{\nabla }}_{K}w=|K{|}^{-1}\sum _{\sigma \in {\mathcal{E}}_K}|\sigma |w_{\sigma }{\mathbf{n}}_{K,\sigma }.$$

$$\forall K\in \mathcal{M},\forall v\in X_{\mathcal{T}},$$

$$\sum _{\sigma \in {\mathcal{E}}_K}{F}_{K,\sigma }^{(n+1)}(v_K-v_{\sigma })={\int }_K\frac{{\mathbf{K}}_K}{\mu (c_K^{(n)})}{\nabla }_H{p}^{(n+1)}(\mathbf{x})\cdot {\nabla }_H\xi (\mathbf{x})d\mathbf{x}$$

$$\sum_{\sigma\in\mathcal{E}_{K}}F^{(n+1)}_{K,\sigma}=\int_{K}q(t^{(n+1)})\quad\mbox{ for all $K\in{\mathcal{M}}$},$$

$$F_{K,\sigma }^{(n+1)}+F_{L,\sigma }^{(n+1)}=0$$

$$F_{K,\sigma }^{(n+1)}=0$$

$$\sum _{K\in \mathcal{M}}|K|p_K^{n+1}=0.$$

$${\int }_{t^{(n)}}^{t^{(n+1)}}{\int }_{\Omega }\left(\varphi \frac{\partial (c\psi )}{\partial t}+\mathbf{D}(\mathbf{x},\mathbf{u})\nabla c\cdot \nabla \psi \right)-{\int }_{t^n}^{t^{n+1}}{\int }_{\Omega }c(\varphi {\psi }_t+\mathbf{u}\cdot \nabla \psi )={\int }_{t^n}^{t^{n+1}}{\int }_{\Omega }{q}_c\psi .$$

$$\frac{d\hat{\mathbf{x}}}{dt}(t)=\frac{{\mathbf{u}}^{(n+1)}(\hat{\mathbf{x}}(t),t)}{\varphi (\hat{\mathbf{x}}(t))},\hat{\mathbf{x}}(t^{(n+1)})=\mathbf{x}.$$

$$\hat{K}=\{\mathbf{x}(t^{(n)}):\mathbf{x}(t^{(n+1)})\in K\}$$

$$\forall K\in \mathcal{M},{\mathbf{U}}_{|K}^{(n+1)}=\frac{{\mathbf{K}}_K}{\mu (c_K^{(n)})}{\overline{\nabla }}_K{p}^{(n+1)}.$$

$$\begin{gathered} \forall K\in \mathcal{M},\forall v\in X_{\mathcal{T}}, \\ \sum _{\sigma \in {\mathcal{E}}_K}{D}_{K,\sigma }^{(n+1)}(v_K-v_{\sigma })={\int }_K\mathbf{D}(\mathbf{x},{\mathbf{U}}_K^{(n+1)}){\nabla }_H{c}^{(n+1)}(\mathbf{x})\cdot {\nabla }_{H}v(\mathbf{x})d\mathbf{x}. \end{gathered}$$

$$\forall K\in \mathcal{M},{\int }_K\varphi c^{(n+1)}-{\int }_{\hat{K}}\varphi c^{(n)}+\Delta t\sum _{\sigma \in {\mathcal{E}}_K}{D}_{K,\sigma }^{(n+1)}={\int }_{t^n}^{t^{n+1}}{\int }_{\Omega }{q}_{c^{(n+1)}}{\psi }_K.$$

$${\int }_K\varphi c^{(n+1)}={\varphi }_K|K|c_K^{(n+1)}\text{mbox}and{\int }_{\hat{K}}\varphi c^{(n)}=\sum _{M\in \mathcal{M}}{\varphi }_M|\hat{K}\cap M|c_M^{(n)},$$

$$\begin{gathered} \forall K\in \mathcal{M}, \\ {\varphi }_K|K|c_K^{(n+1)}+\Delta t\sum _{\sigma \in {\mathcal{E}}_K}{D}_{K,\sigma }^{(n+1)}=\sum _{M\in \mathcal{M}}{\varphi }_M|\hat{K}\cap M|c_M^{(n)}+{\int }_{t^n}^{t^{n+1}}{\int }_{\Omega }{q}_{c^{(n+1)}}{\psi }_K. \end{gathered}$$

$$D_{K,\sigma }^{(n+1)}+D_{L,\sigma }^{(n+1)}=0$$

$$D_{K,\sigma }^{(n+1)}=0$$

$$\forall K\in \mathcal{M},{\varphi }_K|K|{\tilde{c}}_K^{(n+1)}=\sum _{M\in \mathcal{M}}{\varphi }_M|\hat{K}\cap M|c_M^{(n)}+{\int }_{t^n}^{t^{n+1}}{\int }_{\Omega }{q}_{c^{(n+1)}}{\psi }_K.$$

$$\forall K\in \mathcal{M},{\varphi }_K|K|(c_K^{(n+1)}-{\tilde{c}}_K^{(n+1)})+\Delta t\sum _{\sigma \in {\mathcal{E}}_K}{D}_{K,\sigma }^{(n+1)}=0$$

$$\frac{1}{|T_{K,{\sigma }_i}|}\left(F_{{\sigma }_i^{\ast }}-F_{{\sigma }_{i-1}^{\ast }}+F_{K,{\sigma }_i}^{(n+1)}\right)=\frac{1}{|K|}\sum _{\sigma \in {\mathcal{E}}_K}{F}_{K,\sigma }^{(n+1)},\text{ for }i=1,2,\dots ,r.$$

$$\sum _{i=1}^m{\alpha }_i{\mathbf{v}}_i={\mathbf{x}}_K\text{and}\sum _{i=1}^m{\alpha }_i=1.$$

$$\sum _{i=1}^m{\alpha }_i{F}_{{\sigma }_i^{\ast }}=0.$$

$$F_{{\sigma }_i^{\ast }}={\int }_{{\sigma }_i^{\ast }}\mathbf{u}\cdot {\mathbf{n}}_{{\sigma }_i^{\ast }}=(a{\mathbf{x}}_{{\sigma }_i^{\ast }}+\mathbf{b})\cdot {\mathbf{n}}_{{\sigma }_i^{\ast }}|{\sigma }_i^{\ast }|=(a{\mathbf{x}}_{{\sigma }_i^{\ast }}+\mathbf{b})\cdot \mathrm{Rot}({\mathbf{x}}_K-{\mathbf{v}}_i).$$

$$\sum _{i=1}^m{\alpha }_i{F}_{{\sigma }_i^{\ast }}$$

$$=\sum _{i=1}^m{\alpha }_i(a({\mathbf{x}}_{{\sigma }_i^{\ast }}-{\mathbf{x}}_K)+a{\mathbf{x}}_K+\mathbf{b})\cdot \mathrm{Rot}({\mathbf{x}}_K-{\mathbf{v}}_i)$$

Full text

1. Introduction

This study focuses on the recovery of oil in a process known as miscible displacement, in which a solvent, such as a short-chain hydrocarbon or pressurised carbon-dioxide, is injected into the oil reservoir to reduce the viscosity of the resident oil and push it towards a production well. One of the models that describes the said process is the Peaceman model, which was first introduced in [34].

Let $\Omega$ be a bounded domain in $\mathbb{R}^{d}$ and $[0,T]$ be a time interval. Denote by $\mathbf{K}(\mathbf{x})$ and $\phi(\mathbf{x})$ the permeability tensor and the porosity of the medium, respectively. Then, neglecting gravity, the Peaceman model is given by:

[TABLE]

The functions $q^{+}$ and $q^{-}$ represent the injection and production wells respectively, and $\mathbf{D}(\mathbf{x},\mathbf{u})$ denotes the diffusion tensor

[TABLE]

Here, $d_{m}$ is the molecular diffusion coefficient, $d_{l}$ and $d_{t}$ are the longitudinal and transverse dispersion coefficients respectively, and $\mathcal{P}(\mathbf{u})$ is the projection matrix along the direction of $\mathbf{u}$. Also, $\mu(c)=\mu(0)[(1-c)+M^{1/4}c]^{-4}$ is the viscosity of the fluid mixture, where $M=\mu(0)/\mu(1)$ is the mobility ratio of the two fluids. As usually considered in numerical tests, we consider no-flow boundary conditions, and zero initial conditions for the concentration:

[TABLE]

The pressure is fixed by imposing a zero average:

[TABLE]

The purpose of this work is to design and test a numerical scheme for the complete coupled model (1). This scheme is based on an adjoint method of characteristic for the advection in (1b) and a finite volume method for the diffusion terms. The main contributions of our work are:

• •

Usage of the Hybrid Mimetic Mixed (HMM) method for the diffusion terms, that is applicable on meshes with very generic geometries as encountered in real-world applications.

• •

Reconstruction of a velocity field from the numerical Darcy fluxes, that preserves the exact divergence and is locally consistent on $\mathbb{RT}_{0}$ functions; this reconstructed velocity is required to apply the characteristic method.

• •

A new treatment of cells near the injection well that accounts better for the conservativity of the injected fluid.

An early version of the HMM–ELLAM was presented in [10]. This version however did not include the special treatment of cells near the injection well, or the new velocity reconstruction. We show here that the combination of these novel elements lead to improved numerical results compared with this initial version.

A number of various numerical methods have been used to approximate the solutions of (1), from finite difference techniques [32, 17], to finite element schemes [16, 25], to discontinuous Galerkin methods [5, 37], to finite volume methods [7, 8]. Closer to the focus of this article are the Modified Method of Characteristics (MMOC) and the Eulerian-Lagrangian Localized Adjoint Method (ELLAM). These methods, designed to deal with the advective terms, have been applied in conjunction with mixed finite elements (for the pressure equation) and conforming finite elements (for the diffusion terms in the concentration equation) in [9, 24, 41], among others. The combination with FE methods severely restricts the cell geometries that can be managed with such methods. On the contrary, recent finite volume (FV) methods can accommodate very generic mesh geometries, see the review [18] and references therein. Among those, the hybrid mimetic mixed (HMM) method of [21] is a family of numerical schemes for diffusion equations, applicable on generic meshes, which gathers three separately developed numerical methods: hybrid finite volumes [26], mixed-hybrid mimetic finite differences [6], and mixed finite volume [19]. The HMM has been adapted in [7] to the model (1), using an upwind discretisation of the advective term. The drawback of such a discretisation is an additional numerical diffusion, which leads to a widening of the transition layer between the regions $c\approx 1$ and $c\approx 0$.

In this work, we propose to combine the HMM method, for its usability on generic grids, with an ELLAM for the advection, for its reduced numerical diffusion compared to the upwinding method. The HMM method being a FV scheme, using it to discretise the pressure equation (1a) naturally produces numerical fluxes of the Darcy velocity across the mesh faces. The ELLAM however requires a complete velocity field, to track points along it (see equation (8)). To avoid creating artificial wells, this velocity field should have the exact same divergence as prescribed by the numerical Darcy fluxes. To reconstruct such a divergence-preserving field starting from the Darcy fluxes, we use ideas of [31]: the cells are split in simplices (triangles in 2D) and an $\mathbb{RT}_{0}$ velocity field is reconstructed on these simplices such that its divergence on each simplex is equal to the balance of the fluxes on the cell (discrete divergence of the Darcy velocity in the cell). The corresponding system is under-determined and the original method of [31] suggests to select the solution with a certain minimal norm. We propose an alternative approach by using, in the case of dimension 2, a complementary equation which is consistent on functions that are globally $\mathbb{RT}_{0}$ in the entire cell. In other words, if, on a given cell, the numerical Darcy fluxes obtained by solving (1a) are the fluxes of a field of the form $\mathbf{a}+b\mathbf{x}$ on the cell, then the reconstructed velocity is exactly equal to that field in the cell. This construction is performed in dimension $d=2$ here, but already shows promising results; the extension to the three-dimensional case being the purpose of a future work.

As most low order FV methods, the HMM method approximates the concentration by piecewise constant functions in the cells. The test functions are therefore also piecewise constant, and the ELLAM requires to track them along the reconstructed velocity field. This boils down to tracking each cell along the velocity field. This tracking is easier to perform than that of FE basis functions, which are intrinsically more complicated than piecewise constant functions; this is another advantage, in the ELLAM framework, of FV methods over FE methods.

Each cell is tracked by transporting its vertices and edge midpoints (and possibly more points along the edge, depending on the cell regularity) by the reconstructed Darcy velocity; the tracked cell is then approximated by the polygon formed by these tracked points. For cells tracked into the injection well, to mitigate the steepness of the source terms there, we use a different strategy, based on a trace forward algorithm inspired by [4]. The original algorithm of [4] however sometimes lead to an excessive numerical injection of fluid. We therefore modify this specific treatment of cells near the injection well to eliminate this excess. Aside from this, the integral of the source term should be treated properly; otherwise, mass conservation will not be satisfied. To this end, we use a weighted trapezoidal rule, as described in [2].

The paper is organised as follows. The HMM–ELLAM scheme is presented in Section 2, starting with the HMM scheme for the pressure equation (1a). We then proceed, in section 2.2, by presenting an ELLAM scheme for the advective component of concentration equation (1b), whilst still discretising the diffusive component using an HMM scheme. Discretisation of the source term should be treated carefully, and so it will be presented separately in section 2.5. In section 2.3, we present the reconstruction of the velocity, first by using the original method of [31] involving the minimisation of some norm, then by presenting our approach based on a consistent complementary equation. Numerical results are presented, on various types of meshes, in section 3. These results show that the reconstructed velocity based on the consistent complementary equation behaves better, on distorted meshes, than the velocity chosen to have a minimal norm. We also provide qualitative and quantitative comparisons with the results obtained from MFV–upwind and MFEM–ELLAM schemes, outlining the positive and negative qualities of the numerical solutions obtained from each scheme. A few comments about the proper implementation of ELLAM are also made in Section 3, such as the proper number of points to take for tracking each cell, and the impact of choosing the proper quadrature rule for integrating the source term.

2. The HMM–ELLAM

The numerical solution of (1) is usually approximated by using a two-step process. Starting from a known value $c^{(n)}$ of $c$ at time level $n$ (for $n=0$, $c^{(0)}=0$), a numerical solution $p^{(n+1)}$ for $p$ at time level $n+1$ is computed by approximating (1a) with $c=c^{(n)}$. This computation also provides an approximation $\mathbf{u}^{(n+1)}$ of the Darcy velocity at time level $n+1$, and possibly of secondary quantities (e.g., fluxes). The concentration $c^{(n+1)}$ at time level $n+1$ is then computed by approximating (1b) by using $\mathbf{u}=\mathbf{u}^{(n+1)}$ and the aforementioned secondary quantities.

For our discretisation, we will consider polytopal meshes as defined in [20], in dimension $d=2$ (as explained above, this restriction is due to a special consistent reconstruction of the Darcy velocity from the numerical fluxes, whose extension to dimension 3 is the purpose of an ongoing work). Thus, $\mathcal{T}=({\mathcal{M}},{\mathcal{E}},\mathcal{P})$ are the set of cells $K$, edges $\sigma$, and points $\mathbf{x}_{K}$ (one per cell) of our mesh, respectively. The cells can be any star-shaped polygons. For a cell $K\in{\mathcal{M}}$, ${{{\mathcal{E}}}_{K}}\subset{\mathcal{E}}$ denotes the set of edges of the cell $K\in{\mathcal{M}}$. We also denote by $|K|$ the area of a cell $K$, and $|\sigma|$ the length of an edge $\sigma$. As is usual, and not restrictive with respect to applications, we assume that the permeability and porosity are constant in each cell, and we write $\mathbf{K}_{K}$ and $\phi_{K}$ their respective values in the cell $K\in{\mathcal{M}}$.

2.1. Numerical scheme for the pressure equation

The variational formulation for (1a) is given by

[TABLE]

There is an implicit time variable above, which is here just a parameter. The degrees of freedom of our numerical scheme are given by cell and edge values. The space of DOFs is then

[TABLE]

For HMM schemes, a piecewise constant gradient is then defined on a sub-triangulation of cells. Following [21], if $K\in{\mathcal{M}}$ and $(T_{K,\sigma})_{\sigma\in{{{\mathcal{E}}}_{K}}}$ are the triangles defined by $\mathbf{x}_{K}$ and the edges of $K$ (see Figure 1), we set

[TABLE]

Note that $\overline{\nabla}_{K}w$ is a linearly exact reconstruction of the gradient, that is, if $(w_{\sigma})_{\sigma\in{{{\mathcal{E}}}_{K}}}$ interpolate an affine function $A$ at the edge midpoints, then $\overline{\nabla}_{K}w=\nabla A$. In (4), $d_{K,\sigma}$ is the orthogonal distance between $\mathbf{x}_{K}$ and $\sigma$, and the second term can be seen as a stabilisation of the gradient involving a discrete 2nd order Taylor expansion.

Given $p^{(n+1)}\in X_{\mathcal{T}}$, the discrete Darcy fluxes $(F_{K,\sigma}^{(n+1)})_{K\in{\mathcal{M}},\,\sigma\in{{{\mathcal{E}}}_{K}}}$, approximations of $-\int_{\sigma}(\mathbf{K}/\mu(c))\nabla p\cdot\mathbf{n}_{K,\sigma}$, are then defined by

[TABLE]

where $c^{n}_{K}$ is a cell-centred approximate concentration at the previous time step, which is available since $c$ is also approximated in $X_{\mathcal{T}}$.

Fundamental in the formulation of finite volume schemes is the balance and conservativity of fluxes [18], which are given by the equations

[TABLE]

With the definitions (4)–(5), the system (6) provides an approximation $p^{(n+1)}\in X_{\mathcal{T}}$ of $p$ at time $t^{(n+1)}$, as well as approximate Darcy fluxes $(F_{K,\sigma}^{(n+1)})_{K\in{\mathcal{M}}\,,\;\sigma\in{{{\mathcal{E}}}_{K}}}$ at the same time.

2.2. Numerical scheme for the concentration equation

For the concentration equation (1b), we have the following variational formulation:

[TABLE]

A common feature of ELLAM schemes is taking test functions $\psi$ such that $\phi\psi_{t}+\mathbf{u}^{n+1}\cdot\nabla\psi=0$ in the sense of distributions. This equation is completed by a terminal condition at $t=t^{(n+1)}$, whose nature depends on the chosen approximation space for the concentration. Here, the concentration is approximated by piecewise constant functions on the mesh, and we therefore impose $\psi(t^{(n+1)},\cdot)=\mathbf{1}_{K}$, where $K$ is a generic cell and $\mathbf{1}_{K}$ is the characteristic function of $K$. We will therefore denote $\psi=\psi_{K}$.

Computing these tests functions $\psi_{K}$ requires us to use the characteristics of the advective component of (7), that is the solution, for each $\mathbf{x}\in\Omega$, to the ODEs

[TABLE]

This leads to $\psi_{K}(t^{(n)},\mathbf{x})=\mathbf{1}_{K}(\widehat{\mathbf{x}})=\mathbf{1}_{\widehat{K}}(\mathbf{x})$, where $\widehat{K}$ is $K$ tracked back from $t^{(n+1)}$ to $t^{(n)}$ through (8):

[TABLE]

Of course, solving this characteristic equation first requires access to a velocity $\mathbf{u}^{n+1}$. We describe in Section 2.3 how to reconstruct, from the discrete Darcy fluxes obtained in Section 2.1, such a velocity.

The diffusive term is discretised using the HMM method and an implicit time stepping scheme. Fluxes $D_{K,\sigma}^{(n+1)}$ are defined as for the pressure equation (1a), using a piecewise constant Darcy velocity $\mathbf{U}^{(n+1)}$ constructed from the consistent part of the reconstructed pressure gradient (4) and the viscosity at $c^{(n)}$, that is

[TABLE]

Hence, if $c^{(n+1)}\in X_{\mathcal{T}}$ is our discrete concentration sought at time $t^{(n+1)}$, the fluxes are given by

[TABLE]

Note that, on the cell $K$, $\mathbf{D}(\cdot,\mathbf{U}^{(n+1)}_{K})$ is actually constant.

Without discretising the source term yet (this will be discussed in Section 2.5), this leads to the following scheme for (7):

[TABLE]

The integrals are then computed by writing

[TABLE]

which leads to the following discretised form of the concentration equation

[TABLE]

As with the pressure equation, we also have the conservativity of fluxes and no flow boundary conditions given by

[TABLE]

Remark 2.1** (Order of accuracy of the scheme).**

Under regularity assumptions on the mesh and the solutions, error estimates in $H^{1}$ norm for an HMM scheme for an anisotropic diffusion equation similar to (1a) are of order $O(\Delta x)$ [20, Theorem 3.31 and Propositions 13.15, 13.16]. Error estimates for an $\mathbb{RT}_{0}$–ELLAM scheme for an advection–diffusion equation similar to (1b) are obtained in [40], and are of order $O(\Delta x+\Delta t)$. However, for the complete Peaceman model (1), no exact error estimate may be established due to the fact that most of the regularity assumptions violate the physical interpretation of the problem. Nevertheless, convergence of the scheme can still be established, see [11].

Remark 2.2** (Splitting technique).**

The HMM–ELLAM for the concentration equation may also be viewed as an operator splitting technique, which is the usual approach when implementing ELLAM schemes [2, 38]. Introduce the values $\tilde{c}^{(n+1)}=(\tilde{c}^{(n+1)})_{K\in{\mathcal{M}}}$ defined by

[TABLE]

The approximate values $\tilde{c}^{(n+1)}$ correspond to applying one time step of the ELLAM technique on the pure advection equation $\phi\partial_{t}c+\nabla\cdot(\mathbf{u}c)=q_{c}$. Re-starting from $\tilde{c}^{(n+1)}$, construct now $c^{(n+1)}\in X_{\mathcal{T}}$ such that

[TABLE]

with boundary conditions (10). In other words, $c^{(n+1)}$ is obtained from $\tilde{c}^{(n+1)}$ by performing one time step of the HMM scheme for the pure diffusion equation $\phi\partial_{t}c-\nabla\cdot(\mathbf{D}(\mathbf{x},\mathbf{u})\nabla c)=0$. Adding together (11) and (12), we see that $c^{(n+1)}$ is the solution to the complete HMM–ELLAM scheme (9).

2.3. Reconstruction of Darcy velocities

In order to compute the characteristics in (8), two important features of the velocity need to be accounted for: the no-flow boundary conditions $\mathbf{u}^{(n+1)}\cdot\mathbf{n}=0$ on $\partial\Omega$, which ensures that the solutions to (8) do not exit the computational domain, and the preservation of the divergence in (1a), to avoid creating regions with artificial wells or sinks (which lead to non-physical flows), as was mentioned and implemented in [3, 4]. We present this construction in dimension $d=2$, the extension to dimension $d=3$ being the purpose of a future work.

Preserving these features is done by using a technique similar to that of [31]. Each cell $K\in{\mathcal{M}}$ is split into triangles (see Fig. 1) by choosing an interior point $\mathbf{x}_{K}$ and connecting it to the vertices of cell $K$. These vertices are denoted by $(\mathbf{v}_{1},\ldots,\mathbf{v}_{r})$, in counter clockwise order, and with the convention that $\mathbf{v}_{0}=\mathbf{v}_{r}$. We let $\sigma_{i}^{*}$ be the internal edge $[\mathbf{v}_{i},\mathbf{x}_{K}]$ of the triangular subdivision of $K$, and $\sigma_{i}$ be the edge $[\mathbf{v}_{i-1},\mathbf{v}_{i}]$ of $K$. A clockwise oriented interior flux $F_{\sigma_{i}^{*}}$ is then computed on each internal edge $\sigma_{i}^{*}$ of this subdivision. $\mathbf{u}^{(n+1)}$ is the $\mathbb{RT}_{0}$ function reconstructed from these fluxes on the triangular subdivision.

This reconstruction belongs to $H_{\mathop{\rm div}}(\Omega)$ and, to ensure that its divergence is identical to the discrete divergence computed from the Darcy fluxes, the internal fluxes $F_{\sigma_{i}^{*}}$ are constructed so that their balance (along with the fluxes $F_{K,\sigma_{i}}^{(n+1)}$ at the boundary of $K$) in each triangle corresponds to the balance over the cell $K$:

[TABLE]

2.3.1. KR velocity

Although (13) gives us $r$ equations in $r$ unknowns, this system of equations is underdetermined. More specifically, its rank is $r-1$ and it therefore has an infinite number of solutions. The original method of Kuznetzov and Repin [31] consists in selecting the solution that has a minimal $l^{2}$ norm. The velocities reconstructed from these fluxes will be referred to as KR velocities.

Despite its applicability to any space dimension, this choice is questionable from the physical point of view. If, say, the external fluxes $(F_{K,\sigma}^{(n+1)})_{\sigma\in{{{\mathcal{E}}}_{K}}}$ correspond to those of a constant velocity in the cell, this reconstruction does not necessarily recover this natural velocity.

2.3.2. C velocity

We therefore propose a different approach, which consists in adding to the set of equations (13) a complementary equation which selects a unique solution. We form this complementary equation by ensuring that, as (13), it is consistent for internal and external fluxes coming from simple velocities – here, not just constant velocities in the cell, but $\mathbb{RT}_{0}$ velocities in the cell.

Lemma 2.3** (Consistency condition).**

Let $\mathbf{x}_{K}\in K$, $(\mathbf{v}_{1},\mathbf{v}_{2},\dots,\mathbf{v}_{m})$ be some vertices of cell $K$, and $(\alpha_{1},\ldots,\alpha_{m})\in\mathbb{R}^{m}$ be such that

[TABLE]

Let $\mathbf{u}$ be an $\mathbb{RT}_{0}$ function over the cell $K$, that is, $\mathbf{u}(\mathbf{x})=a\mathbf{x}+\mathbf{b}$ for some $a\in\mathbb{R}$ and $\mathbf{b}\in\mathbb{R}^{d}$. Define $F_{\sigma_{i}^{*}}=\int_{\sigma_{i}^{*}}\mathbf{u}\cdot\mathbf{n}_{\sigma_{i}^{*}}$ the flux of $\mathbf{u}$ through the internal edge $\sigma_{i}^{*}$, for all $i=1,\ldots,m$ (here, $\mathbf{n}_{\sigma_{i}^{*}}$ is the clockwise normal to $\sigma_{i}^{*}$, as in Figure 1). Then

[TABLE]

**Proof. ** Denoting by $x_{\sigma_{i}^{*}}$ the center of $\sigma_{i}^{*}$ and $\mathrm{Rot}$ the clockwise rotation by $\pi/2$, write

[TABLE]

Since $(\mathbf{x}_{\sigma_{i}^{*}}-\mathbf{x}_{K})$ is orthogonal to $\mathrm{Rot}(\mathbf{x}_{K}-\mathbf{v}_{i})$, we deduce

[TABLE]

The proof is complete by noticing that $\sum_{i=1}^{m}\alpha_{i}(\mathbf{x}_{K}-\mathbf{v}_{i})=0$ by choice of the coefficients $(\alpha_{i})_{i=1,\ldots,m}$.

This lemma tells us that by choosing 3 or more vertices of our cell $K$ in generic position, we can find an equation that closes the system (13). We only need one such equation, and we do not want to create a bias in constructing it. We therefore use all vertices $(\mathbf{v}_{1},\ldots,\mathbf{v}_{r}$) of $K$ to form this relation. Thus our closing equation will be given by

[TABLE]

where the coefficients $(\alpha_{i})_{i=1,\ldots,r}$ are related to $\mathbf{x}_{K}$ as in Lemma 2.3. Velocities reconstructed from fluxes that satisfy (13)–(14) will be denoted as C velocities (‘C’ for ‘consistent’).

Remark 2.4** (Particular $\mathbf{x}_{K}$ and barycentric combinations).**

If $\mathbf{x}_{K}$ is the iso-barycenter of the vertices of $K$, i.e. $\mathbf{x}_{K}=\frac{1}{r}\sum_{i=1}^{r}v_{i}$, then the consistency relation is simply $\sum_{i=1}^{r}F_{\sigma_{i}^{*}}=0$. If $\mathbf{x}_{K}$ is the center of mass of $K$, then a consistency relation is

[TABLE]

where $T_{K,\sigma_{i-1}}$ is the triangle that shares edge $\sigma_{i-1}^{*}$ with $T_{K,\sigma_{i}}$.

The system (13)–(14) has a unique and explicit solution. Indeed, set

[TABLE]

The system (13) is then

[TABLE]

From these, we easily deduce that

[TABLE]

By noticing that $\sum_{i=1}^{r}a_{i}=0$, we see that (15) also holds for $k=r$. Together with the closing relation (14), and since $\sum_{k=1}^{r}\alpha_{k}=1$, we obtain an explicit expression for $F_{\sigma^{*}_{r}}$, given by

[TABLE]

Equation (16), together with (15), will then give us explicit expressions for $F_{\sigma^{*}_{k}},$ $k=1,2,\dots,r$.

These computations show an advantage of this method over the technique consisting in selecting a minimal norm solution of (13). Here, we do not need to solve any local system, as we have explicit expressions of the solution to (13)–(14). Moreover, the reconstructed velocity is more precise, especially on skewed meshes. As an example, consider a constant velocity field $V=(0,1)$. Pick a cell $K$ and define the fluxes across its boundary edges to be those of $V$. Denote then by $V_{\rm KR}$ and $V_{\rm C}$ the KR and C velocities, respectively, reconstructed as above on a triangular sub-mesh of $K$.

Table 1 represents the relative error obtained between the exact and reconstructed velocities, on a variety of mesh geometries. For square cells (Cartesian mesh), both methods reconstruct the velocity accurately. However, for cells from hexahedral (Fig. 4, right), non-conforming, and Kershaw meshes (Fig. 5), KR velocities noticeably deviate from the actual velocity, by more than $30\%$ on distorted cells. On the other hand, as expected, using the consistency relation as a closure equation enables us to recover the velocity $V$ up to machine precision, regardless of the mesh.

2.4. Approximate traceback region, and tracking points through vertices

In general, we cannot get an exact representation of the region $\widehat{K}$, and thus, for each cell $K$, the traceback region $\widehat{K}$ is approximated in the following manner: we select points $(\mathbf{x}_{i})_{i=1,\ldots,\ell_{K}}$ along the boundary of $K$ (at least all the vertices and edge midpoints are selected), we solve (8) starting from any of these points, thus getting curves $(\widehat{\mathbf{x}}_{i})_{i=1,\ldots,\ell_{K}}$, and we approximate $\widehat{K}$ by the polygon defined by the points $(\widehat{\mathbf{x}}(t^{(n)}))_{i=1,\ldots,\ell_{K}}$.

The reconstructed velocity $\mathbf{u}^{(n+1)}$ is an $\mathbb{RT}_{0}$ function on a triangular sub-mesh. Tracking a point through (8) is naturally done cell-by-cell, using the value of $\mathbf{u}^{(n+1)}$ in a cell $K$ as long as $\widehat{\mathbf{x}}$ stays in $K$, and then, when $\widehat{\mathbf{x}}$ exits $K$ to enter a cell $L$, continuing the tracking by using the value of $\mathbf{u}^{(n+1)}$ in $L$. This type of tracking, with point of exit determined by minimal time of flight, was first implemented by Pollock [35] on meshes characterised by orthogonal grid blocks, e.g. Cartesian meshes. Pollock’s algorithm was then extended to more generic types of cells in [36], and further improved in [30]. Because the fluxes of $\mathbf{u}^{(n+1)}$ are continuous across the edges, this tracking procedure ensures that a point will never do a U-turn, that is, $-\mathbf{u}^{(n+1)}_{|L}$ (we use $-\mathbf{u}^{(n+1)}$ since we are tracking backwards) will not force $\widehat{\mathbf{x}}$ to re-enter $K$ before even entering $L$ (this would in effect freeze $\widehat{\mathbf{x}}$ on the interface between $K$ and $L$).

During this tracking, special care must be taken with points that start or pass through a vertex. An initial position $\mathbf{x}$ corresponding to a vertex is involved in several triangles, and could thus be initially tracked using any of the Darcy velocities in these triangles (see Fig. 2). To avoid non-physical initial tracking, we compute the Darcy velocities at the vertex $\mathbf{x}$ in each of the triangles involved with it. Picking one of these Darcy velocities at random is not acceptable, since if its opposite vector points outside the corresponding triangle, this means that $\mathbf{x}$ would never be tracked back inside this triangle, and that the chosen Darcy velocity is thus not the correct one.

We therefore loop over the triangles $T_{1}$, $T_{2}$, etc. until we reach a triangle $T_{n}$ such that $-\mathbf{u}^{(n+1)}_{|T_{n}}(\mathbf{x})$ points inside $T_{n}$. In Figure 2, this corresponds to triangle $T_{2}$. As a note, regardless of the mesh, our numerical tests suggest that such a triangle always exists. For points that will be tracked into a vertex at some time in $[t^{(n)},t^{(n+1)})$ (see Fig. 3), we consider the negative of the Darcy velocity in the triangle it came from (in this case, $T_{4}$) and determine which triangle it points into (in this case, it points into triangle $T_{2}$); the tracking is then continued based on the reconstructed velocity in this latter triangle.

2.5. Source term

The integral involving the source term should be treated carefully, otherwise the numerical results will feature severe overshoots, especially over the regions around the injection well. For our discretisation, we use a weighted trapezoid rule in time:

[TABLE]

where $q_{c^{(n+1)}}(s)=q^{+}(s)-q^{-}(s)c^{(n+1)}$ (that is, the source term is still fully implicit in $c$).

Note that the left and right quadrature rules correspond to $w=1$ and $w=0$, respectively. To determine the proper weight, we consider an injection cell $K=E$. We mainly focus on injection cells since this is where mass conservation might fail locally. A proper weight that will yield mass conservation has been derived for Cartesian meshes on [2]. This can easily be generalised for arbitrary meshes, and is given by

[TABLE]

Then, for each cell $K$ (injection or not), the integral of the source term is computed using the weighted trapezoid rule (17), where $w$ is obtained as in (18), for some $E$ related to $K$ – see below. We treat the computation of the integral over $\widehat{K}$ in the right hand side of (17) in different manners, depending on whether the cell $K$ is

1. i)

an injection cell, 2. ii)

a cell tracked back into an injection cell (but not an injection cell itself), 3. iii)

or a cell that does not track back into an injection cell.

i) If the cell $K$ is an injection cell $E$, then it tracks back entirely into itself. Hence, over the entire interval $[t^{(n)},t^{(n+1)}]$, $\nabla\cdot\mathbf{u}^{(n+1)}(\widehat{x})=\frac{1}{|E|}\int_{E}q(t^{(n+1)})$ and thus we use, through Liouville’s theorem, an exact computation

[TABLE]

where $\alpha$ is given by (18).

ii) If the cell $K$ is not an injection cell, but is tracked back (at least partially) into an injection cell $E$, then we use a forward tracking algorithm similar to that described in [4]. For these regions, we compute the trace forward region $\widetilde{E}$ of $E$, in which we solve the characteristics (8) with initial condition $\widehat{\mathbf{x}}(t^{(n)})=\mathbf{x}$ instead. The trace forward region $\widetilde{E}$ is then approximated using polygons in a similar manner as the trace back regions in Section 2.4. The integral of the first term on the right hand of (17) is then approximated by

[TABLE]

In physical terms, this means that the volume injected from the well $E$ is transported into each of the cells $K$ proportionally. Note that, on the contrary to [4], only a fraction $(1-e^{-\alpha})$ of $\int_{E}q_{c^{(n+1)}}$ is being spread in the cells $K$ around $E$, since a fraction $e^{-\alpha}$ of $\int_{E}q_{c^{(n+1)}}$ has already been allocated to $E$, as can be seen in (19).

iii) Finally, if a cell $K$ does not track back into an injection cell, then the integral $\int_{\widehat{K}}q_{c^{(n+1)}}$ will be computed using the (approximate) trace back regions as described in Section 2.2. Actually, in that situation, either:

• •

$K$ is not a production cell, $\widehat{K}$ is disjoint from injection cells and production cells (due to $-\mathbf{u}^{(n+1)}$ pointing outside production cells). Therefore, the integrals in (17) are both equal to zero.

• •

or $K$ is a production cell, in which case, by nature of the Darcy flow, it is expected that $K\subset\widehat{K}$, so both integrals in (17) are equal and the value of $w$ is irrelevant.

3. Numerical results

We perform numerical simulations under the following standard test case (see, e.g., [41]):

1. (1)

$\Omega=(0,1000)\times(0,1000)\text{ ft}^{2}$, 2. (2)

injection well at $(1000,1000)$ and production well at $(0,0)$, both with flow rate of $30\text{ft}^{2}/\text{day}$, 3. (3)

constant porosity $\phi=0.1$ and constant permeability tensor $\mathbf{K}=80\mathbf{I}$ md, 4. (4)

oil viscosity $\mu(0)=1.0$ cp and mobility ratio $M=41$, 5. (5)

$\phi d_{m}=0.0\text{ft}^{2}/\text{day}$, $\phi d_{l}=5.0\text{ft}$, and $\phi d_{t}=0.5\text{ft}$

For the time discretisation, we take a time step of $\Delta t=36$ days. These will be simulated on Cartesian type meshes (square cells of dimension $62.5\times 62.5$ ft), hexahedral meshes (see Fig.4), non-conforming meshes, and finally on Kershaw type meshes, as described in [28] (see Fig. 5).

3.1. Effect of the quadrature rule

The following simulations are based on KR velocities (see Section 2.3.1). Figure 6 shows the numerical solution for the concentration at $t=10$ years on a Cartesian mesh using the left and the right rule, respectively. These results support the observations made in [10], i.e., that the left and right hand quadrature rules provide severe underestimates and overshoots, respectively, at the injection well, and are thus not good choices.

Figure 7 (left) shows the numerical solution for the concentration using the proper weight for the trapezoidal rule, and computation of the integrals as described in Section 2.5. This presents a significant improvement from the results obtained through the right and left rule.

The overshoot seen in this figure is at worst around $7\%$, which is commensurate with (or even less than) overshoots already noticed in other other characteristic methods in the absence of specific post-processing [23].

3.2. Comparison with forward tracking in [4]

We now compare our numerical results to the original algorithm of [4]. Instead of implementing i) and ii) as described in Section 2.5,

• •

For an injection cell $E$, $c_{E}^{(n+1)}$ is fixed at $1$, as this is the concentration of the injected solvent.

• •

for cells $K$ tracked back (at least partially) into an injection cell $E$, the following approximation is used:

[TABLE]

where $\widetilde{E}$ is the trace forward region of $E$.

For some discretisation parameters, this implementation might lead to degraded results due to its physical implications. Setting $c_{E}^{(n+1)}=1$ corresponds to distributing a fraction of $\int_{E}q_{c^{(n+1)}}$ into injection cells $E$. In this instance, a good estimate would be given by (19). However, computation of the integral for cells $K$ that track back into an injection cell $E$ by (20) means that we spread the whole of $\int_{E}q_{c^{(n+1)}}$ onto the cells $K$. This then means that at time level $n+1$, an excessive amount of $e^{-\alpha}\int_{E}q_{c^{(n+1)}}$ of fluid has been injected in the cells around $E$. If $\Delta t$ is large enough, then this is a negligible excess as $e^{-\alpha}\approx 0$. However, for moderate to small $\Delta t$, the numerical results do not model the physical phenomenon properly. It is important to note that even though characteristic methods aim for computations using large time steps, we should still have a proper numerical result even when the time steps are small.

Figure 7 (right) shows the numerical solutions obtained at $t=10$ years upon computing the integrals as in [4], with the moderate time step $\Delta t=36$ days. As expected, due to injection of too much fluid, the overshoot at the right of Figure 7 (around 9.9%) is larger than the one on the left (around 7.3%). This particular feature is even more evident if we take smaller time steps. Due to this, we see that the implementation we propose in Section 2.5 is more accurate.

3.3. A criterion for choosing the number of points per edge

In general, a polygon formed by tracking back only vertices and edge midpoints might not give a good approximation to the trace back region $\hat{K}$. The mesh regularity parameter, defined as

[TABLE]

has been used as a criterion for determining the proper number of points to track per edge (see Table 2). To be exact, $\lceil\log_{2}(m_{{\mathcal{M}}\text{reg}})\rceil$ points per edge should be tracked in order to obtain a reasonable numerical solution. These results have been established in [10] for KR velocities. For this paper, we verify that this still holds for C velocities, and on the previously not considered non-conforming mesh.

We start by demonstrating on hexahedral meshes that even for C velocities, tracking only vertices and edge midpoints do not give a good approximation (see Fig. 8 left). As expected, taking 3 points per edge, as suggested in Table 2, then gives a better result, with an overshoot $<7\%$ (see Fig. 8 right).

This heuristic choice of number of points along each edge is further backed up by the numerical solutions for the non-conforming meshes, and also for the very distorted ‘Kershaw’ meshes, see Figures 11 and 12.

We note that not all cells are highly ”irregular”; thus, tracking a lot of points along each edge for the whole mesh introduces unnecessary numerical cost. If the cell $K$ is an injection or production cell, then we still track, as before, $\lceil\log_{2}(m_{{\mathcal{M}}\text{reg}})\rceil$ points along each edge. As an improvement, if $K$ is neither an injection or production cell, we determine the number of points to track along each edge of cell $K$ by measuring instead the cell regularity parameter defined to be

[TABLE]

and track $\lceil\log_{2}(m_{K\text{reg}})\rceil$ points along each edge of cell $K$. By doing so, we reduce the computational cost without degrading the quality of the numerical solutions.

3.4. Comparison of the two reconstructions of the Darcy velocity

In this section, we compare the numerical solutions at $t=10$ years obtained through simulations using KR velocities to those obtained when we do our simulations using C velocities. These will be performed over Cartesian, hexahedral, non-conforming and Kershaw type meshes using the recommended number of points per edge (see Table 2) for tracking.

Figures 9 to 11 compare the numerical solutions at $t=10$ years, on Cartesian, hexahedral, and non-conforming meshes, between KR velocities and C velocities. We notice that for these types of meshes, the numerical solutions obtained using KR velocities are essentially the same as those obtained using C velocities.

Now, looking at the numerical results along Kershaw type meshes on Figure 12, we see that C velocities give us a slightly better numerical result. The maximum overshoot is now around $3.1\%$, as opposed to a maximum overshoot of around $4.5\%$ for KR velocities. Moreover, we observe that the solutions on the Cartesian, hexahedral, and non-conforming meshes are very similar, showing a certain robustness of the method with respect to the choice of mesh. However, the solution on the Kershaw mesh is noticeably different, which signals the presence of a grid effect. In the next section, we consider streamlines to back up our claim that C velocities yield a better numerical approximation than KR velocities. These streamlines also enable us to understand why grid effects are present in Kershaw type meshes.

3.5. Streamlines resulting from the reconstructed Darcy velocities

We start by plotting the streamlines for the velocities reconstructed using both KR and C velocities on Cartesian, hexahedral, and non-conforming meshes, which can be seen in figures 13 to 15. For the following figures, the particles are assumed to have travelled for 3600 days, which is approximately 10 years.

As can be seen in figures 13 to 15, the streamlines along the Cartesian, hexahedral, and non-conforming meshes are quite similar (whether we use KR velocities or C velocities), which explains why the numerical solutions obtained for the concentration are essentially the same, regardless of which velocity we use. We note however, that upon comparing the streamlines of hexahedral type meshes against those of Cartesian or non-conforming meshes, that there is a slight difference in how the fluid travels. In particular, we take note that the streamline arising from the rightmost position of the plots ends near position (770, 60) for the hexahedral meshes, and near position (570, 60) for the other 2 meshes. This particular difference can also be seen in the concentration profiles upon comparing Figure 10 with Figures 9 and 11. Note here that the concentration profile obtained on a hexahedral mesh exhibits a sharper fingering effect along the diagonal, as compared to the concentration profile obtained on the other 2 meshes. This phenomenon is caused by using fluxes generated by the low-order HMM method, which are prone to grid effects; in a future work, we will be exploring the usage, for the pressure equation only, of high-order methods such as the HHO scheme [15] to improve the quality of the numerical fluxes on distorted meshes. The effect is not as severe as the Kershaw type meshes though, since hexahedral meshes are not as distorted.

The streamlines presented in Figure 16 show that, on Kershaw meshes, the flow resulting from the C velocities is more natural than that of KR velocities, which exhibits more staggering in the upper region. This explains why the numerical solution for the concentration obtained using the C velocities is slightly better than that of the KR velocities. Noting that the streamlines obtained along all 4 types of meshes are quite similar, we have still yet to address the grid effects that were present for the Kershaw mesh. We understand this by tracking the particles in the streamline for a shorter time period of 2520 days, or approximately 7 years. For Figures 17 and 18, we will be using C velocities.

As can be seen in Figure 17, the streamlines resulting from both Cartesian and hexahedral meshes are almost symmetric with respect to the line $y=x$. A similar observation can be made for non-conforming meshes (see Fig. 18 left). On the contrary, due to the large distortion of the Kershaw mesh (see Fig. 5 right), the advection field on this mesh is such that particles travelling below the line $y=x$ reach the production well $(0,0)$ faster than those travelling above the line. In particular, upon looking at the plot on the right of Figure 18, we focus on the third streamlines from the right and top boundaries, which should be symmetric about the diagonal $x=y$. The streamline below $y=x$ has travelled near the point (275,100), whereas the streamline above the line $y=x$ has only reached some point near (150,400). Hence, the distorted mesh leads to an advection of the fluid that is skewed towards the lower part of the domain, thus leading to the numerical results that vary from those obtained in the other types of meshes. We may also compare this third streamlines to the streamlines obtained from the other types of meshes. This comparison confirms that, for Kershaw meshes, advection below the line $y=x$ is much faster than expected. Similar observations are obtained upon looking at streamlines obtained from KR velocities.

Remark 3.1**.**

Figures 13 to 18 were obtained from the velocity profile at the first time step; hence, the dependency of the velocity profile on the concentration $c$ due to a high mobility ratio of $M=41$ was not visible. To complete the presentation, we show in Figures 19 and 20 the velocity profile obtained at the final time step. Also, the particles along the streamline are assumed to have traveled 3600 days, or approximately 10 years. Indeed, upon looking at these figures side by side with Figures 9 to 12, we see the dependence of the velocity profile on the concentration, i.e., it tends to flow along the region(s) with high concentration.

3.6. Comparison with other numerical schemes

In this section, we compare the numerical results obtained from HMM–ELLAM to numerical results obtained from other schemes, such as MFV with upwinding [7] and MFEM–ELLAM [41]. This will be done on two test cases. The first test case will be done under the same test data and parameters considered above. The second test case will be done instead with an inhomogeneous permeability tensor $\mathbf{K}=20\mathbf{I}$ md over the region $(200,400)\times(200,400)\cup(200,400)\times(600,800)\cup(600,800)\times(200,400)\cup(600,800)\times(600,800)$ and $\mathbf{K}=80\mathbf{I}$ md elsewhere, while holding all other test data and parameters to be the same as those of the initial test case. This comparison is performed on a Cartesian mesh of size $50\times 50$ ft (so that the discontinuities present in the second test case are aligned with the edges of the cells); other meshes could be considered (triangular for MFEM–ELLAM, and any polytopal mesh for MFV–upwind), with similar conclusions. We perform the second test case with a time step of $\Delta t=2.5$ days, as opposed to the time step of $36$ days for the initial test case. We do this in order to have a fair comparison with the higher order scheme implemented in [1] (which was second order in time, implemented with a time step of $\Delta t=7.2$ days). As a point of reference, we present in Figures 21 and 22 the concentration profile and contour plot at $t=10$ years for the HMM–ELLAM, for the first and second test cases, respectively.

3.6.1. MFEM–ELLAM

As both MFEM–ELLAM and HMM–ELLAM are characteristic methods, tracking is implemented for both schemes for the concentration equation. Typically, HMM–ELLAM schemes only need to track the vertices, together with 1 point per edge for Cartesian type meshes, unless the time step is too large relative to the spacial discretisation (which will result to either a degenerate or self intersecting polygon approximating the trace back region). This can be avoided by reducing the time step or increasing the number of points per edge. However, MFEM–ELLAM schemes need to track a bare minimum of 3-4 points in each cell to get a correct quadrature rule to integrate the basis functions (and much more than 4 points in case these bases functions become too distorted by the tracking velocity [39]); this is why we expect the MFEM–ELLAM to be computationally more expensive than the HMM–ELLAM.

Aside from computational cost, a more important thing to consider would be the quality of the numerical solutions. Figures 23 and 24 give us the numerical solution and contour plot obtained from MFEM–ELLAM at $t=10$ years for the first and second test cases, respectively. These numerical outputs were obtained by a straight application of the MFEM–ELLAM algorithm as presented in [41], with several hundred of quadrature points per cell around the injection well.

The concentration profile and the contour plots obtained from both schemes are quite similar, with overshoots and undershoots that are typical of characteristic-based methods [23]. However, here, the overshoot in the MFEM–ELLAM (around 15% for homogeneous permeability and 20% for inhomogeneous permeability) is higher than that of the HMM–ELLAM (6% and 3% respectively) for both cases. Note that for the source terms, the MFEM–ELLAM integrates a non-constant function through quadrature rules. The main source of error encountered upon computing these integrals arise due to the presence of steep source terms. However, for the HMM–ELLAM, a different treatment of the source terms (see Section 2.5) was implemented. We recall that this can be physically interpreted as spreading the injected fluid around the region surrounding the injection well. Naturally, this resulted to the HMM–ELLAM having a smaller overshoot than the MFEM-ELLAM. We also note that the MFEM–ELLAM exhibits severe undershoots (up to around 11% for homogeneous permeability near the production cell and 90% for inhomogeneous permeability in the low-permeability regions), which is not present in the HMM–ELLAM. This severe undershoot might be due to the fact that conforming FE methods, used for solving the concentration equation, have unknowns on the vertices that sit at the permeability discontinuities. It has been noticed that, for transport of species in heterogeneous domains, schemes with unknowns at the vertices may lead to unacceptable results on coarse meshes, see [27]. On the other hand, the transition layer (from $c\approx 0$ to $c\approx 1$) is thinner for the MFEM–ELLAM than for the HMM–ELLAM.

3.6.2. MFV–upwind

Over each time step, the MFV–upwind requires, for the concentration, the solution of a linear system which has the same sparsity and number of unknowns as the HMM–ELLAM. Moreover, due to the absence of characteristic tracking and computation of integrals over trace-back regions, the computational cost of MFV-upwind scheme is much cheaper than that of the HMM–ELLAM. Next, we compare the quality of the solutions obtained by looking at Figures 25 and 26. It is quite notable that the solution remains bounded between 0 and 1 (actually, no undershoot occurs, and the overshoot is less than 0.01%). However, upwind schemes tend to introduce excess diffusion, and thus the strong viscous fingering effects we expect have been spread out. This can be seen more clearly by looking at the contour plot (Figure 25 right). Upon comparison with contour plots obtained for the HMM–ELLAM and MFEM–ELLAM schemes (Figures 21 and 23, right), we indeed see that the strong viscous fingering expected along the diagonal has been spread out by the upwind scheme. A similar conclusion can be drawn for the inhomogeneous permeability tensors upon comparing Figure 26 to Figures 22 and 24. However, upon comparing Figures 22, 24, and 26 we notice that the concentration profiles obtained on both HMM–ELLAM and MFV-upwind are quite similar. Even the transition layers obtained on the respective contour plots seem to agree with one another. Due to this, we cannot make a definitive conclusion at this point. Another point of comparison will be presented in Section 3.6.3.

3.6.3. Recovered oil

A particular quantity of interest in performing numerical simulations of the model (1) is the amount of oil (percentage of domain) recovered at time $T$, given by $|\Omega|^{-1}\int_{\Omega}\phi c(x,T)$. Figure 27 shows the percentage of domain recovered at $t=10$ years, for each of the three numerical schemes on Cartesian meshes. These results were obtained starting with a very coarse mesh consisting of square cells of dimension $200\times 200$, while refining the spacial discretisation by a factor of 2, leading to a final mesh consisting of square cells of dimension $25\times 25$.

It can be seen here that each of the three schemes seem to yield numerical results that converge to different quantities. The MFV–upwind scheme probably overestimates the amount of oil recovered due to the diffusion it introduces, and thus a ”wider” region has been sweeped. On the other hand, it is likely that the MFEM–ELLAM underestimates the amount of oil recovered due to the presence of undershoots, as we have noted from Figures 23 and 24 (left). Thus, we expect the amount of oil recovered to be somewhere between these two values. Here, the solution obtained from HMM–ELLAM converges to such a value, which seems to sit about 65% for the first test case and 77% for the second test case.

As an element of comparison, we consider the results obtained in [1] on the same model with the Hybrid High-Order (HHO) method. This method is based on polynomials, with an arbitrary chosen degree, in the cells and on the faces, and can theoretically achieve any order of accuracy on any polytopal mesh (at an increased computational cost compared to HMM and ELLAM methods, of course); in practice, though, tests are usually ran up to order 4 or 5; we refer to [15] for the presentation on the pure diffusion equation, and to [13] for advection–diffusion–reaction models. The tests in [1] were performed up to order 4 and seem to indicate that the expected recovery after 10 years is around 65% for the homogeneous test case, and around 75% for the inhomogeneous test case. The low-order, less expensive HMM–ELLAM method seems to provide very similar results, on the contrary to MFV–ELLAM and MFEM–ELLAM. The latter, in particular, only predicts a 50% recovery for the inhomogeneous test case, which is much lower than all other methods; this is naturally expected because of the presence of severe undershoots. The concentration profiles obtained from the HMM–ELLAM are also similar to the one obtained through the HHO in [1].

Moreover, due to the thinner transition layer present in the HMM–ELLAM, this solution is preferred over the MFV-upwind scheme. On the other hand, since the solutions between the HMM–ELLAM and MFV-upwind are quite close to each other for the second test case, one might argue that MFV-upwind would be preferable, due to the cheaper computational cost associated with it. This is not necessarily the case, as we note that the HMM–ELLAM already generates an acceptable concentration profile, as well as a good approximation of the amount of oil recovered upon taking a mesh with square cells of dimension $50\times 50$ ft, whereas MFV-upwind needs to go to the refined mesh with square cells of dimension $25\times 25$ ft (which requires solving a much larger system) before producing such results.

4. Conclusion

Most of the previous work focused on separately developing schemes for either solutions of anisotropic diffusion equations (corresponding to the pressure equation (1a)) or schemes for advection dominated equations (corresponding to the concentration equation (1b)). This work presents a complete scheme for the coupled system, with details on how to combine the scheme for the pressure equation (HMM) with the scheme for the concentration equation (ELLAM). Moreover, this is usable on generic meshes as encountered in real-world applications – with the usual caveats on distorted meshes.

A new method has been introduced to reconstruct Darcy velocities from internal numerical fluxes, based on a consistency relation. This method turns out to be more cost efficient, and more efficient on distorted meshes than the one presented in [31]. Its construction in the 2-dimensional case seems promising enough to envision developing a 3-dimensional version, which is the purpose of ongoing work.

Our analysis also demonstrates the importance of selecting the correct number of approximation points – depending on the regularity of the mesh – to track the cells. An improved treatment of the cells near injection wells has also been introduced, leading to less overshoots and better approximations.

The streamlines in Section 3.5 show that velocities reconstructed from a piecewise constant approximation of the edge fluxes exhibit grid effects, especially on distorted grids. A natural direction for the improvement of the scheme, currently being explored, is to solve the pressure equation (1a) with a slightly higher order method. High order schemes for anisotropic diffusion equations have been developed for generic meshes, and a generic framework which covers a number of these schemes has been presented in [14]. It is therefore possible, on distorted meshes, to implement one of these high order methods for the pressure equation only, with the hope that it provides more accurate fluxes than the HMM to be used in the reconstruction of the Darcy velocity for the ELLAM component on the concentration equation. Since the high order scheme would only be used on the pressure equation, the overall computational cost would be much reduced compared to the full HHO method in [1].

Upon comparison with MFV–upwinding and MFEM–ELLAM schemes, we also got to see that the HMM–ELLAM captures the concentration profile better than upwinding, although as compared to MFEM–ELLAM we have a larger transition layer, as presented in the contour plots in Figures 21 to 26 (right). As for the amount of oil recovered, it seems that the HMM–ELLAM also performs better than both MFEM–ELLAM and MFV–upwind schemes for both test cases. Hence, overall, an HMM–ELLAM scheme is preferred over MFEM–ELLAM or MFV-upwind schemes.

Acknowledgement: this research was supported by the Australian Government through the Australian Research Council’s Discovery Projects funding scheme (project number DP170100605).