### Abstract

A new fully discrete stabilized discontinuous Galerkin method is proposed to solve the incompressible miscible displacement problem. For the pressure equation, we develop a mixed, stabilized, discontinuous Galerkin formulation. We can obtain the optimal priori estimates for both concentration and pressure.

##### Keywords:

Discontinuous Galerkin methods; a priori error estimates; incompressible miscible displacement### 1 Introduction

We consider the problem of miscible displacement which has considerable and practical importance in petroleum engineering. This problem can be considered as the result of advective-diffusive equation for concentrations and the Darcy flow equation. The more popular approach in application so far has been based on the mixed formulation. In a previous work, Douglas and Roberts [1] presented a mixed finite element (MFE) method for the compressible miscible displacement problem. For the Darcy flow, Masud and Hughes [2] introduced a stabilized finite element formulation in which an appropriately weighted residual of the Darcy law is added to the standard mixed formulation. Recently, discontinuous Galerkin for miscible displacement has been investigated by numerical experiments and was reported to exhibit good numerical performance [3,4]. In Hughes-Masud-Wan [5], the method of [2] was extended to the discontinuous Galerkin framework for the Darcy flow. A family of mixed finite element discretizations of the Darcy flow equations using totally discontinuous elements was introduced in [6]. In [7] primal semi-discrete discontinuous Galerkin methods with interior penalty are proposed to solve the coupled system of flow and reactive transport in porous media, which arises from many applications including miscible displacement and acid-stimulated flow. In [8], stable Crank-Nicolson discretization was given for incompressible miscible displacement problem.

The discontinuous Galerkin (DG) method was introduced by Reed and Hill [9], and extended by Cockburn and Shu [10-12] to conservation law and system of conservation laws,respectively. Due to localizability of the discontinuous Galerkin method, it is easy to construct higher order element to obtain higher order accuracy and to derive highly parallel algorithms. Because of these advantages, the discontinuous Galerkin method has become a very active area of research [4-7,13-18]. Most of the literature concerning discontinuous Galerkin methods can be found in [13].

In this paper, we analyze a fully discrete finite element method with the stabilized
mixed discontinuous Galerkin methods for the incompressible miscible displacement
problem in porous media. For the pressure equation, we develop a mixed, stabilized,
discontinuous Galerkin formulation. To some extent, we develop a more general stabilized
formulation and because of the proper choose of the parameters *γ *and *β*, this paper includes the methods of [2,6] and [5]. All the schemes are stable for any combination of discontinuous discrete concentration,
velocity and pressure spaces. Based on our results, we can assert that the mixed stabilized
discontinuous Galerkin formulation of the incompressible miscible displacement problem
is mathematically viable, and we also believe it may be practically useful. It generalizes
and encompasses all the successful elements described in [2,6] and [5]. Optimal error estimate are obtained for the concentration, velocity and pressure.

An outline of the remainder of the paper follows: In Section 2, we describe the modeling equations. The DG schemes for the concentration and some of their properties are introduced in Section 3. Stabilized mixed DG methods are introduced for the velocity and pressure in Section 4. In Section 5, we propose the numerical approximation scheme of incompressible miscible displacement problems with a fully discrete in time, combined with a mixed, stabilized and discontinuous Galerkin method. The boundedness and stability of the finite element formulation are studied in Section 6. Error estimates for the incompressible miscible displacement problem are obtained in Section 7.

Throughout the paper, we denote by *C *a generic positive constant that is independent of *h *and Δ*t*, but might depend on the partial differential equation solution; we denote by *ε *a fixed positive constant that can be chosen arbitrarily small.

### 2 Governing equations

Miscible displacement of one incompressible fluid by another in a porous medium Ω
∈ *R ^{d}*(

*d*= 2, 3) over time interval

*J*= (0,

*T*] is modeled by the system concentration equation:

Pressure equation:

The initial conditions

The no-flow boundary conditions

Dispersion/diffusion tensor

