On the main components of landscape evolution modelling of river systems

Currently, the use of numerical models for reproducing the evolution of river systems and landscapes is part of the day-by-day research activities of fluvial engineers and geomorphologists. However, despite landscape evolution modelling is based on a rather long tradition, and scientists and practitioners are studying how to schematize the processes involved in the evolution of a landscape since decades, there is still the need for improving the knowledge of the physical mechanisms and their numerical coding. Updating past review papers, the present work focuses on the first aspect, discussing six main components of a landscape evolution model, namely continuity of mass, hillslope processes, water flow, erosion and sediment transport, soil properties, vegetation dynamics. The more common schematizations are discussed in a plain language, pointing out the current knowledge and possible open questions to be addressed in the future, towards an improvement of the reliability of such kind of models in describing the evolution of fluvial landscapes and river networks.


Introduction
The present paper reviews the state of the art about modelling of landscape evolution, with a particular focus on the main components typically schematized in numerical codes that can be applied for modelling fluvial terrains shaped by the interaction of internal and external processes, such as precipitation, water flow and sediments, as well as vegetation and soil properties. Nowadays, often the term "model" refers to both the underlying theory and the computer programs used to calculate approximate solutions of the governing equations, involving possible misunderstandings. Numerical models represent an indispensable tool for assisting geomorphologists in reproducing the origins and dynamics of surface landscapes, combining a quantitative characterization of terrain with various theories describing the modification of river system topography by the variety of processes that sculpt it (Mark 1975;Tucker and Hancock 2010;Baas 2017).
In the last decades, the ability of engineers and geomorphologists to measure the topography of river beds and hillslopes has grown tremendously, moving from topography maps, very imprecisely and requiring a massive work to be revised, towards digital elevation models and digital maps, generally having a higher resolution and covering the majority of the emerged landmasses (Gesch et al. 2006), and more easily updatable. In addition, the recent development of high-resolution mapping tools like laser scanners and cameras, as well as satellites, assured a detailed and reliable description of the changes of the Earth's surface, also in regions where the access could be more complicated. At the same time, theories and models of landscape evolution have grown, accounting for more processes and for a more sophisticated description of them. As recently as the second half of the last century, a landscape evolution model was intended as the sequential evolution of a landscape over the geological time. By the end of the century, this term had been associated with a more scientific meaning: a mathematical theory describing how the actions of a multitude of geomorphic processes interact during the time to shape the basin topography.
Aside from their physical meaning, the complexity of the governing equations of landscape evolution requires a numerical solution method to be solved in a closed form.
The growing complexity introduced in landscape evolution models, accompanied by the advances in computing techniques and acquisition of topographic data, has revolutionized the ability of geomorphologists and fluvial engineers to measure and model landforms and their rate of change, as well as to investigate and numerically reproduce how such forms and dynamics arise from the physics of geomorphic processes (Coulthard 2001;Chen et al. 2014;Willgoose 2018).
In geophysics, primary drivers are generally described by means of a system of partial differential equations, which originate from the classical mechanical theory (e.g. equations describing water and sediment transport, rock mechanics, heat transfer, river hydraulics). The scale dependence of such processes complicates the analysis of feedbacks and interactions between them (Bracken et al. 2015). Therefore, depending on the representation scale, models can be computationally intensive and thus typically limited to a small spatial scale of few metres if strictly adherent to the theory, or addressing large scales problems like the evolution of a whole river basin (10 2 -10 5 km 2 ) with a reduced computational effort, but introducing simplified equations (Khosronejad et al. 2014;Larsen et al. 2016). For addressing problems acting at the landscape scale and therefore involving a great number of coupled state variables, a variety of approaches and methodologies is nowadays available. Currently, in fact, landscape evolution models can combine hillslope, channel, tectonic and vegetation processes by linking physically based equations, which represent simplifications of the real world (e.g. De St. Venant equations for the water flow, geomorphic transport laws, erosion/deposition and sediment transport equations), with semiempirical approaches (e.g. organic accumulation, vegetation growth, the presence of external drivers like fire or wind). In addition, these models are able to couple stochastic (e.g. probabilities of sediment detachment, vegetation encroachment and density, distribution of precipitation) and deterministic (e.g. water flow velocities, bank failure) dynamics (Willgoose 2005;Fagherazzi et al. 2012).
Typically, many of the most used landscape evolution models are run over large domains and are therefore computationally intensive (Salles 2016). Consequently, geoscientists are facing the challenge of reducing the computational effort to a minimal level for performing more effectively. As an example, Stark and Passalacqua (2014) developed a simplified landscape evolution model as a set of coupled dynamic systems, aiming to evaluate the changes of biomass and regolith under mass wasting and run-off erosion. An even more simplified approach was proposed by Franzoia and Nones (2017) and tested by Nones et al. (2019), who described the very long-term evolution of a river watershed by applying, at the watershed scale, a 0-D lumped hydromorphological model. The popular, alternative strategy of cellular automata modelling involves describing the physics governing fluid flow or sediment transport by means of discrete rules that control water, air and sediment transport processes on the basis of information from surrounding model grid cells (Willgoose et al. 1991;Liu and Coulthard 2017). This cellular automata strategies have made it possible to reproduce shallow-water flows for hydrological purposes Caviedes-Voullième et al. 2018), but also to simulate the development of braided streams (Murray and Paola 2003), floodplains (Coulthard and Van De Wiel 2006), sand dunes (Zhang et al. 2012), wetland landscape pattern (Williams et al. 2016) and river deltas , for evaluating their response to global changes and human drivers. Aside from cellular automata, precipiton methods can be applied for simulating the evolution of a river landscape, given their ability in mimicking self-organized emerging properties of geomorphological systems, from high-resolution braided patterns to drainage network organization (Davy et al. 2017). Even if based on dissimilar approaches, these two methods have multiple similarities, like the computational effort needed for simulating large domains or the complexity correlated with solving the hydrodynamics in detail.

