Asymptotic problems and numerical schemes for traffic flows with unilateral constraints describing the formation of jams

We discuss numerical strategies to deal with PDE systems describing traffic flows, taking into account a density threshold, which restricts the vehicles density in the situation of congestion. These models are obtained through asymptotic arguments. Hence, we are interested in the simulation of approached models that contain stiff terms and large speeds of propagation. We design schemes intended to apply with relaxed stability conditions.


Introduction
In order to describe traffic flows and to reproduce the formation of congestions, several models based either on Ordinary Differential Equations (ODE) or Partial Differential Equations (PDE) have been proposed. Starting from individual-based "Follow-the-Leader" models [20], a very active stream in the traffic community considers now PDE models. A first example dates back to Lighthill and Whitham in the 50's [25]: the evolution of the density of cars is described by means of a mass conservation equation, where the flux is defined by a prescribed function of the density. In these so-called first-order models, the relation between flux and density is referred to as the fundamental diagram in the traffic flows community. A more accurate description can be expected by considering second-order models where a system of PDE governs the evolution of the density and the speed of cars. A first attempt in this direction is due to Payne [31], strongly inspired by the principles of fluid mechanics. However, Daganzo [15] pointed out the drawbacks of this approach: the Payne-Whitham model may lead to inconsistent behaviors for the flow, such as vehicles going backwards. The model introduced independently by Aw and Rascle [2] and by Zhang [37], which still has the form of a 2 × 2 system of conservation laws, is intended to correct these inconsistencies. In [1], a derivation of the system is proposed from a Follow-the-Leader model. We can also mention that some kinetic models [3,29,30,33,36] are under consideration, after the pioneering work [32].
Further details and references can be found in the survey [4].
This work is concerned with the numerical simulation of certain variants of the Aw-Rascle-Zhang model. Let ρ(x, t) and v(x, t) be the density and the velocity of cars at position x ∈ R and time t > 0, respectively. The Aw-Rascle-Zhang model writes where ρ → p(ρ) plays the role of the pressure in the gas dynamic equations.
In fact, the quantity w = v+p(ρ) describes the desired velocity of the drivers, whereas v corresponds to the actual velocity of the cars. Therefore, the term p(ρ), the velocity offset, stands for the difference between these two velocities, reflecting the fact that the drivers slow down because of the density of cars.
It is convenient to rewrite the equations in the more convenient form of a conservative system; namely (1) is, at least formally, equivalent to    ∂ t ρ + ∂ x (ρv) = 0, ∂ t ρ(v + p(ρ)) + ∂ x ρv(v + p(ρ)) = 0. (2) Of course, a crucial modeling issue relies on the expression of the velocity offset p(ρ). At first glance, again inspired from gas dynamics, we can set p(ρ) = ρ γ for some γ > 1. However, such a model does not permit to impose a priori a limitation to the cars density. Consequently, Berthelin, Degond, Delitala and Rascle proposed in [7] to define the velocity offset as follows: which can be rewritten, for ρ = 0, by p(ρ) , where ρ denotes a maximal value for the density. The velocity offset tends to infinity when ρ → ρ while we get the classical expression p(ρ) ∼ ρ γ when ρ → 0. Such a pressure law also arises in gas dynamics, where it is referred to as the Bethe-Weyl law [5]; for instance it is used to model close-packing effects in multifluid flows, see [8] and the references therein. This expression for p has the role of enforcing that ρ satisfies the constraint 0 ≤ ρ(x, t) ≤ ρ for all x ∈ R and t > 0, as it can be seen from the bounds on the Riemann invariants of the system [7]. Moreover, [7] points out that drivers do not reduce significantly their speed unless they reach a congested region. Accordingly the velocity offset is appropriately rescaled with a small parameter 0 < ε 1, and we use in the Aw-Rascle-Zhang model. We are thus led to the Rescaled Modified Aw-Rascle (RMAR) system In this model, the velocity offset is small unless the density is getting close to the threshold ρ . Finally, [7] studies the limit when ε → 0 in (3). In this regime we obtain (at least formally) the constrained system ∂ t (ρ(v + π)) + ∂ x (ρv(v + π)) = 0, In (4), the limit "pressure" π = lim ε→0 p ε (ρ ε ) appears as the Lagrange multiplier associated to the unilateral constraint 0 ≤ ρ ≤ ρ . In particular, π becomes active only in the congested regions, where the density reaches the threshold ρ . Otherwise, in absence of congestion, the system is reduced to the pressure-less gas dynamic model [11,22]    ∂ t ρ + ∂ x (ρv) = 0, ∂ t (ρv) + ∂ x (ρv 2 ) = 0.
The asymptotic model is further investigated in [7], exhibiting the formation of clusters, and proving the existence of weak solutions to the system (4) through the stability analysis of "sticky blocks" dynamics. It is also worth pointing out the original numerical approach developed in [28] for (4) which uses ideas from the modelling of crowd motion and includes a fine description of the non elastic collision processes.
The asymptotic system (4) is thus specifically intended to describe the formation and the dynamics of jams. In this paper, we are interested in the numerical simulations of the system (4), and in the asymptotic regime ε → 0 in (3). The difficulty is two-fold. On the one hand, in the free flow case, it is well-known that the pressureless gas dynamics system (5) can lead to deltashocks formation, which makes it difficult to treat numerically [10]. On the other hand, with the formation of a congestion, there is no direct access to the limit velocity offset π which is defined in a quite abstract way. Therefore, in order to go beyond the simple particulate approach in [7], we wish to develop numerical simulations of the RMAR model (3) with the velocity offset (VO1) for small values of ε. We are still facing several numerical challenges.
Firstly, the model prohibits that the density exceeds the threshold ρ . Secondly, (one of) the characteristic speeds of the system become very large in congested region, which makes the time step shrink: the smaller ε, the more severe the stability constraint. Therefore, we need to design a scheme which can preserve the natural estimates of the problem, in particular the density limitation. As already observed in [12,13] standard schemes may fail this objective due to the very specific structure of the PDE system. We also refer the reader to [23] for further examples related to fluid mixtures. Moreover, we would like to relax the stability constraints on the time step. To this end, a first attempt would be to adapt the explicit-implicit method pro-posed in [18] for the Euler system with congestion constraint. However, this method applied to (3) is not satisfactory: it produces excessively smoothed density profiles and it overestimates the velocity in congested regions. Thus we design a different splitting strategy, partly inspired from [19]. Beyond the conception of a numerical scheme able to handle the stiffness of the problem, we will also discuss different asymptotic approaches of the constrained problem (4), based on different definitions of the scaled offset velocity in (3) which all lead asymptotically to (4). It is interesting to study how the shape of the pseudo-pressure affects the intermediate states (for not so extreme values of the scaling parameter), and the numerical costs.
The outline of this article is the following. In Section 2, we go back to some properties of the Aw-Rascle-Zhang system and we detail the numerical difficulties we face. Additionally, we propose different velocity offsets and scaling that can be used to recover asymptotically the constrained system (4). Then, in Section 3, we propose a new explicit-implicit scheme based on a splitting strategy. The splitting is constructed to reduce the characteristic speeds in the explicit part so that we can expect to use larger time steps.
Finally, in Section 4, we display some numerical simulations in order to prove the efficiency of the scheme and to compare the behavior of the system when using different velocity offsets.
2 Properties of the Aw-Rascle-Zhang model and numerical difficulties We will describe in this Section the main numerical difficulties we have to deal with, when computing solutions of system (3).