where the unknowns are *p *(the pressure in the fluid mixture), **u **(the Darcy velocity of the mixture, i.e., the volume of fluid flowing cross a unit
across-section per unit time) and *c *(the concentration of the interested species, i.e., the amount of the species per
unit volume of the fluid mixture). *ϕ *= *ϕ*(*x*) is the porosity of the medium, uniformly bounded above and below by positive numbers.
The **E**(**u**) is the tensor that projects onto the **u **direction, whose (*i,j*) component is
*d _{m }*is the molecular diffusivity and assumed to be strictly positive;

*d*and

_{l }*d*are the longitudinal and the transverse dispersivities, respectively, and are assumed to be nonnegative. The imposed external total flow rate

_{t }*q*is sum of sources (injection) and sinks (extraction) and is assumed to be bounded. Concentration

*c** in the source term is the injected concentration

*c*if

_{w }*q*≥ 0 and is the resident concentration

*c*if

*q*< 0. Here, we assume that the

*a*(

*c*) is a globally Lipschitz continuous function of

*c*, and is uniformly symmetric positive definite and bounded.

### 3 Discontinuous Galerkin method for the concentration

#### 3.1 Notation

Let *T _{h }*= (

*K*) be a sequence of finite element partitions of Ω. Let Γ

_{I }denote the set of all interior edges, Γ

_{B }the set of the edges

*e*on

*∂*Ω, and Γ

_{h }= Γ

_{B }+ Γ

_{I}.

*K*

^{+},

*K*

^{- }be two adjacent elements of

*T*; let

_{h}*x*be an arbitrary point of the set

*e*=

*∂K*

^{+ }∩

*∂K*

^{-}, which is assumed to have a nonzero (

*d*- 1) dimensional measure; and let

*n*

^{+},

*n*

^{- }be the corresponding outward unit normals at that point. Let (

**u**,

*p*) be a function smooth inside each element

*K*

^{± }and let us denote by (

**u**

^{±},

*p*

^{±}) the traces of (

**u**,

*p*) on

*e*from the interior of

*K*. Then we define the mean values {{·}} and jumps [[·]] at

^{±}*x*∈ {

*e*} as

For *e *∈ Γ_{B}, the obvious definitions is {{*p*}} = *p*, [[**u**]] = **u**·**n**, with **n **denoting the outward unit normal vector on *∂*Ω. we define the set 〈*K, K'*〉 as

For *s *≥ 0, we define

The usual Sobolev norm on Ω is denoted by ||·||_{m, Ω }[19]. The broken norms are defined, for a positive number *m*, as

The discontinuous finite element space is taken to be

where *P _{r}*(

*K*) denotes the space of polynomials of (total) degree less than or equal to

*r*(

*r*≥ 0) on

*K*. Note that we present error estimators in this paper for the local space

*P*, but the results also apply to the local space

_{r}*Q*(the tensor product of the polynomial spaces of degree less than or equal to

_{r }*r*in each spatial dimension) because

*P*(

_{r}*K*) ⊂

*Q*(

_{r}*K*).

The cut-off operator

where *M *is a large positive constant. By a straightforward argument, we can show that the
cut-off operator

