**Research article**| 04 Jul 2022

# Theoretical and numerical considerations of rivers in a tectonically inactive foreland

Stefan Hergarten

**Stefan Hergarten**Stefan Hergarten

- Institut für Geo- und Umweltnaturwissenschaften, Albert-Ludwigs-Universität Freiburg, Albertstr. 23B, 79104 Freiburg, Germany

- Institut für Geo- und Umweltnaturwissenschaften, Albert-Ludwigs-Universität Freiburg, Albertstr. 23B, 79104 Freiburg, Germany

**Correspondence**: Stefan Hergarten (stefan.hergarten@geologie.uni-freiburg.de)

**Correspondence**: Stefan Hergarten (stefan.hergarten@geologie.uni-freiburg.de)

Received: 24 Feb 2022 – Discussion started: 09 Mar 2022 – Revised: 24 May 2022 – Accepted: 08 Jun 2022 – Published: 04 Jul 2022

Modeling the dynamics of alluvial rivers is theoretically and numerically more challenging than modeling erosion of bedrock channels in active mountain ranges. As a consequence, the majority of the approaches developed in the context of alluvial rivers are one-dimensional. However, recent advances in the numerics of fluvial landform evolution models allow for two-dimensional simulations of erosion and sediment transport over time spans of several million years. This study aims at finding out fundamental properties of rivers in a tectonically inactive foreland of a mountain range by investigating a simple reference scenario theoretically and numerically. This scenario consists of a mountain range and a foreland in a quasi-steady state wherein the material eroded in the mountain range is routed through the foreland. In order to understand the properties of foreland rivers, a subdivision into two classes – carriers and redistributors – is introduced. Carriers originate in the mountain range and are thus responsible for the large-scale sediment transport to the ocean. In turn, redistributors are rivers whose entire catchment is located in the foreland. Using the concept of carriers and redistributors, it is shown that the drainage network in the foreland permanently reorganizes so that a steady state in the strict sense is impossible. However, the longitudinal profiles of carriers are described well by a steady-state approximation. Their concavity index is considerably greater than that of rivers in the mountain range. Carriers predominantly deposit sediment at high rates, while redistributors erode at much lower rates. Despite the low erosion rates, the sediment flux from redistributors into carriers is a major component of the overall sediment budget and finally the main driver of the highly dynamic behavior of the carriers.

Fluvial deposits are among the most important records of Earth's tectonic and climatic history. Numerical models describing the physical processes controlling sediment production, transport, and deposition have become essential tools in this field. Finding out how perturbations in the depositional environment (e.g., changes in sea level) or in the source region (e.g., changes in precipitation or tectonic uplift) propagate through the system has been one of the most important applications of such models (e.g., Armitage et al., 2011, 2013; Mouchené et al., 2017; Yuan et al., 2022).

It seems, however, that our understanding of large-scale and long-term sediment deposition still lags behind our understanding of erosion processes in active mountain ranges. As reviewed by Romans et al. (2016) and Tofelde et al. (2021), large parts of our knowledge about the source-to-sink sediment transfer to the oceans are still on a conceptual level. However, analytical and numerical models have become increasingly important for understanding the general properties of alluvial rivers. Among the recently developed approaches, there are one-dimensional models focusing on individual rivers (e.g., Bolla Pittaluga et al., 2014; Blom et al., 2016, 2017; Malatesta et al., 2017; Wickert and Schildgen, 2019; Braun, 2022) and two-dimensional landform evolution models (e.g., Carretier et al., 2016, 2020; Yuan et al., 2019, 2022). In particular, the latter may improve our understanding of sediment dynamics considerably in the near future. At the moment, however, little is known about the self-organization of alluvial rivers and its effects on rates of erosion and deposition. Even 40 years after Sadler (1981) quantitatively described the decrease in rates of accumulation with increasing time spans of observation, our understanding of this phenomenon is still incomplete.

Modeling sediment transport and deposition seems to be more challenging than modeling erosion. While all models reviewed by Coulthard (2001), Willgoose (2005), and van der Beek (2013) involve a sediment balance, it was already pointed out by Howard (1994) and Kooi and Beaumont (1994) that simulating sediment transport in large rivers requires small time increments. This leads to high computing effort and seriously limits the applicability of the models to large-scale problems with a reasonable spatial resolution. This problem is less severe in the limit of detachment-limited erosion (Howard, 1994) whereby it is assumed that all particles entrained by the river are immediately swept out of the system. Even a fully implicit scheme is available here, which in principle allows for arbitrarily large time increments and thus for large-scale simulations over long time spans (Hergarten and Neugebauer, 2001; Braun and Willett, 2013). Presumably owing to its theoretical and numerical simplicity, the concept of detachment-limited erosion has been applied in numerous studies of landform evolution. However, the central assumption that transported sediment has no effect on landform evolution limits the applicability of detachment-limited erosion to mountain streams.

Concerning the numerics of large-scale models including sediment transport, considerable progress was achieved recently. Yuan et al. (2019) combined the implicit scheme for erosion with a fixed-point iteration for the sediment fluxes. Their scheme achieves high efficiency as long as the conditions are not too close to the transport-limited regime. The concept of transport-limited erosion assumes that the actual sediment flux of a river is always equal to the so-called transport capacity. This means that the rate of erosion or deposition instantaneously adjusts in such a way that the sediment flux equals the transport capacity. In contrast to the iterative scheme proposed by Yuan et al. (2019), the fully implicit scheme introduced by Hergarten (2020) even covers the entire range from detachment-limited to transport-limited erosion (and sediment deposition) at a constant numerical efficiency.

The recent numerical developments allow for large-scale simulations including sediment transport over long time spans and should thus also be able to improve our understanding of sedimentary systems. As a first result, Yuan et al. (2019) observed a permanent reorganization of the drainage pattern in a foreland region without uplift and subsidence over long times even under constant conditions. This reorganization causes an autocyclicity in erosion and deposition of sediments. From a theoretical point of view, the question arises of whether such self-organizing systems are completely irregular or whether their average behavior can be described by simple relations or by one-dimensional models. The central question from a more practical point of view concerns the suitability of sediment archives for recording changes in climate and tectonics. In this context, not only the strength of the autocyclic dynamics compared to allocyclic influences is relevant, but also the respective spatial and temporal scales.

The present study goes a step toward understanding the dynamics of a system that is completely dominated by autocyclic sediment dynamics. In the first part, it will be shown that steady-state topographies in the sense of constant elevation at each point are impossible in absence of uplift and subsidence so that there must be autocyclic aggradation and incision. In the next step, properties of hypothetical steady-state river profiles will be investigated. As a main finding, such rivers capture the properties of large rivers in a regime of permanent reorganization quite well on average, although they cannot be stable over long times. The second part of the paper provides preliminary estimates for rates of erosion and deposition as well as for the timescale of network reorganization.

This study addresses the simplest scenario of rivers in a tectonically inactive foreland. A rectangular domain with a north–south-oriented mountain range at the center and a foreland at each side is considered (Fig. 1). While the mountain range is uniformly uplifted at a constant rate, the foreland regions are neither uplifted nor subsiding. The northern and southern boundaries are periodic, while the eastern and western boundaries are kept at zero elevation and are interpreted as the coast of an ocean.

The open-source landform evolution model OpenLEM
(Hergarten, 2022a) is employed for all numerical
simulations. Since the focus is on a minimum scenario, none of the components
of OpenLEM beyond fluvial erosion and sediment transport,
such as lithospheric flexure
and orographic precipitation (Hergarten and Robl, 2022),
are used.
The fluvial model implemented in OpenLEM is presumably the simplest model of
large-scale
fluvial erosion and sediment transport. Several formulations of this model
were proposed, which are all similar in spirit or even partly
mathematically equivalent: the undercapacity model (Kooi and Beaumont, 1994),
the linear decline model (Whipple and Tucker, 2002), the *ξ*–*q* model
(Davy and Lague, 2009), and the shared stream-power model
(Hergarten, 2020). In this study, the shared stream-power formulation