Different velocity offsets
With the velocity offset (VO1), it is forbidden to produce numerical densities larger than the threshold ρ : if, due to any numerical error, the code returns a density larger than ρ , we cannot update further the system and the simulation breakdowns. To cope with this difficulty, we propose to slightly modify the law, replacing (VO1) by a function ρ →p ε (ρ) which is defined for any positive entry, which behaves like p ε for ρ < ρ , and which blows up In this formula, ρ ε tr is a transition density, which has a modeling nature; it should satisfy ρ ε tr → ρ as ε tends to 0. Beyond the transition,p ε is a second order polynomial, computed so thatp ε remains a C 2 function. We thus set The expected behavior holds for instance with h(ε) = ε, since it satisfies the two following properties when ε → 0: h(ε) → 0 and c 0 ε → +∞. We point out that ρ ε tr is purely a modeling parameter and the model (VO2) leads us to the same difficulties as (VO1) in terms of stability issues, as we will see in the next Section.
Another option is to use the following velocity offset for large values of the exponent γ. In this formula V ref > 0 is a reference velocity, to bear in mind the physical meaning of p (in the numerical simulations below we will simply set V ref = 1). This approach is used in fluid mechanics, for modeling certain free boundary problems where bubbles are immersed in a gas [26]. This function is defined on [0, ∞), it behaves proportionally to the gas-law ρ γ for ρ → 0 and it blows up as γ → +∞ for ρ ≥ ρ . Using (VO3) and the regime of large γ's in traffic flows modeling is quite new; we shall see that this simple law has certain advantages in the numerical simulations of congested situations.

