On the Performance of the Node Control Volume Finite Element Method for Modeling Multi-phase Fluid Flow in Heterogeneous Porous Media

In this paper, we critique the performance of the node control volume finite element (NCVFE) method for modeling multi-phase fluid flow in heterogeneous media. The NCVFE method solves for the pressure at the vertices of elements and a control volume mesh is constructed around them. Material properties are defined on elements, while transport is simulated on the control volumes. These two meshes are not aligned producing inaccurate results and artificial fluid smearing when modeling multi-phase fluid flow in heterogeneous media. We perform numerical tests to quantify and visualize the extent of this artificial fluid smearing in domains with different material properties. The domains are composed of tetrahedron finite elements. Large artificial fluid smearing is observed in coarse meshes; however, it decreases with the increase in mesh resolution. These findings prompt the use of high-resolution meshes for the method and the need for development of novel numerical methods to address this unphysical flow.


Introduction
Finite volume and finite element discretizations have been extensively studied for the past few decades for different flow systems in heterogeneous and fractured porous media (Chen and Ewing 1997;Stefansson et al. 2018;Abushaikha 2018;Abushaikha et al. 2017;Ahmed et al. 2019;Xia and Zhang 2006;Khoei and Haghighat 2011). One important finite element method is the node control volume finite element (NCVFE) method which was developed by Baliga and Patankar (1980) to numerically solve fluid dynamics problems. They subdivided the domain using irregular triangular elements with control volumes surrounding nodes (vertices), Fig. 1. The pressure is solved on the nodes, while 1 3 the velocity components are solved on the elements sequentially. The secondary control volume mesh is necessary to assure local conservation of mass, momentum, and energy since these entities are discontinuous at the element interfaces. Quadrilateral elements were first used by Schneider and Raw (1986), and Fung et al. (1991), Forsyth (1991), while Eymard and Sonier (1994) applied the method to subsurface reservoir engineering problems'. Paluszny et al. (2007) extended the method to domains composed of different types of elements (tetrahedron, hexahedron, prims and pyramids). After that,  implemented discrete fracture models representing subsurface geological fractures using line elements embedded in two-dimensional domains composed of triangle elements, while Monteagudo and Firoozabadi (2004) embedded triangles with tetrahedron elements to solve similar problems.
The usability of the NCVFE method to simulate multi-phase fluid flow in heterogeneous media and complex geological models has increased rapidly (Mello et al. 2009;Sadrnejad et al. 2012;Schmid et al. 2013) thanks to the mesh flexibility of the finite element method and the local conservation characteristics of the control volumes. However, since the geological data (material properties such as porosity and permeability) are assigned on the elements and the dynamic data (transport variables such as saturation) are computed on the control volumes, artificial fluid smearing is observed between regions of different material properties. The effect of this issue can be decreased by refining the mesh or through the use of adaptive spatial meshing techniques as proposed by Jackson et al. (2013). However, the method will always allow the transport variables to be contained in control volumes that represent elements with different material properties, as discussed later. Vermuelen (1973) imposed the properties on the control volumes which decreased the effect of the aforementioned issue; however, the shape of the control volume depends on the finite element mesh, thus the mesh flexibility is considerably decreased. Triangle finite element mesh (dashed lines) with the corresponding node control volume mesh (solid lines) imposed on the vertices of elements. a The material properties (gray color) are defined on the elements, here representing a fracture (white is lower permeability matrix), b The pressure and transport variables are computed on the node control volumes (blue color). The control volume mesh spans both the fracture and matrix material properties promoting unphysical flow across the boundaries of elements Nick andBazar-Afkan and proposed novel numerical approaches to solve the issue of artificial fluid smearing using triangular elements, and Abushaikha et al. (2015) changed the location of the degrees of freedom from the vertices to the interfaces to minimize the degree of smearing using tetrahedron elements. It is important to note that many papers have suggested new schemes that can eliminate the effect of artificial smearing (Abushaikha et al. 2017;Abushaikha 2018;Salinas et al. 2018;Abd et al. 2019;Abd and Abushaikha 2020;Abushaikha and Terekhov 2020;Li et al. 2020;Nardean et al. 2020) even when discrete fracture modeling is implemented .
In this paper, we quantify and numerically assess the artificial smearing produced by the NCVFE for multi-phase fluid flow in heterogeneous media. We use unstructured mesh with tetrahedron elements, while the total fluid mobility is calculated by using an arithmetic average of fluid saturations of the shared control volumes. This allows for the computation of the unphysical error of the artificial smearing for various tests using domains of hugely different material properties. To visualize the artificial smearing better, we plot the results on the actual control volumes contrary to most of the literature Monteagudo and Firoozabadi 2004;Mello et al. 2009;Sadrnejad et al. 2012;Schmid et al. 2013;Jackson et al. 2013;Vermuelen 1973;Nick and Matthai 2011;Bazar-Afkan and Matthai 2011;Edwards 2006) as they plot the solutions on the vertices of elements which tends to underestimate the extent of the problem. The paper is organized as follows. In Sect. 2, we introduce the governing equations followed by a brief description of the NCVFE method in Sect. 3. We then present some numerical tests in Sect. 4 and we end with some conclusions in Sect. 5.

