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

Authors: Tongjun Sun; Keying Ma
Journal: International Journal of Mathematics and Mathematical Sciences
 (2012)

---

## 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 optimalL2-norm error estimates are derived for the scalar
concentration.


---

## Body

## 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)-∇⋅(kμ∇p)≡∇⋅u=q,x∈Ω,t∈J,(1.1b)ϕ∂c∂t+u⋅∇c-∇⋅(D∇c)=(c~-c)q~,x∈Ω,t∈J,(1.1c)u⋅ν=(D(x)∇c)⋅ν=0,x∈∂Ω,t∈J,(1.1d)c(x,0)=c0(x),x∈Ω,where Ω is a bounded domain in ℛ2, J=(0,T], and q~=max{q,0} is nonzero at injection wells only. The variables in (1.1a)–(1.1d) are the pressure p(x,t) in the fluid mixture, the Darcy velocity u=(u1,u2)′, and the relative concentration c(x,t) of the injected fluid. The ν is the unit outward normal vector on boundary ∂Ω.The coefficients and data in (1.1a)–(1.1d) are k(x), the permeability of the porous media; μ(x), the viscosity of the fluid mixture; q(x,t), representing flow rates at wells; ϕ(x), the porosity of the rock; c~(x,t), the injected concentration at injection wells (q>0) and the resident concentration at production wells (q<0); q~=max(q,0). Here, for the diffusion coefficient, we consider a dispersion-free case [5, 6] as follows:
(1.2)D=ϕdmI,
where dm is the molecular diffusivity and I is a 2×2 identity matrix. Furthermore, a compatibility condition ∫Ω‍q(x,t)dx=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 L2-norm error estimates were derived.It should be pointed out that the works mentioned above all gave one order accuracy in time incrementΔt. 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 Δt, symmetric, and unconditionally stable. Optimal error estimates were proved in the framework of L2-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 l∞(L2) error estimates of second order in time increment Δt 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 L2-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:
(2.1)Wpk=W~pk(Ω)={ψ:∂αψ∂xα∈Lp(Ω)for|α|≤k,periodic}
be the periodic Sobolev space on Ω with the usual norm. If p=2, we write
(2.2)Hk=H~k(Ω)=W~2k(Ω)
with norm
(2.3)‖ψ‖k=‖ψ‖H~k(Ω),‖ψ‖=‖ψ‖0=‖ψ‖L2(Ω).Moreover, we adopt some notations for the functional spaces involved, which were used in [12–14]. For a Banach space X and a positive integer m, spaces Cm([0,T],X) and Hm((0,T),X) are abbreviated as Cm(X) and Hm(X), respectively, and endowed with norms
(2.4)‖φ‖Cm(X):=max⁡t∈[0,T]{max⁡j=0,…,m‖φ(j)(t)‖X},‖φ‖Hm(X):=(∫0T‍∑j=0m‍‖φ(j)(t)‖X2dt)1/2,
where φ(j) denotes the jth derivative of φ with respect to time. The Banach space Zm is defined by
(2.5)Zm={f∈Cj(Hm-j(Ω));j=0,…,m},
equipped with the norm
(2.6)‖φ‖Zm:=max⁡{‖φ‖Cj(Hm-j);0≤j≤m}.We also require the following assumptions on the coefficients in (1.1) [3]. Let a*,a*,ϕ*,ϕ*, and K* be positive constants such that
(C)0<a*≤k(x)μ(x)≤a*,0<ϕ*≤ϕ(x)≤ϕ*,|q~(x,t)|+|∂q~∂t(x,t)|≤K*.
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)(kμ∇p,∇χ)=(q(t),χ),χ∈H1(Ω),t∈J.Forhp>0, we discretize (2.7) in space on a quasiuniform triangle mesh 𝒯hp of Ω with diameter of element ≤hp. Let Whp⊂L2(Ω) be a cubic Hermite finite element space for this mesh. The finite element method for the pressure, given at a time t∈J, consists of P∈Whp such that
(2.8)(kμ∇P,∇χ)=(q(t),χ),∀χ∈Whp.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 p (P, resp.) are ensured obviously [15].By the theory of the cubic Hermit element, the finite element spaceWhp possess the following approximation property [15, 16]:
(2.9)inf⁡χ∈Whp‖g-χ‖s≤Khp4-s‖g‖4,s=0,1,∀g∈H4(Ω),
where K is a positive constant independent of hp.Define Π:C2(Ω)→Whp is the interplant operator. Then, we have the following theorem.Theorem 2.1.
Letp and P be the solutions of (1.1a)–(1.1d) and (2.8), respectively, and assume p∈H4(Ω). Then, there exists a positive constant K independent of hp such that
(2.10)‖∇(p-P)‖≤Khp3‖p‖4.Proof.
By Céa lemma, we have(2.11)a‖∇(p-P)‖2≤(kμ∇(p-P),∇(p-P))=(kμ∇(p-P),∇(p-Πp))+(kμ∇(p-P),∇(Πp-P))=(kμ∇(p-P),∇(p-Πp))≤a*‖∇(p-P)‖‖∇(p-ΠP)‖.
Then by (2.9), we can derive
(2.12)‖∇(p-P)‖≤Khp3‖p‖4.


### 2.3. A Second Order Characteristic Method for the Concentration

