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

# Large structure simulation for landscape evolution models

Julien Coatléven and Benoît Chauveau

IFP Énergies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France

Correspondence: Julien Coatléven (julien.coatleven@ifpen.fr)

Received: 7 April 2023 – Discussion started: 26 June 2023

Revised: 5 June 2024 – Accepted: 1 July 2024 – Published: 13 September 2024

**Abstract.** The aim of this paper is to discuss the efficiency of a new methodology to maintain the accuracy of numerical solutions obtained from our landscape evolution model (LEM). As in every LEM, the tricky part is the coupling between water and sediment flows that drives the nonlinear self-amplification mechanisms. But this coupling is also responsible for the emergence and amplification of numerical errors, as we illustrate here. These numerical instabilities being strongly reminiscent of turbulence-induced instabilities in computational fluid dynamics (CFD), we introduce a “large structure simulation” (LSS) approach for LEM, mimicking the large eddy simulation (LESs) used for turbulent CFD. In practice, this treatment consists in a filtering strategy that controls small-scale perturbations in the solution. We demonstrate the accuracy of the LSS approach in the context of our LEM.

## 1 Introduction

Since the pioneering work of Gilbert in the 19th century (Gilbert, 1880), the meaning of the term “landscape evolution model” (LEM) has evolved until reaching its modern definition in the late 20th century. It is now considered to be a numerical application of a mathematical system that seeks to simulate some of the physical processes controlling the landscape dynamic. The capability of LEMs to provide an integrated simulation in which several processes are addressed makes them particularly relevant to tackle a large variety of contexts. The success of those numerical approaches depends on their ability to correctly handle the positive nonlinear feedback between the water flow, sediment erosion, and deposition in a decent computational time. This nonlinear coupling between water and sediments is indeed expected to potentially induce complex water flow networks even on initially small topographic variations, allowing in return the emergence of complex geomorphic landforms. Some algorithms, in particular the family of multiple flow direction (MFD) algorithms, have long been developed for solving surface water flow models in a low computational time. Until very recently, these solvers were not really linked to any physical model, which ruled out the use of an analytic

solution to compare practical numerical results. It was therefore difficult to decipher if the obtained landform results only from physical processes or from the self-amplification of initially small numerical errors. An alternative definition of the specific catchment area often used to model water flow was proposed in Gallant and Hutchinson (2011) and Bonetti et al. (2018), consisting of solving an abstract uniform flow equation in replacement of using one of the MFD algorithms. Independently and following another path, in Coatléven (2020) a first family of MFD algorithms (those for which water is transferred from cell to cell) has been proven to coincide on Cartesian meshes with a classical discretization of the water mass conservation Gauckler–Manning–Strickler model (GMS). The output of the MFD algorithms is exactly a mesh-dependent mean of the water flux associated with the discrete GMS model. This result explains the mesh and numerical dependency since the output of the MFD does not fulfill the consistency criteria, but it also provides a way to correct it in a post-processing step leading to a consistent discrete approximation of the GMS water flux, extended in Coatléven (2020) to general polygonal meshes. As the GMS model can be seen as a generalization of the model proposed in Gallant and Hutchinson (2011) and Bonetti et al. (2018), this finally closes the loop between MFD algorithms and the
specific catchment area defined in Gallant and Hutchinson (2011) and Bonetti et al. (2018) (more details are given in Sect. 2.1). For those reasons, in the present paper we will use a general GMS model to compute our water flow.

This paper has two objectives: (1) to investigate the conditions for which the geomorphic structures simulated from a landscape evolution model derive from numerical instabilities and (2) to introduce a methodology that improves the accuracy of the numerical solution as well as to discuss its potential importance for LEMs. The landscape evolution model used in this paper considers the GMS model for the surface water flow coupled with a representative erosion and deposition sediment flux model detailed in Sect. 2.2 that has been previously used, for instance, in Granjeon (1996), Eymard et al. (2004, 2005), and Peton et al. (2020) and which is a generalization of the models studied in Smith and Bretherton (1972) and Smith et al. (1997). The linear stability analysis of this model highlights the key parameters that control the self-amplification mechanisms of the various water-sediment flow regimes (see Sect. 2.3). To illustrate the related numerical issues, we test the convergence of numerical solutions towards some prescribed analytic solutions for various water-driven and gravity-driven transport coefficients. Comparison between the analytic and numerical solutions leads us to the conclusion that numerical errors must be treated with the greatest care to avoid any misinterpretation of LEM results: the self-amplification processes at the core of the coupling between water flow and sediment evolution can amplify legitimate numerical round-off or solver errors. Thus, estimating the relative impact of numerical errors on the final geomorphologic structures is challenging, making the use of numerical approaches potentially hazardous, in particular those involving implicit time schemes to discuss and quantify the role of self-amplification mechanisms in realistic geodynamic contexts (e.g., valley formation and spacing; Scheingross et al., 2020; Bonetti et al., 2020; Perron et al., 2009; Hooshyar and Porporato, 2021b).

This self-amplification (“butterfly effect”) is very reminiscent of the numerical issues arising in the field of computational fluid dynamics (CFD) for turbulent flows, which prevents the use of direct numerical simulation for high Reynolds numbers unless high-order methods are used over small space scales and timescales (also notice that unmanaged turbulence can sometimes lead to blow-up problems). This comparison with CFD and turbulent flows is not new and was studied in detail, for instance, in Bonetti et al. (2020) and Hooshyar et al. (2020). The modern solution found by the CFD community to achieve reproducible and meaningful simulations is to replace direct numerical simulation (DNS) of the Navier–Stokes equations by large eddy simulation (LES; Berselli et al., 2005). The objective of LES is to obtain a good approximation of local spatial averages of turbulent flows, recovering the correct dynamics only for the organized structures of the flow (the eddies) which are larger than a certain target length scale $\alpha$. Thus, LES chooses

to abandon the idea of resolving all the scales involved in real physical processes, as there is no hope of using a mesh fine enough to resolve the smallest scales correctly. In practice this is done by filtering the solution to distinguish the flow behavior above and below $\alpha$ and obtaining local averages that are smooth and as mesh-independent as possible. To our knowledge, the first attempt at using an LES approach for simulating landscape evolution albeit without explicitly mentioning LES is Perron et al. (2009), where a Laplacian smoothing (equivalent to a mesh-related box filter in the LES terminology) was applied to the topography. More recently, Hooshyar and Porporato (2021a) and Porporato (2022) have used an average in one direction (which is a limit case of filtering) to obtain robust results on channelization statistics and scaling signatures: in other words, they substitute the elevation and the specific drainage area by their mean values in the axial direction of their rectangular simulated domain. In their conclusion they suggest that the use of LES approaches seems a viable avenue for more complex landscape evolution simulations. In line with this observation, we also believe that the success of the attempts of Perron et al. (2009), Hooshyar and Porporato (2021a), and Porporato (2022), as well the numerous analogies between the instabilities arising in landscape evolution models and turbulence reported in Smith and Bretherton (1972), Scheingross et al. (2020), Bonetti et al. (2020), and Hooshyar and Porporato (2021b), and the numerical experiments strongly advocate for the use of some LES technology to overcome the numerical issues arising in the nonlinear coupling of sediment evolution and water flow. Our main contribution is precisely to develop an LES-type methodology for our LEM. We refer to this method by the acronym LSS for “large structure simulation”. Notice that contrary to Hooshyar and Porporato (2021a) and Porporato (2022) and more in line with what is done in the CFD community, we fix a length scale that corresponds to the size of the smallest structures we want to resolve in the problem, quite independently of the domain size. We also consider a more advanced differential filter, namely the Leray-$\alpha$ filter (Cheskidov et al., 2005; Guermond et al., 2003) that is not related to any specific geometric configuration. In this sense, our work can be considered to be a generalization of Hooshyar and Porporato (2021a) and Porporato (2022). We show that when the filter size is correctly defined the results obtained from the LSS are actually free of the nonphysical heterogeneity.

Obtaining a reproducible result that is as error-free and mesh-independent as possible is, of course, what every modeler expects. On the other hand, the emergence of complex geomorphologic structures, which is an objective sought by many LEM users, can require manually introducing relevant physical heterogeneity after handling numerical errors. Several of our simulations are consequently performed using different types of heterogeneity carried by the initial topography or by other physical parameters, such as a variable roughness index or a variable rain map. The emergence of large geomor-
phic structures is discussed by taking into consideration the understanding gained from this work.

The paper is organized as follows. We begin by introducing the water flow and the sediment flow models of the LEM used to perform the simulations discussed in this paper. We then construct analytic solutions and proceed to a comparison with numerical results in the relevant flow regimes. This leads to the first conclusion that for the studied landscape evolution model and the considered classical implicit finite-volume discretization, without any specific treatment, the obtained numerical solutions are potentially controlled by numerical errors. The second step of this work is to introduce and apply the filtering strategy to the water-sediment equation system. The comparison between numerical and analytic solutions clearly shows the crucial role played by this method. Finally, we illustrate the behavior of our LEM in more complex contexts and we test the impact of variable (in space and time) roughness coefficients and rain maps in the final solution.

## 2 Model and notation

Following Smith and Bretherton (1972), we assume that a sedimentary system can be idealized through the following assumptions: (H1) the basin topography can be represented as a mathematical surface, (H2) the principle of mass conservation applies to this surface, and (H3) the sediment flux at any point of the surface is a function of the local slope and the local discharge of water. In other words, using an Eulerian approach (H1) implies that we consider a fixed geographical region over the time period $[0, T]$ mathematically modeled by means of a domain $\Omega \in \mathbb{R}^2$, a function $b: \Omega \times [0, T] \to \mathbb{R}$ describing the basement, i.e., the lower part of the basin in the $z$ direction, and a function $h_s: \Omega \times [0, T] \to \mathbb{R}$ describing the thickness of the sediments (see Fig. 1). Thus, our basin $B: [0, T] \to \mathbb{R}^3$ can be described for almost every (a.e.) $t \in [0, T]$ by

$$ \mathcal{B}(t) = \left\{ \begin{array}{l} (x, y, z) \in \mathbb{R}^3 | (x, y) \in \Omega \text{ and } b(x, y, t) \le z \\ \le b(x, y, t) + h_s(x, y, t) \end{array} \right\}. \quad (1) $$

The evolution of the basement $b$ is governed by several processes, for instance thermal and structural tectonics. In the present paper we assume that the evolution of $b$ is prescribed as input data, and we focus on computing the evolution of the function $h_s$. For the sake of clarity, we give the expression of the mass conservation (H2) equations, neglecting porosity for simplicity:

