Central schemes for networked scalar conservation laws

1.   Introduction

Research on mathematical models on networks understood as directed (one-dimensional) graphs has been successfully conducted over the last decades and we refer to the recent survey ^[9] for details and references. Such models have various applications, such as gas dynamics in pipelines ^[11,39,50], vehicular traffic on road networks ^[25,34], production systems ^[19] and blood flow through systems of blood vessels ^[23] to mention only a few. Apart from modeling questions regarding the partial differential equations on the edges, a major modeling challenge is the description of the dynamics at the network nodes, where adjacent edges connect. Starting with ^[34,39] networked (systems of) conservation or balance laws are defined by (physically induced) coupling conditions, see ^[9] for examples, and e.g, ^[35] for a hierarchical derivation in the case of gas dynamics. Those conditions yield under suitable assumptions boundary conditions. Using wave–front–tracking techniques well–posedness of such coupled problems could be established, ^{[8,16,17,25,34]}. A key analytical concept here is the notion of Riemann solvers ^[24,25] or half–Riemann problems ^[32]. Among others, those concepts require an analytic expression of wave curves of the corresponding models. Here, we are interested in the numerical treatment of networked conservation laws. Numerical methods based on the Riemann solvers have been proposed already in ^[2,26,34] and have recently gained interest in view of high–order methods ^{[1,6,7,10,13,51,54]}, property–preserving schemes ^[49,57] and also for problems where Lax–curves and eigenvalues are not explicit ^[42,52] or not available ^[30,31]. Regarding the development of efficient schemes, most of the higher–order schemes rely on the linearization of the coupling conditions such that the Lax–curves are obtained trivially, see e.g., ^[1,7] for more details. Linearization techniques have also been used to avoid the problem of the explicit computations of eigenvalues and Lax-curves in two–phase problems, see ^[3]. The question of property preserving numerical schemes across networks has been recently investigated in view of well-balanced networked equations ^[49] and entropy-preserving schemes ^[57]. Note that for simplicity we consider here only the case of conservation laws but the numerical schemes directly extend to the case of balance laws. Further, we focus here on finite-volume or discontinuous Galerkin schemes whereas approaches based on finite-element schemes, such as ^[21], require a different treatment of coupling conditions. Similarly, a construction of vanishing viscosity solutions, which has been addressed in the schemes studied in ^[41,55], avoids the use of Riemann problems at the expense of formally treating parabolic systems.

In this manuscript we develop a numerical scheme that does not require a Riemann solver at the junction. To this end, we embed the coupling problem of nonlinear equations in the coupling of linear relaxation systems following the approach developed by Jin-Xin ^[38]. Since for linear systems the Lax curves are multiples of the a priori computable constant eigenvectors of the flux matrix, the computations will be explicit. A similar idea has been used in ^[40] to approximate a scalar conservation law with discontinuous flux. Another related approach has been considered in ^[5], where coupling conditions have been derived from a kinetic relaxation model. Applying an implicit-explicit discretization in time leads further to an explicit scheme contrary to e.g., ^[42]. Following an approach based on hyperbolic relaxation will allow us to define a numerical scheme that does not rely on the solution to the Riemann solvers at the node of the network. It can be extended to higher-order which is also demonstrated here and existing coupling conditions can be embedded in this framework as shown in an example on traffic flow.

2.   Notation and Preliminary Discussion

A network is a directed graph consisting of edges and vertices or junctions. We restrict the discussion to scalar hyperbolic conservation laws posed on each edge and, due to the finite speed of propagation, to the problem at a single coupling node, that we assume at position $x = 0$. The dynamics on the adjacent edge $k$ reads

where the state variable $u^k$ is either given on an incoming edge parameterized by $\mathcal{E}_k = (-\infty, 0)$ if $k\in \delta^- = \{1, \dots, N^- \}$ or on an outgoing edge parameterized by $\mathcal{E}_k = (0, \infty)$ if $k\in\delta^+ = \{N^-+1, \dots, N^-+ N^+ = N\}$. Here, the flux functions are smooth, but not-necessarily convex or concave functions $f_1, \dots, f_{N}:$ $\mathbb{R} \rightarrow \mathbb{R}$. The set of all edges is denoted by $\delta^{\mp} = \delta^- \cup \delta^+$. In addition, we assume given initial data that we denote on each edge $k\in \delta^\mp$ by $u^{k, 0}$.

The coupling is described by a set of conditions on the traces at the coupling node of the form

assuming a mapping $\Psi: \mathbb{R}^{N} \rightarrow \mathbb{R}^\ell$. The number of coupling conditions $\ell$ required to obtain a well-posed problem depends on the number of edges and the choice of flux functions, see e.g., ^[9] and the discussion following below.

The typical (numerical) procedure to obtain conditions on the traces of the solution $u^k$ at $x = 0$ relies on a suitable Riemann solver at the junction ^[25]. Illustrated in the case $|\delta^+| = |\delta^-| = 1$, where Eq (2.1) can be rewritten without edge indices as

the idea is to connect the traces $u_0^-$ and $u_0^+$ on the left and right of the interface to the coupling data $u_L$ and $u_R$ at the interface. The coupling data are determined such that the coupling conditions $\Psi(u_L, u_R) = 0$ are satisfied and the states $u_L$ and $u_R$ can be connected to the traces $u_0^-$ on the left and $u_0^+$ on the right by means of Lax curves corresponding to characteristic fields with negative and positive characteristic speeds, respectively. This problem is formulated using parameterized Lax-curves and reduced to a typically nonlinear but finite-dimensional system of equations, see e.g., ^[9,29,33].

To elaborate on the well-posedness and derivation of coupling conditions we introduce the Riemann problem on the edges of the network, which reads

with initial data given either by

or otherwise by

On incoming edges the right initial data of the Riemann problem at time $t = 0$ denoted by $u^k_R$ is unknown whereas the left data $u^k_0$ is assumed known. Analogously, on outgoing edges the left initial data of (2.4) denoted by $u^k_L$ is unknown whereas the right data $u^{k}_0$ is given. Unknown data is obtained by a Riemann solver.

Definition 1 (Riemann solver for scalar networks). A Riemann solver for problem (2.4) is a mapping that assigns right initial data on incoming edges and left initial data on outgoing edges to given initial data on the respective opposite side, i.e.,