Now, we define the characteristic lines associated with vector fieldu and recall some classical properties satisfied by them. Thus, for given (x,t)∈Ω¯×[0,T], the characteristic line through (x,t) is the vector function Xe(x,t;·) solving the following initial value problem:
(2.13)∂Xe∂τ(x,t;τ)=v(Xe(x,t;τ),τ),Xe(x,t;τ)=x,
where v(Xe(x,t;τ),τ):=u(Xe(x,t;τ),τ)/ϕ.Next, assuming they exist, we denote byFe (by L, resp.) the gradient of Xe (of v, resp.) with respect to the space variable x, that is,
(2.14)(Fe)rs(x,t;τ):=∂(Xe)r∂xs(x,t;τ),Lrs(x,t):=∂vr∂xs(x,t).
We adopted some propositions and lemmas from [13].Proposition 2.2.
Ifv∈C0(Cn(Ω¯)) for n≥1 an integer, then Xe∈C0(Ω¯×[0,T]×[0,T]) and it is Cn with respect to the x variable.In order to compute second order approximations of matricesFe and Fe-1, we need the following equations:
(2.15)∂Fe∂τ(x,t;τ)=L(Xe(x,t;τ),τ)Fe(x,t;τ),∂2Fe∂τ2(x,t;τ)=∇(∂v∂t+Lv)(Xe(x,t;τ),τ)Fe(x,t;τ).Proposition 2.3.
Ifv∈C0(C1(Ω¯)), then
(2.16)‖Fe(x,t;τ)‖≤e‖v‖C0(C1(Ω¯))|τ-t|,∀x∈Ω,t,τ∈[0,T].Proposition 2.4.
Ifv∈C0(C2(Ω¯))∩C1(C1(Ω¯)), then Fe satisfies the Taylor expansions as
(2.17)Fe(x,t;s)=I+(s-t)L(x,t)+∫st‍(τ-s)∇(∂v∂t+Lv)(Xe(x,t;τ),τ)Fe(x,t;τ)dτ,
and its inverse, Fe-1, satisfies the Liouville theorem as
(2.18)Fe-1(x,t;s)=I+(s-t)L(Xe(x,t;s),s)-∫st‍(τ-t)∇(∂v∂t+Lv)(Xe(x,t;τ),τ)Fe(Xe(x,t;s),s;τ)dτ.By using Liouville’s theorem and the chain rule, we obtain(2.19)∂∂τdet⁡Fe-1(x,t;τ)=-det⁡Fe-1(x,t;τ)div⁡v(Xe(x,t;τ),τ),∂2∂τ2det⁡Fe-1(x,t;τ)=-det⁡Fe-1(x,t;τ)((div⁡v)2(Xe(x,t;τ),τ)∂2∂τ2det⁡Fe-1(x,t;τ)×div⁡(∂v∂t+Lv)(Xe(x,t;τ),τ)∂2∂τ2det⁡Fe-1(x,t;τ)-(L⋅LT)(Xe(x,t;τ),τ)).Proposition 2.5.
Ifv∈C0(C1(Ω¯)), then
(2.20)det⁡Fe(x,t;τ)≤e‖v‖C0(C1(Ω¯))|τ-t|,∀x∈Ω,t,τ∈[0,T].Proposition 2.6.
Ifv∈C0(C2(Ω¯))∩C1(C1(Ω¯)), then detFe-1 satisfies
(2.21)det⁡Fe-1(x,t;s)=1-(s-t)div⁡v(x,t)+∫st‍(τ-s)∂2∂τ2(det⁡Fe-1)(x,t;τ)dτ.Variational Formulation
From the definition of the characteristic curves and by using the chain rule, it follows that(2.22)ϕdcdτ(Xe(x,t;τ),τ)=ϕ∂c∂t(Xe(x,t;τ),τ)+u(Xe(x,t;τ),τ)⋅∇c(Xe(x,t;τ),τ).
By writing (1.1b) at point Xe(x,t;τ) and time τ and using (2.22), we have
(2.23)ϕdcdτ(Xe(x,t;τ),τ)-∇⋅(D∇c)(Xe(x,t;τ),τ)=((c~-c)q~)(Xe(x,t;τ),τ).
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.
LetX:Ω¯→X(Ω)¯, X∈C2(Ω¯) be an invertible vector valued function. Let F=∇X and assume that F-1∈C1(Ω¯). Then
(2.24)∫Ω‍div⁡w(X(x))ψ(x)dx=∫Γ‍F-Tn(x)⋅w(X(x))ψ(x)dAx-∫Ω‍F-1w(X(x))⋅∇ψ(x)dx-∫Ω‍div⁡F-T⋅w(X(x))ψ(x)dx,
with w∈H1(X(Ω)) a vector valued function and ψ∈H1(Ω) a scalar function.Now, we can multiply (2.23) by a test function ψ∈H1(Ω), integrate in Ω, and apply the usual Green’s formula and (2.24) with X(x)=Xe(x,t;τ), obtaining
(2.25)∫Ω‍ϕdcdτ(Xe(x,t;τ),τ)ψ(x)dx+∫Ω‍Fe-1(x,t;τ)D∇c(Xe(x,t;τ))⋅∇ψ(x)dxψ(x)+∫Ω‍div⁡Fe-T(x,t;τ)⋅D∇c(Xe(x,t;τ),τ)ψ(x)dx33=∫Ω‍((c~-c)q~)(Xe(x,t;τ),τ)ψ(x)dx.


### 2.4. The Combined Approximation Scheme

