Modelling sediment clasts transport during landscape evolution

. Over thousands to millions of years, the landscape evolution is predicted by models based on ﬂuxes of eroded, transported and deposited material. The laws describing these ﬂuxes, corresponding to averages over many years, are difﬁcult to prove with the available data. On the other hand, sediment dynamics are often tackled by studying the distribution of certain grain properties in the ﬁeld (e.g. heavy metals, detrial zircons, 10Be in gravel, magnetic tracers, etc.). There is a gap 5 between landscape evolution models based on ﬂuxes and these ﬁeld data on individual clasts, which prevent the latter from being used to calibrate the former. Here we propose an algorithm coupling the landscape evolution with mobile clasts. Our landscape evolution model predicts local erosion, deposition and transfer ﬂuxes resulting from hillslope and river processes. Clasts of any size are initially spread in the basement and are detached, moved and deposited according to probabilities using these 10 ﬂuxes. Several river and hillslope laws are studied. Although the resulting mean transport rate of the clasts does not depend on the time step or the model cell size, our approach is limited by the fact that their scattering rate is cell-size dependent. Nevertheless, both their mean transport rate and the shape of the scattering-time curves ﬁt the predictions. Different erosion-transport laws generate different clast movements. These differences show that studying the tracers in the ﬁeld may provide a way 15 to establish these laws on the hillslopes and in the rivers. Possible applications include the interpretation of cosmogenic nuclides in individual gravel deposits, provenance analyses, placers, sediment coarsening or ﬁning, the relationship between magnetic tracers in rivers and the river planform, and the tracing of weathered sediment.


Introduction
Numerical models of landscape evolution have significantly improved our understanding of relief dynamics by recasting competing theories within a general framework (e.g.Kooi and Beaumont, 1994).Within this framework, the relief dynamics are determined by the balance between detached, deposited and transferred sediment fluxes under the influence of tectonics and climate (e.g.Tucker and Hancock, 2010).These fluxes are set by constitutive laws, such as the diffusion law for hillslope creep (Culling, 1960) or the streampower law for sediment river transport (Howard and Kerby, 1983).Nevertheless, it remains a challenge to justify these laws in geomorphology.Different laws produce different topographic evolutions, which have been widely studied (e.g.Howard, 1997;Willgoose et al., 1991;Kooi andBeaumont, 1994, 1996;Whipple and Tucker, 2002;Tucker and Whipple, 2002;Carretier et al., 2009;Davy and Lague, 2009).These laws have been justified or calibrated by mechanics (Whipple et al., 2000) or experiments (Sklar and Dietrich, 2001) or by comparing their predictions with natural river profiles (e.g.Stock and Montgomery, 1999;Lague and Davy, 2003;Carretier et al., 2006;Loget et al., 2006), which has not always Published by Copernicus Publications on behalf of the European Geosciences Union.
resulted in a good agreement (e.g.Tomkin et al., 2003;van der Beek and Bishop, 2003).The erosion-transport laws in landscape evolution models necessarily apply over long time steps ( 1 year) and relatively large model cells ( 10 m).The uncertainty on these laws comes from a lack of methods to directly quantify sediment fluxes over such time and spatial scales.
At the same time, there are many techniques to trace the provenance and transport of rock fragments (clasts) and minerals.For example, detrital zircons, heavy minerals or trace elements in sedimentary rocks and river streams are routinely used to determine sedimentary provenance and/or constrain the exhumation history of orogenic highlands (e.g.Roddaz et al., 2005;Rodríguez et al., 2012).To better understand the dynamics of sediment transport in rivers, other approaches have used painted cobbles (e.g.Church and Hassan, 1992), implanted magnets in cobbles (e.g.Hassan et al., 1991;Haschenburger, 2011), magnetic iron slag as a tracer (Houbrechts et al., 2011) or radio-frequency identification of passive tracers (e.g.Bradley and Tucker, 2012).In parallel, river dynamics have been explored tracing grains in experimental devices (e.g.Lajeunesse et al., 2010;Kasprak et al., 2015).Furthermore, U-series disequilibrium studies in sediment fractions and cosmogenic nuclide measurements in individual clasts provide relevant information about erosion processes and rates (Chabaux et al., 2006;Codilean et al., 2008;Gayer et al., 2008;McPhillips et al., 2014).Many of these methods provide information about individual clasts that are part of bulk material fluxes over geological timescales.Therefore, we hypothesise that this information may be used to derive erosion-transport laws used in landscape evolution models.The link between erosion-transport laws and the spatial evolution of a clast population during landscape evolution requires a model that couples fluxes, topography and clasts.
This type of model would be key to quantitatively link the statistics for provenance tracers with erosion rates in a catchment.For example, Nibourel et al. (2015) recently showed the potential of coupling the Raman spectroscopy analysis on carbonaceous material from detrital sands with a bedrock abrasion model to localise the highest erosion areas in a catchment in New Zealand.More generally, there is a link between zones that erode rapidly, the transport distance of their erosion byproducts, and the topography, climate and base level at a particular period of time.Thus, the coupling of the analysis of tracers ages and distribution such as detrital zircons or apatite with a landscape evolution model could be used to better constrain tectonic or climatic variations.Other issues concerning the relationship between erosion, weathering and CO 2 consumption require a grain-flux model.The efficiency of silicate weathering to regulate atmospheric CO 2 depends on the residence time of the sediment through the mountain-basin system (Anderson et al., 2013;Mudd and Yoo, 2010).Different approaches have led to contradictory mean residence times in the different reservoirs, particularly in the Ganga Plain (Jain and Tandon, 2010;Hoffmann, 2015).A model that couples clasts and landscape evolution may help better understand how sediment mixing by lateral erosion, avulsions, etc. influence the mean residence time of the particles, and may reconcile different estimations (Davy and Lague, 2012).
In order to develop such a model, we couple a landscape evolution model with a clast dispersion model.The landscape evolution model is a modified version of Cidre (e.g.Carretier et al., 2009), which belongs to the family of reducedcomplexity models (Murray, 2007).The clasts move according to probabilities depending on the erosion, deposition and transport fluxes calculated in Cidre.Our main goal here is to show that this algorithm yields a clast population movement that is consistent with the predicted sediment flux in some simple hillslope and river cases.
After briefly reporting previous modelling approaches based on flux-particle duality, we present Cidre and the probabilities used to move clasts.Then we analyse clast movement in the restricted cases of hillslope diffusion and river transport.Finally, we discuss the potential applications of this model.They include the 3-D tracing of weathered material which initially motivated this modelling approach.

