A Direct Mixed–Enriched Galerkin Method on Quadrilaterals for Two-phase Darcy Flow

A Direct Mixed–Enriched Galerkin Method on Quadrilaterals for Two-phase Darcy Flow

Todd Arbogast Todd Arbogast University of Texas at Austin; Department of Mathematics, C1200; Austin, TX 78712-1202 and Institute for Computational Engineering and Sciences, C0200; Austin, TX 78712-1229 22email: arbogast@ices.utexas.eduZhen TaoUniversity of Texas at Austin; Institute for Computational Engineering and Sciences, C0200; Austin, TX 78712–1229 44email: taozhen.cn@gmail.com    Zhen Tao Todd Arbogast University of Texas at Austin; Department of Mathematics, C1200; Austin, TX 78712-1202 and Institute for Computational Engineering and Sciences, C0200; Austin, TX 78712-1229 22email: arbogast@ices.utexas.eduZhen TaoUniversity of Texas at Austin; Institute for Computational Engineering and Sciences, C0200; Austin, TX 78712–1229 44email: taozhen.cn@gmail.com
Abstract

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.

Keywords:
Mixed method Enriched Galerkin AC elements Serendipity elements Capillary flux Entropy viscosity
Msc:
79S05 65M60 65M08

1 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 two-phase 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 , and on hexahedra, only or 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 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) ERW_1984 (); Boffi_Brezzi_Fortin_2013 (), the enhanced velocity method Wheeler_Wheeler_Yotov_2002 (); Thomas_Wheeler_2010 (), the multipoint flux approximation methods (MPFA) and the multipoint flux mixed finite element methods (MFMFE) AEKWY_2007 (); WXY_2012 (); Wheeler_Yotov_2006_multipointFlux (); Ingram_Wheeler_Yotov_2010 (), the mimetic finite difference methods HMSS_2002 (), and discontinuous Galerkin (DG) methods Wheeler_1978 (); BBO_1999 (); RWG_2001 (); Riviere_Wheeler_Girault_1999 (); Houston_Schwab_Suli_2002 (). 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 ABF_2002 (); Bergot_Durufle_2013 (); Siqueira_Devloo_Gomes_2013 (), which have many excess DoFs). Recently, Arbogast and Correa Arbogast_Correa_2016 () 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 Arbogast_Tao_2018_atSpaces () and Cockburn and Fu in Cockburn_Fu_2017_mDecompIII () 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 Arbogast_Tao_2018x_serendipity ().)

Many methods have also been devised for approximation of the transport problem, including some of the methods mentioned above, ELLAM and characteristic methods Douglas_Russell_1982 (); ERW_1984 (); Arbogast_Wheeler_1995 (); WLELQ_2002 (); Arbogast_Huang_2006 (); WZET_2013 (), the standard continuous Galerkin (CG) finite element method, and the enriched Galerkin (EG) method Sun_Liu_2009_eg (); Lee_Lee_Wheeler_2016_eg (). 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 Sun_Liu_2009_eg () 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 Lee_Lee_Wheeler_2016_eg (). The EG method has been applied to miscible displacement problems and two-phase flow in porous media in recent works Lee_Wheeler_2018 (); Lee_Wheeler_2017_adaptive ().

We will base our finite element space for EG on serendipity finite elements Ciarlet_1978 (); Arnold_Awanou_2011 (), 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 Arnold_Awanou_2011 (); Kaliakin_2001 (); Lee_Bathe_1993 (). Recently, the current authors introduced new, direct serendipity finite elements Arbogast_Tao_2018x_serendipity () 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 Hoteit_Firoozabadi_2008 (); Hoteit_Firoozabadi_2008_fractures (), solved using an IMPES operator splitting algorithm CHM_2006 (). 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 Guermond_Pasquetti_Popov_2011_entropyVisc (); Bonito_Guermond_Popov_2014 (); ZGMP_2013 () 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 67, respectively. Numerical tests are given in Section 8, and we summarize and conclude our results in the final section.

2 Direct finite elements on quadrilaterals

In this section, we review the finite element spaces that we will use. Let be a polygonal domain with boundary , and let denote the outward unit normal vector on the boundary. We impose a conforming finite element mesh of quadrilaterals over the domain of maximal diameter . We will assume that the mesh is shape-regular Girault_Raviart_1986 (); Arbogast_Correa_2016 (), which ensures that the mesh elements are not highly elongated nor degenerate nearly to triangles. The outer unit normal to element is . The set of all interior edges or faces of is its skeleton, and it is denoted by . Fix a unit normal vector for each edge .

