ON SOME DIFFERENCE SCHEMES AND ENTROPY CONDITIONS FOR A CLASS OF MULTI-SPECIES KINEMATIC FLOW MODELS WITH DISCONTINUOUS FLUX

. We study a system of conservation laws that describes multispecies kinematic ﬂows with an emphasis on models of multiclass traﬃc ﬂow and of the creaming of oil-in-water dispersions. The ﬂux can have a spatial discontinuity which models abrupt changes of road surface conditions or of the cross-sectional area in a settling vessel. For this system, an entropy inequality is proposed that singles out a relevant solution at the interface. It is shown that “piecewise smooth” limit solutions generated by the semi-discrete version of a numerical scheme the authors recently proposed [R. B¨urger, A. Garc´ıa, K.H. Karlsen and J.D. Towers, J. Engrg. Math. 60:387–425, 2008] satisfy this entropy inequality. We present an improvement to this scheme by means of a special interface ﬂux that is activated only at a few grid points where the discontinuity is located. While an entropy inequality is established for the semi-discrete versions of the scheme only, numerical experiments support that the fully discrete scheme are equally entropy-admissible.


Introduction.
1.1. Scope. This paper concerns the following system of conservation laws that models the one-dimensional kinematic flow of a multi-species mixture: This system is to be solved in Π := {(x, t) | x ∈ R, t ≥ 0}, given the initial data φ i (x, 0) = p i (x), x ∈ R, i = 1, . . . , N. (2)