$$ \left| \begin{array}{ll} \frac{\partial h_s}{\partial t} + \text{div}(\mathbf{J}_s) = S_s & \text{in } \Omega \times ]t_0, T[, \\ -\mathbf{J}_s \cdot \mathbf{n} = B_s & \text{on } \partial\Omega_N \times ]t_0, T[, \\ h_s = 0 & \text{on } \partial\Omega_D \times ]t_0, T[, \\ h_s(t=t_0) = h_{s,0} & \text{in } \Omega, \end{array} \right. \quad (2) $$

where $\mathbf{J}_s$ is the sediment flux, while $S_s$ and $B_s$ are sediment source terms with $S_s$ representing in situ sediment

<figure><img src="image_2.png" /><figcaption>Figure 1. Representation of the two main surfaces considered in a landscape evolution model in the $(x, y, z)$ parameter space, where $z$ is the elevation and $\Omega$ the spatial domain for $(x, y)$ with boundary $\partial\Omega$. The basement $b$ surface represents the bottom part of the simulated block, on which sediments are deposited. The topographic surface is $b+h_s$, where $h_s$ is the sediment thickness. The simulated sedimentary content is denoted as $\mathcal{B}$.</figcaption></figure>

production (or erosion) and $B_s$ boundary sediment supplies. The domain boundary $\partial\Omega$ is divided between $\partial\Omega_N$ where flux (also called Neumann) boundary conditions are imposed and $\partial\Omega_D$ where we enforce fixed elevation (also called Dirichlet) boundary conditions. In the following the $x-y$ coordinates corresponding to the computational domain $\Omega$ will be expressed in kilometers (km), while sediment height $h_s$ and basement $b$ will be expressed in meters (m). Choosing a model corresponds to choosing a specific expression for the sediment flux and the source terms. A common feature of almost all LEMs is that the sediment flux model $\mathbf{J}_s$ and/or the source term $S_s$ depend nonlinearly on the local discharge of water $Q_w$, very often through a power law like $Q_w^{r_s} ||\nabla(h_s+b)||^{p_s+1}$. Self-amplification mechanisms are known to appear at least for $r_s > 1$ (Smith and Bretherton, 1972; Smith et al., 1997).

### 2.1 The water flow model

Landscape evolution models usually requires defining a “local discharge of water” $Q_w$ to account for water-driven erosion and transport. The practical computation of $Q_w$, when not carefully conducted, is the weak point of many models, causing them to lose any hope of consistency in the mathematical sense of the term. It is one of the main reasons why we observe mesh dependency is some LEMs. We re-
call in this section how to define a physically based local discharge of water that maintains consistency. More details can be found in Coatléven (2020), Gallant and Hutchinson (2011), and Bonetti et al. (2018).

Classically, $Q_w$ is computed directly from the so-called drainage or catchment area, CA (also referred to as the contributing area). For a given outlet of the topography (i.e., the top of $B$), it corresponds to the measure of the projection on $\Omega$ of the part of the topography from which the water contributing to this outlet is coming (Maxwell, 1870; Leopold et al., 1964; Bonetti et al., 2018). Despite being a very intuitive notion, it has evaded a precise mathematical definition for a long time. Classical multiple flow direction (MFD) algorithms are intended to provide a practical way to compute a discrete approximation $CA_\epsilon(K)$ of the catchment area, CA, for a mesh cell $K$ (where $\epsilon$ stands for the mesh precision) and in this way a discrete approximation $Q_K$ of $Q_w$ for cell $K$. As is well-documented (Desmet and Govers, 1996; Pelletier, 2010, 2013; Porporato, 2022) the discrete catchment area $CA_\epsilon(K)$ obtained from those algorithms strongly depends on the cell size, geometry, and orientation with respect to the flow. Several attempts can be found in the literature to reduce this mesh dependency, defining the discrete water flow discharge $Q_K$ associated with a mesh cell $K$ as $Q_K = (CA_\epsilon(K)/w(K))$, where $w(K)$ is a normalization factor related to a geometric property of the cell (see Desmet and Govers, 1996) or to an estimate of the flow width (Pelletier, 2010) defining the so-called (discrete approximation of the) specific or unit catchment area (SCA or UCA). A more modern mathematical definition of the specific catchment area $a$ at the continuous level was proposed in Gallant and Hutchinson (2011) and Bonetti et al. (2018), consisting in solving an abstract uniform flow equation:

$$ \left| \begin{array}{l} -\operatorname{div} \left( a \frac{\nabla(h_s+b)}{\|\nabla(h_s+b)\|} \right) = 1 \quad \text{in } \Omega, \\ -a \frac{\nabla(h_s+b)}{\|\nabla(h_s+b)\|} \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega_{\text{in}}, \end{array} \right. \quad (3) $$

where $\partial\Omega_{\text{in}} = \{x \in \partial\Omega | \nabla(h_s+b) \cdot n > 0\}$ is the part of the boundary that is in going and $n$ denotes the outward normal to $\Omega$. Setting $Q_w = a$ at the continuous level, this leads in practice to computing a consistent discrete approximation $a_K$ of $a$ for a mesh cell $K$. This allows reducing the mesh dependency to the usual consistency errors of numerical schemes. At first sight, the model in Eq. (3) could seem very different from MFD algorithms. However, considering, for instance, the classical cell-to-cell algorithms of Freeman (1989, 1991) and Holmgren (1994), one can see that those algorithms act as if we were distributing a fictitious water flow of a mesh cell to the neighboring cells with lower elevation proportionally to a function of the slope, as illustrated in Fig. 2.

One could then legitimately suspect that those MFD algorithms could be related to a discretization of a water flow model. This has been recently demonstrated in Coatléven (2020) for the most classical cell-to-cell MFD algorithms (for

<figure><img src="image_3.png" /><figcaption>Figure 2. Basic principle of the simplest cell-to-cell MFD algorithm: water is distributed to lower neighboring cells proportionally to the slope (reproduced from Coatléven, 2020).</figcaption></figure>

instance, those of Freeman, 1989, 1991; Holmgren, 1994). It became clear that those MFD algorithms are a way of implementing a solver for the following stationary water mass conservation with Gauckler–Manning–Strickler (GMS) flux modeling of surface runoff:

$$ \left| \begin{array}{l} \mathbf{Q}_w = - \frac{k_m h_w \eta_w(h_w)}{s_{\text{ref}}^{p_w}} \|\nabla(h_s + b)\|^{p_w} \nabla(h_s + b), \\ \operatorname{div}(\mathbf{Q}_w) = S_w \quad \text{in } \Omega, \\ \mathbf{Q}_w \cdot \mathbf{n} = B_w \quad \text{on } \partial\Omega_{\text{in}}, \end{array} \right. \quad (4) $$

where $h_w$ is the water height, $s_{\text{ref}} = 1 \text{ m km}^{-1}$ the reference slope, $p_w$ a model parameter, and $\eta_w$ the water mobility function. For simplicity we assume here that the mobility function has no dimension and is a function of $h_w$ only and that the domain source $S_w$ is given in cubic meters per second per kilometer ($\text{m}^3 \text{s}^{-1} \text{km}^{-2}$) such that its integral over a 2D area measured in square kilometers ($\text{km}^2$) coincides with a discharge in cubic meters per second ($\text{m}^3 \text{s}^{-1}$). The boundary influx $B_w$ is measured in cubic meters per second per kilometer ($\text{m}^3 \text{s}^{-1} \text{km}^{-1}$). The coefficient $k_m$ can be thought of as the Strickler coefficient or the inverse of the Gauckler–Manning coefficient up to a change in unit (strictly speaking, this identification is truly valid for channels and if the mobility function $\eta_w$ is equal to a dimensionless hydraulic radius). For this choice of unit for $S_w$, $k_m$ has the unit of a speed, i.e., meters per second ($\text{m s}^{-1}$). Comparing Eq. (4) with Eq. (3),
we see that Eq. (3) corresponds to the particular case where $k_m = 1$, $p_w = -1$, and $a = h_w\eta_w(h_w)$. In this sense the GMS model in Eq. (4) is a generalization of Eq. (3) that allows including the classical ingredients (nonlinear slope dependency and some spatial heterogeneity) of the family of MFD algorithms.

The analysis of Coatléven (2020) explains how the discrete catchment area $CA(O)$ for the outlet of a region $O$ that is computed by MFD algorithms coincides with an intermediate discrete quantity appearing in the most natural discrete solver for Eq. (4). It also allows giving a continuous interpretation and generalization $CA(O)$ for any region $O$ of the discrete $CA_\epsilon$ that is computed by MFD algorithms only for mesh cells:

$$CA(O) = \int_{\partial O} h_w \eta_w(h_w) (-k_m s_{\text{ref}}^{-p_w} ||\nabla(h_s + b)||^{p_w} \nabla(h_s + b) \cdot \mathbf{n})^+, \quad (5)$$

where $h_w$ is the solution of Eq. (4) with $S_w = 1$ and where we have denoted $v^+$ the positive part of $v$ (i.e., $v^+ = \max(0, v)$). Since the model in Eq. (4) describes a water flow, thanks to Coatléven (2020) we can reinterpret the discrete catchment area $CA_\epsilon(K)$ computed through the classical cell-to-cell MFD algorithms as the total flux leaving cell $K$ of a fictitious water flow with a uniform water source $S_w = 1$. Unfortunately, we also see that even at the continuous level, $CA(O)$ strongly depends on the geometry of $O$ and its orientation with respect to the flow. Since it is detailed in Coatléven (2020) that cell-to-cell MFD computations compute in practice the discrete catchment area $CA_\epsilon(K)$ for each cell $K$ of a mesh through a discretized version of Eq. (5), as a result when MFD algorithms are used to estimate the discrete local discharge of water $Q_K$, it produces cell and thus mesh dependency in the simulated surface water distribution. In line with the attempts of Desmet and Govers (1996) and Pelletier (2010) to define a specific catchment area (SCA) by rescaling the CA, the correct scaling would be to set the normalization factor $w$ to the length of the portion of $\partial O$ along which the fictitious water flow is leaving $O$. A corrected definition of the specific catchment at the continuous level in the spirit of Desmet and Govers (1996) and Pelletier (2010, 2013) area would thus be to use

$$SCA(O) = \frac{1}{\int_{\partial O} \chi_{-k_m s_{\text{ref}}^{-p_w} ||\nabla(h_s+b)||^{p_w} \nabla(h_s+b) \cdot \mathbf{n}} > 0} \int_{\partial O} h_w \eta_w(h_w) (-k_m s_{\text{ref}}^{-p_w} ||\nabla(h_s+b)||^{p_w} \nabla(h_s+b) \cdot \mathbf{n})^+, \quad (6)$$

where $\chi$ is the indicator function (i.e., the function with value 1 when the condition is satisfied and 0 otherwise). Depending on the orientation of the flow, such a normalization will sometimes match the choices of Desmet and Govers (1996) and Pelletier (2010, 2013), explaining their partial success.

This continuous SCA scales as an approximation of the continuous water flux magnitude,

$$q_w = |k_m h_w \eta_w(h_w) |s_{\text{ref}}^{-p_w}||\nabla(h_s + b)||^{p_w+1}, \quad (7)$$

(in $m^3 s^{-1} km^{-1}$) but is not equal to it. The SCA defined by Eq. (6) is in fact a mean of $q_w$ along the outflow portion of $\partial O$ and thus still retains some dependency in the geometry of $O$ and its orientation with respect to the flow. Meanwhile, notice that the specific catchment area $a$ of the model in Eq. (3) can be reinterpreted through Eq. (4) as computing $q_w$ since

$$q_w = |k_m h_w \eta_w(h_w) |s_{\text{ref}}^{-p_w}||\nabla(h_s + b)||^{p_w+1} = |a|||\nabla(h_s + b)||^{-1+1} = a,$$

as we have set $a = h_w\eta_w(h_w) \ge 0$, $p_w = -1$, $k_m = 1$, and $s_{\text{ref}} = 1$ to merge Eq. (3) inside Eq. (4). Thus, in view of the success of Bonetti et al. (2018) and within the context of Eq. (4) it seems very natural to choose to set $Q_w = q_w$. One could consider that the equivalence between classical cell-to-cell MFD algorithms established and the consistency correction proposed in Coatléven (2020) that leads to using a discrete version of $q_w$ is another path to recover the conclusions of Eq. (3) and in this sense that $q_w$ is a generalization of $a$ to more complex water flow models.

The consistency correction proposed in Coatléven (2020) for MFD algorithms precisely coincides with the replacement of the computation of $CA_\epsilon(K)$ or $SCA_\epsilon(K)$ for a mesh cell $K$ by a consistent discrete reconstruction $q_K$ of $q_w$ in each cell $K$. Convergence of this discrete version $q_K$ to $q_w$ when the mesh size goes to zero was proven in Coatléven (2020), along with error estimates. Thus, apart from the usual discretization error no anomalous mesh dependency should remain in $q_K$ in practice, contrary to what is observed for $SCA_\epsilon(K)$ given by MFD algorithms. In this sense, $q_K$ can be seen as a consistency correction for $SCA_\epsilon(K)$, as well as a generalization of Eq. (3) to a richer family of flow models. The interpretation of the local water discharge $Q_w$ as being equal to the water flux magnitude $q_w$ given by Eq. (7) from the solution of Eq. (4) is therefore the default configuration chosen in the water flow model used to perform all the simulations we introduce in this paper.

To say that this model corresponds to the GMS model does not necessarily mean that its scope of application is limited to channels: it depends on the specific choice made on the model parameter values. Steady-state analysis (Graf and Altinakar, 2000; Birnir et al., 2001) for channels suggests using values $\eta_w(h_w) = (h_w/h_{\text{ref}})^{1/2}$ and $p_w = -1/2$, while the classical Gauckler–Manning–Strickler formula would coincide with $\eta_w(h_w) = (R_h(h_w)/h_{\text{ref}})^{2/3}$, with $R_h(h_w)$ the hydraulic radius and again $p_w = -1/2$. When applied to large time- and space scales landscape evolution models, these calibrations are no longer valid and at this stage we suggest considering $\eta_w$ and $p_w$ as modeling parameters that can be
tuned for each considered problem. In the following numerical experiments, since we only consider the water flux $q_w$ the choice of the water mobility function has no influence and we set $\eta_w(h_w) = 1$ for simplicity, as well as $p_w = 0$. Notice that the recommendations deduced from the work discussed in this paper would remain valid for more general choices of those parameters.

The application domain of the GMS model is, however, limited by some additional requirements for the topography $h_s+b$. From a purely mathematical point of view, systems in Eq. (3) and Eq. (4) are in fact stationary transport problems for a or $h_w$. Well-posedness, i.e., existence and uniqueness in a suitable function space and continuity with respect to data, can be rigorously established only under some condition on the topography. Many different conditions are possible, all introducing some positivity requirement in the zero-order part of the differential operators applied to a or $h_w$ (see Coatléven, 2020; Bardos, 1970; Veiga, 1987; DiPerna and Lions, 1989; Fernández-Cara et al., 2002; Girault and Tarbar, 2010). In particular, among the possible conditions the simplest ones are undoubtedly

$$-\Delta(h_s + b) > 0 \quad \text{or} \quad -\operatorname{div}(k_m s_{\text{ref}}^{-p_w} ||\nabla(h_s + b)||^{p_w} \nabla(h_s + b)) > 0. \qquad (8)$$

By enforcing a downflow direction to exist everywhere, they both ensure that the model in Eq. (4) is well-posed, at the price of introducing quite stringent restrictions on the admissible topographies. Notice that they are sufficient conditions and not necessary ones: this implies that solutions to Eqs. (3) and (4) can still exist for some topographies not fulfilling one of the sufficient conditions. In particular, saddle-point or valley-like topographies will not easily fulfill those conditions, while it seems reasonable to assume that a solution will exist in such configurations since water can find a downflow direction. This being said, this mathematical requirement that is probably too strong should act as a warning, as it clearly reveals that not all topographies may be admissible for the model in Eq. (3) and its generalization in Eq. (4). More details on the most stringent requirements (no accumulation or flat areas) are given in Sect. 5.2.

## 2.2 The sediment flux model

In the present paper we have chosen to focus on the stratigraphic model that has already been discussed in detail in Granjeon (1996), Eymard et al. (2004, 2005), and Peton et al. (2020) and which is a generalization of the models studied in Smith and Bretherton (1972) and Smith et al. (1997). The corresponding sediment flux $J_s$ takes the following form:

$$J_s = -\eta_s(h_s) s_{\text{ref}}^{-p_s} ||\nabla(h_s + b)||^{p_s} \left( \left( \frac{q_w}{q_{\text{ref}}} \right)^{r_s} \nabla\psi_w(h_s + b) + \nabla\psi_g(h_s + b) \right) \quad \text{in } \Omega \times ]t_0, T[, \qquad (9)$$

where $r_s > 0$ and $p_s > 0$ are model parameters, $q_w$ is the water flux obtained from Eq. (4), $q_{\text{ref}}$ and $s_{\text{ref}}$ are dimensional factors, and $\eta_s$ is a dimensionless sediment mobility function such that

$$0 \le \eta_s(h_s) \le 1 \quad \text{and} \quad \eta_s(0) = 0, \qquad (10)$$

whose main role is to ensure that the sediment height $h_s$ remains positive. In the following we use

$$\eta_s(h_s) = \begin{cases} 1 - \frac{h_c}{h_s + h_c} & \text{if } h_s \ge 0, \\ 0 & \text{otherwise,} \end{cases} \qquad (11)$$

with $h_c$ a parameter. The function with subscript “w” is intended to model the water-driven processes, while the function with subscript “g” models gravity-related processes. We consider here the most common form for functions $\psi_w$ and $\psi_g$ corresponding to

$$\psi_w(h_s + b) = \int_0^{h_s+b} k_w(v)dv \quad \text{and} \quad \psi_g(h_s + b) = \int_0^{h_s+b} k_g(v)dv, \qquad (12)$$

where $k_w$ and $k_g$ are bounded diffusion coefficients depending solely on $h_s$ in such a way that

$$\begin{aligned} \nabla\psi_w(h_s + b) &= k_w(h_s + b)\nabla(h_s + b) \quad \text{and} \\ \nabla\psi_g(h_s + b) &= k_g(h_s + b)\nabla(h_s + b), \end{aligned} \qquad (13)$$

so the sediment flux follows the topographic slope $\nabla(h_s + b)$.

This sediment flux model is implemented in our modeling platform ArcaDES, and all the simulations shown in the following sections are performed using the ArcaDES platform (although ArcaDES is mentioned for the first time in a scientific paper, it has been used since 2015 in the stratigraphic numerical forward model DionisosFlowTM initially developed by Granjeon, 1996). Both soil erosion and sediment deposition are considered. As ArcaDES is designed for large-scale simulations in time and space, we have chosen to express the x-y coordinates in kilometers (km), time in millions of years (Myr), and sediment height $h_s$ and basement $b$ in meters (m). Thus, the unit for sediment sources will be meters per million years (m Myr$^{-1}$). Since we have chosen to use $Q_w = q_w$ with $q_w$ the water flux from Eq. (4), the unit for the water discharge $q_w$ is cubic meters per second per kilometer ($m^3 s^{-1} km^{-1}$), and thus we naturally set $q_{\text{ref}} = 1 m^3 s^{-1} km^{-1}$. The natural unit of coefficients $k_g$ and $k_w$ is square kilometer per million years ($km^2 Myr^{-1}$), with the reference slope again set to $s_{\text{ref}} = 1 m km^{-1}$.

## 2.3 Some insights from perturbation theory

In this subsection, in order to give a feeling of the potential stability issues related to the model in Eqs. (2)-(9)-(4), we will perform a brief analysis of the behavior of solutions under perturbations. Details of the following computation are
relegated to Appendix A. We assume for simplicity that $k_g$ and $k_w$ are constant functions. Let us denote $(h_{s,*}, h_{w,*})$ a reference solution of Eqs. (2)–(9)–(4) with sources $(S_{s,*}, S_{w,*})$, whose stability is to be tested. We denote $(h_{s,\delta}, h_{w,\delta})$ a perturbation of magnitude $\delta$ of this reference solution associated with the source perturbation $(S_{s,\delta}, S_{w,\delta})$ and consider the evolution of $(h_s, h_w) = (h_{s,*} + h_{s,\delta}, h_{w,*} + h_{w,\delta})$ for the perturbed source $(S_{s,*} + S_{s,\delta}, S_{w,*} + S_{w,\delta})$. Since both the perturbed and unperturbed solutions have to satisfy the boundary conditions, we deduce that the perturbation $(h_{s,\delta}, h_{w,\delta})$ itself also satisfies the same boundary conditions. Then, in line with, for instance, the analysis of Smith et al. (1997), injecting $(h_s, h_w)$ into Eqs. (2)–(9), multiplying by $h_{s,\delta}$, and integrating by parts we obtain the equation governing the evolution of the perturbation’s total energy (see Appendix A for details):

$$ \frac{d}{dt} \left( \frac{1}{2} \int_{\Omega} h_{s,\delta}^2 \right) = - \int_{\Omega} j_s (h_s, h_s + b, q_w) ||\nabla h_{s,\delta}||^2 \\ + \int_{\Omega} S_{s,\delta} h_{s,\delta} + \int_{\Omega} (j_s (h_{s,*}, h_{s,*} + b, q_{w,*}) \\ - j_s (h_s, h_s + b, q_w)) \nabla (h_{s,*} + b) \cdot \nabla h_{s,\delta}, \quad (14) $$

where we have denoted

$$ j_s(u, v, w) = \eta_s(u) s_{\text{ref}}^{-p_s} ||\nabla v||^p \left( \left( \frac{w}{q_{\text{ref}}} \right)^{r_s} k_w + k_g \right). $$

The first term of the right-hand side is always negative and thus always contributes to the stability of the system. The second term describes the contribution of potential perturbation sources $S_{s,\delta}$ (other than the initial conditions) to the evolution of the perturbation’s energy. The last term,

$$ A_{\delta} = \int_{\Omega} (j_s(h_{s,*}, h_{s,*} + b, q_{w,*}) - j_s(h_s, h_s + b, q_w)) \\ \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta}, \quad (15) $$

originates partially from the nonlinearity of the sediment transport model but most importantly from the coupling between the flow and the sediment transport. If $A_{\delta}$ is negative or sufficiently small and if the perturbation source is also small enough, then the sediment perturbation energy will decrease with time. In this case, the solution $(h_{s,*}, h_{w,*})$ is said to be stable under perturbation $(h_{s,\delta}, h_{w,\delta})$. However, the sign of $A_{\delta}$ is not always negative and will often take necessarily small positive values. If $A_{\delta}$ is large enough, instead of being diffused by the first term the sediment perturbation energy will grow with time and potentially become as large as the unperturbed solution: the solution $(h_{s,*}, h_{w,*})$ is then unstable under perturbation $(h_{s,\delta}, h_{w,\delta})$. This is a self-amplification mechanism, as the magnitude of $A_{\delta}$ will grow with the perturbation’s magnitude and cancel if the perturbation is zero, and also because of the dependency of the water

flux $q_{w,\delta}$ on the topography perturbation $h_{s,\delta}$. We will say that growing perturbations correspond to the physically unstable regime.

We can anticipate that the relative magnitude of the gravity and water coefficients $k_g$ and $k_w$ will play a key role in the stability of solutions. Indeed denoting $\tau = (k_w q_w^{r_s}) / (k_g q_{\text{ref}}^{r_s})$, if $k_g$ is much larger than $k_w$ and thus $\tau$ is very small we have, assuming for simplicity that $\eta_s = 1$ (see Appendix A for details),

$$ A_{\delta} \approx -k_g s_{\text{ref}}^{-p_s} (p_s ||\nabla(h_{s,*} + b)||^{p_s-2} |\nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta}|^2) \\ + O(\tau) + O(\delta^3) $$

(where we recall that a function $f$ is $O(h)$, there exists a constant $C > 0$ independent of $h$ such that $||f|| \le Ch$ for a suitable norm $||.||$). Then for large values of $k_g$ the term $A_{\delta}$ is always negative and thus stabilizing. On the contrary, if $k_w$ is much larger than $k_g$ then $\tau$ is also very large and we have (see Appendix A for details)

$$ A_{\delta} \approx -k_w s_{\text{ref}}^{-p_s} r_s \frac{q_{w,*}^{r_s-1}}{q_{\text{ref}}^{r_s}} q_{w,\delta} ||\nabla(h_{s,*} + b)||^{p_s} \nabla(h_{s,*} + b) \\ \cdot \nabla h_{s,\delta} - k_w s_{\text{ref}}^{-p_s} \left( \frac{q_{w,*}}{q_{\text{ref}}} \right)^{r_s} (p_s ||\nabla(h_{s,*} + b)||^{p_s-2} |\nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta}|^2) + O(1/\tau) + O(\delta^3). $$

Regions for which $\nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta} < 0$ will amplify the perturbation proportionally to $k_w$ and the power $r_s - 1$ of the water flux. We also see that the term $A_{\delta}$ will behave quite differently if $r_s > 1$ or $r_s < 1$. Indeed, for $r_s > 1$ the water flux will reinforce the amplification term in a kind of positive feedback loop. On the contrary, for $r_s < 1$ the water flux will temper the amplification term, and thus we can anticipate that it will require much larger values of $\tau$ for instability to occur in this situation. For the general case incorporating $\eta_s$ the behavior is roughly speaking the same, with the main difference that the additional term due to $\eta_s$ can on rare occasions also contribute with the wrong sign (see Appendix A for details). The main conclusion to draw from this brief study is that parameter $\tau$ will be the main criterion governing the appearance of instabilities even for our most general model.