Let be the vector space of homogeneous polynomials of exact degree , and let denote the space of polynomials of degree up to . Let be the space of tensor product polynomials of degree up to . Finally, is the space of -dimensional polynomial vectors for which each component is in . We may specify that the domain of definition of the polynomials by writing, e.g., for domain .

At times we will need a reference element, so fix it to be the square . For , the standard bilinear map mapping vertices to vertices is bijective. This gives rise to the map , which maps a function to a function by the rule , where . We also have the Piola transform based on the bilinear map . It maps a vector to a vector , and it preserves the normal components Boffi_Brezzi_Fortin_2013 ().

2.1 The AC mixed finite element spaces

Arbogast and Correa developed two families of mixed finite elements on quadrilateral meshes Arbogast_Correa_2016 () for approximating solving a second order elliptic equation in mixed form

 u=−a∇p,∇⋅u=fin Ω,u⋅ν=0on ∂Ω, (1)

where and the tensor is uniformly positive definite and bounded. For and index , the AC elements are defined in terms of polynomials and a supplemental space of functions . In terms of the vectors , the family of full -approximation elements are

 Vs\rm AC(E)=P2s⊕x~Ps⊕S\rm ACs(E) and Ws% \rm AC(E)=Ps, (2)

while, the family of reduced -approximation elements are, for index ,

 Vs,\rm red\rm AC(E)=P2s⊕S\rm ACs(E) and Ws,\rm red% \rm AC(E)=Ps−1. (3)

The vector elements merge together to form conforming spaces over (i.e., the normal components of the vectors are continuous across each ), while the scalar elements remain discontinuous on . The former family approximates , , and to order , while the latter family approximates to order and and or order , 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 conformity.

