# Numerical convergence of the family of flux-continuous schemes with variable quadrature (qu1,qu2) for single phase flow in porous media

## Mayur Pal1

1Faculty of Mathematics and Natural Sciences, Kaunas University of Technology, Kaunas, Lithuania

Mathematical Models in Engineering, Vol. 8, Issue 2, 2022, p. 26-41. https://doi.org/10.21595/mme.2022.22636
Received 24 April 2022; received in revised form 18 May 2022; accepted 3 June 2022; published 30 June 2022

Copyright © 2022 Mayur Pal. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Editor's Pick
Views 31
CrossRef Citations 0
Abstract.

Finite-volume schemes, which honor pressure and flux-continuity conditions, is developed using double quadrature $\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$, referred as double family scheme. The scheme is applicable to solve the elliptic pressure equation used in reservoir simulation. Schemes are applicable on both regular cartesian and unstructured triangular meshes. The scheme is defined over a control-volume distributed formulation. The scheme can be applied to both diagonal and full permeability tensor elliptic pressure equation with discontinuous coefficients. The scheme removes the first order errors, which are introduced by standard reservoir simulation schemes when applied to full tensor flow. The scheme is quantified with help of a quadrature rule. When the scheme is applied to highly heterogeneous and anisotropic porous media it does not honor maximum principle resulting in unstable solution with oscillatory behavior. The numerical solution is termed non-monotonicity for high anisotropy ratios with results showing oscillations in the numerical pressure solution. In this paper a double $\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ quadrature flux continuous schemes is presented, which with specific choice of quadrature $\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ helps in improved stability of the numerical solutions. Numerical convergence of the scheme is also demonstrated with help of a number of numerical test cases and schemes impact on monotonicity behavior is also demonstrated with numerical examples.

Highlights
• Gridding
• Numerical Solution
• Numerical Convergence

Keywords: finite volume, flux and pressure continuity, monotonicity, stability, permeability anisotropy, elliptic pressure equation.

#### 1. Introduction

Underground porous media reservoirs of hydrogeological flows have complex description in terms of both their physical geometry and geology. Rapid variation in permeability with strong anisotropy is common occurrence in subsurface reservoir simulation, be it for hydrocarbon flows or for underground hydrological studies. Subsurface reservoir simulation requires discretization of continuous flux and pressure equation in order to honor correct local physical interface conditions between the simulation grid blocks, with strong discontinuities and anisotropy in subsurface properties, primarily permeabilities. The derivation of the algebraic flux continuity conditions for full permeability tensor discretization operators has lead to an efficient and robust locally conservative flux-continuous control volume distributed (CVD) finite-volume schemes for determining the discrete pressure and velocity fields in subsurface reservoirs [1-9]. These schemes are also known as multi-point flux approximation or MPFA . Further schemes of this type are presented in [11-13] and more recently in [14, 15]. These numerical schemes could be applied to the diagonal and full tensor elliptic pressure equation with rapid variations in permeability and remove the first order error in numerical approximation in typical reservoir simulation results, which are often based on two-point flux approximation, when applied to full tensor flow approximation. Some other numerical schemes of similar types that also preserve continuity of fluxes have also been developed using mixed finite element methods [16-18] and discontinuous galerkin finite element methods [19, 20]. Numerical convergence of these numerical scheme on unstructured and structured meshes has also been presented earlier in a number of papers [9, 21-23].

When these numerical schemes are applied to strongly anisotropic heterogeneous porous media they fail to satisfy a maximum principle (similar to finite eminent, mixed finite element and other finite-volume methods) and result in loss of solution stability, particularly, for very high anisotropy ratios, which often lead to spurious oscillations in the numerical pressure solution. Some work has been carried out previously with aim of preserving stability of the numerical solutions for highly heterogeneous cases with skewed meshes, please see references [3, 6, 24].

Monotonicity conditions were first derived in  for an M-matrix and later in  for a monotone matrix. Numerical treatment based on mesh/grid optimization methods have been used to help improve stability of the discrete system . Still, all of the published numerical schemes show only a limited impact and their stability is often subject to the order of magnitude of the permeability anisotropy ratios. A more promising strategy based on Flux-splitting techniques has been presented in [7, 8, 27], which helps in locally imposing maximum principle and resulting numerical solutions are oscillation free. Therefore, guaranteeing a more stable numerical pressure solution.

Conditions for the (one parameter) family of flux-continuous schemes to possess symmetric positive definite matrices and for such schemes to have M-matrices with diagonal dominance and negative off-diagonals are presented in [1, 2]. The discretized system obtained for the schemes in the case of a full tensor are shown to be conditionally diagonally dominant and thus the M-matrix property is conditional [1, 28]. For very large anisotropic ratios with skewed meshes, the resulting numerical system for these schemes can be unstable (as in the case of other numerical techniques like FEM & FVM) and the numerical pressure solutions can consequently exhibit oscillatory behavior. This paper is focused on seeking the most robust numerical approximation with help of double parameter $\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ scheme. Numerical convergence of the double parameter schemes is demonstrated with help of series of numerical test cases.

This paper is organized as follows: Description of the single phase flow problem encountered in reservoir simulation with respect to the general tensor pressure equation is given in Section 2. Details of the numerical construction of the double parameter scheme with discretization in physical space are presented in Section 3. Section 4 details the precise conditions required to obtain a stable numerical solution and shows the performance of the scheme with respect to stability conditions. Section 5 presents the numerical convergence results with help of a series of test cases. Section 6 presents numerical examples that demonstrate the impact of scheme quadrature choice on monotonicity behavior with help of some numerical examples. Finally, conclusions are presented in Section 7.

#### 2. Problem description and governing equations

The typical formulation of Darcy’s law states that the problem is to find the pressure $p$ satisfying:

(1)
$-{\int }_{\mathrm{\Omega }}\mathrm{‍}\nabla •\mathbf{K}\left(x,y\right)\nabla pd\tau ={\int }_{\mathrm{\Omega }}\mathrm{‍}qd\tau =\mathbf{M},$

over a given domain $\mathrm{\Omega }$, subjected to a given (Neumann/Dirichlet) boundary conditions on a given boundary $\partial \mathrm{\Omega }$. The right hand side term $\mathbf{M}$ represents a flow rate and . The permeability Matrix $\mathbf{K}$ can either be a diagonal or full cartesian tensor, which has the general form:

(2)
$\mathbf{K}=\left(\begin{array}{ll}{k}_{11}& {k}_{12}\\ {k}_{12}& {k}_{22}\end{array}\right),$

The resulting pressure equation with full permeability matrix is assumed to be elliptic such that:

(3)
${K}_{12}^{2}\le {k}_{11}{k}_{22}.$

The permeability tensor can be discontinuous across internal boundaries of the domain $\mathrm{\Omega }$. The boundary conditions, which are imposed here are of Dirichlet and Neumann type. For incompressible fluid flow in porous media – pressure need to be specified at at-least at one point in the domain to obtain a numerical solutions. For subsurface reservoir simulation, Neumann boundary conditions on $\partial \mathrm{\Omega }$ requires no-flux on the domain boundary such that $\left(K\nabla p\right)\cdot \stackrel{^}{n}=0$, where $\stackrel{^}{n}$ is the outward normal vector to $\partial \mathrm{\Omega }$.

#### 2.1. Elliptic pressure equation

The pressure equation presented in the section above is given in the physical space w.r.t the classical Cartesian coordinate system. In this section we present the elliptic pressure equation in a general curvilinear coordinate system that is defined with respect to a uniform dimensionless transform space with a $\left(\xi ,\eta \right)$ coordinate system. Choosing domain ${\mathrm{\Omega }}_{p}$ to represent a control volume comprised of surfaces that are tangential to constant $\left(\xi ,\eta \right)$ respectively, Eq. (1) is integrated over ${\mathrm{\Omega }}_{p}$ via the Gauss divergence theorem to yield:

(4)
$-{\oint }_{{\partial }_{{\mathrm{\Omega }}_{P}}}\mathrm{‍}\left(\mathbf{K}\nabla p\right)\cdot \stackrel{^}{\mathbf{n}}ds=\mathbf{M},$

where $\partial {\mathrm{\Omega }}_{p}$ is the boundary of ${\mathrm{\Omega }}_{p}$ and $\stackrel{^}{n}$ is the unit outward normal. Spatial derivatives are computed using:

(5)
${p}_{x}=J\left(p,y\right)/J\left(x,y\right),{p}_{y}=J\left(x,p\right)/J\left(x,y\right),$

where $J\left(x,y\right)={x}_{\xi }{y}_{\eta }-{x}_{\eta }{y}_{\xi }$ is the Jacobian. Resolving the $x$, $y$ components of velocity along the unit normals to the curvilinear coordinates $\left(\xi ,\eta \right)$, e.g., for $\xi =$ constant, $\stackrel{^}{\mathbf{n}}ds=\left({y}_{\eta },-{x}_{\eta }\right)d\eta$ gives rise to the general tensor flux components:

(6)

where general tensor $\mathbf{T}$ has elements defined by:

(7)
${T}_{11}=\frac{{k}_{11}{y}_{\eta }^{2}+{k}_{22}{x}_{\eta }^{2}-2{k}_{12}{x}_{\eta }{y}_{\eta }}{J},$
${T}_{22}=\frac{{k}_{11}{y}_{\xi }^{2}+{k}_{22}{x}_{\xi }-2{k}_{12}{x}_{\xi }{y}_{\xi }}{J},$
${T}_{12}=\frac{{k}_{12}\left({x}_{\xi }{y}_{\eta }+{x}_{\eta }{y}_{\xi }\right)-\left({k}_{11}{y}_{\eta }{y}_{\xi }+{k}_{22}{x}_{\eta }{x}_{\xi }\right)}{J},$

and the closed integral can be written as:

(8)
$\int \mathrm{‍}{\int }_{{\mathrm{\Omega }}_{p}}\mathrm{‍}\frac{\left({\partial }_{\xi }\stackrel{~}{F}+{\partial }_{\eta }\stackrel{~}{G}\right)}{J}Jd\xi d\eta ={∆}_{\xi }F+{∆}_{\eta }G=m,$

where e.g. ${∆}_{\xi }F$ is the difference in net flux with respect to $\xi$ and $\stackrel{~}{F}={T}_{11}{p}_{\xi }+{T}_{12}{p}_{\eta }$, $\stackrel{~}{G}={T}_{12}{p}_{\xi }+{T}_{22}{p}_{\eta }$. This condition demonstrates that any such scheme applicable to a full matrix permeability tensor could also be applied to non-K-Orthogonal grids. Note that ${T}_{11}$, ${T}_{22}\ge 0$ and ellipticity of $\mathbf{T}$ follows from Eqs. (3 and 7). Full permeability tensors can arise from upscaling, and local orientation of the grid and the permeability field. For example, by Eq. (7), a diagonal anisotropic Cartesian tensor leads to a full tensor on a curvilinear orthogonal grid  (a grid is called orthogonal if all grid lines intersect at a right angle).

#### 3. Flux-continuous finite volume schemes

