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

HAL
open science

# Well-posedness of a non-local model for material flow on conveyor belts

Elena Rossi, Jennifer Kötz, Paola Goatin, Simone Göttlich

► To cite this version:

Elena Rossi, Jennifer Kötz, Paola Goatin, Simone Göttlich. Well-posedness of a non-local model for material flow on conveyor belts. ESAIM: Mathematical Modelling and Numerical Analysis, 2020, 54 (2), pp.679-704. 10.1051/m2an/2019062. hal-02022654

HAL Id: hal-02022654

https://inria.hal.science/hal-02022654

Submitted on 18 Feb 2019

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.
# Well-posedness of a non-local model for material flow on conveyor belts

Elena Rossi<fnref n="1" />

Jennifer Kötz<fnref n="2" />

Paola Goatin<fnref n="1" />

Simone Göttlich<fnref n="2" />

## Abstract

In this paper, we focus on finite volume approximation schemes to solve a non-local material flow model in two space dimensions. Based on the numerical discretisation with dimensional splitting, we prove the convergence of the approximate solutions, where the main difficulty arises in the treatment of the discontinuity occurring in the flux function. In particular, we compare a Roe-type scheme to the well-established Lax-Friedrichs method and provide a numerical study highlighting the benefits of the Roe discretisation. Besides, we also prove the L<fnref n="1" />-Lipschitz continuous dependence on the initial datum, ensuring the uniqueness of the solution.

Keywords: non-local conservation laws; material flow; Roe scheme; Lax-Friedrichs scheme
2010 Mathematics Subject Classification: 35L65, 65M12

## 1 Introduction

In this paper, we consider the Cauchy problem for a non-local scalar conservation law in two space dimensions, namely

$$ \left\{ \begin{array}{ll} \partial_t \rho + \nabla \cdot (\rho \mathbf{v}^{\text{stat}}(x, y) + \rho \mathbf{v}^{\text{dyn}}(\rho)) = 0, & (t, x, y) \in [0, T] \times \mathbb{R}^2, \\ \rho(0, x, y) = \rho_o(x, y), & (x, y) \in \mathbb{R}^2, \end{array} \right. \qquad (1.1) $$

for any $T > 0$, where for the weighting factor $\varepsilon > 0$

$$ \mathbf{v}^{\text{dyn}}(\rho) = H(\rho - \rho_{\text{max}}) \mathbf{I}(\rho) \quad \text{with} \quad \mathbf{I}(\rho) = -\varepsilon \frac{\nabla(\eta * \rho)}{\sqrt{1 + \|\nabla(\eta * \rho)\|^2}}. \qquad (1.2) $$

Here, $H$ denotes the Heaviside function which becomes active whenever the maximal density $\rho_{\text{max}} > 0$ is exceeded.

The above model was introduced in [9], to describe the flow of objects on a conveyor belt. In particular, the unknown function $\rho = \rho(t, x, y)$ is the density of transported parts, $\mathbf{v}^{\text{stat}}$ is the velocity vector field induced by the conveyor belt, which is constant in time. So, below the maximal density, parts are transported with the velocity of the conveyor belt. On the other

<footnote n="1">¹Inria Sophia Antipolis - Méditerranée, Université Côte d'Azur, Inria, CNRS, LJAD, 2004 route des Lucioles - BP 93, 06902 Sophia Antipolis Cedex, France. E-mail: {elena.rossi, paola.goatin}@inria.fr</footnote>

<footnote n="2">²University of Mannheim, Department of Mathematics, 68131 Mannheim, Germany. Email: jkoetz@mail.uni-mannheim.de, goettlich@uni-mannheim.de</footnote>
hand, the dynamic velocity vector field $v^{dyn}$ is active only at high densities and accounts for colliding parts through the operator $I(\rho)$. The negative gradient of the convolution $\eta * \rho$ pushes the mass towards lower density regions. The particular choice of the collision operator $I(\rho)$ was introduced in [5] to describe crowd dynamics.

Conservation laws with non-local flux function have been recently introduced in the literature to describe transport phenomena accounting for non-local interaction effects among agents, such as road traffic flow [3] or pedestrian dynamics [5]. General well-posedness results have been provided by [2] in the scalar one-dimensional case, while [1] deals with systems of non-local conservation laws in multi-space dimensions. Even if the latter result applies to our problem (1.1), estimates in [1] where obtained for general flux functions using finite volume approximate solutions constructed via Lax-Friedrichs scheme. The aim of the present paper is instead to derive sharp estimates for the Roe scheme, which is known to give less diffusive solutions and is therefore computationally more convenient, especially in the case of non-local problems in multi-D. The same remark holds for the $L^1$-stability estimates, which have been derived from scratch even if more general results are present in the literature [12, 13].

We finally remark that, even if the original model proposes the use of the discontinuous Heaviside function, the stability of the numerical schemes requires a smooth approximation of it. Indeed, it can be noticed that $L^\infty$-norm of the derivative $H'$ appears in the estimates (for instance in the CFL condition (3.9) that guarantees the BV-estimates), which blow up with it.

The paper is organised as follows: for a better overview, we introduce our main results in Section 2 and start then to prove convergence of the approximate solution constructed by the Roe scheme in Section 3. We also add a proof of the Lipschitz continuous dependence on the initial data in Section 4. Sections 5 and 6 are devoted to the comparison of results obtained by the Lax-Friedrichs method. The numerical results emphasise the good performance of the Roe scheme.

## 2 Main results

Throughout the paper, we will denote $\mathcal{I}(r, s) := [\min\{r, s\}, \max\{r, s\}]$, for any $r, s \in \mathbb{R}$. We require the following assumptions to hold:

$$ (v) \quad v^{stat} \in C^2(\mathbb{R}^2; \mathbb{R}^2). $$

(H) The function $H$ is a smooth approximation of the Heaviside function such that, setting

$$ f(r) = r H(r - \rho_{max}), $$

the function $f$ has bounded derivative. In particular, we denote by $L_f$ the Lipschitz constant of the function $f$:

$$ L_f = \|f'\|_{L^\infty(\mathbb{R})}. \qquad (2.1) $$

$$ (\eta) \quad \eta \in (C^3 \cap W^{3,\infty})(\mathbb{R}^2; \mathbb{R}^+). $$

Recall the definition of solution to the Cauchy problem (1.1), see also [1, 2, 5].
**Definition 2.1.** Let $\rho_o \in L^\infty(\mathbb{R}^2; \mathbb{R}^+)$. A map $\rho: [0, T] \to L^\infty(\mathbb{R}^2; \mathbb{R})$ is a solution to (1.1) if it is a Kružkov solution to

$$ \left\{ \begin{array}{ll} \partial_t \rho + \nabla \cdot g(t, x, y, \rho) = 0, & (t, x, y) \in [0, T] \times \mathbb{R}^2, \\ \rho(0, x, y) = \rho_o(x, y), & (x, y) \in \mathbb{R}^2, \end{array} \right. \qquad (2.2) $$

$$ \text{with } g(t, x, y, \rho) = \rho \mathbf{v}^{\text{stat}}(x, y) - \varepsilon \rho H(\rho - \rho_{\text{max}}) \frac{\nabla(\eta * \rho)}{\sqrt{1 + \|\nabla(\eta * \rho)\|^2}}. $$

Above, for the definition of Kružkov solution we refer to [11, Definition 1].

**Theorem 2.2.** Let $\rho_o \in (L^\infty \cap BV)(\mathbb{R}^2; \mathbb{R}^+)$. Let assumptions $(\mathbf{v})$, $(H)$ and $(\eta)$ hold. Then, for all $T > 0$, there exists a unique weak entropy solution $\rho \in (L^\infty \cap BV)([0, T] \times \mathbb{R}^2; \mathbb{R}^+)$ to problem (1.1). Moreover, the following estimates hold: for all $t \in [0, T]$

$$ \begin{array}{l} \|\rho(t)\|_{L^1} = \|\rho_o\|_{L^1}, \\ \|\rho(t)\|_{L^\infty} \le \|\rho_o\|_{L^\infty} e^{C_\infty t}, \\ \text{TV}(\rho(t)) \le \text{TV}(\rho_o) e^{2t\mathcal{K}_1} + \frac{2\mathcal{K}_2}{\mathcal{K}_1} (e^{2t\mathcal{K}_1} - 1), \\ \|\rho(t) - \rho(t - \tau)\|_{L^1} \le 2C_t(t)\tau, \quad \text{for } \tau > 0, \end{array} $$

where $C_\infty$ is defined in (3.12), $\mathcal{K}_1$ is defined in (3.16), $\mathcal{K}_2$ is defined in (3.17) and $C_t(t)$ is as in (3.38).

The proof of Theorem 2.2 is standard, see [1, Theorem 2.3], which refers to [15]. The uniqueness is ensured by Proposition 4.1, which provides the Lipschitz continuous dependence estimate of solutions to (1.1) on the initial data.

The bounds presented in Theorem 2.2 are obtained by passing to the limit in the corresponding discrete bounds.

## 3 Existence

Introduce the uniform mesh of width $\Delta x$ along the $x$-axis and $\Delta y$ along the $y$-axis, and a time step $\Delta t$ subject to a CFL condition, specified later on. For $k \in \mathbb{Z}$ set

$$ \begin{array}{ll} x_k = (k - 1/2)\Delta x, & y_k = (k - 1/2)\Delta y, \\ x_{k+1/2} = k\Delta x, & y_{k+1/2} = k\Delta y, \end{array} $$

where $(x_{i+1/2}, y_j)$ and $(x_i, y_{j+1/2})$ denote the cells interfaces and $(x_i, y_j)$ are the cells centres. Set $N_T = \lfloor T/\Delta t \rfloor$ and let $t^n = n \Delta t$, $n = 0, \dots, N_T$, be the time mesh. Set $\lambda_x = \Delta t/\Delta x$ and $\lambda_y = \Delta t/\Delta y$, and let $\alpha, \beta \ge 1$ be the viscosity coefficients.

For the sake of shortness, sometimes we will use also the notation $x_{i,j} = (x_i, y_j)$.

We approximate the initial datum as follows: for $i, j \in \mathbb{Z}$

$$ \rho_{i,j}^0 = \frac{1}{\Delta x \Delta y} \int_{x_{i-1/2}}^{x_{i+1/2}} \int_{y_{i-1/2}}^{y_{i+1/2}} \rho_o(x, y) \, dx \, dy, $$
and we define a piece-wise constant solution to (1.1) as

$$ \rho_{\Delta}(t, x, y) = \rho_{i,j}^{n} \quad \text{for} \quad \begin{cases} t \in [t^n, t^{n+1}], & n = 0, \dots, N_T - 1, \\ x \in [x_{i-1/2}, x_{i+1/2}], & \text{where } i \in \mathbb{Z}, \\ y \in [y_{j-1/2}, y_{j+1/2}], & \text{where } j \in \mathbb{Z}. \end{cases} \tag{3.1} $$

through a modified Roe-type scheme with dimensional splitting, exactly as in [9]:

**Algorithm 3.1.**

$$ V_1(x, y, u, w) = v_1^{\text{stat}}(x, y) u + \min\{0, v_1^{\text{stat}}(x, y)\}(w - u) \tag{3.2} $$

$$ F(u, w, J(t, x, y)) = J(t, x, y) f(u) + \min\{0, J(t, x, y)\} (f(w) - f(u)) \tag{3.3} $$

$$ V_2(x, y, u, w) = v_2^{\text{stat}}(x, y) u + \min\{0, v_2^{\text{stat}}(x, y)\}(w - u) \tag{3.4} $$

for $n = 0, \dots, N_T - 1$

$$ \rho_{i,j}^{n+1/2} = \rho_{i,j}^{n} - \lambda_x \left[ V_1(x_{i+1/2,j}, \rho_{i,j}^{n}, \rho_{i+i,j}^{n}) - V_1(x_{i-1/2,j}, \rho_{i-1,j}^{n}, \rho_{i,j}^{n}) \right. \\ \left. + F(\rho_{i,j}^{n}, \rho_{i+1,j}^{n}, J_1^n(x_{i+1/2,j}) - F(\rho_{i-1,j}^{n}, \rho_{i,j}^{n}, J_1^n(x_{i-1/2,j})) \right] \tag{3.5} $$

$$ \rho_{i,j}^{n+1} = \rho_{i,j}^{n+1/2} - \lambda_y \left[ V_2(x_{i,j+1/2}, \rho_{i,j}^{n}, \rho_{i,j+1}^{n}) - V_2(x_{i,j-1/2}, \rho_{i,j-1}^{n}, \rho_{i,j}^{n}) \right. \\ \left. + F(\rho_{i,j}^{n}, \rho_{i,j+1}^{n}, J_2^n(x_{i,j+1/2}) - F(\rho_{i,j-1}^{n}, \rho_{i,j}^{n}, J_2^n(x_{i,j-1/2})) \right] \tag{3.6} $$

end

Above, we set $I(\rho^n)(x, y) = (J_1^n(x, y), J_2^n(x, y))$ and the convolution products are computed through the following quadrature formula, for $k = 1, 2$,

$$ (\partial_i \eta * \rho)(x_i, y_j) = \Delta x \Delta y \sum_{k,l \in \mathbb{Z}} \rho_{k,l} \partial_i \eta(x_{i-k}, y_{j-l}), \tag{3.7} $$

where $\partial_1 \eta = \partial_x \eta$ and $\partial_2 \eta = \partial_y \eta$. Remark that the choice of evaluating the numerical flux at $t^n$ for both fractional steps allows to save computational time, because the convolution products (3.7) are computed only once per time step.

Introduce the following notation, which is of use below:

$$ v_{i+1/2} = v_1^{\text{stat}}(x_{i+1/2}). \tag{3.8} $$

## 3.1 Positivity

In the case of positive initial datum, we prove that under a suitable CFL condition the approximate solution to (1.1) constructed via the Algorithm 3.1 preserves the positivity.

**Lemma 3.2. (Positivity)** Let $\rho_o \in L^\infty(\mathbb{R}^2; \mathbb{R}^+)$, Let $(v)$, $(H)$, and $(\eta)$ hold. Assume that

$$ \lambda_x \le \frac{1}{2(\varepsilon + \|v_1^{\text{stat}}\|_{L^\infty})}, \qquad \lambda_y \le \frac{1}{2(\varepsilon + \|v_2^{\text{stat}}\|_{L^\infty})}. \tag{3.9} $$

Then, for all $t > 0$ and $(x, y) \in \mathbb{R}^2$, the piece-wise constant approximate solution $\rho_\Delta$ (3.1) constructed through Algorithm 3.1 is such that $\rho_\Delta(t, x, y) \ge 0$.
**Proof.** Fix $n$ between 0 and $N_T - 1$ and assume that $\rho_{i,j}^n \ge 0$ for all $i, j \in \mathbb{Z}$. Consider (3.5), with the notation (3.2) and (3.3), and observe that:

$$
\begin{align*}
& V_1(x_{i+1/2,j}, \rho_{i,j}^n, \rho_{i+1,j}^n) + F(\rho_{i,j}^n, \rho_{i+1,j}^n, J_1^n(x_{i+1/2,j})) \\
\le {}& v_1^{\text{stat}}(x_{i+1/2,j}) \rho_{i,j}^n + J_1^n(x_{i+1/2,j}) f(\rho_{i,j}^n) \le \left( \|v_1^{\text{stat}}\|_{\mathbf{L}^\infty} + \varepsilon \right) \rho_{i,j}^n, \\
& V_1(x_{i-1/2,j}, \rho_{i-1,j}^n, \rho_{i,j}^n) + F(\rho_{i-1,j}^n, \rho_{i,j}^n, J_1^n(x_{i-1/2,j})) \\
\ge {}& v_1^{\text{stat}}(x_{i-1/2,j}) \rho_{i,j}^n + J_1^n(x_{i-1/2,j}) f(\rho_{i,j}^n) \ge - \left( \|v_1^{\text{stat}}\|_{\mathbf{L}^\infty} + \varepsilon \right) \rho_{i,j}^n.
\end{align*}
$$

Therefore, by (3.5),

$$
\rho_{i,j}^{n+1/2} \ge \rho_{i,j}^{n} - 2 \lambda_x \left( \|v_1^{\text{stat}}\|_{\mathbf{L}^{\infty}} + \varepsilon \right) \rho_{i,j}^{n} \ge 0,
$$

thanks to the CFL condition (3.9). Starting from (3.6), an analogous argument shows that
$\rho_{i,j}^{n+1} \ge 0$, concluding the proof. $\square$

3.2 L¹ bound

The following result on the $\mathbf{L}^1$ bound follows from the conservation property of the Roe scheme.

**Lemma 3.3. (L¹ bound)** Let $\rho_o \in \mathbf{L}^\infty(\mathbb{R}^2; \mathbb{R}^+)$. Let $(v), (H), (\eta)$ and (3.9) hold. Then, for all $t > 0$ and $(x, y) \in \mathbb{R}^2$, $\rho_\Delta$ (3.1) constructed through Algorithm 3.1 satisfies

$$
\| \rho_{\Delta}(t, \cdot, \cdot) \|_{\mathbf{L}^1(\mathbb{R}^2)} = \| \rho_o \|_{\mathbf{L}^1(\mathbb{R}^2)}. \quad (3.10)
$$

3.3 L$^{\infty}$ bound

**Lemma 3.4. (L$^{\infty}$ bound)** Let $\rho_o \in \mathbf{L}^{\infty}(\mathbb{R}^2; \mathbb{R}^+)$. Let $(v), (H), (\eta)$ and (3.9) hold. Then, for all $t > 0$ and $(x, y) \in \mathbb{R}^2$, $\rho_\Delta$ (3.1) constructed through Algorithm 3.1 satisfies

$$
\| \rho_{\Delta}(t, \cdot, \cdot) \|_{\mathbf{L}^{\infty}(\mathbb{R}^2)} \le \| \rho_o \|_{\mathbf{L}^{\infty}} e^{C_{\infty} t}, \quad (3.11)
$$

where

$$
C_{\infty} = \left\| \partial_x v_1^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + \left\| \partial_y v_2^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + 4\varepsilon \left\| \nabla^2 \eta \right\|_{\mathbf{L}^{\infty}} \| \rho_o \|_{\mathbf{L}^1}. \quad (3.12)
$$

**Proof.** Omitting the dependencies on $j$ and exploiting the notation introduced in (3.8), we observe that (3.5) attains its maximum for $v_{i+1/2} < 0$, $v_{i-1/2} \ge 0$, $J_1^n(x_{i+1/2}) < 0$ and $J_1^n(x_{i-1/2}) \ge 0$. In this case

$$
\rho_{i,j}^{n+1/2} \le \rho_i^n - \lambda_x \left( \rho_{i+1}^n v_{i+1/2} + J_1^n(x_{i+1/2}) f(\rho_{i+1}^n) \right) + \lambda_x \left( \rho_{i-1}^n v_{i-1/2} + J_1^n(x_{i-1/2}) f(\rho_{i-1}^n) \right),
$$

where we use the positivity of each $\rho_i^n$ and of the function $f$ and discard all the terms giving
a negative contribution. Moreover, since $v_{i+1/2} < 0$ and $v_{i-1/2} \ge 0$,

$$
\lambda_x \left( -\rho_{i+1}^n v_{i+1/2} + \rho_{i-1}^n v_{i-1/2} \right) \le \lambda_x \|\rho^n\|_{\mathbf{L}^\infty} \left( -v_{i+1/2} + v_{i-1/2} \right)
$$
$$ = \lambda_x \|\rho^n\|_{\mathbf{L}^\infty} (-\Delta x) \partial_x v_1^{\text{stat}}(\hat{x}_i) $$

with $\hat{x}_i \in ]x_{i-1/2}, x_{i+1/2}[$. In a similar way, since $J_1^n(x_{i+1/2}) < 0$ and $J_1^n(x_{i-1/2}) \ge 0$, exploiting also the fact that $f(r) \le r$ for all $r \ge 0$, we get

$$ \begin{aligned} \lambda_x \left( -J_1^n(x_{i+1/2}) f(\rho_{i+1}^n) + J_1^n(x_{i-1/2}) f(\rho_{i-1}^n) \right) &\le \lambda_x \|\rho^n\|_{\mathbf{L}^\infty} \left( -J_1^n(x_{i+1/2}) + J_1^n(x_{i-1/2}) \right) \\ &\le \lambda_x \|\rho^n\|_{\mathbf{L}^\infty} 2 \varepsilon \Delta x \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \|\rho^n\|_{\mathbf{L}^1}, \end{aligned} $$

thanks to (A.2). Therefore,

$$ \rho_{i,j}^{n+1/2} \le \|\rho^n\|_{\mathbf{L}^\infty} \left[ 1 + \Delta t \left( \left\| \partial_x v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2\varepsilon \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \right]. $$

In a similar way we get

$$ \rho_{i,j}^{n+1} \le \|\rho^{n+1/2}\|_{\mathbf{L}^\infty} \left[ 1 + \Delta t \left( \left\| \partial_y v_2^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2\varepsilon \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \right], $$

concluding the proof.

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

## 3.4 BV bound

**Proposition 3.5. (BV estimate in space)** Let $\rho_o \in (\mathbf{L}^\infty \cap \text{BV})(\mathbb{R}^2; \mathbb{R}^+)$, Let $(v), (H)$, and $(\eta)$ hold. Assume that

$$ \lambda_x \le \frac{1}{3(\varepsilon L_f + \|v_1^{\text{stat}}\|_{\mathbf{L}^\infty})}, \qquad \lambda_y \le \frac{1}{3(\varepsilon L_f + \|v_2^{\text{stat}}\|_{\mathbf{L}^\infty})}. \qquad (3.13) $$

Then, for all $t > 0$, $\rho_\Delta$ in (3.1) constructed through Algorithm 3.1 satisfies the following estimate: for all $n = 0, \dots, N_T$,

$$ \sum_{i,j \in \mathbb{Z}} \left( \Delta y \left| \rho_{i+1,j}^n - \rho_{i,j}^n \right| + \Delta x \left| \rho_{i,j+1}^n - \rho_{i,j}^n \right| \right) \le C_x(t^n), \qquad (3.14) $$

where

$$ C_x(t) = e^{2t\mathcal{K}_1} \sum_{i,j \in \mathbb{Z}} \left( \Delta x \left| \rho_{i,j+1}^0 - \rho_{i,j}^0 \right| + \Delta y \left| \rho_{i+1,j}^0 - \rho_{i,j}^0 \right| \right) + \frac{2\mathcal{K}_2}{\mathcal{K}_1} (e^{2t\mathcal{K}_1} - 1), \qquad (3.15) $$

with

$$ \mathcal{K}_1 = 6 \left( \left\| \nabla v^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2\varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right), \qquad (3.16) $$

$$ \mathcal{K}_2 = \left( 4\varepsilon \left( c_1 \|\rho_o\|_{\mathbf{L}^1} + c_2 \|\rho_o\|_{\mathbf{L}^1}^2 \right) + 3 \left\| \nabla^2 v^{\text{stat}} \right\|_{\mathbf{L}^\infty} \right) \|\rho_o\|_{\mathbf{L}^1}, \qquad (3.17) $$

and $c_1, c_2$ are defined in (A.6).

**Remark 3.6.** Observe that the CFL conditions (3.13) are stricter than (3.9).
**Proof.** We follow the idea of [1, Lemma 2.7]. First consider the term

$$ \sum_{i,j \in \mathbb{Z}} \Delta y \left| \rho_{i+1,j}^{n+1/2} - \rho_{i,j}^{n+1/2} \right|. $$

In particular, fixing $i, j \in \mathbb{Z}$ and omitting the dependencies on $y_j$ for the sake of simplicity, by (3.5) we get

$$ \begin{aligned} \rho_{i+1}^{n+1/2} - \rho_i^{n+1/2} &= \rho_{i+1}^n - \rho_i^n - \lambda_x \left[ \begin{aligned} &V_1(x_{i+3/2}, \rho_{i+1}^n, \rho_{i+2}^n) + F(\rho_{i+1}^n, \rho_{i+2}^n, J_1^n(x_{i+3/2})) \\ &- V_1(x_{i+1/2}, \rho_i^n, \rho_{i+1}^n) - F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+1/2})) \\ &- V_1(x_{i+1/2}, \rho_i^n, \rho_{i+1}^n) - F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+1/2})) \\ &+ V_1(x_{i-1/2}, \rho_i^n, \rho_{i-1}^n) + F(\rho_{i-1}^n, \rho_i^n, J_1^n(x_{i-1/2})) \end{aligned} \right] \\ &\quad \pm \lambda_x \left[ \begin{aligned} &V_1(x_{i+3/2}, \rho_i^n, \rho_{i+1}^n) + F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+3/2})) \\ &- V_1(x_{i+1/2}, \rho_i^n, \rho_{i-1}^n) - F(\rho_{i-1}^n, \rho_i^n, J_1^n(x_{i+1/2})) \end{aligned} \right] \\ &= \mathcal{A}_{i,j}^n - \lambda_x \mathcal{B}_{i,j}^n, \end{aligned} $$

