A Second Order Characteristic Method for Approximating Incompressible Miscible Displacement in Porous Media

Abstract

An approximation scheme is defined for incompressible miscible displacement in porous media. This scheme is constructed by two methods. Under the regularity assumption for the pressure, cubic Hermite finite element method is used for the pressure equation, which ensures the approximation of the velocity smooth enough. A second order characteristic finite element method is presented to handle the material derivative term of the concentration equation. It is of second order accuracy in time increment, symmetric, and unconditionally stable. The optimal 𝐿²-norm error estimates are derived for the scalar concentration.

1. Introduction

In this paper, we will consider the following incompressible miscible displacement in porous media, which is governed by a coupled system of partial differential equations with initial and boundary values. The pressure is governed by an elliptic equation and the concentration is governed by a convection-diffusion equation [1–6] as follows:−∇⋅𝑘𝜙𝜕𝑐𝜇∇𝑝≡∇⋅𝑢=𝑞,𝑥∈Ω,𝑡∈𝐽,(1.1a)𝜕𝑡+𝑢𝑢⋅𝜈⋅∇=𝑐(−𝐷𝑐∇((𝑥𝑥,⋅)0(∇𝐷)𝑐=∇)⋅𝑐𝑐𝜈)₀==(0𝑐,̃−𝑥𝑐)∈𝑞̃,𝜕Ω𝑥,∈𝑡∈Ω,𝐽,𝑡∈𝐽,((11..11bc))(𝑥),𝑥∈Ω,(1.1d)where Ω is a bounded domain in ℛ², 𝐽=(0,𝑇], and 𝑞̃=max{𝑞,0} is nonzero at injection wells only. The variables in (1.1a)–(1.1d) are the pressure 𝑝(𝑥,𝑡) in the fluid mixture, the Darcy velocity 𝑢=(𝑢₁,𝑢₂)^, and the relative concentration 𝑐(𝑥,𝑡) of the injected fluid. The 𝜈 is the unit outward normal vector on boundary 𝜕Ω.

The coefficients and data in (1.1a)–(1.1d) are 𝑘(𝑥), the permeability of the porous media; 𝜇(𝑥), the viscosity of the fluid mixture; 𝑞(𝑥,𝑡), representing flow rates at wells; 𝜙(𝑥), the porosity of the rock; 𝑐̃(𝑥,𝑡), the injected concentration at injection wells (𝑞>0) and the resident concentration at production wells (𝑞<0); 𝑞̃=max(𝑞,0). Here, for the diffusion coefficient, we consider a dispersion-free case [5, 6] as follows: 𝐷=𝜙𝑑_𝑚𝐼,(1.2) where 𝑑_𝑚 is the molecular diffusivity and 𝐼 is a 2×2 identity matrix. Furthermore, a compatibility condition ∫_Ω𝑞(𝑥,𝑡)𝑑𝑥=0 must be imposed to determine the pressure.

The pressure equation is elliptic and easily handled, but the concentration equation is parabolic and normally convection dominated. It is well known that standard Galerkin scheme applied to the convection-dominated problems does not work well and produces excessive numerical diffusion or nonphysical oscillation. The characteristic method has been introduced to obtain better approximations for (1.1a)–(1.1d), such as characteristic finite element method [3–7], characteristic finite difference method [8], the modified of characteristic finite element method (MMOC-Galerkin) [9], and the Eulerian-Lagrangian localized adjoint method (ELLAM) [10].

We had considered a combined numerical approximation for (1.1) in [11]. Standard mixed finite element was used for Darcy velocity equation and a characteristics-mixed finite element method was presented for approximating the concentration equation. Characteristic approximation was applied to handle the convection term, and a lowest order mixed finite element spatial approximation was adopted to deal with the diffusion term. Thus, the scalar unknown concentration and the diffusive flux can be approximated simultaneously. This approximation conserves mass globally. The optimal 𝐿²-norm error estimates were derived.

It should be pointed out that the works mentioned above all gave one order accuracy in time increment Δ𝑡. That is to say, the first order characteristic method in time was analyzed. As for higher order characteristic method in time, Rui and Tabata [12] used the second order Runge-Kutta method to approximate the material derivative term for convection-diffusion problems. The scheme presented was of second order accuracy in time increment Δ𝑡, symmetric, and unconditionally stable. Optimal error estimates were proved in the framework of 𝐿²-theory. Numerical analysis of convection-diffusion-reaction problems with higher order characteristic/finite elements were analyzed in [13, 14], which extended the work [12]. The 𝑙^∞(𝐿²) error estimates of second order in time increment Δ𝑡 were obtained.

The goal of this paper is to present a second order characteristic finite element method in time increment to handle the material derivative term of the concentration equation of (1.1a)–(1.1d). It is organized as follows. In Section 2, we formulate the second order characteristic finite element method for the concentration and cubic Hermite finite element method for the pressure, respectively. Then, we present a combined approximation scheme. In Section 3, we analyze the stability of the approximation scheme for the concentration equation. In Section 4, we derive the optimal-order 𝐿²-norm error estimates for the scalar concentration. They are of second order accuracy in time increment, symmetric, and unconditionally stable. We conclude our results in Section 5.

2. Formulation of the Method

2.1. Statements and Assumptions

In this paper, we adopt notations and norms of usual Sobolev spaces. For periodic functions, we use the notation as in [4] as follows. Let: 𝑊^𝑘𝑝=𝑊^𝑘𝑝(Ω)=𝜓∶𝜕^𝛼𝜓𝜕𝑥^𝛼∈𝐿^𝑝(Ω)for|𝛼|≤𝑘,periodic(2.1) be the periodic Sobolev space on Ω with the usual norm. If 𝑝=2, we write 𝐻^𝑘=𝐻^𝑘(Ω)=𝑊^𝑘2(Ω)(2.2) with norm ‖𝜓‖_𝑘=‖𝜓‖𝐻^𝑘_(Ω),‖𝜓‖=‖𝜓‖₀=‖𝜓‖_𝐿²_(Ω).(2.3)

