Julie Binard
email: binard@insa-toulouse.fr

Pierre Degond
email: pierre.degond@math.univ-toulouse.fr

Pascal Noble
email: noble@insa-toulouse.fr

Well-posedness and stability analysis of a landscape evolution model

Keywords: Pattern formation, instabilities, Landscape evolution model, erosion by water, stream incision law, sedimentation Mathematics Subject Classification 35M30, 35Q86, 35B35, 35B36, 35A01, 86-10

In this paper, we study a system of partial differential equations modeling the evolution of a landscape in order to describe the mechanisms of pattern formations. A ground surface is eroded by the flow of water over it, either by sedimentation or dilution. We consider a model, composed of three evolution equations, one on the elevation of the ground surface, one on the fluid height and one on the concentration of sediments in the fluid layer. We first establish the well-posedness of the system in short time, and under the assumption that the initial fluid height does not vanish. Then, we focus on pattern formation in the case of a film flow over an inclined erodible plane. For that purpose, we carry out a spectral stability analysis of constant state solutions in order to determine instability conditions and identify a mechanism for pattern formations. These patterns, which are rills and gullies, are the starting point of the formation of rivers and valleys in landscapes. Finally, we carry out some numerical simulations of the full system in order to validate the spectral instability scenario, and determine the resulting patterns.

Introduction

The modeling of landscape evolution under the effect of water flow has received increasing attention these last decades. Several physical phenomena must be considered in order to describe the erosion and sedimentation processes, which occur when water flows over an erodible surface. The erosion is the removal of sediments from the soil by a fluid, and the sedimentation is the inverse process, when the sediments settle on the surface. Heuristically, the intensity of the erosion process is strongly related to the flow rate of the fluid, the erosion rate is higher when the flow rate is higher. An other process that affects landscape evolution is the creep effect, which tends to smooth the bottom surface on time scales that are much larger than the ones involved in sediment transport. This effect has multiple causes, such as gravitational forces acting on the soil, and corrosion and dilatation caused by physical and chemical factors. It can be described as a simple diffusion process of the soil. The creep effect was introduced in 1892, in [START_REF] Davis | The convex profile of bad-land divides[END_REF]. It was used in [START_REF] Gilbert | The convexity of hilltops[END_REF] to explain the convexity of hilltops profile. Later, in 1963 a derivation of the soil creep effect as a limit of a stochastic process had been done in [START_REF] Culling | Soil creep and the development of hillside slopes[END_REF]. The stochastic effect models soil particles which follow a random motion, under the constrain of gravity. In this article, Culling writes "Soil creep -so gradual as to be imperceptible -appears to be the result of the persistent effect of molecular and macromolecular forces tending to displace the soil particles.", and he describes the soil flow as a quasi-viscous fluid.

The description of the geological processes that occur in landscape evolution was initiated in [START_REF] Gilbert | Report on the Geology of the Henry Mountains[END_REF]. In this book, the author gives the fundamental principles of landscape evolution and explains why the profiles of stream beds are concave upwards. This concavity property had been illustrated later in [START_REF] Culling | Analytical theory of erosion[END_REF], with a mathematical description of the erosion of a slope. The slope evolution is described with a reaction diffusion equation, the exponential source term representing the erosion rate. The landscapes evolve mainly in function of these two competitive factors: creep on the hilltops, and stream incision on the lower slopes. Indeed, in the upper slopes the water flow is weak and dispersive, thus the creep effect predominates. Downward, the effect of shear stress of the flow becomes dominant, and bedload transport process increases. This leads to the formation of stream beds and valley, the profiles being concave. Since then, the complexity of the landscape evolution models increased: see [START_REF] Chen | On the equations of landscape formation[END_REF] and [START_REF] Chen | Landscape evolution models: A review of their fundamental equations[END_REF] for a review of these models. The process of erosion and transport of sediments can be described by two different laws, the alluvial transport law or the stream incision law, which depends on the nature of the soil. When the soil is covered by alluviums, the sediments are directly transported by water flux, and the amount of sediments moved by the water flux in a given time q s follows a law that depends on the water discharge q. This transport discharge law is described, for example, in [START_REF] Smith | Stability and the conservation of mass in drainage basin evolution[END_REF] and is given by

q s = kq m |∇z| n (1)
where k, m, n are constants. This is called the transport limited case.

On the other side, the stream incision law is used when the surface is bedrock, because in this case the sediment transport is limited by the resistance of the bedrock to the shear stress caused by the water flux. This is modeled by an erosion source term in the equation of surface height evolution, as described in [START_REF] Howard | Channel changes in badlands[END_REF]. The sediments removed from the soil by the erosion are supposed to be dissolved in water, thus are transported by the water flux. This stream incision law is described below, in Section 2. This case, that we study in this paper is called the detachment limited case.

In this paper, we study a system of partial differential equations describing the erosion of the soil by water, which had been proposed in [START_REF] Chen | On the equations of landscape formation[END_REF], Section 4. It includes the modeling of the water flow, the erosion of the surface by water, the transport and deposition of sediments, and the creep effect. The characteristic fluid velocity is usually much larger than the erosion rate: in order to describe pattern formation, we only consider large scale fluctuations of the fluid velocity. As a result, we assume that the fluid velocity is proportional to the gradient of the free surface elevation.The other physical principles taken in account are:

• The conservation law for water and sediments dissolved in water,

• The stream incision law: the erosion grows with the water speed and the water height,

• The sedimentation rate is proportional to the concentration of sediments in water,

• The creep effect: the soil is subject to a diffusion process.

Other effects, such as infiltration, vegetation, wind, or ice formation are neglected. Moreover, we suppose that the bottom surface is constituted by one type of sediments. The system studied in this article models the time evolution of the soil and of the fluid. This system is composed by the following three partial differential equations:

       ∂ t h -div(h∇(h + z)) = r(t, x), ∂ t z = K∆z + sc -eh m |v| n , ∂ t (ch) + div(chv) = eh m |v| n -sc.
(

) 2
The first equation is a mass conservation law for the fluid and describes the evolution of the fluid height, denoted by h, and the fluid is transported at speed v := -∇(h + z). The term r is a source term representing an exterior source of water, like the rain. The second equation models the time evolution of the bottom topography (denoted by z). The constant K is the constant of creep whereas e and s are respectively the constants for the incision law and the sedimentation rate. As in [START_REF] Howard | Channel changes in badlands[END_REF], we will suppose that the incision law depends on a power of the norm of the water velocity, and on a power of the water height. Thanks to this hypothesis, as long as the water height and velocity do not vanish, the erosion rate is positive. However, in practice the erosion starts if the water velocity is high enough to break the cohesion of the soil. Thus a threshold effect could be introduced in the model. This effect could be easily added to the numerical scheme of the model, but it would increase significantly the difficulty of its mathematical study. Therefore we choose to ignore the threshold effect in the erosion process. This approximation is justified in the regime that we study, because the water height and water velocity are not close to zero. The last equation is the conservation equation for the sediments dissolved in the fluid. This concentration, denoted by c is the average of the concentration over the height of the fluid, thus is given in gram per square meter. The right hand side is the source term, which represents the exchange of sediments between the ground and the fluid, caused by the erosion and the sedimentation. The model is set in two dimensions, the variables are the time t ∈ R + and the position (x, y), which belongs to a domain Ω ⊂ R 2 .

In [START_REF] Chen | On the equations of landscape formation[END_REF], the authors reviewed various landscape evolution models and focused on system [START_REF] Anand | Linear layout of multiple flow-direction networks for landscape-evolution simulations[END_REF]. They listed some open mathematical problems like local existence in time, regularity of solutions or stability. They performed a numerical study on the erosion of a gaussian shaped hill to demonstrate the ability of the model to exhibit pattern formation and studied the influence of some parameters on the complexity of these patterns. In [START_REF] Lebrun | Numerical simulation of landscape evolution models[END_REF], the authors performed numerical simulations in a more practical context where the initial topography is a real one like the ones found in La Reunion or Madeira islands: they took a particular care of the visualization of landscape evolution through a particular colorization process. The simulations show the capacity of the model to describe a realistic evolution in time of the rivers and gullies on these landscapes, provided its parameters are properly chosen. Note that the detachment-limit hypothesis is valid for bedrock rivers and must be appropriately modified for alluvial rivers, where sediments form a layer of loose material. In this latter case, other classes of models may be considered like the shallow water equations coupled with Exner equations for the transport of sediments: see [START_REF] Fernández-Nieto | Formal deduction of the saint-venant-exner model including arbitrarily sloping sediment beds and associated energy[END_REF] for a formal derivation of these models from bi-layer type models (one layer for the fluid and one layer for the sediment) and [START_REF] Fernandez-Nieto | On the influence of the thickness of the sediment moving layer in the definition of the bedload transport formula in exner systems[END_REF][START_REF] Escalante | Modelling of bedload sediment transport for weak and strong regimes[END_REF] for various numerical strategies to perform simulations of bedload sediment transport. The aim of this paper is to explore the mechanism of pattern formation in the soil caused by water flow, and to precise the form and the frequency of apparition of these patterns. In landscapes, these patterns are channels, bed rivers, valley. We will focus on a simple framework where the initial bottom surface is a tilted plane. This surface is covered by a layer of water, the water flow from the top of the plane. The situation is represented in Figure 1, where the tilted plane is seen from the side. Although it is a less realistic framework than the landscape evolution described in [START_REF] Lebrun | Numerical simulation of landscape evolution models[END_REF] where real initial topographies are considered, it is sufficient to study the emergence of patterns on the surface. Note that experiments of erosion on a flat surfaces were conducted in [START_REF] Guérin | Streamwise dissolution patterns created by a flowing water film[END_REF] with a thin film flow which erodes blocks of plaster or salt.

The main outcomes of this paper are the following:

• We establish the well-posedness of System (2) when the initial fluid height does not vanish.

• We study the spectral stability of constant states when the bottom is an inclined plane and characterize rigorously the spectral instability through a non dimensional number, so called channelization index.

• We perform direct numerical simulations of System (2) for realistic experimental data and on time scales associated to the erosion phenomenon. We tackle the problem of severe CFL restrictions due to very short time scales associated to the fluid flow.

We provide thereafter a short description of our main results. First, we prove the well posedness character of System (2) locally in time under some assumptions on the data and the parameters. We prove the existence of solutions in suitable Sobolev spaces. In order to simplify notations, we will denote for any T > 0 and H a Hilbert space L 2 T (H ) := L 2 ((0, T ); H ); C T (H ) := C((0, T ); H ) the space of functions u : t → u(t, .) ∈ H that are respectively L 2 integrable or continuous in time on the interval (0, T ). We prove the following result: Theorem 1.1. Let m > 0, n > 3 or n = 2, K > 0, and T 0 > 0. Let us fix two constants fluid heights h re f > h min > 0. Suppose that the initial data h 0 , z 0 , c 0 satisfy

h 0 -h re f , z 0 ∈ H k+1 (R 2 ), c 0 ∈ H k , h 0 (x) ≥ 2 h min ∀x ∈ R 2 . with k = 3. Suppose that r ∈ L 2 T 0 (H k ), K h min -||h 0 || 2 L ∞ ≥ 0.
Then there exists 0 < T < T 0 such that System (2) admits a unique solution (h, z, c) with

h -h re f , z ∈ L 2 T (H k+2 ) ∩C T (H k+1 ), c ∈ C T (H k ).
The proof of Theorem 1.1 is based on energy estimates and on a fixed point argument. The problem of the well posedness of this system was raised in the conclusion of [START_REF] Chen | On the equations of landscape formation[END_REF]. We prove that it is well-posed locally in time, and under the assumption that the fluid height does not vanish. With this assumption, the equations on water height and water concentration are parabolic, and this property is used in the proof, for the energy estimates.This validates the model in the situation where the fluid height is close to a constant, hence, far from vanishing. This will be the case in our stability study, and our numerical simulations. The case of vanishing fluid height is more involved as the model becomes degenerated parabolic and the regularizing effect on the fluid height, and thus on the fluid velocity, is lost: one has to consider weaker classes of solutions. More precisely, dropping the z dependence, the equation on h reads

∂ t h = h∆h + |∇h| 2 .
In [START_REF] Brändle | Viscosity solutions for quasilinear degenerate parabolic equations of porous medium type[END_REF], the authors prove the well posedness for this equation within the class of bounded, continuous and non negative viscosity solutions. This is beyond the scope of this paper to couple this type of solutions with the other equations of the model, and we leave this problem open for future studies. An other interesting question is the numerical simulation of the model in the presence of dry areas: this will be carried out in a forthcoming work.

The second part of the paper is dedicated to pattern formation when the initial topography is an inclined plane and the bottom surface is weakly eroded. Our aim is to identify instability mechanisms that could explain the formation of patterns. For that purpose, we linearize System (2) around a constant state and we study the conditions of spectral instability and the nature of the instabilities. We expect that when the system is spectrally unstable for some parameters and wave vectors, a slight perturbation of the system with these unstable modes will grow up and lead to the formation of patterns in the soil. Stability studies have been done previously, for other landscape evolution models. The papers [START_REF] Smith | Stability and the conservation of mass in drainage basin evolution[END_REF] and [START_REF] Loewenherz | Stability and the initiation of channelized surface drainage: a reassessment of the short wavelength limit[END_REF] analyse a model of two equations, where the water is supposed to be at equilibrium. This model has steady solutions for which the soil height can be concave in some areas and convex in other areas, depending on the sediment discharge law. They show that the linearised system is stable in the convex parts and unstable in the concave parts, with a stronger instability in the transverse direction. In their model they use the sediment transport law, thus its not the same framework as in the model [START_REF] Anand | Linear layout of multiple flow-direction networks for landscape-evolution simulations[END_REF]. Note that in a couple of recent papers [START_REF] Anand | Linear layout of multiple flow-direction networks for landscape-evolution simulations[END_REF] and [START_REF] Bonetti | Channelization cascade in landscape evolution[END_REF], the following system of 2 PDEs was considered:

∂ t z = K∆z -e( h H ) m |∇z| n +U, ∂ t h = div(hv 0 ∇z |∇z| ) + R, (3)
where R,U and V 0 are constants. A numerical scheme is designed in [START_REF] Anand | Linear layout of multiple flow-direction networks for landscape-evolution simulations[END_REF] for [START_REF] Andrew | Derivatives of eigenvalues and eigenvectors of matrix functions[END_REF] where the time derivative of the fluid height is neglected with test cases where the initial bottom topography is pyramidal. It is found that a channelization index C I = e m+n K n U 1-n drives the formation of channels: the number of channels and their branching increase with C I . This analysis is completed by a spectral stability analysis of a spatially non homogeneous steady state where the topography is a hillslope which is divided in the middle. It is found numerically that there exists a critical C 0 I such that the steady state is stable if C I ≤ C 0 I and unstable otherwise.