Stability issues
As long as the functions ρ and v are smooth enough, we can rewrite system (1) in the fully non-conservative form The two eigenvalues related to the system are therefore equal to with related eigenvectors The system is strictly hyperbolic, away from the regions where ρ = 0. Let us just note that the largest eigenvalue λ 2 is always linearly degenerate, leading to contact discontinuities and that λ 1 is genuinely non-linear, except for certain forms of the velocity offset we are not considering here. Therefore, the first eigenvalue will admit shocks or rarefaction waves. One of the difficulties of the computations is that vacuum regions may appear. Observe that the information does not travel faster than the actual cars speed v, and that the system preserves the natural properties ρ ≥ 0, v ≥ 0.
We are considering here some Finite Volume (FV) numerical schemes in order to compute the solutions of system (1). Let us denote by ∆t and ∆x the time step and the space step of the method, respectively. We consider the discrete times t n = n∆t, for n ∈ N and the discretization cells , j ∈ Z (neglecting for the time being the issue of the boundary conditions) where x j+1/2 = (j + 1/2)∆x. We go back to the conservative form (2) and we denote the conservative variables. In terms of the conservative variables ρ and y, which recasts as follows, using only the variables ρ and y, The numerical unknown U n j = (ρ n j , y n j ) is thought of as an approximation of the mean value of U (x, t) on the cell C j at time t n . The FV scheme has the following general form which mimics what we obtain by integrating the continuous equation (7) over t n , t n+1 × C j . For the simple schemes we wish to deal with, the numerical flux at the interface x j+1/2 is a function of the neighbouring cells F n j+1/2 = F U n j+1 , U n j . Without entering into the details of the schemes, the numerical stability of such a method relies (at least) on the following inequality, see for example [9,Sect. 2.3.3] or [35], where λ 1 and λ 2 are the two eigenvalues given by (6). It imposes that the state U x j+1/2 , t n+1 on the interface x j+1/2 at time t n+1 = t n + ∆t only depends on the states of the unknown at time t n on the neighbouring cells C j and C j+1 . It can thus be obtained by solving the corresponding Riemann problem with data U n j and U n j+1 .
Let us first consider the case of the pressure (VO1). In case of a congestion formation, ρ ε → ρ but we expect that p ε (ρ ε ) remains bounded and admits the limit π as ε → 0; it leads to the ansatz Accordingly the behavior of the characteristic speeds is given by since v ε should remain bounded when ε → 0. Hence, as ε goes to zero, the time step ∆t shrinks due to the condition (8). The same remarks apply to the velocity offset (VO2), which essentially behaves like (VO1) when ε → 0.
For the velocity offset (VO3), we find max λ which again imposes tiny time steps. This observation motivates the design of a scheme based on splitting strategy so that the fast waves can be treated implicitly.

