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

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


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: where Ω is a bounded domain in R 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 u 1 , u 2 , 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: where d m 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 L 2 -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 convectiondiffusion 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 L 2 -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 ∞ L 2 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 L 2 -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.
International Journal of Mathematics and Mathematical Sciences 3

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: be the periodic Sobolev space on Ω with the usual norm.If p 2, we write 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 C m 0, T , X and H m 0, T , X are abbreviated as C m X and H m X , respectively, and endowed with norms ϕ C m X : max t∈ 0,T max j 0,...,m where ϕ j denotes the jth derivative of ϕ with respect to time.The Banach space Z m is defined by equipped with the norm We also require the following assumptions on the coefficients in 1.1 3 .Let a * , a * , φ * , φ * , and K * be positive constants such that C Other assumptions will be made in individual theorems as necessary.

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 h p > 0, we discretize 2.7 in space on a quasiuniform triangle mesh T h p of Ω with diameter of element ≤ h p .Let W h p ⊂ L 2 Ω 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 ∈ W h p such that k μ ∇P, ∇χ q t , χ , ∀χ ∈ W h p .

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 p P , resp.are ensured obviously 15 .
By the theory of the cubic Hermit element, the finite element space W h p possess the following approximation property 15, 16 : where K is a positive constant independent of h p .Define Π: C 2 Ω → W h p is the interplant operator.Then, we have the following theorem.
Theorem 2.1.Let p and P be the solutions of 1.1a -1.1d and 2.8 , respectively, and assume p ∈ H 4 Ω .Then, there exists a positive constant K independent of h p such that ∇ p − P ≤ Kh 3 p p 4 .

2.10
Proof.By Céa lemma, we have

2.11
Then by 2.9 , we can derive

A Second Order Characteristic Method for the Concentration
Now, we define the characteristic lines associated with vector field u 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 X e x, t; • solving the following initial value problem: where v X e x, t; τ , τ : u X e x, t; τ , τ /φ.Next, assuming they exist, we denote by F e by L, resp. the gradient of X e of v, resp.with respect to the space variable x, that is, F e rs x, t; τ : ∂ X e r ∂x s x, t; τ , L rs x, t : ∂v r ∂x s x, t .

2.14
We adopted some propositions and lemmas from 13 .
In order to compute second order approximations of matrices F e and F −1 e , we need the following equations: Lv X e x, t; τ , τ F e x, t; τ .

2.18
By using Liouville's theorem and the chain rule, we obtain

Variational Formulation
From the definition of the characteristic curves and by using the chain rule, it follows that φ dc dτ X e x, t; τ , τ φ ∂c ∂t X e x, t; τ , τ u X e x, t; τ , τ • ∇c X e x, t; τ , τ .

2.22
By writing 1.1b at point X e x, t; τ and time τ and using 2.22 , we have φ dc dτ X e x, t; τ , τ − ∇ • D∇c X e x, t; τ , τ c − c q X e x, t; τ , τ .

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 X : Ω → X Ω , X ∈ C 2 Ω be an invertible vector valued function.Let F ∇X and assume that

The Combined Approximation Scheme
We now present our sequential time-stepping procedure that combines 2.8 and 2.25 .Part J into 0 The analysis is valid for variable time steps, but we drop the superscript from Δt c for convenience.For functions f on Ω × J, we write f n x for f x, t n .As in 3 , let us part J into pressure time steps 0 Each pressure step is also a concentration step, that is, for each m there exists n such that t m t n , in general, Δt p > Δt c .We may vary Δt p , but except for Δt 1 p we drop the superscript.For functions f on Ω × J, we write f m x for f x, t m ; thus, subscripts refer to pressure steps and superscripts to concentration steps.
If concentration step t n relates to pressure steps by t m−1 < t n ≤ t m , we require a velocity approximation for 2.25 based on U m−1 and earlier values.If m ≥ 2, take the linear extrapolation of U m−1 and U m−2 defined by We retain the superscript on Δt 1 p because EU n is first-order correct in time during the first pressure step and second-order during later steps.

International Journal of Mathematics and Mathematical Sciences
Let X : 0, T → R 2 be a solution of the ordinary differential equation dX dt v X, t , where v X, t : EU X, t φ .

2.28
Then, we can write

2.29
Subject to an initial condition X t n 1 x, we get approximate values of X at t n by the Euler method and the second order Runga-Kutta method, respectively,

2.30
Next, assuming they exist, we denote by F n E resp., by L E the gradient of X n E resp., of v x with respect to the space variable x, that is, 13, 14

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 v C 0 W 1,∞ Ω Δt c < 1/2, we can see that

2.33
Corollary 2.11.Under the assumptions of Lemma 2.10, for all x ∈ Ω, we have

2.34
The time difference for 2.25 will be combined with a standard finite element in the space variables.For h c > 0, we discrete 2.25 in space on a quasi-uniform mesh T h c of Ω with diameter of element ≤ h c .Let M h c ⊂ H 1 Ω be a finite element space.
A characteristic discretization of the weak form 2.25 is given by where C 0 is an initial approximation of exact solution c 0 x into M h c , which will be defined in Section 5.
At each pressure time step t m , we define is the truncation of C to 0, 1 .Then at t m , 2.8 is the following: The steps of calculation are as follows.
Step 1. C 0 known → solve U 0 , P 0 by 2.37 , 2.39a and 2.39b ; Step 2. by 2.35a and 2.35b to solve C 1 → then by 2.35a and 2.35b to solve C 2 ; Step 3. analogously, {C j−1 } n 1 j 1 known → {C j } n 1 j 1 such that t n 1 t 1 ; Step 4. then by 2.37 , 2.39a and 2.39b for U 1 , P 1 ; 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 h c , h p , Δt c , and Δt p , but possibly depending on constants in C .Similarly, ε will denote a generic small positive constant.

