A convergent mass conservative numerical scheme based on mixed finite elements for two-phase flow in porous media

# A convergent mass conservative numerical scheme based on mixed finite elements for two-phase flow in porous media

Florin A. Radu Kundan Kumar Jan M. Nordbotten Iuliu S. Pop
Department of Mathematics, University of Bergen, P. O. Box 7800, N-5020 Bergen, Norway
Department of Mathematics and Computer Science, Eindhoven University of Technology,
P. O. Box 513, 5600 MB Eindhoven, The Netherlands
###### Abstract

Abstract. In this work we present a mass conservative numerical scheme for two-phase flow in porous media. The model for flow consists on two fully coupled, non-linear equations: a degenerate parabolic equation and an elliptic equation. The proposed numerical scheme is based on backward Euler for the temporal discretization and mixed finite element method (MFEM) for the discretization in space. Continuous, semi-discrete (continuous in space) and fully discrete variational formulations are set up and the existence and uniqueness of solutions is discussed. Error estimates are presented to prove the convergence of the scheme. The non-linear systems within each time step are solved by a robust linearization method. This iterative method does not involve any regularization step. The convergence of the linearization scheme is rigorously proved under the assumption of a Lipschitz continuous saturation. Numerical results are presented to sustain the theoretical findings.

Keywords: linearization, two-phase flow, mixed finite element method, convergence analysis, a priori error estimates, porous media, Richards’ equation, degenerate parabolic problems, coupled problems.

## 1 Introduction

Two-phase porous media flow models are widely encountered in real-life applications of utmost societal relevance, including water and soil pollution, oil recovery, geological carbon dioxide sequestration, or nuclear waste management [34, 21]. Such complex problems admit only in very simplified situations analytical solutions, therefore numerical methods for solving multiphase flow in porous media are playing a determining role in understanding and solving the problems. Nevertheless, the design and analysis of robust, accurate and efficient numerical schemes is a very challenging task.

Here we discuss a numerical scheme for a two-phase porous media flow model. The fluids are assumed immiscible and incompressible and the solid matrix is non-deformable. The adopted formulation uses the global pressure and a complementary pressure, obtained by using the Kirchhoff transformation, as primary unknowns (see [11, 3, 12]). This leads to a system of two coupled non-linear partial differential equations, a degenerate elliptic - parabolic one and an elliptic one.

Numerical methods for two-phase flow have been the object of intensive research in the last decades. The major challenge in developing efficient schemes is related to the degenerate nature of the problem. Due to this, the solution typically lacks regularity, which makes lower order finite elements or finite volumes a natural choice for the spatial discretization. In this respect, we refer to [20, 35] for Galerkin finite elements, to [19, 37, 33] for finite volumes, to [16, 13, 14] for methods combining Galerkin finite elements combined with the mixed finite element method (MFEM), and to [17, 52, 22] for the discontinuous Galerkin method. In all cases, the convergence of the numerical schemes is proved rigorously either by compactness arguments, or by obtaining a priori error estimates. A posteriori error estimates are obtained e.g. in . Furthermore, similar issues appear for the Richards equation, which is a simplified model for saturated/unsaturated flow in the case when the pressure of one phase is supposed to be constant. In this context we mention Galerkin finite elements [35, 39], MFEM based works [4, 41, 42, 45, 49, 57, 59], multipointflux approximation (MPFA)  or finite volume - MFEM combined methods . We also refer to [44, 48, 26, 27] for related works concerning reactive transport in porous media.

In this paper we propose a mass conservative scheme based on MFEM (lowest order Raviart-Thomas elements ) and backward Euler for numerical simulation of the two-phase flow in porous media. Continuous, semi-discrete (continuous in space) and fully discrete mixed variational formulations are defined. Existence and uniqueness of solutions is discussed, the equivalence with a conformal formulation being involved in the proof. We show the convergence of the numerical scheme and provide explicit order of convergence estimates. The analysis is inspired by similar results in [4, 14, 41, 45].