Previous models coupling fluxes and particles
Models coupling fluxes and particles have been developed in other scientific fields, in particular in fluid mechanics, and are known as smoothed particle hydrodynamics (SPH) models (Gingold and Monaghan, 1977).The philosophy of these models is to simplify fluids to a limited set of discrete fluid particles moving under a force field.The distribution of the particles at the next time step is then used to estimate the new field of quantities such as forces, temperatures, densities or water flow.In geomorphology, the original landscape evolution model introduced by Mitasova et al. (2004) in the GIS software GRASS belongs to this class of models.In addition, authors used moving grains to predict topographical variations on hillslopes and in rivers.For example, Tucker and Bradley (2010) modelled the morphological evolution of 2-D scarps using moving grains, while Malmaeus and Hassan (2002) and Lajeunesse et al. (2013) among others used grains to model bedload transport.Grain transport modelling has also been carried out to predict the statistical distribution of certain geochemical properties.Within that scope, Repka et al. (1997), Gayer et al. (2008), Codilean et al. (2008), Carretier et al. (2009) and Carretier and Regard (2011) developed approaches in which grains move over a static landscape in order to trace their cosmogenic nuclide concentrations.
The modelling approach presented in this paper is different from these published works in the sense that (1) particles are not used to estimate the water or erosion fields, (2) the topography evolves over time in our model, and (3) our modelling is 3-D instead of 2-D.Nevertheless, our model is inspired by the coupling of the landscape evolution model Eros and sediment particles (Davy and Lague, 2012).

Erosion-sedimentation in Cidre
Cidre is a C++ code modelling the topography dynamics on a regular grid of square cells.At the beginning of a time step, a specified volume of rain falls.Cells are sorted by decreasing elevations.The propagation of water and sediment is proceeded in cascade starting from the highest cell to ensure mass conservation.A multiple-flow algorithm propagates the water flux Q [L 3 /T ] toward all the downstream cells proportionally to the slope in each direction.

Mass balance
The elevation z (river bed or hillslope surface) changes on each cell according to the two following equations (Fig. 1): and we define where is a local erosion (detachment or entrainment) rate [L/T ], D is a local deposition rate [L/T ], U is the uplift or subsidence rate [L/T ], q s is the incoming sediment flux per unit width [L 2 /T ] and L is transport distance.L ≥ dx to keep the deposition flux D dx smaller than the incoming sediment flux q s .For a complete justification of these equations in the case of river transport, see Charru (2006), Davy and Lague (2009) or Lajeunesse et al. (2015).We recall some of the elements in the Appendix.Here, we generalise this approach for both hillslope and river sediment transport processes by specifying and L in both cases.
The length L determines the proportion of incoming sediment flux which is deposited.A large L means that the deposition is small, which is favoured in a real-case scenario by steep slope or high water discharge.The cell outflux per unit width q s results from the sum of the sediment detached from this cell and the sediment eroded from upstream and which then crossed the cell without depositing.The sediment flux q s is thus non-local (e.g.Tucker and Bradley, 2010;Lajeunesse et al., 2015).

