A Broken 𝑃₁-Nonconforming Finite Element Method for Incompressible Miscible Displacement Problem in Porous Media

Abstract

An approximate scheme is defined for incompressible miscible displacement in porous media. This scheme is constructed by using immersed interface finite element method for the pressure equation which is based on the broken 𝑃₁-nonconforming piecewise linear polynomials on interface triangular elements and utilizing finite element method for the concentration equation. Error estimates for pressure in broken 𝐻¹ norm and for concentration in 𝐿² norm are presented.

1. Introduction

Miscible displacement of one incompressible fluid in a porous medium Ω over time interval 𝐽=[0,𝑇] is modeled by the system (a)𝛼(𝑐)(−∇⋅∇𝑝=∇⋅𝑢=𝑞,𝑥∈Ω,𝑡∈𝐽,b)𝜙𝜕𝑐(𝐷(𝑢))=(̃)()(𝜕𝑡+𝑢⋅∇𝑐−∇⋅∇𝑐𝑐−𝑐𝑞=𝑔𝑥,𝑡,𝑐,𝑥∈Ω,𝑡∈𝐽,c)(𝑢)(𝑢⋅𝑛=𝐷∇𝑐⋅𝑛=0,𝑥∈𝜕Ω,𝑡∈𝐽,d)𝑐(𝑥,0)=𝑐₀(𝑥),𝑥∈Ω. (1)

In the formulation, we have utilized the fact that in applications the vertical dimension of a subsurface geological formation is often much smaller than the horizontal dimension. We thus simplify the problem via a vertical average to rewrite the problem as a two-dimensional problem in the horizontal dimension, so the physical domain Ω⊂𝑅². 𝛼(𝑐)=𝛼(𝑥,𝑐)=𝑘(𝑥)/𝜇(𝑐), 𝑘(𝑥) is the permeability tensor of the medium and 𝜇(𝑐) is the viscosity of the fluid mixture that may be discontinuous across some interfaces; 𝜙(𝑥) is the porosity of the medium; 𝑢(𝑥,𝑡) is the Darcy velocity of the mixture; 𝑞(𝑥,𝑡) represents flow rates at wells, commonly a linear combination of Dirac measures; 𝐷(𝑢)=𝐷(𝑥,𝑐,∇𝑝)=𝜙(𝑥){𝑑_𝑚𝐼+|𝑢|(𝑑_𝑙𝐸(𝑢)+𝑑_𝑡𝐸^𝑇(𝑢))} is the diffusion-dispersion tensor, with 𝑑_𝑚, 𝑑_𝑙, and 𝑑_𝑡 being the molecular diffusion, the transverse, and longitudinal dispersivities, respectively, 𝐼 is the identity tensor; 𝑐(̃𝑥,𝑡) is specified at sources and 𝑐(̃𝑥,𝑡)=𝑐(𝑥,𝑡) at sinks. 𝑐₀(𝑥,𝑡) is the initial concentration. The dependent variable is 𝑝(𝑥,𝑡) is the pressure in the fluid mixture, and 𝑐(𝑥,𝑡) is the concentration of a solvent injected into resident reservoir.

For (1) we assume that 𝑞,1=_Ω𝑞(𝑥,𝑡)𝑑𝑥=0,𝑡∈𝐽. (2)

Let 𝑢=𝑢(𝑥,𝑐,∇𝑝)=(𝑢₁(𝑥,𝑐,𝑝_𝑥),𝑢₂(𝑥,𝑐,𝑝_𝑦)); for some 𝜀̂>0 restrict the variable 𝑞₁ to lie between −𝜀̂⩽𝑞₁⩽1+𝜀̂ and 𝑞₂∈𝑅, so we can suppose that the following regularity for 𝛼, 𝐷, 𝜙, and 𝑢 holds: 0<𝛼_∗⩽𝛼𝑥,𝑞₁⩽𝛼^∗,𝐷_∗⩽𝐷(𝑢)⩽𝐷^∗,𝜙_∗⩽𝜙(𝑥)⩽𝜙^∗;𝑢_𝑖𝑥,𝑞₁,𝑞₂⩽𝐾1+𝑞₂,𝑖=1,2,𝑞₂∈𝑅. (3)(4)

