the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
CHONK 1.0: landscape evolution framework: cellular automata meets graph theory
Boris Gailleton
Luca C. Malatesta
Guillaume Cordonnier
Jean Braun
Landscape evolution models (LEMs) are prime tools for simulating the evolution of sourcetosink systems through ranges of spatial and temporal scales. A plethora of various empirical laws have been successfully applied to describe the different parts of these systems: fluvial erosion, sediment transport and deposition, hillslope diffusion, or hydrology. Numerical frameworks exist to facilitate the combination of different subsets of laws, mostly by superposing grids of fluxes calculated independently. However, the exercise becomes increasingly challenging when the different laws are interconnected: for example when a lake breaks the upstream–downstream continuum in the amount of sediment and water it receives and transmits; or when erosional efficiency depends on the lithological composition of the sediment flux. In this contribution, we present a method mixing the advantages of cellular automata and graph theory to address such cases. We demonstrate how the former ensure interoperability of the different fluxes (e.g. water, fluvial sediments, hillslope sediments) independently of the process law implemented in the model, while the latter offers a wide range of tools to process numerical landscapes, including landscapes with closed basins. We provide three scenarios largely benefiting from our method: (i) one where lake systems are primary controls on landscape evolution, (ii) one where sediment provenance is closely monitored through the stratigraphy and (iii) one where heterogeneous provenance influences fluvial incision dynamically. We finally outline the way forward to make this method more generic and flexible.
 Article
(3487 KB)  Fulltext XML
 BibTeX
 EndNote