Governing Equations
In the following, we present the equations of a two-phase immiscible fluid flow of water and oil in heterogeneous porous media. The flow is characterized by the continuity equation and Darcys law (Darcy 1856;Bear 1972), and a slightly compressible rock is considered. The mass balance for the fluid phase is, S is the saturation of the phase, is the porosity of the rock, is the density of the phase, q is the volumetric source-sink rate of the phase, is the Darcy velocity for phase, and t is time. We also assume capillarity and gravity forces are negligible. The Darcy phase velocity is, where P is the phase pressure, K is the absolute rock permeability, and is the phase mobility, where k r is the relative permeability of phase and is fluid viscosity of phase. (1) The relative permeability is saturation dependent, while the fluid viscosity is constant. We assume a closed system with known initial boundary conditions where source/sink terms are represented by wells. Moreover, a pressure equation is written based on the continuity equations for the phases (Eq. 1). The equations of each phase are combined to form a pressure equation where time derivatives of saturations are represented explicitly (Durlofsky 1993), where C r is rock compressibility, q t is the total volumetric source-sink rate of both phases and t is the total mobility given by: Then, Darcy's law, Eq. 2, is used to calculate the phase velocities using the solution of the pressure field, Eq. 4. Moreover, we use Eq. 1 to account for the advection and transport of fluid. The pressure and the advection equations are coupled nonlinearly using the total mobility, Eq. 5. The mobility depends on the saturation that changes over time and space.

Node Control Volume Finite Element (NCVFE) Method
In the NCVFE method, the domain is subdivided into finite elements and a secondary grid is imposed around the nodes (vertices) of the elements forming a control volume mesh, Figs. 1 and 2. The pressure is solved on the vertices of elements using the Galerkin finite element method and the transport variables (i.e., fluid saturation) are solved on the control volumes. The secondary grid is constructed to provide a continuous velocity field between the control volumes, since the velocity field is discontinuous between the elements' interfaces (Fung et al. 1991;). This guarantees mass conservation of the system. Material properties, such as porosity and permeability are still defined on the Fig. 2 a Tetrahedron mesh , b The corresponding node control volume mesh (red lines) imposed on the vertices of the tetrahedron finite elements elements. Some authors refer to this method as the control volume finite element (Baliga and Patankar 1980;Fung et al. 1991;Sadrnejad et al. 2012;Jackson et al. 2013). Here, we include the term "node" as qualifier for the sake of clarity, and to distinguish the method from newly development methods in the literature (Nick and Matthai 2011;Abushaikha et al. 2015).
In the NCFVE method, a multi-phase flow problem is solved in two steps. The first step is to calculate the pressure using finite element method; here we use the well-established Galerkin method (Huyakorn and Pinder 2014). Then, the advection of fluid between the node control volumes is calculated using the finite volume method. Next, we discuss these steps in detail.