such that (a) waves of the solution to (2.4) have negative speed on incoming edges and positive speed on outgoing edges and (b) $\Psi[u^{1}_R, \dots, u^{N^-}_R, u^{N^-+1}_L, \dots, u^{N}_L] = 0$.

A necessary condition for conservation of $u$ at the coupling node is the Kirchhoff condition:

Hence, to conserve the quantity $u$ in the junction of the network, the coupling data obtained by the Riemann solver needs to satisfy the condition

In the scalar network (2.1) condition (2.5) together with admissible boundary data, see ^[20], lead to mass conservation at the coupling node. However, these conditions are not necessarily sufficient for well-posedness of the Riemann solver and the network problem, see ^[9].

3.   Coupled Relaxation System

We follow ^[38], where in addition to the scalar quantity $u$ and the flux function $f$ the auxiliary variable $v \in \mathbb{R}$, the relaxation rate $\varepsilon > 0$ and the relaxation speed $\lambda > 0$ have been introduced and the following system is studied:

The relaxation system is accompanied by initial data given through a scalar function $u^0$ and $v^0 = f(u^0)$. As $\varepsilon\rightarrow0$ the system attains the zero relaxation limit $(u, v) := (u^0, v^0)$ for which Eq (3.1b) necessitates the local equilibrium $\lim_{ \varepsilon\rightarrow 0} v^ \varepsilon = v = f(u)$ and Eq (3.1a) recovers the conservation law

The Chapman-Enskog expansion ^[14] allows for an interpretation of (3.1) as a dissipative equation and shows that $v^ \varepsilon = f(u^ \varepsilon) + O(\epsilon^2)$. The subcharacteristic condition

introduced in ^[48] guarantees that the dissipative approximation is well-posed. The dissipative equation and the conservation law were shown to have the same asymptotic behavior as the relaxation rate goes to zero in ^[15]. An eigenvalue analysis of the relaxation system, see Appendix A, reveals that it is hyperbolic, has the eigenvalues $-\lambda$ and $\lambda$ and can be rewritten in terms of its characteristic variables $w^{ \varepsilon\mp} = v^ \varepsilon/2 \mp \lambda u^ \varepsilon/2$ as

Moreover, the forward and backward Lax-curves of the relaxation system are given by straight lines in the phase plane as follows

3.1.   Relaxation System at 1-to-1 networks

Since the relaxation system describes the behavior of the conservation law in the relaxation limit we use it as a tool to derive suitable coupling data for the 1-to-1 coupling network. To this end we will analyze system (3.1) in the 1-to-1 coupling case in this section and derive admissible coupling data. The relaxation limit will be taken in Section 3.2 in case of an asymptotic preserving discretization of the system, which was shown to converge to the correct limit as the relaxation rate tends to zero, see ^[37] and the references therein.

We consider a 1-to-1 network, which couples two relaxation systems of the form (3.1) on a single incoming and a single outgoing edge at the coupling node and reads

While the scalar quantity $u^ \varepsilon$ is governed by the same equation left and right from the interface the scalar auxiliary variable $v^ \varepsilon$ determining the flux of $u^ \varepsilon$ is governed by the two equations (3.6b) and (3.6c), which account for different flux functions $f_1$ and $f_2$ left and right from the coupling node. Note that we choose the same relaxation rate but allow for different relaxation speeds $\lambda_1, \lambda_2 > 0$ left and right from the node. We assume that subcharacteristic conditions of the form (3.3) hold at both edges. Initial data is given through a smooth and compactly supported function $u^{ \varepsilon, 0}: \mathbb{R} \setminus \{0\}\rightarrow \mathbb{R}$ and the initial condition

for all $x\in \mathbb{R}\setminus \{ 0\}$ with $\chi$ denoting the characteristic function. To close the system coupling conditions of the form

taking into account both variables of the system are required. Since a linear system of two conservation laws is given, two conditions are imposed for well-posedness, i.e., $\Psi: \mathbb{R}^4\rightarrow \mathbb{R}^2$. At a fixed time $t > 0$ we denote the traces left and right from the coupling node by $u_0^{ \varepsilon-}$, $v_0^{ \varepsilon-}$, $u_0^{ \varepsilon+}$ and $u_0^{ \varepsilon-}$. Analogously to Definition 1 we formally define a Riemann solver for the component-wise Riemann problem on the 1-to-1 network assigning coupling data to the traces, i.e.,

In the following, the construction of (3.9) is discussed in detail. To obtain admissible boundary data, $(u_R^{ \varepsilon}, v_R^{ \varepsilon})$ needs to be connected to $(u_0^{ \varepsilon-}, v_0^{ \varepsilon-})$ by a wave with negative velocity, whereas $(u_L^{ \varepsilon}, u_L^{ \varepsilon})$ needs to be connected to $(u_0^{ \varepsilon+}, u_0^{ \varepsilon-})$ by a wave with positive velocity as shown in Figure 1. Thus, by the eigenvalue analysis of the relaxation system, we require the two conditions

Well-posedness of (3.9) is obtained taking into account the two coupling conditions, which must be satisfied by the coupling data obtained by the Riemann solver. As the first coupling condition we impose the Kirchhoff condition (2.5) in Eq (3.6a) and obtain

Similarly, to conserve the mass of the auxiliary variable at the coupling node in the relaxation limit, we impose as second coupling condition (2.5) in Eq (3.6b) and Eq (3.6c) and get the condition

Combining Eq (3.10) and (3.11) a regular linear system is obtained, which determines (3.9). In explicit form the coupling data is given by

Remark 1. Defining $\lambda = \max \{ \lambda_1, \lambda_2\}$ system Eq (3.6) can be rewritten using relaxation speed $\lambda$ in both Eq (3.6b) and Eq (3.6c). In this case we obtain the simplified coupling data

A drawback of this adjustment is the increase of numerical diffusion in the schemes discussed in Section 3.2.

3.2.   Central Scheme for 1-to-1 Networks

In this section we derive the central scheme for 1-to-1 networks of scalar conservation laws. Therefore we will start in Sections 3.2.1 and 3.2.2 from a semi-discretization of the coupled relaxation system making use of the derived coupling data from Section 3.1. A time discretization is introduced in Section 3.2.3 for which the relaxation limit is considered in Section 3.2.4. In Section 3.2.5 we introduce a second order scheme.

3.2.1.   Semi-discrete Scheme

