**Research article**
04 Jun 2020

**Research article** | 04 Jun 2020

# Dimensional analysis of a landscape evolution model with incision threshold

Nikos Theodoratos and James W. Kirchner

^{1}

^{1,2}

**Nikos Theodoratos and James W. Kirchner**Nikos Theodoratos and James W. Kirchner

^{1}

^{1,2}

^{1}Department of Environmental Systems Science, ETH Zurich, Zurich, 8092, Switzerland^{2}Swiss Federal Research Institute WSL, Birmensdorf, 8903, Switzerland

^{1}Department of Environmental Systems Science, ETH Zurich, Zurich, 8092, Switzerland^{2}Swiss Federal Research Institute WSL, Birmensdorf, 8903, Switzerland

**Correspondence**: Nikos Theodoratos (theodoratos@usys.ethz.ch)

**Correspondence**: Nikos Theodoratos (theodoratos@usys.ethz.ch)

Received: 23 Dec 2019 – Discussion started: 09 Jan 2020 – Revised: 27 Apr 2020 – Accepted: 05 May 2020 – Published: 04 Jun 2020

The ability of erosional processes to incise into a topographic surface can be limited by a threshold. Incision thresholds affect the topography of landscapes and their scaling properties and can introduce nonlinear relations between climate and erosion with notable implications for long-term landscape evolution. Despite their potential importance, incision thresholds are often omitted from the incision terms of landscape evolution models (LEMs) to simplify analyses. Here, we present theoretical and numerical results from a dimensional analysis of an LEM that includes terms for threshold-limited stream-power incision, linear diffusion, and uplift. The LEM is parameterized by four parameters (incision coefficient and incision threshold, diffusion coefficient, and uplift rate). The LEM's governing equation can be greatly simplified by recasting it in a dimensionless form that depends on only one dimensionless parameter, the
incision-threshold number *N*_{θ}. This dimensionless parameter is
defined in terms of the incision threshold, the incision coefficient, and
the uplift rate, and it quantifies the reduction in the rate of incision due
to the incision threshold relative to the uplift rate. Being the only
parameter in the dimensionless governing equation, *N*_{θ} is the only
parameter controlling the evolution of landscapes in this LEM. Thus,
landscapes with the same *N*_{θ} will evolve geometrically similarly,
provided that their boundary and initial conditions are normalized according
to appropriate scaling relationships, as we demonstrate using a numerical
experiment. In contrast, landscapes with different *N*_{θ} values will
be influenced to different degrees by their incision thresholds. Using
results from a second set of numerical simulations, each with a different
incision-threshold number, we qualitatively illustrate how the value of *N*_{θ} influences the topography, and we show that relief scales with the quantity *N*_{θ}+1 (except where the incision threshold
reduces the rate of incision to zero).

In the uppermost parts of drainage basins, water is not flowing over the ground surface or is flowing too weakly to incise into it. At least two kinds of limits must typically be overcome for erosion by flowing water to begin. First, sufficient drainage area must be accumulated for overland flow to start; second, this flow must exert sufficient shear stress on the surface to overcome the mechanical resistance of rocks or soils and thus mobilize sediment (e.g., Perron, 2017).

Channel-incision terms in landscape evolution models (LEMs) often capture
both of these limits by including an incision threshold below which no
incision occurs. For instance, if *τ* is the shear stress that water
exerts on the bed and *τ*_{θ} is a critical value of shear stress
(equivalently, *τ* and *τ*_{θ} could refer to stream power),
then the rate of incision is zero for *τ*≤*τ*_{θ}, and it can
be described by a term of the form *k*(*τ*−*τ*_{θ})^{α}, for *τ*>*τ*_{θ}, where *k* and *α* are constants (e.g., Howard, 1994). Including such incision terms in LEMs changes the topographic properties of the landscapes that are synthesized; for example, it leads to decreased drainage densities, more convex hillslopes, and steeper slopes (e.g., Howard, 1994; Tucker and Bras, 1998; Perron et al., 2008).

In addition, incision thresholds can have notable consequences on the relationship between climate and long-term incision rates as described, for example, by Snyder et al. (2003), Tucker (2004), Lague et al. (2005), Perron (2017), and Deal et al. (2018). Specifically, incision thresholds stop smaller events from eroding the surface. In many wet climates, the total annual streamflow is high, but small, frequent events tend to contribute most of this total; in contrast, in many dry climates, a larger fraction of the total annual streamflow tends to be contributed by rare but intense events (e.g., Rossi et al., 2016). Therefore, a sufficiently high incision threshold could render ineffective a larger fraction of the total precipitation in wetter climates than in drier climates. This behavior can lead to a nonlinear dependence of long-term erosion rates on average precipitation; it can even lead to the counterintuitive observation that, in some cases, larger average precipitation corresponds to smaller long-term erosion rates (e.g., DiBiase and Whipple, 2011).

Furthermore, incision thresholds can play a role in setting the smallest
scales of valley dissection, which are among the fundamental scales that
characterize landscapes. For instance, Horton (1945) suggested that valley
dissection stops because further dissection would lead to hillslopes that
are too short to yield flow that can erode the surface. Montgomery and
Dietrich (1992) found that thresholds of the topographic quantity
$A\left(\right|\mathrm{\nabla}z|{)}^{\mathrm{2}}$, where *A* is drainage area and $\left|\mathrm{\nabla}z\right|$ is slope, could define locations of both channel and valley heads, the former being associated with an incision threshold and the latter with the smallest scale of dissection. Perron et al. (2008) studied the spacing of valleys, a scaling property closely related to the smallest scale of dissection. They found that valley spacing is most strongly controlled by the competition between advective and diffusive processes, such as stream incision and soil creep, respectively. However, they found that incision thresholds also control valley spacing by modulating the competition between advection and diffusion.

In Theodoratos et al. (2018), we performed a scaling analysis of an incision–diffusion LEM that did not include an incision threshold. In the present study, we add an incision threshold to that LEM and examine how our analysis needs to be modified to account for this threshold. More specifically, in Theodoratos et al. (2018), we dimensionally analyzed an LEM that includes three parameters – an incision coefficient, a diffusion coefficient, and an uplift rate. For that analysis, we used three characteristic scales (of length, height, and time) that are defined in terms of the three parameters of the LEM. As we explained in detail in Theodoratos et al. (2018), because the characteristic scales depend on the model parameters and because there are three parameters and three characteristic scales, the LEM can be greatly simplified by being recast in a dimensionless form that has no parameters.

Adding an incision threshold to the LEM that we analyzed in Theodoratos et al. (2018) increases the number of its parameters to four (see Eq. 1 below). This leads to the question of whether the LEM with incision threshold can be dimensionally analyzed using the same three characteristic scales that we used to dimensionally analyze the LEM without incision threshold (Theodoratos et al., 2018). Here, we hypothesize that these three scales are reasonable choices even after adding an incision threshold to the LEM, and we test this hypothesis by applying these scales and examining the resulting rescaled equations.

## 2.1 Governing equation

We study an LEM described by the governing equation (e.g., Howard, 1994; Dietrich et al., 2003):

where *z* is elevation at a point with coordinates (*x*, *y*), and *t* is
time; *A*, $\left|\mathrm{\nabla}z\right|$, and ∇^{2}*z* are topographic properties of the point, specifically, drainage area, slope, and Laplacian curvature, respectively; *K* is the incision coefficient, *D* is the diffusion coefficient, *U* is the uplift rate, and *θ* is the incision threshold, a threshold value of the quantity $\sqrt{A}\left|\mathrm{\nabla}z\right|$. The four parameters *K*, *D*, *U*, and *θ* are all assumed to be constant in time and uniform across space. The dimensions of these variables, topographic properties, and parameters are discussed in the following subsection.