Invariant regions
Let us detail another difficulty which is very specific to the traffic flow system (1). The Riemann invariants for the system (1) are given by, see [2], Therefore, the domain is an invariant region for (1): if the initial datum lies in such a region, the solution will still be contained in the same region for all times. However, numerical difficulties arise due to the fact that such domains are non-convex for the conserved quantities ρ, y. This point has been observed in [12] for the traffic flows model, see also [13,23] for similar problems. At first sight, it would be tempting to define the numerical fluxes by using the Godunov scheme, which is a standard for systems of conservation laws. It works into two steps. Owing to the stability condition (8), we solve a set of uncoupled Riemann problems, centered at the interfaces x j+1/2 with data U n j , U n j+1 .
Then, we project the obtained piecewise constant solution to obtain the updated numerical unknown, constant on the cells C j . This projection step does not preserve the invariant region, since the latter is non-convex (for the role of the convexity of the invariant domain we refer the reader to [24] for gas dynamics equations, and more generally to the textbook [9, Prop. 2.11]).
A counter-example is detailed in [12] to explain why the Godunov scheme fails to satisfy the maximum principle for the Riemann invariants of (1), and especially for v. To solve this difficulty, [12] proposes a hybrid scheme which mixes the Glimm scheme to compute the contact discontinuities, and the Godunov procedure to compute 1-shock or 1-rarefaction wave. This hybrid method is well-adapted to handle the specific velocity offsets dealt with in [12,13], see Remark 2 below, which differ from the models we wish to consider here.
We bear in mind that, instead of using the mean of the solutions of the Riemann problems over the cells, the Glimm scheme uses a random sampling strategy in the reconstruction procedure. Hence, by construction, the Glimm scheme preserves the invariant regions, despite the defect of convexity. It is thus well adapted to the simulation of the system (1). Note however the final scheme is non conservative.
Remark 1 Note that, depending on the definition of the numerical fluxes, the stability condition can be even more constrained than (8), for instance in order to fulfill the bound from above on the density with the "close-packinglike" velocity offset (VO1), see e.g. [8].
Remark 2 In [12,13], the discussion focuses on the velocity offset which has some very specific features: • First of all, ρ → p(ρ) is defined on (0, ∞) and non decreasing. The model cannot treat vacuum regions since the velocity offset is not defined for ρ = 0. With this model, the velocity offset is well-defined beyond ρ , and we note that p(ρ) < 0 for 0 < ρ < ρ , which might be physically questionable.
In particular, we get so that the CFL condition depends only on v and it does not shrink as the density approaches ρ .

Description of the scheme
Let us now explain in more details the construction of the scheme that we wish to use for the simulation of (1), with a velocity offset ρ → p(ρ) that introduces some stiffness in order to reproduce the expected behavior of the constrained system (4). In what follows, p thus refers to p ε in (VO1), top ε in (VO2) or to p γ in (VO3).
In order to get rid of the large characteristic speeds, the idea consists in splitting the velocity offset into two parts p = p exp + p imp , so that the system with p exp has a stability CFL condition (8) of order 1 with respect to the scaling parameters. Namely, the eigenvalue λ exp,1 ) is bounded with respect to 0 < ε 1 (resp. γ 1). The corresponding system can thus be treated explicitly by means of the Glimm scheme, which preserves the invariant domains. Next, only the stiff part that involves p imp is treated implicitly. We expect also that the implicit part has a simple structure that can be handled with a not too complicated scheme.
Such splitting approach appeared in [17,18,19] for more standard fluid mechanics systems, for instance as an efficient strategy to handle low Mach number regimes. However, as explained above, the structure of the PDE system (1) significantly differs from the Euler equations, in particular with the lack of convexity of the invariant domains and these methods cannot be directly applied to (1).

Definition of the explicit velocity offset
In the first step of the splitting, we consider the system It has the same structure as the Aw-Rascle-Zhang system (2), just replacing the full velocity offset p, that can be (VO1), (VO2) or (VO3), by p exp . The characteristic speeds are In order to relax the stability constraint, we define p exp so that the characteristic speed does not blow up as the scaling parameters ε goes to 0 or γ goes to ∞. It leads to require that p exp (ρ) is bounded uniformly with respect to ε (resp. γ for the law (VO3)), when ρ lies in a compact set of (0, ∞).
The definition depends on a truncation parameter 0 < ρ num < ρ . Following an idea of [17,18], we set We use a second order polynomial for ρ > ρ num to ensure that p exp is a C 2 function. We will see below that it is important to reach such a regularity.
The transition density ρ num is chosen so that p exp (ρ num ) is bounded when ε → 0 (resp. γ → +∞). If such a condition is satisfied, then the two eigenvalues are bounded with respect to ε (resp. γ) and therefore the stability condition (8) is independent of ε (resp. γ). Let us now explain how to choose ρ num for the velocity offsets under consideration.
a) For the laws (VO1) and (VO2), when ε → 0 and ρ ε → ρ , we expect which leads us to set For the velocity offset (VO2), we point out that ρ num is a truncation parameter of numerical nature, designed to ensure a certain stability property of the scheme, while ρ trans relies on modeling consideration and does not prevent the blow up of the velocity offset. In any case, we have 0 < ρ num < ρ trans < ρ . b) For the law (VO3), we require that p exp (ρ num ) remains bounded when γ → +∞. Denoting again ρ num = ρ (1 − δρ), a simple computation leads to So, we are looking for δρ such that δρ → 0 and both terms in this sum are O(1) when γ → +∞. A simple study of the two sequences of functions (γ exp (−γx)) γ∈N and γ 2 x exp (−γx) γ∈N shows that we should have δρ = O(γ −α ) with α ∈ (0, 1) in order to satisfy the required properties.
In [19] the pressure term is split as p = p exp + p imp with p exp = αp and Here, due to the singularity of p at ρ , this approach does not permit to keep p exp bounded and we have to define p exp and p imp as explained above.

