<figure><img src="image_1.png" /></figure>

HAL
open science

# Existence result for a 2 x 2 system of conservation laws with discontinuous flux and applications

F A Chiarello, S Fagioli, M D Rosini

► To cite this version:

F A Chiarello, S Fagioli, M D Rosini. Existence result for a 2 x 2 system of conservation laws with discontinuous flux and applications. 2024. hal-04798006

HAL Id: hal-04798006

https://hal.science/hal-04798006v1

Preprint submitted on 22 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Existence result for a 2 × 2 system of conservation laws with discontinuous flux and applications

Felisia Angela Chiarello, Simone Fagioli and Massimiliano Daniele Rosini

**Abstract.** This paper analyzes a two-by-two Temple-type system of conservation laws with discontinuous flux, focusing on applications in traffic modeling. We prove the existence of entropy solutions for initial data with sufficiently small total variation. Additionally, we explicitly construct the corresponding Riemann solver and examine its key properties. To illustrate the impact of the discontinuous flux on the solution, we also present numerical simulations.

**Mathematics Subject Classification (2010).** Primary: 35L65; 35B40; Secondary: 35D30; 35Q70; 35L45.

**Keywords.** 2 × 2 systems of conservation law, discontinuous flux, existence result, entropy condition, vehicular traffic, Aw-Rascle-Zhang model, Chiarello-Göttlich-Schillinger-Tosin model, sedimentation model.

## 1. Introduction

Conservation laws with discontinuous flux functions arise in various real-world applications. These applications span multiple fields, including traffic flow, multi-phase fluid dynamics (e.g., two-phase flow models and environmental engineering), material science (e.g., crystallization processes and composite materials), and systems engineering (e.g., acoustic waves and electromagnetic waves). In particular, traffic systems often exhibit sudden changes in road conditions, speed limits, or traffic regulations, which can be effectively modelled using conservation laws with discontinuous flux.

In this paper we study the one-dimensional 2 × 2 system of conservation laws with discontinuous flux

$$ \begin{cases} \partial_t \rho + \partial_x (cV(h) \rho) = 0, \\ \partial_t (\rho(h+p(\rho))) + \partial_x (cV(h) \rho(h+p(\rho))) = 0. \end{cases} \qquad (1) $$

This system finds its relevance in the modeling of vehicular traffic. Indeed, it can be interpreted as a generalisation of two second order vehicular traffic models: the one proposed by Aw, Rascle [5] and Zhang [43], and the one kinetically derived by Chiarello, Göttlich, Schillinger and Tosin in [14]. We refer to the former as the ARZ model, and to the latter as the CGST model. In this context, $t \ge 0$ denotes the time, $x \in \mathbb{R}$ the space, $\rho = \rho(t, x) \ge 0$ the density, $c = c(x) > 0$ is a given discontinuous function that mimics the changes in the road conditions, $p = p(\rho) \ge 0$ is the “pressure” function, and $V \ge 0$ represents the speed law. The interpretation of $h = h(t, x) \ge 0$ depends on the specific model under consideration: it represents the velocity in the ARZ model and the mean headway in the CGST model, see Section 2.2 for further details. It is worth mentioning that in [11], the authors kinetically derive a second order macroscopic traffic model, similar to the CGST model, but with the assumption of a constant headway.

The majority of theoretical advancements concerning conservation laws with discontinuous coefficients are limited to scalar equations, see [3, 35] and the references therein. In particular, in [31] the
authors provide Kružkov-type entropy condition for nonlinear degenerate parabolic equations with discontinuous coefficients, establishing $L^1$-stability, and thus uniqueness, for weak solutions satisfying the entropy condition, provided that the flux function satisfies the so called “crossing condition”. In [32] a convergence proof for the Lax-Friedrichs finite difference scheme for non-convex genuinely non-linear scalar conservation laws with discontinuous coefficients is given. The development of analogous theories for nonlinear systems with discontinuous fluxes, even in the one-dimensional space, remains an open area for further investigation. An exception to this is found in the studies presented in [29, 30], which explore compressible Euler flows within a pipe, in [37] where the ARZ model with seven different types of discontinuous fluxes are analyzed, and in [39] where Riemann solvers are proposed for $2 \times 2$ systems of conservation laws with discontinuous flux, referring to the ARZ model or modeling polymer flooding. Moreover, there are results on networks. For example, in [1] the authors analyze control problems considering a general class of junction distribution controls and inflow controls, establishing the compactness in $L^1$ of a class of flux-traces of solutions. Furthermore, in [16] a single conservation law, defined on a road network that is a collection of roads with junctions, is considered. The authors choose some fixed rules for the distribution of traffic plus optimization criteria for the flux, proving existence of solutions to the Cauchy problem and the Lipschitz continuous dependence by initial data under special assumptions. A generalization can be found in [15].

System (1) falls in the class of Temple systems [41], characterized by coinciding shock and rarefaction Lax curves for the genuinely nonlinear characteristic fields. In fact, only the first characteristic field of (1) is genuinely nonlinear and satisfies this condition, while the second characteristic field is linearly degenerate. Specifically, along the first Lax curve, the solution can exhibit shock waves, rarefaction waves, or mixed waves. In contrast, along the second Lax curve, the solution can contain only contact discontinuities.

The paper is organised as follows. In Section 2, we present our existence result in Theorem 2.1, which applies away from the vacuum. In the following three remarks, we discuss alternative strategies that could potentially extend the existence results to cases including the vacuum. The final remark addresses the entropy condition, which is further examined in details in Section 4. In Section 2.2, we compare some models from the literature with (6), highlighting its possible applications. In Section 2.3, we analyse the eigenvalues and eigenvectors of our system and conduct a study of the Lax curves. In Section 3, we introduce explicitly a Riemann solver and study its main properties. In Section 4 the entropy condition for system (1) is deduced and analysed. As far as we know there is no result in the literature regarding the entropy condition in the case of systems of conservation laws with discontinuous fluxes. In Section 5 some numerical simulations are performed to corroborate the analytical results obtained in the previous Sections. Finally, we defer all technical proofs to Section 6.

## 2. Existence result and applications

In this section, we present our existence result in Theorem 2.1 along with some applications and properties of (1).
## 2.1. Main result

Theorem 2.1 holds away from the vacuum, i.e., for $\rho > 0$, and under the following set of assumptions: for any $\rho > 0$ and $h > 0$ we have

$$p \in \mathbf{C}^1([0, \infty)), \quad p(0) = 0, \quad p'(\rho) > 0, \qquad (2)$$

$$V \in \mathbf{C}^2([0, \infty)), \quad V(0) = 0, \quad V'(h) > 0, \qquad (3)$$

$$V(h) \neq \rho p'(\rho) V'(h), \quad \rho p'(\rho)^2 V''(h) < (2p'(\rho) + \rho p''(\rho)) V'(h), \qquad (4)$$

and there exist $c_{\min}, c_{\infty} > 0$ and $\xi_i, i \in [1, p] := [1, p] \cap \mathbb{Z}$, with $\xi_i < \xi_{i+1}$ such that for almost any $x \in \mathbb{R}$

$$c \in \mathbf{L}_{\text{loc}}^{1}(\mathbb{R}) \cap \mathbf{C}^{2}(\mathbb{R} \setminus \{\xi_{1}, \dots, \xi_{p}\}), \quad c - c_{\infty} \in \mathbf{L}^{1}(\mathbb{R}), \quad c(x) \geq c_{\min}. \qquad (5)$$

Note that $p(\rho) > 0$ for any $\rho > 0$ by (2)$_{2,3}$ and that $V(h) > 0$ for any $h > 0$ by (3)$_{2,3}$.

For $\rho > 0$ system (1) is equivalent to the system of conservation laws

$$\partial_t \mathbf{U} + \partial_x (c F(\mathbf{U})) = \mathbf{0}, \qquad (6)$$

with

$$\mathbf{U} := (\rho, q)^T, \quad F(\mathbf{U}) := V(h(\mathbf{U})) \mathbf{U}, \quad h(\mathbf{U}) := \frac{q}{\rho} - p(\rho), \quad \mathbf{0} := (0, 0)^T, \qquad (7)$$

where $q$ is the generalized momentum and satisfies

$$q = \rho (h + p(\rho)). \qquad (8)$$

The dependent variable $\mathbf{U} = \mathbf{U}(t, x)$ takes values in

$$\mathcal{U} := \{(\rho, q)^T \in \mathbb{R}^2 : \rho \geq 0, q \geq \rho p(\rho) \geq 0\}.$$

Observe that in $\mathcal{U}$ the vacuum $\rho = 0$ corresponds to $\mathbf{U} = \mathbf{0}$ by (8). In this section, even if we do not state it explicitly, we consider states $\mathbf{U}$ away from the vacuum, i.e., states $\mathbf{U}$ in

$$\mathcal{U}_m := \mathcal{U} \setminus \{\mathbf{0}\}. \qquad (9)$$

In the next theorem we give our existence result for the Cauchy problem for the system of conservation laws (6) in

$$\mathcal{U}_e := \{(\rho, q)^T \in \mathbb{R}^2 : \rho > 0, q > \rho p(\rho) > 0\}. \qquad (10)$$

Notice that $\mathcal{U}_e \subset \mathcal{U}_m$ by (9), hence our existence result holds away from the vacuum. The proof is based on [12, Theorem 7.1]. It is worth noticing that we cannot apply [12, Theorem 7.1] directly to (6) because its flux function lacks the necessary regularity since $c$ is discontinuous at $\xi_i, i \in [1, p]$, see (5). To address this issue, as in [39], we introduce a $3 \times 3$ system of conservation laws, see (43).

**Theorem 2.1.** Fix $(\rho_{\infty}, q_{\infty})^T \in \mathcal{U}_e$. Assume that $p, V, c$ are under assumptions (2)-(5). Then, there exists $\delta_0 > TV[c]$ with the following property: if $\bar{\mathbf{U}} := (\bar{\rho}, \bar{q})$ takes values in $\mathcal{U}_e$, $TV[\bar{\mathbf{U}}] \leq \delta_0 - TV[c]$ and $\bar{\rho} - \rho_{\infty}, \bar{q} - q_{\infty} \in \mathbf{L}^1(\mathbb{R})$, then the Cauchy problem