The stream-power incision term $K\left(\sqrt{A}\right|\mathrm{\nabla}z|-\mathit{\theta})$ describes the rate of incision by flowing water. It is a special case of the more general incision term $K\left({A}^{m}\right(\left|\mathrm{\nabla}z\right|{)}^{n}-\mathit{\theta})$, where *m* and *n*, the exponents of drainage area and slope, have values that depend on the
rate law that is assumed to describe incision (such as shear stress or
stream power; e.g., Dietrich et al., 2003). Here, we examine the simplified
case of *m*=0.5 and *n*=1 because it leads to results described by much
simpler formulas; however, these results are valid for generic exponents *m*
and *n* as well but with more complicated formulas (see Appendix A of
Theodoratos et al., 2018, for results pertaining to an LEM without incision
threshold and with generic exponents *m* and *n*). The linear diffusion term
*D*∇^{2}*z* describes the rate of erosion or infilling by hillslope sediment transport processes. Finally, the uplift term *U* gives the rate of tectonic uplift within the model domain (or, equivalently, base level fall at its boundary).

Equation (1) is defined piecewise on two subdomains. The first subdomain,
where $\sqrt{A}\left|\mathrm{\nabla}z\right|\le \mathit{\theta}$, corresponds to areas where the rate of incision is zero because it is fully suppressed by the incision threshold, and, thus, the landscape evolves under the influence of diffusion and uplift only. We refer to these areas as the zones of zero incision, because they tend to form zones along ridges and drainage divides, where drainage area or slope, or both, are small (e.g., see Figs. 3 and 7). The second subdomain, where $\sqrt{A}\left|\mathrm{\nabla}z\right|>\mathit{\theta}$, corresponds to the remaining parts of the landscape. In this subdomain, the incision rate is reduced by a uniform amount *K**θ*, relative to the rate $K\sqrt{A}\left|\mathrm{\nabla}z\right|$ that would prevail with no threshold. Note that the transition between the two subdomains at $\sqrt{A}\left|\mathrm{\nabla}z\right|=\mathit{\theta}$ entails no discontinuity in
incision rates and that if we set *θ*=0, then we obtain the governing
equation of the LEM without an incision threshold (e.g., Howard, 1994),

which is the equation that we dimensionally analyzed in Theodoratos et al. (2018).

## 2.2 Dimensions and characteristic scales

We test whether the characteristic scales defined in Theodoratos et al. (2018) are reasonable choices to analyze the LEM that includes an incision threshold *θ* (Eq. 1). We start by examining the dimensions of the variables and parameters of Eq. (1). The horizontal coordinates (*x*, *y*) have dimensions of length L; elevation *z* has dimensions of height H, and time *t* has dimensions of time T. These fundamental dimensions flow through to the derived topographic quantities *A*, $\left|\mathrm{\nabla}z\right|$, and ∇^{2}*z* and to the parameters *K*, *D*, *U*, and *θ* as seen in Table 1. LEMs typically express length and elevation in units of meters (m), and time in units of years (a). We also use these units in the simulations presented in Sect. 3, and we show incision coefficients *K* in units of per year (a^{−1}),
diffusion coefficients *D* in units of square meters per year (m^{2} a^{−1}), uplift rates *U* in units of meters per year (m a^{−1}), and incision thresholds *θ* in units of meters (e.g., Tables 3 and 4).

Given that all the terms of Eq. (1) have dimensions in H, L, and T, we can dimensionally analyze Eq. (1) using characteristic scales of length, height, and time. In Eqs. (3)–(8) below, we summarize the dimensional analysis of Theodoratos et al. (2018) as necessary background for the new analysis presented here. The dimensional analysis of Theodoratos et al. (2018) is based on a characteristic length that is defined as

a characteristic height that is defined as

and a characteristic time that is defined as

To non-dimensionalize the horizontal coordinates (*x*, *y*), elevation *z*, and time *t*, we divide them by *l*_{c}, *h*_{c}, and *t*_{c}, respectively. Specifically, we define dimensionless coordinates as (*x*^{*}, ${y}^{*}):=(x/{l}_{\mathrm{c}}$, *y*∕*l*_{c}), dimensionless elevation as ${z}^{*}:=z/{h}_{\mathrm{c}}$, and dimensionless time as ${t}^{*}:=t/{t}_{\mathrm{c}}$. We summarize the definitions of characteristic scales and of dimensionless quantities in Table 2.

To non-dimensionalize the rate of elevation change in the left-hand side of
Eq. (1) and the topographic properties in the right-hand side of Eq. (1)
(drainage area *A*, slope $\left|\mathrm{\nabla}z\right|$, and curvature ∇^{2}*z*), we need to divide them by characteristic scales that have the same dimensions as the rate of elevation change and the topographic properties, respectively. We define these characteristic scales by combining the characteristic scales of length, height, and time (*l*_{c}, *h*_{c}, and *t*_{c}) such that we obtain the appropriate dimensions. For instance, the rate of elevation change $\partial z/\partial t$ has dimensions of H T^{−1};
therefore, to non-dimensionalize it, we need to divide it by a characteristic scale with dimensions of height over time. The ratio *h*_{c}∕*t*_{c} has such dimensions. Note that ${h}_{\mathrm{c}}/{t}_{\mathrm{c}}=U$ (see Eqs. 4 and 5). Thus, we can view the uplift rate *U* as a characteristic rate of elevation change and use it to define the dimensionless rate of elevation change as
$\partial {z}^{*}/\partial {t}^{*}:=(\partial z/\partial t)/U$. Likewise, we define a characteristic area with dimensions of L^{2} as

and use it to define the dimensionless drainage area as ${A}^{*}:=A/{A}_{\mathrm{c}}$. Further, we define a characteristic gradient with dimensions of H L^{−1} as

If we divide the slope $\left|\mathrm{\nabla}z\right|$ by the characteristic gradient *G*_{c}, we obtain a dimensionless slope term. We
denote this dimensionless slope by $\left|{\mathrm{\nabla}}^{*}{z}^{*}\right|$ because it is equal to the norm of the gradient of dimensionless elevation *z*^{*} in dimensionless coordinates (*x*^{*}, *y*^{*}). See Table 2 for more details. Finally, we define a characteristic curvature with dimensions of H L^{−2} as

and we use it to define the dimensionless curvature as ${\mathrm{\nabla}}^{{*}^{\mathrm{2}}}{z}^{*}:={\mathrm{\nabla}}^{\mathrm{2}}z/{\mathit{\kappa}}_{\mathrm{c}}$. Note that the
characteristic curvature is opposite to the steady-state curvature at
ridges, drainage divides, and zones of zero incision. Specifically, if
$\partial z/\partial t=\mathrm{0}$ (steady state) and if $\sqrt{A}\left|\mathrm{\nabla}z\right|=\mathrm{0}$ in Eq. (2) or $\sqrt{A}\left|\mathrm{\nabla}z\right|\le \mathit{\theta}$ in Eq. (1), then $D{\mathrm{\nabla}}^{\mathrm{2}}z+U=\mathrm{0}$, which can be rewritten as ${\mathrm{\nabla}}^{\mathrm{2}}z=-U/D=-{\mathit{\kappa}}_{\mathrm{c}}$ (see also Roering et al., 2007; Perron et al., 2009). Note that it can be shown that −*κ*_{c} is also the minimum value of curvature in steady state.

## 2.3 Dimensionless governing equation and the incision-threshold number *N*_{θ}

If we divide all of the terms of the governing equation (Eq. 1) by the
uplift rate *U*, then we obtain an equation that includes only dimensionless
terms. We show how to non-dimensionalize the left-hand side and the
diffusion and uplift terms of Eq. (1) in Table 2. To non-dimensionalize the
incision term, we expand it to $K\sqrt{A}\left|\mathrm{\nabla}z\right|-K\mathit{\theta}$. The first part, $K\sqrt{A}\left|\mathrm{\nabla}z\right|$, corresponds to the value of the incision rate if there were no incision threshold, and the second part, *K**θ*, corresponds to the reduction in the incision rate due to the threshold. If we divide $K\sqrt{A}\left|\mathrm{\nabla}z\right|$ by *U*, then it can be shown that we obtain
the dimensionless product $\sqrt{{A}^{*}}\left|{\mathrm{\nabla}}^{*}{z}^{*}\right|$, and if we divide *K**θ* by *U*, then we obtain *K**θ*∕*U*. This ratio is dimensionless because both *K**θ* and *U* are rates of elevation change, with dimensions of H T^{−1}. We term this dimensionless ratio the incision-threshold number *N*_{θ}:

The incision-threshold number *N*_{θ} quantifies *K**θ*, the
reduction in the rate of incision due to the incision threshold, relative to
the uplift rate *U*. Additionally, *N*_{θ} is the value of the
dimensionless product $\sqrt{{A}^{*}}\left|{\mathrm{\nabla}}^{*}{z}^{*}\right|$ at the transition between the two subdomains of Eq. (1), i.e., at the interface between parts of the landscape where there is no incision and parts of the landscape where incision occurs. Specifically, at that transition, $\sqrt{A}\left|\mathrm{\nabla}z\right|=\mathit{\theta}$; if both sides of this equality are multiplied by *K* and then divided by *U*, then the equality can be shown to become $\sqrt{{A}^{*}}\left|{\mathrm{\nabla}}^{*}{z}^{*}\right|={N}_{\mathit{\theta}}$. Finally, if we rearrange Eq. (9) as ${N}_{\mathit{\theta}}=\mathit{\theta}/(U/K)$, then we see that the incision-threshold number gives the magnitude of the incision threshold *θ* relative to magnitudes of other parameters of the LEM, specifically, relative to the ratio of the uplift rate *U* to the incision coefficient *K*. Note that in the general case in which the drainage area and slope exponents *m* and *n* are not 0.5 and 1, respectively, *θ* and *K* will have different dimensions than in the case of Eq. (1), but their product, *K**θ*, will still have dimensions of H T^{−1}. Thus the ratio *K**θ*∕*U* is dimensionless for any *m* and *n*. Note that *N*_{θ} is equal to the quantity *θ*^{′} in Eq. (19) of Perron et al. (2008).

Bringing together the dimensionless terms derived above, we obtain a dimensionless form of the governing equation (Eq. 1):

Note that the dimensionless quantities that we denote by starred symbols (e.g., *z*^{*}, *A*^{*}, and $\left|{\mathrm{\nabla}}^{*}{z}^{*}\right|$) refer to variables or topographic properties. These quantities vary in space across the landscape and in time as the landscape evolves. By contrast, the incision-threshold number *N*_{θ} depends only on the model parameters *K*, *U*, and *θ*, and thus it plays the role of a parameter in Eq. (10), one that is constant in space and time. The incision-threshold number *N*_{θ} is the only parameter in Eq. (10). Thus, for a given set of boundary and initial
conditions, the value of *N*_{θ} is the only control on the solution of
Eq. (10).

The LEM without incision threshold, which we studied in Theodoratos et al. (2018), has a dimensionless form that does not include any parameters (see Eq. 16 in Theodoratos et al., 2018). Having no parameters to be adjusted, the dimensionless form has a single solution for any given combination of boundary and initial conditions. This implies that landscapes with any parameters but with the same boundary and initial conditions (when
normalized by the characteristic scales *l*_{c} and *h*_{c}) follow geometrically similar evolutionary paths; i.e., they evolve as rescaled copies of each other. We noted that this rescaling property implies that instead of running multiple simulations corresponding to multiple
combinations of parameters, we can explore the entire parameter space of the
LEM by rescaling the results of a single simulation corresponding to just one set of parameters.

In contrast, the dimensionless form of the LEM with an incision threshold,
Eq. (10), includes one parameter, the incision-threshold number *N*_{θ}. Therefore, in general, landscapes with nonzero incision thresholds will
not evolve as rescaled copies. However, Eq. (10) reveals a special case. If
landscapes have the same *N*_{θ}, i.e., if they have incision
thresholds *θ*, incision coefficients *K*, and uplift rates *U* such that they have the same ratios *K**θ*∕*U*, then they will evolve as rescaled
copies of each other, provided that their boundary and initial conditions
are the same when normalized by the characteristic scales of length and
height *l*_{c} and *h*_{c}. In Sect. 3, we numerically demonstrate both the special case of landscapes that have the same *N*_{θ} and evolve geometrically similarly and the general case of landscapes that have different *N*_{θ} and do not follow geometric similarity.

The elimination of three out of four parameter-related degrees of freedom
from the LEM (from the four parameters *K*, *D*, *U*, and *θ* in Eq. (1) to the one parameter *N*_{θ} in Eq. 10) is a substantial simplification. It is a consequence of the fact that we non-dimensionalize Eq. (1) using the
characteristic scales of length, height, and time *l*_{c}, *h*_{c}, and *t*_{c}, which depend on three model parameters (*K*, *D*, and *U*; Eqs. 3–5), and can thus eliminate an equal number of parameter-related degrees of freedom. This simplification validates the hypothesis that *l*_{c}, *h*_{c}, and *t*_{c}, as a group, remain useful in the case of Eq. (1), which includes the incision threshold *θ*. Unfortunately, with only three fundamental dimensions it is not possible to eliminate all four parameters using dimensional analysis, so one dimensionless parameter (in this case *N*_{θ}) must remain.

## 3.1 Special case: landscapes with the same *N*_{θ}

In this section, we numerically demonstrate that landscapes that follow
Eq. (1) but have different parameters will evolve geometrically similarly if
they have equal incision-threshold numbers *N*_{θ}, provided that their
boundary and initial conditions are equivalent when normalized by the
characteristic scales of length and height *l*_{c} and *h*_{c}. Given that we perform numerical simulations on discrete and finite domains, we also normalize the sizes and resolutions of these domains by *l*_{c} (see Sects. 2.2 and 3.2.2 of Theodoratos et al., 2018, for a more detailed discussion regarding the rescaling of domain size and resolution).

In this context, geometric similarity is defined in the following way. Let
the first landscape have characteristic scales *l*_{c} and *h*_{c} and the second have ${l}_{\mathrm{c}}^{\prime}$ and ${h}_{\mathrm{c}}^{\prime}$. The two landscapes are geometrically similar if any point with coordinates (*x*, *y*) and elevation *z* from the first landscape corresponds to a point from the second landscape with coordinates (*x*^{′}, *y*^{′}) and elevation *z*^{′} such that (*x*∕*l*_{c}, $y/{l}_{\mathrm{c}})=({x}^{\prime}/{l}_{\mathrm{c}}^{\prime},{y}^{\prime}/{l}_{\mathrm{c}}^{\prime})$ and $z/{h}_{\mathrm{c}}={z}^{\prime}/{h}_{\mathrm{c}}^{\prime}$. Note that both points correspond to the same point of a dimensionless landscape with coordinates (*x*^{*}, ${y}^{*})=(x/{l}_{\mathrm{c}},y/{l}_{\mathrm{c}})=({x}^{\prime}/{l}_{\mathrm{c}}^{\prime},{y}^{\prime}/{l}_{\mathrm{c}}^{\prime})$ and elevation ${z}^{*}=z/{h}_{\mathrm{c}}={z}^{\prime}/{h}_{\mathrm{c}}^{\prime}$. In other words, the two landscapes are geometrically similar if they correspond to the same
dimensionless landscape. To test whether the two landscapes are geometrically similar during their evolution, we must normalize time by their characteristic timescales *t*_{c} and ${t}_{\mathrm{c}}^{\prime}$. Specifically, we must compare a snapshot of the first landscape at some time *t* to a snapshot of the second landscape at some other time *t*^{′}, such that $t/{t}_{\mathrm{c}}={t}^{\prime}/{t}_{\mathrm{c}}^{\prime}$. Both of these snapshots correspond to
the same snapshot of a dimensionless landscape at a dimensionless time
${t}^{*}=t/{t}_{\mathrm{c}}={t}^{\prime}/{t}_{\mathrm{c}}^{\prime}$.

### 3.1.1 Setup of simulations

We perform numerical simulations using the Channel-Hillslope Integrated Landscape Development (CHILD) model (Tucker et al., 2001). Below, we briefly explain how we set up the simulations, and in Appendix A we present formulas that relate the parameters of CHILD to the parameters of the governing equation (Eq. 1). We refer readers to Theodoratos et al. (2018) for more details about setting up numerical simulations that follow geometric similarity (Sect. 3.1.1 and Appendix C) and about the theory behind such simulations (Appendix B).

