{'fraction_non_alphanumeric': 0.04260645872571087, 'fraction_numerical': 0.022646883182507235, 'mean_word_length': 4.390491519128389, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Maximum-principle-satisfying discontinuous Galerkin methods for incompressible two-phase immiscible flow

This paper proposes a fully implicit numerical scheme for immiscible incompressible two-phase flow in porous media taking into account gravity, capillary effects, and heterogeneity. The objective is to develop a fully implicit stable discontinuous Galerkin (DG) solver for this system that is accurate, bound-preserving, and locally mass conservative. To achieve this, we augment our DG formulation with post-processing flux and slope limiters. The proposed framework is applied to several benchmark problems and the discrete solutions are shown to be accurate, to satisfy the maximum principle and local mass conservation.


Introduction
Multiphase flows in porous media appear in a large number of applications in engineering and sciences, for instance in the environmental clean up of contaminated subsurface or in the energy production of hydrocarbons from reservoirs. This paper introduces a numerical method for solving the immiscible two-phase flow equations, that produces bound-preserving discrete saturations. The proposed method utilizes a fully implicit in time stepping scheme, a discontinuous Galerkin in space discretization and post-processing flux and slope limiters techniques. The resulting numerical saturation is shown to satisfy a maximum principle theoretically and computationally.
The numerical literature for immiscible two-phase flow problems is vast (see Peaceman [2000]; Aziz and Settari [1979]; Chen et al. [2006] and references herein). Suitable numerical methods should be locally mass conservative and should produce bound-preserving discrete saturations. Such methods include finite difference methods and finite volume methods, which are popular methods because of their simplicity and low cost Michel [2003]; Droniou [2014]. However, finite difference methods are not adapted to unstructured meshes and cell-centered finite volume methods suffer from grid distortion and do not easily handle full anisotropy. The class of interior penalty discontinuous Galerkin methods has been applied to model multiphase flows in porous media for more than fifteen years Klieber and Riviere [2006]; Riviere [2007, 2006]; Ern et al. [2010]; Arbogast et al. [2013]; Bastian [2014]; Jamei and Ghafouri [2016] and they have been combined with other locally mass conservative methods like mixed finite element methods in Hoteit and Firoozabadi [2008]; Hou et al. [2016]. DG methods are locally mass conservative, they do not suffer from grid distortion and they are accurate and robust even in the case of anisotropic heterogeneous media. However, it is well known that the DG approximation of the saturation does not satisfy a maximum principle because of local overshoots and undershoots in the neighborhood of the saturation front. While the amount of overshoot and undershoot can be reduced by the choice of implicit time stepping, mesh refinement and appropriate penalty parameters, there is no guarantee that they will completely disappear. The literature on post-processing techniques to reduce or eliminate the amount of overshoots and undershoots for DG methods in general is significant. Slope limiters adjust the gradient of the linear approximation in a heuristic way Burbeau et al. [2001]; Hoteit et al. [2004]; Krivodonova [2007]; Krivodonova et al. [2004]; Kuzmin [2010Kuzmin [ , 2013. Recently, flux limiters related to flux-corrected transport algorithms, were introduced for DG discretizations of conservation laws Frank et al. [2019]; Kuzmin and Gorb [2012].
The main contribution of this paper is the formulation of bound-preserving numerical method for the incompressible two-phase flow problems. Upwind fluxes are employed for the interior penalty discontinuous Galerkin discretization in space. We solve several benchmark problems to investigate the performance of the method and particularly the impact of the limiting techniques on local mass conservation. The numerical method respects maximum principle by limiting the saturation profile to physical upper-and lower-bounds. The violation of maximum principle for the discontinuous approximation of the saturation has been an open problem over the last decade. Our proposed scheme guarantees that the saturation remains bounded in the physical range. In addition, we observe that the monotonicity of the saturation is significantly improved compared to the case of no limiters. Saturation fronts are sharp with minimal numerical diffusion. We present several numerical results that show overshoots and undershoots have been eliminated. We verify that the local mass conservation property is also satisfied. We consider cases where flow is driven by boundary conditions and cases where flow is driven by injection and production wells. In the former case, a theoretical proof of the maximum principle is given.
The content of the paper is as follows. Section 2 describes the mathematical equations. The primary unknowns are the wetting phase pressure and saturation. Section 3 contains the fully implicit numerical scheme, with the construction and analysis of the flux limiters, and the review of slope limiters used in this work. Several numerical results, including benchmark problems and convergence tests, are given in Section 4. Conclusions follow.

GOVERNING EQUATIONS
The incompressible two-phase flow in a porous medium Ω ⊂ R 2 over a time interval [0, T ], is modeled by a system of mass balance equations for each phase, coupled with closure relations.
The other coefficients are the porosity φ, the absolute permeability K, and the gravity vector g. Using (2.2) and (2.3), and choosing for primary unknowns the wetting phase pressure and saturation (P, S) = (P w , S w ), the system of equations reduces to: Let the boundary of the domain be divided into different disjoint sets Dirichlet and Neumann boundary conditions are imposed on parts of the boundary: (2.10) The boundary Γ out,s is referred to as a free boundary because no data is prescribed on that boundary. This means that the surface integrals on this boundary are evaluated in terms of the unknowns. This particular treatment of the outflow boundary has been highlighted in the works of [Papanastasiou et al., 1992;Griffiths, 1997]. In the case of pure homogeneous Neumann boundary conditions (Γ N,s = Γ N,p = ∂Ω) and j p = j s = 0), the flow is driven by injection/production wells (source/sink functions) that depend on the wetting phase saturation as follows: The functionsq and q correspond to the injection and production well rates and s in is the prescribed wetting phase saturation at the injection wells. The fractional flow functions, f α , are the ratios of the phase mobility to the total mobility, f α = λ α /(λ + λ w ).
Finally the model problem is completed by an initial condition on the saturation: S = s 0 .

Numerical method
The domain Ω is decomposed into a non-degenerate partition E h = {E} E consisting of N h triangular or rectangular elements of maximum diameter h. Let Γ h denote the set of all edges and Γ i h denote the set of interior edges. For any e ∈ Γ i h , fix a unit normal vector n e and denote by E + and E − the elements that share the edge e such that the vector n e is directed from E + to E − . We define the jump and average of a scalar function ξ on e as follows: By convention, if e is adjacent to ∂Ω, then the jump and average of ξ on e coincide with the trace of ξ on e and the normal vector n e coincides with the outward normal n. Let P 1 (E) be the space of linear polynomials on an element E. The discontinuous finite element space of order one is: The time interval T is divided into N τ equal subintervals of length τ . Let P n and S n be the numerical solutions at time t n . The proposed discontinuous Galerkin scheme for equations (2.5)-(2.10) reads: The penalty parameter σ is constant on the interior edges and its value is chosen 10 times larger on the Dirichlet boundaries. The quantities (·) ↑v n and (·) ↑v n w denote the upwind values with respect to the vector functions v n and v n w that are scaled quantities of the phase velocities. They depend on the pressure and saturation evaluated at the previous time t n : The definition of the upwind operator with respect to a generic discontinuous vector field v is: At the initial time, the discrete saturation is the L 2 projection of the initial condition.
At each time step, we solve (3.3)-(3.4) together with a Newton solver, followed by flux and slope limiters (see Algorithm 1). Figure 1 is a schematic that describes the actions of both flux and slope limiters on the discrete saturation. The next two sections describe these limiters in detail and contain a proof that the resulting saturation is bound-preserving.
Algorithm 1 DG+FL+SL method Compute initial saturation S 0 for n = 0, . . . , (N τ − 1) do Solve (3.3)-(3.4) with Newton's method Apply flux limiter: S FL n+1 = L flux (S n+1 ) Apply slope limiter to S n+1 = L slope (S FL n+1 ) end for 3.1. Flux limiter. The flux limiter will enforce that the element-wise average of the saturation satisfies the desired physical bounds. We assume that the saturation at the previous time step, t n , satisfies: (3.5) for some constants 0 ≤ s * ≤ s * ≤ 1. The flux limiting is applied to each element E given the element-wise average of the saturation at the previous and current time steps and given a flux function defined on each face e ⊂ ∂E. First we compute the element-wise average at time t n and t n+1 : Next, for a fixed element E, let n E be the unit outward normal vector to E. We define the flux function H n+1 | E = H n+1,E as follows: For an interior face e of the element E, the quantity H n+1,E (e) measures the net mass flux across e into the neighboring element E that also shares the face e. We note that: After application of the flux limiter operator, the limited saturation has a possibly different cellaverage: (3.10) The new cell-average of the saturation is obtained by an iterative process, that takes for input the cell average at the previous time step and the flux function: H n+1 ). Before showing that the limited saturation satisfies (3.5), we describe the algorithm for the operator L avg .
3.1.1. The algorithm for L avg . For a fixed element E, we denote by N E the set of elements that include E and all neighboring elements E that share a face e with E. The algorithm constructs a sequence of flux functions and element-wise averages for E and its neighbors E . While the construction of the element-wise averages are local to E and its neighbors E , the stopping criterion is global to ensure bound-preserving solutions. We first initialize the sequences with the input arguments:S Next, for k ≥ 1, we have the following steps: Step 1 Compute inflow and outflow fluxes: Step 5. Define a global stopping criterion E . Else set k ← k + 1 and go to Step 1. 3.1.2. Bound-preserving solutions. In this section, we show that the solution S FL n+1 obtained in (3.10) has a cell-average that is bound-preserving for the case where flow is driven by boundary conditions only (no wells). Clearly, it suffices to show thatS FL n+1 is bound-preserving. This is done in two steps. First, we show that each iterate in the flux-limiter algorithm is bound preserving. Second, we show that the stopping criterion is reached for some value k 0 .
Lemma 3.1. Let E be a mesh element and let (S belongs to the interval [s * , s * ], it then follows by its definition that α (k−1) E (e) ≥ 0 for all e ⊂ ∂E. We then apply the inequality x ≤ max(0, x) to (3.14) to obtain: Therefore with the definition of Q + E we havē The proof for the lower boundS (k) E ≥ S * follows a similar argument, after applying the identity x ≥ min(0, x) to (3.14).
Lemma 3.2. Assume that the cell averages S n,E at time t n belong to the interval [s * , s * ] for all elements E. Then we have Proof. With Lemma 3.1, it suffices to show that the sequence (H This implies convergence of (H (k) E ) for all elements E, so that there exists k 0 such that 3.2. Slope limiter. The slope limiter operator, denoted by L slope , is applied to the discrete saturation S FL n+1 at each time step. The element-wise mean values of the saturation are left unchanged by this procedure. There is a variety of slope limiters available in the literature. For convenience, we choose a vertex-based slope limiter that is well suited for piecewise linear polynomials Kuzmin [2010] and that consists of two steps.
(i) We first mark the elements in which the maximum principle is not satisfied (i.e., S FL n+1 (x) > s * or S FL n+1 (x) < s * ). We will apply the slope limiter on these marked elements only. (ii) By a Taylor expansion around the centroid c E of element E, the linear saturation takes the form For the marked elements, a slope limiter replaces the local solution S FL n+1 | E by the following linear constrained reconstruction Let v E,i denote the i th vertex of element E. We determine the maximum admissible slope for the constrained reconstruction by choosing values β E ∈ [0, 1] such that boundedness of S n+1 is satisfied at all vertices of E: where S E,i * and S * E,i are defined as maximum and minimum means values of the saturation over all the elements (including E) that contain the vertex v E,i ∈ E.
The bounds of the saturation at all vertices are guaranteed if the correction factor β E is chosen as: (3.21) Using all the previous results, we obtain that the discrete saturation is bound-preserving.
We utilize the MPI-based PETSc library [Balay et al., 2017[Balay et al., , 2018Dalcin et al., 2011] as the linear algebra back-end to solve nonlinear equations (3.3)-(3.4). We use Newton's method with (damped) step line search technique [Crisfield, 1979] and set the relative convergence tolerance to 10 −6 . For the inner linear solve at each Newton iteration, we rely on the MUMPS direct solver [Amestoy et al., 2001[Amestoy et al., , 2019 with relative pivoting threshold of 0.01. MUMPS uses several efficient preordering algorithms to permute the columns of matrix and thereby minimize the fill-in (number of nonzeros in the factorization) in the LU factorization. At each time step, following the Newton solver convergence, we apply flux and slope limiters. Implementation of the flux limiter algorithm discussed in section 3.1 is provided in the module FluxLimiter along with an auxiliary flux wrapper module named Hsign. Global stopping criteria for all problem sets are taken as 1 = 2 = 10 −6 . As for the slope limiter, we use the native VertexBasedLimiter module embedded in the Firedrake project. All simulations are conducted on a single socket Intel i5-8257U node by utilizing a single MPI process.
Codes used to perform experiments in this paper are publicly available at msarrafj/LimitedDG [2021] repository. Firedrake and its component may be obtained from https://www.firedrakeproject.org/. For reproducibility, we also cite archives of the exact software versions used to produce results in this paper. All major Firedrake components have been archived on Zenodo/firedrake [2021]. This record collates DOIs for the components and can be installed following the instructions at https://www.firedrakeproject.org/download.html.

NUMERICAL RESULTS
In this section, several numerical experiments are carried out in following order: (i) We first validate our proposed method on two benchmark problems: one-dimensional Buckley-Leverett problem and two-dimensional Buckley-Leverett problem with gravity. Further, we investigate the convergence rates by using method of manufactured solutions and verify that the flux limiter preserves accuracy. (ii) We then perform various pressure-driven flow problems on structured and (e) ∈ [0, 1] at the discrete level and performing multiple correction cycles (see section 3.1). Note that this procedure is mass conservative and we observe that the decrease of the average value on element E brings about an increase on the average value of neighboring elements E . We subsequently apply a slope limiting procedure to produce pointwise bound-preserving solutions (see section 3.2). The correction factor β E ∈ [0, 1] at each element vertex v E,i takes the average values of neighboring elements as local bounds and determines the maximum admissible slope (figure 3 ○). Note that the local average values are left unchanged during the slope limiting. unstructured meshes, to study the efficacy of limiters on capturing high-accuracy wetting phase saturation profiles. (iii) The robustness of the scheme in the presence of injection and production wells is assessed using the quarter five-spot problem, with homogeneous and discontinuous highly varying permeability fields. For both pressure-driven flow problem and quarter five-spot problems, we examine the element-wise mass balance property associated with the limiters and highlight the capability of the saturation in satisfying the maximum-principle. (iv) Finally, we study the influence of gravity on the flows by testing our scheme with three different gravity numbers.
For all problems, we assume the following parameters unless otherwise specified: The residual saturations imply the physical lower and upper bounds for the saturation:  Buckley et al., 1942], also known as the frontal-advance equation, is a well-known non-linear hyperbolic equation for the description of one-dimensional immiscible displacement in a linear reservoir. Because the problem has a semi-analytical solution, it is widely used to validate numerical methods for two-phase flows in porous media. Since capillary pressure and gravity are neglected, the total velocity of the phases u t = (u t , v t ) T can be written as: (4.1) By substituting ∇P from (4.1) into equation (2.6) and ignoring source/sink terms, we obtain the general form of the Buckley-Leverett equation.
The convection flux f BL = (F (S), G(S)) T reduces in one-dimension to: , and G(S) = 0. (4. 3) The relative permeabilities are chosen as: We take an interval domain Ω = [0, 300] m with uniform mesh, and we fix the following parameters: The Dirichlet boundary condition of g s = 0.85 is weakly prescribed at the left boundary x = 0. We assume outflow boundary at x = 300. This setup gives rise to the classical Buckley-Leverett profile, which consists of a shock wave immediately followed by a rarefaction wave. Both lower and higher DG approximation of solution without any external bound-preserving mechanism do not respect maximum principle [Dawson et al., 2004;Zhang and Shu, 2011]. Here, we employ the first-order implicit DG formulation with our proposed limiter scheme to discretize equations (4.2)-(4.3) in space, and backward Euler scheme is utilized in time. The DG penalty parameter is set to σ = 10 −6 and the flux f BL is approximated with a first-order upwind method [Fambri, 2020; as it provides good results in conjunction with proposed limiters. To implement the flux limiter, the following flux functional H E (e) is adopted on each face e ⊂ ∂E: The final simulation time is T = 800 days, and the saturation profile is depicted in Figure 2 for t = 400 and t = 800 days. We performed a four-step mesh refinement study and linearly refined τ at each step. We observe that the location of the front obtained from the proposed numerical scheme is in good agreement with the location of the front for the analytical solution even for the case of coarse mesh and as we proceed with refinement, the discrete solution converges to the analytical solution. To calculate the semi-analytical solution of Buckley-Leverett equation (i.e., the position of the saturation front), we resorted to Welge graphical method [Welge et al., 1952].  also provides a zoom-in view at the location of front for t = 800 days for better visualization. As expected, the numerical saturation remains within physical bounds and no undershoots and overshoots are observed. The choice of implicit time marching algorithm is shown to have no erroneous smearing effect on the saturation front. 4.1.2. Buckley-Leverett problem with gravity. In this numerical experiment, we study the effect of proposed limiters in the two-dimensional Buckley-Leverett equation that incorporates gravitational effects along the y-axis [Zhang and Tang, 2002;Christov and Popov, 2008;de Loubens, 2007]. Consider equation (4.2) with the following non-convex flux functions in the x-and y-directions: This benchmark problem was solved by finite element method combined with operator-splitting method in [Karlsen et al., 1998]. For comparison purpose, hence, we take u t = v t = 1 m/s, λ w (S) = S 2 , λ (S) = (1 − S 2 ) and φ = 1. We solve (4.2) and (4.5) on the square domain [0, 3] 2 m 2 with structured triangular mesh of size h = 0.03 m subject to the initial condition: s 0 (x, y) = 1, for (x − 1.5) 2 + (y − 1.5) 2 < 0.5 0, otherwise.  Finally, we impose no flow condition u t · n = 0 everywhere on the boundary ∂Ω. Similar to the one-dimensional problem, we use backward Euler time marching and discretize the problem with implicit DG formulation (with σ = 0.1) augmented with the proposed flux and slope limiters scheme. Herein, the flux limiter functional H n+1,E (e) on the interior edges is the same as that of the one-dimensional Buckly-Leverett and on all exterior edges it is set to 0. The simulation runs to T = 0.5 s with 440 time steps. In Figure 3, we show the numerical results at the final time obtained from the implicit DG formulation without limiters (see Figure 3(a)) and with limiters (see Figure  3(b)). We compare the results with the reference solution. Evidently, DG scheme with no limiters produces an oscillatory solution that results in strong violations with respect to maximum principle. However, the application of limiters give rise to bound-preserving solution (i.e., 0 ≤ S n+1 ≤ 1); and undershoots and overshoots are eliminated completely. This result does not exhibit extra numerical diffusion and is consistent with the reference solution shown in Figure 3(c). 4.1.3. Convergence study. We carry out an h-convergence study on two-dimensional structured triangular meshes in order to verify convergence properties of our limiter scheme. The computational domain Ω is the unit square and the exact solutions are: s(x, y, t) = 0.4 + 0.4xy + 0.2 cos(t + x), (4.7a) p(x, y, t) = 2 + x 2 y − y 2 + x 2 sin(y + t) − 1 3 cos(t) + 1 3 cos(t + 1) − 11 6 . (4.7b) Through the method of manufactured solutions, we replace the source/sink terms (i. e., wells flow rates) of equations (2.5)-(2.6) by body force terms obtained from manufactured solutions. Dirichlet boundary conditions are applied on ∂Ω on both saturation and pressure fields. The input parameters are: φ = 0.2, K = 1 m 2 , µ w = µ = 1 Pa · s, s rw = s r = 0, k rw (S) = S 2 , k r (S) = (1 − S) 2 . Table 1. Errors in L 2 and H 1 norms and convergence rates, for τ = h 2 and T = 1 s. Note that flux limiter algorithm preserves the accuracy of the DG discretization. However, slope limiter slightly degrades the rates of convergence. The capillary pressure is based on Brooks-Corey model: where the entry pressure is set to p d = 50 Pa, inhomogeneity characterization parameter is set to θ = 2, and linearization tolerance is set to R = 0.05. The convergence properties are computed by using a time step τ set to h 2 . We note that the admissible global bounds for the flux limiter algorithm are updated throughout the simulation. In other words, at every time step, s * and s * bounds are determined by the maximum and minimum of the exact solution (4.7a), respectively. When using the limiters no upper and lower bound violations are observed in the discrete solution. Table 1 shows the errors in L 2 and H 1 norms evaluated at T = 1 s and the corresponding convergence rates for saturation and pressure. We compare rates for four cases: (i) no limiters (DG), (ii) with both flux and slope limiters (DG+FL+SL), (iii) with only slope limiter (DG+SL), and (iv) with only flux limiter (DG+FL). For both unknowns, DG returns expected optimal convergence rates of 2 in the L 2 norm and 1 in the H 1 norm. However, we observe that applying both limiters results in suboptimal rates. Cases (iii) and (iv) indicate that the application of flux limiters only preserve optimal rates whereas the application of slope limiters only yields a decline in the convergence rates. The slope limiter scheme taken from Kuzmin [2010] is completely independent of the proposed flux limiter in Section 3. Designing a slope limiter that produces optimal rates remains a challenge. We show in Table 2 the errors in the L 2 norm of the cell average for the saturation and the corresponding convergence rates. Optimal rate of 2 is obtained for either DG or DG+FL+SL. This result reiterates that the flux limiter does not reduce the accuracy and the slope limiter does not impact the rates since it does not alter the element-wise averages.
4.2. Two-dimensional pressure-driven flow. We take a computational domain of Ω = [0, 100] 2 m 2 with zero gravity field for all problems in this section. The wetting phase is injected along the left boundary and the non-wetting phase is pushed out through the right boundary. Dirichlet boundary conditions are: P = 3 × 10 6 Pa and S = 0.85 on {0} × (0, 100); and P = 10 6 Pa u ℓ · n = 0 on Γ N,p P = g p = 3 ×  A homogeneous test problem with constant permeability of K = 10 −8 m 2 is examined here, with similar setup and parameters as in the work of Epshteyn and Riviere [2007]. Relative permeabilities and capillary pressure are defined in equations (4.4) and (4.8), respectively, with entry pressure p d = 1000 Pa, θ = 2 and R = 0.05. The viscosities are µ w = 10 −3 Pa · s and µ = 10 −2 Pa · s. Two quadrilateral meshes are considered: (i) a uniform mesh with size of h = 1.25 m and (ii) a non-uniform mesh with 256 elements and with size of h bnd = 1.25 m at the left boundary and h = 6.583 m for the rest of domain (see Figure 5(a)). It is known that slope limiters by design flatten steep slopes near discontinuities (e.g., at left-most elements when simulation starts). Using a mesh with increased density at the (left) boundary reduces the effect of overflattening on the accuracy of solutions [May and Berger, 2013;Giuliani and Krivodonova, 2018]. The time step is chosen as τ = 0.2 s, the final time is T = 300 s, and the penalty parameter is σ = 100. We compare our numerical solutions with a reference unlimited  solution obtained from the fully implicit DG formulation developed by Epshteyn and Riviere [2007] on a quadrilateral mesh with 256 elements. The saturation and pressure profiles obtained with our proposed scheme, along the line y = 50 m are illustrated in Figures 6(a) and 6(b). Numerical solutions, compared to reference solution, are accurate and in good agreement with respect to front location. As expected, the finer mesh tracks the saturation front with more accuracy. It is also evident that our limiting scheme successfully yields pointwise bound-preserving and monotone solutions. However, the reference solution unsurprisingly violates undershoot bound (about 4% right after the saturation front) and produces an oscillatory saturation profile.
To better understand the efficacy of the proposed limiting algorithm (i.e., DG+SL+FL), and distinguish it from the vertex-based slope limiter of Kuzmin [2010] (i.e., DG+SL), we solve the problem again (with same parameters as before) on a crossed structured mesh (shown in figure  5(b)) for total duration of T = 450 s. The initial size of h = 10 m is chosen for this analysis and four-step refinement is performed. Table 3 reports the performance of limiters and compare them with respect to bound-preserving properties, local mass balance violations, and monotonocity. We observe that mesh refinement reduces maximum undershoots of unlimited DG from 47.61 % to 29.21 % and maximum overshoots to less than 0.1 % but does not eliminate violations. The application of slope limiter to DG eliminates undershoots at all time steps and significantly reduces maximum overshoots to 4.79 % for the coarsest mesh and to 3.32 % for the finest mesh. It can be seen that DG+SL falls short to satisfy maximum principle even under excessive mesh refinement. Further, it should be noted that both DG and DG+SL approximations fail to obtain monotone solutions near the saturation front. This means that they are susceptible to local spurious oscillations near the front even when the global bounds are not violated. However, approximations under DG+FL+SL enjoy pointwise maximum principle and the saturation field remains monotone over the entire domain, independently of the mesh size.
As shown in Table 4, only 3 to 4 Newton's iterations are needed at each time step for convergence of either limited DG or unlimited DG approximation. This means that the limiters do not have a significant effect on the number of solver iterations. However, as we refine the mesh, flux limiter algorithm L avg requires more iterations to converge. 4.2.2. Local mass balance. DG methods are known for their local mass conservation properties. [Riviere, 2008;Joshaghani et al., 2019]. In this section, we investigate the effect of the proposed limiters on altering local mass conservation properties. Upon applying element-wise averages and choosing unit test function in (3.4), we obtain the local mass conservation of an element E ∈ E h at (a) Saturation profile (b) Pressure profile Figure 6. Two-dimensional pressure-driven flow in homogeneous domain: This figure exhibits the saturation and pressure profiles obtained from limited DG approximations (with P = 1 and σ = 100). Solutions on uniform and non-uniform meshes are plotted along the line y = 50 m at t = 300 s and are compared with a reference DG solution. Regardless of the mesh size, limiters completely suppress unphysical overshoots and undershoots and accurately predict the location of saturation front which is in good agreement with that of the reference solution. As expected, the uniform finer mesh gives rise to a sharper front. On the other hand, the reference solution is not equipped with any bound-preserving mechanism and thus does not enjoy maximum principle and lower bound violations (S < 0.2) and nonmonotone behavior are captured. Table 3. This table shows the efficacy of the limiters when applied to the pressure-driven flow problem with homogeneous domain. Simulations are carried out for the duration of 450 s on crossed mesh (see Figure 5( (4.9)  We compute the magnitude of mass balance error for the pressure-driven flow problem discussed in the previous section. Table 3 contains the value of maximum error observed throughout the simulation. Evidently, DG+SL scheme is slightly worse than other two schemes with respect to errors, which is consistent for all mesh-sizes. However, values are all very small and below than the solver tolerance. In Figure 7, the values of M(E) are displayed at t = 450 s on a crossed mesh of size h = 2.5 m for three cases of DG, DG+SL, and DG+FL+SL. One can see that applying slope limiter (without flux limiter) instigates an erroneous patch (shown with dark brown color in Figure 7(b)). It is also clear that the proposed numerical scheme (i.e., DG+FL+SL) is locally mass conservative and slightly outperforms DG+SL scheme. 4.2.3. Domain with thin barrier. In this example, the porous medium contains a thin barrier and it is partitioned into an unstructured triangular mesh (see Figure 5(c)). Total time is set to T = 4500 s and the time step is t = 0.5 s. Additionally, noflow boundary conditions are imposed on the barrier edges. All other parameters are the same as in Section 4.2.1. Figure 8 exhibits the saturation profile under limited and unlimited DG at three different time steps. Limited DG, unlike its unlimited version, generates saturation that remains bounded and neither undershoots (blue-colored cells) nor overshoots (red-colored cells) are detected during the simulation. Nonetheless, the saturation front, under both unlimited and limited DG, propagates with the same speed and tends to avoids the barrier as expected. Figure 9 and 10 show the wetting phase pressure contour and velocity field at t = 4500 s. Velocities are computed at time t n , using the formula: u n w = −Kλ w (S n )∇P n . We can see that pressure drops linearly near the top and bottom edges, which confirms that fluid steers clear of the central barrier and flows around it. When no limiter is used, spurious oscillations and erroneous high-velocity regions are visible in velocity solutions. Limited DG, on the other hand, The color mapping for S in [0.2, 0.85] is grayscale, while values below and above bounds are colored blue and red, respectively. As expected, DG approximation with no limiter yields noticeable violations, while limited DG scheme is capable of providing maximumprinciple satisfying results. In spite of this, the front under both unlimited and limited DG, propagates with the same speed.  gives very smooth approximations. From these results we conclude that the proposed numerical scheme is bound-preserving on unstructured meshes. 4.2.4. Non-homogeneous domain. For this problem, permeability is 10 −8 m 2 everywhere except inside a square inclusion of length 20 m located at the center of the domain, where the permeability is 10 4 times smaller. The domain is discretized with a structured rectangular mesh of size h = 1.25 m. Time step is set to τ = 0.5 s and the simulation advances up to T = 650 s. The remaining parameters are the same as in Section 4.2.1. The discrete saturation at different snapshots of t = 350 and t = 650 s are depicted in Figure 11, where saturation values beyond the physical bounds (i.e.,   4.3. Quarter five-spot problem. In this section, the performance and robustness of the limiters are assessed in the presence of wells, for both homogeneous and heterogeneous permeabilities. We employ no flow boundary condition on the entire boundary, as shown in Figure 14; and assume zero capillary pressure. The flow is driven from an injection well at the bottom left corner to a production well at the top right corner. The wells are defined by source/sink terms, which are piecewise constant with compact support. That is,q is nonzero at injection well and q is nonzero at production well. The DG penalty parameters for test problems are taken as σ = 10. 4.3.1. Homogeneous domain. The domain Ω = [0, 100] 2 m 2 is partitioned into a crossed structured mesh of size h = 2.5 m, as depicted in Figure 5(b). The medium is homogeneous with K = 10 −13 m 2 everywhere. We choose Brooks-Corey relative permeabilities as follows: (4.10)  Figure 15 shows the wetting phase saturations at two different times (t = 10 and t = 21 days), for three schemes: DG, DG+SL, DG+FL+SL. The figure shows that violations of the maximum principle for the unlimited DG solution occur in the neighborhood of the injection well and after the front; in addition the DG solution is not monotone before the front. Adding a slope limiter helps with the monotonicity of the solution and with decreasing the number of elements where the maximum principle is not satisfied. The proposed numerical scheme, DG+FL+SL, completely eliminates violation of maximum principle: the solution is monotone and bound-preserving. Figure 16 and 17 show the wetting phase pressure contours and velocity fields at t = 10 days for all three cases. Differences are minimal for the pressure and velocity fields. Finally we display the local mass balance error in Figure 18 for all three cases; the local mass balance error is a piecewise constant field M defined by (4.9). We observe that the local mass balance is negligible (of the order of 10 −11 ) for the DG scheme with or without limiters. 4.3.2. Quarter five-spot problem with heterogeneous domain. We repeat the experiments in Section 4.3.1 with heterogeneous medium of Ω = [0, 1000] 2 m 2 . The permeability fields are discontinuous and values vary over seven orders of magnitude. The permeability data is taken from two layers of the SPE 10 data-set [Christie et al., 2001]; and are scaled to a crossed structured mesh of size h = 20 m (see permeability fields in log-scale in Figure 19). We note that layer 13 varies relatively smoothly, whereas layer 73 contains well-defined channels, which form an additional challenge for any numerical method. We set viscosities to µ w = 5 × 10 −4 Pa · s and µ = 2 × 10 −3 Pa · s and  (4.12) The production and injection wells of size L w = 100 m withq = q = 2.8 × 10 −5 are positioned at opposite corners such that d w = 70 m (see Figure 14). The time step is τ = 4.17 × 10 −3 days and the final time is T = 1.375 days.  We apply our proposed DG scheme with both flux and slope limiters to these porous media. Figure 20 displays the wetting phase saturation contours at different times (t = 0.417, 0.83, 1.375 days) for both layers. As expected the wetting phase floods the domain from the injection well to the production well while avoiding low permeable regions. Because of the location of channels in layer 73, the wetting phase has reached the production well at time t = 1.375 days whereas this is not the case for layer 13. We also observe that the saturation satisfies the maximum principle. Figure 21 shows the magnitude of the wetting phase velocity at the same times. The effect of the heterogeneities can be seen in the velocity fields.

Effect of gravity.
In this section, we examine the success of our limiting scheme in the presence of gravity field and then study the impact of gravity on the pressure-driven flows and quarter five-spot problems. The ratio of gravitational to viscous forces can be represented as a gravity number, Gr. This dimensionless parameter depends on the difference between phase densities; and following the work of Riaz and Tchelepi [2006]; Hassanizadeh and Das [2005]; , can be defined as follows:  Figure 19. Quarter five-spot problem with heterogeneous permeability: This figure illustrates the permeability fields adopted from two horizontal layers of SPE10 model 2 data-set. Layer 13 is taken from relatively smooth Tarbert formation, whereas layer 73 is taken from a highly varying Upper-Ness formation. Values are presented in logarithmic scale.
Pa · s. Here, the characteristic velocity is estimated to be U ≈ 0.1 m/s (using U ≈ −K(P | x=200 − P | x=0 )/µ w L). The wetting phase density is ρ w = 1000 kg/m 3 and the non-wetting phase density takes three different values ρ = 925, 850, 600 kg/m 3 , which yields three values for the gravity number Gr= 0.3, 0.6 and 1.6 respectively. Other parameters and Dirichlet boundary conditions are the same as in Section 4.2.1. The time step is τ = 0.6 s and the final time is T = 600 s. The proposed DG scheme with flux and slope limiters is applied and the penalty parameter is set to σ = 1000. Figure 22 shows the saturation contours at the time t = 600 s. As the gravity number increases, the wetting phase saturation, which is the heaviest, deposits more and more at the bottom of the domain; and the narrow gravity tongue along the bottom edge becomes more pronounced. It should be also noted that similar to earlier problems, the limiting scheme exhibits satisfactory results with respect the maximum principle. This means that for all three cases, solutions always remain between 0.2 and 0.85. Pressure contours and velocity fields are displayed in Figure 23 and 24. Both show the impact of gravity on the solutions. 4.4.2. Quarter five-spot problem. The domain is Ω = [0, 1000] 2 m 2 with permeability of K = 3 × 10 −11 everywhere. Capillary pressure and relative permeabilities are defined in equation (4.8) and (4.10), respectively, with entry pressure P d = 1000 Pa, θ = 2 and R = 0.05. To address wells, we fix the following parameters: L w = 80 m, d w = 80 m,q = q = 9.33 × 10 −6 . The wetting phase density is set to ρ w = 1000 kg/m 3 and the non-wetting phase density takes three different values ρ = 925, 850, 600 kg/m 3 , which yields three values for the gravity number Gr= 0.8, 1.6 and 4.3 respectively. The characteristic velocity in Gr estimation is taken as U ≈ 5.5 × 10 −5 m/s (or 4.8 m/day). Other parameters are the same as in Section 4.3.2. The simulation runs to T = 11 days with 750 time steps. Wetting phase saturation contours, wetting phase pressure contours and wetting phase velocity fields are shown in Figure 25, 26, and 27 respectively. We observe that as the gravity number increases, the inertial forces prevent the saturation to reach the production well. As in the previous section, the discrete solution satisfies the maximum principle. The numerical examples in this section confirm that our proposed numerical method is accurate and robust when gravity dominates.  Figure 20. Quarter five-spot problem with heterogeneous permeability: This figure shows the evolution of the saturation obtained using DG+FL+SL scheme for layer 13 (left) and layer 73 (right). For both cases, the wetting phase moves toward the production well by sweeping the regions with highest permeability values. Another inference is that proposed limiters yield physical values of saturation, without any overshoots and undershoots, even for domains with permeabilities that vary over several orders of magnitudes.

Conclusions
A fully implicit discontinuous Galerkin method is formulated for solving the incompressible two-phase flow equations in porous media. Primary unknowns are the wetting phase pressure and saturation. Nonlinear sytems are solved by Newton's method. Post-processing flux are developed and combined with slope limiters to ensure a bound-preserving saturation at each time step. The  . Quarter five-spot problem with gravity field: This figure depicts saturation contours at t = 11 days for different gravity numbers. DG+FL+SL scheme is applied that leads to satisfactory results with respect to maximum principle. By increasing the difference in phases density, gravitational force dominates the viscous force (from left to right). This results in more wetting phase saturation to be deposited at the bottom of domain and hence less non-wetting phase is recovered at the production well.
numerical method is validated on several benchmark problems and it is applied to problems where permeability fields are highly varying or where gravitational forces are significant. Flooding of the medium is driven by either pressure boundary conditions or by injection and production wells. The various numerical examples show that the scheme is robust and locally mass conservative. The approximation of the saturation is shown to satisfy the maximum principle both theoretically and computationally.
(a) Gr = 0.8 (b) Gr= 1.6 (c) Gr = 4.3 Figure 26. Quarter five-spot problem with gravity field: This figure shows pressure contours at t = 11 days for different gravity numbers. As gravity number increases (from left to right), pressure difference between the injection and the production wells reduces. This is because as the gravitational force dominates, it hinders the wetting phase flow from reaching the production well.
(a) Gr = 0.8 (b) Gr = 1.6 (c) Gr=4.3 Figure 27. Quarter five-spot problem with gravity field: This figure shows velocity field at t = 11 days for different gravity numbers. Increase in gravity number pushes the wetting phase toward the bottom edge and less recovery at the production well, which is also reflected in the decrease in the magnitude of velocities.