For a subclass of model equations, Eqs. (2)–(9)–(4), with $\eta_s = 1$, $k_g = 0$ and $p_w = -1$ and $p_s = 0$, the stability of solutions have been theoretically studied in Smith and Bretherton (1972) and Smith et al. (1997). This would correspond within our notations to the extreme case where $\tau = +\infty$, for which we expect instability to occur. It was, for instance, established in Smith et al. (1997) that if the reference solution is stationary, the second term is negative only if some specific condition on the gradient $\nabla(h_s + b)$ is satisfied on the boundary of the region of interest, here $\Omega$. The linear stability of analytic stationary solutions that are uniform in one direction has also been considered in Smith et al. (1997). Their conclusion
sion is that under periodic perturbations in the transverse direction, for $r_s \le 1$, the linear stability analysis does not reveal any instability, while for $r_s > 1$, the stationary solutions are linearly unstable if the frequency of the periodic perturbation is large enough. This is coherent with the above brief perturbation study. Notice that the case $p_w = -1$ greatly simplifies such studies: the linear stability analysis can be shown to be equivalent to solving a one-dimensional ordinary differential equation.

The studies mentioned above are focused on the stability of physically meaningful solutions. Here, we want to draw attention to the numerical consequences of this self-amplification phenomenon; in this way we focus on the stability of numerical solutions. Let us explain the key idea: assuming that all functions are regular enough, one could consider (for instance in a finite-difference setting) that our numerical solution is roughly speaking a perturbation of the exact continuous solution, where the source terms $S_{s,\delta}$ and $S_{w,\delta}$ represent the unavoidable consistency and solver errors of our solving process. Then the numerical sediment perturbation energy will satisfy Eq. (14) and will self-amplify in the same way that physical perturbations self-amplify. In the unstable regime, this means that the numerical solution can potentially diverge from the exact solution from a large amount up to the point that it cannot be considered a relevant approximation of the continuous solution, even if the numerical perturbation arises from initially small numerical errors. In other words, in the absence of any treatment, numerical errors may dominate the geomorphological responses of the systems. However, since these numerical errors are reworked by a system of physical equations, it is not impossible to obtain good-looking results. This is where the trap lies for the modeler, who might be inclined to interpret a result induced by uncontrolled numerical noise as a physically plausible solution.

### 3 Numerical instabilities arising from the nonlinear coupling of overland flow and sediment dynamic

To illustrate the numerical issues linked to the self-amplification of initially small numerical errors, we consider in this section several situations where we have either full knowledge of the exact solution or a criterion to distinguish it from incorrect solutions. Thanks to this information on the exact solution, we can illustrate the stability issues of simulations using the model in Eqs. (2)–(9)–(4) (discretized by the finite-volume scheme detailed in Appendix C).

#### 3.1 Instabilities for analytic solutions

In this subsection we consider stationary analytic functions of the form

$$ 
\begin{aligned}
h_s^{\text{ex}}(x, y) &= h_{s,x}(x) + \sum_{p=1}^{N_b} g_b \left( \frac{x-x_p}{\delta_x}, \frac{y-y_p}{\delta_y} \right), \\
h_w^{\text{ex}}(x, y) &= h_{w,x}(x),
\end{aligned}
$$

incorporating $N_b$ small smooth bumps randomly positioned at points $(x_p, y_p)$ chosen such that they do interfere with the boundary conditions, with the smooth bump function given by

$$ g_b(x, y) = g_b(r^2) = \begin{cases} H_{\text{pert}} \exp\left(\frac{-\gamma}{1-r^2}\right) \exp(\gamma) & \text{for } r^2 = x^2 + y^2 \le 1, \\ 0 & \text{otherwise,} \end{cases} $$

with in practice $N_b = 5$, $H_{\text{pert}} = 0.03$ m, $\gamma = 10$, and $\delta_x = \delta_y = 0.25$ km. The numerical domain is rectangular and centered at (0, 0) with the dimensions $Lx = 1$ km on the x axis and $Ly = 5$ km on the y axis, and the basement b is set to zero. We impose homogeneous Dirichlet boundary conditions ($h_s = 0$) on the boundaries $y = -Lx/2$ and $x = Lx/2$ and homogeneous Neumann boundary conditions ($\partial_y h_s = 0$) on the boundaries $y = -Ly/2$ and $y = Ly/2$ as illustrated in Fig. 3.

We use for the mono-dimensional functions ($h_{s,x}, h_{w,x}$) the stationary solution of the model in Eqs. (2)–(9)–(4) in the case $\eta_s = 1$ given in Appendix D that satisfies the boundary conditions. For all our simulations, the constant source terms ($S_{s,x}, S_{w,x}$) for the analytic the stationary solution ($h_{s,x}, h_{w,x}$) in the case $\eta_s = 1$ (see Appendix D for details) are always equal to $(10 \text{ m Myr}^{-1}, 1 \text{ m}^3 \text{ s}^{-1} \text{ km}^{-2})$. Injecting ($h_s, h_w$) into Eqs. (2)–(9)–(4), after some straightforward but tedious computations one can derive exact expressions for the corresponding source terms ($S_s^{\text{ex}}, S_w^{\text{ex}}$), making the pair ($h_s, h_w$) an analytic solution of our model for those source terms.