The spectral study carried out in this paper is new and our analysis provides some explanations for the formation of patterns in landscapes. The appearance of channels on the flat plane is indeed the initial stage of development for the formation of valley and rivers in landscapes. The stability of the system depends on the parameters, in particular the constant of creep K plays an important role in this study. We show that there is a critical value K such that if K ≥ K, and if another condition on parameters is satisfied, then the system is spectrally stable at all frequencies. If K < K then there exist some wave numbers and vectors for which the system is spectrally unstable. Moreover, the instabilities grow as the wave vectors of the perturbations points in the direction transverse to the flow, which explains the formation of gullies and channels aligned with the direction of the fluid flow. We then recover qualitatively the results of [START_REF] Bonetti | Channelization cascade in landscape evolution[END_REF].

Finally, our stability study is completed by direct numerical simulations, which illustrate the appearance of patterns for the nonlinear system. The space and time scales of the model can take a large range of values, depending of the environment. On real landscapes, the domain size can be measured in kilometers, with a very slow erosion rate, in the order of 20 -200 millimeters per thousand year. In the experience on salt and plaster made in [START_REF] Guérin | Streamwise dissolution patterns created by a flowing water film[END_REF], the domain has a size of the order of ten centimeters whereas the erosion speed is around one millimeter per hour (so, much faster than in real landscape) and the fluid velocity is 1 meter per second. Parameters chosen in the numerical simulations are based on these experiments. As in [START_REF] Anand | Linear layout of multiple flow-direction networks for landscape-evolution simulations[END_REF][START_REF] Bonetti | Channelization cascade in landscape evolution[END_REF], we have observed that decreasing K (respectively increasing the channelization index C I ) reinforces the channelization process. Note that we have focused here on the formation of channels: unlike simulations made in [START_REF] Lebrun | Numerical simulation of landscape evolution models[END_REF] where the initial state is a matured landscape, we start from a simple state and the landscape evolves by himself in the simulations.

The paper is organised as follow: First, in Section 2, we describe the model and the associated system of equations. Then, Section 3 is devoted to the proof of the well posedness character of the system in short time, (see Theorem 1.1). Next, in Section 4 we carry out a spectral stability analysis of the System (2) linearized about a stationary solution. These stability results are compared to direct numerical simulations of the nonlinear system (2) in Section 4.2. Finally, Section 6 draws a brief conclusion of the paper and provides some future perspectives.

The landscape evolution model

In this section, we introduce the landscape evolution model considered in this paper. This is a system of three partial differential equations for the fluid height h, the bottom topography z and the sediment concentration c.

Evolution of topography. The evolution of the bottom topography z is given by:

∂ t z = K∆z -E + S.
The functions E = E(t, x, y) and S = S(t, x, y) represent the erosion speed of the soil and the sedimentation speed respectively, with t ≥ 0, (x, y) ∈ R 2 .

The parameter K > 0 is a constant, and the term K∆z models the creep of the soil. This phenomenon is a slow diffusive movement of the soil which occurs at large time and space scales. This movement is caused by several processes, such as the gravitational flow of the soil, wind, rain splash, expansions and contractions of the soil due to freeze-thaw, wet-dry and hot-cold cycles, or biological activity. In sufficiently eroded landscape there are generally not many sharp edges, and this creep term, which tends to smooth the bottom surface models this phenomenon. We shall see that this term plays a significant role in the well-posedness of the model. However, the creep effect is not supposed to be relevant in the formation of patterns as it is a short time effect and should be supposed to be small in comparison to the erosion and sedimentation terms.

The erosion of the surface is caused by the shear stress and the friction of the water flow. Assuming that the fluid velocity is constant across the fluid layer, this amounts to consider that the erosion increases with the (norm of the) water velocity, and with the water discharge Q, as in [START_REF] Howard | Channel changes in badlands[END_REF]:

E(t, x, y) = αQ µ v ν .
Consequently, as Q = hv, we will suppose that the erosion speed depends on a power of the norm of the water velocity, and of a power of the water height. Thus we set

E(t, x, y) = e h(t, x, y) H m |v(t, x, y)| V n ,
where e is the erosion speed in the conditions h = H and |v| = V with H,V > 0 that respectively represent a reference fluid height and fluid velocity.

The sedimentation occurs when the concentration of sediments in water is high enough. The sedimentation speed increases with the concentration of sediment in the fluid. For the sake of simplicity, we suppose that this speed is proportional to the concentration and we set:

S(t, x, y) = s c(t, x, y) c sat ,
with s the speed of sedimentation in the reference condition c = c sat . Therefore, the soil elevation evolves according to the equation :

∂ t z = K∆z -e h(t, x, y) H m |v(t, x, y)| V n + s c(t, x, y) c sat . ( 4
)
The landscape evolution model. We complete Equation (4) with two evolution equations for the fluid height h and sediment concentration c. The mass conservation law for the fluid reads

∂ t h + div(hv) = r, (5)
where r is a source term, modeling an incoming flow in a channel or the rain over the bottom.

On the other hand, the mass conservation law for the sediment reads

∂ t (hc) + div(chv) = ρ s (E -S), (6)
where ρ s is the volumetric mass density of the sediments, and is constant. In order to close System (4), ( 5), [START_REF] Chen | On the equations of landscape formation[END_REF], we need to write an equation for the fluid velocity. One possibility would be to write a shallow water type model with an evolution equation for the momentum hv. We rather choose the simpler closure

v = -µ∇(h + z), (7)
where µ > 0 is some characteristic fluid velocity and ∇(h + z) is the gradient of the fluid surface elevation. System (4), ( 5), ( 6), ( 7) is closed and we shall consider its well-posedness in Section 3.

We are also interested in the pattern formation at the surface of the soil. For that purpose, we have chosen to explore the case of water flowing down an inclined plane. This situation was considered experimentally in [START_REF] Guérin | Streamwise dissolution patterns created by a flowing water film[END_REF]. The domain Ω ⊂ R 2 has length L x and width L y :

Ω = [0, L x ] × [0, L y ].
Denote θ the inclination of the plane. We can decompose the bottom topography z as z(t, x, y) = (L xx) tan θ + z(t, x, y) where z is the eroded height of the soil. Thus the fluid velocity is written as v(t, x, y) = µ(tan θ , 0) -µ∇(z + h). Consequently, omitting the tilde over z, System (4), ( 5), ( 6), [START_REF] Chen | Landscape evolution models: A review of their fundamental equations[END_REF] admits the new form:

               ∂ t h + µ tan θ ∂ x h = µ div(h∇(h + z)), h∂ t c + µh tan θ ∂ x c = µh∇(h + z).∇c + ρ s e h(t, x) H m |v(t, x)| V n -ρ s s c(t, x) c sat , ∂ t z = K∆z -e h(t, x) H m |v(t, x)| V n + s c(t, x) c sat . (8a) (8b) (8c)
3 Well-posedness of the landscape evolution model

In this section we study the existence and uniqueness of solutions of the system (4), ( 5), ( 6), [START_REF] Chen | Landscape evolution models: A review of their fundamental equations[END_REF], locally in time.

Hypothesis on the system of equations

We consider System (4), ( 5), ( 6), [START_REF] Chen | Landscape evolution models: A review of their fundamental equations[END_REF] where we set, for simplicity, µ = 1, ρ s = 1, c sat = 1 and

H = V = 1.
The choice of these constants does not change anything in the proof Theorem 1.1, we fix them to simplify the notations. Provided that ∀(t, x) ∈ R + × R 2 , h(t, x) ≥ h min > 0, the equations on h and on z are parabolic equations. As long as h does not vanish, the equation on c can be written as

∂ t c + v.∇c = eh m-1 |v| m -sc/h -rc/h.
Therefore, in order to prove Theorem 1.1, we consider the equivalent following system, composed of two parabolic equations and one transport equation:

       ∂ t h -div(h∇h) -div(h∇z) = r(t, x), ∂ t z = K∆z + sc -eh m |v| n , ∂ t c + v.∇c = eh m-1 |v| n -sc/h -rc/h, (9a) (9b) (9c)
with the initial conditions h(0, x) = h 0 (x), z(0, x) = z 0 (x) and c(0, x) = c 0 (x), and where t ∈ R + , x ∈ R 2 .

In Subsection 3.2, we recall some results concerning Sobolev spaces that will be used to prove Theorem 1.1.

The well posedness of System ( 9) is proved in Subsections 3.3 -3.6 by using a fixed point argument and energy inequalities. We first establish an a priori estimate on the solutions of (9) in Subsection 3.3. Then, we build a sequence of approximate solutions in Subsection 3.4 and provide uniform estimates on these solutions. We show that it forms a Cauchy sequence and converges to a solution of System [START_REF] Culling | Soil creep and the development of hillside slopes[END_REF]. Finally, we prove the uniqueness of solutions in Subsection 3.6. In these sections, we denote by ∇ p f the vector made by the partial derivatives of order p of the function f . We also denote

(∇ f ) 2 = (∂ x i f ∂ x j f , 1 ≤ i ≤ 2, 1 ≤ j ≤ 2).

Sobolev injections

In this section we recall some properties of Sobolev spaces, that will be used in the next sections. The proofs of these properties can be found in [START_REF] Adams | Sobolev spaces[END_REF], and in [START_REF] Evans | Partial differential equations[END_REF] for Proposition 3.5. The first proposition concerns the imbedding of the Hilbert space H k (R 2 ) into Sobolev spaces with smaller derivation index. The symbol → indicates that the injection is continuous.

Proposition 3.1. If 2 < q < +∞ and k ∈ N then H k+1 (R 2 ) → W k,q (R 2 ).
The following proposition gives the imbedding of the Hilbert space H k (R 2 ) into a space of smooth functions.

Proposition 3.2. If k ≥ 2 then H k (R 2 ) → C k-2,α (R 2 ), ∀ 0 < α < 1. In particular, H 2 (R 2 ) ⊂ C 0 b (R 2 ). Proposition 3.3. If k > 1 then H k (R 2 ) → W k-2,∞ (R 2 ).
Finally, Proposition 3.4 provides a bound on the norm of a power of functions, that will be useful to control non linear terms in the equations.

Proposition 3.4. H k (R 2 ) is an algebra for k > 1. Thus if q ∈ N * then || f q || H k (R 2 ) ≤ C|| f || q H k (R 2 ) . If q ∈ R \ N and q > k, one has || f q || H k (R 2 ) ≤ C|| f || q H k (R 2 )
. Finally, Proposition 3.5 concerns spaces involving time, and gives the continuity in time of a function provided this function and its derivative in time have enough regularity.

Proposition 3.5. If f ∈ L 2 T (H k+2 (R 2 )) and ∂ t f ∈ L 2 T (H k (R 2 )) then f ∈ C 0 T (H k+1 (R 2 )).
In what follows, we will denote H k (R 2 ) = H k for the sake of simplicity.

An a priori estimate

In this section we give an a priori estimate on solutions of System ( 9), which will be used in the Section 3.5. We fix k = 3, and assume m > 0, n > 3 or n = 2 as in the hypothesis of Theorem 1.1. [START_REF] Culling | Soil creep and the development of hillside slopes[END_REF]. Assume that the hypothesis of Theorem 1.1 are satisfied. Then there exists T 1 ≤ T such that hh re f and z are bounded in