RAIMUND BÜRGER, KENNETH H. KARLSEN AND JOHN D. TOWERS
The unknowns are the densities φ i = φ i (x, t) of species i, i = 1, . . . , N . The scalar function k(x) is piecewise constant, and is allowed to have a jump at the origin: Our analysis applies to at least two real applications, namely a model of multiclass traffic flow on a highway with possibly heterogeneous road surface conditions (Model 1), and a model of the settling of an oil-in-water dispersion in a vessel with discontinuously varying cross-sectional area (Model 2).
In [10], we proposed simple difference schemes for (1)- (3). In this paper we analyze entropy conditions satisfied by limits of these schemes. We also present an improvement to the scheme by a straightforward modification of the numerical flux across the interface. We show that for discontinuities located away from the interface, limit solutions satisfy a standard integral entropy inequality, that if the system is genuinely nonlinear, any shocks are classical Lax shocks, and that limit solutions Φ satisfy the following interface jump entropy inequality: Here, φ − = φ − (t) and φ + = φ + (t) are the limit values of φ adjacent to x = 0, and φ * is the maximum of f on (0, 1) (precise statements are made later). We prove that for steady solutions, (8) can be obtained by simultaneously smoothing k(x), adding a regularizing diffusion term, and then letting the joint smoothing and regularization parameter tend to zero. In fact, we provide further support for (8) by proving that piecewise smooth limit solutions generated by the semi-discrete version of one of the schemes for (1) introduced in [10] converge to an entropy solution that satisfies (8) (provided that they do converge), and present numerical evidence that the fully discrete scheme produces solutions that satisfy the entropy jump condition to within an error that vanishes as the discretization parameters tend to zero.
1.2. Traffic and dispersion models. We discuss two typical real-world kinematic flow models that give rise to the Cauchy problem (1)- (3). Model 1 is based on the extension of the well-known Lighthill-Whitham-Richards (LWR) kinematic model to multiclass traffic. This extension, first proposed by Benzoni-Gavage and Colombo [7] and Wong and Wong [38] (see also [10,15,18,19,39,40,41]), considers N different species of vehicles, each having its own preferential velocity V i , i = 1, . . . , N , corresponding to a free highway. The function ψ = ψ(φ) describes how drivers adjust their velocities to the global car density φ(x, t) seen at a point (x, t). These assumptions lead to the velocity functions (6). The properties (7) reflect that on a free highway, cars travel at maximum speed (ψ(0) = 1), that there is a gradual increase of hindrance with increasing φ (ψ ′ (φ) < 0 on [0, 1)), and that there is no movement in the "bumper-to-bumper" situation (ψ(1) = 0). For the traffic model, the system (1) under the assumption (3) corresponds to a highway with a change of the road surface conditions at x = 0, which causes an abrupt change of the effective maximal velocity of species i to the left of x = 0, k L V i , to a different value k R V i to the right. We assume the jump in velocities to be "collective" in the sense that k(x) modulates all velocities V 1 , . . . , V N .
Model 2 describes the creaming of small oil droplets in an oil-in-water liquidliquid dispersion [32] with oil droplets of N different volumes (size classes) 0 < V 1 < V 2 < · · · < V N , where z denotes height and φ i = φ i (z, t) is the local volume fraction occupied by size class i. According to [32], the differential motion of the droplets is driven by differences in viscosity. The basic nonlinearity is a viscosity function µ d = µ d (Φ). If µ f denotes the viscosity of pure water, then µ d (Φ) is assumed to satisfy µ d (0, . . . , 0) = µ f , µ d ∈ C 1 (Ω), µ d (Φ) > 0, and ∂µ d /∂φ i > 0 for i = 1, . . . , N . The velocity functions v 1 (Φ), . . . , v N (Φ) are then given by where g, ̺ f and ̺ oil denote the acceleration of gravity, the density of pure water and density of pure oil, respectively. For [10]) the creaming of the dispersion can be described by the system of conservation laws where v i is given by (6) if we define . This application gives rise to a system of the type (1) if we consider a sufficiently large vertical vessel with the variable cross-sectional area S(z). The continuity equations that replace (10) are then Now, let us assume that there exists a constant S 0 > 0 such that S(z) ≥ S 0 throughout, and let us define the new spatial variable x = x(z) := (1/S 0 ) z 0 S(ζ) dζ. In terms of x and t, we may rewrite (11) as If S(z) is constant, but changes discontinuously at z = 0, i.e., S(z) = S L for z < 0 and S(z) = S R for z ≥ 0, then (12) turns into an equation of the type (1), (3) if we define k L := S L /S 0 and k R := S R /S 0 . The role of k(x) as a cross-sectional area produces here a natural "collective" modulation of all velocities.
1.3. Outline of the paper. The remainder of the paper is organized as follows.
In Section 2, we focus on the situation away from x = 0, where k is constant and (1) is a classical hyperbolic system of conservation laws. In Sections 2.1 and 2.2, we recall the definition of weak solutions to (1) and some known properties (see e.g. [39]) of the eigenvalues and eigenvectors of the Jacobian J f (Φ). In Section 2.3 we review the concept of entropy-admissible weak solutions to (1), which involves a convex entropy function E(Φ) with a corresponding entropy flux F (Φ), and is based on a vanishing viscosity approach. In Section 2.4, we review the Lax shock condition for genuinely nonlinear systems of conservation laws, which include the constant-k version of (1) for ψ(φ) = 1 − φ. In Section 2.5 we briefly comment on the relevance of the entropy concepts for our applications. Section 3 deals with the interface entropy jump condition for (1), (3). The basic calculus is presented in Section 3.1. The discussion is simplified under the assumption that φ → f (φ) is unimodal, i.e., has a single maximum. Motivated by our analysis of Model 1 in the scalar case [11], we then propose an interface entropy jump condition that is equivalent to stating that at least on one side of the interface, characteristics impinging on the t-axis must be traceable backwards to the x-axis. We argue that with the exception of one non-generic case, our condition is consistent with that of Zhang et al. [39], which involves one-sided limits of the smallest eigenvalue λ 1 = λ 1 (Φ) of J f (Φ). In Section 3.2 we apply a combined technique of smoothing k(x) and adding viscosity to (1) (the "SVV method") to derive (8). For Model 1 we provide additional justification that is independent of the notion of viscosity. To this end, we consider in Section 3.3 the case ψ(φ) = 1 − φ, for which (1) is genuinely nonlinear. We prove that if (8) rules out a two-state steady solution of (1)-(3) as inadmissible, then this solution is also ruled out by the "speedup impulse" (the condition introduced in [11]).
In Section 4 we consider numerical schemes for (1)-(3). Our analysis is based on a simple fully discrete finite difference scheme (called Scheme 1 herein) taken from [10]. A slight modification of the numerical flux near the interface to avoid spurious oscillations defines Scheme 2; the corresponding semi-discrete versions are termed Schemes 3 and 4. We prove that all schemes preserveΩ as an invariant region under a CFL condition. In Section 5 we prove that if Schemes 3 and 4 converge to a piecewise continuous weak solution Φ taking values inΩ, then Φ satisfies a standard integral entropy inequality on either side of x = 0. Section 6 extends this result to the complete computational domain Π. Precisely, it is shown that if Schemes 3 and 4 converge to a weak solution Φ, then Φ satisfies a standard integral entropy inequality for test functions whose supports do not intersect the interface and the entropy jump conditions described in Section 3.1.
In Section 7 we present numerical examples that illustrate the model, the consequences of choosing either k L < k R or k L > k R in combination with ψ(φ) = 1 − φ or ψ(φ) = (1 − φ) 2 , and the approximate satisfaction of the entropy jump condition by the fully discrete Scheme 2, where we recall that our analysis of Sections 5 and 6 is valid for the semi-discrete variants only (Schemes 3 and 4). Finally, in Section 8 we collect some conclusions of our analysis.
2. Preliminary remarks on the constant-k system.