where we set

$$ \begin{aligned} \mathcal{A}_{i,j}^n &= \rho_{i+1}^n - \rho_i^n - \lambda_x \left[ \begin{aligned} &V_1(x_{i+3/2}, \rho_{i+1}^n, \rho_{i+2}^n) + F(\rho_{i+1}^n, \rho_{i+2}^n, J_1^n(x_{i+3/2})) \\ &- V_1(x_{i+1/2}, \rho_i^n, \rho_{i+1}^n) - F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+1/2})) \\ &+ V_1(x_{i+1/2}, \rho_i^n, \rho_{i-1}^n) + F(\rho_{i-1}^n, \rho_i^n, J_1^n(x_{i+1/2})) \\ &- V_1(x_{i+3/2}, \rho_i^n, \rho_{i+1}^n) - F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+3/2})) \end{aligned} \right], \\ \mathcal{B}_{i,j}^n &= V_1(x_{i+3/2}, \rho_i^n, \rho_{i+1}^n) + F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+3/2})) \end{aligned} $$

$$ \begin{aligned} \mathcal{B}_{i,j}^n &= V_1(x_{i+3/2}, \rho_i^n, \rho_{i+1}^n) + F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+3/2})) \\ &\quad - V_1(x_{i+1/2}, \rho_{i-1}^n, \rho_i^n) - F(\rho_{i-1}^n, \rho_i^n, J_1^n(x_{i+1/2})) \\ &\quad + V_1(x_{i-1/2}, \rho_{i-1}^n, \rho_i^n) + F(\rho_{i-1}^n, \rho_i^n, J_1^n(x_{i-1/2})) \\ &\quad - V_1(x_{i+1/2}, \rho_i^n, \rho_{i+1}^n) - F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+1/2})). \end{aligned} $$

For the sake of shortness, introduce the following notation

$$ H_{k,l}^{n}(u, w) = V_{1}(x_{k,l}, u, w) + F(u, w, J_{1}^{n}(x_{k,l})), \quad (3.18) $$

so that, dropping the $j$ dependencies, $\mathcal{A}_{i,j}^n$ reads

$$ \begin{aligned} \mathcal{A}_{i,j}^n &= \rho_{i+1}^n - \rho_i^n - \lambda_x \left[ \begin{aligned} &H_{i+3/2}^n(\rho_{i+1}^n, \rho_{i+2}^n) - H_{i+1/2}^n(\rho_i^n, \rho_{i+1}^n) \\ &+ H_{i+1/2,j}^n(\rho_{i-1,j}^n, \rho_{i,j}^n) - H_{i+3/2,j}^n(\rho_{i,j}^n, \rho_{i+1,j}^n) \end{aligned} \right] \\ &= \rho_{i+1}^n - \rho_i^n - \lambda_x \frac{H_{i+3/2}^n(\rho_{i+1}^n, \rho_{i+2}^n) - H_{i+3/2}^n(\rho_{i+1}^n, \rho_{i+1}^n)}{\rho_{i+2}^n - \rho_{i+1}^n} (\rho_{i+2}^n - \rho_{i+1}^n) \\ &\quad - \lambda_x \frac{H_{i+3/2}^n(\rho_{i+1}^n, \rho_{i+1}^n) - H_{i+3/2}^n(\rho_i^n, \rho_{i+1}^n)}{\rho_{i+1}^n - \rho_i^n} (\rho_{i+1}^n - \rho_i^n) \\ &\quad + \lambda_x \frac{H_{i+1/2}^n(\rho_i^n, \rho_{i+1}^n) - H_{i+1/2}^n(\rho_i^n, \rho_i^n)}{\rho_{i+1}^n - \rho_i^n} (\rho_{i+1}^n - \rho_i^n) \end{aligned} $$
$$ 
\begin{aligned}
&+ \lambda_x \frac{H_{i+1/2}^n(\rho_i^n, \rho_i^n) - H_{i+1/2}^n(\rho_{i-1}^n, \rho_i^n)}{\rho_i^n - \rho_{i-1}^n} (\rho_i^n - \rho_{i-1}^n) \\
&= \delta_{i+1}^n (\rho_{i+2}^n - \rho_{i+1}^n) + \vartheta_i^n (\rho_i^n - \rho_{i-1}^n) + (1 - \delta_i^n - \vartheta_{i+1}^n) (\rho_{i+1}^n - \rho_i^n),
\end{aligned}
$$

where

$$ \delta_i^n = \begin{cases} -\lambda_x \frac{H_{i+1/2}^n(\rho_i^n, \rho_{i+1}^n) - H_{i+1/2}^n(\rho_i^n, \rho_i^n)}{\rho_{i+1}^n - \rho_i^n} & \text{if } \rho_i^n \neq \rho_{i+1}^n, \\ 0 & \text{if } \rho_i^n = \rho_{i+1}^n, \end{cases} \quad (3.19) $$

$$ \vartheta_i^n = \begin{cases} \lambda_x \frac{H_{i+1/2}^n(\rho_i^n, \rho_i^n) - H_{i+1/2}^n(\rho_{i-1}^n, \rho_i^n)}{\rho_i^n - \rho_{i-1}^n} & \text{if } \rho_i^n \neq \rho_{i-1}^n, \\ 0 & \text{if } \rho_i^n = \rho_{i-1}^n. \end{cases} \quad (3.20) $$

Exploiting (3.18), observe that, whenever $\rho_i^n \neq \rho_{i+1}^n$,