Proposition 3.6. Let (h -h re f , z, c) ∈ L 2 T (H k+2 ) ∩C T (H k+1 ) 2 × C T (H k ) be a solution of System (
C T 1 (H k+1 ) ∩ L 2 T 1 (H k+2 ), c is bounded in C T 1 (H k )
, and they satisfy the estimate:

E (t) ≤ e Ct E (0) +C t 0 r 2 H k , t 0 ∇h 2 H k+1 + ∇z 2 H k+1 ≤ C e Ct E (0) +C t 0 r 2 H k ,
where C is a constant depending on h 0 , z 0 , c 0 , h min and E is defined as:

E (t) = 1 2 h -h re f 2 + ∇h 2 H k + z 2 H k+1 + c 2 H k .
Proof. We first provide Sobolev estimates on the fluid height h. We multiply the equation (9a) by hh re f and integrate it over R 2 . One obtains:

1 2

d dt h -h re f 2 L 2 = R 2 (h -h re f ) div (h∇(h + z)) + R 2 (h -h re f )r. ( 10
)
By integrating by part [START_REF] Culling | Analytical theory of erosion[END_REF], and under the assumption that h ≥ h min , we obtain for any t ≥ 0 fixed:

1 2 d dt h -h re f 2 L 2 + h min ∇h 2 L 2 ≤ r L 2 h -h re f L 2 + h L ∞ ∇z L 2 ∇h L 2 . ( 11
)
Then for all p ∈ {1, . . . , k + 1}, we differentiate p times Equation (9a), multiply it by ∇ p h and integrate it with respect to the space variable:

1 2 d dt ∇ p h 2 L 2 + R 2 ∇ p-1 [div(h∇(h + z))] ∇ p+1 h = - R 2 ∇ p-1 r ∇ p+1 h (12)
We estimate the second term in [START_REF] Escalante | Modelling of bedload sediment transport for weak and strong regimes[END_REF]. As H 1 → L 4 and H 2 → L ∞ , we find:

R 2 ∇ p-1 [div(h∇h)] ∇ p+1 h = R 2 h|∇ p+1 h| 2 + R 2 p-1 ∑ i=1 p i ∇ i h∇ p-i+1 h∇ p+1 h + R 2 ∇h∇ p h∇ p+1 h ≥ h min ||∇ p+1 h|| 2 L 2 -C ∇h 2 H p ∇h H k+1 .
For the third term in [START_REF] Escalante | Modelling of bedload sediment transport for weak and strong regimes[END_REF], we proceed similarly:

R 2 ∇ p-1 div(h∇z)∇ p+1 h = R 2 h∇ p+1 z∇ p+1 h + R 2 p ∑ i=1 p i ∇ i h∇ p-i+1 z∇ p+1 h ≤ ||h|| L ∞ ||∇z|| H k+1 ||∇ p+1 h|| L 2 +C||∇z|| H p ||∇h|| H p ||∇h|| H k+1 .
By inserting these two estimates into [START_REF] Escalante | Modelling of bedload sediment transport for weak and strong regimes[END_REF], one obtains:

1 2 d dt ∇h 2 H k + h min ∇h 2 H k+1 ≤ r H k ∇h H k+1 + h L ∞ ∇z H k+1 ∇h H k+1 +C ( ∇z H k + ∇h H k ) 2 ∇h H k+1 . (13)
Next, we derive an estimate on the bottom topography z. We multiply equation (9b) by z and integrate over space R 2 . We get:

1 2 d dt ||z|| 2 L 2 + K||∇z|| 2 L 2 = s R 2 c z -e R 2 h m |v| n z ≤ s||c|| L 2 ||z|| L 2 + e||h|| m L ∞ ||v|| n L 2n ||z|| L 2 ≤ s||c|| L 2 ||z|| L 2 +Ce h m L ∞ (||∇h|| H 1 + ||∇z|| H 1 ) n ||z|| L 2 . ( 14
)
The last inequality is a consequence of the injection H 1 → L q , with q = 2n ≥ 2. Let us now estimate the derivatives of z of order p ∈ {1, . . . , k + 1}. We differentiate Equation (9b) p times and multiply by ∇ p z.

1 2 d dt ∇ p z 2 L 2 + K ∇ p+1 z 2 L 2 = s R 2 ∇ p c∇ p z -e R 2 ∇ p (h m |v| n ) ∇ p z = -s R 2 ∇ p-1 c∇ p+1 z + e R 2 ∇ p-1 (h m |v| n ) ∇ p+1 z ≤ s||∇ p-1 c|| L 2 ||∇ p+1 z|| L 2 + e||∇ p-1 (h m |v| n )|| L 2 ||∇ p+1 z|| L 2 . ( 15
)
By adding the estimates ( 14) and ( 15) for p = 1, . . . , k + 1, we obtain:

1 2 d dt ∇z 2 H k + K ∇z 2 H k+1 ≤ s c H k ∇z H k+1 + e h m |v| n H k ∇z H k+1 . (16)
There remains to estimate the erosion term h m |v| n in H k norm. We assumed that n > 3 or n = 2, which implies that v → |v| n is C 3 (R 2 ; R). We estimate successively ∇ p (h m |v| n ) for p = 1, 2, 3 = k. By using successively the injections L q → H 1 for q ≥ 2 and H 2 → L ∞ , one finds:

∇(h m |v| n ) L 2 ≤ h m-1 L ∞ ( ∇h H 2 + ∇z H 2 ) n ( h L ∞ + ∇h L 2 ), ∇ 2 (h m |v| n ) L 2 ≤ h m-2 L ∞ ( ∇h H 2 + ∇z H 2 ) n ( h L ∞ + ∇h H 1 ) 2 , ∇ 3 (h m |v| n ) L 2 ≤ h m-3 L ∞ ( ∇h H 3 + ∇z H 3 ) n ( h L ∞ + ∇h H 1 ) 2 ( h L ∞ + ∇h H 3 + ∇z H 3 ) .
We deduce the following estimate on the bottom topography:

1 2 d dt ∇z 2 H k + K ∇z 2 H k+1 ≤ s c H k ∇z H k+1 +C h m-3 L ∞ ( h L ∞ + ∇h H 3 + ∇z H 3 ) n+3 ∇z H k+1 . (17)
Finally, we derive a priori estimate for the sediment concentration c: the basic energy estimates reads:

1 2 d dt c 2 L 2 ≤ ( ∇h H 3 + ∇z H 3 ) c 2 L 2 + e h m-1 L ∞ ( ∇h H 1 + ∇z H 1 ) n c L 2 + (s + r L ∞ ) c 2 L 2 h min . (18)
The first term on the right hand side of ( 18) is related to the advection, the second term to the erosion and the third one to the sedimentation and the source term. Now, for all p ∈ {1, . . . , k}, we differentiate p times Equation (9c) and multiply it by ∇ p c. There is no additional issue with respect to the former computations on z and one finds:

1 2 d dt ||c|| 2 H k ≤ C ( ∇h H k+1 + ∇z H k+1 ) c 2 H k + (s + r H k ) c 2 H k h min P ∇h H k+1 h min +eC h m-4 L ∞ ( h L ∞ + ∇h H 3 + ∇z H 3 ) n+3 . ( 19
)
Here P is some polynomial of degree three with positive coefficients. Let us denote

Q h (a, b) = Ka 2 + h min b 2 -h L ∞ ab.
By combining the estimates ( 11), ( 13), ( 14), ( 17) and ( 19) and the Young inequality on products:

∀(a, b) ∈ R 2 , ∀ε > 0, ab ≤ ε 2 a 2 + 1 2ε b 2 ,
one finds that for any ε > 0, there exists C(ε) such that

d dt E (t) + Q h ( ∇h H k+1 , ∇z H k+1 ) ≤ ε ∇h 2 H k+1 + ∇z 2 H k+1 +C(ε) r 2 H k + F (h, z, c)(t) , (20)
where

F (h, z, c) = h -h re f 2 + h L ∞ ∇z L 2 ∇h L 2 + c L 2 z L 2 + c 2 H k + c 4 H k 1 + (h min ) -2 P 2 ∇h H k+1 h min + h m L ∞ ( ∇h H 1 + ∇z H 1 ) n z L 2 + ∇h H k + ∇z k H 4 + h m-3 2 L ∞ ( h L ∞ + ∇h H k + ∇z H k ) 2(n+1) + h m-4 2 L ∞ ( h L ∞ + ∇h H k + ∇z H k ) (n+1) c H k . (21)
From ( 21), one easily proves that there exists a constant C depending only on h min and α > 0 such that

F (h, z, c)(t) ≤ C (1 +CE (t) α ) E (t)
Note that we have used the estimate:

h L ∞ ≤ h re f + h -h re f L 2 +C ∇h H 1 , h m-i L ∞ ≤ h m L ∞ (h min ) -i , i = 3, 4. Now, the quadratic form Q h is positive semi-definite provided that Kh min > h 2 L ∞ /4. Under the assumption that h 0 L ∞ < √ Kk min , there exists a time T * such that sup (0,T * ) h L ∞ ≤ √ 3Kk min . Thus there exists ε > 0 such that Q h (a, b) ≥ 2ε(a 2 + b 2 ).
We deduce from [START_REF] Lebrun | Numerical simulation of landscape evolution models[END_REF] that for all t ∈ (0, T * ), one has:

d dt E (t) + ε ∇h 2 H k+1 + ∇z 2 H k+1 ≤ C(ε) r 2 H k +C (1 +CE (t) α ) E (t), (22)
Now, we integrate Equation ( 22) with respect to time: one finds

E (t) + ε t 0 ∇h 2 H k+1 + ∇z 2 H k+1 ≤ E (0) +C t 0 r 2 H k +C t 0 (1 +CE (s) α )E (s)ds. (23)
By apply one more time a continuity argument, there exists

T 1 ≤ T * , such that E (t) ≤ 2E (0) for all t ∈ [0, T 1 ].
This implies that

E (t) + ε t 0 ∇h 2 H k+1 + ∇z 2 H k+1 ≤ E (0) +C t 0 r 2 H k +C(1 + 2 α CE (0) α ) t 0 E (s)ds. (24)
By applying Gronwall lemma, one finds

E (t) ≤ e Ct E (0) +C t 0 r 2 H k , t 0 ∇h 2 H k+1 + ∇z 2 H k+1 ≤ C e Ct E (0) +C t 0 r 2 H k .
This concludes the proof of the energy estimate.

The approximate system

In this section, we prove the existence of a sequence of solutions (h i , z i , c i ) i∈N to the following linear system:

           ∂ t h i+1 -div(h i ∇h i+1 ) -div(h i+1 ∇z i ) = r(t, x), ∂ t z i+1 = K∆z i+1 + sc i -eh m i |v i | n , ∂ t c i+1 + v i .∇c i+1 = eh m-1 i |v i | n -sc i /h i -rc i /h i , h i (0, x) = h(0, x), z i (0, x) = z(0, x), c i (0, x) = c(0, x). (25a) (25b) (25c)
with initial conditions h 0 (t, x) = h 0 (x), z 0 (t, x) = z 0 (x) and c 0 (t, x) = c 0 (x). Then we show that these solutions are uniformly bounded with respect to i ∈ N for a suitable Sobolev norm. In what follows, we will denote

h min = in f x∈R 2 (h 0 (x))/2.
The classical theory of parabolic equations and of linear transport equations provides conditions to obtain a well-posed system of equations. For the two parabolic equations (25a) and (25b), we state the following result (see [START_REF] Cherrier | Linear and quasi-linear evolution equations in Hilbert spaces[END_REF] for more details):

Proposition 3.7. Assume that h 0 -h re f , z 0 ∈ H k+1 (R 2 ), r, sc i -eh m i |v i | n ∈ L 2 T (H k ), h i ∈ L ∞ ([0, T ] × R 2 ) with h i ≥ h min and ∇h i , ∇z i ∈ L 1 T (H k ).
Then there exists a unique solution (h i+1 , z i+1 ) to the equations (25a), (25b), and h i+1 -

h re f , z i+1 ∈ L 2 T (H k+2 ) ∩C T (H k+1 (R 2 )).
For the transport equation (25c), we use theorem 7.2.2 in [START_REF] Métivier | Para-differential calculus and applications to the Cauchy problem for nonlinear systems[END_REF] to state the result:

Proposition 3.8. Assume that c 0 ∈ H k (R 2 ), ∇v i ∈ L ∞ T (H k-1 ), eh m-1 i |v i | n -sc i /h i -rc i /h i ∈ L 1 T (H k ) and v i ∈ L 1 T (W 1,∞ ).
Then there exists a unique solution c i ∈ C T (H k ) to the equation (25c).

We prove by induction that for all i ∈ N, System (25) is well-posed. First, for i = 0, the functions of (h i , z i , z i ) are time independent and h 0 ≥ h min . Thus, as T 0 = T is finite, we have sc 0eh m

0 |v 0 | n ∈ L 2 T 0 (H k ), h 0 ∈ L ∞ T 0 (R 2 ) and ∇h 0 , ∇z 0 ∈ L 1 T 0 (H k ).
Then there exists a unique solution

(h 1 , z 1 ) ∈ L 2 T 1 (H k+2 ) ∩ C T 1 (H k+1 ) of (25a) and (25b) with 0 < T 1 ≤ T 0 such that h 1 ≥ h min . Similarly, there exists c 1 ∈ C T 1 (H k ) solution of (25c). Now, if we assume that h i -h re f , z i ∈ L 2 T (H k+2 ) ∩C T i (H k+1 ) and c i ∈ C T i (H k ) then, one has h i ∈ h re f +C T i (H 2 ) ⊂ L ∞ ([0, T ] × R 2 ). Moreover, we have ∇h i , ∇z i ∈ C T i (H k ) ⊂ L 1 T i (H k ) and c i ∈ C T i (H k ) ⊂ L 2 T i (H k ). Finally, one has h m i |v i | n H k ≤ C h m-3 i L ∞ ( h i L ∞ + ∇h i H k + ∇z i H k ) n+3 ∈ L ∞ T i (H k ) ⊂ L 2 T i (H k ).
Thus, the assumptions of Proposition 3.7 are satisfied and there exists a unique solution h i+1 , z i+1 of (25a,25b) such that

h i+1 -h re f , z i+1 ∈ L 2 T i+1 (H k+2 ) ∩ C T i+1 (H k+1 (R 2 )) with 0 < T i+1 ≤ T i such that h i+1 ≥ h min . The existence of a solution c i+1 ∈ C T i+1 (H k ) follows similarly. Proposition 3.9 (Uniform bounds). Denote E i (t) = h i -h re f L 2 + ∇h i H k + z i H k+1 + c i H k .
There exists B > 0 and 0 < T * ≤ T independent of i ∈ N such that, for all i ∈ N

sup t∈[0,T * ] E i (t) ≤ B,
and

h i (t, x) ≥ h min , ∀x ∈ R 2 .
Moreover,

T * 0 ∇h i 2 H k+1 + ∇z i 2 H k+1 ≤ C(T * )B, ∀i ∈ N * .
Proof. We proceed by induction. One has easily E 0 (t) = E (0) ≤ B and h 0 ≥ h min since h 0 , z 0 , c 0 are time independent. Following the strategy used to derive a priori estimates, one can prove that for some ε > 0, for all t ∈ [0, T ],

d dt E 1 (t) + ε ∇h 1 2 H k+1 + ∇z 1 2 H k+1 ≤ C(ε)( r 2 H k + E 0 (t) α )E 1 (t). ( 26
)
for some constant C(ε) depending only on ε and h 0 . By integrating Equation ( 26) with respect to time, one finds:

E 1 (t) + ε t 0 ∇h 1 2 H k+1 + ∇z 1 2 H k+1 ≤ E (0) +C(ε) t 0 r 2 H k + t 0 C(ε)E 0 (s) α E 1 (s)ds. Denote B = 2 E (0) +C(ε) T 0 r 2 H k . We choose T * such that e T * C(ε)B α ≤ 2.
Then, by applying Gronwall's lemma, one obtains:

E 1 (t) ≤ B, ∀t ∈ [0, T * ], T * 0 ∇h 1 2 H k+1 + ∇z 1 2 H k+1 ≤ C(ε) 1 + ln(2) ε B.
Next, we have

h 1 = h 0 + t 0 div(h 0 ∇(h 1 ) + h 1 ∇z 0 )) + r. Thus h 1 -h 0 L ∞ ≤ C( Bt + √ Bt)
for some constant C independent of the problem and related to Sobolev injections whereas B = B + h re f √ B. Then, there exists 0 < T * ≤ T * such that C( Bt + √ Bt) ≤ h min and we deduce that

h 1 ≥ h 0 -h min ≥ h min ∀t ∈ [0, T * ].
This proves the initial step for i = 1. Now assume that E i (t) ≤ B and h i ≥ h min for all t ∈ [0, T * ]. The estimates on h i+1 , z i+1 and c i+1 are a direct consequence of the energy estimate:

d dt E i+1 (t) + ε ∇h i+1 2 H k+1 + ∇z i+1 2
H k+1 ≤ C(ε)( r 2 H k + E i (t) α )E i+1 (t) (27)
which is proved by following the strategy used to derive the a priori estimates. This completes the proof of the proposition.

Proposition 3.10. The sequences

(∂ t h i ) et (∂ t z i ) are uniformly bounded in L 2 T * (H k ). The sequence (∂ t c i ) is uniformly bounded on L 2 T * (H k-1 ).
Proof. For p ∈ {0, . . . , k} we differentiate p times the equations (25a), (25b) (25b) and multiply it by ∇ p ∂ t h i , ∇ p ∂ t z i and ∇ p ∂ t c i respectively, and integrate in space. Using Proposition 3.9, we obtain bounds for