Typical problems involving flow in porous media, like e.g. water and soil pollution or nuclear waste management are spread over decades or even centuries, so that the use of relatively large time steps is a necessity. Due to this, implicit methods are a necessity (our choice here being the first-order backward Euler method, due to the low regularity of the considered problem). Since the original model is non-linear, at each time step one needs to solve non-linear algebraic systems. In this work we propose a robust linearization scheme for the systems appearing at each time step, as a valuable alternative to modified Picard method  or Newton’s method [5, 38, 43, 36, 30, 46, 47] or iterative IMPES [23, 24]. Although the applicability of Newton’s method for parabolic equations is well recognized, its convergence is not straightforward for degenerate equations, where the Jacobian might become singular. A possible way to overcome this is to regularize the problem. However, even in this case convergence is guaranteed only under a severe stability condition for the discretization parameters, see . This has motivated the alternative, robust linearization scheme proposed in this work. The new scheme, called scheme from now on, does not involve the calculations of any derivatives and does not need a regularization step. The -scheme combines the idea of a classical Picard method and the scheme presented in  for MFEM or [58, 54] for Galerkin finite elements. The -scheme was proposed for two-phase flow in combination with the MPFA method in , the proof of convergence there being only sketched and not made completely rigorous. We show here that the -scheme for MFEM based discretizations converges linearly if the time step satisfies a mild condition. This robustness is the main advantage of the scheme when compared to the quadratic, but locally convergent Newton method.

Finally, we mention that the -scheme can be interpreted as a non-linear preconditioner, because the linear systems to be solved within each iteration are much better conditioned than the corresponding systems in the case of modified Picard or Newton’s method. We refer to  for illustrative examples concerning the Richards equation, which is a particular case of the more general model considered in the present work.

To summarize, the main new contributions of this paper are

• We present and analyze a MFEM based numerical scheme for two-phase flow in porous media. Order of convergence estimates are obtained.

• We show the existence and uniqueness of the considered variational formulation. This is based on the equivalence between the conformal and the mixed formulations, which is proved here for the continuous and the time discrete models.

• We present and analyze rigorously a robust, first-order convergent linearization method for MFEM based schemes for two-phase flow in porous media.

The paper is structured as follows. In Section 2, we present the model equations for two-phase flow in porous media and we define the discretization scheme. In Section 3 we analyze the convergence of the discretization scheme based on a priori error estimates. We also prove existence and uniqueness for the problem involved, and give stability estimates. A new MFEM linearization scheme is presented and analyzed in Section 4. Section 5 provides numerical examples confirming the theoretical results. The paper is ending with concluding remarks in Section 6.

## 2 Mathematical model and discretization

In this section we introduce the notations used in this work, the mathematical model and its MFEM/Euler implicit discretization (Problems , and ). A linearization scheme (Problem ) is proposed to solve the non-linear systems appearing at each time step.

The model is defined in the -dimensional bounded domain having a Lipschitz continuous boundary . Further is the final computational time. We use common notations from functional analysis, e.g. is the space of essential bounded functions on , the space of square integrable functions on , or is the subspace of containing functions which have also first order derivatives in . We denote by the space of functions with a vanishing trace on and by its dual. denotes the inner product in , or the duality pairing between and . Further, , and stand for the norms in , , respectively . The functions in are vector valued having a divergence. The norm in is denoted by . denotes the Bochner space of -valued functions defined on , where is a Banach space. Similarly, are -valued functions continuous (w. r. t. norm) on . By we mean a generic positive constant, not depending on the unknowns or the discretization parameters and we denote by the Lipschitz constant of a (Lipschitz continuous) function .

Further, we will denote by an integer giving the time step . For a given , the n time point is . We will also use the following notation for the mean over a time interval. Given the function ( being a Banach space like , or )), its time averaged over the interval is defined as

 ¯gn:=1τ∫tntn−1g(t)dt.

Clearly, this is an element in as well.

The two-phase porous media flow model considered here assumes that the fluids are immiscible and incompressible, and that the solid matrix is non-deformable. By denoting with the wetting and non-wetting phases, the saturation, pressure, flux and density of phase , respectively, the two-phase model under consideration reads (see e.g. [6, 11, 21, 34])

 ∂(ϕραsα)∂t+∇⋅(ραqα) = 0,α=w,n, (1) qα = −kr,αμαk(∇pα−ραg),α=w,n, (2) sw+sn = 1, (3) pn−pw = pcap(sw), (4)