$$\left\{ \begin{array}{l} \partial_t \mathbf{U} + \partial_x (c F(\mathbf{U})) = \mathbf{0} \\ \mathbf{U}(0, x) = \bar{\mathbf{U}}(x) \end{array} \right. \qquad (11)$$

has a weak solution $\mathbf{U} = \mathbf{U}(t, x)$ defined for all $t \geq 0$, that is $\mathbf{U}(0, x) = \bar{\mathbf{U}}(x)$ for a.e. $x \in \mathbb{R}$ and

$$\int_{\mathbb{R}_+} \int_{\mathbb{R}} (\mathbf{U} \varphi_t + c F(\mathbf{U}) \varphi_x) \, dx \, dt = 0$$
for any test function $\varphi \in \mathbf{C}_c^\infty((0, \infty) \times \mathbb{R}; [0, \infty))$.

In addition, if the system of conservation laws in (11) admits a convex entropy pair $(\mathcal{E}, c\mathcal{Q})$ such that

$$ \mathcal{Q}(\mathbf{U}) = D\mathcal{E}(\mathbf{U}) \cdot F(\mathbf{U}), \qquad (12) $$

then one can find a weak solution $\mathbf{U}$ that satisfies the entropy condition

$$ \int_{\mathbb{R}_+} \int_{\mathbb{R}} (\mathcal{E}(\mathbf{U}) \varphi_t + c \mathcal{Q}(\mathbf{U}) \varphi_x) \, dx \, dt \ge 0 \qquad (13) $$

for any test function $\varphi \in \mathbf{C}_c^\infty((0, \infty) \times \mathbb{R}; [0, \infty))$.

We defer the proof to Section 6.1. Some comments on Theorem 2.1 are in order.

**Remark 2.2.** In [26], the authors apply the wave-front tracking method to the ARZ model and obtain an existence result that is valid only away from the vacuum. However, in [28], the authors show that even if the vacuum is not present initially and does not appear immediately, the solution to the ARZ model can reach a vacuum state in finite time. This shows that a theory for the ARZ model away from the vacuum is quite restrictive. We recall that the extension to the case of solutions attaining the vacuum state is proved in [27], [28] and [34], by means of a deterministic particle approximation, BV bounds on the Riemann invariants and the compensated compactness method, respectively. This motivates the extension of our existence result in Theorem 2.1 to the vacuum case in a future paper.

**Remark 2.3.** Away from the vacuum, system (6) can be written in Lagrangian coordinates [42]

$$ \left\{ \begin{array}{l} \partial_t \tilde{\tau} - \partial_{\tilde{x}} (\tilde{c} V (\tilde{w} - p(1/\tilde{\tau}))) = 0, \\ \partial_t \tilde{w} = 0, \end{array} \right. \qquad (14) $$

where we have

$$ t = \tilde{t}, \qquad \partial_{\tilde{x}} x = \tilde{\tau}, \qquad \partial_{\tilde{t}} x = \tilde{c} V \left( \tilde{w} - p \left( \frac{1}{\tilde{\tau}} \right) \right), $$

and

$$ \tilde{\tau}(\tilde{t}, \tilde{x}) = \frac{1}{\rho(t(\tilde{t}), x(\tilde{t}, \tilde{x}))}, \qquad \tilde{w}(\tilde{t}, \tilde{x}) = \frac{q(t(\tilde{t}), x(\tilde{t}, \tilde{x}))}{\rho(t(\tilde{t}), x(\tilde{t}, \tilde{x}))}, \qquad \tilde{c}(\tilde{t}, \tilde{x}) = c(x(\tilde{t}, \tilde{x})). $$

The flux function of system (14) is discontinuous (with respect to both the Lagrangian coordinates $(\tilde{t}, \tilde{x})$), and this prevents from a direct application of [12, Theorem 7.1] to it. However, as in the proof of Theorem 2.1 in Section 6.1, we can introduce a 3 × 3 system of conservation laws, which in this case is

$$ \left\{ \begin{array}{l} \partial_t \tilde{\tau} - \partial_{\tilde{x}} \left( \frac{\tilde{k}}{\tilde{\tau}} V (\tilde{w} - p(\frac{1}{\tilde{\tau}})) \right) = 0, \\ \partial_t \tilde{w} = 0, \\ \partial_t \tilde{k} - \partial_{\tilde{x}} \left( \left( \frac{\tilde{k}}{\tilde{\tau}} \right)^2 V (\tilde{w} - p(\frac{1}{\tilde{\tau}})) \right) = 0, \end{array} \right. $$

where $\tilde{k} := \tilde{c}\tilde{\tau}$. The above system satisfies the hypotheses of [12, Theorem 7.1] and then it is possible to get an existence result analogous to that in Theorem 2.1.

**Remark 2.4.** In [39, Section 4], the author studies the ARZ model with discontinuous flux. Differently from (6), where both the components of the flux are discontinuous due to the presence of the discontinuous coefficient c, in [39, (4.1)-(4.2)] only the second component of the flux is discontinuous (see the discontinuous coefficient k multiplying the pressure term $\rho^\gamma$).
As a first step, the author adds another equation to treat the discontinuity, like we do in the proof of Theorem 2.1, obtains [39, (4.1)-(4.3)] and then rewrites such 3 × 3 system of conservation laws in decoupling coordinates $(\phi, \psi)$, see [36].

We can proceed in an analogous way and obtain the 3 × 3 system of conservation laws

$$ \left\{ \begin{array}{l} \partial_{\psi} \left( \frac{\tilde{\tau}}{\tilde{c}V(\tilde{w}-p(\frac{1}{\tilde{c})))} \right) - \partial_{\phi} \left( \frac{1}{\tilde{c}V(\tilde{w}-p(\frac{1}{\tilde{c})))} \right) = 0, \\ \partial_{\psi} \tilde{w} = 0, \\ \partial_{\phi} \tilde{c} = 0, \end{array} \right. \qquad (15) $$

by introducing in (6) the change of coordinates $(\phi, \psi)$ defined as

$$ \phi_x = -\rho, \qquad \phi_t = c \rho V (h(\rho, q)), \qquad \psi = x, $$

which is well defined when $c \rho V (h(\rho, q)) > 0$, and by letting

$$ \tilde{\tau}(\psi, \phi) = \frac{1}{\rho(t(\psi, \phi), x(\psi))}, \qquad \tilde{w}(\tilde{t}, \tilde{x}) = \frac{q(t(\psi, \phi), x(\psi))}{\rho(t(\psi, \phi), x(\psi))}, \qquad \tilde{c}(\psi) = c(x(\psi)). $$

The last two equations in (15) are decoupled. Moreover, given the initial datum, the values of $(\tilde{w}, \tilde{c})$ for any coordinate point $(\psi, \phi)$ are determined trivially, see [39, Figure 2] for the case of a Riemann datum. Then we can substitute the obtained expressions for $(\tilde{w}, \tilde{c})$ in the first equation and study the resulting scalar conservation law, which has a flux that is in general discontinuous with respect to both the coordinates $(\psi, \phi)$, and can be handled by applying the techniques developed in [17].

**Remark 2.5.** Solutions to system (6) may develop non-classical shocks at the discontinuity points $\xi_i$, $i \in [1, p]$, of the flux $cF$. These discontinuities may fail to satisfy the classical entropy condition (13) for general entropy pairs, not necessarily satisfying (12). At first glance, this seems to contradict the second assertion of Theorem 2.1. To clarify this apparent discrepancy, in Section 4 we construct general entropy pairs, clarifying the role of condition (12). Roughly speaking, we will demonstrate that entropy pairs satisfying (12) cannot select physically admissible solutions. On the other hand, as pointed out in the proof of Theorem 2.1 in Section 6.1, any entropy pair $(\mathcal{E}, c\mathcal{Q})$ satisfies (12) whenever $c$ is not piecewise constant.

A possible way to circumvent this drawback of the entropy condition is to explicitly introduce a (coupling) Riemann solver tailored to (6). As a matter of fact, the discontinuity points of the flux are typically addressed through a Riemann solver, which prescribes the solution's local behaviour at these points $\xi_i$, $i \in [1, p]$, rather than relying solely on the PDE framework. In fact, by prescribing a Riemann solver at $\xi_i$, $i \in [1, p]$, one explicitly encodes the underlying modeling assumptions. This approach is commonly employed when dealing with non-classical solutions to hyperbolic conservation laws, relying on the property of finite propagation speed. We refer to [3] for a comprehensive discussion of discontinuous-flux conservation laws defined via the Riemann solver method. In Section 3, we adapt the general definition of a Riemann solver, as proposed in [37, Section 3.2], to our case.

## 2.2. Applications

We compare now our model (1) with other models from the literature.

**ARZ model.** The ARZ model [5, 43] for vehicular traffic, written in the conservative form, is expressed by the following equations

$$ \left\{ \begin{array}{l} \partial_t \rho + \partial_x (\rho v) = 0, \\ \partial_t (\rho (v + p(\rho))) + \partial_x (\rho v (v + p(\rho))) = 0, \end{array} \right. \qquad (16) $$
where $v$ is the velocity and $p(\rho)$ is a pseudo-pressure function accounting for drivers' anticipation of downstream density changes. This model consists in two equations: the first one expresses the conservation of the density of cars $\rho$, the second one is a balance law for the momentum $\rho v$. System (1) can be rewritten in the form (16) by taking $c \equiv 1$, $h = v$ and $V$ as the identity function. For this reason, (1) can be interpreted as a generalization of (16) to the case of a discontinuous velocity. We can therefore apply Theorem 2.1 to the ARZ model (16) provided that $p$ satisfies (2) and (4). Moreover, Theorem 2.1 also applies to a generalization of the ARZ model to the case of a discontinuous flux, obtained by multiplying the flux by a coefficient function $c$ satisfying (5).

**CGST model.** The vehicular traffic model presented in [14] is represented by the system

$$ \left\{ \begin{array}{l} \partial_t \rho + \partial_x (c V(h) \rho) = 0, \\ \partial_t (\rho h) + \partial_x (c V(h) \rho h) = \kappa \rho^2 \partial_x (c V(h)). \end{array} \right. \qquad (17) $$

Here $h$ is the mean headway, i.e., the mean distance between two consecutive vehicles. This macroscopic model is derived through kinetic techniques and can be interpreted as a generalized version of the ARZ model (16). We can write (17) in the form (1) by taking $p(\rho) = \kappa \rho$, $\kappa > 0$, as pressure function. Indeed, by multiplying the first equation in (17) by $2\kappa \rho$, we formally get

$$ \partial_t (\kappa \rho^2) + 2 \kappa \rho \partial_x (c V(h) \rho) = 0. $$

Let's use the elementary identity $2f(fg)' = (f^2g)' + f^2g'$ for $f = f(x)$, $g = g(x)$, to find that

$$ 2 \rho \partial_x (c V(h) \rho) = \partial_x (c V(h) \rho^2) + \rho^2 \partial_x (c V(h)). $$

Therefore, for $p(\rho) = \kappa \rho$, we reach

$$ \partial_t (\rho p(\rho)) + \partial_x (c V(h) \rho p(\rho)) + \rho p(\rho) \partial_x (c V(h)) = 0. $$

By combining the preceding equation with the latter in (17), specifically

$$ \partial_t (\rho h) + \partial_x (c V(h) \rho h) = \rho p(\rho) \partial_x (c V(h)), $$

we derive the second equation in (1). Hence, Theorem 2.1 applies to the CGST model if the functions $V$ and $c$ are such that (3)-(5) hold.

**Comparison with GARZ and sedimentation.** The $2 \times 2$ system

$$ \left\{ \begin{array}{l} \partial_t \phi + \partial_x (c V(\phi, k) \phi) = 0, \\ \partial_t (k \phi) + \partial_x (c V(\phi, k) k \phi) = 0, \end{array} \right. \qquad (18) $$

can be interpreted as a “discontinuous flux” version of the model proposed in [7]. System (18) reduces to our system (6) by choosing an appropriate velocity function $V$, and then one can apply Theorem 2.1. Recall that in the case $c \equiv 1$, depending on the choice of the velocity function $V$, (18) corresponds to the so called generalized ARZ model (GARZ) proposed in [25], or to the sedimentation model, see [8].

## 2.3. Main properties of system (11)

In this section we analyze the main properties of system (6) in $U_e$ defined in (10). We defer the reader to [24, Chapter 7] for a general introduction to systems of conservation laws in one-space dimension.

The eigenvalues and the corresponding eigenvectors of system (6) are

$$ \lambda_1(\mathbf{U}, x) = c(x) (V(h(\mathbf{U}))) - \rho p'(\rho) V'(h(\mathbf{U}))) , \qquad r_1(\mathbf{U}) = \mathbf{U}, \qquad (19) $$

$$ \lambda_2(\mathbf{U}, x) = c(x) V(h(\mathbf{U})), \qquad r_2(\mathbf{U}) = (\rho, q + \rho^2 p'(\rho))^T, \qquad (20) $$
where $h(\mathbf{U}) := \frac{q}{\rho} - p(\rho)$. System (6) is strictly hyperbolic in $\mathcal{U}_e$ because $\lambda_1(\mathbf{U}, x) < \lambda_2(\mathbf{U}, x)$ for any $\mathbf{U} \in \mathcal{U}_e$ and $x \in \mathbb{R}$ by (2)$_3$, (3)$_3$ and (5)$_3$. Moreover, the first and second characteristic fields are respectively genuinely nonlinear and linearly degenerate because

$$D_{\mathbf{U}}\lambda_1(\mathbf{U}, x) \cdot r_1(\mathbf{U}) = \alpha(\mathbf{U}, x), \qquad D_{\mathbf{U}}\lambda_2(\mathbf{U}, x) \cdot r_2(\mathbf{U}) = 0,$$

where

$$\alpha(\mathbf{U}, x) := c(x) \rho (\rho p'(\rho)^2 V''(h(\mathbf{U})) - (2p'(\rho) + \rho p''(\rho)) V'(h(\mathbf{U})))$$

does not vanish (indeed it is $< 0$) for any $\mathbf{U} \in \mathcal{U}_e$ and $x \in \mathbb{R}$ by (4)$_2$ and (5)$_3$.

The Lagrangian marker $w(\mathbf{U}) := \frac{q}{\rho}$ and the function $h(\mathbf{U}) := \frac{q}{\rho} - p(\rho)$ form a coordinate system of Riemann invariants for (6) by [24, Theorem 7.3.3] because

$$D_{\mathbf{U}}w(\mathbf{U}) \cdot r_2(\mathbf{U}) = \rho p'(\rho), \qquad D_{\mathbf{U}}w(\mathbf{U}) \cdot r_1(\mathbf{U}) = 0,$$

$$D_{\mathbf{U}}h(\mathbf{U}) \cdot r_1(\mathbf{U}) = -\rho p'(\rho), \qquad D_{\mathbf{U}}h(\mathbf{U}) \cdot r_2(\mathbf{U}) = 0.$$

In particular $w$ and $h$ are respectively the first and second Riemann invariants, see [24, Definition 7.3.1]. Notice that the Riemann invariants $w$ and $h$ do not depend on $x$.

As already observed, $p(\rho) > 0$ for any $\rho > 0$, hence $w(\mathbf{U}) > h(\mathbf{U}) > 0$ for any $\mathbf{U} \in \mathcal{U}_e$. Moreover, by definition and the strict monotonicity of the pressure, we have $\rho = p^{-1}(w-h)$ and $q = p^{-1}(w-h)w$. Observe that $\mathcal{U}$ is expressed in the Riemann invariant coordinates as

$$\mathcal{W} := \{(h, w)^T \in \mathbb{R}^2 : w \ge h \ge 0\}. \qquad (21)$$

The vacuum $\rho = 0$ corresponds to $w = h$ and therefore $\mathcal{U}_m$ is expressed in the Riemann invariant coordinates as

$$\mathcal{W}_m := \{(h, w)^T \in \mathbb{R}^2 : w > h \ge 0\}.$$

At last, $\mathcal{U}_e$ becomes

$$\mathcal{W}_e := \{(h, w)^T \in \mathbb{R}^2 : w > h > 0\}.$$

Fix $\mathbf{U}_0 \in \mathcal{U}_e$. The 1-shock and 1-rarefaction Lax curves through $\mathbf{U}_0$ coincide and are described by

$$\mathcal{L}_1^{\mathbf{U}_0} := \{\mathbf{U} \in \mathcal{U}_m : w(\mathbf{U}) = w(\mathbf{U}_0)\}.$$

Indeed, if $\Pi_x$ is the inverse function of $(0, \infty) \ni r \mapsto \lambda_1(r, w(\mathbf{U}_0)r, x)$, which is strictly decreasing by (4)$_2$, then we can parametrize $\mathcal{L}_1^{\mathbf{U}_0}$ with $\phi(\nu) := \Pi_x(\nu)(1, w(\mathbf{U}_0))$, $\nu > \Pi_x^{-1}(\rho_0) = \lambda_1(\mathbf{U}_0, x)$, which is the solution of

$$\begin{cases} \alpha(\phi(\nu)) \phi'(\nu) = r_1(\phi(\nu)), \\ \lambda_1(\phi(\nu), x) = \nu, \end{cases}$$

and therefore $\mathcal{L}_1^{\mathbf{U}_0}$ is the 1-rarefaction Lax curve through $\mathbf{U}_0$. Moreover, $\mathcal{L}_1^{\mathbf{U}_0}$ is also the 1-shock Lax curve through $\mathbf{U}_0$ because, for any $\nu < \Pi_x^{-1}(\rho_0)$, there exists $\sigma = \sigma(\nu, x)$ such that

$$(\phi(\nu) - \mathbf{U}_0) \sigma = c(x) (F(\phi(\nu)) - F(\mathbf{U}_0)).$$

Indeed, it is sufficient to take

$$\sigma(\nu, x) := c(x) \frac{V(w(\mathbf{U}_0) - p(\Pi(\nu))) \Pi(\nu) - V(w(\mathbf{U}_0) - p(\rho_0)) \rho_0}{\Pi(\nu) - \rho_0}.$$

The second characteristic field is linearly degenerate, hence the associated Lax curve is expressed in terms of the first Riemann invariant as follows

$$\mathcal{L}_2^{\mathbf{U}_0} := \{\mathbf{U} \in \mathcal{U}_m : h(\mathbf{U}) = h(\mathbf{U}_0)\}.$$
System (6) is of Temple type [40, 41] because the 1-shock and 1-rarefaction Lax curves coincide and the second characteristic field is linearly degenerate.

**Remark 2.6.** Existence of solution and well-posedness results for Temple systems have been established in several papers, see [2, 6, 9, 13, 18, 38] and the references therein. These results were extended to initial-boundary value problems in [20, 21].

Recently, in [10] an existence result is obtained for Temple systems of nonlinear hyperbolic balance laws on networks via a combination of the wave-front tracking method, the pseudo-polygonal technique and the operator splitting. Such a result can be in principle applied to (6). Indeed, a discontinuity point $\xi_i$ of c can be interpreted as a junction with one incoming road and one outgoing road. However, the hypotheses of Theorem 2.1 differ from those in [10, Theorem 2.7]. Indeed, in [10, Section 2.1] it is assumed that the function c is piecewise constant, V is $C^3$ and $\max_{U_m} \lambda_1 \le \min_{U_m} \lambda_2$, where $\lambda_1$ and $\lambda_2$ are defined in $(19)_1$ and $(20)_1$.

### 3. Riemann solver

In this section we introduce our Riemann solver for (6). To do so, we have to assume that c is a piecewise constant function. It is not limitative to assume that c has a single discontinuity at $x = 0$.

We rewrite (6) in terms of the Riemann invariants and consider the Riemann problem

$$ 
\left\{ 
\begin{array}{l} 
\partial_t \mathbf{U}(\mathbf{W}) + \partial_x (c F(\mathbf{U}(\mathbf{W}))) = \mathbf{0}, \\ 
\mathbf{W}(0, x) = \begin{cases} \mathbf{W}_L & \text{if } x < 0, \\ \mathbf{W}_R & \text{if } x \ge 0, \end{cases} 
\end{array} 
\right. 
\qquad c(x) := \begin{cases} c^- & \text{if } x < 0, \\ c^+ & \text{if } x \ge 0, \end{cases} 
\qquad (22) 
$$

where $\mathbf{W}_L, \mathbf{W}_R \in W$ and $c^-, c^+ > 0$ are given constants, $F$ is defined in $(7)_2$ and $\mathbf{W} := (h, w)^T \in W$ is the dependent variable. With a slight abuse of notation we let

$$ \mathbf{U}(\mathbf{W}) := \begin{pmatrix} \rho(\mathbf{W}) \\ q(\mathbf{W}) \end{pmatrix} := \begin{pmatrix} p^{-1}(w-h) \\ p^{-1}(w-h)w \end{pmatrix}, \qquad \mathbf{W}(\mathbf{U}) := \begin{pmatrix} h(\mathbf{U}) \\ w(\mathbf{U}) \end{pmatrix} := \begin{pmatrix} \frac{q}{\rho} - p(\rho) \\ q/\rho \end{pmatrix}. \qquad (23) $$

Notice that

$$ F(\mathbf{U}(\mathbf{W})) = V(h) p^{-1}(w-h) (1, w)^T. \qquad (24) $$

The motivation to study the system in (22) rather than (1) or (6) is twofold. First, using the Riemann invariants as dependent variables simplifies the forthcoming exposition and analysis. Second, the fluxes in both (1) and (6) are not well-defined at the vacuum, see $(7)_3$, whereas the flux in (22) is, see (24). We recall that, even if the vacuum is not initially present and does not appear immediately, a solution of (6) can, in general, reach a vacuum state in finite time, see Remark 2.2. For this reason, it is preferable to have a framework that can handle the vacuum.

For simplicity in the exposition, we first consider the case $c^- = c^+$ and then the general case.

#### 3.1. The case $c^- = c^+$

Following the argument in [37], we give now the definition of the Riemann solver $RS_{\bar{c}}$ for the Rimeann problem (22) in the case $c^\pm = \bar{c}$ for a given constant $\bar{c} > 0$. We first need to introduce some notation. Let $f: W \to [0, \infty)$ be defined by

$$ f(h, w) := V(h) \rho(h, w). $$
For any $w > 0$, let

$$ 
\begin{aligned}
\lambda_w: [0, w] & \to [V(w), -p^{-1}(w) p'(p^{-1}(w)) V'(0)] \\
h & \mapsto V(h) - \rho(h, w) p'(\rho(h, w)) V'(h)
\end{aligned}
$$

and consider its inverse function

$$ \mathcal{H}_w: [V(w), -p^{-1}(w) p'(p^{-1}(w)) V'(0)] \to [0, w]. $$

By definition we have $\bar{c}\lambda_w(h) = \lambda_1(U(W), x)$, moreover $\lambda_w(w-p(\rho)) = \partial_\rho f(w-p(\rho), w)$ because $f(h, w) = V(h)p^{-1}(w-h)$. Observe that $\rho \mapsto f(w-p(\rho), w)$ is strictly concave by (4)$_2$. Introduce $s: \{(\mathbf{W}_L, \mathbf{W}_R) \in \mathcal{W}^2 : \rho(\mathbf{W}_L) \neq \rho(\mathbf{W}_R)\} \to \mathbb{R}$ defined by

$$ s(\mathbf{W}_L, \mathbf{W}_R) := \frac{f(\mathbf{W}_R) - f(\mathbf{W}_L)}{\rho(\mathbf{W}_R) - \rho(\mathbf{W}_L)}. \quad (25) $$

Observe that if a weak solution has a discontinuity between $\mathbf{W}_L$ and $\mathbf{W}_R$, then it has speed of propagation $\bar{c} s(\mathbf{W}_L, \mathbf{W}_R)$.

**Definition 3.1.** Fix $\bar{c} > 0$. The Riemann solver $\mathcal{R}\mathcal{S}_{\bar{c}}: \mathcal{W}^2 \to \mathbf{BV}(\mathbb{R}; \mathcal{W})$ for the Riemann problem (22) with $c^\pm = \bar{c}$ is defined as follows:

(L.1) If $\mathbf{W}_L, \mathbf{W}_R \in \mathcal{W}$ with $w_L \neq w_R$ and $h_L = h_R$, then

$$ \mathcal{R}\mathcal{S}_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R](\nu) := \begin{cases} \mathbf{W}_L & \text{if } \nu/\bar{c} < V(h_{L,R}), \\ \mathbf{W}_R & \text{if } \nu/\bar{c} \ge V(h_{L,R}). \end{cases} $$

(L.2) If $\mathbf{W}_L, \mathbf{W}_R \in \mathcal{W}$ with $w_L = w_R$ and $h_L > h_R$, then

$$ \mathcal{R}\mathcal{S}_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R](\nu) := \begin{cases} \mathbf{W}_L & \text{if } \nu/\bar{c} < s(\mathbf{W}_L, \mathbf{W}_R), \\ \mathbf{W}_R & \text{if } \nu/\bar{c} \ge s(\mathbf{W}_L, \mathbf{W}_R). \end{cases} $$

(L.3) If $\mathbf{W}_L, \mathbf{W}_R \in \mathcal{W}$ with $w_L = w_R$ and $h_L < h_R$, then

$$ \mathcal{R}\mathcal{S}_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R](\nu) := \begin{cases} \mathbf{W}_L & \text{if } \nu/\bar{c} < \lambda_{w_{L,R}}(h_L), \\ (\mathcal{H}_{w_{L,R}}(\nu), w_{L,R}) & \text{if } \lambda_{w_{L,R}}(h_L) \le \nu/\bar{c} < \lambda_{w_{L,R}}(h_R), \\ \mathbf{W}_R & \text{if } \nu/\bar{c} \ge \lambda_{w_{L,R}}(h_R). \end{cases} $$

(L.4) If $\mathbf{W}_L, \mathbf{W}_R \in \mathcal{W}$ with $w_L \neq w_R$ and $h_L > h_R$, then

$$ \mathcal{R}\mathcal{S}_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R](\nu) := \begin{cases} \mathbf{W}_L & \text{if } \nu/\bar{c} < s(\mathbf{W}_L, \mathbf{W}_M), \\ \mathbf{W}_M & \text{if } s(\mathbf{W}_L, \mathbf{W}_M) \le \nu/\bar{c} < V(h_R), \\ \mathbf{W}_R & \text{if } \nu/\bar{c} \ge V(h_R), \end{cases} $$

where $\mathbf{W}_M := (h_M, w_M)$ with $h_M := h_R$ and $w_M := w_L$.

(L.5) If $\mathbf{W}_L, \mathbf{W}_R \in \mathcal{W}$ with $w_L \neq w_R$ and $h_L < h_R$, then

$$ \mathcal{R}\mathcal{S}_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R](\nu) := \begin{cases} \mathbf{W}_L & \text{if } \nu/\bar{c} < \lambda_{w_L}(h_L), \\ (\mathcal{H}_{w_L}(\nu), w_L) & \text{if } \lambda_{w_L}(h_L) \le \nu/\bar{c} < \lambda_{w_L}(h_M), \\ \mathbf{W}_M & \text{if } \lambda_{w_L}(h_M) \le \nu/\bar{c} < V(h_R), \\ \mathbf{W}_R & \text{if } \nu/\bar{c} \ge V(h_R), \end{cases} $$

where $\mathbf{W}_M := (h_M, w_M)$ with $h_M := \min\{h_R, w_L\}$ and $w_M := w_L$.

(L.6) If $\mathbf{W}_L = \mathbf{W}_R$, then $\mathcal{R}\mathcal{S}_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R] \equiv \mathbf{W}_{L,R}$.
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Construction of $RS_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R]$ in the cases (L.4), left, and (L.5), center and right. Above $\mathbf{W}_R^1$ and $\mathbf{W}_R^2$ represent two possible choices for the right state $\mathbf{W}_R$, moreover $f_L(\rho) = f(w_L - p(\rho), w_L)$.</figcaption></figure>

Some comments on the above definition are in order. In case (L.1) we have that $RS_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R]$ is the 2-contact discontinuity $C_2(\mathbf{W}_L, \mathbf{W}_R)$. In case (L.2) we have that $0 \le \rho(\mathbf{W}_L) < \rho(\mathbf{W}_R) \le p^{-1}(w_{L,R})$ and $RS_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R]$ is the 1-shock $S_1(\mathbf{W}_L, \mathbf{W}_R)$. In case (L.3) we have that $0 \le \rho(\mathbf{W}_R) < \rho(\mathbf{W}_L) \le p^{-1}(w_{L,R})$ and $RS_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R]$ is the 1-rarefaction $R_1(\mathbf{W}_L, \mathbf{W}_R)$. In case (L.4) we have that $0 \le \rho(\mathbf{W}_L) < \rho(\mathbf{W}_M) \le p^{-1}(w_L)$ and $RS_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R]$ is the juxtaposition of $S_1(\mathbf{W}_L, \mathbf{W}_M)$ and $C_2(\mathbf{W}_M, \mathbf{W}_R)$, see Figure 1, left. In case (L.5) we have that $0 \le \rho(\mathbf{W}_M) < \rho(\mathbf{W}_L) \le p^{-1}(w_L)$, with $\rho(\mathbf{W}_M) = 0$ if and only if $h_R \ge w_L$, and $RS_{\bar{c}}[\mathbf{W}_L, \mathbf{W}_R]$ is the juxtaposition of $R_1(\mathbf{W}_L, \mathbf{W}_M)$ and $C_2(\mathbf{W}_M, \mathbf{W}_R)$, see Figure 1, center and right.

