Modelling detrital cosmogenic nuclide concentrations during landscape evolution in Cidre v2.0

The measurement of cosmogenic nuclide (CN) concentrations in riverine sediment has provided breakthroughs in our understanding of landscape evolution. Yet, linking this detrital CN signal and relief evolution is based on hypotheses that are not easy to verify in the field. Models can be used to explore the statistics of CN concentrations in sediment grains. In this work, we present a coupling between the landscape evolution model Cidre and a model of the CN concentration in distinct grains. These grains are exhumed and detached from the bedrock and then transported in the sediment to the catchment outlet with temporary burials and travel according to the erosion–deposition rates calculated spatially in Cidre. The concentrations of various CNs can be tracked in these grains. Because the CN concentrations are calculated in a limited number of grains, they provide an approximation of the whole CN flux. Therefore, this approach is limited by the number of grains that can be handled in a reasonable computing time. Conversely, it becomes possible to record part of the variability in the erosion–deposition processes by tracking the CN concentrations in distinct grains using a Lagrangian approach. We illustrate the robustness and limitations of this approach by deriving the catchment-average erosion rates from the mean ¹⁰Be concentration of grains leaving a synthetic catchment and comparing them with the erosion rates calculated from sediment flux, for different uplift scenarios. We show that the catchment-average erosion rates are approximated to within 5 % uncertainty in most of the cases with a limited number of grains.

1 Introduction

The concentration of cosmogenic nuclides (CN) produced in situ varies according to the depth of the minerals in which they are produced, their altitude, the stable or radioactive nature of the nuclide in question and the magnetic field (Gosse and Phillips, 2001). When a mineral in a grain is exhumed and then transported in rivers, the concentration of CN evolves by integrating stochastic variations of the residence times at different depths and altitudes over time. This integrative characteristic has been used to develop numerous approaches to quantify erosion–deposition processes averaged over catchment areas and millennia (Schaefer et al., 2022). The most widespread approach is one that quantifies the average erosion rate of a catchment from its average CN concentration in a sand sample taken at its outlet (Brown et al., 1995). Since the pioneering work of Repka et al. (1997), other studies have explored the possibility of using the distribution of CN concentrations in distinct grains and possibly grains of different sizes to quantify erosion–deposition processes on hillslopes, in rivers and on alluvial deposits over periods of several thousand to millions of years (Braucher et al., 1998; Dunai et al., 2005; Gayer et al., 2008; Codilean et al., 2008; Carretier et al., 2019). However, it is still difficult to link detrital CN concentration data to specific processes, whether on hillslopes or on a larger scale in river systems (Yanites et al., 2009). This difficulty arises from the stochasticity of grain transport, which can be temporarily stored and then recycled on different scales ranging from a flood event that erodes the banks of an alluvial river to millions of years in the case of exhumation of ancient strata of a foreland basin.

To address this complexity, several models have been developed at the grain and catchment scales (Repka et al., 1997; Niemi et al., 2005; Codilean et al., 2008; Carretier et al., 2009, 2019; Ben-Israel et al., 2022), but without taking the evolution of the relief, and thus of the CN production rate, into account. At the same time, landscape evolution models (LEMs) have made great progress by integrating an increasing number of processes; however, their resolution over long time spans limits taking into account stochastic phenomena that require the intermittent storage of sediment grains on hillslopes (e.g. landslides) and along fluvial systems (e.g. sediment bars and terraces). Carretier et al. (2016) proposed a coupling between the landscape evolution model in Cidre and the transport of tracer grains which are transported stochastically according to simple probability laws depending on the local erosion and deposition rates calculated on each cell of the LEM. This allowed those authors to highlight that some grains may have been stored for a long time before being recycled and evacuated by the rivers (Carretier et al., 2020). In this paper, we present a development in Cidre to calculate the concentration of several CNs within the grains. To our knowledge, the only published model that can jointly model the evolution of the relief and the evolution of the average CN concentration in the sediments is Badlands (Petit et al., 2023). In Badlands, this average concentration is calculated in proportion to the sediment fluxes from the different upstream sources on each grid cell. In Cidre, we adopt a different (i.e. Lagrangian) approach, which enables us to track the full distribution of CN concentrations in a population of grains.

We first present the Cidre model, including grain transport, and then describe the equations used to calculate the CN concentration over time in each grain. We show the accuracy and robustness of the algorithm by comparing the average catchment erosion rate derived from the CN concentrations of outgoing grains with the true rate calculated in several simulations. In the following we call the “true” catchment-average erosion rate the ratio of the sediment “outlux” over the catchment area calculated in Cidre. In particular, we show that this true rate is approximated to within 5 % uncertainty in most of the cases with a limited number of grains.

2 The landscape evolution model Cidre

2.1 Mass balance equation

Cidre is a c++ code that solves the following mass balance on rectangular cells expressed here for the mean elevation z of a cell through time t:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial z}{\partial t}}=-{\mathit{ϵ}}_{\mathrm{r}}-{\mathit{ϵ}}_{\mathrm{h}}+D_{\mathrm{r}}+D_{\mathrm{h}}+U.\end{array}$$

This mass balance refers to the erosion (“ϵ”)–deposition (“D”) model (Davy and Lague, 2009) where the subscript “r” (“river”) refers to erosion driven by flowing water and “h” (“hillslope”) to erosion due only to the topographic gradient or slope S. U is a vertical uplift or subsidence rate. Then, we define a constitutive law for each of these components (Carretier et al., 2016):

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{r}}& {\displaystyle }=K{q}^m{S}^n\text{ for river processes},\text{(3)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{h}}& {\displaystyle }=\mathit{\kappa }S\text{ for hillslope processes},\end{array}$$

where K [L^1−2 m T^m−1] and κ [L T⁻¹] are erodibility parameters; m and n are lithology-dependent (different for bedrock or sediment) erosion parameters; S is the slope (absolute value of the elevation gradient) towards the downstream cell in the steepest direction; and q [L³ T⁻¹] is the water discharge per stream unit width, corresponding to the accumulation of the net precipitation rate (specified in the input or varying dynamically with elevation; Zavala et al., 2020) from the highest to lowest cell according to a multiple flow algorithm, spreading the water discharge of the cell towards all the lower cells among the eight neighbouring cells proportionally to their slope. There is no erosion in a pit cell, only deposition. The deposition rate is