We introduce a uniform grid on the real line by fixing $\Delta x > 0$ and defining the mesh points $x_{j-1/2} = j \, \Delta x$ for any $j \in \mathbb{Z}$. We denote the approximate average of any scalar quantity $q$ in the cell $I_j:= [x_{j}, x_{j+1}]$, which still depends on the time variable, by $q_j$.

We obtain a scheme for the relaxation system in characteristic variables (3.4) by applying the first order upwind discretization, see e.g., ^[46]. Hereby we obtain by the signs of the eigenvalues for any $j\in \mathbb{Z}$

To derive (3.14), we have additionally applied a midpoint discretization to the flux function to approximate $f(q_j) \approx \int_{I_j}f(q(t, x)) \, dx$. Transforming back to the original variables $u_ \varepsilon = (w^{ \varepsilon+} - w^{ \varepsilon-})/\lambda$ and $v_ \varepsilon = w^{ \varepsilon+} + w^{ \varepsilon-}$ we end up with a semi-discrete scheme for system (3.1), that reads for any $j\in \mathbb{Z}$

3.2.2.   Coupled scheme

We consider a discretization of the coupled relaxation system (3.6). Figure 2 shows an illustration of the setting and the space discretization. We denote by $\dots, u_{-2}^ \varepsilon, u_{-1}^ \varepsilon$ and $\dots, v_{-2}^ \varepsilon, v_{-1}^ \varepsilon$ discretizations of $u^{ \varepsilon}$ and $v^{ \varepsilon}$ left from the coupling node and by $u_0^ \varepsilon, u_1^ \varepsilon, \dots$ and $v_0^ \varepsilon, v_1^ \varepsilon, \dots$ discretizations of $u^{ \varepsilon}$ and $v^{ \varepsilon}$ right from the coupling node, respectively. Now we apply scheme (3.15) and use the coupling data derived in Section 3.1 as ghost cell averages beyond the coupling node when approaching $x = 0$ from the left and from the right.

Thus we obtain

for the evolution of the cell average left from the coupling node and

for the evolution of the cell average right from the coupling node. Clearly, (3.16) and (3.17) can be complemented to a scheme over the full real line by additionally considering (3.15) for $j\in \mathbb{Z} \setminus \{-1, 0 \}$ with $\lambda$ and $f$ substituted by $\lambda_1$ and $f_1$ for negative $j$ and $\lambda_2$ and $f_2$ for positive $j$.

In the discretized setting traces are obtained from the cell averages next to the coupling node. Thus we have $u_0^{ \varepsilon-} = u_{-1}^ \varepsilon$, $u_0^{ \varepsilon+} = u_{0}^ \varepsilon$, $v_0^{ \varepsilon-} = v_{-1}^ \varepsilon$ and $v_0^{ \varepsilon+} = v_{0}^ \varepsilon$. Substituting now the coupling data (3.12) into (3.16) and (3.17) we obtain in case of $\lambda = \lambda_1 = \lambda_2$ (cf. Remark 1) left from the coupling node

and right from the coupling node

The corresponding evolution formulas in case of different relaxation speeds $\lambda_1\neq \lambda_2$ are given in Appendix B. The following consistency result follows from (3.18) and (3.19).

Proposition 1 (Consistency of the semi-discrete scheme). We assume $f = f_1 = f_2$ as well as $\lambda_1 = \lambda_2 = \lambda$. Then the semi-discrete scheme for the 1-to-1 network system (3.6), which consists of (3.15) for $j \in \mathbb{Z} \setminus \{-1, 0\}$ supplemented by (3.16), (3.17) and coupling data (3.12), is identical to the semi-discrete scheme for the uncoupled relaxation system (3.1) given by (3.15) for $j \in \mathbb{Z}$.

3.2.3.   Fully discrete scheme

In this section we derive an implicit–explicit scheme and consider to this end a uniform partition of the time line by introducing the time increment $\Delta t > 0$ and setting $t^n = n \Delta t$ for all $n \in \mathbb{N}_0$. The approximate average in the cell $I_j$ of any scalar quantity $q$ at time $t^n$ is denoted by $q_j^n$. For $j \in \mathbb{Z} \setminus \{ -1, 0\}$ an implicit-explicit time discretization of system (3.15) is given by

where $\alpha(j) = 1$ for negative $j$ and $\alpha(j) = 2$ for positive $j$. This scheme was proposed in ^[36] and while it handles most terms explicitly, the stiff relaxation term is treated implicitly for increased stability. Since at each solution update $u_j^{ \varepsilon, n+1}$ can be computed first by the explicit evolution formula (3.20a), it is not necessary to solve a nonlinear system to evaluate Eq (3.20c) afterwards.

We complement this scheme left from the coupling node by

and right from the coupling node by

As (3.21) is a time discretization of (3.16) and (3.17) the coupling data accounting for boundary information in (3.21) is derived from cell averages next to the coupling node by the Riemann solver (3.9) as

Remark 2. In ^[36] the authors further consider an alternative two-stage time discretization, which treats the relaxation implicitly in a similar fashion as scheme (3.20). The limit scheme derived in Section 3.2.4 is also obtained when this alternative time discretization is used.

3.2.4.   Relaxation Limit

In this section we derive the relaxation limit of the scheme for the 1-to-1 relaxation system (3.6) introduced in Section 3.2.3 given by (3.20) and (3.21). Our aim is to obtain in this way a scheme for scalar conservation laws in the 1-to-1 network case (2.3).

In the uncoupled case, scheme Eq (3.20) is asymptotic preserving, see ^[37], and thus attains a suitable limit scheme for the conservation law (3.2) as $\varepsilon$ tends to zero (with discretization parameters fixed). We derive the limit scheme in the coupled case. For the limit process an asymptotic expansion at the relaxation state of the state variables can be considered. As in the scheme the relaxation time appears only in the discretized balance term, the following procedure is equivalent: we keep $\Delta x$ and $\Delta t$ fixed and assume that the magnitude of these quantities as well as of the occurring cell averages are independent of the relaxation time and then consider the limit $\varepsilon\rightarrow 0$ ^[36]. From Eq (3.20c), Eq (3.21b) and Eq (3.21d) we get

for any $n\in \mathbb{N}_0$, where we have used the limit notations $u_j^{n+1} = \lim_{ \varepsilon\rightarrow 0} u_j^{ \varepsilon, n+1}$ and $v_j^n = \lim_{ \varepsilon\rightarrow 0} v_j^{ \varepsilon, n+1}$. Consequently, the evolution formulas for the auxiliary variable can be discarded and we obtain the limit scheme