Finite volume schemes, which honor local pressure and flux continuity conditions and are locally conservative have been presented in [1-6]. The schemes have been developed for different grid types structured and unstructured grids in physical space and transform space. Numerical convergence for a range of quadrature in physical and transform space have been presented in . The schemes have continuity of pressure and flux across primal grid cell control-volume interfaces (Control Volume is defined as the volume over which mass can transfer in and out across its boundary (Pal 2006)). Here we present a short summary of the structured cell centred quadrilateral formulation. The nine point support of the numerical scheme is shown in Fig. 1(a)). The numerical scheme has flow and rock variables (porosity, permeability etc.) defined at the cell centre, with this definition, the discrete approximation points are at then located at the centres of the primal grid cells, which are also the control-volumes. For each group of four nodes, four triangles are drawn as in Fig. 1(a). Each group of nodes are defined within a dual-cell (each group of four cell-centered nodes surrounding a primal grid vertex defines the fundamental corners of a dual cell) which is obtained by joining cell centres with cell edge mid-points as indicated by the dashed contour in Fig. 1(b). The dual cells partition the primal cells (or control-volumes) into subcells (the dual cells partition the primal quadrilateral grid cells (or control volumes) into sub-quadrilateral cells, which are called subcells). Two faces of each subcell also define sub faces of two faces of the parent control-volume.

Fig. 1. a) Flux-continuity at specific north, south, east and west location $\left(n,s,e,w\right)$ on the control-volume sub-cell faces, b) location of quadrature points $\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ on the sub-cell face a) b)

#### 3.1. $\left({\mathbit{q}}_{{\mathbit{u}}_{1}},{\mathbit{q}}_{{\mathbit{u}}_{2}}\right)$ quadrature schemes

Flux continuous schemes are formulated using the principle of pressure and flux continuity conditions on the control-volume subcell faces. The continuity conditions are shown conceptually with help of the Fig. 1(a). It can be seen in the figure that the continuity conditions is met at the four positions $\left(n,s,e,w\right)$, corresponding to north, south, east and west faces. The area is shown with help of shaded triangles and their point of contact in the Fig. 1(b). On each control-volume cell sub-face the point of continuity is parameterized with help of the variable quadrature $q=\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ depending on the location of continuity at the sub-face, please see Fig. 1(b) $\left(0<{q}_{{u}_{1}},{q}_{{u}_{2}}\le 1\right]$. For any given control-volume subcell, the points of continuity of flux and pressure could be located anywhere in the ranges $\left(0<{q}_{{u}_{1}},{q}_{{u}_{2}}\le 1\right]$ on the two faces of each subcell. The precise values of $\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ define the quadrature points, thereby, resulting in the formulation corresponding to families schemes, which are flux-continuous. Fig. 1(b) shows a conceptual plot for quadrature points located at ${q}_{{u}_{1}}=0.01$ and ${q}_{{u}_{2}}=1$. Cell face pressures ${P}_{n}$, ${P}_{e}$, ${P}_{s}$, ${P}_{w}$ are introduced at $n$, $s$, $e$, $w$ locations. Pressure sub-triangles are defined with local triangular support imposed within each quarter (sub-cell) of the dual-cell as shown in Fig. 1(b). Pressure $p$, in local cell coordinates, is piecewise linear over each triangle. The physical space flux-continuity conditions for cells 1 to 4, sharing a common grid vertex inside the dual-cell are expressed as:

(9)
${F}_{N}=-\frac{1}{2}\left({T}_{11}{p}_{\xi }+{T}_{12}{p}_{\eta }\right){|}_{N}^{3}=-\frac{1}{2}\left({T}_{11}{p}_{\xi }+{T}_{12}{p}_{\eta }\right){|}_{N}^{4},$
${F}_{S}=-\frac{1}{2}\left({T}_{11}{p}_{\xi }+{T}_{12}{p}_{\eta }\right){|}_{S}^{1}=-\frac{1}{2}\left({T}_{11}{p}_{\xi }+{T}_{12}{p}_{\eta }\right){|}_{S}^{2},$
${F}_{E}=-\frac{1}{2}\left({T}_{12}{p}_{\xi }+{T}_{22}{p}_{\eta }\right){|}_{E}^{2}=-\frac{1}{2}\left({T}_{12}{p}_{\xi }+{T}_{22}{p}_{\eta }\right){|}_{E}^{3},$
${F}_{W}=-\frac{1}{2}\left({T}_{12}{p}_{\xi }+{T}_{22}{p}_{\eta }\right){|}_{W}^{1}=-\frac{1}{2}\left({T}_{12}{p}_{\xi }+{T}_{22}{p}_{\eta }\right){|}_{W}^{4}.$

The parametric variation in $q=\left({q}_{1},{q}_{2}\right)$ is illustrated further using the sub-cell example of Fig. 1(c), with sub-cell containing sub-triangle $\left(1,S,W\right)$. Let ${r}_{1}=\left({x}_{1},{y}_{1}\right)$ denote the coordinates of the cell-centre and ${r}_{S}=\left({x}_{S},{y}_{S}\right)$,${r}_{W}=\left({x}_{W},{y}_{W}\right)$ denote the local continuity coordinates. Then it is understood that the continuity position is a function of $q$ with ${r}_{S}\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ and ${r}_{W}\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$.

Piecewise constant fluxes are then computed on every pressure sub-triangles belonging to the sub-cells of the dual-cell as shown in Fig. 1(b). The local linear pressure $p$, is expanded in sub-triangle coordinates. The Darcy flux approximation for sub-triangle $\left(1,S,W\right)$ is given below:

(10)
$\left(\begin{array}{l}{p}_{\xi }\\ {p}_{\eta }\end{array}\right)=\left(\begin{array}{c}{p}_{S}-{p}_{1}\\ {p}_{W}-{p}_{1}\end{array}\right),$

and:

(11)

Using Eqs. (10, 11) the discrete Darcy velocity is defined as:

(12)
${v}_{h}=-\mathbf{K}\nabla {p}_{h}=-\mathbf{K}\mathbf{G}\left(q\right)\left(\begin{array}{l}{p}_{\xi }\\ {p}_{\eta }\end{array}\right),$

where $\mathbf{K}$ is the local permeability tensor of cell 1 and dependency of $\nabla {p}_{h}$ on quadrature point arises through:

(13)
$\mathbf{G}\left(q\right)\left(\begin{array}{l}{p}_{\xi }\\ {p}_{\eta }\end{array}\right)=\left(\begin{array}{cc}{y}_{\eta }\left(q\right)& -{y}_{\xi }\left(q\right)\\ -{x}_{\eta }\left(q\right)& {x}_{\xi }\left(q\right)\end{array}\right)1J\left(q\right)\left(\begin{array}{c}{p}_{S}-{p}_{1}\\ {p}_{W}-{p}_{1}\end{array}\right),$

where approximate ${r}_{\xi }\left(q\right)$ and ${r}_{\eta }\left(q\right)$ are defined by Eq. (11). The normal flux at the left hand side of $S$ (Fig. 1(b)) is resolved along the outward normal vector $d{L}_{S}=\left(\mathrm{\Delta }{y}_{{r}_{3},S},-\mathrm{\Delta }{x}_{{r}_{3},S}\right)=\frac{1}{2}\left(\mathrm{\Delta }{y}_{32},-\mathrm{\Delta }{x}_{32}\right)$ (Fig. 1(b) and is defined in terms of the general tensor $T=T\left(q\right)$ as:

(14)
${F}_{S}^{1}={v}_{h}\cdot d{L}_{S}=-\frac{1}{2}\left({T}_{11}^{1}{p}_{\xi }+{T}_{12}^{1}{p}_{\eta }\right){|}_{S}^{1},$

and the resulting coefficients of $-\frac{1}{2}\left({p}_{\xi },{p}_{\eta }\right){|}_{S}^{1}$ denoted by ${T}_{11}{|}_{S}^{1}$ and ${T}_{12}{|}_{S}^{2}$ are sub-cell (physical-space) approximations of the general tensor components (Eq. (9)) at the left hand face at S, and are functions of $q$. A similar expression for flux is obtained at the right hand side of S from cell 2 (Fig. 1(b). Similarly sub-cell fluxes are resolved on the two sides of the other faces at W, N and E. Flux continuity is then imposed across the four cell interfaces at the four positions N, S, E and W (Fig. 1(a) which are specified according to quadrature point $\mathbf{q}$. The local physical space flux continuity conditions are now defined in the dual cell and expressed as follows.

The above system of Eq. (9) is now expressed as:

(15)
$F={A}_{L}{P}_{f}+{B}_{L}{P}_{v}={A}_{R}{P}_{f}+{B}_{R}{P}_{v},$

where $F=\left({F}_{N},{F}_{S},{F}_{E},{F}_{W}{\right)}^{T}$ are the fluxes defined in the dual-cell and ${P}_{f}=\left({p}_{N},{p}_{S},{p}_{E},{p}_{W}{\right)}^{T}$ are the interface pressures. Similarly ${P}_{v}=\left({p}_{1},{p}_{2},{p}_{3},{p}_{4}{\right)}^{T}$ are the cell centered pressures. Thus the four interface pressures are expressed in terms of the four cell centered pressures. Using Eq. (15), ${P}_{f}$ is now expressed in terms of ${P}_{v}$ to obtain the dual-cell flux and coefficient matrix:

(16)
$F=\left({A}_{L}\left({A}_{L}-{A}_{R}{\right)}^{-1}\left({B}_{R}-{B}_{L}\right)+{B}_{L}\right){P}_{v}.$

Thus, the cell-face pressures are eliminated from the flux by being determined locally in terms of the cell centered pressures in a preprocessing step, avoiding introduction of the interface pressure equations into the assembled discretization matrix. The Eq. (16) can also be written as:

(17)
$AF=-\mathrm{\Delta }{P}_{v},$

where the entries of matrix $A$ are accumulated inverse tensor elements and $\mathrm{\Delta }{P}_{v}=\left({p}_{21},{p}_{32},{p}_{34},{p}_{41}{\right)}^{T}$ are the differences of vertex pressures. Consistency of the formulation follows from Eq. (17) which shows that flux is zero for constant potential. The relationship with the mixed method can also be deduced from Eq. (17),  from which it also follows that the above systems Eqs. (16), (17) represent the generalisation of the standard flux with harmonic coefficients to general elements with families of schemes defined by quadrature point $q$, see [2, 4, 9, 21-23] for details. The physical space formulation does not posses a symmetric discretization matrix for arbitrary quadrilaterals, however, transform space (cell and sub-cell) formulations that are symmetric positive definite are presented in [1, 4]. Flux continuity in the case of a general-tensor is obtained while maintaining the standard single degree of freedom per cell. Since the continuity equations depend on both ${p}_{\xi }$ and ${p}_{\eta }$ (unless a diagonal tensor is assumed with cell-face midpoint quadrature resulting in a 2-point flux), the interface pressures ${p}_{f}=\left({p}_{N},{p}_{S},{p}_{E},{p}_{W}{\right)}^{T}$ are locally coupled and each group of four interface pressures is determined simultaneously in terms of the group of four cell centered pressures whose union contains the continuity positions. Finally for a structured grid the scheme is defined by:

(18)
${F}_{i+1/2,j}-{F}_{i-\frac{1}{2},j}+{F}_{i,j+\frac{1}{2}}-{F}_{i,j-\frac{1}{2}}=M,$

where $i$, $j$ are the integer coordinates of the central quadrilateral cell, Fig. 1(a)) and:

(19)
${F}_{i+1/2,j}={F}_{{S}_{i+1/2,j+1/2}}+{F}_{{N}_{i+1/2,j-1/2}},$
${F}_{i,j+1/2}={F}_{{E}_{i-1/2,j+1/2}}+{F}_{{W}_{i+1/2,j+1/2}},$

where $i+1/2$, $j+1/2$ denote the “integer” coordinates of the top right hand side dual-cell, Fig. 1(a)). The unstructured formulation is presented in [22, 23].

#### 4. Monotonicity

Flux-continuous schemes results in a discrete matrix which has five to nine diagonals in two dimensions and three to twenty-seven diagonals in three dimensions on a structured mesh. The discrete matrix can be expressed as:

(20)
$\mathbf{A}p=l,$

where $\mathbf{A}$ is the matrix of discrete operator, $p$ is the unknown pressure and $l$ is the source term. Ideally the discrete system obtained in this way should results in an Eq. (20), which should be a nonsingular matrix and should satisfy maximum principle that is comparable to the continuous counterpart of the discrete problem. The condition ensures that the numerical solution obtained in this way is free from any non-physical oscillations. The maximum principle is satisfied for the discrete system given in 20 if the discrete matrix operator $\mathbf{A}$ is a monotone. The matrix operator, of the discrete system, $\mathbf{A}$ is a non-singular monotone matrix if it obeys the following condition :

(21)
${\mathbf{A}}^{-1}\ge \mathbf{O},$

where $\mathbf{O}$ is a null matrix. A sufficient condition for a maximum principle that can ensure that non-physical oscillations do not occur in the resulting discrete solution is that $\mathbf{A}$ is an M-matrix, i.e. monotone or positive definite with ${a}_{i,j}\le 0$, $i\ne j$.

#### 4.1. M-matrix conditions

The following set of conditions are often easier to verify.

$\mathbf{A}$ is an M-matrix if:

(22)
${a}_{i,i}>0,\forall i,$
${a}_{i,j}\le 0,\forall i,j,i\ne j,$
${\mathrm{\Sigma }}_{j}{a}_{i,j}\ge 0,\forall i.$

In addition, $\mathbf{A}$ must be either strictly diagonally dominant (strict inequality in the latter of Eq. (22) or weakly diagonally dominant with strict inequality for at least one row, $\mathbf{A}$ must also be irreducible. The first M-matrix analysis for schemes of this type is presented in , where conditions for nine-node flux continuous schemes to be an M-Matrix are:

(23)
$min\left({T}_{1,1},{T}_{2,2}\right)\ge \eta \left(q\right)\left({T}_{1,1}+{T}_{2,2}\right)\ge \left|{T}_{1,2}\right|,$

and $\eta \left(q\right)$ is a function of quadrature point. One of the essential conditions here is that:

(24)
$\left|{T}_{1,2}\right|\le min\left({T}_{1,1},{T}_{2,2}\right),$

which is only sufficient for ellipticity  and therefore quite limiting on the range of tensors that are applicable. Tensors that are elliptic with:

(25)
${T}_{1,2}^{2}\le {T}_{1,1}{T}_{2,2},$

and are such that $|{T}_{1,2}|>min\left({T}_{1,1},{T}_{2,2}\right)$ violate the M-Matrix criteria of Eq. (23) and expose the M-Matrix limit.

#### 5. Numerical results: numerical convergence study

While numerical convergence tests for transform and physical space formulations on cartesian and triangular unstructured meshes in two-and three dimensions have previously been performed by [5, 21-23] for the default member of the family of flux-continuous MPFA schemes i.e. for 1.

The aim of the numerical convergence study presented here is to test the effect of variable quadrature points ${q}_{1}$, ${q}_{2}$ on the numerical convergence of the scheme. Although a number of quadrature combinations exist for ${q}_{1}$, ${q}_{2}$, only two specific combination of quadrature points are evaluated ${q}_{{u}_{1}}=$ 0.5, ${q}_{{u}_{2}}=$ 0.1 and ${q}_{{u}_{1}}=$ 0.1, ${q}_{{u}_{2}}=$ 0.5.

In this section numerical convergence results are presented cartesian grids in two dimensions. In all test cases the permeability field remains unchanged under grid refinement, ensuring that each problem is invariant with respect to each grid level for the convergence study. The classical ${L}_{2}$ norm is used to measure errors in fluxes and pressure, which is defined for fluxes and pressure as:

(26)
$||{p}_{h}-p|{|}_{{L}_{2}}={\left(\frac{\sum _{i}\mathrm{‍}\left({V}_{i}\left({p}_{h,i}-{p}_{i}{\right)}^{2}\right)}{\sum _{i}\mathrm{‍}{V}_{i}}\right)}^{\frac{1}{2}},$
(27)
$||{f}_{h}-f|{|}_{{L}_{2}}={\left(\frac{\sum _{j}\mathrm{‍}\left({Q}_{j}\left({f}_{h,j}-{f}_{j}{\right)}^{2}\right)}{\sum _{j}\mathrm{‍}{Q}_{j}}\right)}^{\frac{1}{2}},$

where, $f=-\mathbf{v}\cdot \mathbf{n}$ (where $\mathbf{v}=\mathbf{K}\nabla p$) is the edge normal flow velocity. Subscript $h$ refers to numerical solution. Further ${V}_{i}$ is the volume of the grid cell $i$, and ${Q}_{j}$ is the volume associated with edge $j$ (where two cells are separated by edge $j$). The grid refinement levels used for the ${L}_{2}$ norm calculation where mesh sizes of 16×16, 32×32, 64×64, and 128×128 are used for all test cases in 2D. In each case Dirichlet boundary conditions are prescribed via the exact solutions.

#### 5.1. Numerical test 1: uniformly discontinuous permeability

In this case the pressure field is piecewise quadratically varying over the domain shown in Fig. 2(a). Fig. 2(b) shows analytical pressure results on this test case. Fig. 2(c) shows typical mesh used for the numerical convergence in this case. The domain is discontinuous with discontinuity aligned along the line $x=1/2$, and the analytical solution is given as:

(28)
$p\left(x,y\right)=\left\{\begin{array}{ll}{c}_{l}{x}^{2}+{d}_{l}{y}^{2},& x<1/2,\\ {a}_{r}+{b}_{r}x+{c}_{r}{x}^{2}+{d}_{r}{y}^{2},& x\ge 1/2,\end{array}\right\$
$K=\left\{\begin{array}{ll}\left(\begin{array}{ll}50& 0\\ 0& 1\end{array}\right),& x<1/2,\\ \left(\begin{array}{ll}1& 0\\ 0& 10\end{array}\right),& x\ge 1/2,\end{array}\right\$

The imposed top boundary flux is also discontinuous at $x=1/2$, resulting in a discontinuous tangential flux across the domain. The computed numerical solution and the plots showing ${L}_{2}$ norm of pressure for the varying quadrature ${q}_{1}$ and ${q}_{2}$ are shown in Fig. 3(a) and (b), respectively. ${O}^{{h}^{1.0179}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=\text{0.5}$ and ${q}_{{u}_{2}}=\text{0.1}$ and ${O}^{{h}^{1.0496}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=\text{0.1}$ and ${q}_{{u}_{2}}=\text{0.5}$.

Fig. 2. Case 1: a) Discontinuity of the domain under consideration, b) exact analytical solution for pressure, c) sample mesh used for numerical convergence a) b) c)