The rationale of the research
In the 1950s and 1960s, the discipline of geomorphology turned from a qualitative approach, widely applied at the beginning of the XX century, towards a more quantitative analysis of landscape evolution (Kamp and Owen 2013). Starting from the 1970s, there have been a growing number of excellent review papers covering landscape evolution models and various aspects of geomorphic modelling, offering a very wide vision on the field (see, among many others, Carson and Kirkby 1972;Kirkby 1996;Coulthard 2001;Martin and Church 2004;Khosronejad et al. 2014) and proposing several open questions to be addressed in representing the evolution of the Earth's surface through a mathematical model. As stated by Kamp and Owen (2013), mathematical models were and are still paramount in understanding the landscape processes and the feedback between the involved components, and the recent technological development contributed in increasing their use worldwide.
The aim of the present paper is to discuss six main components usually considered in modelling the evolution of river systems, where the major part of sediments generated on the hillslopes is transported through the drainage network by the water flow. On the one part, river basins are the fundamental geomorphic unit (Chorley 1969), and fluvial landscapes cover most of the Earth's surface. On the other part, being one of the most human-impacted environments, they need particular attention and their future evolution should be addressed with proper methodologies that also account for climate change. For better focusing on river basins, are here not considered additional erosional systems like eolian landscapes, karstified terrains, ocean floors and glacial landscapes. The application of landscape evolution models to river basins can provide additional insights on the physicochemical processes that interact to shape the surface of a fluvial system, transferring the mass from one area to another. Moreover, the opportunity to have a graphical representation of the basin evolution enhances the ability of scientists and non-experts to interpret possible changes of the surface and to quantify the consequences of various hypotheses about fluvial dynamics.
In addition, landscape evolution models can also be applied for evaluating the development of specific features like passive margins (Tucker and Slingerland 1994;Ruetenik et al. 2016;Braun 2018) or mountain chains (Miller and Slingerland 2006). Numerical models can be focused on the long-term evolution of landscape (Bishop 2007) and river systems (Coulthard and Macklin 2001;Di Silvio and Nones 2014;Varrani et al. 2019) or on evaluating tectonics (Beaumont et al. 2000) and other surface processes. Under an environmental point of view, landform evolution models can be used as a methodology for evaluating and managing degraded landscapes such as abandoned mines (Hancock et al. 2000;Hancock and Willgoose 2018) and contaminated sites (Evans 2000), or for projects involving landscapes affected by disturbance of soil and/or vegetation (Coulthard et al. 2002). Because these models allow to visually evaluate the temporal changes of the basin in terms of elevation, catchment size and shape, they can also be applied to support the study of dynamic phenomena like gully network development and valley alluviation or river avulsion, which is generally not possible with fixed-terrain models based on the classic USLE approach (Karydas et al. 2014). In addition, the flexibility of these models permits to evaluate and potentially combine simulations of processes acting at different spatiotemporal scales, spanning from short-term soil loss along single hillslopes (Montgomery and Dietrich 1992;Hancock et al. 2008) to catchment-scale assessments over geologic time (Hancock et al. 2015;Varrani et al. 2019), eventually coupling geomorphic and tectonic models (Beaumont et al. 2000) or accounting for the mobility of the whole river network (Whipple et al. 2017).
Regarding the mathematical description of erosion and transport rate adopted in landscape models, in the past, Carson and Kirkby (1972) discussed a variety of 1-D models of hillslope evolution under different geomorphic scenarios, while Dietrich et al. (2003) provided an overview of rate laws for both hillslope and channel processes. In their reviews on landscape evolution modelling, Coulthard (2001), Chen et al. (2014) and Willgoose (2018) provided a perspective on their strengths and weaknesses, showing that several solutions for reproducing the evolution of terrestrial landscapes exist. Recently, a thorough review on the fundamental equations implied in landscape evolution models was provided by Chen et al. (2014), who pointed out that the numerical implementation is a non-trivial problem, particularly in simulating water flow and sediment transport in an efficient and highly accurate way.
In the past, many contributions were focused on discussing general and philosophical issues relevant to geomorphic modelling. As an example, Carson and Kirkby (1972) firstly and then Kirkby (1996) pointed out the theoretical foundations of the modelling approach adopted in several fluvial landscape evolution models, showing their role in reducing the gap between the theory and the experimental approach. Focusing not only on landscape evolution models but also on the numerical modelling in general, Oreskes et al. (1994) provided a perspective on the codes' structure and problems associated with their verification and validation. Mimicking this approach, other works (Martin and Church 2004;Larsen et al. 2016) showed the limitation of adopting numerical models for describing the nature complexity proposing open questions to be addressed in the future by means of new methodologies.
Moving from the important review proposed by Tucker and Hancock (2010) and maintaining a similar structure because of its effectiveness in driving the message, the present paper summarized their conclusions on the structure and the constitutive equations of landscape evolution models in a relatively simplified way, using a plain language to provide the readers with a general overview rather than with a complex mathematical description of the phenomena involved, as proposed by Chen et al. (2014). Indeed, the review is mostly devoted to students and young researchers that want to understand the basic mechanisms of landscape evolution models, aiming to design their own approach on the problem.
In addition to what was done in previous reviews, this work discusses the importance of soil processes and vegetation dynamics in shaping fluvial landscapes, and the need of a deeper understanding of the feedback between all these processes for adequately implementing them in a numerical code. The next section is divided into six subsections, and each component is analysed individually. In detail, being the basic component of a landscape evolution model, the conservation of mass is firstly reviewed, and then, the hillslope processes are addressed. Two of the main components in landscape evolution models are the water flow and the sediment behaviour (erosion/deposition and transport) and are discussed in the following two subsections. Aside from these four components, which are included in landscape modelling since the very beginning of this research field, for the future, a major effort should be dedicated to evaluating the effects of soil properties and vegetation dynamics in changing the Earth's surface. Final conclusions and open questions for researchers and scholars are summarized at the end of the manuscript, aiming to provide the readers with possible paths to follow in developing landscape evolution models more capable to reproduce the variety and feedback of mechanisms interrelating in nature.

