Well-Posedness and Stability Analysis of a Landscape Evolution Model

Well-Posedness and Stability Analysis of a Landscape Evolution Model Julie Binard1 · Pierre Degond2 · Pascal Noble1 Received: 17 November 2022 / Accepted: 4 November 2023 / Abstract 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, by either 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. Keywords Pattern formation · Instabilities · Landscape evolution model · Erosion by water · Stream incision law · Sedimentation Communicated by Jose A. Carrillo. BJulie Binard [email protected] Pierre Degond [email protected]v-toulouse.fr Pascal Noble [email protected] 1Institut de Mathématiques de Toulouse, UMR 5219, CNRS, INSA, Université de Toulouse, 31077 Toulouse, France 2Institut de Mathématiques de Toulouse, UMR 5219, CNRS UPS, Université de Toulouse, 31062 Toulouse Cedex 9, France 20 Page 2 of 55 Mathematics Subject Classification 35M30 ·35Q86 ·35B35 ·35B36 ·35A01 ·86-10 1 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 timescales 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 Davis (1892). It was used in Gilbert (1909) 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 Culling (1963). 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 Gilbert (1877). 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 Culling (1960), 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 Chen et al. (2014a,b) 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 qsfollows a law that depends on the water discharge q. This transport discharge law is described, for example, in Smith and Bretherton (1972) and is given by qs=kqm|∇z|n(1) Page 3 of 55 20 where k,mand nare 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 Howard and Kerby (1983). The sediments removed from the soil by the erosion are supposed to be dissolved in water and thus are transported by the water flux. This stream incision law is described in Sect.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 Chen et al. (2014a), 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: ⎧ ⎨ ⎩ ∂th−div(h∇(h+z)) =r(t,x), ∂tz=Kz+sc −ehm|v|n, ∂t(ch)+div(chv) =ehm|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 ris 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 Kis the constant of creep, whereas eand sare, respectively, the constants for the incision law and the sedimentation rate. As in Howard and Kerby (1983), we will suppose that the incision law depends on the power of the norm of the water velocity and on the 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 20 Page 4 of 55 Fig. 1 Sketch of the mathematical framework for stability analysis: The water is flowing on a tilted plane 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 and 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 ⊂R2. In Chen et al. (2014a), the authors reviewed various landscape evolution models and focused on system (2). 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 Lebrun et al. (2018), 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 Fernández-Nieto et al. (2017) for a formal derivation of these models from bilayertype models (one layer for the fluid and one layer for the sediment) and Fernandez- Nieto et al. (2014); Escalante et al. (2021) 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 present 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 presented in Fig. 1, where the tilted plane is seen from the side. 123 Page 5 of 55 20 Although it is a less realistic framework than the landscape evolution described in Lebrun et al. (2018) 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 Guérin et al. (2020) 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 nondimensional number, the so-called channelization index. •We perform direct numerical simulations of System (2) for realistic experimental data and on timescales associated with the erosion phenomenon. We tackle the problem of severe CFL restrictions due to very short timescales associated with 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 Ha Hilbert space L2 T(H):= L2((0,T);H);CT(H):= C((0,T);H) the space of functions u:t→ u(t,.) ∈Hthat are, respectively, L2integrable or continuous in time on the interval (0,T). We prove the following result: Theorem 1.1 Let m >0,n >3or n =2,K >0, and T0>0. Let us fix two constants fluid heights href >hmin >0. Suppose that the initial data h0,z 0,c 0satisfy h0−href ,z0∈Hk+1(R2), c0∈Hk,h0(x)≥2hmin ∀x∈R2. with k =3. Suppose that r ∈L2 T0(Hk),Kh min −||h0||2 L∞≥0. Then there exists 0<T<T0such that System (2) admits a unique solution (h,z,c)with h−href ,z∈L2 T(Hk+2)∩CT(Hk+1), c∈CT(Hk). 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 Chen et al. (2014a). 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 20 Page 6 of 55 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 zdependence, the equation on hreads ∂th=hh+|∇h|2. In Brändle and Vázquez (2005), the authors prove the well-posedness for this equation within the class of bounded, continuous and nonnegative 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 Smith and Bretherton (1972) and Loewenherz (1991) analyze 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 linearized 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, it is not the same framework as in the model (2). Note that in a couple of recent papers Anand et al. (2020) and Bonetti et al. (2020), the following system of 2 PDEs was considered: ∂tz=Kz−e(h H)m|∇z|n+U,∂ th=div(hv0∇z |∇z|)+R,(3) where R,Uand V0are constants. A numerical scheme is designed in Anand et al. (2020)for(3) 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 CI=em+n KnU1−n drives the formation of channels: The number of channels and their branching increase with CI. This analysis is completed by a spectral stability analysis of a spatially nonhomogeneous steady state where the topography is a hillslope which is divided in the middle. It is found numerically that there exists a critical C0 Isuch that the steady state is stable if CI≤C0 Iand unstable otherwise. 123 Page 7 of 55 20 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 Kplays an important role in this study. We show that there is a critical value ¯ Ksuch 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 Bonetti et al. (2020). 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 Guérin et al. (2020), 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 m per second. Parameters chosen in the numerical simulations are based on these experiments. As in Anand et al. (2020); Bonetti et al. (2020), we have observed that decreasing K(respectively, increasing the channelization index CI) reinforces the channelization process. Note that we have focused here on the formation of channels: Unlike simulations made in Lebrun et al. (2018) 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 organized as follows: First, in Sect. 2, we describe the model and the associated system of equations. Then, Sect.3is devoted to the proof of the wellposedness character of the system in short time, (see Theorem 1.1). Next, in Sect.4we 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 Sect.4.2. Finally, Sect. 6draws a brief conclusion of the paper and provides some future perspectives. 2 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 zand the sediment concentration c. Evolution of topography. The evolution of the bottom topography zis given by: ∂tz=Kz−E+S. 20 Page 8 of 55 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)∈R2. The parameter K>0 is a constant, and the term Kzmodels 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 with 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 Howard and Kerby (1983): E(t,x,y)=αQμvν. Consequently, as Q=hv, we will suppose that the erosion speed depends on the power of the norm of the water velocity and on the 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 eis the erosion speed in the conditions h=Hand |v|=Vwith 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)=sc(t,x,y) csat , with sthe speed of sedimentation in the reference condition c=csat. Therefore, the soil elevation evolves according to the equation: ∂tz=Kz−eh(t,x,y) Hm|v(t,x,y)| Vn +sc(t,x,y) csat .(4) Page 9 of 55 20 The landscape evolution model. We complete Equation (4) with two evolution equations for the fluid height hand sediment concentration c. The mass conservation law for the fluid reads ∂th+div(hv) =r,(5) where ris 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 ρsis the volumetric mass density of the sediments and is constant. In order to close System (4), (5), (6), 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 Sect.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 Guérin et al. (2020). The domain ⊂R2has length Lxand width Ly:=[0,Lx]×[0,Ly]. Denote θthe inclination of the plane. We can decompose the bottom topography zas z(t,x,y)=(Lx− x)tan θ+˜z(t,x,y)where ˜zis 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), (7) admits the new form: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th+μtan θ∂xh=μdiv(h∇(h+z)), (8a) h∂tc+μhtan θ∂xc=μh∇(h+z).∇c+ρseh(t,x) Hm|v(t,x)| Vn −ρssc(t,x) csat ,(8b) ∂tz=Kz−eh(t,x) Hm|v(t,x)| Vn +sc(t,x) csat .(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), (7), locally in time. 20 Page 16 of 55 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 Cherrier and Milani 2012 for more details): Proposition 3.7 Assume that h0−href ,z0∈Hk+1(R2),r,sc i−ehm i|vi|n∈L2 T(Hk), hi∈L∞([0,T]×R2)with hi≥hmin and ∇hi,∇zi∈L1 T(Hk). Then there exists a unique solution (hi+1,zi+1)to the equations (25a), (25band hi+1−href ,zi+1∈ L2 T(Hk+2)∩CT(Hk+1(R2)). For the transport equation (25c), we use theorem 7.2.2inMétivier(2008) to state the result: Proposition 3.8 Assume that c0∈Hk(R2),∇vi∈L∞ T(Hk−1),ehm−1 i|vi|n−sci/hi− rci/hi∈L1 T(Hk)and vi∈L1 T(W1,∞). Then there exists a unique solution ci∈ CT(Hk)to the equation (25c). We prove by induction that for all i∈N, System (25a) is well posed. First, for i=0, the functions of (hi,zi,zi)are time independent and h0≥hmin. Thus, as T0=Tis finite, we have sc0−ehm 0|v0|n∈L2 T0(Hk),h0∈L∞ T0(R2)and ∇h0,∇z0∈L1 T0(Hk). Then there exists a unique solution (h1,z1)∈L2 T1(Hk+2)∩CT1(Hk+1)of (25a) and (25b) with 0 <T1≤T0such that h1≥hmin. Similarly, there exists c1∈CT1(Hk) solution of (25c). Now, if we assume that hi−href ,zi∈L2 T(Hk+2)∩CTi(Hk+1)and ci∈CTi(Hk)then, one has hi∈href +CTi(H2)⊂L∞([0,T]×R2). Moreover, we have ∇hi,∇zi∈CTi(Hk)⊂L1 Ti(Hk)and ci∈CTi(Hk)⊂L2 Ti(Hk). Finally, one has hm i|vi|nHk≤Chm−3 iL∞hiL∞+∇hiHk+∇ziHkn+3∈L∞ Ti(Hk)⊂L2 Ti(Hk). Thus, the assumptions of Proposition 3.7 are satisfied and there exists a unique solution hi+1,zi+1of (25a,25b) such that hi+1−href ,zi+1∈L2 Ti+1(Hk+2)∩ CTi+1(Hk+1(R2)) with 0 <Ti+1≤Tisuch that hi+1≥hmin. The existence of a solution ci+1∈CTi+1(Hk)follows similarly. Proposition 3.9 (Uniform bounds) Denote Ei(t)=hi−href L2+∇hiHk+ziHk+1+ciHk. There exists B >0and 0<T∗≤T independent of i ∈Nsuch that, for all i ∈N sup t∈[0,T∗] Ei(t)≤B,and hi(t,x)≥hmin,∀x∈R2. Moreover, T∗ 0∇hi2 Hk+1+∇zi2 Hk+1≤C(T∗)B,∀i∈N∗. Page 17 of 55 20 Proof We proceed by induction. One has easily E0(t)=E(0)≤Band h0≥hmin since h0,z0,c0are 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 E1(t)+ε∇h12 Hk+1+∇z12 Hk+1≤C(ε)(r2 Hk+E0(t)α)E1(t), (26) for some constant C(ε) depending only on εand h0. By integrating Equation (26) with respect to time, one finds: E1(t)+εt 0∇h12 Hk+1+∇z12 Hk+1 ≤E(0)+C(ε) t 0r2 Hk+t 0 C(ε)E0(s)αE1(s)ds. Denote B=2E(0)+C(ε) T 0r2 Hk. We choose ˜ T∗such that e˜ T∗C(ε)Bα≤2. Then, by applying Gronwall’s lemma, one obtains: E1(t)≤B,∀t∈[0,˜ T∗], ˜ T∗ 0∇h12 Hk+1+∇z12 Hk+1≤C(ε)1+ln(2) εB. Next, we have h1=h0+t 0 div(h0∇(h1)+h1∇z0)) +r. Thus, h1−h0L∞≤C(˜ Bt +√Bt)for some constant Cindependent of the problem and related to Sobolev injections, whereas ˜ B=B+href √B. Then, there exists 0<T∗≤˜ T∗such that C(˜ Bt +√Bt)≤hmin and we deduce that h1≥h0−hmin ≥hmin ∀t∈[0,T∗]. This proves the initial step for i=1. Now assume that Ei(t)≤Band hi≥hmin for all t∈[0,T∗]. The estimates on hi+1,zi+1and ci+1are a direct consequence of the 20 Page 18 of 55 energy estimate: d dt Ei+1(t)+ε∇hi+12 Hk+1+∇zi+12 Hk+1≤C(ε)(r2 Hk+Ei(t)α)Ei+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 (∂thi)et (∂tzi)are uniformly bounded in L2 T∗(Hk). The sequence (∂tci)is uniformly bounded on L2 T∗(Hk−1). Proof For p∈{0,...,k}, we differentiate ptimes the equations (25a), (25b)(25b) and multiply it by ∇p∂thi,∇p∂tziand ∇p∂tci, respectively, and integrate in space. Using Proposition 3.9, we obtain bounds for ∂thiL2 T∗(Hk),∂tziL2 T∗(Hk)and ∂tciL2 T∗(Hk−1).  3.5 Convergence of the Sequences In what follows, we denote T∗=Tin order to simplify the notations. Proposition 3.11 The sequences (hi)i∈N,(∂thi)i∈Nand (zi)i∈N,(∂tzi)i∈Nare Cauchy sequences in L2 T(L2). The sequences (ci)i∈Nand (∂tci)i∈Nare Cauchy sequences in CT(L2). Proof In the following, we denote the quantities of the form fi+1−fiby δfi. For all i≥1, the equations for δhi,δzi,δciare written as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂tδhi−div(hi∇δhi+δhi−1∇hi)−div(hi+1∇δzi+δhi∇zi−1)=0,(28a) ∂tδzi−Kδzi=sδci−ehm i(|vi|n−|vi−1|n)+|vi−1|n(hm i−hm i−1),(28b) ∂tδci+δvi−1∇ci+1+vi−1δ∇ci =ehm−1 i(|vi|n−|vi−1|n)+|vi−1|n(hm−1 i−hm−1 i−1) −(s+r)ci hi−ci−1 hi−1. (28c) Bounds on δhi,δziand δci: We multiply the equation (28a)byδhiand integrate over space: 1dδhi2 L2+R2(hi∇hi+1−hi−1∇hi)∇δhi+R2(∇zihi+1−∇zi−1hi)∇δhi=0. Page 19 of 55 20 First, on the one hand, we have R2 (hi∇hi+1−hi−1∇hi)∇δhi=R2 hi(∇δhi)2+δhi−1∇hi∇δhi ≥hmin∇δhi2 L2+R2 δhi−1∇hi∇δhi. On the other hand, we have the estimate: R2 (∇zihi+1−∇zi−1hi)∇δhi=1 2R2∇zi∇(δhi)2+R2 δzi−1hi∇δhi =−1 2R2 zi(δhi)2+R2 δzi−1hi∇δhi ≤Cδhi2 L2+hiL∞δzi−1L2∇δhiL2 ≤Cδhi2 L2+1 2hmin δzi−12 L2+hmin 2∇δhi2 L2. Consequently, we obtain: d dtδhi2 L2+hmin∇δhi2 L2≤Cδhi2 L2+Cδhi−12 L2+Cδzi−12 L2.(29) We proceed similarly for δzi. By using Equation (28b), one finds: d dtδzi2 L2+2K∇δzi2≤s(δci2 L2+δzi2 L2)+2eR2(hm i|vi|n−hm i−1|vi−1|n)δzi. In order to bound the right-hand term, we use the inequality ||x|n−|y|n|≤n|x− y|max(|x|,|y|)n−1: hm i|vi|n−hm i−1|vi−1|n=hm i(|vi|n−|vi−1|n)+|vi−1|n(hm i−hm i−1) ≤nhm i|vi−vi−1|max(vi−1L∞,viL∞)n−1 +m hmin |vn i−1|hi−hi−1|max(hiL∞,hi−1L∞)m. Consequently, as hi,hi−1and vi,v i−1are uniformly bounded in L∞([0,T]×R2): R2 (hm i|vi|n−hm i−1|vi−1|n)δzi≤C(vi−vi−1L2+hi−hi−1L2)δziL2 ≤C(∇δhi−1L2+∇δzi−1L2+δhi−1L2)δziL2 ≤C(∇δhi−12 L2+∇δzi−12 L2)+δhi−12 L2+δzi2 L2). Thus, 20 Page 20 of 55 d dtδzi2 L2+2K∇δzi2 L2≤Cδzi2 L2+C(∇δhi−12 L2+∇δzi−12 L2+δhi−12 L2). (30) Finally, for δci, we use the same method as before with the equation (28c) and we obtain: 1 2 d dtδci2 L2=−(δvi−1∇ci+1+vi∇δci)δci +e(hm i|vi|n−hm i−1|vi−1|n)δci−(s+r)ci hi−ci−1 hi−1δci =−δvi−1∇ci+1δci−1 2div(vi)(δci)2 +e(hm i|vi|n−hm i−1|vi−1|n)δci−(s+r)ci hi−ci−1 hi−1δci ≤δvi−1L2∇ci+1L∞δciL2+∇viL∞δci2 L2 +C(∇δhi−12 L2+∇δzi−12 L2+δhi−12 L2) +s+rH2 hmin ciL2+ci−1L2δciL2.(31) Then we add the inequalities (29), (30) and (31) and integrate on [0,t]with 0 ≤ t≤T: δhi2 L2+δzi2 L2+δci2 L2+hmin t 0∇δhi2 L2dt+2Kt 0∇δzi2 L2dt ≤δh02 L2+δz02 L2+δc02 L2+Ct 0δhi2 L2+δzi2 L2+δci2 L2dt +Ct 0δhi−12 L2+δzi−12 L2+δci−12 L2dt +Ct 0 (∇δhi−12 L2+∇δzi−12 L2)dt. Page 21 of 55 20 We apply the Gronwall lemma and obtain, for all t∈[0,T], δhi2 L2+δzi2 L2+δci2 L2+hmin t 0∇δhi2 L2dt +2Kt 0∇δzi2 L2dt≤Cδh02 L2+δz02 L2+δc02 L2 +t 0δhi−12 L2+δzi−12 L2+δci−12 L2dt +t 0 (∇δhi−12 L2+∇δzi−12 L2)dteCt. With this inequality, we deduce by induction on ithat ∀i∈N, δhi2 L2 T(L2)+δzi2 L2 T(L2)+δci2 L∞ T(L2) ≤CTe CTi i!δh02 H1+δz02 H1+δc02 L2. Consequently, the series δhi2 L2 T(L2),δzi2 L2 T(L2)and δci2 L∞ T(L2)converge; thus, the sequences (hi)i∈N,(zi)i∈N,(ci)i∈Nare Cauchy sequences in the required spaces. Bounds on ∂tδhi,∂tδziand ∂tδci: Like the estimates in Proposition 3.10,weuse the system (28a) and Proposition 3.9 to obtain the bounds.  By Proposition 3.11, there exists h−href ,z∈L2 T(L2)such that hi−href converges to h−href and ziconverges to zin L2 T(L2).As(hi−href )and (zi)are uniformly bounded in L2 T(Hk+2), we obtain by interpolation that ∀1/2<θ<1, ∀i>j∈N, T 0hi−hj2 Hθ(k+2)dt≤T 0hi−hj1−θ L2hi−hjθ Hk+2dt ≤hi−hj1−θ CT(L2)T 0hi−hjθ Hk+2dt ≤hi−hj1−θ CT(L2)T2/(2−θ)hi−hjθ L2 T(Hk+2) ≤hi−hj1−θ CT(L2)T2/(2−θ)(2C1)θ−→ n→+∞ 0. Therefore, ∀k/2<s<k, the sequences (hi−href )and (zi)are Cauchy sequences in L2 T(Hs+2); thus, h−href ,z∈L2 T(Hs+2). Moreover, ∂thiconverges to ∂thin L2 T(L2) and (∂thi)is uniformly bounded in L2 T(Hk), similarly for ∂tz. Consequently, ∀s<k, ∂th,∂tz∈L2 T(Hs); thus, by Proposition 3.5,h−href ,z∈CT(Hs+1). Finally the a priori estimates on hand zallow to conclude that h−href ,z∈L2 T(Hk+2)∩CT(Hk+1). In particular, as k+1=4, (hi−href ),(∇hi),(∇2hi),(∂thi)and (zi),(∇zi),(∇2hi), (∂tzi)converges in C([0,T]×R2). 20 Page 22 of 55 Now we consider (ci)i∈N. By Proposition 3.11, there exists c∈CT(L2)limit of (ci)i∈Nin this space. We know that (ci)i∈Nis uniformly bounded in CT(Hk),soby interpolation: ∀0<θ<1, ∀i>j∈N, sup t∈[0,T]ci−cjHθ(k)≤sup t∈[0,T]ci−cj1−θ L2sup t∈[0,T]ci−cjθ Hk ≤ci−cj1−θ CT(L2)(ciCT(Hk)+cjCT(Hk))θ ≤Cci−cj1−θ CT(L2)−→ n,m→+∞ 0. Therefore, (ci)i∈Nis a Cauchy sequence and thus converges to cin the space CT(Hs), for all s<k. Moreover, (∂tci)converges to ∂tcin CT(L2)and is uniformly bounded in CT(Hk−1), so it converges to cin CT(Hs−1). And we conclude by the a priori estimate on cthat c∈CT(Hk). To conclude, (ci)i∈N,(∂tci)i∈Nand (∇ci)i∈Nconverge in CT(R2); thus, we can take the limit in the equations, and (h,z,c)is solutions of System (9a). This concludes the proof of the existence. 3.6 Uniqueness Proposition 3.12 Let (h1,z1,c1)and (h2,z2,c2)be two solutions of System (9a)satisfying the hypotheses of Theorem 1.1. Then ∀t∈[0,T]: h1(t)−h2(t)2 L2+z1(t)−z2(t)2 L2+c1(t)−c2(t)2 L2 ≤h0 1−h0 22 L2+z0 1−z0 22 L2+c0 1−c0 22 L2eCT 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:= h1−h2,δz:= z1−z1and δc:= c1−c2: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂tδh−div(h1∇δh+δh∇h2)−div(h1∇δz+δh∇z2)=0 (32a) ∂tδz−Kδz=sδc−ehm 1(|v1|n−|v2|n)+|v2|n(hm 1−hm 2)(32b) ∂tδc+δv∇c2+v2δ∇c=ehm−1 1(|v1|n−|v2|n)+|v2|n(hm−1 1−hm−1 2) −(s+r)c1 h1−c2 h2(32c) Then the bound is obtained with a similar process as in the proof of Proposition 3.11. 123 Page 23 of 55 20 4 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 (8a) 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. 4.1 Non-dimensionalization and Linearization of the System We write System (8a) in a non-dimensional form in order to identify the important parameters. We introduce several characteristic quantities: Zis a characteristic eroded height, Ha characteristic water height, La characteristic wavelength and Ta characteristic time. We chose Tto be the necessary time to erode the soil of a height Z, with an erosion speed e. Thus, Tverifies eT =Z.Asv=μtan θe1−μ∇(z+h),wefix the characteristic water velocity V=μ. We introduce the dimensionless variables: h:= h H,z=z Z,v := v V,c:= c csat x:= x L,t:= e Zt. In order to simplify the notations, we will assume H=Z. Then, dropping the primes, System (8a) is written as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th+ZV tan θ eL ∂xh=Z2V eL2div(h∇(h+z)), h∂tc+ZV tan θ eL h∂xc=ZV eL2h∇(h+z).∇c+ρs csat hm|tan θe1−Z L∇(h+z)|n−s e ρs csat c, ∂tz=ZK eL2z−hm|tan θe1−Z L∇(h+z)|n+s ec. To simplify the equations, we set Z L=e V, and we define α:= Z L=e V,K:= ZK eL2= K LV . The system reads: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∂th+tan θ∂xh=αdiv(h∇(h+z)), h∂tc+tan θh∂xc=αh∇(h+z).∇c+ρs csat hm|tan θe1−α∇(h+z)|n−ρss csatec, ∂tz=Kz−hm|tan θe1−α∇(h+z))|n+s ec. (33) 20 Page 24 of 55 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 shmtannθ>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: ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th+tan θ∂xh=αh(h+z) ∂tc+tan θ∂xc=ρs csat s emc h2h−αnc htan θ∂x(h+z)−c h ∂tz=Kz−s emc hh−αnc tan θ∂x(h+z)−c(34) We denote f=(h,c,z)T. Then fverifies the equation ∂tf=A0f+A1∂xf+A2f with: A0=⎡ ⎢ ⎢ ⎢ ⎣ 000 amc h−a0 −amc ¯ρs ah ¯ρs 0 ⎤ ⎥ ⎥ ⎥ ⎦,A1=⎡ ⎢ ⎢ ⎢ ⎣ −tan θ00 −αanc tan θ−tan θ−αanc tan θ αanh c ¯ρstan θ0αanh c ¯ρstan θ ⎤ ⎥ ⎥ ⎥ ⎦,A2=⎡ ⎣ αh0αh 000 00K⎤ ⎦, where we have denoted a=sρs ehcsat and ¯ρs=ρs csat . We apply the Fourier transform in space and the equation verified by ˆ f, the Fourier transform in space of f,is: ∂tˆ f=A0+iξA1−(ξ2+η2)A2ˆ 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 Vi(ξ, η). The expressions for λiare not explicit: In the next subsection, we compute numerically the stability of the system. Then Sect. 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. 123 Page 25 of 55 20 Fig. 2 Stability diagrams in the plane (ξ, η) ⊂R2,forK>0. The system is unstable in the red area and stable in the green area (Color figure online) 4.2 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 behavior of the system in the top right quarter of the plane. The choice of parameters for these computations is the same as in Sect.5(unless otherwise specified), see Table 1 for their values. In Fig.2, we have presented an illustration of the stability of the system when K>0, where Ke=5×10−4 3600 m2s−1. The domain represented is a bounded subset of the plane (ξ, η) ⊂R2, and the color represents 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 Kis higher, the stable area is larger. It seems that the unstable area is bounded. This will be confirmed by Proposition 4.2. Then, in Fig.3there is no creep effect: K=0, and we represent various stability diagrams for several values of the ratio n/m.Thevalueofmis the same as in Table 1: m=1.6, and nvaries between m/2 and 10m. We observe that the behavior of the system changes with the ratio m/n, but the unstable area always seems to be unbounded. When nmas in Fig.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,asinFig.3c, the system is stable only for longitudinal perturbations. Finally, when n<mthe system seems to be unstable at all frequencies, as in Fig.3d. 20 Page 32 of 55 Table 1 Parameters of numerical simulations Length of the domain Lx40 cm Width of the domain Ly10 cm Characteristic water speed V=μtan θ1m/s Initial water height h00.5mm Initial sediment concentration c0317 g/m3 Exponent of friction over hm1,6 Exponent of friction over vn3,2 Density of the sediments ρs2.17 ×106g/m3 Concentration of saturation csat 3.17 ×105g/m3 Erosion speed e0.5 mm/hour Sedimentation speed se/2000 Angle of the plane θ39◦ 5 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 Lx,Ly, V,h,ρs,csat,eand θare chosen according to the laboratory experiment (Guérin et al. 2020), which erodes a block of salt. The choice of values for the exponents mand nhas been investigated many times in the literature. The values are chosen between 0 and 3, with an additional relation n=2m, as explained in Chen et al. (2014b) (here the constant ncorresponds to the constant m+nin the literature). We choose n=2m, with nsufficiently 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/eis 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 timescale is one hour and a typical time step would be a minute. However, the numerical time step δtis driven by a CFL: δt≤δx/vwater ≈0.001sif one consider a typical mesh size δx=1mm (for a channel of length 400mm). 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 timescale. A comparison between the two resolution methods is made in Appendix D, and it shows that taking the stationary equation does not affect the results. 123 Page 33 of 55 20 Fig. 7 Perturbation of the initial surface, on the domain [0,Lx]×[0,Ly].The represented quantity is the difference between the flat surface and the perturbed soil, in millimeters. This quantity is zero at yellow points and increases when color is darker (Color figure online) The stationary equations for the water height and concentration of sediments in water are discretized with a finite volume method. The scheme is given in Appendix D. Equation (8b)oncis a linear equation, and we discretize 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: ∂xc+vy vx ∂yc=ρs h(E−S). (37) Equation (8a)onhis nonlinear; thus, it is harder to discretize it. In order to avoid an implicit discretization of this equation, we made the choice to linearize the equation. Denoting by hnthe solution of the equation at time tn, we approximate: div(h∇h)(tn)+div(h∇z)(tn)≈div(hn−1∇hn)+div(hn∇zn). As the solution at the previous time step is known, the right-hand term is linear in hn. Thus, we can discretize 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/Lywhere Lyis the width of the domain. The boundary conditions at the top of the domain are Dirichlet condition for h,cand 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 hand c. In the simulations, the initial surface is a flat tilted plane with a small random perturbation. This surface is presented in Fig.7, and it shows a flat view of the twodimensional plane. The color scale shows the height difference between the actual soil and the flat plane. Thus, the yellow areas are the less dug parts. 20 Page 34 of 55 Journal of Nonlinear Science (2024) 34 :20 Fig. 8 Eroded height of the soil in mm. The darkest areas are the most eroded areas Fig. 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 channels is larger for larger values of K 5.1 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 Fig.8we presents results of the simulation of the system, when K=Ke=5×10−4 3600 m2s−1. In this case, the system is spectrally stable at all frequencies {(ξ, ηn);ξ∈R,n∈ N}. The pictures represent the eroded height z(T)−z0at time T=0.25 hour and T=2 hours. We observe that surface perturbations are quickly smoothed, and tend to disappear. After 15 minutes, we see in Fig. 8a that the amplitude of perturbation has not decreased yet, but these perturbations are smoother than initially. After 2 hours, we see in Fig.8b that the amplitude of the perturbations has decreased. In Figs.9and 10, we present, respectively, the eroded surface and the fluid height when the system is unstable at some frequencies. In Fig.9, we can observe the formation of channels in the soil, in the flow direction. The width of these channels is larger when Kis higher, and they take more time to appear. This can be explained by the 123 Page 35 of 55 20 Fig. 10 Water height on the domain, in millimeters. The water depth is greater in darker areas Fig. 11 Discrete frequency analysis of the result of the simulation, when K=Ke/20 (Fig. 9a) and K=Ke/50 (Fig.9b). The quantity represented is the norm of the discrete Fourier transform of the signal in the transverse direction stability analysis, as discussed below. Figure10 has some similarities with Fig.9since the water tends to fill the eroded channels; the water depth is larger in the channels and smaller between them. 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.InFig.11, we present 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=Lx). The frequencies are the 2πN/Ly,forN∈{0,...,Ny/2}. When K=Ke/20, we can see in Fig. 11 that the dominant frequency is reached when N=2 at the bottom of the tilted plane. When K=Ke/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 Fig.12,we 20 Page 36 of 55 Fig. 12 Maximum of the real part of the three eigenvalues, normalized 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 (Color figure online) present numerical computations of the eigenvalues of the linearized 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/Ly, where Lyis the width of the domain. In Fig.12, the white lines shows these frequencies. For both cases (Fig.12a, b), the system is unstable and the maximum of instability is reached at ξ=0, η1=2π/Ly. Thus, the system destabilizes in the transverse direction, and the most destabilizing 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π/Ly). 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 nonlinear effects of the model that are not taken in account in the stability analysis. 5.2 Simulations Including a Source Term for Water In this part, we present some results of simulations when the source term ris a positive function, to observe the effect of rain in the model. In the following simulations, K=Ke/20. First, we consider a constant in time and uniform in space source term r1= 0.005mm/s. Figure13 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=Lx), than at the top (x=0). As the erosion rate is proportional to the power of the water height h,the surface erodes faster at the bottom of the domain. Then, Fig.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. 123 Page 37 of 55 20 Fig. 13 Results of the simulation when r=r1, after 2000 iterations Fig. 14 Results of the simulation when r=r1, after 16000 iterations Fig. 15 Results of the simulation when r=r2, after 14000 iterations Next, we choose a bigger source term: r2=0.01m/s, and results of this simulation are shown in Fig. 15. We can see that there are bigger perturbations at the bottom of the domain than with the source term r1. These perturbations are also not totally transverse: they undulate a little bit in the longitudinal direction, which is a different behavior from the case without source term (Fig. 9b). Now we take a constant non-uniform source term: r3(x,y)=0.005 exp −60 ∗(y−Ly/2)2 Ly2. 20 Page 38 of 55 Fig. 16 Results of the simulation when r=r3, after 6000 iterations that correspond to a positive longitudinal band of rain in the middle of the domain which decrease exponentially fast when approaching y=0ory=0.1. At time 0.75 hours, the surface is highly eroded at the bottom of the band (Fig. 16). Here, the evolution of the system is driven by the source term, and 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. 6 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 appear 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 123 Page 39 of 55 20 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 nonconstant 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. Acknowledgements 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. Declarations Conflict of interest No funding was received to assist with the preparation of this manuscript. 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, proved in Morris (1962): Proposition A.1 Let P(X)=X3+(a1+ib1)X2+(a2+ib2)X+a3+ib3be 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=−1a1 0b1>0 4= 1a1a2a3 0b1b2b3 01 a1a2 00 b1b2  >0,−6=−  1a1a2a300 0b1b2b300 01 a1a2a30 00 b1b2b30 00 1 a1a2a3 00 0 b1b2b3  >0 This proposition is directly used to prove Theorem 4.1. 20 Page 40 of 55 Proof of Theorem 4.1 The characteristic polynomial of A, denoted by P,is: P(X)= X+iξtan θ+α2h0α2h −amc h+αiξanc tan θX+a+iξtan θαiξanc tan θ amc ¯ρs−αiξanhc ¯ρstan θ−ah ¯ρs X−αiξanhc ¯ρstan θ+2K  = X+iξtan θ+α2h0α2h −amc h+αiξanc tan θX+a+iξtan θαiξanc tan θ 0hX+iξhtan θ ¯ρs X+2K  . 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/tan2θ,M=amc/h. Consequently, the polynomial writes as: P(Y)= Y+¯20¯2 −M+iξNY+ai ¯ ξN 0¯ hY Y +¯2K−i¯ ξ =−¯ hY i¯ ξNY +¯2M+(Y+a)(Y+¯2)(Y−i¯ ξ+¯2¯ K) =Y3+"a+¯2(1+¯ K)−i¯ ξ(1+N¯ h)#Y2 +"¯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 Yhave 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). Page 41 of 55 20 We denote: ¯a1=1+N¯ h, ¯ b1=a+¯2(1+¯ K), ¯a2=a(1+¯ K)−¯ hM +¯2¯ K, ¯ b2=a+¯2. Thus, Q(λ) =λ3−¯ ξ¯a1+i¯ b1)λ2+−¯2¯a2+i¯ ξ¯ b2λ+a¯2(¯ ξ+i¯2¯ K). The first condition: The determinant −2writes as: −2=−1−¯ ξ¯a1 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 4is given by: 4= −¯ b1¯ ξ¯ b2¯4a¯ K 1−¯ ξ¯a1−¯2¯a2 0−¯ b1¯ ξ¯ b2= −¯ b1¯ ξ(¯ b2−¯a1¯ b1)¯2(¯2a¯ K−¯ b1¯a2) 10 0 0−¯ b1¯ ξ¯ b2 =¯2¯ ξ2 ¯2¯ b2¯a1¯ b1−¯ b2+¯ b1¯a2¯ b1−¯2a¯ K. Let t=¯ ξ2 ¯2∈[0,μtan2θ eh ], then one finds that 4>0∀(ξ, η) ∈(R2)∗if and only if t¯ b2¯a1¯ b1−¯ b2+¯ b1¯a2¯ b1−¯2a¯ K>0∀t∈[0,μtan2θ eh ],∀¯2>0. The terms ¯ b1and ¯ b2are positive, and: ¯a1¯ 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,∀(ξ, η) ∈(R2)∗if and only if ¯a2¯ b1−¯2a¯ K≥0∀¯2>0. We compute: ¯a2¯ b1−¯2a¯ K=a(1+¯ K)−¯ hM +¯2¯ Ka+¯2(1+¯ K)−¯2a¯ K =aa(1+¯ K)−¯ hM+¯2(1+¯ K)a(1+¯ K)−¯ hM +¯2¯ K. 20 Page 48 of 55 We can suppose, up to renormalization that V(ξ, η) =O +∞(1). We calculate a Taylor expansion of λ(ξ, η) ξ2+η2=0+1(ξ, η) +2(ξ, η) +o(1), and of V(ξ, η) =V0+V1(ξ, η) +V2(ξ, η) +o(1 ξ2+η2). Order 0: When ξ2+η2→+∞, System (43) becomes A2V0=−0V0. Its solutions are the eigenvalues and eigenvectors of −A2,so1 0=−K,2 0=−αhand 3 0=0. The associated eigenvectors are X1=⎡ ⎣ αh 0 K−αh⎤ ⎦,X2=⎡ ⎣ 1 0 0⎤ ⎦,X3=⎡ ⎣ 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 1and 3 are bifurcating from zero. Consequently, we will study these two cases. C.1TheCaseK>0 In this part, we continue the asymptotic expansion of the eigenvalue 3. The zero-order terms of the eigenvectors are V1 0=X1,V2 0=X2and V3 0=X3. We calculate the first-order term of 3. The first-order terms of System (43) when ξ2+η2→+∞depends on the asymptotic behavior of iξ.Ifξ→+∞, then the first-order terms of the system are: A2V3 1−iξA1 ξ2+η2V3 0=−3 1V3 0.(44) And if ξis bounded, iξA1and A0are of the same order; thus, the first-order terms are: A2V3 1−iξA1+A0 ξ2+η2V3 0=−3 1V3 0.(45) Journal of Nonlinear Science (2024) 34 :20 Page 49 of 55 20 To compute 3 1, we multiply the equations by the left eigenvector of A2associated with the eigenvalue 0, X∗=X3. Then, after some computations: If ξ→+∞, 3 1=−iξtan θ ξ2+η2,V3 1=1 ξ⎡ ⎣ 0 1 0⎤ ⎦. If ξbounded, 3 1=−iξtan θ ξ2+η2−sρs ehcsat 1 ξ2+η2. When ξis bounded, we already computed the second-order term of 3. When ξ→+∞, the second-order terms of the system are: A2V3 2−iξA1 ξ2+η2V3 1−A0 ξ2+η2V3 0=−3 1V3 1−3 2V3 0. Then, one finds: 3 2=− sρs ehcsat 1 ξ2+η2. Finally, the asymptotic expansion of the eigenvalues of A(ξ, η) when ξ2+η2→ +∞and K>0is: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ λ1(ξ, η) =−K(ξ2+η2)+o(ξ2+η2) λ2(ξ, η) =−αh(ξ2+η2)+o(ξ2+η2) λ3(ξ, η) =−iξtan θ−s¯ρs eh +o(1) C.2TheCaseK=0 When K=0, the zero-order terms of the eigenvalues 1and 3vanish. Thus, in order to prove the second point of Proposition 4.2, we need to compute the next order terms of these two eigenvalues. The eigenvectors associated with the two zero eigenvalues of −A2are X1=⎡ ⎣ αh 0 −αh⎤ ⎦,X3=⎡ ⎣ 0 1 0⎤ ⎦. Thus, there exist constants u1,v 1,u3,v 3such that V1 0=u1X1+v1X3and V3 0= u3X1+v3X3. The left eigenvectors associated with the zeros eigenvalues of −Care: X1 ∗=001 ,X3 ∗=010 . 20 Page 50 of 55 If ξ→+∞, then the first-order terms of System (43)are: A2V3 1−iξA1 ξ2+η2V3 0=−3 1V3 0.(46) If ξ= 0 is fixed, the first-order terms are: A2V3 1−iξA1+A0 ξ2+η2V3 0=−3 1V3 0.(47) To compute 1 1and 3 1, we multiply the equations (46) and (47)byX1 ∗and X3 ∗. Then, when ξ→+∞, one has: ⎧ ⎨ ⎩ 1 1=0,V1 0=X1, 3 1=−iξtan θ ξ2+η2,V3 0=X3. 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: A2V2−iξA1 ξ2+η2V1−A0 ξ2+η2V0=−1V1−2V0. Then, one obtains: 1 2=sc eh m−ξ2 2n, 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=0is: ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ λ1(ξ, η) =sc eh m−ξ2 2n+o(1), λ2(ξ, η) =−αh(ξ2+η2)+o(ξ2+η2), λ3(ξ, η) =−iξtan θ−s¯ρs eh +o(1). Finally, when ξ=0 the characteristic polynomial Pcan be factorized: P(X)=XX2+(a+α2h)X+α2ah 1−mc ¯ρs. Page 51 of 55 20 Thus, it has a zero root, λ1=0, and the two other roots are λ2,3=−α2h−a±(a−α2h)2+4α2amhc ¯ρs 2. Consequently, the asymptotic development of the eigenvalues of the system when η2→+∞and ξ=0, K=0is ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ λ1(0,η)=0, λ2(0,η)=−αhη2+o(η2), λ3(0,η)=s¯ρs eh mc ¯ρs−1+o(1). This achieves the proof of Proposition 4.2. 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 cand hare small. In order to justify this assumption, we rescale the equations by introducing the characteristic variables Zthe eroded height, Lthe characteristic length and T=Z/ethe characteristic time. We define the dimensionless variables: h:= h H,z=z Z,v := v V,c:= c csat , x:= x L,t:= e Zt. Therefore, dropping the primes, System (8a) is written as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th+1 div(hv) =r e, ∂t(ch)+1 div(chv) =ρs csat hm|v|n−s ec, ∂tz=$ Kz−hm|v|n+s ec, v=−μ V∇(h+z). Here, we defined := Le HV .If1, as generally ρscsat, we can neglect the terms ∂thand ∂t(ch)in the equation het ch. Thus, assuming enough regularity on the solutions, we can write div(chv) =hv.∇c+cdiv(hv) =hv.∇c. Moreover, we fix 20 Page 52 of 55 V=μfor simplicity. We obtain: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ div(hv) =$r, hv.∇c=ρs csat hm|v|n−s ec, ∂tz=$ Kz−hm|v|n+s ec, v=−∇(h+z). (48) Discretization of the equations The constraint on hat 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(hn−1∇hn)+div(hn∇zn)=−r,(49) where hnrepresents the fluid height at time tn=nδt, whereas the surface height at time n,zn, has been calculated by an explicit Euler method, using the solutions at time n−1, hn−1and cn−1:∀1≤i≤Nx,1≤j≤Ny, zn i,j−zn−1 i,j dt =K zn−1 i+1,j−2zn−1 i−1,j+zn−1 i,j dx +K zn−1 i,j+1−2zn−1 i,j−1+zn−1 i,j dy −(E−S)n i,j. We discretize the equation on hnby a centered finite volume scheme, which writes: −ri,j=hn−1 i+1/2,j(hn i+1,j−hn i,j)−hn−1 i−1/2,j(hn i,j−hn i−1,j) dx2 +hn−1 i,j+1/2(hn i,j+1−hn i,j)−hn−1 i,j−1/2(hn i,j−hn i,j−1) dy2 +hn i+1/2,j(zn i+1,j−zn i,j)−hn i−1/2,j(zn i,j−zn i−1,j) dx2 +hn i,j+1/2(zn i,j+1−zn i,j)−hn i,j−1/2(zn i,j−zn i,j−1) dy2, where hn−1 i+1/2,j=hn−1 i,j+hn−1 i+1,j 2. Finally, Equation (37)oncis discretized by an upwind finite volume scheme, considering the variable xas a time variable. The time discretization is an explicit Euler scheme. ci+1,j−ci,j dt+max(vx[i,j],0) vy[i,j] ci,j−ci,j−1 dy +min(vx[i,j],0) vy[i,j] ci,j−ci,j+1 dy =− ρs hi,jvx[i,j](E−S)n i,j. Page 53 of 55 20 Fig. 17 Eroded height of the soil in mm, for K=Ke,atT=2s Boundary conditions We choose periodic boundary conditions in the ydirection. 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≤Ny, zn 0,j=zn 1,j,zn Nx,j=zn Nx+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. hn 0,j=h0,hn Nx,j=hn Nx+1,j, cn 0,j=c0,cn Nx,j=cn Nx+1,j. 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 (8a), 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: δt≤δx vwater . Therefore, we use a semi-explicit scheme: Linear terms of the water height and sediment concentration equations are implicit, and we explicit a part of the nonlinear terms. The discretization in time of equations on hand cis given by: ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ hk+1−hk δt+μdiv(hh∇hk+1+hk+1∇zk+1)=0, hk+1ck+1−ck δt+hk+1vk+1.∇ck+1=ehk+1mvk+1n−ρssck+1 csat , vk=−μ∇(hk+zk). (50) 20 Page 54 of 55 Discretization 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 nonstationary System (8a)inFig.17a and the stationary System (48)inFig.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 (8a) by System (48) where water is in a stationary regime is valid. References Adams, R.A., Fournier, J.J.: Sobolev Spaces. Elsevier (2003) Anand, S.K., Hooshyar, M., Porporato, A.: Linear layout of multiple flow-direction networks for landscapeevolution simulations. Environ. Model. Softw. 133, 104804 (2020) Andrew, A.L., Chu, K.-W.E., Lancaster, P.: Derivatives of eigenvalues and eigenvectors of matrix functions. SIAM J. Matrix Anal. Appl. 14(4), 903–926 (1993) Bonetti, S., Hooshyar, M., Camporeale, C., Porporato, A.: Channelization cascade in landscape evolution. Proc. Natl. Acad. Sci. 117(3), 1375–1382 (2020) Brändle, C., Vázquez, J.L.: Viscosity solutions for quasilinear degenerate parabolic equations of porous medium type. Indiana Univ. Math. J. 817–860 (2005) Chen, A., Darbon, J., Buttazzo, G., Santambrogio, F., Morel, J.-M.: On the equations of landscape formation. Interfaces Free Bound. 16(1), 105–136 (2014a) Chen, A., Darbon, J., Morel, J.-M.: Landscape evolution models: a review of their fundamental equations. Geomorphology 219, 68–86 (2014b) Cherrier, P., Milani, A.: Linear and quasi-linear evolution equations in Hilbert spaces. American Mathematical Society Providence (2012) Culling, W.: Soil creep and the development of hillside slopes. J. Geol. 71(2), 127–161 (1963) Culling, W.E.H.: Analytical theory of erosion. J. Geol. 68(3), 336–344 (1960) Davis, W.: The convex profile of bad-land divides. Science 508, 245–245 (1892) Escalante, C., Fernández-Nieto, E., Luna, T.M.D., Narbona-Reina, G.: Modelling of bedload sediment transport for weak and strong regimes. Numerical Simulation in Physics and Engineering: Trends and Applications: Lecture Notes of the XVIII ‘Jacques-Louis Lions’ Spanish-French School, pp. 179–189 (2021) Evans, L.C.: Partial Differential Equations, volume 19. American Mathematical Society (2010) Fernández-Nieto, E.D., de Luna, T.M., Narbona-Reina, G., de Dieu Zabsonré, J.: Formal deduction of the saint-venant–exner model including arbitrarily sloping sediment beds and associated energy. ESAIM Math. Model. Numer. Anal. 51(1), 115–145 (2017) Fernandez-Nieto, E.D., Lucas, C., de Luna, T.M., Cordier, S.: On the influence of the thickness of the sediment moving layer in the definition of the bedload transport formula in exner systems. Comput. Fluids 91, 87–106 (2014) Gilbert, G.K.: Report on the Geology of the Henry Mountains. US Government Printing Office (1877) Gilbert, G.K.: The convexity of hilltops. J. Geol. 17(4), 344–350 (1909) Guérin, A., Derr, J., Du Pont, S.C., Berhanu, M.: Streamwise dissolution patterns created by a flowing water film. Phys. Rev. Lett. 125(19), 194502 (2020) Howard, A.D., Kerby, G.: Channel changes in badlands. Geol. Soc. Am. Bull. 94(6), 739–752 (1983) Lebrun, M., Colom, M., Darbon, J., Morel, J.-M.: Numerical simulation of landscape evolution models. Image Process. On Line 8, 219–250 (2018) Loewenherz, D.S.: Stability and the initiation of channelized surface drainage: a reassessment of the short wavelength limit. J. Geophys. Res. Solid Earth 96(B5), 8453–8464 (1991) Métivier, G.: Para-differential calculus and applications to the Cauchy problem for nonlinear systems (2008) Morris, J.: The routh and routh-hurwitz stability criteria: their derivation by a novel method using comparatively elementary algebra. Aircraft engineering and aerospace technology (1962) Smith, T.R., Bretherton, F.P.: Stability and the conservation of mass in drainage basin evolution. Water Resour. Res. 8(6), 1506–1529 (1972) Page 55 of 55 20