Fig. 3. Case 1: a) convergence plots for ${q}_{{u}_{1}}$, ${q}_{{u}_{2}}$. b) exact numerical pressure – physical space a) b)

#### 5.2. Numerical test 2: non-uniform discontinuous permeability field

In this case the medium is divided into two parts as, see Fig. 4(a). The permeability is discontinuous and permeability ratio is 1/10 across the medium discontinuity. The discontinuity is aligned along the line $rx+sy=0$, where $r=\mathrm{t}\mathrm{a}\mathrm{n}\left(\pi /3\right)/\left(1+\mathrm{t}\mathrm{a}\mathrm{n}\left(\pi /3\right)\right)$ and $s=1/\left(1+\mathrm{t}\mathrm{a}\mathrm{n}\left(\pi /3\right)\right)$.The pressure is piecewise linear and varies as:

(29)
$p\left(x,y\right)=\left\{\begin{array}{ll}rx+sy,& rx+sy<0,\\ 10\left(rx+sy\right),& rx+sy>0.\end{array}\right\$

The diagonal permeability tensor $\mathbf{K}=c\mathbf{I}$, where $c=10$ for $rx+sy<0$ and $c=1$ for $rx+sy>0$. The numerical solution shown in Fig. 4(b) was obtained using a grid aligned along the discontinuity. Typical mesh used for the numerical convergence is shown in Fig. 4(c). Numerical convergence for varying ${q}_{{u}_{1}}$ and ${q}_{{u}_{2}}$ quadrature is shown in Fig. 5(a) and the numerical pressure Physical space for the varying quadrature is shown in Fig. 5(b). ${O}^{{h}^{2.0218}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=0.5$ and ${q}_{{u}_{2}}=0.1$ and ${O}^{{h}^{1.9949}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=0.1$ and ${q}_{{u}_{2}}=0.5$.

Fig. 4. Case 1: a) discontinuity of the domain under consideration, b) exact analytical solution for pressure, c) sample mesh used for numerical convergence a) b) c)

