Interactive comment on “ Landscape evolution models using the stream power incision model show unrealistic behavior when m / n equals 0 . 5 ” by

We agree that it is readily seen that there is a slope singularity at the ridge at steady state when drainage area goes to zero. However, without integrating the conservation equation at steady state, the existence of an elevation singularity at the ridges for m/n ≥ 0.5, and its absence for m/n < 0.5 cannot be easily deduced. This is especially true in the 2D model, which has no analytical solution. Within the literature, there has been little discussion specifically oriented to the singular behavior of SPIM in regard to slope at the ridge. To our knowledge, the presence or absence of the elevation singularity (according to the value of m/n) in 2D models has never been shown in the literature.


Introduction
The stream power incision model (SPIM) (e.g., Howard, 1994;Howard et al., 1994) is a commonly used physically based model for bedrock incision.The incision rate, E, can be written as where K is the erodibility coefficient, A is the upslope drainage area, S is the downstream slope, and m and n are exponents.This simple model is thoroughly reviewed in Whipple and Tucker (1999) and Lague (2014), where they hypothesize that m/n is between 0.35 and 0.60.This range is consistent with results inferred from field work and map studies (Flint, 1974;Howard and Kerby, 1983;Tarboton et al., 1989Tarboton et al., , 1991;;Willgoose et al., 1990;Willgoose, 1994;Moglen and Bras, 1995;Snyder et al., 2000).Furthermore, many researchers specifically suggest, or offer as a default, the ratio m/n ∼ 0.5 (Snyder et al., 2000;Banavar et al., 2001;Hobley et al., 2017).The choice of this ratio is paramount in numerical landscape evolution models (LEMs) that utilize SPIM, such as the channel-hillslope integrated landscape development model, CHILD (Tucker et al., 2001).The ratio m/n is also used to describe the relationship between slope and drainage area in describing stream long profiles (Flint, 1974).All models using SPIM, including studies on drainage reorganization and stability (Willett et al., 2014), tectonic histories of landscapes (Goren et al., 2014b;Fox et al., 2014), and persistent drainage migration (Pelletier, 2004), involve the specification of this ratio.In addition, the specific values of m and n are important (Tucker and Whipple, 2002).
Here, however, we focus on the ratio itself, and we show a somewhat unexpected result: when m/n = 0.5, SPIM-based LEMs exhibit elevation solutions that are invariant to shapepreserving stretching of horizontal domain.That is, except for the finest scales on which hillslope diffusion becomes important, the model predicts the same solution for a landscape with a total basin area of 10 km 2 and one with a total basin area of 1000 km 2 under the constraint of identical horizontal basin shape (e.g., square).The extremity of this result underscores a heretofore unrecognized unrealistic aspect of SPIM.
In this paper, we perform a scaling analysis of SPIM.First, we use a 1-D model to analytically derive steady-state river profiles, to illustrate the problem of scale invariance, and to delineate conditions for which elevation singularities occur at the ridge.Then, using a 2-D numerical model, we demonstrate the effects of horizontal scale on the steady-state relief of landscapes and infer the conditions for which elevation singularities occur at ridges.

