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

HAL
open science

# A relaxation scheme for the equations of isentropic gas dynamics on a network with jump transmission conditions

Maya Briani, Roberto Natalini, Magali Ribot

► To cite this version:

Maya Briani, Roberto Natalini, Magali Ribot. A relaxation scheme for the equations of isentropic gas dynamics on a network with jump transmission conditions. 2025. hal-05134865

HAL Id: hal-05134865

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

Preprint submitted on 7 Jul 2025

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.
A RELAXATION SCHEME FOR THE EQUATIONS OF ISENTROPIC GAS DYNAMICS ON A NETWORK WITH JUMP TRANSMISSION CONDITIONS

MAYA BRIANI, ROBERTO NATALINI, AND MAGALI RIBOT

ABSTRACT. In this paper we propose a new numerical scheme of relaxation type to approximate the Euler equations of isentropic gas dynamics on the arcs of a network. At the junction mass conservation and a jump transmission condition on the density are given, and a new solver is introduced to deal with both subsonic and supersonic cases. Consistency properties of the solver are proven and numerical tests are displayed to show its good performance also with respect to other possible solvers.

1. INTRODUCTION

In this article, we consider the Euler equations of isentropic gas dynamics set on the one dimensional arcs of a network. On these arcs the Euler equations are, as usual, given by:

$$ (1) \qquad \begin{cases} \partial_t \rho + \partial_x(\rho u) = 0, \\ \partial_t(\rho u) + \partial_x(\rho u^2 + p(\rho)) = 0, \end{cases} $$

where $x$ is a spatial coordinate, $t > 0$ is the time, $\rho$ is the density of a fluid, so we will assume all the time $\rho > 0$, and $q = \rho u$ its momentum. Also, we shall assume that the pressure $p$ is given by the pressure law for isentropic gases, that is to say

$$ (2) \qquad p(\rho) = p_0 \rho^{\gamma} \quad \text{with } \gamma > 1. $$

Usually, for isentropic fluids, the value of $\gamma$ is taken such that $1 < \gamma < 5/3$. However, here, we have in mind various applications, and in particular to biology, see for instance [26, 50, 46], where the value of $\gamma$ is not clearly defined, so in the following we consider all the values $\gamma > 1$.

While on the external boundary of the network we can assign, for instance, standard no-flux conditions to the gas, at the junction we have to impose some physically motivated transmission conditions. In this paper we address a special class of transmission conditions that, on a 2-arc network – i.e.: a network composed by just two arcs which meet at just one junction point, see Figure1 – at the junction point $x$ reads as:

$$ (3) \qquad (\rho u)_{\ell}(t, x_{\ell}) = (\rho u)_{r}(t, x_{r}) = \kappa(\rho_{\ell}(t, x_{\ell}) - \rho_{r}(t, x_{r})), \quad \text{for all time } t > 0, $$

where the subscripts $\ell$ and $r$ stand respectively for the left and right side of the junction point and $\kappa$ is a given positive constant. In the following we shall refer to these conditions as the Jump Transmission Condition (JTC).

In fact, (JTC) implies that the density flux incoming at the junction is equal to the one outgoing and so the total mass is conserved at the junction. This is similar to the Kirchoff’s law for an electrical circuit. This feature is common to most of the models of flows on networks. However almost all of these models impose continuity condition for the density at the junction, as assumed for instance in [44, 19, 23, 53], both for parabolic and hyperbolic models, see also
[27, 9] for some informed reviews in the hyperbolic case. This is the point where our condition (JTC), given by (3), differs, since it allows jump discontinuities at the junction.

Transmission conditions as those appearing in (3) were introduced, only for diffusive problems, in the simple case of two regions by Kedem and Katchalsky in [41], as permeability conditions in the description of passive transport through a biological membrane, see also [48, 16, 17] for various mathematical and related formulations. More recently they were generalized to hyperbolic problems and proposed as transmission conditions for hyperbolic-parabolic and hyperbolic-elliptic systems describing biological situations such as the movement of microorganisms on networks driven by chemotaxis [32, 33, 34, 35, 36], see also [12, 13, 10] for other models and the numerical treatment of these conditions in related situations.

In general, these density jump conditions like (JTC) seem to be more appropriate than continuity conditions when dealing with mathematical models for movements of individuals, as in traffic flows [28, 11, 29] or in biological phenomena involving bacteria, eukariotic cells, or other microorganism movements [12, 32], where, even if the flux is continuous and the mass is conserved, discontinuities at the inner nodes for the density functions are expected, as a sort of a mathematical counterpart of bottleneck phenomena.

The purpose of this article is to propose a numerical scheme for this system, based on the vector BGK approach [4, 8], which guarantees that condition (3) is fulfilled at each junction of the network without solving the Riemann problem. Actually, we propose a way to solve the (JTC) relation with a scheme which also satisfies a momentum equality junction condition, namely the fact that at the discrete level the equality of the momenta $q = \rho u$ will be enforced both as the unknown of the second equation of system (1) and as the flux of the first equation of system (1).

This approach, as long as it is possible, will be compared with a more classical solver based on the Riemann invariants at the junction, which in the following will be referred as "Riemann invariant junction condition". Unfortunately, as we will discuss in this paper, it is not clear how to deal with this approach in the general "non-subsonic" case, so the comparison will be restricted to the situations where both approaches work.

The article is organized as follows: in Section 2 we study the mathematical formulation of equations (1) on a 2-arc network, with the conditions (3) at the junction. In Section 3, we study how to solve the Riemann problem for (1) on a 2-arc network, with the conditions (3) at the junction using the classical approach with the Riemann invariants. We are able to find a solution in the subsonic case, and for some special initial conditions in the supersonic case.

Then, in Section 4, we introduce our relaxation scheme for the same problem, which can be written for all initial data. In Section 5, we present some numerical tests, and in particular a comparison between the relaxation scheme and the cases where we have an explicit solution for the Riemann problem at the junction. Finally we extend the previous study to the case of a general network in Section 6.

## 2. ISENTROPIC GAS DYNAMICS ON A 2-ARC NETWORK

Consider the 1D system (1) on a simple network, composed of two arcs, as displayed in Figure 1. The quantities linked to the arc on the left, called arc $\ell$ (resp. on the right, called arc $r$) are denoted with subscript $\ell$ (resp. $r$).
<figure><img src="image_2.png" /><figcaption>FIGURE 1. The 2-arc network : the quantities linked to the arc $[0, L_\ell]$ on the left (resp. $[0, L_r]$ on the right) are denoted with subscript $\ell$ (resp. $r$)</figcaption></figure>

We are therefore dealing with the two following systems :

