**Research article**
26 May 2020

**Research article** | 26 May 2020

# Rivers as linear elements in landform evolution models

Stefan Hergarten

**Stefan Hergarten**Stefan Hergarten

- Institut für Geo- und Umweltnaturwissenschaften, Albertstr. 23B, 79104 Freiburg, Germany

- Institut für Geo- und Umweltnaturwissenschaften, Albertstr. 23B, 79104 Freiburg, Germany

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

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

Received: 18 Dec 2019 – Discussion started: 17 Jan 2020 – Revised: 06 Apr 2020 – Accepted: 23 Apr 2020 – Published: 26 May 2020

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.

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({x}_{\mathrm{1}},{x}_{\mathrm{2}},t)$ of a landform evolution model with detachment-limited fluvial erosion reads

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

**that take the nonlinear dependencies of hillslope processes on topography into account (e.g., Andrews and Bucknam, 1987; Howard, 1994; Roering et al., 1999).**

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

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 Tucker, 1999). 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 (Flint, 1974). This relationship is nowadays usually written in the form

where *θ* is the concavity index and *k*_{s} the steepness index.
Assuming that Eq. (3) is the fingerprint of spatially uniform
steady-state conditions, it predicts $\frac{m}{n}=\mathit{\theta}$ 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

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.

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 Mark, 1984) and a fully
implicit scheme (Hergarten and Neugebauer, 2001; Hergarten, 2002) 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.