Fig. 5. Case 1: a) convergence plots for ${q}_{{u}_{1}}$, ${q}_{{u}_{2}}$, b) Exact numerical pressure – physical space a) b)

#### 5.3. Numerical test 3: permeability field with internal discontinuity

This case involves an internal discontinuity. The problem involves a rectangular domain with internal discontinuous permeability variation as indicated in Fig. 6(a). The exact solution in each case takes the form:

(30)
$p\left(r,\theta \right)={r}^{\alpha }\left({a}_{i}\mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \theta \right)+{b}_{i}\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \theta \right)\right).$

The exact analytical pressure solution is shown in Fig. 6(b) and the cartesian mesh used in this case is shown in Fig. 6(c). Numerical convergence for varying ${q}_{{u}_{1}}$ and ${q}_{{u}_{2}}$ quadrature is shown in Fig. 7(a) and the numerical pressure Physical space for the varying quadrature is shown in Fig.7(b). ${O}^{{h}^{0.9955}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=\text{0.5}$ and ${q}_{{u}_{2}}=\text{0.1}$ and ${O}^{{h}^{0.9859}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=\text{0.1}$ and ${q}_{{u}_{2}}=\text{0.5}$.

Fig. 6. a) Internal discontinuity along $\theta =\pi /2$, b) exact analytical pressure, c) sample mesh used for numerical convergence a) b) c)

Fig. 7. Case 1: a) convergence plots for ${q}_{1}$, ${q}_{2}$, b) exact numerical pressure – physical space a) b)

#### 6. Numerical results: test of monotonicity behavior

In this section a series of numerical test are presented to test the quasi-monotonic property of the variable quadrature points ${q}_{{u}_{1}}\text{,}$${q}_{{u}_{2}}$ scheme. The test examples involve subsurface permeability discontinuity with highly heterogeneous and anisotropic permeability tensor.

#### 6.1. Test case 1: v shape permeability distribution

In this test case, a two dimensional domain with a permeability tensor that changes direction in anisotropy halfway across the domain going from $\mathbf{K}$= [1, 0.99, 0.99, 1] to $\mathbf{K}$ = [1, –0.99, –0.99, 1], as shown in Fig. 8, is investigated. A point source is placed in the centre of the 2D domain (of unit length in $x$ and $y$ direction) and Dirichlet conditions are applied. In this test case due to discontinuity in permeability tensor M-matrix conditions are tested under variable coefficient conditions. Fig. 9(a) shows the results with ${q}_{{u}_{1}}={q}_{{u}_{2}}=q=1$ with visible spurious oscillations. In the second case using the different choice of quadrature parameter with ${q}_{{u}_{1}}={q}_{{u}_{2}}=q=$ 0.01 oscillation free results are obtained, see Fig. 9(b).