Motivation
SPIM is a simple model that has been used to gain considerable insight into landscape evolution.Previous studies using SPIM have shown how landscapes respond to tectonic and climate forcing (e.g., Howard, 1994;Howard et al., 1994).Yet like most simple models, SPIM is in some sense an oversimplification.Here, we demonstrate this by showing that it satisfies a curiously unrealistic scale invariance relation.By demonstrating this limitation, we hope to motivate the formulation of models that overcomes it.
The fundamental limitation on SPIM becomes apparent when the ratio m/n = 0.5.Under this condition, SPIM alone will predict the same steady-state relief for a 10 km 2 domain as a 1000 km 2 domain of the same horizontal shape, as illustrated below.LEMs utilizing SPIM often sidestep this problem with the use of a "hillslope diffusion" coefficient (e.g., Passalacqua et al., 2006), a useful but rather poorly constrained parameter that lumps together a wide range of processes (Fernandes and Dietrich, 1997).Alternatively, the problem can be sidestepped with an externally specified "hillslope critical length" (Goren et al., 2014a) that essentially specifies the location of channel heads.For example, the model simulations of Willett et al. (2014) employ the specific value of 500 m for hillslope critical length in their characterization of tendencies for drainage divide migration.The prediction of the hillslope diffusion coefficient and the location of channels are outstanding problems in the field of geomorphology (Montgomery and Dietrich, 1988).The intrinsic nature of the SPIM model, however, is such that scale invariance persists for the case m/n = 0.5 on scales larger than a characteristic hillslope length scale, whether it be externally specified or computed from a diffusion coefficient.
The existence of scale invariance exemplifies an unrealistic aspect of SPIM, which we believe to be associated with its omission of natural processes, such as abrasion due to sediment transport.Gilbert (1877) theorized two roles that sediment moving as bed load could play in bedrock incision: the first as an abrasive agent that incises the bed via collisions and the second as a protector that inhibits collisions of bed load on the bed.These observations have been implemented quantitatively by many modelers (e.g., Sklar and Dietrich, 2001, 2004, 2006;Lamb et al., 2008;Zhang et al., 2015), some of whom have implemented them in LEMs (e.g., Gasparini et al., 2006Gasparini et al., , 2007)).Egholm et al. (2013) have directly compared landscape models using SPIM on the one hand and models using a saltation-abrasion model on the other hand.
Here, we shed light on an unrealistic behavior of SPIM with the goal of motivating the landscape evolution community to develop more advanced treatments that better capture the underlying physics.A further goal is to emphasize the importance of scaling and nondimensionalization in characterizing LEMs.

One-dimensional model: scale invariance and singularities
An LEM can be implemented using the following equation of mass conservation for rock/regolith subject to uplift and denudation: where η is the local landscape elevation, t is time, υ is the rock uplift rate, and D is the hillslope diffusion coefficient.The term D∇ 2 η accounts for hillslope diffusion (Somfai and Sander, 1997;Banavar et al., 2001).The effect of diffusion is commonly neglected at coarse-grained resolution (Somfai and Sander, 1997;Banavar et al., 2001;Passalacqua et al., 2006), at which any resolved channels can be taken to be fluvially dominated bedrock channels (Montgomery and Foufoula-Georgiou,1993).In our analysis, we use Eq. ( 1) to specify the incision term in Eq. ( 2).It should be noted that SPIM refers to the incision in the direction normal to the bed, implying that there are both horizontal and vertical components of incision.In much of the literature using SPIM, however, the horizontal component is neglected in accordance with the original formulation of Howard and Kerby (1983), and incision is assumed to be purely vertical and downward.
Here, we preserve this simplification in order to better understand the overall behavior of SPIM.Last, in correspondence with most 2-D implementations of SPIM within LEM, we neither resolve channels nor compute their hydraulic geometry in our 2-D implementation.The focus of this paper is the most simplified form (e.g., Eq. 1) of SPIM.This way we can analyze the most fundamental behavior of SPIM itself.Equation (2) characterizes landscape evolution in 2-D; i.e., elevation η = η (x, y), where x and y are horizontal coordinates.It is useful for some purposes, however, to simplify Eq. (2) into a 1-D form.Neglecting hillslope diffusion, the 1-D conservation equation is where l is the horizontal stream distance from the ridge, at which l = 0.It should be noted that the negative sign appears in front of the term ∂η/∂l because ∂η/∂l is negative in the downstream direction, so that streambed slope S = −∂η/∂l.In SPIM, slope S is assumed to be positive.In order to solve Eq. ( 3), a relationship between A and l must be established.
Here, we assume a generalized form of Hack's Law (Hack, 1957): where C and h are positive values.Hack's Law assumes that upslope area increases with l h .From empirical data, Hack found the exponent h to be ∼ 1.67 (Hack, 1957).
Previous researchers have presented 1-D analytical solutions for elevation profiles (Chase, 1992;Beaumont et al., 1992, Anderson, 1994;Kooi andBeaumont, 1994, 1996;Tucker and Slingerland, 1994;Densmore et al., 1998;Willett, 1999Willett, , 2010;;Whipple and Tucker, 1999).In their solutions, the effect of the horizontal scale, which in the 1-D model we define as the total length of the stream profile, L 1-D , was neither shown nor discussed.Previous studies that use Eq. ( 4) (Whipple and Tucker, 1999;Willett, 2010) involve nondimensionalization of both the horizontal and vertical coordinates by the total horizontal length of the profile, L 1-D .As we show below, this step obscures the effect of the horizontal scale on the relief of the profile.In our study, we nondimensionalize the vertical coordinate, η, by a combination of υ and the acceleration of gravity, g.Our nondimensionalization of the coordinates is shown below.
Substituting Eqs. ( 4) and ( 5) into Eq.( 3) results in the following dimensionless conservation equation: where the dimensionless number P 1-D , termed the 1-D Pillsbury number herein for convenience, is given by the relation At steady state, Eq. ( 6) becomes From this equation, we see that as we approach the ridge, i.e., l → 0, the slope term −∂ η/∂ l always approaches infinity for positive values of h, m, and n.
The value of the 1-D Pillsbury number P 1-D increases with stream profile length L 1-D when hm/n < 1, is invariant to changes in L 1-D when hm/n = 1, and decreases with L 1-D when hm/n > 1.This can be further illustrated by integrating Eq. ( 8).To solve this first-order differential equation, we need to specify a single boundary condition, shown below.
This boundary condition sets the location and elevation of the outlet, where flow is allowed to exit the system.Integrating Eq. ( 8) yields The steady-state profiles defined by Eq. ( 10) are shown in Fig. 1.Inspecting Eq. ( 10), we see that elevation is infinite at the ridge (l = 0) when hm/n ≥ 1, and elevation is finite when hm/n < 1.In addition, when hm/n = 1, P 1-D is no longer dependent on the horizontal scale, L 1-D , and η is independent of the scale of the basin.Using the empirical value from Hack's original work (1957), i.e., h = 1.67, the ratio m/n must take the value 0.6 for scale invariance.This ratio is within the range reported in the literature (Whipple and Tucker, 1999).