Hillslope
In the following we establish equations for and L first for hillslope processes and then for rivers.As the hillslope model is new, we begin by recalling how long-term ( 1 kyr) hillslope evolution has been previously simulated using a "non-linear" diffusion model (e.g.Carretier et al., 2009;Perron, 2011;Carretier et al., 2014).This "non-linear" diffusion model has been proposed by different authors (e.g.Andrews and Hanks, 1985;Hanks, 1999;Roering et al., 1999) and is supported by 10 Be-derived erosion rates (e.g.Binnie et al., 2007) or experiments (Roering et al., 2001).It is usually presented in the following form (e.g.Roering et al., 1999): where description is thus based on the definition of a flux of transported sediment parallel to the downward slope.When the slope is small, this flux refers to diffusion processes such as soil creep, rain splash or diffuse runoff.This flux increases dramatically when the slope gets closer to the critical slope, simulating on average the cumulative effect of mass wasting events.
Here we use a different approach where the elevation variation results from the difference between a local detachment rate and a deposition rate using Eqs.(1) and (2).In these equations, we specify and the transport length L by where κ [L/T ] is an erodibility coefficient.If the slope S > S c then is set such that S = S c .The detachment rate is proportional to the local gradient.However, the deposition rate (q s /L Eq. 2) depends on the slope and the critical slope: when S S c , most of the sediment entering a cell is deposited on this cell.This is the pure diffusion case.The sediment flux q s does not include sediment eroded from above and is thus "local" (Furbish and Haff, 2010).When S ∼ S c , L becomes infinity and there is no redeposition on the cell.This behaviour corresponds to mass wasting where a grain eroded from a place is able to travel a large distance before stopping.In that case, the flux q s is "non-local" as it incorporates not only sediment that has been detached locally but also sediment in transit that has eroded from above.For an intermediate S, there is a progressive transition between pure creep and "ballistic" transport of material through the hillslope, which seems consistent with experiments (Roering et al., 2001;Gabet and Mendoza, 2012).
Despite these conceptual differences, Eqs. ( 5) and ( 6) predict similar topographic evolutions to the "non-linear" diffusion equations for κ = κ dx.This similarity is illustrated by Fig. 2, which displays the evolution of a transverse hill profile during uplift and relief decline.In order to solve Eqs. ( 3) and ( 4), we use an explicit finite-difference approach along a topographic profile.In order to avoid numerical instability when S is close to S c , we approximate Eq. ( 4) with its linear approximation when S > 0.99S c (e.g.Carretier et al., 2009).Equations ( 5) and ( 6) do not have this stability problem because L is simply set at a huge number when S > S c .Fig. 2 shows that both sets of Eqs.(5, 6) and (3, 4) lead to the same evolution.
The red profiles represent 2-D averaged elevations across 3-D topographies.Note that S c is slightly larger in the "non-linear" model for stability reasons; the non-linear flux in Eq. ( 4) is approximated by the tangent to the curve q s = f (S) when S is close to S c , resulting in a small underestimation of q s .(a) Relaxation from a triangular hill with a slope of 0.5.Profiles every 2.5 Myr.(b) Uplift (0.1 mm yr −1 ) from an initial horizontal surface.Profiles every 0.5 Myr.The elevation is fixed at 0 m at the two boundary conditions in (a) and (b).
They may reconcile the different views concerning the "diffusive" or "non-local" nature of erosion on hillslopes as they predict a progressive transition between local and linear to non-local and non-linear sediment flux when the slope increases (e.g.Roering et al., 2001;van Milligen and Bons, 2002;Roering et al., 2002;Furbish and Haff, 2010;Foufoula-Georgiou et al., 2010;Tucker and Bradley, 2010;Gabet and Mendoza, 2012;Falcini et al., 2013;Furbish and Roering, 2013).These equations are simple to implement in a 3-D model, particularly when there are different sediment and bedrock layers because the erosion of these layers is separated from the sediment passing through the cell or deposited on it.We have not carried out an extensive comparison between the different resolution schemes of Eqs.(3, 4) (for example, see the solutions provided in Perron, 2011) and (5, 6).Nevertheless, the solution of Eqs. ( 5) and ( 6) that we propose remains stable for time steps that are 1 order of magnitude larger than with Eqs. ( 3) and (4) in the examples given in Fig. 2.

Rivers
In the case of river processes, we describe here a simplified version of material detachment (sediment or bedrock), although the detachment threshold and the explicit expression Earth Surf.Dynam., 4, 237-251, 2016 www.earth-surf-dynam.net/4/237/2016/ of bed shear stress in particular can be included (e.g.Tucker, 2004): where K is an erodibility coefficient, q [L 2 /T ] the water discharge per unit flow width on the cell, S is the steepest slope and the exponents are positive.The transport length L comes from the derivation of Davy and Lague (2009) for saltation and ξ [T /L] is a factor depending on particle size and density.This law implies that the deposition rate decreases when the water discharge per unit width q increases.As demonstrated by Charru (2006), Davy and Lague (2009) or Lajeunesse et al. (2015), Eqs. ( 1) and ( 2) are mathematically equivalent to the "under-capacity" or "saturation" transport model (e.g.Beaumont et al., 1992;Whipple and Tucker, 2002;Andreotti et al., 2010), although L has a different physical meaning in both cases.It is a transport length in the Charru (2006) or Davy and Lague (2009) formulation (e.g. a characteristic transport distance of a population of grains entering a river locally by landsliding) and a saturation length to reach a maximum or "equilibrium" sediment flux, also called transport capacity in the Beaumont et al. (1992) formulation (e.g. the distance required for the suspended sediment flux to reach a maximum value downstream an abrupt transition between the bedrock and alluvial river bed).As detailed by Davy and Lague (2009), the notion of transport capacity emerges by rewriting Eqs. ( 1) and (2) as ∂z ∂t = q s − L L , which is the form of the under-capacity model with L defining a transport capacity.Thus, the formulation adopted in Cidre is equivalent to the under-capacity model.
For river processes, the flow width w can be set to the cell width dx or a river width.In the latter case, the flow width w is where k w is a coefficient depending on the lithology and Q [L 3 /T ] is the total water discharge at a river section.In this case, cell erosion in Eq. 1 is recast as w dx (finite-volume approach).
Erosion for sediment is different from that of bedrock (Eqs.5 and 7), and within bedrock, different layers can be defined by their erodibility and detachment or slope thresholds (κ and S c for hillslope processes and K for river processes).On one cell and during a time step dt, different layers are potentially eroded successively.The erosion of each layer consumes part of dt so that less time remains to erode the underlying layer.This time reduction is taken into account by multiplying the time step dt by (1 − volume layer w dx dt ) between the layers.If the surface layer is made of sediment, their erosion takes time and limits the erosion of the bedrock below.This process is known as the "cover" effect for river (e.g.Whipple and Tucker, 2002;Lague, 2010) and thus is taken into account in Cidre.