∂ t h i L 2 T * (H k ) , ∂ t z i L 2 T * (H k ) and ∂ t c i L 2 T * (H k-1 ) .

Convergence of the sequences

In what follows, we denote T * = T in order to simplify the notations.

Proposition 3.11. The sequences

(h i ) i∈N , (∂ t h i ) i∈N and (z i ) i∈N , (∂ t z i ) i∈N are Cauchy sequences in L 2 T (L 2 ). The sequences (c i ) i∈N and (∂ t c i ) i∈N are Cauchy sequences in C T (L 2 ).
Proof. In the following we denote the quantities of the form f i+1f i by δ f i . For all i ≥ 1, the equations for δ h i , δ z i , δ c i are written as:

                 ∂ t δ h i -div(h i ∇δ h i + δ h i-1 ∇h i ) -div(h i+1 ∇δ z i + δ h i ∇z i-1 ) = 0, ∂ t δ z i -Kδ ∆z i = sδ c i -e h m i (|v i | n -|v i-1 | n ) + |v i-1 | n (h m i -h m i-1 ) , ∂ t δ c i + δ v i-1 ∇c i+1 + v i-1 δ ∇c i = e h m-1 i (|v i | n -|v i-1 | n ) + |v i-1 | n (h m-1 i -h m-1 i-1 ) -(s + r) c i h i - c i-1 h i-1 . (28a) (28b) (28c)
Bounds on δ h i , δ z i and δ c i : We multiply the equation (28a) by δ h i , and integrate over space:

1 2 d dt δ h i 2 L 2 + R 2 (h i ∇h i+1 -h i-1 ∇h i )∇δ h i + R 2 (∇z i h i+1 -∇z i-1 h i )∇δ h i = 0.
First, on the one hand, we have

R 2 (h i ∇h i+1 -h i-1 ∇h i )∇δ h i = R 2 h i (∇δ h i ) 2 + δ h i-1 ∇h i ∇δ h i ≥ h min ∇δ h i 2 L 2 + R 2 δ h i-1 ∇h i ∇δ h i .
On the other hand, we have the estimate:

R 2 (∇z i h i+1 -∇z i-1 h i )∇δ h i = 1 2 R 2 ∇z i ∇(δ h i ) 2 + R 2 δ z i-1 h i ∇δ h i = - 1 2 R 2 ∆z i (δ h i ) 2 + R 2 δ z i-1 h i ∇δ h i ≤ C δ h i 2 L 2 + h i L ∞ δ z i-1 L 2 ∇δ h i L 2 ≤ C δ h i 2 L 2 + 1 2h min δ z i-1 2 L 2 + h min 2 ∇δ h i 2 L 2 .
Consequently, we obtain:

d dt δ h i 2 L 2 + h min ∇δ h i 2 L 2 ≤ C δ h i 2 L 2 +C δ h i-1 2 L 2 +C δ z i-1 2 L 2 . ( 29
)
We proceed similarly for δ z i . By using Equation (28b), one finds:

d dt δ z i 2 L 2 + 2K ∇δ z i 2 ≤ s( δ c i 2 L 2 + δ z i 2 L 2 ) + 2e R 2 (h m i |v i | n -h m i-1 |v i-1 | n )δ z i .
In order to bound the right-hand term, we use the inequality

||x| n -|y| n | ≤ n|x -y| max(|x|, |y|) n-1 : h m i |v i | n -h m i-1 |v i-1 | n = h m i (|v i | n -|v i-1 | n ) + |v i-1 | n (h m i -h m i-1 ) ≤ nh m i |v i -v i-1 | max( v i-1 L ∞ , v i L ∞ ) n-1 + m h min |v n i-1 |h i -h i-1 | max( h i L ∞ , h i-1 L ∞ ) m .
Consequently, as

h i , h i-1 and v i , v i-1 are uniformly bounded in L ∞ ([0, T ] × R 2 ): R 2 (h m i |v i | n -h m i-1 |v i-1 | n )δ z i ≤ C( v i -v i-1 L 2 + h i -h i-1 L 2 ) δ z i L 2 ≤ C( ∇δ h i-1 L 2 + ∇δ z i-1 L 2 + δ h i-1 L 2 ) δ z i L 2 ≤ C( ∇δ h i-1 2 L 2 + ∇δ z i-1 2 L 2 ) + δ h i-1 2 L 2 + δ z i 2 L 2 ). Thus, d dt δ z i 2 L 2 + 2K ∇δ z i 2 L 2 ≤ C δ z i 2 L 2 +C( ∇δ h i-1 2 L 2 + ∇δ z i-1 2 L 2 + δ h i-1 2 L 2 ). (30)
Finally, for δ c i , we use the same method as before with the equation (28c) and we obtain:

1 2 d dt δ c i 2 L 2 = -(δ v i-1 ∇c i+1 + v i ∇δ c i ) δ c i + e (h m i |v i | n -h m i-1 |v i-1 | n )δ c i -(s + r) c i h i - c i-1 h i-1 δ c i = - δ v i-1 ∇c i+1 δ c i - 1 2 div(v i )(δ c i ) 2 + e (h m i |v i | n -h m i-1 |v i-1 | n )δ c i -(s + r) c i h i - c i-1 h i-1 δ c i ≤ δ v i-1 L 2 ∇c i+1 L ∞ δ c i L 2 + ∇v i L ∞ δ c i 2 L 2 +C( ∇δ h i-1 2 L 2 + ∇δ z i-1 2 L 2 + δ h i-1 2 L 2 ) + s + r H 2 h min ( c i L 2 + c i-1 L 2 ) δ c i L 2 . ( 31
)
Then we add the inequalities (29), (30) and (31) and integrate on [0,t] with 0 ≤ t ≤ T :

δ h i 2 L 2 + δ z i 2 L 2 + δ c i 2 L 2 + h min t 0 ∇δ h i 2 L 2 dt + 2K t 0 ∇δ z i 2 L 2 dt ≤ δ h 0 2 L 2 + δ z 0 2 L 2 + δ c 0 2 L 2 +C t 0 δ h i 2 L 2 + δ z i 2 L 2 + δ c i 2 L 2 dt +C t 0 δ h i-1 2
L 2 + δ z i-1 2
L 2 + δ c i-1 2 L 2 dt +C t 0 ( ∇δ h i-1 2
L 2 + ∇δ z i-1 2 L 2 )dt.
We apply the Gronwall lemma and obtain, for all t ∈ [0, T ],

δ h i 2 L 2 + δ z i 2 L 2 + δ c i 2 L 2 + h min t 0 ∇δ h i 2 L 2 dt + 2K t 0 ∇δ z i 2 L 2 dt ≤ C δ h 0 2 L 2 + δ z 0 2 L 2 + δ c 0 2 L 2 + t 0 δ h i-1 2
L 2 + δ z i-1 2
L 2 + δ c i-1 2
L 2 dt + t 0 ( ∇δ h i-1 2
L 2 + ∇δ z i-1 2 L 2 )dt e Ct .
With this inequality, we deduce by induction on i that ∀i ∈ N,

δ h i 2 L 2 T (L 2 ) + δ z i 2 L 2 T (L 2 ) + δ c i 2 L ∞ T (L 2 ) ≤ C Te CT i i! δ h 0 2 H 1 + δ z 0 2 H 1 + δ c 0 2 L 2 .
Consequently the series ∑ δ h i

2 L 2 T (L 2 ) , ∑ δ z i 2 L 2 T (L 2 ) and ∑ δ c i 2 L ∞ T (L 2
) converge, thus the sequences (h i ) i∈N , (z i ) i∈N , (c i ) i∈N are Cauchy sequences in the required spaces.

Bounds on ∂ t δ h i , ∂ t δ z i and ∂ t δ c i : Like the estimates in Proposition 3.10, we use the system (28) and Proposition 3.9 to obtain the bounds.

By Proposition 3.11 there exists

h -h re f , z ∈ L 2 T (L 2 ) such that h i -h re f converges to h -h re f and z i converges to z in L 2
T (L 2 ). As (h ih re f ) and (z i ) are uniformly bounded in L 2 T (H k+2 ), we obtain by interpolation that

∀ 1/2 < θ < 1, ∀i > j ∈ N, T 0 h i -h j 2 H θ (k+2) dt ≤ T 0 h i -h j 1-θ L 2 h i -h j θ H k+2 dt ≤ h i -h j 1-θ C T (L 2 ) T 0 h i -h j θ H k+2 dt ≤ h i -h j 1-θ C T (L 2 ) T 2/(2-θ ) h i -h j θ L 2 T (H k+2 ) ≤ h i -h j 1-θ C T (L 2 ) T 2/(2-θ ) (2C 1 ) θ -→ n→+∞ 0.
Therefore ∀k/2 < s < k, the sequences (h ih re f ) and

(z i ) are Cauchy sequences in L 2 T (H s+2 ), thus h - h re f , z ∈ L 2 T (H s+2 ). Moreover ∂ t h i converges to ∂ t h in L 2 T (L 2 ) and (∂ t h i ) is uniformly bounded in L 2 T (H k ), similarly for ∂ t z. Consequently ∀s < k, ∂ t h, ∂ t z ∈ L 2 T (H s ), thus by Proposition 3.5, h -h re f , z ∈ C T (H s+1 ).

Finally the a priori estimates on h and z allow to conclude that

h -h re f , z ∈ L 2 T (H k+2 ) ∩ C T (H k+1 ). In particular, as k + 1 = 4, (h i -h re f ), (∇h i ), (∇ 2 h i ), (∂ t h i ) and (z i ), (∇z i ), (∇ 2 h i ), (∂ t z i ) converges in C([0, T ] × R 2
). Now we consider (c i ) i∈N . By Proposition 3.11, there exists c ∈ C T (L 2 ) limit of (c i ) i∈N in this space. We know that (c i ) i∈N is uniformly bounded in C T (H k ), so by interpolation:

∀ 0 < θ < 1, ∀i > j ∈ N, sup t∈[0,T ] c i -c j H θ (k) ≤ sup t∈[0,T ] c i -c j 1-θ L 2 sup t∈[0,T ] c i -c j θ H k ≤ c i -c j 1-θ C T (L 2 ) ( c i C T (H k ) + c j C T (H k ) ) θ ≤ C c i -c j 1-θ C T (L 2 ) -→ n, m→+∞ 0.
Therefore (c i ) i∈N is a Cauchy sequence, and thus converges to c in the space C T (H s ), for all s < k. Moreover (∂ t c i ) converges to ∂ t c in C T (L 2 ) and is uniformly bounded in C T (H k-1 ), so it converges to c in C T (H s-1 ). And we conclude by the a priori estimate on c that c ∈ C T (H k ). To conclude, (c i ) i∈N , (∂ t c i ) i∈N and (∇c i ) i∈N converge in C T (R 2 ) thus we can take the limit in the equations, and (h, z, c) is solutions of System (9). This concludes the proof of the existence.

Uniqueness

Proposition 3.12. Let (h 1 , z 1 , c 1 ) and (h 2 , z 2 , c 2 ) be two solutions of System (9) satisfying the hypotheses of Theorem 1.1. Then ∀t ∈ [0, T ] :

h 1 (t) -h 2 (t) 2 L 2 + z 1 (t) -z 2 (t) 2 L 2 + c 1 (t) -c 2 (t) 2 L 2 ≤ h 0 1 -h 0 2 2 L 2 + z 0 1 -z 0 2 2 L 2 + c 0 1 -c 0 2 2 L 2 e CT
In particular when the initial conditions are the same for both solutions, these solutions are the same. Consequently this proposition shows the uniqueness of the solution (h, z, c) of the theorem.

Proof. We first write the equations verified by δ h := h 1h 2 , δ z := z 1z 1 and δ c := c 1c 2 :

           ∂ t δ h -div(h 1 ∇δ h + δ h∇h 2 ) -div(h 1 ∇δ z + δ h∇z 2 ) = 0 ∂ t δ z -Kδ ∆z = sδ c -e (h m 1 (|v 1 | n -|v 2 | n ) + |v 2 | n (h m 1 -h m 2 )) ∂ t δ c + δ v∇c 2 + v 2 δ ∇c = e h m-1 1 (|v 1 | n -|v 2 | n ) + |v 2 | n (h m-1 1 -h m-1 2 ) -(s + r) c 1 h 1 - c 2 h 2 (32a) (32b) (32c)
Then the bound is obtained with a similar process as in the proof of Proposition 3.11.

Spectral stability of constant states

In this section, we consider the flow over a topography that is an inclined plane at time t = 0. We assume that r = 0 and we study the spectral stability of constant states. We expect that instability will provide a mechanism for pattern formation. We first write System (8) in a non-dimensional form and then linearize this system around constant states. Then we explore numerically the stability of the system. Finally, we carry out the spectral stability analysis by using Routh-Hurwitz theorem: this provides necessary and sufficient conditions for constant states to be spectrally stable. However, these conditions do not provide any insight on the nature of the instabilities. We complete this analysis by an asymptotic expansion of the spectrum around the origin and in the high frequency regime.

Non-dimensionalization and linearization of the system

We write System (8) in a non-dimensional form in order to identify the important parameters. We introduce several characteristic quantities : Z is a characteristic eroded height, H a characteristic water height, L a characteristic wavelength and T a characteristic time. We chose T to be the necessary time to erode the soil of a height Z, with an erosion speed e. Thus T verifies e T = Z. As v = µ tan θ e 1 -µ∇(z + h), we fix the characteristic water velocity V = µ. We introduce the dimensionless variables :

h := h H , z = z Z , v := v V , c := c c sat x := x L , t := e Z t.
In order to simplify the notations, we will assume H = Z. Then, dropping the primes, System (8) is written as:

               ∂ t h + ZV tan θ eL ∂ x h = Z 2 V eL 2 div(h∇(h + z)), h∂ t c + ZV tan θ eL h∂ x c = ZV eL 2 h∇(h + z).∇c + ρ s c sat h m | tan θ e 1 - Z L ∇(h + z)| n - s e ρ s c sat c, ∂ t z = ZK eL 2 ∆z -h m | tan θ e 1 - Z L ∇(h + z)| n + s e c.
To simplify the equations, we set Z L = e V , and we define α := Z L = e V , K := ZK eL 2 = K LV . The system reads :

         ∂ t h + tan θ ∂ x h = α div(h∇(h + z)), h∂ t c + tan θ h∂ x c = αh∇(h + z).∇c + ρ s c sat h m | tan θ e 1 -α∇(h + z)| n - ρ s s c sat e c, ∂ t z = K∆z -h m | tan θ e 1 -α∇(h + z))| n + s e c. (33)