Weak solutions.
For now, we focus on the behavior of solutions away from the interface. Since k(x) is constant in the region of interest, we simply denote it by the symbol k, and we will discuss the constant-k system This setup has been more extensively studied, and mostly fits within the framework of classical (strictly) hyperbolic systems of conservation laws. Even if the initial data are smooth, solutions of (1) develop discontinuities, and so we seek a weak solution, which is a bounded measurable function Φ satisfying for any smooth test function η(x, t) with compact support contained in Π.

2.3.
Entropy conditions. Consider a discontinuity along a smooth curve x = x(t) that does not intersect x = 0. If the solution Φ is continuous on either side of the curve, the weak formulation (14) implies the Rankine-Hugoniot (RH) condition Here s = dx/dt is the speed of the discontinuity, the shock speed, and Φ + and Φ − are the limiting values of Φ on the right and left sides of the jump, respectively. As is well known, weak solutions of systems of conservation laws are not unique. The RH condition (17) is not sufficient by itself to completely specify a solution. A number of admissibility criteria have been proposed, usually based on the physics being modeled by the system of conservation laws. Suppose that the system (13) admits a strictly convex entropy, meaning that there is a strictly convex function E(Φ) and an entropy flux F (Φ) such that ∇F (Φ) = ∇E(Φ)J f (Φ). We call a weak solution Φ of the constant-k system (13) entropy-admissible [9] if for every smooth nonnegative test function ρ with compact support contained in {(x, t)|x ∈ R, t > 0}. The entropy inequality (18) is justified by considering smooth solutions of the parabolic regularization Multiplying (19) by ∇E yields the scalar equation Letting ε → 0, and assuming that Φ ε converges boundedly a.e. to a limit Φ, we conclude (see e.g. [34]) that Φ satisfies (18) in the sense of distributions. In fact, (18) implies that all discontinuities satisfy the following entropy jump condition: A discontinuity satisfying (21) is called entropy-admissible.
Benzoni-Gavage and Colombo [7] showed that if φ 1 > 0, . . . , φ N > 0, then the constant-k system (1) is symmetrizable, and admits the convex entropy function E(Φ) and the associated entropy flux F (Φ) given by where Ψ is any primitive of ψ: Ψ ′ (z) = ψ(z). In [7] it is also shown that for ψ(z) = 1 − z, the entropy jump condition (21) is equivalent to the condition For both Models 1 and 2, (22) states that for an entropy-admissible shock, the total density (of the traffic or of the oil in the dispersion) to the left (i.e., of the flow approaching the shock) cannot exceed that to the right (i.e., behind the shock).

2.4.
Lax shock condition. In many applications, discontinuities are required to satisfy the Lax shock condition [26], which in addition to satisfaction of the RH condition (17) requires that there should be an index m ∈ {1, . . . , N } for which Such a discontinuity is Lax-admissible, and is called a Lax shock or an m-shock. Condition (23) results from adding to (1) a regularizing dissipative mechanism. The physically relevant solution is the one that results from this regularized system in the limit as the dissipation vanishes. Condition (23) is generally applicable to genuinely nonlinear fields. In the genuinely nonlinear setting, entropy admissibility and Lax admissibility are equivalent for shocks where |Φ − −Φ + | is sufficiently small.

Entropy jump condition in applications.
For a scalar equation (N = 1), condition (22) seems to be generally accepted in the traffic modeling literature, and is called the driver's ride impulse [4]. It is not obvious whether (22) is still relevant for N ≥ 2, but one might argue as follows. In the scalar case, (22) requires that cars slow down as they pass through the shock (since ψ is a decreasing function of φ).
In the system case, it seems reasonable to require that all classes of vehicles slow down as they pass through the shock. Note that in the system case, the velocities are given by kV i ψ(φ), and so (22) is equivalent to this requirement. Figure 1. The fluxes for the associated scalar conservation law.

