# A particle method for stationary transport equations

Rafael Bailo, Julie Binard, Pierre Degond, Pascal Noble

► To cite this version:

Rafael Bailo, Julie Binard, Pierre Degond, Pascal Noble. A particle method for stationary transport equations. 2025. hal-05359514

HAL Id: hal-05359514

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

Preprint submitted on 11 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A particle method for stationary transport equations

Rafael Bailo¹, Julie Binard²,³, Pierre Degond³, and Pascal Noble²

¹Department of Mathematics and Computer Science, Eindhoven University of Technology; Eindhoven, Netherlands, email: r.bailo@tue.nl,

²Institut de Mathématiques de Toulouse; UMR 5219, Université de Toulouse; CNRS, INSA, F-31077 Toulouse, France, email: JB: binard@insa-toulouse.fr, email: PN: noble@insa-toulouse.fr

³Departamento de Matemática Aplicada I, Universidad de Sevilla, E.T.S. Arquitectura. Avda Reina Mercedes, 41012 Sevilla, Spain, email: jbinard@us.es,

⁴Institut de Mathématiques de Toulouse; UMR 5219, Université de Toulouse; CNRS, UPS, F-31062 Toulouse Cedex 9, France, email: pierre.degond@math.univ-toulouse.fr

## Abstract

We present and study a particle method for the stationary solutions of a class of transport equations. This method is inspired by non-stationary particle methods, the time variable being replaced by one spatial variable. Particles trajectories are computed using the “time-dependent” equations, and then the approximation is based on a quadrature method using the particle locations as quadrature points. We prove the convergence of the scheme under suitable regularity assumptions on the data and the solution, together with a “characteristic completeness” assumption (the characteristic curves fullfill the whole computational domain). We also provide an error estimate. The scheme is tested numerically on a two dimensional linear equation and we present a numerical study of convergence. Finally, we use this method to carry out numerical simulations of a landscape evolution model, where an erodible topography evolves under the effects of water erosion and sedimentation. The scheme is then useful to deal with wet/dry areas.

**Keywords** Transport equation, Stationary solutions, Particle method, Landscape evolution model.

**Mathematics Subject Classification** 65N12, 65N75, 35B36, 35M30, 35Q86, 86-10.

## 1 Introduction

In this paper we introduce and analyze a particle method for first order equations in divergence form and with a source term:

$$ \mathrm{div}(au)(x) + a_0(x)u(x) = S(x), \quad \forall x \in \mathbb{R}_+^d := (0, +\infty) \times \mathbb{R}^{d-1}, \quad (1) $$

with $d \ge 2$. We assume that the functions $a, a_0, S$ are defined on an open set $U \subset \mathbb{R}^d$ containing $\overline{\mathbb{R}}_+^d := [0, +\infty) \times \mathbb{R}^{d-1}$, and $a: U \to \mathbb{R}^d \setminus \{0\}$ whereas $a_0, S: U \to \mathbb{R}$. We complete the system (1) with the following
boundary condition:

$$u(0, \xi) = g(\xi), \quad \forall \xi \in \mathbb{R}^{d-1}, \qquad (2)$$

where $x \in \mathbb{R}^d$ is written as $x = (x_1, \xi)$ with $x_1 \in \mathbb{R}$ and $\xi \in \mathbb{R}^{d-1}$. We set $\xi = (x_2, \dots, x_d) = (\xi_1, \dots, \xi_{d-1})$.

The problem originates from the numerical simulation of landscape evolution models: a landscape is modeled by a function $z: x \mapsto z(x)$ and it is eroded by a fluid layer of height $h$. The sediments created by erosion are transported by the fluid, and the sediment density in the fluid is denoted by $c$. The quantities $(z, h, c)$ are functions of the horizontal coordinate $x$ and time $t$. The evolution system for $(z, h, c)$ introduced in [3, 4] and studied in [2] is written:

$$\left\{
\begin{aligned}
\partial_t h + \operatorname{div}(h\mathbf{v}) &= r, \\
\partial_t (hc) + \operatorname{div}(ch\mathbf{v}) &= \rho_s e \left(\frac{h}{H}\right)^m \left(\frac{|\mathbf{v}|}{V}\right)^n - \rho_s s \frac{c}{c_{sat}}, \\
\partial_t z &= K\Delta z - e \left(\frac{h}{H}\right)^m \left(\frac{|\mathbf{v}|}{V}\right)^n + s \frac{c}{c_{sat}}.
\end{aligned}
\right.
\qquad (3)$$

Here, $\rho_s$ represents the sediment volumetric mass, $e$ is an erosion rate, $s$ a sediment redeposition rate whereas $r$ is a source term representing precipitation, evaporation or infiltration. The constant $K$ is a diffusion constant that represents the soil creep effect. The creep effect is a large-scale effect of soil diffusion that was introduced in [6] and used in [12] to explain the convexity of hilltop profiles. The erosion is modeled by a stream incision law, see [15]. This system is written in a closed form with the fluid velocity written as $\mathbf{v} = -V\nabla(h+z)$, with $V \ge 0$ a characteristic speed of the fluid, and where $h+z$ is the surface elevation.

Usually, the fluid speed $|\mathbf{v}|$ is much larger than the erosion rate $|\partial_t z|$, see e.g. the experiment in [14]. In order to focus on the erosion phenomenon, one has to consider large time scales. In this asymptotic regime, the time variations of the fluid height and sediment concentration become negligible: one is left with the equations:

$$\left\{
\begin{aligned}
\operatorname{div}(h\mathbf{v}) &= r, \\
\operatorname{div}(ch\mathbf{v}) &= \rho_s e \left(\frac{h}{H}\right)^m \left(\frac{|\mathbf{v}|}{V}\right)^n - \rho_s s \frac{c}{c_{sat}},
\end{aligned}
\right.
\qquad (4)$$

supplemented with the equation for $z$ which remains unchanged. One problem of interest is the erosion of an inclined plane, which can be carried out experimentally [14]. From a theoretical point of view, this corresponds to a bottom topography $z(x) = -\tan(\theta)x_1 + \tilde{z}(x)$ with $|\tilde{z}| \ll 1$. The fluid height is prescribed on one side of the plane: $h(t, 0, \xi) = h_{\text{ref}}(t, \xi)$ with $\xi \in \mathbb{R}$.

The topography $z$ being fixed, the first equation of (4) is non-linear and could be solved through a fixed point procedure: with $h^0$ an initial guess, one solves

$$\operatorname{div}(\mathbf{a}^n h^{n+1}) = r, \quad \mathbf{a}^n = -V\nabla(h^n + z).$$

We are thus left with a problem similar to the linear problem (1)-(2). By splitting space variables, equation (1) can be written as

$$\frac{\partial}{\partial x_1}(\mathbf{a}_1 u)(x_1, \xi) + \sum_{j=1}^{d-1} \frac{\partial}{\partial \xi_j}(\mathbf{a}_{j+1} u)(x_1, \xi) + a_0(x_1, \xi)u(x_1, \xi) = S(x_1, \xi), \quad \forall x_1 > 0, \forall \xi \in \mathbb{R}^{d-1}. \quad (5)$$
The boundary condition (2) makes sense only if it is not characteristic: in what follows, we will assume that $a_1(0, \xi) > 0, \forall \xi \in \mathbb{R}^{d-1}$. If we further assume that $a_1 > 0$, equation (5) can be seen as a first-order partial differential equation where $x_1$ plays the role of time. Equation (5) can be solved by a classical finite-volume scheme; however, a CFL-type condition should be imposed to preserve the positivity of $u$, which could be restrictive for small slopes. Moreover, it may be difficult to design the mesh of the computational domain in the presence of complex or moving boundaries. In this paper, we consider a particle method.

Particle methods are mesh-free numerical schemes that are mostly used to solve fluid-dynamical problems. The principle of such methods is to use a Lagrangian point of view. The solution is seen as a sum of Dirac distributions carrying some information (density, velocity). These Dirac distributions can be seen as particles advected by the fluid flow and the characteristics of each particle are evolved in time by using the physical laws. These Dirac masses are then approximated by kernel functions $\zeta^h, h > 0$, which then provide a classical and computable approximation of the solution. Smoothed Particle Hydrodynamics (SPH), Particle In Cell (PIC) methods, and vortex methods are examples of such particle methods. The SPH method was initially introduced to solve numerically problems in astrophysics. It was first designed by Gingold and Monaghan in [13], and by Lucy in [18]. The numerical resolution of astrophysical problems in three dimensions, by finite difference schemes requires a very large number of grid points, whereas the SPH method reaches the same accuracy with a comparatively much smaller number of particles. The SPH method is now used in a wider range of applications like oceanography or structural mechanics and high velocity impacts with large deformations or changes of topology. In plasma physics, the Vlasov-Poisson or Vlasov-Maxwell equations are set in a 6 dimensional space and it would be very costly to solve them with a classical scheme. The PIC method [10] is also a particle method which also decreases the computational cost. The vortex method is a type of particle method dedicated to the simulation of incompressible Euler equations [22]: it is based on an alternative formulation of the equations involving the vorticity of the fluid velocity, the fluid velocity being reconstructed with theory of potential. This method was designed and implemented in the 70s in two dimensions [5], and later in three dimensions [16]. The theoretical analysis of the particle methods was initiated by Raviart and coworkers. In [19], Mas-Gallic and Raviart analyzed a particle method for one class of first order linear hyperbolic systems. They proved the convergence of the continuous in time scheme in $L^2$. There is now an important literature on the convergence of the SPH method for various types of equations, both in the linear and in the non-linear case: see e.g. [1] for non-linear conservation laws, more recent references [11] for a systematical derivation and convergence analysis, and [17] for a review of recent results and challenges.

Equation (1) is not, in the strict sense of the term, an evolution equation, and classical particle methods do not apply to it directly. However, the $x_1$ variable could play the role of a time variable, whereas the boundary condition (2) could be seen as an initial condition provided that it not characteristic, which means that $a_1 \neq 0$ on the boundary $\{0\} \times \mathbb{R}^{d-1}$. In this paper we introduce a modified particle method to solve the stationary equation (1), together with the boundary condition (2). This method is based on non-stationary particle methods, and we prove the convergence of the numerical scheme under appropriate assumptions. The idea of the method and a part of the proof of convergence are inspired from [21], which analyses a particle method for linear hyperbolic equations of the first order. The result can be loosely formulated as follows:

**Theorem 1.1.** Let $d \in \mathbb{N}$ such that $d \ge 2$. Assume the data **a**, $a_0, S, g$ are smooth enough and the characteristic curves fullfill the domain $\mathbb{R}_+^d$. Then, for all compact sets $X \subset \mathbb{R}_+^d$, there exist constants $C_X^{(0)}, \dots, C_X^{(3)} > 0$
which do not depend on h, ε, Δs, such that, if $C_X^{(0)} \Delta s \le h$, the approximation $\tilde{u}_h^{\varepsilon, \Delta s}$ given by the scheme (31) verifies the error estimate:

$$ \sup_{x \in X} |u(x) - \tilde{u}_h^{\varepsilon, \Delta s}(x)| \le C_X^{(1)} h^2 + C_X^{(2)} \frac{\varepsilon^d}{h^d} + C_X^{(3)} \frac{\Delta s}{h}, \quad (6) $$

where h is the size of support of the kernel function $\zeta^h$, and where $\varepsilon$ and $\Delta s$ are respectively, the spatial mesh sizes for the transverse variable $\xi \in \mathbb{R}^{d-1}$ and $s \in \mathbb{R}_+$.

The first term in the right hand side of (6) is an interpolation error, stemming from the approximation of u by $u * \zeta^h$ with $\zeta^h$ a kernel function approximating the Dirac mass. The function $u * \zeta^h$ is written in a integral form involving the characteristic curves and the values of the solution along these curves. The second term is a quadrature error, caused by the discretization of the former integral along characteristics. The last term is the discretization error “in time”, associated to the time discretization by an Euler scheme of the ordinary differential equations for the characteristic curves and the values of the solution u along these curves.

This paper is organized as follow. First, in section 2, we introduce a representation formula for the solution u of equation (1)-(2) along the characteristic curves. Then, in section 3, we introduce our numerical method and prove a convergence result on the scheme under suitable assumptions. This result is completed in section 4 by numerical convergence tests for various test cases. Finally, we apply our numerical scheme to carry out numerical simulations of a landscape evolution model (3). Here the landscape is a tilted plane eroded by a fluid.

## 2 A representation formula for the solution using characteristics

In this section, we provide an integral representation of the solution u of (1)-(2) along the characteristic curves. We shall use this representation of the solution u to build our numerical scheme. We suppose that the boundary $\{0\} \times \mathbb{R}^{d-1}$ is non characteristic and that a is bounded on $\{0\} \times \mathbb{R}^{d-1}$: as a consequence, there exist two positive constants $\alpha$ and $\beta$ such that

$$ \forall \xi \in \mathbb{R}^{d-1}, \quad 0 < \alpha \le \mathbf{a}_1(0, \xi), \quad |\mathbf{a}(0, \xi)| \le \beta. \quad (7) $$

The first assumption is consistent with the original problem [2] where water is injected in the domain from the boundary $\{0\} \times \mathbb{R}^{d-1}$ of $\mathbb{R}_+^d$. We further assume that the function $\mathbf{a} \in C^1(U)$. The Cauchy-Lipschitz theorem applies, and we can define the characteristic curves $x_\xi$ as the maximal solutions of

$$ x'_{\xi}(s) = \mathbf{a}(x_{\xi}(s)), \quad \text{with} \quad x_{\xi}(0) = (0, \xi) \in \{0\} \times \mathbb{R}^{d-1}. \quad (8) $$

The characteristics $x_\xi$ are defined on the interval $(T_{\min}(\xi), T_{\max}(\xi))$ with $T_{\min}(\xi) < 0 < T_{\max}(\xi)$ for all $\xi \in \mathbb{R}^{d-1}$. We associate to the characteristic curves $x_\xi$, $\xi \in \mathbb{R}^{d-1}$ the flow function $\Phi$ defined on

$$ \Omega = \{(s, \xi); \xi \in \mathbb{R}^{d-1}, s \in (T_{\min}(\xi), T_{\max}(\xi))\} $$

as

$$ \begin{array}{l} \Phi: \Omega \to \mathbb{R}^d, \\ (s, \xi) \mapsto \Phi(s, \xi) = x_{\xi}(s). \end{array} \quad (9) $$

We then have the following representation formula for the solution u of System (1)-(2) along the characteristics.
<figure><img src="image_2.png" /><figcaption>Figure 1: Domains $\Omega_\ell \subset \Omega$ and $\Phi(\Omega_\ell) = \tilde{U}_\ell \subset U$.</figcaption></figure>

**Lemma 2.1.** Assume $\mathbf{a} \in C^2(U)$, $a_0, S \in C^1(U)$; then, the solution $u \in C^1(\Phi(\Omega))$ of (1)-(2) satisfies

$$ \forall (s, \xi) \in \Omega, \quad u(\Phi(s, \xi)) = \left[ g(\xi) + \int_0^s S(\Phi(t, \xi)) \exp(I(t, \xi)) dt \right] \exp(-I(s, \xi)), \quad (10) $$

$$ \text{where} \quad I(s, \xi) = \int_{0}^{s} \left( \operatorname{div}(\mathbf{a}(\Phi(\tau, \xi))) + a_{0}(\Phi(\tau, \xi)) \right) d\tau. $$

We want to solve the system (1)-(2) on the domain $\mathbb{R}_+^d$. The assumption (7) on $\mathbf{a}_1$ being satisfied, this domain is positively invariant under the flow of the differential system (8):

$$ \Phi(s, \xi) \in \mathbb{R}_+^d, \quad \forall s \ge 0, \quad \forall \xi \in \mathbb{R}^{d-1}. $$

However, this does not ensure that we can define the solution $u$ on the whole domain $\mathbb{R}_+^d$ and we have to make further assumptions on the flow $\Phi$:

**Assumption 2.1.**

(i) There exists $\ell > 0$ such that the flow $\Phi$ is defined on a domain $\Omega_\ell = (-\ell, +\infty) \times \mathbb{R}^{d-1} \subset \Omega$ and $\Phi: \Omega_\ell \to \Phi(\Omega_\ell) := \tilde{U}_\ell \subset U$ is a $C^1$-diffeomorphism (see figure 1):

$$ \forall x \in \tilde{U}_\ell, \quad \exists (s, \xi) \in (-\ell, +\infty) \times \mathbb{R}^{d-1}, \quad x = x_\xi(s). $$

(ii) There exists $\delta > 0$ such that $U_\delta = (-\delta, +\infty) \times \mathbb{R}^{d-1} \subset \tilde{U}_\ell \subset U$. This is illustrated in Figure 2.

These assumptions ensure that we can define the solution $u$ on $\mathbb{R}_+^d$. However, note that we have required a stronger assumption so as to define the solution $u$ on the domain $U_\delta$. Indeed, in the interpolation procedure, we approximate $u$ by $u * \zeta^h$ and we need that the solution $u$ to be defined on a domain $U_h$ to compute $u * \zeta^h$ on $\mathbb{R}_+^d$. Here, we can carry out these computations by assuming $0 < h < \delta$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Characteristic curves on the domain $\tilde{U}_{\ell}$ and definition of $U_{\delta}$.</figcaption></figure>

Note also that our setting does not deal with the problem of “dry areas” when the characteristics curves starting from the boundary $\{0\} \times \mathbb{R}^{d-1}$ do not cover entirely the domain $\mathbb{R}_+^{d-1}$. We will discuss in the numerical part some strategies to extend our framework.

For the reader's convenience, we list the various domains considered in the paper:

• $U$ is the domain of definition of **a**.

• $\Omega$ is the domain of definition of $\Phi$, and $\Phi(\Omega) \subset U$.

• $\Omega_{\ell} = (-\ell, +\infty) \times \mathbb{R}^{d-1} \subset \Omega$.

• $\tilde{U}_{\ell} = \Phi(\Omega_{\ell}) \subset U$.

• $U_{\delta} = (-\delta, +\infty) \times \mathbb{R}^{d-1} \subset \tilde{U}_{\ell}$.

Under suitable assumptions on the vector field **a**, the flow $\Phi$ satisfies the assumptions formulated above. The following proposition is proved in Appendix A.

**Proposition 2.1.** Let $L > 0$, and assume that $U \supset U_L = (-L, +\infty) \times \mathbb{R}^{d-1}$. Suppose that the function **a** is globally Lipschitz continuous on $U$ and satisfies (7). Then,

(i) For all $\xi \in \mathbb{R}^{d-1}$, equation (8) admits a unique maximal solution $x_{\xi}$. This solution is such that: $T_{\min}(\xi) < 0$, $T_{\max}(\xi) = +\infty$. The solution $x_{\xi}$ is $C^1((T_{\min}(\xi), +\infty))$ and the flow $\Phi$ is $C^1(\Omega)$.

(ii) For all $\xi \neq \eta \in \mathbb{R}^{d-1}$, the curves $s \mapsto x_{\xi}(s)$ and $s \mapsto x_{\eta}(s)$ do not intersect each other.

(iii) There exists $\ell > 0$ such that $\Phi$ verifies the first statement of assumption 2.1.

(iv) Assume that $\mathbf{a}_1(x) \ge \alpha > 0, \forall x \in U$. Then $\Phi$ verifies the second statement of assumption 2.1.
We also consider the smoothness of the solution $u$. In order to derive error estimates for our numerical scheme, we need additional regularity on the solution $u$ which is obtained by assuming more regularity on $a$, $a_0$, $S$ and $g$.

**Proposition 2.2.** Let $m \in \mathbb{N}$, and suppose that $a \in C^{m+1}(U)$, $g \in C^m(\mathbb{R}^{d-1})$ whereas $a_0$, $S \in C^m(U)$. Then, one has:

(i) $\Phi : \Omega_\ell \to \tilde{U}_\ell \subset U$ is a $C^{m+1}$ diffeomorphism.

(ii) The system (1)-(2) is well-posed on $\tilde{U}_\ell$ and the solution $u$ belongs to $C^m(\tilde{U}_\ell)$.

See [8, chap V] for the proof of the first item of Proposition 2.2.

# 3 The discretization method

In this section we describe the particle method and provide a proof of convergence of the numerical scheme. We follow the strategy of the SPH method: the solution is written as a sum of Dirac distributions. The Dirac masses are advected by the flow described by the characteristic equations whereas the values of the solution at the Dirac masses locations are evolved by using the partial differential equation. In order to deal with functions, the Dirac masses are replaced by a compactly supported family of functions $\xi^h$, $h > 0$ (so called kernel functions) approximating the Dirac distribution. In our setting, this program can be described as follow: we first replace the solution $u$ by its regularized version $u * \xi^h$. This later function can be written as an integral over $\tilde{U}_\ell$ for $h$ sufficiently small. We then use the flow $\Phi$ to perform a change of variable where the functions are evaluated along characteristics and the integral is taken on $\Omega_\ell$. We use a quadrature formula to discretize this integral with a spatial step size $\varepsilon > 0$ in the transverse direction $\xi \in \mathbb{R}^{d-1}$ and $\Delta s > 0$ in the direction $s > -\ell$. We obtain a fully discretized scheme by using the characteristic equations to compute an approximation of the particles position and an approximation of the solution at these discretization points.

In section 3.1, we introduce the kernel function $\zeta^h$, defined as an approximation of the Dirac delta distribution, and the change of variables formula in integrals defined on $\Phi(\Omega_\ell)$ associated to the characteristic flow $\Phi$. This change of variable is used to write the integral defining $u * \zeta^h$ along the characteristic curves. In section 3.2, we present the quadrature formula, that consists of a semi-discretization of the former integral in the direction $\xi \in \mathbb{R}^{d-1}$. In section 3.3, we introduce the discretization of the characteristic curves which in turn introduces our fully discrete scheme. Finally, the convergence theorem is proved in section 3.4.

## 3.1 The smoothed particle formalism

In what follows, we assume that $x \in \mathbb{R}_+^d := (0, +\infty) \times \mathbb{R}^{d-1}$ and we want to approximate $u$ by a smoother function $u * \zeta^h$ on this domain. For that purpose, we define a kernel function $\zeta^h$, an approximation of the Dirac delta distribution:

$$ \forall x \in \mathbb{R}^d, \quad \zeta^h(x) = \frac{1}{h^d} \zeta\left(\frac{|x|}{h}\right) \qquad (11) $$

with $\zeta: [0, +\infty) \to \mathbb{R}^+$ a function compactly supported in $[0, 1]$ which verifies

$$ \int_{\mathbb{R}^d} \zeta(|x|) dx = \mathrm{meas}(\mathbb{S}^{d-1}) \int_0^1 \zeta(t) t^{d-1} dt = 1, \qquad (12) $$
where meas($\mathbb{S}^{d-1}$) is the $d-1$ dimensional measure of the $d-1$ dimensional sphere. The support of $\zeta^h$ is included into $\bar{B}_d(0, h)$ where, for all $d \in \mathbb{N}$, $x \in \mathbb{R}^d$ and $h > 0$, the sets $B_d(x, h), \bar{B}_d(x, h)$ are given by

$$B_d(x, h) := \{y \in \mathbb{R}^d; |x - y| < h\}, \quad (13)$$

$$\bar{B}_d(x, h) := \{y \in \mathbb{R}^d; |x - y| \le h\}. \quad (14)$$

Then, for any $f \in C^m(U)$ and $x \in \mathbb{R}_+^d$, the function $y \mapsto f(y)\zeta^h(x-y)$ is compactly supported in $U_\delta$ provided that $0 < h < \delta$. This function can be extended by 0 to the domain $\tilde{U}_\ell$ and we get the approximation formula:

$$f(x) \approx f * \zeta^h(x) = \int_{U_\delta} f(y) \zeta^h(x-y) dy = \int_{\tilde{U}_\ell} f(y) \zeta^h(x-y) dy. \quad (15)$$

Next, we use the fact that the flow $\Phi$ is a $C^1$ diffeomorphism so that we can perform a change of variables in the formula (15). Let us denote by $J_\Phi$ the Jacobian matrix of $\Phi$, given by:

$$J_{\Phi}(s, \xi) = \begin{bmatrix} \frac{\partial \Phi_1}{\partial s} & \frac{\partial \Phi_1}{\partial \xi_1} & \cdots & \frac{\partial \Phi_1}{\partial \xi_{d-1}} \\ \vdots & \vdots & & \vdots \\ \frac{\partial \Phi_d}{\partial s} & \frac{\partial \Phi_d}{\partial \xi_1} & \cdots & \frac{\partial \Phi_d}{\partial \xi_{d-1}} \end{bmatrix} \quad (16)$$

As $\Phi$ is invertible, $\det J_\Phi$ does not vanish. Moreover, from the relation $\Phi(0, \xi) = (0, \xi)$, one deduces that

$$\nabla_{\xi} \Phi_1(0, \xi) = 0, \quad \left. \frac{\partial(\Phi_2, \dots, \Phi_d)}{\partial(\xi_1, \dots, \xi_{d-1})} \right|_{(0, \xi)} = \mathrm{Id}_{\mathbb{R}^{d-1}}, \quad \forall \xi \in \mathbb{R}^{d-1}.$$

As a result, for all $\xi \in \mathbb{R}^{d-1}$, one has

$$(\det J_{\Phi})(0, \xi) = \partial_{\xi} \Phi_1(0, \xi) = \mathbf{a}_1(0, \xi) \ge \alpha > 0.$$

Therefore, by using the Liouville formula, one deduces that

$$\forall s \ge -\ell, \forall \xi \in \mathbb{R}^{d-1}, (\det J_{\Phi})(s, \xi) = \exp \left( \int_{0}^{s} \mathrm{div}(\mathbf{a})(\Phi(t, \xi)) dt \right) (\det J_{\Phi})(0, \xi) > 0.$$

Then, for any $f \in C^0(\tilde{U}_\ell)$, by performing the change of variables $\Phi: (-\ell, +\infty) \times \mathbb{R}^{d-1} \to \tilde{U}_\ell$, we get:

$$\int_{\tilde{U}_\ell} f(x) dx = \int_{-\ell}^{+\infty} \int_{\mathbb{R}^{d-1}} f(\Phi(s, \xi)) (\det J_\Phi)(s, \xi) d\xi ds. \quad (17)$$

Consequently, an integral on $\tilde{U}_\ell$ can be converted into an integral on $(-\ell, +\infty) \times \mathbb{R}^{d-1}$, provided that one can determine the function along the characteristic curves. It is easier to perform computations on this domain, and to create a mesh on it.

We can now apply the approximation (15) to the solution $u$ of equation (1) together with the change of variables (17) since we have a representation formula for the solution $u$ along characteristic curves. One finds:

$$\forall x \in \mathbb{R}_+^d, \quad u(x) \approx \int_{\tilde{U}_\ell} u(y) \zeta^h(x-y) dy = \int_{-\ell}^{+\infty} \int_{\mathbb{R}^{d-1}} u(\Phi(s, \xi)) \zeta^h(x - \Phi(s, \xi)) (\det J_\Phi)(s, \xi) d\xi ds. \quad (18)$$

We conclude this section with an error estimate on $u - u * \zeta^h$.
**Lemma 3.1.** Let $u \in C^2(\tilde{U}_\ell)$, $\delta > 0$ such that $\overline{U}_\delta \subset \tilde{U}_\ell$, and $x \in \mathbb{R}_+^d$. Then there exists $C(x, \delta) > 0$ such that, for any $0 < h < \delta$, one has

$$|u(x) - u * \zeta^h(x)| \leq C(x, \delta)h^2. \qquad (19)$$

**Remark 3.1.** This approximation formula is order two, but this could be increased if we remove the hypothesis of positivity on $\zeta$. Actually, if the m first moments of $\zeta$ vanish, then the approximation has order m, see [21]. However, for the applications we have in mind, like the simulation of a landscape evolution model, the unknown u is either a fluid height or a sediment concentration, so we restrict our discussion to $\zeta$ a non-negative function.

**Proof.** Let $x \in \mathbb{R}_+^d$: we expand the function u up order 2 with respect to $z = y - x \to 0$. As u is $C^2(\tilde{U}_\ell)$, for all $y \in \overline{U}_\delta$, one has:

$$u(x) - u(y) = (x - y) \cdot \nabla u(x) + R(x, y),$$

where R is a remainder which verifies:

$$|R(x, y)| \leq |x - y|^2 \sup_{z \in [x, y]} |H(u)(z)|,$$

Where $|H(u)(z)|$ denotes an appropriate matrix norm of the Hessian matrix of u at the point z. Thus, as the integral of $\zeta^h$ is equal to 1, one finds

$$\begin{aligned} u(x) - \int_{\mathbb{R}^d} u(y) \zeta^h(x-y) dy &= \int_{\mathbb{R}^d} (u(x) - u(y)) \zeta^h(x-y) dy \\ &= \int_{\mathbb{R}^d} (x-y) \zeta^h(x-y) dy \cdot \nabla u(x) + \int_{\mathbb{R}^d} R(x,y) \zeta^h(x-y) dy. \end{aligned}$$

As $\zeta^h$ is a radial function, we have

$$\int_{\mathbb{R}^d} (x - y) \zeta^h(x - y) dy = 0.$$

Moreover, the support of the function $y \mapsto \zeta^h(x-y)$ is contained in $\bar{B}_d(x, h) \subset \bar{B}_d(x, \delta) \subset \tilde{U}_\ell$. Consequently, we have

$$|u(x) - u * \zeta^h(x)| \leq \sup_{y \in \bar{B}_d(x, h)} |R(x, y)| \leq C(d)h^2 \sup_{\substack{y \in \bar{B}_d(x, \delta) \\ 1 \leq i, j \leq d}} |\partial_{ij}^2 u(y)| = C(x, \delta)h^2,$$

where $C(d)$ is a constant which only depends on the dimension $d$ and we can set

$$C(x, \delta) = C(d) \sup_{\substack{y \in \bar{B}_d(x, \delta) \\ 1 \leq i, j \leq d}} |\partial_{ij}^2 u(y)| < +\infty, \qquad (20)$$

because $\bar{B}_d(x, \delta)$ is compact. This completes the proof of the lemma. $\square$

**Remark 3.2.** To obtain the bound (19) on $\mathbb{R}_+^d$, we need to suppose that u is defined and smooth on the larger domain $U_\delta$ (see formula (20) for the bound $C(x, \delta)$). Note that with the definition of $C(x, \delta)$, it is a straightforward consequence that for any K compact subset of $\mathbb{R}_+^d$, there exists $C(X, \delta)$ such that we have the uniform estimate:

$$\max_{x \in K} |u(x) - u * \zeta^h(x)| \leq C(X, \delta)h^2.$$
## 3.2 Quadrature formula and the SPH-like discretization

In this section we introduce a discretization of the integral approximation (18) of $u$. For that purpose, we first recall a theorem from [21], which provides a quadrature formula and an error estimate for functions in Sobolev spaces.

**Theorem 3.1.** Let $m$ be an integer such that $m > d$. Then there exists $C > 0$ such that, $\forall f \in W^{m,1}(\mathbb{R}^d)$ and $\forall \varepsilon > 0$:

$$ \left| \int_{\mathbb{R}^d} f(x)dx - \varepsilon^d \sum_{k \in \mathbb{Z}^d} f(x_k) \right| \le C \varepsilon^m \|f\|_{W^{m,1}(\mathbb{R}^d)}, \quad (21) $$

where $(x_k)_i = k_i \varepsilon$, $\forall i = 1, \dots, d$.

**Remark 3.3.** The condition $m > d$ is used to have the continuous embedding $W^{m,1}(\mathbb{R}^d) \hookrightarrow C^0(\mathbb{R}^d)$, which is needed to give a meaning to the point values $f(x_k)$.

Now we discretize the integral (18) with respect to $\xi \in \mathbb{R}^{d-1}$; this yields the following semi-discrete approximation of $u$:

$$ u(x) \approx \varepsilon^{d-1} \sum_{k \in \mathbb{Z}^{d-1}} \int_{-\ell}^{\infty} u(\Phi(s, \xi_k)) \zeta^h(x - \Phi(s, \xi_k))(\det J_{\Phi})(s, \xi_k)ds, \quad (22) $$

with $(\xi_k)_i = k_i \varepsilon$, $\forall i = 1, \dots, d-1$. Introducing $x_k(s) = \Phi(s, \xi_k)$, $\omega_k(s) = \varepsilon^{d-1}(\det J_{\Phi})(s, \xi_k)$ and $u_k(s) = u(x_k(s))$, the approximation (22) reads

$$ u(x) \approx \sum_{k \in \mathbb{Z}^{d-1}} \int_{-\ell}^{\infty} \omega_k(s) u_k(s) \zeta^h(x - x_k(s)) ds, $$

in the spirit of the SPH method.

Now we derive evolution equations for the quantities $x_k$, $u_k$ and $\omega_k$: this will provide a natural way to compute these quantities on a mesh discretizing the interval $(-\ell, +\infty)$.

**Proposition 3.1.** Let $u$ be the solution of equation (1), with boundary condition (2). Let $k \in \mathbb{Z}^{d-1}$. Then, the quantities $x_k(s) = \Phi(s, \xi_k)$, $u_k(s) = u(x_k(s))$ and $\omega_k(s) = \varepsilon^{d-1}(\det J_{\Phi})(s, \xi_k)$ satisfy the differential system

$$ x'_{k}(s) = \mathbf{a}(x_{k}(s)), \quad (23a) $$

$$ \omega'_{k}(s) = (\operatorname{div} \mathbf{a})(x_{k}(s))\omega_{k}(s), \quad \forall s > -\ell \quad (23b) $$

$$ u'_{k}(s) = -u_{k}(s)[(\operatorname{div} \mathbf{a})(x_{k}(s)) + a_{0}(x_{k}(s))] + S(x_{k}(s)), \quad (23c) $$

$$ x_k(0) = (0, \xi_k), \quad \omega_k(0) = \varepsilon^{d-1}\mathbf{a}_1(0, \xi_k), \quad u_k(0) = g(\xi_k). \quad (23d) $$

*Proof.* Equation (23a) follows from the definition of the characteristic system (8). Moreover, equation (23b) is a direct consequence of Liouville's formula:

$$ \frac{d}{ds}(\det J_{\Phi})(s, \xi) = (\operatorname{div} \mathbf{a})(\Phi(s, \xi))(\det J_{\Phi})(s, \xi). $$
Next, we prove (23c). By applying the chain rule, equation (1) and the characteristic system (8), one obtains:

$$
\begin{align*}
\frac{d}{ds}(u(x_k(s))) &= x'_k(s) \cdot \nabla u(x_k(s)) = \mathbf{a}(x_k(s)) \cdot \nabla u(x_k(s)) \\
&= -u(x_k(s)) \left[ (\mathrm{div} \, \mathbf{a})(x_k(s)) + a_0(x_k(s)) \right] + S(x_k(s)).
\end{align*}
$$

This completes the proof of the proposition.

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

We also consider the additional quantity $\rho_k = \omega_k u_k$. It is an easy consequence of proposition (3.1) that these functions satisfy the following result.

**Corollary 3.1.** *If $k \in \mathbb{Z}^{d-1}$, the function $\rho_k = \omega_k u_k$ is the solution of the Cauchy problem:*

$$
\rho'_{k}(s) + a_{0}(x_{k}(s))\rho_{k}(s) = \omega_{k}(s)S(x_{k}(s)), \quad \forall s > -\ell, \quad (24a)
$$

$$
\rho_k(0) = \epsilon^{d-1} \mathbf{a}_1(0, \xi_k) g(\xi_k). \qquad (24b)
$$

Systems (23) and (24) are the analogue of the equations in the SPH method. The difference here is that we replace the time variable by a variable s that parameterizes the space coordinate x₁ = Φ₁(s, ξ) > 0.

We end this section by deriving a priori estimates on the unknowns $\rho_k$ and $\omega_k$ that will be useful in the proof of convergence of our numerical scheme. For that purpose, we make the following additional hypothesis:

**Assumption 3.1.** The data **a**, and **g**, $a_0$, $S$ satisfy the regularity assumptions.

$$
\mathbf{a} \in C^1(\tilde{U}_\ell) \cap W^{2,\infty}(\tilde{U}_\ell), \quad g \in L^\infty(\mathbb{R}^{d-1}), \quad a_0, S \in C^1(\tilde{U}_\ell) \cap W^{1,\infty}(\tilde{U}_\ell).
$$

We obtain the following estimates on the weight $\omega_k$ and on the functions $\rho_k$.

**Lemma 3.2.** Suppose that **a** and **g**, $a_0$, $S$ satisfy the assumption 3.1. Let $(t, T) \in \mathbb{R}^2$ such that $-\ell < t < 0 < T < +\infty$. Then there exist constants $C_0(t, T)$, $C_1(t, T)$, $C_2(t, T)$, $C_3^0(t, T)$, $C_3^2(t, T) > 0$ independent of $\varepsilon$ and $h$ such that for all $k \in \mathbb{Z}^{d-1}$ and for all $s \in [t, T]$, one has:

$$
C_0(t, T) \varepsilon^{d-1} \le \omega_k(s) \le C_1(t, T) \varepsilon^{d-1}, \qquad (25)
$$

$$
|\omega_k''(s)| \le C_2(t, T) \varepsilon^{d-1}, \qquad (26)
$$

$$
|\rho_k(s)| \le C_3^0(t, T) \varepsilon^{d-1}, \quad |\rho_k''(s)| \le C_3^2(t, T) \varepsilon^{d-1}. \qquad (27)
$$

*Proof.* Let us fix $k \in \mathbb{Z}^{d-1}$ and $s \in [t, T]$. Then, by using assumption (7), one finds

$$
\alpha \varepsilon^{d-1} \le \omega_k(0) = \varepsilon^{d-1} \mathbf{a}_1(0, \xi_k) \le \beta \varepsilon^{d-1}.
$$
A direct integration of equation (23b) yields:

$$ \omega_k(s) = \omega_k(0) \exp \left( \int_0^s (\mathrm{div}\,\mathbf{a})(x_k(t)) dt \right) \le \beta \varepsilon^{d-1} \exp \left( (T-t) \|\mathrm{div}\,\mathbf{a}\|_{L^\infty(\tilde{U}_\ell)} \right) := C_1(t,T) \varepsilon^{d-1}. $$

Similarly, one has

$$ \omega_k(s) \ge \alpha \varepsilon^{d-1} \exp \left( (t-T) \|\mathrm{div}\,\mathbf{a}\|_{L^\infty(\tilde{U}_\ell)} \right) := C_0(t,T) \varepsilon^{d-1}. \quad (28) $$

From equations (23b) and (28), one deduces that

$$ \begin{aligned} |\omega_k''(s)| &= |(\mathrm{div}\,\mathbf{a})(x_k(s))^2 \omega_k(s) + \mathbf{a}(x_k(s)) \cdot \nabla(\mathrm{div}\,\mathbf{a})(x_k(s)) \omega_k(s)| \\ &\le 2 \|\mathbf{a}\|_{W^{2,\infty}(\tilde{U}_\ell)^d}^2 \omega_k(s) \le 2 \|\mathbf{a}\|_{W^{2,\infty}(\tilde{U}_\ell)^d}^2 C_1(t,T) := C_2(t,T) \varepsilon^{d-1}. \end{aligned} $$

Now let us show the first estimate of (27). By using (24) and for all $s \in (t, T)$, one finds that:

$$ \begin{aligned} |\rho_k(s)| &= \left| \left[ \varepsilon^{d-1} \mathbf{a}_1(0, \xi_k) g(\xi_k) + \int_0^s \omega_k(\mu) S(x_k(\mu)) \exp \left( \int_0^\mu a_0(x_k(\tau)) d\tau \right) d\mu \right] \exp \left( -\int_0^s a_0(x_k(\mu)) d\mu \right) \right| \\ &\le \varepsilon^{d-1} \exp \left( (T-t) \|a_0\|_{L^\infty(\tilde{U}_\ell)} \right) \left( \|g\|_{L^\infty(\tilde{U}_\ell)} \|a\|_{L^\infty(\tilde{U}_\ell)} + (T-t) C_1(t,T) \|S\|_{L^\infty(\tilde{U}_\ell)} \right) \\ &= C_3^0(t,T) \varepsilon^{d-1}. \end{aligned} $$

Then, from (24) and by using the bounds on $\rho_k$, one obtains:

$$ \begin{aligned} |\rho_k'(s)| &= |-a_0(x_k(s))\rho_k(s) + \omega_k(s)S(x_k(s))| \\ &\le \varepsilon^{d-1} \left( \|a_0\|_{L^\infty(\tilde{U}_\ell)} C_3^0(t,T) + \|S\|_{L^\infty(\tilde{U}_\ell)} C_1(t,T) \right) := C_3^1(t,T) \varepsilon^{d-1}. \end{aligned} $$

Finally, by differentiating (24) with respect to $s$, one deduces that

$$ \begin{aligned} |\rho_k''(s)| &= |x_k'(s) \cdot \nabla a_0(x_k(s)) \rho_k(s) + a_0(x_k(s)) \rho_k'(s) - \omega_k'(s) S(x_k(s)) - \omega_k(s) x_k'(s) \cdot \nabla S(x_k(s))| \\ &\le C_3^2(t,T) \varepsilon^{d-1} \end{aligned} $$

with

$$ C_3^2(t,T) = C(d) \left( \|a_0\|_{W^{1,\infty}(\tilde{U}_\ell)} \left( \|\mathbf{a}\|_{L^\infty(\tilde{U}_\ell)} C_3^0(t,T) + C_3^1(t,T) \right) + C_1(t,T) \|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_\ell)} \|S\|_{W^{1,\infty}(\tilde{U}_\ell)} \right). $$