The stationary states of Equation (33) for a flat surface, denoted by (h, c, z) verify:

∀t ∈ R + , ∀(x, y) ∈ Ω,      h(t, x, y) = h > 0, c(t, x, y) = c = e s h m tan n θ > 0, z(t, x, y) = 0.
This means that the erosion and the deposition process equilibrate each other and the bottom is not eroded, whereas the fluid height is a constant.

Let (h + h, c + c, z), with |h|, |c|, |z| 1, be a small perturbation of the constant state, and solution of System (33). Then, at first order, this solution verifies the following linear system:

             ∂ t h + tan θ ∂ x h = αh ∆(h + z) ∂ t c + tan θ ∂ x c = ρ s c sat s e mc h 2 h - αnc h tan θ ∂ x (h + z) - c h ∂ t z = K∆z - s e mc h h - αnc tan θ ∂ x (h + z) -c (34)
We denote f = (h, c, z) T . Then f verifies the equation

∂ t f = A 0 f + A 1 ∂ x f + A 2 ∆ f with: A 0 =        0 0 0 amc h -a 0 - amc ρs ah ρs 0        , A 1 =       -tan θ 0 0 - αanc tan θ -tan θ - αanc tan θ αanh c ρs tan θ 0 αanh c ρs tan θ       , A 2 =   αh 0 αh 0 0 0 0 0 K   ,
where we have denoted a = sρ s ehc sat and ρs = ρ s c sat .

We apply the Fourier transform in space and the equation verified by f , the Fourier transform in space of f , is:

∂ t f = A 0 + iξ A 1 -(ξ 2 + η 2 )A 2 f := A(ξ , η) f .
Consequently, in order to study the stability of the system, we have to determine the sign of the real part of the eigenvalues of the matrix A(ξ , η). These eigenvalues are denoted by λ i (ξ , η) with i ∈ {1, 2, 3} and the associated eigenvectors are denoted by V i (ξ , η). The expressions for λ i are not explicit: in the next Subsection, we compute numerically the stability of the system. Then Subsection 4.3 gives a stability result on the domain. This result is completed with the asymptotic study of the eigenvalues at low (|ξ | 2 + |η| 2 1) and high (|ξ | 2 + |η| 2 1) frequencies.

Numerical exploration of stability

In this section, we explore numerically the stability of the system (34). For that purpose, we have computed numerically the three eigenvalues of the matrix A(ξ , η), and the system is stable if and only if the real parts of these eigenvalues are negative. Since A(ξ , -η) = A(ξ , η) and A(-ξ , η) = A(ξ , η), the real part of the spectrum remains unchanged under the transformation ξ → -ξ and η → -η. Thus, we examine the behaviour of the system in the top right quarter of the plane. The choice of parameters for these computations is the same as in Section 5 (unless otherwise specified), see Table 1 for their values.

In Figure 2, we have represented an illustration of the stability of the system when K > 0, where

K e = 5 × 10 -4 3600 m 2 s -1 .
The domain represented is a bounded subset of the plane (ξ , η) ⊂ R 2 , and the color represent the stability: the system is stable in the green area, and unstable in the red area. We clearly see the stabilizing effect of the creep effect: when K is higher, the stable area is larger. It seems that the unstable area is bounded. This will be confirmed by Proposition 4.2. Then, in Figure 3 there is no creep effect: K = 0, and we represent various stability diagrams for several values of the ratio n/m. The value of m is the same as in Table 1: m = 1.6, and n vary between m/2 and 10m. We observe that the behaviour of the system changes with the ratio m/n, but the unstable area always seems to be unbounded.

When n m as in Figure 3a, the stable area is bigger than the unstable area, and the system is unstable only for perturbations of transverse dominant direction. When n ≈ m, as in Figure 3c, the system is stable only for longitudinal perturbations. Finally, when n < m the system seems to be unstable at all frequencies, as in Figure 3d.

The stability exploration of the system (34) should be quantified by theoretical results, in order to identify the types of instabilities and the transitions between stability and instability. This is done in the following Section.

Stability analysis

The following theorem provides a necessary and sufficient condition for the stability of the system, depending on the parameters of the model.

If K = 0, System (34) is spectrally unstable. More precisely, for (ξ , η) ∈ R 2 * := R 2 \ {(0, 0)}, the system is stable if and only if m < ρs c and η 2 m < ξ 2 (n -m).
In particular, if n ≤ m then the system is always unstable.

The proof of Theorem 4.1 is given in Appendix A by using Routh-Hurwitz criterion.

This theorem confirms the prominent role of the creep effect that was already observed for the well-posedness property of System (9) since Theorem 1.1 required K to be sufficiently large. Here, it plays a crucial role in the stability of the system. Indeed, the condition K ≥ αmhc/ ρs is verified only if the creep effect on the soil is large enough. Therefore, this effect has a stabilising effect on the system, as observed in the previous section where the unstable region shrinks to 0 as the constant K > 0 increases. Conversely, the two conditions of Theorem 4.1 mean that the constants α, m and n have a destabilising effect on the system. Recall that α measures the ratio between the erosion speed and the water speed. When α is large, the erosion speed is large compared to the fluid velocity which increases the instability of the bottom surface.

When K = 0, the constant m still has a destabilising effect. The other condition implies that arbitrary high frequencies are unstable. When n < m then the system is unstable on the whole domain R 2 . When n ≥ m, the system is stable in an unbounded area of the spectrum :

(ξ , η) ∈ R * × R ; m ≤ ξ 2 ξ 2 + η 2 n .
This area is delimited by the lines of equation

η = ± n -m m ξ ,
thus increases with n: this confirm the observation made in Figure 3. In particular, if n m then the system destabilises only if ξ η, that is for perturbations transverse to the water flux. Moreover, as long as n ≥ m, the system is unstable in the transverse direction (ξ = 0). This may lead to the formation of rills in the direction of the water flux. Theorem 4.1 provides a stability criterion but does not give any insight on the nature of the instabilities. The following propositions determine more precisely the spectrum in the small wavenumber regime ξ 2 + η 2 1 and in the large wavenumber regime ξ 2 + η 2 1. Since the system is spectrally stable if K ≥ αmhc/ ρs , we focus on the case K < αmhc ρs in order to detect low frequency instabilities. The following proposition provides a more complete picture of the stability diagram when K < αmhc/ρ s in the limit ξ 2 +η 2 → 0. The expansion of the eigenvalues depends on the relative position of ξ and η 2 . We split the analysis into two cases: ξ = O(η 2 ) and η 2 = o(ξ ).

Proposition 4.1 (Instabilities at low frequencies). Assume K < αmhc/ ρs . In the limit |ξ | + |η| → 0, System (34) has the following stability properties:

• When ξ = O(η 2 ) then:
under the assumption αmhc/ ρsαh < K < αmhc/ ρs , the system is stable if and only if

ξ 2 < K ρs tan 2 θ (K + αh -αmhc) 2 αmhc -K η 4 .
-If K ≤ αmhc/ ρsαh, the system is unstable for all these frequencies.

• When η 2 = o(ξ ) then the system is stable if and only if

(K ρs + αhc(n -m)) ξ 2 > (αmhc -K ρs ) η 2 .
In particular, if K < α(mn)hc/ρ s the system is unstable at these frequencies. Moreover the most unstable eigenvalue expands as

λ 1 = ξ 2 +η 2 →0 -K + c sat ρ s αnhc - c sat ρ s αmhc ξ 2 -K - c sat ρ s αmhc η 2 + o ξ 2 + η 2 .
Remark 4.1. This proposition means that the system is stable in all directions at low frequencies if and only if K ≥ αmhc/ ρs , and we recover the global criterion.

The consequences of Proposition 4.1 are fully discussed in Section 4.4. However, we can make some preliminary comments. We first note that when K → αmhc/ ρs then

K ρs tan 2 θ (K + αh -αmhc) 2 αmhc -K → +∞ and αmhc -K ρs K ρs + αhc(n -m) → 0.
Consequently, the stable area (the set of frequencies (ξ , η) such that the system is stable at these frequencies) increases with K, until filling all the low frequencies. Conversely, when K → 0 the stable area decreases, up to the area described in the second part of Theorem 4.1, when K = 0. In the unstable regime K < αmhc/ ρs , we can precise the instability scenario. In the case ξ = O(η 2 ), the behaviour of the function f defined as

f (K) := K (αh + K -αmhc/ ρs ) 2
αmhc/ ρs -K is described in Figure 4. On the interval [mhc/ ρsαh, αmhc/ ρs ], the function f is increasing from 0 to +∞, thus for ξ , η fixed there exists a unique K(ξ , η) such that

ξ 2 = 1 ρs tan 2 θ f ( K(ξ , η))η 4 .
If K ≤ K(ξ , η) then the system is unstable at this frequency (ξ , η), and if K > K(ξ , η) then the system is stable at this frequency. In order to prove Proposition 4.1, we compute an asymptotic expansion of the eigenvalues λ i (ξ , η), i = 1, 2 of the matrix A(ξ , η), as |ξ |+|η| goes to 0. Then, we study the sign of their real parts. Note that A(0, 0) = A 0 , which corresponds to homogeneous in space perturbations, admits λ 1 = λ 2 = 0 and λ 3 = -a as eigenvalues. We focus on the expansion of the eigenvalues λ 1 , λ 2 bifurcating from 0 as the third one bifurcates from λ 3 (0, 0) = -a < 0 and its real part remains negative for |ξ | + |η| small enough. The matrix A(0, 0) is diagonalizable and the eigenvalue 0 is semi-simple (its algebraic multiplicity, 2, is equal to its geometric multiplicity). Therefore the eigenvalues λ i (ξ , η) admit a Taylor expansion with respect to the perturbation parameters ξ and η 2 (see [START_REF] Andrew | Derivatives of eigenvalues and eigenvectors of matrix functions[END_REF]). These expansions depends heavily on the ratio ξ /η 2 : we split the analysis between the case |ξ |/η 2 bounded by a constant ,which is studied in Section B.1, and the case |ξ | 2 /η 2 → +∞, which is studied in Section B.2. The details of the proof of Proposition 4.1 can be found in Appendix B.

K f (K) 0 α mhc ρs 0 0 +∞ +∞ (a) mc ≤ ρ s K f (K) 0 K 0 α mhc ρs -αh α mhc ρs 0 0 f (K 0 ) f (K 0 ) 0 0 +∞ +∞ (b) mc > ρ s
The following proposition determine the stability of System (34) in the limit ξ 2 + η 2 → +∞.

Proposition 4.2 (Stability analysis at high frequencies). The stability results for the system at high frequencies are divided in three cases:

• When K > 0 then System (34) is stable at high frequencies. The eigenvalues expand as

       λ 1 (ξ , η) = -K(ξ 2 + η 2 ) + o(ξ 2 + η 2 ) λ 2 (ξ , η) = -αh(ξ 2 + η 2 ) + o(ξ 2 + η 2 ) λ 3 (ξ , η) = -iξ tan(θ ) - s ρs eh + o(1)
.

• If K = 0 and ξ = 0, the Taylor expansion of the eigenvalues when ξ 2 + η 2 → +∞ is given by:

               λ 1 (ξ , η) = sc eh m - ξ 2 ε 2 n + o(1), λ 2 (ξ , η) = -αh(ξ 2 + η 2 ) + o(ξ 2 + η 2 ), λ 3 (ξ , η) = -iξ tan θ - s ρs eh + o(1)
.

• If K = 0 and ξ = 0, the Taylor development of the eigenvalues when η → +∞ is:

           λ 1 (0, η) = 0, λ 2 (0, η) = -αhη 2 + o(η 2 ), λ 3 (0, η) = sc eh m - ρs c + o(1).
The proof of this proposition is postponed in Appendix C.

As a consequence of Proposition 4.2, one finds that when K > 0, the unstable domain in the spectrum is bounded, as observed in Section 6. Thus, there exists a wave vector associated to an unstable eigenvalue with a maximum real part which may provide a description of the pattern geometry. When K = 0, the unstable region is unbounded but there is also a most unstable eigenvalue with its real part bounded by sc eh m: this is obtained when ξ = 0 for pure transverse perturbations.

Discussion on the low frequencies stability analysis: form of the spectrum

The stability result on the whole domain Theorem 4.1 does not give an explicit formula for the localisation of the stable and unstable areas when K > 0. Therefore in this section we interpret the stability analysis at low frequencies Proposition 4.1 to study the form of the limit between stable and unstable areas in the spectrum, at these frequencies, for K > 0. As the real part of the spectrum remains unchanged under the transformations ξ → -ξ and η → -ξ , we examine the behaviour of the system in the top right quarter of the plane. When K ≥ αmhc/ ρs , the system is stable at all low frequencies, thus we focus on the case K < αmhc/ ρs . We first examine the case of the frequencies which verify η 2 = o(ξ ). If K ρs ≤ α(mn)hc, the system is unstable at these frequencies. When α(mn)hc/ ρs ≤ K ≤ αmhc/ ρs , the instability comes from the term in η of the first eigenvalue, and: Consequently the boundary between the stable and unstable areas of the system is close to a straight line, of slope

λ 1 ≥ 0 ⇐⇒ η 2 ξ 2 ≥ K ρs -α(m -n)hc αmhc -Kρ s + o ξ 2 + η 2 .
η ξ = K ρs -α(m -n)hc αmhc -K ρs . ( 35
)
The system is stable below this line, and unstable above. The other case corresponds to the frequency domain ξ = O(η 2 ). If K ≤ αmhc/ ρsαh the system is unstable. When α(mn)hc/ ρs ≤ K ≤ αmhc/ ρs the destabilizing effect is given by the first eigenvalue, and

λ 1 ≥ o ξ 2 + η 2 ⇐⇒ ξ 2 tan 2 θ η 4 ≤ K ρs (K ρs + αh ρs -αmhc) 2 αmhc -K ρs + o ξ 2 + η 2 .
Therefore the boundary between stable and unstable area in this case is close to the curve of equation

η = ρs K αmhc -K ρs (K ρs + αhρ s -αmhc) 2 1/4 tan θ ξ . ( 36
)
Consequently, when min(αmhc/ ρsαh, α(mn)hc/ ρs ) ≤ K ≤ αmhc/ ρs , the system is stable below the line of slope given by (35) and above the curve of Equation (36), and unstable between these curves. This situation is illustrated in Figure 5, for various values of K. The black curves are the boundary curves between stable and unstable areas calculated above, we can see that they fit the calculations. The unstable area are bounded, and this is confirmed by the stability analysis at high frequencies below. Moreover, as K increases, the unstable area decreases.