We now present our sequential time-stepping procedure that combines (2.8) and (2.25). Part J into 0=t0<t1<⋯<tN=T, with Δtcn=tn-tn-1. The analysis is valid for variable time steps, but we drop the superscript from Δtc for convenience. For functions f on Ω×J, we write fn(x) for f(x,tn). As in [3], let us part J into pressure time steps 0=t0<t1<⋯<tM=T, with Δtpm=tm-tm-1. Each pressure step is also a concentration step, that is, for each m there exists n such that tm=tn, in general, Δtp>Δtc. We may vary Δtp, but except for Δtp1 we drop the superscript. For functions f on Ω×J, we write fm(x) for f(x,tm); thus, subscripts refer to pressure steps and superscripts to concentration steps.If concentration steptn relates to pressure steps by tm-1<tn≤tm, we require a velocity approximation for (2.25) based on Um-1 and earlier values. If m≥2, take the linear extrapolation of Um-1 and Um-2 defined by
(2.26)EUn=(1+tn-tm-1tm-1-tm-2)Um-1-tn-tm-1tm-1-tm-2Um-2;
if m=1, set
(2.27)EUn=U0.
We retain the superscript on Δtp1 because EUn is first-order correct in time during the first pressure step and second-order during later steps.LetX:(0,T)→R2 be a solution of the ordinary differential equation
(2.28)dXdt=v^(X,t),wherev^(X,t):=EU(X,t)ϕ.
Then, we can write
(2.29)ϕddtc(X,t)=ϕ∂c∂t+EU(X,t)⋅∇c.
Subject to an initial condition X(tn+1)=x, we get approximate values of X at tn by the Euler method and the second order Runga-Kutta method, respectively,
(2.30)XEn(x)=x-Δtcv^n+1(x),XRKn(x)=x-Δtcv^n+(1/2)(x-Δtc2v^n+1(x)).Next, assuming they exist, we denote byFEn (resp., by LE) the gradient of XEn (resp., of v^(x)) with respect to the space variable x, that is, [13, 14]
(2.31)FEn(x):=∇XEn(x)=I(x)-ΔtcLEn+1(x).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 Hypothesis2.8, if ∥v^∥C0(W1,∞(Ω))Δtc<1/2, we can see that
(2.32)XEn(Ω¯)=XRKn(Ω¯)=Ω¯.Lemma 2.10.
Under Hypothesis2.8, if ∥v^∥C0(W1,∞(Ω))Δtc<1/2, we have
(2.33)(FEn)-1(x)=I+ΔtcLEn+1(x)+(Δtc)2(LEn+1(x))2+O((Δtc)3).Corollary 2.11.
Under the assumptions of Lemma2.10, forallx∈Ω, we have
(2.34)det⁡(FEn)-1(x)=1+Δtcdiv⁡v^n+1(x)+O((Δtc)2),|det⁡(FEn)-1(x)|≤1+Δtc‖v^n+1(x)‖C0(W1,∞(Ω))+O((Δtc)2).The time difference for (2.25) will be combined with a standard finite element in the space variables. For hc>0, we discrete (2.25) in space on a quasi-uniform mesh 𝒯hc of Ω with diameter of element ≤hc. Let Mhc⊂H1(Ω) be a finite element space.A characteristic discretization of the weak form (2.25) is given by {Cn}n=1N∈Mhc such that(2.35a)(ϕCn+1-Cn∘XRKnΔtc,φ)+(D∇Cn+1+(∇Cn)∘XEn2,∇φ)∇Cn+1+Δtc2(D(Ln∇Cn)∘XEn,∇φ)+Δtc2((∇div⁡v^n⋅D∇Cn)∘XEn,φ)∇Cn+1+12(q~n+1Cn+1+(q~nCn)∘XEn,φ)=12(q~n+1c~n+1+(q~nc~n)∘XEn,φ),∀φ∈Mhc,(2.35b)C0=C~0,∀x∈Ω,where g∘XEn and g∘XRKn are compositions
(2.36)(g∘XEn)(x)=g(XEn(x)),(g∘XRKn)(x)=g(XRKn(x)),C~0 is an initial approximation of exact solution c0(x) into Mhc, which will be defined in Section 5.At each pressure time steptm, we define
(2.37)Um=kμ∇Pm,(2.38)C*=min⁡(max⁡(C,0),1)
is the truncation of C to [0,1]. Then at tm, (2.8) is the following: (2.39a)(kμ∇Pm,∇χ)=(qm,χ),∀χ∈Whp,(2.39b)(Pm,1)=0.The steps of calculation are as follows.Step 1.
C
0 known → solve (U0,P0) by (2.37), (2.39a) and (2.39b);Step 2.
by (2.35a) and (2.35b) to solve C1→ then by (2.35a) and (2.35b) to solve C2;Step 3.
analogously,{Cj-1}j=1n1 known →{Cj}j=1n1 such that tn1=t1;Step 4.
then by (2.37), (2.39a) and (2.39b) for (U1,P1);Step 5.
calculate the approximations in turn analogously to get the pressure, velocity, and concentration at other time step, respectively.Throughout the analysis,K will denote a generic positive constant, independent of hc, hp, Δtc, and Δtp, but possibly depending on constants in (C). Similarly, ɛ will denote a generic small positive constant.


## 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:
(2.1)Wpk=W~pk(Ω)={ψ:∂αψ∂xα∈Lp(Ω)for|α|≤k,periodic}
be the periodic Sobolev space on Ω with the usual norm. If p=2, we write
(2.2)Hk=H~k(Ω)=W~2k(Ω)
with norm
(2.3)‖ψ‖k=‖ψ‖H~k(Ω),‖ψ‖=‖ψ‖0=‖ψ‖L2(Ω).Moreover, we adopt some notations for the functional spaces involved, which were used in [12–14]. For a Banach space X and a positive integer m, spaces Cm([0,T],X) and Hm((0,T),X) are abbreviated as Cm(X) and Hm(X), respectively, and endowed with norms
(2.4)‖φ‖Cm(X):=max⁡t∈[0,T]{max⁡j=0,…,m‖φ(j)(t)‖X},‖φ‖Hm(X):=(∫0T‍∑j=0m‍‖φ(j)(t)‖X2dt)1/2,
where φ(j) denotes the jth derivative of φ with respect to time. The Banach space Zm is defined by
(2.5)Zm={f∈Cj(Hm-j(Ω));j=0,…,m},
equipped with the norm
(2.6)‖φ‖Zm:=max⁡{‖φ‖Cj(Hm-j);0≤j≤m}.We also require the following assumptions on the coefficients in (1.1) [3]. Let a*,a*,ϕ*,ϕ*, and K* be positive constants such that
(C)0<a*≤k(x)μ(x)≤a*,0<ϕ*≤ϕ(x)≤ϕ*,|q~(x,t)|+|∂q~∂t(x,t)|≤K*.
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)(kμ∇p,∇χ)=(q(t),χ),χ∈H1(Ω),t∈J.Forhp>0, we discretize (2.7) in space on a quasiuniform triangle mesh 𝒯hp of Ω with diameter of element ≤hp. Let Whp⊂L2(Ω) be a cubic Hermite finite element space for this mesh. The finite element method for the pressure, given at a time t∈J, consists of P∈Whp such that
(2.8)(kμ∇P,∇χ)=(q(t),χ),∀χ∈Whp.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 p (P, resp.) are ensured obviously [15].By the theory of the cubic Hermit element, the finite element spaceWhp possess the following approximation property [15, 16]:
(2.9)inf⁡χ∈Whp‖g-χ‖s≤Khp4-s‖g‖4,s=0,1,∀g∈H4(Ω),
where K is a positive constant independent of hp.Define Π:C2(Ω)→Whp is the interplant operator. Then, we have the following theorem.Theorem 2.1.
Letp and P be the solutions of (1.1a)–(1.1d) and (2.8), respectively, and assume p∈H4(Ω). Then, there exists a positive constant K independent of hp such that
(2.10)‖∇(p-P)‖≤Khp3‖p‖4.Proof.
By Céa lemma, we have(2.11)a‖∇(p-P)‖2≤(kμ∇(p-P),∇(p-P))=(kμ∇(p-P),∇(p-Πp))+(kμ∇(p-P),∇(Πp-P))=(kμ∇(p-P),∇(p-Πp))≤a*‖∇(p-P)‖‖∇(p-ΠP)‖.
Then by (2.9), we can derive
(2.12)‖∇(p-P)‖≤Khp3‖p‖4.


