Abstract.

The focus of the present study is the modified Buckley-Leverett
(MBL) equation describing two-phase flow in porous media. The MBL
equation differs from the classical Buckley-Leverett (BL) equation
by including a balanced diffusive-dispersive combination. The
dispersive term is a third order mixed derivatives term, which
models the dynamic effects in the pressure difference between the
two phases. The classical BL equation gives a monotone water
saturation profile for any Riemann problem; on the contrast, when
the dispersive parameter is large enough, the MBL equation
delivers non-monotone water saturation profile for certain Riemann
problems as suggested by the experimental observations. In this
paper, we first show that the solution of the finite interval
[0,L] boundary value problem converges to that of the half-line
[0,+∞) boundary value problem for the MBL equation as
L→+∞. This result provides a justification for
the use of the finite interval boundary value problem in numerical
studies for the half line problem. Furthermore, we extend the
classical central schemes for the hyperbolic conservation laws to
solve the MBL equation which is of pseudo-parabolic type.
Numerical results confirm the existence of non-monotone water
saturation profiles consisting of constant states separated by
shocks.

2000 Mathematics Subject Classification:

The classical Buckley-Leverett (BL) equation [3] is a simple model for two-phase fluid flow in a porous medium. One application is secondary recovery by water-drive in oil reservoir simulation. In one space dimension the equation has the standard conservation form

ut+(f(u))x=0

in

Q={(x,t):x>0,t>0}

(1.1)

u(x,0)=0

x∈(0,∞)

u(0,t)=uB

t∈[0,∞)

with the flux function f(u) being defined as

(1.2)

f(u)=⎧⎪
⎪
⎪⎨⎪
⎪
⎪⎩0u<0,u2u2+M(1−u)20≤u≤1,1u>1.

In this content, u:¯Q→[0,1] denotes the water saturation (e.g. u=1 means pure water, and u=0 means pure oil), uB is a constant which indicates water saturation at x=0, and M>0 is the water/oil viscosity ratio. The classical BL equation (1.1) is a prototype for conservation laws with convex-concave flux functions.
The graph of f(u) and f′(u) with M=2 is given in Figure 1.1.

(a) f(u)=u2u2+M(1−u)2

(b) f′(u)=2Mu(1−u)(u2+M(1−u)2)2

Figure 1.1. f(u) and f′(u) with M=2.

Due to the possibility of the existence of shocks in the solution of the hyperbolic conservation laws (1.1), the weak solutions are sought. The function u∈L∞(Q) is called a weak solution of the conservation laws (1.1) if

∫Q{u∂ϕ∂t+f(u)∂ϕ∂x}=0for all ϕ∈C∞0(Q).

Notice that the weak solution is not unique. Among the weak
solutions, the entropy solution is physically relevant and unique.
The weak solution that satisfies Oleinik entropy condition [19]

(1.3)

f(u)−f(ul)u−ul≥s≥f(u)−f(ur)u−urforallubetweenulandur

is the entropy solution, where ul, ur are the function values to the left and right of the shock respectively, and the shock speed s satisfies Rankine-Hugoniot jump condition [17, 10]

(1.4)

s=f(ul)−f(ur)ul−ur.

The classical BL equation (1.1) with flux function f(u) as
given in (1.2) has been well studied (see [14]
for an introduction). Let α be the solution of
f′(u)=f(u)u, i.e.,

(1.5)

α=√MM+1.

The entropy solution of the classical BL equation can be classified into two categories:

If 0<uB≤α, the entropy solution has a single shock at xt=f(uB)uB.

If α<uB<1, the entropy solution contains a rarefaction between uB and α for f′(uB)<xt<f′(α) and a shock at xt=f(α)α.

Figure 1.2. The entropy solution of the classical BL equation (M=2,α=√23≈0.8165). \subrefBL1
0<uB=0.7≤α, the solution consists of one shock at
xt=f(uB)uB; \subrefBL2
α<uB=0.98<1, the solution consists of a rarefaction
between uB and α for f′(uB)<xt<f′(α)
and a shock at xt=f(α)α.

In either case, the entropy solution of the classical BL equation (1.1) is a non-increasing function of x at any given time t>0.
However, the experiments of two-phase flow in porous medium reveal complex infiltration profiles, which may involve overshoot, i.e., profiles may not be monotone [7]. This suggests the need of modification to the classical BL equation (1.1).

To better describe the infiltration profiles, we go back to the origins of (1.1). Let Si be the saturation of water/oil (i=w,o) and assume that the medium is completely saturated, i.e. Sw+So=1. The conservation of mass gives