Relief increases with decreasing grid spacing because the smallest
catchment size that can be resolved is ${A}_{min}={\mathit{\delta}}^{\mathrm{2}}$, and the maximum equilibrium
slope is proportional to ${A}_{min}^{-\mathit{\theta}}={\mathit{\delta}}^{-\mathrm{2}\mathit{\theta}}$ 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
*x*_{2} direction) in equilibrium with constant uplift. Linear diffusion,

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*=10^{6} was assumed for each river segment so that the
channel slope should theoretically be $S={\mathrm{10}}^{-\mathrm{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.

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

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

so the fluvial erosion rate required for compensating uplift is higher than it would be without local transport by a factor $\frac{d}{\mathit{\delta}}$. This requires an increase in the channel slope by a factor of ${\left(\frac{d}{\mathit{\delta}}\right)}^{\frac{\mathrm{1}}{n}}$ 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 $\frac{w}{\mathit{\delta}}$, 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 $\frac{w}{\mathit{\delta}}$
yields

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
div** q** by the inverse factor $\frac{\mathit{\delta}}{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

instead of Eq. (6), so that

For *w*≪*d* 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 Tucker, 2015; 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 *A*_{c}, 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.

The simple example considered in the previous section involves a dependence on
grid spacing *δ* according to the factor $\frac{d}{\mathit{\delta}}$ without rescaling
(Eq. 7). Both approaches for rescaling replace the dependence
on *δ* by a dependence on the channel width *w* so that a factor
$\frac{d}{w}$ remains (Eq. 8). This is, however, still a problem
if *w* is not constant. The occurrence of the factor $\frac{d}{w}$ 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 *A*_{c} in such a way that all sites with
*A*≥*A*_{c} are river segments,
while all sites with *A*<*A*_{c} 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.

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

where *A*_{e} 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 *A*_{e} in Eq. (11) is the mean
size, channel steepness will just vary randomly, which is also found
in nature. A systematic dependence of *A*_{e} 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 *A*_{e} depends on *A* and on *A*_{c}.
More precisely, *A*_{e} is the mean size of all hillslope areas
draining into channel sites with a given catchment size *A* at a given fluvial
threshold *A*_{c} (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 *A*_{e} increases
with the fluvial threshold *A*_{c} 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.

The increase in *A*_{e} if *A* approaches *A*_{c} can be explained by distinguishing
between river segments and channel heads. Let us define channel
heads as those sites without any tributary with *A*≥*A*_{c}, i.e.,
as those sites that are only supplied by hillslopes. All other sites
with *A*≥*A*_{c} are considered to be river segments.
All sites with *A*=*A*_{c} are channel heads and thus follow the
relation *A*_{e}=*A* so that all curves start at the dotted line in Fig. 5.
The resulting values *A*_{e} of the river segments (without the channel heads)
are shown by the dashed lines in Fig. 5. The increase in *A*_{e}
if *A* approaches *A*_{c} even turns into a decrease then. This decrease, which arises from the limitation ${A}_{\mathrm{e}}\le A-{A}_{\mathrm{c}}$, holds for all river segments that
have at least one tributary cell contributing at least *A*_{c}. Thus the
contribution of the hillslopes must be small if *A* is only slightly larger
than *A*_{c}. However, the decrease is exaggerated by the logarithmic scale
and concerns only a small number of sites, so it makes sense
to assume that *A*_{e} 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 *A*_{c}. 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
*A*_{c} values. This is, however, not true for the total contributions. Figure 6
shows the ratio of the sum of the *A*_{e} values of all river segments to the sum of
the *A*_{e} 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.

This result suggests that the dependency of *A*_{e} on the threshold *A*_{c} 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*(*A*_{c}) 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*(*A*_{c}) of the considered domain
must erode a given fraction (here about two-thirds)
of the domain, leading to the relation

with $\mathit{\gamma}\approx \frac{\mathrm{2}}{\mathrm{3}}$ for this network.
While *A*_{e} 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 Rinaldo, 1997; Rinaldo et al., 1998; Hergarten and Neugebauer, 2001; Hergarten, 2002; Hergarten et al., 2014, 2016).
It was found that *P*(*A*) follows a power-law distribution,

over a reasonable range, where a range $\mathit{\beta}\in [\mathrm{0.41},\mathrm{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,

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 *A*_{e} only refers to the river segments without the channel
heads so that *P*(*A*_{c}) 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 *A*_{c} values.

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 *A*_{c} from
1 to 10 000, assuming equal errors, so that the large values of *A*_{c} 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 ${A}_{\mathrm{c}}\in [\mathrm{15},\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{000}]$ and less than 1 % for ${A}_{\mathrm{c}}\in [\mathrm{400},\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{000}]$.
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 *A*_{e} on *A*_{c} (Eq. 14) should
be universal. For testing this hypothesis, a set of equilibrium topographies with
$\mathit{\theta}\in \mathit{\{}\mathrm{0.25},\mathrm{0.45},\mathrm{0.5},\mathrm{0.75}\mathit{\}}$ 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; Lague, 2014).
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.

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., Howard, 1990; 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 $\mathit{\theta}=\frac{n-\mathrm{1}}{n+\mathrm{1}}$.
The values of *A*_{e} 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 ${\mathit{\delta}}^{\mathrm{2}}=\frac{{A}_{\mathrm{tot}}}{N}$, where *A*_{tot} 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.

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 *A*_{c} to be a physical (dimensional) area, *A*_{c} has to be
replaced by $\frac{{A}_{\mathrm{c}}}{{\mathit{\delta}}^{\mathrm{2}}}$ in Eq. (14).
Then the fluvial erosion rate (Eq. 11) turns into

so that the fluvial incision term scales like *δ*^{−2β}. For *β*=0.5, the
fluvial term scales like $\frac{\mathrm{1}}{\mathit{\delta}}$. 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 $\mathit{\alpha}\sqrt{{A}_{\mathrm{c}}}$ 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 km^{2} (Montgomery and Foufoula-Georgiou, 1993; Stock and Dietrich, 2003; 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 km^{2}, but it is not clear whether
the erosive action of the resulting small streams follows Flint's law. Reasonable estimates
of *A*_{c} are probably between these two ranges. Assuming a spatial resolution of about
100 m or a bit less, *A*_{c} will be on the order of magnitude of a few to 100 DEM
pixels. As illustrated by the black line in Fig. 8, $\mathit{\alpha}=\sqrt{\mathrm{2}}$
provides a reasonable estimate for this range with simple numbers as
$\mathit{\alpha}{A}_{\mathrm{c}}^{\mathit{\beta}}=\sqrt{\mathrm{2}{A}_{\mathrm{c}}}$. With this estimate,
the scaling factor for the fluvial erosion rate is $\frac{\sqrt{\mathrm{2}{A}_{\mathrm{c}}}}{\mathit{\delta}}$, and
the modified stream-power law for fluvial erosion turns into

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 $\frac{d}{w}$ compared to what is expected from
the erodibility.

It should be noted that this example is not related to the approach to estimate *A*_{e}
from *A*_{c} 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 *A*_{e} does not follow Eq. (14) but is defined by
the geometry as ${A}_{\mathrm{e}}=\frac{d}{\mathit{\delta}}$ (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.

The second example refers to the scenario considered in Fig. 1 but extended
by a fluvial threshold ${A}_{\mathrm{c}}={\mathrm{10}}^{-\mathrm{5}}$ and by linear diffusion
with a diffusivity $D={\mathrm{10}}^{-\mathrm{5}}$. The threshold *A*_{c}
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 *A*≥*A*_{c}, 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 *k*_{s} of the large rivers is plotted as a function of time in
Fig. 10. Large rivers are defined by $A\ge {\mathrm{10}}^{-\mathrm{3}}$ here, which is considerably
larger than *A*_{c} 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 *k*_{s} varies between about 1.6 and 2.0
over the considered range from *N*=10^{5} to *N*=10^{7}. This result does not change
fundamentally if a higher or lower threshold than $A\ge {\mathrm{10}}^{-\mathrm{3}}$ is used for defining
large rivers.

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 *A*_{c}. 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 *A*_{c}. 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 *A*_{c} 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 *A*_{c} 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 $\frac{\mathit{\delta}}{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 *A*_{e} and thus *A*_{c}. 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

for *m*=0.5 and *n*=1 instead of ${l}_{\mathrm{c}}=\sqrt{\frac{D}{K}}$
used by Theodoratos et al. (2018).
This problem also affects the recent extension by an erosion threshold (Theodoratos and Kirchner, 2020).

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
*A*_{c} 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 $\frac{\sqrt{\mathrm{2}{A}_{\mathrm{c}}}}{\mathit{\delta}}$ 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.

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.

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, https://doi.org/10.1029/JB092iB12p12857, 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, https://doi.org/10.5194/esurf-5-47-2017, 2017. a

Culling, W.: Analytical theory of erosion, J. Geol., 68, 336–344, https://doi.org/10.1086/626663, 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, https://doi.org/10.1002/2015JF003618, 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, 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., Willett, S. D., Herman, F., and Braun, J.: Coupled numerical–analytical approach to landscape evolution modeling, Earth Surf. Proc. Land., 39, 522–545, https://doi.org/10.1002/esp.3514, 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, https://doi.org/10.1130/G39820.1, 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., https://doi.org/10.3133/pp294B, 1957. a

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

Hergarten, S.: Self-Organized Criticality in Earth Systems, Springer, Berlin, Heidelberg, New York, https://doi.org/10.1007/978-3-662-04390-5, 2002. a, b

Hergarten, S.: Rivers as linear elements in landform evolution models: codes and data, FreiDok plus, Universitätsbibliothek Freiburg, https://doi.org/10.6094/UNIFR/155182, 2020. a

Hergarten, S. and Neugebauer, H. J.: Self-organized critical drainage networks, Phys. Rev. Lett., 86, 2689–2692, https://doi.org/10.1103/PhysRevLett.86.2689, 2001. a, 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, https://doi.org/10.5194/hess-18-4277-2014, 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, https://doi.org/10.1002/2015WR017530, 2016. a, b, c

Howard, A. D.: Theoretical model of optimal drainage networks, Water Resour. Res., 26, 2107–2117, https://doi.org/10.1029/WR026i009p02107, 1990. a

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

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

Maritan, A., Colaiori, F., Flammini, A., Cieplak, M., and Banavar, J. R.: Universality classes of optimal channel networks, Science, 272, 984–986, https://doi.org/10.1126/science.272.5264.984, 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, https://doi.org/10.1103/PhysRevE.53.1510, 1996b. a, b

Montgomery, D. R. and Foufoula-Georgiou, E.: Channel network source representation using digital elevation models, Water Resour. Res., 29, 3925–3934, https://doi.org/10.1029/93WR02463, 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, https://doi.org/10.1016/S0734-189X(84)80011-0, 1984. a

Pelletier, J. D.: Minimizing the grid-resolution dependence of flow-routing algorithms for geomorphic applications, Geomorphology, 122, 91–98, https://doi.org/10.1016/j.geomorph.2010.06.001, 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, https://doi.org/10.1029/2007JF000977, 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, https://doi.org/10.1029/2019JB018596, 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, https://doi.org/10.1029/92WR00801, 1992. a

Rinaldo, A., Rodriguez-Iturbe, I., and Rigon, R.: Channel networks, Annu. Rev. Earth Pl. Sc., 26, 289–327, https://doi.org/10.1146/annurev.earth.26.1.289, 1998. a, b

Robl, J., Hergarten, S., and Prasicek, G.: The topographic state of fluvially conditioned mountain ranges, Earth Sci. Rev., 168, 290–317, https://doi.org/10.1016/j.earscirev.2017.03.007, 2017. a, b, 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, https://doi.org/10.1029/91WR03033, 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, https://doi.org/10.1029/92GL00938, 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, https://doi.org/10.1029/91WR03034, 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, https://doi.org/10.1029/1998WR900090, 1999. a

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

Stock, J. and Dietrich, W. E.: Valley incision by debris flows: Evidence of a topographic signature, Water Resour. Res., 39, 1089, https://doi.org/10.1029/2001WR001057, 2003. a

Theodoratos, N. and Kirchner, J. W.: Dimensional analysis of a landscape evolution model with incision threshold, Earth Surf. Dynam. Discuss., https://doi.org/10.5194/esurf-2019-80, 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, https://doi.org/10.5194/esurf-6-779-2018, 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, https://doi.org/10.1029/1999JB900120, 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, https://doi.org/10.1016/B978-0-12-374739-6.00226-8, 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, 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 Mars-like rampart crater, Earth Planet. Sc. Lett., 506, 209–220, https://doi.org/10.1016/j.epsl.2018.11.009, 2019. a, b