There exist much literature concerning numerical method and numerical analysis of miscible placement problem in porous media, for example, [1–4], and so on. As we know the porous media equations used to model the interface between oil and an injected fluid in simulations of secondary recovery in oil reservoirs. For the pressure equation (1)(a), the coefficient 𝛼(𝑐) often changes rapidly across fluid interfaces, and this sharp change is accompanied by large changes in the pressure gradient, as a compensatory, yielding a fairly smooth Darcy velocity 𝑢, so we can trade (1)(a) as an interface problem. When the interface is smooth enough, the solution of the interface problem is also very smooth in individual regions where the coefficient is smooth, but due to the jump of the coefficient across the interface, the global regularity is usually very low and has order of 𝐻^1+𝛾, 0<𝛾<1. Due to the low global regularity and the irregular geometry of the interface, it seems to be difficult for the standard finite element method to achieve high accuracy.

For better approximation, the fitted finite element methods whose mesh depends on the smooth interface are developed [5]. However, this method using fitted grids is costly for more complicated time dependent problems in which the interface moves with time and repeated grid generation is called for. Compared with the fitted finite element methods, the immersed interface method proposed by LeVeque and Li [6] allows the mesh to be independent of the interface, such as a Cartesian mesh. In recent years, Li et al. [7] studied an immersed finite element method using uniform grid, and they proved the approximation property of the finite element space of this scheme. On the other hand, Kwak et al. [8] introduced an immersed finite element method based on the broken 𝑃₁-nonconforming piecewise linear polynomials on interface triangular elements; this method uses edge averages as degrees of freedom, and the basis functions are C-R type [9]. Theory analysis and numerical experiments also show the optimal-order convergence of the method.

In this paper, we apply 𝑃₁-nonconforming finite element method (𝑃₁FEM) to pressure equation, while a FEM was used to approximate the concentration equation. Other methods could also be employed for the discretization of the concentration equation, for example, characteristic Galerkin method and so on. However, since the main point of this paper is to show the feasibility of the use of 𝑃₁FEM for pressure, we will discuss the concentration equation in the single case.

The rest of the paper is organized as follows: in Section 2, we first introduce some preliminaries. In Section 3, we briefly describe the 𝑃₁FEM-FEM scheme. In Section 4, we present some projections and lemmas used in the following error analysis. In Section 5, we prove the main error estimate.

2. Preliminaries

In this paper, without loss of generality, we consider the case in which Ω⊂𝑅² is a rectangular domain and the interface Γ is a smooth curve separating Ω into two subdomains Ω⁺­ and Ω⁻ such that Ω=Ω⁺⋃Ω⁻⋃Γ; see Figure 1 for an illustration.

Figure 1 

The proper jump conditions for elliptic (1)(a) are given by 𝑝=0,𝛼(𝑐)∇𝑝⋅𝑛=0acrossΓ, (5) where 𝑝∈𝐻¹(Ω). We assume 𝛼(𝑐) is a positive function bounded below and above by two positive constants.

Let 𝑊^𝑚_𝑞(Ω), 1⩽𝑞⩽+∞, be the Sobolev spaces consisting of functions whose derivatives up to order-𝑚 are 𝑞th integrable on Ω, and 𝐻^𝑚(Ω)∶=𝑊^𝑚₂(Ω) and 𝐻⁰(Ω)=𝐿²(Ω).

For the analysis, we introduce the space 𝐻²(Ω)=𝑝∈𝐻¹(Ω);𝑝∈𝐻²Ω^𝑠,𝑠=+,−, (6) equipped with the norm 𝑝²𝐻²_(Ω)=𝑝²_𝐻2(Ω+)+𝑝²_𝐻2(Ω−),𝑝²𝐻²_(Ω)=𝑝²_𝐻2(Ω+)+𝑝²_𝐻2(Ω−),∀𝑝∈𝐻²(Ω). (7)

We integrate (1)(a) multiplied by any test function 𝑣∈𝐻¹(Ω) over the domain Ω and apply the divergence theorem; then the variational formulation of the pressure equation is as follows Find 𝑝∈𝐻¹(Ω) such that 𝑎𝑝,𝑣=𝑞,𝑣,𝑣∈𝐻¹(Ω), (8) where the bilinear formation 𝑎(𝑝,𝑣)=(𝛼(𝑐)∇𝑝,∇𝑣).

By the Sobolev embedding theorem, for any 𝑝∈𝐻²(Ω), we can get 𝑝∈𝑊¹_𝑠(Ω) for any 𝑠>2. Then we have the following regularity lemma for the weak solution 𝑝 of the variational formulation (8).

