Theory

Basics

The landscape evolution simulations presented on this site have been carried out using the CIDRE code (e.g. Carretier et al., 2016). This website has not been created as an apology for CIDRE. There are other freely usable models with remarkable performances of which a non-exhaustive list is given in the table at the end of this page. All LEMs applied over millions of years have common points (stationarity of water flows in most of the cases, general shape of erosion laws) but also different levels of simplification concerning hydrology and sediment transport which can vary. Some algorithms are frighteningly fast, while others favour the addition of certain processes (e.g. lateral erosion in CIDRE). A feature of CIDRE that exists in few models is its ability to ”spread” water runoff (Multiple flow algorithm instead of steepest descent algorithm), which can therefore diverge and more realistically account for the flow over foothills and plain areas.

In CIDRE, topography is represented by a mesh of square cells (or pixels). Each cell has an altitude and a constant size. Basically, CIDRE is injected with a rain grid or a constant value, an initial topography, a tectonic uplift rate grid or value, an erodibility grid or a value and CIDRE calculates over time the evolution of the topography, the quantity eroded or the quantity of sediments deposited on every cell of the grid.

In practice, over a time step (about 100 years), CIDRE runs the cells in the decreasing direction of the altitudes. It propagates the quantity of water fallen by the rain to calculate a water discharge flowing on each cell, and it also calculates the slopes between this cell and its lower neighbours. From the steepest slope and from the water discharge CIDRE calculates the amount that can be eroded, and the amount of moving sediment from higher up that will be deposited. In addition, the water flow can erode the lateral neighbours cells (bank erosion). Eroded and transiting material is transferred to the lower cells in proportion to the slope. Then CIDRE moves to the cell that is lower in the list and so on until all the cells have been treated. An increment of tectonic uplit is added to elevations. Then a new time step starts and the same process is repeated until the total duration of the model is reached.

Equations

Starting from an initial topography, the modification of the topography proceeds by successive time steps. During a time step, precipitation falls on the grid at a rate P [LT^-1] and a multiple flow algorithm propagates the water flux Q [L³T^-1], toward all downstream cells in proportion to the slope in each direction. Then the elevation z (river bed or hillslope surface) changes on each cell (size dx) according to the balance between erosion ϵ [LT^-1] and deposition D [LT^-1]. The erosion is different for sediment and for bedrock and ϵ is the sum of two values, one corresponding to gravitational processes without involving the runoff, usually dominating on the hillslopes, and another one associated with water discharge, typically dominating in rivers. Water flowing in one direction is also able to detach material from the cells located perpendicular to that direction to simulate river bank erosion. This erosion generates a lateral (bank) sediment discharge Q_sl [L³T^-1] toward the cell where the water is flowing. Finally, elevation changes also by adding an uplift U [LT^-1] (subsidence if negative).

In the following, all the water and sediment fluxes are ”per unit width” [L²T^-1], the width being either a calculated sub-pixel flow width, or the pixel size dx. All the simulations in this website were run using this second option.

The rate of elevation change on a cell is determined by the following mass balance equation (e.g. Davy and Lague, 2009; Carretier et al., 2016; Shobe et al., 2017):

where the subscript ”r” (”river”) denotes rates associated with flowing water and ”h” (”hillslope”) denotes rates that depends only on the topographic gradient or slope S. Then we define a constitutive law for each of these components: (Carretier et al., 2016)

where K [L^1-2mT^m-1], κ [LT^-1] are erodibility parameters, m and n are lithology-dependent (different for bedrock or sediment) erosion parameters, S is the slope, q [L²T^-1] is the water discharge per stream unit width, and

where q_sr and q_sh are the incoming river and hillslope sediment fluxes (total q_s = q_sr + q_sh) per unit width [L²T^-1], ζ is a river transport length parameter [T L^-1] and S_c is a slope threshold. These fluxes are the sum of sediment fluxes leaving upstream neighbour cells while the deposition rates on a cell are a fraction of the incoming sediment.