## 2.3. A Second Order Characteristic Method for the Concentration

Now, we define the characteristic lines associated with vector fieldu and recall some classical properties satisfied by them. Thus, for given (x,t)∈Ω¯×[0,T], the characteristic line through (x,t) is the vector function Xe(x,t;·) solving the following initial value problem:
(2.13)∂Xe∂τ(x,t;τ)=v(Xe(x,t;τ),τ),Xe(x,t;τ)=x,
where v(Xe(x,t;τ),τ):=u(Xe(x,t;τ),τ)/ϕ.Next, assuming they exist, we denote byFe (by L, resp.) the gradient of Xe (of v, resp.) with respect to the space variable x, that is,
(2.14)(Fe)rs(x,t;τ):=∂(Xe)r∂xs(x,t;τ),Lrs(x,t):=∂vr∂xs(x,t).
We adopted some propositions and lemmas from [13].Proposition 2.2.
Ifv∈C0(Cn(Ω¯)) for n≥1 an integer, then Xe∈C0(Ω¯×[0,T]×[0,T]) and it is Cn with respect to the x variable.In order to compute second order approximations of matricesFe and Fe-1, we need the following equations:
(2.15)∂Fe∂τ(x,t;τ)=L(Xe(x,t;τ),τ)Fe(x,t;τ),∂2Fe∂τ2(x,t;τ)=∇(∂v∂t+Lv)(Xe(x,t;τ),τ)Fe(x,t;τ).Proposition 2.3.
Ifv∈C0(C1(Ω¯)), then
(2.16)‖Fe(x,t;τ)‖≤e‖v‖C0(C1(Ω¯))|τ-t|,∀x∈Ω,t,τ∈[0,T].Proposition 2.4.
Ifv∈C0(C2(Ω¯))∩C1(C1(Ω¯)), then Fe satisfies the Taylor expansions as
(2.17)Fe(x,t;s)=I+(s-t)L(x,t)+∫st‍(τ-s)∇(∂v∂t+Lv)(Xe(x,t;τ),τ)Fe(x,t;τ)dτ,
and its inverse, Fe-1, satisfies the Liouville theorem as
(2.18)Fe-1(x,t;s)=I+(s-t)L(Xe(x,t;s),s)-∫st‍(τ-t)∇(∂v∂t+Lv)(Xe(x,t;τ),τ)Fe(Xe(x,t;s),s;τ)dτ.By using Liouville’s theorem and the chain rule, we obtain(2.19)∂∂τdet⁡Fe-1(x,t;τ)=-det⁡Fe-1(x,t;τ)div⁡v(Xe(x,t;τ),τ),∂2∂τ2det⁡Fe-1(x,t;τ)=-det⁡Fe-1(x,t;τ)((div⁡v)2(Xe(x,t;τ),τ)∂2∂τ2det⁡Fe-1(x,t;τ)×div⁡(∂v∂t+Lv)(Xe(x,t;τ),τ)∂2∂τ2det⁡Fe-1(x,t;τ)-(L⋅LT)(Xe(x,t;τ),τ)).Proposition 2.5.
Ifv∈C0(C1(Ω¯)), then
(2.20)det⁡Fe(x,t;τ)≤e‖v‖C0(C1(Ω¯))|τ-t|,∀x∈Ω,t,τ∈[0,T].Proposition 2.6.
Ifv∈C0(C2(Ω¯))∩C1(C1(Ω¯)), then detFe-1 satisfies
(2.21)det⁡Fe-1(x,t;s)=1-(s-t)div⁡v(x,t)+∫st‍(τ-s)∂2∂τ2(det⁡Fe-1)(x,t;τ)dτ.Variational Formulation
From the definition of the characteristic curves and by using the chain rule, it follows that(2.22)ϕdcdτ(Xe(x,t;τ),τ)=ϕ∂c∂t(Xe(x,t;τ),τ)+u(Xe(x,t;τ),τ)⋅∇c(Xe(x,t;τ),τ).
By writing (1.1b) at point Xe(x,t;τ) and time τ and using (2.22), we have
(2.23)ϕdcdτ(Xe(x,t;τ),τ)-∇⋅(D∇c)(Xe(x,t;τ),τ)=((c~-c)q~)(Xe(x,t;τ),τ).
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.
LetX:Ω¯→X(Ω)¯, X∈C2(Ω¯) be an invertible vector valued function. Let F=∇X and assume that F-1∈C1(Ω¯). Then
(2.24)∫Ω‍div⁡w(X(x))ψ(x)dx=∫Γ‍F-Tn(x)⋅w(X(x))ψ(x)dAx-∫Ω‍F-1w(X(x))⋅∇ψ(x)dx-∫Ω‍div⁡F-T⋅w(X(x))ψ(x)dx,
with w∈H1(X(Ω)) a vector valued function and ψ∈H1(Ω) a scalar function.Now, we can multiply (2.23) by a test function ψ∈H1(Ω), integrate in Ω, and apply the usual Green’s formula and (2.24) with X(x)=Xe(x,t;τ), obtaining
(2.25)∫Ω‍ϕdcdτ(Xe(x,t;τ),τ)ψ(x)dx+∫Ω‍Fe-1(x,t;τ)D∇c(Xe(x,t;τ))⋅∇ψ(x)dxψ(x)+∫Ω‍div⁡Fe-T(x,t;τ)⋅D∇c(Xe(x,t;τ),τ)ψ(x)dx33=∫Ω‍((c~-c)q~)(Xe(x,t;τ),τ)ψ(x)dx.


