A Numerical Method Based on the Range-Discrete Grid for One-Dimensional Buckley-Leverett Equation

A numerical method for solving the one-dimensional Buckley-Leverett equation arising in the process of displacement of oil by water is presented. Instead of using the traditional spatial discrete grids, the numerical algorithm is built on a “range discrete” grid, which is obtained by direct subdivisions in the saturation domain. The range discrete grid describes the discontinuities explicitly and completely, and has an adaptive nature in smooth regions. Grid points are divided into two classes: continuous points and discontinuous points. Numerical solution of the Buckley-Leverett equation is achieved by moving continuous points by tracing characteristics and moving discontinuous points by tracking discontinuities. Numerical examples are presented, and the solutions obtained by the proposed method are found of high precision. Especially, shocks are solved with no dissipation, and the sharpness is maintained.


Introduction
The displacement of oil by water in porous media can be described by the Buckley-Leverett (B-L) theory [1], [2], which is an important theory for underground water follow and oil recovery. When capillary effect and gravity is ignored, the B-L equation becomes where S is water saturation and f S is the fractionalflow function, which is typically S-shaped. Eq. (1) is a nonlinear hyperbolic equation and discontinuities may appear in finite time even the initial condition is smooth enough [3], [4]. When discontinuities appear, numerical methods based on traditional spatial discrete grid, i.e. Finite Difference Method (FDM), suffer either massive numerical dissipation or unphysical oscillations [5]- [10]. This is because that the spatial discrete grid describes a discontinuity within a finite-small cell, which, with no doubt, causes numerical problems.
Here, we proposed a numerical method that maintains sharp discontinuities for the one-dimensional B-L equation by introducing a uniformly distributed "rangediscrete" grid. The range-discrete grid shows an adaptive nature in smooth regions. Moreover, discontinuities can be solved with no diffusion. The proposed method combines the method of characteristic and shock tracking method. The numerical method is substantially a moving mesh method and thus abandons the use of traditional spatial discrete grid.  As shown in Fig. 1, a set of intersection points , , ,

The range-discrete grid
exist at a discontinuity. These points are treated as different grid points while sharing the same spatial coordinate value, which is L M Using the rules explained above, a new discrete grid for a general saturation field is constructed. The saturation field is numerically approximated when x j , the spatial coordinate of the grid points 0,1,2, , j A j n n , , with each assigned to a discrete saturation value S i , are known. To differ from the spatially discretized Eulerian grid, we call this grid system a "range-discrete" grid.
Grid points in the range-discrete grid go into two classes intuitively: points in smooth regions are continuous points; while points with a same spatial coordinate are marked as discontinuous points to represent a discontinuity. In this way, the traditional way that describes a discontinuity within limited small space region is abandoned, so the sharpness will be maintained. Furthermore, grid points will gather where saturation various quickly in space and vice versa, which is analogous to the adaptive mesh refinement strategy [11], [12]. While complex refinement strategies have to be processed in these methods, the range-discrete grid is born with an adaptive nature.

Numerical method based on the rangediscrete grid
Using the technique illustrated in section 2, the saturation field at 0 t can be approximated by a range-discrete grid. The movements of the grid points A j can be tracked for a small time step t ' . Once each t t j x ' calculated, approximation of the saturation field at t t ' is obtained. The B-L equation can be solved by repeating this. One should notice that, grid points move with their saturation values unchanged, which is the main difference from the numerical methods based on spatial discrete grids.
In each time step, movements of the grid points can be established explicitly and with little computation cost. However, this may result in multi-valued curves in space due to the nonlinear hyperbolic nature of the B-L equation, so that a correction procedure has to be introduced. Consequently, there are two steps involved in each time step. The first is a prediction step to move these points explicitly and the second is a correction step to solve the multi-value problem.

The prediction step
If marked as a continuous point at t, point A j will remain to be a continuous point and can be traced along its characteristic line, giving j As a set of discontinuous points represents a discontinuity, their movements are determined by the evolution of this discontinuity, which can be established by entropy condition [4]. That is, the points in the part that will change into a shock wave remain to be discontinuous points, moving at the velocity of the shock v shock achieved from Rankine-Hugoniot condition; the points in the part that will change into a rarefaction wave will become continuous points, moving along its characteristic line. Altogether, the spatial coordinates of points A j at t t ' are with if is a discontinuous point after if is a continuous point after