Lemma 1 (see [10]). Assume that 𝑞∈𝐿²(Ω). Then there exists a unique solution 𝑝∈𝐻²(Ω) to the variational formulation (8) such that 𝑝𝐻²_(Ω)⩽𝑞_𝐿2(Ω). (9)

3. Formation of the Method

In this section, firstly, we define a broken 𝑃₁-nonconforming finite element method for pressure equation, and then, we use a finite element method to approximate the concentration equation.

In order to construct a broken 𝑃₁-nonconforming finite element procedure for pressure equation (1)(a), we assume the following situation.

Let 𝑇_ℎ be the usual regular triangulation of the domain Ω such that the elements have diameters bounded by ℎ_𝑝. Without loss of generality, we assume that the triangles in the partition used have the following features.(𝐻₁)If Γ meets one edge at more than two points, then this edge is one part of Γ. (𝐻₂)If Γ meets a triangle at two points, then the two points must be on the different edges of this triangle. (𝐻₃)The interface curve Γ is defined by a piecewise 𝐶² function, and the mesh 𝑇_ℎ is formed such that the subset of Γ in any interface element is 𝐶².

Then we can separate the triangles of partition 𝑇_ℎ into two classes. For an element 𝑇∈𝑇_ℎ, if the interface Γ passes through the interior of 𝑇, we call it an interface element and denote it by 𝑇_𝑚; otherwise, we call it a noninterface element and denote it by 𝑇_𝑛, respectively.

As usual, we will define linear finite element spaces on each element of the partition 𝑇_ℎ, so we have to construct the local basis functions. For a noninterface element 𝑇_𝑛, we can simply use the standard linear shape functions on 𝑇_𝑛 with degrees of freedom at average values along edges 𝑒_𝑗 of 𝑇_𝑛 and construct the linear finite element space 𝑆_ℎ(𝑇_𝑛) as follows: 𝑆_ℎ𝑇_𝑛=span𝜑_𝑖∶𝜑_𝑖∈𝑃₁𝑇_𝑛,1𝑒_𝑗^𝑒_𝑗𝜑_𝑖𝑑𝑠=𝛿_𝑖𝑗,𝑖,𝑗=1,2,3. (10)

For this space, we define the interpolation operator 𝐼_ℎ∶𝐻²(𝑇_𝑛)→𝑆_ℎ(𝑇_𝑛), where the following well-known approximation property is satisfied: 𝑝−𝐼_ℎ𝑝_{𝐿2(𝑇𝑛)}+ℎ_𝑝𝑝−𝐼_ℎ𝑝_{𝐻1(𝑇𝑛)}⩽𝐶ℎ²_𝑝𝑝_{𝐻2(𝑇𝑛)}, (11) and we use 𝑆_ℎ(Ω_𝑛) to denote the space of piecewise 𝑃₁-nonconforming finite element space on the domain Ω_𝑛=⋃𝑇_𝑛.

For interface element 𝑇_𝑚, since the finite element space on a general element can be obtained from the counterpart of a reference element through an affine mapping, we consider a typical interface element 𝑇_𝑚 whose geometric configuration is given in Figure 2 in which the three vertices are given by 𝐴₁=(0,0), 𝐴₂=(1,0), and 𝐴₃=(0,1), and the curve between points 𝐷 and 𝐸 is part of the interface with 𝐷=(0,𝑎) for 0<𝑎⩽1 and 𝐸=(𝑏,0) for 0<𝑏⩽1. Let 𝐷𝐸 denote the line segment connecting points 𝐷 and 𝐸, which divides 𝑇_𝑚 into two parts 𝑇⁺_𝑚 and 𝑇⁻_𝑚 with 𝑇_𝑚=𝑇⁺_𝑚⋃𝑇⁻_𝑚⋃𝐷𝐸.

Figure 2 