Galerkin Finite Element Method
In the finite element method, an integral formulation for the governing flow equation is derived to solve for the pressure at each node in the mesh. Let us consider the differential operator where L[P] is the differential operator defining the pressure, and F is the external source-sink terms. An approximate solution of the pressure in each element is defined as where P (e) is the approximate solution for pressure within element (e), n is the number of nodes within element (e), and P i unknown values of pressure for each node within element (e), and N (e) i is the interpolation function for each node within element (e) and the sum of interpolation functions equals to one at every point within the element In this paper, we use first-order Courant (Courant 1943) (linear) interpolation functions and the derivative of Eq. 7 equals When we substitute the pressure value in Eq. 6 with the approximate solution of pressure, Eq. 7, a residual R occurs at each node in the problem domain, and the pressure equation no longer equals zero Several methods are available to eliminate this error; here we use the method of weighted residuals Istok 2010). In this method, Eq. 10 is multiplied by a weighting function W and the weighted average of the residuals in the domain is forced to equal zero Replacing the differential operator in Eq. 11 with Eq. 4 leads to In the Galerkin finite element method, the weighing function W is identical to the interpolation function N. Based on Eq. 12, the residual for node i in tetrahedron element (e) equals where V (e) is the volume of tetrahedron. The material properties (permeability and porosity) are assumed piece-wise constant within the element Nick and Matthai 2011), and the total fluid mobility is calculated by using an arithmetic average of fluid saturations of the shared control volumes (Baliga and Patankar 1980;Huber and Helmig 1999). For simplicity, we make Eq. 13 equal to where R (e) M,i is the conductance residual, R (e) C,i is the capacitance residual, and R (e) F,i is the force residual.
We integrate by parts the first term of Eq. 14 and replace the derivative of the approximate solution of pressure using Eq. 9 Neglecting the second term in Eq. 15 since no-flow at boundaries and making the conductance matrix by rearranging the conductance residual equations of the nodes of element (e), where [M (e) ] is the conductance matrix, and for tetrahedron element equals, x dS n y , and N (e) n z are node n interpolation function derivatives in the X, Y and Z directions for element (e). They are given in "Appendix A".
For the second term in Eq. 14, we evaluate the time derivative of the approximate solution of pressure over the element volume. In this paper, we use the lumped formulation as it is considered more stable and produces fewer oscillations during multi-phase flow (Eymard and Sonier 1994;Bastian and Helmig 1999). The time derivative of the approximate solution of pressure is defined by, where N * (e) i are the lumped interpolation function for the time derivative of pressure at node i and its product equals, Substituting Eq. 18 in the second term of Eq. 14 gives the capacitance residual for node i in element (e), We make the capacitance matrix by rearranging the capacitance residual equations of the nodes of element (e), where [C e ] is the capacitance matrix of element (e) and applying Eq. 19 for Eq. 20 gives the lumped capacitance matrix for tetrahedron element, For the third term in Eq. 14, where the well is located in element (e), the flow rate is distributed on the elements nodes depending on the well location in the element, where (x 0 , y 0 , z 0 ) are the coordinates of the well location in element (e). (18) Discretizing the time derivative in Eq. 21 using backward Euler method and adding all the matrices in their corresponding global matrices, Eq. 12 becomes the final pressure equation, where t is the time index and we use implicit pressure and explicit saturation (IMPES) (Fanchi 2018). Eq. 24 is an Ax = B equation and is solved for the unknown, next time-step pressure [P] t+1 , using a linear solver. Next, we discuss the construction procedure of a node control volume mesh in a domain of tetrahedron elements.

