Convergence of Fully discrete schemes for diffusive dispersive conservation laws with
discontinuous coefficient
Abstract.
We are concerned with fullydiscrete schemes for the numerical approximation of diffusivedispersive hyperbolic conservation laws with a discontinuous flux function in onespace dimension. More precisely, we show the convergence of approximate solutions, generated by the scheme corresponding to vanishing diffusivedispersive scalar conservation laws with a discontinuous coefficient, to the corresponding scalar conservation law with discontinuous coefficient. Finally, the convergence is illustrated by several examples. In particular, it is delineated that the limiting solutions generated by the scheme need not coincide, depending on the relation between diffusion and the dispersion coefficients, with the classical KružkovOleĭnik entropy solutions, but contain nonclassical undercompressive shock waves.
Key words and phrases:
conservation laws; discontinuous flux; diffusivedispersive approximation; finite difference scheme; convergence; entropy condition; nonclassical shock]rajib.ami@gmail.com ]ujjwal@math.tifrbng.res.in ]deep@math.tifrbng.res.in
1. Introduction
In this paper, we consider a finite difference method for the vanishing diffusivedispersive approximations of scalar conservation laws with a discontinuous flux
(1.1) 
when tends to zero with as . Here is a regularization term, depends upon two parameters and referred to as the diffusion and the dispersion coefficients, motivated by the equations of twophase flow in porous media, is fixed, is the unknown scalar map, the initial data, is a spatially varying (discontinuous) coefficient, and the flux function is a sufficiently smooth scalar function (see Section 2 for the complete list of assumptions).
Motivated by the dynamic capillary pressure [7], we consider in this paper the simplified model
(1.2) 
with a thirdorder mixed derivatives term including one time derivative. Here are fixed parameters. The equation (1.1) along with (1.2) serves as a concrete model of two phase flows in a heterogeneous porous medium.
Moreover, drawing preliminary motivation from phase transition dynamics, we also consider the following specific form of the regularization term
(1.3) 
with fixed.
Furthermore, for the simplicity in the exposition, we assume that the flux function has the following particular form
Note that the flux has a possibly discontinuous spatial dependence through the coefficient , which is allowed to have jump discontinuities.
The scalar conservation laws with a discontinuous flux function
(1.4) 
is a special example of this type of problems, corresponds to the case . A simple physical model corresponding to (1.4) is the Witham model of car traffic flow on a highway (consult the monograph by Leveque [21]), where the spatially varying coefficient corresponds to changing road conditions. Several other models such as two phase flow in a heterogeneous porous medium that arise in petroleum industry, the modeling of the clarifier thickener unit used in waste water treatment plants are also corresponding to (1.4).
Independently of the smoothness of the initial data and , solutions to (1.4) are not necessarily smooth due to the presence of nonlinear flux term in the equation (1.4). Thus, weak solutions must be sought.
Definition 1.1 (Weak solution).
A weak solution of the initial value problem (1.4) is a bounded measurable function satisfying
(1.5) 
for all .
It is well known that (weak) solutions may be discontinuous and they are not uniquely determined by their initial data. Consequently, an entropy condition must be imposed to single out the physically correct solution. If is “smooth”, a weak solution satisfies the entropy condition if for all convex functions
where is defined by .
By standard limiting argument, this implies the Kružkovtype entropy condition
holds for all .
However, the notion of entropy solution described above breaks down when is discontinuous. In view of [14], we use the following notion of entropy solution for (onedimensional) conservation laws with discontinuous flux equations with coefficients that are only spatially dependent. We assume that the spatially varying coefficient is piecewise with finitely many jumps (in and ), located at .
Definition 1.2 (Entropy solution).
A weak solution of the initial value problem (1.4) is called an entropy solution, if the following Kružkovtype entropy inequality holds for all and all test functions .
(1.6)  
The last couple of decades have witnessed remarkable advances in the studies of conservation laws with discontinuous flux function. However, we will not be able to discuss the whole literature here, but only refer to the parts that are pertinent to the current paper. In case of “smooth” , the notion of entropy solution was introduced independently by Kružkov [18] and Vol’pert [31] (the latter author considered the smaller BV class). These authors also proved general existence, uniqueness, and stability results for the entropy solution, see also Oleĭnik [24] for similar results in the convex case .
1.1. Diffusive Dispersive Approximation
It is well known that the conservation law (1.4) is derived by neglecting underlying small scale effects such as diffusion, dispersion, capillarity etc., and may admit physically relevant discontinuous solutions containing shock waves (nonclassical shock) that may depend on underlying smallscale mechanisms. It has been successively recognized that a standard entropy inequality (due to Kružkov, Oleĭnik, and others) does not suffice to single out such a physically relevant solution, and it is important to incorporate these smallscale effects in the entropy condition. In other words, additional admissibility criteria (a kinetic relation) are required in order to characterize these smallscale dependent non classical shock waves uniquely. In [9], the authors developed a framework for the existence and uniqueness of the nonclassical shock waves that arise as limits of diffusivedispersive approximations.
Noting that the solutions to (1.4) can be different due to their explicit dependence on the underlying small scale effects, we focus on a concrete model of two phase flow in porous medium (for a brief derivation of this model consult [4]). The relevant small scale effect is a dynamic capillary pressure term, that was introduced in [7]. Compared to the standard capillary pressure models [1], the addition of the new term resulted in a model that contain higherorder mixed spatiotemporal derivatives (cf. (1.2)).
The diffusivedispersive model has a long tradition, starting with the analysis of linear diffusiondispersion model (1.1). A pioneering study of the effect of vanishing diffusion and dispersion terms in scalar conservation laws, with independent flux function, can be found in Schonbek [26]. The technique of compensated compactness was used to prove convergence toward weak solutions. Kondo and LeFloch [17] studied zero diffusiondispersion limits for independent fluxes under an optimal balance between the sizes of the diffusion and dispersion parameters. LeFloch and Natalini [19] used the concept of measurevalued solution and established convergence results assuming that the diffusion dominates the dispersion. Subsequently, the approach of kinetic decomposition and velocity averaging [25] was introduced by Hwang and Tzavaras [11] to analyze singular limits including nonclassical shock waves. Furthermore, we also mention related works by Wu [32] and Jacobs, McKinney, and Shearer [12] which provides the first existence result of undercompressive shocks for the modified Kortewegde VriesBurgers equation.
It is well known that the relative scaling between and determines the limiting behavior of solutions, and we can distinguish between three cases:

Diffusiondominant regime : The qualitative behavior of solution is same as the solution of the conservation laws.

Dispersiondominant regime : In this case, high oscillations develop (as ) especially in regions of steep gradients of the solutions and only weak convergence is observed.

Balanced diffusiondispersion regime: This typically corresponds to the scenario where . Only mild oscillations are observed near shocks, and the limit solution is a weak solution to conservation laws. Most importantly, in this case, the solution exhibit nonclassical behavior, as they contain undercompressive shocks. However, for the regime , the limit solution coincides with the entropy solution determined by Kružkov theory [18].
1.2. Numerical Schemes
It is well known that standard finite difference, finite volume and finite element methods have been very successful in computing solutions to hyperbolic conservation laws with discontinuous coefficients, including those containing shock waves. However, we mention that most of these wellestablished methods are proven to be not good enough to capture nonclassical shock wave solutions numerically. This well known phenomena has been explained by many authers (Hou and LeFloch [10], Hayes and LeFloch [8], and others) in terms of the equivalent equation associated with discrete schemes through a formal Taylor expansion. The key idea behind capturing nonclassical shocks is to design finite difference schemes whose equivalent equation matches, both, the diffusive and the dispersive terms (cf. (1.3)) in the underlying model. However, these schemes fail to approximate nonclassical solutions with large amplitude, especially strong shocks due to lack of control on higher order error terms present in equivalent equation. A recent work by Ernest et al. [6] has overcome such problems by dominating higher order error terms in amplitude by the leading order terms of the equivalent equation.
In another paper by Chalons and Lefloch [2], the authors introduced a fully discrete scheme for the numerical approximation of diffusivedispersive hyperbolic conservation laws (cf. (1.1)(1.3)) in onespace dimension. An important feature of their scheme is that it satisfies a cell entropy inequality and, as a consequence, the space integral of the entropy is a decreasing function of time. Moreover, they showed that the limiting solutions generated by the scheme contains nonclassical undercompressive shock waves.
On the other hand, there is a sparsity of efficient numerical schemes for (1.1)(1.2) available in the literature. In fact, to the best of our knowledge, this is the first systematic attempt to construct a provably convergent numerical scheme for (1.1)(1.2). Having said this, there are some numerical experiments available in the final section of the recent paper by Coclite et al. [4] without rigorous analysis of the scheme.
1.3. Scope and Outline of the Paper
In view of the above discussion, it is fair to claim that there are no robust and provably stable numerical schemes currently available to simulate the vanishing capillarity approximations of scalar conservation laws equation (1.1)(1.2). In this context, we consider a fullydiscrete (in both space and time) finite difference scheme for (1.1)(1.2) which is provably convergent and able to capture non classical shocks quite well. Since diffusiondispersion model for the conservation laws with discontinuous flux has not been studied in detail, we analyze a fullydiscrete scheme for (1.1)(1.3) as well. While there are several numerical methods which perform well in practice, perhaps better than the one presented here, (see [20] for a recent comparison of diferent numerical methods) we emphasize that we prove the convergence of the schemes proposed in this paper. Here, we mention that a detailed analysis of the scheme introduced by Ernest et al. [6] is beyond the scope of this paper, and will be the topic of an upcoming paper.
To sum up, the schemes in the present paper have the following properties:

Approximate solutions for (1.1)(1.2), generated by the scheme (3.2), converge to the unique entropy solution of (1.4) as long as . A scheme (cf. (7.2)) has been formulated for (1.1)(1.3) and the same techniques can be applied, mutatis mutandis,to prove convergence of approximate solutions to the unique entropy solution of (1.4).
The rest of the paper is organized as follows: In section 2, we present the mathematical framework used in this paper. In particular, we have used a compensated compactness result in the spirit of Tartar [27] but the proof is based on divcurl lemma and does not rely on the Young measure. Section 3 introduces the fullydiscrete finite difference scheme for (1.1)(1.2) . In section 4, we derive a priori estimates for the approximate solutions and a detailed convergence analysis towards weak solutions of (1.4) has been discussed in section 5. Convergence towards the unique entropy solution has been considered in section 6, while a brief discussion on the results for diffusivedispersive approximation (1.1)(1.3) has been addressed in section 7. Finally, numerical results are presented in section 8 to illustrate the performance of the designed schemes.
2. Mathematical Framework
In this section, we list all the assumptions on the data for the problem (1.1), and present relevant mathematical tools to be used in the subsequent analysis. Throughout this paper we use the letters etc. to denote various generic constants independent of approximation parameters, which may change line to line, but the notation is kept unchanged so long as it does not impact the central idea.
The basic assumptions on the data of the problem (1.1) are as follows:

[label=A.0]

For the initial function , we assume that

For the discontinuous coefficient , we assume that

Regarding the flux function , we assume that

Furthermore, we assume that is genuinely nonlinear a.e. in , i.e., , for a.e. .
Remark 2.1.
It is worth mentioning that the Assumption 4 is typically required in the compensated compactness framework. This condition also imposes a condition on the coefficient . In fact, it implies that is genuinely nonlinear (i.e., ) and , for a.e. .
Next, we recapitulate the results required from the compensated compactness method due to Murat and Tartar [23, 27]. For a nice overview of applications of the compensated compactness method to hyperbolic conservation laws, we refer to Chen [3]. Let denote the space of bounded Radon measures on and
If , then
Recall that if and only if for all . We define the norm
The space is a Banach space and it is isometrically isomorphic to the dual space of . Furthermore, we define the space of probablity measures
Before we state the compensated compactness theorem, we shall recall the celebrated divcurl lemma [23].
Lemma 2.1 (divcurl lemma).
Let be a bounded open subset of . Suppose
in as . Furthermore, assume that the two sequences and lie in a (common) compact subset of , where and . Then along a subsequence
Suitably modified for our purpose, we shall use the following compensated compactness result. For a proof, we refer to the paper by Kenneth and Towers [13, Lemma 3.2].
Theorem 2.1.
Assume that 2, 3 and 4 hold. Let be a bounded open set, and assume that is a sequence of uniformly bounded functions such that , for all . Set
(2.1)  
where is an arbitrary constant. If the two sequences
belong to a compact subset of , then there exists a subsequence of that converges a.e. to a function .
We remark that, a feature of the compensated compactness result above is that it avoids the use of the Young measure by following an approach developed by Chen and Lu [3, 22] for the standard scalar conservation law. This is preferable as the fundamental theorem of Young measures applies most easily to functions that are continuous in all variables.
The following compactness interpolation result (known as Murat’s lemma [23]) is useful in obtaining the compactness needed in Theorem 2.1.
Lemma 2.2.
Let be a bounded open subset of . Suppose that the sequence of distributions is bounded in . Suppose also that
where is in a compact subset of and is in a bounded subset of . Then is in a compact subset of .
3. A FullyDiscrete Finite Difference Scheme
We begin by introducing some notation needed to define the fullydiscrete finite difference scheme. Throughout this paper, we reserve to denote small positive numbers that represent the spatial and temporal discretizations parameter of the numerical scheme. Given , we set for , to denote the spatial mesh points. Similarly, we set for , where for some fixed time horizon . Moreover, for any function admitting point values, we write . Furthermore, let us introduce the spatial and spatialtemporal grid cells
where . Let denote the discrete forward and backward differences in space, i.e.,
The discrete Leibnitz rule is given by
while the summationbyparts formula is given by
Furthermore, for any function , using the Taylor expansion on the sequence we obtain
for some between and . In other words, the discrete chain rule is accurate up to an error term of order .
Finally, let denote the discrete forward and backward difference operator in the time, i.e.,
The following identity is readily verified:
(3.1) 
We propose the following fullydiscrete (in space and time) finite difference scheme approximating the limiting solutions generated by the equation (1.1)(1.2)
(3.2)  
(3.3) 
where are fixed parameters, and as . More specifically, we will either use or depending on the quest for the convergence of approximate solution towards a weak solution or the entropy solution, respectively.
Remark 3.1.
Here we used the notation to denote quantities that depend on and are bounded above by , where is a constant independent of . Likewise, we used the notation to denote quantities that depend on and are bounded above by , where is a constant independent of and .
The numerical flux corresponding to the flux function is given by
where is based on a twopoint monotone numerical flux, i.e., nondecreasing with respect to the first argument and nonincreasing with respect to the second argument, and consistent with the actual flux, i.e., . Moreover, in order to maintain monotonicity of the scheme (3.2) without the higher order terms (corresponds to ), the arguments of the numerical flux are transposed when the coefficient is negative. More specifically, we choose
Summing up, the numerical flux is given by
In particular, we focus on EngquistOsher (EO) numerical flux given by
(3.4) 
Remark 3.2.
To this end, observe that the EO flux given by (3.4) is Lipschitz continuous . In fact, for , it has continuous partial derivatives satisfying
(3.5) 
using the conventional notations that and . It is also clear that serves as a Lipschitz constant for EO flux.
For a given initial data , we define the initial grid function by (3.3). Moreover, for the sequence , we associate the function defined by
where denotes the characteristic function of the set . Similarly, for the coefficient approximated at each cell boundary, we associate the function defined by
Note that and are discretized on grids that are staggered with respect to each other. This indeed results in a reduction in complexity, compared with the approach where two discretizations are aligned.
Throughout this paper, we use the notation to denote the functions associated with the sequence . For later use, recall that the discrete , and norms, and BV seminorm for a lattice function are defined respectively as
For the sake of simplified notations, unless specified, we shall use the notation to denote the discrete norm.
4. A Priori Estimates
This section is devoted to the derivation of a priori estimates which turns out to be useful to prove “strong compactness” of the approximate solution . To begin with, following Coclite et al. [4], we assume that the approximate solutions generated by the scheme are uniformly bounded, i.e., . In other words, we assume that