Stability for the Concentration Equation
In this section, we derive the stability for the concentration equation.For a given series of functions {ϕ} N n 0 , we define the following norms and seminorm: q x, t ϕ, ϕ .

3.1
In the following sections, we use positive constants as In our analysis, we need some lemmas.
Lemma 3.1.Under the definitions 2.26 and 2.30 , for n 0, . . ., N and ϕ ∈ L 2 Ω , it holds that Proof.We only need to show a proof in the case i 1.Let J 1 be the Jacobian matrix of the transformation y X n E x x − v n x Δt c as According to the proof of 3.23 in 4 , we have

International Journal of Mathematics and Mathematical Sciences 11
Since Δt c o h c , for h c sufficiently small, we see that Following 3.6 , we have We complete the proof.
Suppose that {C n } N n 0 ⊂ H 1 Ω and c ∈ C 0 L 2 be given.We define linear forms A n 1/2 h and F n 1/2 h on M h c for n 0, . . ., N − 1 by where ϕ ∈ M h c .
Lemma 3.2.Let {C n } N n 0 be a solution of 2.35a and 2.35b .Then it holds that where D Δt c is the forward difference operator defined by

International Journal of Mathematics and Mathematical Sciences
Proof.Substituting ϕ C n 1 into 3.8 , we have

3.14
Then when I 3 ≥ 0 and I 3 < 0, we have

3.15
Similarly, for I 4 , we obtain the estimate

3.16
International Journal of Mathematics and Mathematical Sciences 13 Analogous computations to term I 2 give 3.17 which completes the proof.
From Lemma 3.1, we have

3.19
which completes the proof of the stability by virtue of Gronwall's inequality.

14
International Journal of Mathematics and Mathematical Sciences Theorem 3.3 stability .Let {C n } N n 0 be the solution of 2.35a and 2.35b subject to the initial value C 0 .Then there exists a positive constant c 1 c 1 EU l ∞ W 1,∞ independent of h c and Δt c such that 3.20

Error Estimate Theorem
Now, we turn to derive an optimal priori error estimate in L 2 -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

4.3
If c is the solution of 1.1a -1.1d , we have for n 0, . . ., N − 1 We decompose e as 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 Let the solution of 4.4 satisfy c ∈ Z 3 , ∇c ∈ Z 3 .Then, for each n 0, 1, . . ., N − 1, there exists function Moreover, ξ n 1/2 A ∈ L 2 Ω and the following estimate holds: where c denotes a constant independent of Δt c .
Lemma 4.3 see 13 .Assume the above Hypotheses hold, and that the coefficients of the problem Let the solution of 4.4 satisfy c ∈ Z 3 , ∇c ∈ Z 3 .Then, for each n 0, 1, . . ., N − 1, there exists function ∈ L 2 Ω and the following estimate holds: where c 1 denotes a constant independent of Δt c .

4.10
Proof.From the definition of A n 1/2 h , we have

4.11
Since it holds that we have To estimate I 2 , we divide it into three parts as

4.15
The term I 21 is estimated as

4.16
Using the transformation y X n E x , we have

4.17
It is clear that 18

4.19
Term I 5 can be divided it into three parts like term I 2 as where

18
International Journal of Mathematics and Mathematical Sciences Thus, the same kind of computations used for I 2 leads to the following estimate:
We now turn to estimate e n .From the definition of A n 1/2 h and A n 1/2 , we see

4.23
From Lemmas 4.2, 4.3 and 4.4, we obtain q j 1 e j 1 q j e j • X j E 2 .

4.26
Using the estimates 4.28

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 U. 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 L 2 -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 F e of the gradient of the characteristic line X e and the inverse matrix F −1 e were used to approximate the diffusion term.The properties of these matrices International Journal of Mathematics and Mathematical Sciences 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.

Δt c 2 8 ∇e j 2 c 1 c h 2k c j 2 k 1 4. 27 for j 0 and n 1 ,Theorem 4 . 5
D∇η j , ∇e j ≤ Δt c Δt and Gronwall's inequality, we derive the following error estimate.error estimate .Let {C n } N n 0 be the solution of 2.35a and 2.35b subject to the initial value C 0 .Then there exists a positive constant c 1 c 1 EU C 0 W 1,∞ independent of h c and Δt c such that then F e satisfies the Taylor expansions as t s τ − s ∇ ∂v ∂t Lv X e x, t; τ , τ F e x, t; τ dτ, 2.17 International Journal of Mathematics and Mathematical Sciences and its inverse, F −1 e , satisfies the Liouville theorem as e X e x, t; s , s; τ dτ.
H 1 X Ω a vector valued function and ψ ∈ H 1 Ω a scalar function.Now, we can multiply 2.23 by a test function ψ ∈ H 1 Ω , integrate in Ω, and apply the usual Green's formula and 2.24 with X x X e x, t; τ , obtaining ΩF −1 e x, t; τ D∇c X e x, t; τ • ∇ψ x dx Ω div F −T e x, t; τ • D∇c X e x, t; τ , τ ψ x dx Ω c − c q X e x, t; τ , τ ψ x dx.2.25