Node Control Volume Mesh
A node control volume is a secondary grid constructed around the vertices of the finite element mesh. In this section, we discuss the procedure for constructing a node control volume mesh, and the computation of the area normal vectors for tetrahedral elements. Both the area normal vectors along with node control volume mesh are needed to assure that the fluxes on the interfaces of the control volumes are continuous and perpendicular to ensure local conservation of mass. This procedure is necessary since the velocity field is discontinuous between the elements interfaces.
A node control volume between tetrahedron elements is constructed around the shared node with hexahedron sectors. Each sector has eight vertices. The vertices are located at the shared node, the barycenter of the host element and the barycenter of the three interfaces and the three edges of the host element connected to the node; see Figs. 3 and 4. The area normal vector ( ) of each control volume is calculated where each sector has three faces connected to the flow and each face has two interface line vectors, see Fig. 4. The cross-product of the interface line vectors gives the area normal vectors of each face, The equations are performed for each sector in the node control volume. The pore volume of a node control volume in a tetrahedron element mesh is calculated by, where V (n) is the pore volume of the node control volume (n), and E is the number of elements sharing the control volumes.

Fluid Saturation Calculation in Node Control Volumes
After constructing the node control volume mesh and calculating the area normal vector for each sector in the node control volumes, we integrate the transport equation, Eq. 1, over the node control volume (n). After that, the divergence theorem is applied and the

Fig. 4 Node control volume sector around node I of a tetrahedron (top-left), vertices of sector (top-right), interface line vectors of faces for sector (bottom-left), and area normal vectors of sector (bottom-right) equation is discretized in time using the backward Euler approach (similar to Eq. 24). The next time-step phase saturation in control volume (n) is calculated by,
where SI is the number of faces in node control volume (n), and the flux (n),j) is the flux of face j in node control volume (n) and calculated by, where t w is calculated using the pervious time-step phase saturation of the upstream control volume, and ∇P (e) is calculated from Eq. 9. Equation 24 is calculated for every node control volume in the mesh, and in this paper we calculate the water saturation.

NCVFE Implementation Algorithm
Finally, we present algorithm 1 of the NCVFE implementation . The pressure is implicitly calculated and saturation calculations are explicit (IMPES). The number of nodes in the mesh represent the degrees of freedom of the system.

Numerical Tests
In this section, we perform a numerical study of the NCVFE method. We test its accuracy in modeling two-phase fluid flow in regions with large variations in their material properties, i.e. between matrix regions and sealing/conductive faults (Tests I and II).
A meshing software, GMSH, (Geuzaine and Remacle 2009) is used to create 3-D tetrahedron finite element meshes that are used in the upcoming test. Each test has a mesh of a different element size that is refined by manipulating the element length (h). A quadratic model is used to construct the oil and water relative permeability curves, and a unifrom porosity of 0.2 is used. The water and oil viscosities are 0.4 and 2.5 mPa.s respectively and the rock compressibility equals 4.0 × 10 −10 Pa −1 . At intial conditions, we assume that the domain is fully saturated with oil.
The pore volume injected (PVI) is defined as where V P is the total pore volume of the system. To visualize the tests results better, we plot the solutions on the control volumes.

Buckley-Leverett Validation
In this section, we compare the analytical solution of the Buckley-Leverett function with the saturation profiles produced by the NCVFE method at different mesh resolutions. The solution is one-dimensional with negligible gravity and capillarity effects. The tested domain is a rectangular tube dimensions of 1 m × element length × element length for the tetrahedron elements. The results are plotted in Fig. 5-a where both the analytical and the numerical profiles are shown. The analytical solutions produces a sharper front of water saturation which is expected, and a finer grid is needed (0.001 element length) to achieve a close match to the analytical solution. This shows that the NCVFE method is capable of mapping the water front movement for relatively fine grids. The upcoming tests will discuss in details the effect of the mesh resolution on the water saturation calculations. Furthermore, we computed the error L 2 norm at the different mesh resolutions for tetrahedron elements to measure the rate of convergence using the following equation: The sub-linear convergence rate is around 0.4 and can be noticed in Fig. 5b. The convergence rate is relatively low because of the sharp shock-front of the analytical solution and it is in line with other finite element descritization schemes Abushaikha et al. 2017). The overall results show that the NCVFE method is capable of modeling the nonlinear behavior of the two-phase flow and validated through Buckley-Leverett function.

