1 Introduction

# Bounded domain problem for the modified Buckley-Leverett equation

## 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 boundary value problem converges to that of the half-line boundary value problem for the MBL equation as . 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.

###### Key words and phrases:
conservation laws, dynamic capillarity, two-phase flows, porous media, shock waves, pseudo-parabolic equations, central schemes
###### 2000 Mathematics Subject Classification:
35L65, 35L67, 35K70, 76S05, 65M06, 65M08
This work is supported in part by NSF Grant DMS-0811003 and an Alfred P. Sloan Fellowship.

## 1. Introduction

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 being defined as

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

In this content, denotes the water saturation (e.g. means pure water, and means pure oil), is a constant which indicates water saturation at , and 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 and with is given in Figure 1.1.

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 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−urfor all u between ul and ur

is the entropy solution, where , are the function values to the left and right of the shock respectively, and the shock speed 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 as given in (1.2) has been well studied (see [14] for an introduction). Let be the solution of , i.e.,

 (1.5) α=√MM+1.

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

1. If , the entropy solution has a single shock at .

2. If , the entropy solution contains a rarefaction between and for and a shock at .

These two types of solutions are shown in Figure 1.2 for .

In either case, the entropy solution of the classical BL equation (1.1) is a non-increasing function of at any given time . 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 be the saturation of water/oil () and assume that the medium is completely saturated, i.e. . 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 denotes the discharge of water/oil with , 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 denotes the absolute permeability, is the relative permeability and 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 is the static capillary pressure and is a positive constant, and is the dynamic effects. Using Corey [6, 20] expressions with exponent 2, , rescaling and combining (1.6)-(1.8), the single equation for the water saturation is

 (1.9) ∂u∂t+∂∂x[u2u2+M(1−u)2]=−∂∂x[ϕ2q2k(1−u)2u2μw(1−u)2+μou2∂∂x(pc(u)ϕ−τ∂u∂t)]

where [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 is given as in (1.2). Notice that, if 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 , the MBL equation (1.10) is an extension involving a third order mixed derivative term . 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

 (1.11) ut+(f(u))x = ϵuxx+ϵ2τuxxtinQ={(x,t):x>0,t>0}u(x,0) = u0(x)x∈[0,∞)u(0,t) = gu(t),limx→∞u(x,t) = 0t∈[0,∞)u0(0) = gu(0)compatibility condition

and the other one is finite interval boundary value problem

 (1.12) vt+(f(v))x = ϵvxx+ϵ2τvxxtin˜Q={(x,t):x∈(0,L),t>0}v(x,0) = v0(x)x∈[0,L]v(0,t) = gv(t),v(L,t) = h(t)t∈[0,∞)v0(0) = gv(0),v0(L) = h(0)compatibility condition.

Considering

 (1.13) u0(x)={v0(x)forx∈[0,L]0forx∈[L,+∞),gu(t)=gv(t)≡g(t),h(t)≡0,

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 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.

## 2. The half line problem versus the finite interval problem

Let be the solution to the half line problem (1.11), and let 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 and on the spatial interval at a given finite time . The main result of this section is

###### Theorem 2.1 (The main Theorem).

If satisfies

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

where and are positive constants, then

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

for some , and , 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 , , and 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 -derivative from the -derivative, and formally solving a first order linear ODE in and a second order non-homogeneous ODE in . 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 ). To solve (1.11), we first rewrite (1.11) by separating the -derivative from the -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 to obtain

 (2.3) (I−ϵ2τ∂2∂x2)(u−e−tϵτu0)=∫t0(1ϵτu−(f(u))x)e−t−sϵτds.

Furthermore, we let

 (2.4) A=u−e−tϵτu0,

then (2.3) can be written as

 (2.5) A′′−1ϵ2τA=∫t0(−1ϵ3τ2u+1ϵ2τ(f(u))x)e−t−sϵτds,where  ′=∂∂x.

Notice that (2.5) is a second-order non-homogeneous ODE in -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 and a kernel for the non-homogeneous terms and 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.7) A(x,t)=−1ϵ3τ2∫t0∫+∞0G(x,ξ)u(ξ,s)e−t−sϵτdξds+1ϵ2τ∫t0∫+∞0K(x,ξ)f(u)e−t−sϵτdξds+(g(t)−e−tϵτg(0))e−xϵ√τ

where the Green’s function and the kernel are

 (2.8) G(x,ξ) = ϵ√τ2(e−x+ξϵ√τ−e−|x−ξ|ϵ√τ), (2.9) K(x,ξ) = −∂G(x,ξ)∂ξ = 12(e−x+ξϵ√τ+sgn(x−ξ)e−|x−ξ|ϵ√τ).

To recover the solution for the half line problem (1.11), we refer to the definition of in (2.4). Thus, the implicit solution formula for the half line problem (1.11) is

 (2.10) u(x,t)=−12ϵ2τ√τ∫t0∫+∞0(e−x+ξϵ√τ−e−|x−ξ|ϵ√τ)u(ξ,s)e−t−sϵτdξds+12ϵ2τ∫t0∫+∞0(e−x+ξϵ√τ+sgn(x−ξ)e−|x−ξ|ϵ√τ)f(u)e−t−sϵτdξds+(g(t)−e−tϵτg(0))e−xϵ√τ+e−tϵτu0(x).