$$\begin{array}{}\text{(4)}& {\displaystyle }{D}_{\mathrm{r}}& {\displaystyle }={\displaystyle \frac{q_{\mathrm{sr}}}{\mathit{\zeta }q}}\text{ for river processes},\text{(5)}& {\displaystyle }{D}_{\mathrm{h}}& {\displaystyle }={\displaystyle \frac{q_{\mathrm{sh}}}{\frac{\mathrm{d}x}{\mathrm{1}-(S/S_{\mathrm{c}}{)}^{\mathrm{2}}}}}\text{ for hillslope processes},\end{array}$$

where q_sr and q_sh are the incoming river and hillslope sediment fluxes (total $q_{\mathrm{s}}=q_{\mathrm{sr}}+q_{\mathrm{sh}}$) per unit width [L² T⁻¹], ζ is a river transport length parameter [T L⁻¹] and S_c is a slope threshold. These fluxes are the sum of the sediment fluxes leaving upstream neighbouring cells while the deposition rates on a cell are a fraction of the incoming sediment. Note that the hillslope equations derive from the non-linear diffusion model (Roering et al., 1999) but are written as an erosion–deposition model. Both formulations lead to the similar topographic evolution, but the model used in Cidre is numerically more stable and more adapted to the coupling with grains transport (Carretier et al., 2016).

Equation (1) is solved using the forward finite volume method (z represents the mean elevation of a cell). At each iteration, cells of the model grid are ranked in a decreasing elevation order and then treated successively in that order to ensure that all incoming water and sediment fluxes are known when treating a given cell. For each cell, incoming water fluxes are summed with the local volume of precipitation falling on a cell over the time step. The resulting water flux is then spread among all the downstream cells in proportion to their local topographic gradient (multiple flow), or slope. Then the eroded flux is first calculated on that cell from Eqs. (2) and (3) using the steepest-descent slope among the downstream neighbouring cells. If the erosion potential (the erosion rate multiplied by the time step) is larger than the sediment volume present on that cell, the bedrock is eroded in proportion to the remaining time step ($\mathrm{1}-\text{sediment volume}/\text{erosion potential}$). Then, the deposition flux is calculated using Eqs. (4) and (5). Erosion–deposition is first calculated for the hillslope processes and then for rivers. Once the total eroded and deposited volumes are known, the elevation is modified based on the balance between these two volumes. Sediments that leave the cell are distributed to downstream cells in proportion to the local slope of neighbouring cells. Then, the next cell in the list is treated and so on. When all the cells have been treated, a new iteration begins and the time is incremented by dt.

Other processes, such as lateral erosion, weathering and regolith development, orographic precipitation or the dynamic filling of depressions, are implemented but not described in this paper because the algorithm to calculate the CN concentrations in the grains does not vary according to these processes (see Carretier et al., 2016, 2018; Zavala et al., 2020).

2.2 The grains in Cidre

Grains are clasts of any kind, for example, a single mineral or a pebble composed of multiple minerals, ranging in size from millimetres to decimetres (Carretier et al., 2016). They are localized by the index corresponding to the cell number where they are located (their precise 2D coordinates within a cell are not known) and the depth of their centre beneath the earth's surface. At the beginning of a simulation, their number, radius R, location on the grid and depth z^′ are specified. For example, they can be set randomly on the grid and at depth if the grains are quartz grains and the proportion of quartz is constant in the underlying rock. Alternatively, grains can be set with a higher proportion in some cells or at some depths for which the rock has a higher quartz content.

Grains are moved once the erosion and deposition rates are known on the model grid at the end of a time step. They are passive tracers, they do not influence the erosion and sediment calculated by Cidre and they are all independent from each other. Grains are eroded, transported or deposited according to probabilities depending on the erosion and deposition rates. For a grain on a cell, it is detached if the eroded layer on that time step is thicker than or equal to the depth of the grain's bottom $z_{\mathrm{b}}^{\prime }$. If the grain is not detached, its depth is updated according to the dz calculated on that cell and for this time step. It is possible for a grain to be detached but to not leave the cell to account for the time needed to travel to a larger cell, and for the slower transport of a big grain. Carretier et al. (2016) found that to account for these effects and to predict the correct mean sediment flux, the probability of leaving a cell should be $\mathrm{1.25}\left(\frac{\mathit{ϵ}\cdot \mathrm{d}t}{\mathrm{2}R}\right)(\mathrm{1}-\frac{z_{\mathrm{b}}^{\prime }}{\mathit{ϵ}\cdot \mathrm{d}t})\mathit{\delta }$, where δ=1 if the direction of movement is parallel to rows or columns and $\mathit{\delta }=\mathrm{1}/\sqrt{\mathrm{2}}$ along diagonals. (A longer distance decreases the probability of leaving the cell.) The computation is made successively on a cell for hillslopes (ϵ=ϵ_h) and for river processes (ϵ=ϵ_r). The value of this probability is set to 1 if it exceeds unity, which may occur if the grain is at the surface and the erosion is larger than the grain diameter. If the grain does not leave the cell, its depth is chosen randomly within the layer deposited on that cell during the time step.

If a grain leaves the cell, it goes to one of the downstream neighbouring cells with a probability set as the ratio of the local slope of the considered downstream cell and the sum of the downstream slopes. Once the grain has entered into one of these downstream cells, the probability that it will be deposited is the ratio of the sediment deposition volume in that cell to the sum of the incoming sediment volumes. If not deposited, the grain continues its travel and is exported towards one of the downstream cells, and so on until it is deposited or leaves the model grid. Then a new grain is treated and so on. When all the grains have been treated, a new time step is initiated to calculate the local erosion and deposition on the model grid.

Figure 1 The vertical coordinates of two grains, one staying at depth during a time step and the other moving to a downstream cell where it is deposited. Dashed lines correspond to t−dt and solid lines to t. The upstream cell has a net negative topographic change (erosion) over the time step “dt” (ϵ>D). The downstream cell has a net positive change (sedimentation). The eroded (ϵ⋅dt) and deposited (D ⋅dt) thicknesses are calculated for hillslope and then river processes (see Eqs. 2 and 4). In both cases, a grain is potentially moving if $\mathit{ϵ}\cdot \mathrm{d}t>z_{\mathrm{b}}^{\prime }$. The depth z^′ of a grain's centre is updated according to the net erosion in the upstream cell or as random depth within the deposited layer of the downstream cell.