Simulating the evolution of a river system via a landscape evolution model
To simulate the creation and evolution of a river system, landscape evolution models can be applied. Based on geometrical, morphometric, hydrological and additional (e.g. wind, vegetation, snow, ice, fire, herbivorous) inputs, a landscape evolution model combines such quantities for simulating the future changes of the Earth's surface. These models are based on a system of equations that schematizes the mass continuity, specific geomorphic transport functions for describing the generation and movement of sediments and solutes on the basin hillslopes, as well as to reproduce the erosion phenomenon, the water flow and the transport of water-sediment mixtures along the river network . Depending on the modellers' needs and the structure of the code, a series of numerical methods can be adopted for discretizing the solution in space and time, aiming to obtain more or less approximated solutions of the governing equations.
For adequately reproduce the evolution of a fluvial landscape, the single components should be effectively characterized, as well as the feedback between them. In the following subsections, six main components are described, summarizing the results already available in the literature and pointing out the opportunity to focus the future research on some of them (i.e. vegetation dynamics, soil properties, sediment transport), given that classical mechanisms (i.e. mass continuity, hillslope properties, water flow) are studied since decades (Mark 1975;Tucker and Hancock 2010).

Exchange of mass
Geomorphic systems are rather basic systems, where the mass is conserved absolutely. However, in the case of landscape evolution, there are examples where the mass is not conserved because of sediment detachment, such as models based on stream power erosion formulas (Warren et al. 2019). Moreover, the mass of water is not always conserved, because it also depends on the flow routing assumptions typical of each model. Therefore, in studying the evolution of fluvial landscapes, there are many possible frameworks for addressing the continuity of the mass, depending on the kind of process reproduced, the circumstances under study and the numerical scheme adopted. Given that each of these possible approaches has its own assumptions and limitations, with their pros and cons, system modelling can be considered, to some degree, as a subjective research field. In fact, the included components, the adopted methodology and technology, as well as the modelling schematization, are arbitrarily assumed by the researcher depending on the case study and the research aims, as well as the modeller experience.
In general, the rate of change in a given control volume V can be derived by comparing the mass rate entering the volume with the one going out (Fig. 1). In other words, the process (rate of change) can be computed as a result of the nature and geometry of the idealized model (mass rate difference in-out).
One of the most common continuity expressions in geomorphic models assumes that the control volume can be schematized by means of a vertical column of rock and/ or soil. Starting from this general theory, the modeller can introduce multiple simplifications related to the density or the thickness of the considered material, as well as on its porosity and grain size. Several landscape evolution models consider that all the surface is composed by rock or by a combination of rock and sand. In the first case, assuming that the surface height is a single-valued function of the horizontal position involves the impossibility to describe vertical faces and overhangs. In addition, changes in height due to compaction/expansion of the underlying soil are generally ignored, as well as variations in the thickness of the soil layer. This latter hypothesis means that effects like the dependence of sediment transport rate on the regolith thickness (Carson and Kirkby 1972;Kirkby 1992;Braun et al. 2001;Skinner et al. 2018), the feedbacks between soil water storage capacity, the run-off generation and its effective rate or the weathering and sediment transport processes (Kirkby 1976;Saco et al. 2006;Dochez et al. 2014) are ignored. On the other part, assuming that a contact between the loose, mobile regolith and the underlying rock exists, provides a slightly more complete approach, but still many  (Ahnert 1976(Ahnert , 1987Heimsath et al. 2001;Strudley et al. 2006).
Despite being largely adopted in landscape evolution modelling, the vertical-column continuity approach can sometimes result in being too much simplistic. In fact, there are several landforms that do not fit this framework, such as cliffs, waterfalls and gully headscarps, which have vertical or overhanging faces. Moreover, the 2-D continuity framework commonly implemented in landscape models is often very simple and not able to simulate vertical variations in weathering rates, shallow flows and rock properties. For overcoming these limitations, a fully 3-D approach is therefore necessary, subdividing the column vertically and introducing additional equations to describe the vertical variations of the soil properties (Kirkby 1976;Vanwalleghem et al. 2013). Other approaches included the vertical exchange of mass flux due to the soil strain and the advection of soil layers towards (away from) the eroding (aggrading) land surface, or the consideration of additional drivers like gravitational compaction or changes in mineralogy and geochemistry (Fritsch et al. 2011).

Hillslope erosion processes
Landscape evolution models are based on geomorphic transport functions, which, usually, have a rather general formulation but adopt site-specific parameters, developed for a specific scope (Chen et al. 2014). Moreover, classical approaches to hillslope erosion processes such as the RUSLE (Renard et al. 1997;Nasir and Selvakumar 2018) or the WEPP (Flanagan et al 2012;Nearing et al. 2017) aggregate various geomorphic processes within a given area. (Namely, they can be considered as lumped models.) However, for better addressing the physics of landscape evolution, geomorphic process regimes like sheet flow and gullies should be treated separately (Momm et al. 2018). In the following, a brief overview of the main hillslope erosion processes is reported, showing the need for addressing every single component separately.
Based on long-lasting research, Nearing et al. (2017) defined the critical zone as the region between the top of the forest canopy and the base of the weathering horizon. Processes acting in such zone weaken the rock via mechanisms like mechanical wedging, fracturing, chemical alteration, biological disruption, etc. Although the weathering processes are well studied under a qualitative point of view, their mathematical representation is a rather new field of research, which is mainly focused on predicting rates of change and patterns of rock disintegration induced by specific chemical and physical processes (Cohen et al. 2009;Murphy et al. 2016). Moving from the original intuition proposed by Gilbert (1877), later revised and analytically expressed by Ahnert (1976), one can observe that, assuming a constant rate of regolith production from bare bedrock and a fixed characteristic decay length scale, under quasi-steady conditions (i.e. the regolith thickness varies very slowly with respect to the surface erosion), an inverse relationship between thickness and erosion rate exists. This relationship has been verified against several field data, resulting in being consistent in different environments, ranging from semi-arid and coastal environments to high alpine terrains (Ahnert 1987;McKean et al. 1993). The mathematical formulation of this process can imply several assumptions, but, generally, the obtained decay curve has an exponential trend (Anderson 2002;Saco et al. 2006;Strudley et al. 2006). In the last century, Kirkby (1985) proposed an alternative to the exponential decay models: instead of assuming a sharp contact between bedrock and regolith, he described the transition from rock to regolith as a gradational process having the deficit of soil as a state variable. Such a deficit can be represented by the fraction of unweathered rock remaining at a certain level in the soil profile. He demonstrated that this schematization can be more appropriated in effectively describing a wide interface region between unaltered material and fully weathered soil. Indeed, even if successfully verified in several case studies, the exponential decay rules cannot be considered as a definitive solution for describing the observed regolith thickness patterns, but rather just a relatively simple and appealing method to be used until better approaches will be available. In this sense, further research on the physical mechanics and chemistry of weathering processes is needed for obtaining a better mathematical relationship to correlate rock disintegration rates to factors like subsurface temperature, stresses and mineral alteration (Anderson 1998;Fletcher et al. 2006).
As described in detail in the past (Carson and Kirkby 1972;Tucker and Hancock 2010), to reproduce the longterm phenomenon of soil creep on low-gradient basins, a linear slope-dependent transport function can be adopted (Fernandes and Dietrich 1997), accounting for the convexupward hillslope profiles. Even if widely and successfully applied in many studies since decades, including fault scarps and fluvial systems, as well as marine and lake-shore terraces (Avouac 1993;Arrowsmith and Rhodes 1994;Pelletier et al. 2006), and calibrated against field data derived from cosmogenic nuclide mass-balance measurements (McKean et al. 1993;Heimsath and Ehlers 2005), the calibration constant is still the main source of uncertainty. In fact, an estimate of its magnitude can be obtained from a variety of approaches and for specific processes (Kirkby 1971;Black and Montgomery 1991;Anderson 2002), but, under a general point of view, it should be treated as an empirical parameter. Although a linear relationship, with a constant parameter, provides reliable results, physical considerations suggest that the regolith thickness can influence the soil creep equation. Therefore, many depth-dependent creep functions have been suggested in the literature (Ahnert 1976(Ahnert , 1987Strudley et al. 2006), spanning from very simplified (Rosenbloom and Anderson 1994) to more sophisticated (Braun et al. 2001;Anderson 2002) approaches. However, despite the rising of specific studies, it is evident that more research is needed, especially considering landscapes having steeper gradients, close to the angle of repose for natural soils . A few nonlinear relationships were proposed (Ferdowsi et al. 2018), but, generally, being developed on specific datasets, such formulas result too site-related and therefore not applicable in different contexts. For the future, additional studies should be performed, also involving new technologies like radiometry, laser scanning and micromorphology (Pawlik and Šamonil 2018).
Aside from fluvial processes, transport functions adopted for describing mass movements like shallow, rapid landsliding are more problematic. In fact, there are two general approaches that can be adopted for describing these phenomena: flux-based and event-based models. The first type approximates a series of events in terms of the long-term average rate of mass transfer by means of a transport function (Kirkby 1987). These models are capable to describe the time-averaged sediment transport at a time-scale relevant to landform evolution, but only at the very local spatial scale: the flux in a specific point is represented as a function of the local variables, neglecting the effects of the surroundings and the significance of long-distance transport events on steep slopes (Howard 1994). In the last decades, multiple attempts were made for incorporating the long-distance transport effects into landscape evolution models by accounting for the expected flow paths. However, many uncertainties are still evident, mainly because of the probabilistic behaviour of the sediments and the influence of the soil properties. Indeed, the evident connection between transport statistics, topography and morphological evolution suggests the use of event-based models, but, even if they are more physically rooted being grounded on the current knowledge of landslide triggering and motion (Tucker and Bradley 2010), at the same time, they are computationally not really efficient and therefore not widely applied and tested.

Run-off processes
As well known, the transport of sediment by the flowing water is a fundamental feature of the Earth's processes (Lorang and Hauer 2017). Indeed, geomorphic works in a drainage basin are mostly correlated with surface water, and, therefore, knowing how the water flow is handled in landscape evolution modelling represents a central issue. Despite the possible spatial discretization methods that can be adopted, a common feature between the various models is the need to combine short-time scales (minutes to seasons) associated with hydrologic processes with much longer timescales (years to centuries) that are related to sediment transport and landform evolution.
Typically, in a 2-D model, the flow field is described by means of the De St. Venant (shallow-water) equations, which represent the vertically integrated form of the Navier-Stokes equations for incompressible, free surface flow. They contain a description of the continuity of mass and momentum in the two horizontal dimensions and a friction function to describe the relationship between flow velocity and bed resistance, accounting for four main forces: inertia, gravity, fluid pressure and boundary friction. Being these equations highly complex and generally not analytically solvable, simplified numerical solutions should be applied, accounting for several limitations, as pointed out in the following.
Typically, because channelized flows accelerate only slowly in space (considering a reach-wise averaged velocity), the gradually varied flow approximation is introduced, assuming that the inertial terms in the momentum equation can be neglected. In addition, dropping the time derivative yields the diffusion wave approximation, which is valid in the case of flows mainly driven by pressure and gravity gradients. In the case of small changes of the flow depth in the stream-wise direction with respect to the bed morphology, the gravity represents the main driver of the flow, and the pressure-gradient term can be neglected for obtaining the kinematic-wave equations. In this case, the water gravityrelated acceleration is everywhere balanced by the friction. For gravity-driven (kinematic) flows, the local bed shear stress can be represented as a function of the fluid density (water and sediment mixture), the local water/sediment discharge and the bed slope and friction. In a 2-D schematization, this approximation means that the flow lines follow the surface topography.
Based on this approach, many landscape evolution models use a cellular routing algorithm, imposing that the water flows from a cell to the adjacent one, following the steepest descent. As one can easily figure out, cellular routing algorithms are closely correlated with the spatial discretization of the domain. Indeed, in a numerical model, the continuous landscape surface is typically represented by discrete elements, which can be square cells, leading to pretty simple finite-difference solutions. However, in some cases, such square cells are not flexible enough for representing the computation domains in a proper manner. Therefore, to account for more complex domains, triangular elements associated with a finite-element solution (Maniatis et al. 2009) or triangular irregular (unstructured) cells having the nodes connected using a Delaunay triangulation and the surface nodes area described via Voronoi or Thiessen polygons can be adopted, such as made in common landscape models like the CASCADE and CHILD codes (Forte et al. 2016). Despite Caviedes-Voullième et al. (2012) demonstrated the utility of using triangular unstructured meshes for keeping a low computational cost while ensuring accuracy, only a few models allow for using such a structure (Costabile et al. 2017). The advantages of a cell-routing approach are the simplicity and the speed, but many drawbacks are also present. Firstly, it is hard to handle diverging flows, which is typical of complex river systems. Moreover, the kinematic convergence of the flow depends on the width, which can be supposed equal to the width of the cells, leading to a gridsize dependence of water depth and current velocity (Willgoose et al. 1991). Alternatively, the flow can be assumed as confined to sub-grid cell features, which should have a predetermined (empirical) width (Howard 1994;Tucker and Slingerland 1994). A few landscape models have a less strict approach, relaxing the single-flow-direction assumption by introducing explicit numerical solutions of the steady 2-D kinematic-wave equations (Morgan 1980) or using multiple directions algorithms, which assume that the flow going out from a cell is split among the downslope neighbours, weighted according to the gradient in each direction. Given that this latter type of algorithms provides a better description of shallow overland flow (sheet flow as described by Morgan in 1980) on convex hillslopes and fans (Pelletier 2004), its use is fast spreading. In fact, common and wellknown codes like CAESAR (Van De Wiel et al. 2007) or recent examples such as the Landlab  use multiple flow-direction algorithms, accounting for all the possible directions of flow propagation. These models provide an effective way to approximate time-varying, 2-D flow fields without the computational effort required by the traditional solution of the shallow-water equations, which could be very high especially for large domains (Garcia-Navarro 2016; Shustikova et al. 2019). Commonly, cellbased or kinematic-wave water routing is associated with a steady flow. As an example, the SIBERIA and the DELIM models compute the water discharge as a power function of the watershed area, assuming the local equilibrium between rainfall and run-off (Willgoose et al. 1991;Howard 1994). However, as Sólyom and Tucker (2004) demonstrated, the local geomorphology can highly affect the hydrological behaviour. Moreover, many studies pointed out the importance of the spatial variability of run-off generation, finding out that the run-off excess (saturation) tends to enhance both the hillslope convexity and the hillslope-channel transitions in equilibrium landscapes (Ijjász-Vásquez et al. 1992).
The kinematic-wave theory represents a reliable approximation for channelized flows, but some problems arise in describing the 2-D evolution of landscapes. Indeed, errors can be hindered behind a series of questionable assumptions that lead to the right solutions (Izumi and Parker 2000). On the one part, the problem of flow convergence along valley axes can represent an obstacle to properly capture the transition from distributed to channelized flow, which can be somehow handled only posing major attention on the spatial resolution of the model (Kirkby 1994;Perron et al. 2008). In fact, to overcome such problems, fine-detailed models using the diffusive wave theory can be developed, but there is still the need of employing more powerful computers for evaluating the long-term evolution of large areas, which can hinder their application to real-time forecasts.
Obviously, there are differences in terms of timescale in simulating the run-off which could be observed during a storm event or the long-term evolution of a river watershed. The majority of landscape evolution models deal with this aspect imposing a geomorphically effective run-off event to describe the basin erosion. Namely, such codes assume a single, steady run-off coefficient is equivalent, in terms of geomorphic effectiveness, to a series of run-off events. There are many examples dealing with this approach, spanning from imposing a relationship between a time-averaged sediment transport discharge and the average water flow peak discharge to more complex approaches (Willgoose et al. 1989). However, all these methods assume that the event is somehow stable, while many researchers pointed out the need to consider the role of the discharge variability in time (Lague et al. 2005;Huang and Niemann 2006;Molnar et al. 2006). Regardless of the detail of each method, they commonly agreed that erosion and transport rates increase with the temporal increment in discharge fluctuations, because they depend more than linearly on the water discharge. While there are several examples for which such effective event assumption is reasonable, recent studies proved that the time variability in hydrologic forcing can have a great impact on the landscape dynamics and, therefore, should be incorporated in the landscape evolution modelling, possibly through a stochastic description of both the rainfall and the run-off events (Tucker and Bras 2000;Whipple and Tucker 2002;Armitage et al. 2018).
In the future, many challenges related to modelling the feedback effects between a changing climate, hydrology and landscape evolution in a coupled way should be faced, aiming to account for different spatiotemporal scales and overcome the simplifications generally applied in practice (Sólyom and Tucker 2004;Huang and Niemann 2006;Anders et al. 2008). Moreover, the randomness in the temporal dynamics of run-off processes requires the development of new high-flow statistics for better describing the evolution of landscapes like river floodplains, which are more impacted by extreme flows.

Sediment transport from the hillslopes to the river system
In shaping a river channel, the water flow erodes the bed with a rate limited by the detachment of particles (supply-limited systems) or by the capacity of the flow to transport sediment particles (capacity-limited systems), with a multitude of intermediate behaviours (Carson and Kirkby 1972;Hajigholizadeh et al. 2018;Shobe et al. 2018). Thus, each system needs a different schematization, and the complexity varies depending also on the erosion rate: supply-limited systems result in being the simplest in terms of numerical modelling. In fact, in such systems, the sediment particles disappear as soon as they are eroded (Bagnold 1966;Howard 1994). Therefore, in this case, the erosion rate is assumed as a function of the bed shear stress, which, in its turn, depends on the local slope and discharge, giving rise to the so-called stream power erosion law (Bagnold 1966;Howard and Kerby 1983;Warren et al. 2019). A key property of these systems is the wave-like nature: there is a tendency to form erosional fronts that propagate upstream . In the case of capacity-limited river systems, the erosion rate is a function of the unbalance between sediments entering and going out from the system, assuming a local morphological equilibrium, where the transport rate is everywhere equal to the local carrying capacity. The capacity-based approach is most applicable for describing bedload transport in gravelbed rivers, given that coarser particles have shorter travel distances, so the assumption of immediate adaptation of the transport rate to changes of water discharge or slope is quite reasonable (Einstein 1950). On the other part, in the case of systems mostly driven by the suspended load like sandy rivers, this approach tends to fail because it essentially ignores the time required for sediment grains to settle in the water column in response to transient hydrology. Indeed, the mechanism requires to define an equation representing the mass continuity for sediments in the water column, as well as detachment and settling functions, which are generally correlated with the local shear stress and grain size (Bracken et al. 2015). There are many formulas that can be adopted to describe the erosion and sediment transport phenomena in a river system, but, despite this, they perform in a very similar manner if looking at the long-term longitudinal river-profile evolution under steady conditions Varrani et al. 2019). However, many differences arise in applying the models in transient conditions (Attal et al. 2008;Franzoia and Nones 2017), suggesting the need for using natural experiments to test landscape models (Tucker 2009).
For effectively describing the natural environment and the formation of a river system, a landscape evolution model must correctly reproduce the transition dynamics from the hillslopes to the channel, and the degree to which the surface changes as a function of factors like relief elevation, local climate and river basin lithology (Kirkby 1987;Di Silvio and Nones 2014). The distinction between channels and hillslopes can be explicitly treated, but introducing hardly describable parameters (Willgoose et al. 1991), or representing the channels as sub-grid-scale features where the flow width is prescribed in an empirical way Tucker and Slingerland 1994). Depending on the problem under study, models can be built for representing large-scale mechanisms without requiring a very fine detail (Kooi and Beaumont 1994;Lindim et al. 2016) or to reproduce the evolution of smallscale landforms, implying a grid resolution that can be smaller than the channel width (Perron et al. 2008). For having the order of magnitude of the scales involved, one can consider a regime equation in its original form, correlating the river width with a power of the bankfull discharge (Leopold et al. 1964). One can easily notice that channels are typically some orders of magnitude smaller than the whole basin, meaning that they can be effectively handled as sub-grid features in landscape evolution models. The channel geometry is paramount in defining the volume of sediments available: the narrower is the channel, the more confined is the flow. Assuming that other parameters like the bed slope are constant, a narrow channel means an increase in the bed shear stress and the unit stream power, which translates in a bigger rate of sediment detachment and transported by the current (Lague 2014;Armitage et al. 2018).
The complexity of fluvial morphodynamics needs to be simplified for speeding up the computation, especially in the case of large river basins. However, in doing that, some models lose their physical meaning, imposing that the erosion can be directly computed from the total discharge rather than from the specific one (Willgoose et al. 1989;Kooi and Beaumont 1994) or assuming empirical regime equations (Howard 1994;Tucker and Whipple 2002) hardly verifiable. On the one part, the use of empirical scaling laws has the advantage to explicitly calculate the cross-sectional averaged shear stress and stream power and to permit the application of physically based erosion and transport functions that depend on such quantities. On the other part, relying on simple scaling laws for describing the channel geometry has some drawbacks like the application of an equilibrium assumption to describe nonequilibrium dynamics (Nones and Di Silvio 2016) or the impossibility to describe rivers affected by external forcing factors like tectonics or lithological discontinuities ). In the last years, many models have been developed for reproducing bedrock channel evolution (Stark 2006;Wobus et al. 2008;Langston and Tucker 2018) and changes in channel width (Attal et al. 2008;Nones and Di Silvio 2016), as well as debris and granular flows (Howard 1998;Stock and Dietrich 2006), but there is ample room for improving them towards a more reliable estimate of landscape evolution, accounting for physically based laws, as well as spatially and temporally variable functions, for better incorporating the geological and climatological variability.