3.1.
Interface entropy jump condition. Due to the jump in k at x = 0 (the interface), there will generally be a stationary (s = 0) discontinuity located there. By considering in (14) a test function η with support near the interface, one finds that the RH condition at the interface is If we cancel V i from both sides of each equation, and then add the N equations, we arrive at the following combined RH condition: Thus, (25) is the same RH condition that occurs at the interface for the scalar conservation law (5). In what follows, we assume that ψ(φ) is chosen such that the flux f (φ) = φψ(φ) is unimodal, by which we understand that f satisfies f (φ) > 0 for φ ∈ (0, 1) and that φ → f (φ) has a single maximum φ * ∈ (0, 1). It follows that f is strictly increasing on (0, φ * ) and strictly decreasing on (φ * , 1). The function f is unimodal, for example, if ψ is of the form ψ(φ) = (1 − φ) ν with ν ≥ 1. Note that f is not required to be concave. However, in order for the system to be genuinely nonlinear, we will generally have to assume that ψ(φ) = 1 − φ. When this genuine nonlinearity assumption is in effect, we will state it explicitly.
In [11] we argued that in the scalar case, a generalization of the condition (22) that is meaningful for traffic modeling at the interface is which says that on at least one side of the interface, the characteristics impinging on the t-axis must be traceable backwards in time to the x-axis. If f is unimodal, then the interface condition (26) is equivalent to the condition In the context of the traffic model (Model 1), we referred to our rationale as the speedup impulse [11], which simply states that all cars will increase their speed whenever possible. It seems reasonable to assume that the speedup impulse remains valid in the system setting, and so we take (27) as the analog of (22) that should hold at the interface. For a unimodal flux f , Zhang et al. [39] showed that sgn(φ * − φ) = sgn(λ 1 (Φ)). Thus, our entropy jump condition (27) is equivalent to Note that Zhang et al. [39] use the slightly more restrictive interface condition The one situation that is allowed by our entropy condition that is not allowed by (29) is the one where At least in the scalar case (N = 1), this type of solution is in some sense nongeneric, meaning that it would be unlikely to be observed in an actual traffic flow (see Remark 7 of [11]). For this reason we view our results as consistent with those of [39].
Proof. The equivalence of conditions (26), (27) and (28) follows from the discussion preceding the lemma. It suffices to show that (30) is equivalent to (27). Take the case where k R > k L ; the other case is similar, and we omit it. To show that (30) implies (27), it suffices to show that if (30) is satisfied, we cannot have If (31) holds, then for the left-hand side of (30) we have (where the first equality follows from (25)), which is a contradiction.
To show that (27) implies (30), note that if (27) holds, then at least one of the inequalities φ − ≤ φ * and φ + ≥ φ * is satisfied. Take the case where both are satisfied.
(where the second equality follows from (25)), and thus clearly (30) holds again. If φ − < φ * and φ + < φ * , then where the last equality follows from (25)), and (30) is satisfied. There are two subcases that remain to be checked: The required calculations are similar to those above, and we omit them.

3.2.
Justification of the entropy inequality (30) by smoothing k(x) and vanishing viscosity. In this section, we consider a regularized version of the system (4). Specifically, we add a small viscous term, and simultaneously smooth the discontinuous coefficient k(x). We show that steady state solutions of the system (4) that are obtained as the limit of steady state solutions of the regularized system satisfy the interface entropy condition (30).
To this end, let k δ (x) := (1 − H δ (x))k L + H δ (x)k R denote a regularized approximation to k(x) for δ > 0, where H δ (x) is a C 1 regularization of the Heaviside function H(x) such that dH δ (x)/dx ≥ 0 and H δ → H boundedly a.e.. We then regularize the system (1) by substituting this smoothed version of k(x), and simultaneously including a small artificial viscosity term: In what follows we will employ the following Kružkov entropy/entropy flux pair: Lemma 3.2. Assume that f is unimodal, and that Φ is a steady (Φ t = 0) weak solution of the system (1) which is the limit, boundedly a.e., of steady smooth solutions Φ δ of the regularized system (32). Finally, assume that Then Φ satisfies the interface entropy inequality (30).
Proof. Dividing the ith equation of (32) by V i , then adding all equations, and keeping in mind that Φ t = 0, gives For Moreover, for the left-hand side of this inequality, we compute Combining (35) and (36) yields Now let η : R → R be a smooth nonnegative test function with compact support. Multiplying both sides of (37) by η and integrating by parts gives Letting ε → 0, we obtain by the dominated convergence theorem Next, we apply the identity Another application of the dominated convergence theorem as δ → 0 results in A standard test function argument, which we omit, now completes the proof.