We conclude this section with two remarks.

**Remark 3.2.** According to (22)$_{c^\pm=1}$, if a vehicle is characterized by Lagrangian marker w, then it has maximal speed $V(w)$ and length $1/p^{-1}(w)$. Indeed, if the vehicles have the same Lagrangian marker w and are bumper-to-bumper, then their velocity is zero, $V(h) = 0$; by (3) this corresponds to $h = 0$ and then by (23)$_2$ to the density $\rho = p^{-1}(w)$, hence in any interval with length $L = 1/p^{-1}(w)$ there is $\rho L = 1$ vehicle.

**Remark 3.3.** If $c^- = c^+$, then the fundamental diagrams $\{(\rho, f) : f = c^\pm f(w - p(\rho), w)\}$, $w > 0$, do not intersect away from the vacuum. This property is lost in the case $c^- \ne c^+$.

### 3.2. The case $c^- \ne c^+$

Denote by $RS_\pm, RS: W^2 \to BV(\mathbb{R}; W)$ the Riemann solvers given in Definition 3.1, corresponding to $c = c^\pm$ and $c \equiv 1$, respectively. Before giving a general definition for the Riemann solver for (6), we need to introduce some notation. Let $\mathbf{W}_* := (h_*, w_*) : W \times (0, \infty) \to W$ be defined by