$$ 
\begin{aligned}
\delta_i^n &= - \frac{\lambda_x}{\rho_{i+1}^n - \rho_i^n} \left[ V_1(x_{i+1/2}, \rho_i^n, \rho_{i+1}^n) + F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+1/2})) \right. \\
&\qquad \left. - V_1(x_{i+1/2}, \rho_i^n, \rho_i^n) - F(\rho_i^n, \rho_i^n, J_1^n(x_{i+1/2})) \right] \\
&= - \frac{\lambda_x}{\rho_{i+1}^n - \rho_i^n} \left[ \min \{0, v_1^{\text{stat}}(x_{i+1/2})\} (\rho_{i+1}^n - \rho_i^n) + \min \{0, J_1^n(x_{i+1/2})\} (f(\rho_{i+1}^n) - f(\rho_i^n)) \right] \\
&= - \lambda_x \left( \min \{0, v_1^{\text{stat}}(x_{i+1/2})\} + \min \{0, J_1^n(x_{i+1/2})\} f'(r_{i+1/2}^n) \right),
\end{aligned}
$$

with $r_{i+1/2}^n \in \mathcal{I}(\rho_i^n, \rho_{i+1}^n)$. Since $f'(r) \ge 0$ and by (3.13) we get

$$ \delta_i^n \in [0, \frac{1}{3}]. $$

In a similar way one can prove that $\vartheta_i^n \in [0, 1/3]$. Thus,

$$ \sum_{i,j \in \mathbb{Z}} |A_{i,j}^n| \le \sum_{i,j \in \mathbb{Z}} |\rho_{i+1,j}^n - \rho_{i,j}^n|. \quad (3.21) $$

We pass now to $B_{i,j}^n$. Consider separately the terms involving $V_1$ and those involving $F$. Observe that the maps

$$ x \mapsto \min \{0, v_1^{\text{stat}}(x)\}, \qquad x \mapsto \min \{0, J_1^n(x)\} $$

are Lipschitz continuous, with constant respectively $\|\partial_x v_1^{\text{stat}}\|_{L^\infty}$ and $2\varepsilon\|\nabla^2\eta\|_{L^\infty}\|\rho_o\|_{L^1}$. Exploiting (3.2) we get:

$$ 
\begin{aligned}
& V_1(x_{i+3/2}, \rho_i^n, \rho_{i+1}^n) - V_1(x_{i+1/2}, \rho_{i-1}^n, \rho_i^n) + V_1(x_{i-1/2}, \rho_{i-1}^n, \rho_i^n) - V_1(x_{i+1/2}, \rho_i^n, \rho_{i+1}^n) \\
&= v_1^{\text{stat}}(x_{i+3/2})\rho_i^n + \min \{0, v_1^{\text{stat}}(x_{i+3/2})\} (\rho_{i+1}^n - \rho_i^n) \\
&\quad - v_1^{\text{stat}}(x_{i+1/2})\rho_i^n - \min \{0, v_1^{\text{stat}}(x_{i+1/2})\} (\rho_{i+1}^n - \rho_i^n)
\end{aligned}
$$
$$
\begin{align*}
&+ v_1^{\text{stat}}(x_{i-1/2})\rho_{i-1}^n + \min\{0, v_1^{\text{stat}}(x_{i-1/2})\} (\rho_i^n - \rho_{i-1}^n) \\
&- v_1^{\text{stat}}(x_{i+1/2})\rho_{i-1}^n - \min\{0, v_1^{\text{stat}}(x_{i+1/2})\} (\rho_i^n - \rho_{i-1}^n) \\
&\pm (v_1^{\text{stat}}(x_{i-1/2}) - v_1^{\text{stat}}(x_{i+1/2})) \rho_i^n \\
&= (v_1^{\text{stat}}(x_{i+3/2}) - 2v_1^{\text{stat}}(x_{i+1/2}) + v_1^{\text{stat}}(x_{i-1/2})) \rho_i^n \\
&\quad + (v_1^{\text{stat}}(x_{i-1/2}) - v_1^{\text{stat}}(x_{i+1/2})) (\rho_{i-1}^n - \rho_i^n) \\
&\quad + \left(\min\{0, v_1^{\text{stat}}(x_{i+3/2})\} - \min\{0, v_1^{\text{stat}}(x_{i+1/2})\}\right) (\rho_{i+1}^n - \rho_i^n) \\
&\quad + \left(\min\{0, v_1^{\text{stat}}(x_{i-1/2})\} - \min\{0, v_1^{\text{stat}}(x_{i+1/2})\}\right) (\rho_i^n - \rho_{i-1}^n) \\
&\le 2(\Delta x)^2 \left\|\partial_{xx}^2 v_1^{\text{stat}}\right\|_{\mathbf{L}^\infty} |\rho_i^n| + \Delta x \left\|\partial_x v_1^{\text{stat}}\right\|_{\mathbf{L}^\infty} \left(|\rho_{i+1}^n - \rho_i^n| + 2|\rho_i^n - \rho_{i-1}^n|\right), \tag{3.22}
\end{align*}
$$

since

$$
\begin{align*}
v_1^{\text{stat}}(x_{i+3/2}) - 2 v_1^{\text{stat}}(x_{i+1/2}) + v_1^{\text{stat}}(x_{i-1/2}) &= \Delta x \partial_x v_1^{\text{stat}}(\xi_{i+1}) - \Delta x \partial_x v_1^{\text{stat}}(\xi_i) \\
&= \Delta x (\xi_{i+1} - \xi_i) \partial_{xx}^2 v_1^{\text{stat}}(\zeta_{i+1/2}),
\end{align*}
$$

with $\xi_i \in ]x_{i-1/2}, x_{i+1/2}[ \text{ and } \zeta_{i+1/2} \in ]\xi_i, \xi_{i+1}[.$

Similarly, exploiting (3.3) we obtain

$$
\begin{align*}
& F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+3/2})) - F(\rho_{i-1}^n, \rho_i^n, J_1^n(x_{i+1/2})) \\
&\quad + F(\rho_{i-1}^n, \rho_i^n, J_1^n(x_{i-1/2})) - F(\rho_i^n, \rho_{i+1}^n, J_1^n(x_{i+1/2})) \\
&= J_1^n(x_{i+3/2})f(\rho_i^n) + \min\{0, J_1^n(x_{i+3/2})\} (f(\rho_{i+1}^n) - f(\rho_i^n)) \\
&\quad - J_1^n(x_{i+1/2})f(\rho_i^n) - \min\{0, J_1^n(x_{i+1/2})\} (f(\rho_{i+1}^n) - f(\rho_i^n)) \\
&\quad + J_1^n(x_{i-1/2})f(\rho_i^n) + \min\{0, J_1^n(x_{i-1/2})\} (f(\rho_i^n) - f(\rho_{i-1}^n)) \\
&\quad - J_1^n(x_{i+1/2})f(\rho_i^n) - \min\{0, J_1^n(x_{i+1/2})\} (f(\rho_i^n) - f(\rho_{i-1}^n)) \\
&\quad \pm (J_1^n(x_{i-1/2}) - J_1^n(x_{i+1/2})) f(\rho_i^n) \\
&= (J_1^n(x_{i+3/2}) - 2J_1^n(x_{i+1/2}) + J_1^n(x_{i-1/2})) f(\rho_i^n) \\
&\quad + (J_1^n(x_{i-1/2}) - J_1^n(x_{i+1/2})) (f(\rho_{i-1}^n) - f(\rho_i^n)) \\
&\quad + \left(\min\{0, J_1^n(x_{i+3/2})\} - \min\{0, J_1^n(x_{i+1/2})\}\right) (f(\rho_{i+1}^n) - f(\rho_i^n)) \\
&\quad + \left(\min\{0, J_1^n(x_{i-1/2})\} - \min\{0, J_1^n(x_{i+1/2})\}\right) (f(\rho_i^n) - f(\rho_{i-1}^n)) \\
&\le 2 \varepsilon (\Delta x)^2 C |\rho_i^n| + 2 \varepsilon L_f \Delta x \left\|\nabla^2 \eta\right\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \left(|\rho_{i+1}^n - \rho_i^n| + 2|\rho_i^n - \rho_{i-1}^n|\right), \tag{3.23}
\end{align*}
$$

where we used the fact that $f(r) \le r$, (A.2) and (A.4), with the notation (A.6). Collecting
together (3.22) and (3.23) we therefore obtain

$$
|\mathcal{B}_{i,j}^n| \le 2 (\Delta x)^2 \left( \left\| \partial_{xx}^2 v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + \varepsilon C \right) |\rho_i^n|
$$
$$+ \Delta x \left( \left\| \partial_x v_1^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + 2 \varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^{\infty}} \| \rho_o \|_{\mathbf{L}^1} \right) \left( |\rho_{i+1}^n - \rho_i^n| + 2 |\rho_i^n - \rho_{i-1}^n| \right),$$

so that

$$\sum_{i,j \in \mathbb{Z}} \lambda_x |\mathcal{B}_{i,j}^n| \le 3 \Delta t \left( \left\| \partial_x v_1^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + 2 \varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^{\infty}} \| \rho_o \|_{\mathbf{L}^1} \right) \sum_{i,j \in \mathbb{Z}} |\rho_{i+1}^n - \rho_i^n| \\ + 2 \Delta t \left( \left\| \partial_{xx}^2 v_1^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + \varepsilon C \right) \Delta x \sum_{i,j \in \mathbb{Z}} |\rho_i^n|. \quad (3.24)$$

Therefore, by (3.21) and (3.24), using also Lemma 3.3

$$\sum_{i,j \in \mathbb{Z}} \Delta y |\rho_{i+1,j}^{n+1/2} - \rho_{i,j}^{n+1/2}| \\ \le \sum_{i,j \in \mathbb{Z}} \Delta y \left( |\mathcal{A}_{i,j}^n| + \lambda_x |\mathcal{B}_{i,j}^n| \right) \\ \le \left[ 1 + 3 \Delta t \left( \left\| \partial_x v_1^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + 2 \varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^{\infty}} \| \rho_o \|_{\mathbf{L}^1} \right) \right] \sum_{i,j \in \mathbb{Z}} \Delta y |\rho_{i+1}^n - \rho_i^n| \quad (3.25) \\ + 2 \Delta t \left[ \left\| \partial_{xx}^2 v_1^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + \varepsilon C \right] \| \rho_o \|_{\mathbf{L}^1}.$$

Now pass to the term

$$\sum_{i,j \in \mathbb{Z}} \Delta x |\rho_{i,j+1}^{n+1/2} - \rho_{i,j}^{n+1/2}|.$$

Fix $i, j \in \mathbb{Z}$ and exploit (3.5) again to get

$$\begin{aligned} & \rho_{i,j+1}^{n+1/2} - \rho_{i,j}^{n+1/2} \\ &= \rho_{i,j+1}^{n} - \rho_{i,j}^{n} - \lambda_x \left[ \begin{aligned} & V_1(x_{i+1/2,j+1}, \rho_{i,j+1}^{n}, \rho_{i+1,j+1}^{n}) + F(\rho_{i,j+1}^{n}, \rho_{i+1,j+1}^{n}, J_1^{n}(x_{i+1/2,j+1})) \\ & - V_1(x_{i-1/2,j+1}, \rho_{i-1,j+1}^{n}, \rho_{i,j+1}^{n}) - F(\rho_{i-1,j+1}^{n}, \rho_{i,j+1}^{n}, J_1^{n}(x_{i-1/2,j+1})) \\ & - V_1(x_{i+1/2,j}, \rho_{i,j}^{n}, \rho_{i+1,j}^{n}) - F(\rho_{i,j}^{n}, \rho_{i+1,j}^{n}, J_1^{n}(x_{i+1/2,j})) \\ & + V_1(x_{i-1/2,j}, \rho_{i-1,j}^{n}, \rho_{i,j}^{n}) + F(\rho_{i-1,j}^{n}, \rho_{i,j}^{n}, J_1^{n}(x_{i-1/2,j})) \end{aligned} \right] \\ & \pm \lambda_x \left[ \begin{aligned} & V_1(x_{i+1/2,j+1}, \rho_{i,j}^{n}, \rho_{i+1,j}^{n}) + F(\rho_{i,j}^{n}, \rho_{i+1,j}^{n}, J_1^{n}(x_{i+1/2,j+1})) \\ & - V_1(x_{i-1/2,j+1}, \rho_{i-1,j}^{n}, \rho_{i,j}^{n}) - F(\rho_{i-1,j}^{n}, \rho_{i,j}^{n}, J_1^{n}(x_{i-1/2,j+1})) \end{aligned} \right] \\ &= \mathcal{D}_{i,j}^{n} + \lambda_x \mathcal{E}_{i,j}^{n}, \end{aligned}$$

where we set

$$\begin{aligned} \mathcal{D}_{i,j}^{n} &= \rho_{i,j+1}^{n} - \rho_{i,j}^{n} - \lambda_x \left[ \begin{aligned} & V_1(x_{i+1/2,j+1}, \rho_{i,j+1}^{n}, \rho_{i+1,j+1}^{n}) + F(\rho_{i,j+1}^{n}, \rho_{i+1,j+1}^{n}, J_1^{n}(x_{i+1/2,j+1})) \\ & - V_1(x_{i+1/2,j+1}, \rho_{i,j}^{n}, \rho_{i+1,j}^{n}) - F(\rho_{i,j}^{n}, \rho_{i+1,j}^{n}, J_1^{n}(x_{i+1/2,j+1})) \\ & + V_1(x_{i-1/2,j+1}, \rho_{i-1,j}^{n}, \rho_{i,j}^{n}) + F(\rho_{i-1,j}^{n}, \rho_{i,j}^{n}, J_1^{n}(x_{i-1/2,j+1})) \end{aligned} \right] \end{aligned}$$
$$-V_1(x_{i-1/2,j+1}, \rho_{i-1,j+1}^n, \rho_{i,j+1}^n) - F(\rho_{i-1,j+1}^n, \rho_{i,j+1}^n, J_1^n(x_{i-1/2,j+1}))],$$

$$\begin{aligned}
\mathcal{E}_{i,j}^n = & V_1(x_{i+1/2,j}, \rho_{i,j}^n, \rho_{i+1,j}^n) + F(\rho_{i,j}^n, \rho_{i+1,j}^n, J_1^n(x_{i+1/2,j})) \\
& - V_1(x_{i+1/2,j+1}, \rho_{i,j}^n, \rho_{i+1,j}^n) - F(\rho_{i,j}^n, \rho_{i+1,j}^n, J_1^n(x_{i+1/2,j+1})) \\
& + V_1(x_{i-1/2,j+1}, \rho_{i-1,j}^n, \rho_{i,j}^n) + F(\rho_{i-1,j}^n, \rho_{i,j}^n, J_1^n(x_{i-1/2,j+1})) \\
& - V_1(x_{i-1/2,j}, \rho_{i-1,j}^n, \rho_{i,j}^n) - F(\rho_{i-1,j}^n, \rho_{i,j}^n, J_1^n(x_{i-1/2,j})).
\end{aligned}$$

Similarly as before, rearrange $\mathcal{D}_{i,j}^n$, exploiting the notation (3.18):

$$\begin{aligned}
\mathcal{D}_{i,j}^n = & \rho_{i,j+1}^n - \rho_{i,j}^n - \lambda_x \left[ H_{i+1/2,j+1}^n(\rho_{i,j+1}^n, \rho_{i+1,j+1}^n) - H_{i+1/2,j+1}^n(\rho_{i,j}^n, \rho_{i+1,j}^n) \right. \\
& \left. \qquad + H_{i-1/2,j+1}^n(\rho_{i-1,j}^n, \rho_{i,j}^n) - H_{i-1/2,j+1}^n(\rho_{i-1,j+1}^n, \rho_{i,j+1}^n) \right] \\
& \pm \lambda_x H_{i+1/2,j+1}^n(\rho_{i,j}^n, \rho_{i+1,j+1}^n) \pm \lambda_x H_{i-1/2,j+1}^n(\rho_{i-1,j}^n, \rho_{i,j+1}^n) \\
= & \rho_{i,j+1}^n - \rho_{i,j}^n \\
& - \lambda_x \frac{H_{i+1/2,j+1}^n(\rho_{i,j+1}^n, \rho_{i+1,j+1}^n) - H_{i+1/2,j+1}^n(\rho_{i,j}^n, \rho_{i+1,j+1}^n)}{\rho_{i,j+1}^n - \rho_{i,j}^n} (\rho_{i,j+1}^n - \rho_{i,j}^n) \\
& - \lambda_x \frac{H_{i+1/2,j+1}^n(\rho_{i,j}^n, \rho_{i+1,j+1}^n) - H_{i+1/2,j+1}^n(\rho_{i,j}^n, \rho_{i+1,j}^n)}{\rho_{i+1,j+1}^n - \rho_{i+1,j}^n} (\rho_{i+1,j+1}^n - \rho_{i+1,j}^n) \\
& + \lambda_x \frac{H_{i-1/2,j+1}^n(\rho_{i-1,j}^n, \rho_{i,j+1}^n) - H_{i-1/2,j+1}^n(\rho_{i-1,j}^n, \rho_{i,j}^n)}{\rho_{i,j+1}^n - \rho_{i,j}^n} (\rho_{i,j+1}^n - \rho_{i,j}^n) \\
& + \lambda_x \frac{H_{i-1/2,j+1}^n(\rho_{i-1,j+1}^n, \rho_{i,j+1}^n) - H_{i-1/2,j+1}^n(\rho_{i-1,j}^n, \rho_{i,j+1}^n)}{\rho_{i-1,j+1}^n - \rho_{i-1,j}^n} (\rho_{i-1,j+1}^n - \rho_{i-1,j}^n) \\
= & (1 - \kappa_{i,j}^n - \nu_{i,j}^n)(\rho_{i,j+1}^n - \rho_{i,j}^n) + \nu_{i+1,j}^n(\rho_{i+1,j+1}^n - \rho_{i+1,j}^n) + \kappa_{i-1,j}^n(\rho_{i-1,j+1}^n - \rho_{i-1,j}^n),
\end{aligned}$$