Using the Piola transformation , the space of supplemental vectors on is

 S\rm ACs(E)={span{σs1,σs2},s≥1,span{σ0},s=0, (4)

where , , and

 ⎧⎪ ⎪⎨⎪ ⎪⎩ ^σs1=curl(^xs−1(1−^x2)^y),s≥1,^σs2=curl(^x^ys−1(1−^y2)),s≥1,^σ0=curl(^x^y),s=0, (5)

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 Arnold_Brezzi_1985 (); Boffi_Brezzi_Fortin_2013 ().

2.2 Direct serendipity elements

The current authors defined families of direct serendipity finite elements in a recent paper Arbogast_Tao_2018x_serendipity (). On , the element of index takes the form

 \rm{DS}r(E)=Pr⊕S\rm{DS}r(E), (6)

and each family of elements is determined by the choice of the two supplemental functions spanning . The space can be merged into continuous (i.e., 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\rm{DS}r(E)=dimPr+2=12(r+2)(r+1)+2. (7)

A very general and explicit construction of these supplements is given in Arbogast_Tao_2018x_serendipity (). They can be defined directly on , or they can be defined on and mapped to . In this work, we recall and use one of the simplest families of direct serendipity elements.

When , it is known that the standard space of bilinear polynomials defined on and mapped to an element maintains accuracy Arnold_Awanou_2011 (). These have the form (6), with

 S\rm{DS}1(E)=span{FE(^x^y)}. (8)

Henceforth, we only describe the new direct serendipity finite elements for indices .

For the element , let denote the outer unit normal to edge (denoted ), , and identify the vertices as , , , and (see Fig. 2). We define the linear polynomial giving the distance of to edge  in the normal direction as

 ℓi(x) =−(x−xi)⋅νi,i=1,2,3,4, (9)

where is any point on the edge. The functions are positive over the interior of .

When , we can define the shape functions associated to each vertex as

 ϕv,13(x) =ℓ2(x)ℓ4(x), ϕv,14(x) =ℓ2(x)ℓ3(x), (10) ϕv,23(x) =ℓ1(x)ℓ4(x), ϕv,24(x) =ℓ1(x)ℓ3(x).

Interior shape functions appear when . In this case, we take the shape functions

 ϕE,j⊂ℓ1ℓ2ℓ3ℓ4Pr−4,j=1,…,12(r−2)(r−3) (11)

so that they span the entire space . Since these are internal to , the precise choice of basis is not so important.

The most interesting shape functions are those associated to the edges. Let

 ℓH=ℓ3−ℓ4andℓV=ℓ1−ℓ2, (12)

and define rational functions

 RV(x)=ℓ1(x)−ℓ2(x)ξ−1Vℓ1(x)+η−1Vℓ2(x), (13) RH(x)=ℓ3(x)−ℓ4(x)ξ−1Hℓ3(x)+η−1Hℓ4(x) (14)

(note that the denominators do not vanish on ), where , , and

 ξ−1V=√1−(νH⋅ν1)2, η−1V=√1−(νH⋅ν3)2, (15) ξ−1H=√1−(νV⋅ν2)2, η−1H=√1−(νV⋅ν4)2. (16)

Therefore,

 RV(x)|e1 =−ηV and RV(x)|e2 =ξV, (17) RH(x)|e3 =−ηH and RH(x)|e4 =ξH. (18)

There are shape functions associated to the edges and , and they are

 ϕH,j(x)=ℓ3(x)ℓ4(x)ℓjH(x), j=0,1,…,r−2, (19) ϕH,r−1+j(x)=ℓ3(x)ℓ4(x)to0.0pt$ℓV(x)ℓjH(x),$ (20) j=0,1,…,r−3, ϕH,2r−3(x)=ℓ3(x)ℓ4(x)to0.0pt$RV(x)ℓr−2H(x).$ (21)

In a similar way, shape functions associated with edges and are

 ϕV,j(x)=ℓ1(x)ℓ2(x)ℓjV(x), j=0,1,…,r−2, (22) ϕV,r−1+j(x)=ℓ1(x)ℓ2(x)to0.0pt$ℓH(x)ℓjV(x),$ (23) j=0,1,…,r−3, ϕV,2r−3(x)=ℓ1(x)ℓ2(x)to0.0pt$RH(x)ℓr−2V(x).$ (24)

The edge shape functions are regular polynomials of degree  except the last two functions in each direction, which are rational functions. However, all shape functions restrict to polynomials of degree  on the edges.

Finally,

 S\rm{DS}r(E)=span{ℓ3ℓ4ℓr−2HRV,ℓ1ℓ2ℓr−2VRH}, (25)

and

 \rm{DS}r(E) =span{ϕv,13,ϕv,14,ϕv,23,ϕv,24, (26) ϕH,j,ϕV,j (j=0,1,…,2r−3), ϕE,k (k=1,…,12(r−2)(r−3))} =Pr(E)⊕S\rm{DS}r(E).

The unisolvence of the direct serendipity finite element space is proved generally in Arbogast_Tao_2018x_serendipity (), or specifically for this choice of , , and in Arbogast_Tao_2017_serendipity (). A nodal basis can be constructed using local linear algebra as described in Arbogast_Tao_2017_serendipity (); Arbogast_Tao_2018x_serendipity ().

2.3 Enriched direct serendipity spaces

Let the discontinuous finite element spaces of order over be

 ¯XDSr(Th) ={ϕ∈L2(Ω):ϕ|E∈\rm{DS}r(E),E∈Th}, (27) ¯XQr(Th) ={ϕ∈L2(Ω):ϕ|E=FE(^ϕ), (28) ^ϕ∈Qr(^E),E∈Th}, ¯XPr(Th) ={ϕ∈L2(Ω):ϕ|E∈Pr(E),E∈Th}. (29)

We remark that are the spaces of piecewise constants on each element. The continuous finite element spaces of order over are

 ~XDSr(Th) =¯XDSr(Th)∩C0(Ω), (30) ~XQr(Th) =¯XQr(Th)∩C0(Ω), (31)

where is the set of continuous functions over . Finally, we define the enriched finite element spaces of order over as

 XDSr(Th) =~XDSr(Th)+¯XQ0(Th), (32) XQr(Th) =~XQr(Th)+¯XQ0(Th). (33)

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-, DG-, EG- and EG-, using, respectively, the elements , , , and . In Table 1, we show the number of degrees of freedom of these four methods for an two dimensional quadrilateral mesh. All the methods achieve the same convergence rate. The EG- methods ultilize the fewest number of degrees of freedom.

3 Two-phase flow formulation

Let subscript denote the wetting phase and subscript the non-wetting phase. The pressures, Darcy velocities, and saturations of each phase are denoted , , and , respectively, . The Darcy velocity of each phase satisfies

 uα=−krα(Sw)μαK(∇pα−ραg∇z),α=w,n, (34)

where is the relative permeability of the phase, which depends on the water saturation , is the absolute rock permeability, and are the viscosity and density of the phase (here assumed constant), is the gravitational constant, and defines the direction of gravity. Volume balance requires the algebraic constraint

 Sw+Sn=1, (35)

and conservation of mass requires

 ϕ∂Sα∂t+∇⋅uα=qα,α=w,n, (36)

where is the porosity and are the source or sink (i.e., well) terms. One also needs the capillary pressure relation

 pc(Sw)=pn−pw. (37)

One typical capillary pressure function is given by

 pc(Sw)=−Bclog(Se), (38)

where is a positive parameter inversely proportional to and the normalized saturation is given by

 Se=Sw−Srw1−Srw−Srn, (39)

where are the residual saturations for the wetting and non-wetting phases. An simple example of the relative permeability of each phase is given by

 krw=Sβeandkrn=(1−Se)β, (40)

where is a positive parameter. Other commonly used rock models are those of Brooks-Corey and van Genuchten Bear_1972 (); CHM_2006 ().

For a classical formulation, introduce the phase mobilities

 λα=krαμα,α=w,n, (41)

the total mobility

 λt=λw+λn, (42)

and the fractional flow functions

 fα=λαλt,α=w,n. (43)

Sum the mass conservation equations (36) and combine with the algebraic constraint (35) to obtain the flow equation

 ∇⋅ut=qt, (44)

where is the total flow velocity, and

 ut=−λtK∇pw−λnK∇pc+(ρλ)tgK∇z, (45)

where , . Therefore, the phase velocity and are related to the total velocity by

 uw (46) un (47)

The saturation equation for the water phase becomes

 ϕ∂Sw∂t+∇⋅{Kfw(Sw)λn(Sw)(dpcdSw∇Sw (48) +(ρw−ρn)g∇z)+fw(Sw)u}=qw.

We have three unknowns in the three equations (44), (45), and (48), which together comprise a well defined system (see, e.g., Arbogast_1992a ()).

The classical formulation is challenged when the capillary pressure is degenerate or discontinuous. Suppose the water saturation is near the wetting phase residual saturation . Then the non-wetting phase mobility is about , and at the same time the capillary pressure and its derivative tends to infinity, i.e., and . Thus the term in (45) also tends to infinity. If we choose the non-wetting phase pressure as our primary variable instead of , we can avoid this singularity around . 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 also tends to infinity.

3.1 The Hoteit-Firoozabadi formulation

In Hoteit_Firoozabadi_2008 (); Hoteit_Firoozabadi_2008_fractures (), Hoteit and Firoozabadi presented a new formulation that avoids the drawbacks of the classical one. They use the classic flow potential

 Φα=pα+ραgz,α=w,n, (49)

which leads to the capillary potential

 Φc=Φn−Φw=pc+(ρn−ρw)gz. (50)

The total velocity is then written in terms of two velocity variables and , as follows:

 ut =ua+uc, (51) ua =−λtK∇Φw, (52) uc =−λnK∇Φc. (53)

The velocity variable has the same pressure driving force as the wetting phase velocity , but it has a smoother mobility rather than . Since the total mobility is strictly positive, the value of is bounded. The wetting phase velocity is then simply

 uw=λwλt(−λtK∇Φw)=fwua. (54)

We rewrite the flow equation with the new splitting and the saturation equation in terms of to obtain

 ∇⋅ua=qt−∇⋅uc, (55) ϕ∂Sw∂t=qw−∇⋅(fwua). (56)

4 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 ,

 ua =−λtK∇Φw in Ω, (57) ∇⋅ua =qt−∇⋅uc in Ω, (58) Φw =^ΦB on ΓD, (59) ua⋅ν =uB−uc⋅ν on ΓN, (60)

where has been decomposed into disjoint Neumann and Dirichlet parts of the boundary and and are given. We assume that lies within either or , and denote and . At each time step, we will solve the system for the unknowns and , assuming that and the wetting phase saturation are given.

4.1 A mixed method for the flow equation

For , let denote the mixed finite elements used, which in our case are AC elements of index , as described in Section 2.1 (although one could use any of the direct mixed finite elements defined in Arbogast_Tao_2018x_serendipity ()). If one restricted to rectangular meshes, one could use the classic Raviart-Thomas (RT) elements Raviart_Thomas_1977 (). We solve (57)–(60) using the hybrid form of the mixed method Arnold_Brezzi_1985 (); Boffi_Brezzi_Fortin_2013 (). To this end, we define the Lagrange multiplier space to be the set of piecewise polynomials of degree up to defined on each of the skeleton . Then the global mixed spaces are

 Vh ={vh:vh|E∈Vh(E) ∀E∈Th}, Wh ={wh:wh|E∈Wh(E) ∀E∈Th}, Mh ={μh:μh|γ∈Ps(γ) ∀γ∈Eh,γ∉ΓD}.

Note that we do not enforce conformity on .

The hybrid mixed finite element formulation is: Find , , and such that

 ∫Ω(λtK)−1ua,hvh−∑E∈Th∫EΦw,h∇⋅vh (61) +∑E∈Th∫∂E∖ΓD^Φw,hvh⋅ν∂E=−∫ΓD^ΦBvh⋅ν  ∀vh∈Vh, ∑E∈Th∫E∇⋅ua,hwh=∫Ω(qt−∇⋅uc,h)wh∀wh∈Wh, (62) ∑E∈Th∫∂E∖ΓDua,h⋅ν∂Eμh=∫ΓN(uB−uc⋅ν)μh (63) ∀μh∈Mh.

Since is always bounded, the first term in (61) is well-defined. The approximation spaces and are defined elementwise and the continuity is only enforced through the Lagrange multiplier on the skeleton of the mesh . In terms of the DoFs of the solution, we obtain the following algebraic problem

 ⎛⎜⎝Aλ−BLBT00LT00⎞⎟⎠⎛⎜ ⎜⎝ua,hΦw,h^Φw,h⎞⎟ ⎟⎠=⎛⎜⎝ΦBQUB⎞⎟⎠, (64)

where is symmetric positive definite and block diagonal, with block size equal to . Using the Schur complement technique, we can express

 ua,h=A−1λ(ΦB+BΦw,h−L^Φw,h), (65)

and therefore

 LTua,h =LTA−1λ(ΦB+BΦw,h−L^Φw,h)=Q, (66) BTua,h =BTA−1λ(ΦB+BΦw,h−L^Φw,h)=UB. (67)

If the lowest order RT spaces are used, (65)–(67) is equivalent to equations used in the work of Hoteit and Firoozabadi (Hoteit_Firoozabadi_2008, , equations (17), (20), (25)).

4.2 Construction of the capillary flux

The capillary flux is defined by (53), i.e., . We have assumed the saturation is given, so is known inside each element. However, is discontinuous on , so is not known on the skeleton. Moreover, it vanishes at , so we cannot invert it. A weak form of the equation is to solve it locally on each as follows: Find and such that

 ∫EK−1uc,hvh+∫∂E^Φc,h^λnvh⋅νγ (68) =∫EΦc(Sw)∇⋅(λn(Sw)vh)∀vh∈Vh(E),

and impose continuity of the capillary flux

 ∑E∈Th∫∂Euc,h⋅νγμh=0∀μh∈Mh, (69)

where is the interface value of . Note that in this model, the capillary flux is set to zero on .

Following Hoteit_Firoozabadi_2008 (), the matrix form of the equations can be written as

 (ALλLT0)(uc,h^Φc,h)=(BλΦc0), (70)

and we can eliminate the variable through a Schur complement, i.e.,

 uc,h (71) LTuc,h =LTA−1(BλΦc−Lλ^Φc,h)=0. (72)

Knowing , we solve the latter equation for and construct from the former equation. Equation (72) is equivalent to (Hoteit_Firoozabadi_2008, , equation (31)). The matrix is singular if , and in that case, the linear system (72) cannot be solved. In Hoteit_Firoozabadi_2008 (), an upstream strategy is applied to resolve the singularity. However, if the whole region is filled with the wetting phase (), we can not borrow a nonzero from any neighboring elements.

In order to solve the degenerate equation, we provide an alternate approach. Let , and solve in place of (68): Find and such that

 ∫EK−1uc,hvh+∫∂E^ζhvh⋅νγ (73) =∫EΦc(Sw)∇⋅(λn(Sw)vh)∀vh∈Vh(E).

Now the matrix form of this with (69) is

 (ALLT0)(uc,h^ζh)=(BλΦc0), (74)

and so

 uc,h =A−1(BλΦc−L^ζh), (75) LTuc,h =LTA−1(BλΦc−L^ζh)=0. (76)

The Schur complement 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 is found, we can construct by (75) and finally we obtain the capillary source needed in (62).

We remark that in Lee_Wheeler_2018 () and AJPW_2013 (), the authors enforce the continuity of the capillary pressure (or with the gravity igonored) by penalizing the interface jump in . In this work and in Hoteit_Firoozabadi_2008 (); Hoteit_Firoozabadi_2008_fractures (), weak continuity of the capillary flux is enforced across the face, and no penalty parameter is required.

5 Approximation of transport with entropy stabilization

The saturation equation (56) will be solved assuming that 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, , where the inflow boundary is , and the outflow boundary is