for $j\in \mathbb{Z} \setminus\{-1, 0\}$ in Eq (3.24a) and with $\alpha(j) = 1$ for negative $j$ and $\alpha(j) = 2$ for positive $j$. We emphasize that in general $v^n_R \neq f_1(u^n_R)$ and $v^n_L\neq f_2(u^n_L)$ in Eq (3.24b) and Eq (3.24c) and that, although the scheme approximates the single scalar quantity $u$, coupling data of the auxiliary variable is required for its evaluation. We further note that the limit scheme is fully explicit. Coupling data in Eq (3.24) is determined by the Riemann solver (3.9) as

If we assume $\lambda = \lambda_1 = \lambda_2$ and substitute the coupling data (3.12) derived in Section 3.1 taking into account (3.25), we obtain

which replace Eq (3.24b) and Eq (3.24c) in the scheme (3.24). This makes evident that the limit scheme can be written in the conservative form

using numerical fluxes that depend on the two cell averages next to the cell interface $x_j$ and read

We provide evolution formulas corresponding to (3.26) and (3.27) for differing relaxation speeds in Appendix B. In case $f_1 = f_2$ scheme (3.27) is the first order relaxed scheme from ^[38], which is identical to the Lax–Friedrichs scheme if we additionally assume $\lambda = \frac{\Delta x}{\Delta t}$. Thus the consistency property from Proposition 1 transfers to the limit scheme as follows.

Proposition 2 (Consistency of the limit scheme). Let $f_1 = f_2 = f$ and $\lambda_1 = \lambda_2 = \lambda$. Then the limit scheme (3.24) using coupling data (3.25), (3.12) is identical to the first order relaxed scheme (see ^[38]) for the scalar conservation law.

The relaxed scheme for the uncoupled conservation law was analyzed regarding monotonicity and the result carries over directly to the coupled scheme in the relaxation limit.

Proposition 3 (^[18,38]). Let $f_1 = f_2 = f$, $\lambda_1 = \lambda_2 = \lambda$, $\lambda \Delta t \leq \Delta x$ and assume that the subcharacteristic condition (3.3) holds. Then the limit scheme (3.24) using coupling data (3.25), (3.12) is a consistent and monotone scheme (see e.g., ^[27]) and therefore converges to the entropy solution of (3.2) as both $\Delta t$ and $\Delta x$ tend to zero.

Due to (3.27) and its generalization in Appendix B we state the following result about the conservation of the total mass $\sum_{j\in \mathbb{Z}} u_j^n$ over time. Conservation at the coupling node for general networks is defined in Section 4.3, whereas it follows directly from the conservative form in the 1-to-1 case.

Corollary 1. The limit scheme for the 1-to-1 network (2.3) given by (3.24) and coupling data (3.25), (3.12) is conservative at the coupling node. It is further conservative over the full real line in the case $f_1(0) = f_2(0)$ and over bounded domains if zero-flux boundary conditions are imposed.

We moreover note that scheme Eq (3.27) is in general not positivity-preserving in case of non-negative fluxes and initial data. Other schemes for conservation laws in the 1-to-1 network case have been introduced in the literature. In particular, the work ^[26] analyzes a similar scheme, which can be written as scheme (3.24) with coupling data $u_R^n = u_0^n$, $v_R^n = f_1(u_0^n)$, $u_L^n = u_{-1}^n$ and $v_R^n = f_2(u_{-1}^n)$. Unlike the method we introduce here, this scheme does not admit a conservative form but it has been proven to converge.

3.2.5.   Second-order Extension

In this section we are concerned with a second order scheme for the 1-to-1 scalar network. The scheme is derived by extending the approach in Sections 3.2.1–3.2.4 using piece-wise linear approximations left and right from the coupling node. A similar extension is considered in ^[38] for uncoupled conservation laws.

A second order scheme for system (3.1) is derived by applying the MUSCL scheme ^[56] to the upwind discretization of the relaxation system in characteristic variables (3.14). This way we obtain the semi-discrete scheme

for all $j\in \mathbb{Z}$, which makes use of interface reconstructions. These are obtained by extrapolation from the upwind direction, i.e.,

The slope in each cell is given by the monotonized central-difference limiter

The selected limiter (3.30) was introduced in ^[56] and designed to yield sharp resolutions near discontinuities. In smooth regions it admits the central difference whereas the accuracy at non-sonic critical points is reduced to preserve the monotonicity of the discrete solution. The minmod operator used in its formulation is defined by

As in the derivation of scheme (3.15) we transform (3.28) back to the original variables of the relaxation system, and get

The resulting scheme adds second order extension terms to its first order version, to which we refer in the following by $\mathcal{S}^u_j$ and $\mathcal{S}^v_j$. The included limited slopes (3.30) can be alternatively computed by applying the minmod operator (3.31) to the terms

Analogously to Section 3.2.2, we apply scheme (3.32) to the coupled relaxation system (3.6), see also Figure 2. In the coupled scheme (compare (3.16) and (3.17)) the following second order extensions appear in the evolution formulas of the volumes next to the coupling node

Our coupling approach gives rise to Dirichlet boundary problems on the edges of the network and does not provide any information about the slope beyond the coupling node. Being derived from the Riemann problem coupling data is assumed spatially constant and thus, it is natural to set

We note that by setting $s_j^\mp = 0$ for any $j\in \mathbb{Z}$ in the MUSCL scheme we locally recover the first order scheme. Higher order approximations beyond the coupling node are not considered in this work but have been achieved by transforming spatial to temporal information using an ADER approach as e.g., in ^[1,7]. Moreover, due to (3.5) the coupling data (3.12) satisfies

In fact Eq (3.36) is an equivalent formulation to Eq (3.10) in characteristic variables. Consequently, when using coupling data $u_R^{ \varepsilon}$, $v_R^{ \varepsilon}$, $u_L^{ \varepsilon}$ and $v_L^{ \varepsilon}$ as ghost cell data at the coupling node when computing the linear reconstructions, we get due to Eq (3.30) and Eq (3.31)

Following the steps in Sections 3.2.3 and 3.2.4 we discretize in time and take the relaxation limit. We thus obtain the second order limit scheme

for $j\in \mathbb{Z} \setminus\{-1, 0\}$ in Eq (3.38a) and with $\alpha(j) = 1$ for negative $j$ and $\alpha(j) = 2$ for positive $j$. The time discrete slopes $s^{n\mp}_j$ in (3.38) are obtained by applying the limiter (3.30) to the time discrete characteristic variables