Concerning the river processes, ϵ_r is known as the stream power law and derives from the assumption that ϵ_r is proportional to a power law of the shear stress applied by the flowing water on the river bed (e.g. Whipple et al., 2000; Lague, 2014). The deposition rate D_r is a fraction of the incoming sediment flux and this fraction (ζq) has the dimension of the reverse of a length. We call this length a transport length because it has the physical meaning of a characteristic distance over which a volume of detached material will transit downstream before being deposited. In particular, when the local q is large, few sediment eroded from upstream will deposit on the cell. The transport length depends on ζ, proportional to the reverse of a settling velocity of sediment in water (e.g. Davy and Lague, 2009; Lajeunesse et al., 2013). In instantaneous river models, ζ should be fixed by the grain size of sediment. In landscape evolution models, where the water discharge q averages the periods with and without transport, ζ is an ”apparent” parameter that can take a large range of values in real situations depending on climate variability (Guerit et al., 2019).

Concerning the hillslopes processes, the philosophy is the same, except that the detachment rate ϵ_r and the deposition rate D_r depend only on the slope. The linear slope dependence of ϵ_r describes diffusion processes. D_r depends on a specified critical slope S_c: when the slope is close to S_c, the deposition rate D_r decreases rapidly, simulating in average the onset of shallow landslides. The transport length associated with gravitational processes () is inversely proportional to the probability to deposit sediment on the cell. This erosion-deposition formulation leads to similar solutions as the critical slope-dependent hillslope model studied for example by Roering et al. (1999) (Carretier et al., 2016).

Flowing water in each direction can erode lateral cells perpendicular to that direction. Little is known about the law that describes the widening rate of valley, and establishing a lateral erosion law suitable for landscape evolution models averaging processes over millennia is a challenge (Langston and Tucker, 2018; Langston and Temme, 2019). Here, the lateral sediment flux per unit length q_sl [L²T^-1] eroded from a lateral cell is simply defined as a fraction of the river sediment flux q_sr [L²T^-1] in the considered direction (e.g. Murray and Paola, 1997; Nicholas and Quine, 2007), assuming that lateral mobility of channels, and thus lateral erosion, increases with the flux of river sediment (Bufe et al., 2016, 2019):

where α is a bank erodibility coefficient. α is specified for loose material (sediment) and is implicitly determined for bedrock layers such that the ratio of the lateral erodabilities is equal to the ratio of the fluvial ones (α_loose∕α_bedrock = K_loose∕K_bedrock , with K from Equation 2). If sediment covers the bedrock of a lateral cell, α is weighted by its respective thickness above the target cell.

Finally, the sediment leaving a cell is spread in the same way as water, i.e. proportionally to the downstream slopes. This procedures starts from the most elevated cell and ends with the lowest cell and is repeated in the next time steps until the end of the specified model time (Myr in our case).

Parameter values for the reference simulation

The reference simulation is a constantly uplifting domain grid of 40x40 km² (200x200 cells of size 200 m).

Associating erodibility parameters to a particular rock is challenging. The best that we can do in the current state of our knowledge is to fix the erodibility parameters so that the relief is realistic. We impose here only one bedrock type for which the erodibility is K = 10^-4 m^-0.5 yr^-1, and the water discharge and slope exponent are m = 0.5, n = 1. The bedrock erodibility by pure gravitational processes is κ = 10^-4 m yr^-1. For the sediment (previously deposited but that can be reeroded), we use K = 10^-3 m^-0.5 yr^-1, m = 0.5, n = 1 and κ = 10^-4 m yr^-1.

The transport length parameter ζ is set to 1 yr m^-1, and corresponds to a low value for natural systems (median at 17 yr/m Guerit et al., 2019). It is difficult to link ζ with physical properties of sediment because ζ changes according to the variability of transport periods, but low values seem to correspond to temperate perennial rivers (Guerit et al., 2019).

The lateral erosion parameter α to 0. Finally, the critical slope is S_c =tan(40^o). The northern southern sides are open and fixed to z = 0 m. Periodic boundary conditions are imposed on the two other sides meaning that water and sediment leaving on the one side is reinjected on the other one.

U is fixed to 10^-3 m yr^-1 and P to 1 m yr^-1.

The final time is fixed to 7 Ma for the reference simulation, a duration allowing the topography to reach a dynamic equilibrium, ie. a state where the topography does not change in average and the erosion balances uplift. For some simulations, this time is longer because the topography takes more time to reach a dynamic equilibrium.

Existing LEMs

---

---

References