For our similarity analysis, we simulate nine landscapes, each having a
different combination of the parameters *K*, *D*, and *U* and, thus, a different combination of characteristic scales of length and height *l*_{c} and *h*_{c} (Eqs. 3 and 4). Using Eq. (9), we determine the value of the incision threshold *θ* of each landscape such that the incision-threshold number of all landscapes is *N*_{θ}=1. We show the parameters, characteristic scales, and *θ* and *N*_{θ} values of the nine landscapes in Table 3. The landscapes are named with capital letters, from A to I.

Note that the incision threshold values *θ* of some of the nine landscapes are significantly higher than natural values reported in the
literature (e.g., Prosser and Dietrich, 1995; note the necessary unit
conversions). This is due to the fact that all nine landscapes have
incision-threshold numbers ${N}_{\mathit{\theta}}=\mathit{\theta}/(U/K)=\mathrm{1}$, i.e., due to the fact that each landscape's *θ* value must be equal to the value of its
*U*∕*K* ratio. We chose to use the value *N*_{θ}=1 because it leads to wide zones of zero incision (areas where, according to Eq. 1, there is no
incision, because $\sqrt{A}\left|\mathrm{\nabla}z\right|\le \mathit{\theta}$). These wide zones are readily visible when plotted.

To obtain domains and initial conditions that are equivalent when normalized
by the characteristic scales of length and height *l*_{c} and *h*_{c}, we first synthesize a random triangular irregular network (TIN) in dimensionless space, i.e., a TIN whose vertices have dimensionless horizontal coordinates (*x*^{*}, *y*^{*}) and dimensionless initial elevations *z*^{*}. This TIN's total extent is 60×90, and the average length of its triangle edges is 0.4, resulting in approximately 40 thousand TIN vertices. The initial elevations are a white noise ranging between 0 and 0.1. Second, we multiply (*x*^{*}, *y*^{*}) and *z*^{*} by each landscape's *l*_{c} and *h*_{c}, respectively. Thus we obtain each landscape's dimensional TIN with horizontal coordinates (*x*, *y*) and initial elevations *z*.

Normalizing the initial conditions is necessary for landscapes to evolve
geometrically similarly and to reach geometrically similar steady states.
Specifically, landscapes can be geometrically similar at some time step if
they were geometrically similar at the previous time step. By extension,
landscapes must start from geometrically similar initial conditions. Note
that evolving landscapes must be compared at times that are normalized by
each landscape's characteristic timescale. For example, if two landscapes
have characteristic timescales of *t*_{c} and 2*t*_{c}, then a snapshot of the first landscape with some age *t*_{0} must be compared with a snapshot of the second landscape with age 2*t*_{0}. For more details, see Appendix B in Theodoratos et al. (2018).

Note that landscapes can reach geometrically similar steady states only if
the criteria that define the steady state are normalized by appropriate
characteristic scales, as explained in Sect. 3 of Theodoratos et al. (2018).
In the present study, for instance, we assume that a simulation reaches its
steady state when the absolute rate of elevation change $|\partial z/\partial t|$ falls below a limit *ε* at all points. Given that *ε* is a rate of elevation change, we can normalize it by the uplift rate *U*, which can be viewed as a characteristic rate of elevation change, as we explain in Sect. 2.2. Thus, we set each simulation's limit to *ε*=0.001 *U*.

### 3.1.2 Results: geometric similarity

The nine simulated landscapes are all geometrically similar to each other, both during their evolution and in steady state. In Figs. 1–3, we graphically demonstrate that our simulated landscapes reach geometrically similar steady states. Specifically, we illustrate shaded relief maps in Fig. 1, elevation maps in Fig. 2, and maps of the extents of the zones of zero incision in Fig. 3. In the present study, we illustrate only steady-state results. For examples of graphical demonstrations of geometric similarity during landscape evolution, we refer readers to Figs. 3–5 of Theodoratos et al. (2018). For clarity, we present maps of only four out of the nine landscapes, specifically, of landscapes A–D in Table 3. However, all nine landscapes evolve geometrically similarly.

In Figs. 1–3, the four landscapes are arranged in a 2×2 array, such that the incision threshold *θ* increases from top to bottom and from left to right. The characteristic height *h*_{c} follows the same arrangement as *θ*, because ${h}_{\mathrm{c}}=U/K=\mathit{\theta}/{N}_{\mathit{\theta}}$, and all
landscapes have the same *N*_{θ}. The characteristic length *l*_{c} increases independently of *h*_{c} and *θ*, specifically, from bottom to top and from left to right. The coloring and labeling of Figs. 1–3 highlight both the large differences in scale and the geometric similarity of the four landscapes. Specifically, lengths and elevations on axes and color bars are shown both in units of kilometers or meters using bold fonts and in units of *l*_{c} or *h*_{c} using normal fonts. Further, color scales of elevation maps in Fig. 2 are rescaled by *h*_{c} to assist with comparing the elevations of features. Note that a quantity shown in units of the corresponding characteristic scale has the same numerical value as the dimensionless version of this quantity; e.g., elevation *z* in units of *h*_{c} has the same numerical value as dimensionless elevation *z*^{*} because both values are given by the formula *z*∕*h*_{c}. Therefore, in Figs. 1–3, the values of quantities shown in units of *l*_{c} or *h*_{c} with normal fonts are the same as the values of the corresponding dimensionless quantities.

In the shaded relief maps of Fig. 1, ridges and valleys form identical
plan-view patterns across the four landscapes, illustrating their horizontal
geometric similarity. Note that the characteristic scales of length and
height *l*_{c} and *h*_{c} vary independently, leading to different characteristic gradients *G*_{c} across the landscapes. Therefore, landscape features in different landscapes have different steepness, and, thus, they are shown with different shades of gray.

In the elevation maps of Fig. 2, the spatial pattern of colors is identical
across the four landscapes. This shows that the four landscapes are
geometrically similar both horizontally and vertically, because the color
scales are rescaled by *h*_{c}.

In Fig. 3, we map the zones of zero incision of the four landscapes. To illustrate these zones, we find the Voronoi polygons associated with points for which $\sqrt{A}\left|\mathrm{\nabla}z\right|\le \mathit{\theta}$ and we color them red. Each point of the simulated landscapes is a TIN vertex. The associated Voronoi polygon is the area that is assumed to drain to that point; see Tucker et al. (2001). We observe that the spatial patterns of the red Voronoi polygons in all four maps are geometrically similar. This implies that the zones of zero incision of the four landscapes have geometrically similar horizontal extents.

The landscapes in Figs. 1–3 do not just visually appear to be geometrically
similar. They are in fact geometrically similar. To test this quantitatively, we normalize the elevations *z* of each landscape by its characteristic height *h*_{c} and compare the resulting dimensionless elevations ${z}^{*}=z/{h}_{\mathrm{c}}$ of different landscapes. As we explain further above, the dimensionless elevations *z*^{*} of geometrically similar landscapes must be equal. Indeed, for the nine landscapes of Table 3, we find that the maximum absolute difference between steady-state *z*^{*} values of corresponding points is less than $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{9}}$.

## 3.2 General case: landscapes with different *N*_{θ}

In this subsection, we demonstrate that landscapes with different
incision-threshold numbers *N*_{θ} do not evolve geometrically similarly, even if their domains and initial conditions are rescaled by the
characteristic scales of length and height *l*_{c} and *h*_{c}. Further, we illustrate how the differences in the value of *N*_{θ} are reflected in the topography of these landscapes.

### 3.2.1 Setup of simulations