By supplementing the coupling data (3.25), (3.12) we again obtain the conservative form (3.27a) using the numerical fluxes

At the coupling node the numerical flux is identical to the one of the first order scheme given by Eq (3.27a) and Eq (B.4). In the scheme coupling data (3.25), (3.12) is required to compute $s_{-1}^{n-}$ and $s_0^{n+}$ according to Eq (3.30), which takes into account $w^{n-}_{R} := v_R^n/2 - \lambda_1 u_R^n/2$ and $w^{n+}_{L} := v_L^n/2 + \lambda_1 u_L^n/2$. Alternatively, these slopes can be set to zero for simplicity.

Similar to Proposition 2 consistency to a second order scheme for Eq (3.2) on the real line is given if $f_1 = f_2$, $\lambda_1 = \lambda_2$ and when neglecting the numerical flux at the coupling node. For this reason the following result is deduced.

Proposition 4. Let $f_1 = f_2 = f$, $\lambda_1 = \lambda_2 = \lambda$, $\lambda \Delta t \leq \Delta x/2$ and assume that the subcharacteristic condition (3.3) holds. Then the scheme (3.38) using coupling data (3.25), (3.12) and $s_{-1}^{n-} = s_0^{n+} = 0$ for all $n\in \mathbb{N}_0$ is consistent, $\textrm{total variation diminishing}$ (see e.g., ^[27]) and converges to the weak solution of the conservation law (3.2).

Proof. From a result in ^[38] follows that the scheme is total variation diminishing. In this work the authors considered a scheme of the form (3.38a) (with uniform $\lambda$ and $f$) for $j\in \mathbb{Z}$ but assumed different slopes. Yet, as the proof only requires that the slopes can be written as

which is satisfied by (3.30), (3.31) and

the result transfers to our second order scheme. We neglected time and sign indices above for clarity. Since the numerical fluxes are further consistent, the convergence result follows due to ^[46].

4.   Multiple Incoming and Outgoing Edges

In this section, we generalize the schemes introduced in Sections 3.2.4 and 3.2.5 to scalar conservation laws at nodes with $N^-$ incoming and $N^+$ outgoing edges. To this end we consider the relaxation system introduced in Section 3 in a network setting to derive admissible coupling data. In Section 4.2 we then discuss additional conditions about the flux distribution at the junction, which are required for well posedness of coupling data. In Section 4.3 we eventually present the central schemes for the network problem (2.1).

4.1.   The Relaxation System on Networks

We generalize the approach from Section 3.1 and consider the relaxation system (3.1) on a network that allows for $N^-$ incoming and $N^+$ outgoing edges. Again, we aim to take the relaxation limit in a space discretization of the network system in order to derive a scheme for the scalar network (2.1). By coupling multiple systems of the form (3.1) in the same way we coupled scalar conservation laws in (2.1) we obtain the system

governing the scalar variables $u^k$ and $v^k$, where $k\in\delta^\mp = \delta^- \cup \delta^+ = \{1, \dots, N\}$ indicates the corresponding edge of the network. The system includes a stiff source term that contains the flux functions $f_1, \dots, f_{N}:$ $\mathbb{R} \rightarrow \mathbb{R}$ and the relaxation rate $\varepsilon$ that is chosen uniform over the edges of the network. Compared to the 1-to-1 network (3.6) we do not account for $\varepsilon$ in the variable names here for brevity of presentation. On each edge of the network a relaxation speed $\lambda_k > 0$ is given, which satisfies the subcharacteristic condition, i.e.,

The smooth and compactly supported scalar functions $u^{k, 0}$ determine the initial condition of system (4.1) by

Furthermore, coupling conditions for the system variables at the coupling node have the form

Since the system is linear, $N$ conditions are necessary to obtain a well-posed problem and thus $\Psi: \mathbb{R}^{2N}\rightarrow \mathbb{R}^N$ (compare system (3.6), which can be written as (4.1) with $N = 2$ and requires the two conditions given in (3.11)). At a fixed time $t > 0$ we assume given traces $u_0^{1}, v_0^{1}, \dots, u_0^{N}, v_0^{N}$ at $x = 0$. Unlike in Section 3.1 we do not indicate if a trace is incoming or outgoing by a sign index as this can be seen from the edge index. We formally define a Riemann solver for system (4.1) by

The construction of (4.5) is discussed in the remainder of this section and in Section 4.2. We aim to derive admissible coupling data that verifies the conservation of the system variables in the coupling node. To ensure admissible boundary data, $(u_R^k, v_R^k)$ needs to connect to $(u_0^{k}, v_0^{k})$ by a wave with negative velocity for all $k\in \delta^-$ and $(u_L^{k}, u_L^{k})$ needs to connect to $(u_0^{k}, u_0^{k})$ by a wave with positive velocity for all $k\in \delta^+$. Thus, by the analysis in Appendix A, we get the conditions

Since all Lax-curves are parameterized by a single parameter, Eq (4.6) gives rise to $N$ unknowns, which are to be determined by the coupling conditions. As done in (3.11) for the 1-to-1 network we impose the Kirchhoff condition (2.5) for both system variables and obtain the two conditions

4.2.   Additional Conditions in the General Case

The two conditions given by (4.7) are not sufficient to obtain a well-defined Riemann solver for the network system (4.1). Therefore $N-2$ suitable additional conditions need to be imposed. Here we consider algebraic conditions on the auxiliary variable at the coupling node of the form

with parameters $\beta_k^\ell$, $r^\ell\in \mathbb{R}$ for $k\in\delta^\pm$ and $\ell = 1, \dots, N-2$. In the relaxation limit, the variable $v$ is the flux of the state variable of the conservation law in the relaxation limit and Eq (4.8) can be understood as conditions on the fluxes at the coupling node when the limit scheme is considered. These conditions are motivated by the network (2.1) and in its general form allow for common coupling conditions, e.g., from traffic modeling ^[25]. Motivated by the dissipative interpretation of the relaxation system, parabolic coupling conditions, e.g., ^[22,44], can be alternatively considered. These conditions impose continuity of the state variable at the coupling node in addition to the Kirchhoff condition and require suitable relaxation speeds to satisfy Eq (4.7b).

Taking into account Eq (4.6) we denote by $\sigma_k$ the position of $(u_R^k, v_R^k)$ for $k\in\delta^-$ or of $(u_L^k, v_L^k)$ for $k \in \delta^-$ on the corresponding Lax curve given by (3.5). Then we get from (4.7) and (4.8) the linear system