where denotes the constant gravitational vector. Equation (1) is a mass balance, (2) is the Darcy law, (3) is an algebraic evidence expressing that all pores in the medium are filled by a mixture of the two fluid phases and (4) is the capillary pressure relationship, with supposed to be known. The porosity , permeability , the viscosities are given constants and the relative permeabilities are given functions of . We consider here a scalar permeability, but the results can be easily extended to the case when the permeability is positive-definite tensor.

In this paper we adopt a global/complementary pressure formulation [3, 11, 12]. The global pressure (denoted by ) was introduced in  and the complementary pressure in . They are defined by

 p(x,sw):=pn(x)−∫sw0fw(x,ξ)∂pcap∂ξ(x,ξ)dξ, (5) Θ(x,sw):=−∫sw0fw(x,ξ)λn(x,ξ)∂pcap∂ξ(x,ξ)dξ, (6)

where we denoted by , the phase mobilities and by the fractional flow function. We note the use of Kirchhoff transformation above. In the new unknowns, the resulting system consists of two coupled non-linear partial differential equations, a degenerate parabolic one and an elliptic one. For more details on the modelling we refer to , where the existence and uniqueness of a weak solution is proved for a Galerkin-MFEM formulation. In the new unknowns the system (1)-(4) becomes

 ∂ts(Θ)+∇⋅q = 0, (7) q = −∇Θ+fw(s)u+f1(s), (8) ∇⋅u = f2(s), (9) a(s)u = −∇p−f3(s). (10)

with , , the (wetting) flux, and the total flux. The equations hold true in . The coefficient functions are given and satisfy the assumptions listed below. The system is completed by initial conditions specified below, and by boundary conditions. For simplicity, we restrict our attention to homogeneous Dirichlet boundary conditions, but other kinds of conditions can be considered.

Also note that the results can be extended straightforwardly to the case when a source term is present on the right hand side of (7). For the ease of presentation and in view of the analogy with the model considered in , the source terms are left out here. This simplifies the presentation of the convergence proof in Section 3.

Problem : Continuous mixed variational formulation.
Find , , such that there holds , , and

 ⟨s(Θ(t))−s(Θ0),w⟩+⟨∇⋅∫t0q(y)dy,w⟩ = 0, (11) ⟨∫t0q(y)dy,v⟩−⟨∫t0Θ(y)dy,∇⋅v⟩ −⟨∫t0fw(s(Θ(y)))u(y)dy,v⟩ = ⟨∫t0f1(s(Θ(y)))dy,v⟩, (12) ⟨∇⋅u(t),w⟩ = ⟨f2(s(Θ(t))),w⟩, (13) ⟨a(s(Θ(t)))u(t),v⟩−⟨p(t),∇⋅v⟩+⟨f3(s(Θ(t))),v⟩ = 0 (14)

for all , and , with .

The function is given. For example, if the function is one to one, a natural initial condition is , where is a given initial saturation. By (11), since is continuous in time, it follows that . Then, since and are assumed continuous (see (A1) and (A3) below), from (13) one obtains that . Similarly, is continuous in time as well, so (12)-(14) hold for all .

We now proceed with the time discretization for Problem , which is achieved by the Euler implicit scheme. For a given , we define the time discrete mixed variational problem at time (the time step is denoted by ):

Problem : Semi-discrete variational formulation. Let be given. Find and such that

 ⟨sn−sn−1,w⟩+τ⟨∇⋅qn,w⟩ = 0, (15) ⟨qn,v⟩−⟨Θn,∇⋅v⟩−⟨fw(sn)un,v⟩ = ⟨f1(sn),v⟩, (16) ⟨∇⋅un,w⟩ = ⟨f2(sn),w⟩, (17) ⟨a(sn)un,v⟩−⟨pn,∇⋅v⟩+⟨f3(sn),v⟩ = 0 (18)

for all , and . Initially we take . Throughout this paper stands for , , making the presentation easier.

We can now proceed with the spatial discretization. For this let be a regular decomposition of into closed -simplices; stands for the mesh-size (see ). Here we assume , hence is polygonal. Thus we neglect the errors caused by an approximation of a non-polygonal domain and avoid an excess of technicalities (a complete analysis in this sense can be found in ).

