Articles | Volume 8, issue 2
Research article
26 May 2020
Research article |  | 26 May 2020

Rivers as linear elements in landform evolution models

Stefan Hergarten

Models of detachment-limited fluvial erosion have a long history in landform evolution modeling in mountain ranges. However, they suffer from a scaling problem when coupled to models of hillslope processes due to the flux of material from the hillslopes into the rivers. This scaling problem causes a strong dependence of the resulting topographies on the spatial resolution of the grid. A few attempts based on the river width have been made in order to avoid the scaling problem, but none of them appear to be completely satisfying. Here a new scaling approach is introduced that is based on the size of the hillslope areas in relation to the river network. An analysis of several simulated drainage networks yields a power-law scaling relation for the fluvial incision term involving the threshold catchment size where fluvial erosion starts and the mesh width. The obtained scaling relation is consistent with the concept of the steepness index and does not rely on any specific properties of the model for the hillslope processes.

1 Introduction

Fluvial incision is a major if not dominant component of long-term landform evolution in orogens. When modeling fluvial erosion, restriction to the detachment-limited regime considerably simplifies the equations. Here it is assumed that the erosion rate at any point of a river can be predicted from local properties such as discharge and slope, while sediment transport is not considered. The generic differential equation for the topography H(x1,x2,t) of a landform evolution model with detachment-limited fluvial erosion reads

(1) H t = U - E - div q ,

where U is the uplift rate and E the rate of fluvial incision. The third term describes a local transport process at the hillslopes, where q is the flux density and div the 2-D divergence operator. Linear diffusion is the simplest model here; it was considered in the context of landform evolution by Culling (1960) even before models of fluvial erosion came into play. However, there are also more sophisticated models for q that take the nonlinear dependencies of hillslope processes on topography into account (e.g., Andrews and Bucknam1987; Howard1994; Roering et al.1999).

Concerning the fluvial incision term E, assuming a power-law function of the catchment size A and the channel slope S,

(2) E = K A m S n ,

has become some kind of paradigm. The parameter K is denoted erodibility. 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).

Equation (2) is often called the stream-power approach, 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 Tucker1999). However, the idea behind this approach even dates back to the empirical study of longitudinal channel profiles by Hack (1957). In this study, a power-law relationship between channel slope and drainage area was found, often called Flint's law (Flint1974). This relationship is nowadays usually written in the form

(3) S = k s A - θ ,

where θ is the concavity index and ks the steepness index. Assuming that Eq. (3) is the fingerprint of spatially uniform steady-state conditions, it predicts mn=θ and allows for a convenient interpretation of the erodibility. If local transport (last term in Eq. 1) is neglected, the steepness index follows the relation

(4) k s n = E K .

This relation allows for a simple adjustment of the lumped parameter K in such a way that a given channel steepness is achieved at a given erosion rate.

2 The scaling problem

While widely used and in principle simple, all models of the type described by Eqs. (1) and (2) suffer from a scaling problem. Mathematically, the problem is that catchment sizes are not well-defined in the continuum limit as the catchment of each point degenerates to a line. When considered on a discrete grid, rivers are represented as linear objects with a width of one pixel. Thus, the total surface area of the pixels covering the network of the large rivers decreases with decreasing mesh width.