Given those analytic source terms, initializing the sediment height to the analytic value $h_s(x, y, 0) = h_s^{\text{ex}}(x, y)$ and the water height to the analytic value $h_w(x, y, 0) = h_w^{\text{ex}}(x, y)$ the exact solution of the model in Eqs. (2)–(9)–(4) is of course simply equal to ($h_s^{\text{ex}}, h_w^{\text{ex}}$) for all times. Thus, any reasonable numerical solution should remain a correct approximation of ($h_s^{\text{ex}}, h_w^{\text{ex}}$) for all times. Using the finite-volume discretization described in Appendix C on a Cartesian mesh with square cells for which we denote $\Delta_{xy}$ the size of the edges of the Cartesian cells, we attempt to reproduce the stationary analytic solution by initializing the system to ($h_s(x, y, 0), h_w(x, y, 0) = (h_s^{\text{ex}}(x, y), h_w^{\text{ex}}(x, y))$ and using the analytic source terms ($S_s^{\text{ex}}, S_w^{\text{ex}}$) for various values of the parameters $k_g, k_w, r_s$, and $p_s$. The simulation total time is 0.25 Myr, and we use time steps of maximum length $\Delta t = 0.002$ Myr. The corresponding analytic solutions are
<figure><img src="image_4.png" /><figcaption>Figure 3. Domain configuration for the analytic tests cases.</figcaption></figure>

presented in Figs. 4–8 for the different values of the parameters $k_g$, $k_w$, $r_s$, and $p_s$ we have considered. All those simulations have been performed in parallel on 108 processors using the MPI library. In Fig. 9, we present the obtained convergence curves for all the tested analytic solutions; i.e., we plot the standard $L^2$ error measuring the difference between the simulated sediment height and the exact analytic sediment height. We see in Fig. 9 that for all configurations except the case ($k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$, $k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$), we obtain clean convergences curves, assessing the correctness of our numerical scheme even for the nonlinear couplings. However, for the case ($k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$, $k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$) the two nonlinear couplings ($r_s = 3/2$, $p_s = 1$) and ($r_s = 2$, $p_s = 0$) fail to converge. Looking at Table 1 where we regroup the value of $\tau$ for each test case using our knowledge of the exact solution, we see that convergence problems appear as expected when $\tau$ becomes large. Indeed, since the error increases when we refine the mesh, this error is not a discretization consistency error as all the other test cases validate both our implementation and discretization. On the contrary, it increases with the size of the numerical system, which strongly suggests that it originates from solver (both linear and nonlinear) errors, and this perfectly illustrates the phenomenon of numerical error self-amplification that we have discussed from a theoretical point of view in the Sect. 2.3. Problems are potentially more severe in finer meshes because numerical diffusion that can dissipate residual numerical errors declines with grid spacing. Now, to illustrate how treacherous those numerical solutions are, we present in Figs. 10 and 11 a comparison between the analytic solution and its erroneous numerical counterpart. The erroneous solutions are dangerously “good-looking”: indeed, if only the initial topography and the rain and production data are shown, one could easily be tempted to interpret the quite complex topographies obtained as the realistic self-amplification of the perturbations due to the presence of the bumps. However, since we know the exact solution, we are sure that this is not the case: the appealing numerical solutions are completely wrong. The

overall “geologically realistic” look of the erroneous solution comes from the fact that numerical noise is amplified not by some numerical scheme deficiency but by the capacity of the continuous model to amplify perturbations that we describe in the previous section. In other words, the numerical noise is reworked by the system, giving a “realistic” look to it. This is the reason why we stress that when performing real-life simulations for which of course the correct solution is unknown (otherwise we would not need to simulate anything at all), it can become very hard to decide if the numerical results are correct or blurred by realistic-looking amplified numerical noise. The quality of the numerical scheme, although essential, is not in question: the issue is the self-amplification mechanisms of the continuous model. They are the reason for its physical interest but simultaneously its main issue for performing reliable simulations.

## 3.2 Identifiable instabilities in a non-analytic case

As previously mentioned, having an analytical solution is quite rare when it comes to applying the model in Eq. (2) to realistic cases, and it sometimes becomes difficult to determine whether the numerical solution is correct or not. To illustrate how one can sometimes partially circumvent this difficulty, we consider a simple synthetic topographic surface defined by three constant slope planes. The numerical domain is rectangular with the dimensions $Lx = 400 \text{ km}$ on the x axis and $Ly = 300 \text{ km}$ on the y axis (see Fig. 12a and b). We again use a Cartesian mesh with square cells, the edges of each cell being of length $\Delta_{xy} = 2 \text{ km}$. The gravity diffusion coefficient $k_g$ is equal to $100 \text{ km}^2 \text{ Myr}^{-1}$ in the whole domain, while $k_w = 10 \text{ km}^2 \text{ Myr}^{-1}$ for $h_s + b \ge 0$ and $k_w = 0.1 \text{ km}^2 \text{ Myr}^{-1}$ for $h_s + b < 0$, corresponding to a modulation of the water-induced transport in a fictitious marine domain. Water is supplied by three constant water flux sources located at the domain boundary (black arrows in Fig. 12a), so we call this the “three rivers” test case. Each water source is 12 km large and supplies $1200 \text{ m}^3 \text{ s}^{-1}$ of water.
<table><caption>Table 1. Approximate maximum analytic value of $\tau = \frac{k_w}{k_g} \left( \frac{g_w}{q_{ref}} \right)^{r_s}$ for each convergence test.</caption><thead><tr><th></th><th>(r<sub>s</sub> = 1, p<sub>s</sub> = 0)</th><th>(r<sub>s</sub> = 3/2, p<sub>s</sub> = 1)</th><th>(r<sub>s</sub> = 2, p<sub>s</sub> = 0)</th></tr></thead><tbody><tr><td>(kg, k<sub>w</sub>) = (50, 1) km<sup>2</sup> Myr<sup>-1</sup></td><td>0.01</td><td>0.00353</td><td>0.0025</td></tr><tr><td>(kg, k<sub>w</sub>) = (5, 1) km<sup>2</sup> Myr<sup>-1</sup></td><td>0.1</td><td>0.0353</td><td>0.025</td></tr><tr><td>(kg, k<sub>w</sub>) = (5, 5) km<sup>2</sup> Myr<sup>-1</sup></td><td>0.5</td><td>0.353</td><td>0.25</td></tr><tr><td>(kg, k<sub>w</sub>) = (1, 5) km<sup>2</sup> Myr<sup>-1</sup></td><td>2.5</td><td>1.767</td><td>1.25</td></tr><tr><td>(kg, k<sub>w</sub>) = (1, 50) km<sup>2</sup> Myr<sup>-1</sup></td><td>25</td><td>17.67</td><td>12.5</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4. Sediment height $h_s^{ex}$ of the analytic solution for the case $k_g = 50 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 1 \text{ km}^2 \text{ Myr}^{-1}$.</figcaption></figure>

An essential remark is that the whole configuration is symmetric with respect to the vertical plane $x = Lx/2$. In principle, the equation system consisting of Eqs. (2) and (4), here used with $r_s = 2$, $p_s = 1$, should maintain this symmetry. Since we do not know the exact solution this time, at least we can use symmetry to identify erroneous results that do not fulfill this elementary requirement. Using the finite-volume scheme depicted in Appendix C, we perform a set of three identical simulations in terms of physical parameters but using different numerical settings in order to illustrate the impacts of numerical errors. We perform a sequential computation using GMRES as our linear solver for all systems, its parallel equivalent on four processors, and another sequential simulation using BiCGStab as our linear solver for all systems. The linear solvers are part of the well-known reference PETSc library (Balay et al., 1998) to avoid any potential mistake in their implementation, while the parallelism relies on the Arcane framework (Grospellier and Lelandais, 2009). Final topographies and water flux are shown in the bottom row of Fig. 12. Figure 12c corresponds to sequential GMRES, Fig. 12d to parallel GMRES, and Fig. 12e to sequential BiCGStab.

All the results from these simulations should be almost identical and in any case symmetrical with respect to the vertical plane $x = Lx/2$ in the absence of any spatial heterogeneity in the input data. Clearly, symmetry is lost in the three cases, and what is even more striking is that we get

three very different results. The only difference between the three cases being the numerical solvers, this indicates that it originates from numerical errors. As we are using a decoupled time scheme between water flow and sediment evolution (see Appendix C), one may argue that those instabilities are arising from some violated coupling constraint on the time step. Should this be the case, reducing the time step enough would ultimately lead to clean solutions. However, we have observed the exact opposite: the smaller the time step, the larger the obtained instabilities. The fact that reducing the time step makes things even worse is thus another clear sign that our problems are the result of amplified error accumulation up to the point that it influences flow branching.

## 4 Large structure simulation (LSS): an attempt to get rid of instabilities in LEMs

In this section, we explain how to transpose the ideas underlying the concept of large eddy simulation from the computational fluid dynamics community to our landscape evolution model. The method consists in preventing any self-amplification phenomena that might emerge from the small spatial scales where numerical errors develop. In our opinion, this is a key ingredient for achieving reproducible LEM simulations.
<figure><img src="image_6.png" /><figcaption>Figure 5. Sediment height $h_s^{ex}$ of the analytic solution for the case $k_g = 5 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 1 \text{ km}^2 \text{ Myr}^{-1}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Sediment height $h_s^{ex}$ of the analytic solution for the case $k_g = 5 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 5 \text{ km}^2 \text{ Myr}^{-1}$.</figcaption></figure>

## 4.1 Principles and physical interpretation of filtering

Recall that the main idea of LES is to filter the solution to distinguish between the behavior of the flow above and below a target length scale $\alpha$ to obtain local averages that are smoother and as mesh-independent as possible. This target length scale controls the size of the smallest structures that we will be able to resolve in the problem, quite independently of the domain size. The main practical consequence is that our mesh will have to resolve this length scale; i.e., the mesh size $\epsilon$ will have to be smaller than the chosen length scale.

LES models are probably as numerous as the various authors working on the subject (Berselli et al., 2005); thus, we will be very brief on the subject and refer the reader to the recent review by Zhiyin (2015). The very first LES model is called the Leray-$\alpha$ model. It was used by Leray in 1934 to establish the existence of weak solutions to the Navier-Stokes equations (Leray, 1934). Originally, the filtering in Leray (1934) as well as in many classical LES models was

achieved by using a convolution operator $\mathcal{F}$ defined by

$$ \mathcal{F}(u)(\mathbf{x}) = \int_{\mathbb{R}^d} u(\mathbf{y})g_{\alpha}(\mathbf{x} - \mathbf{y})d\mathbf{y}, \quad \text{where } g_{\alpha}(\mathbf{x}) = \frac{1}{\alpha^d}g\left(\frac{\mathbf{x}}{\alpha}\right), $$

where the filter kernel $g$ satisfies

$$ 0 \le g(\mathbf{x}) \le 1, \quad g(\mathbf{0}) = 1, \quad \int_{\mathbb{R}^d} g(\mathbf{x})d\mathbf{x} = 1. $$

Several kernels are used in the literature, such as a low-pass filter, a box filter, or the very natural Gaussian filter $g(\mathbf{x}) = \pi^{-d/2}e^{-|\mathbf{x}|^2}$. In Fig. 13 we illustrate the smoothing effect of a Gaussian kernel on oscillating data: as expected, it preserves the high-amplitude and low-frequency oscillation while filtering out the high-frequency and low-amplitude oscillations. Such filters might therefore be ideal for our application to landscape evolution models: the small topographic perturbations will be cleaned out such that the flow routing will not be affected by it. Although convolution operators produce averages with the desired properties, they are impractical on bounded domains. The modern way of defining the Leray-$\alpha$ filter for bounded domains consists in using
<figure><img src="image_8.png" /><figcaption>Figure 7. Sediment height $h_s^{ex}$ of the analytic solution for the case $k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 5 \text{ km}^2 \text{ Myr}^{-1}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Sediment height $h_s^{ex}$ of the analytic solution for the case $k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$.</figcaption></figure>

the differential filter $\mathcal{F}_\alpha$ defined by (Cheskidov et al., 2005; Guermond et al., 2003)

$$ \left| \begin{array}{ll} -\alpha^2 \Delta \mathcal{F}_\alpha(u) + \mathcal{F}_\alpha(u) = u & \text{in } \Omega, \\ \nabla \mathcal{F}_\alpha(u) \cdot \mathbf{n} = 0 & \text{on } \partial \Omega_N, \\ \mathcal{F}_\alpha(u) = 0 & \text{on } \partial \Omega_D \end{array} \right. \qquad (16) $$

The filtered result $\mathcal{F}_\alpha(u)$ basically amounts to a convolution of $u$ by Green's function underlying Eq. (16), i.e., the filter applied to the Dirac distribution. Using a finite-volume scheme $\mathcal{F}_\alpha$, this time we can easily obtain a discrete version $\mathcal{F}_{\alpha,h}$, which is one of the main reasons why we have chosen to use this filter, along with its theoretical and practical success for CFD. Notice that contrary to Cheskidov et al. (2005) and Guermond et al. (2003), we use homogeneous Neumann and Dirichlet boundary conditions instead of periodic boundary conditions to simplify the treatment of the boundary. The main drawback of this choice is that our filter does not commute with differential operators. Resorting to only Dirichlet boundary conditions would have solved this issue, but from our numerical experiments we found that

this can create boundary effects unless the chosen Dirichlet boundary condition is adapted to the filtered quantity. The Neumann choice avoids those difficulties without creating any practical issues, which has motivated our choice. For quantities such as the water flux for which Neumann everywhere is a more natural boundary condition, we introduce the alternative filter $\mathcal{F}_\alpha^N$ with only Neumann boundary conditions:

$$ \left| \begin{array}{ll} -\alpha^2 \Delta \mathcal{F}_\alpha^N(u) + \mathcal{F}_\alpha^N(u) = u & \text{in } \Omega, \\ \nabla \mathcal{F}_\alpha^N(u) \cdot \mathbf{n} = 0 & \text{on } \partial \Omega. \end{array} \right. \qquad (17) $$

## 4.2 Leray filtering applied to our landscape evolution model

Our numerical observations suggest that for large values of $\tau$ the model governing the simultaneous evolution of sediment and water is as intractable to solution as the Navier-Stokes system is for large Reynolds numbers. Consequently, mimicking the idea of large eddy simulation (LES) we will now
<figure><img src="image_10.png" /><figcaption>Figure 9. Convergence curves. (a) Case $k_g = 50 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 1 \text{ km}^2 \text{ Myr}^{-1}$. (b) Case $k_g = 5 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 1 \text{ km}^2 \text{ Myr}^{-1}$. (c) Case $k_g = 5 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 5 \text{ km}^2 \text{ Myr}^{-1}$. (d) Case $k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 5 \text{ km}^2 \text{ Myr}^{-1}$. (e) Case $k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$.</figcaption></figure>

apply filtering to key parts of our model problem to obtain a more numerically stable approximate model. This means that the sediment flux used in the mass conservation equations,

$$ 
\left| 
\begin{array}{ll}
\frac{\partial h_s}{\partial t} + \operatorname{div}(\mathbf{J}_s) = S_s & \text{in } \Omega \times ]t_0, T[, \\
-\mathbf{J}_s \cdot \mathbf{n} = B_s & \text{on } \partial\Omega_N \times ]t_0, T[, \\
h_s = 0 & \text{on } \partial\Omega_D \times ]t_0, T[, \\
h_s(t=t_0) = h_{s,0} & \text{in } \Omega,
\end{array}
\right. 
$$

will now be given by

$$ 
\mathbf{J}_s = -\eta_s(h_s)s_{\text{ref}}^{-p_s} \left\| \nabla(h_s+b) \right\|^{p_s} \left( \left( \frac{\mathcal{F}_{\alpha}^{\mathcal{N}}(q_w)}{q_{\text{ref}}} \right)^{r_s} \nabla\psi_w(h_s+b) + \nabla\psi_g(h_s+b) \right) \text{ in } \Omega \times ]t_0, T[, \quad (18)
$$

where we use the filtered water flux magnitude $\mathcal{F}_{\alpha}^{\mathcal{N}}(q_w)$ instead of directly using the water flux $q_w$. In the same way, in the water equations, we will now use the filtered topography
<figure><img src="image_11.png" /><figcaption>Figure 10. Comparison between the sediment height $h_s^{\text{ex}}$ of the analytic solution and numerical solution $h_s$ for the case $k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$, $k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$, $r_s = 3/2$, and $p_s = 1$. (a) Analytic solution $h_s^{\text{ex}}$ and (b) numerical solution $h_s$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Comparison between the sediment height $h_s^{\text{ex}}$ of the analytic solution and numerical solution $h_s$ for the case $k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$, $k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$, $r_s = 2$, and $p_s = 0$. (a) Analytic solution $h_s^{\text{ex}}$ and (b) numerical solution $h_s$.</figcaption></figure>

$\mathcal{F}_{\alpha}(h_s + b)$ instead of the topography $h_s + b$, leading to

$$ 
\left| \begin{array}{l} 
\mathbf{Q}_w = -\frac{k_m h_w \eta_w(h_w)}{s_{\text{ref}}^{p_w}} ||\nabla(\mathcal{F}_{\alpha}(h_s+b))||^{p_w} \nabla(\mathcal{F}_{\alpha}(h_s+b)), \\ 
\text{div}(\mathbf{Q}_w) = S_w \quad \text{in } \Omega, \\ 
\mathbf{Q}_w \cdot \mathbf{n} = B_w \quad \text{on } \partial\Omega_{\text{in}}, 
\end{array} \right. 
\qquad (19) 
$$

with the associated water flux:

$$ q_w = ||k_m h_w \eta_w(h_w) s_{\text{ref}}^{-p_w}||^{p_w} \nabla(\mathcal{F}_{\alpha}(h_s+b)) ||. \qquad (20) $$

Our so-called large structure simulation (LSS) for landscape evolution thus consists in solving Eqs. (2)–(18)–(19)–(20).

The name “large structures” originates from the fact that since we use filtering in the coupling process, the water model does not see anymore topographic details that are smaller than $\alpha$, and in the same way the sediment evolution is no longer influenced by water flow details smaller than $\alpha$. We have thus abandoned the idea of resolving all the scales involved in the landscape evolution problem and will only try to simulate the large sedimentary and water structures, hence the name LSS: only structures several times larger than the filter resolution $\alpha$ will appear in the final result.

## 4.3 Numerical results with filtering

Before turning to numerical experiments, one must choose a value for the filter parameter $\alpha$. Following LES principles, we know that the filter scale $\alpha$ corresponds to the spatial res-
<figure><img src="image_13.png" /><figcaption>**Figure 12.** The “three rivers” test case with $\Delta_{xy} = 2$ km. (a, b) Initial topography; black arrows represent the position of the water inflows. (c–e) Topography and water flux after 6 Myr obtained under different numerical settings. (c) Sequential GMRES, (d) parallel GMRES, and (e) sequential BiCGStab.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Illustration of the effect of the convolution by a Gaussian function.</figcaption></figure>

olution of our continuous approximate model, which in practice one will want to be as small as possible. However, it must naturally be resolved by the grid, meaning we should have at the very least $\Delta_{xy} < \alpha$ for Cartesian grids. Moreover, as we test our numerical solution against an analytic solution for the unfiltered case, we need to make the filter size go to zero at the same speed as the mesh size in order to measure a convergence. For simplicity, we have chosen to use filter parameters $\alpha = \gamma \Delta_{xy}$ with $\gamma > 1$. In Fig. 14 we present the convergence results obtained for the analytic test cases of Sect. 3.1 this time using filters. Convergence is recovered with $\alpha = 1.1\Delta_{xy}$ (i.e., $\gamma = 1.1$) for every case that was already working without a filter, suggesting that the LSS approach at least does not deteriorate correct re-

sults previously obtained. We also see that for the test cases with ($k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$, $k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$) convergence is now obtained for $\alpha = 8\Delta_{xy}$ ($\gamma = 8$). This choice for the ratio $\gamma$ between the filter size $\alpha$ and the mesh size $\Delta_{xy}$ is not random. Indeed, with $\alpha = \gamma \Delta_{xy}$ when $\Delta_{xy}$ tends to zero so does the filter size, and if $\gamma$ is not large enough then the filtering parameter $\alpha$ will no longer be large enough to compensate for solver errors and numerical approximation errors. We illustrate this in Fig. 15. Keeping in mind that we are necessarily using a fixed Newton nonlinear solver tolerance ($1 \times 10^{-6}$ in practice), what we observe on those curves is that when the parameter $\alpha$ becomes smaller than some threshold value that allows controlling the corresponding accumulated solver (and numerical approximation) errors, the obtained solution is no longer correct. Of course, with a larger value of $\gamma$ this threshold is reached for a smaller value of $\Delta_{xy}$, which explains why once $\gamma$ is large enough we can obtain the correct solution along the entire convergence curve. This threshold is likely to depend on $\Delta_{xy}$ in the sense that for finer meshes since the size of the system is larger, so is the solver error. It is also very likely that since we expect larger values of $\tau$ to imply an increase in both the numerical approximation and solver errors, modifying $\tau$ might also influence this threshold value. Results shown in Fig. 15 confirm this behavior. This also explains why we have reported results with $\gamma = 8$ in Fig. 14e: to maintain the convergence over the full range of $\Delta_{xy}$ values used for these simulations despite high $\tau$ values. We nevertheless see in Fig. 15 that for more realistic mesh sizes, smaller values of $\gamma$ will be more than enough to obtain the correct solution and that using filters is not pro-
hibitively costly in realistic configurations. We also observe
that for mesh sizes allowing all the values of the ratio $\gamma$ to
give a correct approximation, the error of course increases
with $\gamma$, which is perfectly expected since $\alpha$ is our largest ap-
proximation parameter.

We finally reproduce the very same experiment that was
performed on the “three rivers” test case with sequential GM-
RES, parallel GMRES, and sequential BiCGStab but using a
filter $\alpha = 2.2$ km for $\Delta_{xy} = 2$ km. Contrary to Fig. 12, the
symmetry is maintained and we obtain almost identical re-
sults for the three configurations (Fig. 16). The expected im-
pact of the filter on the simulated water flow and topogra-
phy is a smoothing effect, which is what is observed when
comparing, for example, the width of the three valleys. How-
ever, the differences remain marginal in this case. Still on our
“three rivers” test case, from our observations on the analytic
cases and as we do not know the exact solution, to assess the
legitimacy of our choice of filter size we analyze the behav-
ior of the solution for various values of the filter parameter $\alpha$
(fixing the grid size to $\Delta_{xy} = 2$ km). Results are displayed in
Fig. 17. We clearly see that symmetric solutions are obtained
for $\alpha > \Delta_{xy}$, while further reducing the filter parameter leads
to behavior similar to the no-filter case. This is first coher-
ent with the principle of LES that the filter should control
what happens below the grid scale, which can only be done
if $\alpha > \Delta_{xy}$, and also a clear sign that our initial choice for
the ratio $\gamma = \alpha/\Delta_{xy}$ belongs to the stable region.

The above results obtained with a filtering strategy repre-
sent, in our opinion, a drastic improvement in the reliability
of these numerical solutions: the anomalous error amplifi-
cation has disappeared, and the results are reproducible and
unaffected by the choice of solver or the number of proces-
sors used for simulation. In the absence of this treatment,
the error amplification phenomenon is most probably a great
source of “anomalous mesh dependency” of the results as
each mesh induces different numerical errors. The positive
impact of our filtering strategy can be seen in the conver-
gence curves (Fig. 15): it allows getting rid of any anoma-
lous mesh dependency and recovering the regular one – that
is to say when refining the mesh the associated sequence of
solutions converges to the correct continuous solution. No-
tice that this is the best kind of “mesh independence” that
one can hope for: in particular, quite large differences will
remain when comparing two simulations defined on two dif-
ferent coarse meshes. Let us avoid any misunderstanding:
the strategy prevents the amplification of numerical errors,
but it does not clean the solution of legitimate numerical
errors. This being said, we believe that following our ap-
proach brings us closer to the best possible mesh indepen-
dence. When LEMs are designed without any filtering strat-
egy, it corresponds to implicitly considering the mesh size
as a cut-off length scale. However, it lacks calibration and
leads to the error self-amplification problems illustrated in
Sect. 3. In our case, we explicitly consider the filter size as
a cut-off length scale, the resolution of the model being then

controlled by the filter size. We will illustrate in the next sec-
tions that the calibration of this length must simply respect
an elementary principle: to be largely lower that the size of
the geomorphic structures the LEM aims to reproduce. Based
on this calibration, the mesh resolution must be chosen so
that the filter is correctly discretized. Thus, we have not de-
teriorated the overall computational situation: we still have
a unique discretization parameter that governs the resolution
and computational cost of the model.

## 4.4 Impacts of filtering on the emergence of geomorphic structures

We now consider two synthetic case studies to observe the
formation of geomorphic features. The idea underlying the
first test case is very simple: we re-use as our initial data
the analytic solution described in Sect. 3.1 in the case ($r_s =
2, p_s = 0$) and ($k_g = 1 \text{ km}^2 \text{ Myr}^{-1}, k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$) as
well as the rectangular domain described in Fig. 3. However,
instead of using the analytic source terms allowing the recov-
ery of the analytic solution for all times, we simply use a con-
stant source term ($S_s = 10 \text{ m Myr}^{-1}$, $S_w = 1 \text{ m}^3 \text{ s}^{-1} \text{ km}^{-2}$),
corresponding to a uniform constant uplift supply and a uni-
form constant rain.

We fix the mesh size to $\Delta_{xy} = 0.005$ km, and we again
perform the simulation over a time period of 0.25 Myr with
maximum time steps of length $\Delta t = 0.002$ Myr. In Fig. 18,
we recall the initial elevation corresponding to our analytic
solution along with the final solution obtained for our now
constant source terms for various values of the filter size
as well as without filters. Since our new source terms and
the analytic ones are of the same magnitude and since ev-
ery other property of the problem is kept the same, we can
guess, using the convergence curves in Fig. 15, which filter
sizes are giving a correct solution (up to the approximation
due to the filtering process itself). Since $\ln \Delta_{xy} \approx -5.298$,
we see in Fig. 15 that for our choice of $\Delta_{xy}$ we can be con-
fident that the filter size $\alpha = 2\Delta_{xy}$ will give us the correct
solution with a small numerical approximation error, and we
use this case as a reference. Thus, the first observation for
the result obtained with $\alpha = 2\Delta_{xy}$ is that the correct solu-
tion this time allows some legitimate geomorphic structures
to appear and self-organize. Those structures originate from
the bumps because if we perform the very same simula-
tion with constant source terms but without bumps, we ob-
tain a clean uniform final state deprived of any geomorphic
complexity. With the larger filter size $\alpha = 4\Delta_{xy}$, we obtain
an averaged version with slightly less geomorphic complex-
ity, illustrating the way the filter only keeps “large” struc-
tures. However for the very large filter $\alpha = 8\Delta_{xy}$ the ap-
proximation for $\Delta_{xy} = 0.005$ km is too crude and we lose
all the geomorphic complexity. We have checked that if we
refine the mesh we recover the correct solution with the ratio
$\alpha = 8\Delta_{xy}$. This confirms that the uniform crude approxima-
tion obtained for $\alpha = 8\Delta_{xy}$ and $\Delta_{xy} = 0.005$ km in Fig. 18
<figure><img src="image_15.png" /><figcaption>Figure 14. Convergence curves with filters. (a) Case $k_g = 50 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 1 \text{ km}^2 \text{ Myr}^{-1}$. (b) Case $k_g = 5 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 1 \text{ km}^2 \text{ Myr}^{-1}$. (c) Case $k_g = 5 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 5 \text{ km}^2 \text{ Myr}^{-1}$. (d) Case $k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 5 \text{ km}^2 \text{ Myr}^{-1}$. (e) Case $k_g = 1 \text{ km}^2 \text{ Myr}^{-1}$ and $k_w = 50 \text{ km}^2 \text{ Myr}^{-1}$.</figcaption></figure>

results, as expected, from an oversizing of $\alpha$. Now, let us consider the final solutions of Fig. 18 for the value $\alpha = 1.1 \Delta_{xy}$ as well as without a filter. Both of those results present more complexity than the reference case $\alpha = 2 \Delta_{xy}$. Using the convergence curves of Fig. 15, we expect the result obtained for $\alpha = 1.1 \Delta_{xy}$ to belong to the hazardous region where the error level starts to increase and this solution, while not completely erroneous, is becoming untrustworthy. However, for the solution without a filter strange small structures appear and the overall topography, despite being the most complex of all, does not have any physical origin.

We now switch to a second synthetic case study. The numerical domain again corresponds to a rectangular grid but this time with dimensions $Lx = 600 \text{ km}$ on the x axis and $Ly = 80 \text{ km}$ on the y axis containing a mesh of resolution $\Delta_{xy} = 0.25 \text{ km}$. The basement $b$ is constant and equal to 0 m, while the sediment thickness $h_s$ is initially given by a uniform in x smooth bump:

$$g(x, y) = \begin{cases} H \exp\left(\frac{-1}{1-r_y^2}\right) & \text{for } r_y = \frac{(y-y_c)}{\delta_y} \le 1, \\ 0 & \text{otherwise,} \end{cases}$$
<figure><img src="image_16.png" /><figcaption>Figure 15. Convergence curves for various values of the ratio α/Δxy: (a) (r_s, p_s) = (3/2, 1) and (b) (r_s, p_s) = (2, 0).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16. The “three rivers” test case with filter α = 2.2 km and Δxy = 2 km. Topography and water flux after 6 Myr. (a) Sequential GMRES, (b) parallel GMRES, and (c) sequential BiCGStab.</figcaption></figure>

with H = 20 m, y_c = 40 km, and δ_y = 20 km. This symmetry in the x direction of the initial topography is then perturbed by N_b = 30 small smooth bumps randomly positioned at points (x_p, y_p):

$$g_{\text{pert}}(x, y) = \begin{cases} H_{\text{pert}} \exp\left(\frac{-1}{1-r^2}\right) & \text{for } r^2 = \frac{(x-x_p)^2}{\delta^2} + \frac{(y-y_p)^2}{\delta^2} \le 1, \\ 0 & \text{otherwise,} \end{cases}$$

with H\_pert = 1 m and δ = 2 km. Rainfall is constant in time and space (3000 mm yr⁻¹) and is the unique water supply for this case. The sediment source (here we simulate a sediment production) goes from S_s = 0 m Myr⁻¹ at y = 0 and y = L_y to S_s = 100 m Myr⁻¹ at y = L_y/2 = y_c. The variation is continuous over the whole domain following

$$S_s(x, y) = \begin{cases} S_{\text{max}} \exp\left(\frac{-1}{1-r_y^2}\right) & \text{for } r_y = \frac{(y-y_c)}{\delta_y} \le 1 \\ 0 & \text{otherwise,} \end{cases}$$

with δ_y = 40 km. Model boundary conditions are fixed elevations on the sides normal to the x axis and a zero gradient

on the sides normal to the y axis. Model parameters controlling the nonlinearity in the water-sediment coupling are set as r_s = 2 and p_s = 0. Simulation takes place over the time period T = 6 Myr.

This second synthetic case study has similarities to the previous one in terms of boundary conditions, but its larger spatial scale makes it relatively close to the case studies published in Perron et al. (2008) and Armitage (2019). We display the initial topography (Fig. 19a) as well as the final topography obtained with and without a filter for k_w = 5 km² Myr⁻¹ and for three different k_g values. The first case considers k_g = 50 km² Myr⁻¹. The relatively high k_g value compared to k_w should not favor the emergence of geomorphic structures. This is, however, not what we observe in the simulation performed without a filter (Fig. 19b). The filter, defined by α = 0.3 km, has a huge impact and no geomorphic structure is produced (Fig. 19c). Refining both α and Δxy for this value of k_g, we always obtain the same uniform solution, which confirms that this is the correct one. An order of magnitude smaller k_g coefficient is used for the second simulation. By decreasing k_g, the emergence of structure may be considered a realistic result. In this case, complex structures
<figure><img src="image_18.png" /><figcaption>Figure 17. The “three rivers” test case with $\Delta_{xy} = 2$ km. Final topography and water flux after 6 Myr obtained with different values of the filter parameter $\alpha$. (a) No filter, (b) $\alpha = 0.2$ km, (c) $\alpha = 1$ km, (d) $\alpha = 2.2$ km, (e) $\alpha = 2.5$ km, and (f) $\alpha = 3$ km.</figcaption></figure>

controlled by at least one wavelength appear in the simulation performed without a filter (Fig. 19d). The effect of the filter, however, indicates the very likely artificial origin of these structures. A residual perturbation can still be observed in the final topography (Fig. 19e), indicating that this $k_g$ and $k_w$ configuration is at the transition between two regimes, the gravity-driven and the water-driven erosion regimes. In our last simulation, we have decreased $k_g$ by a factor of 5 and we indeed observe the emergence of structures even when the filter is active (Fig. 19g). Here again the impact of the filter is important and allows keeping only what we believe to be the correct structures.

This last set of simulations shows the major impact of $k_g$ on the wavelength of the structures that can emerge from our simulation. We have also performed additional simulations (not shown here) using various $k_w$ values for a given $k_g$. The results show that $k_w$ must be high enough to make the structures appear, but they also show that $k_g$ was more important than $k_w$ in the wavelength control. We believe that a dedicated study should be conducted with our model to quantify these effects, but it is beyond the scope of this article. Such a complete study can be found in Perron et al. (2008). Even if it was performed using another LEM, similar conclusions to those drawn from his study are also expected in our case.

## 5 Discussion

We consider this work to belong to the common effort of the scientific community to harmonize landscape evolution models. It is our belief that most of our observations and practical recommendations can also be applied to a wider range of sediment evolution models than the one we use in this study. The implementation of the large structure simulation strategy should be accessible to every LEM satisfying (H1), (H2), and (H3). In particular, we believe that filtering would also be very useful for the models of Perron et al. (2009), Hooshyar and Porporato (2021a), and Porporato (2022) that take the general form

$$ 
\left\{ 
\begin{array}{ll} 
\displaystyle \frac{\partial h_s}{\partial t} + \text{div}(\mathbf{J}_s) = S_s & \text{in } \Omega \times ]t_0, T[, \\ 
-\mathbf{J}_s \cdot \mathbf{n} = B_s & \text{on } \partial\Omega_N \times ]t_0, T[, \\ 
h_s = 0 & \text{on } \partial\Omega_D \times ]t_0, T[, \\ 
h_s(t=t_0) = h_{s,0} & \text{in } \Omega, 
\end{array} 
\right. 
\quad (21) 
$$

with a source given by

$$ S_s = U - \kappa_w s_{\text{ref}}^{-ps,2} \left( \frac{q_w}{q_{\text{ref}}} \right)^{r_s} ||\nabla(h_s + b)||^{ps,2}, $$

with $U$ a sediment source term (or an uplift depending on the interpretation of $b$) and a sediment flux given by

$$ \mathbf{J}_s = -s_{\text{ref}}^{-ps} k_g ||\nabla(h_s + b)||^{ps} \nabla(h_s + b) \text{ in } \Omega \times ]t_0, T[. $$

The behavior of those models is relatively close to the model in Eqs. (2)–(9) that we have studied in detail here, with the
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Results for a mesh size $\Delta_{xy} = 0.005$ km. (a) Initial elevation. Final elevation: (b) no filter, (c) $\alpha = 1.1\Delta_{xy}$, (d) $\alpha = 2\Delta_{xy}$, (e) $\alpha = 4\Delta_{xy}$, (f) $\alpha = 8\Delta_{xy}$.</figcaption></figure>

main difference that the nonlinear term $q_w^{\prime s} ||\nabla(h_s + b)|^{p_s}$ appears as a reaction term rather than a diffusive term. In particular, for $p_w = -1$ the observations of linear stability for the model in Eq. (21) match the conclusion of the linear stability analysis of Smith and Bretherton (1972) and Smith et al. (1997). We can thus expect that the model in Eq. (21) will potentially suffer from similar numerical stability issues to the ones we analyzed in detail for the model in Eqs. (2)–(9), although this certainly requires a dedicated study before drawing conclusions. In particular, several elements can help keep the numerical errors under control: high-order space and time schemes, explicit time schemes, and specific solvers for the

water flow model avoiding inverting a linear system. Nevertheless, an immediate application of the LSS in this context consists of course of replacing $q_w$ by its filtered version $\mathcal{F}_{\alpha}^{\mathcal{N}}(q_w)$ in the second member of Eq. (21) and can only improve the numerical stability. We also believe that the $\xi-q$ model of Davy and Lague (2009) could benefit from a similar filtering strategy.

Correctly using filters requires some understanding of the scales involved in the model. Although this is not such an easy task in general, from generic guidelines concerning the relation between the filter size $\alpha$ and the precision of the results it is clear that the chosen filtering parameter $\alpha$ should
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Final topographies obtained for three different set of diffusive coefficients, systematically tested without filter and with a filter using $\alpha = 0.3$ km. (a) Initial perturbed topography. (b) Solution without filter for $(k_g, k_w) = (50, 5) \text{ km}^2 \text{ Myr}^{-1}$, (c) solution with filter for $(k_g, k_w) = (50, 5) \text{ km}^2 \text{ Myr}^{-1}$, (d) solution without filter for $(k_g, k_w) = (5, 5) \text{ km}^2 \text{ Myr}^{-1}$, (e) solution with filter for $(k_g, k_w) = (5, 5) \text{ km}^2 \text{ Myr}^{-1}$, (f) solution without filter for $(k_g, k_w) = (1, 5) \text{ km}^2 \text{ Myr}^{-1}$, (g) solution with filter for $(k_g, k_w) = (1, 5) \text{ km}^2 \text{ Myr}^{-1}$.</figcaption></figure>

resolve the main sediment structures that one wants to correctly represent in the flow, ideally fulfilling an equivalent of Nyquist's rule. For instance, if an essential valley is 1 km large, then $\alpha$ should be several times smaller (and ideally smaller than 100 m). A good practical test consists in comparing the filtered topography $F_\alpha(h_s + b)$ and the unfiltered one $h_s + b$. The structures of $h_s + b$ that one wants to simulate accurately should be preserved in $F_\alpha(h_s + b)$, of course in a smoother way. For instance, for a given value of $\alpha$ if a small topographic depression in which water could in principle flow is observed on $h_s + b$ but is absent in $F_\alpha(h_s + b)$, then if one really wants to capture water flow inside this “channel” the value of $\alpha$ must be reduced and the mesh refined accordingly if needed. The filter should in any case be able to clean numerical approximation and solver errors, implying that we should at the very least have $\gamma = \alpha/\Delta_{xy} > 1$ to correctly resolve the targeted $\alpha$ spatial scale. To allow the filter to cor-

rectly clean errors that could otherwise have a destabilizing effect on the final configuration, higher values of $\alpha$ should probably be used for increasing values of $\tau$. Nevertheless, our experiments illustrate that even quite small values of $\alpha$ allow cleaning the most relevant geomorphic features.

Notice that in the present paper, we have for simplicity always used uniform meshes with a constant $\Delta_{xy}$, hence obtaining a constant ratio $\gamma = \alpha/\Delta_{xy}$. As an immediate extension, one could resort to adaptive mesh refinement to refine the mesh in areas where $\tau$ becomes large and thus where numerical errors are more likely to be large, mitigating the increase in the system's size and thus the increase in the computational cost. In practice for constant coefficients $k_g$ and $k_w$ this would be equivalent to refining the mesh where water flow occurs. In addition, one could replace the constant parameter $\alpha$ by a space-time-variable coefficient $\alpha(x, y, t)$ in an adaptive filtering strategy, where the filter size could be
chosen in coherence with the local size of the structures one
wants the model to be able to reproduce.

5.1 Recovering realistic landscapes

In principle, the use of filters allows getting rid of the influ-
ence of numerical noise in the solution. An apparent draw-
back is that for unperturbed data, complex topographies will
no longer appear by themselves through the perturbations
induced by either the numerical approximation or the nu-
merical solvers. Moreover, natural landscapes exhibit some
heterogeneity even under a low-τ regime. This suggests an
ingredient is missing, and this ingredient is well-known by
geologists: the heterogeneity. Indeed heterogeneity is every-
where in nature and could be injected in such a model to
make realistic-looking topographies emerge. This idea is of
course not new but we propose investigating the effect of het-
erogeneity in the context of the hydro-sedimentary model we
use for this paper.

We would also like to stress what is in our opinion an es-
sential methodological key point: since natural landscapes
are full of complex heterogeneities that are not contained
in our data, it might be tempting to consider the amplified
numerical noise as a way to recover the complexity lack-
ing in our data. However, the statistical signature of numer-
ical noise will never coincide with physical observations: it
would be like throwing up the pieces of a puzzle and hoping
that they will correctly reconstruct the puzzle while falling
down. The key point is that a numerical simulation is not
supposed to directly reproduce natural observations but only
to compute a correct approximation of a model for a given
dataset. This is the model and its data that should reproduce
nature, and in this case if solved correctly the numerical so-
lution will be a useful approximation. If we do not willingly
add randomness in the data, the numerics should not intro-
duce it out of nowhere and bypass our modeling: it is not
reasonable to rely on numerical hazard to recover the miss-
ing elements in a model or its dataset. Worst of all, numeri-
cal noise lacks two essential modeling requirements: repro-
ducibility and explainability. The first one is missing because
numerical noise depends on the software and algorithms used
and the number of processors, among other factors. The sec-
ond one is missing because it is almost impossible to track
how the numerical errors are generated.

The first heterogeneity we consider here is injected into
the $k_m$ coefficient, reflecting variable soil rugosity. Since ac-
quiring a roughness map adapted to the spatial scales rele-
vant to our approach is difficult and probably not relevant
for a synthetic case study, we resort to an artificial yet effi-
cient trick, namely the Perlin noise (Perlin, 1985) that is often
used in animated movies or video games to produce realistic-
looking mountains or river networks. This type of noise
can easily be used to build isotropic heterogeneity maps
with controlled spatial scales. We thus consider our “three
rivers” test case using variable coefficients $k_m$ in space and

time (Fig. 20). Figure 20b illustrates a typical distribution
in space of the $k_m$ coefficients when using Perlin noise. The
range of values for the $k$ coefficient (from $k_m = 0.01 \text{ m s}^{-1}$
to $k_m = 10 \text{ m s}^{-1}$) is arbitrarily fixed while respecting real-
istic value ranges. Impacted by the heterogeneity in $k_m$, the
water flow is still distributed between neighboring cells ac-
cording to the gradient of the slope, but it also preferentially
chooses to enter the cell with the highest $k_m$, especially when
the slopes become gentle and relatively close between neigh-
bors. The flow then acquires a high degree of complexity de-
spite a filter which, set at $\alpha = 1.1 \Delta_{xy}$, makes it possible to
eliminate numerical errors.

The same approach can be applied to the other synthetic
test case used in Sect. 4.4 using $\alpha = 2\Delta_{xy}$: the simulations
are now performed with spatially and temporally varying
$k_m$ coefficients (the same range of $k_m$ values is used). Figu-
ures 21a and b and 22 show the initial and the final state of the
simulation with a special focus on the geomorphic structures
produced, which are clearly more complex when comparing
to the result shown in Fig. 18d.

Next, we introduce similar heterogeneity in the rain maps.
When we use solely rain heterogeneity incorporating the
same spatial scales as in the $k_m$ maps, the geomorphic struc-
tures produced are very similar to those obtained using only
the heterogeneous $k_m$ coefficients. The most visually satisfy-
ing result is obtained for a simulation using both variable $k_m$
and rain maps (Fig. 21c).

5.2 Overcoming the topographic limitations of the GMS model

The most obvious restriction for the GMS model in Eq. (4)
to be well-posed is that the topography should not have any
flat area; i.e., $||\nabla(h_s + b)||$ should not be identically zero on
any measurable subset of $\Omega$ (subsets with nonzero area). In-
deed, if $||\nabla(h_s + b)|| = 0$ on $\mathcal{O} \subset \Omega$ for $p_w \ge 0$ model in
Eq. (4) becomes identically 0 on $\mathcal{O}$, while for $p_w < 0$ model
in Eq. (4) (and thus also the model in Eq. 3) we have an
under-determined form. In both situations, $h_w$ and $a$ cannot
be computed on $\mathcal{O}$. Arguably, true flat areas are rare for real-
istic topographies, however, a numerical algorithm can pro-
duce a few cells for which the topography is flat. A more
common although less immediately obvious problem arises
from accumulation areas. This is quite easy to understand:
consider a bowl-shaped topography, with a flow of water
coming from the boundary of only one half of the bowl.
From the boundaries with inflow, all water will go straight
down to the bottom of the bowl and stop there, since the flow
can only progress if it finds a downhill direction. Water will
remain stuck at the bottom and will never flow in the sec-
ond half of the bowl. Since the models in Eqs. (4) and (3)
correspond to steady-state water models this implies an infi-
nite value for $h_w$ at the bottom of the bowl. To put this into
more mathematical terms with a very simple example, let us
consider the model in Eq. (4) in the simplified setting where
<figure><img src="image_21.png" /><figcaption>Figure 20. The “three rivers” test case with Perlin-noise-based coefficient $k_m$. (a) Final (at $T = 6$ Myr) elevation and associated water flow with $k_m$ coefficients variable in space and time. (b) $k_m$ coefficients at $T = 6$ Myr.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21. Results with filters and Perlin-noise-based $k_m$ coefficient. (a) Initial elevation, (b) final elevation with variable coefficient $k_m$, and (c) final elevation with variable coefficient and additional Perlin-noise-based perturbation of rainfall.</figcaption></figure>

$k_w = 1$, $p_w = 0$, $S_w = 0$, and $h_s(x, y, t) + b(x, y, t) = x^2 + y^2$ on the unit disk $\Omega = \{(x, y) \in \mathbb{R}^2 | x^2 + y^2 < 1\}$. The model in Eq. (4) becomes

$$ -\operatorname{div}(h_w \nabla (h_s + b)) = -\nabla h_w \cdot \nabla (h_s + b) - h_w \Delta (h_s + b) = -2x \partial_x h_w - 2y \partial_y h_w - 4h_w = 0, $$

leading to solutions of the form $h_w(x, y, t) = \frac{C}{(x+y)^2}$. Assume that the boundary influx is given by $B_w = 1$ for $y \ge 0$ and $B_w = 0$ otherwise. Then, in the half-domain $y \ge 0$ we get $h_w > 0$ with $h_w \to +\infty$ when $(x, y \to (0, 0))$, which is unphysical. In the half-domain $y < 0$, $h_w = 0$ and $q_w = 0$.
<figure><img src="image_23.png" /><figcaption>Figure 23. Comparison of the models in Eqs. (4) and (22) on the “three rivers” test case for α = 2.2 km and Δxy = 2 km.</figcaption></figure>

This illustrates the two problems: infinite values for the water height and a water flux that abruptly stops on the line y = 0.

This is reflected at the discrete level by an abrupt stop of the water flow at the bottom of accumulation areas. Moreover, this prevents recomputing a correct approximation of $h_w$ or a from the intermediate unknown used in MFD algorithms (the total outflow of a cell), since the coefficient relating this intermediate unknown to $h_w$ or a will be zero (see Coatléven, 2020, or Appendix C for details). This coefficient also cancels on flat areas. We can infer that this is a discrete indicator of what could be the weakest theoretical requirements on the topography for the models in Eqs. (3) and (4) to be well-posed: the absence of flat or accumulation areas.

The model in Eq. (4) being in fact a simplification of the shallow-water equation (see Sect. 5.2), this limitation can be seen as the price to pay to simulate the water flow mass conservation with a very low computational expense. At the cost of a higher computational time alternative models also derived from the shallow-water equation can be considered to overcome this limitation (see Sect. 5.2). Notice that in the following numerical experiments, we have been careful to only consider situations for which no well-posedness issues occur.

In the general setting, there is no reason why the sediments should evolve in such a way that one of the sufficient conditions in Eq. (8) is always fulfilled, which can lead to some nonphysical behavior of the GMS model in Eq. (4) and thus also the pure MFD algorithms. This can occur in two obvious situations: in an accumulation area (a topographic depression) or a flat area. In principle, water arriving at an accumulation area should create a “lake” whose bathymetry will be determined by a water balance between incoming flow, infiltration, and evaporation. If the surface reaches the threshold of the lake, then some water leaves the lake and the water flow restarts from the lake threshold. In flat areas, water will spread, diminishing its height until the full area is covered. To reproduce those effects that are not originally considered, implementations of the MFD algorithms all incorporate practical

workarounds. Thanks to our interpretation as the discretization of a continuous model, we can easily propose a generalization of Eq. (4) that overcomes those limitations by noticing that the model in Eq. (4) is in fact a simplification of the shallow-water equations with friction. Indeed, appropriately choosing the friction model and assuming that the mass conservation of water is at steady state, a quite general model arising from applying the hydrostatic approximation to the shallow-water equations would be to consider (see Appendix B)

$$ 
\begin{aligned}
\mathbf{Q}_w &= -\frac{k_m h_w \eta_w(h_w)}{s_{\text{ref}}^{p_w}} ||\nabla(h_w + h_s + b)||^{p_w} \nabla(h_w + h_s + b), \\
\text{div}(\mathbf{Q}_w) &= S_w \quad \text{in } \Omega, \\
\mathbf{Q}_w \cdot \mathbf{n} &= B_w \quad \text{on } \partial\Omega_{\mathcal{N}}, \\
h_w &= 0 \quad \text{on } \partial\Omega_{\mathcal{D}},
\end{aligned} 
\qquad (22) $$

with the associated water flux strength:

$$ q_w = |k_m h_w \eta_w(h_w) |s_{\text{ref}}^{-p_w}||\nabla(h_w + h_s + b)||^{p_w+1}. \qquad (23) $$

This is almost Eq. (4) except that it uses the hydraulic gradient instead of the topographic one. The assumption $\nabla(h_s + b) \approx \nabla(h_w + h_s + b)$, while valid on pronounced slopes, is obviously not valid anymore in accumulation areas (at equilibrium, the hydraulic gradient is almost zero, while the topographic gradient is large) and flat areas (where the topographic gradient is zero and the hydraulic one is not). The nonlinear model in Eq. (22) is thus a natural generalization of the GMS model in Eq. (4) with a built-in handling of accumulation and flat areas which no longer requires practical workarounds. However, the model in Eq. (22) does come with some drawbacks. The first one is that now we must choose the water mobility function $\eta_w$, as we are solving for the water height unknown. This will influence both the water network and the strength of the water flow. In the same way, the absolute value of the coefficient $k_m$ will now impact the strength of the water flux through $h_w$, while only its contrasts were relevant for the GMS model in Eq. (4). Thus, some fine tuning is required for Eq. (22) to produce meaningful results. The last and probably more important drawback is that, Eq. (22) being nonlinear in its unknown $h_w$, its discretization will be more involved and computationally expensive than for Eq. (4). Let us compare the results obtained with the original GMS model in Eq. (4) with the more involved hydrologic model in Eq. (22) on the “three rivers” test case using filters in both cases. The water mobility function $\eta_w$ for Eq. (22) is simply chosen as equal to 1 if $h_w$ is positive and 0 otherwise.

As we can observe in Fig. 23, if the two models of course do not produce exactly the same results, their general behavior remains very similar. Even more close results could certainly be obtained by finely tuning the mobility function. We do not want to explore this any further in the present paper and simply want to illustrate that while suffering from some
limitations, the GMS model in Eq. (4) and thus MFD algorithms remain a very strong and attractive approximation on suitable topographies.

## 6 Conclusions

After illustrating the numerical instabilities arising from the self-amplification phenomenon at the core of coupling overland flow and sediment evolution models, we have proposed mimicking the LES strategy for CFD computation in the context of landscape evolution models, relying on the well-known Leray-$\alpha$ differential filter. Numerical experiments assess whether filtering produces results robust to numerical perturbations. It is our belief that this “large structure simulation” (LSS) approach goes far beyond the specific model considered here and that any LEMs involving a coupling with surface water flow could benefit from it. Indeed, experiments performed without any filtering strategy have shown that it can become extremely difficult to distinguish between the imprint of numerical errors and physical processes. Provided suitable filter parameters and mesh sizes are used, only the nonphysical heterogeneity will disappear. The size of the filter has a real modeling meaning and corresponds to the minimum size of the physical heterogeneity that can be resolved. The choice of its value depends mainly on the accuracy required, which of course also depends on the computational power. Deploying this technique in a geomorphic model that has incorporated the correct water flux expression will make it possible to correct the anomalous mesh dependency so often evoked in papers dealing with LEM behavior. Indeed, if LEMs are designed to reproduce complex emergent phenomena arising from a set of simple physical processes and data, they are in no way intended to give a physical meaning to non-reproducible numerical noise, even if this noise once amplified by physical processes leads to “good-looking” results. The apparently missing visual complexity that previously arose from numerical noise can be physically re-introduced when heterogeneous data are considered. Similarly to LES models, we believe that a mathematical analysis and numerical analysis of the filtered model should be achievable. We hope to be able to publish such an analysis in a future paper. To complete this work, we also plan to use in our next study the full model capacity in building a multi-lithology realistic test case. Finally, pursuing the analogy with LES, an interesting perspective would be to analyze whether it is feasible to develop subfilter models to increase the filtered model accuracy when $\alpha$ is quite large in order to reduce the need for fine $\alpha$ and thus fine meshes and consequently the overall cost of the approach.

## Appendix A: Computational details related to perturbation theory

Recall that we denote ($h_{s,*}$, $h_{w,*}$) a reference solution of Eqs. (2)–(9)–(4) with sources ($S_{s,*}$, $S_{w,*}$) and ($h_{s,\delta}$, $h_{w,\delta}$), a perturbation of magnitude $\delta$ of this reference solution associated with the source perturbation ($S_{s,\delta}$, $S_{w,\delta}$), leading to the overall solution ($h_s$, $h_w$) = ($h_{s,*}$ + $h_{s,\delta}$, $h_{w,*}$ + $h_{w,\delta}$) for the perturbed source ($S_{s,*}$ + $S_{s,\delta}$, $S_{w,*}$ + $S_{w,\delta}$). Since both the perturbed and unperturbed solutions have to satisfy the boundary conditions, we deduce that the perturbation ($h_{s,\delta}$, $h_{w,\delta}$) itself also satisfies the same boundary conditions. Thus, injecting ($h_s$, $h_w$) into Eqs. (2)–(9), multiplying by $h_{s,\delta}$, and integrating by parts we get

$$ 
\begin{aligned}
\frac{d}{dt} \left( \frac{1}{2} \int_{\Omega} h_{s,\delta}^2 \right) &= - \int_{\Omega} \eta_s(h_s) s_{\text{ref}}^{-p_s} ||\nabla(h_s + b)||^{p_s} \\
& \quad \left( \left( \frac{q_w}{q_{\text{ref}}} \right)^{r_s} k_w + k_g \right) ||\nabla h_{s,\delta}||^2 + \int_{\Omega} \eta_s(h_{s,*}) s_{\text{ref}}^{-p_s} \\
& \quad ||\nabla(h_{s,*} + b)||^{p_s} \left( \left( \frac{q_{w,*}}{q_{\text{ref}}} \right)^{r_s} k_w + k_g \right) \nabla(h_{s,*} + b) \\
& \quad \cdot \nabla h_{s,\delta} - \int_{\Omega} \eta_s(h_s) s_{\text{ref}}^{-p_s} ||\nabla(h_s + b)||^{p_s} \left( \left( \frac{q_w}{q_{\text{ref}}} \right)^{r_s} k_w + k_g \right) \\
& \quad k_w + k_g \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta} + \int_{\Omega} S_{s,\delta} h_{s,\delta}.
\end{aligned}
\quad (A1) $$

When denoting

$$ j_s(u, v, w) = \eta_s(u) s_{\text{ref}}^{-p_s} ||\nabla v||^{p_s} \left( \left( \frac{w}{q_{\text{ref}}} \right)^{r_s} k_w + k_g \right), \quad (A2) $$

we recover the claimed equation governing the evolution of the perturbation’s total energy:

$$ 
\begin{aligned}
\frac{d}{dt} \left( \frac{1}{2} \int_{\Omega} h_{s,\delta}^2 \right) &= - \int_{\Omega} j_s(h_s, h_s + b, q_w) ||\nabla h_{s,\delta}||^2 \\
& \quad + \int_{\Omega} S_{s,\delta} h_{s,\delta} + \int_{\Omega} (j_s(h_{s,*}, h_{s,*} + b, q_{w,*}) \\
& \quad - j_s(h_s, h_s + b, q_w)) \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta}.
\end{aligned}
\quad (A3) $$