The discrete subspaces are defined as

 Wh:={p∈L2(Ω)| p is constant on each % element T∈Th},Vh:={q∈H(div;Ω)|q|T(x)=aT+bTx,aT∈R2,bT∈R for all T∈Th}. (19)

So denotes the space of piecewise constant functions, while is the space (see ).

We will use the following projectors (see  and , p. 237):

 Ph:L2(Ω)→Wh,⟨Phw−w,wh⟩=0, (20)

and

 Πh:H(div;Ω)→Vh,⟨∇⋅(Πhv−v),wh⟩=0, (21)

for all , and . For these operators we have

 ∥w−Phw∥≤Ch∥w∥1, respectively ∥v−Πhv∥≤Ch∥v∥1 (22)

for any and .

The fully discrete (non-linear) scheme can now be given. To simplify notations we use in the following the notation , .

Problem : Fully discrete (non-linear) variational formulation. Let , and assume is known. Find and such that there holds

 ⟨snh−sn−1h,wh⟩+τ⟨∇⋅qnh,wh⟩ = 0, (23) ⟨qnh,vh⟩−⟨Θnh,∇⋅vh⟩−⟨fw(snh)unh,vh⟩ = ⟨f1(snh),vh⟩, (24) ⟨∇⋅unh,wh⟩ = ⟨f2(snh),wh⟩, (25) ⟨a(snh)unh,vh⟩−⟨pnh,∇⋅vh⟩+⟨f3(snh),vh⟩ = 0 (26)

for all and all .

The fully discrete scheme (23) – (26) is non-linear, and iterative schemes are required for solving it. Moreover, in (A1) it is only required that is monotone and Lipschitz, so it may have a vanishing derivative, which makes (7) degenerate. This is, indeed, the situation encountered in two-phase porous media flow models. In such cases, usual schemes such as the Newton method may not converge without performing a regularization step. This may affect the mass balance. For the Richards equation, this is proved in . Following the ideas in [40, 54, 58], we propose a robust, first order convergent linearization scheme for solving (23) – (26). The scheme is not requiring any regularization. Similar ideas can be applied in connection with any other spatial discretization method, see e.g.  for multipoint flux approximation. The analysis of the scheme is presented in Section 4.

Let be fixed. Assuming that is Lipschitz continuous, the following iterative scheme can be used to solve the non-linear problem (23) – (26):

Problem : Linearization scheme (-scheme). Let , and let be given. Find and such that

 ⟨L(Θn,ih−Θn,i−1h)+sn,i−1h,wh⟩+τ⟨∇⋅qn,ih,wh⟩ = ⟨sn−1h,wh⟩, (27) ⟨qn,ih,vh⟩−⟨Θn,ih,∇⋅vh⟩−⟨fw(sn,i−1h)un,ih,vh⟩ = ⟨f1(sn,i−1h),vh⟩, (28) ⟨∇⋅un,ih,wh⟩ = ⟨f2(sn,i−1h),wh⟩, (29) = 0 (30)

for all and all . We use the notation , and, as previously, , . For starting the iterations, a natural choice is and, correspondingly, . Note that since we prove below that the iterative scheme is a contraction, this choice is not compulsory for the convergence.

Throughout this paper we make the following assumptions:

• The function is monotone increasing and Lipschitz continuous.

• is Lipschitz continuous and there exists such that for all one has

 0
• and are Lipschitz continuous. Additionally, is uniformly bounded.

• There exits a constant such that , , and for all , the last two being uniformly in and . Here , and are the solution components in Problems P, P and P respectively.

• The function is in .

###### Remark 2.1.

The assumptions are satisfied in most situations of practical interest. Concerning (A4), for this is practically the outcome of the assumptions (A1) and (A3), which guarantee that for every one has , and that the norm is bounded uniformly w.r.t. time. Now, without being rigorous, we observe that by (13) one obtains , where satisfies . Classical regularity theory (see e.g. , Thm. 15.1 in Chapter 3) guarantees that is continuous on the compact , and that the norm can be bounded uniformly in time. For the approximation , one can reason in the same manner, and observe that becomes the projection . Since is bounded uniformly in time, the construction of the projector (see e. g. , Chapter 7.2) guarantees that satisfies the same bounds as . Finally, case of is similar. We also refer to  for a similar situation but in the case of a one phase flow model, where conditions ensuring the validity of (A4) are provided.