Bank or lateral erosion
Flowing water can erode lateral cells (Fig. 1a).These cells are those which are topographically above the considered cell and in a lateral direction perpendicular to each downstream direction.The lateral sediment volumetric flux Q sl [L 3 /T ] is simply a fraction of the sediment flux Q s [L 3 /T ] in the considered direction (e.g.Murray and Paola, 1997;Nicholas and Quine, 2007): where α is a bank erodibility coefficient.α is specified for sediment and is implicitly determined for bedrock layers proportionally to their "fluvial" erodibility such that α sediment /α bedrock = K sediment /K bedrock (K from Eq. 7).If sediment covers the bedrock of a lateral cell, α is weighted by its respective thickness above the target cell.Equation ( 10) is debatable and alternative formulations have been proposed (e.g.Davy and Lague, 2009).The main point is that bank erosion is introduced here to allow clasts from the sides of rivers to be incorporated into the river by lateral transport, as we expect this to control clast spreading.

Sediment spreading
The sediment leaving a cell is spread in the same way as water (Fig. 1a).This rule is justified by the probability of moving proportionally to the slope in the case of diffusive hillslope processes.As water is spread according to this rule, sediment in suspension is spread in the same way.Flume experiments may also support this rule for bedload (Seizilles et al., 2014).
Compared to previous published versions of Cidre (Carretier et al., 2009;Pepin et al., 2010;Carretier et al., 2014), significant modifications have been made: (1) the introduction of the transport length scale L, (2) the new hillslope erosion law and (3) the calculation of erosion using the steepest slope and then the spreading of sediment toward all downstream slopes.In the previous version, a transport capacity was calculated in each downstream direction and then summed to determine the cell erosion.The two approaches are not equivalent, in particular if the river erosion law depends non-linearly on water discharge.The advantage of the new version is that it is more stable and rapid while still allowing flow divergence and spreading.One motivation for these changes was to explicitly separate erosion, deposition and transfer rates on each cell, so that the treatment of clast movement will be straightforward.

Adding and following clasts
Once q s , and D have been calculated over the entire grid during a time step, the clasts are treated one after the other (the order does not matter).The philosophy of our approach is to use the local and D rates calculated by Cidre to determine a probability for a clast to move, cross or deposit onto a cell (Fig. 1b).
A clast has a specified radius R, with no particular limitation, between a small mineral to a large cobble.Its initial position and depth z clast (of its base) are specified.When it is set initially in the bedrock, it corresponds to a piece of rock which will have this size once detached.For example, a clast can represent a quartz mineral individualised by weathering, or a cobble the size of which is determined by a particular fracture spacing in the parent rock.
A clast is detached (eroded) if its depth is shallower than the erosion calculated over the time step on that cell (if z clast ≤ dt).If it is detached, its probability to go in one of the downstream directions is simply the ratio of the local slope and the sum of the downstream slopes, consistent with the multiple-flow algorithm (Fig. 1b).
For a moving clast entering a cell, the probability that it will be deposited is simply the ratio between the volumetric deposition flux and the volumetric incoming flux D w dx / Q s .Therefore, its probability to cross the cell is 1 − D w dx / Q s (Fig. 1b).In practice, D w dx / Q s is compared to a random value between 0 and 1.If this random value is lower than or equal to D w dx / Q s , the clast deposits on the cell.If not, it continues its travel and is transferred towards one of the downstream cells and so on.This cell is chosen following a similar procedure so that the probability to go in one direction is the ratio of the local slope and the sum of the downstream slopes.
A clast may be detached but not leave the cell.This may occur if the clast was at depth.Removing material above the clasts takes part of the time step, so that the remaining time may prevent the clast from leaving the cell.Furthermore, a big clast should have a lower probability to leave the cell than a small one.In order to take these realities into account, a probability to leave the cell is determined by where δ = 1 if the direction of movement is parallel to rows or columns, and δ = 1/ √ 2 along diagonals (a longer distance decreases the probability to leave the cell).The value of probability law Eq.( 11) is set to 1 if it exceeds unity, which may occur if the clast is at the surface and the erosion is larger than the clast diameter.The first term accounts for the inverse relationship between entrainment probability and clast size (e.g.Malmaeus and Hassan, 2002).The numerator dt renders the travel distance independent of the time step.The second term decreases the probability to move if the clasts were at depth.The coefficient 1.25 is justified a posteriori to adjust the mean travel rate of a population of clasts to theoretical predictions (see below).In practice, the decision for a clast to leave or not a cell is determined as follows: a "leave cell" number (between 0 and 1) is calculated using Eq. ( 11).It is compared to a random value between 0 and 1.If the "leave cell" number is larger or equal than this random value, the clast leaves the cell.If not, it remains at the cell surface.
Note that the probability of deposition could also depend on clast size.This is not implemented here but we will return to this point in the discussion.
The erosion-deposition-transfer fluxes calculated from the deterministic rules of Cidre can be viewed as the mean values for the distribution of the clast radius R. Consequently, a 20 cm boulder may travel slowly across a cell because it has a low movement probability (cf.Eq. 11, where the probability to leave the cell is proportional to 1/2R), while the 5 cm sedimentary layer above the bedrock is removed during a time step.This situation can occur because the transport flux calculated by Cidre represents a time average over an implicit distribution of the grain size (this is a volume of transported sediment without specifying its size), whereas the clasts are individuals that move stochastically.This distinction is consistent with what is observed in the field.For example, we can imagine a situation where boulders go down an eroding bare bedrock hillslope.A mean long-term denudation rate for the hillslope can be determined, consistent with the diffusion model, for example, while the boulders go down stochastically.