If K is smaller than αmhc/ ρsαh or α(mn)hc/ ρs then one of the curves (or both) of equation ( 35) and (36) disappear, and the instability area increases. Figure 6 illustrates the case α(mn)hc/ ρs < K < αmhc/ ρsαh, for various values of K. The only boundary between the stable and unstable areas is the line of slope (35). The system is stable below this line and unstable above this line. As K decreases the unstable area increases, and comes closer to the line of slope n-m m which is the boundary when K = 0, given in Theorem 4.1.

Direct numerical simulations

In this section, we present some numerical experiments of the erosion of a tilted plane. The parameters of these experiments come from physical data. The quantities L x , L y , V , h, ρ s , c sat , e and θ are chosen according to the laboratory experiment [START_REF] Guérin | Streamwise dissolution patterns created by a flowing water film[END_REF], which erodes a block of salt. The choice of values for the exponents m and n has been investigated many times in the literature. The values are chosen between 0 and 3, with an additional relation n = 2m, as explained in [START_REF] Chen | Landscape evolution models: A review of their fundamental equations[END_REF] (here the constant n corresponds to the constant m + n in the literature). We choose n = 2m, with n sufficiently large in order to observe the formation of channels in the simulations. Indeed, we found that the effect of digging in depressions is reinforced when these exponents are larger. The choice of parameters is given in Table 1.

Once the ratio between the erosion speed and the water speed has been fixed, we assume that the time variations of the water height and concentration are small, consequently we neglect these variations in the simulation. Indeed, in the simulation, the eroded height of the soil is of the order of a millimeter, thus the characteristic time T = Z/e is of the order of an hour. The water crosses the domain in 0.4 seconds, thus there are four orders of magnitude between the characteristic time of the water flow and that of erosion. A direct numerical simulation of the full system, with time derivatives, would impose a severe CFL restriction: indeed, the fluid velocity is about 1m.s -1 whereas the erosion rate is around 1mm.h -1 . Since we are interested in the erosion process, the natural time scale is one hour and a typical time step would be a minute. However, the numerical time step δt is driven by a CFL: δt ≤ δ x/v water ≈ 0.001s if one consider a typical mesh size ). This increases the numerical cost of the scheme. Instead, we solve stationary problems at each time steps, where the time step is determined by the erosion time scale. A comparison between the two resolution methods is made in Appendix D, it shows that taking the stationary equation does not affect the results. The stationary equations for the water height and concentration of sediments in water are discretised with a finite volume method. The scheme is given in Appendix D. The equation (8b) on c is a linear equation, and we discretise it with an explicit Euler scheme by considering it as an evolution equation with respect to the variable x, as the speed in this direction does not vanish:

∂ x c + v y v x ∂ y c = ρ s h (E -S). (37)
The equation (8a) on h is non linear, thus it is harder to discretise it. In order to avoid an implicit discretisation of this equation, we made the choice to linearise the equation. Denoting by h n the solution of the equation at time t n , we approximate:

div(h∇h)(t n ) + div(h∇z)(t n ) ≈ div(h n-1 ∇h n ) + div(h n ∇z n ).
As the solution at the previous time step is known, the right hand term is linear in h n . Thus we can discretise it with a finite volume scheme in two dimensions. The justification for the quasi-stationary model and the numerical scheme are given in Appendix D.

We have chosen periodic boundary conditions in the transverse direction, therefore in the stability analysis there are only countable frequencies in this direction. The frequencies are the η n = 2πn/L y where L y is the width of the domain. The boundary conditions at the top of the domain are Dirichlet condition for h, c and z.

At the bottom of the domain, z is prescribed by a Dirichlet condition and we suppose that water flows freely. Thus we fix Neumann condition for h and c.

In the simulations, the initial surface is a flat tilted plane with a small random perturbation. This surface is represented in Figure 7, and it shows a flat view of the two dimensional plane. The color scale show the height difference between the actual soil, and the flat plane. Thus the yellow areas are the less dug parts. The represented quantity is the difference between the flat surface and the perturbed soil, in millimeters. This quantity is zero at yellow points, and increase when color is darker.

Simulations without source term for water

In this part, we present some results of simulations when r = 0, as in the spectral stability analysis. These simulations are compared to the theoretical results of stability. First, in Figure 8 we presents results of the simulation of the system, when

K = K e = 5 × 10 -4 3600 m 2 s -1 .
In this case, the system is spectrally stable at all frequencies {(ξ , η n ); ξ ∈ R, n ∈ N}. The pictures represent the eroded height z(T )z 0 at time T = 0.25 hour and T = 2 hours. We observe that surface perturbations are quickly smoothed, and tends to disappear. After 15 minutes, we see in Figure 8a that the amplitude of perturbation has not decreased yet, but these perturbations are smoother than initially. After 2 hours, we see in Figure 8b that the amplitude of the perturbations has decreased. In Figures 9 and10, we represent respectively the eroded surface and the fluid height when the system is unstable at some frequencies. In Figure 9, we can observe the formation of channels in the soil, in the flow direction. The width of these channels are larger when K is higher, and they take more time to appear. This can be explained by the stability analysis, as discussed below. Figure 10 have some similarities with Figure 9 since the water tends to fill the eroded channels ; the water depth is larger in the channels, and smaller between them.

Figure 11: Discrete frequency analysis of the result of the simulation, when K = K e /20 (Figure 9a), and K = K e /50 (Figure 9b). The quantity represented is the norm of the discrete Fourier transform of the signal in the transverse direction.

We quantify the number of channels generated by the simulations in order to validate the stability analysis.

For that purpose, we compute the discrete Fourier transform of the eroded surface at the end of the simulation.

As the main direction of perturbations is transverse to the slope, we have calculated this Fourier transform in this direction, at fixed x. In Figure 11, we represent the norm of the discrete Fourier transform (DFT) of the result of the simulations 9a, 9b in the transverse direction, calculated at the bottom of the plane (for x = L x ). The frequencies are the 2πN/L y , for N ∈ {0, . . . , Ny/2}. When K = K e /20, we can see in Figure 11 that the dominant frequency is reached when N = 2 at the bottom of the tilted plane. When K = K e /50, we have contribution between N = 2 and N = 5 frequency Then, in order to compare this observation with the stability analysis, we compute numerically the frequencies which create the most unstable modes. In Figure 12, we present numerical computations of the eigenvalues of the linearised system, depending on the frequencies. The largest real part of the three eigenvalues is represented, for each frequencies ξ and η in a bounded domain. This quantity control the stability, the system is stable if and only if it is negative. Due to the periodic boundary conditions, the frequencies allowed in the transverse direction are the η n = 2πn/L y , where L y is the width of the domain. In Figure 12, the white lines shows these frequencies. For both cases Figure 12a and 12b, the system is unstable and the maximum of instability is reached at ξ = 0, η 1 = 2π/L y . Thus the system destabilises in the transverse direction, and the most destabilising frequency has period one. This is of the same order of magnitude as computed by the DFT of the simulations, where this frequency is (ξ , η) = (0, 2 × 2π/L y ). The difference between the period of 1 predicted by the stability analysis and the periods of 2 -5 obtained in the simulations could come from the non linear effects of the model, that are not taken in account in the stability analysis.

Simulations including a source term for water

In this part, we present some results of simulations when the source term r is a positive function, to observe the effect of rain in the model. In the following simulations, K = K e /20. First, we consider a constant in time, and uniform in space source term r 1 = 0.005mm/s. Figure 13 shows the surface height (left figure) and water height (right figure) computed by the simulation, at time 0.25 hours.

We can observe that the water height is almost four times higher at the bottom of the domain (x = L x ), than at the top (x = 0). As the erosion rate is proportional to a power of the water height h, the surface erodes faster at the bottom of the domain.

Figure 14: Results of the simulation when r = r 1 , after 16000 iterations.

Then, Figure 14 shows the results of the same simulation, at time 2 hours. The initial perturbations have almost disappeared, and some transverse perturbations have developed at the bottom of the domain. Next, we choose a bigger source term: r 2 = 0.01m/s, and results of this simulation are shown in Figure 15. We can see that there are bigger perturbations at the bottom of the domain than with the source term r 1 . These perturbations are also not totally transverse: they undulate a little bit in the longitudinal direction, which is a different behaviour from the case without source term. Now we take a constant non uniform source term:

r 3 (x, y) = 0.005 exp -60 * (y -Ly/2) 2 Ly 2 ,
that correspond to a positive longitudinal band of rain in the middle of the domain which decrease exponentially fast when approaching y = 0 or y = 0.1. At time 0.75 hours, the surface is highly eroded at the bottom of the band. Here, the evolution of the system is driven by the source term, the erosion landscape depends mainly on this term. These simulation show that the source term may quantitatively influence the results, although qualitatively they look similar to the case with no source term.

Conclusion and perspectives

In this paper, we have considered a model for the evolution of landscape subject to water erosion, in order to study the formation of patterns. This model takes into account the water flow, the dissolved sediments and the main effects of erosion and sedimentation, while remaining simple enough to be studied both theoretically and numerically. We proved that under realistic hypotheses, the system is well posed locally in time. Then a complete spectral analysis of a linearization of the system around stationary solutions has been performed. This analysis highlights the various behaviors of the system depending on the parameters and on the frequencies of the perturbations. A very important parameter of the system is the constant of creep K, which controls the stability. The system can become stable if this constant is large enough, and is unstable if this constant is too small. Moreover, the instabilities in the system appears in the transverse direction, and this leads to the formation of channels parallel to the water flow.

This analysis is a preliminary step in the study of pattern formation on erodible surfaces. We have shown that an instability mechanism can explain the formation of parallel channels at the early stages of erosion. We plan to perform a more complete parametric study of the system, in order to understand the various behaviors of the model. This analysis should be completed by numerical simulation to illustrate these behaviors. There are several interesting direction of research. First, the analysis of this model could be extended to more general landscapes like mountains or valley. In these cases, the well-posedness of the system in short time can be still valid, as long as the water level and the water speed does not vanish. A similar stability analysis could be carried out if steady states exist. One could also consider more involved models. Indeed a lot of factors are not taken into account as weather, vegetation, animal and human activities. It is be quite difficult to include these factors in the model, but a possible amelioration would be to include randomness in the equations. A random term could allow to model these factors which fluctuate over time. Moreover, a non constant source term, that would model the variations of rain in time would make the model more realistic.

Another potential improvement of the model would be to consider more complex laws for the fluid velocity or even consider shallow water type models for the evolution of the layer of fluids.

A Proof of the stability theorem

This section is devoted to the proof of Theorem 4.1. First we state the Routh Hurwitz criteria for complex polynomial of degree 3, proven in [START_REF] Morris | The routh and routh-hurwitz stability criteria: Their derivation by a novel method using comparatively elementary algebra[END_REF]:

Proposition A.1. Let P(X) = X 3 + (a 1 + ib 1 )X 2 + (a 2 + ib 2 )X + a 3 + ib 3 be a polynomial with complex coefficients. Then the roots of P have a positive imaginary part if and only if the three following conditions are satisfied :

                       -∆ 2 = - 1 a 1 0 b 1 > 0 ∆ 4 = 1 a 1 a 2 a 3 0 b 1 b 2 b 3 0 1 a 1 a 2 0 0 b 1 b 2 > 0, -∆ 6 = - 1 a 1 a 2 a 3 0 0 0 b 1 b 2 b 3 0 0 0 1 a 1 a 2 a 3 0 0 0 b 1 b 2 b 3 0 0 0 1 a 1 a 2 a 3 0 0 0 b 1 b 2 b 3 > 0
This proposition is directly used to prove Theorem 4.1.

Proof of Theorem 4.1 . The characteristic polynomial of A, denoted by P, is :

P(X) = X + iξ tan θ + αε 2 h 0 αε 2 h - amc h + αiξ anc tan θ X + a + iξ tan θ αiξ anc tan θ amc ρs -αiξ anhc ρs tan θ - ah ρs X -αiξ anhc ρs tan θ + ε 2 K = X + iξ tan θ + αε 2 h 0 αε 2 h - amc h + αiξ anc tan θ X + a + iξ tan θ αiξ anc tan θ 0 hX + iξ h tan θ ρs X + ε 2 K
.

where we have denoted ε 2 = ξ 2 + η 2 . In order to simplify the computations, we introduce the variables

Y := X + iξ tan θ , ξ = ξ tan θ , ε2 = αhε 2 , K αh = K , ρs h = h.
We also define the constants N = αanc/ tan 2 θ , M = amc/h. Consequently the polynomial writes as:

P(Y ) = Y + ε2 0 ε2 -M + iξ N Y + a i ξ N 0 hY Y + ε2 K -i ξ = -hY i ξ NY + ε2 M + (Y + a)(Y + ε2 )(Y -i ξ + ε2 K) = Y 3 + a + ε2 (1 + K) -i ξ (1 + N h) Y 2 + ε2 (a(1 + K) -hM + K ε2 ) -i ξ (a + ε2 ) Y +a ε2 (ε 2 K -i ξ ). ( 38
)
The system is stable if and only if the roots of X → P(X) have a negative real part. As Y have the same real part as X, this is equivalent to the fact that the roots of Y → P(Y ) have a negative real part. Denoting iλ := Y , the system is stable if and only if the roots of λ → P(iλ ) have a positive imaginary part. Thus we will apply the Routh Hurwitz criteria to the polynomial λ → Q(λ ) := iP(iλ ) where iλ := Y :

Q(λ ) =λ 3 -ξ (1 + N h) + i a + ε2 (1 + K) λ 2 -ε2 (a(1 + K) -hM + K ε2 ) -i ξ (a + ε2 ) λ + a ε2 ( ξ + i ε2 K).
We denote:

ā1 = 1 + N h, b1 = a + ε2 (1 + K), ā2 = a(1 + K) -hM + ε2 K, b2 = a + ε2 . Thus Q(λ ) = λ 3 -ξ ā1 +i b1 ) λ 2 + -ε2 ā2 +i ξ b2 λ + a ε2 ( ξ + i ε2 K).
The first condition : The determinant -∆ 2 writes as :

-∆ 2 = - 1 -ξ ā1 0 -b1 = b1 = a + ε2 (1 + K).
Consequently, as a > 0 and K ≥ 0, the condition -∆ 2 > 0 is always satisfied.

The second condition : The determinant ∆ 4 is given by :