###### Remark 2.2.

One may relax the Lipschitz continuity for the parameter functions to a more general context, where e. g. satisfies a growth condition . We do not exploit this possibility here, but refer to [12, 45, 41] for the details on the procedure.

The following two technical lemmas will be used in Sections 3 and 4. Their proofs can be found e.g. in  and , respectively.

###### Lemma 2.1.

Given a , there exists a such that

 ∇⋅v=wand∥v∥≤CΩ∥w∥,

with not depending on .

###### Lemma 2.2.

Given a , there exits a satisfying

 ∇⋅vh=whand∥vh∥≤CΩ,d∥wh∥,

with not depending on or mesh size.

Also, the following elementary result will be used

###### Proposition 2.1.

Let be a set of vectors. It holds

 N∑n=1⟨an,n∑k=1ak⟩ = 12∥∥ ∥∥N∑n=1an∥∥ ∥∥2+12N∑n=1∥an∥2. (32)

## 3 Analysis of the discretization: existence and uniqueness and a priori error estimates

In this section we analyze the problems , and introduced in Section 2. The existence and uniqueness of a solution will be discussed in Subsection 3.1. For the continuous and semi-discrete cases this will be done by showing an equivalence with conformal variational formulations. The convergence of the numerical scheme will be shown by deriving a priori error estimates. The main convergence result is given in Theorem 3.1. The convergence is established by assuming that the non-linear systems (23) - (26) are solved exactly. We refer to Section 4 for the analysis of the linearization scheme (27)-(30), which was proposed in the previous section to solve these non-linear algebraic systems numerically.

### 3.1 Existence and uniqueness for the variational problems

In this subsection we discuss the existence and uniqueness of the continuous, semi-discrete, and fully discrete variational formulations for the considered model (7)-(10). We establish an equivalence between the continuous mixed formulation and a conformal formulation, which will deliver the existence and uniqueness for the continuous case. The semi-discrete case can be treated analogously. For the fully discrete case we prove below the uniqueness. Existence can be proved by fixed point arguments, using e.g. Lemma 1.4, p. 164 in . We omit the details here as the existence is also a direct consequence of the results in Section 4, where the convergence of a linear iterative scheme is proved. The limit of this iteration is exactly a solution for the fully discrete system.

A conformal variational formulation for the model (7)-(10) reads:

Problem : Continuous conformal variational formulation. Let be given. Find , , such that , , and for all and one has

 ∫T0⟨∂ts(ΘC),v⟩dt+∫T0⟨∇ΘC+fw(sΘC)a(s(ΘC))∇pC,∇v⟩dt +∫T0⟨fw(sΘC)a(s(ΘC))f3(s(ΘC))−f1(s(ΘC)),∇v⟩dt = 0, (33) ⟨1a(s(ΘC))(∇pC+f3(s(ΘC))),∇w⟩ = −⟨f2(s),w⟩. (34)

The existence and uniqueness of a solution for Problem has been studied intensively in the past. Closest to the framework considered here is  (see also ). There the existence and uniqueness is proved, however, for the case that the inverse of is Lipschitz (the so-called slow diffusion case). Here we assume Lipschitz but not necessarily strictly increasing, which is a fast diffusion case. Other relevant references for the existence and uniqueness are [1, 2, 20, 25]. Also, we refer to  for the existence of a solution in heterogeneous media, where the phase pressure differences may become discontinuous at the interface separating two homogeneous blocks.

Having this in mind, one can use the existence and uniqueness results for the conformal formulation to obtain the existence of a solution for Problem (as by-product also establish its regularity). The equivalence is established in Proposition 3.1, whose proof follows the ideas in , Proposition 2.2.

###### Proposition 3.1.

Let be a solution to Problem , define and assume that (A1)-(A5) hold true. Then, a solution to Problem is given by , , and . Conversely, if , are solving Problem , then and is a solution of Problem PC.