## 2.4. The Combined Approximation Scheme

We now present our sequential time-stepping procedure that combines (2.8) and (2.25). Part J into 0=t0<t1<⋯<tN=T, with Δtcn=tn-tn-1. The analysis is valid for variable time steps, but we drop the superscript from Δtc for convenience. For functions f on Ω×J, we write fn(x) for f(x,tn). As in [3], let us part J into pressure time steps 0=t0<t1<⋯<tM=T, with Δtpm=tm-tm-1. Each pressure step is also a concentration step, that is, for each m there exists n such that tm=tn, in general, Δtp>Δtc. We may vary Δtp, but except for Δtp1 we drop the superscript. For functions f on Ω×J, we write fm(x) for f(x,tm); thus, subscripts refer to pressure steps and superscripts to concentration steps.If concentration steptn relates to pressure steps by tm-1<tn≤tm, we require a velocity approximation for (2.25) based on Um-1 and earlier values. If m≥2, take the linear extrapolation of Um-1 and Um-2 defined by
(2.26)EUn=(1+tn-tm-1tm-1-tm-2)Um-1-tn-tm-1tm-1-tm-2Um-2;
if m=1, set
(2.27)EUn=U0.
We retain the superscript on Δtp1 because EUn is first-order correct in time during the first pressure step and second-order during later steps.LetX:(0,T)→R2 be a solution of the ordinary differential equation
(2.28)dXdt=v^(X,t),wherev^(X,t):=EU(X,t)ϕ.
Then, we can write
(2.29)ϕddtc(X,t)=ϕ∂c∂t+EU(X,t)⋅∇c.
Subject to an initial condition X(tn+1)=x, we get approximate values of X at tn by the Euler method and the second order Runga-Kutta method, respectively,
(2.30)XEn(x)=x-Δtcv^n+1(x),XRKn(x)=x-Δtcv^n+(1/2)(x-Δtc2v^n+1(x)).Next, assuming they exist, we denote byFEn (resp., by LE) the gradient of XEn (resp., of v^(x)) with respect to the space variable x, that is, [13, 14]
(2.31)FEn(x):=∇XEn(x)=I(x)-ΔtcLEn+1(x).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 Hypothesis2.8, if ∥v^∥C0(W1,∞(Ω))Δtc<1/2, we can see that
(2.32)XEn(Ω¯)=XRKn(Ω¯)=Ω¯.Lemma 2.10.
Under Hypothesis2.8, if ∥v^∥C0(W1,∞(Ω))Δtc<1/2, we have
(2.33)(FEn)-1(x)=I+ΔtcLEn+1(x)+(Δtc)2(LEn+1(x))2+O((Δtc)3).Corollary 2.11.
Under the assumptions of Lemma2.10, forallx∈Ω, we have
(2.34)det⁡(FEn)-1(x)=1+Δtcdiv⁡v^n+1(x)+O((Δtc)2),|det⁡(FEn)-1(x)|≤1+Δtc‖v^n+1(x)‖C0(W1,∞(Ω))+O((Δtc)2).The time difference for (2.25) will be combined with a standard finite element in the space variables. For hc>0, we discrete (2.25) in space on a quasi-uniform mesh 𝒯hc of Ω with diameter of element ≤hc. Let Mhc⊂H1(Ω) be a finite element space.A characteristic discretization of the weak form (2.25) is given by {Cn}n=1N∈Mhc such that(2.35a)(ϕCn+1-Cn∘XRKnΔtc,φ)+(D∇Cn+1+(∇Cn)∘XEn2,∇φ)∇Cn+1+Δtc2(D(Ln∇Cn)∘XEn,∇φ)+Δtc2((∇div⁡v^n⋅D∇Cn)∘XEn,φ)∇Cn+1+12(q~n+1Cn+1+(q~nCn)∘XEn,φ)=12(q~n+1c~n+1+(q~nc~n)∘XEn,φ),∀φ∈Mhc,(2.35b)C0=C~0,∀x∈Ω,where g∘XEn and g∘XRKn are compositions
(2.36)(g∘XEn)(x)=g(XEn(x)),(g∘XRKn)(x)=g(XRKn(x)),C~0 is an initial approximation of exact solution c0(x) into Mhc, which will be defined in Section 5.At each pressure time steptm, we define
(2.37)Um=kμ∇Pm,(2.38)C*=min⁡(max⁡(C,0),1)
is the truncation of C to [0,1]. Then at tm, (2.8) is the following: (2.39a)(kμ∇Pm,∇χ)=(qm,χ),∀χ∈Whp,(2.39b)(Pm,1)=0.The steps of calculation are as follows.Step 1.
C
0 known → solve (U0,P0) by (2.37), (2.39a) and (2.39b);Step 2.
by (2.35a) and (2.35b) to solve C1→ then by (2.35a) and (2.35b) to solve C2;Step 3.
analogously,{Cj-1}j=1n1 known →{Cj}j=1n1 such that tn1=t1;Step 4.
then by (2.37), (2.39a) and (2.39b) for (U1,P1);Step 5.
calculate the approximations in turn analogously to get the pressure, velocity, and concentration at other time step, respectively.Throughout the analysis,K will denote a generic positive constant, independent of hc, hp, Δtc, and Δtp, 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{φ}n=0N, we define the following norms and seminorm:
(3.1)‖φ‖l∞(L2)=max⁡{‖φn‖;0≤n≤N},‖φn‖l2(L2)={Δtc∑n=0N‍‖φn‖2}1/2,‖φ‖l2(H1)′={Δtc∑n=0N-1‍‖∇φn+1+(∇φn)∘XEn‖D2}1/2,‖φ‖ϕ2=(ϕφ,φ),‖φ‖D2=(Dφ,φ),‖φ‖D2=(q(x,t)φ,φ).
In the following sections, we use positive constants as
(3.2)c0=c0(‖EU‖C0(L∞)),c1=c1(‖EU‖l∞(W(1,∞))),c2=c2(‖EU‖C0(W(2,∞))∩C1(L∞)).
In our analysis, we need some lemmas.Lemma 3.1.
Under the definitions (2.26) and (2.30), for n=0,…,N and φ∈L2(Ω), it holds that
(3.3)‖φ∘Xin‖2≤(1+c1Δtc)‖φ‖2,i=1,2.Proof.
We only need to show a proof in the casei=1. Let J1 be the Jacobian matrix of the transformation y=XEn(x)=x-vn(x)Δtc as
(3.4)J1=(1-∂v1n(x)∂x1Δtc-∂v1n(x)∂x2Δtc-∂v2n(x)∂x1Δtc1-∂v2n(x)∂x2Δtc).
According to the proof of (3.23) in [4], we have
(3.5)|∂vjn(x)∂xiΔtc|≤Kh-1Δtc.
Since Δtc=o(hc), for hc sufficiently small, we see that
(3.6)|det⁡J1-1|≤c1Δtc.
Following (3.6), we have
(3.7)‖φ∘XEn‖2=∫Ω‍φ(XEn(x))2dx=∫Ω‍φ(y)2(det⁡J1)-1dy.
We complete the proof.Suppose that{Cn}n=0N⊂H1(Ω) and c~∈C0(L2) be given. We define linear forms 𝒜hn+1/2 and ℱhn+1/2 on Mhc for n=0,…,N-1 by
(3.8)〈Ahn+1/2C,φ〉≡(ϕCn+1-Cn∘XRKnΔtc,φ)+(D∇Cn+1+(∇Cn)∘XEn2,∇φ)+Δtc2(D(Ln∇Cn)∘XEn,∇φ)+Δtc2((∇div⁡v^n⋅D∇Cn)∘XEn,φ)+12(q~n+1Cn+1+(q~nCn)∘XEn,φ),〈Fhn+1/2,φ〉≡12(q~n+1c~n+1+(q~nc~n)∘XEn,φ),
where φ∈Mhc.Lemma 3.2.
Let{Cn}n=0N be a solution of (2.35a) and (2.35b). Then it holds that
(3.9)〈Ahn+1/2C,Cn+1〉≥DΔtc(12‖Cn‖ϕ2+Δtc4‖∇Cn‖D2+Δtcq*4‖Cn‖2)+12Δtc‖Cn+1-Cn∘XRKn‖ϕ2+14‖∇Cn+1+∇Cn∘XEn‖D2+14q*‖q~n+1Cn+1+(q~nCn)∘XEn‖2-{c12‖Cn‖ϕ2+c1Δtc4‖∇Cn‖D2-q*‖Cn‖2+1+c1Δtcq*‖q~nCn‖2==+c1Δtc2(‖∇Cn+1‖D2+‖∇Cn‖D2+‖Cn+1‖D2)},
where DΔtc is the forward difference operator defined by
(3.10)DΔtcφn=(φn+1-φn)Δtc.Proof.
Substitutingφ=Cn+1 into (3.8), we have
(3.11)〈Ahn+1/2C,Cn+1〉≡(ϕCn+1-Cn∘XRKnΔtc,Cn+1)+(D∇Cn+1+(∇Cn)∘XEn2,∇Cn+1)+Δtc2(D(Ln∇Cn)∘XEn,∇Cn+1)+Δtc2((∇div⁡v^n⋅D∇Cn)∘XEn,Cn+1)+12(q~n+1Cn+1+(q~nCn)∘XEn,Cn+1)=I1+I2+I3+I4+I5.
Lemma 3.1 implies that
(3.12)I1=(ϕCn+1-Cn∘XRKnΔtc,Cn+1+Cn∘XRKn2+Cn+1-Cn∘XRKn2)≥DΔtc(12‖Cn+1‖ϕ2)-c12‖Cn‖ϕ2+12Δtc‖Cn+1-Cn∘XRKn‖ϕ2,(3.13)I2=(D∇Cn+1+(∇Cn)∘XEn2,∇Cn+1+(∇Cn)∘XEn2+∇Cn+1-(∇Cn)∘XEn2)≥DΔtc(Δtc4‖∇Cn‖D2)-c1Δtc4‖∇Cn‖D2+14‖∇Cn+1+(∇Cn)∘XEn‖D2.
Next, by using c1Δtc<1, we obtain
(3.14)I3=Δtc2(D(Ln∇Cn)∘XEn,∇Cn+1)≤Δtc21+c1Δtc‖Ln∇Cn‖D‖∇Cn+1‖D≤c1Δtc2{‖∇Cn‖D2+‖∇Cn+1‖D2}.
Then when I3≥0 and I3<0, we have
(3.15)I3≥-c1Δtc2{‖∇Cn‖D2+‖∇Cn+1‖D2}.
Similarly, for I4, we obtain the estimate
(3.16)I4≥-c1Δtc2{‖∇Cn‖D2+‖Cn+1‖D2}.
Analogous computations to term I2 give
(3.17)I5=(q~n+1Cn+1+(q~nCn)∘X1n2,q~n+1Cn+1+(q~nCn)∘X1n2q~n+1+q~n+1Cn+1-(q~nCn)∘X1n2q~n+1)≥14(q*‖Cn+1‖2-1q*‖(q~nCn)∘X1n‖2)+14q*‖q~n+1Cn+1+(q~nCn)∘X1n‖2≥DΔtc(Δtcq*4‖Cn‖2)+q*4‖Cn‖2-1+c1Δtcq*‖q~nCn‖2+14q*‖q~n+1Cn+1+(q~nCn)∘XEn‖2≥DΔtc(Δtcq*4‖Cn‖2)+q*4‖Cn‖2-2q*‖q~nCn‖2+14q*‖q~n+1Cn+1+(q~nCn)∘XEn‖2,
which completes the proof.From Lemma3.1, we have
(3.18)〈Fhn+1/2,Cn+1〉=12(q~n+1c~n+1+(q~nc~n)∘XEn,Cn+1)≤12ϕ‖Cn+1‖ϕ2+12{‖q~n+1c~n+1‖2+(1+c1Δtc)‖q~nc~n‖2}.
Combining (3.18) with (3.8), we get
(3.19)DΔtc(12‖Cn‖ϕ2+Δtc4‖∇Cn‖D2+Δtcq*4‖Cn‖2)+12Δtc‖Cn+1-Cn∘XRKn‖ϕ2DΔtc+14‖∇Cn+1+∇Cn∘XEn‖D2+14q*‖q~n+1Cn+1+(q~nCn)∘XEn‖2+q*‖Cn‖2D≤c12‖Cn‖ϕ2+c1Δtc4‖∇Cn‖D2+2q*‖q~nCn‖2DΔtc+c1Δtc2(‖∇Cn+1‖D2+‖∇Cn‖D2+‖Cn+1‖D2)+12ϕ‖Cn+1‖ϕ2DΔtc+12{‖q~n+1c~n+1‖2+(1+c1Δtc)‖q~nc~n‖2},
which completes the proof of the stability by virtue of Gronwall’s inequality.Theorem 3.33.3 (stability).
Let{Cn}n=0N be the solution of (2.35a) and (2.35b) subject to the initial value C0. Then there exists a positive constant c1=c1(∥EU∥l∞(W1,∞)) independent of hc and Δtc such that
(3.20)‖C‖l∞(L2)+Δtc‖∇C‖l∞(L2)+Δtc‖C‖l∞(L2)+‖C‖l2(H1)′1234≤c1{‖C0‖+Δtc‖∇C0‖2+‖q~c~‖l2(L2)}.


