Bounded domain problem for the modified Buckley-Leverett equation
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
The classical Buckley-Leverett (BL) equation  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
with the flux function being defined as
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
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 
The entropy solution of the classical BL equation can be classified into two categories:
If , the entropy solution has a single shock at .
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 . 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
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
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
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
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.  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 . The non-monotonicity of the solution profile is consistent with the experimental observations .
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
and the other one is finite interval 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 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 . For , the MBL equation (1.10) gives non-monotone water saturation profiles for certain Riemann problems as suggested by experimental observations . 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).
where and are positive constants, then
for some , and , where
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
By using integrating factor method, we formally integrate (2.2) over to obtain
Furthermore, we let
then (2.3) can be written as
Notice that (2.5) is a second-order non-homogeneous ODE in -variable along with the boundary conditions
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):
where the Green’s function and the kernel are
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
then it satisfies
with the boundary conditions
These boundary conditions affect both the homogeneous solution and the particular solution of (2.12) as follows
where the Green’s function , the kernel and the bases for the homogeneous solutions are
Thus, the implicit solution formula for the finite interval problem (1.12) is
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
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):
In addition, has zero initial condition and boundary conditions at and , i.e.,
Similarly, for , let
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):
Since, in the end, we want to study the difference between and , we define
Now, to estimate , we can estimate and estimate separately. These estimates are done in section 2.3.3.
where and .
if for .
Last but not least, the norm that we will use in Theorem 2.1 and its proof is
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
and establish the following proposition
for some parameter-dependent constants