Let 𝑒_𝑖, 𝑖=1,2,3 be the edges of 𝑇_𝑚, and let 𝜑_𝑒𝑖 denote the average of 𝜑 along 𝑒_𝑖, that is, 𝜑_𝑒𝑖=(1/|𝑒_𝑖|)∫_𝑒𝑖𝜑𝑑𝑠, 𝜑∈𝐻¹(𝑇_𝑚), and then we would like to construct a new function which is linear on 𝑇⁺_𝑚 and 𝑇⁻_𝑚, respectively, and satisfies the jump condition (5) on 𝐷𝐸. For this purpose, we write the modified basis function 𝜑 on an interface element 𝑇_𝑚 as follows: 𝜑𝜑=⁻=𝑎₀+𝑏₀𝑥+𝑐₀𝑦,in𝑇⁻_𝑚,𝜑⁺=𝑎₁+𝑏₁𝑥+𝑐₁𝑦,in𝑇⁺_𝑚, (12) with the following constraints: 𝜑_𝑒𝑖=𝑉_𝑖𝜑,𝑖=1,2,3,₋(𝐷)=𝜑⁺(𝐷),𝜑⁻(𝐸)=𝜑⁺(𝐸),𝛼⁻𝜕𝜑⁻𝜕𝑛_𝐷𝐸=𝛼⁺𝜕𝜑⁺𝜕𝑛_𝐷𝐸, (13) where 𝛼⁻, 𝛼⁺ are averages along 𝐷𝐸 and 𝜌=𝛼⁻/𝛼⁺, 𝑉_𝑖, 𝑖=1,2,3, are given values, and 𝑛_𝐷𝐸 is the unit normal vector on the line 𝐷𝐸, that is, 𝑛_𝐷𝐸=(𝑎,𝑏)/𝑎²+𝑏².

By similar calculation to that given in Theorem 2.2 of [8], we can rewrite (12) and (13) in the matrix form: 𝐴𝑎₀𝑏₀𝑐₀𝑎₁𝑏₁𝑐₁=𝑉₁𝑉₂𝑉₃000, (14) where the coefficient matrix is defined by 1𝐴=00012121𝑎02𝑎²11−𝑎021−𝑎²𝑏12𝑏²101−𝑏21−𝑏²010𝑎−10−𝑎ℎ₂,1𝑏0−1−𝑏00𝜌𝑎𝜌𝑏0−𝑎−𝑏 (15) and the determinant of the matrix 𝐴 is (𝐴)1det=−4𝑎³1𝑏+4𝑎³1𝑏𝜌−4𝑎²𝜌−14𝑎𝑏³+14𝑎𝑏³1𝜌−4𝑏²𝜌=14𝑎²+𝑏²𝜌(𝑎𝑏−1)−𝑎𝑏<0. (16)

Therefore, the coefficients of (12) are uniquely determined by conditions 𝜑_𝑒𝑖=𝑉_𝑖, 𝑖=1,2,3, respectively, and when 𝑉_𝑖, 𝑖=1,2,3, have the same value 𝑉, we can get the piecewise linear function 𝜑=𝑉 by uniqueness.

Therefore we can define the following finite element space on an interface element 𝑇_𝑚: 𝑆_ℎ𝑇_𝑚=span𝜑∶𝜑|_𝑇𝑚iswelldefinedbytheaboveconstruction. (17)

Also we use 𝑆_ℎ(Ω_𝑚) to denote the finite element space defined on the domain Ω_𝑚=⋃𝑇_𝑚.

Remark 2. Although for functions in 𝑆_ℎ(𝑇_𝑚), the flux jump condition is enforced on line segments, they actually satisfy a weak flux jump condition along the interface Γ when 𝛼 is a piecewise constant such that _Γ∩𝑇𝛼⁻∇𝜑⁻−𝛼⁺∇𝜑⁺⋅𝑛_Γ𝑑𝑠=0,∀𝜑∈𝑆_ℎ𝑇_𝑚, (18) which is proved by an application of divergence theorem as in [11].

Combining the definitions of 𝑆_ℎ(Ω_𝑛) and 𝑆_ℎ(Ω_𝑚),we can describe the immersed finite element space 𝑆_ℎ(Ω) on the whole domain Ω, 𝑆_ℎ(Ω)=𝜑∶𝜑|_Ω𝑛∈𝑆_ℎΩ_𝑛,𝜑|_Ω𝑚∈𝑆_ℎΩ_𝑚, (19) endowed with the broken norm, 𝜑²_1,ℎ=_𝑇𝑇∇𝜑²𝑑𝑥,𝜑∈𝑆_ℎ(Ω). (20)

