Consistent MPFA Discretization for Flow in the Presence of Gravity

A standard practice used in the industry to discretizing the gravity term in the two‐phase Darcy flow equations is to apply an upwind strategy. In this paper, we show that this can give a persistent unphysical flux field and an incorrect pressure distribution. As a solution to this problem, we present a new consistent discretization of flow, termed Gravitationally Consistent Multipoint Flux Approximation (GCMPFA), which is valid for both single‐ and two‐phase flows. The discretization is based on the idea that the gravitational term in the flow equations is treated as part of the discrete flux operator and not as a right‐hand side. Here, the traditional formulation representing pressure as a potential is extended to the case including gravity by introducing an additional set of right‐hand side to the local linear system solved in the MPFA construction, thus obtaining an expression of the fluxes in terms of jumps in cell‐center gravities. Numerical examples showing the convergence of the method are provided for both single‐ and two‐phase flows. For two‐phase flow, we show how our new method is capable of eliminating the unphysical fluxes arising when using a standard upwind scheme, thus converging to the correct pressure distribution.

In this work, we discuss MPFA methods. MPFA is a control volume method introduced independently by two different research groups in 1994 (Aavatsmark et al., 1994;Edwards & Rogers, 1994). The two approaches differ on the choice of geometrical points and control volume grids. Here, we only consider the so-called O-method developed by Aavatsmark and coworkers. They first introduced MPFA for general quadrilateral grids in Aavatsmark et al. (1994) and Aavatsmark et al. (1996) and then extended the method to triangular and polygonal grids in Aavatsmark et al. (1998aAavatsmark et al. ( , 1998b. The reader can refer to Aavatsmark (2002) for an excellent review on MPFA methods for quadrilateral grids and to Aavatsmark et al. (2007) for a numerical investigation on its convergence properties. The convergence properties of MPFA have also been investigated for general quadrilateral grids in Eigestad and Klausen (2005), Klausen and Winther (2006a), and Klausen and Winther (2006b) and for unstructured triangular grids in Bause et al. (2010). In particular, using a specific numerical quadrature, the MPFA and MFMFE methods were shown to be equivalent (Klausen & Winther, 2006a). A somewhat simpler MPFA variant is the so-called L-method by Aavatsmark et al. (2008). Finally, when the MPFA method is applied to multiphase flow, a monotone scheme is desirable. Local criteria which ensure monotonicity for general control volume methods on heterogeneous media are given in Nordbotten et al. (2007).
MFMFE methods were introduced for incompressible Darcy flow problems on triangular and convex quadrilaterals in Wheeler and Yotov (2006) using the lowest order Brezzi-Douglas-Marini spaces (Brezzi et al., 1985). Extensions to slightly compressible flow and multiphase flow are presented in Arrarás and Portero (2019) and Wheeler et al. (2012), respectively.
The common feature of all these methods is the treatment of the gravity term in the Darcy flow equations. The traditional approach has been to represent the pressure as a potential and let the discretization consider only deviations from the potential, by ignoring gravity effects in the discretization. Gravity is only considered in the model equations in Arrarás and Portero (2019) and Wheeler et al. (2012), yet it is disregarded in the numerical examples reported therein.
However, this approach is inconsistent when the gravity term in the Darcy flow equation is inhomogeneous, as caused for example by two-phase effects, density variations, stepwise variations of permeability, and certain variants of vertically averaged models for CO 2 storage (Nordbotten & Celia, 2011). While for single-phase flow MPFA can handle discontinuities in the fluid potential caused by, for example, smooth variations of permeability from cell to cell, in two-phase flow discontinuities also arise due to presence and absence of a mobile phase, and this kind of discontinuities can create unphysical fluxes, for example, in the case of a fluid-fluid interface at conditions of vertical equilibrium (see Aavatsmark et al. (1994) for a discussion on this). Here, we show that a standard treatment of the gravity term based on an upwind strategy in the multiphase Darcy flow equations leads to the creation of a persistent unphysical flux field, even in absence of any external forces, and gives an incorrect pressure distribution.
An aim of this paper is therefore to develop a consistent discretization of the porous media flow equations in the presence of gravity which amends this crucial shortcoming. To do so, we treat the gravity term as part of the discrete flux operator and derive an expression of the fluxes in terms of jumps in cell-center gravities. The paper is organized as follows. First, we introduce the governing equations of single-and two-phase flow in porous media in section 2. Then, a standard discretization of the gravity term, and a discussion on its limitations, is presented in section 3. Our new consistent discretization of the flow equations in the presence of gravity is presented for single-and two-phase flow in sections 4 and 5, respectively. Numerical examples are provided in both sections. Finally, concluding remarks are given in section 6.

Single-Phase Flow
Incompressible single-phase Darcy flow in nondeformable porous media is governed by the following equation: where p is pressure, K is the (generally heterogeneous) absolute permeability tensor divided by fluid viscosity, g represents gravitational forces (density times acceleration due to gravity vector), which is a function of space, and is a source term. We emphasize that equation (1) represents incorporation of Darcy's law into a mass conservation equation.

Two-Phase Flow
The Darcy formulation given by equation (1) is extended to incompressible immiscible two-phase flow as follows: where is porosity, t is time, s is the phase saturation associated with phase = 1, 2, and is the phase mobility, which is an increasing function of s . In equation (2), K represents the absolute permeability tensor, as fluid viscosity is incorporated into . Introducing the total quantities Σ = ∑ , and in absence of capillary pressure, that is, p 1 = p 2 = p, summing equation (2) yields where G Σ = ∑ ( g ). The phase fluxes q can then be expressed in terms of the total flux q through the fractional flow function = ∕ in the following manner: Choosing one saturation as primary variable, say, for example, s 2 = s, equation (2) is reformulated in terms of total flux as Equations (3)-(5) form a system of two equations for two unknowns (s and p).

Standard Discretization of Flow
Solution of equation (1) using control volume methods involves the computation of the flux f k through some surface k of the control volume, defined as where n is the unit normal vector to the surface. Likewise, for two-phase flow, from equation (3) one has the total flux

Traditional Potential Formulation for Single-Phase Flow
Representing the pressure as a potential h, and ignoring gravity effects in the discretization, calculation of the flux in equation (6) reduces to the solution of the integral

One-Dimensional Problems
For one-dimensional problems, the flux over the surface between two neighbor cells 1 and 2, f 12 , is approximated by a two-point stencil (see Figure 1 left) as follows: where T 12 is the transmissibility of the surface which is calculated as harmonic average of the cell transmissibilities, that is, The cell transmissibilities are defined as T i = K i ∕Δx i , i = 1, 2, where Δx i is the length of cell i.

Multidimensional Problems
For multidimensional problems, the flux is approximated using the MPFA method as where the coefficients t k,i are called transmissibility coefficients. Calculation of the transmissibility coefficients works as follows. A dual grid is created by connecting the cell centers with the face centers. In this manner, each cell is partitioned into subcells (three and four in triangular and quadrilateral grids and six and eight in tetrahedral and hexahedral grids, respectively) and each face is subdivided into subfaces (two in 2-D grids, four in 3-D grids). Subcells are then grouped together to form an interaction region surrounding each node (see Figure 1 right). Then, the following principles are applied: 1. Potential is assumed to be linear in each subcell. 2. Flux continuity is enforced at the subfaces. 3. Potential continuity is enforced at single points on the subfaces, called continuity points.
There is a whole class of MPFA methods for such grids, depending on the choice of the location of the continuity points. Here we only consider the O-method described by Aavatsmark (2002). Principles (1) and (3) imply that for a subface k with adjacent subcells j 1 and j 2 , one has where ∇h is the subcell gradients and d is the distance between the continuity point and the cell centers. For flux continuity, principle (2) is written where n k is the normal vector of the face and K i is the permeability with respect to cells j 1 and j 2 . Collecting all equations (12)-(13) for each interaction region, a local linear system is recovered of the following form: The first row represents flux balance (13) and involves only the subcell potential gradients. The matrix G contains the discretized Darcy's law, that is, the n·K products on a subcell level. The second row gives pointwise potential continuity (12) over cell faces. Matrix D contains the distances d, while matrix I ± contains ±1 coefficients depending on which side the cell is relative to the face normal vector. The third row together with the right-hand side enforces a unit potential in one cell after another. Equation (14) can be inverted to compute the subcell potential gradients dh as functions of the cell-center potentials, effectively computing basis functions for the discretization. Hence, from solving equation (14) we obtain the transmissibility coefficients of the potential-to flux maps, denoted as k,i , which represent the contribution of cell i to the flux across the sub-face k. To obtain the full face coefficients t k,i , we sum over the subfaces l of face k The discretized flux across the face is then given by equation (11).

Standard Discretization of Single-Phase Flow 3.2.1. One-Dimensional Problems
In one-dimensional problems, we have seen that, in absence of gravity, the face transmissibility is calculated as harmonic average of the transmissibilities of the adjacent cells. If gravity is present, assuming that the pressure p is linear and gravitational forces g are constant within each cell, the flux continuity is given by wherep is the pressure at the interface between the two cells. Introducing the cell transmissibilities, equation (16) can be solved forp to getp Inserting this expression back into, say, the left-hand side of equation (16) gives the flux expression in presence of gravity as Equation (18) shows that for the pressure term the harmonic average of the cell transmissibilities is retrieved, whence the flux due to gravity is given by the product between the harmonic average of the cell transmissibilities times the arithmetic average of the cell gravitational forces.

Multidimensional Problems
Calculation of the full flux in the presence of gravity in equation (6) involves the computation of the term When this term is treated as a right-hand side in equation (1), a standard discretization approach is to extend the result of equation (18) to the multidimensional case and use the harmonic average of the cell permeability tensors. Defining d j = |x k − x j |, the distance between the center of face k, x k , and the center of cell j, x j , where j is either of the two cells j 1 and j 2 with mutual face k, the flux due to gravity is then computed as where the operator ⟨K⟩ k denotes the d-weighted harmonic average of the permeability tensors between the two cells j 1 and j 2 and g k is the weighted arithmetic average of the cell gravity vectors The full flux in the presence of gravity is then given by where the transmissibilities coefficients t k,i are calculated as described in section 3.1.2.

Standard Discretization of Two-Phase Flow 3.3.1. Numerical Method
Solution of equations 3-(5) is done using the Implicit Pressure Explicit Saturation scheme (Chen et al., 2006). Starting from a known saturation s n , the Implicit Pressure Explicit Saturation scheme works as follows: 1. The pressure p n is calculated implicitly by solving (3) and the total flux q n Σ is reconstructed from p n . 2. The saturation is advanced in time explicitly from (5) as Calculation of the face mobilities in equation (24) is done using the method outlined in Moortgat et al. (2011). The method works as follows: 1. First, we pick the phase for which the phase flux has the same sign of the total flux. This is the heaviest phase when (q · n)(Kg · n) > 0 or the lightest phase when (q · n)(Kg · n) < 0. This sign determines the first upwind phase mobility 1 . 2. For the other phase, we have two options. As a first guess, we assume that the second phase has the same sign as the first phase and the total flux and take the upwind mobility 2 accordingly. 3. We can now evaluate the phase fluxes using equation (4) and check consistency. If the guessed sign is retrieved, then the process is complete, otherwise the opposite upwind choice for 2 is made in step (2). When the total flux is zero, the upwind directions can be determined explicitly. A standard discretization of the total flux in equation (7) is then obtained by applying the traditional MPFA construction to the pressure term and an upwind scheme to the gravity term as follows (Enchéry et al., 2002): However, a discretization of such a kind on rough grids is prone to creating unphysical fluxes, as illustrated in the following section. In particular, no flow boundary conditions are assigned to all boundaries, that is, q = q 1 = q 2 = 0. Blue is the heaviest phase (s 2 = 1); green is the lightest phase (s 2 = 0).

Numerical Example
Let us consider the case of a system composed by two incompressible fluids trying to reach vertical equilibrium conditions in absence of any external forces. The two fluids have same viscosity equal to 1.0e −3 Pa·s, the gravity vector g is directed downward along the vertical direction, porosity is equal to 0.2, and the permeability tensor K is heterogenous with four layers of different permeabilities, that is, K = a i KI, with K = 1Da and values of a i reported in Figure 2. No-flow boundary conditions are assigned to all boundaries, that is, q = q 1 = q 2 = 0. Initially, a horizontal interface is considered, with the upper region fully saturated with the heaviest phase ( 2 = 1000 kg/m 3 ) and the lower region fully saturated with the lightest phase ( 1 = 100 kg/m 3 ). Computations are carried out on quadrilateral randomly perturbed grids with five levels of refinement, that is, N = 4, 8, 16, 32, 64 number of cells per side. In virtue of equation (4), counter current flow between the two phases should thus establish, leading eventually to conditions of vertical equilibrium. However, numerical simulations using equation (25) indicate that a persistent spurious flux fields originate  (see Figure 3). After some oscillations, the system yet reaches a stable configuration; however, the computed pressure field is far from vertical equilibrium. This is clearly shown in Figure 4, displaying the time evolution of the errors in pressure and saturation and the maximum total flux. The errors in pressure and saturation are computed using the following L 2 metrics: where is the computed variable, Δ is the element volume, and the exact variables are the ones measured at vertical equilibrium conditions. As Figure 4 clearly shows, the saturation converges to the equilibrium conditions; however, the nonvanishing total flux prevents the pressure to converge to vertical equilibrium conditions. This is because of the inconsistent discretization of the gravity term in the pressure equation, while a standard upwind scheme is sufficient for the transport equation. Nevertheless, the standard upwind scheme converges to the exact pressure field with refinement of the grid (first-order convergence; see Figure 5).
In the following sections, we present a consistent discretization of the flow equations in the presence of gravity, which is capable of eliminating this unphysical flux field and thus gives the correct pressure field.

Numerical Method
For multidimensional problems, a consistent treatment of gravitational forces can be achieved by a more nuanced approach to the local flux balancing within the local construction of the discretization scheme. In the presence of a gravitational field, equation (13) is extended to read where g represents gravitational forces in the cells. We make the observation that jumps in the gravitational forces over the subfaces, [[n k · Kg]] k , will act as a flux imbalance and thus induce an additional pressure gradient in the subcells. To extend the MPFA formulation to equation (27), we introduce an additional set of right-hand side functions, which applies nonzero conditions to the first row of (14). These additional right-hand side functions thus solve We now slightly reformulate equation (11) by considering f k,j , which is the flux in absence of gravity across a subface k as evaluated in cell j (where j is either of the two cells j 1 and j 2 with mutual face k). The extended version of (11) is written for completeness as We make the note that the integral which is approximated is now stated slightly more precisely, in the sense that the integration volume Δ j from which the boundary integral appears is explicit. Also, for the two cells j 1 and j 2 where f k,j is defined, it is clear from equation (14) that k, 1 = k, 2 and therefore k, 1 ,i = k, 2 ,i = k,i . The transmissibility coefficients of the pressure-to-flux maps in absence of gravity, k,i , for the subface k, are then obtained by solving (14).
Similarly, let us denote the coefficients from (28) as k,j,l , which represent the flux across l due to a flux imbalance at k, as evaluated in cell j. We quickly note from equation (28) that for l ≠ k, then as above k, 1 ,i = k, 2 ,i . However, this will not be the case for k = l, due to the flux imbalance, indeed in this case We then obtain the full flux in the presence of gravity as It is noted that due to equation (30), it follows that as expected k, 1 = k, 2 . Thus, the second subscript can be omitted as soon as a convention is chosen for what side the flux evaluation should be considered on. Therefore, one can make equation (31) symmetric by taking the mean of the two sides, that is, Finally, we note that we can represent the K-weighted jump operator over l in terms of vector coefficientš μ l, as and the mean of the cell gravities in terms of the coefficientsμ k, as With this in mind, we obtain the compound coefficients therefore, equation (32) can be written only in terms of cell-center sums as We term this approach Gravitationally Consistent Multipoint Flux Approximation (GCMPFA).

Numerical Examples 4.2.1. Problem Formulation
In these examples, incompressible flow in a unit square domain is considered. The domain has a discontinuity line of equation rx + sy = , where 0 ≤ r, s, ≤ 1 and r + s = 1. Gravitational forces are given as a linear combination of two contributions, namely, a step function across the discontinuity line H(x, y) which is normal to the discontinuity line, and a smooth function P(x, y), as follows: In the latter equation, a 1 and a 2 are two constants and, given the unit vectors e x and e y , H(x, y) and P(x, y) have the following form: with To test the convergence properties of the method, we choose an analytical solution such that so that zero normal flux conditions hold everywhere. The method is tested on different grids, namely, quadrilateral and triangular h-perturbed grids with horizontal discontinuity line (r = 0; see Figures 6a and 6b, respectively) and one regularly perturbed grid with an arbitrary discontinuity line (r = 0.7; see Figure 6c To test the implementation, four cases are considered, depending on the values assigned to the coefficients a i in equation (38) (see Table 1). In the first test, piecewise constant gravitational forces are considered (a 2 = 0). For this test, the GCMPFA method is expected to be exact. In the second test, there is no jump discontinuity, and gravitational forces are represented as a smooth field (a 1 = 0). Tests 3 and 4 have gravitational forces given by linear combination of H(x, y) and P(x, y) with different weighting coefficients. Finally, we make a comparison between our GCMPFA method given by equation (37) and the standard method given by equation (23).

Convergence Results
In the reported results, we consider errors using the following L 2 metrics For the convergence study, all simulations are run on a personal Desktop using Porepy (Keilegavlen et al., 2019), an open-source software framework for flow and transport in deformable fractured porous media  Tables 2 and 3. Table 2 shows the results for the error p and its asymptotic convergence rate O p , while Table 3 shows the results for the error q and its asymptotic convergence rate O q . As expected, when gravitational forces are piecewise constant (Test 1), the GCMPFA method is exact to working precision for both pressure and fluxes, while a standard treatment of the gravity term leads to a discretization error. We remark that the two approaches coincide if the grid is K-orthogonal. It is also noted that for this test, the convergence rate for pressure of the standard method for the h-perturbed grids is generally worse than the second order usually obtained using traditional MPFA methods without gravity. The standard method recovers second-order convergence for Test 2. In this case, the two methods behave similarly. The superiority of the GCMPFA method over the standard method is clearly highlighted when the gravity field is a smooth discontinuous function (Tests 3 and 4). In this case, the GCMPFA method always retains second-order convergence for both pressure and fluxes, independently of the magnitude of the weighting coefficients a i . Conversely, the standard method always shows a reduction in convergence rate for fluxes O q to 1.5 and only achieves second-order convergence for pressure when the magnitude of the smooth part is much greater than that of the discontinuous part, that is, when a 2 ≫ a 1 (Test 4). We summarize the results of Tables 2 and 3 heuristically as follows.
• For h-perturbed grids, the GCMPFA method exhibits a numerical convergence following • For h-perturbed grids, the standard method exhibits a numerical convergence following where is equal to 1 if a 1 = 0 and is equal to 0 otherwise.
We remark that the results presented in equations (44) and (45) are based solely on the tests considered here, as the framework for proving convergence of MPFA methods without assuming smoothness of the permeability coefficient typically does not yield convergence rates, since a priori knowledge of the regularity of the solution cannot be assumed (Agelas & Masson, 2008).

Numerical Method
A consistent discretization of the gravity term for two-phase flow is done by extending the flux formulation given by equation (31) to the total flux in the following manner: In terms of cell-center sums only, equation (46) is written which is the two-phase counterpart of equation (37). The remaining part of the algorithm works as illustrated in section 3.3.1.

Numerical Example
We consider the same example of section 3.3.2, and we test whether the new GCMPFA discretization given by equation (47) is capable of eliminating the spurious flux field arising when using the standard upwind method given by equation (25). Figure 7 shows the time evolution of the cell saturations of the heaviest phase obtained using the GCMPFA method. Comparing Figure 7 to the same figure obtained using the stadard upwind method (see Figure 3), two things can be noted. First, the spurious flux field vanishes once the two fluids approach their equlibrium configuration, that is, for t > 1 × 10 5 s (subfigure (e) onward). Second, countercurrent flow is more uniformly distributed, that is, no oscillating saturations are observed at the near-interface region (compare subfigures (d)-(g)). Figure 8 shows the errors in saturation and pressure and the maximum total flux as a function of time for the two methods for N = 64. As the figure clearly shows, as opposite to the standard method, the GCMPFA method is capable of eliminating the spurious flux field (see Figure 8c), and thus, the pressure converges to conditions of vertical equilibrium (see Figure 8b). The saturation is not substantially affected by the different solution methods for the pressure equation; however, it shows faster convergence with the GCMPFA method (see Figure 8a).

Conclusions
We presented a novel consistent discretization of flow for inhomogeneous gravitational fields, valid for both single-and two-phase flows. The discretization is based on the idea that the gravity term is treated as part of the discrete flux operator and not as a right-hand side. This is achieved by introducing an additional set of right-hand side to the local linear system solved in the MPFA construction, thus obtaining an expression of the fluxes in terms of jumps in cell-center gravities. We provided numerical examples showing the convergence of the method. For single-phase flow, the examples indicate that for rough grids we have a general second-order convergence of the scheme in terms of both pressure and fluxes. This is in contrast to the standard discretization approach for the gravity term using the harmonic average of the cell permeability tensors. For this latter discretization, second-order convergence is reduced when the gravity undergoes stepwise variations from cell to cell. Finally, we provided numerical evidence that, in contrast to the standard upwind strategy used in the industry, the GCMPFA is capable of equilibrating a system of two incompressible fluids in absence of any external forces. This is particularly useful in reservoir simulations applications when vertical equilibrium conditions are sought as initial conditions. We conclude by remarking that, although the numerical examples presented here are two-dimensional, our method is general to multidimensional problems. Besides, extension to three-phase models is straightforward. Extensions to other variants of MPFA methods (such as the L-method) or to slightly compressible flow are also possible, but they are not addressed here.