$$h_*(\mathbf{W}_R, w_L) := \min\{w_L, h_R\}, \qquad w_*(\mathbf{W}_R, w_L) := w_L. \qquad (26)$$

Observe that the maximum of $h \mapsto f(h, w)$, which is

$$F(w) := \max_{h \in [0, w]} f(h, w),$$

is attained at a unique value of $h$, that we denote by $h(w)$ so that

$$f(h(w), w) = F(w).$$
Define $\mathbb{Q}_-: \mathcal{W} \to [0, \infty)$ and $\mathbb{Q}_+: \mathcal{W} \times (0, \infty) \to [0, \infty)$ as follows:

$$ \begin{aligned} \mathbb{Q}_-(\mathbf{W}_L) &:= \max_{h \in [0, w_L]} f(\mathcal{R}\mathcal{S}[\mathbf{W}_L, (h, w_L)](0^{-})) \\ &= \begin{cases} F(w_L) & \text{if } h_L \le h(w_L), \\ f(\mathbf{W}_L) & \text{if } h_L > h(w_L), \end{cases} \end{aligned} \tag{27} $$

$$ \begin{aligned} \mathbb{Q}_+(\mathbf{W}_R, w_L) &:= \max_{h \in [0, w_L]} f(\mathcal{R}\mathcal{S}[(h, w_L), \mathbf{W}_*(\mathbf{W}_R, w_L)](0^+)) \\ &= \begin{cases} f(\mathbf{W}_*(\mathbf{W}_R, w_L)) & \text{if } h_*(\mathbf{W}_R, w_L) < h(w_L), \\ F(w_L) & \text{if } h_*(\mathbf{W}_R, w_L) \ge h(w_L). \end{cases} \end{aligned} \tag{28} $$

Let then $\hat{\mathbf{W}} := (\hat{h}, \hat{w})$, $\check{\mathbf{W}} := (\check{h}, \check{w})$: $\mathcal{W}^2 \to \mathcal{W}$ be defined by

$$ \begin{cases} \hat{h}(\mathbf{W}_L, \mathbf{W}_R) := \min \{h \in [0, h(w_L)] : c^- f(h, w_L) = \mathbb{Q}(\mathbf{W}_L, \mathbf{W}_R)\}, \\ \hat{w}(\mathbf{W}_L, \mathbf{W}_R) := w_L, \end{cases} \tag{29} $$

$$ \begin{cases} \check{h}(\mathbf{W}_L, \mathbf{W}_R) := \max \{h \in [h(w_L), w_L] : c^+ f(h, w_L) = \mathbb{Q}(\mathbf{W}_L, \mathbf{W}_R)\}, \\ \check{w}(\mathbf{W}_L, \mathbf{W}_R) := w_L, \end{cases} \tag{30} $$

where

$$ \mathbb{Q}(\mathbf{W}_L, \mathbf{W}_R) := \min \{c^- \mathbb{Q}_-(\mathbf{W}_L), c^+ \mathbb{Q}_+(\mathbf{W}_R, w_L)\}. \tag{31} $$

Observe that

$$ f(\check{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)) = f(\check{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)) = \mathbb{Q}(\mathbf{W}_L, \mathbf{W}_R). $$

We are now in the position to give the following definition.

**Definition 3.4.** The Riemann solver $\mathcal{R}\mathcal{S}_{c^-, c^+}: \mathcal{W}^2 \to \mathbf{B}\mathbf{V}(\mathbb{R}; \mathcal{W})$ for (22) is defined as follows:

$$ \mathcal{R}\mathcal{S}_{c^-, c^+}[\mathbf{W}_L, \mathbf{W}_R](\nu) := \begin{cases} \mathcal{R}\mathcal{S}_-[\mathbf{W}_L, \hat{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)](\nu) & \text{if } \nu < 0, \\ \mathcal{R}\mathcal{S}_+[\check{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R), \mathbf{W}_R](\nu) & \text{if } \nu \ge 0. \end{cases} \tag{32} $$

In the next proposition we analyze the possible traces of $\mathcal{R}\mathcal{S}_{c^-, c^+}[\mathbf{W}_L, \mathbf{W}_R]$ at $\nu = 0$.

**Proposition 3.5.** Fix $\mathbf{W}_L, \mathbf{W}_R \in \mathcal{W}$ and let $\mathbf{W}_\pm := \mathcal{R}\mathcal{S}_{c^-, c^+}[\mathbf{W}_L, \mathbf{W}_R](0^\pm)$.

* We have $\mathbf{W}_- \in \{\mathbf{W}_L, \hat{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)\}$ and

$$ \mathbf{W}_- = \mathbf{W}_L \neq \hat{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R) \iff h_L > h(w_L) \quad \text{and} \quad c^- f(\mathbf{W}_L) = \mathbb{Q}(\mathbf{W}_L, \mathbf{W}_R). $$

* We have $\mathbf{W}_+ \in \{\mathbf{W}_*(\mathbf{W}_R, w_L), \check{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)\}$ and

$$ \mathbf{W}_+ = \mathbf{W}_* \neq \check{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R) \iff h_* < h(w_L) \quad \text{and} \quad c^+ f(\mathbf{W}_*) = \mathbb{Q}(\mathbf{W}_L, \mathbf{W}_R), $$

where, to simplify the notation, we let $\mathbf{W}_* := (h_*, w_L) := \mathbf{W}_*(\mathbf{W}_R, w_L)$.

*Proof.* By (29) we know that $\hat{w}(\mathbf{W}_L, \mathbf{W}_R) = w_L$. Hence in $\nu < 0$ we have that $\mathcal{R}\mathcal{S}_{c^-, c^+}[\mathbf{W}_L, \mathbf{W}_R]$ is either constant or has a unique 1-wave. In the former case, we have $\mathcal{R}\mathcal{S}_{c^-, c^+}[\mathbf{W}_L, \mathbf{W}_R] \equiv \mathbf{W}_L$ in $\nu < 0$. In the latter case, we have $\mathbf{W}_- = \hat{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)$. Analogously, by (26) and (30) we know that $w_*(\mathbf{W}_R, w_L) = w_L = \check{w}(\mathbf{W}_L, \mathbf{W}_R)$. Hence in $\nu > 0$ we have that $\mathcal{R}\mathcal{S}_{c^-, c^+}[\mathbf{W}_L, \mathbf{W}_R]$ has at most one 1-wave. In this case, we have $\mathbf{W}_+ = \check{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)$, otherwise $\mathbf{W}_+ = \mathbf{W}_*(\mathbf{W}_R, w_L)$.

To complete the proof it is sufficient to observe that by Definition 3.4 we have

$$ c^- f(\mathbf{W}_-) = c^+ f(\mathbf{W}_+) = \mathbb{Q}(\mathbf{W}_L, \mathbf{W}_R) = c^- f(\hat{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)) = c^+ f(\check{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)), $$
and $w_{-} = w_{L} = \hat{w}(\mathbf{W}_{L}, \mathbf{W}_{R}) = \tilde{w}(\mathbf{W}_{L}, \mathbf{W}_{R}) = w_{*}(\mathbf{W}_{R}, w_{L})$. $\square$

Some remarks on Definition 3.4 are in order.

**Remark 3.6.** We stress that $u(t,x) := \mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}](x/t)$ is a weak solution corresponding to the initial datum $(\mathbf{W}_{L}, \mathbf{W}_{R}) \in W^{2}$. This simply follows from the fact that $\mathcal{R}S_{\pm}$ in (32) are Riemann solvers and that by (29), (30) we have

$$c^{-} f(\mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}](0^{-})) = Q(\mathbf{W}_{L}, \mathbf{W}_{R}) = c^{+} f(\mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}](0^{+})).$$

**Remark 3.7.** Definition 3.4 resembles [19, Definition 2.1] for a scalar conservation law with a point constraint on the flow at $x = 0$, which in our case becomes $c^{\pm}f(x = 0^{\mp}) \le Q(\mathbf{W}_{L}, \mathbf{W}_{R})$. The key difference is that we do not check whether the classical solution satisfies the constraint. Instead, we emphasize selecting, among the self-similar weak solutions that meet the constraint, the solution that optimizes the flow.

**Remark 3.8.** For any $(\mathbf{W}_{L}, \mathbf{W}_{R}) \in W^{2}$ such that $w_{L} \ne w_{R}$ we have that the right most wave of $\mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}]$ is a contact discontinuity with speed of propagation $c^{+}V(h_{R})$. We stress that 2- contact discontinuities always separate the two families of vehicles, characterized by the Lagrangian markers $w_{L}$ and $w_{R}$.

**Remark 3.9.** It is easy to see that $\nu \mapsto \mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}](\nu)$ may have total variation greater than that of the initial datum, that is $|h_{L} - h_{R}| + |w_{L} - w_{R}|$; moreover the maximum principle holds for the w-coordinate, but may fail for the h-coordinate. Indeed, if $\mathbf{W}_{L} = \mathbf{W}_{R}$ and $c^{-} \ne c^{+}$, then $\mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}]$ is not constant. We defer the reader to Section 5 for a detailed construction of the solutions corresponding to the constant initial datum (38) and two choices of the function c given in (39), see also Figures 2 and 3.

We conclude this section studying the main properties of the Riemann solver $\mathcal{R}S_{c^{-},c^{+}}$. Specifically, we consider the coherence, the $L_{\text{loc}}^{1}$-continuity with respect to the initial datum and give a sufficient condition on $(\mathbf{W}_{L}, \mathbf{W}_{R})$ to ensure that $\mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}]$ takes values in $W_{e}$. As in [22, Definition 2.2], the coherence of the Riemann solver means that for any $\mathbf{W}_{L}, \mathbf{W}_{R} \in W$ we have