Initial setup
The initial setup of the 3-D landscape model is an initial elevation grid, a grid or uniform value for the uplift-subsidence rates, a grid or uniform value for the precipitation rates, a geological model, and the boundary conditions.The geological model consists of different erodible layers.Their thickness is specified on each pixel.Values for the "erodibility" parameters (Eqs.5, 6 and 7) are attributed to these different materials.Sediment erodibility (resulting from either the deposition of physically detached material or in situ regolith production by bedrock weathering) can be set differently from that of the bedrock.Nevertheless, in the experiments presented in this manuscript, there is only one bedrock type.
In addition, the clasts are initially listed in an input file specifying their initial location in the grid, as well as their depth, radius and mineralogy.There is no limitation to their number except for the one imposed by the computational times.Their initial localisation can be chosen according to a specific goal.For example, they can be grouped within one pixel to follow their transport from a specified source, or spread randomly in the model grid to study the mean transport rate of the sediment particles at the catchment scale.The distribution of the initial clast size can be freely chosen to trace one particular grain size, or a distribution of the grain sizes.

Linear diffusion
One way to validate the above model is to demonstrate that the displacement of a clast population follows predictions in simple cases.We begin with the case of linear diffusion.In Cidre, linear diffusion is obtained by using Eqs.( 5) and ( 6) only, where = κS and L = dx.The resulting linear and local sediment flux per unit width is q s = κ dx S, and the product κ dx defines the diffusivity κ [L 2 /T ].
We consider an inclined plane of slope S. Clasts initially set at the top surface of the plan will travel at a mean velocity ẋ given by the ratio of the sediment flux per unit width q s and a depth corresponding to one clast diameter: Thus, the mean travel distance of the clast population x = κ S 2R t.Furthermore, within the framework of Fickean diffusion theory, the standard divergence of the distances between particles is predicted to vary as σ = √ 2κ t (Einstein, 1905).
These predictions are tested on a bedrock plan of S = 0.1.The same κ is used for the bedrock and sediment (product of erosion) to simplify.Sediment and clasts can only go out on the lowest boundary condition, the elevation of which is fixed at 0 m.Different values of κ, dx and dt are tested in order to evaluate the effect of cell size and time step on the results.The diffusivity is κ = 0.01 m 2 yr −1 in all cases.Initially, 1000 clasts are set at the top of the plan and then followed during their descent.Their radius R = 0.005 m.
Figure 3 illustrates the progressive spreading of the clast population during the plane's erosion.x fits the predicted evolution well and is not sensitive to the cell size or the time step for a given κ .This non-dependence comes from the chosen form of the probability law Eq.11.In this law, the factor 1.25 was added to better adjust the calculated x to κ S 2R t.This coefficient works for other κ, dx, and dt values.It may correspond to the numerical diffusion or to the C++ random function used to move the clasts.
As predicted by the diffusion theory, the scattering of the clasts σ ∝ √ t.It does not depend on dt for a large range of time steps (5 to 50 years in Fig. 3), but σ is significantly lower for larger dt values (250 years).For this time step, the eroded thickness at each step is larger than the diameter of the clasts, so that the "leave cell" probability Eq. 11 is increased.Hence, at each time step, there is less chance that a clast will stay while the other ones move.The clast population moves more homogeneously and, consequently, σ is lower than with a smaller time step.Note that the effect of this time step on σ is consistent with natural processes where erosion occurs by discrete events.A thick erosion event necessarily entrains clasts that have a smaller diameter.For large time steps, the mean travel rate of the clast population remains equal to κ S 2R ; the effect of the time step only affects the scattering.
In all cases, σ strongly overestimates the predicted √ 2κ dx t.However, σ tends towards this predicted value when dx decreases or dt increases.When a clast moves, it travels a distance which is at least dx.Thus, the distance between immobile and moving clasts necessarily depends on dx.This is a limitation of this modelling approach.It is possible to adjust σ to the correct value, independent of dx, by replacing ( dt 2R ) by ( κ dt dx ) in Eq. ( 11).However, the mean travel distance becomes strongly dependent on dx in this case.
In this version, we prefer probability law Eq.( 11), which allows the mean travel distance to be respected, for three reasons.x is directly linked to the sediment mass balance of a portion of the landscape, not σ .In practice, in a realistic topography, convergence zones and canals naturally limit the scattering of the clasts.Finally, many potential applications will be more sensitive to x, which determines the mean residence time, rather than σ .