where

$$\begin{aligned}
\kappa_{i,j}^n = & \begin{cases} \lambda_x \frac{H_{i+1/2,j+1}^n(\rho_{i,j+1}^n, \rho_{i+1,j+1}^n) - H_{i+1/2,j+1}^n(\rho_{i,j}^n, \rho_{i+1,j+1}^n)}{\rho_{i,j+1}^n - \rho_{i,j}^n} & \text{if } \rho_{i,j+1}^n \neq \rho_{i,j}^n, \\ 0 & \text{if } \rho_{i,j+1}^n = \rho_{i,j}^n, \end{cases} \\
\nu_{i,j}^n = & \begin{cases} -\lambda_x \frac{H_{i-1/2,j+1}^n(\rho_{i-1,j}^n, \rho_{i,j+1}^n) - H_{i-1/2,j+1}^n(\rho_{i-1,j}^n, \rho_{i,j}^n)}{\rho_{i,j+1}^n - \rho_{i,j}^n} & \text{if } \rho_{i,j+1}^n \neq \rho_{i,j}^n, \\ 0 & \text{if } \rho_{i,j+1}^n = \rho_{i,j}^n. \end{cases}
\end{aligned}$$

As for $\delta_i^n$ (3.19) and $\vartheta_i^n$ (3.20), it is immediate to prove that $\kappa_{i,j}^n, \nu_{i,j}^n \in [0, \frac{1}{3}]$ for all $i, j \in \mathbb{Z}$.
Hence,

$$\sum_{i,j \in \mathbb{Z}} |\mathcal{D}_{i,j}^n| \le \sum_{i,j \in \mathbb{Z}} |\rho_{i,j+1}^n - \rho_{i,j}^n|. \qquad (3.26)$$

Pass now to $\mathcal{E}_{i,j}^n$: we can proceed analogously to $\mathcal{B}_{i,j}^n$, treating separately the terms involving $V_1$ and those involving $F$. First, by (3.2),

$$V_1(x_{i+1/2,j}, \rho_{i,j}^n, \rho_{i+1,j}^n) - V_1(x_{i+1/2,j+1}, \rho_{i,j}^n, \rho_{i+1,j}^n)$$
$$
\begin{align*}
& + V_1(x_{i-1/2,j+1}, \rho_{i-1,j}^n, \rho_{i,j}^n) - V_1(x_{i-1/2,j}, \rho_{i-1,j}^n, \rho_{i,j}^n) \\
&= v_1^{\text{stat}}(x_{i+1/2,j})\rho_{i,j}^n + \min\{0, v_1^{\text{stat}}(x_{i+1/2,j})\} (\rho_{i+1,j}^n - \rho_{i,j}^n) \\
&\quad - v_1^{\text{stat}}(x_{i+1/2,j+1})\rho_{i,j}^n - \min\{0, v_1^{\text{stat}}(x_{i+1/2,j+1})\} (\rho_{i+1,j}^n - \rho_{i,j}^n) \\
&\quad + v_1^{\text{stat}}(x_{i-1/2,j+1})\rho_{i-1,j}^n + \min\{0, v_1^{\text{stat}}(x_{i-1/2,j+1})\} (\rho_{i,j}^n - \rho_{i-1,j}^n) \\
&\quad - v_1^{\text{stat}}(x_{i-1/2,j})\rho_{i-1,j}^n - \min\{0, v_1^{\text{stat}}(x_{i-1/2,j})\} (\rho_{i,j}^n - \rho_{i-1,j}^n) \\
&\quad \pm (v_1^{\text{stat}}(x_{i-1/2,j+1}) - v_1^{\text{stat}}(x_{i-1/2,j})) \rho_{i,j}^n \\
&= (v_1^{\text{stat}}(x_{i+1/2,j}) - v_1^{\text{stat}}(x_{i+1/2,j+1}) - v_1^{\text{stat}}(x_{i-1/2,j}) + v_1^{\text{stat}}(x_{i-1/2,j+1})) \rho_{i,j}^n \\
&\quad + (v_1^{\text{stat}}(x_{i-1/2,j+1}) - v_1^{\text{stat}}(x_{i-1/2,j})) (\rho_{i-1,j}^n - \rho_{i,j}^n) \\
&\quad + \left(\min\{0, v_1^{\text{stat}}(x_{i+1/2,j})\} - \min\{0, v_1^{\text{stat}}(x_{i+1/2,j+1})\}\right) (\rho_{i+1,j}^n - \rho_{i,j}^n) \\
&\quad + \left(\min\{0, v_1^{\text{stat}}(x_{i-1/2,j+1})\} - \min\{0, v_1^{\text{stat}}(x_{i-1/2,j})\}\right) (\rho_{i,j}^n - \rho_{i-1,j}^n) \\
&\le \Delta x \Delta y \left\|\partial_{xy}^2 v_1^{\text{stat}}\right\|_{\mathbf{L}^\infty} \left|\rho_{i,j}^n\right| + \Delta y \left\|\partial_y v_1^{\text{stat}}\right\|_{\mathbf{L}^\infty} \left(\left|\rho_{i+1,j}^n - \rho_{i,j}^n\right| + 2\left|\rho_{i,j}^n - \rho_{i-1,j}^n\right|\right), \quad (3.27)
\end{align*}
$$

since

$$
\begin{align*}
v_1^{\text{stat}}(x_{i+1/2,j}) - v_1^{\text{stat}}(x_{i+1/2,j+1}) - v_1^{\text{stat}}(x_{i-1/2,j}) + v_1^{\text{stat}}(x_{i-1/2,j+1}) \\
&= \Delta x \partial_x v_1^{\text{stat}}(\xi_i, y_j) - \Delta x \partial_x v_1^{\text{stat}}(\xi_i, y_{j+1}) \\
&= - \Delta x \Delta y \partial_{xy}^2 v_1^{\text{stat}}(\xi_i, \zeta_{j+1/2}),
\end{align*}
$$

with $\xi_i \in ]x_{i-1/2}, x_{i+1/2}]$ and $\zeta_{j+1/2} \in ]y_j, y_{j+1}]$. In a similar way, by (3.3),

$$
\begin{align*}
& F(\rho_{i,j}^n, \rho_{i+1,j}^n, J_1^n(x_{i+1/2,j})) - F(\rho_{i,j}^n, \rho_{i+1,j}^n, J_1^n(x_{i+1/2,j+1})) \\
&\quad + F(\rho_{i-1,j}^n, \rho_{i,j}^n, J_1^n(x_{i-1/2,j+1})) - F(\rho_{i-1,j}^n, \rho_{i,j}^n, J_1^n(x_{i-1/2,j})) \\
&= J_1^n(x_{i+1/2,j}) f(\rho_{i,j}^n) + \min\{0, J_1^n(x_{i+1/2,j})\} (f(\rho_{i+1,j}^n) - f(\rho_{i,j}^n)) \\
&\quad - J_1^n(x_{i+1/2,j+1}) f(\rho_{i,j}^n) - \min\{0, J_1^n(x_{i+1/2,j+1})\} (f(\rho_{i+1,j}^n) - f(\rho_{i,j}^n)) \\
&\quad + J_1^n(x_{i-1/2,j+1}) f(\rho_{i-1,j}^n) + \min\{0, J_1^n(x_{i-1/2,j+1})\} (f(\rho_{i,j}^n) - f(\rho_{i-1,j}^n)) \\
&\quad - J_1^n(x_{i-1/2,j}) f(\rho_{i-1,j}^n) - \min\{0, J_1^n(x_{i-1/2,j})\} (f(\rho_{i,j}^n) - f(\rho_{i-1,j}^n)) \\
&\quad \pm (J_1^n(x_{i-1/2,j+1}) - J_1^n(x_{i-1/2,j})) f(\rho_{i,j}^n) \\
&= (J_1^n(x_{i+1/2,j}) - J_1^n(x_{i+1/2,j+1}) - J_1^n(x_{i-1/2,j}) + J_1^n(x_{i-1/2,j+1})) f(\rho_{i,j}^n) \\
&\quad + (J_1^n(x_{i-1/2,j+1}) - J_1^n(x_{i-1/2,j})) (f(\rho_{i-1,j}^n) - f(\rho_{i,j}^n)) \\
&\quad + \left(\min\{0, J_1^n(x_{i+1/2,j})\} - \min\{0, J_1^n(x_{i+1/2,j+1})\}\right) (f(\rho_{i+1,j}^n) - f(\rho_{i,j}^n)) \\
&\quad + \left(\min\{0, J_1^n(x_{i-1/2,j+1})\} - \min\{0, J_1^n(x_{i-1/2,j})\}\right) (f(\rho_{i,j}^n) - f(\rho_{i-1,j}^n))
\end{align*}
$$
$$ \leq 2 \varepsilon \Delta x \Delta y C |\rho_i^n| + 2 \varepsilon L_f \Delta y \|\nabla^2 \eta\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \left( |\rho_{i+1,j}^n - \rho_{i,j}^n| + 2|\rho_{i,j}^n + \rho_{i-1,j}^n| \right), \quad (3.28) $$

where we used the fact that $f(r) \leq r$, (A.3) and (A.5), with the notation (A.6). Therefore, collecting together (3.27) and (3.28), we get