3 Calculating the concentration for different CNs

3.1 CN evolution

The CN concentration C at the centre of the grain varies based on

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}t}}=-\mathit{\lambda }C+P,\end{array}$$

where λ is the radioactive decay rate [T⁻¹] and P is the CN production rate [𝒩/M/T] at any time t. The solution is calculated as (Carretier et al., 2009)

$$\begin{array}{}\text{(7)}& {\displaystyle }C\left(t\right)=e^{-\mathit{\lambda }t}(C_o+I),\end{array}$$

where C_o is the initial CN concentration (see next section) and

$$\begin{array}{}\text{(8)}& {\displaystyle }I=\sum _n{e}^{\mathit{\lambda }t}P\mathrm{d}t,\end{array}$$

where n is the number of time steps dt since the beginning of the landscape evolution simulation and $t=n\cdot \mathrm{d}t$.

3.2 Production rate

The production rate depends on the nuclide considered. Here, we use a formulation based on three contributions by spallation (subscript “sp”), reactions induced by slow muon capture (subscript “sm”) and by interactions with fast muons (subscript “fm”) (Braucher et al., 2011):

$$\begin{array}{}\text{(9)}& {\displaystyle }P=(P_{\mathrm{sp}}{e}^{-\mathit{\rho }{z}^{\prime }/{\mathrm{\Lambda }}_{\mathrm{sp}}}+P_{\mathrm{sm}}{e}^{-\mathit{\rho }{z}^{\prime }/{\mathrm{\Lambda }}_{\mathrm{sm}}}+P_{\mathrm{fm}}{e}^{-\mathit{\rho }{z}^{\prime }/{\mathrm{\Lambda }}_{\mathrm{fm}}}),\end{array}$$

where z^′ is the depth of the grain's centre; P_sp, P_sm and P_fm are the production rates of a given CN at the earth's surface by spallation, slow muon capture and fast muon interactions, respectively; ρ is the grain density; and Λ_sp, Λ_sm and Λ_fm [M L²] are the respective attenuation factors with depth:

$$\begin{array}{}\text{(10)}& {\displaystyle }{P}_{\mathrm{sp}}& {\displaystyle }=P_{\mathrm{SLHL}}{f}_{\mathrm{sp}}{S}_{\mathrm{sp}},\text{(11)}& {\displaystyle }{P}_{\mathrm{sm}}& {\displaystyle }=P_{\mathrm{SLHL}}{f}_{\mathrm{sm}}{S}_{\mathrm{sm}},\text{(12)}& {\displaystyle }{P}_{\mathrm{fm}}& {\displaystyle }=P_{\mathrm{SLHL}}{f}_{\mathrm{fm}}{S}_{\mathrm{fm}},\end{array}$$

where P_SLHL is the total sea-level and high-latitude production rate of the considered nuclide; f_sp, f_sm and f_fm are the fractions of this production rate due to spallation, slow muon capture and fast muon interactions; and S_sp, S_sm and S_fm are the respective scaling factors depending on latitude and elevation. In the simulations presented here, the model by Stone (2000) is used to calculate the scaling factors and fractions, but it is possible to implement the model by Dunai (2000) or the time-varying model by Lifton et al. (2014) accounting for variations in the magnetic field through geological timescales. Topographic shielding is not taken into account as the slopes in the following simulations are lower than 30^∘ (DiBiase, 2018).

The cosmogenic concentration is calculated for each grain at the end of a time step using Eqs. (7) and (8). For a grain in movement during the time step, the mean elevation and mean depth between its initial and final positions are used to calculate the CN production rate. For a grain leaving the grid, the dt of this iteration in Eq. (8) is decreased to account for the fact that the grain has spent only part of the time step in the grid: dt is multiplied by the ratio between the depth of the grain on the starting cell of this time step and the eroded thickness on that cell during this time step.

Production rates for ¹⁰Be, ²⁶Al, ²¹Ne and ¹⁴C are currently implemented in Cidre, but other nuclides can be easily added with specific production models (e.g. ³⁶Cl, ³He and ${}^{\mathrm{78},\mathrm{80}-\mathrm{82}}$Kr; Dunai et al., 2022; Schaefer et al., 2022).

3.3 Initial CN concentration

At the beginning of a simulation, an initial concentration C_o must be specified for each grain. The choice of C_o depends on the particular situation. For example, if we want to model the post-glacier evolution of a topography that was deeply eroded and previously protected from cosmic rays below a thick glacier, C_o=0 at g⁻¹ may be convenient. Another situation may correspond to a surface eroding at a constant rate ϵ [L T⁻¹] for hundreds of thousands years. In this case, a steady-state solution for C_o is given by

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{C}_o& =P_{\mathrm{sp}}{e}^{-\mathit{\rho }{z}^{\prime }/{\mathrm{\Lambda }}_{\mathrm{sp}}}\cdot {\displaystyle \frac{{\mathrm{\Lambda }}_{\mathrm{sp}}}{\mathit{\rho }\mathit{ϵ}+\mathit{\lambda }{\mathrm{\Lambda }}_{\mathrm{sp}}}}+P_{\mathrm{sm}}{e}^{-\mathit{\rho }{z}^{\prime }/{\mathrm{\Lambda }}_{\mathrm{sm}}}\\ & \cdot {\displaystyle \frac{{\mathrm{\Lambda }}_{\mathrm{sm}}}{\mathit{\rho }\mathit{ϵ}+\mathit{\lambda }{\mathrm{\Lambda }}_{\mathrm{sm}}}}+P_{\mathrm{fm}}{e}^{-\mathit{\rho }{z}^{\prime }/{\mathrm{\Lambda }}_{\mathrm{fm}}}\cdot {\displaystyle \frac{{\mathrm{\Lambda }}_{\mathrm{fm}}}{\mathit{\rho }\mathit{ϵ}+\mathit{\lambda }{\mathrm{\Lambda }}_{\mathrm{fm}}}}.\end{array}\end{array}$$

The deeper the grains are set initially in the rock, the less the choice of C_o matters for the following evolution because the production rate decreases exponentially with depth.

3.4 Grain revival