For these simulations, too, we use CHILD, as described in Appendix A. We
perform nine simulations with incision-threshold numbers *N*_{θ} that
range between 0 and 4. We use a single combination of values for the
incision coefficient *K*, diffusion coefficient *D*, and uplift rate *U*, and we vary the incision-threshold number *N*_{θ} by varying only the incision threshold *θ* (see Eq. 9). Therefore, all nine simulations have the same characteristic scales of length and height (specifically,
*l*_{c}=50 m and *h*_{c}=25 m). Thus, for all nine simulations, we use the same domains and initial conditions. Specifically, we use TINs with total extent of 150*l*_{c}×225*l*_{c} (i.e., 7.5 km × 11.25 km) and average TIN edge length of 0.4*l*_{c} (i.e., 20 m), resulting in approximately 250 thousand TIN vertices. The random initial elevations are drawn from a uniform distribution ranging between 0 and 0.1. The parameters *K*, *D*, and *U* have
values that fall within the typical range seen in the literature (e.g.,
Perron et al., 2008; Tucker, 2009). In contrast, the incision coefficients *θ* that correspond to the highest values of *N*_{θ} that we examine here have values that far exceed real-world incision threshold values typically reported (e.g., Prosser and Dietrich, 1995; note the necessary unit conversions). However, we use these high values to examine
how the LEM behaves when *N*_{θ} is high. The values of *K*, *D*, *U*, and *θ* and the corresponding *N*_{θ} of the nine landscapes are shown in Table 4.

### 3.2.2 Results: lack of geometric similarity and illustration of influence of *N*_{θ} on landscape topography

As we mentioned in the Introduction (Sect. 1), the inclusion of incision
thresholds in LEMs leads to increasing topographic slopes, decreasing
drainage densities, and more convex hillslopes (e.g., Howard, 1994; Tucker
and Bras, 1998; Perron et al., 2008). In Figs. 4–10, we illustrate these
topographic effects using steady-state results of the nine simulations
defined above (Sect. 3.2.1; Table 4). More specifically, we present shaded
relief maps (Figs. 4 and 5), maps of elevation *z* (Fig. 6), maps of the
extents of the zones of zero incision (Fig. 7), maps of curvature ∇^{2}*z* (Fig. 8), and profiles from ridge to outlet along
flow paths (Figs. 9 and 10). We show profiles along each landscape's longest
flow path to make profiles of different landscapes comparable. We mark these
flow paths with blue lines on the maps of Figs. 4–8. The maps in Figs. 4
and 6 show the full extent of the landscapes, which is 7.5 km × 11.25 km (i.e., 150*l*_{c}×225*l*_{c}), whereas the maps in Figs. 5, 7, and 8 show magnified versions of a 5 km × 4 km (i.e., 100*l*_{c}×80*l*_{c}) rectangular region from each map. To make the regions of different
landscapes comparable, we select each region such that it contains the drainage basin of the longest flow path of each landscape. We mark these
regions with blue rectangles in Figs. 4 and 6. Note that in all of these
figures, we show quantities in units of meters (or kilometers in the case of horizontal
lengths) using bold fonts and in units of the corresponding characteristic
scales using normal fonts (which yield the same numerical values as
dimensionless versions of quantities, as we explain in Sect. 3.1.2). Likewise, we show each simulation's incision threshold *θ* (in units
of meters) using bold fonts and the corresponding incision-threshold number *N*_{θ} (dimensionless) using normal fonts.

We observe that landscapes become steeper as *N*_{θ} increases.
Specifically, in the shaded relief maps (Figs. 4 and 5), hillslopes are shown
with darker shades of gray; i.e., they become steeper, and in the profile
plots (Fig. 9), the landscapes' longest flow paths become steeper. Given
that all landscapes have the same horizontal extents, the steepening of
landscapes implies that landscape relief increases. We observe the increase
in relief with *N*_{θ} both in terms of the maximum value of elevation
(see labels at the bottom of elevation maps in Fig. 6) and in terms of the
whole distribution of elevation (see profiles in Fig. 9 and the range of
colors of elevation maps in Fig. 6).

Furthermore, we observe that landscapes become less dissected and appear
smoother in plan view as *N*_{θ} increases. Specifically, in the shaded
relief maps (Figs. 4 and 5), we see that the smooth, undissected areas along
the sides of ridges and interfluves become wider, and the tips of valley
networks move away from the ridges. In the maps of curvature (Fig. 8), we
see that the valley networks become sparser, i.e., the landscapes become
less dissected. For the case of valley heads that fall on the landscapes'
longest flow paths, we see the movement away from the ridges also in the
profile plots of Fig. 9 (see blue circles).

We observe that as *N*_{θ} increases, valleys become deeper (more
concave). Specifically, in the maps of curvature (Fig. 8), the maximum value
of curvature increases with *N*_{θ}, and, thus, the positive values of
curvature become more positive. In other words, concave areas, which can be
defined as valleys (e.g., Howard, 1994), become more concave. For a given
combination of characteristic curvature *κ*_{c} and
incision-threshold number *N*_{θ}, the maximum curvature is mainly
controlled by the size of the domain. Specifically, a larger domain leads to
a larger maximum curvature value. Additionally, in the shaded relief maps
(Figs. 4 and 5), valleys in landscapes with higher *N*_{θ} appear deeper
because their contrast with neighboring hillslopes is higher. Note that the
deepening of valleys is in agreement with the steepening of hillslopes
described above.

Moreover, we observe that as *N*_{θ} increases, the zones of zero
incision (i.e., the areas with $\sqrt{A}\left|\mathrm{\nabla}z\right|\le \mathit{\theta}$; shown with red in Fig. 7) become wider and occupy bigger portions of the hillslopes. We can also observe this in the profile plots of Fig. 9. Specifically, we see that as *N*_{θ} increases, the red dots move away from the ridge and come closer to the blue circles, which implies that the longest flow paths' segments that have zero incision become longer and that they occupy bigger portions of the segments that belong to hillslopes.

Consequently, hillslopes become more convex as *N*_{θ} increases. In
steady state, the curvature in zero-incision zones is equal to −*κ*_{c} (the negative of the characteristic curvature), which is the minimum value of curvature (see Sect. 2.2). Thus, the widening of zero-incision zones implies that bigger portions of hillslopes acquire the minimum curvature, i.e., bigger portions of them become maximally convex. Note, however, that the value of the maximum convexity remains constant as
*N*_{θ} increases, because the minimum curvature remains ${\mathrm{\nabla}}^{\mathrm{2}}z=-{\mathit{\kappa}}_{\mathrm{c}}$. The maps of curvature (Fig. 8) also show that the minimum value of curvature remains constant as *N*_{θ} increases.

Finally, we observe that the widening of the zones of zero incision
eventually leads to a qualitative change in the operation of the laws of the
LEM across the landscapes. Specifically, the zones of zero incision almost
entirely occupy the hillslopes of the landscape with *N*_{θ}=4. We
deduce this by observing in Fig. 7 that the white areas (i.e., areas with
$\sqrt{A}\left|\mathrm{\nabla}z\right|>\mathit{\theta}$, where incision does operate) follow the pattern of the dendritic valley network. The almost complete occupation of hillslopes by the zones of zero incision implies that incision operates almost exclusively in valleys, which is a qualitative change. The governing equation without incision threshold (Eq. 2) is based on the fundamental assumption that all of its processes (incision, diffusion, and uplift) operate everywhere across a landscape (e.g., Howard, 1994). By including the incision threshold *θ*, the governing equation (Eq. 1) becomes piecewise, with a first subdomain with $\sqrt{A}\left|\mathrm{\nabla}z\right|\le \mathit{\theta}$ where only diffusion and uplift operate and a second subdomain with $\sqrt{A}\left|\mathrm{\nabla}z\right|>\mathit{\theta}$ where all three processes operate. This formulation does not exclude incision from hillslopes in principle. In effect, however, incision is excluded from hillslopes for high values of the
incision-threshold number *N*_{θ}, as revealed by the white dendritic
patterns in Fig. 7. Thus, for *N*_{θ}=4 the governing equation (Eq. 1)
is, in effect, reminiscent of LEMs that explicitly define distinct laws for
hillslopes and valleys (e.g., Willgoose et al., 1991; Goren et al., 2014).
Note that increasing *N*_{θ} beyond the value of 4 would not lead to
the additional qualitative change in zero-incision zones starting to occupy
valleys, because zero-incision zones have negative curvature
(${\mathrm{\nabla}}^{\mathrm{2}}z=-{\mathit{\kappa}}_{\mathrm{c}}$; see Sect. 2.2). Note that *N*_{θ}=4 is the value for which hillslopes are completely occupied by zero-incision zones in the landscapes that we examine here. However, in landscapes with different boundary and initial conditions, the qualitative change described here could occur at different values of the incision-threshold number *N*_{θ}.