3.3.
Interface entropy condition for the traffic model (Model 1) justified via speedup impulse. Lemma 3.2 partially justifies the entropy inequality (8). For the traffic model, however, we also would like to have a justification that does not rely on artificial viscosity, but is rather related to driver/vehicle behavior. In [11], we proposed a criterion, the speedup impulse, to justify our entropy jump conditions for the scalar version (5) of (1). The speedup impulse says that, if possible, the solution will evolve to one with a higher velocity on the left side of the interface. This is explained by the desire of individual drivers to increase their speeds, or equivalently to reach their destinations as quickly as possible. If (1) is genuinely nonlinear, (i.e., we use ψ(φ) = 1 − φ), then the speedup impulse can also be used in the present context. We concentrate on the admissibility for two-state solutions of the form (1) is genuinely nonlinear. If the solution (38) is inadmissible under any of the equivalent entropy conditions (26), (27), (28), (30), and Φ 0 − , Φ 0 + ∈ Ω, then there is a similarity solution Φ 1 (x, t) of (4) such that Φ 1 (x, 0) = Φ 0 (x) and for t > 0, φ 1 (0−, t) < φ 0 − . Moreover, any jumps in Φ 1 (x, t) away from the interface are classical Lax shocks of the form (23).

Remark 1.
In [11] we originally formulated the "speedup impulse" as follows: A two-state solution like (38) has "speedup potential" if when this two-state solution is used as initial data, it can evolve (under the conservation law) to a new solution where away from the interface all discontinuities satisfy the usual entropy condition, while at the left of the interface the velocity (density) is larger (smaller) than in the initial data. The principle of "speedup impulse" says that no two-state solution should have "speedup potential". In the scalar case (with concave flux) we showed that our entropy principle is the same as the "speedup impulse". More specifically, a two-state solution satisfies our entropy condition if and only if it has no "speedup potential". In the present paper, we show it in one direction. Namely, that all two-state solutions ruled out by our entropy principle have "speedup potential". In other words, such solutions are ruled out by the "speedup impulse" principle. More precisely, Lemma 3.3 says that if the solution (38) is inadmissible under our entropy condition, then it is also ruled out by the speedup impulse. Specifically, the speedup impulse will cause (38) to evolve into a solution with higher velocity for cars entering the interface. Here we are using that the speed of the ith class entering the interface is k L V i ψ(φ − ), and ψ is strictly decreasing, so a decrease in φ is equivalent to an increase in speed for every class. We state our speedup criterion in terms of speed rather than φ because relative speed is more readily observable by a driver than density.