**Lemma 3.1 **[7] (Property of operator
*The cut-off operator *
*defined as in Equation 3.4 is uniformly Lipschitz continuous with a Lipschitz constant
one, that is*

*We shall also use the following inverse inequalities, which can be derived using the
method in *[20]. *Let K *∈ *T _{h}*,

*v*∈

*P*(

_{r}*K*)

*and h*

_{K }is the diameter of K. Then there exists a constant C independent of v and h_{K}, such that

#### 3.2 Discontinuous Galerkin schemes

Let ∇_{h }· *v *and ∇* _{h}v *be the functions whose restriction to each element

*v*, ∇

*v*, respectively. We introduce the bilinear form

*B*(

*c, w*;

**u**) and the linear functional

*L*(

*w*;

**u**,

*c*)

with

here *c*_{11 }> 0 is a constant independent of the meshsize.

We now define the weak formulation on which our mixed discontinuous method is based

Let *N *be a positive integer,
*t _{m }*=

*m*Δ

*t*for

*m*= 0, 1, ...,

*N*. The approximation of

*c*at

_{t }*t*=

*t*

^{n+1 }can be discreted by the forward difference. The DG schemes for approximating concentration are as follows. We seek

*c*∈

_{h }*W*

^{1,∞}(0,

*T*;

*D*

_{k-1}(

*T*)) satisfying

_{h}

where
**u**_{h }defined below

### 4 A stabilized mixed DG method for the velocity and pressure

#### 4.1 Elimination for the flux variable *u*

Letting *α*(*c*) = *a*(*c*)^{-1}. For the velocity and pressure, we define the following forms

The discrete problem for the velocity and pressure can be written as: find **u**_{h }∈ (*D*_{l-2}(*T _{h}*))

^{d}, (

*l*≥ 2),

*p*∈

_{h }*D*

_{l-1}(

*T*) such as

_{h}

In order to eliminate the flux variable, we first recall a useful identity, that holds
for vectors **u **and scalars *ψ *piecewise smooth on *T _{h}*:

Using (4.4) we have

Substituting (4.5) in the first equation of (4.3) we obtain

We introduce the lift operator *R*:*L*^{1}(∪*∂K*) → (*D*_{l-2}(*T _{h}*))

^{d }defined by

From (4.6) and (4.7) we have

We also introduce the *L*^{2}-projection *π *onto (*D*_{l-2}(*T _{h}*))

^{d}

Equation 4.8 gives now

Noting that ∇_{h}*D*_{l-1}(*T _{h}*) ⊂ (

*D*

_{l-2}(

*T*))

_{h}^{d}, we have

*π*∇

*≡ ∇*

_{h}p_{h }*for all*

_{h}p_{h }*p*∈

_{h }*D*

_{l-1}(

*T*). The Equation 4.10 gives

_{h}

Using (4.5) and the lifting operator *R *defined in (4.7) we have

Substituting (4.12) in the second equation of (4.3) and using (4.11) we have

For future reference, it is convenient to rewrite (4.13) as follows

where

#### 4.2 Stabilization of formulation (4.3)

We write first (4.3) in the equivalent form: find (**u**_{h}, *p _{h}*) ∈ (

*D*

_{l-2}(

*T*))

_{h}^{d }×

*D*

_{l-1}(

*T*) such that

_{h}

where

In a sense, (4.16) can be seen as a Darcy problem. The usual way to stabilized it
is to introduce penalty terms on the jumps of *p *and/or on the jumps of *u*. In [2], Masud and Hughes introduced a stabilized finite element formulation in which an
appropriately weighted residual of the Darcy law is added to the standard mixed formulation.
In Hughes-Masud-Wan [5], the method was extend within the discontinuous Galerkin framework. A family of mixed
finite element discretizations of the Darcy flow equations using totally discontinuous
elements was introduced in [6]. In this paper, we consider the following stabilized formulation which includes the
methods of [2,6] and [5].

The stabilized formulation of (4.16) is

where

where *γ *and *β *are chosen as the following (i) *γ *= 1, *β *= 1. (ii) *γ *= 0, *β *= 1, *δ *could assume either the value +1 or the value -1. The definition of *θ *will be given in the following content.

### 5 A mixed stabilized DG method for the incompressible miscible displacement problem

By combining (3.8) with (4.18), we have the stabilized DG for the approximating (2.1)-(2.5):
seek *c _{h }*∈

*W*

^{1,∞}(0,

*T*;

*D*

_{k-1}(

*T*)) =:

_{h}*W*,

_{h}*p*∈

_{h }*W*

^{1,∞}(0,

*T*;

*D*

_{l-1}(

*T*)) =:

_{h}*Q*and

_{h }**u**

_{h }∈ (

*W*

^{1,∞}(0,

*T*;

*D*

_{l-2}(

*T*)))

_{h}^{d }=:

*V*satisfying

_{h }

We define the "stability norm" by

where

### 6 Stability and consistency

From [6], we can state the following results.

**Lemma 6.1 **[6]*There exist two positive constants C*_{1 }*and C*_{2}, *depending only on the minimum angle of the decomposition and on the polynomial degree*

**Lemma 6.2 **[6]*There exists two positive constants C*_{1 }*and C*_{2}, *depending only on the minimum angle of the decomposition such that*

**Lemma 6.3 **[6]*Let *
*be a Hilbert spaces, and λ and μ positive constants. Then, for every ξ and η in *
*we have*

**Theorem 6.1 **(Stability) *For δ *= 1, *problem (4.18) is stable for all θ *∈ (0,1).

*Proof *Consider first the case *γ *= 1, *β *= 1. From the definition of *A*_{stab}(·,·;·,·;·), we have

We remark that (6.4) can be rewritten as

and the stability in the norm (5.2) follows from

Consider now the case *γ *= 0, *β *=1. Using the equivalent expressions (4.11) and (4.12) for the first and second equation
of (4.3), respectively, the problem (4.18) for *γ *= 0 can be rewritten as: find **u**_{h }∈ (*D*_{l-2}(*T _{h}*))

^{d},

*p*∈

_{h }*D*

_{l-1}(

*T*) such that

_{h}

From the first equation in (6.6) and (4.9) we have

Substituting the expression (6.7) in the second equation of (6.6) for *δ *= 1, we have

Denote by *B _{1h}*(·,·) the bilinear form (6.8), we have

and the stability in the norm (5.3) follows from Lemma 6.1. This completes the proof. □

**Theorem 6.2 **For *δ *= -1, problem (4.18) is stable for all *θ *< 0.

*Proof *Consider first the case *γ *= 1, *β *= 1. The problem (4.18) for *δ *= -1 reads

Using the arithmetic-geometric mean inequality, we have

and since *θ *< 0 the result follows.

Consider now the case *γ *= 0, *β *= 1. From (6.7) the second equation of (6.6) for *δ *= -1 can be written as

We remark that formulation (6.12) can be rewritten as

where *A _{BO}*(

*p*,

_{h}*ψ*) is introduced by Baumann and Oden [14], and given by

Denote by *B*_{2h}(·,·) the bilinear form (6.13), we have

and since *θ *< 0 the result follows again from Lemma 6.3 and 6.1. □

**Theorem 6.3 **(Consistency) *If p,c and ***u ***are the solution of (2.1)-(2.5) and are essentially bounded, then*

*provided that the constant M for the cut-off operator is sufficiently large*.

To summarize, for all the bilinear forms in (6.4), (6.10), (6.8) or (6.13) we have:
∃*C *> 0 such that

and ∃*C *> 0 such that

where (6.17) clearly holds for every *θ *∈ (0,1) for the case ((6.4), (6.8)), and for every *θ *< 0 for the case ((6.10), (6.13)). On the other hand, since ∇_{h}*D*_{l-1}(*T _{h}*) ⊂ (

*D*

_{l-2}(

*T*))

_{h}^{d }holds, boundedness of the bilinear form in (6.8) and (6.13) follows directly from the boundedness of the bilinear forms

*A*

_{BR }and

*A*

_{BO}, as proved in [13], thanks to the equivalence of the norms (6.1) and (6.2). Thus, we have: ∃

*C*> 0 such that

### 7 Error estimates

Let
**u**, *p*, *c*) such that