The correction step
As shown in Fig. 2(a) impossible. The real saturation curve is spatially discontinuous [1], [13]. Considering that f S typically have inflection points, we combine the material balance method suggested by Buckley and Leverett with entropy condition to solve the multi-value problem after the prediction step.

Correction procedure of two continuous points
As a numerical method, it is quite straightforward to mark the two points in Fig. 2(a) as discontinuous points with their spatial coordinates assigned to a same value x s , representing the discontinuity that will be introduced to correct the multi-valued portion to a locally single-valued one, as shown in dash in Fig. 2(a).
We should mention that only discontinuities that satisfy entropy condition can be produced in this correction step, so that an entropy checking is needed. It is easy to prove that if one continuous point goes ahead another one, that is point 2 advances faster than point 3 in this situtation, the entropy condition is always checked. So, these two points are admitted to form a new discontinuity here.
After this, x s can be calculated using the material balance principle: where x * and x represent the spatial coordinates before and after the correction procedure. Eq. (4) maintains the conservation of mass in the correction step. The integral in Eq. (4) can be calculated by linear interpolation or a higher resolution can be expected with the interpolation introduced by Forjoun etc. [14].
where the subscript a and b represent two adjacent points. Thus, with the position of the new discontinuity calculated, the unphysical multi-value problem is solved.

Correction procedures when discontinuous points involved
After the prediction step, the multi-valued problem happens not only between continuous points(see Fig.  2(a)), but also exits when discontinuous points occurred (see Fig. 2(b) , 2(c) and 2(d)). Actually, as a set of discontinuous points represents a single discontinuity, it is reasonable to treat the set of discontinuous points as a whole here. If we do so, things in Fig. 2(b), 2(c) and 2(d) become exactly the same as the situation in Fig. 2(a). Consequently, the correction procedures when discontinuous points involved are actually the same as the procedure illustrated in 3.2.1, with only one more thing to notice.
That is, if we mark all the points in the multi-valued portion of the curve as discontinuous points, it may produce a discontinuity that does not satisfy the entropy condition. So, from the most upstream point in the multi-valued portion to the most downstream point, an entropy checking is processing after marking them as discontinuous points: take the discrete saturation values of each point and the most downstream point in the multivalued portion as the upstream and downstream value of a discontinuity respectively. If this discontinuity does not satisfy the entropy condition, this point will be marked as a continuous point until an entropy satisfied discontinuity occurs. In this way, the resulting discontinuity is enforced to satisfy the entropy condition.
If the time step t ' is small enough, the multi-value problem only exists between adjacent points so that all possible situations are shown in Fig. 2. So, with the technique illustrated above, the multi-valued curve after the prediction step can be corrected to a single-valued saturation curve, giving the approximation of the saturation curve at 0 t t ' . With mass balance principle and entropy condition applied explicitly in the correction step, the solution achieved is thought to be mass conserved and satisfy the entropy condition.  Here, the fractional-flow function is set as . The initial condition is set to with no discontinuities as The solutions obtained by FDM and the proposed method at different time are shown in Fig. 3, where the reference solution is obtained by the FDM with a refined grid. As we can see, the solution develops a moving shock. The shock is solved by the proposed method with no diffusion. Using equal number of grid points, first order FDM gives almost a smooth solution. We should mention that, as no linear equations has to be solved, the proposed method cost even much less computation effort than the first order FDM.

Two initial discontinuities
In this example, f S is designed to be a third order polynomial as 2 3 The initial saturation field consists of two discontinuities given by  Figure 4. Solutions at different time from the two initial discontinuities. (a) t=0.15 (b) t=0. 6. Grid points of the proposed method at a discontinuity is represented by two grid points at both ends for simplicity.
At t>0, each of these two discontinuities will change into a shock followed by a rarefaction and will interact with each other. The solutions at t=0.15 and t=0.6 are shown in Fig. 4. While the FDM gives rather wrong answers, the proposed method solves the position and strength of the shocks accurately with even fewer discrete points, with sharpness of the shocks kept and grid points located adaptively.

Conclusions
A range-discrete grid system and a numerical method based on the grid for the one-dimensional nonlinear hyperbolic B-L equation is presented. The numerical method solves the discontinuities sharply, and has an adaptive nature in smooth regions. With mass balance principle and the entropy condition employed explicitly, the solution is mass conserved and satisfies the entropy condition.
Our future work is to expand this method to higher dimensions, as it is straightforward that, in two dimensional problems, the range-discrete grid becomes contour lines in the saturation map.