4.
The difference schemes. Let ∆x > 0 denote the spatial discretization parameter. The spatial domain R is divided into cells I j := [x j − ∆x/2, x j + ∆x/2) with centers at the points x j = j∆x for j ∈ Z. Let χ j (x) be the characteristic function for the interval I j . We use the symbols Φ j (t) and φ i,j (t) to denote the approximate solution at x j : We will often suppress the dependence on t. We discretize (2) according to If p(x) = (p 1 (x), . . . , p N (x)) T ∈Ω for all x ∈ R, then Φ 0 j ∈Ω for all j ∈ Z. We discretize k(x) by k j = k(x − j ), so that k j = k L for j ≤ 0 and k j = k R for j ≥ 1. There are two versions of the scheme for variable k, denoted Scheme 1 and Scheme 2. Scheme 1 has no special processing for the interface (this is the original scheme appearing in [10]). Scheme 2 has a special interface flux, which mostly gets rid of spurious bumps that appear for certain initial data. We refer to the semi-discrete version of Scheme 1 as Scheme 3, i.e., We are using the dot to denote d/dt. The flux in (41) is defined by Like the Lax-Friedrichs scheme, our algorithms do not require knowledge of the eigenvectors and eigenvalues of the system. This is of particular interest for the system (1) because computing the eigenvectors and eigenvalues is not straightforward (but see [18]). Our numerical experiments in [10] indicate that Scheme 1 is robust, and is easily extended to a formally second order version via a standard MUSCL/Runge-Kutta extension (which we, however, do not pursue herein).
In certain situations, Scheme 3 (defined by (41)) causes a spurious overshoot at the interface. We dealt with this in [11] for the scalar version of the problem by limiting the magnitude of the flux at the interface. Our modified flux for (1) is a generalization of that idea. Let f max = max w∈[0,1] f (w). Due to our discretization of k(x), the interface is located between the grid points x 0 and x 1 . If we are at the interface, the components of the unmodified numerical flux are The modified interface flux is defined by Schemes 2 and 4 are the versions of the discrete Scheme 1 and the semi-discrete Scheme 3, respectively, defined by using (43) instead of (42) for j = 0. With the definition (43), it is easy to check that To motivate this condition, we note that the combined RH condition (25) requires a similar inequality for piecewise continuous solutions of (1). Indeed, the quantity on the left side of (44) should approximate k L f (φ − ) and k R f (φ + ), and these quantities satisfy Our entropy results for the two variants of the difference scheme are valid for the semi-discrete setup above (Schemes 3 and 4), but we will also consider the fully discrete versions (Schemes 1 and 2). To this end, time is discretized via t n = n∆t for n ∈ Z 0 + := {0, 1, . . .} (Z + := {1, 2, . . .}), resulting in time strips [t n , t n+1 ). We use the notation Φ n j and φ n i,j to denote the fully discrete approximations: With λ = ∆t/∆x, the fully discrete versions, Schemes 1 and 2, take the form We always assume that the CFL condition is in effect for the fully discrete Schemes 1 and 2. Note that the CFL condition does not require explicit knowledge of the eigenvalues of the Jacobian.
The following theorem establishes that our difference schemes preserve the invariant regionΩ. For the original scheme (without the interface flux), this result basically follows from Theorem 3.1 of [10]. The novelty is that this invariance also holds for the version that incorporates the interface fix.
Theorem 4.1. The fully discrete scheme (46), with or without the interface fix (Schemes 1 and 2, respectively), preserves the invariant regionΩ. Similarly, the semi-discrete scheme (41), with or without the interface fix (Schemes 3 and 4, respectively), preserves the invariant regionΩ. In other words, if the initial data lies in the invariant region, Φ 0 (x) ∈Ω all x ∈ R, then for Schemes 1 and 2, at each time level n ≥ 0 we will have Φ n j ∈Ω for all j ∈ Z. For Schemes 3 and 4, for each t ≥ 0, we will have Φ j (t) ∈Ω for all j ∈ Z.
Proof. For our specific numerical fluxes, we may rewrite Schemes 1 and 2 as Here θ n i,j+1/2 ∈ [0, 1] is a parameter that accounts for whether or not the interface fix is in effect. Suppose that Φ n j ∈Ω for all j. Then (48) implies that In view of (47), we obtain from the last inequality that φ n+1 i,j ≥ 0. To show that φ n+1 j ≤ 1, note that from (48) we have Summing this over i and using the assumption that φ n j−1 ≤ 1, we obtain φ n+1 . Since ψ(1) = 0, we have that G(1) = 1, and that G ′ (φ n j ) = 1 + λ max{k L , k R }V N ψ ′ (φ n j ). Thanks to the CFL condition (47), 1]. Combining this with G(1) = 1, we have G(φ n j ) ≤ 1, and thus φ n+1 j = G(φ n j ) ≤ 1. We have proven that both variants of the fully discrete scheme preserve the invariant region.
For the semi-discrete variants, Schemes 3 and 4, we keep ∆x fixed, and let ∆t ↓ 0 in the fully discrete variants (Schemes 1 and 2). The approximations converge to the solutions of the system of differential equations (41) defining the semi-discrete approximations. This follows from a standard fact concerning the convergence of Euler's method when approximating systems of ordinary differential equations [21]. This result is applicable here since the numerical flux (both with and without the interface fix) is Lipschitz continuous in all arguments. The limit solutions (the semi-discrete approximations) then also remain inΩ.