Soil properties
For terrestrial life, the soil represents one of the most important substances, supporting both the life (Lin 2011) and being a medium for transport and storage of water and gases (Strahler and Strahler 2006). Indeed, hydrological and morphological processes are a function of the soil characteristics (Bryan 2000), but also depend on the ratio between soil and rock coverage (Poesen and Lavee 1994). Therefore, the understanding of the formation, global distribution and functional properties of the soil is paramount in catching the mechanisms driving the landscape dynamics.
Aiming to link the small scale of soil characteristics to the large scale of landscape evolution, in the past years, many statistical methods have been developed for determining and mapping different soil properties depending on other soil characteristics and the basin geomorphology (Behrens and Scholten 2006). However, one of the shortcomings of such an approach is the need of having a very large and detailed dataset of soil attributes, such as the particle size distribution or the amount of organic matter, which is needed for predicting hardly measurable soil properties like the water content. In fact, even if applicable at the small scale, analyses at large (basin-wide) scale require distributed samples, which can be prohibitively expensive and definitely time-consuming (Scull et al. 2003).
While the spatial mapping of soil properties is important, understanding the evolution of these properties and processes at the required scale is also fundamental. For quantifying such processes and predicting the time evolution of the soil characteristics, modellers can apply processbased models (Hoosbeek and Bryant 1992;Minasny et al. 2008;Schoorl and Veldkamp 2016). On the other hand, the most tested but out-of-date process-based models cannot be applied to large domains due to an excessive need of computational resources; therefore, new methods based on state-space matrix methodology were recently introduced (Cohen et al. 2009), also accounting for multiple soil layers (Welivitiya et al. 2016). These models are able to adequately predict the soil properties of an individual pixel, but failed in modelling the spatial interconnectivity between the various parts of the soil catena that result from transport-limited erosion and deposition. To correctly predict the temporal changes of the spatially distributed soil properties, since the end of the last century, many researchers tried to couple the soil profile evolution with the landform evolution. As an example, Minasny and McBratney (2001) modelled the influence of soil and weathering processes on landform evolution using a single layer of soil, while Vanwalleghem et al. (2013) developed the MILESD code, which accounts for four layers (the bottommost bedrock layer and three soil layers above it) for reproducing the evolution of landforms, with a particular focus on Australia. Recently, the soil evolution module adopted in MILESD has been modified adding additional layers (Temme and Vanwalleghem 2016). Indeed, limiting the description of the topsoil to only three layers can hinder the importance of soil characteristics like the particle size distribution, which can be an index of various soil attributes such as the soil moisture content (Schaap et al. 2001;Minasny et al. 2015). As a matter of fact, future landscape evolution models should consider soil characteristics like depth and water holding capacity explicitly, preferably via a physically based approach, because they constrain the grading and amount of material eroded across the river basin.