∆ 4 = -b1 ξ b2 ε4 a K 1 -ξ ā1 -ε2 ā2 0 -b1 ξ b2 = -b1 ξ ( b2 -ā1 b1 ) ε2 (ε 2 a K -b1 ā2 ) 1 0 0 0 -b1 ξ b2 = ε2 ξ 2 ε2 b2 ā1 b1 -b2 + b1 ā2 b1 -ε2 a K . Let t = ξ 2 ε2 ∈ [0, µ tan 2 θ eh ], then one finds that ∆ 4 > 0 ∀ (ξ , η) ∈ (R 2 ) * if and only if t b2 ā1 b1 -b2 + b1 ā2 b1 -ε2 a K > 0 ∀t ∈ [0, µ tan 2 θ eh ], ∀ ε2 > 0.
The terms b1 and b2 are positive, and :

ā1 b1 -b2 = a + ε2 (1 + K) + N h(a + ε2 (1 + K)) -a -ε2 = ε2 K + N h a + ε2 (1 + K) > 0.
Consequently, one has

∆ 4 > 0, ∀ (ξ , η) ∈ (R 2 ) * if and only if ā2 b1 -ε2 a K ≥ 0 ∀ ε2 > 0.
We compute:

ā2 b1 -ε2 a K = a(1 + K) -hM + ε2 K a + ε2 (1 + K) -ε2 a K = a a(1 + K) -hM + ε2 (1 + K) a(1 + K) -hM + ε2 K As a consequence, one has ∆ 4 > 0 ∀ (ξ , η) ∈ (R 2 ) * if and only if a(1 + K) -hM ≥ 0. This condition is equivalent to: ρs c + K ρs αhc -m ≥ 0.
The third condition : The determinant ∆ 6 is given by:

∆ 6 = -b1 ξ ( b2 -ā1 b1 ) ε2 (ε 2 a K -ā2 b1 ) ξ ε2 a b1 0 -b1 ξ b2 ε4 a K 0 1 -ξ ā1 -ε2 ā2 ξ ε2 a 0 -b1 ξ b2 ε4 a K = -b1 2 ε2 (ε 2 a K -ā2 b1 ) + ξ 2 ā1 ( b2 -ā1 b1 ) ξ ε2 a b1 + ξ ε2 ā2 ( b2 -ā1 b1 ) -ξ 2 ε2 a( b2 -ā1 b1 ) ξ ( b2 -ā1 b1 ) ε2 (ε 2 a K -ā2 b1 ) ξ ε2 a b1 -b1 ξ b2 ε4 a K .
We can factorize out the term ε2 in the third column. Then we develop this expression with respect to the third row. One has:

∆ 6 = -ε 6 b 2 1 T 0 + tT 1 + t 2 T 2 with T 0 = ε2 a K ε2 a K -ā2 b1 2 > 0, T 1 = ε2 a K ε2 a K ā1 b2 -ā1 b1 + ā1 b1 -b2 a b1 + ā2 b2 -a b1 b2 +a b1 ā1 b1 -b2 + a ā2 b1 2 b2 -a 2 b1 3 = ε2 a K ā1 b1 -b2 2a b1 + ā2 b2 -ε2 a ā1 K -a b1 b2 + a b1 2 ā2 b2 -a b1 , T 2 = a b2 2 ā1 b1 -b2 > 0.
We thus have to determine the sign of T 1 . First, one finds that

ā1 b2 -a b1 = a ( ¡ 1 + K) -hM + ε2 K a + ε2 -a ¡ a + ε2 ( ¡ 1 + K) = a K -hM + ε2 K a + ε2 -ε2 a K = a a K -hM + ε2 a K -hM + ε4 K. The constant term in T 1 is a 4 a K -hM . When ε2 = o(t), that is ε 4 = o(ξ 2
), the dominant term in ∆6 is the constant term of T 1 times t, that is a 4 a K -hM t. Thus a necessary condition for the positivity of ∆6 is a K ≥ hM. We can write -∆ 6 b1 2 ε6 := -∆6 as a polynomial of degree two in a 1 , with a negative coefficient in front of a 2 1 :

-∆6 = -ε 4 a 2 K 2 b1 t ā1 2 + ε2 a K b1 t 2 ā1 b1 + ā2 b2 ) + ε2 a K b2 + a b1 b2 2 t 2 ā1 + a b1 2 ā2 b2 -a b1 t -ε2 a K b2 2a b1 + ā2 b2 + a b1 b2 t -a b2 3 t 2 + ε2 aK ā2 b1 -ε2 aK 2 .
The term ā1 writes as ā1 = 1 + N h ≥ 1. When N = 0 and a K ≥ hM, we can verify that ∀(ξ , η) ∈ (R 2 ) * , -∆6 > 0. Consequently when a K ≥ hM the polynomial -∆6 (ā 1 ) has two roots, the first one x 1 < 1 and the second one x 2 > 1. Now we suppose that a K ≥ hM, and we define

δ = inf{x 2 (ξ , η)); (ξ , η) ∈ (R 2 ) * } ≥ 1.
According to the asymptotic calculus of stability, when a K ≥ hM the system is stable near 0 and +∞. Consequently δ is an infimum of a continuous function on a compact set of R 2 , thus it admits a minimum.

To conclude, ∆ 6 > 0 ∀(ξ , η) ∈ (R 2 ) * if and only if ā1 < δ and a K ≥ hM. That is, denoting γ = tan 2 θ (δ -1)/(αnhc), this condition reduces to n < γ and αmhc ≤ ρs K.

The third condition when K = 0: In this case, the computations on ∆ 6 are fully explicit. The third determinant reads:

- ∆ 6 b1 2 ε6 = ta b1 2 ā2 b2 -a b1 + t 2 a b2 2 ā1 b1 -b2 = ta h a + ε2 3 (-M + tN) .
Thus the stability is directly related to the sign of tN -M.

As t = ξ 2 ε2 = tan 2 θ αh ξ 2 ε 2 , and M = amc h , N = α anc tan 2 θ , then this condition reads m < ξ 2 ε 2 n.
This concludes the proof of Theorem 4.1.

B Proof of the stability results at low frequencies

In this section we prove Proposition 4.1.

B.1 The first case

In this part, we consider the case ξ → 0, η → 0 with ξ = O(η 2 ). We write the characteristic polynomial of the matrix A(ξ , η), and use the fact that λ i (0, 0) = 0 (i = 1, 2) to calculate the two roots of this polynomial bifurcating from zero. As in the proof of Theorem 4.3, we define the variables Y := X + iξ tan θ , ξ = ξ tan θ , ε2 = αhε 2 , K αh = K, and ρs h = h. We also define the constants N = αanc/ tan 2 θ , M = amc/h. We recall that the characteristic polynomial of A, denoted by P is given by Equation (38) :

P(Y ) =Y 3 + a + ε2 (1 + K) -i ξ (1 + N h) Y 2 + ε2 (a(1 + K) -hM + K ε2 ) -i ξ (a + ε2 ) Y + a ε2 (ε 2 K -i ξ ).
The zero order term of the eigenvalues λ 1 and λ 2 is zero, thus λ i = O(η 2 ) for i = 1, 2. Consequently λ 1 and λ 2 are solution of the approximate equation:

Y 2 + -i ξ + η2 (1 + K -hM/a) Y + η2 ( η2 K -i ξ ) = O(η 6 ),
where η2 = αhη 2 . Therefore, one has

           λ 1 (ξ , η) = - i ξ + η2 (1 + K -hM/a) 2 + 1 2 √ ∆ + O(η 3 ), λ 2 (ξ , η) = - i ξ η2 (1 + K -hM/a) 2 - 1 2 √ ∆ + O(η 3 ), ∆ = -ξ 2 + 2i ξ η2 (1 + K -hM/a) + η4 ((1 + K -hM/a) 2 -4 K) + O(η 6 ).
When ξ is of the order of η 2 , we can write ξ = q η2 with q = 0. Thus,

∆ = -q 2 + 2iq(1 + K -hM/a) + ((1 + K -hM/a) 2 -4 K) η4 + O(η 6 ).
Therefore we cannot compute explicitly the expression of its square roots, neither directly determine the sign of the real part of the eigenvalues. Consequently, we compute their sign. First, the real part of the square root of the discriminant ∆ is given by:

ℜ( √ ∆) = Re(∆) + |∆| 2
, where:

|∆| 2 = ξ 4 + 2 ξ 2 η4 2(1 -K + hM/a) 2 -(1 + K -hM/a) 2 + 4 K + η8 (1 + K -hM/a) 2 -4 K 2 + O(η 12 ).
We study the sign of ℜ(λ 1 ): since ℜ(λ 1 ) ≥ ℜ(λ 2 ), this will determine the spectral stability of the system. One has:

Re(λ 1 ) = - η2 (1 + K -hM/a) 2 + 1 2
Re(∆) + |∆| 2 + O(η 3 ) > 0 ⇔ Re(∆) + |∆| 2 > η2 (1 + K -hM/a) + O(η 3 ).
We then study two cases, depending on the sign of 1 + K -hM/a :

• If 0 ≤ K ≤ hM/a -1, that vis 1 + K -hM/a ≤ 0, then Re(λ 1 ) + O(η 3 ) > 0. This situation can occurs only if hM ≥ a.

• If K > hM/a -1, then ℜ(λ 1 ) > 0 ⇔ Re(∆) + |∆| > 2 η4 (1 + K -hM/a) 2 + O(η 6 ) ⇔ -ξ 2 -4 η4 K + |∆| > η4 (1 + K -hM/a) 2 + O(η 6 ) ⇔ |∆| 2 > ξ 2 + η4 ((1 + K -hM/a) 2 + 4 K) 2 + O(η 12 ) ⇔ ξ 2 ( hM/a -K) > η4 K(1 + K -hM/a) 2 + O(η 8 ).
Consequently, if K ≥ hM/a then Re(λ 1 ) + O(η 3 ) < 0, and if hM/a -1 ≤ K < hM/a then Re(λ

1 ) + O(η 3 ) > 0 ⇔ ξ 2 η4 + O(η 4 ) > f ( K), with f ( K) = K (1 + K -hM/a) 2 hM/a - K . (39)

B.2 The second case

In this section, we consider values of ξ and η which go to zeros, with η 2 = o(ξ ). The matrix A has no zero order terms in ξ and η, and no first order terms in η so we can write the Taylor expansion of the eigenvalues as λ (ξ , η)

= λ 0 + λ 1 (ξ ) + λ 2 (ξ 2 , η 2 ) + o(ξ 2 , η 2
). The method of the previous section gives only the first order terms of the expansion, but in this case they are imaginary. Consequently, we need to calculate the second order terms and the previous method is much more complicated in this case, so we proceed differently. The equation satisfied by λ and V is

∀ξ , η ∈ R, A(ξ , η)V (ξ , η) = λ (ξ , η)V (ξ , η).
We compute the terms of the Taylor expansion of λ and V step by step using this equation and identifying the terms of same order. The Taylor expansion of V is written as

V (ξ , η) = V 0 +V 1 (ξ ) +V 2 (ξ , η) + o(ξ 2 + η 2 )
(up to a multiplicative constant). For the following calculations we will solve equations of the form A 0 X = Y where Y ∈ Im(A 0 ). The matrix A 0 has a one dimensional image, generated by the vector X 3 (which is defined below). Consequently, the solutions X are of the form αX 3 + K, where K ∈ ker A 0 . We also need to compute the eigenvalues and the eigenvectors of A 0 . The eigenvalues are 0, 0, -a, and the associated eigenvectors are

X 1 =   0 0 1   , X 2 =   h mc 0   , X 3 =   0 ρs -h   .
Moreover, the two left eigenvectors of A 0 associated to the eigenvalue 0 are used to cancel some terms in the calculus. These vectors are:

X * 1 = [1, 0, 0] , X * 2 = [0, h, ρs ] .
Order 0 : When ξ = η = 0 the system AV = λV reduces to A 0 V 0 = λ 0 V 0 , so λ and V are respectively the eigenvalues and the eigenvectors of A 0 . Consequently

λ 1 0 = λ 2 0 = 0 and V 1 0 , V 2 0 ∈ vect(X 1 , X 2 ).
Order 1 : We identify the first order terms of the equation AV = λV :

A 0 V 1 + iξ A 1 V 0 = λ 1 V 0 , (40)
because λ 0 = 0. In order to compute λ 1 1 and λ 2 1 we multiply the system (40) by X * 1 and X * 2 . We obtain the following system, denoting V 0 = uX 1 + vX 2 :

0 λ 1 + iξ tan θ λ 1 mc(λ 1 + α iξ tan θ ) u vh = 0.
This system admits two solutions up to a multiplying factor :

λ 1 1 = 0, V 1 0 =   0 0 1   and λ 2 1 = -iξ tan θ , V 2 0 =   h mc 0   .
Then, using Equation (40) we compute V 2 1 and V 2 1 :

V

1 1 = -iξ αnc ρs tan θ X 3 + K 1 , V 2 1 = -iξ αnhc ρs tan θ X 3 + K 2 , with K 1 , K 2 ∈ ker A 0 = vect(X 1 , X 2 ).
As X 1 appears in the leading order of V 1 , we can suppose (up to renormalisation of

V 1 ) that K 1 = k 1 X 2 with k 1 ∈ R. Consequently, V 1 1 = -iξ αnc ρs tan θ X 3 + k 1 X 2 .
Similarly, as V 2 = X 2 +V 1 2 + o(ξ , η) we can write

V 2 1 = -iξ αnhc ρs tan θ X 3 + k 2 X 1 , with k 2 ∈ R.
Order 2 : The second order terms of the system AV i = λ i V i for i ∈ {1, 2} are:

A 0 V 2 + iξ A 1 V 1 -(ξ 2 + η 2 )A 2 V 0 = λ 1 V 1 + λ 2 V 0 . (41)
For i = 1, we multiply System (41) by the two left eigenvectors of A 0 associated to the eigenvalue 0 and we obtain the system :

iξ tan θ h 0 iξ tan θ mhc ρs k 1 λ 2 1 = -(ξ 2 + η 2 )αh -ξ 2 αnhc -(ξ 2 + η 2 )K ρs .
Consequently λ 1 2 = -(K + α nhc ρsα mhc ρs )ξ 2 -(Kα mhc ρs )η 2 . For the second eigenvalue we multiply System (41) by X * 1 and X * 2 , and we obtain :

0 αh -iξ ρs mhc k 2 λ 2 2 = -(ξ 2 + η 2 )h 2 -ξ 2 αnhc .
One finds λ 2 2 = -αh(ξ 2 + η 2 ).