Now, assume that the approximation of concentration 𝐶 is known. Then, the pressure 𝑃∈𝑆_ℎ(Ω) can be determined by the following system: 𝑎_ℎ𝑃,𝜑=𝑞,𝜑,𝜑∈𝑆_ℎ(Ω), (21) where 𝑎_ℎ(𝑃,𝜑)=∑_{𝑇∈𝑇ℎ}∫_𝑇𝛼(𝐶)∇𝑃∇𝜑𝑑𝑥 is the bilinear formulation defined on 𝐻_ℎ(Ω)×𝐻_ℎ(Ω) with 𝐻_ℎ(Ω)=𝐻¹(Ω)⊕𝑆_ℎ(Ω).

The question at hand is to discretize the concentration equation.

Multiplying the concentration equation (1)(b) by a test function 𝑣∈𝐻¹(Ω) and integrating over Ω which leads to the following weak formulation for concentration 𝑐: 𝜙𝜕𝑐+𝑢+𝐷=𝑔(𝑐)𝜕𝑡,𝑣𝑐,∇𝑝⋅∇𝑐,𝑣𝑐,∇𝑝∇𝑐,∇𝑣,𝑣,𝑣∈𝐻¹(Ω). (22)

Let 𝑀_ℎ denote the standard piecewise linear finite element space associated with 𝑇_ℎ, and ℎ_𝑐 denote the bound of elements diameters. Then the discrete procedure of (22) is to find 𝐶∶𝐽→𝑀_ℎ such that 𝜙𝜕𝐶𝜕𝑡,𝑣_ℎ+𝑢(𝐶,∇𝑃)⋅∇𝐶,𝑣_ℎ+𝐷(𝐶,∇𝑃)∇𝐶,∇𝑣_ℎ=𝑔(𝐶),𝑣_ℎ,𝑣_ℎ∈𝑀_ℎ,𝐶(𝑥,0)=𝑐₀̃(𝑥),𝑥∈Ω, (23) where 𝑐₀̃(𝑥) is the elliptic projection of 𝑐₀(𝑥).

4. Projection and Some Lemmas

For convergence analysis conducted in the following section, we need an interpolation operator 𝐼_ℎ∶𝐻²(𝑇)→𝑆_ℎ(𝑇) using the average of 𝑝 on each edge by 𝐼_ℎ𝑝_𝑒𝑖=𝑝_𝑒𝑖,𝑖=1,2,3, (24) and when 𝑝∈𝐻²(Ω), we can define 𝐼_ℎ𝑝 by (𝐼_ℎ𝑝)|_𝑇=𝐼_ℎ(𝑝|_𝑇). Then the following interpolation estimate hold for the 𝑃₁FEM approximation [8] 𝑝−𝐼_ℎ𝑝_𝐿2(Ω)+ℎ_𝑝𝑝−𝐼_ℎ𝑝_1,ℎ⩽𝐶ℎ²_𝑝𝑝𝐻²_(Ω),∀𝑝∈𝐻²(Ω). (25)

Next, let 𝐶∶𝐽→𝑀_ℎ be the elliptic projection of 𝑐 defined by 𝑢𝑐,∇𝑝⋅∇𝐶−𝑐,𝑣_ℎ+𝐷𝑐,∇𝑝∇𝐶−𝑐,∇𝑣_ℎ+𝜆𝐶−𝑐,𝑣_ℎ=0,∀𝑣_ℎ∈𝑀_ℎ. (26)

The function 𝜆 will be chosen to assure coercivity of the form. Then, the following estimates hold [12]: 𝑐−𝐶_𝐿2(Ω)+ℎ_𝑐𝑐−𝐶_𝐻1(Ω)⩽𝐾ℎ²_𝑐‖𝑐‖_𝐻²_(Ω),𝜕(𝑐−𝐶)𝜕𝑡_𝐿2(Ω)⩽𝐾ℎ²_𝑐‖𝑐‖_𝐻²_(Ω)+𝜕𝑐𝜕𝑡_𝐻²_(Ω). (27)

For convergence analysis conducted in Section 5, we need to apply quote some lemmas from [9, 13] and references therein.

Lemma 3 (see [9] (the second Strang lemma)). Let 𝑝 and 𝑃 be the solution of (8) and (21). Then there exists a constant 𝐶>0 such that 𝑝−𝑃_1,ℎ⩽𝐶inf_𝑤ℎ∈𝑆_ℎ(Ω)𝑝−𝑤_ℎ1,ℎ𝜑+sup_∈𝑆_ℎ(Ω)𝑎_ℎ𝑝,𝜑−𝑞,𝜑𝜑_1,ℎ. (28)