Vegetation dynamics
Both the soil formation and the establishment of vegetation are paramount in changing the hydrological fluxes by accommodating soil moisture and facilitating the formation of sub-surface flow paths, affecting the form and the magnitude of erosion, sediment transport (Ebel et al. 2007) and deposition (Molina et al. 2009). Therefore, numerical models devoted to predicting the evolution of landscape features are highly dependent on vegetation properties (Casadei et al. 2003;Collins et al. 2004;Istanbulluoglu and Bras 2005;Yetemen et al. 2010). In addition, plants convert solar energy into geomorphic forces with very significant impacts from the regional to the global scales (Phillips 2009). Amundson et al. (2015) summarized the importance of vegetation in landscape evolution in a few points: (1) if water is available, a world without plants would likely have little or no soil on hillslopes; (2) plants may control the soil thickness; (3) soil production rates may be very high with respect to outcrop erosion rates (around one order of magnitude); (4) given that the soil residence times are constrained within a broad window of nutrient sufficiency/optimization, environments characterized by high weathering and low denudation rates can suffer from a deficit of rock-derived elements; (5) local feedbacks between plants, nutrients and soil thickness are possible; (6) at the very long-term (millennia), the vegetation evolution can impact (and be impacted by) geomorphic conditions.
As for the soil, the inclusion of vegetation dynamics in landscape evolution models is related to the spatial scale. Indeed, for adequately simulating how soils and biota interact with climate and bedrock (Corenblit and Steiger 2009), and, consequently, for modulating the geomorphic response at the catchment scale, it is necessary to collect detailed data spatially distributed across the basin. To reduce the computational effort in creating such a dataset, generally spatial analyses of remotely sensed images and digital maps of elevation, geology, soils and vegetation in relation to the local climate and sediment yield are developed (Newton et al. 2009). Aside from creating new databases, the attention of geomorphologists and soil scientists is now focused on qualitatively and quantitatively understand the effects of vegetation on the landscape physical processes, aiming to provide a more reliable schematization of them to be included in numerical codes. In this sense, the present major challenge is how to explicitly and quantitatively account for the role of biota in the production of soil from bedrock and its transport downslope, investigating the combined evolution and feedback of soil, plants, hydrology and climate.