Consequently, at low frequencies when η 2 = o(ξ ), the eigenvalues λ 2 and λ 3 have a negative real part. Thus the instability is given by λ 1 , and depends on the parameters. When K ρs > mhc then λ 1 has a negative real part so System (34) is stable around the stationary solution (h, c, z). When hc(mn) < K ρs < mhc then the system is stable in the longitudinal direction but not in the transverse direction. In the last case, when K ρs < hc(mn) then the system is unstable.

C Proof of the stability results at high frequencies

In this Appendix, we prove the stability result at high frequencies, Proposition 4.2. We use the same method as in Subsection B.2.

Proof of Proposition 4.2. The system AV = λV is written as

A 2 - iξ A 1 + A 0 ξ 2 + η 2 V (ξ , η) = - λ (ξ , η) ξ 2 + η 2 V (ξ , η). (43)
We can suppose, up to renormalisation that V (ξ , η) = O +∞ [START_REF] Adams | Sobolev spaces[END_REF]. We calculate a Taylor expansion of λ (ξ , η) ξ 2 + η 2 = Λ 0 + Λ 1 (ξ , η) + Λ 2 (ξ , η) + o(1), and of V (ξ , η) = V 0 +V 1 (ξ , η) +V 2 (ξ , η) + o( 1

ξ 2 + η 2 ).
Order 0 : When ξ 2 + η 2 → +∞, System (43) becomes A 2 V 0 = -Λ 0 V 0 . Its solutions are the eigenvalues and eigenvectors of -A 2 , so Λ 1 0 = -K, Λ 2 0 = -αh and Λ 3 0 = 0. The associated eigenvectors are

X 1 =   αh 0 K -αh   , X 2 =   1 0 0   , X 3 =   0 1 0   .
When K > 0 there is one eigenvalue bifurcating from zero, Λ 3 , and two other eigenvalues have negative real parts. When K = 0 the two eigenvalues Λ 1 and Λ 3 are bifurcating from zero. Consequently we will study these two cases.

The eigenvectors associated to the two zero eigenvalues of -A 2 are

X 1 =   αh 0 -αh   , X 3 =   0 1 0   .
Thus, there exist constants u 1 , v 1 , u 3 , v 3 such that V 1 0 = u 1 X 1 + v 1 X 3 and V 3 0 = u 3 X 1 + v 3 X 3 . The left eigenvectors associated to the zeros eigenvalues of -C are:

X 1 * = 0 0 1 , X 3 * = 0 1 0 . If ξ → +∞ then the first order terms of System (43) are:

A 2 V 3 1 - iξ A 1 ξ 2 + η 2 V 3 0 = -Λ 3 1 V 3 0 . ( 46
)
If ξ = 0 is fixed, the first order terms are:

A 2 V 3 1 - iξ A 1 + A 0 ξ 2 + η 2 V 3 0 = -Λ 3 1 V 3 0 . ( 47
)
To compute Λ 1 1 and Λ 3 1 , we multiply the equations ( 46) and (47) by X 1 * and X 3 * . Then, when ξ → +∞, one has:

   Λ 1 1 = 0, V 1 0 = X 1 , Λ 3 1 = - iξ tan θ ξ 2 + η 2 , V 3 0 = X 3 .
And when ξ = 0 is fixed,

     Λ 1 1 = scm eh , Λ 3 1 = -iξ tan θ + s ρs eh 1 ξ 2 + η 2 .
When ξ → +∞, the second order terms of System (43) are:

A 2 V 2 - iξ A 1 ξ 2 + η 2 V 1 - A 0 ξ 2 + η 2 V 0 = -Λ 1 V 1 -Λ 2 V 0 .
Then, one obtains:

Λ 1 2 = sc eh m - ξ 2 ε 2 n , Λ 3 2 = - s ρs eh 1 ξ 2 + η 2 .
As a conclusion, the asymptotic expansion of the eigenvalues of A(ξ , η), when ξ 2 + η 2 → +∞ and ξ = 0, K = 0 is: Finally, when ξ = 0 the characteristic polynomial P can be factorised: P(X) = X X 2 + (a + αε 2 h)X + αε 2 ah 1 -mc ρs .

              
Thus it has a zero root, λ 1 = 0, and the two other roots are λ 2,3 = -αε 2 ha ± (aαε 2 h) 2 + 4αε 2 amhc ρs 2 .

Consequently the asymptotic development of the eigenvalues of the system when η 2 → +∞ and ξ = 0, K = 0 are This achieves the proof of Proposition 4.2.

           λ 1 (0, η) = 0,

D Numerical scheme

In this appendix, we give the numerical scheme used for the simulations of the system.

Stationary regime for the water

The ratio between the erosion speed and the fluid velocity is small and we will suppose that the time derivatives of c and h are small. In order to justify this assumption, we rescale the equations by introducing the characteristic variables Z the eroded height, L the characteristic length, and T = Z/e the characteristic time. We define the dimensionless variables:

h := h H , z = z Z , v := v V , c := c c sat , x := x L , t := e Z t.
Therefore, dropping the primes, System (8) is written as:

                 ∂ t h + 1 ε div(hv) = r e , ∂ t (ch) + 1 ε div(chv) = ρ s c sat h m |v| n - s e c , ∂ t z = K∆z -h m |v| n + s e c, v = - µ V ∇(h + z).
Here, we defined ε := Le HV . If ε 1, as generally ρ s c sat , we can neglect the terms ∂ t h and ∂ t (ch) in the equation h et ch. Thus, assuming enough regularity on the solutions, we can write div(chv) = hv.∇c + c div(hv) = hv.∇c. Moreover we fix V = µ for simplicity. We obtain: (48)

         div(hv) =
Discretisation of the equations The constraint on h at each time div(h∇(h + z)) = 0 leads to a fully nonlinear problem that may be hard to solve numerically. Instead, we discretize this equation as div(h n-1 ∇h n ) + div(h

n ∇z n ) = -r, (49)
where h n represents the fluid height at time t n = nδt whereas the surface height at time n, z n , has been calculated by an explicit Euler method, using the solutions at time n -1, h n-1 and c n-1 : ∀1 ≤ i ≤ N x , 1 ≤ j ≤ N y , z n i, jz n-1 i, j dt = K z n-1 i+1, j -2z n-1 i-1, j + z n-1 i, j dx + K z n-1 i, j+1 -2z n-1 i, j-1 + z n-1 i, j dy -(E -S) n i, j .

We discretise the equation on h n by a centered finite volume scheme, which writes:

-r i, j = h n-1 i+1/2, j (h n i+1, jh n i, j )h n-1 i-1/2, j (h n i, jh n i-1, j ) dx 2 + h n-1 i, j+1/2 (h n i, j+1h n i, j )h n-1 i, j-1/2 (h n i, jh n i, j-1 ) dy 2 + h n i+1/2, j (z n i+1, jz n i, j )h n i-1/2, j (z n i, jz n i-1, j ) dx 2 + h n i, j+1/2 (z n i, j+1z n i, j )h n i, j-1/2 (z n i, jz n i, j-1 ) dy 2 , where h n-1 i+1/2, j = h n-1 i, j +h n-1 i+1, j

2

. Finally, Equation (37) on c is discretised by an upwind finite volume scheme, considering the variable x as a time variable. The time discretisation is an explicit Euler scheme. c i+1, jc i, j dt + max(v x [i, j], 0) v y [i, j] c i, jc i, j-1 dy + min(v x [i, j], 0) v y [i, j] c i, jc i, j+1 dy = -ρ s h i, j v x [i, j] (E -S) n i, j .

Boundary conditions We choose periodic boundary conditions in the y direction. There remains two boundaries, the top and the bottom of the tilted plane. The boundary conditions for the soil height are Neumann conditions: ∀1 ≤ j ≤ N y , z n 0, j = z n 1, j , z n N x , j = z n N x+1 , j .

The boundary conditions for the water height and for the concentration of sediments are a Dirichlet condition at the top because the incoming flow is fixed, and a free flow Neumann condition at the bottom.

h n 0, j = h 0 , h n N x , j = h n N x+1 , j , c n 0, j = c 0 , c n N x , j = c n N x+1 , j .

Figure 1 :

 1 Figure 1: Sketch of the mathematical framework for stability analysis: the water is flowing on a tilted plane.

Figure 2 :

 2 Figure 2: Stability diagrams in the plane (ξ , η) ⊂ R 2 , for K > 0. The system is unstable in the red area, and stable in the green area.

2 Figure 3 :

 23 Figure 3: Stability diagrams in the plane (ξ , η) ⊂ R 2 , for K = 0 and various values of n/m. The system is stable in the green area and unstable in the red area.

Theorem 4 . 1 .

 41 If K > 0, there exists a constant γ > 0 independent of n such that System (34) is spectrally stable at all frequencies (ξ , η) ∈ R 2 * if and only if K ≥ αmhc/ ρs and n < γ.

Figure 4 :

 4 Figure 4: Variations of the function f

Figure 5 :

 5 Figure 5: Stability of the system at low frequencies in the plane (ξ , η), when K > αmhc/ ρsαh. The system is stable in the green area and unstable in the red area. The black curves are the line of slope (35) and the curve of Equation (36) .

Figure 6 :

 6 Figure 6: Stability of the system at low frequencies in the plane (ξ , η), when K < αmhc/ ρsαh. The system is stable in the green area and unstable in the red area. The black line is the line of slope (35), the blue line is the line of slope n-m m (the boundary line when K = 0).

Figure 7 :

 7 Figure 7: Perturbation of the initial surface, on the domain [0, L x ] × [0, L y ]. The represented quantity is the difference between the flat surface and the perturbed soil, in millimeters. This quantity is zero at yellow points, and increase when color is darker.

2 Figure 8 :

 28 Figure 8: Eroded height of the soil in mm. The darkest areas are the most eroded areas.

Figure 9 :Figure 10 :

 910 Figure 9: Eroded height of the soil in mm with respect to the flat surface. The darkest areas are the most eroded areas. The width of the channel are larger for larger values of K

Figure 12 :

 12 Figure 12: Maximum of the real part of the three eigenvalues, normalised by the higher one, for (ξ , η) ∈ [0, 20] × [0, 300]. The most unstable areas are in yellow, and in the black parts the three eigenvalues have a negative real part, thus the system is stable.

Figure 13 :

 13 Figure 13: Results of the simulation when r = r 1 , after 2000 iterations.

Figure 15 :

 15 Figure 15: Results of the simulation when r = r 1 , after 14000 iterations.

Figure 16 :

 16 Figure 16: Results of the simulation when r = r 3 , after 6000 iterations.

40Conclusion:s αmhc η 2 + o ξ 2 + η 2 ,λ 2 =

 2222 Finally, the three eigenvalues of A(ξ , η) when ξ , η → 0 and ξ /η 2 → +∞ are : -iξ tan θαh ξ 2 + η 2 + o ξ 2 + η 2 , λ 3 = -sρ s ehc sat + o(1).

λ 1 (ξ , η) = sc eh m - ξ 2 ε 2 n

 12 + o(1), λ 2 (ξ , η) = -αh(ξ 2 + η 2 ) + o(ξ 2 + η 2 ), λ 3 (ξ , η) = -iξ tan θ -s ρs eh + o(1).

λ 2

 2 (0, η) = -αhη 2 + o(η 2 ), λ 3 (0, η) = s ρs eh mc ρs -1 + o(1).

  r, hv.∇c = ρ s c sat h m |v| ns e c , ∂ t z = K∆zh m |v| n + s e c, v = -∇(h + z).

Table 1 :

 1 Parameters of numerical simulations δ x = 1mm (for a channel of length 400mm

	Length of the domain	L x	40 cm
	Width of the domain	L y	10 cm
	Characteristic water speed	V = µ tan θ 1 m/s
	Initial water height	h 0	0.5 mm
	Initial sediment concentration c 0	317 g/m 3
	Exponent of friction over h	m	1, 6
	Exponent of friction over v	n	3, 2
	Density of the sediments	ρ s	2.17 × 10 6 g/m 3
	Concentration of saturation	c sat	3.17 × 10 5 g/m 3
	Erosion speed	e	0.5 mm/hour
	Sedimentation speed	s	e/2000
	Angle of the plane	θ	39 •

Acknowledgments PD holds a visiting professor association with the Department of Mathematics, Imperial College London.

Availability of Data and Materials

The code of numerical simulations generated during the current study is available at https://github.com/juliebinard/landscape_evol_finite_volum.

C.1 The case K > 0

In this part, we continue the asymptotic expansion of the eigenvalue Λ 3 . The zero order terms of the eigenvectors are V 1 0 = X 1 , V 2 0 = X 2 and V 3 0 = X 3 .

We calculate the first order term of Λ 3 . The first order terms of System (43) when ξ 2 + η 2 → +∞ depends on the asymptotic behaviour of iξ . If ξ → +∞ then the first order terms of the system are:

And if ξ is bounded, iξ A 1 and A 0 are of the same order, thus the first order terms are:

To compute Λ 3 1 we multiply the equations by the left eigenvector of A 2 associated to the eigenvalue 0, X * = X 3 . Then, after some computations:

When ξ is bounded we already computed the second order term of Λ 3 . When ξ → +∞, the second order terms of the system are:

Then, one finds:

Finally, the asymptotic expansion of the eigenvalues of A(ξ , η) when ξ 2 + η 2 → +∞ and K > 0 is:

C. [START_REF] Anand | Linear layout of multiple flow-direction networks for landscape-evolution simulations[END_REF] The case K = 0

When K = 0, the zero order terms of the eigenvalues Λ 1 and Λ 3 vanish. Thus, in order to prove the second point of Proposition 4.2, we need to compute the next order terms of these two eigenvalues.

Comparison with a scheme that solves the non stationary system In this paragraph, we drop the stationary assumption for the water height and sediment concentration, and compare the results to those where the stationary assumption is made. We solve System (8), with a finite volume scheme. A fully explicit scheme fails to solve the system, because the computed solution quickly blows up, even with a time step smaller than the one given by the CFL condition:

Therefore, we uses a semi-explicit scheme: linear terms of the water height and sediment concentration equations are implicit, and we explicit a part of the non linear terms. The discretisation in time of equations on h and c is given by:

Discretisation in space is done by a finite volume scheme, the same as for the stationary System (48). The numerical parameters are given by Table 1, and δt = 0.1s. Figure 17 shows the results of two simulations with the same parameters, the non stationary System (8) in Figure 17a and the stationary System (48) in Figure 17b. At T = 0.25 hours, we can observe that these two simulations are almost the same. This result supports the fact that the approximation of the complete System (8) by System (48), where water is in a stationary regime is valid.

Declarations

Competing interests No funding was received to assist with the preparation of this manuscript.