With the above observations in mind, we can explain the observation that
landscapes become steeper as *N*_{θ} increases in two distinct ways,
one referring to areas outside zero-incision zones and one referring to
areas within them. First, channels become steeper to compensate for the
reduction in the strength of incision by the incision threshold.
Equation (1) shows that incision operates in areas with $\sqrt{A}\left|\mathrm{\nabla}z\right|>\mathit{\theta}$, but the rate of incision is reduced by
the quantity *K**θ* relative to $K\sqrt{A}\left|\mathrm{\nabla}z\right|$, which is the rate of incision in a landscape without incision threshold. Therefore, for a given drainage area *A*, the landscape must have steeper slope $\left|\mathrm{\nabla}z\right|$ to achieve the same incision rate and thus balance the other processes and reach equilibrium. This effect becomes stronger as *N*_{θ} increases. Second, for purely geometrical reasons, the fact that hillslopes become more convex as *N*_{θ} increases implies that they also become steeper. Typically, the more negative the Laplacian curvature ∇^{2}*z* of an
area, the faster is the increase in slope over a given flow path length.
Exceptions can be areas with negative contour curvature but positive
profile curvature, where slope decreases along flow paths, e.g., wind gaps
(see also Fig. 2c in Mitášová and Hofíerka, 1993.) Therefore, as *N*_{θ} increases and hillslopes become more convex, the slope at a
given distance from the ridge becomes steeper.

In an alternative interpretation, one could potentially view the quantity *K**θ* not as a reduction of the rate of incision but rather as a virtual source term, i.e., as a virtual increase in the uplift rate *U*. Thus the observed increase in relief would be interpreted as resulting from the virtual increase in the uplift rate because, all else remaining equal,
higher uplift rates lead to higher reliefs (e.g., Tucker and Whipple, 2002;
Theodoratos et al., 2018). However, this view is not meaningful in the zones
of zero incision, because in the first subdomain of Eq. (1) the quantity *K**θ* does not appear, and, thus, *U* is the only source term (this is also reflected in the fact that ridgelines do not become more sharply convex as they would if the uplift rate were actually increased; rather, the curvature of ridgelines remains equal to $-{\mathit{\kappa}}_{\mathrm{c}}=-U/D)$. To quantify how the uplift rate's virtual increase depends on the incision-threshold number *N*_{θ}, we can rearrange the right-hand side of the second subdomain of the governing equation (Eq. 1). We take the quantity *K**θ* from the incision term and we group it with the uplift rate *U*. Thus, we form the virtual uplift rate *K**θ*+*U*, which we rewrite as

Because Eq. (11) does not apply within the zones of zero incision, treating *K**θ* as a virtual increase in the uplift rate implies that one must also treat the landscape as having two distinct uplift rates, (*N*_{θ}+1)*U* outside the zones of zero incision and *U* within them.

Equation (11) suggests that the quantity *N*_{θ}+1 can predict how the
relief of a landscape (outside the zones of zero incision) depends on the
value of the incision-threshold number *N*_{θ}. All else being equal,
relief is proportional to the uplift rate (e.g., see definition of the
uplift erosion number *N*_{E} in Tucker and Whipple, 2002, or interpretations of our characteristic height *h*_{c} in Theodoratos et al., 2018). Therefore, Eq. (11) suggests that relief (outside zero-incision zones) is proportional to *N*_{θ}+1 (because the virtually increased uplift rate is proportional to *N*_{θ}+1), implying that elevations (outside zero-incision zones) in landscapes that differ only in their *N*_{θ} values would be equal when normalized by *N*_{θ}+1.

We can test this hypothesis using the profiles of Fig. 9, since they belong
to landscapes that have different incision-threshold numbers *N*_{θ}
but the same parameters, characteristic scales, domains, and initial
conditions (see Table 4). Specifically, we divide elevations along each
profile of Fig. 9 by *N*_{θ}+1, and we plot the resulting normalized
profiles in Fig. 10. The hypothesis will not be rejected if the normalized
profiles have the same normalized elevations outside the zones of zero
incision. Indeed, we observe that the normalized elevations are nearly equal, especially in those reaches of each profile that are not near its zone of zero incision. This suggests that away from the zero-incision zones, landscape relief nearly scales with *N*_{θ}+1.

In Fig. 10, we observe that the elevations of the normalized profiles deviate systematically from one another. Specifically, we observe that, whereas the reliefs of the original (un-normalized) profiles grow as *N*_{θ} increases, the reliefs of the normalized profiles decrease as *N*_{θ} increases. In Table 5 we show
an example of this reversal using the original and normalized elevations of the profiles at a distance of 0.5 km from the ridge, which falls outside the
zones of zero incision of all profiles; see arrows in Fig. 10. This
reversal implies that normalizing elevations by *N*_{θ}+1 is an
overshoot, as it lowers the profiles by a larger factor than what would make
them equal to each other. In other words, as *N*_{θ} increases, the
elevation of the original profiles is increased less than proportionally to
*N*_{θ}+1, i.e., less than what is predicted by viewing the quantity *K**θ* as a virtual increase in the uplift rate. This suggests that the incision threshold could be resulting in additional effects, which oppose the virtual increase in the uplift rate, and that these effects depend on the value of *N*_{θ}. Future work can study such effects. For example, it is known that incision's competition with diffusion for the propagation of elevation perturbations can be influenced by the incision threshold (e.g., see relationship between the Péclet number and the incision threshold in Perron et al., 2008). Thus it may be productive to examine whether changes in the competition between incision and diffusion alter how the incision threshold affects the rate of incision.

## 4.1 On the definition of zones of zero incision

Unlike the LEM studied here, other LEMs, such as those of Tucker (2004) or Deal et al. (2018), do not define zones of zero incision, i.e., areas where incision never operates, because those LEMs define incision terms based on conceptually different temporal averaging of rainfall events, in comparison to the LEM examined here.

Specifically, those other LEMs derive long-term incision rates by integrating stochastic rainfall over time, assuming that incision occurs when the shear stress (or, equivalently, the stream power) exceeds a threshold value. Given that the value of shear stress depends on discharge and slope, points with different drainage areas or slopes will experience different shear stress values during any given event. Therefore, any given combination of drainage area and slope corresponds to a critical rainfall intensity that is sufficient to generate a shear stress that equals the threshold shear stress. Long-term incision rates can be derived by integrating over the rainfall events that exceed this critical rainfall intensity. This approach implies that, in theory, any point with nonzero drainage area and slope can experience incision with a nonzero probability (provided that rainfall can theoretically become sufficiently intense). Therefore, in LEMs that follow this approach, zero-incision zones are not defined. Note, however, that in those LEMs one can define zones of probability of exceedance of the critical rainfall intensity, i.e., of probability of incision.

In contrast, the LEM studied here assumes constant, uniform rainfall, which
leads to constant stream power for any given combination of drainage area
and slope (i.e., for any given value of $\sqrt{A}\left|\mathrm{\nabla}z\right|$). Therefore, instead of explicitly including a stream-power threshold, the LEM's governing equation (Eq. 1) uses a threshold of the quantity $\sqrt{A}\left|\mathrm{\nabla}z\right|$ itself, namely, the incision threshold *θ* (see the relationship between *θ* and the threshold of stream power in Appendix A). This formulation of Eq. (1) has the advantage of being much simpler than those of LEMs that use stochastic rainfall and shear-stress (or stream-power) thresholds. However, Eq. (1) has the disadvantage of being unable to explore the nonlinear relationship between average precipitation and long-term incision rates that we describe in the Introduction (Sect. 1).

## 4.2 On the choice of characteristic scales