where $\nu_k$ denotes the sign of the edge defined as

When assigning parameters for Eq (4.8) care must be taken that the system (4.9) has full rank. For example in case of a 2-to-1 network and $\lambda_1 = \lambda_2 = \lambda_3$ this is satisfied whenever $\beta_1^1 \neq \beta_2^1$. We introduce the following vector notations for traces and coupling data

Moreover, we introduce the $N \times N$ diagonal matrices $\Lambda$ and $\mathcal{N}$ with diagonal entries $\lambda_1, \dots, \lambda_N$ and $\nu_1, \dots, \nu_N$, respectively. This allows us to express the coupling data by means of the linear systems

which hence define Eq (4.5). Here $\mathcal{A}$ and ${\bf b}$ denote the system matrix and the right-hand side of (4.9), i.e.,

In the rest of this section we give general suggestions for the parameters in Eq (4.8). However, we note that they are an issue of modeling and can be freely adapted to the problem at hand. We will assume non-negativity of the flux functions $f_1, \dots, f_N$ and the trace data $v_0^1, \dots, v_0^N$. In the relaxation limit, the latter follows directly from the non-negativity of the fluxes.

4.2.1.   Incoming edges

In case of multiple incoming edges a reasonable assumption motivated from free flow cases is that the relation between the incoming flux from a given edge to the total incoming flux in the coupling node is inherited from the traces. Since the variable $v$ determines the flux of the primary variable $u$ this is stated as

provided that both $\sum_{k\in\delta^-} v_0^k \neq 0$ and $\sum_{k\in\delta^-} v_R^k \neq 0$. This can be rewritten as the following conditions on the coupling data

It is sufficient to impose Eq (4.14) for all but one incoming edges, then still, by a summation argument, it follows for all incoming edges. Thus to account for assumption (4.14), the parameters for $\ell = 1, \dots, N^–1$ in Eq (4.9c) can be chosen as $\beta_\ell^\ell = \sum_{k\in\delta^-\setminus\{\ell\}} v_0^k$, $\beta_k^\ell = -v_0^m$ for $k \in \delta^-\setminus\{\ell\}$, $\beta_k^\ell = 0$ for $k \in \delta^+$ and $r^\ell = 0$. These trace dependent parameters lead to a nonlinear Riemann solver (4.5), whose evaluation only requires the assembly of (4.13) and the solution of the systems (4.12). In case $\sum_{k\in\delta^-} v_0^k \neq 0$ the system (4.9) (ignoring $\ell\geq N^-$) has maximal rank. To handle the case $\sum_{k\in\delta^-} v_0^k = 0$, we can regularize Eq (4.14) by adding a small constant $\epsilon > 0$ to both denominators. To account for this regularization in the parameters, we add $\epsilon$ to $\beta_\ell^\ell$ and set $r^\ell = \epsilon v^\ell$ for $\ell = 1, \dots, N^–1$. This way the Riemann solver can only assign nonzero incoming fluxes in the edge $k = N^-$ if all corresponding incoming traces are zero.

4.2.2.   Outgoing edges

The distribution of the fluxes among the outgoing edges in a network is often described by a distribution matrix, see e.g., ^[25]. This matrix is comprised of rates that describe how the flux from any incoming edge is distributed to the outgoing edges. Similarly, we impose by

that the flux that enters the coupling node from edge $k\in\delta^-$ is distributed to edge $m\in \delta^+$ with rate $\alpha^m_k$. The rates are chosen such that $0 \leq \alpha^m_k\leq 1$ and

Similar to condition (4.14) in case of incoming edges, condition (4.16) only needs to be imposed to all but one outgoing edges, then it also follows for the remaining outgoing edge due to Eq (4.7a) and Eq (4.17). Thus, if we set the parameters in Eq (4.9) for all $\ell = N^-, \dots, N-2$ to $\beta_k^\ell = \alpha^\ell_k$ for $k \in \delta^-$, $\beta^\ell_{\ell+1} = - 1$, $\beta_k^\ell = 0$ for $k\in\delta^+\setminus\{\ell+1\}$ and $r^\ell = 0$, condition Eq (4.16) holds for $m \in \delta^+$. When we further use the suggested parameters accounting for the incoming edges as discussed above, system (4.9) has full rank and the Riemann solver (4.5) is well-defined.

4.3.   Central Scheme for Networks of Scalar Conservation Laws

A scheme for the scalar network (2.1) is obtained by a straightforward generalization of the steps in Sections 3.2.1–3.2.5 and taking into account the generalized Riemann solver at the coupling node introduced in Sections 4.1 and 4.2. Unlike in Section 3.2 we do not substitute explicit formulas of coupling data into the scheme in the general network case as was done in (3.26). Instead each edge in the derivation is considered separately.

Hence, after discretizing system Eq (4.1) in space by scheme (3.15), discretizing in time as in (3.20) and taking the relaxation limit we obtain a network generalization of scheme (3.24) that can be written in conservative form as

where $\mathcal{I}_k = \mathbb{Z}^-$ if $k\in \delta^-$ or $\mathcal{I}_k = \mathbb{Z}^+_0$ if $k\in \delta^+$. The numerical fluxes are given by

If the included high order extension terms ${\mathcal{S}}^{k, n}_{j-1/2}$ are set to zero, the first order scheme is obtained. To get the second order scheme, as derived in Section 3.2.5 for the 1-to-1 network, we set

where the slopes of the linear reconstructions are given by

$w^{k, n\mp}_{j} = \frac 1 2 f_k(u^{k, n}_j) \mp \frac {\lambda_k}{2} u^{k, n}_j$ and the minmod operator (3.31). Coupling data is obtained by applying the Riemann solver (4.5) and setting

The form (4.18a) does not imply mass conservation at the coupling node, which we instead analyze by considering the incoming and outgoing numerical fluxes and the coupling data. Taking into account (3.5) in (4.18b) for $j = 0$, it follows that the numerical fluxes (4.18) coincide with the fluxes obtained by the Riemann solver, i.e.,

Thus the following result is a consequence of Eq (4.7a) and Eq (4.22).

Proposition 5 (Conservation property of the central scheme). The central scheme for the scalar network (2.1) given by (4.18) and coupling data (4.21), (4.5) is conservative in the coupling node, i.e., it holds