Test I: Barriers Case
In our first test, we incorporate 20 barriers into a domain to test the sealing capability of the NCVFE method. The tested domain is a square of a side length of 1 m and a depth of 0.05 m where barriers aperture (width) is 0.04 m . The injector is located at the left boundary of the domain, where water is injected at a constant rate. The resolution of the mesh is varied at the injector side as well. Figure 6 and Table 1 show the mesh properties along with error in the saturation calculation (we will discuss this error in details later). Figure 7 shows the water saturation employing the three different meshes from Table 1 at different injected pore volume. In the coarse mesh, as the water propagates inside the domain, it fully invades the barriers and does not honour the zero permeability property of the elements. The method allows the water to move in a control volume that contains elements with zero permeability; the correct solution should have no flow in such elements. The degree of nonphysical water invasion decreases as the mesh is refined. This is a classic case of the artificial fluid smearing produced by the NCVFE method between regions of different of material properties. To measure and quantify this behavior, we calculate the error in water saturation in the barriers using   Table 1 1 3 where V Z is the total pore volume of the barriers, S ref w is the reference solution for water saturation in the zone which equals the initial water saturation. Figure 8 and Table 1 show this error. The error is large and decreases as the mesh is refined. We also calculate the convergence rate of the NCVFE method for this test (Janicke and Kost 1996), where N DF is the number of degrees of freedom, p is the order of convergence and d is the spatial dimension. For this case, p ≡ 1 (where d = 3 ) which is in agreement with the linear tetrahedron elements employed in the method (Janicke and Kost 1996).
Sealing faults or non-permeable layers are heavily present in subsurface reservoirs. They divide a field in to different reservoirs or compartments where the fluids are stored. Therefore modeling nonphysical communication between the reservoirs or permeability zones will produce inaccurate flow behaviors and predictions of recovery. High mesh resolutions are required to minimize this behavior, as we saw in this test.
To summarize , we used the NCVFE method to model water invasion into the nonpermeable regions of a porous domain. This invasion caused the appearance of artificial smearing problems inside the region. This nonphysical flow is attributed to the procedure of constructing the control volume around the vertices of the elements. These elements with different material properties share interfaces which promotes nonphysical communication, however this behavior becomes less significant as the mesh of the domain is refined.