Lemma 4 (see [13]). Let 𝑒 be an edge of  𝑇. Then there exists a constant 𝐶>0 such that for all 𝜑,𝑣∈𝐻¹(𝑇): _𝑒𝜑𝑣−𝑣_𝑒𝑑𝑠⩽𝐶ℎ_𝑝𝜑_1,𝑇|𝑣|_1,𝑇, (29) where 𝑣_𝑒∶=(1/|𝑒|)∫_𝑒𝑣𝑑𝑠.

5. Convergence Analysis

In this section, we will present convergence analysis for the pressure and concentration.

Theorem 5. Let 𝑝∈𝐻²(Ω), 𝑃∈𝑆_ℎ(Ω) be the solution of (8) and (21); then there exists a constant 𝐾>0 independent of ℎ_𝑝 and the location of the interface, such that the following result holds: 𝑝−𝑃_1,ℎ⩽𝐾ℎ_𝑝𝑝𝐻²_(Ω)+𝐾‖𝑐−𝐶‖_𝐿²_(Ω). (30)

Proof. Since the immersed finite element formulation (21) is nonconforming, we can use Lemma 3 to prove the error bound. The first part in (28) is an approximation error, which can be estimated simply: inf_𝑤ℎ∈𝑆_ℎ(Ω)𝑝−𝑤_ℎ1,ℎ⩽𝐾ℎ_𝑝𝑝𝐻²_(Ω). (31)
For the second part, it is a consistency error estimate. By the definition of 𝑎_ℎ(⋅,⋅) and Green’s formula, we can get 𝑎_ℎ𝜑−𝜑=𝑝,𝑞,_{𝑇∈𝑇ℎ𝑇}𝛼(𝐶)∇𝑝∇𝜑𝑑𝑥−_Ω𝑞=𝜑𝑑𝑥_{𝑇∈𝑇ℎ𝑇}𝛼(𝐶)∇𝑝∇𝜑𝑑𝑥−_{𝑇∈𝑇ℎ𝑇}𝛼(𝑐)∇𝑝∇𝜑𝑑𝑥+_{𝑇∈𝑇ℎ}𝛼(𝑐)𝜕𝑝,𝜑𝜕𝑛_𝜕𝑇⩽_{𝑇∈𝑇ℎ𝑇}(𝛼(𝐶)(𝑐+−𝛼))∇𝑝∇𝜑𝑑𝑥_{𝑇∈𝑇ℎ}𝛼(𝑐)𝜕𝑝,𝜑𝜕𝑛_𝜕𝑇⩽∇𝑝_𝐿∞(Ω)‖𝑐−𝐶‖_𝐿²_(Ω)𝜑_1,ℎ+_{𝑇∈𝑇ℎ}𝛼(𝑐)𝜕𝑝,𝜑𝜕𝑛_𝜕𝑇⩽𝐾‖𝑐−𝐶‖_𝐿²_(Ω)𝜑_1,ℎ+_{𝑇∈𝑇ℎ}𝛼(𝑐)𝜕𝑝,𝜑𝜕𝑛_𝜕𝑇, (32) where 𝑛 is the unit outward normal vector on each 𝜕𝑇, 𝜑∈𝑆_ℎ(Ω), and by the construction of the space 𝑆_ℎ(Ω) we have a well-defined property on the interior edges, that is, ∫_𝑒𝜑|_𝑇1=∫_𝑒𝜑|_𝑇2 for 𝑒 is the common edge of adjacent element 𝑇₁ and 𝑇₂. By boundary condition, we note that 𝛼(𝑐)(𝜕𝑝/𝜕𝑛) vanishing average on the boundary. So we rewrite ∑_{𝑇∈𝑇ℎ}<𝛼(𝑐)(𝜕𝑝/𝜕𝑛),𝜑>_𝜕𝑇 and use Lemma 4 to derive that _{𝑇∈𝑇ℎ}𝛼(𝑐)𝜕𝑝,𝜑𝜕𝑛_𝜕𝑇=_{𝑇∈𝑇ℎ𝑒∈𝜕𝑇}𝛼(𝑐)𝜕𝑝−𝜕𝑛𝛼(𝑐)𝜕𝑝,𝜑𝜕𝑛_𝑒⩽_{𝑇∈𝑇ℎ}𝐶ℎ_𝑝𝛼(𝑐)𝜕𝑝𝜕𝑛_1,𝑇𝜑_1,𝑇⩽𝐶ℎ_𝑝𝑝𝐻²_(Ω)𝜑_1,ℎ. (33)
Combining (28) with (31)–(33) yields the theorem result.