Then, if $k_g$ is much larger than $k_w$ and thus $\tau$ is very small, we have, assuming for simplicity that $\eta_s = 1$ from Taylor’s expansion,

$$ 
\begin{aligned}
A_{\delta} & \approx k_g s_{\text{ref}}^{-p_s} (||\nabla(h_{s,*} + b)||^{p_s} - ||\nabla(h_s + b)||^{p_s}) \\
& \quad \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta} + O(\tau) \approx -k_g s_{\text{ref}}^{-p_s} \\
& \quad (p_s ||\nabla(h_{s,*} + b)||^{p_s-2} ||\nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta}||^2) \\
& \quad + O(\tau) + O(\delta^3).
\end{aligned}
\quad (A4) $$
In the same way, if $k_w$ is much larger than $k_g$ then $\tau$ is also very large and we have, again using Taylor's expansion,

$$
\begin{aligned}
A_{\delta} & \approx -k_w s_{\text{ref}}^{-p_s} q_{\text{ref}}^{-r_s} (q_{w,*}^{r_s} ||\nabla(h_{s,*} + b)||^{p_s} - q_w^{r_s}||\nabla(h_{s,*} + b)||^{p_s} \\
& \qquad \nabla(h_s + b)||^{p_s} \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta} + O(1/\tau) \\
& \approx -k_w s_{\text{ref}}^{-p_s} r_s \frac{q_{w,*}^{r_s-1}}{q_{\text{ref}}^{r_s}} q_{w,\delta} ||\nabla(h_{s,*} + b)||^{p_s} \nabla(h_{s,*} + b) \\
& \qquad \cdot \nabla h_{s,\delta} - k_w s_{\text{ref}}^{-p_s} \left(\frac{q_{w,*}}{q_{\text{ref}}}\right)^{r_s} (p_s ||\nabla(h_{s,*} + b)||^{p_s-2} \\
& \qquad \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta}|^2) + O(1/\tau) + O(\delta^3).
\end{aligned}
\quad \text{(A5)}
$$