Let us define interpolation errors, finite element solution errors and auxiliary errors

It was proven in [18] that

hold for all *t *∈ *J *with the constant *C *independent only on bounds for the coefficient *α*(*c*), but not on *c *itself.

**Theorem 7.1 **(Error estimate for the velocity and pressure) *Let (***u**, *p, c) be the solution to (2.1)-(2.5), and assume p *∈ *L*^{2}(0, *T*; *H ^{l}*(

*T*)),

_{h}**u**∈ (

*L*

^{2}(0,

*T*;

*H*

^{l-1}(

*T*)))

_{h}*∈*

^{d }and c*L*

^{2}(0,

*T*;

*H*(

^{k}*T*)).

_{h}*We further assume that p*, ∇

*p*,

*c and*∇

*c are essentially bounded. If the constant M for the cut-off operator is sufficiently large, then there exists a constant C independent of h such that*

*Proof *For the sake of brevity we will assume
*γ *= 1, *β *= 1. From the second equation of (5.1) and (6.16) we have

That is

Choosing *v *= *ξ*_{1}, *ψ *= *η*_{1 }and splitting *e _{p }*according

*e*=

_{p }*η*

_{1}-

*η*

_{2}, from (7.1) and we obtain

Let us first consider the left side of error equation (7.5)

