the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Transportlimited fluvial erosion – simple formulation and efficient numerical treatment
Stefan Hergarten
Most of the recent studies modeling fluvial erosion in the context of tectonic geomorphology focus on the detachmentlimited regime. One reason for this simplification is the simple relationship of the constitutive law used here – often called streampower law – to empirical results on longitudinal river profiles. Another no less important reason lies in the numerical effort that is much higher for transportlimited models than for detachmentlimited models. This study proposes a formulation of transportlimited erosion where the relationship to empirical results on river profiles is almost as simple as it is for the streampower law. As a central point, a direct solver for the fully implicit scheme is presented. This solver requires no iteration for the linear version of the model, allows for arbitrarily large time increments, and is almost as efficient as the established implicit solver for detachmentlimited erosion. The numerical scheme can also be applied to linear hybrid models that cover the range between the two endmembers of detachmentlimited and transportlimited erosion.
Rivers play a major if not dominant part in largescale landform evolution. If horizontal displacement of the crust is not taken into account, models describing the evolution of a topography $H({x}_{\mathrm{1}},{x}_{\mathrm{2}},t)$ are typically written in the form
where U and E are the rates of uplift of crustal material relative to a given datum and of erosion, respectively.
Two endmembers – detachmentlimited and transportlimited erosion – are widely considered in the context of fluvial landform evolution. The term “detachmentlimited erosion” was presumably coined by Howard (1994). The idea behind this concept is that all particles entrained by the river are immediately removed from the system. The erosion rate E can be considered as a function of local properties at each point. In the simplest approach, these are catchment size and channel slope (slope in the direction of steepest descent), while all other influences are subsumed in a lumped parameter often called erodibility.
In all scenarios other than the detachmentlimited case, a sediment balance must be considered. If no material is directly removed, the erosion rate is
where q is the sediment flux per unit width (volume per time and cross section length) and div is the 2D divergence operator. It is usually assumed that q follows the direction of the channel slope, so only its absolute value q varies between different models.
The concept of transportlimited erosion assumes that the rate of bed erosion is limited by the ability of the flow to transport the eroded material, rather than by the availability of potentially mobile sediment. The implementation of this concept in fluvial landform evolution models presumably dates back to Willgoose et al. (1991b). Transportlimited models directly define the sediment flux per unit width q instead of the erosion rate E at each point as a function of local properties such as catchment size and channel slope.
Mathematically, both concepts differ fundamentally. Equation (1) only involves derivatives of first order with regard to time and with regard to the spatial coordinates (arising from the channel slope) in the detachmentlimited scenario. So it is a hyperbolic differential equation of the advection type. Propagation of information in one direction only – upstream here – is a characteristic property of this type. Anything that happens at a given point and a given time only affects the region upstream of this point in the future. In contrast, Eq. (1) contains spatial derivatives of second order in the transportlimited regime since q inside the divergence operator depends on the channel slope. Equation (1) combined with Eq. (2) is a parabolic differential equation of the diffusion type then, where information propagates in both upstream and downstream direction.
Several comprehensive numerical models of fluvial landform evolution have been developed since the 1990s. All models reviewed by Coulthard (2001), Willgoose (2005), and van der Beek (2013) involve a sediment balance. In recent years, however, using the detachmentlimited model has become a popular choice, although the idea that all particles are immediately excavated is limited and has been questioned (e.g., Turowski, 2012). All types of bedload transport are obviously not captured by this concept. Nevertheless, even some recent studies using models that are able to simulate sediment transport focus on the detachmentlimited case (e.g., Duvall and Tucker, 2015; Theodoratos et al., 2018; Eizenhöfer et al., 2019).
At least three aspects make the detachmentlimited approach appealing. First, the relationship to empirical studies of longitudinal channel profiles is particularly simple here. Hack (1957) observed a powerlaw relationship between channel slope S and upstream catchment size A in several rivers. This relationship is nowadays often called Flint's law (Flint, 1974) and written in the form
where θ is the concavity index and k_{s} the steepness index. Assuming that Eq. (3) is the fingerprint of a spatially constant erosion rate under uniform conditions, it can be assumed that
where f is an arbitrary function. A powerlaw function
where the parameter K is denoted erodibility, is a common choice in this context. The fluvial erosion rate is often written in the form
with m=θn. Equation (6) is often called streampower law or streampower incision model since it can be interpreted in terms of energy dissipation of the water per channel bed area if an empirical relationship between channel width and catchment size is used (e.g., Whipple and Tucker, 1999).
The concavity index $\mathit{\theta}=\frac{m}{n}$ appears to be well constrained, so most modeling studies either use the value θ=0.5 originally found by Hack (1957) or a slightly lower reference value θ=0.45 (e.g., Whipple et al., 2013; Lague, 2014). In turn, the value of the exponent n is less well constrained since it cannot be determined from the shape of equilibrium profiles under uniform conditions. The model is linear with regard to H (if the flow pattern is given) for n=1, which simplifies both theoretical considerations and the numerical implementation. Thus, the remaining uncertainty in the effective value of n often serves as a reason for choosing n=1, although the results compiled by Lague (2014) rather suggest n>1, for example. If θ is well constrained and n=1 is accepted as a convenient choice, the erodibility K remains as the only parameter. It is a lumped parameter subsuming all influences on erosion other than channel slope and catchment size. So it is not only a property of the rock, but also depends on climate in a nontrivial way (e.g., Ferrier et al., 2013; Harel et al., 2016). However, it just defines how steep rivers will become at a given uplift rate, so reasonable values can be found, e.g., by analyzing river profiles in situations where estimates of the uplift rate are available.
Constitutive laws based on powerlaw relations, however, have not only been employed in detachmentlimited models. Even the earliest numerical model of transportlimited erosion (Willgoose et al., 1991b) used a power law for the sediment flux density based on physical relations for the shear stress at the bed. The empirical results on real rivers represented by Eq. (3) were also used to constrain the parameter values before the detachmentlimited concept became popular (Willgoose et al., 1991a). However, the transportlimited approach never reached the simplicity of the detachmentlimited approach with regard to the small number of parameters and their quite direct relation to the properties of real river profiles.
The simplicity of the differential equation itself serves as a second argument in favor of the detachmentlimited approach. In the linear case (n=1), Eq. (1) combined with Eq. (6) can be solved analytically for any given uplift pattern and history. The term KA^{θ} defines the velocity of advection then, so disturbances propagate upstream at this velocity. The treatment can be simplified by the χ transform introduced by Perron and Royden (2013). It transforms the upstream coordinate x to a new coordinate,
where A_{0} is an arbitrary reference catchment size and the integration starts from an arbitrary reference point. This transformation eliminates the inherent curvature of river profiles arising from the upstream decrease in catchment size, so equilibrium profiles under spatially uniform conditions turn into straight lines. The solutions of this equation and their potential for unraveling the uplift and erosion history were investigated by Royden and Perron (2013), and a formal inversion procedure for the linear case (n=1) was presented by Goren et al. (2014). So the detachmentlimited model can be reconciled with real river profiles not only under steadystate conditions, but also in the context of temporal changes.
As a third point, detachmentlimited erosion can be implemented in numerical models more efficiently than transportlimited erosion. Here, even a fully implicit scheme that allows for arbitrary time increments with linear time complexity, also known as O(N), is available. This means that the computing effort increases only linearly with the total number of nodes N. The scheme was introduced in the context of fluvial erosion by Hergarten and Neugebauer (2001), described in detail for n=1 and n=2 by Hergarten (2002), and made popular by Braun and Willett (2013).
So far there is no comparable implementation for transportlimited erosion. As mentioned above, transportlimited erosion corresponds to a diffusiontype equation. The challenge is that the diffusivity depends on the catchment size and thus varies over several orders of magnitude. Multigrid methods (e.g., Hackbusch, 1985) are still the only schemes for the diffusion equation in more than one dimension with linear time complexity. However, convergence breaks down if the diffusivity varies by some orders of magnitude, so multigrid methods have not been applied in the context of fluvial erosion. So far none of the existing landform evolution models treat the transportlimited case with a fully implicit scheme that allows for arbitrarily large time increments.
The advantage of the detachmentlimited model concerning the numerical complexity persists if explicit schemes are used here, too. The main reason for using explicit schemes for detachmentlimited erosion is the artificial smoothing of knickpoints by the implicit discretization, while explicit schemes that preserve the shape of knickpoints better are available. A comparison was given by Campforts et al. (2017). As already pointed out by Howard (1994), explicit schemes for the transportlimited case typically require time steps 3 to 4 orders of magnitude shorter than for the detachmentlimited case.
Howard (1994) already developed an approximation that makes the explicit scheme for the transport term numerically more stable. Kooi and Beaumont (1994) proposed an approach that increases stability and also allows for a physical interpretation, often called an undercapacity model or – in a more general context – linear decline model (Whipple and Tucker, 2002). It defines an equilibrium flux per unit width q_{e} from local properties (channel slope, catchment size, etc.) and assumes that the erosion rate is
The parameter l defines a length scale and can be seen as inertia of sediment detachment and deposition against changes in fluvial conditions.
An alternative physical interpretation of the linear decline model was developed by Davy and Lague (2009). The detachmentlimited model (Eq. 6) was extended by a sediment deposition term proportional to the actual sediment flux. As a main point, Davy and Lague (2009) found an expression for the rate of deposition that keeps equilibrium river profiles consistent with Eq. (3), which is not the case for the original undercapacity model (Whipple and Tucker, 2002).
Yuan et al. (2019) implemented an implicit numerical scheme for this model based on a Gauss–Seidel iteration in the upstream direction. The rate of convergence was found to be independent of the size of the grid, so the scheme is indeed of linear time complexity. However, the convergence slows down strongly for faster deposition, i.e., when approaching the transportlimited endmember. The scheme of Yuan et al. (2019) is therefore presumably the most efficient implementation of sediment transport in largescale fluvial erosion models but achieves its full power only if we do not come too close to the transportlimited endmember.
In the following section, a formulation of transportlimited erosion is proposed that can be directly reconciled with the concept of the erodibility. Then, Sect. 3 presents a fully implicit, direct scheme for solving the equation numerically. After presenting a numerical example in Sect. 4, Sect. 6 provides a discussion of several versions of the linear decline model and an extension of the numerical scheme for this class of models.
Let us start from the interpretation of Hack's empirical relation (Eq. 3) as the fingerprint of uniform erosion under spatially constant conditions, regardless of the mechanism of erosion. This implies that the erosion rate is a function of the steepness index (Eq. 4). Then the sediment flux Q (volume per time, not per unit width) through any cross section of a river is
where the integral extends over the upstream catchment. For uniform erosion, the integral reduces to the product of the erosion rate and the catchment size,
If a powerlaw function (Eq. 6) is used as it is in the detachmentlimited model, the sediment flux becomes
In contrast to the more common formalism based on the flux per unit width q (Eq. 2), these relations use the total sediment flux Q (volume per time) passing the entire cross section of a channel segment. This total flux cannot be inserted formally into the divergence operator in Eq. (2) to form a continuous differential equation. Practically, however, this is not a problem for a discrete channel network. If any pixel of the considered topography has a unique drainage direction towards a single neighbor and sediment transport follows flow direction, the respective discrete version of the divergence operator at the node i is
where Q_{i} is the flux from the node i to its flow target. The sum extends over all neighbors that deliver their discharge und thus their sediment flux to the node i, called donors in the following. Finally, s_{i} is the area of the considered node, i.e., the pixel size for a regular mesh or the area of the respective cell in a general finitevolume discretization. As the model describes the total sediment flux and not flux per unit width, an integration over the edges of the cell is not necessary.
Inserting Eqs. (2) and (12) into Eq. (1) then yields the simplest form of a transportlimited fluvial erosion model,
where Q_{i} is defined by Eqs. (10) or (11).
As mentioned above, using powerlaw functions for sediment transport is not new. In combination with empirical relations for the channel width, physicallybased relations for the sediment flux density (e.g., Willgoose et al., 1991b) support the hypothesis of a powerlaw dependence of Q on A and S (Eq. 11). However, the relations where never written in such a simple form as in Eq. (11) with parameters that are related so closely to the concepts of concavity index and steepness index (Eq. 3). Equation (11) was discussed in the literature (e.g., Whipple and Tucker, 2002) in the context of equilibrium river profiles, but apparently never used directly for defining a transportlimited erosion model. In view of Hack's findings this is, however, as straightforward as describing detachmentlimited erosion by Eqs. (4) or (6). Even the physical unit of the erodibility K is the same in both models, and the same values of K yield the same erosion rate at the same topography for spatially uniform erosion.
The two models are, however, not equivalent for nonuniform erosion. According to Eqs. (4) and (5), the steepness index k_{s} directly reflects the erosion rate at the considered point in the form
for the detachmentlimited model. In turn, Eq. (11) of the transportlimited model can be combined with Eq. (9) to become
This relation is basically the same as Eq. (14) except that the righthand side is the mean erosion rate of the upstream catchment instead of the local erosion rate at the considered point. So channel steepness directly reflects the local erosion rate in the detachmentlimited model, but the mean erosion rate of the catchment for the transportlimited model.
The same holds for the interpretation of the erodibility K. In the detachmentlimited model, it describes how much material is eroded at the considered location for a given steepness index. In turn, it describes how much material is eroded on average in the upstream catchment in the transportlimited model. From a processoriented point of view, K would rather be considered a transport coefficient than an erodibility here. However, this is just a matter of terminology.
The model proposed in the previous section can be treated with an efficient, fully implicit numerical scheme in the linear case (n=1). The reason why this is possible in contrast to the 2D diffusion equation lies in the tree structure of the flow and sediment transport pattern.
The fully implicit discretization of Eq. (13) reads
where the time step extends from t_{0} to t and $\mathit{\delta}t=t{t}_{\mathrm{0}}$. The solution at t_{0} is known, and the solution at t is computed.
Let the node b be the flow target of the node i, so H_{b} serves as a base level for the node i. As the entire problem is linear, the height H_{i} responds linearly to base level changes. Figure 1 illustrates this behavior in a simple numerical example of a river (solid lines) with one tributary (dashed lines). The initial state (t=0, blue) was a steady state under constant uplift. The red curves (t=1) show the result of an implicit time step without uplift and with the same base level (H_{b}=0) as for t=0, while the orange curves correspond to different base levels H_{b}. The four red and orange curves of each river are equidistant at each point x for equal increments in H_{b}, so the change in height H_{i}(t) due to changes in base level is proportional to the change in base level H_{b}(t).
Due to the linearity, the sediment flux Q_{i} to the node b also responds linearly to base level changes and can therefore be written in the form
Here, ${Q}_{i}^{\mathrm{0}}$ is the flux that occurs if the base level H_{b} remains constant (H_{b}(t)=H_{b}(t_{0})), and ${Q}_{i}^{\prime}$ is the derivative of Q_{i}(t) with regard to base level changes. Inserting Eq. (17) for the donors into Eq. (16) yields
and thus
with the terms
introduced in order to keep the equations short. Similarly to the detachmentlimited model, nodes without any donors act as boundaries within the domain. These nodes do not require any specific treatment except that the respective sums in Eqs. (18) and (20) are empty.
The channel slope at the node i is
where d_{i} is the distance between the nodes i and b. So the sediment flux is
according to Eq. (11) for n=1. This leads to
Inserting this relation into Eq. (19) yields
which can be rearranged in the form
Comparing this expression with Eq. (17) yields
and
Equations (26) and (27) allow for the computation of ${Q}_{i}^{\mathrm{0}}$ and ${Q}_{i}^{\prime}$ from the respective values of the donors (because α_{i} and β_{i} depend on these) and from known elevation values at the time t_{0}. All values ${Q}_{i}^{\mathrm{0}}$ and ${Q}_{i}^{\prime}$ can thus be computed successively in the downstream direction. As the required order of the nodes is the same as for computing the catchment sizes A_{i}, it is most efficient to calculate ${Q}_{i}^{\mathrm{0}}$ and ${Q}_{i}^{\prime}$ in the same sweep over the nodes where the catchment sizes are computed.
Once the values ${Q}_{i}^{\mathrm{0}}$ and ${Q}_{i}^{\prime}$ have been computed for all nodes, the sediment flux Q_{i}(t) can be computed using Eq. (17). This sediment flux is then used for computing the elevation H_{i}(t) from Eq. (23). As these steps require the elevation of the flow target H_{b}(t), they have to be performed successively in upstream order. This order is the same as that used in the implicit scheme for detachmentlimited erosion.
So the numerical scheme consists of three sweeps over the grid:
 Sweep 1:

Compute the flow directions b of all nodes. The nodes can be processed in any order.
 Sweep 2:

Compute the catchment size A and the properties Q^{0} (Eq. 26) and Q^{′} (Eq. 27) of all nodes. The nodes have to be processed in downstream order. This is implemented most conveniently in a recursive scheme with a function that computes the three above properties for each node. Before computing these values, the function checks which of the donors have already been treated and invokes itself for those donors that have not been considered before.
 Sweep 3:

Compute Q(t) according to Eq. (17) and H(t) from Eq. (23) for all nodes. The nodes must be processed in upstream order, which is also performed conveniently by a recursive implementation. The principle is the same as in sweep 2 except that the flow target has to be considered instead of the donors.
The scheme is a direct scheme without any iterative component. The derivatives Q^{′} are always negative (lower base level leads to a higher sediment flux), so that the properties α and thus the denominator in Eqs. (26) and (27) are always positive. So the scheme is unconditionally stable, and its time complexity is linear (O(N)) under all conditions.
The workflow with the three sweeps is basically the same as in the implicit scheme for detachmentlimited erosion. The structure is the same without any extra loops, conditions, or functions to be invoked. Additional effort only arises from floatingpoint operations. Table 1 provides an estimate of the time complexity compared to detachmentlimited erosion. All results were obtained using the landform evolution model OpenLEM that was used in some previous studies (e.g., Robl et al., 2017; Wulf et al., 2019; Hergarten, 2020a) but has not been published explicitly. A regular 5000×5000 grid was used, and the results of several runs involving 100 to 1000 time steps were checked for consistency. The CPU time was normalized to the total effort of one time step for detachmentlimited erosion. The difference in time complexity between both models is marginal.
With regard to memory complexity, the scheme presented here requires two additional variables per node, Q^{0} and Q^{′}. When performing the third sweep, one of them can be recycled for storing the original surface height H(t_{0}) that is needed later when Eq. (17) is applied to the donors. The remaining variable can be used for storing the actual sediment flux Q(t) in case it is needed later.
The transportlimited model proposed in Sect. 2 is equivalent to the detachmentlimited model only for uniform erosion. Transient states are typically characterized by spatially variable erosion, so the two endmembers cannot yield the same transient behavior. This result is, however, already clear from more general arguments since both endmembers are described by differential equations of different types (parabolic vs. hyperbolic) as discussed in Sect. 1.
This section presents a numerical example showing that nonuniform conditions result in strong differences between the two models even in a steady state. The example uses a square domain of 5000×5000 nodes. The northern and southern boundaries are kept at H=0, while the two other boundaries are periodic. All horizontal lengths and areas are measured in terms of pixels. An exponent m=0.5 was assumed, so that equilibrium rivers have a concavity index of θ=0.5 for the linear model (n=1). The erodibility was set to K=1.
An equilibrium topography obtained for uniform uplift U=1 was used as a reference. This topography (Fig. 2, left) was generated by starting from a flat initial topography with a small random disturbance. As the transportlimited and the detachmentlimited models are equivalent for uniform erosion, this topography is an equilibrium topography for both models.
As a simple nonuniform uplift pattern, tentshaped uplift is considered. The maximum uplift rate U=1 is achieved here in the middle between the northern and southern boundary (x_{2}=2500) and decreases linearly to zero towards the boundaries. In order to get similar flow patterns (Fig. 2), the equilibrium topography corresponding to constant uplift was used as an initial condition.
All equilibrium topographies were computed by starting with small increments δt that are increased through time when the number of changes in the flow direction per time step is sufficiently small. This procedure is useful for generating steadystate topographies with similar largescale flow patterns at a reasonable number of time steps. At large δt, smaller random values of δt were used in each second step in order to avoid periodic oscillations between topographies with different flow patterns that prevent the topography from reaching a steady state.
The tentshaped uplift pattern causes an overall increase in uplift in upstream directions, at least for large rivers. This increase results in an upstream increase in steepness. As the steepness reflects the mean erosion rate of the upstream catchment (Eq. 15) instead of the local erosion rate for transportlimited erosion, it varies more gently with the uplift rate here than for detachmentlimited erosion.
Figure 3 shows swath profiles through the three topographies. The maximum surface height (uppermost curve of the respective color) is dominated by the steep slopes at small catchment sizes. Since these depend on the local uplift rate in equilibrium, the maximum elevation roughly follows the tentshaped uplift pattern with minor differences between transportlimited and detachmentlimited erosion. The absolute difference between the two models is similar for maximum, mean, and minimum elevation, so it can be attributed to the different heights of large valleys, while local relief is similar.
The profiles of the large rivers marked in Fig. 2 are shown in Fig. 4. For a clearer representation, the longitudinal coordinate was χ transformed according to Eq. (7) with A_{0}=1. With the value K=1 used here, equilibrium profiles follow a straight line H=χ at a uniform uplift rate U=1. In turn, χtransformed equilibrium profiles are concave if the uplift rate increases in the upstream direction. This concavity is weaker for the transportlimited model than for the detachmentlimited model as the local slope reflects the mean erosion rate of the upstream catchment, while it reflects the local erosion rate for detachmentlimited erosion. In the upper part of the catchment, however, both turn into parallel straight lines. In the lower part of the catchment, the river profiles of the transportlimited model are steeper than those of the detachmentlimited model because the river also has to carry away the material from the upper part with high erosion rates.
While the χtransformed river profiles of the transportlimited model are more straight than for detachmentlimited erosion, local collinearity of tributaries is lost. For detachmentlimited erosion, profiles of tributaries start with the same slope as the trunk stream and deviate more and more with increasing distance. In contrast, tributaries and the trunk stream may contribute different amounts of sediment per catchment size due to different mean erosion rates in their upstream catchments, which leads to different slopes immediately above the point of confluence in the transportlimited model. As a consequence, the capture of tributaries leads to stable knickpoints in the trunk stream for transportlimited erosion.
The most important lesson to be learned from this simple example concerns the estimation of the exponent n. Figure 5 shows the relation between the steepness index k_{s} and the erosion rate E, which is the same as the uplift rate U in a steady state. According to Eq. (14), the erosion rate is proportional to ${k}_{\mathrm{s}}^{n}$ in the detachmentlimited regime. So comparing k_{s} at different locations exposed to different erosion is a common approach to estimate n (e.g., Lague, 2014).
The detachmentlimited model with K=1 and n=1 reproduces the expected relation E=k_{s}, while this is not the case for the transportlimited model. Here the curve for the trunk stream (blue line) rather looks like a straight line with an offset. If we analyzed this curve without knowing that it originated from the transportlimited model, we would find that erosion starts at a threshold steepness index k_{s}≈0.7. If only a few points from this line were available, we would arrive at a nonlinear relation with n>1. This is, however, an extreme example as $E=U=\mathrm{0}$ at the boundary, while the sediment flux from the domain requires a nonzero channel steepness. Qualitatively, the result would be similar in any situation where the uplift rate increases in an upstream direction. If we interpret a long transportlimited river profile in terms of detachmentlimited erosion, the exponent n would be systematically overestimated.
In turn, comparing the three tributaries that are predominantly oriented in the east–west direction would yield an exponent close to the correct exponent n=1 used here. The reason is that these tributaries are not subject to strong variations in uplift rate. All estimates reviewed by Lague (2014) suggest n>1 except for one data set. This data set describes strikeparallel tributaries originally investigated by Kirby and Whipple (2001) where a reanalysis by Wobus et al. (2006) resulted in n≈1. This finding sheds new light on the apparent evidence for exponents n>1 obtained from analyzing river profiles under nonuniform conditions. An unrecognized contribution of sediment transport may result in an overestimation of n here. This problem makes estimating the effective values of n even more difficult and deserves further consideration in the future.
While fully implicit schemes for diffusiontype equations are unconditionally stable for arbitrary δt, their numerical error increases linearly with δt. This error is a systematic error in the sense that the response of the system to temporal changes (e.g., in uplift here) is always too slow.
In order to assess this error and to estimate a reasonable maximum δt, the response of the steadystate topography with unit uplift considered in Sect. 4 to changes in uplift rate is investigated in the following. Due to the linearity of the model, the response to changes in U is also linear as long as the flow directions do not change, which is the case if the change in U is sufficiently small. Technically, the uplift rate can simply be set to U=0 at t=0, and the simulation is performed without recomputing the flow directions.
Figure 6 shows the results for different values of δt in terms of the total sediment flux Q out of the domain or, more precisely, in terms of the change in this flux per change in uplift rate, $\frac{\mathit{\delta}Q}{\mathit{\delta}U}$. The numerical error is in the order of magnitude of some percent for δt=1. It should, however, be emphasized that this error is just some kind of time lag, while the sediment balance itself is satisfied exactly. The curves for δt≤0.1 are hardly distinguishable in the diagram. If we assume, for example, K=2.5 Myr^{−1} (Robl et al., 2017), a unit of nondimensional time corresponds to 400 000 years. So time increments in the order of magnitude of 100 000 years should be no problem concerning the accuracy of the implicit scheme even for large grids (5000×5000 nodes in this example).
While the numerical error of the implicit scheme is even smaller than for the detachmentlimited model (Fig. 6), it must be taken into account that all models in this field compute flow directions and changes in topography in separate steps. Changes in the flow directions introduce an additional numerical error. In the worst case, the implicit scheme acts on a deprecated flow pattern over almost the entire time step. So this error is formally also linear in δt, but strongly depends on the topography. Regions with small relief are particularly susceptible to artifacts. Under erosion, unreasonable river networks may deeply incise at large δt, and reorganization may take a long time afterwards. In the aggradational regime, large rivers may even turn into weird ridges within a single, large time step.
In many situations, the limitation of the maximum δt arising from changes in the flow directions is more severe than the numerical error of the implicit scheme itself. However, as it strongly depends on the topography, it is difficult to provide an estimate for a reasonable δt then. Practically, tracking the number of changes in the flow direction and adjusting δt so that the number of changes per time step does not exceed a given threshold provides a feasible criterion.
Detachmentlimited erosion and transportlimited erosion can be seen as endmembers of a more general framework. In particular, the extension of the detachmentlimited model by sediment transport proposed by Davy and Lague (2009) is receiving growing interest in this context. Recently, Guerit et al. (2019) derived estimates of the sediment deposition parameter occurring in this model from analyzing natural and experimental topographies. The authors concluded that “natural landscapes seem to describe a continuum between the two modes with a preference for TL mode” (transportlimited mode) as already suggested by Davy and Lague (2009).
Whipple and Tucker (2002) already proposed the generic form of the model of Davy and Lague (2009) and coined the term “linear decline model”. This concept starts from the detachmentlimited model and assumes that the sediment flux reduces the ability of the river to erode. Assuming that the decrease in erosion rate is linear, this leads to the expression
(from Eq. 4), where ψ is an arbitrary function.
In addition to Eq. (28), the sediment balance (Eq. 2 or the respective discrete form, Eq. 12) must be satisfied. Inserting Eq. (28) and the sediment balance into Eq. (1) yields a system of two coupled partial differential equations for the surface height H and the sediment flux Q.
The sediment balance can be written conveniently in integral form (Eq. 9). Combining this expression with Eqs. (1) and (28) yields a single integrodifferential equation for the surface height,
instead of a system of two differential equations.
The detachmentlimited endmember corresponds to ψ(A)=0. In this case, the two differential equations are decoupled, so the equation for H can be solved without computing Q. Approaching the detachmentlimited endmember is mathematically more complicated. This can be achieved by increasing f and ψ in such a way that f→∞ and ψ→∞, while the ratio $\frac{f}{\mathit{\psi}}$ converges to a finite, nonzero value. Then, Eq. (28) turns into
The resulting sediment flux Q defines the transport capacity.
If we request that equilibrium river profiles under uniform conditions are still consistent with Hack's findings (Eq. 3), the entire erosion rate defined by Eq. (28) must be a function of the product A^{θ}S. Inserting Eq. (10) into Eq. (28) yields
So ψ(A) must be inversely proportional to A,
with a nondimensional constant G. This is exactly the relation proposed by Davy and Lague (2009) (with Θ instead of G there). The second term in Eq. (28) turns into $\frac{GQ}{A}$, which was interpreted as deposition of sediments by Davy and Lague (2009).
The equilibrium erosion rate under uniform conditions is
in this model. So sediment deposition effectively reduces the erosion rate by a factor of 1+G under uniform equilibrium conditions, which makes equilibrium profiles by a factor of 1+G steeper in the linear model (n=1) as already stated by Yuan et al. (2019). This effect can be compensated by rescaling the erodibility K by a factor of 1+G, which modifies the model of Davy and Lague (2009) to
Both versions differ concerning the interpretation of the erodibility K. While it characterizes the process of detachment in the original model, it is interpreted as the fingerprint of spatially uniform erosion including sediment transport in the rescaled version. In contrast to the original version, the rescaled version also captures the transportlimited endmember for G→∞ since
The linear decline model can be interpreted in several ways. If we define an equilibrium sediment flux by
Eq. (28) turns into
This is the undercapacity model (Kooi and Beaumont, 1994) written in terms of sediment flux instead of flux per unit width (Eq. 8).
The formulation of transportlimited erosion proposed in Sect. 2 allows for an alternative definition of a hybrid model that can also be interpreted as a linear decline model. Let us write the detachmentlimited endmember (Eq. 6) in the form
and the transportlimited endmember (Eq. 11) in the form
In contrast to the previous considerations, different symbols K_{d} and K_{t} are used here. As discussed above, their meaning is in principle the same, but there is no reason why the values should be the same under all conditions.
The simplest combination of both endmembers is assuming that the property A^{m}S^{n} that is responsible for both detachment and transport is shared among the two processes, i.e.,
This model approaches the detachmentlimited regime for zero sediment flux and the transportlimited regime for high sediment flux.
Equation (40) can also be written in the form
so
and
The formulation defined by Eq. (40) could be called “shared stream power model”. Compared to the concepts of detachment and deposition and the undercapacity model, this is rather a generic model. In turn, the formulation in terms of K_{d} and K_{t} may help to understand rivers passing different lithologies. Here we could expect that K_{d} shows a stronger variation than K_{t}, although the ability to transport material also depends on the characteristics of the sediments at the river bed.
The numerical scheme described in Sect. 3 can be extended towards the linear version of the linear decline model, i.e., if the first term in Eq. (28) is also linear in channel slope S (n=1 in Eq. 6). The general form of this model reads
with any functions ϕ and ψ. For the rescaled version of the model proposed by Davy and Lague (2009) and the shared stream power version, the two functions are
while the term 1+G would not occur in the original version.
Inserting Eq. (44) into the general landform evolution model (Eq. 1) yields
and after inserting difference quotients for time derivative and channel slope
This equation can be rearranged in the form
Plugging this result into Eq. (19), rearranging the resulting equation to yield Q_{i}(t), and comparing the obtained expression to Eq. (17) finally yields
and
The scheme is very similar to that presented in Sect. 3 for transportlimited erosion. Equations (50) and (51) have to be used in sweep 2. Sweep 3 is now based on Eq. (49), while Eq. (17) is still used in its original form.
Preliminary numerical tests revealed that the time complexity of this version is very close to the transportlimited case, while that of the iterative scheme proposed by Yuan et al. (2019) is close to the detachmentlimited case in each iteration step. In the first iteration, sweep 2 computes the catchment sizes here, while it integrates the upstream erosion rate to yield the sediment flux (Eq. 9) in subsequent iterations. Taking the values from Table 1, there would be a slight advantage of the iterative scheme (100 % vs. 112 %) if the iterative scheme could be applied with a single iteration step. This is, however, not possible if the flow direction of any node has changed because the sediment flux Q is not available for the actual flow pattern then. In the best case, the iterative scheme requires two steps. According to the numerical tests of Yuan et al. (2019), this is achieved for small values of G in the order of magnitude of 0.01 at n=1. This yields 112 % vs. 162 % effort, so the direct scheme is at least 30 % faster than the iterative procedure. The advantage of the direct scheme rapidly grows with increasing sediment transport. The data repository of the recent study of Guerit et al. (2019) found a median value of G=1.6 for n=1 by analyzing several natural river profiles. The iterative scheme requires about eight iterations at this value, resulting in an effort of 112 % vs. 534 %. So the direct scheme is almost 5 times faster under these conditions. Guerit et al. (2019) also reported on higher values of G, where the convergence of the iterative scheme would become very slow. In addition, the direct scheme has the advantage of an exact solution without the need for checking convergence.
7.1 Adding transportlimited and detachmentlimited erosion
The models of the linear decline type discussed in Sect. 6 enforce a strict balance for the sediment flux for any nonzero function ψ(A). However, we could also assume that a part of the eroded material is immediately excavated, while the rest is transported. This could be seen as a first step towards considering different particle sizes where one class of particles it so finegrained that these will not be deposited. For simplicity, this version is elaborated only as an extension of the transportlimited model here, although a combination with the linear decline model is also possible. The linear version of this model reads
or inserted into Eq. (1) and discretized in fully implicit form
Here, Γ is any function that describes the excavation of material. Similarly to the functions ϕ and ψ used in the previous section, Γ may in principle depend on all properties except for surface heights in order to maintain the linearity. It may be tempting to use the expression $\mathrm{\Gamma}=\stackrel{\mathrm{\u0303}}{K}{A}^{m}$ from the detachmentlimited model, where $\stackrel{\mathrm{\u0303}}{K}$ has the same meaning and physical unit as K. However, Eq. (52) combines a sediment balance with immediate excavation, which causes a scaling problem if rivers are considered as linear objects (Howard, 1994; Perron et al., 2008; Pelletier, 2010; Hergarten, 2020a). As a consequence, an additional rescaling factor depending on the pixel size must be introduced in the definition of Γ in order to avoid an artificial dependence of the results on the spatial resolution. Different approaches for this scaling factor are discussed in the above references.
Apart from this scaling problem, the numerical implementation is straightforward. Using Eq. (11), the last term in Eq. (53) can be expressed as
This results in a factor $\mathrm{1}+\frac{{s}_{i}\mathrm{\Gamma}}{K{A}_{i}^{m+\mathrm{1}}}$ in front of Q_{i}(r) in Eq. (19). This factor propagates to the denominator of Eqs. (26) and (27), so that we finally arrive at
and
7.2 Hillslope diffusion
Linear diffusion (e.g., Culling, 1960) as the simplest model of hillslope erosion can be implemented more efficiently than in the detachmentlimited model because the flux components in the direction of the channel network can be integrated into the implicit scheme. If d_{ij} is the distance between the node i and a neighbor j and l_{ij} is the length of the respective edge in a finitevolume representation, the diffusive flux in this direction is
where D is the diffusivity. In contrast to the discrete divergence of the fluvial sediment flux density (Eq. 12), this simple expression is only valid if the edge is normal to the connecting line, as is the case, for example, in a Voronoi discretization.
An implicit scheme for the flux components in the direction of the channel network combined with an explicit scheme for the other components requires an additional variable B_{i} (where practically either ${Q}_{i}^{\mathrm{0}}$ or ${Q}_{i}^{\prime}$ might be used) for the balance of the diffusive fluxes. For each node i, a loop over all neighbors j with H_{j}<H_{i} except for the flow target b is employed. The respective values Q_{ij} are added to B_{j} and subtracted from B_{i}. After considering all nodes i, the values $\frac{{B}_{i}}{{s}_{i}}$ are added to the uplift rates U_{i}. This part of the scheme captures the diffusive fluxes except for those in the direction of the channel network. Comparing Eq. (57) to Eq. (22), it is easily recognized that the diffusive flux from each node i to its flow target b can be included by replacing the term $K{A}_{i}^{\mathit{\theta}+\mathrm{1}}$ by $K{A}_{i}^{\mathit{\theta}+\mathrm{1}}+D{l}_{ib}$ throughout the calculations of Sect. 3.
As this scheme is not fully implicit, the maximum time increment is still limited. However, as the flux component in the flow direction is the largest among all, its partly implicit treatment improves the stability of the diffusion term and thus increases the maximum possible time increment.
While the approach presented here is efficient and can be applied to a large class of problems, some limitations should also be mentioned.
First, any kind of sediment transport that transfers material from one site to more than one target site destroys the treelike topology of sediment fluxes. Such processes are thus not compatible with the implicit scheme presented here. This applies to hillslope processes as well as to fluvial processes with multiple flow directions as implemented, e.g., in the model TTLEM (Campforts et al., 2017). However, the implicit scheme for detachmentlimited erosion is subject to the same limitation.
Concerning numerics, nonlinearity is the only point where the approach suggested here falls behind the implicit scheme for detachmentlimited erosion. The latter can be solved directly for n=1 and for n=2 (Hergarten, 2002) but can be treated by finding the roots of a scalar nonlinear equation at each point for any value of n. In contrast, nonlinearity can only be included in the approach proposed here either by treating the nonlinear terms in an explicit manner or to employ an iteration.
Finally, the treatment of lakes, i.e., local depressions in the topography, is a problem. In the detachmentlimited model, local depressions result in negative channel slopes and thus in negative erosion rates without any specific treatment. However, these negative erosion rates can be cut off easily in the implicit scheme. In the transportlimited model, local depressions result in a sediment flux opposite to the flow direction. Erosion of dams may be too fast then, so that the lifetime of lakes may be too short. This effect cannot be fixed easily in the fully implicit scheme.
This study proposes a simple formulation of transportlimited fluvial erosion. This formulation can be immediately reconciled with the empirical results of Hack (1957) on longitudinal river profiles. The interpretation of Hack's findings as the fingerprint of spatially uniform erosion is equivalent for transportlimited erosion and for detachmentlimited erosion where it has been widely used. In turn, the behavior of both models differs if erosion is nonuniform.
As a main point, a new numerical scheme for treating transportlimited erosion with a fully implicit discretization in time was presented. It is a direct solver without any iteration and is unconditionally stable for arbitrarily large time increments. It is of linear time complexity (O(N)) where the computing effort is marginally higher than for detachmentlimited erosion. The scheme can also be applied to combined linear models of detachmentlimited erosion and sediment transport such as the linear decline model. Here it also allows for approaching the transportlimited endmember without any loss of performance and provides a numerical efficiency that is better than the iterative scheme suggested by Yuan et al. (2019).
All codes and computed data can be downloaded from the FreiDok data repository (Hergarten, 2020b) (https://doi.org/10.6094/UNIFR/166660). The author is happy to assist interested readers in reproducing the results and performing subsequent research.
The author declares that there is no conflict of interest.
The author would like to thank Xiaoping Yuan, Wolfgang Schwanghart, Jean Braun, and an anonymous reviewer for their constructive and encouraging comments.
The article processing charge was funded by the BadenWuerttemberg Ministry of Science, Research and Art and the University of Freiburg in the funding programme Open Access Publishing.
This paper was edited by Jean Braun and reviewed by Wolfgang Schwanghart, Xiaoping Yuan, and one anonymous referee.
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
Campforts, B., Schwanghart, W., and Govers, G.: Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model, Earth Surf. Dynam., 5, 47–66, https://doi.org/10.5194/esurf5472017, 2017. a, b
Coulthard, T. J.: Landscape evolution models: a software review, Hydrol. Process., 15, 165–173, https://doi.org/10.1002/hyp.426, 2001. a
Culling, W.: Analytical theory of erosion, J. Geol., 68, 336–344, https://doi.org/10.1086/626663, 1960. 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, c, d, e, f, g, h, i
Duvall, A. R. and Tucker, G. E.: Dynamic ridges and valleys in a strikeslip environment, J. Geophys. Res.Earth, 120, 2016–2026, https://doi.org/10.1002/2015JF003618, 2015. a
Eizenhöfer, P. R., McQuarrie, N., Shelef, E., and Ehlers, T. A.: Landscape response to lateral advection in convergent orogens over geologic time scales, J. Geophys. Res.Earth, 124, 2056–2078, https://doi.org/10.1029/2019JF005100, 2019. a
Ferrier, K. L., Perron, J. T., Mukhopadhyay, S., Rosener, M., Stock, J. D., Huppert, K. L., and Slosberg, M.: Covariation of climate and longterm erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua'i, GSA Bull., 125, 1146–1163, https://doi.org/10.1130/B30726.1, 2013. a
Flint, J. J.: Stream gradient as a function of order, magnitude, and discharge, Water Resour. Res., 10, 969–973, https://doi.org/10.1029/WR010i005p00969, 1974. a
Goren, L., Fox, M., and Willett, S. D.: Tectonics from fluvial topography using formal linear inversion: Theory and applications to the Inyo Mountains, California, J. Geophys. Res.Earth, 119, 1651–1681, https://doi.org/10.1002/2014JF003079, 2014. a
Guerit, L., Yuan, X. P., Carretier, S., Bonnet, S., Rohais, S., Braun, J., and Rouby, D.: Fluvial landscape evolution controlled by the sediment deposition coefficient: Estimation from experimental and natural landscapes, Geology, 47, 853–856, https://doi.org/10.1130/G46356.1, 2019. a, b, c
Hack, J. T.: Studies of longitudinal profiles in Virginia and Maryland, no. 294B in US Geol. Survey Prof. Papers, US Government Printing Office, Washington D.C., https://doi.org/10.3133/pp294B, 1957. a, b, c
Hackbusch, W.: MultiGrid Methods and Applications, Springer, Berlin, Heidelberg, New York, https://doi.org/10.1007/9783662024270, 1985. a
Harel, M.A., Mudd, S. M., and Attal, M.: Global analysis of the stream power law parameters based on worldwide 10Be denudation rates, Geomorphology, 268, 184–196, https://doi.org/10.1016/j.geomorph.2016.05.035, 2016. a
Hergarten, S.: SelfOrganized Criticality in Earth Systems, Springer, Berlin, Heidelberg, New York, https://doi.org/10.1007/9783662043905, 2002. a, b
Hergarten, S.: Rivers as linear elements in landform evolution models, Earth Surf. Dynam., 8, 367–377, https://doi.org/10.5194/esurf83672020, 2020a. a, b
Hergarten, S.: Transportlimited fluvial erosion – simple formulation and efficient numerical treatment: codes and data, FreiDok plus, Universitätsbibliothek Freiburg, https://doi.org/10.6094/UNIFR/166660, 2020b. a
Hergarten, S. and Neugebauer, H. J.: Selforganized critical drainage networks, Phys. Rev. Lett., 86, 2689–2692, https://doi.org/10.1103/PhysRevLett.86.2689, 2001. a
Howard, A. D.: A detachmentlimited model for drainage basin evolution, Water Resour. Res., 30, 2261–2285, https://doi.org/10.1029/94WR00757, 1994. a, b, c, d
Kirby, E. and Whipple, K. X.: Quantifying differential rock uplift rates via stream profile analysis, Geology, 29, 415–418, https://doi.org/10.1130/00917613(2001)029<0415:QDRURV>2.0.CO;2, 2001. a
Kooi, H. and Beaumont, C.: Escarpment evolution on highelevation rifted margins: insights derived from a surface process model that combines diffusion, advection and reaction, J. Geophys. Res., 99, 12191–12209, 1994. a, b
Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Process. Land., 39, 38–61, https://doi.org/10.1002/esp.3462, 2014. a, b, c, d
Pelletier, J. D.: Minimizing the gridresolution dependence of flowrouting algorithms for geomorphic applications, Geomorphology, 122, 91–98, https://doi.org/10.1016/j.geomorph.2010.06.001, 2010. a
Perron, J. T. and Royden, L.: An integral approach to bedrock river profile analysis, Earth Surf. Process. Land., 38, 570–576, https://doi.org/10.1002/esp.3302, 2013. a
Perron, J. T., Dietrich, W. E., and Kirchner, J. W.: Controls on the spacing of firstorder valleys, J. Geophys. Res.Earth, 113, F04016, https://doi.org/10.1029/2007JF000977, 2008. a
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
Royden, L. and Perron, J. T.: Solutions of the stream power equation and application to the evolution of river longitudinal profiles, J. Geophys. Res.Earth, 118, 497–518, https://doi.org/10.1002/jgrf.20031, 2013. a
Theodoratos, N., Seybold, H., and Kirchner, J. W.: Scaling and similarity of a streampower incision and linear diffusion landscape evolution model, Earth Surf. Dynam., 6, 779–808, https://doi.org/10.5194/esurf67792018, 2018. a
Turowski, J. M.: Semialluvial channels and sedimentfluxdriven bedrock erosion, in: Gravel‐Bed Rivers, edited by: Church, M., Biron, P., and Roy, A., chap. 29, John Wiley & Sons, Ltd, 399–418, https://doi.org/10.1002/9781119952497.ch29, 2012. a
van der Beek, P.: Modelling landscape evolution, in: Environmental Modelling: Finding Simplicity in Complexity, edited by: Wainwright, J. and Mulligan, M., 2 edn., WileyBlackwell, Chichester, 309–331, 2013. a
Whipple, K. X. and Tucker, G. E.: Dynamics of the stream power river incision model: Implications for height limits of mountain ranges, landscape response time scales and research needs, J. Geophys. Res., 104, 17661–17674, https://doi.org/10.1029/1999JB900120, 1999. a
Whipple, K. X. and Tucker, G. E.: Implications of sedimentfluxdependent river incision models for landscape evolution, J. Geophys. Res., 107, 2039, https://doi.org/10.1029/2000JB000044, 2002. a, b, c, d
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/B9780123747396.002268, 2013. a
Willgoose, G.: Mathematical modeling of whole landscape evolution, Annu. Rev. Earth Planet. Sci., 33, 443–459, https://doi.org/10.1146/annurev.earth.33.092203.122610, 2005. a
Willgoose, G., Bras, R. L., and RodriguezIturbe, I.: A physical explanation of an observed link areaslope relationship, Water Resour. Res., 27, 1697–1702, https://doi.org/10.1029/91WR00937, 1991a. a
Willgoose, G., Bras, R. L., and RodriguezIturbe, I.: Results from a new model of river basin evolution, Earth Surf. Proc. Land., 16, 237–254, https://doi.org/10.1002/esp.3290160305, 1991b. a, b, c
Wobus, C., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., Crosby, B., and Sheehan, D.: Tectonics from topography: Procedures, promise, and pitfalls, in: Tectonics, Climate, and Landscape Evolution, edited by: Willett, S. D., Hovius, N., Brandon, M. T., and Fisher, D. M., vol. 398 of GSA Special Papers, Geological Society of America, Boulder, Washington, D.C., 55–74, https://doi.org/10.1130/2006.2398(04), 2006. a
Wulf, G., Hergarten, S., and Kenkmann, T.: Combined remote sensing analyses and landform evolution modeling reveal the terrestrial Bosumtwi impact structure as a Marslike rampart crater, Earth Planet. Sc. Lett., 506, 209–220, https://doi.org/10.1016/j.epsl.2018.11.009, 2019. 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, https://doi.org/10.1029/2018JF004867, 2019. a, b, c, d, e, f
 Abstract
 Introduction
 Simple formulation of transportlimited erosion
 A fully implicit numerical algorithm for transportlimited erosion
 A numerical example
 Numerical accuracy
 Extension towards the linear decline model
 Further extensions
 Limitations
 Conclusions
 Code and data availability
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Simple formulation of transportlimited erosion
 A fully implicit numerical algorithm for transportlimited erosion
 A numerical example
 Numerical accuracy
 Extension towards the linear decline model
 Further extensions
 Limitations
 Conclusions
 Code and data availability
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References