"Non-linear" diffusion with a slope threshold
Here we illustrate the scattering of the clasts in the case of "non-linear" diffusion.We use Eqs.( 5) and ( 6) with S c = www.earth-surf-dynam.net/4/237/2016/Earth Surf.Dynam., 4, 237-251, 2016 0.105, a value just above the plan slope S = 0.1.This S c value is much lower than published values ∼ 0.6 (e.g.Roering et al., 2001), but our goal here is to compare this with the previous linear diffusion case.We also use coarser clasts with R = 0.05 m to obtain similar mean transport rates as in the previous case.The other parameter values are the same as in the linear diffusion case.
Figure 4 shows the pattern of clasts at different times.Compared to the linear diffusion case, once detached, clasts move a longer distance because L dx when S ∼ S c .Consequently, the clast population presents a much more elongated distribution along the hillslope.This pattern seems consistent with colluvial sediment covering steep hillslopes in the real world.The mean travel distance x does not depend on dx nor dt, while σ does depend on dx.Note that σ varies as √ t but with much larger values representing the downstream spreading of the clasts.The different patterns obtained in the linear and non-linear diffusion cases illustrate the progressive transition from diffusion for S S c to ballistic movements for S ∼ S c .

River transport
River transport usually implies the formation of incisions, local depositions and lateral movements of sediment by bank erosion.The movement of clasts associated with this dynamic is an active research field (e.g.Seizilles et al., 2014).In order to verify that x is consistent with the description of the river process using erosion and deposition rates, we simplify the problem by using a planar and stable portion of a river bed.This bed corresponds to the same bedrock plan with the slope of 0.1 used above.Detachment is determined by = K q S. The same K is used for both bedrock and sediment.The flow width w is set at cell size dx.Thus, corresponds to the "stream power law" (e.g.Whipple and Tucker, 1999).The transport length L = ξ q is close to dx, by setting ξ to a low value 10 −6 (L cannot be smaller than dx).Consequently, the sediment flux per unit width q s is close to dx ∝ q.The linear relationship between q s and q prevents incisions from forming (Smith and Bretherton, 1972), which is what we want here.A precipitation rate of P = 1 m yr −1 is imposed on the grid.Thus, the water discharge per unit width q = P x grows linearly downstream, as does the probability of clast entrainment.
The mean travel velocity of a clast ẋ can be written as , where x o is the initial distance.For a clast population, x must follow this solution.
We use larger clasts with a radius of 0.05 m in order to have a transport rate of the same order of magnitude as in the diffusion case.Figure 5 shows that clasts are spread from the injection point, forming a downstream plume.This pattern results from the downstream increase in the entrainment probabilities.x fits the above solution acceptably well and does not depend on dx or dt.The scattering σ increases more than linearly because the water discharge increases downstream.
Note that, despite similar plan view patterns, σ evolves very differently in the non-linear diffusion and river cases.This difference illustrates the interest of coupling clasts and flux in order to identify erosion-transport processes.

A more general case
We present here an illustration of clasts moving in a mountain-foreland context.Our goal is not to precisely analyse the clast dynamics in that case but to qualitatively describe a possible situation in a real-case scenario.
The mountain-foreland system consists of an uplifting block (the mountain) and a stable domain where sediment eroded from the block is deposited or in transit.The grid is 200 × 200 cells, the cell size is 500 m, and the mountain is 37.5 km wide, which could correspond to the Pyrenees or certain Andean foothills in Argentina, for example.A constant uplift rate is fixed at 0.5 mm yr −1 .The precipitation rate is homogeneous and constant at 1 m yr −1 .Sediment, water and clasts can escape the model grid on the southern boundary, where the elevation is fixed at 0 m.The east and west boundaries are periodic, so that sediment, water and clasts can go out from one side and re-enter on the other side.No sediment, water or clasts can leave the northern side.Erosion parameters are chosen so that the maximum elevation is 3000 m when the topography reaches a rough steady state at around 3 Myr (see Fig. 6 caption).Additional figures in the Supplement provide more information concerning the slope-area distribution of the uplifting block at steady state (Fig. S1 and S2) and the scaling between valley width and water discharge (Fig. S3).
The experiment begins from a flat topography with small random elevations (< 1 m) (Fig. 6).The sediment eroded from the uplifting block is propagated through the foreland, forming alluvial fans and rivers.An upstream drainage develops in the mountain, and a divergent river network progrades over the alluvial apron.Initially, 2000 clasts were seeded at two different places in the mountain at depths of 500 to 600 m (green clasts) and 400 to 500 m (red clasts), as if they were belonging to two intrusive bodies.When erosion reaches this depth, clasts begin to exhume, move and deposit in the foreland.They stop because they are progressively buried under sediment.The red clasts exhume first and are spread within an alluvial fan.Deposition of the green clasts ends 0.5 Myr after the deposition of the red clasts.The complete deposition of both populations lasts 1 Myr.Both populations are mainly included in a sediment layer that is thinner than 100 m.They are mixed at similar depths at the transition zone between the two alluvial fans.This situation illustrates a lateral diachronism between close alluvial aprons which may have been erroneously correlated in the field.
If these grains were detrital zircons of different ages, the analysis of their age distribution at different places would allow the mixed zone to be mapped and, thus, the dynamics of the lateral alluvial fan to be reconstructed.