We know that (7.2) and quasi-regularity that
*L*^{∞}(Ω). So the right side of the error equation (7.5) can be bounded from below. Noting
that

The second and the third terms of the right side of the error equation (7.5) can be bounded using Cauchy-Schwartz inequality and approximation results,

The fourth term can be bounded in a similar way as that for the first term

The last two terms can be bounded as follows

Substituting all these inequalities into Equation 7.5, we have

Using the facts

The theorem follows from the triangle inequality.

Now consider the case *γ *= 0, *β *= 1. The bilinear form (6.8) from the second equation of (5.1) for

where

Replacing (6.8) with *p _{h }*=

*p*and subtracting it from (7.13) we finally obtain

Choosing *ψ *= *η*_{1}, we have

Let us first estimate the left side of (7.15). From (6.17) and using the fact

The first and the second terms of the right side of (7.15) can be bounded using Lemma 6.1 and (3.5)

Note that
*L*^{∞}(Ω) and

Substituting all these inequalities into the (7.15) and using the triangle inequality we have

We easily deduce,using (7.19)

which completes the proof. □.

From [7], we state two lemmas for the properties of the dispersion-diffusion tensor, which will be used to prove error estimates for the concentration.

**Lemma 7.1 **[7] (Uniform positive definiteness of **D**(**u**)) *Let ***D**(**u**) *defined as in Equation 2.6, where ϕd _{m}*(

*x*) ≥ 0,

*d*(

_{l}*x*) ≥ 0

*and d*(

_{t}*x*) ≥ 0

*are non-negative functions of x*∈ Ω.

*Then*

*If, in addition, φd _{m}*(

*x*) ≥

*d*

_{m,* }> 0

*uniformly in the domain*Ω,

*then*

**D**(

**u**)

*is uniformly positive definite in*Ω:

**Lemma 7.2 **[7] (Uniform Lipschitz continuity of **D(u)**) *Let ***D(u) ***defined as in Equation 2.6, where d _{m}*(

*x*) ≥ 0,

*d*(

_{l}*x*) ≥ 0

*and d*(

_{t}*x*) ≥ 0

*are non-negative functions of x*∈ Ω,

*and the dispersivities d*(

_{l }and d_{t }is uniformly bounded, i.e., d_{l}*x*) ≤

*d*

_{l}^{* }

*and d*(

_{t}*x*) ≤

*d*

_{t}^{*}.

*Then*

*where *
*is a fixed number (d = 2 or 3 is the dimension of domain *Ω*)*.

**Theorem 7.2 **(Error estimate for concentration) *Let (***u**, *p*, *c) be the solution to (2.1)-(2.5), and assume p *∈ *L*^{2}(0, *T*; *H ^{l}*(

*T*)),

_{h}**u**∈ (

*L*

^{2}(0,

*T*;

*H*

^{l-1}(

*T*)))

_{h}*∈*

^{d }and c*L*

^{2}(0,

*T*;

*H*(

^{k}*T*)).

_{h}*We further assume that p*, ∇

*p*,

*c and*∇

*c are essentially bounded. If the constant M for the cut-off operator is sufficiently large, then there exists a constant C independent of h and*Δ

*t such that*

*Proof *The first equation of (5.1) is

It can be written as

Subtracting the DG scheme equation from the weak formulation, we have for any *w *∈ *D*_{k-1}(*T _{h}*)

that is

Choosing *w *= *τ*_{1}^{n+1}, we obtain

Let us first consider the left side of the error equation (7.25). The first term can be bounded as

The second term of Equation 7.25 is

The second term of *B*(·,·;·) can be estimated using the boundedness of **u**_{M }and

Thus

Let us bound the right side of the error equation (7.25).

Using Taylor series expansion, we have

The fourth term in the right side of the error equation (7.25) is