Proof. Clearly, and defined above have the regularity required in Problem . Furthermore, and are elements of . Recalling that is Lipschitz continuous, one immediately obtains that . Since this shows that . With and arbitrary chosen, taking now in (33) and using the definition of gives

 ⟨s(ΘC)−s(Θ0),ϕ⟩−⟨∫t0q(y)dy,∇ϕ⟩=0, (35)

for all . In other words, in distributional sense. The regularity of mentioned above implies that, actually, lies in as well, and that . Moreover, by density arguments (11) holds for any . Also, from (35) one gets .

In a similar way, using (34) and the definition of one obtains

 ⟨−u,∇w⟩=⟨f2(sC),w⟩,

for all , so for a. e. . In view of the regularity of and of the assumptions on , this means that and that (13) holds for every .

To obtain (14) one uses (34), the definition of and that has a vanishing trace on . This gives

for a. e. and for all . Recalling the continuity in time for and , it follows that (14) is valid for every . Finally, one can use similar ideas to show that (12) holds true as well.
Q. E. D.

Let now be the solution of Problem . We have to show that is the solution of , i.e. to show that the functions are actually in and that they satisfy (33)-(34). Further, since in (A1) is assumed Lipschitz, it follows that as well. Clearly, (11) gives for a.e. . Since one gets .

Taking arbitrary as test function in (14) we get

 ⟨a(s(Θ(y)))u,v⟩+⟨∇p,v⟩+⟨f3(s(Θ(y))),v⟩=0, (36)

which implies

 ∇p=−a(s(Θ(y)))u−f3(s(Θ(y))) (37)

first in a distributional sense, and then by using the regularity of and in sense. It follows that . It remains to verify that has a vanishing trace on the boundary of . Taking now and using the regularity of and (14) one obtains

 ⟨∇p,v⟩=−⟨a(s(Θ(y)))u,v⟩−⟨f3(s(Θ(y))),v⟩=−⟨p,∇⋅v⟩, (38)

for all . By using now the Green theorem for , see , pg. 91

 ∫Γpv⋅nds=(∇p,v)+(p,∇⋅v)=0.

It follows immediately that . From (37) and (13) one gets that satisfies (34). In a similar manner one can show that and it satisfies (33).
Q. E. D.

We can now state an analogous result for the semi-discrete case.

Problem : Semi-discrete conformal variational formulation. Let and be given. Find , , such that for all one has

 ⟨s(ΘnC)−s(Θn−1C),v⟩+⟨∇ΘnC+fw(sΘnC)a(s(ΘnC))∇pnC,∇v⟩ +⟨fw(sΘnC)a(s(ΘnC))f3(s(ΘnC))−f1(s(ΘnC)),∇v⟩ = 0, (39) ⟨1a(s(ΘnC))(∇pnC+f3(s(ΘnC))),∇w⟩ = −⟨f2(s(ΘnC)),w⟩. (40)
###### Proposition 3.2.

Let , be given and assume that (A1)-(A5) hold true. If is a solution to Problem then, a solution to Problem is given by , , and . Conversely, if , are solving Problem , then and is a solution of Problem .

The proof of Proposition 3.2 is similar with the one of Proposition 3.1 (also see , Prop. 2.3) and it will be skipped here. Then the existence and uniqueness of a solution for the semi-discrete variational formulation (15)–(18) is a direct consequence of the equivalence result.

The following proposition establish the uniqueness for the fully-discrete variational problem (23)-(26). As explained in the introduction of Sec. 3.1, the existence of a solution will follow from the convergence of the linear iteration scheme.

###### Proposition 3.3.

Let be fixed. If (A1)-(A5) hold true and the time step is sufficiently small, then the problem (23)-(26) has at most one solution.

Proof. Let us assume that there exists two solutions of (23)-(26), with . Further, we let stand for the two saturations. These solutions are then satisfying

 ⟨snh,i−sn−1h,wh⟩+τ⟨∇⋅qnh,i,wh⟩ = 0, (41) ⟨qnh,i,vh⟩−⟨Θnh,i,∇⋅vh⟩−⟨fw(snh,i)unh,i,vh⟩ = ⟨f1(snh,i),vh⟩, (42) ⟨∇⋅unh,i,wh⟩ = ⟨f