Fig. 8. Change in direction of anisotropy in Permeability Tensor at $y=\text{0.5}$ Fig. 9. a) Numerical pressure contours with instabilities like oscillations ${q}_{{u}_{1}}={q}_{{u}_{2}}=q=\text{1}$, b) numerical pressure contours free of instabilities e.g., no oscillation. ${q}_{{u}_{1}}={q}_{{u}_{2}}=q=\text{0.01}$ a) b)

Results using the double family with ${q}_{{u}_{1}}=$ 0.005025125, ${q}_{{u}_{2}}=\text{1}$ are shown in Fig. 10. In these cases the M-matrix conditions are satisfied. Using a reduced support scheme (7-point – constant tensor) motivated by the above analysis yields the result in Fig. 11(a). If the alternative orientation is selected the scheme still has reduced support and the scheme does not posses an M-matrix and the result obtained is shown in Fig. 11(b). There are clear oscillations and spread in the solution in this case, which serves as a test of the angled condition.

#### 6.2. Test case 2: strong anisotropic permeability distribution

In this case the same boundary conditions apply as in case 2. The permeability tensor changes direction in anisotropy halfway across the domain as before, now with stronger tensor field, going from $\mathbf{K}=$ [750.25 432.58 432.58 250.75] to [750.25 –432.58 –432.58 250.75] as indicated in Fig. 8. In this case the tensor has a principal anisotropy ratio of 1/1000, while the tensor is elliptic, the sufficient condition for ellipticity, is clearly violated. Consequently, the conditions for an M-matrix are also violated.

Fig. 10. Numerical pressure contours free of instabilities like oscillations for double quadrature ${q}_{{u}_{1}}=$ 0.005025125, ${q}_{{u}_{2}}=$ 1 Fig. 11. a) Numerical pressure 7-pt scheme tensor not aligned with scheme stencil, b) numerical pressure contours 7-pt scheme tensor aligned with scheme stencil a) b)

Fig. 12. a) Numerical pressure contours 9-pt scheme tough tensor, $q=\text{1}$ Results are presented for ${q}_{{u}_{1}}=\text{1}$ and ${q}_{{u}_{2}}=\text{1}$ Fig. 12 for a 64×64 grid. There are strong oscillations in the solution showing clear violation of the maximum principle. Next, we turn to the reduced support scheme with angled approximation with orientation of opposite sign to the cross terms, results are in Fig. 13(a) which show severe oscillations. Finally, we employ the reduced support scheme with angled approximation and orientation that matches the sign of the cross terms, consistent with the M-matrix analysis, for the same 64×64 grid Fig. 13(b). While oscillations are present, they are seen to be reduced. The results clearly show a reduction in oscillations in each case, however at each grid level and as convergence is approached, the consistent reduced support scheme yields almost oscillation free results or quasi-monotonic.

Fig. 13. a) Numerical pressure contours 7-pt scheme tough tensor not aligned with scheme stencil, b) numerical pressure 7-pt scheme tough tensor aligned with scheme stencil a) b)

#### 7. Conclusions

A two parameter $\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ scheme is presented, which is locally conservative and grantees continuity of flux and pressure. Numerical convergence of the scheme is also presented, which shows numerical accuracy of order ${O}^{h}-{O}^{{h}^{2}}$ on the tested examples.

When a classical scheme is used that violates the governing conditions oscillations appear in the solutions. In this paper results are presented for highly anisotropic full-tensor fields that verify that double family scheme offers quasi-monotonic behavior. When the governing conditions are satisfied the for given choice of quadrature the numerical pressure solution is free oscillations.

The variable support $\left({q}_{{u}_{1}},{q}_{{u}_{2}}\right)$ schemes are able to compute numerical pressure solutions with minimal or no oscillations and has much improved resolution compared to schemes relying on fixed quadrature rules. The scheme also proves effective for high anisotropy ratios.