Going back to the general case for $\eta_s$, we have for small values of $\tau$

$$
\begin{aligned}
A_{\delta} & \approx -k_g s_{\text{ref}}^{-p_s} (\eta'_{s}(h_{s,*}) ||\nabla(h_{s,*} + b)||^{p_s} h_{s,\delta} \\
& \qquad \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta} + p_s \eta_s(h_{s,*}) ||\nabla(h_{s,*} + b)||^{p_s-2} \\
& \qquad \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta}|^2) + O(\tau) + O(\delta^3),
\end{aligned}
\quad \text{(A6)}
$$

while for large values of $\tau$ we have

$$
\begin{aligned}
A_{\delta} & \approx -k_w s_{\text{ref}}^{-p_s} r_s \frac{q_{w,*}^{r_s-1}}{q_{\text{ref}}^{r_s}} q_{w,\delta} \eta_s (h_{s,*}) ||\nabla(h_{s,*} + b)||^{p_s} \\
& \qquad \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta} - k_w s_{\text{ref}}^{-p_s} \left(\frac{q_{w,*}}{q_{\text{ref}}}\right)^{r_s} (\eta'_{s}(h_{s,*}) \\
& \qquad \qquad ||\nabla(h_{s,*} + b)||^{p_s} h_{s,\delta} \nabla(h_{s,*} + b) \cdot \nabla h_{s,\delta} \\
& \qquad \qquad + p_s \eta_s(h_{s,*}) ||\nabla(h_{s,*} + b)||^{p_s-2} |\nabla(h_{s,*} + b) \\
& \qquad \qquad \qquad \cdot \nabla h_{s,\delta}|^2) + O(1/\tau) + O(\delta^3).
\end{aligned}
\quad \text{(A7)}
$$

The behavior is roughly speaking the same, with the main difference that the additional term in $\eta'_s(h_{s,*})$ can also contribute with the wrong sign. Since $\eta'_s(h_{s,*})$ will be almost zero as soon as $h_{s,*}$ is large enough (see Eq. 11), this can only happen in regions where $h_{s,*}$ is close to zero (in particular near Dirichlet boundaries). In this case, the potential contribution to the instabilities is controlled by the magnitude of $|\eta'_s(h_{s,*})h_{s,\delta}| \le |h_{s,\delta}|/h_c$. If the perturbation is not amplified by other engines, which will be the case if $\tau$ is small, and if the parameter $h_c$ is not chosen to be too small (a typical valid value is 20 cm), then no severe instability can occur through this additional term. Thus, we can be confident that parameter $\tau$ will be the main criterion governing the appearance of instabilities even for our most general model.

## Appendix B: From shallow-water model to the steady-state hydrologic model in Eq. (22)

Recall that the shallow-water system is given by (see Birnin et al., 2001; Peton et al., 2020)

$$
\left|
\begin{array}{l}
\displaystyle \frac{\partial h_w}{\partial t} + \operatorname{div}(h_w \mathbf{u}_w) = 0, \\[1.5ex]
\displaystyle \frac{\partial}{\partial t} (h_w \mathbf{u}_w) + \operatorname{div}(h_w \mathbf{u}_w \otimes \mathbf{u}_w) + g h_w \nabla(h_s + b + h_w) + \\[1.5ex]
\displaystyle = -\kappa_w (h_w, ||\nabla(h_w + h_s + b)||) |\mathbf{u}_w|^{r_w} \mathbf{u}_w,
\end{array}
\right.
$$

(B1)

where $\mathbf{u}_w$ denotes the water speed, $g$ is the acceleration due to gravity, and $\kappa_w$ is the friction coefficient. Then, following Peton et al. (2020) and defining $H_{s,c}$ to be the characteristic sediment height, $H_{w,c}$ the characteristic water height, $L_c$ the characteristic domain length, and $T_c$ the characteristic time, as well as defining the nondimensional variables,

$$
\begin{gathered}
\hat{h}_s = \frac{h_s}{H_{s,c}}, \quad \hat{b}_s = \frac{b}{H_{s,c}}, \quad \hat{h}_w = \frac{h_w}{H_{w,c}}, \quad \hat{u}_w = \frac{T_c \mathbf{u}_w}{L_c}, \\
\hat{x} = \frac{x}{L_c}, \quad \hat{y} = \frac{y}{L_c}, \quad \hat{t} = \frac{t}{T_c},
\end{gathered}
$$

we see that Eq. (B1) is equivalent to

$$
\left|
\begin{aligned}
& \frac{\partial \hat{h}_w}{\partial t} + \operatorname{div}(\hat{h}_w \hat{u}_w) = 0, \\
& \frac{\partial}{\partial t}(\hat{h}_w \hat{u}_w) + \operatorname{div}(\hat{h}_w \hat{u}_w \otimes \hat{u}_w) \\
& \qquad + g \frac{H_{s,c} T_c^2}{L_c^2} \hat{h}_w \hat{\nabla}(\hat{h}_s + \hat{b}) + g \frac{H_{w,c} T_c^2}{L_c^2} \hat{h}_w \hat{\nabla}(\hat{h}_w), \\
& = -\kappa_w (h_w, ||\nabla(h_w + h_s + b)||) \frac{L_c}{H_{w,c}} \left(\frac{L_c}{T_c}\right)^{r_w-1} |\hat{u}_w|^{r_w} \hat{u}_w.
\end{aligned}
\right.
$$

The “shallow” hypothesis corresponds to assuming that $L_c/H_{w,c} \gg 1$, while the two numbers

$$F_{r,w} = \frac{L_c}{\sqrt{g H_{w,c} T_c}} \quad \text{and} \quad F_{r,s} = \frac{L_c}{\sqrt{g H_{s,c} T_c}}$$

are equivalent to Froude numbers for the water and sediment flows. For long-term sediment evolution, it is reasonable to assume that $F_{r,w} \ll 1$ and $F_{r,s} \ll 1$, i.e., that gravity is the dominant phenomenon. Combined with the shallow-water assumption this suggests neglecting the inertia terms in the nondimensional momentum balance, leading to the hydrostatic assumption:

$$
gh_w \nabla(h_s + b + h_w) = -\kappa_w (h_w, ||\nabla(h_w + h_s + b)||) |\mathbf{u}_w|^{r_w} \mathbf{u}_w. \quad \text{(B2)}
$$

Inverting Eq. (B2) we obtain the following expression for the water speed:

$$\mathbf{u}_w = -\mu_w (h_w, ||\nabla(h_w + h_s + b)||) \nabla(h_s + b + h_w), \quad \text{(B3)}$$

where

$$
\begin{aligned}
& \mu_w (h_w, ||\nabla(h_w + h_s + b)||) \\
& = \frac{g^{\frac{1}{r_w+1}} h_w^{\frac{1}{r_w+1}}}{\kappa_w (h_w, ||\nabla(h_w + h_s + b)||)^{\frac{1}{r_w+1}}} \\
& \qquad ||\nabla(h_s + b + h_w)||^{-\frac{r_w}{r_w+1}}.
\end{aligned}
\quad \text{(B4)}
$$
Thus, appropriately choosing the friction model, for instance by setting $r_w = 0$ and

$$ \kappa_w(h_w, ||\nabla(h_w + h_s + b)||) = \frac{gh_w}{k_m \eta_w(h_w) s_{\text{ref}}^{-p_w} ||\nabla(h_w + h_s + b)||^{p_w}}, \quad (B5) $$

and assuming that the mass conservation of water is at steady state, we obtain the following quite general hydrostatic approximation to the shallow-water equations:

$$ \left| \begin{array}{l} \mathbf{Q}_w = -\frac{k_m h_w \eta_w(h_w)}{s_{\text{ref}}^{p_w}} ||\nabla(h_w + h_s + b)||^{p_w} \nabla(h_w + h_s + b), \\ \text{div}(\mathbf{Q}_w) = S_w \quad \text{in } \Omega, \\ \mathbf{Q}_w \cdot n = B_w \quad \text{on } \partial\Omega_{\mathcal{N}}, \\ h_w = 0 \quad \text{on } \partial\Omega_{\mathcal{D}}, \end{array} \right. $$