$$\mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}](0^{\pm}) = \mathbf{W}_{\pm} \quad \implies \quad \mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{-}, \mathbf{W}_{+}](0^{\pm}) = \mathbf{W}_{\pm}.$$

Observe that coherence can be thought as a stability property: the ordered pair of traces is, in a sense, a fixed point for $\mathcal{R}S_{c^{-},c^{+}}$. Incoherence leads to the instability of the numerical solution, see for instance [23, Figure 5-center].

**Proposition 3.10.** The Riemann solver $\mathcal{R}S_{c^{-},c^{+}}: W^{2} \to BV(\mathbb{R}; W)$ satisfies the following properties:

1) It is coherent.

2) It is $L_{\text{loc}}^{1}$-continuous with respect to the initial datum.

3) If $\mathbf{W}_{L}, \mathbf{W}_{R} \in W_{e}$ and $w_{L} > h_{R}$, then $\mathcal{R}S_{c^{-},c^{+}}[\mathbf{W}_{L}, \mathbf{W}_{R}]$ takes values in $W_{e}$.

It is worth underlying that no condition on $c^{-}$ and $c^{+}$ are required in the above proposition. The proof is deferred to Section 6.2.

## 4. On the entropy condition

In Definition 4.1, we introduce the concept of weak solutions. In general, weak solutions may include non-physical shocks. This motivates the introduction of entropy solutions in Definition 4.2 to select
physically admissible shocks. To achieve this, we first construct a one-parameter family of entropy pairs, following the approach in [4], and then provide the definition of entropy solutions.

We begin by rewriting system (6) in Riemann invariant coordinates, as shown in (33). Furthermore, for simplicity, we assume that $c$ is piecewise constant; additional considerations regarding the general case are provided in Remark 4.3. Consider the Cauchy problem

$$ \begin{cases} \partial_t \mathbf{U}(\mathbf{W}) + \partial_x (c F(\mathbf{U}(\mathbf{W}))) = \mathbf{0}, \\ \mathbf{W}(0, x) = \bar{\mathbf{W}}(x), \end{cases} \qquad (33) $$

where $F$ is defined in (7)$_2$, $\mathbf{W} := (h, w)^T$ is the dependent variable expressed in terms of Riemann invariants, and the initial datum $\bar{\mathbf{W}} := (\bar{h}, \bar{w})^T$ is assumed to be in $\mathbf{L}^\infty(\mathbb{R}; W)$, where $W$ is defined in (21). With a slight abuse of notation, we define

$$ \mathbf{U}(\mathbf{W}) := \begin{pmatrix} \rho(\mathbf{W}) \\ q(\mathbf{W}) \end{pmatrix} := \begin{pmatrix} p^{-1}(w-h) \\ p^{-1}(w-h)w \end{pmatrix}, \qquad \mathbf{W}(\mathbf{U}) := \begin{pmatrix} h(\mathbf{U}) \\ w(\mathbf{U}) \end{pmatrix} := \begin{pmatrix} \frac{q}{\rho} - p(\rho) \\ q/\rho \end{pmatrix}. \qquad (34) $$

Notice that

$$ F(\mathbf{U}(\mathbf{W})) = V(h) p^{-1}(w-h) (1, w)^T. $$

The motivation to study (33) rather than (1) or (6) is twofold. First, the use of the Riemann invariants as independent variables eases the forthcoming analysis. Second, the entropy pairs are well defined at the vacuum in the Riemann invariant coordinates $\mathbf{W}$, whereas they are multi-valued in the conserved coordinates $\mathbf{U}$. Observe that, in general, even if the vacuum is not present initially and does not appear immediately, a solution of (6) can reach a vacuum state in finite time, see Remark 2.2. For this reason, a definition of solution possibly taking into account vacuum states is of interest. As first step to cope with the vacuum, we extend the flux $F$ to the whole $\mathcal{U}$ by taking $F(0, 0) := 0$. This choice is motivated by the observation that the flux $F(\mathbf{U})$ is the product of $\mathbf{U}$ and the uniformly bounded $V(h)$, where $0 \le h \le w \le \|\bar{w}\|_\infty$.

**Definition 4.1.** A weak solution of the Cauchy problem (33) with $\bar{\mathbf{W}} \in \mathbf{L}^\infty(\mathbb{R}; W)$ is a function $\mathbf{W}$ in $\mathbf{L}^\infty(\mathbb{R}_+ \times \mathbb{R}; W) \cap \mathbf{C}^0(\mathbb{R}_+; \mathbf{L}_{\text{loc}}^1(\mathbb{R}; W))$ that satisfies the initial condition a.e. in $\mathbb{R}$ and for any test function $\varphi \in \mathbf{C}_c^\infty((0, \infty) \times \mathbb{R}; [0, \infty))$

$$ \int_{\mathbb{R}_+} \int_{\mathbb{R}} p^{-1}(w-h) (\varphi_t + cV(h)\varphi_x) \, dx \, dt = 0, $$

$$ \int_{\mathbb{R}_+} \int_{\mathbb{R}} p^{-1}(w-h) w (\varphi_t + cV(h)\varphi_x) \, dx \, dt = 0. $$

Before giving the definition of entropy solution, we first construct a one parameter family of entropy pairs. Away from the vacuum, we can apply [24, (7.4.12)] and obtain that an entropy pair $(\mathcal{E}, c\mathcal{Q})$ expressed in terms of the Riemann invariants satisfies the system

$$ \partial_h (c(x) \mathcal{Q}) = \lambda_1 (\mathbf{U}(\mathbf{W}), x) \partial_h \mathcal{E}, $$

$$ \partial_w (c(x) \mathcal{Q}) = \lambda_2 (\mathbf{U}(\mathbf{W}), x) \partial_w \mathcal{E}, $$

where $\lambda_1$ and $\lambda_2$ are respectively defined in (19)$_1$ and (20)$_1$, together with the integrability condition [24, (7.4.13)]

$$ \rho(h, w) p'(\rho(h, w)) \partial_h \partial_w \mathcal{E} + \left( 1 + \frac{\rho(h, w) p''(\rho(h, w))}{p'(\rho(h, w))} \right) \partial_h \mathcal{E} + \partial_w \mathcal{E} = 0, $$
where $\rho(\mathbf{W})$ is defined in (34)$_1$. The general solutions of the above system are

$$ \mathcal{E}(h, w) := \left( \mathfrak{f}(w) - \int_{0}^{h} \frac{\zeta'(w - \nu) \mathfrak{g}(\nu)}{c V'(\nu)} d\nu \right) p^{-1}(w - h), \quad (35) $$

$$ \mathcal{Q}(h, w) := V(h) \mathcal{E}(h, w) - \int_{0}^{h} \mathfrak{g}(\nu) d\nu + \mathfrak{b}, $$

where $\mathfrak{f}$ and $\mathfrak{g}$ are sufficiently smooth functions, $\mathfrak{b}$ is an arbitrary constant and

$$ \zeta(\nu) := 1/p^{-1}(\nu). $$

In analogy to [4, (2.13)], we consider the one parameter family of entropy pairs obtained from (35) by taking $\mathfrak{f} \equiv 0$, $\mathfrak{g}(\nu) := V'(\nu) \operatorname{sgn}_+(\nu - k)$ and $\mathfrak{b} = 0$ and defined for any $k \ge 0$ as

$$ \mathcal{E}_k(h, w) := \begin{cases} 1 - \frac{p^{-1}(w-h)}{p^{-1}(w-k)} & \text{if } h \in (k, w], \\ 0 & \text{otherwise,} \end{cases} \quad \mathcal{Q}_k(h, w) := \begin{cases} V(k) - V(h) \frac{p^{-1}(w-h)}{p^{-1}(w-k)} & \text{if } h \in (k, w], \\ 0 & \text{otherwise.} \end{cases} \quad (36) $$

We are now in the position to give the definition of entropy solution analogous to [37, Definition 1.1].

**Definition 4.2.** A weak solution $\mathbf{W}$ of the Cauchy problem (33) with $\overline{\mathbf{W}} \in \mathbf{L}^{\infty}(\mathbb{R}; \mathcal{W})$ is an entropy solution if for any test function $\varphi \in \mathbf{C}_c^{\infty}((0, \infty) \times (\mathbb{R} \setminus \{\xi_1, \dots, \xi_p\}); [0, \infty))$ and for any $k > 0$

$$ \int_{\mathbb{R}_+} \int_{\mathbb{R}} (\mathcal{E}_k(\mathbf{W})) \partial_t \varphi + c \mathcal{Q}_k(\mathbf{W}) \partial_x \varphi \, dx \, dt \ge 0. \quad (37) $$

Furthermore, if $\mathbf{W}$ has a discontinuity along the curve $x = s(t)$, then its traces $\mathbf{W}^-(t) := \mathbf{W}(t, s(t)^-)$ and $\mathbf{W}^+(t) := \mathbf{W}(t, s(t)^+)$ satisfy the following conditions:

(A) For a.e. $t > 0$, if the discontinuity $x = s(t)$ occurs away from $x = \xi_i$, $i \in [1, p]$, and involves a vacuum state, then

$$ \mathcal{R} \mathcal{S}_{c(s(t))} [\mathbf{W}^{-}(t), \mathbf{W}^{+}(t)](\nu) = \begin{cases} \mathbf{W}^{-}(t) & \text{if } \nu < s(t), \\ \mathbf{W}^{+}(t) & \text{if } \nu \ge s(t). \end{cases} $$

(B) For a.e. $t > 0$, if the discontinuity $x = s(t)$ occurs at $x = \xi_i$, $i \in [1, p]$, then

$$ \mathcal{R} \mathcal{S}_{c(\xi_i^-), c(\xi_i^+)} [\mathbf{W}^{-}(t), \mathbf{W}^{+}(t)](\nu) = \begin{cases} \mathbf{W}^{-}(t) & \text{if } \nu < s(t), \\ \mathbf{W}^{+}(t) & \text{if } \nu \ge s(t). \end{cases} $$

Some comments on the above definition are in order. Differently from the test functions considered for the entropy condition (13) in Theorem 2.1, those in (37) are supported outside the set of discontinuities $\{\xi_1, \dots, \xi_p\}$ of the flux $cF$. Hence, condition (37) deals with discontinuity of $\mathbf{W}$ away from the discontinuity points $\{\xi_1, \dots, \xi_p\}$; the discontinuities of $\mathbf{W}$ at such points are addressed in (B). Furthermore, in (37) we exploit the entropy pairs (36), which do not select a physically reasonable solution if a vacuum state is involved, see [4, Remark 2.1]; this motivates condition (A).

**Remark 4.3.** Assume now that $c$ satisfies (5) and is not piecewise constant. In this case, we cannot introduce a (self-similar) Riemann solver. For this reason it would be desirable to rely again on condition (37) and the entropy pairs in (36) to select physically admissible shocks, at least away from the vacuum. However, this is not possible. Indeed, on one hand we can understand if a shock between $\mathbf{W}^- := (h_L, \bar{w})$ and $\mathbf{W}^+ := (h_R, \bar{w})$ is physically admissible by considering entropy pairs given in (36) with $k$ between $h_L$ and $h_R$. On the other hand, as underlined in the proof of Theorem 2.1, see
Section 6.1, if c is not piecewise constant, then the entropy pairs have to satisfy also condition (48) in Section 6, which is equivalent to (12). Unfortunately, we have that:

(A) The entropy pair given by (36) satisfies (12) if and only if $k \in \{0, w\}$.

This has a severe drawback because, as we mentioned before, the entropy pairs $(\mathcal{E}_w, \mathcal{Q}_0)$ and $(\mathcal{E}_w, \mathcal{Q}_w)$ alone are not sufficient to select physically reasonable shocks. We have in fact that:

(B) A discontinuity away from $\xi_i$, $i \in [1, p]$, between $\mathbf{W}^- := (h_L, \bar{w})$ and $\mathbf{W}^+ := (h_R, \bar{w})$ with $h_L < h_R$ satisfies the entropy condition (37) with $k \in \{0, w\}$.

Observe that the shock described in (B) is not physically reasonable.

We prove now (A). The case $k = w$ is indeed trivial because $(\mathcal{E}_w, \mathcal{Q}_w) = (0, 0)$. Assume now that $k < w$. If we rewrite the entropy pair given in (36) in the conservative coordinates, then for any $h \in (k, w]$ we get