5.
Discrete entropy inequality away from the interface. In this section we show that for grid points away from the interface, the semi-discrete difference scheme satisfies a discrete entropy inequality, which along with the conservation form of the scheme guarantees that if the scheme converges to a piecewise continuous solution, any discontinuities that do not intersect the interface are entropy-admissible. This means that in the genuinely nonlinear case, i.e., ψ(φ) = 1 − φ, all discontinuities with |Φ − − Φ + | sufficiently small are Lax shocks.
We concentrate on the scheme (41), (42) with j = 0, 1 (i.e., away from the interface). Thus, the original scheme coincides with the one with the interface flux, (43). Also note that k j = k j+1 since we are away from the interface.
Let D + denote the divided spatial difference operator, i.e., D + z j = (z j+1 − z j )/∆x. Employing an identity by Osher [28], one can readily show that for Φ j , Φ j+1 ∈ Ω, our semi-discrete scheme satisfies Here Γ j+1/2 is any piecewise smooth curve in Ω connecting Φ j to Φ j+1 , andF j−1/2 is the numerical entropy flux Following [28], it is straightforward to prove for any vectors Φ j , Φ j+1 ∈ Ω that the right-hand side of (49) is non-positive, and thereby to prove that the semi-discrete scheme satisfies the discrete entropy inequalitẏ The proof of [28] involves the evaluation of ∇E(Φ) and ∇ 2 E(Φ), which are not defined for vectors with zero components. However, below we prove by an approximation argument that (50) remains valid for vectors Φ j , Φ j+1 ∈Ω.
Then the semi-discrete scheme satisfies the discrete entropy inequality (50).
Summarizing, we have shown that the semi-discrete scheme satisfies the approximate discrete entropy inequality Taking the limit µ ↓ 0 proves the result.
Theorem 5.2. Assume that f is unimodal, and that the approximate solutions generated by the semi-discrete scheme (Schemes 3 and 4) converge boundedly a.e. to a piecewise continuous weak solution Φ ∈Ω. If the limit solution Φ has a discontinuity described by a smooth curve x = x(t) that does not intersect the interface x = 0, then that discontinuity is entropy-admissible, i.e., the entropy jump condition (21) is satisfied. If in addition, we choose ψ(φ) = 1 − φ, so that the system is genuinely nonlinear, then the discontinuity is a Lax shock, i.e., the Lax condition (23) for shocks is satisfied, as well as the condition (22) that the density increases across the shock.
Proof. Recall that the approximate (and thus limit) solutions remain inΩ, consult Theorem 4.1. Let ρ ≥ 0 be a smooth test function with compact support in R × R + . Assume that the support of ρ does not intersect the interface x = 0. We multiply the discrete entropy inequality (5.1) by ρ(x j , t), sum over j and integrate over t > 0, and then proceed as in the proof of the Lax-Wendroff theorem to conclude that the limit solution satisfies the entropy inequality (18). A standard test function argument then yields the jump condition (21). The remaining assertions of the proof follow from the discussion in Section 2.4.