The timescale of sediment transport along the sourcetosink system from upstream erosion to downstream deposition is relatively short compared to the timescales of other geological processes. However, its large spatial extent on areas and the sometimes great intermittence of activity make it difficult to measure and observe directly (e.g. Sadler, 1981; Jerolmack and Sadler, 2007; Ganti et al., 2011; Schumer et al., 2017). Various models, analogue and numerical, help explore sourcetosink systems at different temporal and spatial scales that complement field observations. Analogue models offer a time compression in scaled experiments to rapidly simulate long time spans with complex physics but relatively simple environmental forcing (Babault et al., 2005; Paola et al., 2009; Guerit et al., 2014). Alternatively, numerical landscape evolution models (LEMs) have the advantage of giving complete control of the simulation. However, they rely mostly on empirical laws and are often limited to specific geoscience problems. For example, the evolution of surface topography over millions of years can be efficiently explored with erosion laws that only indirectly consider sediment transport in their numerical scheme (e.g. Yuan et al., 2019; Hergarten, 2020) or even completely ignore them (e.g. Braun and Willett, 2013) to the benefit of numerical performance. On the contrary, bedrock incision can be advantageously ignored when the focus is a highresolution modelling of sediment redistribution at very short timescales (e.g. Sklar and Dietrich, 1998; Croissant et al., 2017; Coulthard et al., 2013; Roelvink and Van Banning, 1995). A major challenge therefore lies in finding the right combination of laws to best address a given problem (Barnhart et al., 2019). That is complicated by the significant impacts of any change in numerical or physical parameters both in terms of quantitative results and computation time (e.g. Campforts et al., 2017; Armitage, 2019; Grieve et al., 2016).
While some process laws are implemented in standalone models (e.g. Hergarten, 2020; Braun and Sambridge, 1997; Coulthard et al., 2013), mature frameworks exist to facilitate the combination of different LEM components and to explore their results (e.g. Barnhart et al., 2020; Bovy, 2019; Schwanghart and Scherler, 2014; Mudd et al., 2019). However, these frameworks and models are designed to combine process laws at grid level where, for example local minima, flow routing or river incision are successively and serially solved. This can be a problem when studying more complex sourcetosink systems with multiple processes that are interconnected. Let us picture a situation where a lake acts as a local sink. Its sediment and water budget depends on all processes involving sediments or water upstream (Fig. 1). Its filling will, in turn, impact the behaviour of the process law downstream by modifying the amount of water and sediments they will transmit. This is incompatible with grids solved independently, or it requires exchange data structures that increase the runtime exponentially. Beside the question of local minima as in the lake example, the role of sediment fluxes is perhaps the most representative example of interconnectivity. Sediment starving or an abundance of clasts in a river will impede bedrock erosion by a lack of tools or an excess cover (Sklar and Dietrich, 2004; Finnegan et al., 2007; Geurts et al., 2018). The relative strengths of erosive clasts and erodible bedrock can also significantly enhance or diminish the erosion efficiency of a river and trigger local and nonlocal consequences (e.g. Gailleton et al., 2021; Sklar and Dietrich, 2001, respectively). Alluvial dynamics matter for sourcetosink studies: aggradation–incision cycles in mountain piedmonts can delay sediment delivery to the ocean by >10 kyr, a climatically relevant timescale, and recycle old signals in the sediment stream (e.g. Clift and Giosan, 2014; Malatesta et al., 2018; Dingle et al., 2020). Modulation of sediment fluxes also leads to prominent alluvial terraces and surfaces that are key for landscape interpretation (Bufe et al., 2017; Tofelde et al., 2017; Malatesta and Avouac, 2018). Increasingly fine resolution in stratigraphic studies warrants a renewed attention to the trajectory of sediment tracers across the landscape (Tofelde et al., 2021). Novel radiometric methods enable the exploitation of new sedimentary signatures that require a precise understanding of the rate and path of transport of sediment across landscapes (e.g. Lupker et al., 2017). Modelling sediment fluxes at the level of details that field and analytical studies can now attain benefits from the holistic approach presented here rather than the independent implementation of individual processes.
In this contribution, we propose a novel methodology, CHONK, to develop frameworks that include finegrained modularity in a cellbased referential, so as to ensure interconnectivity between LEM properties. CHONK is built to guarantee unconditional access to a common numerical toolkit regardless of the type of geomorphological laws employed. The cellbased referential allows one to track the parameters and/or to explore the dynamic feedback within the different fluxes transported from one cell to another. We demonstrate the potential of integrating cellular automata elements with graphbased finite difference methods to resolve sedimentary dynamics necessary for sedimentological studies of landscape evolution. This contribution presents the core architecture of CHONK, while several collaborative projects employ and apply the framework to address sedimentological and geomorphological challenges. These projects concurrently inform the development of a userfriendly platform to be progressively released in the coming months and years.
First, we concisely present and explain the new method. We then detail the model structure, its different algorithms, and the process laws we picked for the demonstration. Finally, we present and discuss different scenarios demonstrating the capabilities of this new method.
The new formulation we introduce in this contribution combines the advantages offered by the cellular automata methods (von Neumann, 1951; Wolfram, 1984) and graphbased finite difference methods commonly used in LEMs and frameworks (e.g. Bovy, 2019; Barnhart et al., 2020; GarciaCastellanos and JiménezMunt, 2015; Braun and Willett, 2013). We first briefly define and review the existing methods and framework to explain our rationale for creating a new one.
2.1 Graphbased frameworks and methods
LEM frameworks typically solve the different components of landscape evolution modelling independently following a graphbased logic applied on data grids. In other words, fluxes and other quantities (e.g. elevation, erosion, water) are discretized on 2D arrays that are calculated and combined successively. Geomorphological processes typically require downstream transfers of fluxes (e.g. drainage area, water or sediments); upstream propagation of numerical schemes (e.g. Braun and Willett, 2013; Campforts et al., 2017); or even successive iterations of both (e.g. Yuan et al., 2019; Hergarten, 2020). LEMs compute those by building a directed acyclic graph (DAG) where each discretized location defines a node (or vertex) with directed link (or edge) towards its downstream neighbour(s). This data structure enables operations such as graph traversals or topological sorting which ensure the required downstream or upstream processing of nodes. LEMs can integrate the graph structure explicitly (i.e. computing the vertex and edge structures), taking advantage of graph theory algorithms such as topological sorting (e.g. Braun and Willett, 2013; Anand et al., 2020) or more sophisticated localminima processing methods (e.g. Cordonnier et al., 2019; Barnes et al., 2020); they can also use intermediate data structures such as priority queues to navigate in complex, depressionbearing landscapes without having to store edges (e.g. Barnes et al., 2014b); or simply sort nodes by decreasing elevation after eventually processing local minima (e.g. Braun and Sambridge, 1997; Carretier et al., 2016; Hergarten, 2020).
A typical graphbase LEM flow can be illustrated with the streampower incision model (Howard and Kerby, 1983), a widely used fluvial incision equation where erosion rate is defined as a function of slope and drainage area. LEMs first compute a graph structure to calculate drainage area and weigh it by water influx to have a proxy for water discharge. At this stage, local minima (i.e. lakes, endorheic basins or noise) are filled with water or carved out at this stage in order to ensure flow continuity (e.g. Braun and Willett, 2013; Cordonnier et al., 2019; Barnes et al., 2014b; Salles, 2019). Then the models compute topographic slopes from the previously filled surfaces and finally combine both data grids with an erodibility parameter to calculate fluvial incision rates in each cell. Such gridbased formulation is very flexible and allows for modifications in order to adapt the model to a geoscientific problem. For example adding hillslope diffusion to the workflow (e.g. Roering et al., 1999) would consist in calculating it after fluvial incision and combining the grids of elevation changes at the end of the process. Implementing alternative methods is also straightforward: for example switching from the stream power incision model to a transportlimited one (e.g. Sklar and Dietrich, 1998) only requires replacing the section of the sequential process that deals with the fluvial process.
However, the combination of independently calculated grids reaches its limitations when processes are interdependent (e.g. calculating fluvial incision function of the nature of its upstream sediment input combining all the processes and what could have been stored in potential lakes). Let us consider an example where lakes are of great importance for landscape evolution. They act as intermediate traps in the domain and the amount of sediments and water they may transmit downstream as a function of the overall amount they receive. The influx can only be known if all the processes happening upstream of the lake have been processed. This is not compatible with a sequential treatment of processes and requires specific implementations (e.g. GarciaCastellanos and JiménezMunt, 2015; Salles, 2019) where any modification in the methods or processes requires significant work to redesign its whole implementation. This is not always straightforward, is often highly dependent on the actual numerical format of the LEM and is accompanied by an unavoidable loss of flexibility and modularity.
2.2 Cellular automata method
Cellular automata models are reducedcomplexity models designed to tackle discretized problems on networks of connected cells (von Neumann, 1951; Wolfram, 1984). The cells have given properties and states which evolve as a function of the states of their neighbours according to a set of rules. Road traffic modelling by Nagel and Schreckenberg (1992) is a good illustration of cellular automata logic. Cells represent a stretch of road and their properties include, for example the presence or absence of cars, their velocity, or whether they are Honda Jazz. Cells evolve as a function of the presence of cars in their linked counterparts following simple rules to simulate road traffic. Cellular automata methods can also include more sophisticated equations and processes and have been utilized for modelling elements of landscape evolution such as water and sediment fluxes (Coulthard et al., 2013); tracking particles in flow (e.g. Tucker et al., 2016); hillslope evolution (e.g. Tucker and Bradley, 2010; Jyotsna and Haff, 1997); soil erosion (e.g. D'Ambrosio et al., 2001); sediment transport and channel morphology (e.g. Salles et al., 2007). Frameworks exist to take advantage of cellular automata (e.g. Barnhart et al., 2020, partially implemented in Tucker et al., 2016). It is important to note the cells are processed in no particular order and cannot propagate nonlocal fluxes such as drainage area within a single time step.
2.3 Hybrid solution: cellular automata on a graph
We developed a new formulation combining the advantages of the graphbased and cellular automata methods. The aim is to make generic interactivity between fluxes and processes an intrinsic feature of LEM design. Building on the lake example we used in Sect. 2.1, if a simulation requires the sediment and water budget of a lake, one should be able to edit the processes affecting water and sediment fluxes without having to modify the entire workflow. This requires a number of numerical constraints: (i) fluxes should be defined separately from processes in order to let a theoretically unconditional number of processes affect the fluxes; (ii) the processing of the graph should be as independent as possible of the processes – meaning that the resolving of local minima should not imply they are systematically filled; and (iii) every interconnected processes should be processed simultaneously within a single location before transmitting fluxes to the next. In this contribution, we present CHONK 1.0, a proofofconcept for this modelling design with its functional workflow.
To do so, we process a topographic grid from which we calculate a depressionaware graph of downstream/upstream connectivity. The latter does not assume the depression systems will be systematically filled, but instead preprocesses a data structure allowing for different scenarios (i.e. different levels of details in the processing). Every node (i.e. every discretized location as described in Sect. 2.1) on the grid is then treated as a cell only processed in the downstream direction. The properties of the cells are the quantities and fluxes needed by the equations implemented in the model – e.g. elevation, bedrock incision, water, sediment fluxes (Fig. 2). Like cellular automata methods, all the processes affecting fluxes and quantities are processed before transfer to downstream cells. The combination of both methods ensures that when a cell belongs to a topographic depression, the definitive fluxes are known and this information can be used to fill that depression. The cells are processed in a specific order following a topology dictated by the topographic graph. Lakes are just one example benefiting from this method, and we will demonstrate that this method is particularly adapted to solving problems and challenges involving interdependent feedbacks that are difficult to tackle otherwise.
The prototype we developed for this contribution is a first step towards a fully fledged, generic and dedicated framework like Barnhart et al. (2020) or Bovy (2019). We implemented a specific set of process laws in order to test and demonstrate the advantages of this method.
2.4 Comparison with existing models or frameworks
A number of numerical tools already exists for landscape evolution modelling. It is not in the scope of this paper to provide an exhaustive review of all of them, but it is important to mention those most relevant to our goals.
Cellular automata models have been utilized for LEMs at basin scale. Coulthard et al. (2013) developed CAESARLISFLOOD, a cellular automaton model approximating the shallowwater equations (Bates et al., 2010) and designed to explore fluvial sediment transport and bedrock erosion over the timescale of few thousand years. Like other LEMs solving similar equations (e.g. Davy et al., 2017; Adams et al., 2017), this family of methods is not designed for geological and/or mountainrange scale because they are (i) numerically limited by the short time steps required to keep the finite difference scheme stable and (ii) philosophically limited by the amount of external constraints required (highresolution precipitation patterns, for example). CAESARLISFLOOD also processes all the cells in any arbitrary order (or even in parallel) and only transfer sediments and water from a cell to its immediate neighbour within a single numerical time step – as opposed to a full landscape traversal per time steps for longerterm LEMs.
The landscape evolution modelling community benefits from wellestablished frameworks to develop and design LEMs and topographic analysis tools (Barnhart et al., 2020; Schwanghart and Kuhn, 2010; Mudd et al., 2019; Bovy, 2019). They all rely on routines manipulating a topographic grid and building a graph of node connectivity for it. However, existing frameworks are primarily designed to be solved on grids alone and they inherit their limitations (Sect. 2.1).
Accounting for lakes is one of the main features of our method and we note that different approaches already exist. A common family of methods consists in preprocessing local minima by directly altering the topography in order to force an outflow by either carving a way to an output (Lindsay, 2016) or filling them (Wang and Liu, 2006; Barnes et al., 2014a). Both force the local minima to connect to the rest of the landscapes and the flow to escape via the edges. Bovy (2019) utilizes an alternative method by Cordonnier et al. (2019) leveraging graph theory to simulate carving/filling without affecting topography. It is worth noting that some algorithms have been specifically develop to process, calculate and fill depressions with an arbitrarily given amount of water. Among these, the closest to our aim are the developments by Callaghan and Wickert (2019); Barnes et al. (2020, 2021). They designed a set of methods to (i) identify, (ii) hierarchize and (iii) fill the depression with a particular focus on numerical efficiency. It is worth noting that the model developed by Callaghan and Wickert (2019) is a cellular automaton. While we partially built our numerical method to manage lakes on these previous developments, there are a couple of differences, most of them related to our need to integrate the lake solver into a preexisting multipleflow graph for node connectivity. More detailed differences are outlined in Sect. 3.3.2. Geurts et al. (2018), for example utilized the method of Braun and Sambridge (1997) to simulate lake filling by stopping flow at the lake bottom and only connecting the lake to the rest of the landscapes once filled with fluvial sediments. Campforts et al. (2020) or Yuan et al. (2019) enhanced fluvial deposition in lake areas in order to roughly approximate lake deposition. These methods acknowledge the importance of lakes in the landscapes, but do not treat them as separate domains with dedicated processes. Salles (2019) characterize lakes by first filling the topography with the approach of Barnes et al. (2014a) and identifying areas of topographic change. Their model then traps all the sediment carried in these domains, transmitting potential excess to the downstream landscape. This method is close to what we aim to achieve but considers lake as unconditionally filled and outletting, thus not designed for endorheic basins. TISC (initially described by GarciaCastellanos and JiménezMunt, 2015) is a pioneer in terms of integrating endorheism to LEMs and recognizing its impact on landscape evolution (GarciaCastellanos et al., 2003; GarciaCastellanos, 2006; Struth et al., 2021). TISC calculates the topography of the depressions and fills them gradually with the available sediment and water. Excess material is only transmitted to the outlet and downstream landscapes if available, successfully simulating closed lakes and endorheism when lake evaporation or infiltration balances precipitation. The implementation of TISC, however, is not compatible with our design as water fluxes are calculated separately from the rest of the processes. Runoff is first calculated and the lakes are gradually filled, dynamically accounting for evaporation and lake spilling. Other processes are only calculated after the water flux is defined whereas CHONK calculates all the fluxes simultaneously for each separate cell, allowing interconnectivity between their properties.
Finally, tracking sediment provenance in LEMs has been done in different ways. Carretier et al. (2016) add discrete Lagrangian particles on top of Eulerian grids. They postprocess erosion, entrainment and deposition of sediment fluxes to determine the movement of these particles with a probabilistic approach. Sharman et al. (2019) integrate the erosion field to backcalculate provenance from labelled areas. These existing methods have in common that they are postprocessing the tracking, i.e. calculating the proportion of sediment provenance after the calculation of the surface process laws. We aim in this contribution to embed the trackers into the model in order to make possible their integration directly in the process laws, for example adjusting fluvial erosivity to the proportion of a certain rock type in the sediments relative to the local bedrock type.
3.1 Generic numerical structure
Before describing the technical details inherent to CHONK 1.0, we provide a generic description of the modelling design – outlining the steps required for a general implementation following the same principles.
As illustrated in Fig. 2, the first step is to build the cellular structure by determining the fluxes and properties needed. These can be spatial data (e.g. precipitations, elevation, sediment provenance in stratigraphy), fluxes (e.g. fluvial sediments, water) or process parameters (e.g. erodibility). The second step defines the processes, i.e. the laws defining the interactions between all the fluxes and processes (e.g. fluvial incision, hillslope diffusion). Processes and/or fluxes can be domainspecific (e.g. marine, fluvial, lake, glacier). Finally, a graph structure providing the order of processing for the nodes needs to be determined. The graph structure has to be processagnostic and capable of acknowledging domains of different topology. Algorithm 1 presents a simplified simulation.
An ideal numerical implementation of this principle should numerically separate fluxes, properties, processes and graphs. While complicated, numerical designs such as loose coupling can achieve this and ensure the different elements of the framework do not require presence or awareness of the others. In other words, it enables the addition, replacement or removal of processes affecting the same fluxes without needing to modify the rest of the model. For example, the model could change from a regular grid to a 1D profile or a Voronoi grid by simply “switching” the graph module. Some existing frameworks (e.g. Barnhart et al., 2020; Bovy, 2019) follow similar numerical design, but not in a cellular referential as we advocate in this contribution.
3.2 Building a directed acyclic graph
The first step consists in building a graph of connectivity on the landscape in order to determine a processing order for the cells that takes into account the topography of endorheic basins. Here, we use a regular rectangular grid to discretize topographic elevation, z. Each individual location, i, is a node from the point of view of the graph and holds a cellular automata cell. Each cell is connected to adjacent neighbours with the D8 direction, i.e. encompassing all the cells in diagonals and side directions. This defines the node graph, where for any given cell i we call all the connected cells with lower z receivers r_{i} and all the connected cells with higher z donors d_{i}. Cells with no donors are referred to as “source cells” while cells with no receivers are “pits” (if internal) or “edges” (if located on a matrix boundary). We implemented two types of boundary conditions at the edges: (i) open boundaries, where fluxes can escape the model and cells have no receivers, and (ii) periodic boundaries, where the fluxes communicate with the opposite cells (e.g. a cell at the eastern boundary is linked to its opposite at the western boundary). The graph hence created is a directed acyclic graph (DAG): each cell is linked to one or several receivers and cannot cycle back to itself. In graph theory, setting up a DAG allows for the use of a wide range of dedicated algorithms for topological ordering or graph traversals. The type of flow emulated by this DAG is called a “multiple flow direction” (Schwanghart and Scherler, 2014) as one cell can be linked to multiple receivers.
Note that the case of numerically flat surfaces, i.e. a node surrounded by others with the exact same elevation at numerical precision, needs particular care. In such situations, neighbours of i can end up being neither a receiver nor a donor and can generate cycles. Methods exist to process the flat surfaces (e.g. Barnes et al., 2014a). We use the carving algorithm by Cordonnier et al. (2019) to approximate an acyclic flow direction on these flat surfaces; the algorithm is detailed in the next subsection.
3.3 Computing a depressionaware topological order
Once the connection between cells is established – i.e. the receivers and donors of each cell are determined by the topography or by the rerouting algorithm on flat surfaces – we compute the topological order. It is a crucial step for any LEM: it determines the order in which cells need to be processed starting from the source nodes and finishing with the model edges. Alternative methods exist: it is possible, for example to utilize an iterative method accumulating fluxes progressively (e.g. Braun and Sambridge, 1997); solving large sparse matrices (e.g. Perron, 2011); using priority queue data structures to traverse the graph of cells dynamically (e.g. Barnes et al., 2014b, 2020).
Our implementation of an algorithm calculating a lakeaware topological order needs to satisfy a number of conditions: (i) conservation of the original topography of the depression in order to take its characteristics into account, and (ii) respect of the notion of upstream and downstream including potential lake and depression systems. We implemented two different algorithms to incorporate local minima in the model. First, a topological order can “passively” reroute local minima and approximate the flow path as if depressions were recognized but assumed to be entirely filled up to the elevation of the outlet. Second, an algorithm accounts for the volume of potential lakes, and uses separate dedicated processes within them. Both of the algorithms modify the DAG in order to emulate a notion of upstream/downstream by linking the pit nodes of the different depressions to an adequate outlet node. Finally, we apply a topological sorting algorithm on the modified DAG to calculate the depressionaware topological order.
We use the same topological sorting algorithm to calculate the topological order for both scenarios (detailed in Sect. 3.3.1 and 3.3.2). The algorithm is a modified implementation of the one in fastscape (Bovy, 2019) and very similar to the one described by Anand et al. (2020). It is O(n) in complexity, with n being the total number of links between the cells and their receivers in the graph. In short, a queue is initialized with the source cells. In turn, these are popped out of the queue, pushed into a stack array and their receivers are visited. An array tracks the number of times each cell is visited. If the number of visits equals the number of donors of a given node, it is saved into the stack and the process continues. Once the queue is emptied, all the cells have come through and the stack array contains all the node indices ordered. This stack array can be traversed in normal or reverse order to respectively process upstream or downstream cells first and is illustrated in Fig. 3. This process is equivalent to the steepest descent alternative of Braun and Willett (2013).
3.3.1 Topological order for landscapes with passive lakes
The solver for passive lakes is designed for cases where depressions are a secondary feature of the landscape evolution study. It ensures flow continuity through the landscape and conservation of original topography by connecting the pit of each depression to an outlet that will eventually reach the model edge (Cordonnier et al., 2019). The solver bypasses the computational expense of considering the exact geometry of the depressions while still accounting for their existence. Our method is adapted from the work of Cordonnier et al. (2019), where the steepest descent graphs reroute local minima towards model edges. We only modify it to be compatible with multiple flow directions. The algorithm first links every node to either a single edge or a pit using a steepest descent route to define basins. It then links pairs of adjacent drainage basins using their lowest connections from the most internal to the most external one. This defines a receiver cell for the internal pit of each internal basin in order to ultimately drain to the edge.
Our implementation adds a couple of extra steps. First, our algorithm actually carves the surface in the case of flat surfaces in order to avoid 0slopes (see Algorithm 3 in Cordonnier et al. (2019) which inverts the nodetonode steepest descent connections from the sill to the pit). We make sure to reassess all the potential multipleflow links impacted by this singleflow rerouting (e.g. cells from the target basin partly flowing to the source basin). After this step, the topological order can be computed and will route flow through depression.
This method has the advantage of speed, versatility and stability as demonstrated by the benchmarks of Cordonnier et al. (2019). However, the links between basins are estimated with a steepest descent algorithm, which might shift the location of the geometrical outlet of the depression by a few pixels. It also maintains unconditional connectivity between local minima and their outlets, ignoring endorheism.
3.3.2 Topological order for depressionaware simulations
The depressionaware solver fully embraces the topographic complexity of depression systems. It does not assume the lakes outflow and treats them as separate domains. The geometry of depression systems can be convoluted with multiple levels of subdepressions (Fig. 4). To deal with this complexity, we build a binary trees for each depression system with a principle adapted from Barnes et al. (2020): each vertex of the binary tree can only have up to two children, one parent and one twin. It is built with a “vertical” logic illustrated in Fig. 4 where each vertex corresponds to a spatially identifiable domain made of a single or multiple merged depressions. Building such trees ensures efficient operations to numerically navigate through individual depression systems (Barnes et al., 2020). We refer to the work of Barnes et al. for a detailed description of this binary tree and how to efficiently build it. This data structure has been utilized for flooding landscapes with finite amount of water (Barnes et al., 2021).
Our implementation differs from Barnes et al. (2020) for a couple of important points. For this contribution, the binary tree needs to be fully integrated within the topographic graph in order to enable topological sorting and downstream traversals. We therefore sacrificed some of the computational efficiency (in terms of memory and CPU) to store more information for communication between the topographic DAG and the binary trees of depression. Barnes et al. (2020) build a forest of binary trees connected to each other and to the “ocean” that Barnes et al. (2021) use to iteratively flood the landscape from one depression to another. Water flows through the whole landscape to the depression bottom and is then redistributed from one depression to another until all the water is used or all depressions are filled. Instead, we build independent local trees that are only connected to their surrounding DAG. Numerically, we only label nodes belonging to a depression inside the corresponding depression system instead of labelling and linking all nodes across the landscape to a depression like Barnes et al. (2020). We also store a lot of information on a cell basis, for example which cell belongs to which depression sorted by elevation. We opted for this heavier solution because, contrary to Barnes et al. (2021), we do not fill the depressions iteratively and only visit each cell once, as explained in Sect. 2.3. The most significant difference is perhaps the flow topology: we built CHONK to be compatible with multiple flow directions while the Barnes et al. (2020) and Barnes et al. (2021) counterparts are singleflow oriented. Thus a depression system can be linked to multiple others, whereby a steepest descent route can only link one depression system to another at a time. This point makes our algorithm significantly more convoluted, especially in the presence of complex systems of nested depressions (e.g. white noise). We note that all these modifications added to fit our needs complicate the original algorithm and can make it significantly slower in many cases. We are not presenting a version with better computing speed or accuracy compared to the work of Barnes et al. (2020) and Barnes et al. (2021). We adapt its use to our prototype to explore the consequences of explicitly considering lakes in LEMs. A cleaner, performanceoriented solution could benefit from being entirely based on Barnes et al. (2020) and Barnes et al. (2021); however, using their version out of the box would require significant work to achieve all the features we need for CHONK.
We make heavy use of priority queuebased algorithms to build the graph (see Barnes et al., 2014a and Barnes et al., 2020 for full details about this data structure). This enables the dynamic sorting of selected cells as a function of their elevation. First, we place each internal pit cell in individual priority queues as a starting point for all the base depressions and we label the cell with a unique depression ID (black dots in Fig. 4a). We process each priority queue until it is empty by popping out the lowest elevation cell and checking all of its neighbours. If the neighbour is higher in elevation, it is placed in the queue for later processing or labelling. This process runs until the cell being processed is already labelled as belonging to another depression. In this case both are registered as twins. Each twin records the connecting cell as their tipping node (e.g. depression 2 and 3 in Fig. 4a). If one of the neighbours has a lower elevation, that cell is labelled as outlet and this depression is placed at the top of its tree – or remains an outlet as long as it is not labelled as a twin by another priority queue. The trees are complete once all priority queues are empty. Note that while we do not detail each and every one of them for the sake of conciseness, the algorithm needs to potentially manage a lot of specific edge cases.
The data structure allows us to process the following metrics for each depression:

a depression level, which represents the maximum distance in the tree from a base depression. Each base depression is at level 0, and each parent's level is equal to the maximum level of their children plus 1,

the minimum volume of a depression (0 if base depression, the minimum volume to fill all the children and “reach” the depression in the tree),

the volume of the depression V_{total} if filled, note that it includes the volume of their children if any,

the maximum elevation of the depression if filled,

the tipping node of the depression, which represents either the outlet of the whole subsystem, or the node joining two twins.
In addition to the depressionspecific information, the model stores a number of internal structures to navigate between the topographic graph and the depression tree. Note that the maximum volume of water can account for potential evaporation if it is enabled in the model.
Our depression tree relies on the principle of uniqueness of the tipping points which can be invalidated by numerically flat surfaces or if depression borders have equal elevation. To prevent this, we add minute numerical noise between ${\mathrm{10}}^{\mathrm{6}}$ and 10^{−6} m at each time step and carve depressions with insignificant volumes using Algorithm 3 of Cordonnier et al. (2019).
After building the depression tree, we can finally calculate the topological order for the depressionaware lake solver. This is achieved in the DAG by temporarily linking the pit cell of each base depression to the cells that lie downstream of the outlet of the above depression in each system (Fig. 4a). These links ensure that any lake system will be processed before its downstream counterparts and is cancelled after the calculation of a topological order.
3.4 Cellular automata structure
3.4.1 Properties, parameterization and tracking
Once the DAG is built, the model skeleton is ready and a cell is attributed to each node. The information held by each cell can be adjusted and expanded on a needs basis. In the current implementation, cells have the following properties updated at each time step (illustrated in Fig. 2):

topographic elevation (in m)

thickness of the immobile sediment layer (in m)

volumetric water flux Q_{w} in m^{3} yr^{−1} traversing the cell

volumetric sediment flux Q_{s} in m^{3} yr^{−1} traversing the cell (mobile sediments)

proportion of sediment flux from the river or the hillslope systems

a list of the downstream cells receiving either sediment or water in transit, calculated from the graph and from the process law implemented in the model

lists of weights describing the proportions of sediment and water transmitted to each downstream receiving cell

erosion, sediment entrainment and deposition fluxes

tracking information if activated (e.g. proportion of the sediment flux coming from a given source area).
Three kinds of parameter inputs are currently available. First, external parameters which can be single values (e.g. dx, dy, dt), or global arrays (e.g. 2D matrices of precipitation or uplift), varying in space and/or time. Second, parameters that are labeldependent: a 2D matrix of labels defines discrete spatial areas and each label has a set of distinct parameters, for example different rock type can be associated with different erodibility and diffusivity (Gailleton, 2021). And third, parameters that are fully dynamic: they are interdependent of each other and defined by a function rather than a given value. Examples of the latter are detailed in Sect. 4.4.
The tracking capabilities of the method also rely on the labels. While the numerical implementation is tedious, its principle is simple and powerful: any material eroded by any process from any location keeps track of its label when it is incorporated in the mobile sediment flux. The mobile sediment flux is partitioned to the receivers alongside the water flux, adjusted for local erosion and deposition processes. In the stratigraphy, a dynamic sparse matrix of cells stacks “containers” of sediments and keeps track of label proportions to guarantee tracking if reeroded.
It is worth noting, however, that this cellular automata structure has some numerical limitations. To maintain all the advances detailed in this contribution, all the calculations need to be processed from ridge to outlet, which is not necessarily compatible with all numerical laws. For example, solving stream powerlike equations implicitly necessitates multiple graph traversals in the upstream and downstream directions, therefore limiting the amount of upstream information a cell can use in the processes (e.g. Braun and Willett, 2013; Campforts et al., 2017; Hergarten, 2020). However, they are not fully incompatible: one could imagine calculating a “static” erosion field with one of these implicit schemes and postprocess them using the cellular automata method to integrate upstream information (e.g. provenance), only sacrificing the dynamic adjustment capabilities.
3.4.2 Cellprocessing order for local minima
The model processes the cells following the upstreamtodownstream topological order first assuming that there are no lakes. Before their turn, unprocessed cells receive water and sediments from upstream neighbours. When a cell is next in the topological order, the model applies external flux modifiers on it: precipitation, infiltration or any related process law affecting the water or sediment flux by addition or removal from external sources. Then the process laws affecting the cell are executed in the following order: water routing; fluvial incision, deposition and sediment entrainment; and/or hillslope diffusion following equations described in Sect. 4.1. At that stage, the model calculates weights for the distribution of sediment and water across the cell receivers. Finally, transport processes transmit water and sediment to the unprocessed receiver cells, along with the proportion of sediment fluxes respectively belonging to hillslope and fluvial domains.
When the solver for passive lakes is activated, cells in the area affected by the local minima are processed like any other cells. However, flow is rerouted from the pit cell to the lake outlet and this effectively reduces the topographic gradient, enhances deposition processes, and reduces erosion processes.
Solving lakes is done in multiple steps. Cells are processed normally, i.e. with fluvial and hillslope processes, following the downstream topological order. By definition, all the pit cells of a given depression system are processed before any section of the landscapes downstream of the lake. If the processed pit cell is the last of its depression system (in the case of a simple lake there is only one pit, but nested depression systems can have multiple pits) we can use the full volume of sediment and water in these cells to fill the lake(s) using the precomputed depression trees.
The first step consists in calculating the total amount of water entering the full depression system by summing Q_{w} for each pit node of base depressions in the system. The tree is traversed from bottom to top, propagating the water from children to parents. In the end, the following volume of water V_{w in} enters each depression:
where i_{pit} is the cell index of every pit cell downstream of a given depression.
The second step determines whether the depression system needs breaking into subtrees: the full tree is assessed from the top depression down. If the sum of available water is more than what the top depression can store, the whole lake system fills with water and will outflow. Otherwise if the minimum amount of water storable in the top depression is less than V_{w in}, the lake does not outflow but all the child depressions will be filled. Finally, if the minimum amount of water storable in the top depression is greater than V_{w in}, the local tree is divided in two and the assessment is reiterated until all the subtrees are filled. Note that all the water entering can also evaporate, in which case no lake is created.
The third step consists in calculating the elevation of the lake (z_{w}). Because our current implementation solves explicit finite difference schemes, we assume that, within a time step, the volume of water in the lake solely determines z_{w}. Elevation changes due to lake sediment deposition are only applied at the end of the time step. If the lake outflows, h_{w} equals the elevation of the outlet cell. Underfilled depressions lead to more complications (GarciaCastellanos and JiménezMunt, 2015). In these cases, the model calculates a balance between lake evaporation and the available amount of water. Using a priorityqueuebased graph traversal (see Sect. 3.3.2), we traverse the depression cells in increasing elevation order. Cells are included one by one and contribute in turn to storing the available amount of water V_{w avail} while giving their elevation to h_{w}:
where N_{lake} tracks the number of cells already in the lake, Q_{evap} is water lost to evaporation, and Δz the elevation difference between the current h_{w} and the elevation of the next node in the priority queue. The final h_{w} is calculated once ${V}_{\text{w avail.}}<{N}_{\mathrm{lake}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}x\phantom{\rule{0.125em}{0ex}}\mathrm{d}y\phantom{\rule{0.125em}{0ex}}(\mathrm{\Delta}z+{Q}_{\mathrm{evap}})$.
3.4.3 Water and sediment fluxes into and across lakes
Once the water elevation is determined, the model backcalculates sediments. All the cells below water are “deprocessed” from continental processes: fluvial and hillslope processes are reversed with adequate correction on cells sediments and water contents. The volume of sediment stored in the lake, V_{w evap}, can now be stored in the lake in a straightforward manner as its final volume is known. Any excess is transmitted to the outlet cell. As noted by GarciaCastellanos (2006) and GarciaCastellanos and JiménezMunt (2015), the outlet of the lake needs particular care as its behaviour through time ultimately controls its draining. The deprocessing of the outlet is only partial as it gives sediments to the lake and to the downstream landscape. Only the part of the fluxes going into the lake needs to be cancelled and the other parts need to be recalculated with the new amount of sediment and water. The latter are determined by subtracting the incoming V_{w} and V_{s} by what has been stored. Additional care is needed to consider water and sediment coming to the outlet from its nonlacustrine upstream neighbours.
We cannot stress enough how convoluted this deprocessing can be, numerically speaking. One needs to account carefully for all the neighbouring cells of the outlet and not remove/readd fluxes multiple times. Given the critical nature of this task, we make sure our model is not plagued by uncovered edge cases and we implemented massbalance checkers making sure no water or sediments are lost due to the transfer processes. Massbalance for a transferable flux can simply be defined as follows:
where Q_{in} encompasses any fluxes being added to the system and Q_{out} any fluxes leaving the system. For water flux, Q_{w in} includes effective precipitation as well any water stored in a lake at the previous time step (when using the depressionaware lake solver). Q_{w out} encompasses any water stored in a lake at the current time step, evaporation and water leaving the system via the model edges. For sediment fluxes, Q_{w evap} includes any process eroding material and putting it in transported flux (e.g. incision, entrainment, diffusion) and Q_{s out} any processes depositing sediment from this flux (e.g. fluvial deposition, lake deposition) or exiting the model via the edges. The mass balance is respected if M=0 (plus or minus numerical precision errors). Note that a simpler alternative would be to process the water fluxes separately to avoid the deprocessing. However, this contribution aims to develop a method able to keep unconditional interoperability between all the fluxes and processes, which would be broken by such sequential separation.
Finally, once all cells have been processed, each cell updates the topography and the sediment layer with its erosion and deposition fields. The model also calculates and formats data to monitor the direct model outputs (e.g. maps of erosion, water fluxes, sediment thickness) and indirect outputs such as the sum of the sediment fluxes outletting the model versus the sum of sediment fluxes stored in sediment layers.
We demonstrate the capabilities of the method with three fields of applications. First we test the effect of considering lakes in a tectonically active range with an internal basin. We then illustrate the tracking capabilities of the model by monitoring the sediments flux coming from a magmatic pluton. Finally, we explore the dynamic parametrization feature with the previous pluton settings, and adapt parameters in function of sediment flux composition. All the models start from the same nearsteadystate landscape obtained after running a simulation until drainage stabilization with block uplift and nonsubsiding foreland. The model has been tested on a computer with an Intel i910980HK and 32 GB or RAM on both MacOS 11.7, Windows 10 and Linux Ubuntu 22.04.
4.1 Process laws
To test the framework, we implemented a set of process laws that simulate longterm hydrology, fluvial, and hillslope processes.
4.1.1 Hydrology
Hydrology in longterm LEMs is usually approximated by a flowrouting algorithm distributing weighted drainage area from source nodes to outlets in the downstream direction. The weights represent the spatial variation of precipitation rates (see Leonard and Whipple, 2021, for a comprehensive review on the subject). First, the local effective precipitation is added to the water discharge in the considered cell i:
where P_{i} is the local effective precipitation weight factor that can include infiltration.
The second step is the routing to receivers. It can follow the steepest descent singleflow direction (e.g. O'Callaghan and Mark, 1984; Braun and Willett, 2013) or a multipleflow direction (e.g. Tarboton, 1997; Schwanghart and Heckmann, 2012; Armitage, 2019). We implemented an adaptive algorithm routing water with multiple flow following the method of Bovy (2019). We added an optional parameter to enable dynamic switching to singleflow routing after an arbitrary threshold of discharge, in order to roughly simulate a transition from hillslope to fluvial domains. Note that we use relatively large cell sizes (dx=200 m) and this parameter is optional. In the multipleflow domain, water is split according to the local slope. Following Bovy (2019), an exponent p_{r} is calculated for each receiver r of a cell i:
and then normalized to satisfy $\sum {p}_{\mathrm{r}}=\mathrm{1}$ and conserve mass balance. The water flux is then transmitted to each receiver with
4.1.2 Fluvial erosion and deposition
We simulate fluvial erosion and deposition using the SPACE model (Shobe et al., 2017), a hybrid law allowing for simultaneous treatment of the detachmentlimited and transportlimited portion of the rivers based on Davy and Lague (2009). SPACE can process all kinds of landscapes, whether sediments are absent or they saturate the system. The process law separates sediment entrainment ^{f}E_{s}, bedrock incision ^{f}E_{r} and sediment deposition ^{f}D_{s} into three equations solved simultaneously:
where K_{s} is the sediment entrainment coefficient regulating the ease with which sediment cover can be mobilized; K_{r} is the erodibility coefficient ultimately controlling local rock strength (proxying various factors such as weathering or fracturing); m and n are exponents regulating the relative importance of topographic gradient and water flux in the stream power law (e.g. Harel et al., 2016); H is the sediment height; H^{*} is the bed roughness index linked to the proportion of bedrock not covered by sediment; V is a dimensionless settling velocity coefficient encompassing information about the turbulence and composition of the suspended load. Details on all these parameters can be found in the original paper by Shobe et al. (2017). In the case of multiple flow departing from a single cell, the process is simply summed for each receiver: S, Q_{w} and Q_{s} being different for each one.
4.1.3 Hillslope diffusion
Following the same philosophy, we implemented the nonlinear hillslope diffusion of Carretier et al. (2016). This law separates sediment entrainment from deposition (i) enabling greater numerical stability than the purely nonlinear explicit scheme (Roering et al., 1999) while (ii) keeping the nonlocal, nonlinear aspect of the diffusion process. This law is versatile and collapses to both linear and nonlinear end members under different contexts as demonstrated in Carretier et al. (2016). Material entrainment follows a local, straightforward linear diffusion scheme which is defined as
where E_{rock} and E_{soil} are the entrainment rate in [L T^{−1}] for bedrock and sediment respectively; and κ_{rock} and κ_{soil} are modulating parameters as a function of the physical characteristic of the substrate and soil. Note that it is possible to disable bedrock diffusion to consider soil movements. In the case of multiple flow, we respect the numerical implementation of Carretier et al. (2016) considering that the steepest slope is the main driver to calculate $\frac{\mathrm{d}z}{\mathrm{d}x}$. If both E_{rock} and E_{soil} are active, E_{soil} is applied first. If E_{soil}×dt is greater than the soil thickness, the remaining E_{rock} is applied proportionally to the remaining fraction of bedrock. For example, if ${E}_{\mathrm{soil}}\times \mathrm{d}t=\mathrm{0.2}$ m but soil thickness is 0.1 m, then E_{rock} is applied at 50 %. Deposition of sediment by hillslope processes is nonlocal and relies on a transport length approach based on Davy and Lague (2009):
where S_{c} is a critical slope parameter (Roering et al., 1999). If $\mathrm{d}z/\mathrm{d}x\ll {S}_{\mathrm{c}}$, most of the sediments are deposited and we approach the linear side of the equation. When $\mathrm{d}z/\mathrm{d}x\to {S}_{\mathrm{c}}$, most of the sediments go to the receivers as predicted by the nonlinear diffusion. In the case of $(\mathrm{d}z/\mathrm{d}x)>{S}_{\mathrm{c}}$, the process recasts the slope to S_{c}, adding any excess material to Q_{s}. A conceptual difference with Carretier et al. (2016) is that we express volumetric flux rather than flux by unit width. This does not affect the physical behaviour of the process but is more consistent with the rest of our implementation. Q_{s} is modified according to E_{rock}, E_{soil} and D_{hill} and fluxes are distributed to multiple receivers proportionally to the slope.
4.1.4 K and κ coefficients for erosion and sediment transport
The coefficients for hillslope and fluvial erosion or sediment transport – κ_{s}, κ_{r}, K_{r} and K_{s} – are empirical and their value can greatly vary from one site to another (e.g. Harel et al., 2016; Carretier et al., 2016). In streampowerlike models, they are roughly a function of m, n and local conditions. In hillslope diffusion, they are a function of local soil and lithologically driven heterogeneity (Carretier et al., 2018). Because both of these empirical coefficients encompass many processes (e.g. Tucker and Slingerland, 1996; Whipple et al., 2013), we use a common base value for each parameter across the whole landscape, or parts of it. These values can be modulated by local or global heterogeneities. The base values can be estimated with sensitivity analyses of spatially variable weighting coefficients and obtain relevant elevations.
4.1.5 Lacustrine sedimentation
Lake deposition is approximated with a simple draping algorithm. Once the final state of a lake is known (see Sect. 3.3), we calculate the proportion of the lake that can be filled with incoming sediment in each pixel: ${V}_{\mathrm{slake}}/\left({V}_{\mathrm{totlake}}\phantom{\rule{0.125em}{0ex}}{h}_{\mathrm{lake}}\right)$. While simplistic, it serves the purpose of this contribution to be a proof of concept in treating lakes as separate entities and paves the way to more detailed lacustrine processes.
4.2 Application I: considering lakes in longterm landscape evolution
In this first set of experiments, we assess the role of lakes and closed basins in longterm landscape evolution. Earlier work by GarciaCastellanos (2006) and GarciaCastellanos and JiménezMunt (2015) (1D and 2D respectively) already noted that endorheism in LEMs was a function of complex relationships between climate (precipitation, evaporation), tectonics and surface processes. Their experiments highlighted the potential importance of integrating endorheism in longterm and largescale landscape evolution studies. Here, we exploit our the capacity of our method to process lakes in order to assess how it could impact simulation results for a given setting.
We ran three simulations for 10 Myr in an idealized mountain range with a frontal thrust, a foreland, and a normal fault in its hinterland (Fig. 5). A uniform, semiarid, yearly precipitation rate was set at 700 mm. Scenario 1 uses the passive lake solver (Sect. 3.3.1), scenario 2a runs with the depressionaware lake solver (Sect. 3.3.2) and scenario 2b has the depressionaware lake solver with lake evaporation. Figure 5 displays a snapshot of the landscape after 5 Myr as a N–S median profile of median and minimum elevation and median sediment thickness and a hillshaded mapview of Q_{w} and lake extents. Figure 6 shows time series of sediment fluxes escaping the southern border of the model as well as the total volume of deposited sediment over the landscapes.
In scenario 1, an unrealistically deep (−500 m after 5 Myr), underfilled and subsiding basin has formed on the footwall of the normal fault. The main E–W drainage divide migrates significantly to the south (Fig. 5a). Over a total 10 Myr long evolution, the two basins store 4×10^{11} m^{3} of sediments while the exported sediment flux is only mildly impacted by the onset of the normal fault (Fig. 6) and shows steady increase after 2 Myr. In scenario 2a, using a depressionaware lake solver, sedimentation in the internal basin balances off subsidence. A longlived very shallow lake is continuously connected to the foreland via a single river (Fig. 5b). The sediment export through time initially is nearly halved from 4 to 2.5×10^{5} m^{3} yrs^{−1} as the depression grows and stabilizes beyond 3 Myr (Fig. 6). Finally, in scenario 2b, a closed basin forms on the hanging wall of the normal fault (Fig. 5c). Its elevation increases through time and it traps all the incoming sediments and water. It quickly becomes disconnected from the foreland (Fig. 5c). The exported sediment flux is halved from 4 to 2×10^{5} m^{3} yr^{−1}, and increases again slightly and steadily through time (Fig. 6).
In scenario 1, the internal depression is unconditionally connected to the rest of the outlet by the passive lake solver and only fluvial deposition can fill the basin. The subsiding basin surface in Fig. 5a demonstrates fluvial deposition is not efficient enough to balance the subsidence. If the topographic signature of the normal fault is exaggerated, Fig. 6 shows that its sediment flux signature is greatly attenuated (minor drop for 2 Myr). More striking is the steady increase of sediment export, which can be explained by the constant lowering of the internal base level and the steepening of the internal basin. The steeper slopes erode faster and the ever greater volume of sediment is exported to the foreland due to the unconditional rerouting. Ultimately, if scenario 1 ran for longer it would display a meaningless landscape inversion draining to the depocenter of the internal basin and “teleporting” sediments to the model edge.
In scenario 2a, the lake almost constantly outflows, maintaining connectivity to the foreland the whole time. This is due to the large amount of water coming from the basin flanks compared to the accommodation space offered by the lake. The erosion of the outlet is barely impacted by the relatively small amount of water stored in the lake (this point is discussed in greater extent later in the Discussion). The balance between maintaining the connection to the foreland base level but maintaining the ability to trap sediment explains the stability of the basin elevation (Fig. 5b) and sediment export through time (Fig. 6) in equilibrium with the tectonic conditions. This results in very low actual lake depth (<1 m most of the time); however, this needs to be interpreted bearing in mind the time step of our simulation is 1000 years and represents an average of processes in that time span. In reality, this could be translated into patches of migrating but more realistically deep lakes. More sophisticated acknowledgement of lake sediment dynamics such as compaction could also enhance the creation of more realistic lakes (Håkanson, 1982).
Finally, scenario 2b is the only one breaking the connectivity to the rest of the landscapes, effectively simulating a closed basin. Lake evaporation balances water input in the lake and enables a decoupling where the wouldbe outlet of the lake does not receive any water or sediments from the lake, inhibiting its erosion compared to scenario 2a. The absence of outlet for the depression means all sediments are trapped within, explaining the highest volume of sediments stored and the lowest export to the model edges. The elevation of the overall model also rises, and if run for longer, the model would probably reach a steady state where the basin would be eventually captured by a river draining externally. The globally higher elevation and the increase of erosion export through time (Fig. 6) result from the increasing elevation of the internal basin and steepening of the landscape.
4.3 Application II: monitoring the sourcetosink system
This case demonstrates the ability of the CHONK framework to provide finegrained detailed information about provenance in the stratigraphy. LEMs have been widely used to investigate the sourcetosink systems (e.g. Guerit et al., 2019; Yuan et al., 2019; Sharman et al., 2019). One particular need in this context is tracking the provenance and destination of material during their erosion, transport and sedimentation processes. This can be done by tracking discrete individual particles (Carretier et al., 2016), or with a bulk approach (Sharman et al., 2019). The latter usually postprocesses this information by integrating the erosion and sedimentation field. Our approach enables an easy embedding of such information within the cell. Provenance tracking is built in and straightforward. Provenance can be tracked within the stratigraphy and reutilized in later time steps without information loss. We demonstrate the model capabilities with a run similar to scenario 2 from Sect. 4.2, but exhuming a simple plutonlike body of harder rock type in the range. Greater rock strength was simulated with a decrease in erodibility. We refer to the harder rock type as “granite” and the background rock type as “substrate” for simplicity.
We ran the simulation for 10 Myr. Figure 7a illustrates highresolution monitoring of sediments with a granite provenance. Thanks to the 3D cellular system storing this information, it can be retrieved with different resolutions, for example in the full sediment column or in the first 10 m as illustrated in the left and right parts of Fig. 7a. Figure 7b and c display this information in crosssection views which highlight largescale stratigraphic structures. Note that here, the provenance data are displayed as relative proportion instead of absolute volume, but both options are possible. The E–W and N–S cross sections in Fig. 7 illustrate the irregularity in the stratigraphic patterns of deposition as the distributary system sweeps across the foreland.
4.4 Application III: erosivity and erodibility captured by dynamic parameters
While the tracking capabilities open many options to monitor the sourcetosink system, they can also be used to integrate feedbacks between processes and characteristics of the sediment flux. Because tracking is dynamic, the state of the fluxes is always known and it can be used to directly influence the process laws. In the following example, we alter the K coefficients of erosion efficiency in Eqs. (9) and (7) to incorporate a notion of relative strengths between sediment and substrate. We assume that harder tools (e.g. granite) impacting softer bedrock (e.g. mudstone) yield greater river incision than softer tools (e.g. schist) on harder material (e.g. Sklar and Dietrich, 1998; Sklar, 2001; Sklar and Dietrich, 2004). We take advantage of the dynamic parametrization of CHONK to implement a firstorder tool strength principle:
where K_{eff} is the effective erodibility used in the equation, K_{r} is the bedrock erodibility, K_{sed} is the erodibility of the mobile sediment, s is an exponent regulating the sensitivity of the system and K_{i} is the local erodibility factor. K_{sed} is a weighted average proportional to the content of each lithology in the model. We store the proportion of each lithology as detailed in Sect. 3.4.1. K_{eff} encompasses nonlocal effects that K_{r} and K_{s} cannot express. The latter are simply linked to the local condition of the cell and have no information about upstream conditions. This interdependence between the nature of nonlocal sediment flux and local erodibility would not be possible without an integrated approach like the one offered by CHONK.
We ran a modified simulation with an uplifting range and a static foreland, essentially Fig. 7a without the normal fault. We start from steadystate conditions and exhume a simplified granitoid. Figure 8 shows the profile of the main river draining through the granitoid at t=0 and t=3 Myr for an unmodified simulation using Eq. (9) and the tool effect simulation using Eq. (14) to highlight nonlinear and nonlocal effects. The lower effective erodibility of the granite traversed by weaker bedload leads to a steeper stream. With hard tools enhancing incision, the slope of the area downstream of the harder rocks is reduced, which drops the base level and propagates knickpoints up all tributaries. Because this enhanced incision is a function of the quantity of granite in the mobile sediments, its effect fades downstream as more softer sediment is added in the mobile flux, affecting the concavity of the river profile nonlinearly.
In this contribution we explored the potential of a modelling framework that separates landscape topology managed by a processagnostic graph on the one hand, from the processes and fluxes managed by a cellular automata numerical structure on the other. We illustrated how this method is particularly suited to tackle research questions involving multiple interconnected processes in complex environments, for example cases where lakes disturb the fluxes of sediment and water independently.
Our approach is built upon existing contributions (e.g. GarciaCastellanos et al., 2003; Tucker and Hancock, 2010; Braun and Willett, 2013; GarciaCastellanos and JiménezMunt, 2015; Carretier et al., 2016; Anand et al., 2020; Barnes et al., 2021), and we show, with three case examples, how we can address scientific questions that were not straightforward or impossible to answer with previous methods without significant refactoring. The main advantages of our modelling design are that (i) it is built for interoperability between fluxes and parameters and (ii) it allows for finegrained monitoring of fluxes independently of surface laws, making it a prime tool for sourcetosink and other sedimentological or stratigraphic studies. We illustrated this interoperability with the simulation of a simple “tool” effect where upstream sediment nature and provenance (from any processes) influence fluvial erosivity. Crossing it with graph theories enables full and efficient control of topology independently of the process and fluxes simulated, even in regions where imbrications of local minima complicate it significantly. Whether lakes, endorheic basins or insignificant noise, our method can process local minima with a lot of flexibility depending on the case study. They can be treated as fully separated domains with dedicated process laws and/or act as partial or full trap in the sourcetosink sediment and water routine. Local minima can also be simply rerouted to ensure flow continuity without affecting computing performance or requiring dedicated processes. Local minima are often overlooked or bypassed, and we demonstrated that the way they are integrated into the model (i) significantly impacts the simulated landscape evolution and (ii) can be fully separated from the surface processes implemented in the LEM. The main breakthrough is the generic processing of these closed domains independent of process laws encouraging seamless integration within LEMs.
The dynamic nature of the model also enables advanced monitoring of fluxes. We illustrated how the pointtracking of sediment provenance and storage in the stratigraphy can inform process laws, whereas existing models commonly postprocess that information from an erosion field. While we focused on the provenance, this opens a wide range of possibilities linked to any information that can be tracked in the cells. For example, one could extend this provenance information to geochemical tracers, or detrital thermochronometers and cosmonuclides (e.g. Petit et al., 2023). In the end, a tracker only needs to be associated with its transporting flux whether hillslope or fluvial sediment or water. Another field of possibility is the tracking of more indirect properties, such as residence time, which are crucial to model luminescence and cosmogenic signals.
The model prototype is opensource, and updated information is available in a github repository (https://github.com/bgailleton/CHONK, last access: 20 December 2023) with instructions on usage, installation and updates. The exact version used in this paper has been archived in https://doi.org/10.5281/zenodo.7746465 (Gailleton, 2023). No data sets were used in this article.
BG was responsible for the software development, BG and LM designed and tested the method and GC and JB provided advice on the numerical methods. BG wrote the manuscript with the help of LM. JB and GC greatly contributed to improve the final version of the text.
The contact author has declared that none of the authors has any competing interests.
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.
We thank Benoit Bovy for his helpful comments on the method. We thank Sébastien Carretier, Kerry L. Callaghan and Andrew D. Wickert for their constructive and helpful reviews that greatly improved the manuscript.
The article processing charges for this openaccess publication were covered by the Helmholtz Centre Potsdam – GFZ German Research Centre for Geosciences.
This paper was edited by Andrew Wickert and reviewed by Sebastien Carretier and Kerry Callaghan.
Adams, J. M., Gasparini, N. M., Hobley, D. E. J., Tucker, G. E., Hutton, E. W. H., Nudurupati, S. S., and Istanbulluoglu, E.: The Landlab v1.0 OverlandFlow component: a Python tool for computing shallowwater flow across watersheds, Geosci. Model Dev., 10, 1645–1663, https://doi.org/10.5194/gmd1016452017, 2017. a
Anand, S. K., Hooshyar, M., and Porporato, A.: Linear layout of multiple flowdirection networks for landscapeevolution simulations, Environ. Modell. Softw., 133, 104804, https://doi.org/10.1016/j.envsoft.2020.104804, 2020. a, b, c
Armitage, J. J.: Short communication: flow as distributed lines within the landscape, Earth Surf. Dynam., 7, 67–75, https://doi.org/10.5194/esurf7672019, 2019. a, b
Babault, J., Bonnet, S., Crave, A., and Van Den Driessche, J.: Influence of piedmont sedimentation on erosion dynamics of an uplifting landscape: An experimental approach, Geology, 33, 301–304, https://doi.org/10.1130/G21095.1, 2005. a
Barnes, R., Lehman, C., and Mulla, D.: An efficient assignment of drainage direction over flat surfaces in raster digital elevation models, Comput. Geosci., 62, 128–135, https://doi.org/10.1016/j.cageo.2013.01.009, 2014a. a, b, c, d
Barnes, R., Lehman, C., and Mulla, D.: Priorityflood: An optimal depressionfilling and watershedlabeling algorithm for digital elevation models, Comput. Geosci., 62, 117–127, https://doi.org/10.1016/j.cageo.2013.04.024, 2014b. a, b, c
Barnes, R., Callaghan, K. L., and Wickert, A. D.: Computing water flow through complex landscapes – Part 2: Finding hierarchies in depressions and morphological segmentations, Earth Surf. Dynam., 8, 431–445, https://doi.org/10.5194/esurf84312020, 2020. a, b, c, d, e, f, g, h, i, j, k, l
Barnes, R., Callaghan, K. L., and Wickert, A. D.: Computing water flow through complex landscapes – Part 3: Fill–Spill–Merge: flow routing in depression hierarchies, Earth Surf. Dynam., 9, 105–121, https://doi.org/10.5194/esurf91052021, 2021. a, b, c, d, e, f, g, h
Barnhart, K. R., Glade, R. C., Shobe, C. M., and Tucker, G. E.: Terrainbento 1.0: a Python package for multimodel analysis in longterm drainage basin evolution, Geosci. Model Dev., 12, 1267–1297, https://doi.org/10.5194/gmd1212672019, 2019. a
Barnhart, K. R., Hutton, E. W. H., Tucker, G. E., Gasparini, N. M., Istanbulluoglu, E., Hobley, D. E. J., Lyons, N. J., Mouchene, M., Nudurupati, S. S., Adams, J. M., and Bandaragoda, C.: Short communication: Landlab v2.0: a software package for Earth surface dynamics, Earth Surf. Dynam., 8, 379–397, https://doi.org/10.5194/esurf83792020, 2020. a, b, c, d, e, f
Bates, P. D., Horritt, M. S., and Fewtrell, T. J.: A simple inertial formulation of the shallow water equations for efficient twodimensional flood inundation modelling, J. Hydrol., 387, 33–45, https://doi.org/10.1016/j.jhydrol.2010.03.027, 2010. a
Bovy, B.: fastscapelem/fastscape: v0.1.0alpha, Zenodo [data set], https://doi.org/10.5281/zenodo.3479426, 2019. a, b, c, d, e, f, g, h, i
Braun, J. and Sambridge, M.: Modelling landscape evolution on geological time scales: a new method based on irregular spatial discretization, Basin Res., 9, 27–52, https://doi.org/10.1046/j.13652117.1997.00030.x, 1997. a, b, c, d
Braun, J. and Willett, S. D.: A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution, Geomorphology, 180–181, 170–179, https://doi.org/10.1016/j.geomorph.2012.10.008, 2013. a, b, c, d, e, f, g, h, i
Bufe, A., Burbank, D. W., Liu, L., Bookhagen, B., Qin, J., Chen, J., Li, T., Thompson Jobe, J. A., and Yang, H.: Variations of Lateral Bedrock Erosion Rates Control Planation of Uplifting Folds in the Foreland of the Tian Shan, NW China, J. Geophys. Res., 122, 2431–2467, 2017. a
Campforts, B., Schwanghart, W., and Govers, G.: Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model, Earth Surf. Dynam., 5, 47–66, https://doi.org/10.5194/esurf5472017, 2017. a, b, c
Campforts, B., Shobe, C. M., Steer, P., Vanmaercke, M., Lague, D., and Braun, J.: HyLands 1.0: a hybrid landscape evolution model to simulate the impact of landslides and landslidederived sediment on landscape evolution, Geosci. Model Dev., 13, 3863–3886, https://doi.org/10.5194/gmd1338632020, 2020. a
Carretier, S., Martinod, P., Reich, M., and Godderis, Y.: Modelling sediment clasts transport during landscape evolution, Earth Surf. Dynam., 4, 237–251, https://doi.org/10.5194/esurf42372016, 2016. a, b, c, d, e, f, g, h, i
Carretier, S., Goddéris, Y., Martinez, J., Reich, M., and Martinod, P.: Colluvial deposits as a possible weathering reservoir in uplifting mountains, Earth Surf. Dynam., 6, 217–237, https://doi.org/10.5194/esurf62172018, 2018. a
Clift, P. D. and Giosan, L.: Sediment fluxes and buffering in the postglacial Indus Basin, Basin Res., 26, 369–386, 2014. a
Cordonnier, G., Bovy, B., and Braun, J.: A versatile, linear complexity algorithm for flow routing in topographies with depressions, Earth Surf. Dynam., 7, 549–562, https://doi.org/10.5194/esurf75492019, 2019. a, b, c, d, e, f, g, h, i, j
Coulthard, T., Neal, J., Bates, P., Ramirez, J., De Almeida, G., and Hancock, G.: Integrating the LISFLOODFP 2D hydrodynamic model with the CAESAR model: implications for modelling landscape evolution, Earth Surf. Proc. Land., 38, 1897–1906, https://doi.org/10.1002/esp.3478, 2013. a, b, c, d
Croissant, T., Lague, D., Steer, P., and Davy, P.: Rapid postseismic landslide evacuation boosted by dynamic river width, Nat. Geosci., 10, 680–684, https://doi.org/10.1038/ngeo3005, 2017. a
D'Ambrosio, D., Di Gregorio, S., Gabriele, S., and Gaudio, R.: A Cellular Automata model for soil erosion by water, Phys. Chem. Earth Pt. B, 26, 33–39, https://doi.org/10.1016/S14641909(01)850115, 2001. a
Davy, P. and Lague, D.: Fluvial erosion/transport equation of landscape evolution models revisited, J. Geophys. Res.Sol. Ea., 114,, F03007, https://doi.org/10.1029/2008JF001146, 2009. a, b
Davy, P., Croissant, T., and Lague, D.: A precipiton method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities, J. Geophys. Res.Earth, 122, 1491–1512, https://doi.org/10.1002/2016JF004156, 2017. a
Dingle, E. H., Sinclair, H. D., Venditti, J. G., Attal, M., Kinnaird, T. C., Creed, M., Quick, L., Nittrouer, J. A., and Gautam, D.: Sediment dynamics across gravelsand transitions: Implications for river stability and floodplain recycling, Geology, 48, 468–472, https://doi.org/10.1130/G46909.1, 2020. a
Finnegan, N. J., Sklar, L. S., and Fuller, T. K.: Interplay of sediment supply, river incision, and channel morphology revealed by the transient evolution of an experimental bedrock channel, J. Geophys. Res.Earth, 112, F03S11, https://doi.org/10.1029/2006JF000569, 2007. a
Gailleton, B.: fastscapelem/fastscapelitho: fastscapelitho 0.0.1, Zenodo [code], https://doi.org/10.5281/zenodo.4773791, 2021. a
Gailleton, B.: CHONK 1.0, Zenodo [code], https://doi.org/10.5281/zenodo.7746465, 2023. a
Gailleton, B., Sinclair, H. D., Mudd, S. M., Graf, E. L. S., and Mațenco, L. C.: Isolating Lithologic Versus Tectonic Signals of River Profiles to Test Orogenic Models for the Eastern and Southeastern Carpathians, J. Geophys. Res.Earth, 126, e2020JF005970, https://doi.org/10.1029/2020JF005970, 2021. a
Ganti, V., Straub, K. M., FoufoulaGeorgiou, E., and Paola, C.: Spacetime dynamics of depositional systems: Experimental evidence and theoretical modeling of heavytailed statistics, J. Geophys. Res.Earth, 116, https://doi.org/10.1029/2010JF001893, 2011. a
GarciaCastellanos, D.: Longterm evolution of tectonic lakes: Climatic controls on the development of internally drained basins, Geological Society of America, ISBN 9780813723983, https://doi.org/10.1130/2006.2398(17), 2006. a, b, c
GarciaCastellanos, D. and JiménezMunt, I.: Topographic Evolution and Climate Aridification during Continental Collision: Insights from Computer Simulations, PLOS ONE, 10, e0132252, https://doi.org/10.1371/journal.pone.0132252, 2015. a, b, c, d, e, f, g
GarciaCastellanos, D., Vergés, J., GasparEscribano, J., and Cloetingh, S.: Interplay between tectonics, climate, and fluvial transport during the Cenozoic evolution of the Ebro Basin (NE Iberia), J. Geophys. Res.Sol. Ea., 108, 2347, https://doi.org/10.1029/2002JB002073, 2003. a, b
Geurts, A. H., Cowie, P. A., Duclaux, G., Gawthorpe, R. L., Huismans, R. S., Pedersen, V. K., and Wedmore, L. N. J.: Drainage integration and sediment dispersal in active continental rifts: A numerical modelling study of the central Italian Apennines, Basin Res., 30, 965–989, https://doi.org/10.1111/bre.12289, 2018. a, b
Grieve, S. W. D., Mudd, S. M., Milodowski, D. T., Clubb, F. J., and Furbish, D. J.: How does gridresolution modulate the topographic expression of geomorphic processes?, Earth Surf. Dynam., 4, 627–653, https://doi.org/10.5194/esurf46272016, 2016. a
Guerit, L., Métivier, F., Devauchelle, O., Lajeunesse, E., and Barrier, L.: Laboratory alluvial fans in one dimension, Phys. Rev. E, 90, 022203, https://doi.org/10.1103/PhysRevE.90.022203, 2014. a
Guerit, L., Yuan, X.P., Carretier, S., Bonnet, S., Rohais, S., Braun, J., and Rouby, D.: Fluvial landscape evolution controlled by the sediment deposition coefficient: Estimation from experimental and natural landscapes, Geology, 47, 853–856, https://doi.org/10.1130/g46356.1, 2019. a
Håkanson, L.: Bottom dynamics in lakes, in: Sediment/Freshwater Interaction, edited by: Sly, P. G., Developments in Hydrobiology, Springer Netherlands, Dordrecht, 9–22, ISBN 9789400980099, https://doi.org/10.1007/9789400980099_2, 1982. a
Harel, M. A., Mudd, S. M., and Attal, M.: Global analysis of the stream power law parameters based on worldwide ^{10}Be denudation rates, Geomorphology, 268, 184–196, https://doi.org/10.1016/j.geomorph.2016.05.035, 2016. a, b
Hergarten, S.: Transportlimited fluvial erosion – simple formulation and efficient numerical treatment, Earth Surf. Dynam., 8, 841–854, https://doi.org/10.5194/esurf88412020, 2020. a, b, c, d, e
Howard, A. D. and Kerby, G.: Channel changes in badlands, GSA Bulletin, 94, 739–752, https://doi.org/10.1130/00167606(1983)94<739:CCIB>2.0.CO;2, 1983. a
Jerolmack, D. J. and Sadler, P.: Transience and persistence in the depositional record of continental margins, J. Geophys. Res.Earth, 112, F03S13, https://doi.org/10.1029/2006JF000555, 2007. a
Jyotsna, R. and Haff, P. K.: Microtopography as an indicator of modern hillslope diffusivity in arid terrain, Geology, 25, 695–698, https://doi.org/10.1130/00917613(1997)025<0695:MAAIOM>2.3.CO;2, 1997. a
Callaghan, K. L. and Wickert, A. D.: Computing water flow through complex landscapes – Part 1: Incorporating depressions in flow routing using FlowFill, Earth Surf. Dynam., 7, 737–753, https://doi.org/10.5194/esurf77372019, 2019. a, b
Leonard, J. S. and Whipple, K. X.: Influence of Spatial Rainfall Gradients on River Longitudinal Profiles and the Topographic Expression of Spatially and Temporally Variable Climates in Mountain Landscapes, J. Geophys. Res.Earth, 126, e2021JF006183, https://doi.org/10.1029/2021JF006183, 2021. a
Lindsay, J. B.: Efficient hybrid breachingfilling sink removal methods for flow path enforcement in digital elevation models, Hydrol. Process., 30, 846–857, https://doi.org/10.1002/hyp.10648, 2016. a
Lupker, M., Lavé, J., FranceLanord, C., Christl, M., Bourlès, D., Carcaillet, J., Maden, C., Wieler, R., Rahman, M., Bezbaruah, D., and Xiaohan, L.: ^{10}Be systematics in the TsangpoBrahmaputra catchment: the cosmogenic nuclide legacy of the eastern Himalayan syntaxis, Earth Surf. Dynam., 5, 429–449, https://doi.org/10.5194/esurf54292017, 2017. a
Malatesta, L. C. and Avouac, J.P.: Contrasting river incision in north and south Tian Shan piedmonts due to variable glacial imprint in mountain valleys, Geology, 46, 659–662, https://doi.org/10.1130/G40320.1, 2018. a
Malatesta, L. C., Avouac, J.P., Brown, N. D., Breitenbach, S. F. M., Pan, J., Chevalier, M.L., Rhodes, E., SaintCarlier, D., Zhang, W., Charreau, J., Lavé, J., and Blard, P.H.: Lag and mixing during sediment transfer across the Tian Shan piedmont caused by climatedriven aggradation–incision cycles, Basin Res., 30, 613–635, https://doi.org/10.1111/bre.12267, 2018. a
Mudd, S. M., Clubb, F. J., Grieve, S. W. D., Milodowski, D. T., Hurst, M. D., Gailleton, B., and Valters, D. A.: LSDTopoTools2, Zenodo [data set], https://doi.org/10.5281/zenodo.3245041, 2019. a, b
Nagel, K. and Schreckenberg, M.: A cellular automaton model for freeway traffic, J. Phys. I, 2, 2221–2229, https://doi.org/10.1051/jp1:1992277, 1992. a
O'Callaghan, J. F. and Mark, D. M.: The extraction of drainage networks from digital elevation data, Computer Vision, Graphics, & Image Processing, 28, 323–344, https://doi.org/10.1016/S0734189X(84)800110, 1984. a
Paola, C., Straub, K., Mohrig, D., and Reinhardt, L.: The “unreasonable effectiveness” of stratigraphic and geomorphic experiments, EarthSci. Rev., 97, 1–43, https://doi.org/10.1016/j.earscirev.2009.05.003, 2009. a
Perron, J. T.: Numerical methods for nonlinear hillslope transport laws, J. Geophys. Res.Earth, 116, F02021, https://doi.org/10.1029/2010JF001801, 2011. a
Petit, C., Salles, T., Godard, V., Rolland, Y., and Audin, L.: River incision, ^{10}Be production and transport in a sourcetosink sediment system (Var catchment, SW Alps), Earth Surf. Dynam., 11, 183–201, https://doi.org/10.5194/esurf111832023, 2023. a
Roelvink, J. and Van Banning, G.: Design and development of DELFT3D and application to coastal morphodynamics, Oceanographic Literature Review, 11, 925, https://www.infona.pl/resource/bwmeta1.element.elsevier1ca19bb625b93bf5bfe9e96a7027c553 (last access: 2 January 2024), 1995. a
Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resour. Res., 35, 853–870, https://doi.org/10.1029/1998WR900090, 1999. a, b, c
Sadler, P. M.: Sediment Accumulation Rates and the Completeness of Stratigraphic Sections, J. Geol., 89, 569–584, https://doi.org/10.1086/628623, 1981. a
Salles, T.: eSCAPE: Regional to Global Scale Landscape Evolution Model v2.0, Geosci. Model Dev., 12, 4165–4184, https://doi.org/10.5194/gmd1241652019, 2019. a, b, c
Salles, T., Lopez, S., Cacas, M. C., and Mulder, T.: Cellular automata model of density currents, Geomorphology, 88, 1–20, https://doi.org/10.1016/j.geomorph.2006.10.016, 2007. a
Schumer, R., Taloni, A., and Furbish, D. J.: Theory connecting nonlocal sediment transport, earth surface roughness, and the Sadler effect, Geophys. Res. Lett., 44, 2281–2289, https://doi.org/10.1002/2016GL072134, 2017. a
Schwanghart, W. and Heckmann, T.: Fuzzy delineation of drainage basins through probabilistic interpretation of diverging flow algorithms, Environ. Modell. Softw., 33, 106–113, https://doi.org/10.1016/j.envsoft.2012.01.016, 2012. a
Schwanghart, W. and Kuhn, N. J.: TopoToolbox: A set of Matlab functions for topographic analysis, Environ. Modell. Softw., 25, 770–781, https://doi.org/10.1016/j.envsoft.2009.12.002, 2010. a
Schwanghart, W. and Scherler, D.: Short Communication: TopoToolbox 2 – MATLABbased software for topographic analysis and modeling in Earth surface sciences, Earth Surf. Dynam., 2, 1–7, https://doi.org/10.5194/esurf212014, 2014. a, b
Sharman, G. R., Sylvester, Z., and Covault, J. A.: Conversion of tectonic and climatic forcings into records of sediment supply and provenance, Scientific Reports, 9, 4115, https://doi.org/10.1038/s41598019397546, 2019. a, b, c
Shobe, C. M., Tucker, G. E., and Barnhart, K. R.: The SPACE 1.0 model: a Landlab component for 2D calculation of sediment transport, bedrock erosion, and landscape evolution, Geosci. Model Dev., 10, 4577–4604, https://doi.org/10.5194/gmd1045772017, 2017. a, b
Sklar, D. W.: Sediment and rock strength controls on river incision into bedrock, Geology, 29, 1087–1090, https://doi.org/10.1130/00917613(2001)029<1087:SARSCO>2.0.CO;2, 2001. a
Sklar, L. and Dietrich, W.: River longitudinal profiles and bedrock incision models: Stream power and the influence of sediment supply, Rivers over rock: fluvial processes in bedrock channels, American Geophysical Union, 237–260, Print ISBN 9780875900902, Online ISBN 9781118664292, https://doi.org/10.1029/GM107p0237, 1998. a, b, c
Sklar, L. S. and Dietrich, W. E.: Sediment and rock strength controls on river incision into bedrock, Geology, 29, 1087–1090, 2001. a
Sklar, L. S. and Dietrich, W. E.: A mechanistic model for river incision into bedrock by saltating bed load, Water Resour. Res., 40, W06301, https://doi.org/10.1029/2003WR002496, 2004. a, b
Struth, L., GarcíaCastellanos, D., RodríguezRodríguez, L., ViaplanaMuzas, M., Vergés, J., and JiménezDíaz, A.: Topographic, lithospheric and lithologic controls on the transient landscape evolution after the opening of internallydrained basins. Modelling the North Iberian Neogene drainage, BSGF – Earth Sciences Bulletin, 192, 45, https://doi.org/10.1051/bsgf/2021036, 2021. a
Tarboton, D. G.: A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resour. Res., 33, 309–319, https://doi.org/10.1029/96WR03137, 1997. a
Tofelde, S., Schildgen, T. F., Savi, S., Pingel, H., Wickert, A. D., Bookhagen, B., Wittmann, H., Alonso, R. N., Cottle, J., and Strecker, M. R.: 100 kyr fluvial cutandfill terrace cycles since the Middle Pleistocene in the southern Central Andes, NW Argentina, Earth Planet. Sc. Lett., 473, 141–153, 2017. a
Tofelde, S., Bernhardt, A., Guerit, L., and Romans, B. W.: Times Associated With SourcetoSink Propagation of Environmental Signals During Landscape Transience, Frontiers in Earth Science, 9, 628315, https://doi.org/10.3389/feart.2021.628315, 2021. a
Tucker, G. E. and Bradley, D. N.: Trouble with diffusion: Reassessing hillslope erosion laws with a particlebased model, J. Geophys. Res.Earth, 115, F00A10, https://doi.org/10.1029/2009JF001264, 2010. a
Tucker, G. E. and Hancock, G. R.: Modelling landscape evolution, Earth Surf. Proc. Land., 35, 28–50, https://doi.org/10.1002/esp.1952, 2010. a
Tucker, G. E. and Slingerland, R.: Predicting sediment flux from fold and thrust belts, Basin Res., 8, 329–349, https://doi.org/10.1046/j.13652117.1996.00238.x, 1996. a
Tucker, G. E., Hobley, D. E. J., Hutton, E., Gasparini, N. M., Istanbulluoglu, E., Adams, J. M., and Nudurupati, S. S.: CellLabCTS 2015: continuoustime stochastic cellular automaton modeling using Landlab, Geosci. Model Dev., 9, 823–839, https://doi.org/10.5194/gmd98232016, 2016. a, b
Von Neumann, J.: The general and logical theory of automata, in: Systems Research for Behavioral Science, Routledge, 97–107, 1968. a, b
Wang, L. and Liu, H.: An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling, Int. J. Geogr. Inf. Sci., 20, 193–213, https://doi.org/10.1080/13658810500433453, 2006. a
Whipple, K. X., DiBiase, R. A., and Crosby, B. T.: Bedrock Rivers, in: Treatise on Geomorphology, Fluvial Geomorphology, 9, 550–573, https://doi.org/10.1016/B9780123747396.002542, 2013. a
Wolfram, S.: Cellular automata as models of complexity, Nature, 311, 419–424, https://doi.org/10.1038/311419a0, 1984. a, b
Yuan, X. P., Braun, J., Guerit, L., Rouby, D., and Cordonnier, G.: A New Efficient Method to Solve the Stream Power Law Model Taking Into Account Sediment Deposition, J. Geophys. Res.Earth, https://doi.org/10.1029/2018JF004867, 2019. a, b, c, d