$$ \begin{aligned} |\mathcal{E}_{i,j}^n| & \leq \Delta x \Delta y \left( \left\| \partial_{xy} v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2 \varepsilon C \right) |\rho_i^n| \\ & \quad + \Delta y \left( \left\| \partial_y v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2 \varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \left( |\rho_{i+1,j}^n - \rho_{i,j}^n| + 2|\rho_{i,j}^n + \rho_{i-1,j}^n| \right), \end{aligned} $$

so that

$$ \sum_{i,j \in \mathbb{Z}} \lambda_x |\mathcal{E}_{i,j}^n| \leq 3 \lambda_x \Delta y \left( \left\| \partial_y v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2 \varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \sum_{i,j \in \mathbb{Z}} |\rho_{i+1,j}^n - \rho_{i,j}^n| \\ + \Delta t \left( \left\| \partial_{xy} v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2 \varepsilon C \right) \Delta y \sum_{i,j \in \mathbb{Z}} \rho_{i,j}^n. \quad (3.29) $$

Hence, by (3.26) and (3.29), using also Lemma 3.3, we obtain

$$ \begin{aligned} & \sum_{i,j \in \mathbb{Z}} \Delta x |\rho_{i,j+1}^{n+1/2} - \rho_{i,j}^{n+1/2}| \\ & \leq \sum_{i,j \in \mathbb{Z}} \Delta x \left( \left| \mathcal{D}_{i,j}^n \right| + \lambda_x |\mathcal{E}_{i,j}^n| \right) \\ & \leq \sum_{i,j \in \mathbb{Z}} \Delta x |\rho_{i,j+1}^n - \rho_{i,j}^n| \\ & \quad + 3 \Delta t \left( \left\| \partial_y v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2 \varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \sum_{i,j \in \mathbb{Z}} \Delta y |\rho_{i+1,j}^n - \rho_{i,j}^n| \quad (3.30) \\ & \quad + \Delta t \left( \left\| \partial_{xy} v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 2 \varepsilon C \right) \|\rho_o\|_{\mathbf{L}^1}. \end{aligned} $$

Setting

$$ K_1 = 3 \left( \left\| \partial_x v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + \left\| \partial_y v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + 4 \varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right), \quad (3.31) $$

$$ K_2 = \left( 4 \varepsilon C + 2 \left\| \partial_{xx}^2 v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} + \left\| \partial_{xy}^2 v_1^{\text{stat}} \right\|_{\mathbf{L}^\infty} \right) \|\rho_o\|_{\mathbf{L}^1}, \quad (3.32) $$

by (3.25) and (3.30) we conclude

$$ \begin{aligned} & \sum_{i,j \in \mathbb{Z}} \left( \Delta y |\rho_{i+1,j}^{n+1/2} - \rho_{i,j}^{n+1/2}| + \Delta x |\rho_{i,j+1}^{n+1/2} - \rho_{i,j}^{n+1/2}| \right) \\ & \leq (1 + \Delta t K_1) \sum_{i,j \in \mathbb{Z}} \left( \Delta x |\rho_{i,j+1}^n - \rho_{i,j}^n| + \Delta y |\rho_{i+1,j}^n - \rho_{i,j}^n| \right) + \Delta t K_2. \end{aligned} $$
Analogous computations yield

$$ \sum_{i,j \in \mathbb{Z}} \left( \Delta y \left| \rho_{i+1,j}^{n+1} - \rho_{i,j}^{n+1} \right| + \Delta x \left| \rho_{i,j+1}^{n+1} - \rho_{i,j}^{n+1} \right| \right) \\ \leq (1 + \Delta t K_3) \sum_{i,j \in \mathbb{Z}} \left( \Delta x \left| \rho_{i,j+1}^{n+1/2} - \rho_{i,j}^{n+1/2} \right| + \Delta y \left| \rho_{i+1,j}^{n+1/2} - \rho_{i,j}^{n+1/2} \right| \right) + \Delta t K_4, $$

where

$$ K_3 = 3 \left( \left\| \partial_x v_2^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + \left\| \partial_y v_2^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + 4 \varepsilon L_f \left\| \nabla^2 \eta \right\|_{\mathbf{L}^{\infty}} \| \rho_o \|_{\mathbf{L}^1} \right), \quad (3.33) $$

$$ K_4 = \left( 4 \varepsilon C + 2 \left\| \partial_{yy} v_2^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} + \left\| \partial_{xy}^2 v_2^{\text{stat}} \right\|_{\mathbf{L}^{\infty}} \right) \| \rho_o \|_{\mathbf{L}^1}. \quad (3.34) $$

Observe that, using the notation (3.16) and (3.17),

$$ K_1, K_3 \leq K_1, \quad K_2, K_4 \leq K_2. $$

A recursive argument yields the desired result:

$$ \sum_{i,j \in \mathbb{Z}} \left( \Delta y \left| \rho_{i+1,j}^n - \rho_{i,j}^n \right| + \Delta x \left| \rho_{i,j+1}^n - \rho_{i,j}^n \right| \right) \\ \leq e^{2n \Delta t K_1} \sum_{i,j \in \mathbb{Z}} \left( \Delta x \left| \rho_{i,j+1}^0 - \rho_{i,j}^0 \right| + \Delta y \left| \rho_{i+1,j}^0 - \rho_{i,j}^0 \right| \right) + \frac{2 K_2}{K_1} (e^{2n \Delta t K_1} - 1). $$

□

**Corollary 3.7. (BV estimate in space and time)** Let $\rho_o \in (\mathbf{L}^\infty \cap \mathbf{BV})(\mathbb{R}^2; \mathbb{R}^+)$, Let $(v)$, $(H)$, $(\eta)$, (3.13) hold. Then, for all $t > 0$, $\rho_\Delta$ in (3.1) constructed through Algorithm 3.1 satisfies the following estimate: for all $n = 1, \dots, N_T$,

$$ \sum_{m=0}^{n-1} \sum_{i,j \in \mathbb{Z}} \Delta t \left( \Delta y \left| \rho_{i+1,j}^m - \rho_{i,j}^m \right| + \Delta x \left| \rho_{i,j+1}^m - \rho_{i,j}^m \right| \right) + \sum_{m=0}^{n-1} \sum_{i,j \in \mathbb{Z}} \Delta x \Delta y \left| \rho_{i,j}^{m+1} - \rho_{i,j}^m \right| \leq C_{xt}(t^n), \quad (3.35) $$

where

$$ C_{xt}(t) = t(C_x(t) + 2C_t(t)), \quad (3.36) $$

with $C_x$ as in (3.15) and $C_t$ as in (3.38).

**Proof.** By Proposition 3.5 we have

$$ \sum_{m=0}^{n-1} \sum_{i,j \in \mathbb{Z}} \Delta t \left( \Delta y \left| \rho_{i+1,j}^m - \rho_{i,j}^m \right| + \Delta x \left| \rho_{i,j+1}^m - \rho_{i,j}^m \right| \right) \leq n \Delta t C_x(n \Delta t). \quad (3.37) $$

Since

$$ \left| \rho_{i,j}^{m+1} - \rho_{i,j}^{m} \right| \leq \left| \rho_{i,j}^{m+1} - \rho_{i+1,j}^{m+1/2} \right| + \left| \rho_{i,j}^{m+1/2} - \rho_{i,j}^{m} \right|, $$
we focus first on

$$ \sum_{i,j \in \mathbb{Z}} \Delta x \Delta y |\rho_{i,j}^{m+1/2} - \rho_{i,j}^{m}|. $$

By the scheme (3.5), we have, using the notation (3.8),

$$
\begin{align*}
\rho_{i,j}^{m+1/2} - \rho_{i,j}^{m} &\le \lambda_x \left[ V_1(x_{i+1/2,j}, \rho_{i,j}^{m}, \rho_{i+1,j}^{m}) + F(\rho_{i,j}^{m}, \rho_{i+1,j}^{m}, J_1^m(x_{i+1/2,j})) \right. \\
&\qquad \left. - V_1(x_{i-1/2,j}, \rho_{i-1,j}^{m}, \rho_{i,j}^{m}) - F(\rho_{i-1,j}^{m}, \rho_{i,j}^{m}, J_1^m(x_{i-1/2,j})) \right] \\
&= \lambda_x \left[ v_{i+1/2,j} \rho_{i,j}^{m} + \min\{0, v_{i+1/2,j}\} (\rho_{i+1,j}^{m} - \rho_{i,j}^{m}) \right. \\
&\qquad \left. - v_{i-1/2,j} \rho_{i-1,j}^{m} - \min\{0, v_{i-1/2,j}\} (\rho_{i,j}^{m} - \rho_{i-1,j}^{m}) \right. \\
&\qquad \left. + J_1^m(x_{i+1/2,j}) f(\rho_{i,j}^{m}) + \min\{0, J_1^m(x_{i+1/2,j})\} (f(\rho_{i+1,j}^{m}) - f(\rho_{i,j}^{m})) \right. \\
&\qquad \left. - J_1^m(x_{i-1/2,j}) f(\rho_{i-1,j}^{m}) - \min\{0, J_1^m(x_{i-1/2,j})\} (f(\rho_{i,j}^{m}) - f(\rho_{i-1,j}^{m})) \right. \\
&\qquad \left. \pm v_{i-1/2,j} \rho_{i,j}^{m} \pm J_1^m(x_{i-1/2,j}) f(\rho_{i,j}^{m}) \right] \\
&= \lambda_x \left[ \Delta x \partial_x v_1^{\text{stat}}(\xi_i, y_j) \rho_{i,j}^{m} + \left( v_{i-1/2,j} - \min\{0, v_{i-1/2,j}\} \right) (\rho_{i,j}^{m} - \rho_{i-1,j}^{m}) \right. \\
&\qquad \left. + \min\{0, v_{i+1/2,j}\} (\rho_{i+1,j}^{m} - \rho_{i,j}^{m}) + \left( J_1^m(x_{i+1/2,j}) - J_1^m(x_{i-1/2,j}) \right) f(\rho_{i,j}^{m}) \right. \\
&\qquad \left. + \left( J_1^m(x_{i-1/2,j}) - \min\{0, J_1^m(x_{i-1/2,j})\} \right) (f(\rho_{i,j}^{m}) - f(\rho_{i-1,j}^{m})) \right. \\
&\qquad \left. + \min\{0, J_1^m(x_{i+1/2,j})\} (f(\rho_{i+1,j}^{m}) - f(\rho_{i,j}^{m})) \right] \\
&\le \Delta t \left( \|\partial_x v_1^{\text{stat}}\|_{\mathbf{L}^\infty} + 2\varepsilon \|\nabla^2 \eta\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \rho_{i,j}^{m} \\
&\quad + \lambda_x \left( \|v_1^{\text{stat}}\|_{\mathbf{L}^\infty} + \varepsilon L_f \right) \left( |\rho_{i,j}^{m} - \rho_{i-1,j}^{m}| + |\rho_{i+1,j}^{m} - \rho_{i,j}^{m}| \right),
\end{align*}
$$

where $\xi_i \in ]x_{i-1/2}, x_{i+1/2}[ \text{ and we used } f(r) \le r, \text{(A.1) and (A.2). Therefore,}$

$$
\begin{align*}
\sum_{i,j \in \mathbb{Z}} \Delta x \Delta y |\rho_{i,j}^{m+1/2} - \rho_{i,j}^{m}| &\le 2 \Delta t \left( \|v_1^{\text{stat}}\|_{\mathbf{L}^\infty} + \varepsilon L_f \right) \sum_{i,j \in \mathbb{Z}} \Delta y |\rho_{i+1,j}^{m} - \rho_{i,j}^{m}| \\
&\quad + \Delta t \left( \|\partial_x v_1^{\text{stat}}\|_{\mathbf{L}^\infty} + 2\varepsilon \|\nabla^2 \eta\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \|\rho_o\|_{\mathbf{L}^1} \\
&\le 2 \Delta t \left( \|v_1^{\text{stat}}\|_{\mathbf{L}^\infty} + \varepsilon L_f \right) C_x(m \Delta t) \\
&\quad + \Delta t \left( \|\partial_x v_1^{\text{stat}}\|_{\mathbf{L}^\infty} + 2\varepsilon \|\nabla^2 \eta\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \|\rho_o\|_{\mathbf{L}^1} \\
&\le \Delta t C_t(m \Delta t),
\end{align*}
$$

where we set

$$ C_t(s) = 2 \left( \|v^{\text{stat}}\|_{\mathbf{L}^\infty} + \varepsilon L_f \right) C_x(s) + \left( \|\nabla v^{\text{stat}}\|_{\mathbf{L}^\infty} + 2\varepsilon \|\nabla^2 \eta\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1} \right) \|\rho_o\|_{\mathbf{L}^1}. \quad (3.38) $$
Analogously, we get

$$
\begin{align*}
\sum_{i,j \in \mathbb{Z}} \Delta x \Delta y |\rho_{i,j}^{m+1} - \rho_{i,j}^{m+1/2}| &\le 2 \Delta t \left( \|v_2^{\text{stat}}\|_{\mathbf{L}^{\infty}} + \varepsilon L_f \right) \sum_{i,j \in \mathbb{Z}} \Delta x |\rho_{i,j+1}^{m+1/2} - \rho_{i,j}^{m+1/2}| \\
&\quad + \Delta t \left( \|\partial_y v_2^{\text{stat}}\|_{\mathbf{L}^{\infty}} + 2\varepsilon \|\nabla^2 \eta\|_{\mathbf{L}^{\infty}} \|\rho_o\|_{\mathbf{L}^1} \right) \|\rho_o\|_{\mathbf{L}^1} \\
&\le 2 \Delta t \left( \|v_2^{\text{stat}}\|_{\mathbf{L}^{\infty}} + \varepsilon L_f \right) C_x(m \Delta t) \\
&\quad + \Delta t \left( \|\partial_y v_2^{\text{stat}}\|_{\mathbf{L}^{\infty}} + 2\varepsilon \|\nabla^2 \eta\|_{\mathbf{L}^{\infty}} \|\rho_o\|_{\mathbf{L}^1} \right) \|\rho_o\|_{\mathbf{L}^1} \\
&\le \Delta t C_t(m \Delta t).
\end{align*}
$$

Hence

$$
\sum_{m=0}^{n-1} \sum_{i,j \in \mathbb{Z}} \Delta x \Delta y |\rho_{i,j}^{m+1} - \rho_{i,j}^{m}| \le 2 n \Delta t C_t(n \Delta t), \quad (3.39)
$$

which together with (3.37) completes the proof.

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

3.5 Discrete entropy inequalities

Following [1], see also [7, 8], introduce the following notation: for $i, j \in \mathbb{Z}$, $n = 0, \dots, N_T - 1$ and $\kappa \in \mathbb{R}$,

$$
\Phi_{i+1/2,j}^n(u,v) = V_1(x_{i+1/2,j}, u \lor \kappa, v \lor \kappa) + F(u \lor \kappa, v \lor \kappa, J_1^n(x_{i+1/2,j})) \\
\qquad - V_1(x_{i+1/2,j}, u \land \kappa, v \land \kappa) - F(u \land \kappa, v \land \kappa, J_1^n(x_{i+1/2,j})),
$$

$$
\Gamma_{i,j+1/2}^{n}(u,v) = V_2(x_{i,j+1/2}, u \lor \kappa, v \lor \kappa) + F(u \lor \kappa, v \lor \kappa, J_2^n(x_{i,j+1/2})) \\
\qquad - V_2(x_{i,j+1/2}, u \land \kappa, v \land \kappa) - F(u \land \kappa, v \land \kappa, J_2^n(x_{i,j+1/2})),
$$

with $V_1$, $V_2$ and $F$ defined as in (3.2), (3.4) and (3.3) respectively.

**Lemma 3.8. (Discrete entropy condition)** Fix $\rho_o \in (\mathbf{L}^\infty \cap \mathbf{BV})(\mathbb{R}^2; \mathbb{R}^+)$.
Let $(v)$, $(H)$,
$(\eta)$, (3.13) hold. Then, the solution $\rho_\Delta$ in (3.1) constructed through Algorithm 3.1 satisfies
the following discrete entropy inequality: for $i, j \in \mathbb{Z}$, for $n = 0, \dots, N_T - 1$ and $\kappa \in \mathbb{R}$,

$$
\begin{align*}
\left| \rho_{i,j}^{n+1} - \kappa \right| - \left| \rho_{i,j}^n - \kappa \right| &+ \lambda_x \left( \Phi_{i+1/2,j}^n(\rho_{i,j}^n, \rho_{i+1,j}^n) - \Phi_{i-1/2,j}^n(\rho_{i-1,j}^n, \rho_{i,j}^n) \right) \\
&+ \lambda_x \operatorname{sgn}(\rho_{i,j}^{n+1/2} - \kappa) \left( v_1^{\text{stat}}(x_{i+1/2,j}) - v_1^{\text{stat}}(x_{i-1/2,j}) \right) \kappa \\
&+ \lambda_x \operatorname{sgn}(\rho_{i,j}^{n+1/2} - \kappa) \left( J_1^n(x_{i+1/2,j}) - J_1^n(x_{i-1/2,j}) \right) f(\kappa) \\
&+ \lambda_y \left( \Gamma_{i,j+1/2}^n(\rho_{i,j}^{n+1/2}, \rho_{i,j+1}^{n+1/2}) - \Gamma_{i,j-1/2}^n(\rho_{i,j-1}^{n+1/2}, \rho_{i,j}^{n+1/2}) \right) \\
&+ \lambda_y \operatorname{sgn}(\rho_{i,j}^{n+1} - \kappa) \left( v_2^{\text{stat}}(x_{i,j+1/2}) - v_2^{\text{stat}}(x_{i,j-1/2}) \right) \kappa \\
&+ \lambda_y \operatorname{sgn}(\rho_{i,j}^{n+1} - \kappa) \left( J_2^n(x_{i,j+1/2}) - J_2^n(x_{i,j-1/2}) \right) f(\kappa) \le 0.
\end{align*}
$$

The proof is omitted, being entirely analogous to that of [2, Proposition 2.8], see also [1, Lemma 2.8].
## 4 Lipschitz continuous dependence on initial data

**Proposition 4.1.** Fix $T > 0$. Let $(v)$, $(H)$ and $(\eta)$ hold. Let $\rho_o$, $\sigma_o \in (\mathbf{L}^\infty \cap \mathbf{BV})(\mathbb{R}^2; \mathbb{R}^+)$.
Call $\rho$ and $\sigma$ the corresponding solutions to (1.1). Then the following estimate holds:

$$ \|\rho(t) - \sigma(t)\|_{\mathbf{L}^1(\mathbb{R}^2)} \le \|\rho_o - \sigma_o\|_{\mathbf{L}^1(\mathbb{R}^2)} e^{tA(t)}. $$

with $A(t)$ defined in (4.2).

**Proof.** In the rest of the proof, to avoid heavy notation, we will denote pairs in $\mathbb{R}^2$ by $x$ or $y$. Introduce the following notation:

$$ \mathbf{R}(t, x) = (\mathbf{I}(\rho(t)))(x), \qquad \mathbf{S}(t, x) = (\mathbf{I}(\sigma(t)))(x). \qquad (4.1) $$

The idea is to apply the *doubling of variables method* introduced by Kružkov in [11], exploiting in particular the proof of [4, Lemma 4]. There, a flux of the form $f(t, x, \rho) V(t, x)$ is taken into account, with $x \in \mathbb{R}$, the proof being valid also in the multidimensional case, i.e. $x \in \mathbb{R}^n$. Therefore we are going to use this result for what concerns the part of the flux of type $f(\rho) \mathbf{R}(t, x)$.

For the sake of completeness, we recall that a flux function of type $l(x) g(\rho)$ is considered in [10], with $x \in \mathbb{R}^n$, and the proof of [4, Lemma 4] follows the lines of that of [10, Theorem 1.3]. Thus, here we are adding the dependence on time to the function $l(x)$ considered in [10].

Let $\varphi \in \mathbf{C}_c^1[0, T[\times \mathbb{R}^2; \mathbb{R}^+)$ be a test function as in the definition of solution by Kružkov. Let $Y \in \mathbf{C}_c^\infty(\mathbb{R}; \mathbb{R}^+)$ be such that

$$ Y(z) = Y(-z), \qquad Y(z) = 0 \text{ for } |z| \ge 1, \qquad \int_{\mathbb{R}} Y(z) dz = 1. $$

Define, for $h > 0$, $Y_h(z) = \frac{1}{h}Y(\frac{z}{h})$. Clearly, $Y_h \in \mathbf{C}_c^\infty(\mathbb{R}; \mathbb{R}^+)$, $Y_h(z) = Y_h(-z)$, $Y_h(z) = 0$ for $|z| \ge h$, $\int_{\mathbb{R}} Y_h(z) dz = 1$ and $Y_h \to \delta_0$ as $h \to 0$, $\delta_0$ being the Dirac delta in 0. Define moreover

$$ \psi_h(t, \mathbf{x}, s, \mathbf{y}) = \varphi \left( \frac{t+s}{2}, \frac{\mathbf{x}+\mathbf{y}}{2} \right) Y_h(t-s) \prod_{i=1}^{2} Y_h(x_i - y_i). $$

Introduce the space $\Pi_T = ]0, T[\times \mathbb{R}^2$ and, from the definition of solution, derive the following entropy inequalities for $\rho = \rho(t, \mathbf{x})$ and $\sigma = \sigma(s, \mathbf{y})$:

$$ \iiint_{\Pi_T \times \Pi_T} \left\{ \begin{aligned} & |\rho - \sigma| \partial_t \psi_h + |\rho - \sigma| \nabla_x \mathbf{v}^{\text{stat}}(\mathbf{x}) \cdot \nabla_x \psi_h + \text{sgn}(\rho - \sigma) (f(\rho) - f(\sigma)) \mathbf{R}(t, \mathbf{x}) \cdot \nabla_x \psi_h \\ & - \text{sgn}(\rho - \sigma) \sigma \text{div}_x \mathbf{v}^{\text{stat}}(\mathbf{x}) \psi_h - \text{sgn}(\rho - \sigma) f(\sigma) \text{div}_x \mathbf{R}(t, \mathbf{x}) \psi_h \end{aligned} \right\} d\mathbf{x} dt d\mathbf{y} ds \ge 0, $$

$$ \iiint_{\Pi_T \times \Pi_T} \left\{ \begin{aligned} & |\sigma - \rho| \partial_s \psi_h + |\sigma - \rho| \nabla_y \mathbf{v}^{\text{stat}}(\mathbf{y}) \cdot \nabla_y \psi_h + \text{sgn}(\sigma - \rho) (f(\sigma) - f(\rho)) \mathbf{S}(s, \mathbf{y}) \cdot \nabla_y \psi_h \\ & - \text{sgn}(\sigma - \rho) \rho \text{div}_y \mathbf{v}^{\text{stat}}(\mathbf{y}) \psi_h - \text{sgn}(\sigma - \rho) f(\rho) \text{div}_y \mathbf{S}(s, \mathbf{y}) \psi_h \end{aligned} \right\} d\mathbf{x} dt d\mathbf{y} ds \ge 0. $$

Sum the two inequalities above and rearrange the terms therein, following the proof of [11, Theorem 1] for what concerns the linear part of the flux and the proof of [4, Lemma 4], see also [10, Theorem 1.3], for the other part:

$$ \iiint_{\Pi_T \times \Pi_T} \left\{ |\rho - \sigma|(\partial_t \psi_h + \partial_s \psi_h) + \text{sgn}(\rho - \sigma)(\rho \mathbf{v}^{\text{stat}}(\mathbf{x}) - \sigma \mathbf{v}^{\text{stat}}(\mathbf{y})) \cdot (\nabla_x \psi_h + \nabla_y \psi_h) \right\} d\mathbf{x} dt d\mathbf{y} ds \ge 0. $$
$$
\begin{aligned}
& + \operatorname{sgn}(\rho - \sigma) \sigma \left[ (\mathbf{v}^{\text{stat}}(\mathbf{y}) - \mathbf{v}^{\text{stat}}(\mathbf{x})) \cdot \nabla_x \psi_h - \operatorname{div}_x \mathbf{v}^{\text{stat}}(\mathbf{x}) \psi_h \right] \\
& + \operatorname{sgn}(\rho - \sigma) \rho \left[ (\mathbf{v}^{\text{stat}}(\mathbf{y}) - \mathbf{v}^{\text{stat}}(\mathbf{x})) \cdot \nabla_y \psi_h + \operatorname{div}_y \mathbf{v}^{\text{stat}}(\mathbf{y}) \psi_h \right] \\
& + \operatorname{sgn}(\rho - \sigma) (f(\rho) \mathbf{R}(t, \mathbf{x}) - f(\sigma) \mathbf{S}(s, \mathbf{y})) \cdot (\nabla_x \psi_h + \nabla_y \psi_h) \\
& + \operatorname{sgn}(\rho - \sigma) f(\sigma) \left[ (\mathbf{S}(s, \mathbf{y}) - \mathbf{R}(t, \mathbf{x})) \cdot \nabla_x \psi - \operatorname{div}_x \mathbf{R}(t, \mathbf{x}) \psi_h \right] \\
& + \operatorname{sgn}(\rho - \sigma) f(\rho) \left[ (\mathbf{S}(s, \mathbf{y}) - \mathbf{R}(t, \mathbf{x})) \cdot \nabla_y \psi + \operatorname{div}_y \mathbf{S}(s, \mathbf{y}) \psi_h \right] \right\} d\mathbf{x} \, dt \, d\mathbf{y} \, ds \geq 0.
\end{aligned}
$$

Let $h \to 0$, which gives

$$
\iint_{\Pi_T} \left\{ \begin{aligned}
& |\rho - \sigma| \partial_t \varphi_h + \operatorname{sgn}(\rho - \sigma) \left[ (\rho - \sigma) \mathbf{v}^{\text{stat}}(\mathbf{x}) + \mathbf{S}(t, \mathbf{x}) (f(\rho) - f(\sigma)) \right] \cdot \nabla_x \varphi \\
& \qquad + \operatorname{sgn}(\rho - \sigma) f(\rho) \operatorname{div}_x (\mathbf{S}(t, \mathbf{x}) - \mathbf{R}(t, \mathbf{x})) \varphi \\
& \qquad + \operatorname{sgn}(\rho - \sigma) f'(\rho) (\mathbf{S}(t, \mathbf{x}) - \mathbf{R}(t, \mathbf{x})) \partial_x \rho(t, \mathbf{x}) \varphi
\end{aligned} \right\} d\mathbf{x} \, dt \, d\mathbf{y} \, ds \ge 0.
$$

Choosing a suitable test function $\varphi$ leads to

$$
\int_{\mathbb{R}^2} |\rho(t, \mathbf{x}) - \sigma(t, \mathbf{x})| d\mathbf{x} - \int_{\mathbb{R}^2} |\rho(\tau, \mathbf{x}) - \sigma(\tau, \mathbf{x})| d\mathbf{x} \\
+ \int_{\tau}^{t} \int_{\mathbb{R}^2} \left| \operatorname{div} (\mathbf{S}(s, \mathbf{x}) - \mathbf{R}(s, \mathbf{x})) \right| f(\rho(s, \mathbf{x})) d\mathbf{x} \, ds \\
+ \int_{\tau}^{t} L_f \left\| \mathbf{S}(s) - \mathbf{R}(s) \right\|_{\mathbf{L}^{\infty}(\mathbb{R}^2)} \operatorname{TV}(\rho(s)) \, ds \geq 0.
$$

Observe that, following [6, Lemma 4.1], the following bounds hold

$$
\begin{aligned}
\| \mathbf{S}(s) - \mathbf{R}(s) \|_{\mathbf{L}^{\infty}(\mathbb{R}^2)} &\le 2\varepsilon \| \nabla \eta \|_{\mathbf{L}^{\infty}} \| \rho(s) - \sigma(s) \|_{\mathbf{L}^1(\mathbb{R}^2)}, \\
\| \operatorname{div} (\mathbf{S}(s) - \mathbf{R}(s)) \|_{\mathbf{L}^{\infty}(\mathbb{R}^2)} &\le \varepsilon \| \rho(s) - \sigma(s) \|_{\mathbf{L}^1(\mathbb{R}^2)} \| \Delta \eta \|_{\mathbf{L}^{\infty}} (1 + \| \sigma(s) \|_{\mathbf{L}^1(\mathbb{R}^2)} \| \nabla \eta \|_{\mathbf{L}^1}).
\end{aligned}
$$

Thus, letting $\tau \to 0$ and exploiting the bounds on $\rho$ and $\sigma$ given by Theorem 2.2, as well as $f(r) \le r$, we get

$$
\begin{aligned}
& \int_{\mathbb{R}^2} |\rho(t, \mathbf{x}) - \sigma(t, \mathbf{x})| d\mathbf{x} \\
\le & \int_{\mathbb{R}^2} |\rho_o(\mathbf{x}) - \sigma_o(\mathbf{x})| d\mathbf{x} + 2\varepsilon \|\nabla\eta\|_{\mathbf{L}^\infty} \|\rho_o\|_{\mathbf{L}^1(\mathbb{R}^2)} \int_0^t \int_{\mathbb{R}^2} |\rho(s, \mathbf{x}) - \sigma(s, \mathbf{x})| d\mathbf{x} \, ds \\
& + \varepsilon L_f \|\Delta\eta\|_{\mathbf{L}^\infty} (1 + \|\sigma_o\|_{\mathbf{L}^1(\mathbb{R}^2)} \|\nabla\eta\|_{\mathbf{L}^1}) \int_0^t \operatorname{TV}(\rho(s)) (\int_{\mathbb{R}^2} |\rho(s, \mathbf{x}) - \sigma(s, \mathbf{x})| d\mathbf{x}) \, ds \\
= & \int_{\mathbb{R}^2} |\rho_o(\mathbf{x}) - \sigma_o(\mathbf{x})| d\mathbf{x} + \int_0^t A(s) (\int_{\mathbb{R}^2} |\rho(s, \mathbf{x}) - \sigma(s, \mathbf{x})| d\mathbf{x}) \, ds,
\end{aligned}
$$

with

$$
A(s) = 2\varepsilon \|\nabla \eta\|_{\mathbf{L}^{\infty}} \|\rho_o\|_{\mathbf{L}^1(\mathbb{R}^2)} + \varepsilon L_f \|\Delta \eta\|_{\mathbf{L}^{\infty}} (1 + \|\sigma_o\|_{\mathbf{L}^1(\mathbb{R}^2)} \|\nabla \eta\|_{\mathbf{L}^1}) \operatorname{TV}(\rho(s)). \quad (4.2)
$$
An application of Gronwall Lemma, together with

$$ \int_{0}^{t} A(s) \exp \left( \int_{s}^{t} A(r) dr \right) ds = -1 + \exp \left( \int_{0}^{t} A(s) ds \right), $$

yields the desired estimate

$$ \|\rho(t) - \sigma(t)\|_{\mathbf{L}^1(\mathbb{R}^2)} \le \|\rho_o - \sigma_o\|_{\mathbf{L}^1(\mathbb{R}^2)} e^{tA(t)}. $$

□

**Remark 4.2.** We can interpret $\rho$ and $\sigma$ as solutions of the following Cauchy problems:

$$ \left\{ \begin{array}{l} \partial_t \rho + \nabla \cdot g(t, \mathbf{x}, \rho) = 0, \\ \rho(0, \mathbf{x}) = \rho_o(\mathbf{x}), \end{array} \right. $$

$$ \left\{ \begin{array}{l} \partial_t \sigma + \nabla \cdot h(t, \mathbf{x}, \sigma) = 0, \\ \sigma(0, \mathbf{x}) = \sigma_o(\mathbf{x}), \end{array} \right. \quad (t, \mathbf{x}) \in [0, T] \times \mathbb{R}^2, \quad \mathbf{x} \in \mathbb{R}^2, $$

where

$$ g(t, \mathbf{x}, r) = r v^{\text{stat}}(\mathbf{x}) + f(r) R(t, \mathbf{x}), \quad h(t, \mathbf{x}, r) = r v^{\text{stat}}(\mathbf{x}) + f(r) S(t, \mathbf{x}), $$

so that the $\mathbf{L}^1$ distance between the solutions at time $t > 0$ can be estimated by [12, Proposition 2.10], see also the refinement in [13, Proposition 2.9]. However, making use of the explicit expression of the flux in the present case, one may see that the bound provided by Proposition 4.1 is sharper than that coming from [13, Proposition 2.9].

# 5 Lax-Friedrichs scheme

It is also possible to consider a piece-wise constant solution $\rho_\Delta$ to (1.1) as in (3.1) defined through a Lax-Friedrichs type finite volume scheme with dimensional splitting. The algorithm reads as follows

**Algorithm 5.1.**

$$ \text{for } n = 0, \dots, N_T - 1 \quad F^n(x, y, u, w) = \frac{1}{2} \left[ v_1^{\text{stat}}(x, y)(u+w) + J_1^n(x, y) (f(u) + f(w)) \right] - \frac{\alpha}{2} (w-u) \quad (5.1) $$

$$ G^n(x, y, u, w) = \frac{1}{2} \left[ v_2^{\text{stat}}(x, y)(u+w) + J_2^n(x, y) (f(u) + f(w)) \right] - \frac{\beta}{2} (w-u) \quad (5.2) $$

$$ \rho_{i,j}^{n+1/2} = \rho_{i,j}^{n} - \lambda_x \left[ F^n(x_{i+1/2,j}, \rho_{i,j}^{n}, \rho_{i+1,j}^{n}) - F^n(x_{i-1/2,j}, \rho_{i-1,j}^{n}, \rho_{i,j}^{n}) \right] \quad (5.3) $$

$$ \rho_{i,j}^{n+1} = \rho_{i,j}^{n+1/2} - \lambda_y \left[ G^n(x_{i,j+1/2}, \rho_{i,j}^{n}, \rho_{i,j+1}^{n}) - G^n(x_{i,j-1/2}, \rho_{i,j-1}^{n}, \rho_{i,j}^{n}) \right] \quad (5.4) $$

end

The algorithm is close to that studied in [1], except that in the present case to compute $\rho^{n+1}$ the flux is evaluated at $\rho^n$, instead of $\rho^{n+1/2}$.

Following closely the proofs presented in [1], it is possible to recover also for Algorithm 5.1 the bounds on the approximate solution necessary to prove the convergence. Below, we state only the final results, omitting the computations.
**Lemma 5.2. (Positivity)** Let $\rho_o \in L^\infty(\mathbb{R}^2; \mathbb{R}^+)$. Let assumptions $(v)$, $(H)$ and $(\eta)$ hold. Assume that

$$ \alpha \ge \|v_1^{\text{stat}}\|_{L^\infty} + \varepsilon L_f, \qquad \lambda_x \le \frac{1}{3} \min \left\{ \frac{1}{\alpha}, \frac{1}{2\varepsilon L_f + \Delta x \|v_1^{\text{stat}}\|_{L^\infty}} \right\}, \qquad (5.5) $$

$$ \beta \ge \|v_2^{\text{stat}}\|_{L^\infty} + \varepsilon L_f, \qquad \lambda_y \le \frac{1}{3} \min \left\{ \frac{1}{\beta}, \frac{1}{2\varepsilon L_f + \Delta x \|v_2^{\text{stat}}\|_{L^\infty}} \right\}. \qquad (5.6) $$

Then, for all $t > 0$ and $(x, y) \in \mathbb{R}^2$, the piece-wise constant approximate solution $\rho_\Delta$ (3.1) constructed through Algorithm 5.1 is such that $\rho_\Delta(t, x, y) \ge 0$.

**Lemma 5.3. (L¹ bound)** Let $\rho_o \in L^\infty(\mathbb{R}^2; \mathbb{R}^+)$. Let $(v)$, $(H)$, $(\eta)$, (5.5) and (5.6) hold. Then, for all $t > 0$, $\rho_\Delta$ in (3.1) constructed through Algorithm 5.1 satisfies (3.10), that is

$$ \|\rho_\Delta(t)\|_{L^1(\mathbb{R}^2)} = \|\rho_o\|_{L^1(\mathbb{R}^2)}. $$

**Lemma 5.4. (L^∞ bound)** Let $\rho_o \in L^\infty(\mathbb{R}^2; \mathbb{R}^+)$. Let $(v)$, $(H)$, $(\eta)$, (5.5) and (5.6) hold. Then, for all $t > 0$, $\rho_\Delta$ in (3.1) constructed through Algorithm 5.1 satisfies

$$ \|\rho_\Delta(t)\|_{L^\infty(\mathbb{R}^2)} \le \|\rho_o\|_{L^\infty} e^{\tilde{C}_\infty t}, $$

where

$$ \tilde{C}_\infty = \|\partial_x v_1^{\text{stat}}\|_{L^\infty} + \|\partial_y v_2^{\text{stat}}\|_{L^\infty} + 4\varepsilon L_f \|\nabla^2 \eta\|_{L^\infty} \|\rho_o\|_{L^1}. $$

**Remark 5.5.** Compare the $L^\infty$ estimate obtained in Lemma 5.4 using the Lax-Friedrichs scheme with that in Lemma 3.4, given by the Roe scheme. Although they look very similar, the constants appearing in the exponent are actually different: when comparing $\tilde{C}_\infty$ above to $C_\infty$ as in (3.12), we see that in $\tilde{C}_\infty$ the last addend is multiplied by $L_f$.

**Proposition 5.6. (BV estimate in space)** Let $\rho_o \in (L^\infty \cap \text{BV})(\mathbb{R}^2; \mathbb{R}^+)$. Let $(v)$, $(H)$, $(\eta)$, (5.5) and (5.6) hold. Then, for all $t > 0$, $\rho_\Delta$ in (3.1) constructed through Algorithm 5.1 satisfies the following estimate: for all $n = 0, \dots, N_T$,

$$ \sum_{i,j \in \mathbb{Z}} \left( \Delta y \left| \rho_{i+1,j}^n - \rho_{i,j}^n \right| + \Delta x \left| \rho_{i,j+1}^n - \rho_{i,j}^n \right| \right) \le \tilde{C}_x (t^n), $$

where

$$ \tilde{C}_x(t) = e^{2t\tilde{K}_1} \sum_{i,j \in \mathbb{Z}} \left( \Delta x \left| \rho_{i,j+1}^0 - \rho_{i,j}^0 \right| + \Delta y \left| \rho_{i+1,j}^0 - \rho_{i,j}^0 \right| \right) + \frac{2\tilde{K}_2}{\tilde{K}_1} (e^{2t\tilde{K}_1} - 1), \qquad (5.7) $$

with

$$ \tilde{K}_1 = 2 \|\nabla v^{\text{stat}}\|_{L^\infty} + 4\varepsilon L_f \|\nabla^2 \eta\|_{L^\infty} \|\rho_o\|_{L^1}, $$

$\tilde{K}_2$ as in (3.17) and $c_1, c_2$ are defined in (A.6).

**Remark 5.7.** Observe that $\tilde{K}_1 < \tilde{K}_2$ in (3.16).
**Corollary 5.8. (BV estimate in space and time)** Let $\rho_o \in (L^\infty \cap BV)(\mathbb{R}^2; \mathbb{R}^+)$, Let $(v)$, $(H)$, $(\eta)$, (5.5) and (5.6) hold. Then, for all $t > 0$, $\rho_\Delta$ in (3.1) constructed through Algorithm 5.1 satisfies the following estimate: for all $n = 1, \dots, N_T$,

$$ \sum_{m=0}^{n-1} \sum_{i,j \in \mathbb{Z}} \Delta t \left( \Delta y \left| \rho_{i+1,j}^m - \rho_{i,j}^m \right| + \Delta x \left| \rho_{i,j+1}^m - \rho_{i,j}^m \right| \right) + \sum_{m=0}^{n-1} \sum_{i,j \in \mathbb{Z}} \Delta x \Delta y \left| \rho_{i,j}^{m+1} - \rho_{i,j}^m \right| \le \tilde{C}_{xt}(t^n), $$

where

$$ \tilde{C}_{xt}(t) = t \left( \tilde{C}_x(t) + 2 \tilde{C}_t(t) \right), $$

with $\tilde{C}_x$ as in (5.7) and

$$ \tilde{C}_t(s) = 2 \left( \|v^{\text{stat}}\|_{L^\infty} + \varepsilon L_f \right) \tilde{C}_x(s) + \left( \|\nabla v^{\text{stat}}\|_{L^\infty} + \varepsilon \|\nabla^2 \eta\|_{L^\infty} \|\rho_o\|_{L^1} \right) \|\rho_o\|_{L^1}. $$

# 6 Numerical results

We consider the test setting given in [9] to compare the results of the Roe scheme, cf. Algorithm 3.1, to the results of the Lax-Friedrichs type scheme, cf. Algorithm 5.1.

## 6.1 Test setting

A total number of $N = 192$ parts in the shape of metal cylinders are transported on a conveyor belt moving with speed $v_T = 0.42$ m/s and are redirected by a diverter. The diverter is positioned at an angle of $\vartheta = 45$ degree with respect to the border of the conveyor belt. Figure 1 illustrates the static velocity field of the conveyor belt. Parts are transported with velocity $v^{\text{stat}} = (v_T, 0)$ in Region A and the diverter redirects the parts (Region C). The area behind the diverter (Region B) is modelled in a way that should prohibit parts from passing through the diverter, see [9]. The point $(x_d, y_d)$ marks the end of the diverter.

<figure><img src="image_4.png" /><figcaption>**Figure 1** – Schematic view of the static field of the conveyor belt.</figcaption></figure>

## 6.2 Discretisation and solution properties

To numerically model the setting, we introduce a uniform grid $\Delta x = \Delta y$ on the selected area of the conveyor belt. Initial conditions for the density at time $t = 0$ are given by the
experimental data and normalised so that $\rho_{max} = 1$. The mollifier $\eta$ which is used in the operator $I(\rho)$ (1.2) is chosen as follows

$$ \eta(x) = \frac{\sigma}{2\pi} e^{-\frac{1}{2}\sigma \|x\|_2^2}, $$

with $\sigma = 10000$.

In the original model formulation [9], the Heaviside function was introduced to avoid densities larger than $\rho_{max}$. In this work, we investigate the performances of two numerical schemes with two types of smooth approximations of the Heaviside function, one sensibly closer than the other. The former is the approximation $H_t$ (atan), obtained using the inverse tangent

$$ H_t(u) = \frac{\arctan(50(u-1))}{\pi} + 0.5, \qquad (6.1) $$

while the latter is denoted $H_p$ (polynomial) and it is obtained by cubic spline interpolation with the following conditions

$$ \begin{array}{llllll} H_p(u) = 0 & \forall u \le d_l < 1, & H_p(u) = 1 & \forall u \ge d_r > 1, & & \\ H_p(d_l) = 0, & H_p(1) = \frac{1}{2}, & H_p(d_r) = 1, & H'_p(d_l) = 0, & H'_p(d_r) = 0. & \end{array} $$

The approximation $H_p$ for $d_l = \frac{1}{2}$ and $d_r = \frac{8}{5}$ together with the inverse tangent approximation are depicted in Figure 2.

<figure><img src="image_5.png" /><figcaption>Figure 2 – The Heaviside function and approximations $H_t$ (atan) and $H_p$ (polynomial).</figcaption></figure>

Using the inverse tangent approximation corresponds to activating the collision operator $I(\rho)$ very close to $\rho_{max}$. On the other hand, with the polynomial approximation, the collision operator starts activating at $\frac{1}{2}\rho_{max}$, which implies that clusters with densities values between $\frac{1}{2}\rho_{max}$ and $\rho_{max}$ are already dispersed to some extent. Numerically, this means that densities above the maximum one are more likely to appear when using the inverse tangent approximation, while exploiting the polynomial approximation prevents from reaching such high values of the density.

Clearly, different approximations of the Heaviside function lead to different Lipschitz constants $L_f$ (2.1) and therefore influence the CFL time steps of the Roe scheme (3.13) and of the
Lax-Friedrichs scheme (5.5)–(5.6). Moreover, the constants $\alpha$ and $\beta$ given by the CFL condition for the Lax-Friedrichs scheme depend on the Lipschitz constant $L_f$. Larger Lipschitz constants, and therefore higher viscosity coefficients $\alpha$ and $\beta$, add additional viscosity to the Lax-Friedrichs scheme and therefore more diffusion, as shown in [14]. Note that in general, the Lax-Friedrichs scheme is more diffusive than the Roe scheme. To ensure conservation of mass within the given area of the conveyor belt, we impose zero-flux-conditions at the boundaries of the conveyor belt for the Lax-Friedrichs scheme. Therefore, at the boundary, the flux that would exit the domain is set to zero.

The Lipschitz constants of the approximations of the Heaviside functions depicted in Figure 2, as well as the corresponding CFL conditions, are displayed in Table 1, for a fixed space step size $\Delta x$. The inverse tangent approximation has a greater Lipschitz constant, leading to small CFL time steps and thus to an increased computational effort.

<table><caption>**Table 1** – CFL time steps for different Heaviside approximations and fixed $\Delta x, \Delta y$.</caption><thead><tr><th>Approximation</th><th>L<sub>f</sub></th><th>Δx = Δy [m]</th><th>CFL time step Roe [s]</th><th>CFL time step LxF [s]</th></tr></thead><tbody><tr><td>H<sub>t</sub> (atan)</td><td>16.42</td><td>1 · 10<sup>-2</sup></td><td>2.37 · 10<sup>-4</sup></td><td>1.21 · 10<sup>-4</sup></td></tr><tr><td>H<sub>p</sub> (polynomial)</td><td>2.09</td><td>1 · 10<sup>-2</sup></td><td>1.63 · 10<sup>-3</sup></td><td>9.50 · 10<sup>-4</sup></td></tr></tbody></table>

We analyze the amount of parts that pass the obstacle, i.e. the outflow at the end of the obstacle ($x_d, y_d$). The time-dependent mass function $U(t)$ counts the measured parts that are located in the region $\Omega_0$

$$U(t) = \frac{1}{N} \sum_{i=1}^{N} \chi_{\Omega_0}(x^{(i)}(t), y^{(i)}(t)) \qquad \chi_{\Omega_0}(x, y) = \begin{cases} 1, & (x, y) \in \Omega_0 \\ 0, & \text{otherwise} \end{cases} \qquad (6.2)$$

where $\Omega_0 = \{(x, y) \in \mathbb{R}^2 | x \le x_d\}$ is the left sided region upstream the obstacle and $(x^{(i)}(t), y^{(i)}(t))$ is the position of part $i$, $i \in 1, \dots, N$, at time $t$. The time-dependent mass function describing the outflow to the solution of the conservation law is given by

$$U_{\rho}(t) = \frac{1}{\int_{\Omega_0} \rho(x, 0) \, dx} \int_{\Omega_0} \rho(x, t) \, dx. \qquad (6.3)$$

The outflow curves obtained using Roe scheme, Lax-Friedrichs scheme and the outflow measured experimentally are shown in Figure 3. The parameters chosen for each scheme are those given in Table 1. Figure 4 displays the $L^\infty$ norms of the solution over time.
<figure><img src="image_6.png" /><figcaption>Figure 3 – Outflow.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4 – Time evolution of the $L^\infty$ norm of the approximate solutions.</figcaption></figure>

The outflow curve given by Roe scheme for both approximations of the Heaviside function is closer to the experimental data, due to the fact that the scheme captures more congestion, as indicated also by the $L^\infty$ norm. With Roe scheme, as expected, the density piles up even more when using the inverse tangent approximation. We observe that a maximum principle is not verified. On the contrary, with the polynomial approximation, higher densities are avoided, since the collision operator is activated earlier. Due to the influence of the viscosity coefficients, an opposite behaviour is observable with the Lax-Friedrichs scheme. Results for the $L^\infty$ norm of the solution are quite promising using the polynomial approximation, whereas the viscosity of the scheme is too large in the case of the inverse tangent approximation. The $L^\infty$ norm of the solution is constantly decreasing over time because of diffusion.

Figure 5 displays the parts' positions in the experiment and the density distribution computed with Roe and Lax-Friedrichs scheme using the polynomial Heaviside approximation at time $t = 1.5$ s. The density plot of the Roe scheme matches the experimental data quite well: regions with higher densities mostly coincide with regions in the experiments, where the parts are side by side. In contrast, the Lax-Friedrichs scheme produces a more widely spread density distribution. Even the parts on the upper section of the belt, which are transported to the right with the velocity of the conveyor belt, are not correctly portrayed.

<figure><img src="image_8.png" /><figcaption>Figure 5 – Experimental data (left), results of the Roe scheme (middle) and the Lax-Friedrichs scheme (right) at $t = 1.5$ s.</figcaption></figure>

Since the Roe scheme using the sharper approximation of the Heaviside function provides
the best result in comparison to the experimental data, we analyse its behaviour for $\Delta x, \Delta y \to 0$. Table 2 shows the error norms of the outflow of the simulations with the Roe scheme and the inverse tangent approximation of the Heaviside function compared to the outflow given by the experimental data.

<table><caption>**Table 2** – Error norms of Roe scheme, with the inverse tangent approximation of the Heaviside function (6.1), against experimental data.</caption><thead><tr><th>Δx = Δy</th><th>CFL time step</th><th>L<sup>1</sup>-error</th><th>L<sup>2</sup>-error</th><th>L<sup>∞</sup>-error</th></tr></thead><tbody><tr><td>4 · 10<sup>-2</sup></td><td>9.47 · 10<sup>-4</sup></td><td>0.42</td><td>0.26</td><td>0.20</td></tr><tr><td>2 · 10<sup>-2</sup></td><td>4.73 · 10<sup>-4</sup></td><td>0.16</td><td>0.10</td><td>0.10</td></tr><tr><td>1 · 10<sup>-2</sup></td><td>2.37 · 10<sup>-4</sup></td><td>0.09</td><td>0.07</td><td>0.09</td></tr><tr><td>5 · 10<sup>-3</sup></td><td>1.18 · 10<sup>-4</sup></td><td>0.07</td><td>0.05</td><td>0.07</td></tr></tbody></table>

The scheme is evaluated for different space step sizes and their corresponding CFL time steps. We observe that the error of the Roe scheme decreases as the space step decreases, suggesting the convergence of the outflow to the experimental data, compare also Figure 6.

<figure><img src="image_9.png" /><figcaption>**Figure 6** – Outflow computed by Roe scheme with different space step sizes using the inverse tangent approximation of the Heaviside function (6.1).</figcaption></figure>

## Acknowledgments

This work was financially supported by the project “Non-local conservation laws for engineering applications” co-funded by DAAD (Project-ID 57445223) and the PHC Procope (Project no. 42503RM) and the DFG project GO 1920/7-1. Elena Rossi was partially supported by University of Mannheim and acknowledges its hospitality during the accomplishment of this work.
Appendix A Technical Lemma

**Lemma A.1.** Let $\eta \in (C^3 \cap W^{3,\infty})(\mathbb{R}^2; \mathbb{R})$. Then, for $n = 0, \dots, N_T$, for $i,j \in \mathbb{Z}$, the following estimates hold:

$$ \|J_k^n\|_{L^\infty} \le \varepsilon \quad \text{for } k = 1, 2, \qquad (A.1) $$

$$ \left| J_1^n(x_{i+1/2,j}) - J_1^n(x_{i-1/2,j}) \right| \le 2 \varepsilon \Delta x \left\| \nabla^2 \eta \right\|_{L^\infty} \|\rho^n\|_{L^1}, \qquad (A.2) $$

$$ \left| J_1^n(x_{i+1/2,j}) - J_1^n(x_{i+1/2,j+1}) \right| \le 2 \varepsilon \Delta y \left\| \nabla^2 \eta \right\|_{L^\infty} \|\rho^n\|_{L^1}, \qquad (A.3) $$

$$ \left| J_2^n(x_{i,j+1/2}) - J_2^n(x_{i,j-1/2}) \right| \le 2 \varepsilon \Delta y \left\| \nabla^2 \eta \right\|_{L^\infty} \|\rho^n\|_{L^1}, $$

$$ \left| J_2^n(x_{i+1,j+1/2}) - J_2^n(x_{i,j+1/2}) \right| \le 2 \varepsilon \Delta x \left\| \nabla^2 \eta \right\|_{L^\infty} \|\rho^n\|_{L^1}, $$

$$ \left| J_1^n(x_{i+3/2,j}) - 2 J_1^n(x_{i+1/2,j}) - J_1^n(x_{i-1/2,j}) \right| \le 2 \varepsilon (\Delta x)^2 (c_1 \|\rho^n\|_{L^1} + c_2 \|\rho^n\|_{L^1}^2), \qquad (A.4) $$

$$ \left| J_2^n(x_{i,j+3/2}) - 2 J_2^n(x_{i,j+1/2}) - J_2^n(x_{i,j-1/2}) \right| \le 2 \varepsilon (\Delta y)^2 (c_1 \|\rho^n\|_{L^1} + c_2 \|\rho^n\|_{L^1}^2), $$

$$ \left| J_1^n(x_{i+1/2,j}) - J_1^n(x_{i+1/2,j+1}) - J_1^n(x_{i-1/2,j}) - J_1^n(x_{i-1/2,j+1}) \right| \le 2 \varepsilon \Delta x \Delta y C, \qquad (A.5) $$

$$ \left| J_2^n(x_{i,j+1/2}) - J_2^n(x_{i+1,j+1/2}) - J_2^n(x_{i,j-1/2}) - J_2^n(x_{i+1,j-1/2}) \right| \le 2 \varepsilon \Delta x \Delta y C, $$

where we set

$$ C = c_1 \|\rho^n\|_{L^1} + c_2 \|\rho^n\|_{L^1}^2, \quad c_1 = 2 \|\nabla^3 \eta\|_{L^\infty}, \quad c_2 = 3 \|\nabla^2 \eta\|_{L^\infty}^2. \qquad (A.6) $$

**Proof.** The proof of (A.1) is immediate.

Pass now to (A.2). For the sake of simplicity, introduce the following notation:

$$ D_+ = \sqrt{1 + \left\| (\nabla \eta * \rho^n)(x_{i+1/2}, y_j) \right\|^2}, \quad D_- = \sqrt{1 + \left\| (\nabla \eta * \rho^n)(x_{i-1/2}, y_j) \right\|^2}. $$

Hence,

$$ \left| J_1^n(x_{i+1/2}, y_j) - J_1^n(x_{i-1/2}, y_j) \right| \qquad (A.7) $$

$$ = \varepsilon \left| \frac{\Delta x \Delta y}{D_+} \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \partial_1 \eta(x_{i+1/2-k}, y_{j-\ell}) - \frac{\Delta x \Delta y}{D_-} \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \partial_1 \eta(x_{i-1/2-k}, y_{j-\ell}) \right| $$

$$ \le \varepsilon \left| \frac{\Delta x \Delta y}{D_+} \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \left( \partial_1 \eta(x_{i+1/2-k}, y_{j-\ell}) - \partial_1 \eta(x_{i-1/2-k}, y_{j-\ell}) \right) \right| \qquad (A.8) $$

$$ + \varepsilon \Delta x \Delta y \left| \frac{1}{D_+} - \frac{1}{D_-} \right| \sum_{k, \ell \in \mathbb{Z}} \left| \rho_{k, \ell}^n \right| \left| \partial_1 \eta(x_{i-1/2-k}, y_{j-\ell}) \right|. \qquad (A.9) $$

Consider (A.8): since $D_+ \ge 1$ and

$$ \left| \partial_1 \eta(x_{i+1/2-k}, y_{j-\ell}) - \partial_1 \eta(x_{i-1/2-k}, y_{j-\ell}) \right| \le \int_{x_{i-1/2-k}}^{x_{i+1/2-k}} \left| \partial_{11}^2 \eta(x, y_{j-\ell}) \right| dx, $$
we obtain

$$ [(A.8)] \leq \varepsilon \Delta x \left\| \partial_{11}^2 \eta \right\|_{\mathbf{L}^{\infty}} \| \rho^n \|_{\mathbf{L}^1}. \quad (A.10) $$

On the other hand, to estimate (A.9), compute

$$ \left| \frac{1}{D_+} - \frac{1}{D_-} \right| = \frac{|D_+ - D_-|}{D_+ D_-}. $$

Introduce $a(x) = \nabla \eta * \rho^n(x)$ and $b(z) = (1 + \|z\|^2)^{1/2}$, for $z \in \mathbb{R}^2$. In particular compute
$b'(z) = \frac{\|z\|}{(1 + \|z\|^2)^{1/2}}$ and observe that $|b'(z)| \le 1$. Then

$$ \begin{aligned} |D_+ - D_-| &= |b(a(x_{i+1/2})) - b(a(x_{i-1/2}))| = |b'(a(\tilde{x}_i)) a'(\tilde{x}_i) (x_{i+1/2} - x_{i-1/2})| \\ &= \left| \frac{a(\tilde{x}_i)}{(1 + a(\tilde{x}_i)^2)^{1/2}} (\partial_x \nabla \eta * \rho^n)(\tilde{x}_i) \Delta x \right| \\ &\le \Delta x \|\rho^n\|_{\mathbf{L}^1} \|\nabla^2 \eta\|_{\mathbf{L}^\infty}. \end{aligned} \quad (A.11) $$

Therefore,

$$ \varepsilon \Delta x \Delta y \left| \frac{1}{D_+} - \frac{1}{D_-} \right| \sum_{k, \ell \in \mathbb{Z}} \left| \rho_{k, \ell}^n \right| \left\| \partial_1 \eta(x_{i-1/2-k}, y_{j-\ell}) \right\| \le \varepsilon \Delta x \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty} \| \rho^n \|_{\mathbf{L}^1}. \quad (A.12) $$

Inserting (A.10) and (A.12) into the estimate of (A.7) yields the desired result.

Consider now (A.4). Introduce the following notation: for $\mu \in \{-1; 1; 3\}$ set

$$ D_{\mu} = \sqrt{1 + \left\| (\nabla \eta * \rho^n)(x_{i+\mu/2}, y_j) \right\|^2}. $$

Thus

$$ \begin{aligned} & J_1^n(x_{i+3/2,j}) - 2 J_1^n(x_{i+1/2,j}) + J_1^n(x_{i-1/2,j}) \\ &= - \varepsilon \left( \frac{(\partial_1 \eta * \rho^n)(x_{i+3/2,j})}{D_3} - 2 \frac{(\partial_1 \eta * \rho^n)(x_{i+1/2,j})}{D_1} + \frac{(\partial_1 \eta * \rho^n)(x_{i-1/2,j})}{D_{-1}} \right. \\ & \qquad \left. \pm \frac{(\partial_1 \eta * \rho^n)(x_{i+3/2,j})}{D_1} \pm \frac{(\partial_1 \eta * \rho^n)(x_{i-1/2,j})}{D_1} \right) \\ &= - \varepsilon \left( \left( \frac{1}{D_3} - \frac{1}{D_1} \right) (\partial_1 \eta * \rho^n)(x_{i+3/2,j}) + \frac{1}{D_1} \left( (\partial_1 \eta * \rho^n)(x_{i+3/2,j}) - (\partial_1 \eta * \rho^n)(x_{i+1/2,j}) \right) \right. \\ & \qquad \left. + \frac{1}{D_1} \left( (\partial_1 \eta * \rho^n)(x_{i-1/2,j}) - (\partial_1 \eta * \rho^n)(x_{i+1/2,j}) \right) + \left( \frac{1}{D_{-1}} - \frac{1}{D_1} \right) (\partial_1 \eta * \rho^n)(x_{i-1/2,j}) \right). \end{aligned} $$

Consider the terms separately, forgetting for a moment the $\varepsilon$ in front of everything. Focus first on the terms with common denominator $D_1$:

$$ \frac{1}{D_1} \left( (\partial_1 \eta * \rho^n)(x_{i+3/2,j}) - (\partial_1 \eta * \rho^n)(x_{i+1/2,j}) + (\partial_1 \eta * \rho^n)(x_{i-1/2,j}) - (\partial_1 \eta * \rho^n)(x_{i+1/2,j}) \right) $$
$$
\begin{align*}
&= \frac{\Delta x \Delta y}{D_1} \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \left( \partial_1 \eta(x_{i+3/2-k}, y_{j-\ell}) - \partial_1 \eta(x_{i+1/2-k}, y_{j-\ell}) \right. \\
&\qquad \left. + \partial_1 \eta(x_{i-1/2-k}, y_{j-\ell}) - \partial_1 \eta(x_{i+1/2-k}, y_{j-\ell}) \right) \\
&= \frac{\Delta x \Delta y}{D_1} \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \Delta x \left( \partial_{11}^2 \eta(\hat{x}_{i+1-k}, y_{j-\ell}) - \partial_{11}^2 \eta(\hat{x}_{i-k}, y_{j-\ell}) \right) \\
&= \frac{\Delta x \Delta y}{D_1} \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \Delta x \int_{\hat{x}_{i-k}}^{\hat{x}_{i+1-k}} \partial_{111}^3 \eta(x, y_{j-\ell}) \, dx \\
&\le 2 (\Delta x)^2 \left\| \partial_{111}^3 \eta \right\|_{\mathbf{L}^\infty} \| \rho^n \|_{\mathbf{L}^1}, \tag{A.13}
\end{align*}
$$

with $\hat{x}_{i-k} \in ]x_{i-1/2-k}, x_{i+1/2-k}[$. We are left with

$$
\left( \frac{1}{D_3} - \frac{1}{D_1} \right) (\partial_1 \eta * \rho^n)(x_{i+3/2,j}) + \left( \frac{1}{D_{-1}} - \frac{1}{D_1} \right) (\partial_1 \eta * \rho^n)(x_{i-1/2,j}). \quad (A.14)
$$

Add and subtract to (A.14)

$$
\left( \frac{1}{D_{-1}} - \frac{1}{D_1} \right) (\partial_1 \eta * \rho^n)(x_{i+3/2,j}).
$$

Hence,

$$
\left( \frac{1}{D_3} - 2 \frac{1}{D_1} + \frac{1}{D_{-1}} \right) (\partial_1 \eta * \rho^n)(x_{i+3/2,j}) \quad (A.15)
$$

$$
+ \left( \frac{1}{D_{-1}} - \frac{1}{D_1} \right) \left( (\partial_1 \eta * \rho^n)(x_{i-1/2,j}) - (\partial_1 \eta * \rho^n)(x_{i+3/2,j}) \right). \quad (A.16)
$$

Consider first (A.16): exploiting also (A.11), we obtain

$$
\begin{align*}
\text{[(A.16)]} &= \frac{D_1 - D_{-1}}{D_1 D_{-1}} \Delta x \Delta y \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \left( \partial_1 \eta(x_{i-1/2-k}, y_{j-\ell}) - \partial_1 \eta(x_{i+3/2-k}, y_{j-\ell}) \right) \\
&= \frac{D_1 - D_{-1}}{D_1 D_{-1}} \Delta x \Delta y \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \int_{x_{i+3/2-k}}^{x_{i-1/2-k}} \partial_{11}^2 \eta(x, y_{j-\ell}) \, dx \\
&\le 2 (\Delta x)^2 \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty}^2 \| \rho^n \|_{\mathbf{L}^1}^2. \tag{A.17}
\end{align*}
$$

As far as (A.15) is concerned, focus on the terms in the brackets:

$$
\begin{align*}
\frac{1}{D_3} - 2 \frac{1}{D_1} + \frac{1}{D_{-1}} &= \frac{D_1 D_{-1} - 2 D_3 D_{-1} + D_3 D_1}{D_3 D_1 D_{-1}} \\
&= \frac{D_{-1}(D_1 - D_3) - D_3(D_{-1} - D_1) \pm D_3(D_1 - D_3)}{D_3 D_1 D_{-1}} \\
&= \frac{(D_{-1} - D_3)(D_1 - D_3)}{D_3 D_1 D_{-1}} - \frac{D_{-1} - 2 D_1 + D_3}{D_1 D_{-1}}. \tag{A.18}
\end{align*}
$$

Inserting the first addend of (A.18) back into (A.15) yields

$$
\frac{(D_{-1} - D_3)(D_1 - D_3)}{D_3 D_1 D_{-1}} (\partial_1 \eta * \rho^n)(x_{i+3/2,j}) \le 2 (\Delta x)^2 \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty}^2 \| \rho^n \|_{\mathbf{L}^1}^2, \quad (A.19)
$$
where we exploit (A.11) twice and use the fact that $\frac{(\partial_1 \eta * \rho^n)(x_{i+3/2,j})}{D_3} \le 1$. Concerning the second addend of (A.18), focus on its numerator: with the notation introduced before (A.11),

$$
\begin{align*}
D_{-1} - 2D_1 + D_3 &= b(a(x_{i+3/2})) - 2b(a(x_{i+1/2})) + b(a(x_{-1/2})) \\
&= b'(a(\tilde{x}_{i+1})) a'(\tilde{x}_{i+1}) (x_{i+3/2} - x_{i+1/2}) - b'(a(\tilde{x}_i)) a'(\tilde{x}_i) (x_{i+1/2} - x_{i-1/2}) \\
&= \Delta x (b'(a(\tilde{x}_{i+1})) (\partial_1 \nabla \eta * \rho^n)(\tilde{x}_{i+1}) - b'(a(\tilde{x}_i)) (\partial_1 \nabla \eta * \rho^n)(\tilde{x}_i)) \\
&\quad \pm \Delta x b'(a(\tilde{x}_i)) (\partial_1 \nabla \eta * \rho^n)(\tilde{x}_{i+1}) \\
&= \Delta x [b'(a(\tilde{x}_{i+1})) - b'(a(\tilde{x}_i))] (\partial_1 \nabla \eta * \rho^n)(\tilde{x}_{i+1}) \\
&\quad + \Delta x b'(a(\tilde{x}_i)) [(\partial_1 \nabla \eta * \rho^n)(\tilde{x}_{i+1}) - (\partial_1 \nabla \eta * \rho^n)(\tilde{x}_i)] \\
&= \Delta x b''(a(\bar{x}_{i+1/2})) a'(\bar{x}_{i+1/2}) (\tilde{x}_{i+1} - \tilde{x}_i) (\partial_1 \nabla \eta * \rho^n)(\tilde{x}_{i+1}) \\
&\quad + \Delta x b'(a(\tilde{x}_i)) \Delta x \Delta y \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n (\partial_1 \nabla \eta(\tilde{x}_{i+1-k}, y_{j-\ell}) - \partial_1 \nabla \eta(\tilde{x}_{i-k}, y_{j-\ell})) \\
&= \Delta x b''(a(\bar{x}_{i+1/2})) a'(\bar{x}_{i+1/2}) (\tilde{x}_{i+1} - \tilde{x}_i) (\partial_1 \nabla \eta * \rho^n)(\tilde{x}_{i+1}) \\
&\quad + \Delta x b'(a(\tilde{x}_i)) \Delta x \Delta y \sum_{k, \ell \in \mathbb{Z}} \rho_{k, \ell}^n \int_{\tilde{x}_{i-k}}^{\tilde{x}_{i+1-k}} \partial_{11}^2 \nabla \eta(x, y_{j-\ell}) dx
\end{align*}
$$

where $\tilde{x}_i \in ]x_{i-1/2}, x_{i+1/2}[$ and $\bar{x}_{i+1/2} \in ]\tilde{x}_i, \tilde{x}_{i+1}[$. Now insert this estimate back into (A.18) and (A.15): since $|b''(z)| \le 1$,

$$
\left| \frac{D_{-1} - 2D_1 + D_3}{D_1 D_{-1}} (\partial_1 \eta * \rho^n)(x_{i+3/2,j}) \right| \le 2 (\Delta x)^2 \left[ \| \partial_1 \nabla \eta \|_{\mathbf{L}^\infty}^2 \| \rho^n \|_{\mathbf{L}^1}^2 + \| \partial_{11}^2 \nabla \eta \|_{\mathbf{L}^\infty} \| \rho^n \|_{\mathbf{L}^1} \right]. \quad (A.20)
$$

Collecting together (A.13), (A.17), (A.19) and (A.20) yields

$$
\left| J_1^n(x_{i+3/2}, y_j) - 2J_1^n(x_{i+1/2}, y_j) + J_1^n(x_{i-1/2}, y_j) \right| \\
\le 2\epsilon (\Delta x)^2 \left( 2 \left\| \nabla^3 \eta \right\|_{\mathbf{L}^\infty} \| \rho^n \|_{\mathbf{L}^1} + 3 \left\| \nabla^2 \eta \right\|_{\mathbf{L}^\infty}^2 \| \rho^n \|_{\mathbf{L}^1}^2 \right).
$$

□

References

[1] A. Aggarwal, R. M. Colombo, and P. Goatin. Nonlocal systems of conservation laws in several space dimensions. *SIAM J. Numer. Anal.*, 53(2):963–983, 2015.

[2] P. Amorim, R. M. Colombo, and A. Teixeira. On the numerical integration of scalar nonlocal conservation laws. *ESAIM Math. Model. Numer. Anal.*, 49(1):19–37, 2015.

[3] S. Blandin and P. Goatin. Well-posedness of a conservation law with non-local flux arising in traffic flow modeling. *Numer. Math.*, 132(2):217–241, 2016.

[4] F. A. Chiarello, P. Goatin, and E. Rossi. Stability estimates for non-local scalar conservation laws. *Nonlinear Anal. Real World Appl.*, 45:668–687, 2019.
[5] R. M. Colombo, M. Garavello, and M. Lécureux-Mercier. A class of nonlocal models for pedestrian traffic. *Math. Models Methods Appl. Sci.*, 22(4):1150023, 34, 2012.

[6] R. M. Colombo and E. Rossi. Hyperbolic predators vs. parabolic prey. *Commun. Math. Sci.*, 13(2):369–400, 2015.

[7] M. Crandall and A. Majda. The method of fractional steps for conservation laws. *Numer. Math.*, 34(3):285–314, 1980.

[8] M. G. Crandall and A. Majda. Monotone difference approximations for scalar conservation laws. *Math. Comp.*, 34(149):1–21, 1980.

[9] S. Göttlich, S. Hoher, P. Schindler, V. Schleper, and A. Verl. Modeling, simulation and validation of material flow on conveyor belts. *Appl. Math. Model.*, 38(13):3295–3313, 2014.

[10] K. H. Karlsen and N. H. Risebro. On the uniqueness and stability of entropy solutions of nonlinear degenerate parabolic equations with rough coefficients. *Discrete Contin. Dyn. Syst.*, 9(5):1081–1104, 2003.

[11] S. N. Kružhkov. First order quasilinear equations with several independent variables. *Mat. Sb.* (N.S.), 81 (123):228–255, 1970.

[12] M. Lécureux-Mercier. Improved stability estimates for general scalar conservation laws. *J. Hyperbolic Differ. Equ.*, 8(4):727–757, 2011.

[13] M. Lécureux-Mercier. Improved stability estimates on general scalar balance laws. *ArXiv e-prints*, July 2013.

[14] R. J. LeVeque. *Numerical methods for conservation laws*. Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, second edition, 1992.

[15] R. Sanders. On convergence of monotone finite difference schemes with variable spatial differencing. *Math. Comp.*, 40(161):91–106, 1983.