Perspectives
The modelling approach described in this paper and its developments may have different applications, which we propose in the following.

River-erosion law
River dynamics involve processes acting on a large range of time periods, from hours in the case of catastrophic flooding to thousands of years to transport huge volumes of glacier sediment, for example.Determining simplified laws to predict this complexity remains a challenge (Murray, 2007).Recently, Kasprak et al. (2015) showed a tight relationship between the inter-bar distance and travel distance of marked grains in experiments on braided rivers.At the same time, Davy and Lague (2009) evidenced a relationship between this inter-bar distance and L, a fundamental transport-length parameter used in the Eros landscape evolution model and which has been also introduced in Cidre.The link between both observations can be made using the modelling approach proposed in this paper.By using a try-and-guess approach, investigating different erosion and deposition laws of different complexities, it may be possible to identify which laws predict the distribution of the clasts and the river bed forms observed in experiments (Kasprak et al., 2015) or rivers (Bradley et al., 2010;Bradley and Tucker, 2012).As far as we know, the study that was carried out by introducing moving clasts in Eros is the only study of this kind (Davy and Lague, 2012).

Grain sorting, coarsening
Fining of the grain size by attrition in very steep catchments and by selective deposition in most of the other catchments is a well-observed phenomenon.(Paola et al., 1992;Paola and Seal, 1995;Attal and Lave, 2006;Fedele and Paola, 2007;Miller et al., 2014).Selective deposition could be simply included by allowing the deposition probability to depend on clast size.Physical attrition could also be included by decreasing the clast size according to the travel distance.Local surface coarsening could be also studied, as clast entrainment depends on their size.This model may represent an alterna- tive to models requiring a depth discretisation of sediment layers containing different grain sizes (Cohen et al., 2009).

Grain-size-dependent landscape evolution
In the real world, material erosion and transport depend on the clast size and the grain size evolves spatially.The feedback of grain size on the landscape dynamics has been little explored (Sklar et al., 2008).In Cidre, entrainment and deposition parameters may vary according to the mean radius of the grains on each cell.The number of model parameters would be kept small (K, L specified initially for the median size, for example), while accounting for complex feedbacks.This could be an alternative to other approaches which incorporate explicitly different classes of grain sizes to simulate fluvial dynamics (e.g.Coulthard et al., 2002).

Provenance problems and alluvial or placer-type deposits
Provenance studies on detrital grains help constrain the chronology of the exhumation of the sediment source (e.g.Roddaz et al., 2005;Rodríguez et al., 2012).Grains of a particular lithology can be used to identify the provenance of certain material of economic interest, like distant grav-els eroded from a hidden porphyry copper deposit (Sillitoe, 2010).These provenance studies come up against uncertainties regarding the landscape evolution.The coupled modelling of landscape evolution and clast movement may help constrain this link in specific cases.Placer-type deposits are secondary ores that can contain free particles with very fine gold and other native metals (e.g.platinum-group elements).The occurrence of gold grains in supergene environments, such as soil, sediments and placers in rivers, is controlled by physical (as well as bio-geochemical) processes of redistribution from a distant gold-quartz vein (Reich and Vasconcelos, 2015).Gold grains between 0.1 and 4 mm in diameter constitute most of the alluvial and eluvial gold, and represent economically important ores, such as the Witwatersrand palaeoplacer gold deposits in South Africa (Mossman et al., 1999;Fairbrother et al., 2013).Measurements of geometric grain parameters such as roundness, degree of bending, and grain flatness have been used to constrain the source type of the gold and transportation distances (e.g.Townley et al., 2003;Mudaliar et al., 2007).Therefore, an approach that involves the modelling of gold grain transport during the landscape evolution might provide an independent test to check the validity of transport distance estimations based on grain morphology.

Cosmogenic nuclides
The analysis and modelling of cosmogenic nuclide concentration in individual clasts give quantitative information about the erosion-transport processes at the landscape scale (Repka et al., 1997;Codilean et al., 2008;Gayer et al., 2008;McPhillips et al., 2014;Carretier et al., 2009;Carretier and Regard, 2011;Vassallo et al., 2011).It would be straightforward to compute the cosmogenic nuclide concentration evolution in clasts in Cidre, using, for example, the model of Carretier and Regard (2011).This may be particularly useful to link cosmogenic nuclides from buried sediment in basins to quantify the erosion-sediment history of mountainforeland systems.