[label=B.0,resume]

For almost every , , for some fixed constants ;
Remark 4.1.
It is worth mentioning that the above assumption on the approximate solutions is the manifestation of the specific structure of the flux function (depends explicitly on space variable). In fact, to obtain bound on the solution, one requires a priori bound on the solution (cf. (4.4)). This assumption can be toppled by replacing the “space dependent flux function” to a flux function which depends explicitly only on the solution. In such a scenario, one can use framework of the compensated compactness result [5], to reproduce all the results in this paper.
To proceed further, we first collect all the available estimates on the approximate solutions in the following lemma.
Lemma 4.1.
Let be a sequence of approximations generated by the scheme (3.2). Moreover, assume that the initial data lies in . Then the following estimate holds
(4.1)  
provided and satisfies the following CFL condition
(4.2) 
with . Here and the constant is independent of .
In particular, the estimate (4.1) guarantees following spacetime estimates:
(4.3a)  
(4.3b)  
(4.3c)  
(4.3d) 
Remark 4.2.
In light of the CFL condition (4.2), we want to emphasize that if (required to prove convergence towards a weak solution, cf. Theorem 5.1), then we need to be kept fixed. On the other hand, if (required to prove convergence towards the entropy solution, cf. Theorem 6.1), then we need , with , to be kept fixed. To sum up, we need a stronger CFL condition to prove convergence of approximate solutions towards the unique entropy solution of (1.1). To this end, we mention that in the subsequent analysis the CFL condition (4.2) is always assumed to hold.
Proof.
To start with, we multiply the scheme (3.2) by and subsequently sum over . Then, using summationbyparts formula and the identity (3.1), we obtain
Note that, the identity (3.1) also implies that
Next, we move on to estimate the term . Using summationbyparts formula, we obtain
where is the primitive of , i.e., . A first order Taylor’s expansion together with the monotonicity of the numerical flux function gives us an estimate of . To see this, notice that
where lies in between and . To proceed further, we consider the following two cases:
Case 1: Assume that , then by the definition of numerical flux
. If , then
On the other hand, if , then
Thus, in any case we have .
Case 2: Assume that , i.e., . If , then
On the other hand, if , then
Therefore, in this case also, we have . Having this in mind and making use of the Assumption 1, we conclude that
(4.4) 
where is a constant independent of . Finally, combining all the above estimates, we obtain
(4.5)  
Next, we multiply the scheme (3.2) by and sum over to obtain,
(4.6)  
(4.7) 
Again, use of the identity (3.1) reveals that