The fate of a grain is to leave the model grid at one of its outlets, where it then becomes a “dead” grain. Every single grain may leave the model grid before the end of a simulation. For some applications, such as the study of the riverine detrital ¹⁰Be evolution at mountain outlets, the flux of grains reaching the outlet must be continuous. One approach is to populate the initial bedrock with a large number of grains at great depth so that there are always grains exhuming during the whole simulation. This is possible but may require long computational times. One alternative is to reuse the grains leaving the model grid at each time step. We set them back to their initial cell, at a random depth between two specified values. Their initial concentration C_o can be assumed to correspond to the steady-state value given by Eq. (13) with the erosion rate corresponding to the erosion rate on that cell calculated at the last time step. This approach allows us to handle a limited number of grains.

The other advantage of reviving grains at their initial location also deals with the theory of the ¹⁰Be-derived catchment mean erosion rate, which we illustrate in the next section. This theory requires that the mean ¹⁰Be concentration of grains at the outlet of a catchment reflect the ratio between the flux of ¹⁰Be atoms and the flux of quartz (Brown et al., 1995). Because the revival of grains occurs more frequently where cells erode faster, there are proportionally more grains coming from these cells that reach the outlet compared with those coming from cells that are eroding more slowly. Statistically, this ensures that the mean ¹⁰Be concentration at the catchment outlet reflects the ratio between the ¹⁰Be flux and quartz flux.

Using the steady-state C_o when a grain is reintroduced in the grid is only an estimation of the true value that should integrate previous variations in the erosion rate, but if the depth at which the grains are set back is deep enough, this approach provides a good compromise between computing time and precision, as shown in the following examples. An appropriated revival depth is below the attenuation length $\frac{{\mathrm{\Lambda }}_{\mathrm{sp}}}{\mathit{\rho }}$ of CN production by spallation. For a granitoid rock with a density of 2.7 g cm⁻², this attenuation length is approximately 65 cm. The production rate decreases exponentially by a factor of 2.7 between the surface and 65 cm. At 1 m the production rate is only roughly 21 % of that at the surface. Deeper down (> 3 m), the production by muons dominates, but the production rate by muons adapts itself more slowly to variations in erosion rates (Braucher et al., 2003). Thus, setting back grains at depths deeper than 65 cm attenuates the error associated with the steady-state assumption for defining C_o if there are high frequency variations in the erosion rate. Conversely, the record of a previously higher erosion rate in the past can be lost in the case of a topography which has evolved very slowly and therefore has undergone a very strong decrease in erosion rate. Again, the deeper the grains are set back, the smaller the bias but also the smaller the number of outgoing grains at each time step.

3.5 Pseudo code for the Cidre erosion–deposition and grain transport algorithm presented in this study

Read the input parameters, initial elevations and grains list
While the specified final time is not reached:
time = time + dt
Rank the cells in the order of decreasing elevation
Do for each cell in this order:
Calculate the slopes from the cell towards all the downstream directions
Calculate the water discharge by summing the incoming water flux and the local precipitation
Calculate the outgoing water flux in each downstream cell direction in proportion to the local slope
(the water discharge qdx is now known on the grid)
Do for each cell in the same order:
Calculate the potential eroded volume of sediment by hillslope processes (Eq. 2)
If the erosion is larger than the available volume of sediment on the cell:
Erode a volume of bedrock according to Eq. 2 weighted by (1-sediment volume/potential erosion of sediment)
Calculate the potential eroded volume of sediment by river processes (Eq. 2)
Calculate the deposited volume of sediment by hillslope processes (Eq. 4)
Calculate the deposited volume of sediment by river processes (Eq. 4)
Calculate the balance between incoming, deposited and eroded volumes and spread the outgoing sediment volume among the downstream cells in proportion to their slope
Add the net elevation change from the difference between erosion and deposition (Eq. 1)
Add the uplift to the elevation
Do for each grain in the grid until the grain is deposited or leaves the model grid:
If the grain does not move:
Update its depth
Else If it moves:
Draw the next cell of the grain among all downstream cells with a probability proportional to the slope
If the grain is not deposited on the next cell:
Continue to move the grain to the next cell
If the grain left the model grid definitively:
If the revival option is true:
Set the grain back to its initial cell at a random depth between specified values
Do for each grain in the grid:
Calculate its CN concentration (Eqs. 7 to 10) in selected nuclides using the mean elevation and depth of its travel during the time step
Save the results if the time fits the user-defined output time
This concludes the code.

4 Example: deriving the catchment mean concentration from ¹⁰Be in uplifting and down-wearing landscapes

4.1 Steady state

We carry out a first test of the algorithm by comparing the true mean erosion rate calculated by Cidre (Eq. 1) in a catchment with the ¹⁰Be-derived mean erosion rate inferred from the mean ¹⁰Be concentration of grains leaving the model grid at each time step. For this, we design a reference simulation using a grid of 100×100 cells with a size of 100 m (10 km²), with only one imposed outlet in a corner fixed at an elevation of 0 m, and starting from a nearly horizontal topography with variations in elevation obeying a Gaussian distribution with a standard deviation of 0.5 m. We impose a constant precipitation rate of 1 m yr⁻¹ on the grid. The computation time step dt is 100 years. The other model parameters are given in Table 1. We impose a first 10 Myr period with a constant uplift rate of 1 mm yr⁻¹ so that a dynamic equilibrium is reached, with a steady topography and a mean erosion rate on the grid equal to the uplift rate (Kooi and Beaumont, 1996). Then, we put 100 000 grains in the bedrock at a depth between 10 and 100 m with a ¹⁰Be concentration C_o of 0 at g⁻¹. These grains are progressively exhumed during 0.1 Myr and then they leave the model grid. Once they leave the model, they are considered “dead” for the rest of the simulation and we do not recycle them. We record the ¹⁰Be concentrations of the grains leaving the model grid at each time step and we calculate their mean ¹⁰Be concentration. In order to infer the mean erosion rate, we use the classic steady-state ¹⁰Be concentration model for a constant erosion rate given by (Lal, 1991; Braucher et al., 2011)

$$\begin{array}{}\text{(14)}& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{ϵ}}{\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }\stackrel{\mathrm{‾}}{C_o}}}(\stackrel{\mathrm{‾}}{P_{\mathrm{sp}}}{\mathrm{\Lambda }}_{\mathrm{sp}}+\stackrel{\mathrm{‾}}{P_{\mathrm{sm}}}{\mathrm{\Lambda }}_{\mathrm{sm}}+\stackrel{\mathrm{‾}}{P_{\mathrm{fm}}}{\mathrm{\Lambda }}_{\mathrm{fm}}),\end{array}$$