Test II: Fractures Case
In this case, we model the fluid flow for a fractured network based on an outcrop sample taken from the Lias formation (limestone) from Bristol Channel, UK (Abushaikha 2013). We model the fractures explicitly using tetrahedron elements with different aperture sizes for the fractures ranging from 0.4 to 3 cm. The domain has dimensions of 68.4 × 48.6 × 1 cm and has 14 matrix blocks surrounded by conductive fractures; see Fig. 9. The injector is Fig. 8 The error of Test I resulting from estimating the water saturation in the using NCVFE method at 0.5 and 1 pore volume injected as the element length varies located at the right boundary where water is injected at a constant rate in a similar fashion to the schematic shown in Fig. 6b. The mesh properties are shown in Table 2.
To understand the implications of the flow attributed to the NCVFE method, we change the permeability of the matrix zones to negligible values. The results are considered as a reference solution that represents the water flow through the fractured zones without invading the matrices.The outcrop is shown in Fig. 9, where the permeability of the fractures is 10,000 mD. Figure 10 shows the water saturation in three different meshes from Table 2 at different injected pore volumes. For the coarse mesh, the water heavily invades the matrix regions delaying the flow in the fractured network. As the mesh is refined, the water propagates faster inside the fracture regions with fewer invasion of the matrix. To analyze this behavior, we calculate the water cut at the right hand boundary as a function of pore volumes injected to measure the breakthrough time. The water cut is given by where q t = q w + q 0 , q w and q o is the water and oil production flow rate, respectively. Figure 11a shows the water breakthrough times at the right hand boundary for the meshes of Table 2 using Eq. 35. As observed in Fig. 10, the water reaches the right hand boundary faster when the mesh is refined. This behavior is in contradiction with common mesh convergence studies since the opposite is the norm (Fanchi 2018). To distinguish which type of dispersion is causing this behavior, physical or numerical, we mesh the fracture network without the matrix regions and estimate the water breakthrough times at the right hand boundary for the element lengths of Table 2 in Fig. 11b. We notice that as the fracture network is refined, the water breakthrough time is delayed; opposite to Fig. 11a. Therefore, the artificial smearing produced by NCVFE is causing this behavior where the physical dispersion is eclipsing the truncated numerical error of the mesh refinements promoting the delay of breakthrough times for the coarse meshes.
To measure this artificial smearing in the matrix regions, we use Eq. 33 to define the error. Figure 12 and Table 2 show the error produced, and as the mesh is refined the (35) WC = q w q t Fig. 9 a Outcrop sample from Lias formation (limestone) from Bristol Channel 2, b discrete fracture models of the outcrop case error decreases. Similarly, using Eq. 34, a first-order convergence rate is noticed where d = 3 . It is worth mentioning that the first half of the convergence rate slope is also first order because the element length of the first two meshes is larger than the depth of the domain; effectively we have d = 2.
In this test, a delay of flow in the fracture region and breakthrough time was noticed due to the numerical extension of the fracture width. This can be attributed to artificial smearing that causes inaccurate predictions of water saturation calculations. The role of fractures is very important in understanding the flow in reservoirs and optimizing recovery factors, and their properties (dimensions and aperture size) greatly influence the kinetics and the recovery of hydrocarbons. Therefore, high-resolution meshes are  recommended for NCVFE when modeling highly conductive zones especially between regions with large variations in their material properties, as we saw in this test and the previous test. Simulations such as these could be used to determine empirical fracturematrix transfer functions (Abushaikha 2008) that encapsulate the average recovery in geological realistic fracture networks to be used in full field case studies.

Summary and Conclusions
In this paper, we have investigated the accuracy of the NCVFE method to model two-phase fluid flow in heterogeneous and fractured porous media. In modeling sealing faults, the method performed poorly as the barriers were fully penetrated by the displacing water. The water did not respect the zero-permeabilities defined on the elements and flowed through them, because the NCVFE encompasses elements with both zero and finite permeabilities. In representing conductive faults (fractures), the NCVFE method did not perform any better. Very high mesh Fig. 11 a The water cut at the right-hand boundary for Test II versus pore volumes injected for the meshes without the matrix regions. b The water cut at the right hand boundary for Test II versus pore volumes injected Fig. 12 The error of Test II resulting from estimating the water saturation in the using NCVFE method at 0.1 and 0.45 pore volume injected as the element length varies resolutions were needed to obtain a reasonable representation of the fluid flow inside the fracture region. The low resolution mesh suffered a large amount of unphysical flow to the matrix regions and a major delay in breakthrough time compared to the high-resolution meshes. We conclude that NCVFE requires high mesh resolutions to model multi-phase flow in heterogeneous systems reasonably. The critical drawback of the NCVFE method is that the material properties are defined on the elements, while the fluid saturation is calculated on the vertices of elements. Essentially, the saturation is updated on a coarse mesh that is not aligned with the assignment of rock properties. = (x i , y i , z i ) = (x j , y j , z j ) = (x k , y k , z k ) = (x l , y l , z l )