## 4. Error Estimate Theorem

Now, we turn to derive an optimal priori error estimate inL2-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] Πh:C0(Ω¯)→Mhc.Lemma 4.1.
There exists a positive integerk such that
(4.1)‖Πhc-c‖s≤Khk+1-s‖c‖k+1,s=0,1,∀c∈Hk+1(Ω)∩C0(Ω¯).
Let en=Cn-Πhcn and ηn=cn-Πhcn. By Lemma 4.1, it holds that
(4.2)‖ηn‖1≤Khk‖cn‖k+1,‖DΔtcηn‖≤KhkΔtc‖∂c∂t‖L2((tn,tn+1);Hk).Letc∈C1(L2)∩C0(H2), u∈C0(L∞), and q~∈C0(L2) be given. Corresponding to 𝒜hn+1/2 and ℱhn+1/2, we introduce linear forms 𝒜n+1/2 and ℱn+1/2 on (H1(Ω))′ for n=0,…,N-1 by (2.25) as follows:
(4.3)〈An+1/2c,φ〉≡(ϕ(dcdt)n+(1/2)∘Xen+(1/2),φ)+((Fen+(1/2))-1(D∇cn+(1/2))∘Xen+(1/2),∇φ)+(div⁡(Fen+(1/2))-T(D∇cn+(1/2))∘Xen+(1/2),φ)+((q~n+(1/2)cn+(1/2))∘Xen+(1/2),φ),〈Fn+1/2,φ〉≡((q~n+(1/2)c~n+(1/2))∘Xen+(1/2),φ),∀φ∈H1(Ω).
If c is the solution of (1.1a)–(1.1d), we have for n=0,…,N-1(4.4)(An+(1/2)c,φ)=(Fn+(1/2),φ),∀φ∈H1(Ω).We decomposee as
(4.5)Ahn+(1/2)e=(Fhn+(1/2)-Fn+(1/2))+(An+(1/2)-Ahn+(1/2))c+Ahn+(1/2)η.
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 v∈C0(W3,∞(Ω))∩C1(W2,∞(Ω)), q∈W2,∞(Ω), and ∥v^∥C0(W1,∞(Ω))Δtc<1/2. Let the solution of (4.4) satisfy c∈Z3, ∇c∈Z3. Then, for each n=0,1,…,N-1, there exists function ξAn+(1/2):Ω→R, such that
(4.6)〈(An+(1/2)-Ahn+(1/2))c,φ〉=(ξAn+(1/2),φ),φ∈H1(Ω).
Moreover, ξAn+(1/2)∈L2(Ω) and the following estimate holds:
(4.7)‖ξAn+(1/2)‖0≤c~1(Δtc)2(‖c‖Z3+‖∇c‖Z3),
where c~ denotes a constant independent of Δtc.Lemma 4.3 (see [13]).
Assume the above Hypotheses hold, and that the coefficients of the problem (1.1a)–(1.1d) satisfy v∈C0(W2,∞(Ω))∩C1(W1,∞(Ω)), q∈W2,∞(Ω), and ∥v^∥C0(W1,∞(Ω))Δtc<1/2. Let the solution of (4.4) satisfy c∈Z3, ∇c∈Z3. Then, for each n=0,1,…,N-1, there exists function ξfn+(1/2):Ω→R, such that
(4.8)〈(Fn+(1/2)-Fhn+(1/2),φ〉=(ξfn+(1/2),φ),φ∈H1(Ω).
Moreover, ξfn+(1/2)∈L2(Ω) and the following estimate holds:
(4.9)‖ξfn+(1/2)‖0≤c~1(Δtc)2‖q~c~‖Z2,
where c~1 denotes a constant independent of Δtc.Lemma 4.4.
Letη=c-Πhc. There exists
(4.10)〈Ahn+1/2η,en+1〉≤18‖∇en+1+∇en∘XEn‖2+12{(D∇ηn+1,∇en+1)-(D∇ηn,∇en)}+c1{h2k(1Δtc‖∂c∂t‖L2(tn,tn+1;Hk)2+‖c‖C0(Hk+1)2)+c1+1+1+h2k+2‖c‖C0(Hk+1)2+Δtc(‖∇en‖2+‖∇en+1‖2)‖∂c∂t‖L2(tn,tn+1;Hk)21Δtc}.Proof.
From the definition of𝒜hn+1/2, we have
(4.11)〈Ahn+1/2η,en+1〉≡(ϕηn+1-ηn∘XRKnΔtc,en+1)+(D∇ηn+1+(∇ηn)∘XEn2,∇en+1)+Δtc2(D(Ln∇ηn)∘XEn,∇en+1)+Δtc2((∇div⁡vn⋅D∇ηn)∘XEn,en+1)+12(q~n+1ηn+1+(q~nηn)∘XEn,en+1)=I1+I2+I3+I4+I5.
Since it holds that
(4.12)‖ϕηn+1-ηn∘XRKnΔtc‖≤c1Δtc{‖∂η∂t‖L2(tn,tn+1;L2)+‖η‖L2(tn,tn+1;H1)},
we have
(4.13)I1≤(1+c1Δtc)(1+2c1)2K2h2k2Δtc(‖∂c∂t‖L2(tn,tn+1;Hk)2+‖c‖L2(tn,tn+1;Hk+1)2)+12‖en+1‖2.
To estimate I2, we divide it into three parts as
(4.14)I2=12[(D∇ηn+1,∇en+1)-(D∇ηn,∇en)]+I21+I22,
where
(4.15)I21=12[((D∇ηn)∘XEn,∇en+1+∇en∘XEn)],I22=12[(D∇ηn,∇en)-((D∇ηn)∘XEn,∇en∘XEn)].
The term I21 is estimated as
(4.16)I21≤18‖∇en+1+∇en∘XEn‖2+2‖(D∇ηn)∘XEn‖2≤18‖∇en+1+∇en∘XEn‖2+(1+c1Δtc)K(d*)2h2k‖cn‖Hk+12.
Using the transformation y=XEn(x), we have
(4.17)I22=12∫ΩD∇ηn⋅∇en{1-det⁡(∂y∂x)-1dy}dx≤c1(d*)2Δtc4(K2h2k‖cn‖Hk+12+‖∇en‖2).
It is clear that
(4.18)I3≤(1+c1Δtc)c12(d*)2ΔtcK2h2k4‖cn‖Hk+12+14‖∇en+1‖2,(4.19)I4≤(1+c1Δtc)c12(d*)2ΔtcK2h2k4‖cn‖Hk+12+14‖∇en+1‖2.
TermI5 can be divided it into three parts like term I2 as
(4.20)I5=I51+I52+I53,
where
(4.21)I51=12[(q~n+1ηn+1,q~n+1en+1)-(q~nηn,q~nen)],I52=12[((q~nηn)∘XEn,q~n∘XEn(en+1+en∘XEn))],I53=12[(q~nηn,q~nen)-((q~nηn)∘XEn,(q~nen)∘XEn)].
Thus, the same kind of computations used for I2 leads to the following estimate:
(4.22)I5≤I51+(c1c3Δtc4+(1+c1Δtc)c3)K2h2k‖cn‖Hk+12+(2c1+c12Δtc)c3Δtc8(‖q~nηn‖2+‖q~n+1ηn+1‖2)+18‖q~n+1ηn+1+(q~nηn)∘XEn‖2.
Gathering (4.11)–(4.19) together, we complete the proof.We now turn to estimateen. From the definition of 𝒜hn+(1/2) and 𝒜n+(1/2), we see
(4.23)〈Ahn+(1/2)e,en+1〉=〈Ahn+(1/2)η,en+1〉+〈(An+(1/2)-Ahn+(1/2))c,en+1〉+〈(Fhn+(1/2)-Fn+(1/2)),en+1〉.
From Lemmas 4.2, 4.3 and 4.4, we obtain
(4.24)DΔtc(12‖en+1‖ϕ2+Δtc4‖∇en‖D2+Δtc4‖en‖q2)+12Δtc‖en+1-en∘XRKn‖2DΔtc+18‖∇en+1+∇en∘XEn‖D2+14‖q~n+1en+1+q~nen∘XEn‖2=≤c1{‖en+1‖2+Δtc‖∇en+1‖D2+‖en‖2+Δtc‖∇en‖D2+Δtc‖en‖q2}DΔtc+12{(D∇ηn+1,∇en+1)-(D∇ηn,∇en)}DΔtc+c2{h2k(1Δtc‖∂c∂t‖L2(tn,tn+1;Hk)2+‖c‖C0(Hk+1)2)+h2k+2‖c‖C0(Hk+1)2DΔtc+c22+c2=+(Δtc)4(‖|∇c|‖32+‖|c|‖32)‖∂c∂t‖L2(tn,tn+1;Hk)2}.
Summing up the above equation about time from t=0 to tn, we get
(4.25)an+1+Δtc∑j=0n‍bj≤c1Δtc∑j=0n‍aj+Δtc2{(D∇ηn+1,∇en+1)-(D∇η0,∇e0)}+c2{h2k(1Δtc‖∂c∂t‖L2(tn,tn+1;Hk)2+‖c‖C0(Hk+1)2)+h2k+2‖c‖C0(Hk+1)2+c2=+c2=+(Δtc)4(‖|∇c|‖32+‖|c|‖32)‖∂c∂t‖L2(tn,tn+1;Hk)2},
where
(4.26)aj=12‖ej‖ϕ2+Δtc4‖∇ej‖D2+Δtc4‖ej‖q2,bj=12Δtc‖ej+1-ej∘XRKj‖2+18‖∇ej+1+∇ej∘XEj‖D2+14‖q~j+1ej+1+q~jej∘XEj‖2.
Using the estimates
(4.27)Δtc2(D∇ηj,∇ej)≤Δtc8‖∇ej‖2+c1Δtch2k‖cj‖k+12
for j=0 and n+1, and Gronwall’s inequality, we derive the following error estimate.Theorem 4.5 (error estimate).
Let{Cn}n=0N be the solution of (2.35a) and (2.35b) subject to the initial value C0. Then there exists a positive constant c1=c1(∥EU∥C0(W1,∞)) independent of hc and Δtc such that
(4.28)‖c-C‖l∞(L2)+Δtc‖∇(c-C)‖l∞(L2)+Δtc‖c-C‖l∞(L2)+‖c-C‖l2(H1)′‖|Δc|‖32≤c2{hk(‖∂c∂t‖L2(hk)+‖c‖C0(Hk+1))+(Δtc)2(‖|Δc|‖32+‖|c|‖3)}.


## 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 velocityU. 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 L2-norm error estimates for the scalar concentration. They are of second order accuracy in time increment, symmetric, and unconditionally stable.In the paper, the matrixFe of the gradient of the characteristic line Xe and the inverse matrix Fe-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 u only in this paper. We will consider the more actual model in petroleum applications later.



---

*Source: 870402-2012-08-30.xml*