where $\stackrel{\mathrm{‾}}{P_{\mathrm{sp}}}$, $\stackrel{\mathrm{‾}}{P_{\mathrm{sm}}}$ and $\stackrel{\mathrm{‾}}{P_{\mathrm{fm}}}$ are the spatially averaged production rates over the model grid. As the CN production rates depend on elevation, we use the topography corresponding to the studied model time to calculate these values.

This model is based on Eq. (13) but neglects the radioactive decay so that an analytical solution exists. At each time step, we calculate the ¹⁰Be-derived erosion rate if there are more than 10 grains that leave the model grid during the time step. At the end of the simulation we can also calculate the average ¹⁰Be concentration of all the grains that went out during the 0.1 Myr steady-state period to calculate a mean ¹⁰Be-derived erosion rate that should equal the constant erosion rate of 1 mm yr⁻¹. We carry out other simulations with different uplift rates (0.01, 0.1, 0.5 and 1.5 mm yr⁻¹) setting the initial depths of the grains so that there is the same number (∼100) of outgoing grains at each time step for all these simulations on average (e.g. between 10 and 160 m for U=1.5 mm yr⁻¹).

Figure 2 shows the comparison between the true mean erosion rate (calculated from the landscape evolution model in Cidre) and the ¹⁰Be-derived erosion rate in each case. When averaged over the whole steady-state period, there is a good fit between the ¹⁰Be-derived and true erosion rates. The ¹⁰Be-derived erosion rate overestimates the true erosion rate for low erosion rates (0.01 mm yr⁻¹ and 0.1 mm yr⁻¹) by about 10 % because the radioactive decay is neglected in order to calculate the ¹⁰Be-derived erosion rate, which is well known (Balco et al., 2008). From one time step to another, there is a variability that depends on the number of outgoing grains as well as on the magnitude of local erosion on the grid. The variability is higher for larger erosion rates (±20 % for U=1.5 mm yr⁻¹), which is consistent with erosion that is dominated by hillslope processes with slopes near the critical slope S_c (Eq. 4). At each time step, a thick layer of tens of centimetres can be removed, including one grain on average in each cell, at different depths depending on the cell and, thus, with very different ¹⁰Be concentrations.

Figure 2 Comparison between the true (calculated from the landscape evolution model in Cidre) black and ¹⁰Be-derived erosion rates at dynamic equilibrium ($U=\stackrel{\mathrm{‾}}{\mathit{ϵ}}$) for different values of uplift rates (i.e. different simulations). The overestimation (shift) of the ¹⁰Be-derived erosion rates for a low erosion rate (10⁻⁵ and 10⁻⁴ m yr⁻¹) comes from the absence of radioactive decay in Eq. (14) used to infer the ¹⁰Be-derived erosion rates, whereas radioactive decay is taken into account in Cidre. Radioactive decay slightly decreases the mean ¹⁰Be concentration. The apparent inferred erosion rate is inversely proportional to the ¹⁰Be concentration (Eq. 14), but because it is calculated by neglecting radioactive decay, the apparent inferred erosion rate is slightly overestimated.

4.2 Transient erosion rate

Here, we test the consistency between the true and ¹⁰Be-derived erosion rates during the transient stage of the topography adaptation to uplift. We use the same parameters as in the reference simulation described in the previous section. We impose a first 5 Myr period with a constant uplift rate of 1 mm yr⁻¹ and a second 5 Myr-period with a constant uplift rate of 0.1 mm yr⁻¹ which is 10 times lower. The (true) evolution of the mean erosion rate is composed of a transient period and then a dynamic equilibrium where ϵ=U on each cell (Fig. 3a). During the first period, the transient period begins with the establishment of a drainage network and then an increase in the slopes leading to an increase in the mean erosion rate with a classic convex curve, before reaching the dynamic equilibrium (Fig. 3a) (Bonnet and Crave, 2003; Carretier et al., 2009). The maximum dynamic equilibrium elevation reaches 1800 m. In the second period, the mean erosion rate decreases to match the lower uplift rate value at the new dynamic equilibrium. The maximum elevation is 340 m during this new equilibrium period.

Table 1 Model parameters used in the reference simulation. Parameters written in bold are those that are varied in the other simulations.

At the beginning of the simulation, we spread 10 000 grains of quartz with a radius of 1 mm (1 per cell on average) at a randomly chosen location and depth, where the depth is between 3 and 6 m. Their initial C_o=0 at g⁻¹. Contrary to previous simulations, where grains left definitively the model grid, we recycle them each time they leave the model grid. The grains are exhumed and transported, and when they leave the model grid, they are set back to their initial cell at random depths between 1 and 2 m. A grain is set back with a ¹⁰Be concentration C_o corresponding to the steady-state concentration (Eq. 13) calculated using the current erosion rate of the cell where the grain is set back. The minimum depth of 1 m ensures that most of the ¹⁰Be acquisition (79 %) occurs later. This depth will be evaluated later.

Figure 3b shows three snapshots of the topography and grains at 0.35, 5 and 10 Myr. In Fig. 3b, the larger cubes show grains that are in movement during the last time step, whereas the smaller ones show grains still in the bedrock at depth. The first snapshot corresponds to the period of drainage network growth (see Fig. 3a). Not all cells are connected to the imposed outlet, and therefore, the divide of the associated catchment is expanding towards the boundary of the model grid. The 5 Myr step is the topography at dynamic equilibrium with a homogeneous and constant erosion rate of 1 mm yr⁻¹. The 10 Myr snapshot corresponds to the dynamic equilibrium of the second period with a final mean denudation rate of 0.1 mm yr⁻¹, and consequently a lower elevation (see Fig. 3a).

Figure 3 Reference simulation. Panel (a)shows the evolution of the mean erosion rate evolution with two periods of constant uplift rates. Panel (b) is a snapshot of three stages in the topographic evolution. The number of displayed grains is divided by 10 for easier viewing.

Figure 4 Reference simulation. Panel (a) shows the evolution of the number of outgoing grains during a time step (100 years) and the mean ¹⁰Be concentration averaged over the outgoing grains. Panel (b) shows the true mean erosion rate calculated by Cidre and the mean ¹⁰Be-derived erosion rate calculated from the mean ¹⁰Be concentration of outgoing grains.