If local transport is not considered, the scaling problem leads to a canyon-like topography, where the width of the valleys decreases with mesh width. This behavior is illustrated in Figs. 1 and 2, where two steady-state topographies with mesh widths of δ=0.01 (100×100 nodes) and δ=0.002 (500×500 nodes) are considered. All parameter values are set to unity except for m=0.5 so that θ=0.5. The northern and southern boundaries are held at zero elevation, while the western and eastern boundaries are periodic. The topographies were obtained from the landform evolution model OpenLEM that was used in some previous studies (e.g., Robl et al.2017; Wulf et al.2019) but has not been published explicitly. It uses the D8 flow-routing scheme (O'Callaghan and Mark1984) and a fully implicit scheme (Hergarten and Neugebauer2001; Hergarten2002) so that large time steps can be performed in order to ensure that a steady state is achieved. The simulation on the fine grid was started from a flat topography with a small random disturbance, while the simulation on the coarse grid was started from a downsampled version of the finer topography.

Figure 1Fluvial equilibrium topographies computed for identical parameter values on grids with different spacing (δ=0.01, 100×100 nodes, and δ=0.002, 500×500 nodes). The horizontal lines refer to the profiles analyzed in Fig. 2 and the rectangle to the region shown in Fig. 4.


Figure 2Profiles through the topographies shown in Fig. 1.


Relief increases with decreasing grid spacing because the smallest catchment size that can be resolved is Amin=δ2, and the maximum equilibrium slope is proportional to Amin-θ=δ-2θ according to Eq. (3). As nodes with small catchment sizes can drain directly into large rivers, this increase is not restricted to major drainage divides but also results in steep valley flanks. The heights of the valley floors are, however, hardly affected by the spatial resolution. Catchment sizes of large rivers even converge in the limit δ→0 so that longitudinal profiles of large rivers become stable for δ→0 according to Eq. (3). Thus, relief and also mean elevation depend on the spatial resolution for the simplest model without local transport, while large rivers are hardly affected.

The independence of river steepness of resolution is, however, lost as soon as local transport comes into play. Figure 3 shows the example of short, parallel river segments with unit spacing (periodic in x2 direction) in equilibrium with constant uplift. Linear diffusion,

(5) q = - D H ,

was assumed as the simplest model for local transport. As in the previous example, all parameters except for m=0.5 were set to unity. A catchment size of A=106 was assumed for each river segment so that the channel slope should theoretically be S=10-3 in equilibrium with U=1. While the topography of the hillslopes is in principle independent of the grid spacing δ, the river segment becomes steeper if δ decreases.

Figure 3River segments in equilibrium with uplift for different mesh widths δ.


The reason for the increasing channel steepness is that the local transport is conservative, so the river not only has to incise into the rock at its bed but also has to remove the material coming from the hillslopes. Regardless of the model used for local transport, a flux of (dδ)U per river length enters the site that contains the river in equilibrium, where d is the valley spacing. Then the discretized divergence of the flux density is

(6) div q = - ( d - δ ) U δ .

Inserting this result into the steady-state version of Eq. (1) yields

(7) E = U - div q = d δ U ,

so the fluvial erosion rate required for compensating uplift is higher than it would be without local transport by a factor dδ. This requires an increase in the channel slope by a factor of dδ1n according to Eq. (2).

This scaling issue has been known for more than 25 years, and two approaches have been suggested to overcome the problem. Howard (1994) suggested a subpixel representation of the rivers, where a river segment only covers a fraction of a grid cell. It was assumed that this fraction is wδ, where w is the river width, and then the fluvial incision term E was multiplied by this factor. Perron et al. (2008) transferred this concept to the detachment-limited case. According to Eq. (7), rescaling E by the factor wδ yields

(8) E = d w U ,

so the dependency on δ indeed vanishes.

While straightforward at first sight, this scaling approach is not free of problems. The channel width in general increases in the downstream direction so that equilibrium river profiles are no longer consistent with Eq. (3). Perron et al. (2008) avoided this problem by assuming a constant channel width and postponing it to subsequent studies. As discussed by Pelletier (2010), taking an increase in channel width in the downstream direction into account would require a reduction of the exponent m in Eq. (2) in order to keep it consistent with Eq. (3). However, the unit and meaning of the erodibility K would change then.

In order to overcome this problem, Pelletier (2010) suggested leaving the fluvial incision term as is and rescaling the local transport term divq by the inverse factor δw at sites containing rivers. Practically, this rescaling means that the flux of material coming from the hillslopes is not distributed over the entire grid cell but only over the part of the area covered by the river. Thus it can be seen as the inverse of the subpixel approach of Howard (1994) and Perron et al. (2008) applied to the local transport instead of the fluvial erosion. For the steady-state example considered above, this rescaling leads to

(9) div q = - ( d - δ ) U w ,

instead of Eq. (6), so that

(10) E = U - div q = ( d + w - δ ) w U .

For wd and δd, however, this relation approaches Eq. (8), so this concept suffers from the same problem as the approach of Howard (1994) and Perron et al. (2008).

Thus there seems to be no completely satisfactory solution of the scaling problem so far. Several contemporary modeling studies (e.g., Duvall and Tucker2015; Gray et al.2018; Wulf et al.2019; Reitman et al.2019) use neither of the two approaches but implement Eq. (1) as is without taking its dependence on the grid scale into account. This is not a crucial problem as long as simulations with different spatial resolutions are not compared and as long as we are aware that the erodibility K has a limited meaning. As soon as the relevance of fluvial erosion and hillslope processes is assessed quantitatively or scaling relations are developed (e.g., Theodoratos et al.2018), the problem may become crucial. A further discussion is given in Sect. 5.

Other recent approaches navigate around the scaling problem by neglecting the flux of material from the hillslopes into the rivers. The recently presented landform evolution model TTLEM (Campforts et al.2017) makes a distinction by catchment size in such a way that fluvial erosion only acts on sites with a catchment size above a given threshold Ac, while hillslope processes only act at smaller catchment sizes. It is assumed that all hillslope material entering the rivers is immediately excavated without any further effect so that fluxes from hillslopes into rivers can be disregarded and the scaling problem does not occur. This approach reduces the interaction between rivers and hillslopes to a one-way coupling, where only the rivers have an influence on the evolution of the hillslopes and can be seen as an implementation of bedrock incision in the strict sense. While it seems that the terms detachment-limited erosion and bedrock incision are sometimes used synonymously, it should be clarified that the applicability of the concept of pure bedrock incision is probably much narrower than that of detachment-limited erosion, in particular if highly resistant material is brought into the channels (Shobe et al.2016). The same in principle holds for the model most widely used in the context of drainage divide migration (Goren et al.2014), where analytical solutions for hillslope processes are used on the subpixel scale.

3 A new scaling approach

The simple example considered in the previous section involves a dependence on grid spacing δ according to the factor dδ without rescaling (Eq. 7). Both approaches for rescaling replace the dependence on δ by a dependence on the channel width w so that a factor dw remains (Eq. 8). This is, however, still a problem if w is not constant. The occurrence of the factor dw suggests that the valley spacing d would be a more suitable characteristic length scale for rescaling than w if we want to preserve the form of the erosion law (Eq. 2) without changing the exponents m and n. In the following, a concept that generalizes the simple example of parallel rivers to dendritic networks is developed.

Let us start from the simplest approach to distinguish channel sites from hillslopes by defining a threshold catchment size Ac in such a way that all sites with AAc are river segments, while all sites with A<Ac belong to hillslopes. As local transport is conservative, all material eroded anywhere has to be removed by the river sites so that we have to determine how much material each river site receives from the hillslopes. The area of the respective hillslopes can be determined for a given topography without any specific assumptions on the transport process except for the direction of transport. The simplest model is to assume that local transport follows the hypothetic channel network at the hillslopes, i.e., the direction of steepest descent on a purely fluvial topography. Figure 4 illustrates this concept. Each colored area consists of one channel site and the hillslope area that delivers its eroded material to this site, i.e, of those sites that drain into the considered site without passing any other upstream channel site.

Figure 4Flow pattern of the central region of Fig. 1. Black lines show rivers with AAc for Ac=100 pixels. Gray lines are channels with A<Ac, considered to be hillslope sites. Each colored area consists of one channel site plus the hillslope area that drains into this site without passing another upstream channel site.


If the size of this area was the same for each river site, rescaling the fluvial erosion rate (Eq. 2) according to

(11) E = A e K A m S n ,

where Ae is the size of this area measured in DEM (digital elevation model) pixels (i.e., the number of sites), would already solve the scaling problem. However, it is immediately recognized in Fig. 4 that the sizes of these areas are highly variable. A random variation in these sizes is not a problem. If Ae in Eq. (11) is the mean size, channel steepness will just vary randomly, which is also found in nature. A systematic dependence of Ae on catchment size would, however, be a problem. In this case, equilibrium river profiles would be no longer consistent with Eq. (3), so the problem would be basically the same as in the previous approach for a non-constant channel width.

In the following, numerically obtained equilibrium drainage networks are analyzed in order to find out how Ae depends on A and on Ac. More precisely, Ae is the mean size of all hillslope areas draining into channel sites with a given catchment size A at a given fluvial threshold Ac (plus the respective channel site). For simplicity, all areas are measured in DEM pixels in the following considerations, i.e., as a number of sites. The starting point of the analysis is the drainage network of a fluvial equilibrium topography on a square L×L grid with L=10 000. Boundary conditions and parameter values except for the grid size are the same as those in the smaller examples shown in Fig. 1.

Figure 5 reveals that the eroded area Ae increases with the fluvial threshold Ac but becomes independent of A if the catchment size A is sufficiently large. This means that the hillslopes draining into large rivers are not systematically larger than those draining into small rivers. It is the reason why we will arrive at a scaling relation that preserves the form of Eq. (2) and avoids the problem occurring if the river width is used for scaling.

Figure 5Eroded area Ae as a function of the catchment size A for different fluvial thresholds Ac. Raw data were used for the catchment sizes that occurred at least 1000 times on the grid. Otherwise, data were binned dynamically so that there are at least 1000 points in each bin.


The increase in Ae if A approaches Ac can be explained by distinguishing between river segments and channel heads. Let us define channel heads as those sites without any tributary with AAc, i.e., as those sites that are only supplied by hillslopes. All other sites with AAc are considered to be river segments. All sites with A=Ac are channel heads and thus follow the relation Ae=A so that all curves start at the dotted line in Fig. 5. The resulting values Ae of the river segments (without the channel heads) are shown by the dashed lines in Fig. 5. The increase in Ae if A approaches Ac even turns into a decrease then. This decrease, which arises from the limitation AeA-Ac, holds for all river segments that have at least one tributary cell contributing at least Ac. Thus the contribution of the hillslopes must be small if A is only slightly larger than Ac. However, the decrease is exaggerated by the logarithmic scale and concerns only a small number of sites, so it makes sense to assume that Ae is independent of A for river segments.

Both the number of river segment sites and the number of channel head sites decrease with an increasing threshold Ac. The decrease in the latter is faster so that the ratio of the numbers of head sites to river sites converges to zero for large Ac values. This is, however, not true for the total contributions. Figure 6 shows the ratio of the sum of the Ae values of all river segments to the sum of the Ae values of the channel heads. It can also be interpreted as the ratio of the total area that must be eroded by the river segments to the total area that must be eroded by the channel heads. The results shown for different grid sizes shown in Fig. 6 suggest that this ratio becomes constant in the limit of large grid sizes. It apparently approaches a value of about 2 here, which means that the river segments contribute about two-thirds, and the channel heads one-third, to total fluvial erosion.

Figure 6Ratio of total area eroded by all river segments to total area eroded by all channel head sites as a function of the fluvial threshold Ac.


This result suggests that the dependency of Ae on the threshold Ac is determined by the cumulative distribution P(A) of the catchment sizes in the drainage network. This distribution describes the probability that a randomly selected site has a catchment size A. The probability P(Ac) evaluated at the fluvial threshold is the ratio of the area covered by all channel pixels to the total area. It can be interpreted as a drainage density (river length per total area) on a discrete grid. Then a fraction P(Ac) of the considered domain must erode a given fraction (here about two-thirds) of the domain, leading to the relation

(12) A e = γ P ( A c ) ,

with γ23 for this network. While Ae can be measured directly for the considered drainage network, its relation to P(A) (Eq. 12) is useful, as this distribution has already been investigated in several studies on natural and modeled drainage networks (Rodriguez-Iturbe et al.1992a; Maritan et al.1996b; Rodriguez-Iturbe and Rinaldo1997; Rinaldo et al.1998; Hergarten and Neugebauer2001; Hergarten2002; Hergarten et al.2014, 2016). It was found that P(A) follows a power-law distribution,

(13) P ( A ) A - β ,

over a reasonable range, where a range β[0.41,0.46] was found except for the two latest studies. In these studies, larger networks were considered to be making use of increasing data availability and computing capacities. An exponent very close to 0.5 was found for both optimal channel networks (OCNs; see below) (Hergarten et al.2014) and a real river pattern at the continental scale (Hergarten et al.2016).

Equations (12) and (13) suggest a power-law relation,

(14) A e = α A c β ,

between the eroded area and the fluvial threshold. The validity of Eqs. (12), (13), and (14) is investigated in Fig. 7. Comparing the two solid curves reveals that Eq. (12) does not hold exactly, since the curves come closer to each other for decreasing catchment sizes. The reason for this is that Ae only refers to the river segments without the channel heads so that P(Ac) in Eq. (12) should also exclude the channel head sites. The dashed colored line in Fig. 7 showing the accordingly reduced distribution P(A) illustrates that Eq. (12) indeed holds then and that the effect vanishes for large Ac values.

Figure 7Black axes: eroded area as a function of the fluvial threshold. Colored axes: cumulative distribution of the catchment sizes.


The black dashed line in Fig. 7 refers to the best-fit power-law relation according to Eq. (14). It is based on all integer values of Ac from 1 to 10 000, assuming equal errors, so that the large values of Ac practically have a high weight in the fit. The power law with the obtained values α=1.360 and β=0.465 fits the data well, with a relative error of less than 5 % for Ac[15,10000] and less than 1 % for Ac[400,10000]. The deviations are larger for smaller fluvial thresholds due to the fact that dendritic networks cannot be represented well on a regular lattice at small scales.

The relation to the catchment-size distribution (Eqs. 12 and 13) suggests that the power-law dependency of Ae on Ac (Eq. 14) should be universal. For testing this hypothesis, a set of equilibrium topographies with θ{0.25,0.45,0.5,0.75} was analyzed. These values cover the range that has been found so far under relatively homogeneous conditions (e.g., Robl et al.2017). The value θ=0.45 was added, as it is often used as a reference value instead of θ=0.5 (e.g., Whipple et al.2013; Lague2014). Parameter values and boundary conditions are the same as in the previous example. Since the exponent n has no immediate effect on equilibrium topographies, values n≠1 were not considered.

The power-law parameters α and β obtained from equilibrium topographies on different lattice sizes L are given in Table 1. In addition, the original data for the largest grids are shown in Fig. 8. The results are overall similar, with a tendency to lower exponents β for increasing θ. A notable deviation is only found for the very high concavity index θ=0.75. Here the slopes become very steep at small catchment sizes, resulting in a slower migration of drainage divides during the simulation (Robl et al.2017). As a result, the topography reaches a steady state quite soon so that there is finally less reorganization in the drainage network with regard to the initial random pattern. In this sense, the lower exponents found for θ=0.75 can be seen as a fingerprint of poorly organized river patterns but are probably not relevant for the rivers that were the empirical basis of the stream-power law. These findings confirm that the concavity index θ has a minor effect on the topology of the drainage networks, although it strongly affects the shape of longitudinal river profiles and thus the topography.

Table 1Parameter values of the power-law relation between eroded area and fluvial threshold (Eq. 14) obtained from different simulated drainage networks on regular lattices with L×L nodes.

Download Print Version | Download XLSX

Figure 8Eroded area Ae as a function of the fluvial threshold Ac for the considered drainage networks. For clarity, only the results obtained from the largest domains are plotted.


In addition, Table 1 and Fig. 8 also contain results obtained from optimal channel networks (OCNs) on a grid with L=4096. Optimal channel networks are derived from the principle of minimum energy dissipation and have been widely used in the context of river networks (e.g., Howard1990; Rodriguez-Iturbe et al.1992c, b; Rinaldo et al.1992, 1998; Maritan et al.1996a, b). The networks considered here are those shown in Fig. 1 of Hergarten et al. (2014), where θ is related to the parameter n used there by θ=n-1n+1. The values of Ae of OCNs are overall slightly higher than those of the equilibrium topographies, and the variation with θ is lower. As OCNs are organized more strongly than drainage networks of arbitrary equilibrium topographies, the lower variability among OCNs is not surprising.

Table 2 provides additional results obtained from steady-state topographies on triangulated irregular networks (TINs). Numbers of neighbors, distances to neighbors, and areas of pixels are variable here. Areas of pixels are defined by the Voronoi diagram. Nondimensional areas (in DEM pixels) are normalized to the mean pixel size given by δ2=AtotN, where Atot is total area and N the number of nodes. The values listed in Table 2 and the respective curve in Fig. 8 show that the results obtained from TINs are close to those obtained from regular meshes.

Table 2Parameter values of the power-law relation between eroded area and fluvial threshold (Eq. 14) obtained from different simulated drainage networks on triangular lattices with N nodes for θ=0.5.

Download Print Version | Download XLSX

These results suggest defining the values α=1.508 and β=0.478 obtained from the OCN with θ=0.5 as reference values. The question is, however, whether such precision is useful for applications. In particular, β=0.5 would be more convenient than lower values. In the considerations made above, all areas are measured in DEM pixels and are thus nondimensional properties. Considering Ac to be a physical (dimensional) area, Ac has to be replaced by Acδ2 in Eq. (14). Then the fluvial erosion rate (Eq. 11) turns into

(15) E = α A c δ 2 β K A m S n ,

so that the fluvial incision term scales like δ−2β. For β=0.5, the fluvial term scales like 1δ. This is not only convenient but also leads to basically the same scaling relation assumed by Perron et al. (2008). The only difference is that the term αAc occurring here was interpreted as a channel width w and then assumed to be constant for all rivers so that it lost its physical meaning. Thus the new formulation of the fluvial incision term also fixes the concern raised by Pelletier (2010) that led to the alternative formulation where the hillslope transport term was rescaled.

In order to estimate α for β=0.5, it is helpful to know which region of Fig. 8 is occupied by typical model applications. A breakdown of Flint's law (Eq. 3) was reported at catchment sizes between between about 0.1 and 5 km2 (Montgomery and Foufoula-Georgiou1993; Stock and Dietrich2003; Wobus et al.2006). However, channel steepness declines at small catchment sizes, so this breakdown implies that other erosion processes come into play rather than that fluvial erosion is no longer active. In turn, many small springs in mountain regions have discharges on the order of magnitude of 0.1 L s−1 (e.g., Hergarten et al.2016), corresponding to catchment sizes A<0.01 km2, but it is not clear whether the erosive action of the resulting small streams follows Flint's law. Reasonable estimates of Ac are probably between these two ranges. Assuming a spatial resolution of about 100 m or a bit less, Ac will be on the order of magnitude of a few to 100 DEM pixels. As illustrated by the black line in Fig. 8, α=2 provides a reasonable estimate for this range with simple numbers as αAcβ=2Ac. With this estimate, the scaling factor for the fluvial erosion rate is 2Acδ, and the modified stream-power law for fluvial erosion turns into

(16) E = 2 A c δ K A m S n .
4 Numerical examples

Let us first return to the example of parallel rivers considered in Fig. 3. It was found in Sect. 2 that the topography of the hillslopes was robust against the spatial resolution, while the channel slope increases with decreasing grid spacing δ. Both approaches previously published fix this problem, but the channel slopes are too steep by a factor of dw compared to what is expected from the erodibility.

It should be noted that this example is not related to the approach to estimate Ae from Ac for dendritic networks (Eqs. 15 and 16) but can only test the validity of the principal scaling approach (Eq. 11). The size of the area Ae does not follow Eq. (14) but is defined by the geometry as Ae=dδ (measured in DEM pixels). Figure 9 shows the numerical results for the parameter values used in Fig. 3 for different values of δ. The simulation was started from a flat topography where the flow paths of the parallel rivers are predefined. As the problem is linear for n=1, this example can also be seen as the change in the river profile over time if uplift suddenly increases at t=0, while the base level remains constant. The results show that the equilibrium profile achieved for long times is reproduced correctly and that the time-dependent behavior is also robust against the resolution. This means that the scaling approach itself (Eq. 11) yields both the correct equilibrium behavior and the correct timescale.

Figure 9Numerical results for the scenario considered in Fig. 3. The river profiles obtained for δ=0.025 and δ=0.01 cannot be distinguished visually.


The second example refers to the scenario considered in Fig. 1 but extended by a fluvial threshold Ac=10-5 and by linear diffusion with a diffusivity D=10-5. The threshold Ac is a property of the fluvial erosion process, while the diffusive hillslope process is not related to it. It is thus assumed that fluvial erosion acts only at sites where AAc, while diffusion is active everywhere. A TIN representation is used in order to avoid artifacts from the combination of the eight-neighbor (D8) flow-routing scheme with the standard four-neighbor diffusion scheme on a regular mesh. The simulations are started from an almost flat topography with unit uplift. Uplift is switched off at t=50 in order to observe the decay of the topography.

The mean steepness index ks of the large rivers is plotted as a function of time in Fig. 10. Large rivers are defined by A10-3 here, which is considerably larger than Ac but much smaller than the domain. As expected, the simulations performed without any rescaling of the erodibility (dashed lines) are strongly affected by the spatial resolution. The steepness index increases with an increasing number of nodes N, i.e., with decreasing pixel size. In turn, the results obtained using the simple scaling relation (Eq. 16; solid lines) have a much weaker dependence on resolution. There is, however, a residual variation in channel steepness. The mean value of ks varies between about 1.6 and 2.0 over the considered range from N=105 to N=107. This result does not change fundamentally if a higher or lower threshold than A10-3 is used for defining large rivers.

Figure 10Mean steepness index ks of the large rivers obtained from simulations on TINs with different resolutions, defined by the total number of nodes N. Solid lines refer to the simplified scaling approach suggested in this paper (Eq. 16), while dashed lines refer to simulations performed without any rescaling. The latter are plotted only for N≤106.


5 Discussion

It may be surprising that the example of fluvial incision and hillslope diffusion considered in the previous section yields a mean steepness index greater than 1, although the scaling concept was developed in order to preserve channel steepness. The concept is, however, based on a generic hillslope process where the direction of transport follows a hypothetic fluvial equilibrium pattern and turns into fluvial erosion at a given threshold catchment size Ac. It is questionable whether any hillslope process occurring in nature comes close to this simple model. In the example considered here, the diffusion process is characterized by a diffusivity D and is not related to Ac. The fluvial domain is affected by diffusion more and more with increasing diffusivity. As a consequence, slopes of small channels decrease so that they erode less efficiently. This has to be compensated by the larger rivers so that they become steeper.

This is, however, a real property of the hillslope process here, and it is not the goal of the scaling approach to remove it. The concept presented here aims at removing the dependence on the resolution and providing the way in which values of the erodibility should be interpreted. Here it is suggested that they should be considered in combination with a fluvial threshold Ac in such a way that they would yield the expected channel steepness if the generic hillslope model were valid.

In turn, the residual dependence of channel steepness on resolution is a problem, in particular because it is not clear whether it converges in the limit δ→0 (N→∞). The problem arises from network reorganization, which also affects the fluvial region. Diffusion disturbs the dendritic topology towards parallel flow where the model based on Hack's findings (Eq. 2) is not valid. Using an improved flow-routing scheme that is able to distinguish channelized flow from parallel flow as suggested by Pelletier (2010) and letting Ac self-adjust might reduce the problem. However, the aim of this study is to develop a simple, quite universal rescaling approach that avoids or at least reduces the dependence on resolution without modifying the applied model seriously. In this sense, Eq. (16) should be a good trade-off.

Nevertheless it is important to keep the difference between detachment-limited erosion and pure bedrock incision in mind. Here it is assumed that the ability of the river to take up particles and carry them away concerns both the riverbed and material coming from adjacent hillslopes. If we, conversely, assume that all material coming from the hillslopes is instantaneously removed by the river without any consequences, there is no feedback of the hillslopes to the rivers, and Eq. (1) does not require any rescaling.

The results of this study have consequences for scaling relations in coupled models of rivers and hillslopes. Theodoratos et al. (2018) conducted a comprehensive analysis of the problem with linear diffusion without rescaling. The parameters they used were the same as in the previous example (Fig. 10), so it is immediately clear that their numerical results strongly depend on resolution. The authors argued that, following the approach of Pelletier (2010), both grid spacing and channel width are rescaled so that the ratio δw remains constant, and the scaling issue is consistent throughout all scales. However, the results presented here show that the property relevant for compensating δ is not channel width but Ae and thus Ac. These parameters are, however, physical properties of the erosion process, so they do not scale with the size of the domain. As a consequence, the characteristic horizontal length scale of the coupled system should rather be

(17) l c = D A c K ,

for m=0.5 and n=1 instead of lc=DK used by Theodoratos et al. (2018). This problem also affects the recent extension by an erosion threshold (Theodoratos and Kirchner2020).

6 Conclusions

This study presents a simple scaling relation for the fluvial incision term in landform evolution models involving detachment-limited fluvial erosion and hillslope processes. In order to avoid a dependence of the simulated topographies on the spatial resolution of the grid, the fluvial incision term must be multiplied by a scaling factor depending on the ratio of the threshold catchment size Ac where fluvial erosion starts and the pixel size δ2 of the grid. The analysis of several simulated drainage networks yields a power-law dependence of the scaling factor in Eq. (15) with an exponent slightly lower than 0.5. However, for application in numerical models, a simpler approximation where the fluvial erosion rate is rescaled by a factor 2Acδ is suggested. As this relation assumes a simple, generic hillslope process, it cannot provide an exact solution for all types of hillslope processes. In combination with such processes, e.g., diffusion, the dependence on the spatial resolution is not completely removed. Nevertheless, the simple scaling relation appears to be a reasonable trade-off between accuracy and simplicity.

Code and data availability

All codes and computed data can be downloaded from the FreiDok data repository 155182 (Hergarten2020). The author is happy to support interested readers in reproducing the results and performing subsequent research.

Competing interests

The author declares that there is no conflict of interest.


The author would like to thank the two anonymous reviewers for their thorough consideration and for their very constructive suggestions to improve the readability of the paper. The author would also like to thank Wolfgang Schwanghart for the editorial handling.

Review statement

This paper was edited by Wolfgang Schwanghart and reviewed by Taylor Perron and two anonymous referees.


Andrews, D. J. and Bucknam, R. C.: Fitting degradation of shoreline scarps by a nonlinear diffusion model, J. Geophys. Res., 92, 12857–12867,, 1987. 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,, 2017. a

Culling, W.: Analytical theory of erosion, J. Geol., 68, 336–344,, 1960. a

Duvall, A. R. and Tucker, G. E.: Dynamic ridges and valleys in a strike-slip environment, J. Geophys. Res.-Earth, 120, 2016–2026,, 2015. a

Ferrier, K. L., Perron, J. T., Mukhopadhyay, S., Rosener, M., Stock, J. D., Huppert, K. L., and Slosberg, M.: Covariation of climate and long-term erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua´i, GSA Bull., 125, 1146–1163,, 2013. a

Flint, J. J.: Stream gradient as a function of order, magnitude, and discharge, Water Resour. Res., 10, 969–973,, 1974. a

Goren, L., Willett, S. D., Herman, F., and Braun, J.: Coupled numerical–analytical approach to landscape evolution modeling, Earth Surf. Proc. Land., 39, 522–545,, 2014. a

Gray, H. J., Shobe, C. M., Hobley, D. E. J., Tucker, G. E., Duvall, A. R., Harbert, S. A., and Owen, L. A.: Off-fault deformation rate along the southern San Andreas fault at Mecca Hills, southern California, inferred from landscape modeling of curved drainages, Geology, 46, 59–62,, 2018. 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.,, 1957. 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,, 2016. a

Hergarten, S.: Self-Organized Criticality in Earth Systems, Springer, Berlin, Heidelberg, New York,, 2002. a, b

Hergarten, S.: Rivers as linear elements in landform evolution models: codes and data, FreiDok plus, Universitätsbibliothek Freiburg,, 2020. a

Hergarten, S. and Neugebauer, H. J.: Self-organized critical drainage networks, Phys. Rev. Lett., 86, 2689–2692,, 2001. a, b

Hergarten, S., Winkler, G., and Birk, S.: Transferring the concept of minimum energy dissipation from river networks to subsurface flow patterns, Hydrol. Earth Syst. Sci., 18, 4277–4288,, 2014. a, b, c

Hergarten, S., Winkler, G., and Birk, S.: Scale invariance of subsurface flow patterns and its limitation, Water Resour. Res., 52, 3881–3887,, 2016. a, b, c

Howard, A. D.: Theoretical model of optimal drainage networks, Water Resour. Res., 26, 2107–2117,, 1990. a

Howard, A. D.: A detachment-limited model for drainage basin evolution, Water Resour. Res., 30, 2261–2285,, 1994. a, b, c, d

Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Proc. Land., 39, 38–61,, 2014. a

Maritan, A., Colaiori, F., Flammini, A., Cieplak, M., and Banavar, J. R.: Universality classes of optimal channel networks, Science, 272, 984–986,, 1996a. a

Maritan, A., Rinaldo, A., Rigon, R., Giacometti, A., and Rodriguez-Iturbe, I.: Scaling laws for river networks, Phys. Rev. E, 53, 1510–1515,, 1996b. a, b

Montgomery, D. R. and Foufoula-Georgiou, E.: Channel network source representation using digital elevation models, Water Resour. Res., 29, 3925–3934,, 1993. a

O'Callaghan, J. F. and Mark, D. M.: The extraction of drainage networks from digital elevation data, Comput. Vision Graph., 28, 323–344,, 1984. a

Pelletier, J. D.: Minimizing the grid-resolution dependence of flow-routing algorithms for geomorphic applications, Geomorphology, 122, 91–98,, 2010. a, b, c, d, e

Perron, J. T., Dietrich, W. E., and Kirchner, J. W.: Controls on the spacing of first-order valleys, J. Geophys. Res.-Earth, 113, F04016,, 2008. a, b, c, d, e

Reitman, N. G., Mueller, K. J., Tucker, G. E., Gold, R. D., Briggs, R. W., and Barnhart, K. R.: Offset channels may not accurately record strike-slip fault displacement: Evidence from landscape evolution models, J. Geophys. Res.-Sol. Ea., 124, 13427–13451,, 2019. a

Rinaldo, A., Rodriguez-Iturbe, I., Bras, R. L., Ijjasz-Vasquez, E., and Marani, A.: Minimum energy and fractal structures of drainage networks, Water Resour. Res., 28, 2181–2195,, 1992. a

Rinaldo, A., Rodriguez-Iturbe, I., and Rigon, R.: Channel networks, Annu. Rev. Earth Pl. Sc., 26, 289–327,, 1998. a, b

Robl, J., Hergarten, S., and Prasicek, G.: The topographic state of fluvially conditioned mountain ranges, Earth Sci. Rev., 168, 290–317,, 2017. a, b, c

Rodriguez-Iturbe, I. and Rinaldo, A.: Fractal River Basins. Chance and Self-Organization, Cambridge University Press, Cambridge, UK, New York, USA, Melbourne, Australia, 1997. a

Rodriguez-Iturbe, I., Ijjasz-Vasquez, E., Bras, R. L., and Tarboton, D. G.: Power law distribution of mass and energy in river basins, Water Resour. Res., 28, 1089–1093,, 1992a. a

Rodriguez-Iturbe, I., Rinaldo, A., Rigon, R., Bras, R. L., Ijjasz-Vasquez, E., and Marani, A.: Fractal structures as least energy patterns: The case of river networks, Geophys. Res. Lett., 19, 889–892,, 1992b. a

Rodriguez-Iturbe, I., Rinaldo, A., Rigon, R., Bras, R. L., Marani, A., and Ijjasz-Vasquez, E.: Energy dissipation, runoff production, and the three-dimensional structure of river basins, Water Resour. Res., 28, 1095–1103,, 1992c. a

Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resour. Res., 35, 853–870,, 1999. a

Shobe, C. M., Tucker, G. E., and Anderson, R. S.: Hillslope-derived blocks retard river incision, Geophys. Res. Lett., 43, 5070–5078,, 2016. a

Stock, J. and Dietrich, W. E.: Valley incision by debris flows: Evidence of a topographic signature, Water Resour. Res., 39, 1089,, 2003. a

Theodoratos, N. and Kirchner, J. W.: Dimensional analysis of a landscape evolution model with incision threshold, Earth Surf. Dynam. Discuss.,, in review, 2020. a

Theodoratos, N., Seybold, H., and Kirchner, J. W.: Scaling and similarity of a stream-power incision and linear diffusion landscape evolution model, Earth Surf. Dynam., 6, 779–808,, 2018. a, b, c

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,, 1999. 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 Treatise on Geomorphology, Academic Press, San Diego, CA, USA, 550–573,, 2013. a

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., USA, 55–74,, 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 Mars-like rampart crater, Earth Planet. Sc. Lett., 506, 209–220,, 2019. a, b

Short summary
Models of fluvial erosion have a long history in landform evolution modeling. Interactions between rivers and processes acting at hillslopes (e.g., landslides) are receiving growing interest in this context. While present-day computer capacities allow for applying such coupled models, there is still a scaling problem when considering rivers to be linear elements on a topography. Based on a reinterpretation of old empirical results, this study presents a new approach to overcome this problem.