with the associated water flux strength:

$$ q_w = |k_m h_w \eta_w(h_w) |s_{\text{ref}}^{-p_w} ||\nabla(h_w + h_s + b)||^{p_w+1}. $$

**Remark B.1.** The friction model in Eq. (B5) becomes singular when $||\nabla(h_w + h_s + b)|| = 0$. Thus, an alternative choice would be to use something like

$$ \kappa_w(h_w, ||\nabla(h_w + h_s + b)||) = \frac{gh_w}{k_m \eta_w(h_w) (\beta + s_{\text{ref}}^{-p_w} ||\nabla(h_w + h_s + b)||^{p_w})} $$

for some $\beta > 0$ (the same holds for function $\eta_w$ such that $\eta(0) = 0$). This alternative choice is probably more physical, as the term in $s_{\text{ref}}^{-p_w} ||\nabla(h_w + h_s + b)||^{p_w}$ can be interpreted as modeling some deceleration in accumulation areas. We have chosen to use Eq. (B5) to be as close as possible to the MFD algorithms of the literature.

## Appendix C: Finite-volume discretization

In this section we describe the full finite-volume discretization of the system in Eqs. (2)–(18)–(19)–(20). Let $\Omega$ be a bounded polyhedral connected domain of $\mathbb{R}^2$, whose boundary is denoted as $\partial\Omega = \Omega \setminus \Omega$. We recall the usual finite-volume notations describing a mesh $\mathcal{M} = (\mathcal{T}, \mathcal{F})$ of $\Omega$. The set of the cells of the mesh $\mathcal{T}$ is a finite family of connected open disjoint polygonal subsets of $\Omega$ such that $\bar{\Omega} = \cup_{K \in \mathcal{T}} \bar{K}$. For any $K \in \mathcal{T}$, we denote $|K|$ the measure of $K$, $\partial K = \bar{K} \setminus K$ the boundary of $K$, $\rho_K$ its diameter, and $\mathbf{x}_K$ its barycenter. The set of faces of the mesh $\mathcal{F}$ is a finite family of disjoint subsets of $\mathbb{R}^2$ included in $\bar{\Omega}$ such that, for all $\sigma \in \mathcal{F}$, its measure is denoted $|\sigma|$, its diameter $h_\sigma$, and its barycenter $\mathbf{x}_\sigma$. For any $K \in \mathcal{T}$, the faces of cells $K$ correspond to the subset $\mathcal{F}_K$ of $\mathcal{F}$ such that $\partial K = \cup_{\sigma \in \mathcal{F}_K} \bar{\sigma}$. Then, for any face $\sigma \in \mathcal{F}$, we denote $\mathcal{T}_\sigma = \{K \in \mathcal{T} | \sigma \in \mathcal{F}_K\}$ the cells of which $\sigma$ is a face. Next, for all cell $K \in \mathcal{T}$ and

all face $\sigma \in \mathcal{F}_K$ of cell $K$, we denote $\mathbf{n}_{K,\sigma}$ the unit normal vector to $\sigma$ outward to $K$, and $d_{K,\sigma} = |\mathbf{x}_\sigma - \mathbf{x}_K|$. The set of boundary faces is denoted $\mathcal{F}_{\text{ext}}$, while interior faces are denoted $\mathcal{F}_{\text{int}}$. Finally for any $\sigma \in \mathcal{F}_{\text{int}}$, whenever the context is clear we will denote $K$ and $L$ the two cells forming $\mathcal{T}_\sigma = \{K, L\}$, as well as $d_{KL} = |\mathbf{x}_K - \mathbf{x}_L|$. This, for instance, allows when looping over the faces $\sigma$ of cell $K$ denoting $L$ the other face of $\sigma$ without resorting to a notation that is too heavy. To avoid any confusion with water and sediment heights, $\epsilon = \max_{K \in \mathcal{T}} \rho_K$ will denote the mesh size. For any continuous quantity $u$, its discrete counterpart will be denoted $u_T = ((u_K)_{K \in \mathcal{T}}, (u_\sigma)_{\sigma \in \mathcal{F}_{\text{ext}}})$, where for any $K \in \mathcal{T}$ $u_K$ is the constant approximation of $u$ in cell $K$, while for any $\sigma \in \mathcal{F}_{\text{ext}}$ $u_\sigma$ is the constant approximation of $u$ over face $\sigma$.

In the following we will assume that the mesh is orthogonal, i.e., there exists a family of centroids $(\bar{\mathbf{x}}_K)_{K \in \mathcal{T}}$ such that

$$ \bar{\mathbf{x}}_K \in \mathcal{K} \quad \forall K \in \mathcal{T} \quad \text{and} \\ \frac{\bar{\mathbf{x}}_L - \bar{\mathbf{x}}_K}{|\bar{\mathbf{x}}_L - \bar{\mathbf{x}}_K|} = \mathbf{n}_{K,\sigma} \quad \text{for } \sigma \in \mathcal{F}_{\text{int}}, \sigma = \{K, L\}, $$

and let us denote $\bar{\mathbf{x}}_\sigma$ as the orthogonal projection of $\bar{\mathbf{x}}_K$ to the hyperplane containing $\sigma$ for any $\sigma \in \mathcal{F}_K$ and any $K \in \mathcal{T}$ with $\bar{d}_{K,\sigma} = |\bar{\mathbf{x}}_K - \bar{\mathbf{x}}_\sigma|$, as well as $\bar{d}_{KL} = |\bar{\mathbf{x}}_K - \bar{\mathbf{x}}_L|$. Then, one can use a two-point finite-volume scheme to discretize diffusion operators with scalar diffusion coefficients (no tensors).

We also assume that the mesh is compatible with the boundary decomposition, i.e., there exist subsets $(\mathcal{F}_{\text{ext}}^N)$ and $(\mathcal{F}_{\text{ext}}^D)$ such that

$$ \bar{\partial \Omega}_{\mathcal{N}} = \bigcup_{\sigma \in \mathcal{F}_{\text{ext}}^{\mathcal{N}}} \bar{\sigma} \quad \text{and} \quad \bar{\partial \Omega}_{\mathcal{D}} = \bigcup_{\sigma \in \mathcal{F}_{\text{ext}}^{\mathcal{D}}} \bar{\sigma}. $$

Notice that all our simulations without filters employ the same numerical schemes but of course replacing the filtered values by the original ones.

## C1 Leray-$\alpha$ filtering equation

Using the two-point flux approximation (TPFA) the approximate filter $\mathcal{F}_{\alpha,h}$ is defined for

$$ u_T = ((u_K)_{K \in \mathcal{T}}, (u_\sigma)_{\sigma \in \mathcal{F}_{\text{ext}}}) $$

by

$$ \mathcal{F}_{\alpha,h}(u_T) = ((\mathcal{F}_{\alpha,K}(u_T))_{K \in \mathcal{T}}, (\mathcal{F}_{\alpha,\sigma}(u_T))_{\sigma \in \mathcal{F}_{\text{ext}}}), $$
where

$$
\left\{
\begin{aligned}
& \alpha^2 \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}} \frac{|\sigma|}{\bar{d}_{KL}} (\mathcal{F}_{\alpha,K}(u_T) - \mathcal{F}_{\alpha,L}(u_T)) \\
& + |K| \mathcal{F}_{\alpha,K}(u_T) = |K| u_K \quad \text{for all } K \in \mathcal{T}, \\
& \mathcal{F}_{\alpha,\sigma}(u_T) = \mathcal{F}_{\alpha,K}(u_T) \quad \text{for all } K \in \mathcal{T} \\
& \text{and all } \sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}^{\mathcal{N}}, \\
& \mathcal{F}_{\alpha,\sigma}(u_T) = 0 \quad \text{for all } K \in \mathcal{T} \text{ and all } \sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}^{\mathcal{D}}.
\end{aligned}
\right.
\tag{C1}
$$

The discrete Neumann filter $\mathcal{F}_{\alpha,h}^{\mathcal{N}}$ of course satisfies Eq. (C1) but with Neumann boundary conditions on every $\sigma \in \mathcal{F}_{\text{ext}}$.

C2 Sediment mass conservation equations

We now assume that the time interval $[0, T]$ is subdivided into $N_T$ subintervals $[t_n, t_{n+1}]$, where $t_0 = 0$ and $t_{N_T+1} = T$. We denote $\Delta t^n = t_{n+1} - t_n$. The discrete quantities associated with time $t_n$ will be denoted as usual with a superscript $n$. The TPFA finite-volume scheme for the mass conservation of sediments in Eq. (2) for the flux in Eq. (18) is given by

$$
\left|
\begin{aligned}
& \frac{|K|}{\Delta t^n} (h_{s,K}^{n+1} - h_{s,K}^n) + \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}} \frac{|\sigma|}{\bar{d}_{KL} s_{\text{ref}}^{p_w}} \eta_{s,\sigma}^{n+1} \Delta \Psi_{KL}^{n,n+1} \\
& + \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}^{\mathcal{D}}} \frac{|\sigma|}{\bar{d}_{K\sigma} s_{\text{ref}}^{p_w}} \eta_{s,\sigma}^{n+1} \Delta \Psi_{K\sigma}^{n,n+1} \\
& - \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}^{\mathcal{N}}} |\sigma| B_{s,\sigma}^{n+1} = |K| S_{s,K}^n \quad \text{for all } K \in \mathcal{T}, \\
& h_{s,\sigma}^{n+1} + b_{\sigma}^{n+1} = h_{s,K}^{n+1} + b_{K}^{n+1} + \mathbf{G}_{s,K}^{n+1} \cdot (\bar{\mathbf{x}}_{\sigma} - \bar{\mathbf{x}}_{K}) \\
& \text{for all } K \in \mathcal{T} \text{ and all } \sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}^{\mathcal{N}}, \\
& h_{s,\sigma}^{n+1} = 0 \quad \text{for all } \sigma \in \mathcal{F}_{\text{ext}}^{\mathcal{D}},
\end{aligned}
\right.
\tag{C2}
$$

where

$$
\Delta \Psi_{KL}^{n,n+1} = (q_{w,\sigma}^{n+1})^{r_s} \left( \begin{aligned} & ||\mathbf{G}_{s,\sigma}^{\dagger,n+1}||^{p_s} (\psi_w(h_{s,K} + b_K) \\ & - \psi_w(h_{s,L} + b_L)) + ||\mathbf{G}_{s,\sigma}^{\dagger,n+1}||^{p_s} \\ & (\psi_g(h_{s,K} + b_K) - \psi_g(h_{s,L} + b_L)) \end{aligned} \right) \quad \text{(C3)}
$$

and

$$
\Delta \Psi_{K\sigma}^{n,n+1} = (q_{w,\sigma}^{n+1})^{r_s} \left( \begin{aligned} & ||\mathbf{G}_{s,\sigma}^{\dagger,n+1}||^{p_s} (\psi_w(h_{s,K} + b_K) \\ & - \psi_w(h_{s,\sigma} + b_\sigma)) + ||\mathbf{G}_{s,\sigma}^{\dagger,n+1}||^{p_s} \\ & (\psi_g(h_{s,K} + b_K) - \psi_g(h_{s,\sigma} + b_\sigma)), \end{aligned} \right) \quad \text{(C4)}
$$

where the mobility $\eta_{s,\sigma}^{n+1}$ taken upwind with respect to $\Delta\Psi_{KL}^{n,n+1}$ for $\sigma \in \mathcal{F}_{\text{int}}$,

$$
\eta_{s,\sigma}^{n+1} = \begin{cases} \eta_s (h_{s,K}^{n+1}) & \text{if } \Delta\Psi_{KL}^{n,n+1} \ge 0, \\ \eta_s (h_{s,L}^{n+1}) & \text{if } \Delta\Psi_{KL}^{n,n+1} < 0, \end{cases} \quad \text{(C5)}
$$

and using $\Delta\Psi_{K\sigma}^{n,n+1}$ for $\sigma \in \mathcal{F}_{\text{ext}}^{\mathcal{D}}$,

$$
\eta_{s,\sigma}^{n+1} = \begin{cases} \eta_s (h_{s,K}^{n+1}) & \text{if } \Delta\Psi_{K\sigma}^{n,n+1} \ge 0, \\ \eta_s (h_{s,\sigma}^{n+1}) & \text{if } \Delta\Psi_{K\sigma}^{n,n+1} < 0, \end{cases} \quad \text{(C6)}
$$

where the filtered water flux magnitude is approximated by the harmonic mean whenever possible and the mean value otherwise:

• if $\sigma \in \mathcal{F}_{\text{ext}}^{\mathcal{D}}$:

$$
q_{w,\sigma}^{n+1} = \mathcal{F}_{\alpha,K}^{\mathcal{N}}(q_{w,\tau}^{n+1})
$$

• if $\sigma \in \mathcal{F}_{\text{int}}$ and $\mathcal{F}_{\alpha,K}^{\mathcal{N}}(q_{w,\tau}^{n+1}) > 0$ and $\mathcal{F}_{\alpha,L}^{\mathcal{N}}(q_{w,\tau}^{n+1}) > 0$

$$
q_{w,\sigma}^{n+1} = \frac{\bar{d}_{KL} \mathcal{F}_{\alpha,K}^{\mathcal{N}}(q_{w,\tau}^{n+1}) \mathcal{F}_{\alpha,L}^{\mathcal{N}}(q_{w,\tau}^{n+1})}{\mathcal{F}_{\alpha,K}^{\mathcal{N}}(q_{w,\tau}^{n+1}) \bar{d}_{L\sigma} + \mathcal{F}_{\alpha,L}^{\mathcal{N}}(q_{w,\tau}^{n+1}) \bar{d}_{K\sigma}}
$$

• if $\sigma \in \mathcal{F}_{\text{int}}$ and $\mathcal{F}_{\alpha,K}^{\mathcal{N}}(q_{w,\tau}^{n+1}) = 0$ or $\mathcal{F}_{\alpha,L}^{\mathcal{N}}(q_{w,\tau}^{n+1}) = 0$

$$
q_{w,\sigma}^{n+1} = \frac{1}{2} (\mathcal{F}_{\alpha,K}^{\mathcal{N}}(q_{w,\tau}^{n+1}) + \mathcal{F}_{\alpha,L}^{\mathcal{N}}(q_{w,\tau}^{n+1}))
$$

The discrete full topographic gradient is given for any cell $K \in \mathcal{T}$ by

$$
\mathbf{G}_{s,K}^{n} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}} \frac{|\sigma|}{\bar{d}_{KL}} (h_{s,L}^{n} + b_{L}^{n} - h_{s,K}^{n} - b_{K}^{n}) \\
(\mathbf{x}_{\sigma} - \mathbf{x}_{K}) + \frac{1}{|K|} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} \frac{|\sigma|}{\bar{d}_{K\sigma}} (h_{s,\sigma}^{n} + b_{\sigma}^{n} \\
-h_{s,K}^{n} - b_{K}^{n}) (\mathbf{x}_{\sigma} - \mathbf{x}_{K}),
$$

while its stabilized version $\mathbf{G}_{s,\sigma}^{\dagger,n}$ is given by $\mathbf{G}_{s,\sigma}^{\dagger,n} = \mathbf{G}_{s,\sigma}^{n} + \mathbf{R}_{s,\sigma}^{n}$ with

$$
\mathbf{G}_{s,\sigma}^{n} = \begin{cases} \frac{1}{2} (\mathbf{G}_{s,K}^{n} + \mathbf{G}_{s,L}^{n}) & \text{if } \mathcal{T}_{\sigma} = \{K, L\}, \\ \mathbf{G}_{s,K}^{n} & \text{if } \mathcal{T}_{\sigma} = \{K\}, \end{cases} \quad \text{(C7)}
$$

as well as

• If $\mathcal{T}_{\sigma} = \{K, L\}$

$$
\mathbf{R}_{s,\sigma}^{n} =
$$

$$
\frac{1}{\bar{d}_{KL}^2} (h_{s,L}^n + b_L^n - h_{s,K}^n - b_K^n - \mathbf{G}_{s,\sigma}^n \cdot (\bar{\mathbf{x}}_L - \bar{\mathbf{x}}_K)) (\bar{\mathbf{x}}_L - \bar{\mathbf{x}}_K)
$$

• If $\mathcal{T}_{\sigma} = \{K\}$

$$
\mathbf{R}_{s,\sigma}^{n} =
$$

$$
\frac{1}{\bar{d}_{K\sigma}^2} (h_{s,\sigma}^n + b_\sigma^n - h_{s,K}^n - b_K^n - \mathbf{G}_{s,\sigma}^n \cdot (\bar{\mathbf{x}}_\sigma - \bar{\mathbf{x}}_K)) (\bar{\mathbf{x}}_\sigma - \bar{\mathbf{x}}_K)
$$
## C3 Water equations

The finite-volume scheme for the water equations in Eqs. (19)–(20) is simply obtained by applying the corrected MFD algorithm of Coatléven (2020) to the filtered topography and reconstructing a consistent water flux by setting $q_K^{n+1} = ||Q_K^{n+1}||$ with

$$ 
\begin{aligned}
Q_K^{n+1} = & \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}, \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) > \mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n)} \\
& \frac{\tau_{KL}^{n,n+1} \tilde{q}_K^{n+1}}{|K|s_K^{n,n+1}} (\mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) - \mathcal{F}_{\alpha,L}) \\
& \qquad (h_{s,\tau}^n + b_{\tau}^n)) (\mathbf{x}_{\sigma} - \mathbf{x}_K) \\
& - \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}, \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) < \mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n)} \\
& \frac{\tau_{KL}^{n,n+1} \tilde{q}_L^{n+1}}{|K|s_L^{n,n+1}} (\mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n) \\
& - \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n)) (\mathbf{x}_{\sigma} - \mathbf{x}_K) \\
& - \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\sigma| B_{w,\sigma}^{n+1}
\end{aligned}
\quad (C8) $$

and

