{'fraction_non_alphanumeric': 0.042360449854595056, 'fraction_numerical': 0.01492652546516715, 'mean_word_length': 4.167974464689453, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Direct Mixed-Enriched Galerkin Method on Quadrilaterals for Two-phase Darcy Flow

We develop a locally conservative, finite element method for simulation of two-phase flow on quadrilateral meshes that minimize the number of degrees of freedom (DoFs) subject to accuracy requirements and the DoF continuity constraints. We use a mixed finite element method (MFEM) for the flow problem and an enriched Galerkin method (EG) for the transport, stabilized with an entropy viscosity. Standard elements for MFEM lose accuracy on quadrilaterals, so we use the newly developed AC elements which have our desired properties. Standard tensor product spaces used in EG have many excess DoFs, so we would like to use the minimal DoF serendipity elements. However, the standard elements lose accuracy on quadrilaterals, so we use the newly developed direct serendipity elements. We use the Hoteit-Firoozabadi formulation, which requires a capillary flux. We compute this in a novel way that does not break down when one of the saturations degenerate to its residual value. Extension to three dimensions is described. Numerical tests show that accurate results are obtained.


Introduction
For many years now, finite element methods for Darcy flow and transport have been developed on rectangular and simplicial meshes. However, many methods lose accuracy when posed on distorted rectangular meshes. In this work, we develop an accurate and efficient numerical method for twophase flow in porous media on meshes of quadrilateral elements. We also discuss extensions to meshes of cuboidal hexahedra for problems posed in three dimensions.
We are interested in non-rectangular meshes for at least three reasons: Firstly, space discretization on quadrilaterals uses half the number of elements compared to discretization on triangles with the same scale h, and on hexahedra, only 1/5 or 1/6 compared to tetrahedra. Secondly, in most geoscience applications, discretization is determined by natural geologic layering. Lastly, problems with slightly deformable porous media require distorted meshes.
A classic procedure in the finite element method is to define appropriate finite elements on squares or cubes, which is relatively easy to do, and then to map these elements to quadrilaterals or hexahedra. This works well under affine mappings. However, a general quadrilateral or hexahedron is not affine equivalent to a reference square or cube. Rather, the simplest map is bi-or tri-linear, which leads to inaccuracies in the approximation properties of the finite element spaces. The fundamental idea we use to circumvent this problem, used also by other authors, is to define the finite element space directly on the physical element instead of mapping it from a reference element. The direct finite element will contain polynomials, and so accurate approximation will result. Two-phase flow in porous media is governed by a system of two partial differential equations (PDEs) combined with two equality constraints. The PDEs can be formulated arXiv:1811.01099v1 [math.NA] 2 Nov 2018 as a nearly elliptic or parabolic flow problem and a nearly hyperbolic transport problem.
Many locally conservative methods have been developed for elliptic flow problems. We mention just a few of these, the mixed finite element method (MFE) [23,17], the enhanced velocity method [46,43], the multipoint flux approximation methods (MPFA) and the multipoint flux mixed finite element methods (MFMFE) [1,48,49,30], the mimetic finite difference methods [29], and discontinuous Galerkin (DG) methods [47,13,38,39,28]. We will approximate the flow problem using mixed finite element methods, which are locally conservative and give very accurate velocities. However, as is well-known, classic mixed finite elements defined on a square or cube and mapped to a general convex quadrilateral or cuboidal hexahedron perform poorly; in fact, they fail to approximate the divergence in an optimal way (except the ABF and Devloo et al. spaces [11,16,40], which have many excess DoFs). Recently, Arbogast and Correa [3] resolved the problem on quadrilaterals. They defined the AC spaces, two families of mixed finite elements that achieve optimal convergence properties on quadrilaterals, while maintaining efficiency by using the minimal number of Degrees of freedom (DoFs) possible. The current authors in [8] and Cockburn and Fu in [21] defined H(div)-conforming mixed finite elements on cuboidal hexahedra that maintain accuracy and use the minimal number of DoFs on general hexahedra. (The current authors also defined some new mixed elements on quadrilaterals similar to the AC spaces in [7].) Many methods have also been devised for approximation of the transport problem, including some of the methods mentioned above, ELLAM and characteristic methods [22,23,9,44,4,45], the standard continuous Galerkin (CG) finite element method, and the enriched Galerkin (EG) method [41,34]. The DG method utilizes discontinuous piecewise polynomial finite element spaces to approximate the solution and weakly enforces interelement continuity by penalty terms. However, the method uses a very high number of DoFs. The CG method uses many fewer DoFs, but it fails to provide locally conservative saturations, which can lead to non-physical results. We will use an EG method, which resolves this deficiency of the CG method. Sun and Liu [41] defined the EG method by enriching the approximation space of the CG method with elementwise constant functions and using it in the DG formulation. The EG method significantly reduces the global number of DoFs compared to DG, and, in fact, there is an efficient solver for elliptic and parabolic problems with EG approximations [34]. The EG method has been applied to miscible displacement problems and twophase flow in porous media in recent works [36,35].
We will base our finite element space for EG on serendipity finite elements [20,10], since these use a fewer number of DoFs than full tensor product Lagrange finite elements. However, it is well known that the accuracy of serendipity finite elements degrades when mapped from a reference element to a quadrilateral or hexahedron (at least for elements higher order than bilinear), while Lagrange finite elements maintain accuracy [10,31,33]. Recently, the current authors introduced new, direct serendipity finite elements [7] that have the same number of degrees of freedom as the classic serendipity elements but maintain accuracy on general non-degenerate convex quadrilaterals. We will use these elements.
There are many ways to separate the PDE system into flow and transport problems. We use the approach promoted by Hoteit and Firoozabadi [27,26], solved using an IMPES operator splitting algorithm [19]. The formulation requires construction of the divergence of the capillary flux, and we will provide a novel implementation that does not break down when the system degenerates (i.e., one of the saturations tends to its residual value). To avoid spurious oscillations due to sharp gradients in the solution of the transport problem, the stabilization technique of Guermond, Popov, and collaborators [25,18,50] will be used.
To simplify the presentation, we concentrate on applying our new method to problems on quadrilateral meshes. In the next section we define the finite elements that we use on quadrilaterals. The equations governing two-phase flow are given in Section 3. We give the finite element method for flow in Section 4 and for transport in Section 5. The IMPES coupling and extension to three dimensions is described briefly in Sections 6-7, respectively. Numerical tests are given in Section 8, and we summarize and conclude our results in the final section.

Direct finite elements on quadrilaterals
In this section, we review the finite element spaces that we will use. Let Ω ⊂ R 2 be a polygonal domain with boundary ∂ Ω , and let ν ν ν denote the outward unit normal vector on the boundary. We impose a conforming finite element mesh T h of quadrilaterals over the domain Ω of maximal diameter h. We will assume that the mesh T h is shape-regular [24,3], which ensures that the mesh elements are not highly elongated nor degenerate nearly to triangles. The outer unit normal to element E ∈ T h is ν ν ν ∂ E . The set of all interior edges or faces of T h is its skeleton, and it is denoted by E h . Fix a unit normal vector ν ν ν γ for each edge γ ∈ E h . LetP k = span{x i y k−i : i = 0, 1, . . . , k} be the vector space of homogeneous polynomials of exact degree k ≥ 0, and let P k = k n=0Pn = span{x i y j : i, j = 0, 1, . . . , k; i + j ≤ k} denote the space of polynomials of degree up to k ≥ 0. Let Q k = span{x i y j : i, j = 0, 1, . . . , k} be the space of tensor product polynomials of degree up to k ≥ 0. Finally, P 2 k is the space of 2-dimensional polynomial vectors for which each component is in P k . We may specify that the domain of definition of the polynomials by writing, e.g., P k (E) for domain E.
At times we will need a reference element, so fix it to be the squareÊ = [−1, 1] 2 . For E ∈ T h , the standard bilinear map F E :Ê → E mapping vertices to vertices is bijective. This gives rise to the map F E , which maps a functionf : We also have the Piola transform P E based on the bilinear map F E . It maps a vectorv :Ê → R 2 to a vector v : E → R 2 , and it preserves the normal componentsv ·ν ν ν [17].

The AC mixed finite element spaces
Arbogast and Correa developed two families of mixed finite elements on quadrilateral meshes [3] for approximating (u, p) solving a second order elliptic equation in mixed form where f ∈ L 2 (Ω ) and the tensor a is uniformly positive definite and bounded. For E ∈ T h and index s ≥ 0, the AC elements are defined in terms of polynomials and a supplemental space of functions S AC s (E). In terms of the vectors xP s = {(x, y) p(x, y) : p ∈P s }, the family of full H(div)approximation elements are while, the family of reduced H(div)-approximation elements are, for index s ≥ 1, The vector elements merge together to form H(div) conforming spaces over Ω (i.e., the normal components of the vectors are continuous across each γ ∈ E h ), while the scalar elements remain discontinuous on Ω . The former family approximates u, ∇ · u, and p to order O(h s+1 ), while the latter family approximates u to order O(h s+1 ) and ∇ · u and p or order O(h s ), as long as the mesh is shape regular. These are the optimal convergence rates. Moreover, the elements have the minimal number of DoFs possible for these rates under the restriction of H(div) conformity. Using the Piola transformation P E , the space of supplemental vectors on E is where σ s i = P Eσ s i , i = 0, 1, 2, and are defined on the reference elementÊ. Different from the implementation of classical mixed finite elements, the Piola mapping is only used to define the supplemental vectors. It is easy to implement these elements using the hybrid mixed finite element formulation [12,17].  The current authors defined families of direct serendipity finite elements in a recent paper [7]. On E ∈ T h , the element of index r ≥ 1 takes the form and each family of elements is determined by the choice of the two supplemental functions spanning S DS r (E). The space can be merged into continuous (i.e., H 1 conforming) finite elements, and it has a minimal number of degrees of freedom (DoFs), which are depicted in Fig. 1. Its local dimension is dim DS r (E) = dim P r + 2 = 1 2 (r + 2)(r + 1) + 2.
A very general and explicit construction of these supplements is given in [7]. They can be defined directly on E, or they can be defined onÊ and mapped to E. In this work, we recall and use one of the simplest families of direct serendipity elements. When r = 1, it is known that the standard space of bilinear polynomials Q 1 defined onÊ and mapped to an element E ∈ T h maintains accuracy [10]. These have the form (6), with Henceforth, we only describe the new direct serendipity finite elements for indices r ≥ 2.
x v,24 = e 2 ∩ e 4 (see Fig. 2). We define the linear polynomial i (x) giving the distance of x ∈ R 2 to edge e i in the normal direction as where x i ∈ e i is any point on the edge. The functions are positive over the interior of E.
When r ≥ 2, we can define the shape functions associated to each vertex as Interior shape functions appear when r ≥ 4. In this case, we take the shape functions so that they span the entire space 1 2 3 4 P r−4 . Since these are internal to E, the precise choice of basis is not so important. The most interesting shape functions are those associated to the edges. Let and define rational functions (note that the denominators do not vanish on Therefore, There are 2(r − 1) shape functions associated to the edges e 1 and e 2 , and they are In a similar way, shape functions associated with edges e 3 and e 4 are The edge shape functions are regular polynomials of degree r except the last two functions in each direction, which are rational functions. However, all shape functions restrict to polynomials of degree r on the edges. Finally, and The unisolvence of the direct serendipity finite element space is proved generally in [7], or specifically for this choice of R V , R H , V and H in [6]. A nodal basis can be constructed using local linear algebra as described in [6,7].

Enriched direct serendipity spaces
Let the discontinuous finite element spaces of order r over T h bē We remark thatX Q 0 =X P 0 are the spaces of piecewise constants on each element. The continuous finite element spaces of order r over T h arẽ where C 0 (Ω ) is the set of continuous functions over Ω . Finally, we define the enriched finite element spaces of order r over T h as When used in a standard discontinuous Galerkin weak formulation, we can obtain locally mass conservative, discontinous or enriched Galerkin approximations. We denote these methods by DG-Q r , DG-P r , EG-Q r and EG-DS r , using, respectively, the elementsX Q In Table 1, we show the number of degrees of freedom of these four methods for an n × n two dimensional quadrilateral mesh. All the methods achieve the same convergence rate. The EG-DS r methods ultilize the fewest number of degrees of freedom. Table 1: A comparison of the global number of DoFs of DG-Q r , DG-P r , EG-Q r , and EG-DS r on an n × n mesh.
3 Two-phase flow formulation Let subscript w denote the wetting phase and subscript n the non-wetting phase. The pressures, Darcy velocities, and saturations of each phase are denoted p α , u α , and S α , respectively, α = w, n. The Darcy velocity of each phase satisfies where k rα is the relative permeability of the α phase, which depends on the water saturation S w , K is the absolute rock permeability, µ α and ρ α are the viscosity and density of the α phase (here assumed constant), g is the gravitational constant, and ∇z defines the direction of gravity. Volume balance requires the algebraic constraint and conservation of mass requires where φ is the porosity and q α are the source or sink (i.e., well) terms. One also needs the capillary pressure relation One typical capillary pressure function is given by where B c is a positive parameter inversely proportional to √ k and the normalized saturation S e is given by where S rα are the residual saturations for the wetting and non-wetting phases. An simple example of the relative permeability of each phase is given by where β is a positive parameter. Other commonly used rock models are those of Brooks-Corey and van Genuchten [15,19]. For a classical formulation, introduce the phase mobilities the total mobility and the fractional flow functions Sum the mass conservation equations (36) and combine with the algebraic constraint (35) to obtain the flow equation where u t = u w + u n is the total flow velocity, and where q t = q w + q n , (ρλ ) t = ρ w λ w + ρ n λ n . Therefore, the phase velocity u w and u n are related to the total velocity by The saturation equation for the water phase becomes We have three unknowns {p w , u t , S w } in the three equations (44), (45), and (48), which together comprise a well defined system (see, e.g., [2]). The classical formulation is challenged when the capillary pressure is degenerate or discontinuous. Suppose the water saturation is near the wetting phase residual saturation S rw . Then the non-wetting phase mobility is about 1/µ n , and at the same time the capillary pressure and its derivative tends to infinity, i.e., p c (S rw ) → ∞ and p c (S rw ) → ∞. Thus the term λ n K∇p c in (45) also tends to infinity. If we choose the non-wetting phase pressure p n as our primary variable instead of p w , we can avoid this singularity around S rw . However, if the derivative of the capillary pressure also tends to infinity at the non-wetting phase residual saturation, e.g., in a van Genuchten model, we will meet the same problem even with the non-wetting phase pressure as the primary variable. Moreover, when there are different rock types present, i.e., capillarity is heterogeneous as in (38), the spatial gradient of the capillary pressure ∇p c also tends to infinity.

The Hoteit-Firoozabadi formulation
In [27,26], Hoteit and Firoozabadi presented a new formulation that avoids the drawbacks of the classical one. They use the classic flow potential which leads to the capillary potential The total velocity u t is then written in terms of two velocity variables u a and u c , as follows: The velocity variable u a has the same pressure driving force as the wetting phase velocity u w , but it has a smoother mobility λ t rather than λ w . Since the total mobility λ t is strictly positive, the value of λ −1 t is bounded. The wetting phase velocity is then simply We rewrite the flow equation with the new splitting and the saturation equation in terms of u a to obtain

Approximation of flow and construction of the capillary flux
For ease of exposition, we collect the equations from the previous section that govern flow, and we add appropriate boundary conditions. At a fixed time t, where ∂ Ω = Γ N ∪ Γ D has been decomposed into disjoint Neumann and Dirichlet parts of the boundary andΦ B and u B are given. We assume that γ ∈ E h lies within either Γ D or At each time step, we will solve the system for the unknowns u a and Φ w , assuming that u c and the wetting phase saturation S w are given.

A mixed method for the flow equation
For denote the mixed finite elements used, which in our case are AC elements of index s ≥ 0, as described in Section 2.1 (although one could use any of the direct mixed finite elements defined in [7]). If one restricted to rectangular meshes, one could use the classic Raviart-Thomas (RT) elements [37]. We solve (57)-(60) using the hybrid form of the mixed method [12,17]. To this end, we define the Lagrange multiplier space M h to be the set of piecewise polynomials of degree up to s defined on each of the skeleton E h . Then the global mixed spaces are Note that we do not enforce H(div) conformity on V h . The hybrid mixed finite element formulation is: Find Since λ −1 t is always bounded, the first term in (61) is welldefined. The approximation spaces V h and W h are defined elementwise and the continuity is only enforced through the Lagrange multiplierΦ w on the skeleton of the mesh T h . In terms of the DoFs of the solution, we obtain the following algebraic problem where A λ is symmetric positive definite and block diagonal, with block size equal to dim V h (E). Using the Schur complement technique, we can express and therefore If the lowest order RT spaces are used, (65)-(67) is equivalent to equations used in the work of Hoteit and Firoozabadi [27, equations (17), (20), (25)].

Construction of the capillary flux
The capillary flux u c is defined by (53), i.e., u c = −λ n K∇Φ c . We have assumed the saturation S w is given, so Φ c (S w ) is known inside each element. However, S w is discontinuous on γ ∈ E h , so λ n is not known on the skeleton. Moreover, it vanishes at S w = 1 − S rn , so we cannot invert it. A weak form of the equation is to solve it locally on each E ∈ T h as follows: Find u c,h ∈ V h andΦ c,h ∈ M h such that and impose continuity of the capillary flux whereλ n is the interface value of λ n (S w ). Note that in this model, the capillary flux is set to zero on ∂ Ω . Following [27], the matrix form of the equations can be written as and we can eliminate the variable u c,h through a Schur complement, i.e., Knowing Φ c , we solve the latter equation forΦ c,h and construct u c,h from the former equation. Equation (72) is equivalent to [27, equation (31)]. The matrix L T A −1 L λ is singular if λ n = 0, and in that case, the linear system (72) cannot be solved. In [27], an upstream strategy is applied to resolve the singularity. However, if the whole region is filled with the wetting phase (S e = 1), we can not borrow a nonzero λ n from any neighboring elements. In order to solve the degenerate equation, we provide an alternate approach. Letζ h =λ nΦc,h , and solve in place of (68):

Now the matrix form of this with (69) is
and so The Schur complement L T A −1 L is now symmetric and does not depend on the value of saturation. The solution can be shifted by a constant and the differential equation still holds, so we need to set a constraint on the average or assign a Dirichlet boundary condition on a small portion of the boundary. Then the Schur complement system is non-singular and can always be solved. Onceζ h is found, we can construct u c,h by (75) and finally we obtain the capillary source B T u c,h needed in (62). We remark that in [36] and [5], the authors enforce the continuity of the capillary pressure p c (or Φ c with the gravity igonored) by penalizing the interface jump in p c . In this work and in [27,26], weak continuity of the capillary flux u c is enforced across the face, and no penalty parameter is required.

Approximation of transport with entropy stabilization
The saturation equation (56) will be solved assuming that u a is known. In that case, it is of hyperbolic type for the saturation. For simplicity, we impose Neumann conditions on the boundary ∂ Ω . To this end, decompose it into two disjoint parts, ∂ Ω = Γ in ∪ Γ out , where the inflow boundary is Γ in = {x ∈ ∂ Ω : u a ·ν ν ν < 0}, and the outflow boundary is Γ out = {x ∈ ∂ Ω : u a ·ν ν ν ≥ 0}. We impose the inflow boundary condition S w = S B on Γ in , where S B is given.
Because (56) is hyperbolic, numerical schemes need to be stabilized. We choose to stabilize our hyperbolic equation by adding an entropy viscosity µ h . The stabilized transport problem, posed for t > 0, becomes The inflow condition (78) corresponds to S w = S B for the unstabilized equation (56) (i.e., µ h = 0).

Enriched Galerkin for the saturation equation
Let subscript m = 0, 1, . . . denote the time step index. Within a transport time step (t m ,t m+1 ], we assume that we know the velocity u m a and the wetting phase saturation S m w at the earlier time t m .
We assume that the mesh skeleton E h is aligned with the interfaces between Γ in and Γ out . Let the set of boundary faces in E h be denoted by Γ h,in and Γ h,out . Suppose two elements E i and E j in T h are neighbors, so γ = ∂ E i ∩ ∂ E j ∈ E h is nonempty, and that the unit normal ν ν ν γ points from E i to E j . For a function φ which is piecewise continuous on each E ∈ T h , we define the average and jump on γ as Let the upwind value of the interface saturation be Define the interior penalty term J σ mapping to R as where σ γ > 0 is the penalty parameter for face γ ∈ E h , r is the degree of the polynomials used in the finite element space, d = 2 is the dimension of the domain, and |γ| is the measure of γ. Let where µ m h will be defined in the next subsection and s form ∈ {−1, 0, 1} determines whether one chooses to use the nonsymmetric (NIPG, s form = −1), incomplete (IIPG, s form = 0), or symmetric (SIPG, s form = 1) interior penalty Galerkin formulation. Let the linear functional L m DG (w) for the time step m be Finally, after applying backward Euler for the time discretization, we have the fully discrete Galerkin formulation: Given where δt m+1 = t m+1 − t m and X h (T h ) can be any of the discontinuous Galerkin finite element spaces defined in (27)- (29) or any of the enriched spaces defined in (32)-(33).

Entropy viscosity
We review the entropy viscosity method of Guermond, Pasquetti, and Popov [25] for the conservation law An entropy pair is a set of functions E(c) and F(c) such that E(c) is convex and F (c) = E (c)F (c). The function E(c) is called an entropy, and F(c) is called the associated entropy flux. A commonly used choice is E(c) = 1 2 c 2 . Then the entropy solution of (89) satisfies with equality holding except at the shocks. Given a numerical approximation c h (·,t) at time t, we define the entropy residual R h (·,t) = ∂ t E(c h (·,t)) + ∇ · F(c h (·,t)). (91) Following [25,50,18,35], we use this residual to define a viscosity µ E on each E ∈ T h as where λ E is a tunable parameter, h E = diam(E), · ∞ is the maximum norm, the denominator is a normalization constant, andĒ(c h ) = 1 |Ω | Ω E(c h ). However, in our applications, the (true) saturation is always between 0 and 1, and we intend to take the entropy E(c) = 1 2 c 2 , so the normalization E(c h ) −Ē(c h ) ∞,Ω ≤ 1/2. Therefore we can simplify the entropy viscosity (92) tõ We introduce an upper bound for some tunable parameter λ max , and set the entropy viscosity as The original problem is modifed to Now for our equation (56), we adapt the ideas as follows.
We define the fully discrete entropy residual locally on an element E ∈ T h (so the term ∇S m w,h is welldefined). In fact, we drop the final term, since ∇ · u m a = 0 except at wells. We defineμ m h using (93), evaluating the maximum of |R m h | by considering only the quadrature points used in the numerical integration formulas. Similarly, we evaluate (94), which is now and then we can set µ m h as in (95).

IMPES coupling of flow and transport
We summarize the solution procedure, which is of the implicit pressure, explicit saturation (IMPES) type [19]. The wetting phase saturation S 0 w,h ∈ X h (T h ) is given as an approximation of the initial condition S 0 in (80), perhaps as the L 2 -projection as stated in (88). We set S −1 w,h = S 0 w,h to start the algorithm.   [34].

Extension to three dimensions
Our results extend to three dimensions, that is, to meshes of cuboidal hexahedra. New families of mixed finite elements that approximate the velocity and pressure accurately were defined by the current authors in [8], and they have the minimal number of DoFs for general hexahedra. One could also use the macro-elements of Cockburn and Fu [21]. Unfortunately, at this writing, it is not known how to define serendipity elements that maintain accuracy on general hexahedra. However, direct serendipity elements were defined in the Ph.D. dissertation of Tao [42] for special hexahedra, namely, truncated cubes. A truncated cube is a cuboidal hexahedron with two pairs of parallel faces (so the cross-section is parallelogram).
The dimension of the direct serendipity element of index r on a truncated cube is dim P r + 3(r + 1). Table 2 shows the number of DoFs for an n × n × n mesh. Compared with using fully discontinuous spaces or standard mapped tensor product elements, the reduction in the number of degrees of freedom is much more significant than in two dimensions. The enhanced Galerkin method is also dramatically better on hexahedra than using tetrahedra, since it takes at least five tetrahedra to fill a fixed hexahedron. Table 2: A comparison of the global number of DoFs of DG-Q r , DG-P r , EG-Q r , and EG-DS r on an n × n × n mesh.

Numerical tests
In this work, the absolute permeability tensor is chosen as K = kI, where I is the identity matrix and k > 0 is possibly heterogeneous over the domain. We consider three numerical test cases. They all use s form = −1, i.e., the NIPG variant of the method. We implemented the numerical methods within the deal.II framework [14].

Coupled flow and transport
We begin with a simple test of the numerical procedure, taken from [41], except that we use quadrilateral or hexahedral mesh elements instead of triangles. Let λ t = 1, f w (S w ) = S w , and ignore capillary and gravity forces. The model reduces to a coupled linear flow and transport model that is used to represent single phase flow of the concentration c = S w of a dilute solute. That is, we have We take φ = 1, set k = 1, and assume that there are no injection nor extraction sources within the domain (i.e., q = q c = 0).

Two dimensions
The problem is solved in a unit square with a circular hole in the middle. The flow equation (99)-(100) is given a Dirichlet boundary condition p = 1.5 at the left side and p = 0.5 at the right, and the top and bottom sides, as well as the circular hole, are given a no flow Neumann boundary condition (u · ν ν ν = 0). The transport equation (101) is given an inflow boundary condition c = 1 on the right side (i.e., on Γ in ), and the initial condition is c = 0. Figure 3 is an example of the mesh at refinement level 2 with 128 convex quadrilateral elements. The actual simulation is conducted with a similar mesh at refinement level 5 with 8192 convex quadrilateral elements. The velocity u is approximated using the first order full H(div)-approximation Arbogast-Correa space (u h ∈ AC 1 ) [3], which means that the pressure is approximated by a discontinuous polynomial of degree 1 (i.e., p h ∈X P 1 ). In Figure 4, we show the pressure and velocity for the flow problem (99)-(100).
The transport equation (101) is solved using the EG-DS 2 method, i.e., we take X h (T h ) = X DS 2 (T h ) in (87). The time step is δt = 0.1h min = 0.00105998, where h min is the minimal cell diameter of all elements in the mesh T h . The simplified version of entropy viscosity is applied. The results are shown in Figure 5, where the left column show the concentration values at time step m = 200, 500, 800, and 1000 when only the maximum stabilization is applied (µ h = µ max ). The right column of Figure 5 show the concentration values at the same time steps with the entropy stabilization (µ h = min(µ max ,μ E )). It is obvious that the concentration values with the entropy stabilization show a sharper front.
Our numerical results compare favorably in accuracy to those of Sun and Liu [41]. The advantage to our approach is that we use many fewer DoFs. We achieve this by using quadrilateral elements and serendipity-based enriched Galerkin finite element spaces rather than triangular meshes.  This test was also performed in [34], and again we see comparable accuracy. However, that work was constrained to use rectangular meshes for accuracy, and so the internal hole was modified to be a rectangle. We were not subject to this restriction, since our methods are accurate on quadrilateral meshes and we can follow a curved interface with them (recall Figure 3). with the entropy stabilization (µ h = min(µ max ,μ E )). We can observe that the concentration values with the entropy stabilization show sharper transport fronts and less numerical diffusion.

Example 1: Effect of capillary pressure
In this example, we show the effect of capillary pressure in a heterogeneous medium. The test is modified from that given in [27] in that we use quadrilateral mesh elements instead of rectangles.
This is a nondimensionalized model problem, so all physical units are ignored in the following description. The computational domain is Ω = (0, 1.25) × (0, 0.875). It is composed of layers of alternate permeabilities k = 0.5 and k = 0.01, as shown in Figure 9, and rock porosity φ = 0.2. The fluid viscosities are µ w = 1 and µ n = 0.45. The relative permeabilities are set as in (40) with β = 2, the residual saturations are S rw = 0 and S rn = 0 (so S e = S w ) and the capillary pressure is  for the non-zero capillary pressure simulation, otherwise, it is set to be zero. We ignore gravity in this test, and therefore the densities are not required. There are no injection and extraction sources inside the domain (q w = q n = 0). For the flow equation, we set a Dirichlet boundary condition Φ w = 1 on the left hand side and Φ w = 0 on the right hand side, and a no flow Neumann boundary condition on the top and bottom sides of the domain. The inflow boundary for the saturation equation is the left hand side, and S B = 1 there. The domain is initially saturated with oil (the non-wetting phase), so S 0 = 0.
The computational mesh T h consists of randomly distorted quadrilaterals, as shown in Figure 9. The finite element approximation spaces used are u a,h ∈ AC 1 , Φ w,h ∈ X P 1 , and S w,h ∈ X DS 2 . The time step is δt = 0.001, and we use NIPG for the entropy stabilization (σ γ = 1). The entropy function E(S w ) = 1 2 S 2 w and F (S w ) = S w f w (S w )u a , and λ max = 0.2 and λ E = 1.
In Figure 10, we compare the simulated wetting phase saturation with and without capillary pressure at the same time steps. We observe that the front of the saturation profile with capillary pressure slows down and is more diffuse compared to the zero capillary pressure case. Moreover, the velocity is much greater near the interface between different rock types due to the presence of heterogeneous capillarity, as similar numerical tests in [36,27,5,32] show.
Our results compare favorably to those in [27], where rectangular elements were used. The test shows that grid distortion does not degrade the accuracy of the results.

Example 2: Heterogeneous capillary pressure
In this second example of two-phase flow, we show the difference of homogeneous and heterogeneous capillary pressure in a heterogeneous medium. Similar to the previous example, this is a nondimensionalized model problem so that all physical units are omitted. The computational domain is Ω = (0, 1) 2 . A random permeability field is distributed over Ω as depicted in Figure 11. The rock porosity φ , the fluid viscosities µ w and µ n , the relative permeabilities, the residual saturations, and initial saturation S 0 are the same as in Example 1, and gravity is ignored here as well. The heterogeneous capillary pressure is given in (102), and the homogeneous capillary pressure is set as which does not depend on the rock permeability. There are no injection and extraction sources inside the domain (q w = q n = 0). The water (wetting phase) injection boundary is {0} × (0, 0.05), where Φ w = 1, and the extraction boundary is {1} × (0.95, 1), where Φ w = 0. The rest of the flow boundary is set with a no flow Neumann boundary condition. The inflow boundary is the injection boundary, and S B = 1 there.
The computational mesh has 100 × 100 distorted quadrilateral elements (as in the coarse mesh of Figure 12). The time step, approximation spaces, entropy stabilization and the NIPG formulation are the same as Example 1, except λ E = 1.5 in this example. Figure 13 shows the profiles of the wetting phase saturation with homogeneous and heterogeneous capillary pressures at the same time steps. At time step m = 1000, there are two low permeability regions that are not invaded by the wetting phase if a homogeneous capillary pressure is applied. On the other hand, with a heterogeneous capillary pressure, the wetting phase invades these regions. At time step m = 2000, the water front has almost reached the extraction boundary. The homogeneous capillary pressure case shows an unrecovered region near the bottom right of the domain, while that region is much better saturated if the heterogeneous capillary pressure is applied.
Being able to use quadrilateral meshes greatly aids the set-up of this problem. It allows us to accurately model the shape of the low permeability inclusions. We maintain logically rectangular mesh indexing, and we use many fewer DoFs than triangles would provide.

Summary and conclusions
We described recently developed mixed and direct serendipity finite elements [3,7]. These elements provide accurate approximations on meshes of quadrilaterals and maintain the minimal number of degrees of freedom (DoFs) needed to maintain finite element conformity. Elements that approximate to any order are available. We developed a new enriched Galerkin method [41] on quadrilaterals based on the direct serendipity finite elements.
We applied the new elements to computational simulation of two-phase flow in porous media. We used the Hoteit-Firoozabadi formulation [27], which separates the system in   terms of the advective and capillary velocities u a and u c into an elliptic pressure equation for the flow velocity u a and a hyperbolic saturation equation for S w . A standard IMPES solution procedure allows the two parts of the system to be solved independently.
The pressure equation was discretized on quadrilateral meshes using a mixed finite element approximation and the accurate Arbogast-Correa spaces [3]. However, the formulation requires construction of the divergence of the capillary flux. We provided a novel implementation that does not break down when the system degenerates (i.e., one of the saturations tends to the residual value).
The saturation equation was also discretized accurately on quadrilateral meshes using the new enriched Galerkin method employing the direct serendipity spaces of the current authors [7]. Being hyperbolic, the system needs to be stabilized, and we used the entropy stabilization procedure of Guermond, Pasquetti, and Popov [25], adapted to our specific equation.
Extension to three space dimensions is straightforward, up to the definition of the minimal DoF finite element spaces.
The mixed spaces have been defined in [8,21], and the direct serendipity spaces were defined in [42] for truncated cubes.
Numerical results showed that accurate results are obtained by our numerical method, in both two and three dimensions. The ability to use quadrilateral and hexahedral meshes greatly increases one's ability to set up appropriate meshes for the problem to be solved. Our method uses a minimal number of DoFs, and so is quite efficient. The low number of DoFs is due to two sources. First, quadrilateral and hexahedral meshes use far fewer elements than ones based on triangles or simplices. Second, we use minimal DoF finite element spaces, which greatly reduce the number of DoFs compared to, say ABF or Devloo et al. mixed finite element spaces [11,16,40] and standard enriched Galerkin methods based on tensor product polynomials. Moreover, the method is efficient because, in many problems, we can maintain a logically rectangular mesh indexing.
Our numerical results also showed that our novel construction of the divergence of the capillary flux captures well the effects of capillary pressure.