The number of grains leaving the model grid at each time step increases through time in the first period (Fig. 4a) because the mean erosion rate increases (Fig. 4b). During the second period, the mean erosion rate decreases and so does the number of grains leaving the model grid (Fig. 4). The mean ¹⁰Be concentration increases in the first hundreds of thousands of years because, during the progressive establishment of the drainage network, many grains have stayed in the former surface for a long time before being exported out of the catchment. During the uplift of this surface, the ¹⁰Be production rate increases. The long residence time at shallow depth and increasing ¹⁰Be production rate have produced grains with high ¹⁰Be concentrations, which explains the increase in the mean ¹⁰Be concentration during network growth. Once the drainage network is completely established, the ¹⁰Be concentration decreases slightly and stabilizes to a nearly constant value, although the erosion rate increases rapidly: the increase in the ¹⁰Be production rate due to the increasing elevation is compensated by a decrease in the residence time of the grains as the erosion rate increases. This evolution illustrates that a record showing a constant ¹⁰Be concentration may not be a diagnostic for a constant erosion rate. During the second period, the mean concentration in ¹⁰Be increases because, although the elevation and thus the ¹⁰Be production rate decrease, the ¹⁰Be mean concentration is dominated by the longer residence time of grains following the decrease in erosion rate.

Figure 4a shows that the mean ¹⁰Be concentration of grains leaving the outlet does not depend on the number of grains, and high frequency variations in this number from one time step to another generate smaller variations in the mean ¹⁰Be concentration. The continuous and smooth evolution of the mean ¹⁰Be concentration shows that the outgoing grains correctly sample the model grid although they are 14–100 times less numerous than the model cells, even during the transient periods where the erosion rate is heterogeneous.

We now use the mean ¹⁰Be concentration $\stackrel{\mathrm{‾}}{C}$ of the grains leaving the catchment during a time step to calculate the catchment-averaged erosion rate $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ through time using Eq. (14). $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ is not calculated if the number of grains is smaller than 20. Furthermore, $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ is not calculated during the development of the drainage network for practical reasons as the catchment is smaller than the whole model grid. Each ¹⁰Be-derived $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ value corresponds to the average of outgoing grains over a time step of 100 years. Figure 4c shows the good match between the true mean erosion rate and $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ through time. The variation in ¹⁰Be-derived $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ is less than 5 % around the true value for the two periods.

Doubling the initial number of grains decreases the scattering to 1 % around the true value, while halving the number of grains increases it to 5 % (Fig. 5).

Decreasing the time step to 20 years increases the scattering to 5 % because there are fewer grains leaving the model grid during a time step (Fig. 6). There is also a very slight (∼ 1 %) overestimation of the mean erosion rate on average when it becomes larger than about 0.8 mm yr⁻¹ because a small time step increases the probability for a grain to be temporarily stored in a cell. Grains coming from cells far from the outlet may take several time steps to reach the outlet once they are detached from the bedrock. A grain close to the outlet may take less time. Consequently, there are proportionately more grains coming from cells close to the outlet. As these grains are located at lower elevations, with smaller ¹⁰Be production rates and thus smaller ¹⁰Be concentrations, the mean ¹⁰Be concentration of the outgoing grains is slightly underestimated. In turn, the ¹⁰Be-derived erosion rate is slightly overestimated.

Figure 5 Effect of seeding the model with a different number of grains on the mean ¹⁰Be-derived erosion rate calculated from the mean ¹⁰Be concentration of outgoing grains. Panel (a) represents 20 000 grains (two per cell on average). Panel (b)represents 5000 grains (one out of every two cells on average).

Figure 6 Effect of dividing the calculation time step by 10 on the mean ¹⁰Be-derived erosion rate calculated from the mean ¹⁰Be concentration of outgoing grains. The variability in the mean ¹⁰Be-derived erosion rate increases because there is a smaller number of outgoing grains over a smaller time step.

As LEMs can be sensitive to cell size, we tested the result of decreasing the cell size to 25 m, and the time step to 50 years does not change the goodness of fit between the ¹⁰Be-derived and true $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ (Fig. 7).

Figure 7 Effect of decreasing the cell size to 25 m (keeping the same number of cells) on the mean ¹⁰Be-derived erosion rate calculated from the mean ¹⁰Be concentration of outgoing grains. The different shape of the denudation curve compared with the reference simulation in Fig. 4 comes from the increase in slopes in this simulation due to the smaller cell size, leading to a larger contribution of hillslope processes when the slope approaches the critical slope S_c in Eq. (4).

In order to test the effect of catchment size on the calculation of ¹⁰Be-derived $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$, we use a model grid with 200×200 cells, which is four times larger than the reference simulation, but we leave the other parameters as in the reference simulation. We seed the model grid with 40 000 grains, i.e. one per cell on average, as in the reference simulation. Figure 8 shows the good match between the ¹⁰Be-derived $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ and the true value in this case too. The smaller variability in the ¹⁰Be-derived $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ around the true value is due to the fact that there are four times more grains that leave the grid at each time step.

Figure 8 Effect of multiplying the model grid size by four on the mean ¹⁰Be-derived erosion rate calculated from the mean ¹⁰Be concentration of outgoing grains. The number of grains is also multiplied by four to have one grain per cell on average as in the reference simulation.

We now rerun the reference simulation but divide the uplift rate by 10 in the two periods. Figure 9a shows that the ¹⁰Be-derived $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ is overestimated by approximately 2 % during the dynamic equilibrium of the first period (4–5 Myr). This is due to the model in Eq. (14) used to calculate $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ that neglects the radioactive decay. Here, radioactive decay influences the ¹⁰Be concentration because the clast residence time in the bedrock is longer for small uplift rates (Lal, 1991). In the second period, the ¹⁰Be-derived erosion rate is overestimated by roughly 20 % during the transient adjustment to a smaller uplift rate (6–8 Myr). This overestimation comes from the delayed response of CN for a low erosion rate. The response time to adjust to a new erosion rate ϵ is about four times $\mathrm{1}/(\mathit{\lambda }+\frac{\mathit{\rho }}{{\mathrm{\Lambda }}_{\mathrm{sp}}}\mathit{ϵ})$ when only considering CN produced by spallation (Lal, 1991). With ϵ between 0.1 and 0.01 mm yr⁻¹, the response time is between 60 000 and 600 000 years. Consequently, the CN concentration is out of phase during the rapid transient decrease in the erosion rate for the second period. This delay was not observed in the reference simulation because the response times were 10 times shorter. Furthermore, the variability in ¹⁰Be-derived $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ tends to increase for smaller uplift rates because the number of outgoing grains is smaller at each time step. In this simulation, the number of outgoing grains varies between 100 and 10. It is remarkable that with such a small number of grains (particularly during the second period with a very low uplift rate), the ¹⁰Be-derived $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ remains a good estimate of the true $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$. When the number of grains is multiplied by four, this decreases the variability (Fig. 9b).