1. M. G. Edwards and C. F. Rogers, “Finite volume discretization with imposed flux continuity for the general tensor pressure equation,” Computational Geosciences, Vol. 2, No. 4, pp. 259–290, 1998, https://doi.org/10.1023/a:1011510505406 [Publisher]
2. M. G. Edwards, “Unstructured, control-volume distributed, full-tensor finite-volume schemes with flow based grids,” Computational Geosciences, Vol. 6, No. 3/4, pp. 433–452, 2002, https://doi.org/10.1023/a:1021243231313 [Publisher]
3. M. G. Edwards, “M-Matrix flux splitting for general full tensor discretization operators on structured and unstructured grids,” Journal of Computational Physics, Vol. 160, No. 1, pp. 1–28, May 2000, https://doi.org/10.1006/jcph.2000.6418 [Publisher]
4. M. G. Edwards and M. Pal, “Symmetric positive definite general tensor discretization: a family of sub-cell flux continuous schemes on cell centred quadrilateral grids,” (in press). [Search CrossRef]
5. M. Pal, M. G. Edwards, and A. R. Lamb, “Convergence study of a family of flux‐continuous, finite‐volume schemes for the general tensor pressure equation,” (in Press), International Journal for Numerical Methods in Fluids, Vol. 51, No. 9-10, pp. 1177–1203, Jul. 2006, https://doi.org/10.1002/fld.1211 [Publisher]
6. Mayur Pal and Michael G. Edwards, “Effective upscaling using a family of flux-continuous finite-volume schemes for the pressure equation,” in ACME UK, Apr. 2006. [Search CrossRef]
7. M. Pal and M. G. Edwards, “Family of flux-continuous finite-volume schemes with improved monotonicity,” in Proceedings of the 10th European Conference on the Mathematics of Oil Recovery, 2006. [Search CrossRef]
8. M. Pal and M. G. Edwards, “Flux-Splitting Schemes for improved montonicity of discrete solution of elliptic equation with highly anisotropic coefficients,” in ECCOMAS CFD 2006: Proceedings of the European Conference on Computational Fluid Dynamics, 2006. [Search CrossRef]
9. M. Pal and M. G. Edwards, “Non‐linear flux‐splitting schemes with imposed discrete maximum principle for elliptic equations with highly anisotropic coefficients,” International Journal for Numerical Methods in Fluids, Vol. 66, No. 3, pp. 299–323, May 2011, https://doi.org/10.1002/fld.2258 [Publisher]
10. I. Aavatsmark, “Introduction to multipoint flux approximation for quadrilateral grids,” Computational Geosciences, Vol. 6, No. 3/4, pp. 405–432, 2002, https://doi.org/10.1023/a:1021291114475 [Publisher]
11. S. H. Lee, P. Jenny, and H. A. Tchelepi, “A finite-volume method with hexahedral multiblock grids for modeling flow in porous media,” Computational Geosciences, Vol. 6, No. 3/4, pp. 353–379, 2002, https://doi.org/10.1023/a:1021287013566 [Publisher]
12. S. Verma, “Flexible grids for reservoir simulation,” Ph.D. Thesis, Stanford University, 1996. [Search CrossRef]
13. S. K. Khattri, “Analyzing finite-volume for single-phase flow in porous media,” Journal of Porous Media, Vol. 10, No. 2, pp. 109–123, 2007. [Publisher]
14. O. Iliev and I. Rybak, On Approximation Property of Multipoint Flux Approximation Method. 2007. [Search CrossRef]
15. A. A. Lukyanov, “A meshless multi-point flux approximation method,” in Numerical Analysis and Applied Mathematics ICNAAM 2012: International Conference of Numerical Analysis and Applied Mathematics, 2012, https://doi.org/10.1063/1.4756377 [Publisher]
16. L. J. Durlofsky, “A Triangle Based Mixed Finite Element-Finite Volume Technique for Modeling Two Phase Flow through Porous Media,” Journal of Computational Physics, Vol. 105, No. 2, pp. 252–266, Apr. 1993, https://doi.org/10.1006/jcph.1993.1072 [Publisher]
17. A. S. Abushaikha, D. V. Voskov, and H. A. Tchelepi, “Fully implicit mixed-hybrid finite-element discretization for general purpose subsurface reservoir simulation,” Journal of Computational Physics, Vol. 346, pp. 514–538, Oct. 2017, https://doi.org/10.1016/j.jcp.2017.06.034 [Publisher]
18. A. S. Abushaikha and K. M. Terekhov, “A fully implicit mimetic finite difference scheme for general purpose subsurface reservoir simulation with full tensor permeability,” Journal of Computational Physics, Vol. 406, p. 109194, Apr. 2020, https://doi.org/10.1016/j.jcp.2019.109194 [Publisher]
19. B. Rivière, M. F. Wheeler, and K. Banaś, “Discontinuous Glaerkin method applied to a single phase flow in porous media,” Computational Geosciences, Vol. 4, No. 4, pp. 337–349, 2000, https://doi.org/10.1023/a:1011546411957 [Publisher]
20. T. F. Russell and M. Wheeler, “Mathematics of reservoir simulation: finite element and finite difference methods for continuous flows in porous media,” in Frontiers in Applied Mathematics SIAM, 1983, p. 35. [Publisher]
21. M. Pal and M. G. Edwards, “The competing effects of discretization and upscaling – a study using the q-family of CVD-MPFA,” 11th European Conference on the Mathematics of Oil Recovery, Vol. 68, No. 1, Sep. 2008, https://doi.org/10.3997/2214-4609.20146372 [Publisher]
22. M. Pal and M. G. Edwards, “A family of multi-point flux approximation schemes for general element types in two and three dimensions with convergence performance,” International Journal for Numerical Methods in Fluids, Vol. 69, No. 11, pp. 1797–1817, Aug. 2012, https://doi.org/10.1002/fld.2665 [Publisher]
23. Mayur Pal, “Families of control-volume distributed CVD(MPFA) finite volume schemes for the porous medium pressure equation on structured and unstructured grids,” Ph.D. thesis, Swansea University, 2007. [Search CrossRef]
24. J. M. Nordbotten and G. T. Eigestad, “Discretization on quadrilateral grids with improved monotonicity properties,” Journal of Computational Physics, Vol. 203, No. 2, pp. 744–760, Mar. 2005, https://doi.org/10.1016/j.jcp.2004.10.002 [Publisher]
25. J. Nordbotten and I. Aavatsmark, “Monotonicity conditions for control volume methods on uniform parallelogram grids in homogeneous media,” Computational Geosciences, Vol. 9, pp. 61–72, 2005. [Publisher]
26. M. J. Mlacnik and L. J. Durlofsky, “Unstructured grid optimization for improved monotonicity of discrete solutions of elliptic equations with highly anisotropic coefficients,” Journal of Computational Physics, Vol. 216, No. 1, pp. 337–361, Jul. 2006, https://doi.org/10.1016/j.jcp.2005.12.007 [Publisher]
27. T. Aadland, I. Aavatsmark, and H. K. Dahle, “Performance of a flux splitting scheme when solving the single-phase pressure equation discretized by MPFA,” Computational Geosciences, Vol. 8, No. 4, pp. 325–340, Dec. 2004, https://doi.org/10.1007/s10596-004-3774-y [Publisher]
28. G. T. Eigestad, I. Aavatsmark, and M. Espedal, “Symmetry and M-Matrix issues for the O-Method on an unstructured grid,” Computational Geosciences, Vol. 6, No. 3/4, pp. 381–404, 2002, https://doi.org/10.1023/a:1021239130404 [Publisher]
29. O. Axelsson, Iterative Solution Methods. Cambridge University Press, 1994, https://doi.org/10.1017/cbo9780511624100 [Publisher]
30. I. Aavatsmark, G. T. Eigestad, B. T. Mallison, and J. M. Nordbotten, “A compact multipoint flux approximation method with improved robustness,” Numerical Methods for Partial Differential Equations, Vol. 24, No. 5, pp. 1329–1360, Sep. 2008, https://doi.org/10.1002/num.20320 [Publisher]