$$ 
\begin{aligned}
\tilde{q}_K^{n+1} - & \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}, \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) < \mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n)} \\
& \frac{\tau_{KL}^{n,n+1} \tilde{q}_L^{n+1}}{s_L^{n,n+1}} (\mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n) - \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n)) \\
& - \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\sigma| B_{w,\sigma}^{n+1} = |K| S_{w,K}^n \quad \text{for all } K \in \mathcal{T}, \\
S_K^{n,n+1} = & \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}, \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) \ge \mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n)} \\
& \tau_{KL}^{n,n+1} (\mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) - \mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n)) \\
& \tau_{KL}^{n,n+1} = \frac{|\sigma| k_{m,\sigma}^{n+1}}{\bar{d}_{KL} s_{ref}^{p_w}} ||\mathbf{G}_{\mathcal{F},s,\sigma}^n||^{p_w},
\end{aligned}
\quad (C9) $$

where

$$ \mathbf{G}_{\mathcal{F},s,\sigma}^n = \begin{cases} \frac{1}{2}(\mathbf{G}_{\mathcal{F},s,K}^n + \mathbf{G}_{\mathcal{F},s,L}^n) & \text{if } \mathcal{T}_{\sigma} = \{K, L\}, \\ \mathbf{G}_{\mathcal{F},s,K}^n & \text{if } \mathcal{T}_{\sigma} = \{K\}, \end{cases} \quad (C10) $$

and the gradient of the filtered topography is of course given by

$$ 
\begin{aligned}
\tilde{q}_K^{n+1} - & \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}, \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) < \mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n)} \\
& \frac{\tau_{KL}^{n,n+1} \tilde{q}_L^{n+1}}{s_L^{n,n+1}} (\mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n) - \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n)) \\
& - \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} |\sigma| B_{w,\sigma}^{n+1} = |K| S_{w,K}^n \quad \text{for all } K \in \mathcal{T}, \\
S_K^{n,n+1} = & \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{int}}, \mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) \ge \mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n)} \\
& \tau_{KL}^{n,n+1} (\mathcal{F}_{\alpha,K} (h_{s,\tau}^n + b_{\tau}^n) - \mathcal{F}_{\alpha,L} (h_{s,\tau}^n + b_{\tau}^n)) \\
& \tau_{KL}^{n,n+1} = \frac{|\sigma| k_{m,\sigma}^{n+1}}{\bar{d}_{KL} s_{ref}^{p_w}} ||\mathbf{G}_{\mathcal{F},s,\sigma}^n||^{p_w},
\end{aligned}
$$

## Appendix D: Derivation of analytic solutions

For simplicity, we consider in this section the special case where $b=0$, $k_w$ and $k_g$ are constants, and the water mobility function and coefficient $k_m$ are both equal to 1 ($\eta_w(h_w) = 1$ and $k_m = 1$). To ease the reading, we will not write the dimension constants $s_{ref}$ and $q_{ref}$, as they are both equal to 1 in the chosen unit system. The sediment flux simplifies into

$$ J_s = -\eta_s(h_s) ||\nabla(h_s + b)||^{p_s} (q_w^r k_w \nabla(h_s + b) + k_g \nabla(h_s + b)) \quad \text{in } \Omega \times ]t_0, T[. \quad (D1) $$

We consider the simplified setting where $\eta_s(h_s) = 1$. This setting corresponds to the analytic steady-state solutions studied in Smith et al. (1997). Since $\eta_s(h_s) \approx 1$ as soon as $h_s$ is large enough, we label those solutions as “quasi-steady state”. We seek quasi-steady-state solutions that are moreover uniform in the y variable $h_s(x, y, t) = h_{s,x}(x)$ and symmetric with respect to the axis $x=0$, and we consider only the interval $]0, Lx/2[$. We finally assume that $S_s$ and $S_w$ are equal to two constants $S_{s,x}$ and $S_{w,x}$. We consequently have $\nabla(h_s + b) = \partial_x h_{s,x} e_x$, and the water equation reduces to

$$ -\partial_x (h_{w,x} |\partial_x h_{s,x}|^{p_w} \partial_x h_{s,x}) = S_{w,x}. \quad (D2) $$

Assuming $-\partial_x h_{s,w} > 0$ (the solution is decreasing from the center of the domain to its boundary) this leads to $\partial_x q_{w,x} = S_{w,x}$, $q_{w,x} = -h_{w,x} |\partial_x h_{s,x}|^{p_w} \partial_x h_{s,x}$, and finally $q_{w,x} = q_w(0) + S_{w,x}x$. In the same way, the conservation of sediments after integrating $x$ and again using our hypothesis on the sign of $\partial_x h_{s,x}$ reduces to

$$ \partial_x h_{s,x} = - \frac{(S_{s,x}x + \gamma)^{\frac{1}{p_s+1}}}{(k_g + k_w (q_w(0) + S_{w,x}x)^{\frac{r_s}{p_s+1}})} \quad (D3) $$
To ensure the continuity of the derivatives at $x = 0$, let us assume that $\partial_x h_{s,x}(0) = 0$ and thus $\gamma = 0$, and consequently $q_w(0) = 0$. The above relation simplifies into

$$ \partial_x h_{s,x} = -(S_{s,x}x)^{\frac{1}{p_s+1}} (k_g + k_w (S_{w,x}x)^{r_s})^{-\frac{1}{p_s+1}}. \quad (D4) $$

Notice that this is coherent with our assumption $-\partial_x h_{s,x} > 0$. The water height $h_{w,x}$ can then be obtained by setting

$$ h_{w,x} = (-1)^{p_w+1} \frac{S_{w,x}x}{\partial_x h_{s,x}^{p_w+1}} = (S_{w,x}x)(S_{s,x}x)^{\frac{-(p_w+1)}{p_s+1}} \\ (k_g + k_w (S_{w,x}x)^{r_s})^{\frac{p_w+1}{p_s+1}}, \quad (D5) $$

which is positive as expected. At this stage, integration for $h_{s,x}$ was simpler in Smith et al. (1997) because of the absence of $k_g$. Indeed, for $k_g = 0$ we immediately have

$$ h_{s,x} = h_{s,x}(0) - \frac{1+p_s}{2+p_s-r_s} (S_{s,x}k_w^{-1}S_{w,x}^{-r_s})^{\frac{1}{p_s+1}} x^{2+p_s-r_s}. \quad (D6) $$

Conversely, if $k_w = 0$ (no coupling between water and sediments), we get

$$ h_{s,x} = h_{s,x}(0) - \frac{p_s+1}{p_s+2} (k_g^{-1} S_{s,x})^{\frac{1}{p_s+1}} x^{\frac{p_s+2}{p_s+1}}. \quad (D7) $$

In the general case using the variable change $v = u^{r_s}$, $u = v^{1/r_s}$, and $du = \frac{1}{r_s}v^{(1-r_s)/r_s}dv$ we need to compute

$$ h_{s,x} = h_{s,x}(0) - \frac{1}{r_s} S_{s,x}^{\frac{1}{p_s+1}} \int_0^{x^{r_s}} v^{\frac{(1-r_s)(p_s+1)+1}{r_s(p_s+1)}} \\ (k_g + k_w S_{w,x}^{r_s} v)^{-\frac{1}{p_s+1}} dv, \quad (D8) $$

which will lead to easily computable analytic solutions, in particular for the special combinations of values of $r_s$ and $p_s$ that satisfy $(1-r_s)(p_s+1)+1 = 0$ and cancel the exponent $\frac{(1-r_s)(p_s+1)+1}{r_s(p_s+1)}$. In the special case $p_s = 0$ and $r_s = 2$, we have

$$ h_{s,x} = h_{s,x}(0) - \frac{S_{s,x}}{2k_w S_{w,x}^2} \left( \ln(k_g + k_w S_{w,x}^{r_s} x^2) - \ln(k_g) \right). \quad (D9) $$

In the other cases for which $(1-r_s)(p_s+1)+1 = 0$, this leads to

$$ h_{s,x} = h_{s,x}(0) - \frac{p_s+1}{p_s r_s} \frac{S_{s,x}^{\frac{1}{p_s+1}}}{k_w S_{w,x}^{r_s}} \\ \left( (k_g + k_w S_{w,x}^{r_s} x^{r_s})^{p_s/(p_s+1)} - k_g^{p_s/(p_s+1)} \right). \quad (D10) $$

Apart from those cases that cancel the exponent appearing in the integral, another interesting special case is the linear case

$p_s = 0$ and $r_s = 1$ for which we have

$$ h_{s,x} = h_{s,x}(0) - S_{s,x} \left( \frac{x}{k_w S_{w,x}} - \frac{k_g}{k_w^2 S_{w,x}^2} \ln |k_g| + k_w S_{w,x} x \right) + \frac{k_g}{k_w^2 S_{w,x}^2} \ln |k_g| \quad (D11) $$

It is then easy to choose the value for $h_{s,x}(0)$ such that $h_{s,x}(Lx/2) = 0$.

**Code availability.** All the numerical schemes used in this paper are fully described in Appendix C. Implementation was performed in code ArcaDES, which is available through the commercial simulator DionisosFlow™.

**Author contributions.** JC: conceptualization, writing, software, simulations. BC: writing, software, simulations.

**Competing interests.** Both the authors are core developers of the ArcaDES simulator supporting DionisosFlow™, a commercial stratigraphic simulator.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

**Acknowledgements.** The authors would like to thank John J. Armitage and Didier Granjeon for their careful reading of the present paper, as well as the referees for their numerous useful comments.

**Review statement.** This paper was edited by Wolfgang Schwanghart and reviewed by Stefan Hergarten and three anonymous referees.

## References

Armitage, J. J.: Short communication: flow as distributed lines within the landscape, Earth Surf. Dynam., 7, 67–75, https://doi.org/10.5194/esurf-7-67-2019, 2019.

Balay, S., Gropp, W., McInnes, L. C., and Smith, B. F.: PETSc, the portable, extensible toolkit for scientific computation, Argonne National Laboratory, 1998.

Bardos, C.: Problèmes aux limites pour les équations aux dérivées partielles du premier ordre à coefficients réels, Théorèmes d’approximation, Application à l’équation de transport, Ann. Sci. Ec. Norm. Sup. Ser. 4, 3, 185–233, 1970.
Berselli, L. C., Iliescu, T., and Layton, W. J.: Mathematics of Large Eddy Simulation of Turbulent Flows, Springer, Berlin, Heidelberg, ISBN 10 3642065791,, 2005.

Birnir, B., Smith, T. R., and Merchant, G. E.: The scaling of fluvial landscapes, Comput. Geosci., 27, 1189-1216, https://doi.org/10.1016/S0098-3004(01)00022-X, 2001.

Bonetti, S., Bragg, A. D., and Porporato, A.: On the theory of drainage area for regular and non-regular points, P. Roy. Soc. A, 474, 20170693, https://doi.org/10.1098/rspa.2017.0693, 2018.

Bonetti, S., Hooshyar, M., Camporeale, C., and Porporato, A.: Channelization cascade in landscape evolution, P. Natl. Acad. Sci. USA, 117, 1375-1382, 2020.

Cheskidov, A., Olson, E., Holm, D., and Titi, E.: On a Leray-$\alpha$ model of turbulence, P. Roy. Soc. Lond. A, 146, 1-21, 2005.

Coatléven, J.: Some multiple flow direction algorithms for overland flow on general meshes, ESAIM: Mathematical Modelling and Numerical Analysis, 54, 1917-1949, 2020.

Davy, P. and Lague, D.: Fluvial erosion/transport equation of landscape evolution models revisited, J. Geophys. Res.-Earth, 114, F03007, https://doi.org/10.1029/2008JF001146, 2009.

Desmet, P. J. J. and Govers, G.: Comparison of routing algorithms for digital elevation models and their implication for predicting ephemeral gullies, Int. J. Geo. Inf. Syst., 10, 311-331, 1996.

DiPerna, R. and Lions, P.: Ordinary differential equations, transport theory and Sobolev spaces, Invent. Math. 98, 511-547, 1989.

Eymard, R., Gallouët, T., Gervais, V., and Masson, R.: Existence and uniqueness of a weak solution to a stratigraphic model, edited by: Feistauer, M., Dolejši, V., Knobloch, P., and Najzar, K., Springer, Berlin, 278-287, https://doi.org/10.1007/978-3-642-18775-9 25, 2004.

Eymard, R., Gallouët, T., Gervais, V., and Masson, R.: Convergence of a numerical scheme for stratigraphic modeling, SIAM J. Numer. Anal., 43, 474-501, 2005.

Fernández-Cara, E., Guillén, F., and Ortega, R.: Mathematical modeling and analysis of visco-elastic fluids of the Oldroyd kind, in: Numerical Methods for Fluids, Part 2, Handbook of Numerical Analysis, vol. VIII, edited by: Ciarlet, P. G. and Lions, J. L., North-Holland, Amsterdam, Elsevier, 543-661, https://doi.org/10.1016/S1570-8659(02)08005-5, 2002.

Freeman, T. G.: Drainage with divergent flow over a regular grid, in: Proc. 8th Biennial Conf. Simulation Society of Australia, Canberra, 160-165, 1989.

Freeman, T. G.: Calculating catchment area with divergent flow based on a regular grid, Comput. Geosci., 17, 413-422, 1991.

Gallant, J. C. and Hutchinson, M. F.: A differential equation for specific catchment area, Water Resour. Res., 47, W05535, https://doi.org/10.1029/2009WR008540, 2011.

Gilbert, G.: Geology of the Henry Mountains, US Geological and Geological Survey, Washington, D.C., https://doi.org/10.3133/70038096, 1880.

Girault, V. and Tartar, L.: $L^p$ and $W^{1,p}$ regularity of the solution of a steady transport equation, C. R. Acad. Sci. Paris Ser. I, 348, 885-890, 2010.

Graf, W. H. and Altinakar, M. S.: Hydraulique fluviale: Ecoulement et phénomènes de transport dans les canaux à géométrie simple, Traité de génie civil, vol. 16, Presses polytechniques et universitaires romandes, ISBN-10 2880748127,, 2000.

Granjeon, D.: Modélisation stratigraphique déterministe: Conception et applications d'un modèle diffusif 3-d multilithologique, PhD thesis, Université de Rennes I, 1996.

Grospellier, G. and Lelandais, B.: The Arcane development framework, POOSC 09: Proceedings of the 8th workshop on Parallel/High-Performance Object-Oriented Scientific Computing, 1-11, https://doi.org/10.1145/1595655.1595659, 2009.

Guermond, J.-L., Oden, J., and Prudhomme, S.: An interpretation of the Navier-Stokes-$\alpha$ model as a frame-indifferent Leray regularization, Physica D, 177, 23-30, 2003.

Holmgren, P.: Multiple flow direction algorithms for runoff modelling in grid based elevation models: an empirical evaluation, Hydrol. Process., 8, 327-334, 1994.

Hooshyar, M. and Porporato, A.: Mean dynamics and elevation-contributing area covariance in landscape evolution models, Water Resour. Res., 57, e2021WR029727, https://doi.org/10.1029/2021WR029727, 2021a.

Hooshyar, M. and Porporato, A.: Spectral signature of landscape channelization, Geophys. Res. Lett., 48, e2020GL091015, https://doi.org/10.1029/2020GL091015, 2021b.

Hooshyar, M., Bonetti, S., Singh, A., Foufoula-Georgioui, E., and Porporato, A.: From turbulence to landscapes: Logarithmic mean profiles in bounded complex systems, Phys. Rev. E, 102, 033107, https://doi.org/10.1103/PhysRevE.102.033107, 2020.

Leopold, L. B., Wolman, M. G., and Miller, J. P.: Fluvial Processes in Geomorphology, W. H. Freeman, San Francisco, California, ISBN-10 0486845524,, 1964.

Leray, J.: Sur le mouvement d'un fluide visqueux emplissant l'espace, Acta Math., 63, 193-248, 1934.

Maxwell, J. C.: On hills and dales, Philos. Mag. J. Sci., 4/40, 421-427, 1870.

Pelletier, J.: 2.3 Fundamental Principles and Techniques of Landscape Evolution Modeling, in: Treatise on Geomorphology, edited by: Schroder, J. F., Academic Press, San Diego, 29-43, https://doi.org/10.1016/B978-0-12-818234-5.01001-4, 2013.

Pelletier, J. D.: Minimizing the grid-resolution dependence of flow-routing algorithms for geomorphic applications, Geomorphology, 122, 91-98, 2010.

Perlin, K.: An image synthesizer, ACM SIGGRAPH Comput. Graph., 19, 287-296, https://doi.org/10.1145/325165.325247, 1985.

Perron, J. T., Dietrich, W. E., and Kirchner, J. W.: Controls on the spacing of first-order valleys, J. Geophys. Res., 113, F04016, https://doi.org/10.1029/2007JF000977, 2008.

Perron, J. T., Kirchner, J. W., and Dietrich, W. E.: Formation of evenly spaced ridges and valleys, Nature, 460, 502-505, 2009.

Peton, N., Cancès, C., Granjeon, D., Tran, Q.-H., and Wolf, S.: Numerical scheme for a water flow-driven forward stratigraphic model, Comput. Geosci., 24, 37-60, 2020.

Porporato, A.: Hydrology without dimensions, Hydrol. Earth Syst. Sci., 26, 355-374, https://doi.org/10.5194/hess-26-355-2022, 2022.

Scheingross, J. S., Limaye, A. B., McCoy, S. W., and Whittaker, A. C.: The shaping of erosional landscapes by internal dynamics, Nat. Rev. Earth Environ., 1, 661-676, 2020.

Smith, T. R. and Bretherton, F. P.: Stability and the Conservation of Mass in Drainage Basin Evolution, Water Resour. Res., 8, W03417, https://doi.org/10.1029/WR008i006p01506, 1972.
Smith, T. R., Birnir, B., and Merchant, G. E.: Towards an elementary theory of drainage basin evolution: I The theoretical basis, Comput. Geosci., 23, 811–822, 1997.

Veiga, H. B. D.: Existence results in Sobolev spaces for a stationary transport equation, Ricerche Mat. Suppl., XXXVI, 173–184, 1987.

Zhiyin, Y.: Large-eddy simulation: Past, present and the future, Chinese J. Aeronaut., 28, 11–24, 2015.