A time-splitting scheme
As said above, it is convenient to work on the conservative form (7) of the system (1), dealing with the unknowns U = ρ, y where y = ρ(v + p(ρ)).
We use now a time-splitting scheme. Knowing some approximate values U n = (ρ n , y n ) at time t n , we proceed as follows.
• Step 1: Solve with an explicit scheme the system of conservation laws As said above, this system has the same structure as the original problem (7). In particular the invariant domains are non-convex. It can be solved with the Glimm scheme adapted for the pressure p exp . More details will be given in Section 3.3. It defines some intermediate values • Step 2: Solve implicitly the system Note that the system has a simple structure and the two equations decouple. The first equation is a non linear scalar conservation law for the density ρ, the second is a linear transport equation for y where the velocity −p imp (ρ) can be considered as given. The numerical method to solve the system (12) is thus not that complicated. More details will be given in Section 3.4.

A few words about the Glimm scheme for (1)
In order to use a Glimm scheme, we need to know the Riemann solutions of the problem. This computation has already been done in [2] and in [7, Section 6] where all the details can be found. We refer the reader to some classical books [34,35] for general discussions about the role of Riemann problems in the theory of conservation laws and to [2,7] for the specific case of the traffic flow system. We recap here only the Riemann solutions, omitting the details on the elementary waves and on the admissibility of solutions. (1) Let us just recall that the second eigenvalue λ 2 given by (6) of system (1) is always linearly degenerate, leading to contact discontinuities and that λ 1 is genuinely non-linear, leading to shocks or rarefaction waves. One of the difficulties of the computations is that vacuum regions may appear. Therefore, the Riemann solution of system (1), with an initial datum

A brief overview on the Riemann problem for
can be computed according to the five following cases: • if ρ L > 0, ρ R > 0 and 0 ≤ v R ≤ v L , the solution consists of a 1-shock that connects ρ L , v L to the intermediate state (ρ * , v * ) and a contact discontinuity between (ρ * , v * ) and ρ R , v R ; , the solution consists of a 1-rarefaction wave that connects ρ L , v L to (ρ * , v * ) and a contact discontinuity between (ρ * , v * ) and ρ R , v R ; • if ρ L > 0, ρ R > 0 and v L + p ρ L < v R , a vacuum region appears; the solution consists of a 1-rarefaction wave that connects ρ L , v L to (0, v * ), then a vacuum region between (0, v * ) and (0, v R ) and a contact discontinuity between 0, v R and ρ R , v R ; • if ρ L > 0, ρ R = 0, the solution is only a 1-rarefaction wave connecting ρ L , v L and 0, v R ; • if ρ L = 0, ρ R > 0, the solution is only a 2-contact discontinuity con- The intermediate state (ρ * , v * ) of the Riemann solution is computed with the Riemann invariants, that is to say In the case of a shock, the speed of the shock between (ρ * , v * ) and ρ L , v L is given by In the case of a rarefaction wave, the self-similar solution (ρ, v)(ξ) with ξ = x/t is given by the following formulae which apply for ξ ∈ λ 1 ρ L , v L , λ 1 ρ R , v R .

Remark 3
In practice, we compute the self-similar solutions of equation (13), by using the Newton algorithm. The method requires that p has the C 2 regularity. This remark explains the construction of the explicit part of the velocity offset (11) (we have observed bad behaviors of the scheme when p has jumps).