Tracing weathering material
In the 1990s, Raymo and Ruddiman (1992) proposed that the uplift of Tibet has driven the post-Eocene global atmospheric cooling by accelerating the chemical weathering of silicate rocks and the associated consumption of atmospheric CO 2 .However, recent studies at the foot of the Himalayas and Andes suggest that the weathering flux may be larger in the foreland basin than in the mountain range itself (Lupker et al., 2012).A key issue is to understand the distribution of the residence times in the different reservoirs (Mudd and Yoo, 2010;Anderson, 2015;Hoffmann, 2015).These studies call for the development of 3-D models able to trace the weathered material from source to sink (Anderson et al., 2012;Vanwalleghem et al., 2013;Minasny et al., 2015).The clast-flux modelling approach may be developed within that scope, by allowing clasts to weather.This will be the aim of a future paper.

Conclusions
The algorithm predicts a consistent clast velocity and surface erosion rate in simple cases.The mean travel distance of the clasts does not depend on the model cell size or time step.The scattering of the clasts depends on the cell size and is overestimated.Nevertheless, decreasing the cell size decreases the overestimation.This model has numerous potential applications allowing field data on distinct grains to be linked to a large-scale landscape evolution.The differences between the simple river and hillslope cases illustrated here (e.g.Figures 4 and 5) show that the movement of tracers in the field or in experiments may provide a way to establish local erosion-transport laws at spatial and temporal scales adapted to landscape evolution models.

Figure 1 .
Figure 1.Illustration of erosion-deposition processes in Cidre.(a) Multiple flow and the different calculated fluxes: q s [L 2 /T ] incoming sediment flux per unit width.T i [L/T ] transferred sediment rate in each downward direction.[L/T ] detachment rate of material (sediment or bedrock) from the cell calculated in the steepest direction.D [L/T ] deposited sediment rate on the cell.Qs L [L 3 /T] lateral volumetric sediment flux deposited on the cell.S i slope in each downward direction.w = flow width.dx = cell width.The incoming water flow is propagated towards lower cells proportionally to the local slope in each direction (not illustrated).(b) Probabilities that clasts will sediment, transfer or detach using fluxes calculated by Cidre.

Figure 3 .
Figure 3. Model test in the linear diffusion case ( = κ S and L = dx).The initial plan slope is S = 0.1.1000 clasts are initially set at the surface and top of this plan.The diffusivity κ = κ dx = 0.01 m 2 yr −1 in all the cases.(a) Example of clasts spreading in the case κ = 0.0001, dx = 100 m and dt = 50 years.(b) Evolution of the mean travel distances of the clasts for the different dx, dt and κ values.Note that the mean transport distance fits well with the prediction and does not significantly depend on dt or dx.(c) Evolution of the standard deviation of the travel distances of the clasts σ .The black dashed curve is a √ t fit to the reference (greenred) model curves, showing that σ ∝ √ t, consistent with the diffusion theory.Nevertheless, this fit overestimates the predicted curve (σ = √2κ dx t) by more than 1 order of magnitude, which should correspond to a value of 0.036 km at t = 50 kyr for the green and red curves.

Figure 4 .
Figure 4. Model test in the non-linear diffusion case ( = κ S and L = dx 1−(S/S c ) 2 ).The initial plan slope is S = 0.1 and S c is set close at 0.105.The mean sediment flux per unit width q s = dx is the same in all cases (same κ = κ dx).One thousand clasts are initially set at the surface and top of this plan.(a) Example of clasts spreading in the case κ = 0.0001, dx = 100 m and dt = 50 years.(b) Evolution of the mean travel distances of the clasts for the different dx, dt and κ values.Note that the mean transport distance fits well with the prediction and does not depend significantly on dt or dx.(c) Evolution of the standard deviation of the travel distances of the clasts σ .σ ∝ √ t, consistent with the diffusion theory.

Figure 5 .
Figure 5. Model test in the river case ( = K q S and L = dx).The initial plan slope is S = 0.1 and the precipitation rate is 1 m yr −1 .Initially, 1000 clasts are set at the surface and top of this plan.The mean sediment flux per unit width q s = dx is the same in all cases by keeping K dx = 310 −5 constant.(a) Example of clasts spreading in the case K = 310 −7 m −1 , dx = 100 m and dt = 50 m.(b) Evolution of the mean travel distances of the clasts for the different dx, dt and K.The blue colouring indicates an increase in the downstream water discharge.Note that the mean transport distance fits well with the prediction and does not significantly depend on dt or dx.(c) Evolution of the standard deviation for the travel distances of the clasts σ .Note that σ increases exponentially downstream, differently from the diffusion cases.

Figure 6 .
Figure 6.Example of clasts exhumation from two locations (intrusive body) located initially at a depth of 0.5 km (green clasts) and 0.4 km (red clasts).The final maximum elevation for the mountain is 3000 m.The domain size is 100 km × 100 km.There is constant rainfall over the whole grid (1 m yr −1 ).The northern block is continuously uplifting at 0.5 mm yr −1 .The elevation of the clasts is increased so that they are visible.Note that the green grains stop depositing 1 Myr after the red ones.Nevertheless, both populations are mixed between the two alluvial fans at similar depth, showing that 1 Myr of difference may separate two adjacent sediment layers that are apparently synchronous.K sediment = 0.6 × 10 −3 m yr −1 , K bedrock = 0.3 × 10 −3 m yr −1 , m = 0.5, n = 1, ξ = 0.1 year m −1 , α sediment = 0.01, κ sediment = κ bedrock = 10 −4 m yr −1 , Sc sediment = Sc bedrock = 0.8, dx = 500 m, dt = 100 year.