Two-dimensional model: scale invariance
In 2-D, the conservation equation using SPIM and neglecting hillslope diffusion can be written as To understand the behavior of Eq. ( 11) in response to scale, we need to use a dimensionless formulation in a fashion similar to the previous 1-D analysis.Here, L 2-D denotes the horizontal length of the entire domain, which is taken to be square for convenience.For the 2-D analysis, our nondimen- The form of Eq. ( 11) in which x, y, and A have been made dimensionless using the definitions shown in Eq. ( 12) is where the dimensionless number P 2-D , termed the 2-D Pillsbury number, is given as At steady state, Eq. ( 13) becomes The form of the parameter P 2-D specified by Eq. ( 14) is similar to the 1-D form, Eq. ( 7), but is different due to the different dimensionality.The parameter, P 2-D , scales with the relief of the landscape; as it increases, the slope term on the RHS (right-hand side) of Eq. ( 15) also increases.
The value of P 2-D increases with L 2-D for m/n < 0.5, remains constant with L 2-D for m/n = 0.5, and decreases with L 2-D for m/n > 0.5.For the ratio m/n = 0.5, the exponent to which L 2-D is raised in Eq. ( 14) becomes 0 and the relief of the landscape becomes invariant to horizontal scale.When m/n = 0.5, the same steady-state solution to Eq. ( 15) prevails regardless of the value of L 2-D .We note here that this scale invariance, which is the key result of this paper, is intrinsic to the model itself and is not a function of the discretization scheme used in implementing numerical solutions.
Our 2-D model was solved using the following boundary conditions: The bottom (outlet) side of the domain presented in Fig. 2 is fixed at the base level η = 0 m, corresponding to an open boundary where flow can exit the system while satisfying Eq. ( 16).The top side of the domain is designated as an impermeable boundary to flow, i.e., the drainage divide satisfies Eq. ( 17).Periodic boundary conditions satisfying Eq. ( 18) are applied at the left and right boundaries.Flow, slope, and drainage area are determined using the D8 flow algorithm, where flow follows the route of steepest descent (O'Callaghan and Mark, 1984).The initial condition is a gently sloped plane oriented towards the outlet with small random elevation perturbations.
For the results of Fig. 2, we use regular grids that contain 100 2 cells.The number of cells is constant, regardless of the value of L 2-D .This is in contrast to holding cell size constant and instead increasing the number of cells with L 2-D .We argue that the former shows the fundamental behavior of SPIM, while the latter obscures this behavior due to the existence of slope and elevation singularities near the ridges in the landscape.The next sections show this singular behavior in the 2-D numerical model.
Figure 2a shows steady-state solutions for m/n = 0.5 and two values of L 2-D using the same initial condition.At each corresponding grid cell between the two solutions, the slope, S, decreases as L 2-D increases.However, the relief structures of each landscape are identical.By relief structure, we are describing the elevation value at each corresponding grid cell in the two steady-state solutions.This is confirmed by nondimensionalizing the horizontal scale of landscape without adjusting the vertical scale (Fig. 2b).Using the same numerical methods and the parameters from Fig. 2a, the results of a similar analysis using different ratios m/n = 0.4, 0.5, and 0.6 are shown in Fig. 2c.
In Fig. 2c, the case of scale invariance can be seen when m/n = 0.5.For m/n = 0.4, the relief of the entire landscape increases with increasing L 2-D , and for m/n = 0.6, the relief decreases with increasing L 2-D .When m/n = 0.5, the landscapes do not exhibit scale invariance.However, the overall planform drainage network structure shows resemblance across scales.That is, the location of the major streams and rivers in the numerical grid are similarly organized.It should be noted that the landscapes are not identical.When the landscapes are shown in dimensional space, as shown in Fig. 2a, the landscapes appear to be quite different.In the case of Fig. 2b, however, the smaller landscape can be stretched horizontally to be precisely identical to the large one.The drainage network structure described above persists in each simulation due to the imprinting of the initial condition, which always consists of the same randomized perturbations.

Two-dimensional model: quasi-theoretical analysis of singular behavior
Like the 1-D model of Eq. ( 8), the 2-D model, Eq. ( 15), has slope, S, approaching infinity as area, A, approaches 0 at steady state.In contrast to the 1-D model, however, general steady-state solutions for elevation in the 2-D model, Eq. ( 15), cannot be determined analytically.However, the ratio m/n for which elevation singularities occur can be determined by analyzing the behavior of the 2-D numerical model in close proximity to a ridge.Here, we first develop a quasitheoretical treatment to study near-ridge behavior, and we then use it to infer singular behavior in the numerical model.Converting the coordinate system from Cartesian to a system that follows the stream-wise direction, we rewrite Eq. ( 11) as Earth Surf.Dynam., 5, 807-820, 2017 www.earth-surf-dynam.net/5/807/2017/where s is the distance along the path of steepest descent away from the ridge.From dimensional considerations, A [L 2 ] must scale with s 2 [L 2 ] near the ridge (s = 0), and therefore, where β is the scaling factor.For this analysis, our nondimensionalization is where L R is the horizontal ridge scale.Near the ridge, Eq. ( 19) can be nondimensionalized into where P R is another dimensionless Pillsbury number, here denoted as At steady state (∂η/∂t = 0), Eq. ( 22) becomes From Eq. ( 24), we see that at the ridge (ŝ = 0), there is a singularity in slope, i.e., the slope (−∂ η/∂ ŝ) goes to infinity.Integration of Eq. ( 24) using the downstream boundary condition, η ŝ=1 = 0, allows for the delineation of the conditions for elevation singularities in the 2-D model.The profile is given as Instead of the elevation singularity occurring when hm/n ≥ 1 as seen in the 1-D model, Eq. ( 10), this analysis for the 2-D model shows an elevation singularity at the ridge when m/n ≥ 0.5.= 40 2 , 80 2 , and 160 2 .To make the relief of the landscapes comparable, the 2-D Pillsbury number, P 2-D , is set to 3.10 × 10 23 for solutions of all m/n ratios with L 2-D = 10 km.To achieve this for υ = 1 mm yr −1 , K = 6.31 × 10 −12 m 0.2 s −1 , 1.00 × 10 −12 s −1 , and 1.58 × 10 −13 m −0.2 s −1 for m/n = 0.4, 0.5, and 0.6, respectively.Relief increases with the number of cells because the ridge singularity is resolved at finer resolution.

Two-dimensional model: numerical analysis of singular behavior
In Figs. 3 and 4 we present results which serve to distinguish the fundamental behavior of SPIM from the numerical behavior associated with a varying density of discretization.Figures 3 and 4 each show nine steady-state simulations, each using three values of M 2 and three values of m/n, i.e., 0.4, 0.5, and 0.6.In both figures, the number of cells is quadrupled from column to column.The leftmost column contains 40 2 cells, the middle column contains 80 2 cells, and the rightmost column contains 160 2 cells.Figure 3 shows simulations where the horizontal length scale, L 2-D , is held constant in all simulations.By increasing the number of cells, the grid size decreases.In all cases of m/n, the maximum relief increases with the number of cells.However, our quasi-theoretical analysis predicted the absence of an elevation singularity at the ridge for m/n < 0.5.To illustrate this point, we take a different approach, shown later in this section.Figure 4 contains simulations where grid size is held constant at 125 m.Here, the horizontal length scale, L 2-D , increases with the number of cells.In Fig. 4, the leftmost column contains 40 2 cells with L 2-D = 5 km, the middle column contains 80 2 cells with L 2-D = 10 km, and the rightmost column contains 160 2 cells with L 2-D = 20 km.Regardless of the m/n ratio and whether L 2-D or grid size is kept constant, the maximum relief of the landscape increases as the number of cells increases.Relief increases in both sets of simulations because with more grid cells, we are numerically sampling closer to ridges, and by sampling closer to ridges, we are resolving the ridge singularity on a finer scale.We emphasize, however, that the issue of dependence of the solution = 40 2 , 80 2 , and 160 2 .To make the relief of the landscapes comparable, the 2-D Pillsbury number, P 2-D , is set to 3.10 × 10 23 for solutions of all m/n ratios with L 2-D = 10 km.To achieve this for υ = 1 mm yr −1 , K = 6.31 × 10 −12 m 0.2 s −1 , 1.00 × 10 −12 s −1 , and 1.58 × 10 −13 m −0.2 s −1 for m/n = 0.4, 0.5, and 0.6, respectively.Like Fig. 3, relief increases with the number of cells because the ridge singularity is resolved at a finer resolution.on grid size is separate from the issue of scale invariance for m/n = 0.5, the latter result being deduced from the governing equation itself (Eq.15) before any discretization is implemented and illustrated in Fig. 2c.
Our quasi-theoretical analysis infers the conditions for singular behavior in the 2-D model.If elevation singularities exist, the model will not satisfy grid invariance, causing the relief between the ridge and outlet to increase indefinitely as grid size decreases.In contrast, in simulations where singularities do not exist, the relief between the ridge and outlet can be expected to converge as the grid size decreases.In both cases, understanding ridge behavior in the 2-D model requires studying solution behavior as grid size approaches 0.
We do this by extracting river profiles from 13 landscape simulations of different scales for each of three values of m/n, i.e., 0.4, 0.5, and 0.6.The largest simulation is for L 2 2-D = 10 6 km 2 ; simulations were also performed at progressively 1 order of magnitude less in area down to L 2 2-D = 10 −6 km 2 .The number of grid cells, M 2 , is held constant at 25 2 .In each simulation, then, the closest distance to the ridge that can be resolved is one grid cell, given by From each of the simulations, we construct two synthetic river profiles: one that intersects the highest point of the basin divide (high profile) and one that intersects the lowest point of the basin divide (low profile).The choice of these two www.earth-surf-dynam.net/5/807/2017/Earth Surf.Dynam., 5, 807-820, 2017 elevations was made so as to bracket the possible range of behavior; analogous results would be obtained from starting points along the basin divide at intermediate elevations.We use these synthetic profiles to characterize whether or not the numerical model is tending toward a singularity near ridges.We do this because the numerical model itself cannot directly capture singular behavior.We outline the details of the methodology for the high profile only, as the case of the low profile involves a transparent extension.
The 13 simulations result in 13 elevation profiles η i , where i = 1, 2... 13, each extending from l i (i.e., one grid point from the divide) to a downstream value l Di that is somewhat larger than the value 10 3−(i/2) km (because the downchannel path of steepest descent does not follow a straight line).We assemble a synthetic channel profile, η S (l), from these as follows.The first leg of η S (l) is identical to η 1 (l) and extends from l = l 1 to l D1 .We extend the synthetic profile by translating the second profile upward until its elevation at its downstream point l D2 matches with η S (l D2 ), as shown in Fig. 5a.The profile η S (l) now extends from l 2 to l D1 .As shown in Fig. 5a, we repeat this process until all 13 profiles have been used to assemble the synthetic profile, which now extends from l 13 to l D1 .
This procedure results in a high synthetic profile encompassing all 13 profiles (circles) and in a low synthetic profile (crosses) (Fig. 5b).One-dimensional analytical solutions, Eq. ( 10), are then fitted to the profiles of the 2-D simulations using the 1-D Pillsbury number, P 1-D , as a fitting parameter.To account for the difference in dimensionality, the 1-D steady-state profiles with hm/n = 0.8, 1.0, and 1.2 are fitted to the 2-D data for m/n = 0.4, 0.5, and 0.6, respectively.The scatter in the synthetic profile is due to the randomness in the pathway, as dictated by the initial conditions.
Figure 5b shows good fit between the 2-D results and the corresponding 1-D steady-state profiles.This allows us to make inferences concerning asymptotic behavior at a ridge.The analytical curves for elevation that best fit the 2-D data for m/n < 0.5 converge to finite values as l approaches 0 and infinity for m/n ≥ 0.5.While these results do not constitute analytical proof of this asymptotic behavior, they provide compelling evidence for it.

Scale behavior in other landscape evolution models
We offer here an example of a landscape model that does not necessarily satisfy horizontal-scale invariance, i.e., that of Gasparini et al. (2007).They incorporate the formulation of Sklar and Dietrich (2004) for bedrock abrasion due to wear in their model.The rate of erosion E is given as where K GA is the abrasion coefficient, Q s is the bed load sediment flux, W is the channel width, and Q t is the bed load transport capacity.Gasparini et al. (2007) use the following relation for Q t : where K t is a transport constant and m t and n t are exponents.At steady state, the total sediment flux at any point in the landscape must equal the production rate of sediment due to rock uplift: where K B is the fraction of sediment produced that contributes to bed load (the remainder being moved out of the system as wash load).For channel width, they use a relation of the form where Q is the water flow discharge, k w is the hydraulic geometry constant, b is the hydraulic geometry exponent (e.g., Finnegan et al., 2005).The value of b has been found to vary between 0.3 and 0.5 for bedrock rivers (Whipple, 2004); Gasparini et al. ( 2007) use b = 0.5 in their model.They also estimate discharge as an effective precipitation rate, k q , multiplied by a drainage area to the power of c, where c ≤ 1: The resulting relation for steady-state slope is Using the nondimensionalization terms from Eq. ( 12), we nondimensionalize Eq. ( 32) to where P G1 and P G2 are two dimensionless Pillsbury numbers: Horizontal-scale invariance results only when both dimensionless numbers are independent of the horizontal length scale, L 2-D .Gasparini et al. (2007) use m t = 1.5 and n t = 1.0.This parameter does indeed make the exponent 2 − 2m t + n t , equal to 0, so that P G1 is independent of L 2-D .The parameter P G2 is invariant to the horizontal scale when the product of b and c is equal to 1.However, realistic values of b are between 0.3 and 0.5 (Whipple, 2004) and the value of c is less than or equal to 1.This means that the maximum value of bc is 0.5.It follows that P G2 is not independent of the horizontal scale and that the model of Gasparini et al. (2007) does not satisfy horizontal-scale invariance.

Sensitivity of relief to hillslope length and profile length
In the river profiles of Figs. 1 and 5b, we see that a sizable proportion of the relief is confined to the headwaters, i.e., near a ridge.In our 1-D model, for hm/n ≥ 1, ridge elevation is infinite, thus formally implying infinite relief.This problem has been sidestepped by introducing a critical hillslope length l c , upstream of which it is assumed that there is no channel (e.g., Goren et al., 2014a).This point may be thought of as loosely corresponding to the channel-hillslope transition in the slope-area relation discussed by Montgomery andDietrich (1988, 1992).Here, then, we let the hillslope zone cover the range 0 ≤ l ≤ l c , where l c is an appropriately small fraction of profile length L 1-D .Modifying Eq. ( 10) accordingly, we can determine the total relief, R, of the channel profile as follows: where We remind the reader that according to Eq. ( 7), We now consider the scale-invariant case, hm/n = 1, and inquire as to how the relief of the basin might change.Increasing L 1-D does not increase relief because the parameter P 1-D ∼ L 0 1-D .It is thus seen from Eqs. ( 36) and ( 37) that relief can be increased only by decreasing lc .But from Eq. ( 36), R → ∞ as lc → 0. It follows that relief is extremely sensitive to the choice of lc .Based on our previous analysis, we expect that this result carries over to the case m/n = 0.5 for the 2-D model.We next provide an example illustrating the dependence of relief on hillslope length and profile length when hm/n = 1.Specifically, we consider the case hm/n = 0.9, with a dimensionless hillslope length lc = 0.01.According to Eq. ( 36), a halving of lc to 0.005 increases the relief by 11.4 %.In order to achieve the same increase in relief by changing profile length L 1-D while holding lc constant, L 1-D would have to be increased by 196 %.It is thus seen that relief of the channel profile can be more sensitive to a relative change in dimensionless critical channel length than it is to a relative change in horizontal scale.

Discussion and conclusion
Our 1-D analytical solutions, Eq. ( 10) and Fig. 1, characterize the scale behavior of 1-D SPIM, with horizontal-scale invariance satisfied when hm/n = 1.0.Our 2-D numerical solutions shown in Fig. 2 illustrate our analytical result that 2-D SPIM shows horizontal-scale invariance when m/n = 0.5.That is, 2-D models using SPIM with m/n = 0.5 show the same relief structure regardless of the horizontal scale.This scale invariance has been previously demonstrated for neither the 1-D nor the 2-D SPIM model.Our result calls into question the common usage of the ratio m/n = 0.5 in landscape evolution models (Gasparini et al., 2006).For example, the Python-based landscape modeling environment, Landlab (Hobley et al., 2017) offers a default m/n ratio of 0.5.Our result also motivates further investigation as to why analysis of field data commonly yields values of m/n ∼ 0.5 (e.g., Snyder et al., 2000).It should be noted that local empirical measurements indicating m/n = 0.5 do not necessarily mean m/n = 0.5 should be used as a universal ratio in SPIM.Gasparini and Brandon (2011) used multiple incision laws, other than SPIM, to simulate steady-state landscapes and were able to fit E, A, and S in Eq. ( 1) to find empirical values of m and n (prime denotes an empirical value).They found that the ratio of m /n was sensitive to the incision model's parameters as well as the rock uplift pattern in each landscape.This implies that both m and n have dependency on landscapes properties and are not universal from landscape to landscape.
In addition to the horizontal-scale invariant case m/n = 0.5 for the 2-D SPIM model, we also emphasize the relationship between the steady-state landscape relief and horizontal when m/n = 0.5.Equations ( 14) and ( 15) and the results in Fig. 2c show that the relief structure of the landscape scales with P 2-D .Within P 2-D , the horizontal length scale term is L 1−2m/n 2-D .For the m/n ratio range 0.35 to 0.6 (Whipple and Tucker, 1999), the corresponding exponent range in the horizontal length scale term is −0.2 to 0.3.This means that over the stated range of m/n, the relief structure has a weak dependence on the horizontal length scale.For m/n < 0.5, relief weakly increases with horizontal scale.
For m/n > 0.5, relief weakly but unrealistically decreases with horizontal scale.The underlying physics of channel and hillslope processes that might dictate such behavior are, at present, unhelpfully opaque.In natural systems, larger landscapes would yield longer rivers.Since elevation monotonically increases with upstream distance, one would expect relief to increase with horizontal scale.The results of SPIM, where m/n ≥ 0.5, clearly contradict this intuitive understanding.
Our work neglects the effect of hillslope diffusion because our intent is to study the behavior of SPIM itself.Without hillslope diffusion, SPIM causes singular behavior at ridges in both the 1-D and 2-D formulation.Indeed, both the 1-D and 2-D models exhibit singularities in slope at ridges for all hm/n ratios (1-D) and all m/n ratios (2-D).For hm/n ≥ 1 (1-D) and m/n ≥ 0.5 (2-D), the models exhibit singular behavior in elevation at ridges as well.When relief is limited by a hillslope length l c , elevation and slope do indeed reach finite values at the channel heads, but the effects of the singularity persist.For example, for the case hm/n ≥ 1 in the 1-D model, relief approaches infinity as hillslope length approaches 0. Our analysis of ridge singularities in SPIM shows that the choice of hillslope parameterization plays a key role in determining the relief of natural landscapes.
Numerical solutions of the 2-D model indicate that it cannot be grid-invariant for m/n ≥ 0.5.In the absence of hillslope diffusion, ridges reach infinite elevation as grid size becomes vanishingly small.This result underlines the critical role of hillslope diffusion in obtaining meaningful results from the 2-D model.Field estimates of hillslope diffusion have been obtained on the hillslope scale, but there are unanswered questions about their application to large-scale models (Fernandes and Dietrich, 1997).Our results suggest that for the ratio m/n < 0.5, there are steady-state grid-invariant solutions.However, the grid size below which grid invariance is realized may be so small, e.g., sub-meter scale, that the validity of Eq. ( 1) is called into question.Issues with SPIM when used on large scales include the following.Studies commonly neglect the effect of hillslope diffusion when the scale of the grid is larger than the hillslope scale (Somfai and Sander, 1997;Banavar et al., 2001;Passalacqua et al., 2006).On coarse-grained scales, increasing the size of the numerical domain, while keeping the number of cells constant, will result in the behavior shown in Fig. 2. In Fig. 4 we see that adding more cells to compensate for the increase in size of the domain, such that the grid size remains constant, produces heavily biased (i.e., ever more singular) behavior near the ridges.
Our analysis illustrates that SPIM has two important limitations: (a) unrealistic scale invariance when m/n takes the commonly used value 0.5, so that a 10 km 2 basin has identical relief to a 1000 km 2 basin, and (b) singular behavior near the ridges for m/n ≥ 0.5 that makes maximum relief entirely and unrealistically dependent on grid size.SPIM has been shown to be of considerable use in the study of the general Earth Surf.Dynam., 5, 807-820, 2017 www.earth-surf-dynam.net/5/807/2017/

Figure 3 .
Figure 3. Nine 2-D numerical simulations at steady state for three different values of M 2 and three different values of m/n.In this figure, the horizontal scale is kept constant: L 2-D = 10 km for all solutions.The value of K has been chosen to be different for each value of m/n for clarity in the figures.From left to right, the number of cells M 2= 40 2 , 80 2 , and 160 2 .To make the relief of the landscapes comparable, the 2-D Pillsbury number, P 2-D , is set to 3.10 × 10 23 for solutions of all m/n ratios with L 2-D = 10 km.To achieve this for υ = 1 mm yr −1 , K = 6.31 × 10 −12 m 0.2 s −1 , 1.00 × 10 −12 s −1 , and 1.58 × 10 −13 m −0.2 s −1 for m/n = 0.4, 0.5, and 0.6, respectively.Relief increases with the number of cells because the ridge singularity is resolved at finer resolution.

Figure 4 .
Figure 4. Nine 2-D numerical simulations at steady state for three different values of M 2 and three different values of m/n.In this figure, the grid size, L 2-D /M = 125 m, is used in all the solutions.The value of K has been chosen to be different for each value of m/n for clarity in the figures.From left to right, the number of cells M 2= 40 2 , 80 2 , and 160 2 .To make the relief of the landscapes comparable, the 2-D Pillsbury number, P 2-D , is set to 3.10 × 10 23 for solutions of all m/n ratios with L 2-D = 10 km.To achieve this for υ = 1 mm yr −1 , K = 6.31 × 10 −12 m 0.2 s −1 , 1.00 × 10 −12 s −1 , and 1.58 × 10 −13 m −0.2 s −1 for m/n = 0.4, 0.5, and 0.6, respectively.Like Fig.3, relief increases with the number of cells because the ridge singularity is resolved at a finer resolution.