6.
Discrete entropy inequality at the interface.
Theorem 6.1. Assume that f is unimodal, and that the approximate solutions generated by the semi-discrete schemes (Schemes 3 and 4) converge boundedly a.e. to a piecewise continuous weak solution Φ. Fix a time t 0 > 0, and assume that in a neighborhood N α := (−α, α) × (t 0 − α, t 0 + α) of the point (0, t 0 ), the solution is piecewise continuous, with the only discontinuity in N α being a single jump along the line x = 0. Letφ ∆ i := j∈Z χ j (x)φ i,j (t), where χ j is the characteristic function for I j , and assume thatφ Then at t 0 the limit solution Φ satisfies the equivalent entropy jump conditions (26), (27), (28), and (30). In particular, this result applies when Φ is a steady (Φ t = 0) weak solution of the system.
Proof. Recall that the approximate (and thus limit) solutions remain inΩ, consult Theorem 4.1. For now, assume that we are dealing with Scheme 3 (the semidiscrete scheme without the special interface flux). As in the proof of Lemma 5.1, we adapt the semi-discrete entropy inequality of Osher [28].
In either of these last two cases, (67) follows immediately from either of the two possibilities that h int 1/2 = h 1/2 or h int 1/2 = min{k L , k R }f (φ * ). We now have (for both Scheme 3 and Scheme 4) Let ρ ≥ 0 be a test function with supp(ρ) ⊆ N α . Multiplying both sides of (72) by ρ(x j , t), and then rearranging, we get Let ρ ∆ (x, t) = j∈Z χ j (x)ρ(x j , t). We sum (73) over j, integrate over t > 0, and use that k j+1 = k j for j = 0 to get On the left-hand side, we sum by parts, which yields By applying the dominated convergence theorem, along with the assumption (63), we obtain If we now substitute into (75) test functions of the form ρ ε (x, t) = η(t)σ ε (x) satisfying η(t) ≥ 0, σ ε (x) ≥ 0, σ ε (x) = 0 for |x| > ε, σ ε (0) = 1 and η(t) = 0 for |t − t 0 | > α and then take the limit as ε ↓ 0, the second integral on the right side of (75) vanishes since ∂ t z ∈ L 1 (N α ), and we obtain Fixing t = t 0 , and dividing both sides by η(t 0 ), we have the entropy inequality (30), and thus also (26), (27) and (28). To put the effect of a discontinuous coefficient k(x) (shown in Examples 2 to 5) into the proper context, we first present an example for the constant-k case, i.e., k ≡ 1, and ψ(φ) = 1 − φ. Figure 2 shows the numerical solution. The discontinuity initially located at x = −1 evolves into a series of kinematic shocks. The t = 8  Table 1. Indicator d(∆t) of approximate satisfaction of the interface entropy jump condition. The final time is T = 40. plot of Figure 2 shows that from the left to the right, a zone forms in which only the slowest species 1 is present, followed by a zone containing only species 1 and 2, and so on. The remaining plots show that species 1 is eventually left behind. Of course, the other species will equally segregate as a consequence of the ordering V 1 < V 2 < · · · < V N . Note that all discontinuities produced by the numerical solution, appearing here as vertical segments of profiles, are associated with a sudden increase of φ, in agreement with the entropy jump condition (22).
In Examples 2 and 3 we choose ψ(φ) = 1 − φ. In Example 2 we set k L = 1 and k R = 0.5, and in Example 3, k L = 0.5 and k R = 1. Figure 3 shows the corresponding numerical solutions obtained by Scheme 2. Note that there is a decreasing stationary jump at x = 0. For Examples 2 and 3 we illustrate in Figure 4 the convergence properties of the scheme. Figure 4 shows portions of the total density φ of the numerical solution with ∆x = 1/2000 for Examples 2 and 3, along with numerical values for ∆x = 1/J with J = 100, J = 200 and J = 400. We observe that standard traveling discontinuities are approximated in "smeared" fashion, while the stationary discontinuity located at x = 0 is sharply resolved. In Examples 4 and 5, we choose k L and k R as in Examples 2 and 3, respectively, but now utilize ψ(φ) = (1 − φ) 2 . Figure 5 shows the numerical solutions. While the overall behaviour is similar to that of Examples 2 and 3, we observe some features that are not possible with ψ(φ) = 1 − φ; for example, the first two plots of Figure 5 indicate a kinematic shock to the left of x = 0 that is moving to the left.
Furthermore, for Examples 2 to 5 we also test the approximate satisfaction of the interface jump entropy condition (8). To this end, we fix λ and determine the  8. Conclusions. The present analysis is based on a system of conservation laws with discontinuous flux in the very particular form (1), (3), which is, however, relevant to several real-world applications. In particular, these applications give rise to systems with arbitrarily large numbers of species N . The structural properties of the flux vector and the fact that the nonlinearity is defined by the single scalar function ψ = ψ(φ) admit an entropy criterion (30) that only appeals to limits of φ, but not to the composition Φ. This property is well established for standard shocks (those that occur in the constant-k case).
Inequality (8) is essentially consistent with the jump condition by Zhang et al. [39]. In that paper it is mentioned that it is usually assumed that f (ϕ) is concave (f ′′ (ϕ) < 0 for 0 ≤ ϕ ≤ 1), but the properties that are indeed required for the analysis are satisfied by functions f that are just unimodal. Since we focus on the properties of difference schemes, our analysis complements the results of [39,40], which are based on the discussion of eigenvalues, characteristics and waves.
Concerning the convergence analysis, condition (63) may seem rather restrictive. However, from the point of view of computations, it is not unusual to observe approximations that are consistent with this condition.
For Model 1 we do not consider individual jumps in the maximum velocity of each species, but this should provide an interesting extension of the present analysis. Likewise, we have started to analyze systems of conservation laws modelling the sedimentation of polydisperse suspensions (see e.g. [8,10,15]), which are similar to (1) with the same motivation of a multiplicative flux discontinuity as that given in Section 1.2 for Model 2, but where v i depends on Φ in a more involved way that has so far precluded the determination of an entropy-entropy flux pair.