is used, where *E* is the erosion rate, *Q* the sediment flux (volume per time),
*A* the upstream catchment size, and *S* the channel slope. The model
involves four parameters *K*_{d}, *K*_{t}, *m*, and *n*.
The term *A*^{m}*S*^{n} is often called the stream-power term, and the model
implements the idea that this term is used jointly by erosion and
sediment transport.

While the equation for the change in surface elevation *H* at a given
uplift rate *U* is straightforward,

an additional balance equation for the sediment fluxes is required in order
to obtain a closed system of equations.
Assuming that each node *i* of a discrete grid delivers its entire sediment flux
*Q*_{i} to a single neighbor, this sediment balance equation reads

where *s*_{i} is the size (area) of the respective grid cell. The right-hand
side of Eq. (3) is a discrete representation of the divergence
operator; the sum extends over all neighbors that deliver their sediment
flux to the cell *i*. This sum is empty for nodes that do not receive sediment,
which may be seen as internal boundaries. For consistency,
the same routing scheme should be used for water (catchment size *A* and direction
of channel slope *S*) and sediment flux. The implementation in OpenLEM
uses the D8 routing scheme (O'Callaghan and Mark, 1984) on a regular grid.

In the absence of transported sediment (*Q*=0), the shared stream-power model reduces to

which is the stream-power incision model widely used in the context
of detachment-limited erosion. The parameter *K*_{d} is called
erodibility; the subscript emphasizes the relation
to detachment-limited erosion. In turn, neither erosion nor deposition
takes place (*E*=0) if the sediment flux is

So this term defines the transport capacity, wherein the transport
coefficient *K*_{t} describes the ability to transport sediment
at given *A* and *S*. The shared stream-power model turns into the
stream-power incision model for *K*_{t}→∞ and into
a transport-limited model for *K*_{d}→∞. In the latter
limit, the sediment flux is always equal to the transport capacity defined
by Eq. (5). As discussed by Hergarten (2020),
the system of equations reduces to a single diffusion equation then.

For spatially uniform erosion, the sediment flux is *Q*=EA, and Eq. (1)
collapses to a form analogous to the stream-power incision model (Eq. 4) with an
effective erodibility *K* according to

River profiles follow the relation

with $\mathit{\theta}=\frac{m}{n}$ then. The exponent *θ* in Eq. (7) is called
concavity index in the context of analyzing river profiles.
Concavity indexes of real rivers have been investigated in numerous studies,
starting from the seminal work of Hack (1957). Values *θ*≈0.5
are typically found for rivers at uniform erosion;
either *θ*=0.45 or *θ*=0.5 is often used as a reference value
(e.g., Whipple et al., 2013; Lague, 2014). So the ratio of the exponents *m* and *n* is
constrained quite well by the concavity of real-world rivers.

The absolute values of the exponents *m* and *n* are, however, more uncertain
than their ratio since they cannot be determined from the shape
of river profiles under uniform erosion. Assuming *n*=1 simplifies both theoretical
considerations and the numerical implementation since the model is linear
with regard to the topography then.
In turn, the results compiled by Lague (2014) as well as some recent
studies (Harel et al., 2016; Hilley et al., 2019; Adams et al., 2020) rather suggest *n*>1.

In this study, the linear version of the shared stream-power model
(*n*=1) is used for numerical reasons. The fully implicit scheme
introduced by Hergarten (2020) is not only stable at
arbitrary time increments, but also avoids oscillations
in elevation or rates of erosion and deposition, e.g., if a river is suddenly exposed to a
large sediment flux after an avulsion event. It is, however, restricted
to *n*=1. The implementation in OpenLEM
contains a semi-implicit extension for *n*>1, which is still stable at large
time increments but not able to avoid oscillations completely.
While these oscillations are not a problem in many applications, they
affect the short-term rates
of erosion investigated in Sect. 8 and may even cause an artificially
increased frequency of avulsion events. While some simulations
were also performed for *n*=2, the results are not included in this paper.

Since the foreland region
will be covered by alluvial deposits, different lithologies are assumed
in the two domains, in contrast to the recent study of Yuan et al. (2022). As the simplest
approach, it is assumed that *K*_{t} is the same in both domains, while
the values of *K*_{d} differ. The ratio of *K*_{d} and *K*_{t}
is mathematically equivalent to the parameter Θ or *G* (depending on the
notation) in the *ξ*–*q* model. The results obtained by Davy and Lague (2009) and
Guerit et al. (2019) suggest *K*_{d} *⪆* *K*_{t} for *n*=1.
For simplicity, *K*_{d}=*K*_{t} is assumed in the mountain range,
which was assumed by Yuan et al. (2022) for the entire domain.

We will see in Sect. 3 that the choice *K*_{d}=*K*_{t} for the
mountain range has a minor effect on the rivers in the foreland.
In turn, the parameter choice for the
foreland is more critical. As discussed by Hergarten (2021), *K*_{d}
refers to the properties of the actual riverbed. As a consequence, *K*_{d}
should be much larger than *K*_{d} of the bedrock if a previously deposited
thick alluvial cover is eroded. In an environment of deposition (*E*<0),
Eq. (1) should even be replaced by the transport-limited version
(*K*_{d}→∞). Otherwise, assuming a lower erodibility
*K*_{d} would reduce the rate of deposition if all other parameters
remain constant, which would not make much sense.
In order to keep the model as simple as possible, the transport-limited
end-member (*K*_{d}→∞)
is used in the entire foreland region, no matter whether
the rivers are actually depositing sediments or eroding.
An alternative scenario will be considered in Sect. 10.

The implicit scheme inhibits the formation of lakes in individual streams. However, lakes may occur at confluences. If a large amount of sediment is deposited in a large river, its elevation may exceed that of its tributaries. Then the lowermost ranges of the tributaries turn into lakes. Formally, this leads to an upstream sediment flux so that the lowermost ranges of the tributaries serve as accommodation area for the sediments. Since this phenomenon is not unrealistic and does not affect the stability of the scheme, there is no need to inhibit these backward sediment fluxes. It should, however, be noted that this effect does not add floodplains to the model in a realistic way since it is only based on the topography of the bed and does not take into account the water level.

In order not to introduce unnecessary constraints, a nondimensional representation
is used in all simulations and theoretical considerations.
Assuming *n*=1, the choice *m*=0.5 (so *θ*=0.5) is convenient since it
avoids odd physical units of *K*_{d} and *K*_{t} and simplifies the
scaling from nondimensional to real-world properties. For this choice of
*m* and *n*, the dimensions of *K*_{d}, *K*_{t}, and *K* (Eq. 6)
are inverse time so that each of them can be used for defining the timescale independently
of the spatial scales. Since *K*_{t} is constant in the entire domain,
the nondimensional representation is based on *K*_{t}=1.
Then one nondimensional time unit
corresponds to an absolute time span of $T=\frac{\mathrm{1}}{{K}_{\mathrm{t}}}$ (with the real-world
value of *K*_{t}). However, real-world estimates are typically obtained
from the steepness of rivers at a given erosion rate and thus refer to
*K* (Eq. 6) instead of *K*_{d} or *K*_{t}. Assuming
*K*_{d}=*K*_{t} leads to

With the value *K*=2.5 Myr^{−1} used by Robl et al. (2017),
this would yield a timescale
*T*=200 000 years. So some 100 000 years should be a realistic magnitude
for one nondimensional time unit.

The vertical length scale *L*_{v} is defined by the uplift rate
*U* and the timescale *T* in the form

Using a nondimensional uplift rate *U*=1, a convenient vertical
length scale of *L*_{v}=100 m would emerge at a real
uplift rate of 0.5 mm yr^{−1}
for *T*=200 000 years.

The horizontal length scale is arbitrary and independent of all model parameters and of the other scales for $\frac{m}{n}=\mathrm{0.5}$. In principle, this allows for an arbitrary scaling of the horizontal coordinates. However, the pixel size of the model imposes two limitations. First, hillslope processes are not taken into account so that the shared stream-power model is assumed to hold even for single-pixel catchments. Therefore, the pixel size should not be too small. More important for this study, the pixel size also defines the area over which deposited sediments are distributed since the model does not explicitly include floodplains.

All simulations were performed on a grid with 8192 rows and 5632 columns
(Fig. 1), starting from a flat topography superimposed by a small random disturbance
with a range of 10^{−4}. The mountain range has a width of 512 cells. So each side
of the foreland is 5 times as wide as the mountain range. This choice ensures
that a large part of the foreland consists of an alluvial plain that is not
affected strongly by the properties of the mountain range and its alluvial
fans. In turn, the width of the mountain range is large enough to generate rivers
with a wide spectrum of catchment sizes in the mountain range (up to about
100 000 pixels).

The choice of the time increment *δ**t* is not trivial, although
the fully implicit scheme implemented in OpenLEM is unconditionally
stable for any *δ**t*. As discussed by Hergarten (2020),
the limitation arises from changes in flow direction rather than from the
accuracy of the scheme itself. Since each cell can change its flow direction
only once per time step, large time increments slow down the dynamics of the
system. Preliminary tests with a similar setup revealed that this effect
becomes dominant for *δ* *⪆* 10^{−2}. In order to stay clearly
below this limit, $\mathit{\delta}t={\mathrm{2}}^{-\mathrm{10}}\approx {\mathrm{10}}^{-\mathrm{3}}$ is used in this
study, although the results are by far not independent of *δ**t*
then (see also Sects. 3 and 9). According to
the timescale estimated above, *δ**t* is in an order of magnitude
of some 100 years.
Using much smaller values of *δ**t* would not only increase the
numerical effort, but would also question the applicability of the model since
we would probably have to take into account the effects of individual
flood events on river avulsions. However, simulating individual events
is outside the scope of the stream-power concept.

A strong reorganization of the drainage pattern occurs
during the first phase of the simulation. This reorganization slows down in the tectonically
active region with increasing incision of the rivers and has
almost ceased at *t*≈10. Afterwards, there is little
reorganization in the mountain range, although the topography is still
far from equilibrium. Only 32 changes in flow direction occur
at points with a catchment size between 1000 and 10 000 pixels
(the largest catchment in the mountain range is about 100 000 pixels large)
from *t*=10 to the end of the simulation at *t*=500, which means less
than 1 change per 15 time units.

In turn, the reorganization of the drainage network
continues in the foreland. On average, 2.5 % of all sites
change their flow direction in each time step, while avulsions of large
rivers take place more frequently. Even 24 % of all sites with
*A*≥100 000 change their flow direction in each time step on
average. However, the respective frequencies of avulsion depend
on the time increment *δ**t*. Performing a shorter simulation
with an 8 times smaller $\mathit{\delta}t={\mathrm{2}}^{-\mathrm{13}}$ yielded percentages
of 0.51 % and 9.9 %, respectively. So the mean rate of avulsions taken over
all sites increases moderately by a factor of 1.6 if *δ**t*
is reduced by a factor of 8. For the large rivers, however, it
is a factor of about 3. However, we will see in
Sect. 9 that the effect of *δ**t* on
the long-term reorganization of the drainage pattern is weaker
than it seems here. Nevertheless, we should keep in mind that
the frequency of avulsions of large rivers in particular depends on the time increment
and that the frequency of large floods may play an important part
in nature.

The highest mean and peak elevations
are reached in the mountain range at *t*≈80. Afterwards,
the topography is still not constant but decreases very slowly
due to the slow network reorganization. This effect was described
for the stream-power incision model by Robl et al. (2017), but
here it is not relevant since it is much slower than the ongoing network
reorganization in the foreland.
Finally, the time span from *t*=100 to *t*=500 was used
in the following analyses; most of the results were obtained
by averaging over 401 equally spaced snapshots (Δ*t*=1).

Figure 2 shows a swath profile of the mean topography
over the considered time span. Minimum and maximum values are not
taken over *x*_{2} and *t* but only over *x*_{2} and then averaged
over all snapshots.
The foreland topography becomes increasingly steep
close to the mountain range and reaches a mean value of about
5.6, corresponding to about one-fifth of the peak elevation and one-third
of the mean elevation along the center of the mountain range.

The topography obtained under transport-limited conditions in the mountain
range (*K*_{t}=1 and *K*_{d}=∞ as in the foreland)
is also shown in Fig. 2 for comparison (however, based on only 160 snapshots). This change mainly affects
the topography in the mountain range. While the effective erodibility
was *K*=0.5 according to Eq. (6) before, it is twice as high (*K*=1) for
transport-limited conditions. As a consequence, equilibrium channel slopes and thus
relief are 2 times lower. As discussed by Hergarten (2021),
*K*_{d} determines the speed at which knickpoints migrate upstream
and thus the response of the sediment flux from the mountain range
to changes in topography in the foreland. However, the feedback
on the topography in the foreland appears to be minor. In this context, we
should keep in mind that assuming transport-limited erosion in the mountain region
is an unrealistic extreme scenario, although the ratio of *K*_{d} and
*K*_{t} is still not constrained well from real-world data. So these results
suggest that a moderate deviation from *K*_{d}=*K*_{t}
will have an effect on the height of the mountain range compared to the topography
of the foreland but will not affect the properties of the foreland seriously.

Figure 3 provides a more detailed analysis of the mean topography
(averaged over both sides).
The relief (maximum minus minimum elevation taken in the *x*_{2} direction)
is greater than the mean elevation in the mountain range, corresponding
to deeply incised valleys and high peaks. This relation is inverted
in the foreland, corresponding to the rapid decline in maximum elevation
visible in Fig. 2. In the orange-colored region (see
Fig. 1), relief is even less than one-third of the mean
elevation. So the topography is quite smooth here.
This decrease in relief goes along with a distinct minimum
in the correlation of the topography with the topography
at the edge of the mountain range ($\left|{x}_{\mathrm{1}}\right|=\mathrm{256}$).
Surprisingly, the correlation
slightly recovers in the
green and turquoise domains, and a weak positive correlation persists
in the entire domain.

Further away from the mountain range, both mean elevation and relief decrease toward the ocean. The decrease in relief is, however, slower than the decrease in mean elevation. So the topography becomes smoother toward the ocean on an absolute scale but rougher in relation to the mean elevation.

Figure 4 gives a first insight into the dynamics of the
rivers in the foreland. The river considered here is not the biggest
river overall (which is not the same at all times) but the river
with the largest catchment at the edge of the foreland (orange-colored
line in Fig. 1). Since
network reorganization in the mountain range is weak, this river
remains the same for all times *t*≥100 considered here;
changes take place only in the uppermost part of its catchment.

After leaving the mountain range, the area covered by the river widens
moderately. As discussed above, the rivers are confined in narrow gorges
in the mountain range, but their relief rapidly decreases in the foreland,
allowing for wider valleys. At some point, however, the behavior changes
abruptly to form a large alluvial fan. Directions of flow vary by more
than 90^{∘} here, while the river itself is more or less straight.

In contrast to typical alluvial fans in real-world topographies, the alluvial fan observed here is not sharply bounded downstream but rather dissolves by systematically diverting the river toward the ocean. This difference is related to the topography already being in a quasi-steady state wherein the fan is no longer growing. Systematically analyzing the predicted sizes of the alluvial fans and comparing them to real-world data would be interesting. The data reviewed by Blair and McPherson (2009) suggest a linear relation between the areas of the fans and the respective upstream catchment sizes. Braun (2022) reproduced such a linear relation by assuming a specific relation between catchment size and river length in the foreland. In turn, the behavior of the numerical model considered here suggests that the distance of the river to the nearest river of similar size may be the primary control rather than the absolute catchment size. While the results may finally be similar, a thorough analysis could be a starting point of a subsequent study.

In our scenario of an active mountain range and entirely passive foreland regions, it makes sense to classify rivers into two categories. Let us define carriers as rivers that receive discharge (and thus sediment flux) from the mountain range. In turn, rivers whose entire catchment is located in the passive foreland are called redistributors in the following.

The majority of the large rivers are carriers. In Fig. 1, 39 out of the 50 biggest rivers are carriers. However, it is visible that the sources of large redistributors are either close to carriers or to valleys in the mountain range. This observation suggests that large redistributors were either carriers in the past and were disconnected by avulsion events or will turn into carriers in the future.

Carriers and redistributors differ fundamentally in their properties since carriers not only receive discharge from the mountain range, but also a sediment flux. In turn, the sediment flux of redistributors arises solely from erosion in the tectonically inactive foreland. Since the foreland topography is typically not very steep, erosion rates and sediment fluxes of redistributors are rather low. As long as the base level of a redistributor (either at the ocean or at the confluence with a carrier) remains constant, the topography of its catchment decays. So redistributors are predominantly erosive with moderate rates.

Conversely, carriers must deposit sediment on average since the erosion of the redistributors would result in an ongoing decrease in topography in the foreland otherwise. Furthermore, carriers must be steeper than redistributors on average owing to their higher sediment flux. This difference is responsible for the more elongated shapes of catchments in the foreland compared to the mountain range, which is immediately recognized in Fig. 1. In the mountain range, large rivers are less steep than small rivers so that small rivers tend to flow into large rivers instead of flowing directly toward the edge of the mountain range. In contrast, the largest rivers in the foreland are carriers. Since these are rather steep, the tendency of redistributors to flow toward large carriers instead of draining into the ocean is much weaker than in the mountain range, which results in strongly elongated catchments.

Real-world topographies are typically not in a steady state. The event- or threshold-based characteristics of at least some involved processes are primary reasons for the absence of steady states. As an example, an individual big flood may contribute to landform evolution. As a second aspect, approaching a steady state may take a long time, while the tectonic and climatic conditions are typically not constant over sufficiently long time spans.

Using the *ξ*–*q* model, Yuan et al. (2019) already observed a permanent
reorganization of the rivers in a passive foreland, although
all conditions (including the uplift rate in the mountain range)
were constant. In contrast,
constant uplift rates typically result in steady-state topographies.
So there seems to be a fundamental difference in the properties of the
model between active and passive regions. This difference is not only
relevant for our understanding of landform evolution, but also for
the question of whether a record that suggests non-steady conditions
is necessarily related to changes in tectonic or climatic conditions
or whether it may be the result of self-organization.

Using the properties of carriers and redistributors described in the previous section, it can easily be shown that steady-state topographies are indeed impossible in a passive foreland for models of the type considered here. If the topography was in a steady state, the topography of all redistributors would be flat. This would imply that the catchments of all redistributors that drain directly into the ocean would be at sea level, while those that drain into a carrier would be exactly at the elevation of the point of confluence. This situation is illustrated in Fig. 5 for the drainage pattern from Fig. 1. It would immediately result in inconsistent flow directions at drainage divides. This inconsistency is obvious at the borders of the catchments of redistributors that drain directly into the ocean (dark green areas). These drainage divides would move into the adjacent higher-elevation catchments. However, even the catchments of redistributors draining into the same carrier would be inconsistent since the catchments that drain into the carrier more downstream are lower in elevation. So there would be a shift in the drainage divides, with the redistributors playing an essential part. These arguments are not restricted to the specific model considered here but only rely on the property that the catchments of redistributors become flat through time.

In the previous section, we have seen that the drainage network and thus also the topography cannot reach a steady state in the foreland even under constant conditions. Nevertheless, it is useful to investigate the properties of the rivers in a hypothetical steady state in order to find out whether they have anything in common with the rivers in a state of permanent reorganization. So let us assume that the drainage network was frozen and consider the properties of the carriers for a steady-state topography like the one shown in Fig. 5, although the flow directions are not consistent with the topography at the drainage divides.

According to Eq. (5), the channel slope in a steady state
at zero uplift (*E*=0) follows the relation

If *Q* was constant along the river (the sediment flux from the mountain range), the
concavity of the river profile would arise from the downstream increase
in discharge alone according to

with $\mathit{\theta}=\frac{m+\mathrm{1}}{n}$. Then the concavity index *θ* would
be by $\frac{\mathrm{1}}{n}$ greater than the concavity index at uniform
erosion ($\mathit{\theta}=\frac{m}{n}$). This would, however, only be
true if all tributaries were redistributors, which would not contribute
sediment in a steady state. Confluences with other
carriers lead to a downstream increase in sediment flux and thus
to a weaker concavity.

Let us consider a cross section in the *x*_{2} direction (parallel to
the mountain range). If *d* is the mean spacing of the carriers crossing
this line, each carrier has to accommodate a sediment flux of

on average, where *w* is the width of the mountain range
(so $\frac{w}{\mathrm{2}}d$ is half of the area of the mountain range)
and *U* the uplift
rate.

In order to relate the equilibrium sediment flux *Q* to the catchment size
*A*, the relation between *d* and *A* was investigated numerically.
For this purpose,
the number of carriers and their mean catchment size were measured
for each line of the
grid over the 401 snapshots. Note that the mean
catchment size cannot be estimated from the total area upstream of the
line and the respective number of carriers alone since a part of the
domain is directly drained into the ocean by redistributors.

The results shown in Fig. 6 suggest that
there are two regimes with simple scaling relations between *d* and
*A*. A linear relation *d*∝*A* is found close to the mountain range.
More important, the relation turns into a power law *d*∝*A*^{γ} with
*γ*=0.53 at greater distances, starting from about half the width of the mountain
range, so for mean catchment sizes *A* *⪆* 50 000. The question of whether
it is a fractal relation with *γ*>0.5 or *γ*=0.5
($d\propto \sqrt{A}$) is not important here.

Inserting this result into Eq. (12) yields

This relation can be reconciled with the analysis of real-world drainage networks by Prasicek et al. (2020), although that study refers to mountain catchments. It was found there that confluences of rivers of the same stream order (so of similar sizes) and capture of smaller tributaries contribute almost equally to the downstream increase in catchment size. Since the largest rivers are typically carriers, while redistributors are smaller, this means that confluence with other carriers and capture of redistributors contribute almost equally to the downstream increase in catchment size. As redistributors do not deliver sediments in a hypothetical steady state, this implies

for the downstream increases *δ**Q* and *δ**A*, which is
equivalent to Eq. (13) for *γ*≈0.5. This agreement
with the findings of Prasicek et al. (2020) suggests that Eq. (13)
is not specific to the linear version (*n*=1) of the model used here
but reflects a general property of drainage networks. So Eq. (13)
may also be helpful for developing one-dimensional versions of more
complex models.

It should, however, be kept in mind that water and sediments are routed toward
a single neighbor in the model used here. So Eq. (13) and the
following considerations are only valid for dendritic channel networks. Finding
out whether bifurcations in the network only affect the factor of
proportionality or also the exponent *γ* and thus also the concavity index
would require further research.

Inserting Eq. (13) into Eq. (10) yields

So the fundamental relation for the concavity of rivers (Eq. 11) also holds for carriers in a steady state but with a concavity index

So the concavity index of carriers in a steady state is by about
$\frac{\mathrm{1}}{\mathrm{2}n}$ higher than at uniform erosion.
For the parameter values considered here (*m*=0.5, *n*=1),
the concavity index is *θ*≈1.

Since it was shown in the previous section that steady states cannot exist
in the foreland, the question arises of whether the properties of steady-state
carriers obtained above are relevant at all. Figure 7 shows an
example of the five biggest rivers from Fig. 1 (*t*=128, solid lines)
and Fig. 5 (steady state, dashed lines).

Since the uppermost parts of the rivers are located in the mountain range, they follow straight lines with a negative slope of $\mathit{\theta}=\frac{m}{n}=\mathrm{0.5}$ in the double-logarithmic plot in equilibrium (dashed lines). The respective curves from the snapshot (solid lines) show the same overall behavior in the mountain range but with distinct local deviations in channel slope. These are disturbances propagating upstream, so they are mobile knickpoints. While these mobile knickpoints originate from changes in foreland topography, their feedback on the rivers in the foreland by means of changes in sediment flux is small.

In the foreland region,
steady-state river segments with a concavity index $\mathit{\theta}=\frac{m+\mathrm{1}}{n}=\mathrm{1.5}$
are found, as predicted for carriers that do not receive sediment flux
from their tributaries. These segments are displaced horizontally and vertically
at confluences with other carriers. As illustrated by the line $S\propto {A}^{-\mathrm{1}}$,
the overall decrease in *S* with *A* follows the prediction (Eq. 16) quite well not only for
the steady-state profiles, but also for the carriers obtained from the snapshot.
As a main difference, non-steady river segments between confluences of carriers do not
follow the steady-state relation $S\propto {A}^{-\mathrm{1.5}}$ exactly.
However, this difference apparently
has no effect on the overall decline in steepness, $S\propto {A}^{-\mathrm{1}}$; it even
seems that the rivers from the snapshot are slightly closer to this relation than the
steady-state carriers.

Actual and steady-state elevations are compared in Fig. 8. The results
were again obtained from the 401 snapshots. In order
to quantify the average behavior, mean elevations along lines in the
*x*_{2} direction are considered instead of individual data points.

If the average over the entire area is considered, which means over all carriers and redistributors
at a given *x*_{1} value, the equilibrium elevations are much lower than the actual
elevations. This finding is not surprising since almost the entire foreland
area is covered by redistributors and their catchments. These are flat in
equilibrium, wherein the catchments of redistributors draining directly into the ocean
are even at zero elevation.

In turn, the mean elevation along the actual carriers is quite close to the respective equilibrium elevation. The maximum relative deviation occurs far away from the mountain range, where the mean elevation is about 80 % of the equilibrium elevation (dark blue domain). However, absolute elevations are small here so that the deviation is also small on an absolute scale. The ratio approaches 1 toward the mountain range, where it is even between about 0.99 and 1.01 in the green, orange, and yellow regions. This result confirms the observation that the largest rivers are described reasonably well by the concept of carriers in equilibrium (Fig. 7).

The deviation in elevation between the actual carriers and the respective equilibrium profiles seems to be unimportant at first. However, it implies that carriers are less steep than expected close to the ocean and thus deliver less sediment to the ocean in total than they receive from the mountain range. This deficit is related to the drainage pattern, with a considerable part of the foreland directly draining into the ocean by redistributors. So a part of the sediment flux from the mountain range is deposited by the carriers and transported into the ocean by redistributors after the network has reorganized.

Figure 9 shows the mean sediment budget obtained from the 401 snapshots, in which all discharges (blue arrows, equivalent to catchment sizes) are expressed as percentages of the overall discharge and all sediment fluxes (orange-colored arrows) as percentages of the sediment flux from the mountain range. One-third of the total domain (including the mountain range) is drained directly into the ocean by redistributors. However, the contribution of these rivers to the total sediment delivery is only 9 % since their catchments are rather flat. So the carriers deliver more than 90 % of the sediment flux from the mountain range to the ocean.

However, this result does not imply that the carriers are routing more than 90 % of their sediment flux to the ocean and depositing less than 10 % to be cleaned up later by redistributors. Here, the sediment flux from the redistributors into the carriers also plays an important part. This flux amounts to 33 % of the influx from the mountain range. So the total sediment input to the carriers is in fact 133 %. Therefore, the 91 % delivered to the ocean is less than 70 % of the total sediment input, while about 30 % (42 out of 133) is deposited.

So the simple concept of depositing sediment by the carriers and cleaning up the deposits later by the redistributors contains an important internal component. In sum, the redistributors erode as much material as deposited by the carriers on average, equivalent to 42 % of the sediment flux from the mountain range. Only a quite small fraction of this sediment flux is delivered to the ocean (here about $\frac{\mathrm{9}}{\mathrm{42}}\approx \mathrm{21}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$), while the majority arrives in the actual carriers. A considerable part of this material is deposited further downstream and waits there to be eroded by a new generation of redistributors after the network has reorganized.

It should, however, be emphasized that these numbers not only depend
on the considered model, but also on the size of the foreland in relation
to the size of the mountain range. In order to obtain a sufficiently large
region far away from the mountain range, the foreland is in total
10 times as large as the mountain range. For a smaller foreland region,
the fluxes from the redistributors into the carriers and the respective
balances (±42 *%* of the sediment flux from the mountain range)
would be smaller.

The quite large sediment flux from redistributors into carriers also explains the high rates of sediment deposition in carriers, which will be investigated in the following section. As found in Sect. 6, longitudinal profiles of carriers are described well by equilibrium profiles on average. This means that they are able to transport almost all of the material eroded in the mountain range on average. Only in combination with the additional sediment supply from the redistributors does the total amount considerably exceed the transport capacity, which enforces rapid deposition.

Figure 10 shows the cumulative distribution of
the erosion rates of all redistributors, with negative rates
referring to deposition. These rates were evaluated in each step of the
simulation from *t*=100 to 500 (so not only at the snapshots).
So the rates are average rates over time intervals
of $\mathit{\delta}t={\mathrm{2}}^{-\mathrm{10}}\approx {\mathrm{10}}^{-\mathrm{3}}$ (typically some 100 years).

It is immediately recognized from the cumulative probabilities at zero erosion rate that the vast majority of all redistributors is indeed eroding as discussed earlier. The fraction of eroding redistributor sites is greater than 90 % everywhere and even greater than 97 % close to the mountain range (yellow curve).

The rates are overall rather low compared to the mountain range ($E=U=\mathrm{1}$). In the four domains close to the mountain range (yellow to green curves), the net mean erosion rate (deposition contributes negatively) is about 0.1, which is about 10 % of the erosion rate in the mountain range. Further away from the mountain range (blue curves), the rates decrease rapidly. The average erosion rate over all redistributors is 0.042. Taking into account that the foreland area is 10 times as large as the mountain range and covered almost entirely by redistributor sites, this result could also be obtained directly from the balance of the redistributors shown in Fig. 9 (total erosion amounting to 42 % of the sediment flux from the mountain range).

The respective rates for the carriers are shown in Fig. 11. As carriers predominantly deposit sediments, rates of deposition are shown here, with negative rates describing erosion. The rates are overall high compared to the redistributors and also compared to the mountain range. The overall mean rate is 4.8, which is almost 5 times higher than the erosion rate in the mountain range. Since the carrier sites cover only a small part of the area, this result cannot be derived directly from the overall balance (Fig. 9).

In contrast to the redistributors, the lowest mean net rates occur close to the mountain range. The yellow domain is the only region where the mean net rate of deposition (≈0.7) is lower than the erosion rate in the mountain range. As shown in Fig. 4, the lateral mobility is also small in this region, at least for the largest river.

The rates increase rapidly with increasing distance from the mountain range and reach a mean net rate greater than 10 in the turquoise domain. This domain extends to a distance of 512 pixels from the mountain range, which is equal to the full width of the mountain range. Both deposition and erosion rates are quite high here. Almost 50 % of the carriers in this domain deposit at rates higher than the erosion rate in the mountain range, and 30 % of them erode at such rates. The rate of deposition is even higher than 60 times the erosion rate in the mountain range at more than 5 % of the carrier sites here. With regard to these high rates of deposition, we must keep in mind that the model does not explicitly include floodplains. Since sediments are deposited pixel by pixel, wider floodplains can only be filled by rapid avulsions. This may result in an overestimation of rates of deposition.

Further away from the mountain range, the rates decrease (blue domains). However, this decrease is much slower than for the redistributors, and the mean rates of deposition stay clearly above the erosion rate in the mountain range. The occurrence of a maximum in deposition rate at some distance from the mountain range is related to the drainage pattern. Close to the mountain range, almost all rivers are carriers. So a rather small sediment flux from the redistributors is distributed among a large number of carriers and thus has little effect compared to the sediment flux from the mountain range.

The strong autocyclicity in the sediment deposition by carriers raises the
question of whether it may even obscure effects of climatic variations.
Recently, Yuan et al. (2022) investigated the effect of periodic oscillations
in precipitation with a similar model and found a clear
link between modeled climate oscillations and sediment signals.
While a variation of some tens of meters in the mean elevation
of the riverbed was found at the edge of the mountain range,
the respective rates of aggradation are rather low.
Assuming a relative variation in precipitation by ±50 *%*, the maximum
rate of aggradation is about
0.25 mm yr^{−1} at an uplift rate of 1 mm yr^{−1} (independent
of the period, estimated from their Fig. 2). This rate is already lower
than the mean rate found here for the very proximal region of the alluvial
fans and more than a decade lower than the mean rates in the more distal
regions. So autocyclicity will probably dominate the record at short timescales
and small spatial scales except for the very proximal region of the alluvial
fans. However, finding out which scales are needed to overcome the dominance of
autocyclic aggradation requires further investigations.

While the interplay of carriers and redistributors is the key to understanding the long-term reorganization of the drainage network, the respective timescale cannot be determined directly from the rates of erosion and deposition. In this section, a simple concept is used for quantifying the reorganization of the drainage network and its potential fingerprint in the deposits. Knowing that carriers predominantly deposit sediment and redistributors predominantly erode, the time span since the previous flooding by a carrier is investigated. Figure 12 shows a map of this time span for the snapshot from Fig. 1. The spatial pattern is irregular. While dark blue areas depict the mountain range and the actual carriers, there are more or less continuous blueish areas depicting large regions that were recently flooded by carriers. In turn, there are also yellow areas indicating that some regions were not flooded by carriers for more than 3 time units. Some of these yellow areas are intersected by darker lines, which means that individual carriers crossed the area later without affecting the area as a whole.

The cumulative statistical distribution of the time span since the previous flooding by a carrier is shown in Fig. 13. The distributions were derived from the 401 snapshots. If flooding by carriers was a random process wherein the probability is independent of the time span since the previous event, the time spans would follow an exponential distribution (a straight line in the plot). The results suggest an exponential distribution at large time spans.

The time spans since the previous flooding by a carrier increase from the mountain range toward the ocean. While the mean time span is about 0.2 in the yellow domain, it increases to about 1.25 in the two outermost domains. So we can expect older deposits at the surface further away from the mountain range.

As it was found in Sect. 3 that the frequency of avulsions of large
rivers depends strongly on the time increment *δ**t* used in the simulation,
the timescales of network reorganization also depend on *δ**t*. However,
this dependence is weaker than expected from the frequency of avulsions.
While reducing *δ**t* from 2^{−10} to 2^{−13} increases
the frequency of avulsions of large rivers by a factor of almost 3, the mean time
spans since the previous flooding by a carrier decrease only by about 20 % to 30 %.
So a large fraction of the observed avulsions seem to be related to either
back-and-forth oscillations in flow direction owing to the single-neighbor (D8) flow-routing
scheme or to filling floodplains pixel by pixel. In particular,
the increase in the time span since the previous flooding by a carrier
toward the ocean appears to be robust.

The increase in the time span since the previous flooding by a carrier might be interpreted as a decreasing rate of river avulsion toward the ocean, similarly to the decreasing rates of deposition. However, we also have to take into account the spacing of the carriers here, which also increases toward the ocean as recognized Fig. 6. The increase in spacing is even stronger than the increase in the time span since the previous flooding. So the individual rivers do not become less active toward the ocean, and the longer time spans arise from a smaller number of carriers sweeping over the area.

The distributions shown in Fig. 13 deviate from exponential distributions at short time spans. The rapid decline of the distribution in the yellow domain indicates a clustering of events in the sense that there is an increased probability that the river will return to a location where it was recently. A random walk is the simplest process with this property, which here could be diverting the river randomly in the regime of the alluvial fans close to the mountain range. Far away from the mountain range, the behavior is opposite. Here, the probability of flooding by a carrier increases if the previous flooding was long ago. This effect is presumably related to the topography. Since carriers predominantly deposit sediment and redistributors erode, large channels that have not been carriers for a long time are rather flat and are thus favored candidates for becoming carriers in the future.

Using cosmogenic Ne isotopes, Sinclair et al. (2019) obtained a strong variation in the residence times of pebbles about 1000 km downstream of the source region located in the Rocky Mountains. It was found that some of the investigated clasts had even spent some million years in the river. In the same study, a conceptual model for the recycling of pebbles by tributaries was proposed. In principle, this concept already implicitly contains the idea of carriers and redistributors.

Concerning the order of magnitude, the timescale of network reorganization obtained in this study aligns well with the transit times obtained by Sinclair et al. (2019). However, we should keep in mind that the time span since the previous flooding by a carrier is not equivalent to transit times of sediment clasts. The generic form of the shared stream-power model (Eq. 1) predicts only a net rate but not the exchange of particles between the river and an alluvial cover. Therefore, tracking sediment clasts in OpenLEM would not only require some technical effort, but also an extension of the shared stream-power model itself.

Carretier et al. (2016) already added a component for tracking clasts to the model CIDRE, which explicitly defines erosion and deposition rates. Simulating the distribution of clasts in a smaller region (equivalent to the yellow, orange, and green domains in this study), a variation of transit times over several decades was found (Carretier et al., 2020). The obtained distribution may even be heavy-tailed with a power-law tail, consistent with results obtained from a canyon in the central Andes (Carretier et al., 2019). In turn, the exponential distribution of the time span since the previous flooding by a carrier suggests that the reorganization of the drainage network alone would not produce heavy-tailed distributions. So finding out whether transit-time distributions of sediment clasts are indeed heavy-tailed and, if so, explaining their origin requires further research.

So far, a transport-limited model has been used for the foreland region. However, we have seen in the previous section that a considerable part of the area may be covered by deposits older than some 100 000 years; the question arises of whether eroding such deposits is described well by a transport-limited model.

In order to investigate the effect of
a finite erodibility *K*_{d} in the erosive regime, an
extension of the numerical scheme proposed by Hergarten (2020)
was developed and implemented in OpenLEM. This extension
switches between the shared stream-power model with a finite erodibility
*K*_{d} and the transport-limited end-member
(*K*_{d}→∞) at each node. By integrating the decision
into the scheme, the fully implicit character of the scheme can be
almost completely preserved. Only the base level of each node
(the elevation of the flow target) has to be adopted from the beginning
of the respective time step, while the actual values of all other properties
can be included in the decision.

As an extreme scenario, the same parameters as in the mountain range were assumed for the erosive regime in the foreland (${K}_{\mathrm{d}}={K}_{\mathrm{t}}=\mathrm{1}$). This would be an instantaneous consolidation of all deposits to a rock with the same properties as the bedrock in the mountain range. Although unrealistic, this extreme scenario is useful for investigating the effect of not fully transport-limited conditions in the erosive regime.

It was already recognized in Sect. 7 that the sediment fluxes from
the redistributors into the carriers are very important for the high rates
of sediment deposition in the carriers. The mean erosion rates of the redistributors
indeed decrease strongly; a more than fivefold decrease was found except
for the two regions closest to the mountain range
(yellow and orange). This may be surprising
at first because the effective erodibility (Eq. 6) is reduced only by a
factor of 2 (from 1 to 0.5). However, it was already shown in the context of knickpoint
migration that disturbances propagate upstream at a velocity defined by *K*_{d},
which explains the big difference between *K*_{d}→∞ and
*K*_{d}=1 in the transient behavior.

The reduction in mean erosion rate equivalently reduces the sediment fluxes from the redistributors into the carriers and into the ocean. As shown in Fig. 14, this reduction strongly affects the sediment balance of the carriers. Now the total sediment input is only 107 % instead of 133 % (Fig. 9), while 98 % is delivered to the ocean instead of 91 %. So 92 % of the total sediment input is delivered to the ocean and only 8 % is deposited. The respective ratio was about 30 % in the transport-limited model (Fig. 9). So moving from the transport-limited model to the shared stream-power model with ${K}_{\mathrm{d}}={K}_{\mathrm{t}}=\mathrm{1}$ also considerably reduces the rates of sediment deposition in the carriers, although the decrease is not as strong as in the erosion rates of the redistributors. In agreement with the results of the sediment balance, the decrease is by a factor of 3.5 to 4.8, except for the first and last domain (yellow and dark blue), for which it is less than 3.

The increase in the rates of erosion and deposition also slows down the dynamics of network reorganization. While the distributions of the time since the previous flooding by a carrier were found to be qualitatively similar to the distribution shown in Fig. 13, the timescale is stretched. The mean time since the previous flooding by a carrier increases by a factor of 4.6 to 5.8 for the individual domains. These factors follow the decrease in erosion rates in the redistributors rather than the decrease in deposition rates in the carriers. This finding emphasizes the relevance of the erosion by redistributors and the resulting sediment flux to the carriers for the dynamics of the foreland rivers.

This paper is intended to define some kind of reference scenario for fluvial landform evolution in a tectonically inactive foreland of a mountain range. Additionally, it can also be seen as a starting point for further studies.

The considered scenario combines one of the simplest models of large-scale fluvial erosion and sediment transport – the shared stream-power model – with a simple geometry consisting of a mountain range and an inactive foreland. In order to understand the behavior of the model, the foreland rivers were subdivided into two classes – carriers and redistributors. Carriers originate in the mountain range and are thus responsible for the large-scale sediment transport to the ocean. In turn, redistributors are rivers whose entire catchment is located in the foreland.

Using the concept of carriers and redistributors,
it was shown that a steady-state topography in the strict sense is impossible
in the foreland even under constant conditions. Although the topography becomes
more or less constant on average over long times, the drainage network in the
foreland permanently reorganizes. On the other hand, longitudinal profiles
of carriers are described well by a hypothetical steady state on average,
with the sediment flux from the mountain range just routed to the ocean.
The concavity index of carriers is typically greater than the concavity index
$\mathit{\theta}=\frac{m}{n}$ of rivers in a mountain range at uniform erosion but smaller than
the value $\mathit{\theta}=\frac{m+\mathrm{1}}{n}$ expected for carriers without sediment supply
by tributaries. By analyzing the topology of the drainage network, it was found that
*θ* is in the middle between these two values, so *θ*≈1 for
the linear version of the shared stream-power model (*n*=1).

It was found that redistributors predominantly erode at rates lower than the erosion rate in the mountain range. In turn, carriers predominantly deposit sediments, with the rates typically much higher than the erosion rate in the mountain range. As a major result, the sediment flux from the redistributors into the carriers plays a central part for the deposition of sediments and for the reorganization of the drainage network. While the erosion rates of the redistributors are rather low, the respective areas are large, generating considerable sediment input in total. As a consequence, the assumptions on the erosion in the foreland are more important for the dynamics of the rivers than it may seem at first.

While these results might be fundamental, there are two major limitations, which are related to each other. First, there is the finite length of the time increment used in the numerical simulation. Concerning the accuracy, much smaller time increments than used in this study would be desirable. This would, however, not only increase the numerical effort, but would also shift the timescale into a range in which individual flood events become relevant in reality. As a second point, floodplains are not explicitly taken into account. Valleys are filled with sediments pixel by pixel in combination with a high frequency of avulsions. Finding an approach to extend deposition laterally over floodplains would probably also help to solve the issue with the time increment.

There are several further aspects for which subsequent studies should go deeper. This also includes the consideration of transient states and the comparison to real-world topographies. In addition, rates of sediment deposition were only investigated at a given timescale. Given the numerical efficiency of the model OpenLEM used here, the scaling properties could be investigated over a range from some 100 years up to millions of years. As a next step, tracking ages of deposits in different depths would be interesting and would allow for a validation by real-world data, although this would be technically more challenging.

In addition, the nonlinear version of the shared stream-power model (so with exponents
*n* > 1) should also be investigated in subsequent studies.
From a theoretical point of view, the still existing uncertainty concerning the
exponent *n* in the stream-power formulation is still a challenge. The theoretical
considerations of the concavity index of carriers suggest that
the exponent *n* has an effect in the foreland even under spatially uniform
conditions, in contrast to active mountain ranges. So a more thorough investigation
of the influence of the exponent in combination with real-world river profiles may
considerably contribute to our knowledge of the value of this exponent.

All code is available in a Zenodo repository at https://doi.org/10.5281/zenodo.6770230 (Hergarten, 2022b). This repository also contains a part of the simulated data. The full set of simulated data will be stored for a limited time and is available upon request. Users who are interested in using the landform evolution model OpenLEM in their own research are advised to download the most recent version from http://hergarten.at/openlem (Hergarten, 2022a). The author is happy to assist interested readers in reproducing the results and performing subsequent research.

The author has declared that there are no competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The author would like to thank Jean Braun and two anonymous reviewers for their constructive comments and Sagy Cohen for the editorial handling.

This open-access publication was funded by the University of Freiburg.

This paper was edited by Sagy Cohen and reviewed by Jean Braun and two anonymous referees.

Adams, B. A., Whipple, K. X., Forte, A. M., Heimsath, M., and Hodges, K. V.: Climate controls on erosion in tectonically active landscapes, Sci. Adv., 6, eaaz3166, https://doi.org/10.1126/sciadv.aaz3166, 2020. a

Armitage, J. J., Duller, R. A., Whittaker, A. C., and Allen, P. A.: Transformation of tectonic and climatic signals from source to sedimentary archive, Nat. Geosci., 4, 231–235, https://doi.org/10.1038/ngeo1087, 2011. a

Armitage, J. J., Dunkley Jones, T., Duller, R. A., Whittaker, A. C., and Allen, P. A.: Temporal buffering of climate-driven sediment flux cycles by transient catchment response, Earth Planet. Sc. Lett., 369–370, 200–210, https://doi.org/10.1016/j.epsl.2013.03.020, 2013. a

Blair, T. C. and McPherson, J. G.: Processes and forms of alluvial fans, in: Geomorphology of Desert Environments, edited by: Parsons, A. J. and Abrahams, A. D., Springer, Dordrecht, 413–467, https://doi.org/10.1007/978-1-4020-5719-9_14, 2009. a

Blom, A., Viparelli, E., and Chavarrias, V.: The graded alluvial river: Profile concavity and downstream fining, Geophys. Res. Lett., 43, 6285–6293, https://doi.org/10.1002/2016GL068898, 2016. a

Blom, A., Arkesteijn, L., Chavarrias, V., and Viparelli, E.: The equilibrium alluvial river under variable flow and its channel-forming discharge, J. Geophys. Res.-Earth, 122, 1924–1948, https://doi.org/10.1002/2017JF004213, 2017. a

Bolla Pittaluga, M., Luchi, R., and Seminara, G.: On the equilibrium profile of river beds, J. Geophys. Res. -Earth, 119, 317–332, https://doi.org/10.1002/2013JF002806, 2014. a

Braun, J.: Comparing the transport-limited and ξ–q models for sediment transport, Earth Surf. Dynam., 10, 301–327, https://doi.org/10.5194/esurf-10-301-2022, 2022. a, b

Braun, J. and Willett, S. D.: A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution, Geomorphology, 180–181, 170–179, https://doi.org/10.1016/j.geomorph.2012.10.008, 2013. a

Carretier, S., Martinod, P., Reich, M., and Godderis, Y.: Modelling sediment clasts transport during landscape evolution, Earth Surf. Dynam., 4, 237–251, https://doi.org/10.5194/esurf-4-237-2016, 2016. a, b

Carretier, S., Regard, V., Leanni, L., and Farias, M.: Long-term dispersion of
river gravel in a canyon in the Atacama Desert, Central Andes,
deduced from their ^{10}Be concentrations, Sci. Rep.-UK, 9, 17763,
https://doi.org/10.1038/s41598-019-53806-x, 2019. a

Carretier, S., Guerit, L., Harries, R., Regard, V., Maffre, P., and Bonnet, S.: The distribution of sediment residence times at the foot of mountains and its implications for proxies recorded in sedimentary basins, Earth Planet. Sc. Lett., 546, 116448, https://doi.org/10.1016/j.epsl.2020.116448, 2020. a, b

Coulthard, T. J.: Landscape evolution models: a software review, Hydrol. Process., 15, 165–173, https://doi.org/10.1002/hyp.426, 2001. a

Davy, P. and Lague, D.: Fluvial erosion/transport equation of landscape evolution models revisited, J. Geophys. Res.-Earth, 114, F03007, https://doi.org/10.1029/2008JF001146, 2009. a, b

Guerit, L., Yuan, X. P., Carretier, S., Bonnet, S., Rohais, S., Braun, J., and Rouby, D.: Fluvial landscape evolution controlled by the sediment deposition coefficient: Estimation from experimental and natural landscapes, Geology, 47, 853–856, https://doi.org/10.1130/G46356.1, 2019. a

Hack, J. T.: Studies of longitudinal profiles in Virginia and Maryland, no. 294-B in US Geol. Survey Prof. Papers, US Government Printing Office, Washington D.C., https://doi.org/10.3133/pp294B, 1957. a

Harel, M.-A., Mudd, S. M., and Attal, M.: Global analysis of the stream power
law parameters based on worldwide ^{10}Be denudation rates,
Geomorphology, 268, 184–196, https://doi.org/10.1016/j.geomorph.2016.05.035, 2016. a

Hergarten, S.: Transport-limited fluvial erosion – simple formulation and efficient numerical treatment, Earth Surf. Dynam., 8, 841–854, https://doi.org/10.5194/esurf-8-841-2020, 2020. a, b, c, d, e, f

Hergarten, S.: The influence of sediment transport on stationary and mobile knickpoints in river profiles, J. Geophys. Res.-Earth, 126, e2021JF006218, https://doi.org/10.1029/2021JF006218, 2021. a, b

Hergarten, S.: OpenLEM, http://hergarten.at/openlem, last access: 9 June 2022a. a, b

Hergarten, S.: Rivers in a tectonically inactive foreland, Zenodo [code] and [data set], https://doi.org/10.5281/zenodo.6770230, 2022b. a

Hergarten, S. and Neugebauer, H. J.: Self-organized critical drainage networks, Phys. Rev. Lett., 86, 2689–2692, https://doi.org/10.1103/PhysRevLett.86.2689, 2001. a

Hergarten, S. and Robl, J.: The linear feedback precipitation model (LFPM 1.0) – a simple and efficient model for orographic precipitation in the context of landform evolution modeling, Geosci. Model Dev., 15, 2063–2084, https://doi.org/10.5194/gmd-15-2063-2022, 2022. a

Hilley, G. E., Porder, S., Aron, F., Baden, C. W., Johnstone, S. A., Liu, F., Sare, R., Steelquist, A., and Young, H. H.: Earth’s topographic relief potentially limited by an upper bound on channel steepness, Nat. Geosci., 12, 828–832, https://doi.org/10.1038/s41561-019-0442-3, 2019. a

Howard, A. D.: A detachment-limited model for drainage basin evolution, Water Resour. Res., 30, 2261–2285, https://doi.org/10.1029/94WR00757, 1994. a, b

Kooi, H. and Beaumont, C.: Escarpment evolution on high-elevation rifted margins: insights derived from a surface process model that combines diffusion, advection and reaction, J. Geophys. Res., 99, 12191–12209, https://doi.org/10.1029/94JB00047, 1994. a, b

Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Proc. Land., 39, 38–61, https://doi.org/10.1002/esp.3462, 2014. a, b

Malatesta, L. C., Prancevic, J., and Avouac, J.-P.: Autogenic entrenchment patterns and terraces due to coupling with lateral erosion in incising alluvial channels, J. Geophys. Res.-Earth, 121, 335–355, https://doi.org/10.1002/2015JF003797, 2017. a

Mouchené, M., van der Beek, P., Carretier, S., and Mouthereau, F.: Autogenic versus allogenic controls on the evolution of a coupled fluvial megafan–mountainous catchment system: numerical modelling and comparison with the Lannemezan megafan system (northern Pyrenees, France), Earth Surf. Dynam., 5, 125–143, https://doi.org/10.5194/esurf-5-125-2017, 2017. a

O'Callaghan, J. F. and Mark, D. M.: The extraction of drainage networks from digital elevation data, Computer Vision, Graphics, and Image Processing, 28, 323–344, https://doi.org/10.1016/S0734-189X(84)80011-0, 1984. a

Prasicek, G., Hergarten, S., Deal, E., Herman, F., and Robl, J.: A glacial buzzsaw effect generated by efficient erosion of temperate glaciers in a steady state model, Earth Planet. Sc. Lett., 543, 116350, https://doi.org/10.1016/j.epsl.2020.116350, 2020. a, b

Robl, J., Hergarten, S., and Prasicek, G.: The topographic state of fluvially conditioned mountain ranges, Earth Sci. Rev., 168, 290–317, https://doi.org/10.1016/j.earscirev.2017.03.007, 2017. a, b

Romans, B. W., Castelltort, S., Covault, J. A., Fildani, A., and Walsh, J. P.: Environmental signal propagation in sedimentary systems across timescales, Earth Sci. Rev., 153, 7–29, https://doi.org/10.1016/j.earscirev.2015.07.012, 2016. a

Sadler, P. M.: Sediment accumulation rates and the completeness of stratigraphic sections, J. Geol., 89, 569–584, 1981. a

Sinclair, H. D., Stuart, F. M., Mudd, S., McCann, L., and Tao, Z.: Detrital
cosmogenic ^{21}Ne records decoupling of source-to-sink signals by
sediment storage and recycling in Miocene to present rivers of the Great
Plains, Nebraska, USA, Geology, 47, 3–6, https://doi.org/10.1130/G45391.1, 2019. a, b

Tofelde, S., Bernhardt, A., Guerit, L., and Romans, B. W.: Times associated with source-to-sink propagation of environmental signals during landscape transience, Front. Earth Sci., 9, 628315, https://doi.org/10.3389/feart.2021.628315, 2021. a

van der Beek, P.: Modelling landscape evolution, in: Environmental Modelling: Finding Simplicity in Complexity, edited by: Wainwright, J. and Mulligan, M., Wiley-Blackwell, Chichester, 2 edn., 309–331, ISBN 978-0-470-74911-1, 2013. a

Whipple, K. X. and Tucker, G. E.: Implications of sediment-flux-dependent river incision models for landscape evolution, J. Geophys. Res., 107, 2039, https://doi.org/10.1029/2000JB000044, 2002. a

Whipple, K. X., DiBiase, R. A., and Crosby, B. T.: Bedrock rivers, in: Fluvial Geomorphology, edited by: Shroder, J. and Wohl, E., vol. 9 of Treatise on Geomorphology, Academic Press, San Diego, CA, 550–573, https://doi.org/10.1016/B978-0-12-374739-6.00254-2, 2013. a

Wickert, A. D. and Schildgen, T. F.: Long-profile evolution of transport-limited gravel-bed rivers, Earth Surf. Dynam., 7, 17–43, https://doi.org/10.5194/esurf-7-17-2019, 2019. a

Willgoose, G.: Mathematical modeling of whole landscape evolution, Annu. Rev. Earth Pl. Sc., 33, 443–459, https://doi.org/10.1146/annurev.earth.33.092203.122610, 2005. a

Yuan, X. P., Braun, J., Guerit, L., Rouby, D., and Cordonnier, G.: A new efficient method to solve the stream power law model taking into account sediment deposition, J. Geophys. Res.-Earth, 124, 1346–1365, https://doi.org/10.1029/2018JF004867, 2019. a, b, c, d, e

Yuan, X. P., Guerit, L., Braun, J., Rouby, D., , and Shobe, C. M.: Thickness of fluvial deposits records climate oscillations, J. Geophys. Res.-Sol. Ea., 127, e2021JB023510, https://doi.org/10.1029/2021JB023510, 2022. a, b, c, d, e

- Abstract
- Introduction
- Model setup
- First results
- The concept of carriers and redistributors
- The impossibility of steady-state solutions
- Properties of carriers in a steady state
- Source-to-sink considerations
- Rates of erosion and deposition
- The timescale of network reorganization
- The effect of consolidation
- Conclusions and outlook
- Code and data availability
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Model setup
- First results
- The concept of carriers and redistributors
- The impossibility of steady-state solutions
- Properties of carriers in a steady state
- Source-to-sink considerations
- Rates of erosion and deposition
- The timescale of network reorganization
- The effect of consolidation
- Conclusions and outlook
- Code and data availability
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References