$$ 
\begin{aligned}
D_{\mathbf{U}} \mathcal{E}_k (\mathbf{W}(\mathbf{U})) \cdot F(\mathbf{U}) &= \frac{V(h(\mathbf{U}))}{p^{-1}\left(\frac{q}{\rho}-k\right)^2} \left( -p^{-1}\left(\frac{q}{\rho}-k\right) - \frac{\rho p'(p^{-1}\left(\frac{q}{\rho}-k\right))}{1-p^{-1}\left(\frac{q}{\rho}-k\right)} \right) \cdot \begin{pmatrix} \rho \\ q \end{pmatrix} \\
&= -V(h(\mathbf{U})) \frac{\rho}{p^{-1}\left(\frac{q}{\rho}-k\right)}, \\
\mathcal{Q}_k (\mathbf{W}(\mathbf{U})) &= V(k) - V(h(\mathbf{U})) \frac{\rho}{p^{-1}\left(\frac{q}{\rho}-k\right)},
\end{aligned}
$$

and by imposing condition (12) we obtain that $V(k) = 0$, namely $k = 0$. This ends the proof of (A). We prove now (B), namely that if $k \in \{0, w\}$ then

$$ (\mathcal{E}_k(\mathbf{W}_R) - \mathcal{E}_k(\mathbf{W}_L)) c(x) s(\mathbf{W}_L, \mathbf{W}_R) - c(x) (\mathcal{Q}_k(\mathbf{W}_R) - \mathcal{Q}_k(\mathbf{W}_L)) = 0, $$

where $s$ is defined in (25) and $x \in \mathbb{R} \setminus \{\xi_1, \dots, \xi_p\}$ is the position where the discontinuity occurs. The case $k = w$ is indeed trivial because $(\mathcal{E}_w, \mathcal{Q}_w) = (0, 0)$. If $k = 0$, then the above condition writes

$$ \left( \frac{p^{-1}(w - h_L)}{p^{-1}(w)} - \frac{p^{-1}(w - h_R)}{p^{-1}(w)} \right) \frac{V(h_R) p^{-1}(w - h_R) - V(h_L) p^{-1}(w - h_L)}{p^{-1}(w - h_R) - p^{-1}(w - h_L)} \\ - \left( -V(h_R) \frac{p^{-1}(w - h_R)}{p^{-1}(w)} + V(h_L) \frac{p^{-1}(w - h_L)}{p^{-1}(w)} \right) = 0, $$

which is clearly satisfied for any $w > 0$. This completes the proof of (B).

## 5. Numerical Simulations

In this Section we perform some numerical simulations. For simplicity, below we always take

$$ V(h) := h, \qquad p(\rho) := \rho^2. $$

We discretize system (1) by using the classical Lax-Friedrichs numerical scheme [33, §4.6] adapted to deal with discontinuous fluxes. Below we choose a space mesh with length $\Delta x = 10^{-4}$ and a time step $\Delta t := 0.2 \Delta x$.

In Figure 2 we show two simulations corresponding to the constant initial datum

$$ \bar{\rho} = 0.4, \qquad \bar{q} = 0.4, \qquad \bar{h} = 0.84, \qquad \bar{w} = 1, \qquad (38) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Numerical simulations for the initial datum as in (38) and two possible choices for $c$ given in (39). Final time $T = 0.5$.</figcaption></figure>

and the two piecewise constant functions

$$ \text{left: } c(x) = \begin{cases} 1 & \text{if } x < 0, \\ 0.5 & \text{if } x \ge 0, \end{cases} \qquad \text{right: } c(x) = \begin{cases} 0.5 & \text{if } x < 0, \\ 1 & \text{if } x \ge 0. \end{cases} \qquad (39) $$

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Representation of the solutions to the Riemann problems corresponding to (38), (39). Above $W_L = W_R = (\bar{h}, \bar{w})$, $\tilde{W} = \dot{W}(W_L, W_R)$ and $\tilde{W} = \dot{W}(W_L, W_R)$.</figcaption></figure>

In these cases, the exact solutions are obtained by applying $RS_{1,0.5}$ and $RS_{0.5,1}$ given in Definition 3.4 to $W_L = (\bar{h}, \bar{w}) = W_R$. By construction, $RS_{1,0.5}[W_L, W_R]$ is the juxtaposition of the 1-shock $S_1(W_L, \dot{W}(W_L, W_R))$, the non-classical shock $N(\tilde{W}(W_L, W_R), \dot{W}(W_L, W_R))$ and the 1-rarefaction $R_1(\tilde{W}(W_L, W_R), W_R)$, see Figure 3, left. On the other hand, $RS_{0.5,1}[W_L, W_R]$ is the juxtaposition of the non-classical shock $N(W_L, \dot{W}(W_L, W_R))$ and the 1-shock $S_1(\dot{W}(W_L, W_R), W_R)$, see Figure 3, right. Observe that such behaviours are captured by the numerical simulations in Figure 2.
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Numerical simulations for the initial datum as in (38) and two possible choices for c given in (40). Final time $T = 0.5$.</figcaption></figure>

In Figure 4 we show simulations considering the same constant initial datum (38), but taking c periodic, with period $T = 0.5$, and defined for $x \in [0, 0.5)$ by

left: $c(x) = 1 - x$,

right: $c(x) = 0.5 + x$. (40)

The above simulations show the relevance and the effects of the discontinuous coefficient c on the solution. Indeed, in each of the both examples, the solution is not constant even if initial datum is so. In particular, the first two examples show that a single jump of c can produce a solution with two discontinuities. Moreover, the periodic behaviour of the solutions of the last two examples is due only to the periodic behaviour of c, as the initial datum is constant.

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Numerical simulations for the initial datum as in (41) and two possible choices for c given in (42). Final time $T = 0.5$.</figcaption></figure>
In Figure 5 we show numerical simulations with the following initial datum

$$ 
\begin{align} 
\bar{\rho}(x) &= \begin{cases} 0.97 & \text{if } x < 0, \\ 0.5 & \text{if } x \ge 0, \end{cases} \qquad \bar{q}(x) &= \begin{cases} 0.97 & \text{if } x < 0, \\ 0.75 & \text{if } x \ge 0, \end{cases} \tag{41} \\ 
\bar{h}(x) &= \begin{cases} 0.0591 & \text{if } x < 0, \\ 1.25 & \text{if } x \ge 0, \end{cases} \qquad \bar{w}(x) = \begin{cases} 1 & \text{if } x < 0, \\ 1.5 & \text{if } x \ge 0, \end{cases} \nonumber 
\end{align} 
$$

and c defined by

$$ 
\text{left: } c(x) = \begin{cases} 0.5 & \text{if } x < 0, \\ 1 & \text{if } x \ge 0, \end{cases} \qquad \text{right: } c(x) = \begin{cases} 0.5 & \text{if } x < -\epsilon, \\ \frac{x+\epsilon}{\epsilon} + 0.5 & \text{if } -\epsilon \le x < 0, \\ 1 & \text{if } x \ge 0, \end{cases} \tag{42} 
$$

with $\epsilon = 0.1$.

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Representation of the solution to the Riemann problem corresponding to (41) and (42)$_1$, left. Above $\tilde{W} = \tilde{W}(W_L, W_R)$ and $\tilde{W} = \tilde{W}(W_L, W_R)$.</figcaption></figure>

If c is given by (42), left, then the exact solution is obtained by computing $RS_{0.5,1}[W_L, W_R]$ with $W_L = (0.0591, 1)$ and $W_R = (1.25, 1.5)$. By construction, $RS_{0.5,1}[W_L, W_R]$ is the juxtaposition of the 1-rarefaction $R_1(W_L, \tilde{W}(W_L, W_R))$, the non-classical shock $N(\tilde{W}(W_L, W_R), \tilde{W}(W_L, W_R))$, the 1-rarefaction $R_1(\tilde{W}(W_L, W_R), (w_L, w_L))$ and the 2-contact discontinuity $C_2((w_L, w_L), W_R)$, see Figure 6.

Observe that the behaviour of the exact solution is captured by the numerical simulation in Figure 5, left. In particular, the numerical solution “reaches” the vacuum state both in the case with a discontinuous c, left, and in the case of a continuous c, right. However, it is worth pointing out that the w-component of the exact solution has a single discontinuity, which moves together with the 2-contact discontinuity $C_2((w_L, w_L), W_R)$. This is not respected by the behaviour of the numerically computed profile of w. A possible explanation is that the Lax-Friedrichs numerical scheme is diffusive. Moreover, handling the vacuum state numerically is generally challenging, as the expressions for both w and h involve division by $\rho$. This suggests the need for an ad hoc scheme. However, a detailed numerical study of system (1) is beyond the scope of this paper, where simulations are used solely to corroborate the analytical results.
# 6. Technical details

## 6.1. Proof of Theorem 2.1

Consider the Cauchy problem for a 3 × 3 system of conservation laws

$$ \left\{ \begin{array}{l} \partial_t \mathbf{X} + \partial_x G(\mathbf{X}) = \mathbf{0}, \\ \mathbf{X}(0, x) = \bar{\mathbf{X}}(x), \end{array} \right. \qquad (43) $$

where both the conserved quantity $\mathbf{X} := (\rho, q, \zeta)^T$ and the initial datum $\bar{\mathbf{X}}$ take values in

$$ \Omega_e := \{(\rho, q, \zeta)^T \in \mathbb{R}^3 : (\rho + \rho_\infty, q + q_\infty) \in \mathcal{U}_e, \zeta + c_\infty > 0\}, \qquad (44) $$

where $\mathcal{U}_e$ is defined in (10), the flux $G: \Omega_e \to \mathbb{R}^3$ is defined by

$$ G(\rho, q, \zeta) := \begin{pmatrix} (\zeta + c_{\infty}) F(\rho + \rho_{\infty}, q + q_{\infty}) \\ 0 \end{pmatrix}, \qquad (45) $$

and, with a slight abuse of notation, $\mathbf{0} := (0, 0, 0)^T$. Observe that $\Omega_e$ is open and $\mathbf{0} \in \Omega_e$ because by hypothesis $(\rho_\infty, q_\infty)^T \in \mathcal{U}_e$ and $c_\infty > 0$. Moreover, $G$ is $\mathbf{C}^1$ in $\Omega_e$ by (2)$_1$ and (3)$_1$. The eigenvalues and the corresponding eigenvectors of the system of conservation laws in (43) are

$$ \begin{align*} \lambda_0(\mathbf{X}) = 0, \qquad r_0(\mathbf{X}) &= \begin{pmatrix} V(h_\infty(\mathbf{X})) & (\rho + \rho_\infty) \\ V(h_\infty(\mathbf{X})) & (q + q_\infty) \\ ((\zeta + c_\infty) ((\rho + \rho_\infty) p'(\rho + \rho_\infty) V'(h_\infty(\mathbf{X}))) - V(h_\infty(\mathbf{X}))) & 0 \end{pmatrix}, \\ \lambda_1(\mathbf{X}) &= (\zeta + c_\infty) (V(h_\infty(\mathbf{X}))) - (\rho + \rho_\infty) p'(\rho + \rho_\infty) V'(h_\infty(\mathbf{X})), \qquad r_1(\mathbf{X}) = \begin{pmatrix} \rho + \rho_\infty \\ q + q_\infty \\ 0 \end{pmatrix}, \\ \lambda_2(\mathbf{X}) &= (\zeta + c_\infty) V(h_\infty(\mathbf{X})), \qquad r_2(\mathbf{X}) = \begin{pmatrix} \rho + \rho_\infty \\ q + q_\infty + (\rho + \rho_\infty)^2 p'(\rho + \rho_\infty) \\ 0 \end{pmatrix}, \end{align*} $$

where

$$ h_\infty(\mathbf{X}) := \frac{q + q_\infty}{\rho + \rho_\infty} - p(\rho + \rho_\infty). $$

Observe that $h_\infty(\mathbf{X}) > 0$ for any $\mathbf{X} \in \Omega_e$ by (44). System (43) is strictly hyperbolic in $\Omega_e$ because $\lambda_1(\mathbf{X}) < \lambda_2(\mathbf{X})$ by (2)$_3$ and (3)$_3$, $\lambda_1(\mathbf{X}) \neq 0 = \lambda_0(\mathbf{X})$ by (4)$_1$, and $\lambda_2(\mathbf{X}) \neq 0 = \lambda_0(\mathbf{X})$ by (3)$_{2,3}$. The first and the last characteristic fields are linearly degenerate because $\nabla\lambda_0 \cdot r_0 \equiv 0$ and $\nabla\lambda_2 \cdot r_2 \equiv 0$. The second characteristic field is genuinely non-linear because by (4)$_2$ we have

$$ \begin{align*} \nabla\lambda_1(\mathbf{X}) \cdot r_1(\mathbf{X}) &= \\ &= (\zeta + c_\infty)(\rho + \rho_\infty)((\rho + \rho_\infty)p'(\rho + \rho_\infty)^2 V''(h_\infty(\mathbf{X})) \\ &\quad - (2p'(\rho + \rho_\infty) + (\rho + \rho_\infty)p''(\rho + \rho_\infty))V'(h_\infty(\mathbf{X}))) \neq 0. \end{align*} $$

By [12, Theorem 7.1] there exists $\delta_0 > 0$ such that for every initial condition $\bar{\mathbf{X}}$ in $\mathbf{L}^1$ with $\text{TV}[\bar{\mathbf{X}}] \le \delta_0$, the Cauchy problem (43) has a weak solution $\mathbf{X}$ defined for all $t \ge 0$. In addition, if system (43) admits a convex entropy $\eta$, with entropy flux $\mu$, then for any test function $\varphi \in \mathbf{C}_c^\infty((0, \infty) \times \mathbb{R}; (0, \infty))$

$$ \int_{\mathbb{R}_+} \int_{\mathbb{R}} (\eta(\mathbf{X}) \varphi_t + \mu(\mathbf{X}) \varphi_x) \, dx \, dt \ge 0. \qquad (46) $$

We take now $\bar{\mathbf{X}}(x) := (\bar{\rho}(x) - \rho_\infty, \bar{q}(x) - q_\infty, c(x) - c_\infty)^T$. By assumption $\bar{\mathbf{X}} \in \mathbf{L}^1(\mathbb{R})$ and $\text{TV}[\bar{\mathbf{X}}] \le \delta_0$. Furthermore, $\bar{\mathbf{X}}$ takes values in $\Omega_e$ because $(\bar{\rho}, \bar{q})$ takes values in $\mathcal{U}_e$ and by (5)$_3$. The Cauchy
problem (43) corresponding to such initial datum has then a weak solution $\mathbf{X} := (\rho, q, \zeta)$ defined for all $t \ge 0$, namely $\mathbf{X}(0, x) = \bar{\mathbf{X}}(x)$ for a.e. $x \in \mathbb{R}$, and for any test function $\varphi \in \mathbf{C}_c^\infty((0, \infty) \times \mathbb{R}; \mathbb{R})$

$$ 
\begin{align*}
\int_{\mathbb{R}_+} \int_{\mathbb{R}} \left( \rho \varphi_t + (\zeta + c_\infty) V \left( \frac{q + q_\infty}{\rho + \rho_\infty} - p(\rho + \rho_\infty) \right) (\rho + \rho_\infty) \varphi_x \right) dx dt &= 0, \\
\int_{\mathbb{R}_+} \int_{\mathbb{R}} \left( q \varphi_t + (\zeta + c_\infty) V \left( \frac{q + q_\infty}{\rho + \rho_\infty} - p(\rho + \rho_\infty) \right) (q + q_\infty) \varphi_x \right) dx dt &= 0, \\
\int_{\mathbb{R}_+} \int_{\mathbb{R}} \zeta \varphi_t dx dt &= 0.
\end{align*}
$$

This implies that $\zeta = c - c_\infty$ a.e. in $\mathbb{R}$ and therefore $\mathbf{U} := (\rho + \rho_\infty, q + q_\infty)$ is a weak solution of (11) taking values in $\mathcal{U}_e$.

Let $(\mathcal{E}, c, \mathcal{Q})$ be a convex entropy pair for the system of conservation laws in (11) satisfying (12). Since $(\mathcal{E}, c, \mathcal{Q})$ is a convex entropy pair, for a.e. $x \in \mathbb{R}$ we have that

$$ D_{\mathbf{U}}\mathcal{E}(\mathbf{U}) \cdot D_{\mathbf{U}}(c(x) F(\mathbf{U})) = D_{\mathbf{U}}(c(x) \mathcal{Q}(\mathbf{U})), \quad (47) $$

$$ D_{\mathbf{U}}\mathcal{E}(\mathbf{U}) \cdot \partial_x (c(x) F(\mathbf{U})) = \partial_x (c(x) \mathcal{Q}(\mathbf{U})), \quad (48) $$

and $D_{\mathbf{U}}^2\mathcal{E}(\mathbf{U}) \ge 0$. Observe that (12) implies (48), but the converse is not true if $c$ is piecewise constant. By (5)$_3$, equation (47) is equivalent to

$$ D_{\mathbf{U}}\mathcal{E}(\mathbf{U}) \cdot D_{\mathbf{U}}F(\mathbf{U}) = D_{\mathbf{U}}\mathcal{Q}(\mathbf{U}). \quad (49) $$

Define

$$ \eta(\rho, q, \zeta) := \mathcal{E}(\rho + \rho_\infty, q + q_\infty), \quad \mu(\rho, q, \zeta) := (\zeta + c_\infty) \mathcal{Q}(\rho + \rho_\infty, q + q_\infty). $$

We have that $(\eta, \mu)$ is an entropy pair for the system of conservation laws in (43) because, by (45), (49), and (12), it holds that

$$ 
\begin{align*}
D_{\mathbf{X}}\eta(\mathbf{X}) \cdot D_{\mathbf{X}}G(\mathbf{X}) &= (D_{\mathbf{U}}\mathcal{E}, 0) \cdot \begin{pmatrix} ((\zeta + c_\infty) D_{\mathbf{U}}F & F \\ \mathbf{0} & 0 \end{pmatrix} \\
&= ((\zeta + c_\infty) D_{\mathbf{U}}\mathcal{E} \cdot D_{\mathbf{U}}F, D_{\mathbf{U}}\mathcal{E} \cdot F) = ((\zeta + c_\infty) D_{\mathbf{U}}\mathcal{Q}, \mathcal{Q}) = D_{\mathbf{X}}\mu(\mathbf{X}).
\end{align*}
$$

Moreover $\eta$ is convex because

$$ D_{\mathbf{X}}^2\eta(\rho, q, \zeta) = \begin{pmatrix} D_{\mathbf{U}}^2\mathcal{E}(\rho + \rho_\infty, q + q_\infty) & 0 \\ 0 & 0 \end{pmatrix} \ge 0. $$

With the above choice for $(\eta, \mu)$ and recalling that $\zeta = c - c_\infty$ a.e. in $\mathbb{R}$, the entropy condition (46) becomes (13). This ends the proof of Theorem 2.1.

## 6.2. Proof of Proposition 3.10

In this section we need to represent the construction of the solutions to some Riemann data in the $(\rho, f)$-plane. Observe that $\mathbf{W} = (h, w)^T \in \mathcal{W}$ is represented in the $(\rho, f)$-plane by the two points $(\rho, f^-) = (p^{-1}(w-h), c^- f(\mathbf{W}))$ and $(\rho, f^+) = (p^{-1}(w-h), c^+ f(\mathbf{W}))$. This motivates the introduction of the following notation. For a given $\mathbf{W}_L \in \mathcal{W}$, take

$$ F_L^- := (p^{-1}(w_L - h_L), c^- f(\mathbf{W}_L)) \qquad F_L^+ := (p^{-1}(w_L - h_L), c^+ f(\mathbf{W}_L)). $$

Analogous notation is introduced for $F_R^\pm, F_*^\pm$, and so on.
1) **Coherence.** We prove now that the Riemann solver $RS_{c^-,c^+}$ is coherent in $W^2$. Observe that to do so it is sufficient to prove

$$Q(\mathbf{W}_{-}, \mathbf{W}_{+}) = Q(\mathbf{W}_{L}, \mathbf{W}_{R}). \qquad (50)$$

Indeed, (50) implies that either $\mathbf{W}_{-} = \hat{\mathbf{W}}(\mathbf{W}_{-}, \mathbf{W}_{+})$ or the 1-shock $S_1(\mathbf{W}_{-}, \hat{\mathbf{W}}(\mathbf{W}_{-}, \mathbf{W}_{+}))$ has zero speed of propagation, as well as either $\mathbf{W}_{+} = \tilde{\mathbf{W}}(\mathbf{W}_{-}, \mathbf{W}_{+})$ or the 1-shock $S_1(\tilde{\mathbf{W}}(\mathbf{W}_{-}, \mathbf{W}_{+}), \mathbf{W}_{+})$ has zero speed of propagation. In any of these cases, we have $RS_{c^-,c^+}[\mathbf{W}_{-}, \mathbf{W}_{+}](0^{\pm}) = \mathbf{W}_{\pm}$.

Assume by contradiction that there exists $(\mathbf{W}_L, \mathbf{W}_R) \in W^2$ such that (50) is not satisfied, namely

$$Q(\mathbf{W}_{-}, \mathbf{W}_{+}) \neq Q(\mathbf{W}_{L}, \mathbf{W}_{R}).$$

By (26), (29) and (30) we have $w_*(\mathbf{W}_R, w_L) = \hat{w}(\mathbf{W}_L, \mathbf{W}_R) = \tilde{w}(\mathbf{W}_L, \mathbf{W}_R) = w_L$, hence by Proposition 3.5

$$w_{-} = w_{L} = w_{+}. \qquad (51)$$

This implies that $\mathbf{W}_*(\mathbf{W}_+, w_-) = \mathbf{W}_+$. By (27), (28) and (31) we have then

$$Q(\mathbf{W}_{-}, \mathbf{W}_{+}) > c^{\pm} f(\mathbf{W}_{\pm}) = Q(\mathbf{W}_{L}, \mathbf{W}_{R}).$$

As a consequence, by (51) we have

$$Q(\mathbf{W}_{-}, \mathbf{W}_{+}) = \min\{c^{-}, c^{+}\} F(w_{L}) > c^{\pm} f(\mathbf{W}_{\pm}) = Q(\mathbf{W}_{L}, \mathbf{W}_{R}),$$

and therefore

$$Q(\mathbf{W}_L, \mathbf{W}_R) < \min\{c^-, c^+\} F(w_L), \quad h_- < h(w_L), \quad h_+ = h_*(\mathbf{W}_+, w_-) > h(w_L).$$

This implies that

$$\mathbf{W}_{-} = \hat{\mathbf{W}}(\mathbf{W}_{L}, \mathbf{W}_{R}), \qquad \mathbf{W}_{+} = \tilde{\mathbf{W}}(\mathbf{W}_{L}, \mathbf{W}_{R}) = \mathbf{W}_{*}(\mathbf{W}_{+}, w_{-}).$$

The above considerations imply that $Q_-(\mathbf{W}_L) \ge Q(\mathbf{W}_L, \mathbf{W}_R)$, with the equality holding if and only if $\mathbf{W}_L = \mathbf{W}_- = \hat{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)$. However, this means that

$$Q(\mathbf{W}_L, \mathbf{W}_R) = c^+ f(\mathbf{W}_*(\mathbf{W}_R, w_L)), \qquad h_*(\mathbf{W}_R, w_L) < h(w_L).$$

On the other hand, the above conditions imply that $\mathbf{W}_+ = \mathbf{W}_*(\mathbf{W}_R, w_L)$ and therefore

$$h_+ = h_*(\mathbf{W}_R, w_L) < h(w_L) < h_*(\mathbf{W}_+, w_-) = h_+.$$

This gives a contradiction and therefore $RS_{c^-,c^+}$ is coherent.

2) **$L^1_{\text{loc}}$-continuity.** We prove now that the Riemann solver $RS_{c^-,c^+}$ is $L^1_{\text{loc}}$-continuous with respect to the initial datum in $W^2$. Observe first that $Q$ is continuous because $Q_-$ and $Q_+$ are so. By the continuity of the speed of propagation of the waves, it is sufficient to consider initial data for which traces at $\nu = 0$ are discontinuous. More precisely, assume that $(\mathbf{W}_L^\epsilon, \mathbf{W}_R^\epsilon)$ converges to $(\mathbf{W}_L, \mathbf{W}_R)$ (as $\epsilon$ tends to zero), but the traces $\mathbf{W}_\pm := RS_{c^-,c^+}[\mathbf{W}_L^\epsilon, \mathbf{W}_R^\epsilon](0^\pm)$ do not converge to $\mathbf{W}_\pm := RS_{c^-,c^+}[\mathbf{W}_L, \mathbf{W}_R](0^\pm)$. A case by case study shows that $RS_{c^-,c^+}[\mathbf{W}_L^\epsilon, \mathbf{W}_R^\epsilon]$ has then a unique 1-wave, either in $\nu < 0$ or in $\nu > 0$, and a non-classical shock at $\nu = 0$. In both cases, the $L^1_{\text{loc}}$-continuity follows by observing that the speed of propagation of such 1-wave is continuous and goes to zero as $\epsilon$ tends to zero. Here we omit the details of the proof and defer the reader to Figure 7 for the case $c_- > c_+$, the case $c_- < c_+$ is analogous. Observe that in the case as in Figure 7, left, $RS_{c^-,c^+}[\mathbf{W}_L^\epsilon, \mathbf{W}_R]$ and $RS_{c^-,c^+}[\mathbf{W}_L^\epsilon, \mathbf{W}_R]$ differ only in $\nu \le 0$: the former is the juxtaposition of the 1-shock $S_1(\mathbf{W}_L^\epsilon, \hat{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R))$ and the non-classical shock $N(\hat{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R), \mathbf{W}_*(\mathbf{W}_R, w_L))$, whereas the latter is just the non-classical shock $N(\mathbf{W}_L, \mathbf{W}_*(\mathbf{W}_R, w_L))$. In this case, the $L^1_{\text{loc}}$-convergence of $RS_{c^-,c^+}[\mathbf{W}_L^\epsilon, \mathbf{W}_R]$ to $RS_{c^-,c^+}[\mathbf{W}_L, \mathbf{W}_R]$
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Two representative examples of discontinuous traces. We let $\mathbf{W}_* = \mathbf{W}_*(\mathbf{W}_R, w_L)$, $\tilde{\mathbf{W}} = \tilde{\mathbf{W}}(\mathbf{W}_L, \mathbf{W}_R)$ and $\tilde{\mathbf{W}}^\varepsilon = \tilde{\mathbf{W}}(\mathbf{W}_L^\varepsilon, \mathbf{W}_R)$. We omit the representation of $\mathbf{F}_R$ as it is not strictly necessary.</figcaption></figure>

holds because the speed of propagation of the 1-shock goes to zero as $\varepsilon$ tends to zero. The case as in Figure 7, right, is analogous.

3) $\mathbf{W}_e$. Observe that $RS_{c^-,c^+}[\mathbf{W}_L, \mathbf{W}_R]$ involves a vacuum state if and only if $\mathbf{W}_*(\mathbf{W}_R, w_L)$ is a vacuum state. By (26) this can happen if and only if $\mathbf{W}_*(\mathbf{W}_R, w_L) = (w_L, w_L)$, namely $w_L < h_R$. Furthermore, if $h_L, h_R > 0$ then by construction we have $Q(\mathbf{W}_L, \mathbf{W}_R) > 0$, hence the $h$-component of $RS_{c^-,c^+}[\mathbf{W}_L, \mathbf{W}_R]$ does not vanish.

### Acknowledgment

F.A.C., S.F. and M.D.R. are members of Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM). F.A.C. and S.F. are partially supported by the Ministry of University and Research (MUR), Italy, under the grant PRIN 2020 - Project N. 20204NT8W4, "Non-linear evolution PDEs, fluid dynamics and transport equations: theoretical foundations and applications". F.A.C. recognizes financial support from the Project "Sistemi iperbolici: modelli di traffico veicolare e modelli idrodinamici di tipo flocking", Progetti di Ateneo per Avvio alla Ricerca 2023, Universiy of L'Aquila. S.F. is partially supported by the Italian "National Centre for HPC, Big Data and Quantum Computing" - Spoke 5 "Environment and Natural Disasters", by the InterMaths Network, www.intermaths.eu. M.D.R. recognizes financial support from the PRIN 2022 Project "Modeling, Control and Games through Partial Differential Equations", CUP D53D23005620006, funded by the European Union-Next Generation EU, and from "INDAM - GNAMPA Project", CUP E53C23001670001.