(1.6)

ϕ∂Si∂t+∂qi∂x=0

where ϕ is the porosity of the medium (relative volume occupied by the pores) and qi denotes the discharge of
water/oil with qw+qo=q, which is assumed to be a constant in space due to the complete saturation assumption.
Throughout of this work, we consider it constant in time as well.
By Darcy’s law

(1.7)

qi=−kkri(Si)μi∂Pi∂x,i=w,o

where k denotes the absolute permeability, kri is the relative permeability and μi is the viscosity of water/oil.
Instead of considering constant capillary pressure as adopted by the classical BL equation (1.1),
Hassanizadeh and Gray [8, 9] have defined the dynamic capillary pressure as

(1.8)

Pc=Po−Pw=pc(Sw)−ϕτ∂Sw∂t

where pc(Sw) is the static capillary pressure and
τ is a positive constant, and ∂Sw∂t is the dynamic effects. Using Corey [6, 20]
expressions with exponent 2, krw(Sw)=S2w,kro(So)=S2o, rescaling xϕq→x and combining
(1.6)-(1.8), the single equation for the
water saturation u=Sw is

where M=μwμo[22]. Linearizing the right
hand side of (1.9) and rescaling the equation as in
[21, 20],
the modified Buckley-Leverett equation (MBL) is derived as

(1.10)

∂u∂t+∂f(u)∂x=ϵ∂2u∂x2+ϵ2τ∂3u∂x2∂t

where the water fractional flow function f(u) is given as in
(1.2). Notice that, if Pc in (1.8) is taken to
be constant, then (1.9) gives the classical BL
equation; while if the dispersive parameter τ is taken to be
zero, then (1.10) gives the viscous BL equation, which still
displays monotone water saturation profile. Thus, in addition to
the classical second order viscous term ϵuxx, the MBL
equation (1.10) is an extension involving a third order mixed derivative term
ϵ2τuxxt. Van Dujin et al. [21] showed
that the value τ is critical in determining the type of the
solution profile. In particular, for certain Riemann problems, the solution
profile of (1.10) is not monotone when τ is larger than
the threshold value τ∗, where τ∗ was numerically
determined to be 0.61 [21]. The non-monotonicity of the
solution profile is consistent with the experimental observations
[7].

The classical BL equation (1.1) is hyperbolic, and the
numerical schemes for hyperbolic equations have been well
developed (e.g. [14, 15, 4, 5, 18, 12]
). The MBL
equation (1.10), however, is pseudo-parabolic, we will
illustrate how to extend the central schemes [18, 12, 13] to solve
(1.10) numerically. Unlike the finite domain of dependence
for the classical BL equation (1.1), the domain of dependence for the MBL equation (1.10) is infinite.
This naturally raises the question for the choice of
computational domain. To answer this question, we will first study the
MBL equation equipped with two types of domains and corresponding
boundary conditions. One is the half line boundary value problem

we will show the relation between the solutions of problems
(1.11) and (1.12). To the best knowledge of
the authors, there is no such study for MBL equation (1.10).
Similar questions were answered for BBM equation [1, 2].

The organization of this paper is as follows. Section
2 will bring forward the exact theory comparing the
solutions of (1.11) and (1.12). The
difference between the solutions of these two types of problems
decays exponentially with respect to the length of the interval
L for practically interesting initial profiles. This provides a
theoretical justification for the choice of the computational
domain. In section 3, high order central schemes will
be developed for MBL equation in finite interval domain. We
provide a detailed derivation on how to extend the central schemes
[18, 12] for conservation laws to solve the MBL equation
(1.10). The idea of adopting numerical schemes originally
designed for hyperbolic equations to pseudo-parabolic equations is
not restricted to central type schemes only ([23, 24]).
The numerical results in section 4 show that the
water saturation profile strongly depends on the dispersive
parameter τ value as studied in [21]. For
τ>τ∗, the MBL equation (1.10) gives non-monotone
water saturation profiles for certain Riemann problems as
suggested by experimental observations [7]. Section
5 gives the conclusion of the paper and the
possible future directions.

Let u(x,t) be the solution to the half line
problem (1.11), and let v(x,t) be the solution
to the finite interval problem (1.12). We
consider the natural assumptions (1.13). The goal of
this section is to develop an estimate of the difference between
u and v on the spatial interval [0,L] at a given finite time
t. The main result of this section is

Theorem 2.1 (The main Theorem).

If u0(x) satisfies