It is further conservative on the full network in the case $f_1(0) = f_2(0) = \dots = f_N(0) = 0$ and on bounded networks, where all edges that connect to the node are bounded and zero-flux boundary conditions are imposed.

5.   Numerical Experiments

In this section we apply the derived schemes in various numerical experiments to demonstrate their capabilities and performance. We focus on the limit schemes for system (2.1) and do not take into account discretizations of the relaxation system. The experiments consider scalar conservation laws on 1-to-1 and 2-to-1 networks.

We assume that the network nodes are bounded and parameterized by $(-1, 0)$ if they are incoming or by $(0, 1)$ if they are outgoing. Each edge is discretized over $m$ uniform mesh cells of size $\Delta x = 1/m$. We use a fixed relaxation speed $\lambda$ over the full network and take time increments as

Details on the used number of mesh cells, Courant number and boundary conditions are provided in the individual experiment descriptions. The employed computer programs are implemented in the Julia programming language ^[4] and publicly available from ^[43].

5.1.   Inviscid Burgers' equation

In the first numerical experiment we consider the inviscid Burgers' equation on a 1-to-1 network. To test the accuracy of the schemes we impose identical fluxes on both edges of the network given by (2.3) and $f_1(u) = f_2(u) = 1/2 \, u^2$. We consider an experiment adapted from ^[45] and employ the initial data

and periodic boundary conditions. While the solution of this problem is smooth at small times it later develops a shock discontinuity. This can be seen in Figure 3, which shows the numerical solutions computed by the first order scheme (3.27) and the MUSCL scheme (3.27a), (3.40) in the two top rows. In the computations $m = 50$ mesh cells were used left and right from the coupling node, the relaxation speed was assumed $\lambda = 1$ and Courant numbers were chosen $0.9$ in the first order and $0.2$ in the MUSCL scheme. The computations show that the higher order scheme achieves a sharper resolution of the developing shock. Furthermore, the reduction in accuracy in the coupling node ($x = 0$), which this scheme experiences, is visible at time instance $t = 0.2$.

To further investigate the impact of the coupling node we computed the $L^1$ and $L^\infty$ errors after mesh refinement as well as the corresponding experimental order of convergence (EOC). Besides the first order scheme (central scheme) we distinguish three variants of the MUSCL scheme: the first (central MUSCL) makes use of the coupling data in the computation of the slopes $s_{-1}^{n-}$ and $s_0^{n+}$, see Section 3.2.5, the second (central MUSCL TVD) sets these slopes to zero and is TVD according to Proposition 4 and the third (uncoupled MUSCL) applies the scheme to the uncoupled case, where the inviscid Burgers' equation is considered on the domain $(-1, 1)$ discretized over $2m$ cells. Solutions of all four schemes are presented in Figure 3. To see how the error depends on $\Delta x$ we computed different mesh solutions with time increments chosen according to Eq (5.1) and $\text{CFL} = 0.49$ for the first order scheme and fixed $\Delta t = 2 \times 10^{-6}$ for the MUSCL schemes. We computed the errors at time $t = 0.5$ when the problem still admits a smooth solution.

The computed $L^1$ errors shown in Table 1 confirm the expected orders of convergence, i.e., first order in the central scheme and second order in the MUSCL variants. The presence of the coupling node in the schemes central MUSCL and central MUSCL TVD only slightly reduces the accuracy when compared to the uncoupled MUSCL scheme and does not interfere in the convergence order in $L^1$. Thereby the central MUSCL scheme, which was also TVD in all numerical tests, yields a higher accuracy than the central MUSCL TVD scheme. The errors in $L^\infty$ shown in Table 2 behave differently. While the central scheme and the uncoupled MUSCL scheme yield the expected first and second experimental order, respectively, the experimental orders of the schemes central MUSCL and central MUSCL TVD are significantly reduced due to the handling of the coupling node. Still, these two schemes achieve high accuracy in $L^\infty$ similar to the uncoupled scheme.

5.2.   Traffic Flow

The second numerical experiment is concerned with a traffic scenario and imposes the Lighthill-Whitham-Richards (LWR) model ^[47,53]

on the edges of a 2-to-1 network. In more details, we consider the scalar problem Eq (2.1) with two incoming edges ($N^- = 2$) and one outgoing edge ($N^+ = 1$). We assume that the outgoing edge has larger capacity than the incoming ones and allows for higher traffic densities. This is reflected in the flux functions, which we set

Coupling conditions for traffic models have been a topic of high interest, see e.g., ^[25,28]. We are interested how the common flow maximization approach for coupling of (5.3) on networks compares to the coupling model implied by the presented scheme.

We sketch how coupling data is obtained according to flow maximization on 2-to-1 networks. Hereby we focus on the computation of the fluxes $v^1_R = f_1(u^1_R)$, $v^2_R = f_2(u^2_R)$ and $v^3_L = f_3(u^3_L)$. To constitute admissible boundary data these fluxes need to satisfy the demand and supply conditions, which impose upper bounds of the form

for details see ^[25]. Coupling fluxes are then chosen maximal under constraints given by Eq (5.4) and the Kirchhoff conditon (4.7a). Two cases can occur. In the so called free flow case, where $d_1(u_0^1) + d_2(u_0^2) \leq s_3(u_0^3)$, we take $v^1_R = d_1(u_0^1)$, $v^2_R = d_2(u_0^2)$ and $v^3_L = d_1(u_0^1) + d_2(u_0^2)$. The complimentary case is referred to as congestion. Here we set $v^3_L = s_3(u_0^3)$ and employ the right of way parameter $0\leq \beta \leq 1$ to set $v_R^1 = \beta s_3(u_0^3)$ and $v_R^2 = (1-\beta) s_3(u_0^3)$. If this leads to a violation of Eq (5.4) for either $v^1_R$ or $v^2_R$ the affected coupling flux is chosen as the respective upper bound and its counterpart is computed from the Kirchhoff condition. This procedure defines a Riemann solver, which can be applied to cell averages of a numerical scheme. To account for flow maximization in numerical simulations we employed scheme (4.18) and replaced the fluxes at the coupling node by

where the coupling fluxes were computed from the trace data $u_0^{1, n}$, $u_0^{2, n}$ and $u_0^{3, n}$ as above. In numerical experiments this choice of nodal flux has required smaller Courant numbers to prevent oscillations compared to the central scheme in Section 5.1. Therefore we fixed $\text{CFL} = 0.2$ for all schemes in our numerical comparison.