### References

[1] F. Ancona, A. Cesaroni, G. M. Coclite, and M. Garavello. On the optimization of conservation law models at a junction with inflow and flow distribution controls. *SIAM Journal on Control and Optimization*, 56(5):3370–3403, Jan. 2018.

[2] F. Ancona and P. Goatin. Uniqueness and stability of $L^\infty$ solutions for Temple class systems with boundary and properties of the attainable sets. *SIAM J. Math. Anal.*, 34(1):28–63, 2002.

[3] B. Andreianov. New approaches to describing admissibility of solutions of scalar conservation laws with discontinuous flux. In *CANUM 2014—42e Congrès National d'Analyse Numérique*, volume 50 of *ESAIM Proc. Surveys*, pages 40–65. EDP Sci., Les Ulis, 2015.
[4] B. Andreianov, C. Donadello, and M. D. Rosini. A second-order model for vehicular traffics with local point constraints on the flow. *Mathematical Models and Methods in Applied Sciences*, 26(04):751-802, 2016.

[5] A. Aw and M. Rascle. Resurrection of “second order” models of traffic flow. *SIAM Journal on Applied Mathematics*, 60(3):916-938, 2000.

[6] P. Baiti and A. Bressan. The semigroup generated by a Temple class system with large data. *Differential Integral Equations*, 10(3):401-418, 1997.