Conclusions and open questions
Using a plain language, the present review proposes a short summary about landscape evolution modelling and the main components of such codes, showing that, even if characterized by a quite long history, this research field is still very active and several improvements are forecasted for the future for answering to a series of open questions towards a more reliable representation of the Earth's surface (Willgoose 2018). In fact, the understanding of landscape dynamics requires a deeper knowledge of the recursive, multi-scale interactions among abiotic and biotic states and processes (Phillips and Van Dyke 2017).
As for the geomorphic transport functions, in the future, additional research should be performed towards a better evaluation of the dynamics associated with sediment characteristics such as a varying grain size distributions, the role of the basin lithology and the horizontal movement of geomorphic features due to processes like scarp retreat and tectonic displacement.
Most importantly, for obtaining a consistent schematization of the natural landscape, the dynamics of soils must be incorporated in new numerical models, overcoming limitations shown by past codes. In fact, the soil represents one of the most important substances found on the Earth, given that, covering its uppermost layer, it provides the support for all the terrestrial organisms and guarantees the terrestrial life (Lin 2011). In the last decade, many tentative were made for incorporating the soil behaviour into landscape evolution models, as well as in combining soilscape and landscape modelling (Ebel et al. 2007;Welivitiya et al. 2019). Indeed, on the one part, the soil controls the interaction between vegetation and water, as well as the atmosphere in terms of carbon and nutrient cycling. On the other part, in combination with the vegetation, the soil determines the rate of erosion and deposition and therefore cannot be ignored in adequately simulating the long-term dynamics of fluvial landscapes (Wilkes et al. 2019).
Aiming to draw a more complete picture of the evolution of a natural landscape, research should be directed also in including the role of biota, the dynamics of streamchannel adjustment, the erosion and transport of sediments and material by means of woody and debris flows or other mass movements and the formation and evolution of the critical zone . In summary, there is an evident need for a better understanding of mechanics and feedback of the physical, chemical and biological controls that can have a role in shaping the landscape forms. Even if the importance of the sediment dynamics is well recognized (Schumm and Lichty 1965;Robinson and Slingerland 1998), a major focus should be posed towards the study of the soil properties, improving the understanding on the influence of grain size, transport and sorting in shaping river systems, given that such aspects received some attention only recently (Gasparini et al. 2004;Di Silvio and Nones 2014;Sklar et al. 2017). Moreover, further research shall be devoted in better interpreting and modelling the links between local climate, relief and grain size delivery to sedimentary basins, aiming to obtain a most reliable estimate of the processes acting at the watershed scale. In fact, fluvial transport capacity and competence are highly sensitive to grain size composition, and, consequently, phenomena like abrasion, weathering and armouring can have a significant impact on the transport mechanisms (Gasparini et al. 2004;Attal and Lave 2006), pinpointing the opportunity to account for them in landscape models at the small scale (Willgoose and Sharmeen 2006;Román-Sánchez et al. 2019) or to justify their absence in the case of basin-wide simplified approaches Varrani et al. 2019).
The geometry of river represents another challenge for landscape evolution modelling given that size and shape of a channel controls (and is indirectly controlled by) the distribution of friction and energy dissipation across its wetted perimeter. There are an increasing number of numerical and experimental studies focused on evaluating how the channel geometry follows the changes of base-level controls, tectonic tilting and water and sediment supply (Stark 2006;Finnegan et al. 2007;Wobus et al. 2008;Davy et al. 2017;O'Hara et al. 2019), but these relationships are far to be fully understood. Even if landscape evolution models that couple surface changes with vegetation dynamics have begun to appear in the last decades and now are becoming a major research field (Murray and Paola 2003;Nones and Di Silvio 2016), the scientific community has just started to find a consensus on the quantitative relationships between hydrological, biological and geomorphic processes. In addition, challenges posed by modelling such mechanisms are related also to their spatiotemporal scales, which could be significantly different.
Aside from the theoretical and physical understanding of the mechanisms involved, as well as their mathematical schematization and the associated computing challenges, one of the greatest limitations in widely applying landscape evolution models is the evident lack of data and methods to test them. As stated by Mark (1975), the reliability of a landscape model can be assessed through four ways, depending on the process under evaluation. First, in the case of rapid landform development measurable in timescales of months or years, such as in the case of gully formation and postmining landscape (Hancock 2006;Hancock et al. 2017), the model predictions can be tested against direct observations. However, because of uncertainties in understanding delayed effects of processes like vegetation encroachment and weathering, problems can arise in extrapolating information from such newly created landscapes to the long term (Moliere et al. 2002). Second, there are situations where real-time measurements of sediment and solute fluxes can provide a useful basis for evaluating the model performance, even in the case of a slow rate of landform variation (Montgomery and Dietrich 1992). Third, the development of landscapes can be evaluated by means of scaled experiments, where the involved process can be adequately measured in a controlled environment. Since decades, laboratory experiments are very helpful in addressing specific issues focusing on a few geomorphic features (Hasbargen and Paola 2000;Pelletier 2008), but, generally, they perform well only under a qualitative point of view. The recent boom in high-speed computing resources and digital photogrammetry permits to overcome the operative limitations present in the past towards a more quantitative estimation of the landscape evolution at the laboratory scale, which could ultimately support numerical models. Indeed, because of many limitations correlated with laboratory tests and parameters uncertainty (Skinner et al. 2018), physical experiments and numerical model should be combined to adequately reproduce complex landscape changes. Fourth, landscape models can be tested by comparison with natural experiments, which are case studies having sufficiently constraints to allow for a quantitative comparison between field observations and numerical outputs (Montgomery and Dietrich 1992;Tucker 2009). While, since years, there is a great potential in combining natural experiments and modelling runs (Hancock and Willgoose 2001), the need to develop probabilistic frameworks and robust statistical measures for discriminating between dissimilar natural landscapes and scales remains (Schumer et al. 2017).
As visible from the present review, landscape evolution models can be extremely complicated, depending on the processes considered. In fact, the basin surface is not simply shaped by the interaction between water and sediment, but also many other factors can play a major role. As an example, in different parts of the world, the landscape is shaped by wind (Okin and Gillette 2001), snow (Liston and Elder 2006), ice (Ugelvig et al. 2016) or the fire (Scott 2018) since millennia. However, the description of such phenomena is outside the scope of the present work and is therefore not addressed here.
Human-induced alterations of river channels and basins as well as the climate change are actually having a major role in reshaping geomorphic systems, altering the natural relationship between the components and posing additional challenges to river modellers. As suggested by recent studies, topography and climate are generally coupled, and precipitation increases because of orographic effects during the uplift of a high mountain till an elevation of about 1000-2000 m (Bookhagen and Burbank 2006). Therefore, altered climatic conditions can drive to a change in precipitation patterns, with consequences on the watershed topography. However, the landscape sensitivity to the timescale of climatic variations is not yet completely understood and represents a hot research topic for the future (Moussirou and Bonnet 2018). In addition to the large-scale effects of climate change, local changes on the landscape can be also caused by the presence of animals like herbivorous. Nowadays, the correlation between animals and landscape evolution is generally studied only at the local level (Butler et al. 2007), and a general framework is still missing because of the high complexity of considering the presence of human beings and animals.
The present review showed that substantial progress has been made in quantitative modelling the evolution of the Earth's surface in terms of water-sediment-vegetation interactions, but much still remains to be accomplished and there are many open questions that can be addressed in the future. On the one part, there is the need for refining and testing landscape evolution models in a larger variety of cases to cover a multitude of spatial and temporal scales, by means of new and improved computing techniques. On the other part, one of the major challenges lies in developing experimental and field-based datasets for testing and validating numerical models across a wide range of spatiotemporal scales and covering different geomorphic environments (Rixhon et al. 2017).
Funding The work was also supported within statutory activities No. 3841/E-41/S/2019 of the Ministry of Science and Higher Education of Poland.

Compliance with ethical standards
Conflict of interest The author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.