This completes the proof of the lemma. $\square$
### 3.3 Discretization of the particle trajectories: the fully discrete scheme

Now we introduce the discretization with respect to $s \in (-\ell, +\infty)$. Indeed, we want to approximate $u * \zeta^h(x)$ for all $x \in \mathbb{R}_+^d$. For that purpose, we need to introduce a discretization of $\Omega_\ell = (-\ell, +\infty) \times \mathbb{R}^{d-1}$. Let us denote $n_\ell \in \mathbb{Z}$ such that $(n_\ell - 1)\Delta s < -\ell \le n_\ell \Delta s$, where $\Delta s > 0$. We will approximate $u$ with the formula:

$$ 
\begin{aligned}
u(x) &\approx \sum_{k \in \mathbb{Z}^{d-1}} \int_{-\ell}^{+\infty} \omega_k(s) u_k(s) \zeta^h(x - x_k(s)) ds, \\
&\approx \Delta s \sum_{i=n_\ell}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \omega_k(s_i) u_k(s_i) \zeta^h(x - x_k(s_i)), \\
&= \sum_{i=n_\ell}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \tilde{\omega}_k(s_i) u_k(s_i) \zeta^h(x - x_k(s_i)), \\
&= \sum_{i=n_\ell}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \tilde{\rho}_k(s_i) \zeta^h(x - x_k(s_i)) := \Pi_h^{\varepsilon, \Delta s}(u)(x),
\end{aligned} 
\quad (29) $$

where $s_i = i\Delta s$ and $\forall i \ge n_\ell$, $\tilde{\omega}_k(s_i) = \Delta s\omega_k(s_i)$ and $\tilde{\rho}_k(s_i) = \Delta s\rho_k(s_i)$. Note that $\tilde{\rho}_k$, $\tilde{\omega}_k$ and $\rho_k$, $\omega_k$ satisfy the same differential system: the only difference lies in the initial conditions that are multiplied by a factor $\Delta s$. In what follows, we remove the tilde from $\tilde{\rho}_k$ and $\tilde{\omega}_k$ and simply assume that $\rho_k$, $\omega_k$ now satisfy (23b),(24a) with the modified initial conditions

$$ \omega_k(0) = \Delta s \varepsilon^{d-1} \mathbf{a}_1(0, \xi_k), \quad \rho_k(0) = \Delta s \varepsilon^{d-1} \mathbf{a}_1(0, \xi_k) g(\xi_k). $$

Next we compute an approximation of $x_k(s_i)$, $u(x_k(s_i))$, $\omega_k(s_i)$ and $\rho_k(s_i)$ through the time discretization of the differential systems (23), (24) with an explicit Euler scheme:

$$ x_k^{j+1} = x_k^j + \Delta s \mathbf{a}(x_k^j), \quad (30a) $$

$$ \omega_k^{j+1} = \omega_k^j (1 + \Delta s (\operatorname{div} \mathbf{a})(x_k^j)), \quad (30b) $$

$$ \rho_k^{j+1} = \rho_k^j - \Delta s a_0 (x_k^j) \rho_k^j + \Delta s \omega_k^j S(x_k^j), \quad (30c) $$

$$ x_k^0 = (0, \xi_k), \quad u_k^0 = g(\xi_k), \quad \omega_k^0 = \Delta s \varepsilon^{d-1} \mathbf{a}_1(0, \xi_k), \quad \rho_k^0 = u_k^0 \omega_k^0. \quad (30d) $$

In order to simplify the discussion, we suppose for simplicity that the function $\operatorname{div} \mathbf{a}$ can be computed exactly. Now, we can define the fully discrete approximation of the function $u$:

$$ \tilde{u}_h^{\varepsilon, \Delta s}(x) := \sum_{j=n_\ell}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \omega_k^j u_k^j \zeta^h (x - x_k^j) = \sum_{j=n_\ell}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \rho_k^j \zeta^h (x - x_k^j). \quad (31) $$

Note that the sum is finite and is only composed of the indices $(j, k) \in [n_\ell, +\infty[ \times \mathbb{Z}^{d-1}$ such that $x_k^j \in \bar{B}_d(x, h)$.

We conclude this section by computing an error estimate associated to the Euler scheme for $x_k$, $\omega_k$, $\rho_k$. In the sequel, we will denote respectively $\lfloor \cdot \rfloor$ the floor function and $\lceil \cdot \rceil$ the ceiling function.
**Lemma 3.3.** Let $-l < t < 0 < T$ and set $n_t := \lfloor \frac{t}{\Delta s} \rfloor$, $N_T := \lfloor \frac{T}{\Delta s} \rfloor$. Assume $\mathbf{a} \in C^1(\tilde{U}_l) \cap W^{1,\infty}(\tilde{U}_l)$. Then there exists a constant $C_\Phi(t,T) > 0$ such that:

$$ \forall k \in \mathbb{Z}^{d-1}, \quad \max_{i \in [n_t, n_T]} |x_k(s_i) - x_k^i| \le C_\Phi(t, T) \Delta s. \qquad (32) $$

If, in addition, $\mathbf{a}$, $a_0$ and $S$ satisfy Assumption 3.1, then there exist constants $C_\omega(t,T) > 0$, $C_\rho(t,T) > 0$ independent of $\Delta s$, $\varepsilon$ and $h$ such that:

$$ \forall i \in [n_t, n_T], \quad \forall k \in \mathbb{Z}^{d-1}, \quad \begin{cases} |\omega_k(s_i) - \omega_k^i| \le C_\omega(t, T) \Delta s^2 \varepsilon^{d-1}, \\ |\rho_k(s_i) - \rho_k^i| \le C_\rho(t, T) \Delta s^2 \varepsilon^{d-1}. \end{cases} \qquad (33) $$

Moreover, there exist constants $C_\omega^0(t,T)$ and $C_\rho^0(t,T)$ such that for all $k \in \mathbb{Z}^{d-1}$

$$ \max_{i \in [n_t, n_T]} |\omega_k^i| \le C_\omega^0(t, T) \Delta s \varepsilon^{d-1}, \qquad (34) $$

$$ \max_{i \in [n_t, n_T]} |\rho_k^i| \le C_\rho^0(t, T) \Delta s \varepsilon^{d-1}. \qquad (35) $$

*Proof.* Suppose that $\mathbf{a} \in C^1(\tilde{U}_l) \cap W^{1,\infty}(\tilde{U}_l)$ and let $k \in \mathbb{Z}^{d-1}$. From equation (8), one deduces that

$$ x_k''(s) = (\mathbf{a} \cdot \nabla) \mathbf{a}(x_k(s)), \quad |x_k''(s)| \le \|\mathbf{a}\|_{L^\infty(\tilde{U}_l)} \|\nabla \mathbf{a}\|_{L^\infty(\tilde{U}_l)}, $$

where $\mathbf{a} \cdot \nabla = \sum_{i=1}^d \mathbf{a}_i \partial_{x_i}$. By computing a Taylor's expansion of $x_k$ on the interval $[s_j, s_{j+1}]$, one gets:

$$ x_k(s_{j+1}) = x_k(s_j) + \Delta s \mathbf{a}(x_k(s_j)) + \int_{s_j}^{s_{j+1}} x_k''(t)(s_{j+1} - t)dt. $$

One then obtains the estimate:

$$ |x_k(s_{j+1}) - x_k^{j+1}| \le (1 + \Delta s C(d) \|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_l)}) |x_k(s_j) - x_k^j| + C(d) \|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_l)}^2 \frac{\Delta s^2}{2}. $$