We now are in the position to prove the error estimate for concentration.

Theorem 6. Assume that the true solution (𝑐,𝑝) of (1) satisfies 𝑐∈𝐿^∞(𝐻²)∩𝐿²(𝐻²), 𝑝∈𝐿²(𝐻²). Let 𝐶 be the projection of 𝑐 satisfying (26) and let 𝐶 be the solution of (23). Suppose that the mesh parameters satisfy ℎ²_𝑐=𝑂(ℎ_𝑝). Then the following result holds: 𝐶−𝐶_{𝐿∞(𝐿2)}+∇𝐶−𝐶_𝐿2(𝐿2)⩽𝐾ℎ²_𝑐+ℎ_𝑝. (34)

Proof. We decompose the error 𝑐−𝐶 to 𝜉=𝐶−𝐶, and 𝜂=𝑐−𝐶. Then, combining (22), (23), and (26) results in the following error equation for the concentration: 𝜙𝜕𝜉𝜕𝑡,𝑣_ℎ+𝐷(𝐶,∇𝑃)∇𝜉,∇𝑣_ℎ=𝜙𝜕𝜂𝜕𝑡,𝑣_ℎ−𝜆𝜂,𝑣_ℎ+𝑔(𝐶)−𝑔(𝑐),𝑣_ℎ−𝑢(𝐶,∇𝑃)∇𝜉,𝑣_ℎ+𝐷𝑐,∇𝑝−𝐷(𝐶,∇𝑃)∇𝐶,∇𝑣_ℎ+𝑢𝑐,∇𝑝−𝑢(𝐶,∇𝑃)∇𝐶,𝑣_ℎ. (35)
We choose the test function as 𝑣_ℎ=𝜉 in (35) and integrate from 0 to 𝑡 on both sides to derive ^𝑡₀𝜙𝜕𝜉𝜕𝑡,𝜉𝑑𝜏+^𝑡₀(𝐷())=𝐶,∇𝑃∇𝜉,∇𝜉𝑑𝜏_𝑡0𝜙𝜕𝜂𝜕𝑡,𝜉𝑑𝜏−^𝑡₀+𝜆𝜂,𝜉𝑑𝜏_𝑡0𝑔(𝐶)(𝑐)−−𝑔,𝜉𝑑𝜏_𝑡0(𝑢())+𝐶,∇𝑃∇𝜉,𝜉𝑑𝜏_𝑡0𝐷𝑐,∇𝑝−𝐷(𝐶,∇𝑃)∇𝐶,∇𝜉𝑑𝜏+^𝑡₀𝑢𝑐,∇𝑝−𝑢(𝐶,∇𝑃)∇𝐶,𝜉𝑑𝜏. (36)
The terms on the left side can be easily bounded by ^𝑡₀𝜙𝜕𝜉𝜕𝑡,𝜉𝑑𝜏+^𝑡₀(𝐷())⩾1𝐶,∇𝑃∇𝜉,∇𝜉𝑑𝜏2𝜙^∗𝜉²+𝐷_∗^𝑡₀∇𝜉²𝑑𝜏. (37)
Next, we should estimate terms on the right side one by one.
By (27), it is trivial that ^𝑡₀𝜙𝜕𝜂ℎ𝜕𝑡,𝜉𝑑𝜏⩽𝐾⁴_𝑐+^𝑡₀𝜉²𝑑𝜏. (38)
Similarly, we can bound the second and third terms as follows: 𝜆𝜂,𝜉𝑑𝜏⩽𝐾ℎ^𝑡₀⁴_𝑐+^𝑡₀𝜉²𝑑𝜏_𝑡0,𝑔(𝐶)−𝑔(𝑐),𝜉𝑑𝜏⩽𝐾ℎ⁴_𝑐+^𝑡₀𝜉²𝑑𝜏, (39) provided that 𝑔 is Lipschitz continuous.
For the fifth term, we can derive the estimate by using the boundness of 𝐶 and (27), as well as (30), ^𝑡₀𝐷()∇∇𝐶𝑐,∇𝑝−𝐷𝐶,∇𝑃𝐶,∇𝜉𝑑𝜏⩽𝐾_{𝐿∞(𝐿∞)}^𝑡₀𝜉+𝜂+𝑝−𝑃_1,ℎ∇𝜉𝑑𝜏⩽𝐾^𝑡₀𝜉+ℎ²_𝑐+ℎ_𝑝∇𝜉𝑑𝜏⩽𝐾^𝑡₀𝜉²𝑑𝜏+ℎ⁴_𝑐+ℎ²_𝑝+𝐷_∗4^𝑡₀∇𝜉²𝑑𝜏. (40)
Similarly, we can derive ^𝑡₀𝑢()∇ℎ𝑐,∇𝑝−𝑢𝐶,∇𝑃𝐶,𝜉𝑑𝜏⩽𝐾⁴_𝑐+ℎ²_𝑝+^𝑡₀𝜉²𝑑𝜏. (41)
In order to prove the forth term, we need the following induction hypothesis: ‖∇𝑃‖_𝐿^∞_(𝐿^∞₎⩽𝐾^∗. (42)
As the statement in [1], we assume that for some 𝜀>0,𝜀⩽𝜀̂, −𝜀≤𝐶⩽1+𝜀 (43) and with this hypothesis and (4), we can get ^𝑡₀(𝑢())‖𝐶,∇𝑃∇𝜉,𝜉𝑑𝜏⩽𝐾‖∇𝑃_𝐿^∞_(𝐿^∞₎+1²×^𝑡₀𝜉²𝐷𝑑𝜏+^∗4^𝑡₀∇𝜉²𝑑𝜏⩽𝐾^𝑡₀𝜉²𝐷𝑑𝜏+^∗4^𝑡₀∇𝜉²𝑑𝜏. (44)
Therefore, collecting all the bounds derived above and using Gronwall inequality leads to 𝜉²_{𝐿∞(𝐿2)}+∇𝜉²_𝐿2(𝐿2)⩽𝐾ℎ⁴_𝑐+ℎ²_𝑝. (45)
To complete the argument, we have to check hypothesis (42). We use (25), (27), (30), (45), inverse inequality, and boundness of 𝐼_ℎ𝑝 to get that ‖∇𝑃‖_𝐿^∞_(𝐿^∞₎⩽∇𝑃−𝐼_ℎ𝑝_{𝐿∞(𝐿∞(𝑇))}∇𝐼_ℎ𝑝_{𝐿∞(𝐿∞(𝑇))}⩽𝐾+𝐾ℎ_𝑝⁻¹_{𝑇∈𝑇ℎ}∇𝑃−𝐼_ℎ𝑝_{𝐿∞(𝐿2(𝑇))}⩽𝐾+𝐾ℎ_𝑝⁻¹𝑃−𝑝_{𝐿∞(1,ℎ)}+𝑝−𝐼_ℎ𝑝_{𝐿∞(1,ℎ)}⩽𝐾+𝐾ℎ_𝑝⁻¹ℎ_𝑝+ℎ²_𝑐⩽𝐾^∗. (46)

Combining (27) with (30) and (45), we deduce the error estimates for the pressure and the concentration.

Theorem 7. Assume that the true solution (𝑐,𝑝) of (1) satisfies 𝑐∈𝐿^∞(𝐻²)∩𝐿²(𝐻²), 𝑝∈𝐿²(𝐻²). Let 𝐶, 𝑃 be the solution of (23) and (21). Suppose that the mesh parameters satisfy ℎ²_𝑐=𝑂(ℎ_𝑝). Then the following result holds: 𝑝−𝑃_{𝐿∞(1,ℎ)}+‖𝑐−𝐶‖_𝐿^∞_(𝐿²₎⩽𝐾ℎ²_𝑐+ℎ_𝑝,‖∇(𝑐−𝐶)‖_𝐿²_(𝐿²₎⩽𝐾ℎ_𝑐+ℎ_𝑝. (47)

Remark 8. In this paper, we just get a priori error estimates for the coefficient 𝐷=𝐷(𝑥,𝑐,∇𝑝). However, some mathematical models for miscible displacements, which are currently being used by oil companies, make the assumption that the coefficient 𝐷=𝐷(𝑥,𝑐) and it does not depend on ∇𝑝; with this assumption, we can prove the optimal order convergence similarly.

Acknowledgments

This work is supported in part by the National Natural Science Foundation of China (no. 10971254), National Natural Science Foundation of Shandong Province (no. Y2007A14), and Young Scientists Fund of Shandong Province (no. 2008BS01008).