Coupling data in the central scheme does not require the parameter $\beta$ and was obtained from the linear systems (4.12). The parameters in system (4.9) were chosen so that Eq (4.14) was satisfied, i.e., $\beta^1_1 = v^2_0$, $\beta^1_2 = -v_0^1$ and $\beta^1_3 = r^1 = 0$. In the numerical computations for both coupling models we used $m = 200$ cells on each edge and assumed $\lambda = 1$. We further imposed zero-flux boundary conditions on the incoming edges and homogeneous Neumann boundary conditions on the outgoing edge.

To investigate the case of free flow, we first consider a numerical experiment with constant initial data on the edges chosen as $u^{0, 1} = 0.07$, $u^{0, 2} = 0.15$ and $u^{0, 3} = 0.2$. In both coupling models the traffic freely propagates from the incoming edges to the outgoing edge and eventually out of the network. The numerical solutions computed by the first order central scheme are presented in Figure 4 and show that the central scheme reproduces the dynamics of flow maximization in this case.

In a second numerical experiment we consider another Riemann problem using the modified initial data $u^{0, 1} = 0.6$, $u^{0, 2} = 0.35$ and $u^{0, 3} = 0.35$. The larger traffic densities lead to congestion at the coupling node. The numerical results by the first order scheme depicted in Figure 5 exhibit backward moving waves in the first incoming edge for both the flow maximization and the central approach. The behavior of the solution in the flow maximization case depends on the right of way parameter $\beta$. While smaller $\beta$ (solution for $\beta = 0.2$ shown) lead to backward moving waves only in the first incoming edge, larger $\beta$ (solution for $\beta = 0.5$ shown) lead to waves in both incoming edges. The solution of the central approach on the incoming edges is qualitatively similar to the case of flow maximization and a small right of way parameter. On the outgoing edge the central approach introduces a layer next to the coupling node connecting the sum of both incoming traffic densities near the coupling node to a decaying profile that is also obtained by the flow maximization approach (the numerical solution has been independent of $\beta$ here), see Figure 6. This might be a consequence of condition (4.7b). Computations on various meshes reveal that the layer is mesh dependent and decreases as the mesh is refined. The layer in the outgoing edge has not occured in case of free flow.

Moreover, Figure 7 shows that the central MUSCL scheme recovers the same dynamics as the first order central scheme but yields a smaller layer at the coupling node and higher resolution of the discontinuities. For an insight in the role of the relaxation speed we present analogue solutions with larger $\lambda = 10$ and $\lambda = 50$ in the same figure. On the incoming edge the solutions exhibit qualitatively similar behaviour compared to the solution assuming $\lambda = 1$ except for an increase in viscosity. On the outgoing edge the layer close to the coupling node increases in width with the relaxation speed and the density profile falls to a slightly lower baseline on the left part of the edge.

5.3.   Coupled Two-Phase Flow Model

Lastly, we apply our approach to the Buckley–Leverett equation ^[12], a simple model of two-phase flow. Given a mixture of water and oil in a tube of porous media, the water fraction $u$ in the model is governed by a scalar conservation law with the non-convex flux function

Again, we consider a 2-to-1 network and impose the model on its edges by taking the flux functions $f_1 = f_2 = f_3 = f$. In a numerical experiment we reproduce a scenario, in which water is pumped into two tubes in order to displace oil and enforce its outflow through a third tube. To this end we use the initial data

and homogeneous Neumann boundary conditions at the edges. The numerical solution computed by the central MUSCL scheme employing $m = 300$ cells per edge, $\text{CFL = 0.49}$ and $\lambda = 2.5$ is shown in Figure 8. In the first incoming edge a shock wave is formed that is followed by a rarefaction wave and passes through the coupling node to the outgoing edge, where it interacts with a second shock wave originating from the second incoming edge. We emphasize that our numerical approach resolved these network dynamics without analysis of the underlying complex (due to the non-convexity of the flux function) Riemann problem.

Acknowledgment

The authors thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for the financial support through 320021702/GRK2326, 333849990/IRTG-2379, CRC1481, HE5386/18-1, 19-2, 22-1, 23-1, ERS SFDdM035 and under Germany's Excellence Strategy EXC-2023 Internet of Production 390621612 and under the Excellence Strategy of the Federal Government and the Länder. Support through the EU ITN DATAHYKING is also acknowledged. The authors acknowledge the support of the Banff International Research Station (BIRS) for the Focused Research Group [22frg198] "Novel perspectives in kinetic equations for emerging phenomena", July 17-24, 2022, where part of this work was done.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A: Eigenvalue analysis of the relaxation system

In this appendix we derive eigenvalues, eigenvectors and characteristic variables for the relaxation system (3.1) with given relaxation rate and speed $\varepsilon, \lambda > 0$. We rewrite the system in vector form as

using the notation ${\bf q}^ \varepsilon = (u^ \varepsilon, q^ \varepsilon)$. This is a linear system of balance laws and by diagonalizing the flux matrix we obtain

which reveals $-\lambda$ and $\lambda$ as eigenvalues of the system whereas its eigenvectors are given by ${\bf r}_1 = (- 1/ \lambda, 1)^T$ and ${\bf r}_2 = (1/ \lambda, 1)^T$. The characteristic variables consequently read

Appendix B: Schemes for differing relaxation speeds

In this appendix we provide evolution formulas and schemes for differing relaxation speeds in the coupled relaxation system (3.6). In Section 3.2, the corresponding formulas are, for brevity, only discussed in the simplified case of equal relaxation speeds left and right from the coupling node (see Remark 1).

The coupled semi-discrete scheme Analogously to (3.18) and (3.19) in Section 3.2.2 we obtain by substituting the coupling data (3.12) into (3.16) and (3.17) in the general case

for the evolution of the volumes left to the coupling node and

for the evolution of the volumes right to the coupling node. Consistency in the case $f_1 = f_2$ to the scheme (3.15) in the uncoupled case is only given if also $\lambda_1 = \lambda_2$, see Proposition 1.

The limit scheme If we allow for differing relaxation speeds when substituting the coupling data given by (3.12) and (3.25) into Eq (3.24b) and (3.24c), in analogy to the derivation of (3.26) in Section 3.2.4, we obtain

These formulas then replace (3.24b) and (3.24c) in scheme (3.24). Also in this more general case the limit scheme can be rewritten in the conservative form (3.27a) using modified numerical fluxes given by