$$ (4) \qquad \left\{ \begin{array}{l} \partial_t \rho_\ell(t, x) + \partial_x(\rho_\ell u_\ell)(t, x) = 0, \\ \partial_t(\rho_\ell u_\ell)(t, x) + \partial_x(\rho_\ell u_\ell^2 + p(\rho_\ell))(t, x) = 0, \quad \text{for } x \in [0, L_\ell], \end{array} \right. $$

and

$$ (5) \qquad \left\{ \begin{array}{l} \partial_t \rho_r(t, x) + \partial_x(\rho_r u_r)(t, x) = 0, \\ \partial_t(\rho_r u_r)(t, x) + \partial_x(\rho_r u_r^2 + p(\rho_r))(t, x) = 0, \quad \text{for } x \in [0, L_r], \end{array} \right. $$

complemented with no-flux boundary conditions, at least for smooth solutions, at the external points of the two intervals, namely, always setting for the flux $q = \rho u$ in the following :

$$ (6) \qquad \left\{ \begin{array}{l} \partial_x \rho_\ell(t, 0) = 0, \quad q_\ell(t, 0) = 0, \quad \text{for all } t \ge 0 \text{ on the left arc,} \\ \partial_x \rho_r(t, L_r) = 0, \quad q_r(t, L_r) = 0, \quad \text{for all } t \ge 0 \text{ on the right arc.} \end{array} \right. $$

Clearly the boundary conditions, as usual for hyperbolic problems, have to be intended in a proper sense, since sometimes, due to the direction of the characteristics with respect to the boundary, they need to be interpreted, see [52] and references therein for more details.

Finally, we define transmission conditions at the junction, which couple the unknowns of both systems. For that purpose, we consider the equality of fluxes in order to have the mass conservation at the junction, i.e. :

$$ (7) \qquad q_\ell(t, L_\ell) = q_r(t, 0), \quad \text{for all } t \ge 0. $$

We supplement this condition with a jump transmission condition (JTC) for the density at the junction, that is to say :

$$ (8) \qquad q_\ell(t, L_\ell) = q_r(t, 0) = \kappa(\rho_\ell(t, L_\ell) - \rho_r(t, 0)), \quad \text{for all } t \ge 0, $$

where $\kappa$ is a permeability coefficient that we suppose to be constant. These previous conditions represent some transport conditions through a jump (given by a membrane or a change of domain), and the fluxes are linked with the difference of densities on both sides of the jump.

This condition (JTC) is the hyperbolic analogous of the parabolic Kedem-Katchalsky transmission condition, see for instance [41, 48, 16, 17]. The connection between condition (8) and the Kedem-Katchalsky condition can be seen in the diffusive limit of the damped version of equation (1). By scaling this equation as in [25], we obtain the perturbed system

$$ \left\{ \begin{array}{l} \partial_t \rho + \partial_x(\rho u) = 0, \\ \partial_t(\epsilon^2 \rho u) + \partial_x(\epsilon^2 \rho u^2 + p(\rho)) = -\rho u, \end{array} \right. $$

where the inertial terms $\partial_t(\rho u) + \partial_x(\rho u^2)$ are supposed to be negligible.

So, as $\epsilon$ tends to zero, the momentum $\rho u$ approaches $-\partial_x p(\rho)$ and, as shown in [25], we obtain the porous medium equation

$$ \partial_t \rho - \partial_{xx} p(\rho) = 0, $$
see also [47] in the framework of asymptotic preserving numerical schemes for isentropic Euler equations coupled with chemotaxis.

Replacing $q = \rho u$ by $-\partial_x p(\rho)$ in (8), we obtain the condition

$$-\partial_x p(\rho_\ell)(t, L_\ell) = -\partial_x p(\rho_r)(t, 0) = \kappa(\rho_\ell(t, L_\ell) - \rho_r(t, 0))$$

which is in the form of the diffusive Kedem-Katchalsky condition treated for instance in [30, 18].

Clearly, condition (7) together with condition (6) ensure the mass conservation of the system, that is to say, working with smooth solutions :

$$(9) \qquad \frac{d}{dt}m(t) = 0, \text{ with } m(t) = \int_{0}^{L_{\ell}} \rho_{\ell}(t, x) \, dx + \int_{0}^{L_{r}} \rho_{r}(t, x) \, dx.$$

## 2.1. Entropy dissipation.
Let us prove now that, at least for smooth solutions, system (4)-(8) is endowed with an entropy dissipation property at the continuous level under the condition that the density jump is small enough at the junction.

First, on an interval, there exists at least an entropy-entropy flux pair $(\eta, G)$ for system (1), see [24] p.212, namely the functions

$$(10) \qquad \eta(\rho, \rho u) = \frac{1}{2} \rho u^2 + \frac{p_0}{\gamma - 1} \rho^\gamma \quad \text{and} \quad G(\rho, \rho u) = \left( \frac{1}{2} \rho u^2 + \frac{p_0 \gamma}{\gamma - 1} \rho^\gamma \right) u.$$

The function $\eta$ is convex and, for smooth solutions $(\rho, q = \rho u)$ to system (1) we have

$$\partial_t \eta(\rho, q) + \partial_x G(\rho, q) = 0.$$

For weak solutions, the admissible (entropy) solutions have to verify (in the weak sense) the inequality

$$(11) \qquad \partial_t \eta(\rho, q) + \partial_x G(\rho, q) \le 0,$$

and we can say that the admissible solutions dissipate the entropy (see again [24]). In general, for one single interval $[0, L]$ with no flux conditions, the quantity

$$\int_0^L \eta(\rho(t, x), q(t, x)) dx$$

is decreasing in time.

Now, we want to investigate if condition (JTC) in (8) preserves this dissipation property, at least in some range of solutions. Let us consider smooth solutions to equation (11) on the 2 arcs of the previous network and integrate it with respect to time and space :

$$\begin{align*}
S(t) &:= \int_{0}^{L_{\ell}} \eta(\rho_{\ell}(t, x), q_{\ell}(t, x)) dx + \int_{0}^{L_{r}} \eta(\rho_{r}(t, x), q_{r}(t, x)) dx \\
&\le \int_{0}^{L_{\ell}} \eta(\rho_{\ell}(0, x), q_{\ell}(0, x)) dx + \int_{0}^{L_{r}} \eta(\rho_{r}(0, x), q_{r}(0, x)) dx \\
&\quad - \int_{0}^{t} (G(\rho_{\ell}(s, L_{\ell}), q_{\ell}(s, L_{\ell})) - G(\rho_{\ell}(s, 0), q_{\ell}(s, 0))) ds \\
&\quad - \int_{0}^{t} (G(\rho_{r}(s, L_{r}), q_{r}(s, L_{r})) - G(\rho_{r}(s, 0), q_{r}(s, 0))) ds.
\end{align*}$$

Thanks to boundary conditions on the external domain, we find that since $\gamma > 1$, we have

$$G(\rho_r(s, L_r), q_r(s, L_r)) = G(\rho_\ell(s, 0), q_\ell(s, 0)) = 0,$$
which leads to the inequality :

$$
\begin{equation} \tag{12}
\begin{split}
S(t) = \int_0^{L_\ell} & \eta(\rho_\ell(t, x), q_\ell(t, x)) dx + \int_0^{L_r} \eta(\rho_r(t, x), q_r(t, x)) dx \\
& \le \int_0^{L_\ell} \eta(\rho_\ell(0, x), q_\ell(0, x)) dx + \int_0^{L_r} \eta(\rho_r(0, x), q_r(0, x)) dx \\
& \phantom{\le} - \int_0^t (G(\rho_\ell(s, L_\ell), q_\ell(s, L_\ell)) - G(\rho_r(s, 0), q_r(s, 0))) ds.
\end{split}
\end{equation}
$$

Therefore, if the quantity $\Delta G(t) := G(\rho_{\ell}(t, L_{\ell}), q_{\ell}(t, L_{\ell})) - G(\rho_{r}(t, 0), q_{r}(t, 0)) \ge 0$ for all $t > 0$, we have entropy dissipation, that is to say for all $t > 0$ it holds

$$
S(t) \leq S(0)
$$

for all $t > 0$.

Now, using the conditions at the junction, we can compute $\Delta G$ as follows. Using that $\rho > 0$, we rewrite $G(\rho, q)$ as :

$$
G(\rho, q) = G(\rho, \rho u) = \frac{1}{2} \rho u^3 + \frac{p_0 \gamma}{\gamma - 1} \rho^{\gamma} u = \frac{1}{2} \frac{q^3}{\rho^2} + \frac{p_0 \gamma}{\gamma - 1} \rho^{\gamma-1} q
$$

and therefore

$$
\begin{align*}
\Delta G(t) &= G(\rho_{\ell}(t, L_{\ell}), q_{\ell}(t, L_{\ell})) - G(\rho_{r}(t, 0), q_{r}(t, 0)) \\
&= \frac{1}{2} \left( \frac{(q_{\ell}(t, L_{\ell}))^3}{(\rho_{\ell}(t, L_{\ell}))^2} - \frac{(q_{r}(t, 0))^3}{(\rho_{r}(t, 0))^2} \right) + \frac{p_0 \gamma}{\gamma - 1} \left( (\rho_{\ell}(t, L_{\ell}))^{\gamma-1} q_{\ell}(t, L_{\ell}) - (\rho_{r}(t, 0))^{\gamma-1} q_{r}(t, 0) \right)
\end{align*}
$$

Using the conditions (8) at the junction, the previous expression can be simplified as :

$$
\begin{align*}
\Delta G(t) &= G(\rho_{\ell}(t, L_{\ell}), q_{\ell}(t, L_{\ell})) - G(\rho_{r}(t, 0), q_{r}(t, 0)) \\
&= -\frac{1}{2}\kappa^3(\rho_{\ell}(t, L_{\ell}) - \rho_{r}(t, 0))^4 \frac{\rho_{r}(t, 0) + \rho_{\ell}(t, L_{\ell})}{(\rho_{r}(t, 0))^2(\rho_{\ell}(t, L_{\ell}))^2} \\
&\quad + \frac{p_0\gamma}{\gamma-1}\kappa(\rho_{\ell}(t, L_{\ell}) - \rho_{r}(t, 0))((\rho_{\ell}(t, L_{\ell}))^{\gamma-1} - (\rho_{r}(t, 0))^{\gamma-1}).
\end{align*}
$$

The first term on the right hand-side is negative whereas the second one is positive. We remark that as long as the difference $\rho_{\ell}(t, L_{\ell}) - \rho_{r}(t, 0)$ is small,

$$
\Delta G(t) \sim \frac{p_0 \gamma}{\gamma - 1} \kappa (\rho_{\ell}(t, L_{\ell}) - \rho_{r}(t, 0)) ((\rho_{\ell}(t, L_{\ell}))^{\gamma-1} - (\rho_{r}(t, 0))^{\gamma-1}) \geq 0
$$

and the dissipation entropy property holds, whereas as the difference $\rho_{\ell}(t, L_{\ell}) - \rho_{r}(t, 0)$ becomes large, we have

$$
\Delta G(t) \sim -\frac{1}{2} \kappa^3 (\rho_{\ell}(t, L_{\ell}) - \rho_{r}(t, 0))^4 \frac{\rho_{r}(t, 0) + \rho_{\ell}(t, L_{\ell})}{(\rho_{r}(t, 0))^2 (\rho_{\ell}(t, L_{\ell}))^2} \le 0,
$$

and then, a priori, the dissipation property is no longer true. In Section 5 we will test this dissipation property more closely, for various initial data, since in some cases the dissipation inside the arcs can compensate the growth of the entropy at the junction.

3. THE RIEMANN PROBLEM AT THE JUNCTION - THE CASE OF A 2-ARC NETWORK

We want to solve the Riemann problem for (4)–(8) at the junction using the classical approach with Riemann invariants.
### 3.1. The Riemann problem on a single interval.
Let us first describe Riemann invariants on an interval before generalizing it to the junction case, see for example [24] for more details.
Let consider isentropic gas dynamics given by Eq.(1) on a single arc, under the form

$$ (13) \qquad \partial_t U + \partial_x F(U) = 0, $$

where $U \in \mathbb{R}^2$ is the vector of unknowns $U = \begin{pmatrix} \rho \\ q = \rho u \end{pmatrix}$, and the function $F: \mathbb{R}^2 \to \mathbb{R}^2$ is the flux of the system defined by $F(U) = \begin{pmatrix} \rho u \\ \rho u^2 + p(\rho) \end{pmatrix}$.

The eigenvalues of the system are

$$ \mu_1(\rho, q) = \frac{q}{\rho} - \sqrt{p'(\rho)}, \quad \mu_2(\rho, q) = \frac{q}{\rho} + \sqrt{p'(\rho)}, $$

where $q = \rho u$, and the eigenvectors are

$$ r_1(\rho, q) = \begin{pmatrix} \rho \\ q - \rho \sqrt{p'(\rho)} \end{pmatrix}, \quad r_2(\rho, q) = \begin{pmatrix} \rho \\ q + \rho \sqrt{p'(\rho)} \end{pmatrix}. $$

The system is strictly hyperbolic in $\mathcal{D} = \{(\rho, q) : \rho > 0 \text{ or } q = 0\}$ with $\mu_1 < \mu_2$. It is genuinely nonlinear in both characteristic fields, i.e. $\nabla \mu_i \cdot r_i \neq 0$, $i=1,2$. The two functions $z_1 = z_1(\rho, q)$, $z_2 = z_2(\rho, q)$ are called the *Riemann invariants* of system (1) corresponding to $\mu_1$ and $\mu_2$ if they satisfy the equations

$$ (14) \qquad \nabla z_i \cdot r_i = 0, \quad i = 1, 2. $$

One solution of (14), when $p$ follows the pressure law of isentropic gases (2) is

$$ (15) \qquad z_1(\rho, q) = \frac{q}{\rho} + \alpha_{\gamma} \sqrt{p'(\rho)}, \quad z_2(\rho, q) = \frac{q}{\rho} - \alpha_{\gamma} \sqrt{p'(\rho)}, \quad \text{with } \alpha_{\gamma} = \frac{2}{\gamma - 1}. $$

Consider now the Riemann problem for (13), that is to say system (13) complemented with piecewise constant initial data such that

$$ U(x, 0) = \begin{cases} U^{-} & \text{if } x < 0, \\ U^{+} & \text{if } x > 0, \end{cases} $$

for two given constant left and right states $U^{-} = (\rho^{-}, q^{-})$ and $U^{+} = (\rho^{+}, q^{+})$. The solution to the Riemann problem is then explicitly known and defined by an intermediate state $\hat{U} = (\hat{\rho}, \hat{q})$ verifying the following conditions:

i) The left state $U^{-}$ is connected to $\hat{U}$ by a 1-shock or a 1-rarefaction wave such that

$$ z_1(\rho^{-}, q^{-}) = z_1(\hat{\rho}, \hat{q}). $$

Moreover, shock waves satisfy $\mu_1(U^{-}) > \mu_1(\hat{U})$ (Lax admissibility condition) while rarefaction waves satisfy $\mu_1(U^{-}) < \mu_1(\hat{U})$. We then have a 1-shock if $\hat{\rho} > \rho^{-}$ and $z_1(\rho^{-}, q^{-}) = z_1(\hat{\rho}, \hat{q})$ or a 1-rarefaction if $\hat{\rho} < \rho^{-}$ and $z_1(\rho^{-}, q^{-}) = z_1(\hat{\rho}, \hat{q})$.

ii) The intermediate state $\hat{U}$ is connected to the right state $U^{+}$ by a 2-shock if $\hat{\rho} > \rho^{+}$ and $z_2(\hat{\rho}, \hat{q}) = z_2(\rho^{+}, q^{+})$ or by a 2-rarefaction if $\hat{\rho} < \rho^{+}$ and $z_2(\hat{\rho}, \hat{q}) = z_2(\rho^{+}, q^{+})$.

The intermediate state $\hat{U}$ is then uniquely defined if, in the $(\rho, u)$-plane, the following two curves $u_1(\rho) = u^{-} + \alpha_{\gamma}(\sqrt{p'(\rho^{-})} - \sqrt{p'(\rho)})$ and $u_2(\rho) = u^{+} - \alpha_{\gamma}(\sqrt{p'(\rho^{+})} - \sqrt{p'(\rho)})$ intersect. Since $u_1(\rho)$ is strictly decreasing, unbounded, and starts at $u^{-} + \alpha_{\gamma}\sqrt{p'(\rho^{-})}$ and
$u_2(\rho)$ is strictly increasing, unbounded, with minimum $u^+ - \alpha_\gamma \sqrt{p'(\rho^+)}$, the Riemann problem has a unique solution in the region where

$$u^{-} + \alpha_{\gamma} \sqrt{p'(\rho^{-})} \geq u^{+} - \alpha_{\gamma} \sqrt{p'(\rho^{+})}$$

or equivalently

$$(16) \qquad z_1(\rho^{-}, q^{-}) \geq z_2(\rho^{+}, q^{+}).$$

Furthermore, we say that a state $(\rho, q) \in D\{(0,0)\}$ is subsonic if $\mu_1(\rho, q) < 0$ or $\mu_2(\rho, q) < 0$ or $\mu_1(\rho, q) > 0$ or $\mu_2(\rho, q) > 0$.

**3.2. The Riemann problem at the junction: the subsonic case.** Now, we shall consider a simple junction connecting two arcs as sketched in Figure 1. We use the same notation as in Section 2 and in particular, we shall use the subscript $\ell$ for the left arc and $r$ for the right one.

**Notations 3.1.** From now on, we shall use the following notations:

* $U_{\ell}^{-} = (\rho_{\ell}^{-}, q_{\ell}^{-})$ and $U_{r}^{+} = (\rho_{r}^{+}, q_{r}^{+})$ denote the constant states on the left arc and on the right arc;

* $U_{\ell}^{*} = (\rho_{\ell}^{*}, q_{\ell}^{*})$ and $U_{r}^{*} = (\rho_{r}^{*}, q_{r}^{*})$ represent the unknown values of the left and right states at the junction.

Following what has been done previously, we define the two sonic curves that determine the boundary between subsonic and supersonic regions, see the black curves on Figure 2, thanks to the following functions

$$(17) \qquad \begin{cases} q_s^+(\rho) = \rho\sqrt{p'(\rho)}, \\ q_s^-(\rho) = -\rho\sqrt{p'(\rho)}. \end{cases}$$

We define the constant data $U_{\ell}^{-} = (\rho_{\ell}^{-}, q_{\ell}^{-})$ and $U_{r}^{+} = (\rho_{r}^{+}, q_{r}^{+})$ given on the left and right arc respectively. We will assume that we are in the subsonic case at the junction, that is to say that $q_{s}^{-}(\rho_{\ell}^{-}) < q_{\ell}^{-} < q_{s}^{+}(\rho_{\ell}^{-})$ and $q_{s}^{-}(\rho_{r}^{+}) < q_{r}^{+} < q_{s}^{+}(\rho_{r}^{+})$.

Now, in order to solve the Riemann problem, the unknown states $U_{\ell}^{*}$ and $U_{r}^{*}$ at the junction have to lie on the Riemann invariant curves defined from $U_{\ell}^{-}$ and $U_{r}^{+}$. Considering the jump condition (JTC) in (8), and assuming that the states are all subsonic, we then get the following system of equations

$$(18) \qquad \begin{cases} z_1(\rho_{\ell}^*, q_{\ell}^*) = z_1(\rho_{\ell}^{-}, q_{\ell}^{-}), \\ z_2(\rho_r^*, q_r^*) = z_2(\rho_r^+, q_r^+), \\ q_{\ell}^* = q_r^* = \kappa(\rho_{\ell}^* - \rho_r^*), \end{cases}$$

subject to the constraints

$$(19) \qquad \begin{cases} q_s^-(\rho_{\ell}^*) < q_{\ell}^* < q_s^+(\rho_{\ell}^*), \\ q_s^-(\rho_r^*) < q_r^* < q_s^+(\rho_r^*). \end{cases}$$

Following Eq. (15), we define the Riemann invariant curves

$$(20) \qquad \begin{cases} q_{\ell}^{*}(\rho) = z_{1}(\rho_{\ell}^{-}, q_{\ell}^{-})\rho - \alpha_{\gamma}\rho\sqrt{p'(\rho)}, \\ q_{r}^{*}(\rho) = z_{2}(\rho_{r}^{+}, q_{r}^{+})\rho + \alpha_{\gamma}\rho\sqrt{p'(\rho)}. \end{cases}$$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Riemann problem in the subsonic case with $p(\rho) = p_0\rho^2$ (i.e. $\gamma = 2$): sonic curves $q_s^-$ and $q_s^+$ defined at Eq.(17) delimiting the subsonic region and Riemann invariant curves $q_\ell^*$ and $q_r^*$ defined at Eq.(20)</figcaption></figure>

We notice that in the case where $1 < \gamma \le 3$ and in the subsonic case, then $\mu_1(\rho_\ell^-, q_\ell^-) < 0 < \mu_2(\rho_r^+, q_r^+)$, implies $z_1(\rho_\ell^-, q_\ell^-) > 0 > z_2(\rho_r^+, q_r^+)$. Moreover, $q_\ell^*$ is a strictly concave function of $\rho$ and $q_r^*$ is strictly convex. Thus condition (16) holds and there exists a unique intersection point, apart from 0, for the two curves of functions $q_\ell^*$ and $q_r^*$. From now on, we shall call this point $U_I = (\rho_I, q_I)$ and it satisfies $q_I = q_\ell^*(\rho_I) = q_r^*(\rho_I)$.

Now, in the case of the conditions (18) we consider here, we may allow a discontinuity on densities, i.e. $\rho_\ell^* \neq \rho_r^*$. In the next lemma, we shall show that system (18), subject to (19), admits a unique solution, which lies therefore in the subsonic region.

**Lemma 3.1.** We assume $p$ of the form (2) with $1 < \gamma \le 3$. Let be given two points $U_\ell^- = (\rho_\ell^-, q_\ell^-)$ and $U_r^+ = (\rho_r^+, q_r^+)$ such that $q_s^-(\rho_\ell^-) < q_\ell^- < q_s^+(\rho_\ell^-)$ and $q_s^-(\rho_r^+) < q_r^+ < q_s^+(\rho_r^+)$, then there exists a unique couple of points $U_\ell^*, U_r^*$ verifying system (18) subject to constraints (19).

*Proof.* First of all, note that the solution $q^* = q_\ell^*(\rho_\ell^*) = q_r^*(\rho_r^*)$ to system (18) should be sought, see Figure2, in the region defined by

$$ \{q : q_{min}^* = \min_{\rho>0} q_r^*(\rho) \le q \le \max_{\rho>0} q_\ell^*(\rho) = q_{max}^*\}. $$

Note that the intersection point $(\rho_I, q_I)$ satisfies $q_{min}^* < q_I < q_{max}^*$.

Then, it is easy to check that the intersection point, apart from 0, between the sonic curve defined by $q_s^+$ and the curve defined by $q_\ell^*$ coincides with the maximum point for $q_\ell^*$. Indeed,

$$ q_s^+(\rho) = q_\ell^*(\rho) \text{ and } \rho > 0 \iff \sqrt{p'(\rho)} = \frac{\gamma - 1}{\gamma + 1} z_1(\rho_\ell^-, q_\ell^-) $$

and

$$ (q_\ell^*)'(\rho) = 0 \iff \sqrt{p'(\rho)} = \frac{\gamma - 1}{\gamma + 1} z_1(\rho_\ell^-, q_\ell^-). $$
Analogously, the minimum point for $q_r^*$ coincides with the intersection point, apart from 0, between the sonic curve defined by $q_s^-$ and the curve defined by $q_r^*$.

Therefore, in the subsonic region, $q_\ell^*$ (resp. $q_r^*$) is a continuous and strictly decreasing (resp. strictly increasing) function of $\rho$ and we can therefore define its inverse function $\phi_\ell$ (resp. $\phi_r$), which is also a continuous and strictly decreasing (resp. strictly increasing) function of $q$.

We can now prove that for all $\kappa \in [0, +\infty)$ there exists a unique solution to system (18) subject to subsonic conditions (19).

Such a solution satisfies $\rho_\ell^* = \phi_\ell(q^*)$, $\rho_r^* = \phi_r(q^*)$ and $q^* = \kappa(\rho_\ell^* - \rho_r^*) = \kappa(\phi_\ell(q^*) - \phi_r(q^*))$ and the existence and unicity of the solution follows from the fact that the function $K$ defined by

$$K(q) = \frac{q}{\phi_\ell(q) - \phi_r(q)}$$

is continuous and strictly monotone, i.e. we have that

(i) if $q_I > 0$ then $K: [0, q_I] \to [0, +\infty)$ is continuous and strictly increasing;

(ii) if $q_I < 0$ then $K: [q_I, 0] \to [0, +\infty)$ is continuous and strictly decreasing.

□

Staying in the subsonic region is equivalent to limiting the amplitude of the jump, specifically for $\kappa \to 0$ then $q^* \to 0$ which belongs to the subsonic region and for $\kappa \to +\infty$ then $q^* \to q_I$ and the jump amplitude goes to zero, i.e. $\rho_\ell^* - \rho_r^* \to 0$.

**Remark 3.1.** When $\gamma > 3$, if we impose condition (16) for the left and right states, the intersection point $q_I$ exists and system (18) admits a solution, but we cannot guarantee that this solution is subsonic.

### 3.3. Riemann invariants and the junction - the supersonic case.
In the supersonic case, the Riemann solver leading to system (18) does not apply. A different formulation is needed, which in general is still an open problem. A number of possible approaches have been proposed in the literature, see for instance [39], but here we follow the approach proposed in [14]. The idea stated in this paper can be resumed as follows:

* Assume to be on a left arc and to have to the left of the junction a supersonic state $U_\ell^-$ with $\mu_1(\rho_\ell^-, q_\ell^-) > 0$, namely $q_\ell^- > q_s^+(\rho_\ell^-)$. Introducing the point $\tilde{U}_\ell^- = (\tilde{\rho}_\ell^-, \tilde{q}_\ell^-)$ such that $\tilde{q}_\ell^- = q_\ell^-$, $z_1(\tilde{\rho}_\ell^-, q_\ell^-) = z_1(\rho_\ell^-, q_\ell^-)$ and $\mu_1(\tilde{U}_\ell^-) < 0$, the region of possible junction values $U_\ell^*$ is defined as the set of points verifying

$$z_1(\rho_\ell^*, q_\ell^*) = z_1(\tilde{\rho}_\ell^-, q_\ell^-) \quad \text{with} \quad \rho_\ell^* > \tilde{\rho}_\ell^-.$$

* Assume to be on a right arc and to have to the right of the junction a supersonic state $U_r^+$ with $\mu_2(\rho_r^+, q_r^+) < 0$, namely $q_r^+ < q_s^-(\rho_r^+)$. Introducing the point $\tilde{U}_r^+ = (\tilde{\rho}_r^+, \tilde{q}_r^+)$ such that $\tilde{q}_r^+ = q_r^+$, $z_2(\tilde{\rho}_r^+, q_r^+) = z_2(\rho_r^+, q_r^+)$ and $\mu_2(\tilde{U}_r^+) > 0$, the region of possible junction values $U_r^*$ is defined as the set of points verifying

$$z_2(\rho_r^*, q_r^*) = z_2(\tilde{\rho}_r^+, q_r^+) \quad \text{with} \quad \rho_r^* < \tilde{\rho}_r^+.$$

To summarize, on the one hand, if we are in the supersonic region with $\mu_1(\rho_\ell^-, q_\ell^-) > 0$ (resp. $\mu_2(\rho_r^+, q_r^+) < 0$) the idea is to connect first the left (resp. right) state $U_\ell^-$ (resp. $U_r^+$) with a shock with zero speed to the subsonic state $\tilde{U}_\ell^-$ (resp. $\tilde{U}_r^+$), i.e. $\frac{q_\ell^- - \tilde{q}_\ell^-}{\rho_\ell^- - \tilde{\rho}_\ell^-} = 0$ (resp.
$\frac{q_r^+ - \tilde{q}_r^+}{\rho_r^+ - \tilde{\rho}_r^+} = 0$), and then the latter to $U_\ell^*$ (resp. $U_r^*$) with a shock with negative (resp. positive) speed.

On the other hand, if we are in the supersonic region with $\mu_1(\rho_\ell^-, q_\ell^-) < 0$, that is $q_\ell^-$ is negative and $q_\ell^- < q_s^- (\rho_\ell^-)$ (resp. $\mu_2(\rho_r^+, q_r^+) > 0$, namely $q_r^+$ is positive and $q_r^+ > q_s^+(\rho_r^+)$), see Figure 3, an admissible solution in the sense of the definition of Riemann solver given in [14] is to stay on the invariant curve $z_1(\rho_\ell^*, q_\ell^*) = z_1(\rho_\ell^-, q_\ell^-)$ (resp. $z_2(\rho_r^*, q_r^*) = z_2(\rho_r^+, q_r^+)$) with the constraint $\mu_2(\rho_\ell^*, q_\ell^*) < 0$ (resp. $\mu_1(\rho_r^*, q_r^*) > 0$).

Moreover, in both cases, the values $U_\ell^*$ and $U_r^*$ have to verify the jump condition (8). It is clear that this procedure does not always lead to a solution. However, as stated in [14], the region of admissible junction solutions may be enlarged but losing the uniqueness, and additional conditions have to be added to define the solution at a junction. This discussion is beyond the scope of our work, and we will simply compare simulations in cases where the above procedure is well defined.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Riemann problem in the supersonic case with $p(\rho) = p_0\rho^2$: on the left, right supersonic case with $\mu_2(\rho_r^+, q_r^+) < 0$, i.e. $q_r^+ < q_s^-(\rho_r^+)$, and construction of the point $\tilde{U}_r^+ = (\tilde{\rho}_r^+, \tilde{q}_r^+)$. On the right, left supersonic case with $\mu_1(\rho_\ell^-, q_\ell^-) < 0$, i.e. $q_\ell^- < q_s^-(\rho_\ell^-)$; here, we stay on the invariant curve $q_\ell^*$.</figcaption></figure>

## 4. A RELAXATION SCHEME FOR THE JUMP JUNCTION CONDITION - CASE OF A 2-ARC NETWORK

For the sake of simplicity, we explain in this section the numerical scheme and the numerical jump junction condition we develop in the case of a simple network, composed of two arcs. An extension to more general networks will be detailed in Section 6.
4.1. **A relaxation scheme on a single interval.** Now, let us describe the numerical scheme we use to approximate (4)-(8). To begin with, we use a relaxation scheme, first introduced in [40], see also [3, 45, 4, 8], which we will see, away from the junction, is just another interpretation of the HLL scheme introduced in [37].

We rewrite our equations in the vector form already presented in Section 3, the same as equation (13):

$$ (21) \qquad \partial_t U + \partial_x F(U) = 0, $$

where $U \in \mathbb{R}^2$ is the vector of unknowns, and the function $F: \mathbb{R}^2 \to \mathbb{R}^2$ is the flux of the system.

Now, we introduce its relaxation approximation, namely a semilinear hyperbolic system with a singular perturbation parameter $\epsilon$:

$$ (22) \qquad \left\{ \begin{array}{l} \partial_t f_1^\epsilon + c_1 \partial_x f_1^\epsilon = \frac{1}{\epsilon} (M_1(U^\epsilon) - f_1^\epsilon), \\ \\ \partial_t f_2^\epsilon + c_2 \partial_x f_2^\epsilon = \frac{1}{\epsilon} (M_2(U^\epsilon) - f_2^\epsilon). \end{array} \right. $$

Here $c_1, c_2$ are real constants such that $c_1 < c_2, f_i^\epsilon \in \mathbb{R}^2$, for $i = 1, 2$, and $U^\epsilon := f_1^\epsilon + f_2^\epsilon$. The functions $M_i(U): \mathbb{R}^2 \to \mathbb{R}^2$, for $i = 1, 2$, are called Maxwellians for the system and are vector functions satisfying the two following relations:

$$ (23) \qquad \left\{ \begin{array}{l} M_1(U) + M_2(U) = U, \\ c_1 M_1(U) + c_2 M_2(U) = F(U). \end{array} \right. $$

From (23), we can obtain explicit expressions for $M_1(U)$ and $M_2(U)$ as

$$ (24) \qquad M_1(U) = \frac{c_2 U - F(U)}{c_2 - c_1} \quad \text{and} \quad M_2(U) = \frac{c_1 U - F(U)}{c_1 - c_2}. $$

Notice that, setting $V^\epsilon = c_1 f_1^\epsilon + c_2 f_2^\epsilon$, from system (22) we easily obtain

$$ \left\{ \begin{array}{l} \partial_t U^\epsilon + \partial_x V^\epsilon = 0, \\ \\ \partial_t V^\epsilon + \partial_x ((c_1 + c_2)V^\epsilon - c_1c_2U^\epsilon) = \frac{1}{\epsilon} (F(U^\epsilon) - V^\epsilon). \end{array} \right. $$

So, as $\epsilon$ goes to zero, we expect $U^\epsilon$ converges to the solution to (21), as shown for the scalar case in [45], see also [51, 7, 6] for the system case in one space dimension. Clearly, when $c_1 = -\lambda < 0 < c_2 = \lambda$, for some real constant $\lambda > 0$, we recover the classical Jin-Xin model of relaxation [40].

Notice that, as shown in [7] (see also [4]), a necessary condition to have the stability for this approximation is that the Jacobians of the Maxwellian functions $M_i(U)$, are strictly positive matrices for $i = 1, 2$ and all values of $U$ under consideration. Using the explicit expression (24), it is easy to see that the so-called subcharacteristic condition has to be verified. Namely, if $\mu(U)$ is a real eigenvalue of the Jacobian matrix $F'(U)$, we need to have that

$$ (25) \qquad c_1 \le \mu(U) \le c_2 \quad \text{for all } U \text{ under consideration.} $$

This condition will be required not only to assess the consistency of the correspondent scheme, but also its stability.
The strategy to use this relaxation approach to create a numerical scheme for equation (21), is in general to solve numerically at each time step the homogeneous part of system (22)

$$ \left\{ \begin{array}{l} \partial_t f_1 + c_1 \partial_x f_1 = 0, \\ \partial_t f_2 + c_2 \partial_x f_2 = 0, \end{array} \right. $$

and then to reconstruct the solution to (21) from a suitable projection of the functions ($f_1, f_2$) on the Maxwellian functions $M_1(U)$, $M_2(U)$. This strategy has been fully implemented in [4] in the framework of Finite Difference schemes, see also the recent developments in the framework of Lattice Boltzmann schemes in [2]. In the following, as already done in [46], we prefer to follow the Finite Volume approach of [8]. So, let us introduce some notations to present our scheme.

We begin by describing the discretization process on a single interval, which will be the numerical strategy that we will use both on arc $\ell$ and on arc $r$.

Let us denote by $\Delta x$ the spatial discretization step. We discretize the interval using $J$ cells

$$ C_j = [x_{j-1/2}, x_{j+1/2}], \quad 1 \le j \le J $$

centered at the nodes $x_j = (j - \frac{1}{2})\Delta x$, $1 \le j \le J$.

We also denote by $\Delta t$ the time step. Discrete times are denoted by $t_n = n\Delta t, n \in \mathbb{N}$. The time index will be omitted when possible for the sake of simplicity.

Consider system (1) under the form (21). We define $U_j^n$ to be the approximation of the mean value of vector $U$ on cell $C_j$ at time $t_n$, that is to say $U_j^n \sim \frac{1}{\Delta x} \int_{C_j} U(t^n, x)dx$. According to the finite volume framework, an explicit conservative finite volume scheme for discretizing equation (21) can be expressed as follows

$$ (26) \qquad U_j^{n+1} = U_j^n - \frac{\Delta t}{\Delta x} \left( \mathcal{F}_{j+1/2}^n - \mathcal{F}_{j-1/2}^n \right), $$

with $\mathcal{F}_{j+1/2}^n = \mathcal{F}(U_j^n, U_{j+1}^n)$ for an explicit three point scheme. The scheme is therefore defined by the function $\mathcal{F}$ such that $\mathcal{F}(U^{-}, U^{+}) = \begin{pmatrix} \mathcal{F}^{\rho}(U^{-}, U^{+}) \\ \mathcal{F}^{q}(U^{-}, U^{+}) \end{pmatrix} \in \mathbb{R}^2$, which is called the numerical flux.

We now describe the scheme we will use on each arc, which will be equivalent to the HLL scheme. Here we follow the presentation given in [8].

Consider again the homogeneous counterpart to the relaxation system (22), with $f_i = f_i(x, t) \in \mathbb{R}^2$ for $i = 1, 2$,

$$ (27) \qquad \left\{ \begin{array}{l} \partial_t f_1 + c_1 \partial_x f_1 = 0, \\ \partial_t f_2 + c_2 \partial_x f_2 = 0, \end{array} \right. $$

where $c_1 < c_2$ are two speeds to be chosen later on.

Noticing that system (27) can be solved explicitly, we obtain an approximate solution to system (21) as a function of the solution to system (27). Indeed, an exact Riemann solver for
(27) gives

$$ (28) \qquad \mathcal{R}(x/t, f^-, f^+) = \begin{cases} (f_1^-, f_2^-) & \text{if } x/t < c_1, \\ (f_1^+, f_2^-) & \text{if } c_1 < x/t < c_2, \\ (f_1^+, f_2^+) & \text{if } c_2 < x/t, \end{cases} $$

with $f_i^\pm \sim M_i(U^\pm)$, $i=1,2$. Thus we get the following approximate HLL numerical flux for the scheme in conservative form (26)

$$ (29) \qquad \mathcal{F}(U^-, U^+) = \begin{cases} c_1 M_1(U^-) + c_2 M_2(U^-) = F(U^-), & \text{if } 0 < c_1, \\ c_1 M_1(U^+) + c_2 M_2(U^-) = \frac{c_2 F(U^-) - c_1 F(U^+)}{c_2 - c_1} + \frac{c_2 c_1}{c_2 - c_1} (U^+ - U^-), & \text{if } c_1 < 0 < c_2, \\ c_1 M_1(U^+) + c_2 M_2(U^+) = F(U^+), & \text{if } c_2 < 0. \end{cases} $$

Note that for stability reasons and to have some entropy dissipation, $c_1$ and $c_2$ have to be chosen accordingly to the inequality (25).

From now on, to deal with the transmission condition (JTC), we always stay in the case $c_1 < 0 < c_2$. For simplicity reasons, we set $c_1 = -c_2 = \lambda$ with $\lambda > \max_U |\mu(U)| > 0$, where $\mu(U)$ is the maximal eigenvalue of $F'(U)$, even if the same considerations could be made in the general case. The numerical flux (29) therefore reduces to

$$ (30) \qquad \mathcal{F}(U^-, U^+) = \frac{F(U^-) + F(U^+)}{2} - \frac{\lambda}{2}(U^+ - U^-) $$

and the Maxwellian functions (24) to

$$ (31) \qquad M_1(U) = \frac{1}{2} \left( U - \frac{F(U)}{\lambda} \right) \quad \text{and} \quad M_2(U) = \frac{1}{2} \left( U + \frac{F(U)}{\lambda} \right). $$

Finally, let us describe the discrete numerical boundary conditions in the case of a single interval. We assume that system (21) is complemented with Neumann boundary conditions on both extremities, which guarantee the mass conservation, that is to say $\partial_x \rho(t, 0) = q(t, 0) = 0$ and $\partial_x \rho(t, L) = q(t, L) = 0$ for all $t \ge 0$.

In the discrete setting, we have to define the quantities in the external ghost cells of the network, $C_0$ and $C_{J+1}$. To ensure the mass conservation also at the discrete level, we impose on both boundaries

$$ \left\{ \begin{aligned} \mathcal{F}^{\rho}(U_0, U_1) &= \frac{1}{2}(q_0 + q_1) - \frac{\lambda}{2}(\rho_1 - \rho_0) = 0, \\ \mathcal{F}^{\rho}(U_J, U_{J+1}) &= \frac{1}{2}(q_J + q_{J+1}) - \frac{\lambda}{2}(\rho_{J+1} - \rho_J) = 0, \end{aligned} \right. $$

see Eq. (30), which is satisfied if

$$ (32) \qquad \rho_0 = \rho_1 \text{ and } q_0 = -q_1 $$

and

$$ (33) \qquad \rho_J = \rho_{J+1} \text{ and } q_J = -q_{J+1}. $$
4.2. **The relaxation scheme on a 2-arc network.** Let us now go back to the case of a simple network with two arcs sketched in Figure 1.

We still denote by $\Delta x$ the spatial discretization step, that we assume to be the same for both arcs for simplicity reasons and by $\Delta t$ the time step. We discretize the left arc (resp. right arc) using $J_\ell$ (resp. $J_r$) cells

$$C_{j,\ell} = [x_{j-1/2,\ell}, x_{j+1/2,\ell}], 1 \le j \le J_{\ell} \quad (\text{resp. } C_{j,r} = [x_{j-1/2,r}, x_{j+1/2,r}], 1 \le j \le J_r)$$

centered at the nodes $x_{j,\ell} = (j - \frac{1}{2})\Delta x$, $1 \le j \le J_\ell$ (resp. $x_{j,r} = (j - \frac{1}{2})\Delta x$, $1 \le j \le J_r$).

On each arc we solve the system (21) by the approximation (26)-(30) with $\lambda = \lambda_\ell$ and $\lambda = \lambda_r$ for arc $\ell$ and arc $r$ respectively. Again, more general speeds could be considered, but we limit this presentation to the simplest case.

Boundary conditions at the extremities of the network are simply deduced from (32)-(33), that is to say

$$(34) \qquad \rho_{0,\ell} = \rho_{1,\ell} \quad \text{and} \quad q_{0,\ell} = -q_{1,\ell}.$$

and

$$(35) \qquad \rho_{J_r,r} = \rho_{J_r+1,r} \quad \text{and} \quad q_{J_r,r} = -q_{J_r+1,r}.$$

4.3. **Numerical treatment of the junction.** The current objective is to establish precise and consistent conditions at the junction of the two arcs, that is to say consistent values for the ghost cells at the junction, $C_{J_\ell+1,\ell}$ and $C_{0,r}$.

**Notations 4.1.** From now on we shall use the following notations, see notations in blue in Figure 4:

* $U_\ell^* = (\rho_\ell^*, q_\ell^*)$ and $U_r^* = (\rho_r^*, q_r^*)$ represent the unknown values at the junction on the ghost cell $C_{J_\ell+1,\ell}$ on arc $\ell$ and on the cell $C_{0,r}$ on arc $r$, respectively;

* $U_\ell^- = (\rho_\ell^-, q_\ell^-) = (\rho_{J_\ell,\ell}, q_{J_\ell,\ell})$ denotes the values on the last cell $C_{J_\ell,\ell}$ on arc $\ell$;

* $U_r^+ = (\rho_r^+, q_r^+) = (\rho_{1,r}, q_{1,r})$ stands for values on the first cell $C_{1,r}$ on arc $r$,

* the notation $v_{j,a}^n$ will be employed to denote variable $v$ ($v = \rho$ or $q$) with $j$ the cell index, $a$ denotes the arc and can take the values $a = \ell$ or $a = r$ and $n$ denotes the time index;

* in the relaxation setting, $f_{i,\ell}^-$ and $f_{i,r}^+$, $i = 1,2$, denote the left and right states and $f_{i,\ell}^*$ and $f_{i,r}^*$, $i = 1,2$, denote the left and right intermediate states, see Figure 4.

Therefore, to complement the scheme described in Section 4.2, we need now to find some values for the four quantities at the junction, i.e. $q_\ell^*, q_r^*, \rho_\ell^*$ and $\rho_r^*$. As stated above, we need to impose some transmission conditions at the junction. We consider a discrete version of the equality of the fluxes (7) in order to have the mass conservation at the junction and the (JTC) transmission conditions (8) linking the fluxes with the difference of densities. Namely, with the notations introduced above,

$$q_{\ell}^{*} = q_{r}^{*} = \kappa(\rho_{\ell}^{*} - \rho_{r}^{*}).$$

Since we have a total of 4 unknowns to be found at the junction, we need to impose two supplementary conditions. Let us now consider the relaxation system (27) on both arcs $\ell$ and $r$ with $c_{1,\ell} = -c_{2,\ell} = \lambda_\ell > 0$ and $c_{1,r} = -c_{2,r} = \lambda_r > 0$ respectively. By the exact Riemann solver (28) at the junction and considering that the junction stands for a shock with 0-velocity,
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Notations for the quantities at the junction : in blue, $(\rho_{\ell}^{-}, q_{\ell}^{-})$ and $(\rho_{r}^{+}, q_{r}^{+})$ are the left and right states and $(\rho_{\ell}^{*}, q_{\ell}^{*})$ and $(\rho_{r}^{*}, q_{r}^{*})$ the two intermediate states. In red, the two corresponding vectors $f_{1}$ and $f_{2}$ for the relaxation approximation.</figcaption></figure>

we get

$$ (36) \qquad f_{2,\ell}^{*} = f_{2,\ell}^{-} \quad \text{on arc } \ell, $$

$$ (37) \qquad f_{1,r}^{*} = f_{1,r}^{+} \quad \text{on arc } r, $$

where here $f_{i,\ell}^{*}$ (resp. $f_{i,r}^{*}$), $i = 1,2$ represents the (unknown) value of the solution for $-\lambda_{l} < x/t < 0$ (resp. for $0 < x/t < \lambda_{r}$), $f_{i,\ell}^{-}$ is the solution for $x/t < -\lambda_{\ell}$ on arc $\ell$ and $f_{i,r}^{+}$ the solution for $x/t > \lambda_{r}$ on arc $r$.

Switching to the first components of maxwellian functions given by (31), we can rewrite (36) and (37) as

$$ \left\{ \begin{aligned} \frac{1}{2} \left( \rho_{\ell}^{*} + \frac{q_{\ell}^{*}}{\lambda_{\ell}} \right) &= \frac{1}{2} \left( \rho_{\ell}^{-} + \frac{q_{\ell}^{-}}{\lambda_{\ell}} \right), \\ \frac{1}{2} \left( \rho_{r}^{*} - \frac{q_{r}^{*}}{\lambda_{r}} \right) &= \frac{1}{2} \left( \rho_{r}^{-} + \frac{q_{r}^{-}}{\lambda_{r}} \right), \end{aligned} \right. $$

which gives us the two extra conditions at the junction we were looking for.

In summary, at the junction, we need to solve the following linear system for the unknowns $q^{*} = q_{\ell}^{*} = q_{r}^{*}$, $\rho_{\ell}^{*}$ and $\rho_{r}^{*}$:

$$ (38) \qquad \left\{ \begin{aligned} q^{*} &= \kappa(\rho_{\ell}^{*} - \rho_{r}^{*}), \\ \rho_{\ell}^{*} + \frac{q_{\ell}^{*}}{\lambda_{\ell}} &= \rho_{\ell}^{-} + \frac{q_{\ell}^{-}}{\lambda_{\ell}}, \\ \rho_{r}^{*} - \frac{q_{r}^{*}}{\lambda_{r}} &= \rho_{r}^{-} + \frac{q_{r}^{-}}{\lambda_{r}}. \end{aligned} \right. $$
This system can be solved explicitely as :

$$ (39) \qquad \left\{ \begin{array}{l} q^* = q_\ell^* = q_r^* = \frac{\kappa}{\kappa(\lambda_\ell + \lambda_r) + \lambda_\ell \lambda_r} (\lambda_r q_\ell^- + \lambda_\ell q_r^+ + \lambda_\ell \lambda_r (\rho_\ell^- - \rho_r^+)), \\ \rho_r^* = \rho_r^+ + \frac{1}{\lambda_r} (q^* - q_r^+), \\ \rho_\ell^* = \rho_\ell^- - \frac{1}{\lambda_\ell} (q^* - q_\ell^-). \end{array} \right. $$

The solution can then be used as boundary values for the arcs on the ghost cells at the junction.

Notice that, the system given by (38) describing the quantities at the junction, can be recovered looking at the numerical flux $\mathcal{F}$ by considering directly the discretization of system (4)-(8).

We still use here the previous notations, namely, we denote by $U_\ell^- = (\rho_\ell^-, q_\ell^-) = (\rho_{J_\ell,\ell}, q_{J_\ell,\ell})$ (resp. $U_r^+ = (\rho_r^+, q_r^+) = (\rho_{1,r}, q_{1,r})$) the values on the last cell $C_{J_\ell,\ell}$ on arc $\ell$ (resp. on the first cell $C_{1,r}$ on arc $r$).

At the junction, the mass conservation at the discrete level gives the condition

$$ (40) \qquad \mathcal{F}^\rho(U_\ell^-, U_\ell^*) = \mathcal{F}^\rho(U_r^*, U_r^+) $$

that is coupled with the following relations, coming from the discretization of conditions (8) :

$$ \mathcal{F}^\rho(U_\ell^-, U_\ell^*) = \mathcal{F}^\rho(U_r^*, U_r^+) = \kappa(\rho_\ell^* - \rho_r^*). $$

Again, we have 4 unknowns at the junction and we need two supplementary conditions. We find them by imposing the equality between the numerical flux for the mass equation $F^\rho = \rho u = q$ and the numerical discretization for the momentum $q = \rho u$, namely :

$$ (41) \qquad \mathcal{F}^\rho(U_\ell^-, U_\ell^*) = q_\ell^*, \quad \mathcal{F}^\rho(U_r^*, U_r^+) = q_r^*. $$

We obtain therefore a linear system of 4 equations with 4 unknowns :

$$ (42) \qquad \left\{ \begin{array}{l} \mathcal{F}^\rho(U_\ell^-, U_\ell^*) = \mathcal{F}^\rho(U_r^*, U_r^+) = \kappa(\rho_\ell^* - \rho_r^*), \\ \mathcal{F}^\rho(U_\ell^-, U_\ell^*) = q_\ell^*, \\ \mathcal{F}^\rho(U_r^*, U_r^+) = q_r^*. \end{array} \right. $$

Note that system (42) holds for any numerical flux $\mathcal{F}$ and not only for the HLL flux.

In the case of the HLL numerical flux, the previous system resumes to

$$ (43) \qquad \left\{ \begin{array}{l} q^* = q_\ell^* = q_r^* = \kappa(\rho_\ell^* - \rho_r^*), \\ \frac{q_\ell^- + q^*}{2} - \frac{\lambda_\ell}{2}(\rho_\ell^* - \rho_\ell^-) = q^*, \\ \frac{q^* + q_r^+}{2} - \frac{\lambda_r}{2}(\rho_r^+ - \rho_r^*) = q^*. \end{array} \right. $$

which is the same as system (38).

**Remark 4.1.** Let us briefly compare the two junction conditions that have been considered in this paper, i.e. system (38) and system (18).

* Note that system (38) admits always a solution, both for subsonic and supersonic states, unlike system (18). It can therefore be seen as an extension of system (18) that works in both cases.

* Then, Riemann conditions leading to system (18) should not preserve the mass conservation at the discrete level, although conditions in system (38) are designed to do so.
* Finally, the reasoning for Riemann conditions at Subsection 3.2 only holds for a power $\gamma$ in Eq.(2) such that $1 < \gamma \le 3$, whereas system (38) holds a solution for every $\gamma > 1$. Note that the parameter $\gamma$ is only present in Eq.(38) through the values of $\lambda_{\ell}$ and $\lambda_{r}$.

**Remark 4.2.** In the subsonic region, the two last conditions given in (38) may be viewed as a local linearization of relations (20). In fact, substituting (20) with the tangent lines of the two Riemann invariants in the points $(\rho_{\ell}, q_{\ell}^{-})$ and $(\rho_{r}^{+}, q_{r}^{+})$ respectively, we get

$$ (44) \qquad \left\{ \begin{array}{l} q_{\ell}^{*}(\rho) = q_{\ell}^{-} + \mu_{1}(\rho_{\ell}^{-}, q_{\ell}^{-})(\rho - \rho_{\ell}^{-}), \\ q_{r}^{*}(\rho) = q_{r}^{+} + \mu_{2}(\rho_{r}^{+}, q_{r}^{+})(\rho - \rho_{r}^{+}). \end{array} \right. $$

On the other hand, the two last equations in (38) can be rewritten as

$$ (45) \qquad \left\{ \begin{array}{l} q_{\ell}^{*} = q_{\ell}^{-} - \lambda_{\ell}(\rho_{\ell}^{*} - \rho_{\ell}^{-}), \\ q_{r}^{*} = q_{r}^{+} + \lambda_{r}(\rho_{r}^{*} - \rho_{r}^{+}). \end{array} \right. $$

Then, the two unknown values $U_{\ell}^{*}$ and $U_{r}^{*}$ belong to the decreasing line $q(\rho) = q_{\ell}^{-} - \lambda_{\ell}(\rho - \rho_{\ell}^{-})$ and to the increasing line $q(\rho) = q_{r}^{+} + \lambda_{r}(\rho - \rho_{r}^{+})$ respectively. Since we are assuming $\mu_{1}(\rho_{\ell}^{-}, q_{\ell}^{-}) < 0$ and $\mu_{2}(\rho_{r}^{+}, q_{r}^{+}) > 0$, (44) reads as

$$ (46) \qquad \left\{ \begin{array}{l} q_{\ell}^{*} = q_{\ell}^{-} - |\mu_{1}(\rho_{\ell}^{-}, q_{\ell}^{-})|(\rho_{\ell}^{*} - \rho_{\ell}^{-}), \\ q_{r}^{*} = q_{r}^{+} + \mu_{2}(\rho_{r}^{+}, q_{r}^{+})(\rho_{r}^{*} - \rho_{r}^{+}). \end{array} \right. $$

The two systems (45) and (46) differ only in the choice of the angular coefficient of the two lines.

### 4.4. Consistency properties of the numerical scheme.

#### 4.4.1. Mass conservation.
First we prove the following property, which is equivalent to the continuous mass preservation (9).

**Proposition 4.1.** Let us consider system (4)-(8) set on a 2-arc network displayed in Figure1 with boundary conditions (6) on the outer nodes and junction conditions (8). The numerical scheme (26) with discrete boundary conditions (34)-(35) on the outer nodes and discrete junction conditions (42) is mass-preserving, that is to say the discrete mass

$$ m^n = h \sum_{j=1}^{J_{\ell}} \rho_{j,\ell}^n + h \sum_{j=1}^{J_{r}} \rho_{j,r}^n \text{ is independent of } n \in \mathbb{N}. $$

*Proof.* We consider the first component of equation (26). The conservation of the discrete mass is given by the boundary conditions for the outer nodes (34)-(35) and by condition (40) at the junction. $\square$

#### 4.4.2. Positivity of the solution.
Now, we prove that the positivity of the solution is preserved in the particular case of HLL numerical flux:

**Proposition 4.2.** Let us consider scheme (26) with flux (30), boundary conditions (34) – (35) and condition (43) at the junction in order to discretize system (4)-(8) on a 2-arc network. We define $\lambda_{\ell}^{n} > 0$ (resp. $\lambda_{r}^{n} > 0$) such that

$$ (47) \qquad \lambda_{\ell}^{n} \ge \max_{1 \le j \le J_{\ell}} |u_{j,\ell}^{n}| + \sqrt{p'(\rho_{j,\ell}^{n})} \quad (\text{resp. } \lambda_{r}^{n} \ge \max_{1 \le j \le J_{r}} |u_{j,r}^{n}| + \sqrt{p'(\rho_{j,r}^{n})}). $$
If the initial condition is positive, that is to say

$$ \rho_{j,l}^0 \ge 0 \text{ for all } 1 \le j \le J_l \text{ (resp. } \rho_{j,r}^0 \ge 0 \text{ for all } 1 \le j \le J_r), $$
then the solution remains positive in time, that is to say

$$ \rho_{j,l}^n \ge 0 \text{ for all } 1 \le j \le J_l \text{ (resp. } \rho_{j,r}^n \ge 0 \text{ for all } 1 \le j \le J_r), \text{ for all } n \in \mathbb{N}. $$

Notice that condition (47) is just a formulation in this case of condition (25), which is necessary to obtain the stability of the numerical scheme.

*Proof.* We prove the proposition by induction.

Assume that $\rho_{j,l}^n \ge 0$ for all $1 \le j \le J_l$ (resp. $\rho_{j,r}^n \ge 0$ for all $1 \le j \le J_r$). Since $\lambda_l^n \ge |u_{j,l}^n|$ for all $1 \le j \le J_l$ (resp. $\lambda_r^n \ge |u_{j,r}^n|$ for all $1 \le j \le J_r$), we use the explicit expressions obtained at Eq.(39) to prove that $\rho_r^{*,n} \ge 0$ and $\rho_l^{*,n} \ge 0$. From these two inequalities, we demonstrate that $\rho_{j,l}^{n+1} \ge 0$ for all $1 \le j \le J_l$ (resp. $\rho_{j,r}^{n+1} \ge 0$ for all $1 \le j \le J_r$). $\square$

## 5. NUMERICAL TESTS IN THE CASE OF A 2-ARC NETWORK

In this section, we examine the characteristics of the numerical method proposed in Section 4 and we evaluate its accuracy by examining a simple network consisting of two arcs. The pressure is given by the pressure law given at Eq.(2) with $\gamma = 2$, unless another value is specified. We therefore consider the numerical scheme defined by Eq.(26) with flux (29) on the two arcs, conditions (34) - (35) at the extremities of the network and conditions (38) at the junction; this scheme will be referred as the HLL-(38) scheme in the following.

We consider a uniform spatial grid with a fixed spatial step $\Delta x = 0.05$ on the two arcs with equal length of $L_l = L_r = 2$. Then, the number of cells verify $J_l = J_r = J$. The time step is dynamically computed at each time iteration $n$ by the CFL condition

$$ \Delta t^n \le \frac{\Delta x}{\max\{\lambda_l^n, \lambda_r^n\}}, $$

where $\lambda_l^n$ and $\lambda_r^n$ have been defined in Proposition 4.2, such that the subscharacteristic condition (47) is verified.

### 5.1. Numerical convergence results for the numerical scheme HLL-(38). Let the following Gaussian function be given

$$ g(x; \sigma, c) = e^{-(x-c)^2/2\sigma^2}. $$

To study the numerical order of convergence of the scheme applied to the whole network, we set as initial data (at time $t = 0$), on the left arc, for $x \in [0, L_l]$,

$$ (48) \qquad \rho_{\ell}(0, x) = 1.5 + g(x, 0.2, 0.8) \text{ and } q_{\ell}(0, x) = 0.5, $$

and on the right-arc, for $x \in [0, L_r]$,

$$ (49) \qquad \rho_r(0, x) = 1 + g(x, 0.1, 0.5) \text{ and } q_r(0, x) = 0.5. $$

This initial condition is subsonic all over the network and verifies at the junction the continuity of the flux and the (JTC) condition (8) with $\kappa = 1$, i.e.

$$ q_{\ell}(0, L_{\ell}) = q_{r}(0, 0) = \rho_{\ell}(0, L_{\ell}) - \rho_{r}(0, 0) = 0.5. $$
The accuracy is measured by the $L^1$ norm of the error at the final time $T = 1$, defined as

$$ \text{err}_{\rho,\ell} = \Delta x \sum_{j=1}^{J_{\ell}} |\rho_{j,\ell}^{N} - \rho_{j,\ell}^{\text{ref}}|, \quad \text{err}_{q,\ell} = \Delta x \sum_{j=1}^{J_{\ell}} |q_{j,\ell}^{N} - q_{j,\ell}^{\text{ref}}|, $$

with the same definition for arc $r$, where $(\rho^{\text{ref}}, q^{\text{ref}})$ is a reference solution obtained with a small $\Delta x \sim 10^{-5}$ on both arcs $\ell$ and $r$. Considering the system (38) at the junction, Table 1 shows that the conditions imposed at the junction do not affect the convergence order of the scheme, which remains first order.

Furthermore, we show that the system (38) does not alter the asymptotic behavior of the solution along the entire 2-arc network. We consider again the initial data (48)–(49) and calculate the rate of convergence of the solution to the expected asymptotic state

$$ (50) \qquad \bar{\rho}_{\text{asympt}} = \frac{1}{L_{\ell} + L_{r}} \int_{0}^{L_{\ell}+L_{r}} (\rho_{\ell}(0,x) + \rho_{r}(0,x)) dx, \quad \bar{q}_{\text{asympt}} = 0. $$

Given N data points $(t_i, e_{\rho,\ell}(t_i))_{i=1,\dots,N}$, we shall define $\gamma$ and $C$ as the solution of the following least square problem,

$$ (51) \qquad \min_{C, \gamma} \sum_{i=1}^{N} |\ln(e_{\rho,\ell}(t_i)) - \ln(C t_i^{-\gamma})|^2, $$

where $e_{\rho,\ell}(t_i) = \|\rho_{\ell}(t_i, \cdot) - \bar{\rho}_{\text{asympt}}\|_{L^1,L^{\infty}}$, with the same definition for $e_{\rho,r}$, $e_{q,\ell}$ and $e_{q,r}$. Table 2 presents the optimal values of $\gamma$ and $C$. Once again, the system imposed at the junction does not alter the asymptotic behavior in time of the solution, see [31], where it is proved that the solution decays as $t^{-1}$, in the case of a single interval.

<table><caption>TABLE 1. $L^1$ error and order of convergence of the numerical scheme HLL-(38) on both arcs starting from the regular initial data (48)–(49).</caption><thead><tr><th></th><th colspan="2">ρ<sub>ℓ</sub></th><th colspan="2">ρ<sub>r</sub></th><th colspan="2">q<sub>ℓ</sub></th><th colspan="2">q<sub>r</sub></th></tr><tr><th>Δx</th><th>err<sub>ρ,ℓ</sub></th><th>order</th><th>err<sub>ρ,r</sub></th><th>order</th><th>err<sub>q,ℓ</sub></th><th>order</th><th>err<sub>q,r</sub></th><th>order</th></tr></thead><tbody><tr><td>0.0063</td><td>0.02651</td><td>0.70286</td><td>0.04090</td><td>0.81492</td><td>0.05104</td><td>0.77079</td><td>0.06378</td><td>0.63251</td></tr><tr><td>0.0031</td><td>0.01495</td><td>0.82633</td><td>0.02151</td><td>0.92712</td><td>0.02727</td><td>0.90406</td><td>0.03558</td><td>0.84207</td></tr><tr><td>0.0016</td><td>0.00743</td><td>1.00890</td><td>0.01030</td><td>1.06292</td><td>0.01407</td><td>0.95527</td><td>0.01742</td><td>1.02998</td></tr><tr><td>0.0008</td><td>0.00334</td><td>1.15396</td><td>0.00453</td><td>1.18482</td><td>0.00655</td><td>1.10264</td><td>0.00775</td><td>1.16828</td></tr></tbody></table>

<table><caption>TABLE 2. Optimal values of $\gamma$ and $C$ given by the least square problem (51) with $e_{\rho,\ell}(t_i) = \|\rho_{\ell}(t_i, \cdot) - \bar{\rho}_{\text{asympt}}\|_{L^1,L^{\infty}}$ for $i = 1, \dots, N$ computed in $L^1$ and $L^{\infty}$ norms, same definition for $e_{\rho,r}$, $e_{q,\ell}$ and $e_{q,r}$.</caption><thead><tr><th></th><th colspan="2">ρ<sub>ℓ</sub></th><th colspan="2">ρ<sub>r</sub></th><th colspan="2">q<sub>ℓ</sub></th><th colspan="2">q<sub>r</sub></th></tr><tr><th></th><th>C</th><th>γ</th><th>C</th><th>γ</th><th>C</th><th>γ</th><th>C</th><th>γ</th></tr></thead><tbody><tr><td>L<sup>1</sup></td><td>0.2356</td><td>1.0248</td><td>0.2301</td><td>1.0274</td><td>0.5642</td><td>0.9802</td><td>0.5193</td><td>0.9544</td></tr><tr><td>L<sup>∞</sup></td><td>0.2295</td><td>1.0558</td><td>0.2295</td><td>1.0558</td><td>0.5662</td><td>1.0493</td><td>0.5446</td><td>1.0371</td></tr></tbody></table>
5.2. **Comparison between the two numerical schemes HLL-(38) and HLL-(18).** Now, we also consider the so-called HLL-(18) scheme, which is defined as the previous HLL-(38) scheme, unlike the junction condition which is taken here as the solution of system (18) defined at Sec.3 thanks to Riemann invariant considerations.

In this section we compare the dynamics obtained by the two schemes HLL-(38) and HLL-(18), examining the subsonic and supersonic cases. In all tests, we start from data verifying condition (8) with $\kappa = 1$ at the junction, i.e.

$$q_{\ell}(0, L_{\ell}) = q_{r}(0, 0) = \rho_{\ell}(0, L_{\ell}) - \rho_{r}(0, 0).$$

In the following figures, these initial data are drawn in dashed black lines. We consider the following test problems:

**Test C1:** Case of a subsonic initial datum, constant on both arcs,

$$\begin{cases} \rho_{\ell}(0, x) = 4.5, & q_{\ell}(0, x) = 0.5, & \text{for all } x \in [0, L_{\ell}], \\ \rho_{r}(0, x) = 4, & q_{r}(0, x) = 0.5, & \text{for all } x \in [0, L_{r}]. \end{cases}$$

As can be seen in Figure 5, the two schemes HLL-(38) and HLL-(18) give exactly the same results on both the arcs and the junction. We also notice that solutions are converging in large times to the asymptotic constant state described at Eq.(50).

<figure><img src="image_6.png" /><figcaption>FIGURE 5. **Test C1.** Comparison of the solutions given by scheme HLL-(38) (in solid red line) and by scheme HLL-(18) (in dashed blue line), starting from a constant subsonic initial datum on both arcs. The blue and red dots represent the values at the junction. On the left, function $\rho$ and on the right, function $q$ for different times $t = 0.25, 0.5$ and $0.75$ with final time $T = 1$.</figcaption></figure>

**Test C2:** Case of a non constant initial datum given by

$$\begin{cases} \rho_{\ell}(0, x) = 0.4 + g(x, 0.2, 0.8) & \text{and} & q_{\ell}(0, x) = 0.1, & \text{for all } x \in [0, L_{\ell}], \\ \rho_{r}(0, x) = 0.3 + g(x, 0.1, 0.5) & \text{and} & q_{r}(0, x) = 0.1, & \text{for all } x \in [0, L_{r}], \end{cases}$$

which is subsonic at the junction. As can be seen in Figure 6, the two approaches show some differences in the calculated values at the junction (look particularly at Figure 6-(a) at time $T/4$), but these do not propagate along the arcs and vanish over time. Here again the solutions converge to the expected asymptotic profile for large times.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Test C2. Comparison of the solutions given by scheme HLL-(38) (in solid red line) and by scheme HLL-(18) (in dashed blue line), starting from a subsonic non constant initial datum (in dashed black line). The blue and red dots represent the values at the junction. On the left, function $\rho$ and on the right, function $q$ for different times $t = 0.25, 0.5$ and $0.75$ with final time $T = 1$.</figcaption></figure>

**Test C3:** Case of a constant initial datum supersonic to the right of the junction (right-arc) and subsonic on the left,

$$ \begin{cases} \rho_{\ell}(0, x) = 2.5, & q_{\ell}(0, x) = 0.5, & \text{for all } x \in [0, L_{\ell}], \\ \rho_{r}(0, x) = 2, & q_{r}(0, x) = 0.5, & \text{for all } x \in [0, L_{r}]. \end{cases} $$

As can be seen in Figure 7, on the right arc the solution at the junction remains supersonic and the two algorithms give the same results. This is an example where the datum on the right is supersonic but with positive velocity, namely $q_r(0, L_r) > 0$ and $q_r(0, L_r) > q_s^+(\rho_r(0, L_r))$. In this case, the Riemann based HLL-(18) scheme selects the solution on the right invariant curve, see Subsec.3.3, and so does the HLL-(38) scheme.

**Test C4:** Case of a non constant initial datum given by

$$ \begin{cases} \rho_{\ell}(0, x) = 0.1 + g(x, 0.2, 0.8) & \text{and} & q_{\ell}(0, x) = -0.2, & \text{for all } x \in [0, L_{\ell}], \\ \rho_{r}(0, x) = 0.3 + g(x, 0.1, 0.5) & \text{and} & q_{r}(0, x) = -0.2, & \text{for all } x \in [0, L_{r}], \end{cases} $$

which is supersonic to the left of the junction (left-arc) and subsonic to the right of the junction. As can be seen in Figure 8, the two methods produce different results, particularly on the left arc. This is even more visible in Figure 9, where we plot the evolution in time of the four junction values $\rho_{\ell}^*, \rho_r^*, q_{\ell}^*$ and $q_r^*$. We observe, on the left arc the flow transition from the supersonic region to the subsonic one, and that the two solvers, HLL-(38) and HLL-(18), produce some slightly different results, especially before the transition. This is an example where the datum on the left is supersonic with negative velocity, namely $q_{\ell}(0, L_{\ell}) < 0$ and $q_{\ell}(0, L_{\ell}) < q_s^-(\rho_{\ell}(0, L_{\ell}))$. In this case the Riemann based algorithm selects the solution on the right invariant curve, see Subsection 3.3, while the relaxation one (38) gives a result which is not quite the same.
<figure><img src="image_8.png" /><figcaption>FIGURE 7. **Test C3.** Comparison of the solutions given by scheme HLL-(38) (in solid red line) and by scheme HLL-(18) (in dashed blue line), starting from a constant initial datum supersonic on the right-arc and subsonic on the left-arc. The blue and red dots represent the values at the junction and the dashed black line the initial value for $q$. On the left, function $\rho$ and on the right, function $q$ for different times $t = 0.25, 0.5$ and $0.75$ with final time $T = 1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 8. **Test C4.** Comparison of the solutions given by scheme HLL-(38) (in solid red line) and by scheme HLL-(18) (in dashed blue line), starting from a non-constant initial datum supersonic on the left-arc and subsonic on the left-arc (in dashed black line). The blue and red dots represent the values at the junction. On the left, function $\rho$ and on the right, function $q$ for different times $t = 0.25, 0.5$ and $0.75$ with final time $T = 1$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Test C4. Evolution with respect to time of the 4 junction values: $\rho_{*,\ell}$ on the top left subfigure, $\rho_{*,r}$ on the top right, $q_{*,\ell}$ on the bottom left and $q_{*,r}$ on the bottom right, computed by scheme HLL-(38) (in solid red line) and by scheme HLL-(18) (in dashed blue line). The two dotted lines in the figures on the bottom delineate the boundary between the subsonic and supersonic regions (the sonic curve given for the corresponding value of $\rho$), in pink for the HLL-(38) scheme and in light blue for the HLL-(18) scheme. The regime is subsonic inside these two curves and supersonic outside.</figcaption></figure>

5.3. **Study of different values for parameter $\kappa$.** In this test case, we show the numerical dynamics obtained by scheme HLL-(38) for two different values of parameter $\kappa$. We assume a constant initial datum given by

$$ \begin{cases} \rho_{\ell}(0, x) = 2.5, & \text{for all } x \in [0, L_{\ell}], \\ \rho_{r}(0, x) = 0.5, & \text{for all } x \in [0, L_{r}], \end{cases} $$

and

$$ q_{\ell}(0, x) = q_{r}(0, x) = \kappa(\rho_{\ell}(0, x) - \rho_{r}(0, x)) = 2\kappa. $$

In Figure 10, we show the density evolution obtained with the two values $\kappa = 100$ and $\kappa = 0.1$.

On the one hand, as expected, for $\kappa = 100$ the density values at the junction start out different and then tend to converge until they coincide. This is coherent since for large values for $\kappa$, we expect the system to behave similarly to the system with the continuity of densities at the junction, that is to say the junction condition that is usually considered in other articles. On the other hand, for small value of $\kappa = 0.1$, the density remains discontinuous at the interface.

Note that the Riemann-based HLL-(18) scheme does not provide a solution to test case with $\kappa = 100$, whereas scheme HLL-(38) does; it would be necessary to expand the region of admissible solutions and add further conditions to uniquely define the unknown values at the junction.

5.4. **Entropy dissipation at the discrete level.** In this subsection, following the discussion in Subsec. 2.1, we investigate the evolution of the total entropy in the network,

$$ (52) \qquad S(t) = \int_{0}^{L_{\ell}+L_{r}} \eta(\rho(t, x), q(t, x))dx \approx S^{n} = \Delta x \sum_{j=0}^{J_{\ell}+J_{r}} \eta_{j}^{n}, $$
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Different values for parameter $κ$. Density solution obtained with the approximation HLL-(38), starting from a constant initial datum on both arcs (in solid black line), for two different values of parameter $κ$: $κ = 100$ (dotted orange line) and $κ = 0.1$ (dotted green line), for different times $t = 0.5, 1$ and $1.5$ with final time $T = 2$.</figcaption></figure>

where $η$ is the entropy function defined in (10) and $η_j^n = η(ρ_j^n, q_j^n)$ for $j$ running along the nodes of the two arcs forming the network. We assume two different constant initial data which verify (8):

$$q_{\ell}(0, x) = q_{r}(0, x) = \kappa(\rho_{\ell}(0, x) - \rho_{r}(0, x)),$$

with $κ = 1$, and which differ in the amplitude of the density jump at the junction:

(E1) $ρ_ℓ(0, x) = 1.5$ and $ρ_r(0, x) = 1$ such that $q_ℓ(0, x) = q_r(0, x) = ρ_ℓ(0, x) - ρ_r(0, x) = 0.5$.

(E2) $ρ_ℓ(0, x) = 6$ and $ρ_r(0, x) = 1$ such that $q_ℓ(0, x) = q_r(0, x) = ρ_ℓ(0, x) - ρ_r(0, x) = 5$.

We recall here the entropy dissipation inequality (12): $S(t) \le S(0) - \int_0^t \Delta G(s)ds$.

As can be seen in Figure 11-(right) and Figure 12-(right), for both initial data the total entropy (52) is decreasing in time and numerically verify a discrete version of the inequality (12), i.e.

$$(53) \quad S^n = \Delta x \sum_{j=0}^{J_{\ell}+J_r} \eta_j^n \le \Delta x \sum_{j=0}^{J_{\ell}+J_r} \eta_j^0 - \Delta t \sum_{k=0}^{n} (G_{\ell}^{*,k} - G_{r}^{*,k}) = S^0 - \Delta t \sum_{k=0}^{n} \Delta G^{*,k},$$

where $G_{\ell}^{*,n} = G(ρ_{\ell}^{*,n}, q_{\ell}^{*,n})$ and $G_{r}^{*,n} = G(ρ_{r}^{*,n}, q_{r}^{*,n})$ represent the values of the entropy flux, as defined in (10), computed at the left and right junction states obtained from the solution of system (38) at time step $n$. As in Subsec. 2.1, we also define the quantity $\Delta G^{*,n} := G_{\ell}^{*,n} - G_{r}^{*,n}$, which is a discrete version at time $t_n$ of quantity $\Delta G$.

Specifically, for the small jump assumed in test case (E1), the quantity $\Delta G$, displayed in Figure 11-(a), continuous line, is positive for all time $t > 0$ as expected from the discussion in Subsec. 2.1, thus ensuring the dissipation of entropy with respect to time. We notice in Figure
11-(b) that the estimate $S(0) - \int_{0}^{t} \Delta G(s)ds$ is above the entropy for all times, as expected, and is equal up to a constant to the integral displayed on Figure 11-(a). We also observe that in that case the estimate is sharp.

However, for big density jumps at the junction, like in test case (E2), the variation $\Delta G(t)$ is negative for all times, as displayed in Figure 12-(a), continuous line. But, even if the right hand part of the inequality (53) is increasing, and the bound from above $S(0) - \int_{0}^{t} \Delta G(s)ds$ is therefore not accurate, the total entropy remains decreasing in time, see Figure 12-(b). We can interpret this as follows : the dissipation of entropy along the arcs is large enough to compensate the small enough non-dissipation of entropy at the junction.

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Entropy dissipation, test case (E1), small jump at the junction.</figcaption></figure>

(a) Plot of $\Delta G$ with respect to time (solide black line and range of values on the left) and of $-\int_{0}^{t} \Delta G(s) ds$ (dotted grey line and range of values on the right).

(b) Plot of the total entropy $S$ (solid black line) as a function of time and of the right hand side $S(0) - \int_{0}^{t} \Delta G(s)ds$ of inequality (53) (dotted grey line).

5.5. Simulations in the case $\gamma = 4/3$. To conclude this part of the numerical tests, let us consider the problem with a different value of $\gamma$ in the pressure law (2). We take $\gamma = 4/3$ and the constant initial datum (same as (E2)) given by

$$ \begin{cases} \rho_{\ell}(0, x) = 6, & q_{\ell}(0, x) = 5, & \text{for all } x \in [0, L_{\ell}], \\ \rho_{r}(0, x) = 1, & q_{r}(0, x) = 5, & \text{for all } x \in [0, L_{r}]. \end{cases} $$

Note that here the simulations are performed using only the scheme HLL-(38), since the scheme HLL-(18) is not able to handle this case (there is no simple solution to the corresponding Riemann problem).

The results are shown in Figure 13. We notice first that the evolution is much slower than in the case $\gamma = 2$ and that function $q$ is decreasing towards 0 but very slowly. We also observe that before going asymptotically towards a constant state, the density $\rho$ tends to form a bump like solution on the right arc of the network, see the function $\rho$ at different times on the left subfigure.
<figure><img src="image_13.png" /><figcaption>FIGURE 12. **Entropy dissipation**, test case (E2), large jump at the junction.</figcaption></figure>

(a) Plot of $\Delta G$ with respect to time (solide black line and range of values on the left) and of $-\int_0^t \Delta G(s) ds$ (dotted grey line and range of values on the right).

(b) Plot of the total entropy $S$ (solid black line) as a function of time and of the right hand side $S(0) - \int_0^t \Delta G(s) ds$ of inequality (53) (dotted grey line).

<figure><img src="image_14.png" /><figcaption>FIGURE 13. **Simulation for $\gamma = 4/3$**. Solutions given by scheme HLL-(38) (in solid red line), starting from a constant initial datum (in dashed black line) on both arcs. The red dots represent the values at the junction. On the left, function $\rho$ and on the right, function $q$ for different times $t = 0.25, 0.5$ and $0.75$ with final time $T = 1$.</figcaption></figure>

## 6. EXTENSIONS TO A GENERAL NETWORK

To begin with, we consider an extension of what has been previously done to more complicated network, i.e. networks with any number of arcs where nodes can be connected to any number of arcs.
Let us define a network or a connected graph $G = (\mathcal{N}, \mathcal{A})$, as composed of two finite sets, a set of nodes (or vertices) $\mathcal{N}$ and a set of $N$ arcs $\mathcal{A}$, such that an arc connects a pair of nodes. Since arcs are bidirectional the graph is non-oriented, but we need to fix an artificial orientation in order to fix a sign to the velocities. The network is therefore composed of "oriented" arcs and there are two different types of intervals at a node $p \in \mathcal{N}$: incoming ones – the set of these intervals is denoted by $\mathcal{I}_p$ – and outgoing ones – whose set is denoted by $\mathcal{O}_p$. See, for example, Figure 14 for the representation of a 4-arc and 1node-network, with 2 incoming arcs and 2 outgoing arcs, where $1, 2 \in \mathcal{I}_1$ and $3, 4 \in \mathcal{O}_1$.

**6.1. Extension of junction conditions (43).** Let us consider now a node $p \in \mathcal{N}$, which is the junction of $N_p$ arcs, denoted by $i$, with $1 \le i \le N_p$. We denote by $\mathcal{I}_p$ the set of incoming arcs and by $\mathcal{O}_p$ the set of outgoing arcs.

<figure><img src="image_15.png" /><figcaption>FIGURE 14. An example of a 4 arc network with 2 incoming arcs and 2 outgoing arcs</figcaption></figure>

We will first generalize at the continuous level conditions (6) for the outer nodes and then conditions (8) for junction nodes. Then, we give a numerical discretization of these conditions, namely a generalization of discrete conditions (34)-(35) on outer nodes and conditions (43) at junction nodes. Note that the "discrete equality junction conditions" (43) are straightforward to generalize unlike Riemann invariant junction conditions (18), which would be much more involved to consider for general networks.

Therefore, the generalization of conditions (6) at the outer nodes reads as :

$$ (54) \qquad \left\{ \begin{array}{ll} \partial_x \rho_i(\cdot, L_i) = 0 \text{ and } q_i(\cdot, L_i) = 0, & \text{if } i \in \mathcal{O}_p, \\ \partial_x \rho_i(\cdot, 0) = 0 \text{ and } q_i(\cdot, 0) = 0, & \text{if } i \in \mathcal{I}_p, \end{array} \right. $$

We now consider the following generalization of conditions (8) at the junction node $p \in \mathcal{N}$, which says that the flux is equal to a linear combination of the different differences of densities
on the various arcs :

$$
(55) \quad \left\{
\begin{array}{ll}
-q_i(\cdot, 0) = \sum_{j \in \mathcal{I}_p} \kappa_{i,j}(\rho_i(\cdot, 0) - \rho_j(\cdot, L_j)) + \sum_{j \in \mathcal{O}_p} \kappa_{i,j}(\rho_i(\cdot, 0) - \rho_j(\cdot, 0)) & \text{if } i \in \mathcal{O}_p, \\
q_i(\cdot, L_i) = \sum_{j \in \mathcal{I}_p} \kappa_{i,j}(\rho_i(\cdot, L_i) - \rho_j(\cdot, L_j)) + \sum_{j \in \mathcal{O}_p} \kappa_{i,j}(\rho_i(\cdot, L_i) - \rho_j(\cdot, 0)) & \text{if } i \in \mathcal{I}_p,
\end{array}
\right.
$$

where $\kappa_{i,j} \ge 0$ are generalizations of permeability coefficients between two arcs.

In order to have the mass conservation property at the junction $p \in \mathcal{N}$, we also impose a
generalization of condition (7), which reads as :

$$
(56) \qquad \sum_{i \in \mathcal{O}_p} -q_i(\cdot, 0) + \sum_{i \in \mathcal{I}_p} q_i(\cdot, L_i) = 0.
$$

If we impose that the transmission coefficients are symmetric, namely

$\kappa_{i,j} = \kappa_{j,i},$

condition (56) is automatically satisfied. This condition ensures that the total mass of the
whole system is conserved, namely

$$
(57) \qquad \frac{d}{dt} m(t) = 0, \text{ with } m(t) = \sum_{i \in \mathcal{A}} \int_{0}^{L_i} \rho_i(t, x) \, dx.
$$

Now, let us consider the discretization of the previous conditions on the outer nodes and at
the junction, following the ideas of Sec.4.

For the outer nodes, we set, in the spirit of conditions (34)-(35) :

$$
(58) \quad \begin{cases} \rho_{0,i} = \rho_{1,i} \quad \text{and} \quad q_{0,i} = -q_{1,i} & \text{if } i \in \mathcal{O}_p, \\ \rho_{J_i,i} = \rho_{J_{i-1},i} \quad \text{and} \quad q_{J_i,i} = -q_{J_{i-1},i} & \text{if } i \in \mathcal{I}_p. \end{cases}
$$

Now, at each node $p \in \mathcal{N}$, we will need to solve a linear system. Namely, at node $p \in \mathcal{N}$, considering the junction conditions with $N_p$ arcs, we need to compute $2N_p$ values, that we denote $\rho_i^*$ and $q_i^*$ for $1 \le i \le N_p$. Equations (55) give at the discrete level the following $N_p$ equations :

$$
(59) \quad \left\{
\begin{array}{ll}
-q_i^* = \sum_{j \in \mathcal{I}_p} \kappa_{i,j} (\rho_i^* - \rho_j^*) + \sum_{j \in \mathcal{O}_p} \kappa_{i,j} (\rho_i^* - \rho_j^*) & \text{if } i \in \mathcal{O}_p, \\
q_i^* = \sum_{j \in \mathcal{I}_p} \kappa_{i,j} (\rho_i^* - \rho_j^*) + \sum_{j \in \mathcal{O}_p} \kappa_{i,j} (\rho_i^* - \rho_j^*) & \text{if } i \in \mathcal{I}_p.
\end{array}
\right.
$$

We obtain $N_p$ other equations by imposing a generalization of equation (41) :

$$
(60) \quad \left\{
\begin{array}{ll}
\mathcal{F}^{\rho}(U_{1,i}, U_i^*) = q_i^*, & \text{if } i \in \mathcal{O}_p, \\
\mathcal{F}^{\rho}(U_{J_i,i}, U_i^*) = q_i^*, & \text{if } i \in \mathcal{I}_p.
\end{array}
\right.
$$

In the case of the HLL flux, the system (60) becomes

$$
\left\{
\begin{array}{ll}
\displaystyle \frac{q_i^* + q_{1,i}}{2} + \frac{\lambda_i}{2}(\rho_i^* - \rho_{1,i}) = q_i^*, & \text{if } i \in \mathcal{O}_p, \\
\\
\displaystyle \frac{q_i^* + q_{J_i,i}}{2} - \frac{\lambda_i}{2}(\rho_i^* - \rho_{J_i,i}) = q_i^*, & \text{if } i \in \mathcal{I}_p.
\end{array}
\right.
$$

or equivalently

$$
(61) \quad \left\{
\begin{array}{ll}
q_i^* = q_{1,i} + \lambda_i (\rho_i^* - \rho_{1,i}), & \text{if } i \in \mathcal{O}_p, \\
q_i^* = q_{J_i,i} - \mu_i (\rho_i^* - \rho_{J_i,i}), & \text{if } i \in \mathcal{I}_p.
\end{array}
\right.
$$
Combining Eq.(61) and (59), we obtain the following system composed of $N_p$ equations:

$$ (62) \quad \left\{ \begin{array}{ll} \sum_{j \in \mathcal{I}_p} \kappa_{i,j} (\rho_i^* - \rho_j^*) + \sum_{j \in \mathcal{O}_p} \kappa_{i,j} (\rho_i^* - \rho_j^*) + \mu_i \rho_i^* = -q_{1,i} + \mu_i \rho_{1,i} & \text{if } i \in \mathcal{O}_p, \\ \sum_{j \in \mathcal{I}_p} \kappa_{i,j} (\rho_i^* - \rho_j^*) + \sum_{j \in \mathcal{O}_p} \kappa_{i,j} (\rho_i^* - \rho_j^*) + \mu_i \rho_i^* = q_{J_i,i} + \mu_i \rho_{J_i,i} & \text{if } i \in \mathcal{I}_p. \end{array} \right. $$

Therefore, at node $p \in \mathcal{N}$, the vector $(\rho_1^*, \dots, \rho_{N_p}^*)$ is solution to a linear system, which is invertible, as soon as $\kappa_{i,j} \ge 0$ and $\mu_i > 0$, since the matrix

$$ (63) \quad \left( \begin{array}{cccc} \sum_{j \neq 1} \kappa_{1,j} + \mu_1 & & & \\ & \ddots & -\kappa_{i,j} & \\ & & \ddots & \\ -\kappa_{i,j} & & & \sum_{j \neq N_p} \kappa_{N_p,j} + \mu_{N_p} \end{array} \right), $$

is a strictly diagonally dominant matrix. Values for momenta $q_i^*$, $1 \le i \le N_p$ are then given by Eq.(61).

**6.2. Properties of the numerical scheme.** It is now easy to generalize the mass conservation and the positivity of solution, already demonstrated at Sec.4.4 in the case of a 2 arc network, to the case of a $N$ arc network. However, the entropy dissipation property is far more complicated and will not be treated here. First, we prove the following property, which is equivalent to the continuous mass preservation (57).

**Proposition 6.1.** Let us consider system (4) set on a $N$-arc network with boundary conditions (54) on the outer nodes and junction conditions (55) with symmetric coefficients, that is to say $\kappa_{i,j} = \kappa_{j,i}$. The numerical scheme (26) with discrete boundary conditions (58) on the outer nodes and discrete junction conditions (59) is mass-preserving, that is to say the discrete mass

$$ m_{tot}^n = h \sum_{i=1}^{N} \sum_{j=1}^{J_i} \rho_{j,i}^n \text{ is independent of } n \in \mathbb{N}. $$

*Proof.* We consider the first component of equation (26). The conservation of mass is given by Eq.(58) and (59) with $\kappa_{i,j} = \kappa_{j,i}$, which leads to

$$ \sum_{i \in \mathcal{O}_p} -q_i^* + \sum_{i \in \mathcal{I}_p} q_i^* = 0. $$

□

Now, we prove the generalization of the positivity of the solution in the particular case of HLL numerical flux:

**Proposition 6.2.** Let us consider scheme (26) with flux (30) on a $N$ arc network, boundary conditions (58) and conditions (59)-(60) at the junction in order to discretize system (4)-(54)-(55) on a general network. We define $\lambda_i^n > 0$ such that

$$ \lambda_i^n \ge \max_{1 \le j \le J_i} |u_{j,i}^n| + \sqrt{p'(\rho_{j,i}^n)}, \text{ for } 1 \le i \le N. $$
If the initial condition is positive, that is to say

$\rho_{j,i}^0 \ge 0$ for all $1 \le i \le N$, for all $1 \le j \le J_i$,

then the solution remains positive in time, that is to say

$\rho_{j,i}^n \ge 0$ for all $1 \le i \le N$, for all $1 \le j \le J_i$, for all $n \in \mathbb{N}$.

*Proof.* We prove the proposition by induction, following the proof of Prop.4.2. Assume that $\rho_{j,i}^n \ge 0$ for all $1 \le i \le N$, $1 \le j \le J_i$.

We consider node $p \in \mathcal{N}$. Since $\lambda_i^n \ge |u_{j,i}^n|$ for all $1 \le j \le J_i$, we can show that the components of the second member of linear system (62) are nonnegative. Moreover, since $\kappa_{i,j} \ge 0$ and $\mu_i^n > 0$, it can be proved easily that the matrix (63) of linear system (62) is a monotone matrix. Consequently the solution $(\rho_1^*, \dots, \rho_{N_p}^*)$ has positive components and therefore, we can deduce that $\rho_{j,i}^{n+1} \ge 0$ for all $1 \le i \le N$, $1 \le j \le J_i$. $\square$

**6.3. Numerical tests.** We present now some numerical tests for two different general networks.

**Network 1:** The first test is performed on a 12-arc network with 12 nodes, see [12] : 4 junction nodes and 8 outer nodes, see Figure15. The 4 junction nodes are forming a square and each of them is connected to 2 outer nodes. Each junction node connects therefore 4 arcs together and the permeability coefficients $\kappa_{i,j}$ are given in the following matrix, the same for each junction node, note that $\kappa_{i,j} = \kappa_{j,i}$ :

$$K = \begin{pmatrix} 0 & 0.3 & 0.2 & 0.5 \\ 0.3 & 0 & 0.2 & 0.1 \\ 0.2 & 0.2 & 0 & 0.2 \\ 0.5 & 0.1 & 0.2 & 0 \end{pmatrix}.$$

Here, all the arcs have the same length.

We start initially with density $\rho$ and momentum $q$ null on each arc except the first arc, where $\rho$ is a small perturbation of a constant state equal to 100 and the velocity is constant equal to 10. We let the system evolves until $T = 50$.

At time $T = 10$, we see the propagation of the density along the network : the system is nearly at equilibrium but we notice some larger values for the density on the arcs which are connected to the first arc, especially the "neighboring" one which is also connected to an outer node. The value for $q$ has also strongly decreased in the first arc and the arcs where the modulus of the momentum is the highest are the ones linked to node 11.

We notice that at time $T = 50$, the system has reached an asymptotic state with some null momentum $q$ on each arc of the network and a constant density $\rho$ which is the same on each arc of the network; the constant density can be computed thanks to mass conservation and the initial mass on the whole network, see Figure16.

<figure><img src="image_16.png" /><figcaption>FIGURE 15. Network 1 configuration with 12 arcs and 12 nodes</figcaption></figure>
<figure><img src="image_17.png" /></figure>

(a) initial condition for $\rho$

(b) initial condition for $q$

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

(c) Function $\rho$ at time $T = 10$

(d) Function $q$ at time $T = 10$

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

(e) Function $\rho$ at time $T = 50$

(f) Function $q$ at time $T = 50$

FIGURE 16. Network 1: on top, initial condition; in the middle, solution at time $T = 10$; on bottom, solution at time $T = 50$. On the left, function $\rho$; on the right, function $q$.

**Network 2:** We consider now a grid-shaped network with 26 arcs and 18 nodes, see [13]: 16 junction nodes and 2 outer nodes, see Figure17. The junction nodes are connecting to 2, 3 or 4 arcs, namely nodes inside the grid are connected to 4 nodes, nodes on the arcs of the grid to 3 nodes and nodes on the corners of the grids to 2 or 3 other nodes. Permeability coefficients depend on the number of arcs at the node, more precisely if we denote by $N_p$ the number of arcs at node $p$, then $\kappa_{i,j} = 1/N_p$ for all $i \neq j$ and 0 otherwise. The lengths of the arcs are also different and are equal to $L = 10$ except $L_1 = L_5 = L_9 = L_{10} = L_{14} = L_{21} = L_{25} = L_{26} = 0.5$.

The initial density $\rho$ is a small perturbation on each arc of a constant equal to 0.45 and the initial momentum $q$ is constant everywhere equal to 0.1.

At time $T = 10$, we notice that positive values for $q$ are concentrated around the nodes of the network, whereas null or positive values of the momentum take place on the arcs. Regarding the density, we remark the highest values on the arcs around one outer node of the network and the smallest values around the other outer node.

At time $T = 50$, the values for the density follow the same pattern as for $T = 10$, but the values are more homogeneous all around the network. As for the momentum, it is null on each arc and $q$ have positive and negative values around the nodes of the networks, where the main changes occur.

Finally at time $T = 150$, once again, the system converges in large time towards a constant state for $\rho$ and $q$, but the convergence time is much longer than for the previous test cases, since the network is bigger see Figure18. The value of the density is slightly higher on the arcs with smaller lengths, see the red path on the figure.

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

FIGURE 17. Network 2 configuration with 26 arcs and 18 nodes; some of the arcs are smaller
network4_rhoinit-eps-converted-to.pdf network4_qinit-eps-converted-to.pdf

(a) initial condition for $\rho$

(b) initial condition for $q$

network4_rho_T10-eps-converted-to.pdf network4_q_T10-eps-converted-to.pdf

(c) Function $\rho$ at time $T = 10$

(d) Function $q$ at time $T = 10$

network4_rho-eps-converted-to.pdf

network4_q-eps-converted-to.pdf

(e) Function $\rho$ at time $T = 50$

(f) Function $q$ at time $T = 50$

network4_rho_T150-eps-converted-to.pdf network4_q_T150-eps-converted-to.pdf

(g) Function $\rho$ at time $T = 150$

(h) Function $q$ at time $T = 150$

FIGURE 18. Network 2: on top, initial condition; in the second row, solution at time $T = 10$; in the third row, solution at time $T = 50$; on bottom, solution at time $T = 150$. On the left, function $\rho$; on the right, function $q$.

REFERENCES

[1] B.Andreianov, G.Coclite, C. Donadello, *Well-posedness for vanishing viscosity solutions of scalar conservation laws on a network*, Discrete and Continuous Dynamical Systems, Volume 37, Number 11, 5913-5942 (2017).

[2] D. Aregba-Driollet, *Convergence of lattice Boltzmann methods with overrelaxation for a nonlinear conservation law*, ESAIM, Math. Model. Numer. Anal. 58, No. 5, 1935-1958 (2024).

[3] D. Aregba-Driollet, R. Natalini, *Convergence of relaxation schemes for conservation laws*, Appl. Anal. 61, No. 1-2, 163-190 (1996).

[4] D. Aregba-Driollet, R. Natalini, *Discrete kinetic schemes for multidimensional systems of conservation laws*, SIAM J. Numer. Anal. 37, No. 6, 1973-2004 (2000).

[5] F. Bellamoli, L.Müller, E.Toro, *A numerical method for junctions in networks of shallow-water channels*, Appl. Math. Comput., 337, 190-213 (2018).

[6] S. Bianchini, *Hyperbolic limit of the Jin-Xin relaxation model*, Commun. Pure Appl. Math. 59, No. 5, 688-753 (2006).

[7] F. Bouchut, *Construction of BGK models with a family of kinetic entropies for a given system of conservation laws*, J. Stat. Phys. 95, No. 1-2, 113-170 (1999).

[8] F. Bouchut, Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources. Frontiers in Mathematics. Basel: Birkhäuser (ISBN 3-7643-6665-6/pbk). viii, 135 p. (2004).

[9] A. Bressan; S. Canic; M. Garavello; M. Herty; B. Piccoli, *Flows on networks: Recent results and perspectives*, EMS Surv. Math. Sci., 1, 47-111 (2014).

[10] E. Braun, G. Bretti, and R. Natalini, *Mass-preserving approximation of a chemotaxis multi-domain transmission model for microfluidic chips*, Mathematics, 9(6), 2021.

[11] G. Bretti, R. Natalini, and B. Piccoli, *Numerical approximations of a traffic flow model on networks*, Netw. Heterog. Media, 1(1):57-84, 2006.

[12] G. Bretti, R. Natalini, and M. Ribot, *A hyperbolic model of chemotaxis on a network: a numerical study*, ESAIM Math. Model. Numer. Anal., 48(1):231-258, 2014.

[13] G. Bretti, R. Natalini, *Numerical approximation of nonhomogeneous boundary conditions on networks for a hyperbolic system of chemotaxis modeling the Physarum dynamics*, J. Comput. Methods Sci. Eng., 18(1):85-115, 2018.

[14] M. Briani, B. Piccoli, *Fluvial to torrential phase transition in open canals*, Netw. Heterog. Media 13, No. 4, 663-690 (2018).
[15] M. Briani, G.Puppo, M.Ribot, Angle dependence in coupling conditions for shallow water equations at channel junctions, Comput. Math. Appl. 108, 49-65 (2022).

[16] F. Calabrò, P. Zunino, Analysis of parabolic problems on partitioned domains with nonlinear conditions at the interface: application to mass transfer through semi-permeable membranes, Mathematical models and methods in applied sciences 16, 4, 479-501 (2006).

[17] A. Cangiani, R. Natalini, A spatial model of cellular molecular tracking including active transport along microtubules, Journal of theoretical biology 267, 4, 614-625 (2010).

[18] G. Ciavolella, N. David, A. Poulain, Effective interface conditions for a porous medium type problem, Interfaces Free Bound. 26, No. 2, 161-188 (2024).

[19] G. Coclite, M. Garavello, Vanishing viscosity for traffic on networks, SIAM Journal on Mathematical Analisys, 42(4), 1761-1783 (2010).

[20] G. Coclite, C. Donadello, Vanishing viscosity on a star-shaped graph under general transmission conditions at the node, Netw. Heterog. Media, 15(2), 197-213 (2020).

[21] R. Colombo, M. Garavello, A well posed Riemann problem for the p-system at a junction, Netw. Heterog. Media 1, No. 3, 495-511 (2006).

[22] R. Colombo, M. Herty, V. Sachers, On 2×2 conservation laws at a junction, SIAM Journal on Mathematical Analysis, 40, 605-622 (2008)

[23] R. Dager and E. Zuazua, Wave Propagation, Observation and Control in 1-D Flexible Multi-Structures, Math. Appl. 50, Springer-Verlag, Berlin, 2006.

[24] C.Dafermos, Hyperbolic Conservation Laws in Continuum Physics, Springer, Heidelberg, 2016.

[25] M. Di Francesco, P. Marcati, Singular convergence to nonlinear diffusion waves for solutions to the Cauchy problem for the compressible Euler equations with damping, Math. Models Methods Appl. Sci. 12, No. 9, 1317-1336 (2002).

[26] A. Gamba, D. Ambrosi, A. Coniglio, A. de Candia, S. Di Talia, E. Giraudo, G. Serini, L. Preziosi, and F. Bussolino, Percolation, morphogenesis, and Burgers dynamics in blood vessels formation, Phys. Rev. Lett., 90(11), 118101, 2003.

[27] M. Garavello, A review of conservation laws on networks. Netw. Heterog. Media 5, No.3, 565-58 (2010).

[28] M. Garavello, K. Han, B. Piccoli, Models for vehicular traffic on networks. AIMS Series on Applied Mathematics, 9. American Institute of Mathematical Sciences (AIMS), Springfield, MO, 2016.

[29] M. Garavello, R. Natalini, B. Piccoli, A. Terracina, Conservation laws with discontinuous flux, Netw. Heterog. Media 2, No. 1, 159-179 (2007).

[30] M. Chaplain, C.Giverso, T.Lorenzi, L.Preziosi, Derivation and application of effective interface conditions for continuum mechanical models of cell invasion through thin membranes, SIAM Journal on Applied Mathematics, 79(5), 2011-2031 (2019).

[31] J. Glimm, P.D. Lax, Decay of solutions of systems of nonlinear hyperbolic conservation laws, Mem. Am. Math. Soc. 101, 112 p. (1970).

[32] F.R.Guarguaglini, R.Natalini, Global smooth solutions for a hyperbolic chemotaxis model on a network, SIAM J.Math.Anal. vol. 47, No. 6, 4652-4671 (2015).

[33] F.R.Guarguaglini, R.Natalini, Vanishing viscosity approximation for linear transport equations on finite star-shaped networks, J. Evol. Equ. 21, No. 2, 2413-2447 (2021).

[34] F.R. Guarguaglini, Stationary solutions and asymptotic behaviour for a chemotaxis hyperbolic model on a network, Netw. Heterog. Media, Vol. 13, No. 1, 47-67 (2018).

[35] F. R. Guarguaglini, Global solutions for a chemotaxis hyperbolic-parabolic system on network with nonhomogeneous boundary conditions, Commun. Pure Appl. Anal., vol. 19 n.2, 1057-1087 (2020).

[36] F.R. Guarguaglini, M. Papi, F. Smarrazzo, Local and global solutions for a hyperbolic-elliptic model of chemotaxis on a network, Math. Models Methods Appl. Sci., vol.29 n.8, 1465-1509 (2019).

[37] A. Harten, P. Lax, B. Van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Review 25 (1983), 35-61.

[38] M. Herty, J. Mohring, V. Sachers, A new model for gas flow in pipe networks, Math. Methods Appl. Sci. 33, No. 7, 845-855 (2010).

[39] Y. Holle, M. Herty, M. Westdickenberg, New coupling conditions for isentropic flow on networks, Netw. Heterog. Media 15, No. 4, 605-631 (2020).

[40] S. Jin, Z. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions, Comm. Pure Appl. Math., Vol. 48, 235-277 (1995).
[41] O. Kedem, A. Katchalsky, *Thermodynamic analysis of the permeability of biological membranes to nonelectrolytes*, Biochim. Biophys. Acta 27, 229-246 (1958).

[42] P. Lefloch, M. Thanh, *The Riemann problem for fluid flows in a nozzle with discontinuous cross-section*, Commun. Math. Sci., 1, 4, 763-797 (2003).

[43] A. Marigo, *Entropic solutions for irrigation networks*, SIAM Journal on Applied Mathematics, 70, 1711-1735 (2010).

[44] D. Mugnolo, Semigroup methods for evolution equations on networks. Understanding Complex Systems. Springer, Cham, 2014.

[45] R. Natalini, *A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws*, J. Differ. Equations 148, No. 2, 292-317 (1998).

[46] R. Natalini, M. Ribot, M. Twarogowska, *A well-balanced numerical scheme for a one-dimensional quasilinear hyperbolic model of chemotaxis*, Commun. Math. Sci. 12, No. 1, 13-39 (2014).

[47] R. Natalini, M. Ribot, M. Twarogowska, *A numerical comparison between degenerate parabolic and quasilinear hyperbolic models of cell movements under chemotaxis*, J. Sci. Comput. 63, No. 3, 654-677 (2015).

[48] A. Quarteroni, A. Veneziani, P. Zunino, *Mathematical and numerical modeling of solute dynamics in blood flow and arterial walls*, SIAM Journal on Numerical Analysis 39, 5, 1488-1511 (2002).

[49] A. Roggensack, *A kinetic scheme for the one-dimensional open channel flow equations with applications on networks*, Calcolo, 50, 4, 255-282 (2013).

[50] G. Serini, D. Ambrosi, E. Giraudo, A. Gamba, L. Preziosi, and F. Bussolino, *Modeling the early stages of vascular network assembly*, The EMBO J., 22(8), 1771-1779, April 2003.

[51] D. Serre, *Semi-linear and kinetic relaxation for systems of conservation laws. (Relaxations semi-linéaire et cinétique des systèmes de lois de conservation.)* Ann. Inst. Henri Poincaré, Anal. Non Linéaire 17, No. 2, 169-192 (2000).

[52] D. Serre, Systems of conservation laws 2. Geometric structures, oscillations, and initial-boundary value problems. Translated from the French by I. N. Sneddon. Cambridge: Cambridge University Press. xi, 269 p. (2000).

[53] J. Valein and E. Zuazua, *Stabilization of the wave equation on 1-D networks*, SIAM J. Control Optim., 48, 2771-2797 (2009).

ISTITUTO PER LE APPLICAZIONI DEL CALCOLO, CONSIGLIO NAZIONALE DELLE RICERCHE

Email address: maya.briani@cnr.it

ISTITUTO PER LE APPLICAZIONI DEL CALCOLO, CONSIGLIO NAZIONALE DELLE RICERCHE

Email address: roberto.natalini@cnr.it

INSTITUT DENIS POISSON, UNIVERSITÉ D'ORLÉANS, CNRS, UNIVERSITÉ DE TOURS & UNIVERSITÉ PARIS-SACLAY, MAIAGE, INRAE JOUY-EN-JOSAS

Email address: magali.ribot@univ-orleans.fr