Assume that $0 \le j \le n_T$ and $x_k^0 = x_k(0)$. By applying a discrete Gronwall argument, one finds that

$$ \max_{0 \le j \le n_T} |x_k(s_j) - x_k^j| \le C(d) \|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_l)}^2 \frac{e^{C(d)\|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_l)}^T} - 1}{2} \Delta s. $$

By using a similar argument for $n_t \le j \le 0$, we get the error estimate (32) if we define $C_\Phi(t,T)$ as:

$$ C_\Phi(t, T) = C(d) \|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_l)}^2 \frac{e^{C(d)\|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_l)}(T-t)} - 1}{2}. \qquad (36) $$

Now we suppose that $\mathbf{a}, a_0, S$ satisfy the assumption 3.1 and we compute an error estimate on $e_k^j := \omega_k(s_j) - \omega_k^j$, for $n_t \le j \le n_T$. From equation (23b) and by applying Taylor's theorem, we get:

$$ \omega_k(s_j) = \omega_k(s_{j-1})(1 + \Delta s (\text{div } \mathbf{a})(x_k(s_{j-1}))) + \Delta s^2 \int_0^1 (1-t) \omega_k''(s_{j-1} + t\Delta s) dt. \qquad (37) $$
By applying lemma 3.2, one obtains

$$ \left| \int_{0}^{1} (1-t) \omega_k''(s_{j-1} + t \Delta s) dt \right| := \left| E_k^{j-1} \right| \leq \sup_{u \in [t,T]} |\omega_k''(u)| \leq C_2(t,T) \Delta s \varepsilon^{d-1}. $$

Then, by combining equation (30b) and (37), one obtains

$$ \begin{aligned} e_k^j &= e_k^{j-1} + \Delta s \left( (\operatorname{div} \mathbf{a})(x_k(s_{j-1})) \omega_k(s_{j-1}) - (\operatorname{div} \mathbf{a})(x_k^{j-1}) \omega_k^{j-1} \right) + \Delta s^2 E_k^{j-1} \\ &= e_k^{j-1} + \Delta s \left[ \left( (\operatorname{div} \mathbf{a})(x_k(s_{j-1})) - (\operatorname{div} \mathbf{a})(x_k^{j-1}) \right) \omega_k(s_{j-1}) - (\operatorname{div} \mathbf{a})(x_k^{j-1}) (\omega_k(s_{j-1}) - \omega_k^{j-1}) \right] \\ &\quad + \Delta s^2 E_k^{j-1} \\ &= e_k^{j-1} (1 - \Delta s (\operatorname{div} \mathbf{a})(x_k^{j-1})) + \Delta s \left( (\operatorname{div} \mathbf{a})(x_k(s_{j-1})) - (\operatorname{div} \mathbf{a})(x_k^{j-1}) \right) \omega_k(s_{j-1}) + \Delta s^2 E_k^{j-1}. \end{aligned} $$

As a consequence, by applying Lemma 3.2, one deduces that:

$$ |e_k^j| \leq |e_k^{j-1}| (1 + \Delta s \|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_\ell)}) + (C_\Phi(t,T)C_1(t,T)\|\mathbf{a}\|_{W^{2,\infty}(\tilde{U}_\ell)} + C_2(t,T)) \Delta s^3 \varepsilon^{d-1}. $$

Consequently, by using a discrete Gronwall argument and $e_k^0$, one obtains:

$$ |e_k^j| \leq (C_\Phi(t,T)C_1(t,T)\|\mathbf{a}\|_{W^{2,\infty}(\tilde{U}_\ell)} + C_2(t,T)) \Delta s^2 \varepsilon^{d-1} \frac{e^{\|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_\ell)} T} - 1}{\|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_\ell)}}. $$

One proceeds similarly for $n_t \le j \le 0$ and gets the first estimate of (33) by setting

$$ C_{\omega}(t,T) = (C_{\Phi}(t,T)C_{1}(t,T)\|\mathbf{a}\|_{W^{2,\infty}(\tilde{U}_{\ell})} + C_{2}(t,T)) \frac{e^{\|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_{\ell})}(T-t)} - 1}{\|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_{\ell})}}. $$

Now we consider the error estimate $l_k^j := \rho_k(s_j) - \rho_k^j$, for $n_t \le j \le n_T$. By using equation (24a) and computing a Taylor expansion of $\rho_k$ on the interval $[s_j, s_{j+1}]$, one finds:

$$ \begin{aligned} \rho_k(s_j) = \rho_k(s_{j-1}) &- \Delta s a_0(x_k(s_{j-1})) \rho_k(s_{j-1}) + \Delta s \omega(x_k(s_{j-1})) S(x_k(s_{j-1})) \\ &+ \Delta s^2 \int_0^1 (1-t) \rho_k''(s_{j-1} + t \Delta s) dt. \end{aligned} \quad (38) $$

From Lemma 3.2, one deduces that

$$ \left| \int_{0}^{1} \rho_k''(s_{j-1} + t \Delta s) dt \right| := L_k^{j-1} \leq \sup_{s \in [t,T]} |\rho_k''(s)| \leq C_3^2(t,T) \Delta s \varepsilon^{d-1}. $$

Then, by combining equation (30c) and (38), one has:

$$ \begin{aligned} |l_k^j| &\leq |l_k^{j-1} - \Delta s (a_0(x_k(s_{j-1}))\rho_k(s_{j-1}) - a_0(x_k^{j-1})\rho_k^{j-1})| \\ &\quad + |\Delta s (S(x_k(s_{j-1}))\omega_k(s_{j-1}) - S(x_k^{j-1})\omega_k^{j-1})| + \Delta s^2 L_k^{j-1} \\ &\leq |l_k^{j-1}| (1 + \Delta s \|a_0\|_{L^\infty(\tilde{U}_\ell)}) + \Delta s^2 C_\Phi(t,T) \|\nabla a_0\|_{L^\infty(\tilde{U}_\ell)} |\rho_k(s_{j-1})| \\ &\quad + \Delta s^2 C_\Phi(t,T) |\omega_k(s_{j-1})| \|\nabla S\|_{L^\infty(\tilde{U}_\ell)} + \Delta s |\omega_k(s_{j-1}) - \omega_k^{j-1}| \|S\|_{L^\infty(\tilde{U}_\ell)} + C_3^2(t,T) \Delta s^3 \varepsilon^{d-1} \\ &\leq (1 + \Delta s \|a_0\|_{L^\infty(\tilde{U}_\ell)}) |l_k^{j-1}| + C_4(t,T) \Delta s^3 \varepsilon^{d-1}, \end{aligned} $$
with

$$C_4(t, T) = \|a_0\|_{W^{1,\infty}(\tilde{U}_\ell)} C_\Phi(t, T) C_3^0(t, T) + \|S\|_{W^{1,\infty}(\tilde{U}_\ell)} (C_\Phi(t, T) C_1(t, T) + C_\omega(t, T)) + C_3^2(t, T).$$

Consequently, by the same computations as above, and by $l_k^0 = 0$, one obtains

$$|l_k^j| \le C_4(t, T) \frac{e^{\|a_0\|_{L^\infty(\tilde{U}_\ell)}(T-t)} - 1}{\|a_0\|_{L^\infty(\tilde{U}_\ell)}} \Delta s^2 \varepsilon^{d-1} := C_\rho(t, T) \Delta s^2 \varepsilon^{d-1}. \quad (39)$$

Finally, we derive estimates for $\omega_k^j$ and $\rho_k^j$ with $k \in \mathbb{Z}^{d-1}$ and $j \in [n_t, n_T]$: it is a straightforward computation to show that

$$|\omega_k^{j+1}| \le (1 + \Delta s \|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_\ell)}) |\omega_k^j|.$$

As a result, one obtains:

$$|\omega_k^j| \le \exp\left((T-t)\|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_\ell)}\right) \omega_k^0 \\ \le \exp\left((T-t)\|\mathbf{a}\|_{W^{1,\infty}(\tilde{U}_\ell)}\right) \|\mathbf{a}\|_{L^\infty(\tilde{U}_\ell)} \Delta s \varepsilon^{d-1} := C_\omega^0(t, T) \Delta s \varepsilon^{d-1}.$$

We finish the proof of this proposition by deriving an estimate on $\rho_k^j$. One has

$$\rho_k^{j+1} \le (1 + \Delta s \|a_0\|_{L^\infty(\tilde{U}_\ell)}) \rho_k^j + \|S\|_{L^\infty(\tilde{U}_\ell)} C_\omega^0(t, T) \Delta s^2 \varepsilon^{d-1}.$$

By a discrete Gronwall argument, one obtains $\rho_k^j \le C_\rho^0(t, T) \Delta s \varepsilon^{d-1}$ with

$$C_\rho^0(t, T) = \|S\|_{L^\infty(\tilde{U}_\ell)} C_\omega(t, T) \frac{e^{\|a_0\|_{L^\infty(\tilde{U}_\ell)}(T-t)} - 1}{\|a_0\|_{L^\infty(\tilde{U}_\ell)}}. \quad (40)$$

□

## 3.4 The main theorem

In this section, we prove the convergence of our numerical scheme.

**Theorem 3.2.** Let $d \in \mathbb{N}$ such that $d \ge 2$. Suppose that the following conditions are satisfied:

1. Smoothness of data: $\mathbf{a}$ and $a_0, S, g$ satisfy assumption 3.1 and we suppose that

$$\mathbf{a} \in C^{d+1}(\tilde{U}_\ell), \quad g \in C^d(\tilde{U}_\ell), \quad a_0, S \in C^d(\tilde{U}_\ell), \quad \zeta \in C^d(\mathbb{R}).$$

2. Characteristic Flow: the flow $\Phi$ satisfies Assumption 2.1. Moreover, we suppose that

$$\nabla (\Phi^{-1}) \in L^\infty(\tilde{U}_\ell).$$
3. Numerical parameters: $0 \le h \le \delta$ and $\varepsilon, \Delta s > 0$. Moreover, we assume

$$ \Delta s = \mathcal{O}_{\varepsilon \to 0}(\varepsilon), \quad \varepsilon = \mathcal{O}_{h \to 0}(h). $$

Then, for all compact sets $X \subset \mathbb{R}_+^d$, there exists $t_X < 0 < T_X$ and a compact set $\Xi_X \subset \mathbb{R}^{d-1}$ such that

$$ \bigcup_{x \in X} \bar{B}_d(x, \delta) \subset \Phi((t_X, T_X) \times \Xi_X). $$

Moreover, there exist constants $C_X^{(0)}, \dots, C_X^{(3)} > 0$ which do not depend on $h$, $\varepsilon$, $\Delta s$ such that, if $C_X^{(0)} \Delta s \le h$, the approximation $\tilde{u}_h^{\varepsilon, \Delta s}$ defined by equation (31) verifies:

$$ \sup_{x \in X} |u(x) - \tilde{u}_h^{\varepsilon, \Delta s}(x)| \le C_X^{(1)} h^2 + C_X^{(2)} \frac{\varepsilon^d}{h^d} + C_X^{(3)} \frac{\Delta s}{h}. \quad (41) $$

**Remark 3.4.** The error term is split into three contributions:

(i) The interpolation error: $u(x) - u * \zeta^h(x)$.

(ii) The quadrature error $u * \zeta^h(x) - \Pi_h^{\varepsilon, \Delta s}(u)(x)$.

(iii) The “time” discretisation error $\Pi_h^{\varepsilon, \Delta s}(u)(x) - \tilde{u}_h^{\varepsilon, \Delta s}(x)$ (stemming from the discretization of the characteristic equations).

The first assumption of the theorem implies that the flow $\Phi \in C^d(\Omega_\ell)$ and that the solution $u \in C^d(\tilde{U}_\ell)$. Now, we consider the first statement of the theorem. Let us first prove the following lemma.

**Lemma 3.4.** Let $x \in \mathbb{R}_+^d$, and suppose that $U_\delta = (-\delta, \infty) \times \mathbb{R}^{d-1} \subset \Phi(\Omega_\ell)$. Then there exists $-\ell < t_x < T_x$, and $\xi \in \mathbb{R}^{d-1}$, $\gamma_x > 0$ such that

$$ \bar{B}_d(x, \delta) \subset \Phi([t_x, T_x] \times \bar{B}_{d-1}(\xi_x, \gamma_x)), $$

(see Figure 3 for an illustration).

*Proof.* The set $\bar{B}_d(x, \delta)$ is a compact subset of $\mathbb{R}^d$ included in $\tilde{U}_\ell = \Phi(\Omega_\ell)$, and $\Phi : \Omega_\ell \to \tilde{U}_\ell$ is continuous. Thus $\Phi^{-1}(\bar{B}_d(x, \delta))$ is a compact subset of $\mathbb{R}^d$ included in $\Omega_\ell$. As a consequence there exists $-\ell < t_x < T_x$, and $\xi_x \in \mathbb{R}^{d-1}$, $\gamma_x > 0$ such that:

$$ \Phi^{-1}(\bar{B}_d(x, \delta)) \subset [t_x, T_x] \times \bar{B}_{d-1}(\xi, \gamma). $$

Since $\Phi$ is bijective, this shows the result. $\square$
<figure><img src="image_5.png" /><figcaption>Figure 3: Illustration of Lemma 3.4. The set $\bar{B}_d(x, \delta)$ is represented by the green circle. The set $\Phi([t_x, T_x] \times \bar{B}_{d-1}(\xi_x, \gamma_x))$ is in grey.</figcaption></figure>

The proof of the first statement of the theorem is then straightforward: let $K \subset \mathbb{R}_+^d$ be a compact set. For all $x \in K$, there exist $-\ell < t_x < T_x$, $\xi_x \in \mathbb{R}^{d-1}$, $\gamma_x > 0$ such that $\bar{B}_d(x, \delta) \subset \Phi((t_x, T_x) \times B_{d-1}(\xi_x, \gamma_x))$. We then choose

$$t_X = \min_{x \in K}(t_x), \quad T_X = \max_{x \in K}(T_x), \quad \Xi_X = \bigcup_{x \in K} B_{d-1}(\xi_x, \gamma_x)$$

and we get

$$\bigcup_{x \in X} \bar{B}_d(x, \delta) \subset \Phi((t_X, T_X) \times \Xi_X).$$

Now we prove the convergence of our numerical scheme.

*Proof.* Let $x \in \mathbb{R}_+^d$. We first split the difference $u(x) - \tilde{u}_h^{\varepsilon, \Delta s}(x)$ into four terms:

$$\begin{aligned}
u(x) - \tilde{u}_h^{\varepsilon, \Delta s}(x) = & u(x) - (u * \zeta^h)(x) \\
& + (u * \zeta^h)(x) - \sum_{j=n_\ell}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \rho_k(s_j) \zeta^h(x - x_k(s_j)) \\
& + \sum_{j=n_\ell}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} (\rho_k(s_j) - \rho_k^j) \zeta^h(x - x_k(s_j)) \\
& + \sum_{j=n_\ell}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \rho_k^j (\zeta^h(x - x_k(s_j)) - \zeta^h(x - x_k^j)) \\
& := I_h(x) + Q_{h,\varepsilon,\Delta s}(x) + F_{h,\Delta s}^{\rho} + F_{h,\Delta s}^{\zeta}(x).
\end{aligned}$$