[7] F. Betancourt, R. Bürger, C. Chalons, S. Diehl, and S. Farås. A random sampling method for a family of temple-class systems of conservation laws. *Numerische Mathematik*, 138(1):37-73, July 2017.

[8] F. Betancourt, R. Bürger, S. Diehl, and S. Farås. Modeling and controlling clarifier-thickeners fed by suspensions with time-dependent properties. *Minerals Engineering*, 62:91-101, July 2014.

[9] S. Bianchini. The semigroup generated by a Temple class system with non-convex flux function. *Differential Integral Equations*, 13(10-12):1529-1550, 2000.

[10] R. Borsche, M. Garavello, and B. Gunarso. Temple system on networks. *J. Hyperbolic Differ. Equ.*, 20(3):669-705, 2023.

[11] R. Borsche and A. Klar. A nonlinear discrete velocity relaxation model for traffic flow. *SIAM Journal on Applied Mathematics*, 78(5):2891-2917, Jan. 2018.

[12] A. Bressan. *Hyperbolic Systems of Conservation Laws. The one-dimensional Cauchy problem.*, volume 20. Oxford Lecture Series in Mathematics and its Applications. Oxford University Press, 2000.

[13] A. Bressan and P. Goatin. Stability of $L^\infty$ solutions of Temple class systems. *Differential Integral Equations*, 13(10-12):1503-1528, 2000.

[14] F. Chiarello, S. Göttlich, T. Schillinger, and A. Tosin. Hydrodynamic traffic flow models including random accidents: A kinetic derivation. *Commun. Math. Sci.*, 2023.

[15] F. A. Chiarello. *Conservation laws and vanishing viscosity on network*. Master thesis, University of Bari (Italy) and Université de Franche-Comté (France), March 2017.

[16] G. M. Coclite, M. Garavello, and B. Piccoli. Traffic flow on a road network. *SIAM Journal on Mathematical Analysis*, 36(6):1862-1886, Jan. 2005.

[17] G. M. Coclite and N. H. Risebro. Conservation laws with time dependent discontinuous coefficients. *SIAM Journal on Mathematical Analysis*, 36(4):1293-1309, 2005.

[18] R. M. Colombo and A. Corli. On a class of hyperbolic balance laws. *J. Hyperbolic Differ. Equ.*, 1(4):725-745, 2004.

[19] R. M. Colombo and P. Goatin. A well posed conservation law with a variable unilateral constraint. *J. Differential Equations*, 234(2):654-675, 2007.

[20] R. M. Colombo and M. D. Rosini. Well posedness of balance laws with boundary. *J. Math. Anal. Appl.*, 311(2):683-702, 2005.

[21] R. M. Colombo and M. D. Rosini. Well posedness of balance laws with non-characteristic boundary. *Boll. Unione Mat. Ital. Sez. B Artic. Ric. Mat.* (8), 10(3):875-894, 2007.

[22] A. Corli, M. Figiel, A. Futa, and M. D. Rosini. Coupling conditions for isothermal gas flow and applications to valves. *Nonlinear Anal. Real World Appl.*, 40:403-427, 2018.

[23] A. Corli, U. Razafison, and M. D. Rosini. Coherence and flow-maximization of a one-way valve. *ESAIM Math. Model. Numer. Anal.*, 56(5):1715-1739, 2022.

[24] C. Dafermos. *Hyperbolic conservation laws in continuum physics*, volume 325 of *Grundlehren der mathematischen Wissenschaften*. Springer-Verlag, Berlin, third edition, 2010.

[25] S. Fan, M. Herty, and B. Seibold. Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model. *Networks and Heterogeneous Media*, 9(2):239-268, 2014.
[26] R. E. Ferreira and C. I. Kondo. Glimm method and wave-front tracking for the Aw-Rascle traffic flow model. *Far East J. Math. Sci. (FJMS)*, 43(2):203-223, 2010.

[27] M. D. Francesco, S. Fagioli, and M. D. Rosini. Many particle approximation of the Aw-Rascle-Zhang second order model for vehicular traffic. *Mathematical Biosciences and Engineering*, 14(1):127-141, 2017.

[28] M. Godvik and H. Hanche-Olsen. Existence of solutions for the Aw-Rascle traffic flow model with vacuum. *J. Hyperbolic Differ. Equ.*, 5(1):45-63, 2008.

[29] H. Holden, N. H. Risebro, and H. Sande. The solution of the Cauchy problem with large data for a model of a mixture of gases. *J. Hyperbolic Differ. Equ.*, 6(1):25-106, 2009.

[30] H. Holden, N. H. Risebro, and H. Sande. Front tracking for a model of immiscible gas flow with large data. *BIT*, 50(2):331-376, 2010.

[31] K. H. Karlsen, N. Risebro, and J. D. Towers. $L^1$ stability for entropy solutions of nonlinear degenerate parabolic convection-diffusion equations with discontinuous coefficients. *Skr. K. Nor. Vidensk. Selsk*, 3:1-49, 2003.

[32] K. H. Karlsen and J. D. Towers. Convergence of the Lax-Friedrichs scheme and stability for conservation laws with a discontinuous space-time dependent flux. *Chinese Annals of Mathematics*, 25(03):287-318, jul 2004.

[33] R. J. LeVeque. *Finite volume methods for hyperbolic problems*. Cambridge Texts in Applied Mathematics. Cambridge University Press, Cambridge, 2002.

[34] Y.-G. Lu. Existence of global bounded weak solutions to nonsymmetric systems of Keyfitz-Kranzer type. *J. Funct. Anal.*, 261(10):2797-2815, 2011.

[35] S. Mishra. Numerical methods for conservation laws with discontinuous coefficients. In *Handbook of numerical methods for hyperbolic problems*, volume 18 of *Handb. Numer. Anal.*, pages 479-506. Elsevier/North-Holland, Amsterdam, 2017.

[36] A. P. Pires, P. G. Bedrikovetsky, and A. A. Shapiro. A splitting technique for analytical modelling of two-phase multicomponent flow in porous media. *Journal of Petroleum Science and Engineering*, 51(1):54-67, 2006. Special Issue on Petroleum Production Research in Brazil.

[37] M. D. Rosini. Systems of conservation laws with discontinuous fluxes and applications to traffic. *Ann. Univ. Mariae Curie-Skłodowska Sect. A*, 73(2):135-173, 2019.

[38] D. Serre. Richness and the classification of quasilinear hyperbolic systems. In *Multidimensional hyperbolic problems and computations (Minneapolis, MN, 1989)*, volume 29 of *IMA Vol. Math. Appl.*, pages 315-333. Springer, New York, 1991.

[39] W. Shen. Global Riemann solvers for several 3 × 3 systems of conservation laws with degeneracies. *Math. Models Methods Appl. Sci.*, 28(8):1599-1626, 2018.

[40] B. Temple. Systems of conservation laws with coinciding shock and rarefaction curves. *Contemp. Math.*, 17(143):151, 1983.

[41] B. Temple. Systems of conservation laws with invariant submanifolds. *Trans. Amer. Math. Soc.*, 280(2):781-795, 1983.

[42] D. Wagner. Equivalence of the Euler and Lagrangian equations of gas dynamics for weak solutions. *J. Differential Equations*, 68(1):118-136, 1987.

[43] H. Zhang. A non-equilibrium traffic model devoid of gas-like behavior. *Transportation Res. Part B*, 36(3):275-290, 2002.

Felisia Angela Chiarello
University of L'Aquila,
Department of Engineering and Information Science and Mathematics,
Via Vetoio, Ed. Coppito 1, L'Aquila, 67100, Italy
e-mail: felisiaangela.chiarello@univaq.it
Simone Fagioli
University of L'Aquila,
Department of Engineering and Information Science and Mathematics,
Via Vetoio, Ed. Coppito 1, L'Aquila, 67100, Italy
e-mail: simone.fagioli@univaq.it

Massimiliano Daniele Rosini
Maria Curie-Skłodowska University,
Plac Marii Curie-Skłodowskiej 1,
Lublin, 20-031, Poland,
Department of Management and Business Administration,
University "G. d'Annunzio" of Chieti-Pescara,
viale Pindaro, 42,
Pescara, 65127, Italy
e-mail: massimiliano.rosini@unich.it