As there is a rapid decrease in the erosion rate, but a long residence time of grains in the last metres because the erosion rates are low in these last simulations, one may question the chosen revival depths of the grains because they may be too shallow to record the previous decrease in erosion. We test this by changing the revival depth of the grains to between 2 and 4 m. Figure 9c shows that the delay in the ¹⁰Be-derived erosion rate between 6 and 8 Myr is similar, with a larger variability due to the smaller number of outgoing grains at each model time.

Figure 9 Effect of an uplift rate that is 10 times smaller. In (a) the number of grains is 10 000 as in the reference simulation. In (b) the number of grains is 40 000. In (c) the number of grains is 10 000, but the (revival) depth at which the grains are set back to the model grid is between 2 and 4 m instead of 1 and 2 m. In the three cases, the slight 1 %–2 % overestimation of the ¹⁰Be-derived erosion rate during the first dynamic equilibrium stage (4–5 Myr) comes from the model used to derive erosion rates that neglect the ¹⁰Be radioactive decay, which only matters for low erosion rates. During the second period (6–8 Myr), there is a significant overestimation of the true erosion by approximately 20 %. This overestimation comes mainly from the delayed response of the CN concentration for a low erosion rate. In (c) the variability is larger because there are fewer grains going out at each model time step.

In order to further test if a revival depth between 1 and 2 m gives good results, we carried out three final simulations using the parameters of the reference simulation but imposing a constant uplift rate of 0.1 mm yr⁻¹ and an oscillatory precipitation rate between 0.5 and 1 m yr⁻¹ during 10 Myr. The three simulations correspond to three different periods of oscillation constituting the Milankovitch cycles: 23, 41 and 100 kyr. It is predicted that the ¹⁰Be-derived erosion rate should be shifted with a lag that increases with the period of oscillation (Schaller and Ehlers, 2006). Figure 10 illustrates several cycles in the last 200 kyr of these simulations. It shows that the lag between the ¹⁰Be-derived and true erosion rate signals increases with the oscillation period. The lags are very similar to the ones found by Schaller and Ehlers (2006) in the case of one-point source simulations for the same oscillation periods, magnitude and true erosion rates (see Fig. 5d, e and f in Schaller and Ehlers, 2006). We also observe that the amplitude of the ¹⁰Be-derived erosion rate decreases when the frequency of the true erosion rate variation increases, consistent with the results of Schaller and Ehlers (2006). Indeed, when a grain reaches shallow depths (< 1 m) during a low erosion rate period, its ¹⁰Be concentration is relatively high. If the grain is then rapidly exhumed, it will reach the surface with a concentration that is too high compared with what it would have been with a high rate of erosion. Once detached, if we use this concentration to determine an erosion rate, we underestimate the erosion rate. This memory effect causes the cosmogenic signal to be damped out.

Figure 10 Simulations using the parameters of the reference simulation but a constant uplift rate of 0.1 mm yr⁻¹ and oscillations of the precipitation rates between 0.5 and 1 m yr⁻¹ with different periods: (a) 23 kyr, (b) 41 kyr and (c) 100 kyr. The observed lags ϕ are similar to the ones calculated by Schaller and Ehlers (2006) for similar conditions but for a theoretical one-point source catchment. This consistency shows that the revival procedures of grains at depths between 1 and 2 m give robust results.

5 Discussion

5.1 Limitations and advantages

The main limitation of this approach deals with the computational time when the number of grains is large. Some study cases may require a huge number of grains. In practice, the feasibility of these simulations will depend on the computation facilities. To give an indication, the reference simulation with 100×100 cells and 10 000 grains runs in half an hour on a personal computer. The simulation with a grid and number of grains four times larger runs in 2.5 h. Parallelism of the grains' calculation is straightforward as they are independent of each other and would decrease the running time. As the precision of the results depends on the number of grains, the model can be adapted according to the studied process. For example, a higher density of grains is required to study the CN detrital signal associated with climate variations for the 20 kyr period compared with the density required for variations in the 100 kyr period. It is possible to determine a simple quality control of the results by setting a minimum number of grains leaving the model grid to interpret the CN concentrations. Beyond the number of grains, the choice of time step and revival depths have an influence on the resulting cosmogenic concentrations, but the simulations presented here show that these choices do not imply variations in inferred erosion rates larger than 10 % around the true value (when the CN and the erosion rate signals are in phase). This error remains smaller than or equal to the uncertainty of the CN-derived erosion rate in most cases, given the external uncertainties on the CN production rates (Balco et al., 2008).

In the current implementation, grains do not physically erode during their transport (but they can weather and decrease in size chemically; Carretier et al., 2018). In the real world, attrition can significantly change the CN concentration within and at the surface of decimetric pebbles. In the model, this would require tracking the CN concentration of the layers within the pebbles. This was done in the model of Carretier and Regard (2011), which served as a basis for the grain model in Cidre; therefore, it would be straightforward to implement it in Cidre as well.

Modelling deep landslides (Campforts et al., 2020) would be a limitation to the grain revival process. If the specified maximum depth at which a grain is set back to the grid is much smaller than the thickness of the landslides, there is a risk that the mean CN concentration of grains in the landslide would overestimate the mean value of the landslide. In such cases, the maximum revival depth of the grains must be adequately adapted to be greater than the maximum landslide thickness.