In this study, we have examined whether the characteristic scales of length,
height, and time (*l*_{c}, *h*_{c}, and *t*_{c}; Eqs. 3–5), which we introduced in Theodoratos et al. (2018), remain useful after the inclusion of an incision threshold in the LEM, and we find that they do. Furthermore, while non-dimensionalizing Eq. (1) using this group of characteristic scales, we obtain the dimensionless incision-threshold number *N*_{θ}, which has useful properties. These results, however, do not imply that *l*_{c}, *h*_{c}, and *t*_{c} are the only possible choices of characteristic scales or even that they are the best choices for all problems. For any given model, different characteristic scales may be more appropriate for different applications.

Dimensional analysis can ensure that a set of characteristic scales is
dimensionally consistent and can provide the number of degrees of freedom
that can be eliminated from a model (e.g., Buckingham, 1914), but it cannot
show a priori which characteristic scales should be used. For example, in
the case of Eq. (1), if we assume that length L and height H are distinct
dimensions, then together with time T they form a group of three dimensions, and dimensional analysis will show that any manipulation of Eq. (1) can eliminate at most 3 degrees of freedom. Because the characteristic scales *l*_{c}, *h*_{c}, and *t*_{c} are defined by the parameters *K*, *D*, and *U*, eliminating 3 degrees of freedom eliminates these three parameters. If, instead, one defined characteristic scales that depended, for example, on the measurements of the domain (e.g., Perron et al., 2008), the corresponding degrees of freedom that could be eliminated would be related to these domain scales. Such an approach might be more appropriate for characterizing extensive properties of a landscape as a whole (e.g., Perron et al., 2012), whereas the approach that we use here may be more appropriate for characterizing processes and intensive properties that vary across a landscape (e.g., Theodoratos et al., 2018). It may be difficult to predict a priori which choices of characteristic scales will be better for a given problem, and the only way to find out may be to try several different alternatives. In general, dimensional analysis can be used to simplify governing equations, and it can point to useful numerical, field, or lab experiments, but it cannot fully substitute for the information contained in empirical results (e.g., Huntley, 1967).

A future study could examine the utility of characteristic scales that are defined to depend on model parameters (Eqs. 3–5) in the case of parameters that vary in space or time. In such a case, the characteristic scales would also vary. We expect that if the parameter variation is gradual and follows a systematic pattern (e.g., the differential uplift across a fold described by Kirby and Whipple, 2001), then the resulting variable characteristic scales could be useful. For example, designing a lab-scale sandbox landscape that models differential uplift might benefit from these nonuniform characteristic scales. However, if the parameters were randomly heterogeneous, or they varied greatly over distances much smaller than typical landscape units, then the resulting “characteristic” scales might not be characteristic of any landscape properties, and thus they might lose their explanatory power.

In this study, we perform a dimensional analysis of an LEM that includes
terms describing stream-power incision, linear diffusion, and uplift
(Eq. 1). The LEM assumes that incision is limited by a threshold,
specifically, that there is no incision at points with drainage area *A* and slope $\left|\mathrm{\nabla}z\right|$ such that the quantity $\sqrt{A}\left|\mathrm{\nabla}z\right|$ is below a threshold value *θ* and that this threshold also reduces incision at points with $\sqrt{A}\left|\mathrm{\nabla}z\right|>\mathit{\theta}$.

Our dimensional analysis is based on characteristic scales of length, height, and time (*l*_{c}, *h*_{c}, and *t*_{c}) that depend only on parameters of the LEM (specifically, on the incision coefficient *K*, diffusion coefficient *D*, and uplift rate *U*; Eqs. 3–5). We introduced these scales as a group in Theodoratos et al. (2018), where we analyzed a related LEM that did not include the incision threshold (reproduced here as Eq. 2). The distinction between *l*_{c} and *h*_{c} is based on the assumption that horizontal lengths and vertical heights are dimensionally distinct.

In Sect. 2.3, using the characteristic scales *l*_{c}, *h*_{c}, and *t*_{c}, we derive Eq. (10), a dimensionless form of the governing equation of the LEM that includes only one parameter, the incision-threshold number ${N}_{\mathit{\theta}}=K\mathit{\theta}/U=\mathit{\theta}/(U/K)$ (Eq. 9). This dimensionless parameter quantifies the value of *K**θ*, which is the reduction in the rate of incision due to the incision threshold, relative to the uplift rate *U* or, equivalently, the relative magnitude of the incision threshold *θ* versus the ratio *U*∕*K*. The original, dimensional LEM (Eq. 1) includes four parameters (*K*, *D*, *U*, and *θ*). Because the three characteristic scales (*l*_{c}, *h*_{c}, and *t*_{c}) depend on three model parameters (*K*, *D*, and *U*), in deriving the dimensionless Eq. (10) we can eliminate three out of four parameter-related degrees of freedom, which is a notable simplification. This suggests that this group of characteristic scales is useful in the case of the LEM that includes an incision threshold.

The fact that the incision-threshold number *N*_{θ} is the only
parameter in the dimensionless governing equation (Eq. 10) implies that it
is the only control on this equation, for any given set of boundary and
initial conditions. As a consequence, the evolution of all landscapes with a
given *N*_{θ} value will be geometrically and temporally similar,
provided that their domains, boundary conditions, and initial conditions are
rescaled by *l*_{c} and *h*_{c} (see Theodoratos et al., 2018, for more detailed theoretical exposition of these arguments). In Sect. 3.1, we present numerical simulations of landscapes with different parameters but
equal incision-threshold numbers *N*_{θ}. In Figs. 1–3, we demonstrate
that these landscapes indeed evolve geometrically similarly. In contrast,
landscapes with different *N*_{θ} values evolve without geometric
similarity, as we show with a second set of numerical simulations in Sect. 3.2.

In Sect. 3.2.2, we explore how these different *N*_{θ} values influence
steady-state topographic properties of the resulting landscapes. We illustrate the topographic influence of *N*_{θ} in Figs. 4–10. As *N*_{θ} increases, we find that valleys become sparser and deeper (Figs. 4, 5, and 8), that hillslopes become more convex (Figs. 7 and 8), and that relief increases (Figs. 6 and 9) and, thus, hillslopes and channels become steeper (Figs. 4, 5, and 9).

Finally, we derive a quantitative prediction of the increase in relief with
the incision-threshold numbers *N*_{θ}. Specifically, we show that in
areas with $\sqrt{A}\left|\mathrm{\nabla}z\right|>\mathit{\theta}$, where incision operates, relief tends to scale with the quantity *N*_{θ}+1, and thus elevations tend to become equal if they are normalized by *N*_{θ}+1 (Fig. 10). Our simulation results show deviations from this prediction, but we observe that these deviations are systematic (Sect. 3.2.2, Table 5), and we posit that the incision threshold causes additional effects, which can be the focus of future work.

To implement the governing equation of the LEM (Eq. 1) with CHILD, we use
CHILD's detachment-limited module, and we set the parameter DETACHMENT_LAW equal to 0. Furthermore, we use constant, uniform, and continuous precipitation, define infiltration to be 0, and set the hydraulic geometry scaling exponents *ω*_{b} and *ω*_{s} to be equal to 0.5 and the detachment capacity exponents *m*_{b}, *n*_{b}, and *P*_{b} to be equal to 1 (see Tucker et al., 2001, and Tucker, 2010, for definitions of CHILD's assumptions, modules, and
parameters).

For this choice of exponents, CHILD uses the following equations to calculate the rate of elevation change due to incision (in CHILD notation):

where *D*_{c} is the maximum detachment capacity in meters per year (m a^{−1}), *τ*_{0} is stream power per unit bed area in watts per square meter (W m^{−2}), *τ*_{c} is the threshold value of *τ*_{0}, below which there is no incision, also in watts per square meter (W m^{−2}), *k*_{b} is a detachment rate coefficient (units: m a^{−1} (W m${}^{-\mathrm{2}}{)}^{-\mathrm{1}}$) (i.e., *k*_{b} is the rate of elevation change per each unit of stream power per unit bed area), *k*_{t} is the specific weight of water in newtons per cubic meter (N m^{−3}), *P* is the precipitation intensity in meters per year (m a^{−1}), *k*_{w} is bankfull width per unit scaled discharge (units: s^{0.5} m^{−0.5}), and *S* is slope (Tucker et al., 2001; Tucker, 2010).