As explained in Remark 3.4, $I_h$ is the interpolation error, $Q_{h,\varepsilon,\Delta s}(x)$ is the quadrature error, and $F_{h,\Delta s}^{\rho}$ and $F_{h,\Delta s}^{\zeta}$ are related to the “time” discretization errors of the characteristic differential equations.
We first consider the interpolation error $I_h(x)$. By applying Lemma 3.1, one has

$$|I_h(x)| = |u(x) - u * \zeta^h(x)| \leq C(x, \delta)h^2, \quad C(x, \delta) = C(d) \|u\|_{W^{2,\infty}(B_d(x, \delta))}.$$

As a result, we obtain

$$\max_{x \in X} |u(x) - u * \zeta^h(x)| \leq C_X^{(1)} h^2, \quad C_X^{(1)} = \max_{x \in X} C(x, \delta)$$

We get $C_X^{(1)} < +\infty$ from the smoothness assumption.

Next, we consider the quadrature error $Q_{h,\varepsilon,\Delta s}(x)$. By using the expression (18) and performing the change of variables $s = \frac{\Delta s}{\varepsilon}t$, one finds:

$$\begin{aligned} \left(u * \zeta^h\right)(x) &= \int_{-\ell}^{+\infty} \int_{\mathbb{R}^{d-1}} u(\Phi(s, \xi)) (\det J_{\Phi})(s, \xi) \zeta^h(x - \Phi(s, \xi)) d\xi ds \\ &= \int_{-\frac{\varepsilon \ell}{\Delta s}}^{+\infty} \int_{\mathbb{R}^{d-1}} f \circ g(t, \xi) d\xi dt. \end{aligned}$$

with

$$g(t, \xi) = \left(\frac{\Delta s}{\varepsilon} t, \xi\right), \quad f(s, \xi) = \frac{\Delta s}{\varepsilon} u(\Phi(s, \xi)) (\det J_{\Phi})(s, \xi) \zeta^h(x - \Phi(s, \xi)).$$

The term $Q_{h,\varepsilon,\Delta s}(x)$ is the error caused by the approximation of this integral in time by

$$\sum_{j=n_{\ell}}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \rho_k(s_j) \zeta^h(x - x_k(s_j)).$$

We apply Theorem 3.1 to the function $(t, \xi) \mapsto f \circ g(t, \xi)$ with a mesh of size $\varepsilon$. This function is $C^d(\Omega_\ell)$ (a is $C^{d+1}(U)$, thus $\Phi \in C^{d+1}(\Omega_\ell)$ and $\det J_\Phi \in C^d(\Omega_\ell)$). As a consequence, one finds:

$$|Q_{h,\varepsilon,\Delta s}(x)| \leq C(d) \varepsilon^d \|f\|_{W^{d,1}(\Omega_\ell)}.$$

We introduce $g(t, \xi) = (\frac{\Delta s}{\varepsilon}t, \xi)$, and $K_x := g^{-1}(\Phi^{-1}(\bar{B}_d(x, \delta)))$. For $\alpha \in \mathbb{N}^d$ we introduce the differential

$$D_{t, \xi}^{\alpha} := \left(\frac{\partial}{\partial t}\right)^{\alpha_1} \left(\frac{\partial}{\partial \xi_1}\right)^{\alpha_2} \cdots \left(\frac{\partial}{\partial \xi_{d-1}}\right)^{\alpha_d}.$$

As $\text{supp}(f) \subset K_x$, and by the change of variables $s = \frac{\Delta s}{\varepsilon}t$ we get

$$\begin{aligned} \|f \circ g\|_{W^{d,1}(\Omega_\ell)} &= \|f \circ g\|_{W^{d,1}(K_x)} = \sum_{p=0}^{d} \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\alpha|=p}} \int_{K_x} |D_{t, \xi}^{\alpha} f \circ g(t, \xi)| d\xi dt \\ &= \frac{\Delta s}{\varepsilon} \sum_{p=0}^{d} \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\alpha|=p}} \int_{\Phi^{-1}(\bar{B}_d(x, \delta))} \left(\frac{\Delta s}{\varepsilon}\right)^{\alpha_1} |D_{s, \xi}^{\alpha} (u(\Phi(s, \xi))(\det J_{\Phi})(s, \xi) \zeta^h(x - \Phi(s, \xi)))| \frac{\varepsilon}{\Delta s} d\xi ds \\ &\leq \|\det J_{\Phi}\|_{W^{d,\infty}(\Phi^{-1}(\bar{B}_d(x, \delta)))} \sum_{p=0}^{d} \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\alpha|=p}} \int_{\Phi^{-1}(\bar{B}_d(x, \delta))} \left(\frac{\Delta s}{\varepsilon}\right)^{\alpha_1} |D_{t, \xi}^{\alpha} (u(\Phi(s, \xi)) \zeta^h(x - \Phi(s, \xi)))| d\xi ds. \end{aligned}$$
Then, by performing the change of variable associated to $\Phi^{-1}$, one obtains:

$$ \begin{aligned} \mathcal{X} := & \sum_{p=0}^{d} \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\alpha|=p}} \int_{\Phi^{-1}(\bar{B}_d(x, \delta))} \left(\frac{\Delta s}{\varepsilon}\right)^{\alpha_1} |D_{t, \xi}^{\alpha} (u(\Phi(s, \xi)) \zeta^h(x - \Phi(s, \xi)))| \, d\xi \, ds \\ \le & \sum_{p=0}^{d} \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\alpha|=p}} \left(\frac{\Delta s}{\varepsilon}\right)^{\alpha_1} \sum_{\substack{\beta \in \mathbb{N}^d \\ |\beta| \le |\alpha|}} \int_{\bar{B}_d(x, \delta)} |D_y^\beta (u(y) \zeta^h(x-y))| P_{\alpha\beta}(|D\Phi|, \dots, |D^{|\alpha|}\Phi|) (\Phi^{-1}(y)) |\det J_\phi^{-1}(y)| \, dy \end{aligned} \quad (42) $$

where $P_{\alpha\beta}$ are polynomials with positive coefficients, and $D^k\Phi(y)$ is a vector containing all the partial derivatives $D^\alpha\Phi(y)$ with $|\alpha| = k$. As $\Delta s = \mathcal{O}_{\varepsilon\to 0}(\varepsilon)$, there exists $A > 0$ such that $\Delta s \le A\varepsilon$ for $\varepsilon$ sufficiently small. Then, by rearranging the terms, we get:

$$ \begin{aligned} \mathcal{X} \le & \| \det J_{\Phi}^{-1} \|_{L^{\infty}(\bar{B}_d(x, \delta))} \sum_{p=0}^{d} \sum_{\substack{\beta \in \mathbb{N}^d \\ |\beta|=p}} \int_{\bar{B}_d(x, \delta)} |D_y^\beta (u(y) \zeta^h(x-y))| \\ & \times \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\beta| \le |\alpha| \le p}} A^{\alpha_1} P_{\alpha\beta} (|D\Phi|, \dots, |D^{|\alpha|}\Phi|) (\Phi^{-1}(y)) \, dy \\ \le & \| \det J_{\Phi}^{-1} \|_{L^{\infty}(\bar{B}_d(x, \delta))} C \left( \|\Phi\|_{W^{d,\infty}(\bar{U}_\ell)} \right) \|u\|_{W^{d,\infty}(\bar{B}_d(x, \delta))} \|\zeta^h\|_{W^{d,1}(\bar{B}_d(x, \delta))}, \end{aligned} $$

where $C$ is a bounded function on any compact set of $[0, +\infty)$. By definition of $\zeta^h$, for any $p \in \mathbb{N}$ we have:

$$ \| (\zeta^h)^{(p)} \|_{L^1(\mathbb{R}^d)} = \frac{\mathrm{meas}(\mathbb{S}^{d-1}) \| \zeta^{(p)} \|_{L^1([0, +\infty), \mathcal{T}^{d-1}dt)}}{hp} := \frac{C_p}{hp}. \quad (43) $$

We assumed that $h \le \delta$ and thus we get $1/h^p \le \delta^{d-p}/h^d$. This leads to

$$ |Q_{h,\varepsilon,\Delta s}(x)| \le C_x \frac{\varepsilon^d}{h^d}, \quad (43) $$

where $C_x$ is a constant independent of $\varepsilon$, $\Delta s$, $h$. So, we get the second error term with $C_x^{(2)} = \max_{x \in \bar{X}} C_x$, which is finite with the smoothness assumption of the theorem 3.2.

Finally, we focus on the error terms $F_{h,\Delta s}^\rho(x)$ and $F_{h,\Delta s}^\zeta(x)$. By Lemma 3.4, there exists $n_x < N_x \in \mathbb{N}$ and $\xi_x \in \mathbb{R}^d$, $\gamma_x > 0$ such that

$$ \bar{B}_d(x, \delta) \subset \Phi([n_x \Delta s, N_x \Delta s] \times \bar{B}_{d-1}(\xi_x, \gamma_x)). $$

We note $t_x = n_x \Delta s$ and $T_x = N_x \Delta s$. Using Lemma 3.3, we have for $k \in \mathbb{Z}^{d-1}$, $n_x \le j \le N_x$:

$$ |\rho_k(s_j) - \rho_k^j| \le C_p (t_x, T_x) \Delta s^2 \varepsilon^{d-1}. $$

Moreover, by Lemma 3.2 (applied to $\tilde{\omega}_k = \Delta s \omega_k$ that we now denote by $\omega_k$),

$$ \omega_k(s_j) \ge C_0(t_x, T_x) \Delta s \varepsilon^{d-1}. $$
Consequently,

$$
\begin{align}
|F_{h,\Delta s}^{\rho}(x)| &= \left| \sum_{j=n_{\ell}}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \frac{\rho_k(s_j) - \rho_k^j}{\omega_k(s_j)} \omega_k(s_j) \zeta^h(x-x_k(s_j)) \right| \nonumber \\
&\le \frac{C_{\rho}(t_x, T_x)}{C_0(t_x, T_x)} \Delta s \sum_{j=n_{\ell}}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \omega_k(s_j) \zeta^h(x-x_k(s_j)). \tag{44}
\end{align}
$$

Now, we can write:

$$
\begin{align}
& \sum_{j=n_{\ell}}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \omega_k(s_j) \zeta^h(x-x_k(s_j)) \nonumber \\
&= \int_U \zeta^h(x-y) dy + \left| \int_U \zeta^h(x-y) dy - \sum_{j=n_{\ell}}^{+\infty} \sum_{k \in \mathbb{Z}^{d-1}} \omega_k(s_j) \zeta^h(x-x_k(s_j)) \right| \nonumber \\
&\le 1 + C(d) \|\zeta\|_{W^{d,1}([0, +\infty), t^{d-1}dt)} \frac{\varepsilon^d}{h^d} \le 1 + C(d) B^d \|\zeta\|_{W^{d,1}([0, +\infty), t^{d-1}dt)}. \tag{45}
\end{align}
$$

The first inequality is obtained by applying the same strategy as for the estimation of $\mathcal{Q}$ above, namely using Theorem 3.1. The second inequality follows from the fifth condition in Theorem 3.2. By combining estimates (44) and (45), one obtains:

$$
|F_{h,\Delta s}^{\rho}(x)| \leq C_x \Delta s \leq C_x \delta \frac{\Delta s}{h},
$$

for some constant $C_x$ independent of $\Delta s, \varepsilon$ and $h$. The second inequality is obtained by using $h \le \delta$. Next, we define $C_X^{(3,1)} = \max_{x \in X} \delta C_x$ which is finite thanks to the expressions (28) and (39) of $C_0(t_x, T_x)$ and $C_\rho(t_x, T_x)$.
We deduce that

$$
\sup_{x \in X} |F_{h, \Delta s}^{\rho}(x)| \le C_X^{(3,1)} \frac{\Delta s}{h}. \quad (46)
$$

Next, we compute the estimate of $F_{h,\Delta s}^{\zeta}(x)$. For all $n_x \le j \le N_x$, using Lemma 3.3, we have:

$$
|x_k(s_j) - x_k^j| \le C_{\Phi}(t_x, T_x) \Delta s.
$$

Then, one finds:

$$
\begin{align}
\left| \zeta^h(x - x_k(s_j)) - \zeta^h(x - x_k^j) \right| &\le \frac{|x_k(s_j) - x_k^j|}{h^{d+1}} \sup_{y \in [0,1]} |\zeta'(y)| \nonumber \\
&\le C_{\Phi}(t_x, T_x) \|\zeta'\|_{L^{\infty}([0,1])} \frac{\Delta s}{h^{d+1}} := C_{\zeta}(x) \frac{\Delta s}{h^{d+1}}. \tag{47}
\end{align}
$$

By using lemma 3.3 to bound $\rho_k^j$ and that $\zeta^h(x - \cdot)$ is compactly supported in $\bar{B}_d(x, h)$, one finds:

$$
|F_{h,\Delta s}^{\zeta}(x)| = \left| \sum_{(k,j) \in E_h(x)} \rho_k^j \left( \zeta^h(x - x_k(s_j)) - \zeta^h(x - x_k^j) \right) \right| \le C_{\zeta}(x) C_{\rho}^0(t_x, T_x) \frac{\Delta s^2}{h^{d+1}} \varepsilon^{d-1} |E_h(x)|, \quad (48)
$$
where

$$
E_h(x) = \{ (k, j) \in \mathbb{Z}^{d-1} \times \{n_l, \dots, +\infty\} ; x_k(s_j) \text{ or } x_k^j \in \bar{B}_d(x, h) \},
$$

and $|E_h(x)|$ denotes the cardinal of $E_h(x)$. In order to compute the cardinal of $E_h(x)$, first we introduce the set $E_{2h}^1(x)$, consisting of pairs $(k, j)$ such that $x_k(s_j)$ belongs to the larger set $\bar{B}_d(x, 2h)$:

$$
E_{2h}^{1}(x) = \{ (k, j) \in \mathbb{Z}^{d-1} \times \{n_{l}, \dots, +\infty\} ; x_{k}(s_{j}) \in \bar{B}_{d}(x, 2h) \}.
$$

We show that

$$
\Phi^{-1}(\bar{B}_d(x, 2h)) \subset \bar{B}_d \left( \Phi^{-1}(x), 2h \| (\nabla \Phi)^{-1} \|_{L^{\infty}(\tilde{U}_t)} \right).
$$

Let us note $(\sigma, \xi) := \Phi^{-1}(x)$, and let $(\sigma', \xi') \in \Phi^{-1}(\bar{B}_d(x, 2h))$. Then, $\nabla(\Phi^{-1})$ being bounded, one finds:

$$
\begin{aligned}
| (\sigma', \xi') - (\sigma, \xi) | &= | \Phi^{-1}(\Phi(\sigma', \xi')) - \Phi^{-1}(x) | \\
&\leq | \Phi(\sigma', \xi') - x | \| (\nabla \Phi)^{-1} \|_{L^{\infty}(\tilde{U}_t)} \leq 2h \| (\nabla \Phi)^{-1} \|_{L^{\infty}(\tilde{U}_t)},
\end{aligned}
$$

since $\Phi(\sigma', \xi') \in \bar{B}_d(x, 2h)$. Now, we define by $F_{2h}(x) \supset E_{2h}^1(x)$ the set:

$$
F_{2h}(x) = \{ (k, j) \in \mathbb{Z}^{d-1} \times \{-n_l, \dots, +\infty\} ; (s_j, \xi_k) \in \bar{B}_d \left( (\sigma, \xi), 2h \|(\nabla \Phi)^{-1}\|_{L^{\infty}(\tilde{U}_t)} \right) \},
$$

with $\xi_k$ such that $\Phi(s_j, \xi_k) = x_k(s_j)$. We have:

$$
\begin{align*}
F_{2h}(x) \subset G_{2h}(x) &:= \left[ \left\lfloor \frac{\sigma - 2h \|(\nabla \Phi)^{-1}\|_{L^{\infty}(\tilde{U}_t)}}{\Delta s} \right\rfloor, \left\lfloor \frac{\sigma + 2h \|(\nabla \Phi)^{-1}\|_{L^{\infty}(\tilde{U}_t)}}{\Delta s} \right\rfloor \right] \\
&\quad \times \prod_{1 \le i \le d-1} \left[ \left\lfloor \frac{\xi_i - 2h \|(\nabla \Phi)^{-1}\|_{L^{\infty}(\tilde{U}_t)}}{\varepsilon} \right\rfloor, \left\lfloor \frac{\xi_i + 2h \|(\nabla \Phi)^{-1}\|_{L^{\infty}(\tilde{U}_t)}}{\varepsilon} \right\rfloor \right], \\
|G_{2h}(x)| &\le \left( \frac{4h \|(\nabla \Phi)^{-1}\|_{L^{\infty}(\tilde{U}_t)}}{\varepsilon} + 2 \right)^{d-1} \left( \frac{4h \|(\nabla \Phi)^{-1}\|_{L^{\infty}(\tilde{U}_t)}}{\Delta s} + 2 \right) \\
&\le \left( 4 \|(\nabla \Phi)^{-1}\|_{L^{\infty}(\tilde{U}_t)} + 2B \max(A, 1) \right)^d \frac{h^d}{\Delta s \varepsilon^{d-1}} = C(x) \frac{h^d}{\Delta s \varepsilon^{d-1}}.
\end{align*}
$$

As a consequence:

$$
|E_{2h}^{1}(x)| \leq |F_{2h}(x)| \leq |G_{2h}(x)| \leq C(x) \frac{h^d}{\Delta s \varepsilon^{d-1}}. \quad (49)
$$

Let $k \in \mathbb{Z}^{d-1}$, $n_x \le j \le N_x$, one has:

$$
|x_k^j - x| \ge |x_k(s_j) - x| - |x_k^j - x_k(s_j)| \ge |x_k(s_j) - x| - C_{\Phi}(t_x, T_x) \Delta s.
$$
Now, we suppose that $C_{\Phi}(t_x, T_x)\Delta s \le h$. Then, if $x_k(s_j) \notin \bar{B}_d(x, 2h)$, we have $x_k^j \notin \bar{B}_d(x, h)$. Indeed:

$$|x_k^j - x| \ge 2h - C_x(t_x, T_x)\Delta s \ge h.$$

As a result, $E_h(x) \subset E_{2h}^1(x)$. Thus, we apply equation (49) and deduce that

$$|E_h(x)| \le C(x) \frac{h^d}{\Delta s \epsilon^{d-1}}.$$

We define $C_X^{(0)} = \max_{x \in X} C_{\Phi}(t_x, T_x)$ and $C_X^{(3,2)} = \max_{x \in X} C(d)C_{\zeta}(x)C_{\rho}^0(t_x, T_x)$. These quantities are finite by using estimates (36), (47), (40). Then, if $C_X^{(0)}\Delta s \le h$ and by using (48), one obtains

$$\sup_{x \in X} |F_{h, \Delta s}^{\zeta}(x)| \le C_X^{(3,2)} \frac{\Delta s}{h}. \qquad (50)$$

By combining (46) and (50), and setting $C_X^{(3)} = C_X^{(3,1)} + C_X^{(3,2)}$, we obtain

$$\sup_{x \in X} |F_{h, \Delta s}^{\rho}(x)| + \sup_{x \in X} |F_{h, \Delta s}^{\zeta}(x)| \le C_X^{(3)} \frac{\Delta s}{h},$$

which leads to the last error term in Theorem 3.2. This completes the proof of the theorem. $\square$

# 4 Numerical tests

In this section we present various numerical simulations to test our scheme. Equation (1) is solved on a bounded domain of $\mathbb{R}^2$, $\Omega := [0, L_x] \times [0, L_y]$. We first provide some details on the implementation of the scheme: choice of the kernel function $\eta^h$, boundary conditions and computational costs. Then we present numerical tests on the convergence of the scheme in the linear case. Finally, we use our scheme to carry out numerical simulation of a non linear landscape evolution model.

## 4.1 Practical implementation of the scheme

We provide some details on the implementation of our numerical scheme in two dimensions.

**The kernel function $\zeta^h$** We choose the cubic spline kernel, commonly used in the SPH method, as in [20]. For $x \in \mathbb{R}^2$, $h > 0$, this function is given by:

$$\zeta^h(x) = \frac{10}{7\pi h^2} \begin{cases} 1 - \frac{3}{2}q^2 + \frac{3}{4}q^3 & \text{if } 0 \le q \le 1, \\ \frac{1}{4}(2 - q^3) & \text{if } 1 \le q \le 2, \\ 0 & \text{otherwise,} \end{cases} \qquad (51)$$

where $q = |x|/(4h)$. This kernel function has compact support, is positive, and belongs to $C^2(\mathbb{R}^2)$.
**Boundary conditions** In order to simplify the discussion on the treatment of boundary conditions, we assume that the solution u of (1) is periodic in the transverse direction, that is

$$ \forall (x, y) \in \Omega, u(x, y + L_y) = u(x, y). $$

Therefore, $\Omega$ has two boundaries: $\Gamma_1 := \{0\} \times [0, L_y]$ and $\Gamma_2 := \{L_x\} \times [0, L_y]$. The equation (1) is a first-order partial differential equation: we only need to fix one boundary condition. To initialize the scheme, the particles $(\xi_k)_{k \in \{1, \dots, n\}}$ are located at the non-characteristic boundary $\Gamma_1$. We first compute the particle positions $x_k^j$ for $-1 \le j \ge m_k^0$ with the system (30), where:

$$ m_k^0 := \max \left\{ j \le -1; \left(x_k^j\right)_1 < -h \right\} + 1, $$

that is, the position of particles outside $\Omega$ that are located at a distance smaller than $h$ to $\Gamma_1$. Then, we compute the position of the particles inside the domain $\Omega$. These are the particles $x_k^j$ such that $0 \le j \le m_k^1$ with

$$ m_k^1 := \min \left\{ j \ge 1; \left(x_k^j\right)_1 > L_x \right\} - 1. $$

Finally, we compute the position of the particles $x_k^j$ such that $m_k^1 \le j \le m_k^2$ in the domain $]L_x, L_x+h] \times [0, L_y]$, located to the right of the right boundary $\Gamma_2$ with $m_k^2$ defined by:

$$ m_k^2 := \min \left\{ j \ge 1; \left(x_k^j\right)_1 > L_x + h \right\} - 1. $$

We refer to Fig 4 for an illustration of the particles positions.

**Computation of the approximate solution** Formula (31) and system (30) are used to compute the approximation $\tilde{u}_h^{\epsilon, \Delta s}$ of u in the domain $\Omega$, and in its extension $[-h, 0[ \times [0, L_y]$ and $]L_x, L_x+h] \times [0, L_y]$. We denote by $\tilde{\Omega} := [-h, L_x+h] \times [0, L_y]$. We compute the approximate solution at points x on a grid fine enough to provide a good insight into the variations of the approximate solution $\tilde{u}_h^{\epsilon, \Delta s}$. One has to compute the weights $\omega_k^j$ and $\rho_k^j$. For $1 \le k \le n$, $\omega_k^j$ verifies (30b). The function div a is not analytically known, so it needs to be approximated numerically. Consequently, we create a rectangular mesh of $\tilde{\Omega}$, and compute the derivatives of a by a finite difference scheme. Then, we use an linear interpolation to determine the value of div a at the particle positions. The computation of $\rho_k^j$ follows easily from (30c).

**Reduction of the computational cost** For any $x \in \Omega$, the approximate solution $\tilde{u}_h^{\epsilon, \Delta s}$ is given by:

$$ \tilde{u}_h^{\epsilon, \Delta s}(x) = \sum_{k=1}^{n} \sum_{j=m_k^0}^{m_k^2} \rho_k^j \zeta^h (x - x_k^j). \qquad (52) $$

In this sum, only a few terms are not zero, as the support of $\zeta^h(x - \cdot)$ is restricted to $\bar{B}_d(x, h)$. Therefore, computing every quantity $\zeta^h(x - x_k^j)$ for all values of $k$ and $j$ is unnecessary, and can be improved. For that purpose, we used the cell linked-lists method [9]. This method consists of two steps:
<figure><img src="image_6.png" /><figcaption>Figure 4: Illustration of particle positions in the scheme. Particles are generated along the time-discretization of the blue dotted curve (a characteristic curve). The red points are the extremal particles inside the domain, and the green one are the extremal particles outside it.</figcaption></figure>

1. First, a mesh grid with a spatial step h is created in the domain, together with an array representing this grid. Each particle is positioned in its corresponding cell inside this mesh grid.

2. If $x \in \Omega$, the particles that interact with this point (i.e., those corresponding to a non-zero term in the sum (52)) are located in the cell containing $x$ and its eight neighbours. Therefore, the sum is computed only on the particles that belongs to these nine cells, using the previous array.

The array of particle positions in the grid needs to be computed only once. Figure 5 represents the grid and the particles in the domain. The red particles interact with the point $x$, the green domain is composed of the nine neighboring cells. Green particles add unnecessary computations, but are not too many compared to the total number of particles in the domain. A direct implementation of (52) would require $N$ operations, where $N$ is the number of particles on the domain. Using the cell linked-list method, the computation of the sum requires approximately $9Nh^2 \ll N$ operations (equals to the number of particles in the nine neighboring cells to the point $x$). The approximate solution $\tilde{u}_h^{\epsilon, \Delta}$ is evaluated on a grid with $N_x = L_x/\Delta x$ and $N_y = L_y/\Delta y$ in the streamwise and cross-stream directions. The total computational cost is then approximately $9Nh^2N_xN_y$ with the cell linked list.

## 4.2 Numerical study of convergence

In this section we present numerical experiments using the particle scheme in the domain $\Omega$ and study the convergence. With a view to landscape evolution models simulation, we consider the following linear system:

$$ \left\{ \begin{array}{l} \mathrm{div}(au)(x,y) = r(x,y), \quad \forall (x,y) \in \tilde{\Omega} \\ u(0,y) = u_0(y), \end{array} \right. \qquad (53) $$
<figure><img src="image_7.png" /><figcaption>Figure 5: Grid of the linked-list method, of size h. The particles in the nine green neighboring cells of the blue point x are in color. The red particles are located in the support of y ↦ ζ^h(x - y) (represented by the red circle), the green particles are inside the green cells but outside this support.</figcaption></figure>

We define $\mathbf{a} := (\mathbf{a}_x, \mathbf{a}_y)$ with

$$ \left\{ \begin{array}{l} \mathbf{a}_x(x,y) = \tan \theta - \partial_x z(x,y), \\ \mathbf{a}_y(x,y) = -\partial_y z(x,y). \end{array} \right. \qquad (54) $$

For $(x, y) \in \tilde{\Omega} \setminus \Omega$, the vector field $\mathbf{a}$ is given by:

$$ \left\{ \begin{array}{l} \mathbf{a}_x(x,y) = \tan \theta, \\ \mathbf{a}_y(x,y) = 0. \end{array} \right. \qquad (55) $$

The function $(x, y) \mapsto \tan \theta(L_x - x) + z(x, y)$ will represent the bottom surface height in the landscape evolution model. When $z = 0$, this surface is a tilted plane, $\theta$ being the angle formed by this plane with the horizontal plane. This system describes the stationary distribution of fluid height $u$ on a given topography when a constant fluid flow is injected at the boundary $x = 0$.

In what follows, we solve system (53) with our particle method for various values of the parameters. The modeling parameters are given in Table 1, taken or estimated from an experiment [14] with a block of plaster set on a small tilted plane and eroded by a layer of water with a constant height. The numerical parameters are given in Table 2. We consider that equation (53) is written in a non-dimensional form, so that these parameters have no physical units.

In order to test the convergence of the scheme, three numerical parameters are tuned: $h$, $\varepsilon$ and $\Delta s$. The weights $w_0$ at $x = 0$ are the product of the distance between two particles, $\varepsilon$, and the distance traveled by a
<table><caption>Table 1: Parameters of numerical simulations</caption><tr><td>Length of the domain</td><td>L<sub>x</sub></td><td>1</td></tr><tr><td>Width of the domain</td><td>L<sub>y</sub></td><td>0.25</td></tr><tr><td>Plane inclination</td><td>θ</td><td>39°</td></tr><tr><td>Water height at the boundary</td><td>u<sub>0</sub>(y)</td><td>10<sup>-3</sup></td></tr><tr><td>Source term</td><td>r(x, y)</td><td>Variable</td></tr><tr><td>Surface perturbation</td><td>z(x, y)</td><td>Variable</td></tr></table>

<table><caption>Table 2: Parameters of the scheme.</caption><tr><td>Diameter of the kernel function</td><td>h</td><td>Variable</td></tr><tr><td>Distance between particles at the inbound boundary</td><td>ε</td><td>Variable</td></tr><tr><td>Discretization step for variable s</td><td>Δs</td><td>Variable</td></tr><tr><td>Weight of particles at the boundary</td><td>w<sub>0</sub></td><td>tan θ Δs ε</td></tr></table>

particle during one time step: a(0, y) Δs = tan θ Δs. This gives w₀ := tan θ Δs ε.

We consider the following values for h, ε and Δs:

$$h_{char} = \frac{1}{40}, \quad \varepsilon_{char} = \frac{h_{char}}{10}, \quad \Delta s_{char} = \frac{1}{200}.$$

The values of h, ε and Δs can vary in the following range:

$$h_{min} := 2^{-4}h_{char} \le h \le h_{max} := 2^{4}h_{char},$$

$$\varepsilon_{min} := 2^{-2}\varepsilon_{char} \le \varepsilon \le \varepsilon_{max} := 2^{3}\varepsilon_{char},$$

$$\Delta s_{min} := 2^{-5}\Delta s_{char} \le \Delta s \le \Delta s_{max} := 2^{5}\Delta s_{char}.$$

These bounds are chosen to ensure that the support of kernel is smaller than the domain, and that there is an average of more than one particle in the kernel function support. We present some error curves, depending on parameters. The error is the $\ell^\infty$ relative error, defined by:

$$\text{Error} = \frac{\max_{1 \le i \le N_x, 1 \le j \le N_y} \left( \tilde{u}_h^{\varepsilon, \Delta s}(idx, jdy) - u(idx, jdy) \right)}{\max_{1 \le i \le N_x, 1 \le j \le N_y} u(idx, jdy)}, \quad (56)$$

where $dx, dy > 0$ are given by $dxN_x = L_x$ and $dyN_y = L_y$.

### 4.2.1 Numerical convergence of the particle method

**Test on a flat surface.** First we study the convergence of the scheme in the simplest framework $z = 0$. The upstream water height is constant, of value one and the right hand side $r$ is set to 0. Thus the exact solution is a constant water height, equal to one. The numerical results performed with $h = h_{char}$, $ε = ε_{char}$ and $Δs = Δs_{char}$ produce a relative error with the exact solution of about 0.14%. As $u$ and $a$ are constants,
there is neither convolution error nor error in the computation of the $\omega_k^j$ and $\rho_k^j$. The only error comes from the quadrature error and according to Theorem 3.2, we have:

$$ \text{Error} \leq \frac{C\varepsilon^2}{h^2}. \qquad (57) $$

Thus the error grows with $\varepsilon/h$.

<figure><img src="image_8.png" /><figcaption>Figure 6: Flat surface test case: graphs of the $\ell^\infty$ relative errors (56) as functions of $h$, $\varepsilon$ and $\Delta s$, in log-log scale (using decimal logarithms). $z = 0, r = 0$.</figcaption></figure>

In figure 6a, we vary $h$ between $h_{min}$ and $h_{max}$, and the values of $\varepsilon$ and $\Delta s$ are varied accordingly with the formula: $\varepsilon = h/10$ and $\Delta s = 4\varepsilon$. Therefore, the number of particles in $\bar{B}_d(x, h)$ does not change. The computation times for the different values of $h$ remain almost constant. This is due to the linked-list method which keeps constant the number of computations when the number of particles in $\bar{B}_d(x, h)$ is constant. The only computation which takes more operations when $h$ becomes smaller is the computation of the linked-list grid.

For values of $h$ greater than 1/40 the error grows with $h$. Indeed, a more precise formula for (57) would
be

$$ \text{Error} \leq C\varepsilon^2 \left(1 + \frac{1}{h} + \frac{1}{h^2}\right), \qquad (58) $$

(this can be seen by ignoring the simplification made just before (43) in the proof of Theorem 3.2). The term $C(1+1/h)$ is responsible for this growth. However, this effect only appears when $h$ is large enough (see below). In practice the value of $h$ is small and this effect does not appear. Indeed, for values of $h$ smaller than 1/40 the error is constant, of about $10^{-2.9}$, which is consistent with the error bound (57).

Figure 6b shows the relative error depending on $h$, when $\varepsilon$ and $\Delta s$ are fixed. In accordance with equation (57) this error grows when $h$ becomes smaller, because the number of particles in the kernel diminishes with $h$. Figure 6c shows the error when $h = h_{char}$ and $\varepsilon$ varies with $\Delta s$. It shows that the error decreases if $\varepsilon$ and $\Delta s$ are decreased in agreement with equation (57). In figure 6d, $h = h_{char}$, $\varepsilon = \varepsilon_{char}$ and only $\Delta s$ varies. The error decreases with $\Delta s$, and reaches a minimum of about $10^{-5}$, when $\Delta s \le 2\varepsilon$. This is due to a simplification made in the proof of the theorem 3.2 where we used the inequality $\Delta s \le A\varepsilon$ in the proof of the quadrature error, after equation (42). Therefore a more precise formula for (57) contains terms that grow with $\Delta s/\varepsilon$, which is consistent with figure 6d. This shows that the spacing between quadrature points along $s$ and along $\xi$ should be roughly equal.

In figure 6c, we observe a minimum relative error on the order of $10^{-5.3}$, but which could be made smaller with smaller values of $\varepsilon$. Setting $h = h_{char}$, $\varepsilon = h/40$ and $\Delta s = 2\varepsilon$, the $\ell^\infty$ relative error is $7 \times 10^{-7}$. Of course, the computation time dramatically increases with such a small value of $\varepsilon$.

We have also considered a non-zero constant source term $r = 1/20$. The exact solution can be computed explicitly:

$$ \forall(x,y) \in \Omega, \quad u(x,y) = u_0 + \frac{rx}{\tan\theta}. $$

We found that the $\ell^\infty$ relative error when $h = h_{char}$, $\varepsilon = \varepsilon_{char}$, and $\Delta s = \Delta s_{char}$ is $2.9 \times 10^{-3}$.

**Test on a non-flat surface, transverse perturbation.** We consider the system with a surface varying in the transverse direction $y$:

$$ \forall(x,y) \in \Omega, \quad z_1(x,y) = \frac{u_0}{2} \cos \left( \frac{2\pi}{L_y} \left( y - \frac{L_y}{2} \right) \right). $$

Figure 7a shows a heat map of this function. Figure 7b presents the water height computed by the particle scheme, and particle positions are showed in Figure 7c. We can see that, when $x$ increases, particles tend towards the sides $y = 0$ and $y = 0.25$, where the surface height is smaller.

Next, we present some error curves. The exact solution of System (53) cannot be computed. Therefore we compute a reference solution with a finite-volume method as in [2], and a discretisation of 400 × 100 points; the error is computed using this reference solution. First, in Figure 8a, $h$ varies with $\varepsilon/h$ and $\Delta s/h$ is kept constant. For $1/3 \ge h \ge 1/20$, the error decreases with $h$, due to the first term of the error bound (41). To decrease the error, $\varepsilon$ or $\Delta s$ should be taken small compared to $h$. Then, Figure 8b shows the relative error depending on $h$, when $\varepsilon$ and $\Delta s$ are fixed. This error grows when $h$ becomes smaller, because of the terms in powers of $1/h$ in Equation (41). Figure 8c shows the error when $h = h_{ref}$ is fixed and $\varepsilon$, $\Delta s$ vary, and in Figure 8d only $\Delta s$ varies. In both Figures, the error decreases when $\varepsilon$ and $\Delta s$ diminish, because of the last two terms of Equation (41).
<figure><img src="image_9.png" /><figcaption>Figure 7: Non-flat surface, transverse perturbation: surface height (a), water height (b) and particle positions (c) in the domain $\Omega$, with $z = z_1$, $r = 0$ and $h = 4h_{char}$, $\epsilon = \epsilon_{char}$, $\Delta s = 2\Delta s_{char}$.</figcaption></figure>

**Test on a non-flat surface, longitudinal perturbation.** Now, we consider the system with a surface varying in $x$, in the longitudinal direction:

$$ \forall (x, y) \in \Omega, \quad z_2(x, y) = 20u_0 \cos\left(\frac{8\pi}{L_x}x\right). \qquad (59) $$

The profile of this surface is shown in Figure 9a. We can see that the function $x \mapsto z_2(x, 0)$ is strictly decreasing. Thus the vector field **a** does not vanish. Figure 9b presents the water height computed by the particle scheme. The quantitative behavior of the error curves is similar to the previous example 7, but the errors are approximately fifty times bigger. It can be explained by the fact that the minimum of the function $\det J_{\Phi}$ is quite small, and this term appears in the error estimate (41). As a comparison, the relative error is ten times smaller if the function $z$ is divided by two compared with (59), when $h = h_{char}$, $\epsilon = \epsilon_{char}$, and $\Delta s = 2\Delta s_{char}$.

### 4.2.2 Numerical simulations with dry areas

We have designed a numerical scheme under the assumption that the characteristic curves cover entirely the computational domain. Using the terminology of landscape evolution models, this means that there are no "dry areas" in the domain. Here, we discuss an extension of the numerical scheme in the presence of "dry areas".

**Definition 4.1.** A dry area is a connected component of $\tilde{U}_l \setminus \overline{\Phi(\Omega_l)}$. That is, it is a set of points that do not belong to a characteristic curve issued from the boundary $x_1 = 0$. We denote by $D_{\tilde{U}_l}$ the union of all the dry areas of $\tilde{U}_l$, that is:

$$ D_{\tilde{U}_l} = \tilde{U}_l \setminus \overline{\Phi(\Omega_l)}. $$

In the case of dry areas, we may construct a solution on the whole domain if we assume that the source term $S$ of (1) vanishes on $D_{\tilde{U}_l}$. We then consider the solution, given by:

* Equation (10) in $\mathrm{Int}(\Phi(\Omega_l))$,

* $u = 0$ in $D_{\tilde{U}_l}$.
<figure><img src="image_10.png" /><figcaption>Figure 8: Non-flat surface, transverse perturbation: graphs of the $\ell^\infty$ relative errors (56) as functions of $h$, $\epsilon$ and $\Delta s$, in log-log scale. $z = z_1$, $r = 0$.</figcaption></figure>

The solution may be discontinuous on $\partial(\Phi(\Omega_\ell)) = \partial(D_{\tilde{U}_\ell})$ and we do not define it on this set. See figure 10 for an example of such a situation: the dashed area is not reached by characteristic curves. It can be seen as a hill, whose top is the red cross. The left pointed edge of the dashed area is a mountain pass, and a vanishes at this point.

We have considered from a numerical point of view and chosen a surface with more variations than in previous section, see figure 11a. For $(x, y) \in \Omega$, let

$$z_3(x, y) = 2u_0 \cos \left( \frac{4\pi}{L_y} \left( y - \frac{L_y}{2} \right) \right).$$

The solution of system (53) presents some dry areas, on the highest areas of the surface. Figure 11b presents the water height computed from the particle scheme, and figure 11c the position of particles on the domain. We observe that they regroup in the lowest parts of the bottom surface, the so-called thalweg in hydrology terms. Between these thalwegs there are no particles, which means that the water height vanishes.
<figure><img src="image_11.png" /><figcaption>Figure 9: Non-flat surface, longitudinal perturbation: surface profile (a) and water height (b). $z = z_2$, $r = 0$, $h = 4h_{char}$, $\epsilon = \epsilon_{char}$, $\Delta s = 2\Delta s_{char}$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Characteristic curves on the domain $\mathbb{R}^{d-1} \times \mathbb{R}^+$. The vector field $\mathbf{a}$ vanishes at points indicated by red crosses.</figcaption></figure>

## 4.3 Application to a landscape evolution model

In this section we present numerical simulations of a landscape evolution model studied in [2]. The model consists of three partial differential equations on the bottom topography $z$, the fluid height $h$ and the sediment concentration $c$. Focusing on topography erosion, the time derivatives of $h, c$ are negligible and the system is given by: $\forall t \ge 0, \forall x \in \Omega \subset \mathbb{R}^2$,

$$ \text{div}(hv) = r(t,x), \qquad (60a) $$

$$ \text{div}(chv) = \rho_s e \left(\frac{h(t,x)}{H}\right)^m \left(\frac{|v(t,x)|}{V}\right)^n - \rho_s s \frac{c(t,x)}{c_{sat}}, \qquad (60b) $$

$$ \partial_t z = K \Delta z(t,x) - e \left(\frac{h(t,x)}{H}\right)^m \left(\frac{|v(t,x)|}{V}\right)^n + s \frac{c(t,x)}{c_{sat}}, \qquad (60c) $$

$$ v(t,x) = V(\tan \theta, 0) - V \nabla (h+z)(t,x). \qquad (60d) $$
<figure><img src="image_13.png" /><figcaption>Figure 11: Test case with dry areas: surface height (a), water height (b) and positions of the particles (c). $z = z_3$, $r = 0$ and $h = 4h_{char}$, $\epsilon = \epsilon_{char}$, $\Delta s = 2\Delta s_{char}$.</figcaption></figure>

The function $r$ represents the water source term, $\rho_s$ is the volumetric mass density of the sediments, $c_{sat}$ is the concentration of saturation for sediments, $e$ is the erosion rate, $s$ is the sedimentation rate, and $K$ is the constant of creep. The constants $H$ and $V$ are characteristic values for water height and water velocity.

Equation (60c) is discretized in time by an explicit Euler scheme. Denoting by $dt$ the time step of this scheme, we define $t_i := i dt$. Then, due to the expression of water velocity, both equations (60a) and (60b) are fully non linear. Therefore, these equations cannot be solved directly by the previous particle scheme, at $t \ge 0$ fixed. However, at each time step $t_i$, the quantities $z_i := z(t_i, \cdot)$ and $v_i := v(t_i, \cdot)$ differ only by an order $O(dt)$ term from those at the previous time step $t_{i-1}$. Consequently, we can propose a consistent in time discretization of these equations by evaluating these quantities at the previous time step in the expressions of the water velocity and of the source terms. We consider the following time discretization of (60). It leads to a linear system for the unknowns $h_{i+1}$, $c_{i+1}$, $z_{i+1}: \forall i \in \mathbb{N}$, $\forall x \in \Omega \subset \mathbb{R}^2$,

$$ \text{div}(h_{i+1}v_i)(x) = r_{i+1}(x), \qquad (61a) $$

$$ \text{div}(c_{i+1}h_{i+1}v_i)(x) = \rho_s e \left(\frac{h_i(x)}{H}\right)^m \left(\frac{|v_i(x)|}{V}\right)^n - \rho_s s \frac{c_i(x)}{c_{sat}}, \qquad (61b) $$

$$ z_{i+1}(x) = z_i(x) + dtK\Delta z_i(x) - e dt \left(\frac{h_i(x)}{H}\right)^m \left(\frac{|v_i(x)|}{V}\right)^n + s dt \frac{c_i(x)}{c_{sat}}, \qquad (61c) $$

$$ v_i(x) = V(\tan \theta, 0) - V \nabla (h_i + z_i)(x), $$

where $r_i(x) \approx r(t_i, x)$, $h_i(x) \approx h(t_i, x)$, and $c_i(x) \approx c(t_i, x)$. The first two equations of this system, equation. (61a) and equation (61b) can be solved by the particle method. The scheme is initialized with a constant elevation, that is $h_0(x, y) + z_0(x, y) = h_0$. Therefore, the initial water velocity is constant: $v_0(x, y) = V(\tan \theta, 0)$. The initial concentration $c$ is defined to be constant: $c_0(x, y) = c_0$.

**Choice of parameters** The choice of parameters is the same as in [2], where numerical simulations of system (60) were carried out with a finite-volume method. It is given in Table 3. The choice of numerical parameters is given in Table 4. The computation of $\nabla z$ and $\Delta z$ is performed with a finite difference discretization on a grid of size $N_x \times N_y$.

In the numerical simulations, the initial surface is a flat tilted plane with a small random perturbation, in the form of channels dug at random positions. This surface, seen from above is depicted in Figure 12. The
<table><caption>Table 3: Parameters of the model.</caption><tr><td>Length of the domain</td><td>L<sub>x</sub></td><td>40 cm</td></tr><tr><td>Width of the domain</td><td>L<sub>y</sub></td><td>10 cm</td></tr><tr><td>Characteristic water speed</td><td>V</td><td>1 m/s</td></tr><tr><td>Exponent of friction over h</td><td>m</td><td>1, 6</td></tr><tr><td>Exponent of friction over v</td><td>n</td><td>3, 2</td></tr><tr><td>Density of the sediments</td><td>ρ<sub>s</sub></td><td>2, 17.10<sup>6</sup> g/m<sup>3</sup></td></tr><tr><td>Concentration at saturation</td><td>c<sub>sat</sub></td><td>3, 17.10<sup>5</sup> g/m<sup>3</sup></td></tr><tr><td>Erosion speed</td><td>e</td><td>0, 5 mm/hour</td></tr><tr><td>Sedimentation speed</td><td>s</td><td>e/2000</td></tr><tr><td>Angle of the plane</td><td>θ</td><td>39°</td></tr><tr><td>Water height at the inbound boundary</td><td>h<sub>0</sub></td><td>0.5 mm</td></tr><tr><td>Concentration of sediments at the inbound boundary</td><td>c<sub>0</sub></td><td>317 g/m<sup>3</sup></td></tr></table>

<table><caption>Table 4: Parameters of numerical simulations</caption><tr><td>Number of points of the grid in the x direction</td><td>N<sub>x</sub></td><td>400</td></tr><tr><td>Number of points of the grid in the y direction</td><td>N<sub>y</sub></td><td>100</td></tr><tr><td>Time step for the equation on z</td><td>dt</td><td>4.5 s</td></tr><tr><td>Kernel diameter</td><td>h</td><td>L<sub>y</sub>/4</td></tr><tr><td>Distance between particles at the boundary</td><td>ε</td><td>h/10</td></tr><tr><td>Discretisation of variable s</td><td>Δs</td><td>0.005 s</td></tr></table>

computation time is about three hours for 1000 time steps, that is a final time of 1.25 hours in the simulation.

**Results of numerical simulations** The constant of creep *K* plays an important role in the behavior of the system. If *K* is large enough, the system is stable around the stationary solution of the flat plane, and when it is smaller the system is unstable, and channels are dug on the surface. The stability analysis of the system is performed in [2]. First, we show some numerical results with the creep parameter *K*:

$$K = K_e = \frac{5 \times 10^{-4}}{3600} m^2 s^{-1},$$

which correspond to a stable case.

Figures 13a and 13b show the surface height *z*, after 200 and 800 time steps. We observe that the initial perturbations have been flattened. Indeed, after 800 time steps the amplitude of the variations of the surface height is 2% of the initial amplitude. Then, Figure 13c shows the concentration of sediments in water, after 800 time steps. The concentration at the left of the domain is the same order as $c_0 = 317g/m^3$, and is larger on the right due to erosion. Figure 13d shows the water height after 800 time steps of the scheme. Its value is almost constant, around 0.5mm.

In Figure 14a and 14b we display the surface height, when $K = K_e/20$ and $K = K_e/50$. In these cases the system is unstable at some wavelengths, and we expect to observe a transverse instability. We can observe
<figure><img src="image_14.png" /><figcaption>Figure 12: Heat maps of the initial surface (A) and initial water height (B), in the domain $[0, L_x] \times [0, L_y]$.</figcaption></figure>

the formation of channels in the soil, in the flow direction. The amplitude of these channels is about 32% of the amplitude of the initial perturbations in the first case, and a bit more than 57% in the second case. Moreover, the depth of the channels is bigger when $K$ is smaller, in accordance with the stability results of [2]. After a sufficiently large number of time iterations with $K = K_e/50$, channels are deeper and water is concentrated inside these channels. Water height vanishes in highest areas, as in figure 15b. We observe in figure 15c that these areas do not contain particles. Particles gather inside the channels, so we could consider to merge some particles when they are too close, to improve the method. The surface height is represented in figure 15a, and has variations of 4 millimeters, which is 57 times bigger than the amplitude of the initial random perturbations.

In [2], system (3) was solved with a finite-volume method. This scheme can only be used with a positive water height. If the water height approaches zero at some points in the domain, some negative values of $h$ can appear in the computed solution. A finite-volume scheme that can handle dry areas could be designed, as in [7], but require a more sophisticated method. Our particle method is able to cope with the appearance of regions of vanishing water height (these areas simply do not contain any particles). Consequently, this particle method can be used on more general landscapes than the finite-volume method of [2], in particular, landscapes which may contain valleys and mountains.

**Convergence tests** Figures 16 and 17 show results of the simulations with $K = K_e/50$ at time 3.5 hours, for different values of the numerical parameters. In Figure 16, the time step of the landscape evolution model vary. We observe similar behavior between the initial figure 16a, when the time step is reduced 16b and when it is increased 16c.

Then, when the number of particles on the domain is reduced in the transverse or the longitudinal direction ($\epsilon$ or $\Delta$s bigger) as shown in Figures 17b and 17d the channels are shallower. On the other hand, if the number of particles on the domain is increased ($\epsilon$ or $\Delta$s smaller) in Figures 17c and 17e, the solution is quite similar to the reference solution 17a.
<figure><img src="image_15.png" /><figcaption>Figure 13: Heat maps of the surface height in millimeters $K = K_e$ after 200 time steps (a) and 800 time steps (b), sediment concentration in grams per cubic meters (c) and water height in millimeters (d) both after 800 time steps, in the domain $[0, L_x] \times [0, L_y]$, for $K = K_e$. Scales are given by the color bar to the right of each picture.</figcaption></figure>

# 5 Conclusion and perspectives

Motivated by applications in the simulation of landscape evolution models, we have proposed a new method to solve scalar equations in a divergence form that preserves the positivity of the solutions. This is a particle method and, to our knowledge, this is the first time it is applied to solve stationary problems. We have proven the convergence of the scheme under some regularity assumptions on the data (which in turn guarantee the smoothness of the solution), under some (classical) restrictions on the numerical parameter and under the assumption that the characteristic curves fulfill the whole computational domain. This later assumption is restrictive, but we observe that the numerical solution is well behaved even in the presence of dry areas. When the source term vanishes in these areas, one can extend the solution by assuming it vanishes also. We have performed several convergence tests which confirm our convergence result. We have also applied our numerical scheme in the context of landscape evolution models when fluid height vanishes and standard finite-volume methods may be numerically unstable.

Regarding applications to the simulation of landscape evolution model, some further developments are needed to make our scheme more efficient or applicable to more general landscapes. First, this scheme is
<figure><img src="image_16.png" /><figcaption>Figure 14: Heat maps of the surface height in millimeters for $K = K_e/20$ (A) and $K = K_e/50$ (B) after 1600 time steps, in the domain $[0, L_x] \times [0, L_y]$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: Heat maps of the surface height (A) water height (B) in millimeters and particle positions (C) for $K = K_e/50$ after 6400 time steps, in the domain $[0, L_x] \times [0, L_y]$.</figcaption></figure>

very costly compared to a finite-volume scheme. For the simulation of the landscape evolution model, the computation time of one time step is 200 times shorter with the finite-volume scheme than with the particle scheme, to obtain a similar precision on the solutions. This is due to the non-linear terms depending on $\nabla h$ which have to be computed by a sum over the neighboring particles. However, note that the computation of the particle method can be parallelized, as the computations for each trajectory can be run on independent processors. From a numerical point of view, we can also improve the accuracy of the scheme in the non-linear case: here, we proposed a linearized version of the problem at each time step which may generate a loss of accuracy. One may consider also solving directly the non-linear problem at each time step by a Newton method or a fixed point procedure. Note also that particles computed as solutions to the linear system tend to gather in the lowest parts of the surface. If they are too close, then the gradient of $h$ can become singular, which leads to considerable errors in the computation of the solution. Consequently, the time steps have to be taken small enough to obtain good convergence. One could also consider merging particles that are too close to prevent the formation of singularities.
<figure><img src="image_18.png" /><figcaption>Figure 16: Heat maps of the surface height in millimeters for $K = K_e/50$ at time 3.5 hours, in the domain $[0, L_x] \times [0, L_y]$. Parameters of table 4 (A), smaller time step (B), bigger time step (C).</figcaption></figure>

A Proof of proposition 2.1

*Proof.* Items (i) and (ii) of the theorem are proved by a direct application of the Cauchy-Lipschitz theorem. Then we prove Item (iii). First, as **a** is P-Lipschitz continuous, by (7) and the fact that $\Phi(0, \xi) = (0, \xi)$, we have for $t \ge 0$:

$$ \begin{aligned} |\Phi(t, \xi) - \Phi(0, \xi)| &\le \int_0^t |\mathbf{a}(\Phi(s, \xi))| ds \le t\beta + \int_0^t |\mathbf{a}(\Phi(s, \xi)) - \mathbf{a}(\Phi(0, \xi))| ds \\ &\le t\beta + P \int_0^t |\Phi(s, \xi) - \Phi(0, \xi)| ds. \end{aligned} $$

A similar inequality can be shown for $t \le 0$. Then, by Gronwall's lemma, denoting by $\Phi_1(t, \xi)$ the first component of $\Phi(t, \xi)$, and using that $\Phi_1(0, \xi) = 0$ we get:

$$ |\Phi_1(t, \xi)| \le \frac{\beta}{P} (e^{P|t|} - 1). $$
<figure><img src="image_19.png" /><figcaption>(a) $dt = 4.5$ s, $\epsilon = h/10$ cm and $\Delta s = 0.005$ s.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) $dt = 4.5$ s, $\epsilon = h/5$ cm and $\Delta s = 0.005$ s.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(c) $dt = 4.5$ s, $\epsilon = h/20$ cm and $\Delta s = 0.005$ s.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(d) $dt = 4.5$ s, $\epsilon = h/10$ cm and $\Delta s = 0.01$ s. Different color bar scale.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(e) $dt = 4.5$ s, $\epsilon = h/10$ cm and $\Delta s = 0.0025$ s.</figcaption></figure>

Figure 17: Heat maps of the surface height in millimeters for $K = K_e/50$ at time 3.5 hours (2800 time steps), in the domain $[0, L_x] \times [0, L_y]$. Scales are the same for all pictures except (D) which is larger. Parameters of table 4 (A), fewer/more particles in the transverse direction (B)/(C), fewer/more particles in the longitudinal direction (D)/(E).
Thus for all $|t| \le \ell := \frac{1}{\bar{p}} \ln(1 + \frac{LP}{\bar{p}})$, we have $|\Phi_1(t, \xi)| \le L$. This implies that $\Phi_1(t, \xi) \ge -L$ for $t \in (-\ell, 0]$. Then for $t > 0$ we show that $\Phi_1(t, \xi) > 0$. Indeed, because of (7) we have $\Phi_1(t, \xi) > 0$ for $t \in (0, \varepsilon_0]$ with $\varepsilon_0 > 0$ small. By way of contradiction, suppose that $t_1 > 0$ is the smallest time such that $\Phi_1(t_1, \xi) = 0$. So we have $\Phi_1(t, \xi) > 0$ for $t \in [t_1 - \varepsilon_1, t_1)$ with $\varepsilon_1 > 0$ small enough. This contradicts the fact that

$$ \partial_t \Phi_1(t_1, \xi) = \mathbf{a}_1(\Phi(t_1, \xi)) = \mathbf{a}_1(0, \Phi_2(t_1, \xi), \dots, \Phi_d(t_1, \xi)) > 0, $$

by (7). Hence $\Phi_1(t, \xi) > -L$, $\forall (t, \xi) \in (-\ell, +\infty) \times \mathbb{R}^{d-1}$. This shows the item (iii) (which is also Assumption 2.1 (i)).

Finally we prove Item (iv). We first show that with the sole property (7), there exists $\delta > 0$ such that $(-\delta, \delta) \times \mathbb{R}^{d-1} \subset \Phi(\Omega_\ell)$. Indeed, since $\mathbf{a}$ is $P$-Lipschitz continuous and with (7), setting $\delta = \frac{\alpha}{2} \min(\frac{1}{\bar{p}}, \ell)$ we have:

$$ \mathbf{a}_1(x_1, \eta) \ge \alpha - P|x_1| > \frac{\alpha}{2}, \quad \forall (x_1, \eta) \in (-\delta, \delta) \times \mathbb{R}^{d-1}. $$

Now, for such $(x_1, \eta)$, define $\chi(t)$ as the solution of

$$ \chi'(t) = \mathbf{a}(\chi(t)), \quad \chi(0) = (x_1, \eta). \qquad (62) $$

Without loss of generality, suppose $x_1 \in (-\delta, 0)$. We have

$$ \chi_1(t_1) - \chi_1(t_0) = \int_{t_0}^{t_1} \mathbf{a}_1(\chi(s))ds \ge \frac{\alpha}{2}(t_1 - t_0) > 0, \qquad (63) $$

for all $t_1 > t_0 > 0$ provided that $\chi_1(s) \in (-\delta, \delta)$, $\forall s \in [0, t_1]$. Therefore $\chi_1$ is strictly increasing on some interval $(0, T_1]$, with $T_1 > 0$. Since $0 \in (-\delta, \delta)$, there exists a time $t_2 > 0$ such that $0 < t_2 < T_1$ and $\chi_1(t_2) = 0$. Furthermore, using (63) with $(t_0, t_1) = (0, t_2)$, $t_0$ satisfies

$$ \delta > -x_1 = \int_0^{t_2} \mathbf{a}_1(\chi(s))ds > \frac{\alpha t_2}{2}, $$

hence $t_2 < 2\delta/\alpha \le \ell$. We denote $\chi(t_2) = (0, \xi)$. Then, by the uniqueness in the Cauchy-Lipschitz theorem we have $\Phi(t - t_2, \xi) = \chi(t)$ for all $t \in (0, t_2)$. Thus for $t = 0$, one has: $\Phi(-t_2, \xi) = \chi(0) = (x_1, \eta)$. This shows that $(-\delta, \delta) \times \mathbb{R}^{d-1} \subset \Phi(\Omega_\ell)$.

Now we assume that $\mathbf{a}_1(x_1, \xi) \ge \alpha > 0$ for all $(x_1, \xi) \in U$ and consider $(x_1, \eta) \in [\delta, +\infty) \times \mathbb{R}^{d-1}$. Let us consider the solution $\chi(t)$ of (62). For $t > 0$, we have:

$$ |\chi(-t) - \chi(0)| = |\int_{0}^{t} (\mathbf{a}(x_1, \eta) + \mathbf{a}(\chi(-s)) - \mathbf{a}(\chi(0))) ds| \le |\mathbf{a}(x_1, \eta)|t + P \int_{0}^{t} |\chi(-s) - \chi(0)|ds. $$

Hence, by Gronwall's lemma,

$$ |\chi(-t) - \chi(0)| \le \frac{|\mathbf{a}(x_1, \eta)|}{P} (e^{Pt} - 1). $$

This shows that $\chi(-t)$ is defined for all $t > 0$ as long as $\chi(-t) \ge -L$. Thus, there exists $t_0 > 0$ such that $\chi_1(-t_0) = 0$. Setting $\chi(-t_0) = (0, \xi)$, we have $\Phi(t_0, \xi) = (x_1, \eta)$ which finishes the proof of Assumption (2.1) (ii). $\square$
References

[1] B Ben Moussa and JP Vila. Convergence of SPH method for scalar nonlinear conservation laws. *SIAM Journal on Numerical Analysis*, 37(3):863-887, 2000.

[2] Julie Binard, Pierre Degond, and Pascal Noble. Well-posedness and stability analysis of a landscape evolution model. *Journal of Nonlinear Science*, 34(1):20, 2024.

[3] Alex Chen, Jérôme Darbon, Giuseppe Buttazzo, Filippo Santambrogio, and Jean-Michel Morel. On the equations of landscape formation. *Interfaces and Free Boundaries*, 16(1):105-136, 2014.

[4] Alex Chen, Jérôme Darbon, and Jean-Michel Morel. Landscape evolution models: A review of their fundamental equations. *Geomorphology*, 219:68-86, 2014.

[5] Alexandre Joel Chorin and Peter S Bernard. Discretization of a vortex sheet, with an example of roll-up. *Journal of Computational Physics*, 13(3):423-429, 1973.

[6] WM Davis. The convex profile of bad-land divides. *Science*, (508):245-245, 1892.

[7] Olivier Delestre. *Simulation du ruissellement d'eau de pluie sur des surfaces agricoles*. PhD thesis, Université d'Orléans, 2010.

[8] Jean Dieudonné. *Eléments d'analyse: tome I*. Gauthier-Villars, 1969.

[9] JM Domínguez, AJC Crespo, M Gómez-Gesteira, and JC2869628 Marongiu. Neighbour lists in smoothed particle hydrodynamics. *International Journal for Numerical Methods in Fluids*, 67(12):2026-2042, 2011.

[10] Martha W Evans, Francis Harvey Harlow, and Eleazer Bromberg. The particle-in-cell method for hydrodynamic calculations. 1957.

[11] Joep HM Evers, Iason A Zisis, Bas J van der Linden, and Manh Hong Duong. From continuum mechanics to sph particle systems and back: Systematic derivation and convergence. *ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik*, 98(1):106-133, 2018.

[12] Grove Karl Gilbert. The convexity of hilltops. *The Journal of Geology*, 17(4):344-350, 1909.

[13] Robert A Gingold and Joseph J Monaghan. Smoothed particle hydrodynamics: theory and application to non-spherical stars. *Monthly notices of the royal astronomical society*, 181(3):375-389, 1977.

[14] Adrien Guérin, Julien Derr, Sylvain Courrech Du Pont, and Michael Berhanu. Streamwise dissolution patterns created by a flowing water film. *Physical Review Letters*, 125(19):194502, 2020.

[15] Alan D Howard and Gordon Kerby. Channel changes in badlands. *Geological Society of America Bulletin*, 94(6):739-752, 1983.

[16] Anthony Leonard. Computing three-dimensional incompressible flows with vortex elements. *Annual Review of Fluid Mechanics*, 17(1):523-559, 1985.
[17] Steven J Lind, Benedict D Rogers, and Peter K Stansby. Review of smoothed particle hydrodynamics: towards converged lagrangian flow modelling. *Proceedings of the royal society A*, 476(2241):20190801, 2020.

[18] Leon B Lucy. A numerical approach to the testing of the fission hypothesis. *Astronomical Journal*, vol. 82, Dec. 1977, p. 1013-1024.

[19] S Mas-Gallic and PA Raviart. A particle method for first-order symmetric systems. *Numerische Mathematik*, 51(3):323-352, 1987.

[20] Joe J Monaghan. Smoothed particle hydrodynamics. *Annual review of astronomy and astrophysics*, 30(1):543-574, 1992.

[21] Pierre-Arnaud Raviart. An analysis of particle methods. In *Numerical Methods in Fluid Dynamics: Lectures given at the 3rd 1983 Session of the Centro Internationale Matematico Estivo (CIME) held at Como, Italy, July 7-15, 1983*, pages 243-324. Springer, 2006.

[22] Louis Rosenhead. The formation of vortices from a surface of discontinuity. *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character*, 134(823):170-192, 1931.