### 2.2. Finite interval problem

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

 (2.11) AL=v−e−tϵτv0,

then it satisfies

 (2.12) (AL)′′−1ϵ2τAL=∫t0(−1ϵ3τ2v+1ϵ2τ(f(v))x)e−t−sϵτdswhere′=∂∂x

with the boundary conditions

 AL(0,t)=v(0,t)−e−tϵτv0(0)=g(t)−e−tϵτg(0),AL(L,t)=v(L,t)−e−tϵτv0(L)=h(t)−e−tϵτh(0).

These boundary conditions affect both the homogeneous solution and the particular solution of (2.12) as follows

 (2.13) AL(x,t)=−1ϵ3τ2∫t0∫L0GL(x,ξ)v(ξ,s)e−t−sϵτdξds+1ϵ2τ∫t0∫L0KL(x,ξ)f(v)e−t−sϵτdξds+c1(t)ϕ1(x)+c2(t)ϕ2(x)

where the Green’s function , the kernel and the bases for the homogeneous solutions are

 (2.14) GL(x,ξ)=ϵ√τ2(e2Lϵ√τ−1)(ex+ξϵ√τ+e2L−(x+ξ)ϵ√τ−e|x−ξ|ϵ√τ−e2L−|x−ξ|ϵ√τ),
 (2.15) KL(x,ξ)=−12(e2Lϵ√τ−1)(ex+ξϵ√τ−e2L−(x+ξ)ϵ√τ+sgn(x−ξ)e|x−ξ|ϵ√τ−sgn(x−ξ)e2L−|x−ξ|ϵ√τ),
 (2.16) c1(t)=g(t)−e−tϵτg(0),c2(t)=h(t)−e−tϵτh(0), (2.17) ϕ1(x)=eL−xϵ√τ−e−L+xϵ√τeLϵ√τ−e−Lϵ√τ,andϕ2(x)=exϵ√τ−e−xϵ√τeLϵ√τ−e−Lϵ√τ.

Thus, the implicit solution formula for the finite interval problem (1.12) is

 (2.18) v(x,t)=−12ϵ2τ√τ(e2Lϵ√τ−1)∫t0∫L0(ex+ξϵ√τ+e2L−(x+ξ)ϵ√τ−e|x−ξ|ϵ√τ−e2L−|x−ξ|ϵ√τ)v(ξ,s)e−t−sϵτdξds−12ϵ2τ(e2Lϵ√τ−1)∫t0∫L0(ex+ξϵ√τ−e2L−(x+ξ)ϵ√τ+sgn(x−ξ)e|x−ξ|ϵ√τ−sgn(x−ξ)e2L−|x−ξ|ϵ√τ)f(v)e−t−sϵτdξds+c1(t)ϕ1(x)+c2(t)ϕ2(x)+e−tϵτv0(x).

### 2.3. Comparisons

In this section, we will prove that the solution to the half line problem can be approximated as accurately as one wants by the solution 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 ( respectively) into two parts: and ( and respectively), such that ( respectively) enjoys zero initial condition and boundary conditions at and . We estimate the difference between and by estimating the differences between and , and , 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 as sum of two terms and , 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 and are given in (2.16) and (2.17) respectively. With this definition, takes care of the initial condition and boundary conditions at and for . Then satisfies an equation slightly different from the equation satisfies in (1.11):

 (2.20) Ut−ϵUxx−ϵ2τUxxt=(ut−ϵuxx−ϵ2τuxxt)−((uL)t−ϵ(uL)xx−ϵ2τ(uL)xxt)=−(f(u))x+1ϵτuL(x,t)

In addition, has zero initial condition and boundary conditions at and , i.e.,

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

Similarly, for , 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 , and , are given in (2.16) and (2.17) respectively. With this definition, takes care of the initial condition and boundary conditions and at and for . Then satisfies an equation slightly different from the equation 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 and , we define

 W(x,t)=V(x,t)−U(x,t)forx∈[0,L].

Because of (2.20) and (2.23), we have

 (2.25) Wt−ϵWxx−ϵ2τWxxt=−(f(v)−f(u))x+1ϵτ(vL−uL).

In lieu of (2.21) and (2.24), also has zero initial condition and boundary conditions at and , i.e.,

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

Now, to estimate , we can estimate and estimate 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 and satisfies

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

where and are positive constants. Notice that the constraint is crucial in Lemmas 2.3, 2.4.

where  and .

1. .

2. .

1. .

2. .

3. .

###### Lemma 2.5.
1. .

2.  for  .

3.  if   for  .

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 in section 2.3.3. By comparing and given in (2.19) and (2.22) respectively, it is clear that the coefficient for appeared in (2.19) needs to be compared with the corresponding coefficient for appeared in (2.22). We thus define a space-dependent function

 (2.29) Uc2(x,t)=u(x,t)−c1(t)e−xϵ√τ−e−tϵτu0(x)

and establish the following proposition

###### Proposition 2.6.
 (2.30) ∣∣Uc2(L,t)∣∣≤aτ(t)ebτtϵτe−λLϵ√τ+cτtϵτe(bτ−1)tϵτe−λ(L−L0)ϵ√τ

for some parameter-dependent constants