Glimm scheme
Hence, we have at hand formula to compute the solution of the Riemann problems, which are the elementary brick of the Glimm's scheme (like for Godunov's scheme). This scheme has been introduced for theoretical purposes [21], and its implementation for hyperbolic systems is further discussed in [12,14,27]. Let U n j = ρ n j , y n j be the approximated mean value on the cell C j of U = (ρ, y) at time t n . We proceed as follows to update the numerical unknown: • We solve the associated Riemann problem at each interface x j+1/2 , namely all the Riemann problems with U L = U n j , U R = U n j+1 .
• Let a n be a number picked randomly in [0, 1]. We define the value U n+1 j of the numerical unknown in the cell C j at time t n+1 as to be the solution of the Riemann problem evaluated at the point x j− 1 2 +a n ∆x ∈ C j . The scheme does not use any averaging or projection procedure and the obtained solution, by construction, remains in the invariant region of the PDE system. In practice, we use the Van Der Corput quasirandom sequence (a n ) n∈N (see [14]) defined by a n =

Treatment of the implicit part
Let us now discuss how we handle the system (12) where we remind the reader that p imp contains the stiff part of the velocity offset. As said above, the system decouples and it has a very simple structure. Let us set Φ : . We solve the scalar conservation law for ρ with the classical Engquist-Osher scheme. The Engquist-Osher flux is defined by where [X] + = max(X, 0) and [X] − = min(X, 0). The implicit scheme takes the following form where ρ n+1/2 j is the result from the first (explicit) step of the scheme.
Here, the velocity field Φ (ρ) = − d dρ ρp imp (ρ) = −p imp (ρ) − ρp imp (ρ) of this scalar equation is always non positive. Accordingly, the numerical flux reduces to a mere function of the right density Consequently, the non-linear equation (14) becomes It forms a triangular system of non linear scalar equations. If J stands for the number of grid points, we have to solve J non linear scalar equations, which means J executions of the scalar Newton algorithm to update the density. Then a mere linear system defines the updated y. Indeed, once ρ n+1 is determined, we solve the transport equation for updating y. To this end, we use the standard implicit upwind scheme The specific case p imp ≥ 0, yields It forms a triangular linear system of equations that can be solved by backward substitution, leading to the straightforward formula

Simulations
As indicated in the introduction, it is far from clear how to design a "natural" scheme for a direct simulation of the constrained model (4). We can only mention the recent approach for crowd dynamics proposed in [28]; here we are rather motivated by the asymptotic issues. Our aim is two-fold. On the one hand we wish to discuss the asymptotic behavior of the different models (VO1), (VO2) and (VO3) for the velocity offset, which are all expected to capture asymptotically (for ε → 0 or γ → ∞) the features of the limit system (4). On the other hand, we shall discuss the numerical difficulties and the ability of the time-splitting strategy, which will be compared to the standard Glimm scheme with a small enough time step, in handling the asymptotic behaviour.
The simulations presented below are thought of as Riemann problems and we impose boundary conditions that maintain constant the inflow conditions.
Of course, the method can be adapted to treat further boundary conditions.
In particular imposing zero-influx produces vacuum regions, a numerical difficulty that our method is able to handle, as shown with the decongestion case below.

Case of simple transport
To begin with, we test the case of a simple transport: the computational domain is the interval [0, 1] and for the initial data we set see Figure 1 (cyan curves). We compare the six following situations: • system (1) with pressure (VO1), using the Glimm scheme, • system (1) with pressure (VO1), using the scheme presented in Section 3, • system (1) with pressure (VO2), using the Glimm scheme, • system (1) with pressure (VO2), using the scheme presented in Section 3, • system (1) with pressure (VO3), using the Glimm scheme, • system (1) with pressure (VO3), using the scheme presented in Sec- The space step is equal to ∆x = 10 −3 and the time step is computed in order to satisfy the stability condition (evaluated with the full p for the Glimm scheme, and with p exp for the implicit-explicit method). The parameters are taken as follows: • Pressure (VO1) with γ = 2, ε = 10 −3 , ρ = 1. For the explicit-implicit scheme, the numerical threshold is chosen as ρ num = ρ 1 − 1 5 ε 1/(γ+1) .
For the explicit-implicit scheme, the numerical threshold is chosen as

Case of decongestion
Next, we study the case of a decongestion in the traffic. The data are defined The initial density is close to the threshold. Since the vehicles ahead are  In the former case, the limit behavior is captured with γ = 100 and for the latter, we get a satisfactory result with ε = 10 −5 if γ = 2, and ε = 10 −7 if γ = 3. We notice also at Figure 4 that if we take γ = 3 for pressure (VO2) instead of γ = 2, we need to take a smaller value of ε, namely ε = 10 −7 instead of ε = 10 −5 . Note that the results for the original model (VO1) and for the modified model (VO2) are totally equivalent.   The simulations are performed with the implicit-explicit scheme and the initial condition is plotted in cyan.