Moreover, we adopt some notations for the functional spaces involved, which were used in [12–14]. For a Banach space 𝑋 and a positive integer 𝑚, spaces 𝐶^𝑚([0,𝑇],𝑋) and 𝐻^𝑚((0,𝑇),𝑋) are abbreviated as 𝐶^𝑚(𝑋) and 𝐻^𝑚(𝑋), respectively, and endowed with norms ‖𝜑‖_𝐶^𝑚_(𝑋)∶=max_{[]𝑡∈0,𝑇}max_{𝑗=0,…,𝑚}‖‖𝜑^(𝑗)(𝑡)‖‖_𝑋,‖𝜑‖^{𝐻𝑚(𝑋)}∶=^𝑇0𝑚_𝑗=0‖‖𝜑^(𝑗)(𝑡)‖‖^2𝑋𝑑𝑡^1/2,(2.4) where 𝜑^(𝑗) denotes the 𝑗th derivative of 𝜑 with respect to time. The Banach space 𝑍^𝑚 is defined by 𝑍^𝑚=𝑓∈𝐶^𝑗𝐻^𝑚−𝑗(Ω);𝑗=0,…,𝑚,(2.5) equipped with the norm ‖𝜑‖_𝑍^𝑚∶=max‖𝜑‖_𝐶^𝑗_(𝐻^𝑚−𝑗₎;0≤𝑗≤𝑚.(2.6)

We also require the following assumptions on the coefficients in (1.1) [3]. Let 𝑎_∗,𝑎^∗,𝜙_∗,𝜙^∗, and 𝐾^∗ be positive constants such that 0<𝑎_∗≤𝑘(𝑥)𝜇(𝑥)≤𝑎^∗||𝑞̃(𝑥,𝑡)||+|||𝜕,𝑞̃0<𝜙_∗≤𝜙(𝑥)≤𝜙^∗,𝜕𝑡(𝑥,𝑡)|||≤𝐾^∗.(C) Other assumptions will be made in individual theorems as necessary.

2.2. A Cubic Hermite Element Method for the Pressure

The variational form for the pressure equation (1.1a) is equal to the following: 𝜇𝑘∇𝑝,∇𝜒=(𝑞(𝑡),𝜒),𝜒∈𝐻¹(Ω),𝑡∈𝐽.(2.7)

For ℎ_𝑝>0, we discretize (2.7) in space on a quasiuniform triangle mesh 𝒯_ℎ𝑝 of Ω with diameter of element ≤ℎ_𝑝. Let 𝑊_ℎ𝑝⊂𝐿²(Ω) be a cubic Hermite finite element space for this mesh. The finite element method for the pressure, given at a time 𝑡∈𝐽, consists of 𝑃∈𝑊_ℎ𝑝 such that 𝜇𝑘∇𝑃,∇𝜒=(𝑞(𝑡),𝜒),∀𝜒∈𝑊^ℎ𝑝.(2.8)

Since the left-hand side of (2.7) ((2.8), resp.) is a continuous and coercive bilinear form and the right-hand side of (2.7) ((2.8), resp.) is a continuous linear functional, existence and uniqueness of 𝑝 (𝑃, resp.) are ensured obviously [15].

By the theory of the cubic Hermit element, the finite element space 𝑊_ℎ𝑝 possess the following approximation property [15, 16]: inf_{𝜒∈𝑊ℎ𝑝}‖𝑔−𝜒‖_𝑠≤𝐾ℎ^𝑝4−𝑠‖𝑔‖₄,𝑠=0,1,∀𝑔∈𝐻⁴(Ω),(2.9) where 𝐾 is a positive constant independent of ℎ_𝑝.

Define Π: 𝐶²(Ω)→𝑊_ℎ𝑝 is the interplant operator. Then, we have the following theorem.

Theorem 2.1. Let 𝑝 and 𝑃 be the solutions of (1.1a)–(1.1d) and (2.8), respectively, and assume 𝑝∈𝐻⁴(Ω). Then, there exists a positive constant 𝐾 independent of ℎ_𝑝 such that ‖∇(𝑝−𝑃)‖≤𝐾ℎ^3𝑝‖𝑝‖₄.(2.10)

Proof. By Céa lemma, we have 𝑎‖∇(𝑝−𝑃)‖²≤𝑘=𝑘𝜇∇(𝑝−𝑃),∇(𝑝−𝑃)=𝜇𝑘∇(𝑝−𝑃),∇(𝑝−Π𝑝)+𝜇𝑘∇(𝑝−𝑃),∇(Π𝑝−𝑃)≤𝑎_∗𝜇∇(𝑝−𝑃),∇(𝑝−Π𝑝)‖∇(𝑝−𝑃)‖‖∇(𝑝−Π𝑃)‖.(2.11) Then by (2.9), we can derive ‖∇(𝑝−𝑃)‖≤𝐾ℎ^3𝑝‖𝑝‖₄.(2.12)

2.3. A Second Order Characteristic Method for the Concentration

Now, we define the characteristic lines associated with vector field 𝑢 and recall some classical properties satisfied by them. Thus, for given (𝑥,𝑡)∈Ω×[0,𝑇], the characteristic line through (𝑥,𝑡) is the vector function 𝑋_𝑒(𝑥,𝑡;⋅) solving the following initial value problem: 𝜕𝑋_𝑒𝜕𝜏(𝑥,𝑡;𝜏)=𝑣𝑋^𝑒(𝑥,𝑡;𝜏),𝜏,𝑋_𝑒(𝑥,𝑡;𝜏)=𝑥,(2.13) where 𝑣(𝑋_𝑒(𝑥,𝑡;𝜏),𝜏)∶=𝑢(𝑋_𝑒(𝑥,𝑡;𝜏),𝜏)/𝜙.

Next, assuming they exist, we denote by 𝐹_𝑒 (by 𝐿, resp.) the gradient of 𝑋_𝑒 (of 𝑣, resp.) with respect to the space variable 𝑥, that is, 𝐹_𝑒_𝑟𝑠(𝑥,𝑡;𝜏)∶=𝜕𝑋^𝑒_𝑟𝜕𝑥_𝑠(𝑥,𝑡;𝜏),𝐿_𝑟𝑠(𝑥,𝑡)∶=𝜕𝑣^𝑟𝜕𝑥_𝑠(𝑥,𝑡).(2.14) We adopted some propositions and lemmas from [13].

Proposition 2.2. If 𝑣∈𝐶⁰(𝐶^𝑛(Ω)) for 𝑛≥1 an integer, then 𝑋_𝑒∈𝐶⁰(Ω×[0,𝑇]×[0,𝑇]) and it is 𝐶^𝑛 with respect to the 𝑥 variable.

In order to compute second order approximations of matrices 𝐹_𝑒 and 𝐹^𝑒−1, we need the following equations: 𝜕𝐹_𝑒𝜕𝜏(𝑥,𝑡;𝜏)=𝐿𝑋^𝑒𝜕(𝑥,𝑡;𝜏),𝜏𝐹_𝑒(𝑥,𝑡;𝜏),₂𝐹_𝑒𝜕𝜏²(𝑥,𝑡;𝜏)=∇𝜕𝑣𝜕𝑡+𝐿𝑣𝑋^𝑒(𝑥,𝑡;𝜏),𝜏𝐹_𝑒(𝑥,𝑡;𝜏).(2.15)

Proposition 2.3. If 𝑣∈𝐶⁰(𝐶¹(Ω)), then ‖‖𝐹_𝑒(𝑥,𝑡;𝜏)‖‖≤𝑒^{‖𝑣‖𝐶01((𝐶Ω))|𝜏−𝑡|},∀𝑥∈Ω,𝑡,𝜏∈[0,𝑇].(2.16)

Proposition 2.4. If 𝑣∈𝐶⁰(𝐶²(Ω))∩𝐶¹(𝐶¹(Ω)), then 𝐹_𝑒 satisfies the Taylor expansions as 𝐹_𝑒(𝑥,𝑡;𝑠)=𝐼++(𝑠−𝑡)𝐿(𝑥,𝑡)_𝑡𝑠(𝜏−𝑠)∇𝜕𝑣𝜕𝑡+𝐿𝑣𝑋^𝑒(𝑥,𝑡;𝜏),𝜏𝐹_𝑒(𝑥,𝑡;𝜏)𝑑𝜏,(2.17) and its inverse, 𝐹^𝑒−1, satisfies the Liouville theorem as 𝐹^𝑒−1(𝑥,𝑡;𝑠)=𝐼−+(𝑠−𝑡)𝐿𝑋_𝑒(𝑥,𝑡;𝑠),𝑠_𝑡𝑠(𝜏−𝑡)∇𝜕𝑣𝜕𝑡+𝐿𝑣𝑋^𝑒(𝑥,𝑡;𝜏),𝜏𝐹_𝑒𝑋_𝑒(𝑥,𝑡;𝑠),𝑠;𝜏𝑑𝜏.(2.18)

By using Liouville’s theorem and the chain rule, we obtain 𝜕𝜕𝜏det𝐹^𝑒−1𝜕(𝑥,𝑡;𝜏)=−det𝐹^𝑒−1(𝑥,𝑡;𝜏)div𝑣𝑋_𝑒(𝑥,𝑡;𝜏),𝜏₂,𝜕𝜏²det𝐹^𝑒−1(𝑥,𝑡;𝜏)=−×ddeivt𝐹𝜕^𝑒−1𝑣(𝑥,𝑡;𝜏)(div𝑣)²𝑋_𝑒(𝑥,𝑡;𝜏),𝜏𝜕𝑡+𝐿𝑣𝑋^𝑒−𝐿⋅𝐿_𝑇(𝑥,𝑡;𝜏),𝜏𝑋_𝑒(𝑥,𝑡;𝜏),𝜏.(2.19)

Proposition 2.5. If 𝑣∈𝐶⁰(𝐶¹(Ω)), then det𝐹_𝑒(𝑥,𝑡;𝜏)≤𝑒^{‖𝑣‖𝐶01((𝐶Ω))|𝜏−𝑡|},∀𝑥∈Ω,𝑡,𝜏∈[0,𝑇].(2.20)

Proposition 2.6. If 𝑣∈𝐶⁰(𝐶²(Ω))∩𝐶¹(𝐶¹(Ω)), then det𝐹^𝑒−1 satisfies det𝐹^𝑒−1(𝑥,𝑡;𝑠)=1−(𝑠−𝑡)div𝑣(𝑥,𝑡)+^𝑡𝑠(𝜏−𝑠)𝜕²𝜕𝜏²det𝐹^𝑒−1(𝑥,𝑡;𝜏)𝑑𝜏.(2.21)

Variational Formulation
From the definition of the characteristic curves and by using the chain rule, it follows that 𝜙𝑑𝑐𝑑𝜏𝑋^𝑒(𝑥,𝑡;𝜏),𝜏=𝜙𝜕𝑐𝜕𝑡𝑋^𝑒(𝑥,𝑡;𝜏),𝜏+𝑢𝑋_𝑒(𝑥,𝑡;𝜏),𝜏⋅∇𝑐𝑋_𝑒(𝑥,𝑡;𝜏),𝜏.(2.22) By writing (1.1b) at point 𝑋_𝑒(𝑥,𝑡;𝜏) and time 𝜏 and using (2.22), we have 𝜙𝑑𝑐𝑑𝜏𝑋^𝑒(𝑥,𝑡;𝜏),𝜏−∇⋅(𝐷∇𝑐)𝑋_𝑒(𝑥,𝑡;𝜏),𝜏=((𝑐̃−𝑐)𝑞̃)𝑋_𝑒(𝑥,𝑡;𝜏),𝜏.(2.23)
Before giving a week formulation of (2.23), we adopted a lemma from [13], which can be considered as Green’s formula.

Lemma 2.7. Let 𝑋∶Ω→𝑋(Ω), 𝑋∈𝐶²(Ω) be an invertible vector valued function. Let 𝐹=∇𝑋 and assume that 𝐹⁻¹∈𝐶¹(Ω). Then _Ωdiv𝐰(𝑋(𝑥))𝜓(𝑥)𝑑𝑥=_Γ𝐹^−𝑇𝑛(𝑥)⋅𝐰(𝑋(𝑥))𝜓(𝑥)𝑑𝐀_𝑥−_Ω𝐹⁻¹𝐰(𝑋(𝑥))⋅∇𝜓(𝑥)𝑑𝑥−_Ωdiv𝐹^−𝑇⋅𝐰(𝑋(𝑥))𝜓(𝑥)𝑑𝑥,(2.24) with 𝐰∈𝐻¹(𝑋(Ω)) a vector valued function and 𝜓∈𝐻¹(Ω) a scalar function.

Now, we can multiply (2.23) by a test function 𝜓∈𝐻¹(Ω), integrate in Ω, and apply the usual Green’s formula and (2.24) with 𝑋(𝑥)=𝑋_𝑒(𝑥,𝑡;𝜏), obtaining _Ω𝜙𝑑𝑐𝑑𝜏𝑋^𝑒+(𝑥,𝑡;𝜏),𝜏𝜓(𝑥)𝑑𝑥+_Ω𝐹^𝑒−1(𝑥,𝑡;𝜏)𝐷∇𝑐𝑋_𝑒(𝑥,𝑡;𝜏)⋅∇𝜓(𝑥)𝑑𝑥_Ω=div𝐹^𝑒−𝑇(𝑥,𝑡;𝜏)⋅𝐷∇𝑐𝑋_𝑒(𝑥,𝑡;𝜏),𝜏𝜓(𝑥)𝑑𝑥_Ω((𝑐̃−𝑐)𝑞̃)𝑋_𝑒(𝑥,𝑡;𝜏),𝜏𝜓(𝑥)𝑑𝑥.(2.25)

2.4. The Combined Approximation Scheme

We now present our sequential time-stepping procedure that combines (2.8) and (2.25). Part 𝐽 into 0=𝑡⁰<𝑡¹<⋯<𝑡^𝑁=𝑇, with Δ𝑡^𝑛𝑐=𝑡^𝑛−𝑡^𝑛−1. The analysis is valid for variable time steps, but we drop the superscript from Δ𝑡_𝑐 for convenience. For functions 𝑓 on Ω×𝐽, we write 𝑓^𝑛(𝑥) for 𝑓(𝑥,𝑡^𝑛). As in [3], let us part 𝐽 into pressure time steps 0=𝑡₀<𝑡₁<⋯<𝑡_𝑀=𝑇, with Δ𝑡^𝑚𝑝=𝑡_𝑚−𝑡_𝑚−1. Each pressure step is also a concentration step, that is, for each 𝑚 there exists 𝑛 such that 𝑡_𝑚=𝑡^𝑛, in general, Δ𝑡_𝑝>Δ𝑡_𝑐. We may vary Δ𝑡_𝑝, but except for Δ𝑡^1𝑝 we drop the superscript. For functions 𝑓 on Ω×𝐽, we write 𝑓_𝑚(𝑥) for 𝑓(𝑥,𝑡_𝑚); thus, subscripts refer to pressure steps and superscripts to concentration steps.

If concentration step 𝑡^𝑛 relates to pressure steps by 𝑡_𝑚−1<𝑡^𝑛≤𝑡_𝑚, we require a velocity approximation for (2.25) based on 𝑈_𝑚−1 and earlier values. If 𝑚≥2, take the linear extrapolation of 𝑈_𝑚−1 and 𝑈_𝑚−2 defined by 𝐸𝑈^𝑛=1+𝑡^𝑛−𝑡_𝑚−1𝑡_𝑚−1−𝑡_𝑚−2𝑈^𝑚−1−𝑡^𝑛−𝑡_𝑚−1𝑡_𝑚−1−𝑡_𝑚−2𝑈_𝑚−2;(2.26) if 𝑚=1, set 𝐸𝑈^𝑛=𝑈₀.(2.27) We retain the superscript on Δ𝑡^1𝑝 because 𝐸𝑈^𝑛 is first-order correct in time during the first pressure step and second-order during later steps.

Let 𝑋∶(0,𝑇)→𝑅² be a solution of the ordinary differential equation 𝑑𝑋𝑑𝑡=𝑣̂(𝑋,𝑡),where𝑣̂(𝑋,𝑡)∶=𝐸𝑈𝜙(𝑋,𝑡).(2.28) Then, we can write 𝜙𝑑𝑑𝑡𝑐(𝑋,𝑡)=𝜙𝜕𝜕𝑐𝑡+𝐸𝑈(𝑋,𝑡)⋅∇𝑐.(2.29) Subject to an initial condition 𝑋(𝑡^𝑛+1)=𝑥, we get approximate values of 𝑋 at 𝑡^𝑛 by the Euler method and the second order Runga-Kutta method, respectively, 𝑋𝑋^𝑛𝐸(𝑥)=𝑥−Δ𝑡_𝑐𝑣̂^𝑛+1(𝑥),_𝑛𝑅𝐾(𝑥)=𝑥−Δ𝑡_𝑐𝑣̂^𝑛+(1/2)𝑥−Δ2𝑡^𝑐𝑣̂^𝑛+1(𝑥).(2.30)

Next, assuming they exist, we denote by 𝐹^𝑛𝐸 (resp., by 𝐿_𝐸) the gradient of 𝑋^𝑛𝐸 (resp., of 𝑣(̂𝑥)) with respect to the space variable 𝑥, that is, [13, 14] 𝐹^𝑛𝐸(𝑥)∶=∇𝑋^𝑛𝐸(𝑥)=𝐼(𝑥)−Δ𝑡_𝑐𝐿^𝐸𝑛+1(𝑥).(2.31)

Hypothesis 2.8. For convenience, we assume that (1.1a)–(1.1d) is Ω-periodic ([3]), that is, all functions will be assumed to be spatially Ω-periodic throughout the rest of this paper.

This is physically reasonable, because no-flow conditions (1.1c) are generally treated by reflection, and, in general, interior flow patterns are much more important than boundary effects in reservoir simulation. Thus, the boundary conditions (1.1c) can be dropped.

Lemma 2.9. Under Hypothesis 2.8, if ‖𝑣‖̂_𝐶⁰_(𝑊^1,∞_(Ω))Δ𝑡_𝑐<1/2, we can see that 𝑋^𝑛𝐸Ω=𝑋^𝑛𝑅𝐾Ω=Ω.(2.32)

Lemma 2.10. Under Hypothesis 2.8, if ‖𝑣‖̂_𝐶⁰_(𝑊^1,∞_(Ω))Δ𝑡_𝑐<1/2, we have 𝐹^𝑛𝐸⁻¹(𝑥)=𝐼+Δ𝑡_𝑐𝐿^𝐸𝑛+1(𝑥)+Δ𝑡_𝑐²𝐿^𝐸𝑛+1(𝑥)²+𝑂Δ𝑡_𝑐³.(2.33)

Corollary 2.11. Under the assumptions of Lemma 2.10, forall𝑥∈Ω, we have |||detd𝐹e_𝑛𝐸t𝐹^𝑛𝐸⁻¹(𝑥)=1+Δ𝑡_𝑐div𝑣̂^𝑛+1(𝑥)+𝑂Δ𝑡_𝑐²,⁻¹(𝑥)|||≤1+Δ𝑡_𝑐‖‖𝑣̂^𝑛+1(𝑥)‖‖_{𝐶0(𝑊1,∞(Ω))}+𝑂Δ𝑡_𝑐².(2.34)

The time difference for (2.25) will be combined with a standard finite element in the space variables. For ℎ_𝑐>0, we discrete (2.25) in space on a quasi-uniform mesh 𝒯_ℎ𝑐 of Ω with diameter of element ≤ℎ_𝑐. Let 𝑀_ℎ𝑐⊂𝐻¹(Ω) be a finite element space.

A characteristic discretization of the weak form (2.25) is given by {𝐶^𝑛}^𝑁𝑛=1∈𝑀_ℎ𝑐 such that𝜙𝐶^𝑛+1−𝐶^𝑛∘𝑋^𝑛𝑅𝐾Δ𝑡_𝑐,𝜑+𝐷∇𝐶^𝑛+1+(∇𝐶^𝑛)∘𝑋^𝑛𝐸+Δ𝑡_𝑐2,∇𝜑2𝐷(𝐿^𝑛∇𝐶^𝑛)∘𝑋^𝑛𝐸,∇𝜑+Δ𝑡^𝑐2∇div𝑣^𝑛̂+1⋅𝐷∇𝐶^𝑛∘𝑋^𝑛𝐸,𝜑2𝑞̃^𝑛+1𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛𝐸,𝜑=12𝑞̃^𝑛+1𝐶𝑐̃^𝑛+1+(𝑞̃^𝑛𝑐̃^𝑛)∘𝑋^𝑛𝐸,𝜑,∀𝜑∈𝑀_ℎ𝑐,(2.35a)₀=𝐶⁰,∀𝑥∈Ω,(2.35b)where 𝑔∘𝑋^𝑛𝐸 and 𝑔∘𝑋^𝑛𝑅𝐾 are compositions 𝑔∘𝑋^𝑛𝐸(𝑥)=𝑔𝑋^𝑛𝐸(𝑥),𝑔∘𝑋^𝑛𝑅𝐾(𝑥)=𝑔𝑋^𝑛𝑅𝐾(𝑥),(2.36)𝐶⁰ is an initial approximation of exact solution 𝑐₀(𝑥) into 𝑀_ℎ𝑐, which will be defined in Section 5.

At each pressure time step 𝑡_𝑚, we define 𝑈_𝑚=𝑘𝜇∇𝑃^𝑚𝐶,(2.37)_∗=min(max(𝐶,0),1)(2.38) is the truncation of 𝐶 to [0,1]. Then at 𝑡_𝑚, (2.8) is the following: 𝜇𝑘∇𝑃^𝑚,∇𝜒𝑃=𝑞_𝑚,𝜒,∀𝜒∈𝑊_ℎ𝑝,(2.39a)_𝑚,1=0.(2.39b)The steps of calculation are as follows.

Step 1. 𝐶⁰ known → solve (𝑈₀,𝑃₀) by (2.37), (2.39a) and (2.39b);

Step 2. by (2.35a) and (2.35b) to solve 𝐶¹ → then by (2.35a) and (2.35b) to solve 𝐶²;

Step 3. analogously, {𝐶^𝑗−1}^𝑛1_𝑗=1 known → {𝐶^𝑗}^𝑛1_𝑗=1 such that 𝑡^𝑛1=𝑡₁;

Step 4. then by (2.37), (2.39a) and (2.39b) for (𝑈₁,𝑃₁);

Step 5. calculate the approximations in turn analogously to get the pressure, velocity, and concentration at other time step, respectively.

Throughout the analysis, 𝐾 will denote a generic positive constant, independent of ℎ_𝑐, ℎ_𝑝, Δ𝑡_𝑐, and Δ𝑡_𝑝, but possibly depending on constants in (C). Similarly, 𝜀 will denote a generic small positive constant.

3. Stability for the Concentration Equation

In this section, we derive the stability for the concentration equation. For a given series of functions {𝜑}^𝑁𝑛=0, we define the following norms and seminorm: ‖𝜑‖‖𝜑_𝑙^∞_(𝐿²_)𝑛=max{‖𝜑^𝑛‖;0≤𝑛≤𝑁},‖_𝑙²_(𝐿²₎=Δ𝑡_𝑐𝑁_𝑛=0‖𝜑^𝑛‖²^1/2‖𝜑‖_{𝑙2𝐻1},=Δ𝑡_𝑐𝑁−1_𝑛=0‖‖∇𝜑^𝑛+1+(∇𝜑^𝑛)∘𝑋^𝑛𝐸‖‖^2𝐷^1/2‖𝜑‖_2𝜙,=(𝜙𝜑,𝜑),‖𝜑‖^2𝐷=(𝐷𝜑,𝜑),‖𝜑‖^2𝐷=(𝑞(𝑥,𝑡)𝜑,𝜑).(3.1) In the following sections, we use positive constants as 𝑐₀=𝑐₀‖𝐸𝑈‖_𝐶⁰_(𝐿^∞₎,𝑐₁=𝑐₁‖𝐸𝑈‖_𝑙^∞_(𝑊^(1,∞)₎,𝑐₂=𝑐₂‖𝐸𝑈‖_𝐶⁰_(𝑊^(2,∞)_)∩𝐶¹_(𝐿^∞₎.(3.2) In our analysis, we need some lemmas.

Lemma 3.1. Under the definitions (2.26) and (2.30), for 𝑛=0,…,𝑁 and 𝜑∈𝐿²(Ω), it holds that ‖‖𝜑∘𝑋^𝑛𝑖‖‖²≤1+𝑐₁Δ𝑡_𝑐‖𝜑‖²,𝑖=1,2.(3.3)

Proof. We only need to show a proof in the case 𝑖=1. Let 𝐽₁ be the Jacobian matrix of the transformation 𝑦=𝑋^𝑛𝐸(𝑥)=𝑥−𝑣^𝑛(𝑥)Δ𝑡_𝑐 as 𝐽₁=⎛⎜⎜⎜⎜⎜⎝1−𝜕𝑣^𝑛1(𝑥)𝜕𝑥₁Δ𝑡_𝑐−𝜕𝑣^𝑛1(𝑥)𝜕𝑥₂Δ𝑡_𝑐−𝜕𝑣^𝑛2(𝑥)𝜕𝑥₁Δ𝑡_𝑐1−𝜕𝑣^𝑛2(𝑥)𝜕𝑥₂Δ𝑡_𝑐⎞⎟⎟⎟⎟⎟⎠.(3.4) According to the proof of (3.23) in [4], we have ||||𝜕𝑣^𝑛𝑗(𝑥)𝜕𝑥_𝑖Δ𝑡_𝑐||||≤𝐾ℎ₋₁Δ𝑡_𝑐.(3.5) Since Δ𝑡_𝑐=𝑜(ℎ_𝑐), for ℎ_𝑐 sufficiently small, we see that ||det𝐽₁−1||≤𝑐₁Δ𝑡_𝑐.(3.6)
Following (3.6), we have ‖‖𝜑∘𝑋^𝑛𝐸‖‖²=_Ω𝜑𝑋^𝑛𝐸(𝑥)²𝑑𝑥=_Ω𝜑(𝑦)²det𝐽₁⁻¹𝑑𝑦.(3.7) We complete the proof.

Suppose that {𝐶^𝑛}^𝑁𝑛=0⊂𝐻¹(Ω) and 𝑐̃∈𝐶⁰(𝐿²) be given. We define linear forms 𝒜^ℎ𝑛+1/2 and ℱ^ℎ𝑛+1/2 on 𝑀_ℎ𝑐 for 𝑛=0,…,𝑁−1 by 𝒜^ℎ𝑛+1/2𝐶,𝜑≡𝜙𝐶^𝑛+1−𝐶^𝑛∘𝑋^𝑛𝑅𝐾Δ𝑡_𝑐,𝜑+𝐷∇𝐶^𝑛+1+(∇𝐶^𝑛)∘𝑋^𝑛𝐸+Δ𝑡_𝑐2,∇𝜑2𝐷(𝐿^𝑛∇𝐶^𝑛)∘𝑋^𝑛𝐸,∇𝜑+Δ𝑡^𝑐2∇div𝑣^𝑛̂+1⋅𝐷∇𝐶^𝑛∘𝑋^𝑛𝐸,𝜑2𝑞̃^𝑛+1ℱ𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛𝐸,𝜑_ℎ𝑛+1/2,,𝜑≡12𝑞̃^𝑛+1𝑐̃^𝑛+1+(𝑞̃^𝑛𝑐̃^𝑛)∘𝑋^𝑛𝐸,𝜑,(3.8) where 𝜑∈𝑀_ℎ𝑐.

Lemma 3.2. Let {𝐶^𝑛}^𝑁𝑛=0 be a solution of (2.35a) and (2.35b). Then it holds that 𝒜^ℎ𝑛+1/2𝐶,𝐶^𝑛+1≥𝐷_Δ𝑡𝑐12‖𝐶^𝑛‖^2𝜙+Δ4𝑡^𝑐‖∇𝐶^𝑛‖^2𝐷+Δ𝑡^𝑐𝑞_∗4‖𝐶^𝑛‖²+12Δ𝑡_𝑐‖‖𝐶^𝑛+1−𝐶^𝑛∘𝑋^𝑛𝑅𝐾‖‖^2𝜙+14‖‖∇𝐶^𝑛+1+∇𝐶^𝑛∘𝑋^𝑛𝐸‖‖^2𝐷+14𝑞_∗‖‖𝑞̃^𝑛+1𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛𝐸‖‖²−𝑐¹2‖𝐶^𝑛‖^2𝜙+𝑐¹Δ𝑡_𝑐4‖∇𝐶^𝑛‖^2𝐷−𝑞_∗‖𝐶^𝑛‖²+1+𝑐¹Δ𝑡_𝑐𝑞_∗‖𝑞̃^𝑛𝐶^𝑛‖²+𝑐¹Δ𝑡_𝑐2‖‖∇𝐶^𝑛+1‖‖^2𝐷+‖∇𝐶^𝑛‖^2𝐷+‖‖𝐶^𝑛+1‖‖^2𝐷,(3.9) where 𝐷_Δ𝑡𝑐 is the forward difference operator defined by 𝐷_Δ𝑡𝑐𝜑^𝑛=𝜑^𝑛+1−𝜑^𝑛Δ𝑡_𝑐.(3.10)

Proof. Substituting 𝜑=𝐶^𝑛+1 into (3.8), we have 𝒜^ℎ𝑛+1/2𝐶,𝐶^𝑛+1≡𝜙𝐶^𝑛+1−𝐶^𝑛∘𝑋^𝑛𝑅𝐾Δ𝑡_𝑐,𝐶^𝑛+1+𝐷∇𝐶^𝑛+1+(∇𝐶^𝑛)∘𝑋^𝑛𝐸+Δ𝑡_𝑐2,∇𝐶^𝑛+12𝐷(𝐿^𝑛∇𝐶^𝑛)∘𝑋^𝑛𝐸,∇𝐶^𝑛+1+Δ𝑡^𝑐2∇div𝑣^𝑛̂+1⋅𝐷∇𝐶^𝑛∘𝑋^𝑛𝐸,𝐶^𝑛+12𝑞̃^𝑛+1=𝐼₁𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛𝐸,𝐶^𝑛+1+𝐼₂+𝐼₃+𝐼₄+𝐼₅.(3.11) Lemma 3.1 implies that 𝐼₁=𝜙𝐶^𝑛+1−𝐶^𝑛∘𝑋^𝑛𝑅𝐾Δ𝑡_𝑐,𝐶^𝑛+1+𝐶^𝑛∘𝑋^𝑛𝑅𝐾2+𝐶^𝑛+1−𝐶^𝑛∘𝑋^𝑛𝑅𝐾≥𝐷_Δ𝑡𝑐221‖‖𝐶^𝑛+1‖‖^2𝜙−𝑐2¹‖𝐶^𝑛‖^2𝜙+2Δ1𝑡_𝑐𝐼‖‖𝐶^𝑛+1−𝐶^𝑛∘𝑋^𝑛𝑅𝐾‖‖^2𝜙,(3.12)₂=𝐷∇𝐶^𝑛+1+(∇𝐶^𝑛)∘𝑋^𝑛𝐸2,∇𝐶^𝑛+1+(∇𝐶^𝑛)∘𝑋^𝑛𝐸2+∇𝐶^𝑛+1−(∇𝐶^𝑛)∘𝑋^𝑛𝐸≥𝐷_Δ𝑡𝑐2Δ4𝑡^𝑐‖∇𝐶^𝑛‖^2𝐷−𝑐¹Δ𝑡_𝑐4‖∇𝐶^𝑛‖^2𝐷+14‖‖∇𝐶^𝑛+1+(∇𝐶^𝑛)∘𝑋^𝑛𝐸‖‖^2𝐷.(3.13) Next, by using 𝑐₁Δ𝑡_𝑐<1, we obtain 𝐼₃=Δ𝑡^𝑐2𝐷(𝐿^𝑛≤Δ𝑡_𝑐∇𝐶^𝑛)∘𝑋^𝑛𝐸,∇𝐶^𝑛+12√1+𝑐₁Δ𝑡_𝑐‖𝐿^𝑛∇𝐶^𝑛‖_𝐷‖‖∇𝐶^𝑛+1‖‖_𝐷≤𝑐¹Δ𝑡_𝑐2‖∇𝐶^𝑛‖^2𝐷+‖‖∇𝐶^𝑛+1‖‖^2𝐷.(3.14) Then when 𝐼₃≥0 and 𝐼₃<0, we have 𝐼₃≥−𝑐¹Δ𝑡_𝑐2‖∇𝐶^𝑛‖^2𝐷+‖‖∇𝐶^𝑛+1‖‖^2𝐷.(3.15) Similarly, for 𝐼₄, we obtain the estimate 𝐼₄≥−𝑐¹Δ𝑡_𝑐2‖∇𝐶^𝑛‖^2𝐷+‖‖𝐶^𝑛+1‖‖^2𝐷.(3.16) Analogous computations to term 𝐼₂ give 𝐼₅=𝑞̃^𝑛+1𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛12,𝑞̃^𝑛+1𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛12𝑞̃^𝑛+1+𝑞̃^𝑛+1𝐶^𝑛+1−(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛1≥12𝑞̃^𝑛+14𝑞^∗‖‖𝐶^𝑛+1‖‖²−1𝑞_∗‖‖(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛1‖‖²+41𝑞_∗‖‖𝑞̃^𝑛+1𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛1‖‖²≥𝐷_Δ𝑡𝑐Δ𝑡^𝑐𝑞_∗4‖𝐶^𝑛‖²+𝑞4^∗‖𝐶^𝑛‖²−1+𝑐¹Δ𝑡_𝑐𝑞_∗‖𝑞̃^𝑛𝐶^𝑛‖²+14𝑞_∗‖‖𝑞̃^𝑛+1𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛𝐸‖‖²≥𝐷_Δ𝑡𝑐Δ𝑡^𝑐𝑞_∗4‖𝐶^𝑛‖²+𝑞4^∗‖𝐶^𝑛‖²−2𝑞_∗‖𝑞̃^𝑛𝐶^𝑛‖²+14𝑞_∗‖‖𝑞̃^𝑛+1𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛𝐸‖‖²,(3.17) which completes the proof.

From Lemma 3.1, we have ℱ^ℎ𝑛+1/2,𝐶^𝑛+1≤=121𝑞̃^𝑛+1𝑐̃^𝑛+1+(𝑞̃^𝑛𝑐̃^𝑛)∘𝑋^𝑛𝐸,𝐶^𝑛+12𝜙‖‖𝐶^𝑛+1‖‖^2𝜙+12‖‖𝑞̃^𝑛+1𝑐̃^𝑛+1‖‖²+1+𝑐₁Δ𝑡_𝑐‖𝑞̃^𝑛𝑐̃^𝑛‖².(3.18) Combining (3.18) with (3.8), we get 𝐷_Δ𝑡𝑐12‖𝐶^𝑛‖^2𝜙+Δ𝑡^𝑐4‖∇𝐶^𝑛‖^2𝐷+Δ𝑡^𝑐𝑞_∗4‖𝐶^𝑛‖²+2Δ1𝑡_𝑐‖‖𝐶^𝑛+1−𝐶^𝑛∘𝑋^𝑛𝑅𝐾‖‖^2𝜙+14‖‖∇𝐶^𝑛+1+∇𝐶^𝑛∘𝑋^𝑛𝐸‖‖^2𝐷+14𝑞_∗‖‖𝑞̃^𝑛+1𝐶^𝑛+1+(𝑞̃^𝑛𝐶^𝑛)∘𝑋^𝑛𝐸‖‖²+𝑞_∗‖𝐶^𝑛‖²≤𝑐¹2‖𝐶^𝑛‖^2𝜙+𝑐¹Δ𝑡_𝑐4‖∇𝐶^𝑛‖^2𝐷+2𝑞_∗‖𝑞̃^𝑛𝐶^𝑛‖²+𝑐¹Δ𝑡_𝑐2‖‖∇𝐶^𝑛+1‖‖^2𝐷+‖∇𝐶^𝑛‖^2𝐷+‖‖𝐶^𝑛+1‖‖^2𝐷+21𝜙‖‖𝐶^𝑛+1‖‖^2𝜙+12‖‖𝑞̃^𝑛+1𝑐̃^𝑛+1‖‖²+1+𝑐₁Δ𝑡_𝑐‖𝑞̃^𝑛𝑐̃^𝑛‖²,(3.19) which completes the proof of the stability by virtue of Gronwall’s inequality.

Theorem 3.3 3.3 (stability). Let {𝐶^𝑛}^𝑁𝑛=0 be the solution of (2.35a) and (2.35b) subject to the initial value 𝐶⁰. Then there exists a positive constant 𝑐₁=𝑐₁(‖𝐸𝑈‖_𝑙^∞_(𝑊^1,∞₎) independent of ℎ_𝑐 and Δ𝑡_𝑐 such that ‖𝐶‖_𝑙^∞_(𝐿²₎+√Δ𝑡_𝑐‖∇𝐶‖_𝑙^∞_(𝐿²₎+√Δ𝑡_𝑐‖𝐶‖_𝑙^∞_(𝐿²₎+‖𝐶‖^𝑙_2𝐻1≤𝑐₁‖‖𝐶⁰‖‖+√Δ𝑡_𝑐‖‖∇𝐶⁰‖‖²+‖𝑞̃𝑐̃‖_𝑙²_(𝐿²₎.(3.20)

4. Error Estimate Theorem

Now, we turn to derive an optimal priori error estimate in 𝐿²-norm for the concentration of approximation (2.35a) and (2.35b). In order to state error estimates, we need the following Lagrange interpolation operator [15] Π_ℎ∶𝐶⁰(Ω)→𝑀_ℎ𝑐.

Lemma 4.1. There exists a positive integer 𝑘 such that ‖‖Π_ℎ𝑐−𝑐‖‖_𝑠≤𝐾ℎ^{𝑘+1−𝑠}‖𝑐‖_𝑘+1,𝑠=0,1,∀𝑐∈𝐻^𝑘+1(Ω)∩𝐶⁰Ω.(4.1) Let 𝑒^𝑛=𝐶^𝑛−Π_ℎ𝑐^𝑛 and 𝜂^𝑛=𝑐^𝑛−Π_ℎ𝑐^𝑛. By Lemma 4.1, it holds that ‖𝜂^𝑛‖₁≤𝐾ℎ^𝑘‖𝑐^𝑛‖_𝑘+1,‖‖𝐷_Δ𝑡𝑐𝜂^𝑛‖‖≤𝐾ℎ^𝑘√Δ𝑡_𝑐‖‖‖𝜕𝑐𝜕𝑡‖‖‖_𝐿²_((𝑡^𝑛_,𝑡^𝑛+1_);𝐻^𝑘₎.(4.2)

Let 𝑐∈𝐶¹(𝐿²)∩𝐶⁰(𝐻²), 𝑢∈𝐶⁰(𝐿^∞), and 𝑞̃∈𝐶⁰(𝐿²) be given. Corresponding to 𝒜^ℎ𝑛+1/2 and ℱ^ℎ𝑛+1/2, we introduce linear forms 𝒜^𝑛+1/2 and ℱ^𝑛+1/2 on (𝐻¹(Ω))^ for 𝑛=0,…,𝑁−1 by (2.25) as follows: 𝒜^𝑛+1/2𝑐,𝜑≡+𝜙d𝑑𝑑iv𝑐𝑡𝐹^𝑛+(1/2)∘_{𝑒𝑛+(1/2)}𝑋^{𝑒𝑛+(1/2)},𝜑+𝐹^{𝑒𝑛+(1/2)}⁻¹𝐷∇𝑐^𝑛+(1/2)∘𝑋^{𝑒𝑛+(1/2)},∇𝜑+𝑞̃_𝑛+(1/2)^−𝑇𝐷∇𝑐^𝑛+(1/2)∘𝑋^{𝑒𝑛+(1/2)},𝜑ℱ𝑐^𝑛+(1/2)∘𝑋^{𝑒𝑛+(1/2)},𝜑_𝑛+1/2,,𝜑≡𝑞̃^𝑛+(1/2)𝑐̃^𝑛+(1/2)∘𝑋^{𝑒𝑛+(1/2)},𝜑,∀𝜑∈𝐻¹(Ω).(4.3) If 𝑐 is the solution of (1.1a)–(1.1d), we have for 𝑛=0,…,𝑁−1𝒜^𝑛+(1/2)𝑐,𝜑=ℱ^𝑛+(1/2),𝜑,∀𝜑∈𝐻¹(Ω).(4.4)

We decompose 𝑒 as 𝒜^{ℎ𝑛+(1/2)}𝑒=ℱ^{ℎ𝑛+(1/2)}−ℱ^𝑛+(1/2)+𝒜^𝑛+(1/2)−𝒜^{ℎ𝑛+(1/2)}𝑐+𝒜^{ℎ𝑛+(1/2)}𝜂.(4.5) In order to estimate the terms of the right-hand side of (4.5), we need the following lemmas.

Lemma 4.2 (see [13]). Assume the above Hypotheses hold, and that the coefficients of the problem (1.1a)–(1.1d) satisfy 𝑣∈𝐶⁰(𝑊^3,∞(Ω))∩𝐶¹(𝑊^2,∞(Ω)), 𝑞∈𝑊^2,∞(Ω), and ‖𝑣‖̂_𝐶⁰_(𝑊^1,∞_(Ω))Δ𝑡_𝑐<1/2. Let the solution of (4.4) satisfy 𝑐∈𝑍³, ∇𝑐∈𝑍³. Then, for each 𝑛=0,1,…,𝑁−1, there exists function 𝜉^{𝐴𝑛+(1/2)}∶Ω→𝑅, such that 𝒜^𝑛+(1/2)−𝒜^{ℎ𝑛+(1/2)}𝑐,𝜑=𝜉^{𝐴𝑛+(1/2)},𝜑,𝜑∈𝐻¹(Ω).(4.6) Moreover, 𝜉^{𝐴𝑛+(1/2)}∈𝐿²(Ω) and the following estimate holds: ‖‖𝜉^{𝐴𝑛+(1/2)}‖‖₀≤𝑐̃₁Δ𝑡_𝑐²‖𝑐‖_𝑍³+‖∇𝑐‖_𝑍³,(4.7) where 𝑐̃ denotes a constant independent of Δ𝑡_𝑐.

Lemma 4.3 (see [13]). Assume the above Hypotheses hold, and that the coefficients of the problem (1.1a)–(1.1d) satisfy 𝑣∈𝐶⁰(𝑊^2,∞(Ω))∩𝐶¹(𝑊^1,∞(Ω)), 𝑞∈𝑊^2,∞(Ω), and ‖𝑣‖̂_𝐶⁰_(𝑊^1,∞_(Ω))Δ𝑡_𝑐<1/2. Let the solution of (4.4) satisfy 𝑐∈𝑍³, ∇𝑐∈𝑍³. Then, for each 𝑛=0,1,…,𝑁−1, there exists function 𝜉^{𝑓𝑛+(1/2)}∶Ω→𝑅, such that 𝐹^𝑛+(1/2)−ℱ^{ℎ𝑛+(1/2)},𝜑=𝜉^{𝑓𝑛+(1/2)},𝜑,𝜑∈𝐻¹(Ω).(4.8) Moreover, 𝜉^{𝑓𝑛+(1/2)}∈𝐿²(Ω) and the following estimate holds: ‖‖𝜉^{𝑓𝑛+(1/2)}‖‖₀≤𝑐̃₁Δ𝑡_𝑐²‖𝑞̃𝑐̃‖_𝑍²,(4.9) where 𝑐̃₁ denotes a constant independent of Δ𝑡_𝑐.

Lemma 4.4. Let 𝜂=𝑐−Π_ℎ𝑐. There exists 𝒜^ℎ𝑛+1/2𝜂,𝑒^𝑛+1≤81‖‖∇𝑒^𝑛+1+∇𝑒^𝑛∘𝑋^𝑛𝐸‖‖²+12𝐷∇𝜂^𝑛+1+𝑐₁,∇𝑒^𝑛+1−(𝐷∇𝜂^𝑛,∇𝑒^𝑛)ℎ^2𝑘Δ1𝑡_𝑐‖‖‖𝜕𝑐𝜕𝑡‖‖‖^2𝐿2_𝑡^𝑛_,𝑡^𝑛+1_;𝐻^𝑘_+‖𝑐‖^2𝐶_{0𝐻𝑘+1}+ℎ_2𝑘+2‖𝑐‖^2𝐶_{0𝐻𝑘+1}+Δ𝑡_𝑐‖∇𝑒^𝑛‖²+‖‖∇𝑒^𝑛+1‖‖².(4.10)

Proof. From the definition of 𝒜^ℎ𝑛+1/2, we have 𝒜^ℎ𝑛+1/2𝜂,𝑒^𝑛+1≡𝜙𝜂^𝑛+1−𝜂^𝑛∘𝑋^𝑛𝑅𝐾Δ𝑡_𝑐,𝑒^𝑛+1+𝐷∇𝜂^𝑛+1+(∇𝜂^𝑛)∘𝑋^𝑛𝐸+Δ𝑡_𝑐2,∇𝑒^𝑛+12𝐷(𝐿^𝑛∇𝜂^𝑛)∘𝑋^𝑛𝐸,∇𝑒^𝑛+1+Δ𝑡^𝑐2(∇div𝑣^𝑛+1⋅𝐷∇𝜂^𝑛)∘𝑋^𝑛𝐸,𝑒^𝑛+12𝑞̃^𝑛+1=𝐼₁𝜂^𝑛+1+(𝑞̃^𝑛𝜂^𝑛)∘𝑋^𝑛𝐸,𝑒^𝑛+1+𝐼₂+𝐼₃+𝐼₄+𝐼₅.(4.11) Since it holds that ‖‖‖𝜙𝜂^𝑛+1−𝜂^𝑛∘𝑋^𝑛𝑅𝐾Δ𝑡_𝑐‖‖‖≤𝑐₁√Δ𝑡_𝑐‖‖‖𝜕𝜕𝜂𝑡‖‖‖_𝐿²_(𝑡^𝑛_,𝑡^𝑛+1_;𝐿²₎+‖𝜂‖_𝐿²_{(𝑡𝑛,𝑡𝑛+1;𝐻}¹₎,(4.12) we have 𝐼₁≤1+𝑐¹Δ𝑡_𝑐1+2𝑐₁²𝐾²ℎ^2𝑘2Δ𝑡_𝑐‖‖‖𝜕𝜕𝑐𝑡‖‖‖^2𝐿2_𝑡^𝑛_,𝑡^𝑛+1_;𝐻^𝑘_+‖𝑐‖^2𝐿_{2𝑡𝑛,𝑡𝑛+1;𝐻𝑘+1}+21‖‖𝑒^𝑛+1‖‖².(4.13) To estimate 𝐼₂, we divide it into three parts as 𝐼₂=12𝐷∇𝜂^𝑛+1,∇𝑒^𝑛+1−(𝐷∇𝜂^𝑛,∇𝑒^𝑛)+𝐼₂₁+𝐼₂₂,(4.14) where 𝐼₂₁=12(𝐷∇𝜂^𝑛𝐼)∘𝑋^𝑛𝐸,∇𝑒^𝑛+1+∇𝑒^𝑛∘𝑋^𝑛𝐸₂₂,=12(𝐷∇𝜂^𝑛,∇𝑒^𝑛)−(𝐷∇𝜂^𝑛)∘𝑋^𝑛𝐸,∇𝑒^𝑛∘𝑋^𝑛𝐸.(4.15) The term 𝐼₂₁ is estimated as 𝐼₂₁≤18‖‖∇𝑒^𝑛+1+∇𝑒^𝑛∘𝑋^𝑛𝐸‖‖²+2‖‖(𝐷∇𝜂^𝑛)∘𝑋^𝑛𝐸‖‖²≤18‖‖∇𝑒^𝑛+1+∇𝑒^𝑛∘𝑋^𝑛𝐸‖‖²+1+𝑐₁Δ𝑡_𝑐𝐾𝑑^∗²ℎ^2𝑘‖𝑐^𝑛‖^2𝐻_𝑘+1.(4.16) Using the transformation 𝑦=𝑋^𝑛𝐸(𝑥), we have 𝐼₂₂=12_Ω𝐷∇𝜂^𝑛⋅∇𝑒^𝑛1−det𝜕𝜕𝑥𝑦⁻¹≤𝑐𝑑𝑦𝑑𝑥₁𝑑^∗²Δ𝑡_𝑐4𝐾²ℎ^2𝑘‖𝑐^𝑛‖^2𝐻_𝑘+1+‖∇𝑒^𝑛‖².(4.17) It is clear that 𝐼₃≤1+𝑐¹Δ𝑡_𝑐𝑐²¹𝑑^∗²Δ𝑡_𝑐𝐾²ℎ^2𝑘4‖𝑐^𝑛‖^2𝐻_𝑘+1+14‖‖∇𝑒^𝑛+1𝐼‖‖²,(4.18)₄≤1+𝑐¹Δ𝑡_𝑐𝑐²¹𝑑^∗²Δ𝑡_𝑐𝐾²ℎ^2𝑘4‖𝑐^𝑛‖^2𝐻_𝑘+1+14‖‖∇𝑒^𝑛+1‖‖².(4.19)
Term 𝐼₅ can be divided it into three parts like term 𝐼₂ as 𝐼₅=𝐼₅₁+𝐼₅₂+𝐼₅₃,(4.20) where 𝐼₅₁=12𝑞̃^𝑛+1𝜂^𝑛+1,𝑞̃^𝑛+1𝑒^𝑛+1−√𝐼₅₂𝑞̃^𝑛𝜂^𝑛,√𝑞̃^𝑛𝑒^𝑛,=12√𝐼𝑞̃^𝑛𝜂^𝑛∘𝑋^𝑛𝐸,√𝑞̃^𝑛∘𝑋^𝑛𝐸𝑒^𝑛+1+𝑒^𝑛∘𝑋^𝑛𝐸₅₃,=12√𝑞̃^𝑛𝜂^𝑛,√𝑞̃^𝑛𝑒^𝑛−√𝑞̃^𝑛𝜂^𝑛∘𝑋^𝑛𝐸,√𝑞̃^𝑛𝑒^𝑛∘𝑋^𝑛𝐸.(4.21) Thus, the same kind of computations used for 𝐼₂ leads to the following estimate: 𝐼₅≤𝐼₅₁+𝑐¹𝑐₃Δ𝑡_𝑐4+1+𝑐¹Δ𝑡_𝑐𝑐₃𝐾²ℎ^2𝑘‖𝑐^𝑛‖^2𝐻_𝑘+1+2𝑐¹+𝑐²¹Δ𝑡_𝑐𝑐₃Δ𝑡_𝑐8‖‖√𝑞̃^𝑛𝜂^𝑛‖‖²+‖‖‖𝑞̃^𝑛+1𝜂^𝑛+1‖‖‖²+18‖‖‖𝑞̃^𝑛+1𝜂^𝑛+1+√𝑞̃^𝑛𝜂^𝑛∘𝑋^𝑛𝐸‖‖‖².(4.22) Gathering (4.11)–(4.19) together, we complete the proof.

We now turn to estimate 𝑒^𝑛. From the definition of 𝒜^{ℎ𝑛+(1/2)} and 𝒜^𝑛+(1/2), we see 𝒜^{ℎ𝑛+(1/2)}𝑒,𝑒^𝑛+1=+𝒜ℱ^{ℎ𝑛+(1/2)}𝜂,𝑒^𝑛+1_{ℎ𝑛+(1/2)}+𝒜^𝑛+(1/2)−𝒜^{ℎ𝑛+(1/2)}𝑐,𝑒^𝑛+1−ℱ^𝑛+(1/2),𝑒^𝑛+1.(4.23) From Lemmas 4.2, 4.3 and 4.4, we obtain 𝐷_Δ𝑡𝑐12‖‖𝑒^𝑛+1‖‖^2𝜙+Δ𝑡^𝑐4‖∇𝑒^𝑛‖^2𝐷+Δ𝑡^𝑐4‖𝑒^𝑛‖^2𝑞+2Δ1𝑡_𝑐‖‖𝑒^𝑛+1−𝑒^𝑛∘𝑋^𝑛𝑅𝐾‖‖²+18‖‖∇𝑒^𝑛+1+∇𝑒^𝑛∘𝑋^𝑛𝐸‖‖^2𝐷+14‖‖𝑞̃^𝑛+1𝑒^𝑛+1+𝑞̃^𝑛𝑒^𝑛∘𝑋^𝑛𝐸‖‖²+≤1𝑐₁‖‖𝑒^𝑛+1‖‖²+Δ𝑡_𝑐‖‖∇𝑒^𝑛+1‖‖^2𝐷+‖𝑒^𝑛‖²+Δ𝑡_𝑐‖∇𝑒^𝑛‖^2𝐷+Δ𝑡_𝑐‖𝑒^𝑛‖^2𝑞2𝐷∇𝜂^𝑛+1+𝑐₂,∇𝑒^𝑛+1−(𝐷∇𝜂^𝑛,∇𝑒^𝑛)ℎ^2𝑘Δ1𝑡_𝑐‖‖‖𝜕𝑐𝜕𝑡‖‖‖^2𝐿2_𝑡^𝑛_,𝑡^𝑛+1_;𝐻^𝑘_+‖𝑐‖^2𝐶_{0𝐻𝑘+1}+ℎ^2𝑘+2‖𝑐‖^2𝐶_{0𝐻𝑘+1}+Δ𝑡_𝑐⁴‖‖||∇𝑐||‖‖²³+‖|𝑐|‖²³.(4.24) Summing up the above equation about time from 𝑡=0 to 𝑡^𝑛, we get 𝑎_𝑛+1+Δ𝑡_𝑐𝑛_𝑗=0𝑏_𝑗≤𝑐₁Δ𝑡_𝑐𝑛_𝑗=0𝑎_𝑗+Δ𝑡^𝑐2𝐷∇𝜂^𝑛+1+𝑐₂,∇𝑒^𝑛+1−𝐷∇𝜂⁰,∇𝑒⁰ℎ^2𝑘Δ1𝑡_𝑐‖‖‖𝜕𝑐𝜕𝑡‖‖‖^2𝐿2_𝑡^𝑛_,𝑡^𝑛+1_;𝐻^𝑘_+‖𝑐‖^2𝐶_{0𝐻𝑘+1}+ℎ^2𝑘+2‖𝑐‖^2𝐶_{0𝐻𝑘+1}+Δ𝑡_𝑐⁴‖‖||∇𝑐||‖‖²³+‖|𝑐|‖²³,(4.25) where 𝑎_𝑗=12‖‖𝑒^𝑗‖‖^2𝜙+Δ𝑡^𝑐4‖‖∇𝑒^𝑗‖‖^2𝐷+Δ𝑡^𝑐4‖‖𝑒^𝑗𝑏_𝑗‖‖^2𝑞,=12Δ𝑡_𝑐‖‖𝑒^𝑗+1−𝑒^𝑗∘𝑋^𝑗𝑅𝐾‖‖²+18‖‖∇𝑒^𝑗+1+∇𝑒^𝑗∘𝑋^𝑗𝐸‖‖^2𝐷+14‖‖𝑞̃^𝑗+1𝑒^𝑗+1+𝑞̃^𝑗𝑒^𝑗∘𝑋^𝑗𝐸‖‖².(4.26) Using the estimates Δ𝑡_𝑐2𝐷∇𝜂^𝑗,∇𝑒^𝑗≤Δ𝑡^𝑐8‖‖∇𝑒^𝑗‖‖²+𝑐₁Δ𝑡_𝑐ℎ^2𝑘‖‖𝑐^𝑗‖‖^2𝑘+1(4.27) for 𝑗=0 and 𝑛+1, and Gronwall’s inequality, we derive the following error estimate.

Theorem 4.5 (error estimate). Let {𝐶^𝑛}^𝑁𝑛=0 be the solution of (2.35a) and (2.35b) subject to the initial value 𝐶⁰. Then there exists a positive constant 𝑐₁=𝑐₁(‖𝐸𝑈‖_𝐶⁰_(𝑊^1,∞₎) independent of ℎ_𝑐 and Δ𝑡_𝑐 such that ‖𝑐−𝐶‖_𝑙^∞_(𝐿²₎+√Δ𝑡_𝑐‖∇(𝑐−𝐶)‖_𝑙^∞_(𝐿²₎+√Δ𝑡_𝑐‖𝑐−𝐶‖_𝑙^∞_(𝐿²₎+‖𝑐−𝐶‖^𝑙_2𝐻1≤𝑐₂ℎ^𝑘‖‖‖𝜕𝜕𝑐𝑡‖‖‖_𝐿²_(ℎ^𝑘₎+‖𝑐‖_𝐶⁰_(𝐻^𝑘+1₎+Δ𝑡^𝑐²‖‖||Δ𝑐||‖‖²³+‖|𝑐|‖₃.(4.28)

5. Conclusions

We have presented an approximation scheme for incompressible miscible displacement in porous media. This scheme was constructed by two methods. This paper is our sequential research work. Cubic Hermite finite element method for the pressure equation was used to ensure the higher regularity of the approximation velocity 𝑈. A second order characteristic finite element method was presented to handle the material derivative term of the concentration equation. We analyzed the stability of the approximation scheme and derived the optimal-order 𝐿²-norm error estimates for the scalar concentration. They are of second order accuracy in time increment, symmetric, and unconditionally stable.

In the paper, the matrix 𝐹_𝑒 of the gradient of the characteristic line 𝑋_𝑒 and the inverse matrix 𝐹^𝑒−1 were used to approximate the diffusion term. The properties of these matrices were very important and derived by the complicated theoretical analysis. This paper is the first step of our sequential research work. At current stage, we consider a simple case only for the theoretical aim, which was not related with actual petroleum applications. So, we consider the diffusion coefficient independent of the velocity 𝑢 only in this paper. We will consider the more actual model in petroleum applications later.

Acknowledgments

This paper supported by the Major State Basic Research Development Program of China (Grant no. G19990328), the Scientific Research Award Fund for Excellent Middle-Aged and Young Scientists of Shandong Province (nos. BS2009NJ003 and BS2009HZ015), and the NSF of Shandong Province (no. ZR2010AQ019), China.