(2.1)

u0(x)={Cux∈[0,L0]0x>L0

where L0<L and Cu, are positive constants, then

∥u(⋅,t)−v(⋅,t)∥H1L,ϵ,τ≤D1;ϵ,τ(t)e−λLϵ√τ+D2;ϵ,τ(t)e−λ(L−L0)ϵ√τ

for some 0<λ<1,
D1;ϵ,τ(t)>0 and D2;ϵ,τ(t)>0, where

∥Y(⋅,t)∥H1L,ϵ,τ:=√∫L0Y(x,t)2+(ϵ√τYx(x,t))2dx

Notice that
the initial condition (2.1) we considered is the Riemann
problem. Theorem 2.1 shows that the solution to the
half line problem (1.11) can be approximated as
accurately as one wants by the solution to the finite interval
problem (1.12) in the sense that
D1;ϵ,τ(t), D2;ϵ,τ(t), λLϵ√τ and λ(L−L0)ϵ√τ can be controlled.

To prove theorem 2.1, we first derive the implicit solution formulae for the half line problem and the finite interval problem in section 2.1 and section 2.2 respectively. The implicit solution formulae are in integral form, which are derived by separating the x-derivative from the t-derivative, and formally solving a first order linear ODE in t and a second order non-homogeneous ODE in x. In section 2.3, we use Gronwall’s inequality multiple times to obtain the desired result in theorem 2.1.

2.1. Half line problem

In this section, we derive the implicit solution
formula for the half line problem (1.11) (with gu(t)=g(t)).
To solve (1.11), we first rewrite (1.11) by separating the x-derivative from the
t-derivative,

(2.2)

(I−ϵ2τ∂2∂x2)(ut+1ϵτu)=1ϵτu−(f(u))x.

By using integrating factor method, we formally integrate (2.2) over [0,t] to obtain

Notice that (2.5) is a second-order non-homogeneous ODE in
x-variable along with the boundary conditions

(2.6)

A(0,t)=u(0,t)−e−tϵτu0(0)=g(t)−e−tϵτg(0),A(∞,t)=u(∞,t)−e−tϵτu0(∞)=0.

To solve (2.5), we first solve the corresponding linear
homogeneous equation with the non-zero boundary conditions
(2.6). We then find a particular solution for the
non-homogeneous equation with zero boundary conditions by
introducing a Green’s function G(x,ξ) and a kernel K(x,ξ)
for the non-homogeneous terms u and
(f(u))x respectively. Combining the solutions for the two
non-homogeneous terms and the homogeneous part with boundary
conditions, we get the solution for equation (2.5)
satisfying the boundary conditions (2.6):

2.2. Finite interval problem

The implicit solution for the finite interval problem (1.12) (with gv(t)=g(t)) can be solved in a similar way.
The only difference is that the additional boundary condition h(t) at x=L in (1.12) gives different boundary conditions for the non-homogeneous ODE in x-variable. Denote

2.3. Comparisons

In this section, we will prove that the solution u(x,t) to the half line problem can
be approximated as accurately as one wants by the solution
v(x,t) to the finite interval problem as stated
in Theorem 2.1.

Due to the difference in the integration domains, we do not use
(2.10) and (2.18) directly for the comparison.
Instead, we decompose
u(x,t) (v(x,t) respectively) into two parts: U(x,t) and uL(x,t)
(V(x,t) and vL(x,t) respectively),
such that U(x,t) (V(x,t) respectively) enjoys zero
initial condition and boundary conditions at x=0 and x=L.
We estimate the difference between u(⋅,t) and v(⋅,t) by estimating the differences between
uL(⋅,t) and vL(⋅,t), U(⋅,t) and V(⋅,t), then applying the triangle inequality.

Definitions and lemmas

To assist the proof of Theorem 2.1 in section 2.3.3, we introduce some new notations in this section. We first decompose u(x,t) as sum of two terms U(x,t) and uL(x,t), such that

u(x,t)=U(x,t)+uL(x,t)x∈[0,+∞)

where

(2.19)

uL=e−tϵτu0(x)+c1(t)e−xϵ√τ+(u(L,t)−c1(t)e−Lϵ√τ−e−tϵτu0(L))ϕ2(x)

and c1(t) and ϕ2(x) are given in (2.16) and
(2.17) respectively. With this definition, uL takes care
of the initial condition u0(x) and boundary conditions g(t)
at x=0 and x=L for u(x,t).
Then U satisfies an equation slightly different from the
equation u satisfies in (1.11):

In addition, U(x,t) has zero initial condition and boundary conditions at x=0 and x=L, i.e.,

(2.21)

U(x,0)=0,U(0,t)=0,U(L,t)=0.

Similarly, for v(x,t), let

v(x,t)=V(x,t)+vL(x,t)x∈[0,L]

where

(2.22)

vL=e−tϵτv0(x)+c1(t)ϕ1(x)+c2(t)ϕ2(x)

and c1(t), c2(t) and ϕ1(x), ϕ2(x) are given in
(2.16) and (2.17) respectively. With this definition,
vL takes care of the initial condition v0(x) and boundary
conditions g(t) and h(t) at x=0 and x=L for v(x,t).
Then V satisfies an equation slightly different from the
equation v satisfies in (1.12):

(2.23)

Vt−ϵVxx−ϵ2τVxxt=−(f(v))x+1ϵτvL(x,t)

with

(2.24)

V(x,0)=0,V(0,t)=0,V(L,t)=0.

Since, in the end, we want to study the difference between U(x,t) and V(x,t), we define

In lieu of (2.21) and (2.24), W(x,t) also has zero initial condition and boundary conditions at x=0 and x=L, i.e.,

(2.26)

W(x,0)=0,W(0,t)=0,W(L,t)=0.

Now, to estimate ∥u−v∥, we can estimate
∥W∥=∥V−U∥ and estimate ∥uL−vL∥ separately.
These estimates are done in section 2.3.3.

Next, we state the lemmas needed in the proof of Theorem
2.1. The proof of the lemmas can be found in the
appendix A and [22]. In all the lemmas, we
assume 0<λ<1 and u0(x) satisfies

(2.27)

u0(x)={Cux∈[0,L0]0x>L0

where L0<L and Cu are positive constants. Notice that the constraint
λ∈(0,1) is crucial in Lemmas 2.3,
2.4.

Lemma 2.2.

f(u)=u2u2+M(1−u)2≤Du where
D=f(α)α and α=√MM+1.

Lemma 2.3.

∫+∞0∣∣∣e−x+ξϵ√τ−e−|x−ξ|ϵ√τ∣∣∣eλx−λξϵ√τdξ≤2ϵ√τ1−λ2.

∫+∞0∣∣∣e−x+ξϵ√τ−e−|x−ξ|ϵ√τ∣∣∣eλx−ξϵ√τdξ≤ϵ√τe(1−λ) .

∫+∞0∣∣∣e−x+ξϵ√τ−e−|x−ξ|ϵ√τ∣∣∣eλxϵ√τ|u0(ξ)|dξ≤2Cuϵ√τeλL0ϵ√τ
.

Lemma 2.4.

∫+∞0∣∣∣e−x+ξϵ√τ+sgn(x−ξ)e−|x−ξ|ϵ√τ∣∣∣eλx−λξϵ√τdξ≤2ϵ√τ1−λ2
.

∫+∞0∣∣∣e−x+ξϵ√τ+sgn(x−ξ)e−|x−ξ|ϵ√τ∣∣∣eλx−ξϵ√τdξ≤ϵ√τ+ϵ√τe(1−λ) .

∫+∞0∣∣∣e−x+ξϵ√τ+sgn(x−ξ)e−|x−ξ|ϵ√τ∣∣∣eλxϵ√τ|u0(ξ)|dξ≤2Cuϵ√τeλL0ϵ√τ .

Lemma 2.5.

∣∣∣ϕ1(x)−e−xϵ√τ∣∣∣=e−Lϵ√τ|ϕ2(x)| .

|ϕ2(x)|≤1 for x∈[0,L] .

∣∣ϕ′2(x)∣∣≤2ϵ√τ if ϵ≪1 for x∈[0,L] .

Last but not least, the norm that we will use in Theorem 2.1 and its proof is

(2.28)

∥Y(⋅,t)∥H1L,ϵ,τ:=√∫L0Y(x,t)2+(ϵ√τYx(x,t))2dx.

A proposition

In this section, we will give a critical estimate, which is essential in the calculation of maximum difference ∥uL(⋅,t)−vL(⋅,t)∥∞ in section 2.3.3.
By comparing uL(x,t) and vL(x,t) given in (2.19) and (2.22) respectively, it is clear that the coefficient u(L,t)−c1(t)e−Lϵ√τ−e−tϵτu0(L) for ϕ2(x) appeared in (2.19) needs to be compared with the corresponding coefficient c2(t) for ϕ2(x) appeared in (2.22).
We thus define a space-dependent function