A jump of velocity creating a congestion
Finally, we turn to the simulation of a congestion in the traffic. The initial conditions are given by The density is initially close to the threshold; since the cars ahead are slower, a congestion might occur and the Lagrange multiplier becomes active to prevent an excess of vehicles density.
Indeed, discontinuous solutions are characterized by the Rankine-Hugoniot conditions: with t → s(t) the speed of the discontinuity curve, we havė We can check that and a Lagrange multiplier active only in the congestion domain is solution of (4). The presence of slow vehicles ahead of the fast ones instantaneously creates a congestion behind the velocity jump: the slow vehicles ahead make the faster ones behind brake. This is typical of the Follow-the-Leader approach, which has led to a derivation of the Aw-Rascle-Zhang system [1,20]. However, it is likely that solutions of the constrained model (4) are not uniquely defined for such data; we refer the reader to [6] for such considerations.
The parameters are defined as in Section 4.1 and we show the solutions obtained at the final time T = 0.01 in Figure 5. We observe exactly the same behaviors between the Glimm scheme and the implicit-explicit scheme with (VO1) or (VO2). We observe that with these parameters, the three models do not find the solution (18). The time steps for the Glimm scheme are smaller than with the explicit-implicit scheme, but, quite surprisingly, by a factor 3 or 4 only.
We observe that these parameters provide a result closer to the explicit solution (ρ 1 , v 1 ). The results for (VO3) with different values of γ are given at Figure 7 and for (VO2) with different values of γ and ε at Figure 8.
The two schemes behave equivalently in that case, with some advantage in terms of time step for the implicit-explicit method, see Table 1 (the smaller ε, resp. the larger γ, the more important the gain). Note that in the case when γ = 3, the results are exactly the same because the density is below the numerical threshold; consequently, the implicit part of the scheme is not used and the final result is the same as the one given by the explicit step, that is to say the final result is the same as the one given by the Glimm scheme.
As expected, making ε → 0 for (VO2) or γ → +∞ for (VO3) allows us to obtain a result compatible with (ρ 1 , v 1 ). In particular, we point out that the approach of (4) as an asymptotic model from (3)

A fast cluster reaching a slower one
We now consider the situation of a slow car cluster reached by a faster cluster of vehicles. The initial conditions are 0 otherwise. and the expected solution at T = 0.3 is      When using the Glimm scheme, we obtain the expected limit solution as we make the parameters vary: for (VO1) as ε goes to 0, see Figure 9, for (VO2) as ε goes to 0, see Figure 10, Here again, the quantity of mass is not conserved: the difference between the initial density and the final one is around 23%. We point out that this simulation is quite tough; the results obtained with the Glimm scheme are neater but at the price of a significative numerical cost.

Riemann problems: comparison with [16]
In this section we present the numerical results obtained with the Glimm scheme and the implicit-explicit schemes of Section 3 for the sub-cases AI and AIII of [7,16]. So we consider in the following a road [0, 1] with the initial data    and where we use the same labels as in [7,16]. A solution for (AI) is and is a solution for (AIII) (note that v does not make sense in the vacuum region x ∈ [0.5 + 0.1t, 0.5 + 0.5t]). In order to compare our scheme with the results of [16] we take the same parameters, namely for the velocity offset   According to Figure 13, all the different numerical approaches give similar results which coincide with the solution given by (22).

Conclusion
The model (4) is intended to describe the formation and the dynamics of traffic jams, through a Lagrange multiplier that accounts for a density threshold.
This model can be motivated, at least formally, through asymptotic arguments from the Aw-Rascle-Zhang system with a rescaled velocity-offset. It raises the question of simulating efficiently the Aw-Rascle-Zhang system with potentially stiff velocity offsets. Depending on the values of the parameters it can be seen either as the simulation of a model for traffic flows with stiff parameters or as a way to access the limiting behavior described by (4), alternative for instance to the approach of [28]. However the scaling induces fast propagation waves and, in turn, severe stability conditions. In this paper, we propose several approaches to obtain asymptotically (4) and we introduce an implicit-explicit method in order to cope with the large characteristic speeds of the system.
This study exhibits numerical difficulties, related to both the lack of convexity of the invariant domains of (1) and the large characteristic speeds.
We have proposed a time-splitting method, based on a decomposition of the