The revival of dead grains is a useful approach for the statistics of CN at the catchment outlet, but this approach is an approximation. It is fundamental that the depth at which the grains are set back to the model grid be much larger than the attenuation length of production by spallation ($\frac{{\mathrm{\Lambda }}_{\mathrm{sp}}}{\mathit{\rho }}\sim$ 65 cm in granitic rocks). If not, there is a risk that the revival CN concentration C_o of the grain would be a very poor approximation of the erosion history of its cell if the erosion rate has decreased very quickly and with a huge magnitude. Figure 10 shows that setting back the outgoing clasts at depths larger than 1 m with the ¹⁰Be concentration corresponding to the equilibrium value with the local erosion rate gives robust results.

In a Lagrangian formulation of CN concentration evolution, the approach by discrete grains has advantages. The alternative Eulerian approach would require calculating and tracking the concentration of CNs at different depths below each cell at each time step (Niemi et al., 2005). Petit et al. (2023) present a simplified approach implemented in the LEM Badlands model based on the computation of the CN concentration at the earth's surface only without calculating the CN concentration profile at depth. This level of simplification is required in this case for practical reasons of computation in an LEM. Tracking the whole depth profile evolution in the whole model grid cell where erosion and deposition can alternate at each iteration would be prohibitive in terms of computational time and would be difficult to implement (Petit et al., 2023). One advantage of the Lagrangian approach is that it can tackle more complex erosion–deposition scenarios compared with the Eulerian approach (Knudsen et al., 2019). For example, it is simpler to calculate the evolution of a given grain experiencing a complex erosion–deposition history of the surface above it rather than the evolution of a (deep) depth profile of the concentration for each time step below a surface that is constantly changing in terms of elevation. From the grain's point of view, it is easy to calculate the evolution of CN in the grain when it is stored, buried and eroded again stochastically, or when the soil above it is alternatively eroded or buried, as only the grain's depth has to be adapted. (Note that the difference in density between rock and sediment is not taken into account when calculating the CN production rate in the simulations presented in this work, but it is possible to simply implement this.) The Eulerian approach requires tracking the sources, and averaging the concentration from these different sources in cells during the transport of sediment, assuming perfect mixing. This process loses some information about the distribution of CN for a population of grains. On the contrary, the distribution of the CN concentration is fully conserved by treating the grains separately. One drawback to the grain-by-grain approach is that the average CN concentration of the transported sediment is not precisely known, with an uncertainty that decreases with the number of grains. On the other hand, when the full distribution of the CN concentrations of a population of grains is known, then it becomes possible to study part of the stochasticity of erosion–deposition processes on hillslopes and in rivers, including the long-term temporary storage of some grains that may contribute significantly to the CN concentration average (Carretier et al., 2019). The potential applications of this are described below.

5.2 Potential applications

The method using CN in riverine sediment to quantify modern and palaeo catchment-averaged denudation rates assumes that the CN concentrations of any grain have been entirely acquired on the hillslopes. Yet, the temporary storage and recycling of sediment grains between the eroding sources and sedimentary basins may change their CN concentration depending on the considered nuclide or the climatic context (Wittmann and VonBlanckenburg, 2009; Carretier et al., 2019). Storage and recycling can delay or obscure the transmission of an erosion signal from source to sink (Carretier et al., 2020; Tofelde et al., 2021). Detrital signals in the CN concentrations can also be used positively to study the erosion–deposition processes on hillslopes (Slosson et al., 2022), burial histories in basins (Balco et al., 2013; Sanchez et al., 2021), the distribution of residence times and transport lengths in rivers (Carretier et al., 2019) or to infer the palaeo-denudation rates of mountains (Charreau et al., 2011). Nevertheless, linking a CN detrital signal with landscape evolution is not straightforward. A forward model, such as Cidre, which simulates CN in distinct grains that move stochastically, would help. For example, the conclusions of several studies on palaeo-denudation rates over the Pliocene–Pleistocene periods suggest either constant or modest variations in the denudation rate in Asia or in the Andes (Charreau et al., 2011; Madella et al., 2018; Lenard et al., 2020; Charreau et al., 2021), or a significant change recorded in the last glacial cycle in river terraces or at the outlet of small mountainous catchments (Schaller et al., 2002; Mariotti et al., 2021). Some factors that should be explored numerically with a numerical model, such as Cidre, to better understand the CN detrital signals (Petit et al., 2023) include the following: the effect of catchment size, uplift rate or grain size; whether a flood plain is present or not; the frequency of climatic variations; and variations in the elevation of the eroded sources.

Cidre includes the possibility to produce a regolith corresponding to the weathering of the underlying bedrock (Carretier et al., 2014, 2018). Once a grain is located in the regolith during its exhumation on the hillslopes, it continues to be exhumed towards the surface until it is detached. Alternatively, its depth could be chosen randomly at each time step to simulate bioturbation or physical creep within the soil. This option is easy to implement and would make it possible to analyse the effect of different soil processes on the riverine detrital CN signal within the framework of reservoir theories (Mudd and Yoo, 2010).

One of the advantages of modelling grains is the simplicity in tracking the evolution of distinctive cosmogenic nuclides, taking advantage of their different radioactive decay rates. When grains experience complex temporary burial and recycling histories, their initial concentration ratio on hillslopes varies downstream. It would be also quite simple to implement the evolution of the meteoritic ¹⁰Be concentration at the surface of the grains, or their optically stimulated luminescence (OSL) dosimetry based on Guyez et al. (2023). The distributions of various CN concentrations in a river sample, e.g. ¹⁰Be, ²⁶Al and ¹⁴C, contain information about possible past changes in the catchment-averaged erosion rates and residence times in fluvial systems (Repasch et al., 2020; Ben-Israel et al., 2022). We need to establish a link between erosion rate changes and the detrital CN concentrations in different nuclides and potentially other properties such as OSL dosimetry. Doing so may allow us to address recent changes over the past few centuries associated with natural and anthropic modifications of the landscape.

6 Conclusions

We present a new coupling of the landscape evolution model Cidre with a model of CN concentrations in individual grains. The algorithm is tested by comparing the catchment-averaged erosion rate derived from the ¹⁰Be concentration of grains leaving an uplifting catchment and the true catchment-averaged erosion rate calculated by Cidre. The main limitation is the number of grains that has to be set in a simulation to achieve the desired precision. The catchment-average erosion rates are approximated to within 5 % uncertainty in most of the cases with a limited number of grains. This Lagrangian approach allows us to fill the gap that exists between landscape modelling, which is used to help understand variations in the elevation and erosion of landscapes, and field data, which often correspond to the CN concentrations of grains in a soil or river sample.