Equating $K\left(\sqrt{A}\right|\mathrm{\nabla}z|-\mathit{\theta})$, i.e., the incision term of Eq. (1), to *D*_{c} of Eq. (A1a) we can relate the incision coefficient *K* and the incision threshold *θ* of Eq. (1) with CHILD's parameters according to

Equations (A2) and (A3) include the unit conversion factor $\sqrt{\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{a}}/\sqrt{\mathrm{31557600}\phantom{\rule{0.125em}{0ex}}\mathrm{s}}$ because the input files of CHILD include variables with units of both years and seconds.

In Eqs. (A2) and (A3), we assume constant values of *k*_{t}=9810 N m^{−3}, *P*≈1.31 m a^{−1}, and
*k*_{w}=10 s^{0.5} m^{−0.5}, and we obtain the desired values of *K* and *θ* by entering the appropriate values of *k*_{b} and *τ*_{c} into CHILD's input files.

All data used in this study were synthesized using the CHILD model (Tucker et al., 2001). The input files needed to reproduce this data are available from the corresponding author upon request.

NT derived analytical results, and NT and JWK interpreted them. NT designed, performed, and analyzed numerical simulations, and NT and JWK interpreted them. NT drafted the paper, and NT and JWK edited it.

The authors declare that they have no conflict of interest.

This study was made possible by financial support from ETH Zurich.

This paper was edited by Jean Braun and reviewed by Wolfgang Schwanghart, Eric Deal, and one anonymous referee.

Buckingham, E.: On physically similar systems; illustrations of the use of dimensional equations, Phys. Rev., 4, 345–376, 1914.

Deal, E., Braun, J., and Botter, G.: Understanding the role of rainfall and hydrology in determining fluvial erosion efficiency, J. Geophys. Res.-Earth, 123, 744–778, https://doi.org/10.1002/2017JF004393, 2018.

DiBiase, R. A. and Whipple, K. X.: The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate, J. Geophys. Res., 116, F04036, https://doi.org/10.1029/2011JF002095, 2011.

Dietrich, W. E., Bellugi, D. G., Sklar, L. S., Stock, J. D., Heimsath, A. M., and Roering, J. J.: Geomorphic Transport Laws for Predicting Landscape form and Dynamics, in: Prediction in Geomorphology, Geophysical Monograph Series 135, edited by: Wilcock, P. R. and Iverson, R. M., American Geophysical Union, Washington, D.C., USA, 103–132, https://doi.org/10.1029/135GM09, 2003.

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.

Horton, R. E.: Erosional development of streams and their drainage basins; hydrophysical approach to quantitative morphology, Geol. Soc. Am. Bull., 56, 275–370, 1945.

Howard, A. D.: A detachment-limited model of drainage basin evolution, Water Resour. Res., 30, 2261–2285, 1994.

Huntley, H. E.: Dimensional analysis, Dover Publications, New York, USA, 1967.

Kirby, E. and Whipple, K.: Quantifying differential rock-uplift rates via stream profile analysis, Geology, 29, 415–418, 2001.

Lague, D., Hovius, N., and Davy, P.: Discharge, discharge variability, and the bedrock channel profile, J. Geophys. Res., 110, F04006, https://doi.org/10.1029/2004JF000259, 2005.

Mitášová, H. and Hofierka, J.: Interpolation by regularized spline with tension: II. Application to terrain modeling and surface geometry analysis, Math. Geol., 25, 657–669, https://doi.org/10.1007/BF00893172, 1993.

Montgomery, D. R. and Dietrich, W. E.: Channel Initiation and the Problem of Landscape Scale, Science, 255, 826–830, 1992.

Perron, J. T.: Climate and the pace of erosional landscape evolution, Ann. Rev. Earth Plan. Sc., 45, 561–591, https://doi.org/10.1146/annurev-earth-060614-105405, 2017.

Perron, J. T., Dietrich, W. E., and Kirchner, J. W.: Controls on the spacing of first-order valleys, J. Geophys. Res., 113, F04016, https://doi.org/10.1029/2007JF000977, 2008.

Perron, J. T., Kirchner, J. W., and Dietrich, W. E.: Formation of evenly spaced ridges and valleys, Nature, 460, 502–505, https://doi.org/10.1038/nature08174, 2009.

Perron, J. T., Richardson, P. W., Ferrier, K. L., and Lapôtre, M.: The root of branching river networks, Nature, 492, 100–103, https://doi.org/10.1038/nature11672, 2012.

Prosser, I. P. and Dietrich, W. E.: Field experiments on erosion by overland flow and their implication for a digital terrain model of channel initiation, Water Resour. Res., 31, 2867–2876, 1995.

Roering, J. J., Perron, J. T., and Kirchner, J. W.: Functional relationships between denudation and hillslope form and relief, Earth Planet. Sc. Lett., 264, 245–258, https://doi.org/10.1016/j.epsl.2007.09.035, 2007.

Rossi, M. W., Whipple, K. X., and Vivoni, E. R.: Precipitation and evapotranspiration controls on daily runoff variability in the contiguous United States and Puerto Rico, J. Geophys. Res.-Earth, 121, 128–145, https://doi.org/10.1002/2015JF003446, 2016.

Snyder, N. P., Whipple, K. X., Tucker, G. E., and Merritts, D. J.: Importance of a stochastic distribution of floods and erosion thresholds in the bedrock river incision problem, J. Geophys. Res., 108, 2117, https://doi.org/10.1029/2001JB001655, 2003.

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.

Tucker, G. E.: Drainage basin sensitivity to tectonic and climatic forcing: Implications of a stochastic model for the role of entrainment and erosion thresholds, Earth Surf. Proc. Land., 29, 185–205, https://doi.org/10.1002/esp.1020, 2004.

Tucker, G. E.: Natural experiments in landscape evolution, Earth Surf. Proc. Land., 34, 1450–1460, https://doi.org/10.1002/esp.1833, 2009.

Tucker, G. E.: CHILD Users Guide for version R9.4.1, Coop. Inst. Res. Environ. Sci. (CIRES) and Dep. Geol. Sci., Univ. Colo., Boulder, CO, USA, available at: http://csdms.colorado.edu/mediawiki/images/Child_users_guide.pdf (last access: 30 March 2018), 2010.

Tucker, G. E. and Bras, R. L.: Hillslope processes, drainage density, and landscape morphology, Water Resour. Res., 34, 2751–2764, 1998.

Tucker, G. E. and Whipple, K. X.: Topographic outcomes predicted by stream erosion models: Sensitivity analysis and intermodal comparison, J. Geophys. Res., 107, 2179, https://doi.org/10.1029/2001JB000162, 2002.

Tucker, G. E., Lancaster, S., Gasparini, N., and Bras, R.: The channel-hillslope integrated landscape development model (CHILD), in: Landscape erosion and evolution modeling, III, edited by: Harmon, R. S. and Doe, W. W., Kluwer Academic/Plenum Publishers, New York, USA, 349–388, 2001.

Willgoose, G., Bras, R. L., and Rodriguez-Iturbe, I.: A coupled channel network growth and hillslope evolution model: 1. Theory, Water Resour. Res., 27, 1671–1684, https://doi.org/10.1029/91WR00935, 1991.

- Abstract
- Introduction
- Dimensional analysis of LEM that includes incision threshold
- Numerical simulations
- Discussion
- Summary and conclusions
- Appendix A: Implementation of governing equation with CHILD
- Data availability
- Author contributions
- Competing interests
- Acknowledgements
- Review statement
- References

- Abstract
- Introduction
- Dimensional analysis of LEM that includes incision threshold
- Numerical simulations
- Discussion
- Summary and conclusions
- Appendix A: Implementation of governing equation with CHILD
- Data availability
- Author contributions
- Competing interests
- Acknowledgements
- Review statement
- References