THE SCALING OF FLUVIAL LANDSCAPES

Björn Birnir ^{‡} Terence R. Smith ^{†}^{∗} George E. Merchant ^{†}

‡ Department of Mathematics, University of California at Santa Barbara, California 93106, USA and University of Iceland, Science Institute, 3 Dunhaga, 107 Reykjavík † Department of Geography, University of California at Santa Barbara, California 93106, USA ∗ Department of Computer Science, University of California at Santa Barbara, California 93106, USA

ABSTRACT

The analysis of a family of physically-based landscape models leads to the analysis of two stochastic processes that seem to determine the shape and structure of river basins. The partial differential equation determine the scaling invariances of the landscape through these processes. The models bridge the gap between the stochastic and deterministic approach to landscape evolution because they produce noise by sediment divergences seeded by instabilities in the water flow. The first process is a channelization process corresponding to Brownian motion of the initial slopes. It is driven by white noise and characterized by the spatial roughness coefficient of 0.5. The second process, driven by colored noise, is a maturation process where the landscape moves closer to a mature landscape determined by separable solutions. This process is characterized by the spatial roughness coefficient of 0.75 and is analogous to an interface driven through random media with quenched noise. The values of the two scaling exponents, which are interpreted as reflecting universal, but distinct, physical mechanisms involving diffusion driven by noise, correspond well with field measurements from areas for which the advective sediment transport processes of our models are applicable. Various other scaling laws, such as Hack’s law and the Law of Exceedence Probabilities, are shown to result from the two scalings, and Horton’s Laws for a river network are derived from the first one.

1 Introduction

In three preceding papers [Smith et al., 1997a, Smith et al., 1997b, Smith et al., 2000], we showed that a family of landscape evolution models capture, at least to a first approximation, the emergence and development of stable, dendritic patterns of valleys and ridges. In this paper we demonstrate the manner in which these models also capture the effects of random influences in driving the processes of landscape evolution. In particular, our results provide a physical basis for explaining various fundamental scaling relationships that characterize fluvial landscapes and suggest a bridge between deterministic and stochastic theories of drainage basin evolution.

An important criticism of the discrete stochastic models is that, despite their impressive simulations of many planar features of fluvial landscapes, they contain no acceptable bridge to the more physically acceptable theories based on PDEs. The complexity of fluvial processes and landforms, however, has long provided a justification for models that incorporate stochastic elements. The growing body of research focused on scaling relationships in landscapes, recently summarized in [Rodriguez-Iturbe and Rinaldo, 1997], emphasizes the importance of integrating deterministic and stochastic representations of fluvial processes into a tractable and useful theory. A general goal of this paper, therefore, is to show that by interpreting PDE-based models of landscape evolution as driven by random influences, one obtains significant insights not only into the nature of landscape scaling laws, but also into the manner in which the deterministic and stochastic lines of investigation may be bridged.

In relation to our modeling context, we note that the models described in Section 2 below focus on the advective entrainment and transport of sediment in transport limited conditions as the dominant process in drainage basin evolution. While other processes, such as detachment limited entrainment and gravity driven diffusion, are important processes in many geomorphic contexts, the modeling results presented in [Smith et al., 1997a, Smith et al., 1997b] indicate that advective entrainment and transport processes are in some sense the fundamental processes underlying drainage basin evolution. We believe that isolating and understanding the effects of these processes is a necessary and critical step in understanding the effects of many interacting processes on landscape development.

In relation to our analytical context, we apply the methods developed in a large and rapidly growing body of research focused on explaining the dynamics and forms of complex surfaces and interfaces evolving in a variety of physical contexts. Important goals of this literature are to characterize the qualitative behavior of high- and infinite-dimensional, nonlinear systems that are driven by noise and to explain the origin of temporal and spatial scaling behavior in a wide variety of phenomena [Schroeder, 1991, Zipf, 1949, Willis, 1924, Mandelbrot, 1983]. It is now known that the evolution of surfaces whose dynamics are driven by various forms of noise are frequently characterized by scaling laws. Such laws often, but not always, indicate that the system is insensitive to differences in the details of the underlying mechanisms and processes. This insensitivity leads to the observation that seemingly different systems exhibit qualitatively similar, or universality classes of, behaviors (see, for example, Sneppen [Sneppen, 1998].) This fact is of great value insofar as it allows one to extend results that are known about the response of one system to noisy inputs to other systems in the same universality class. In the following analysis, we make use of such extensions in understanding how the two types of noise (homogeneous and quenched) lead to distinct emergent scalings in the eroding surfaces described by our models.

We employ several specific techniques from the emerging theory of complex surface evolution in investigating the models discussed in [Smith et al., 1997a, Smith et al., 1997b] as systems driven by noise or stochastic processes. First, we characterize the statistical structure of eroding surfaces and flows in terms of various structure functions that represent the statistical correlation structure of complex surfaces. Second, we apply known results from this theory concerning the form of scalings that emerge from appropriate universality classes of PDEs when subjected to random driving forces of a specific form. As noted above, the rationale for such application is that systems belonging to the same universality class manifest qualitatively similar behaviors. In order to take advantage of these known results we employ suitable approximations of our models in conjunction with specific models of noise that we justify in terms of the basic properties of our models. We also connect part of our theory to the concept of self-organized criticality (SOC) as proposed by Bak et al. [Bak et al., 1988, Bak and Paczuski, 1995, Paczuski, 1995, Sneppen, 1998]. Although this concept has yet to find an appropriate mathematical formulation, it is based on the underlying idea of an attractor poised at criticality to which a system is driven by random perturbations. Such systems are characterized by essential relationships between the spatial and the temporal scalings of events [Sapozhnikov and FouFoula-Georgiou, 1996] and by specific asymptotic forms of the structure functions^{1}.

In relation to the statistical structure of our eroding surfaces and flows, our main results are the discovery of two fundamental scaling relationships that are represented in terms of the asymptotic forms of a structure function characterizing the width, or height-height correlation function (see section 2.3 below), of eroding surfaces. A first scaling relation, with an exponent of 0.5 characterizing the depth of water flow and the slope of the water surface, emerges early. This scaling is fundamentally associated with the development of channelized flows on unstable surfaces and indicates the importance of stochastic influences in determining the configuration of emerging channel networks. A second scaling, with an exponent of approximately 0.75, is associated with self-similar separable solutions which are associated with the mature stages of development and the emergence of a variational principle. This scaling characterizes the form of the land surfaces and its slope. The differences between the two scalings suggest that they are determined by different mechanisms.

1We believe, in fact, that our theory may provide some of the foundation for a mathematical theory of self-organized criticality.

We note with great interest that field values of the same height-height correlation function obtained for landscapes in Ethiopia, Somalia, and Saudi Arabia by Weissel and Pratson [Weissel and Pratson, 1994], fall neatly into the range 0.5-0.7. Although empirical scaling results of this form tend to cluster around a relatively small set of values, these results provide some validation for our theory, particularly since our models may be interpreted as representing conditions under which our advective models of sediment entrainment and transport apply, since the diffusive mass-wasting processes of humid landscapes are relatively unimportant. We may therefore interpret our models as providing physically-based explanations for such empirical scaling relationships.

In relation to the application of known results from the theory of complex surfaces to explain the physical mechanisms underlying these scalings, we first construct appropriate linearizations of our models in the two relevant regimes: the early channel initiation regime and the mature, separable surface regime. We show that, when viewed at their fast time scales, these linearized models possess noise generating mechanisms. These mechanisms act as seeds of the noise that drives the erosion equations at their long time scales and thereby produce the observed scalings. We may interpret the initial, small, random variations as characteristic of the eroding elevations and of the rainfall process^{2}.

The mechanisms that generate this driving noise are important since they provide a linkage between the observed scalings and the physical mechanisms underlying the evolution of the eroding surfaces. We identify two plausible mechanisms of noise generation that stem from instabilities in the flow of water over the eroding surfaces and that take place on the short time-scales associated with such flows. A mechanism that we associate with the regime of early channel initiation involves the formation of shocks or hydrolic jumps ^{3} in the flow of water over locally convex portions of landscapes in which channels are first emerging. These shock formations are of two types. The first is associated with nonlinearities in the water flow, the second with both nonlinearities in the sediment flow and the linear propagation of singularities in the water flow. As we show below, noise generated by shocks of both types may be viewed as seeds for the instabilities in the sediment flow that drive the process of channelization. A second mechanism of noise generation that we associate with the regime of mature, separable surfaces again takes place on a short time-scale and involves the formation of shocks on characteristically concave surfaces at the points where the derivatives of the surface elevation function are singular, and hydrolic jumps on the lower reaches of the slopes and in the valleys. While we emphasize that these mechanisms of nonlinear noise generation are specific to our simple advective models, mechanisms involving analogous effects almost certainly exist in more complex models. The details of how the noise generating mechnisms make the water flow turbulent and lead to divergences in the sediment flow are not completely understood. These mechanism and how to model them mathematically are the topic of current investigations. We will elaborate below on the intuitive reasoning that instabilities in the water flow act as seeds to the divergences in the sediment flow and that the resulting mathematical model for the sediment transport must be a noise-driven gradient flow.

2In computing solutions to our basic model, we simulate such noise as a combination of perturbations to the initial eroding surface and the round-off error in the numerical computations.

3A shock is a discontinuity in a wave-like form that occurs when the higher parts of the wave propagate at a higher velocity than the lower part (see, for example, Smoller [Smoller, 1994]).

The cumulative effects of the noise generated in the flow of water at short time-scales impact the sediment flows at the longer time-scale of erosion and lead to the observed scalings. Hence we view the first scaling relation, which arises from an erosion process driven by noise, generated by water-seeded divergences in the sediment flow, as resulting from universal phenomena that are modeled as diffusive (random) processes driven by homogeneous noise. Since this scaling is associated with the development of channelized flows on unstable surfaces, it indicates the importance of stochastic influences in determining the configuration of emerging channel networks. In particular, we show that such configurations may be viewed as arising from a random walk. The second scaling relation, which arises from erosion processes being seeded by instabilities in the water flow but over separable surfaces with singularities [Smith et al., 1997a], is also a universal phenomenon generated by a noise-driven diffusive process. Whereas the noise generated by shocks (hydrolic jumps) over the initial channelizing surfaces is spatially homogeneous, the noise seeded by the shock-inducing singularities of the separable surfaces is highly inhomogeneous in space due to nonlinear quenching. Thus the sediment divergences are seeded by homogenous instabilities in the water flow in case of the first process but colored instabilites in case of the second process. In addition the sediment divergences themselves differ in the two cases. Hence, although the two processes are similar in both being diffusive and noise-driven, the differences in forms of the associated noise result in differing scaling exponents.

We show that various other scaling relations follow from our two basic scaling relations. In particular, we show that the second scaling leads to a simple derivation of the established scaling relationship between stream length and drainage area known as Hack’s Law. Furthermore, the noise that leads to the first scaling also drives the initial growth of channels, allowing us to characterize, in terms of random walks, the nature of the process by which channels emerge from unstructured surfaces. We emphasize that our contributions in relation to our scaling results are not simply in showing how one particular model gives rise to (self-similar) scaling, but to show that the scaling of surfaces and processes in a large class of relevant models may be viewed as being caused by universal physical processes.

Many details of this theory remain for investigation and significant efforts will be required to generalize the theory to cases involving other geomorphic processes. Specific questions of interest include whether the two suggested mechanisms of noise generation can be generalized to other, more comprehensive models; whether the scaling results can be shown to follow in the full nonlinear regime of the equations; and whether the two scalings are sufficient to determine all of the other scaling relations of fluvial landscapes, or whether there are other scalings that we have either yet to find or that arise as a result of processes not encompassed within our advective models. Despite these daunting research problems, we believe that the results presented in this paper provide a reasonable qualitative model of the emergence of form and scalings over a broad range of fluvial contexts and that such results are likely to be independent of many specific modeling details.

The paper is structured as follows. In Section 2 we discuss the mathematical models, the nonlinear shock formation that seeds the noise, and the structure functions that give rise to the scaling relations. The results of numerical simulations of the mathematical models, which are presented in Section 3, provide a basis for the empirical and analytical investigations presented in the remaining sections of the paper. In Section 4, we provide a theoretical explanation for the first scaling relation that is associated with the emergence of channelization on unstructured surfaces. In Section 5, we provide a theoretical explanation for the second scaling relation that is associated with the separable landscapes. We discuss the implications of our results for the understanding of geomorphic phenomena in Section 6, while Section 7 contains a summary of our results. A derivation of the formation of hydrolic jumps (shocks) in the water flow over concave surfaces in contained in Appendix A. A proof of the formation of shocks from surface singularities and their linear propagation, and an explanation for the observed scaling of the gradient of the water surface are contained in Appendices B and C respectively and a derivation of a relationship between Horton’s laws is contained in Appendix D.

2 Modeling Eroding Surfaces and their Scaling Laws

Because the surfaces generated by fluvial erosion are the result of highly nonlinear processes driven by noisy inputs, they are complex and difficult to represent in traditional terms. We now discuss the models of nonlinear erosion processes from which our results are generated, and briefly indicate the nature and origin of the random effects that drive the equations. We then describe the manner in which such surfaces may be represented in terms of structure functions, the asymptotic forms of such structure functions, and the associated scaling laws.

2.1 A Family of Landscape Evolution Models

The scaling results presented in this paper are derived from the same family of models studied analytically and numerically by Smith, Birnir, and Merchant [Smith et al., 1997a], Smith, Merchant, and Birnir [Smith et al., 1997b], and Birnir, Smith, and Merchant [Birnir et al., 1998]. These models represent the advective entrainment and transport of sediment in transport limited conditions [Howard, 1994]. Based on our previous analyses and on the results reported in this paper, we believe that they capture essential aspects of fluvial erosion at small to medium scales of spatial resolution (see, for example, Howard, Dietrich, and Seidl [Howard et al., 1994]) ^{4}. While it is straightforward to extend these models to represent other processes that are significant in the evolution of fluvial landscapes, such as gravitydriven diffusion processes on hillslopes, we believe that their effects would mask the scaling relations associated with the advective processes of sediment transport in channel and overland flows^{5}.

We now provide a summary description of the family of models. A more detailed description of their derivation is provided in [Smith et al., 1997a]. We also provide a brief discussion of those aspects of the characteristic time scales of the models that are relevant to the present context. The models are based on conditions describing the conservation of water and sediment fluxes over a continuous, erodible surface z = z(x,y,t)

− R + \frac{∂h}{∂t} = ∇ · (u_{w} q_{w}),  (1)

\frac{∂z}{∂t} = ∇ · (u_{w} q_{s})  (2)

in which h = h(x,y,t) is the depth of water varying continuously over the landsurface, −u_{w} = −\frac{∇H}{|∇H|} is a unit vector in the direction of both water and the advected sediment flows, H = H(x,y,t) = z(x,y,t) + h(x,y,t) is a free water surface, q_{w} represents the flux of water per unit width, and q_{s} represents the advected flux of sediment per unit width.

There are three natural time scales that characterize the dynamics of equations (1), (2). These time scales may be derived by transforming the variables of equations (1), (2) to dimensionless form, using relations v = [v]v^{∗} in which [v] is a scale and v^{∗} a dimensionless variable. Given values for the scaling parameters [H],[h],[q_{w}],[q_{s}],[R] and [x] = [y] that are characteristic of the variables in equations (1), (2),^{6} we may then define representative values of the scaling parameter for the time variable ([t]) by setting the values of specific dimensionless parameters of the problem to unity.

4Whether they capture essential aspects at large scales of spatial resolution is an open issue.

5We plan to explore the effects of such processes on scaling in later papers.

6A useful set of such values characterize the representative scales when the basic model (1, 2) is used to represent the erosion of a small, linear ridge.

A short time scale [t] for t may be defined by the relation [t][q_{w}]/[h][x] = 1, which may be interpreted to imply that [t] is the time scale at which the volume of water flowing over a lateral cross section of the ridge is approximately the same as the volume of water on the surface. Applying this scaling to dimensionless versions of equations (1), (2) and dropping the asterisks denoting dimensionless variables, we obtain

\frac{∂h}{∂t} = ∇ · (u_{w} q_{w})+ \frac{[x][R]}{[q_{w}]} R
\frac{∂H}{∂t} − \frac{[h]}{[H]}\frac{∂h}{∂t} = \frac{[h][q_{s}]}{[H][q_{w}]} ∇ · (u_{w} q_{s}).

It is natural to choose a scaling that makes the dimensionless parameter [x][R]/[q_{w}] = 1 in the water flow equation and interpret all terms in this equation as being O(1) (see, for example, [Fowler, 1997]). The significance of this fact, as we show in Section 2.4 and in Appendix A, is that the flows described by this equation are characterized by shock waves that we may interpret as short term noise. The order of the dimensionless parameter [h][q_{s}]/[H][q_{w}] in the erosion equation is much smaller than unity and smaller than [h]/[H]. Hence we obtain the (approximate) erosion equation

\frac{∂H}{∂t} ≈ \frac{[h]}{[H]}\frac{∂h}{∂t}

which indicates that variations in the water surface H are small at this time scale and driven by short term variations in the flow depth h described by the water flow equation. This equation also implies that erosion of the surface is negligible at this time scale.

A time scale [t] that is intermediate between the previous short scale and longer scales may be defined by the relation [t][q_{s}]/[h][x] = 1, which implies that we may interpret [t] as a time scale at which the volume of sediment flowing over a lateral cross section of the ridge is approximately equal to the volume of water on the ridge at any time. On applying this scaling to equations (1), (2) we obtain

\frac{[q_{s}]}{[q_{w}]}\frac{∂h}{∂t} = ∇ · (u_{w} q_{w})+ R
\frac{∂H}{∂t} = \frac{[h]}{[H]}∇ · (u_{w} q_{s})+ \frac{[h]}{[H]}\frac{∂h}{∂t}

where we have once more used the scaling relation [x][R]/[q_{w}] = 1. We note that the dimensionless parameter [q_{s}]/[q_{w}] that multiplies the term \frac{∂h}{∂t} has a magnitude that is greater than its magnitude at the short time scale. Since the evolution of the water surface H is determined by two terms that we may assume to be O(1), see [Fowler, 1997], and since both are multiplied by the same small dimensionless parameter [h]/[H], we may interpret this to mean that changes in the surface H are determined to a significant degree by small fluctuations in the depth of water h, indicated by the water flow equation, as well as by the erosion of the z surface.

A long time scale [t] may be defined by the relation [t][q_{s}]/[H][x] = 1 and we may interpret

[t] as the time scale at which the volume of sediment flowing over a lateral cross section of the ridge represents a significant proportion of the volume of the ridge. On applying this scaling to equations (1), (2) and again using the scaling relation [x][R]/[q_{w}] = 1, we obtain

\frac{[h][q_{s}]}{[H][q_{w}]} ; \frac{∂h}{∂t} = ∇ · (u_{w} q_{w})+ \frac{[x][R]}{[q_{w}]} R
\frac{∂H}{∂t} = ∇ · (u_{w} q_{s})+ \frac{[h]}{[H]}\frac{∂h}{∂t}.

The small size of the dimensionless parameter on the LHS of the water flow equation, [h][q_{s}]/[H][q_{w}] << 1, suggests that the basic flow is essentially an equilibrium flow down the surface gradient. It is not unreasonable, therefore, to adopt the following approximation to the water flow relation

0 ≈ ∇ · (u_{w} q_{w})+ \frac{[x][R]}{[q_{w}]} R  (3)

as was done, for example, in the analyses presented in [Smith and Bretherton, 1972, Loewenherz-Lawrence, 1991, Smith et al., 1997a, Smith et al., 1997b]. While we may interpret the dominant balance in the erosion equation as being between the first two terms, we may interpret the relatively small term [h]/[H]∂h/∂tas representing random variations that drive the sediment flow.

2.2 Previous Results for Long Time Scale Approximations

Smith, Birnir, and Merchant [Smith et al., 1997a, Smith et al., 1997b] analyzed a subfamily of the models (1), (2) that was obtained with the use of the following assumption. First, the long time scale approximation (3) to the basic model was assumed, together with the additional scaling relationship [h]/[H] = [q_{s}]/[q_{w}] ≡ η. Second, a Manning-type constitutive relation

q_{w} = nh^{\frac{5}{3}}|∇H |^{\frac{1}{2}}  (4)

where n, the inverse of the channel roughness, is a constant, was adopted to represent the flux of water. Third, a power law relation

q_{s} = kh^{\frac{5}{3}γ}|∇H |^{\frac{β}{2}+δ}  (5)

was employed in representing the flux of sediment q_{s} [Julien and Simons, 1985]. These assumptions lead to the following dimensionless equations:

−R = ∇ · [\frac{∇H}{|∇H |} h^{\frac{5}{3}} |∇H |^{\frac{1}{2}}]  (6)

\frac{∂H}{∂t} − η\frac{∂h}{∂t} = ∇ · [\frac{∇H}{|∇H |} h^{\frac{5}{3}γ}|∇H |^{\frac{γ}{2}+δ} ].  (7)

The initial and boundary conditions employed in these analyses modeled a linear ridge undergoing erosion as a result of a uniform rainfall (see equations (9), (10) below.)

Applying the numerical methods described in [Smith et al., 1997b] to equations (6)-(7), it was shown [Smith et al., 1997b] that initially unstructured but randomly perturbed twodimensional surfaces are unstable when eroded by water. Channelized flows develop on the surfaces with a region of maximum channelization first emerging towards the lower boundary of surfaces with an initially planar configuration^{7}. This process of channelization involves the merging of “rivulets” and the formation of larger “channels”. The longitudinal profiles of the ridge simultaneously develop a concavity that also emerges initially near the base of the ridge under the influence of the fixed lower boundary. The region characterized by the maximum rate of channelization and by the emergence of the concavity gradually moves towards the upper boundary of the ridge as the surface erodes.

A pattern of dendritic valleys separated by ridges emerges and comes to dominate the eroding surface. Such a pattern is illustrated in Figures 1 and 2, which represent the surface after 10% of an initially planar surface has been eroded. After a characteristic period of erosion, the surfaces evolve towards stable landscapes of concave valleys and ridges that decay slowly in a self-similar manner. Such a pattern is illustrated in Figures 3 and 4, which represent the surface after 60% of an initial planar surface has been eroded^{8}.

The convergence to these forms occurs from initial surfaces that are both planar and nonplanar. The characteristic period of erosion that precedes this convergence, which we measure in terms of the percentage of the original surface eroded, depends on the nature of the initial surface, the rainfall rate, and the parameters as γ and δ. These mature landscapes, consisting of stable patterns of concavities, valleys, ridges, and associated flows, are welldescribed by a class of solutions to the nonlinear PDEs (6)-(7) for cases in which η \frac{∂h}{∂t} may be ignored in equation (7) relative to \frac{∂H}{∂t} [Smith et al., 1997a]. Such solutions depend on the characteristics of the water and sediment transport laws and on the boundary conditions for the PDEs; they are separable in time and space

7A variety of numerical experiments have shown that channelization occurs on surfaces with a large array of initial configurations.

8We note, in particular, the sharpness of the longitudinal ridges, or divides, separating the valleys.

H(x,y,t) = T_{H}(t)H_{o}(x,y), h(x,y,t) = T_{h}(t)h_{o}(x,y);  (8)

and they are stable [Smith et al., 1997a]. They are also characterized in terms of a variational principle [Smith et al., 1997a], by which a simple function of the sediment transport over the surface is minimized, subject to constraints involving the conservation of water flow and the elevations of the initial surface.

2.3 The Models Underlying the Present Analysis

The analysis presented in the current paper is similarly based on the use of equations (1), (2), the Manning-type constitutive relation (4) for the flux of water, a power law relation (5) for the flux of sediment, and the initial and boundary conditions prescribed in equations (11), (12). For reasons that we discuss below, however, we no longer assume that the term \frac{∂h}{∂t} in equation (1) is always negligible. As we show this term leads to significant disturbances in the flow of water at the short time scales discussed in Section 2.1. We interpret these disturbances as noise, whose cumulative effects, as we discuss below, are significant in seeding the instabilities in the flow of sediment at the longer time scales.

We therefore focus our attention on the model

η^{2} \frac{∂h}{∂t} = ∇ · [\frac{∇H}{|∇H |} h^{\frac{5}{3}} |∇H |^{\frac{1}{2}}] + R,  (9)

\frac{∂H}{∂t} − η\frac{∂h}{∂t} = ∇ · [\frac{∇H}{|∇H |} h^{\frac{5}{3}γ}|∇H |^{\frac{γ}{2}+δ} ]  (10)

in which we have adopted the scaling relation [h]/[H] = [q_{w}]/[q_{s}] ≡ η. In relation to the three time scales discussed above, we note that there is no loss of generality in this representation. In particular, we may obtain the short time scale version of the equations by applying the transformation τ = t/η^{2} and the intermediate time scale version by applying the the transformation τ = t/η.

We employ the same boundary and initial conditions as [Smith et al., 1997a, Smith et al., 1997b] to model a linear ridge extending uniformly in the lateral (x)-direction and defined over a rectangular domain of length L and width W,

D = {(x,y) ∈ IR^{2}|0 ≤ x ≤ L, 0 ≤ y ≤ W },

with initial conditions

h(x,y,0) = d(y), d(0) = h_{o}, d(W) = 0,
H(x,y,0) = cy + h_{o}, 0 ≤ y ≤ W  (11)

and boundary conditions

h(x,W,t) = 0,
H(x,0,t) = h_{0} = h(x,0,t)  (12)

corresponding to a water depth of zero at the top of the ridge and an absorbing body of water at the base of the ridge. While the water surface must be considered to be a free surface at the top of the ridge, it may be viewed as consisting of finitely many smooth curves that are solutions of a nonlinear ODE (the PDE restricted to the boundary). These curves are joined in a continuous, but not smooth, moving boundary (see, for example, Figure 3.) The upper boundary is characterized by the additional conditions

q_{w} = q_{s} = 0,

indicating the absence of any flux of water or sediment over this boundary. Since q_{w} and q_{s} are expressed as powers of h and ∇H in the constitutive relations, these conditions imply that the vanishing of the water depth dominates the blow-up of the gradient of H (in q_{w} and q_{s}) and that the normal derivative of H may become infinite at the upper boundary

n · ∇H(x,w,t) = ∞.

The boundary conditions on the lateral boundaries of the ridge at x = 0 and x = L are taken to be periodic, modeling a linear ridge of infinite extent.

2.4 The Role of Noise and Its Generation

The reasons for retaining the term η^{2} \frac{∂h}{∂t} in the water flow equation (9) is related to the issue of noise, its effect in driving the erosion process, and its relation to the scaling relations that characterize eroding land surfaces. Many investigations of PDEs that model the evolution of complex surfaces and interfaces have shown the significant role played by noisy inputs in determining the behavior of such surfaces. As shown, for example, by Sneppen [Sneppen, 1998] the nature of the noise driving the evolution of complex surfaces and interfaces is often a major determinant of the structure, and hence the scaling properties, in both emergent and mature systems.

It is especially essential to consider the effects of noise when modeling the evolution of geomorphic phenomena. The complexities of geomorphic processes and surfaces make it virtually impossible to model many aspects of landscape evolution without representing various classes of model inputs in stochastic terms. The nonlinear properties of such systems may lead to invalid analyses if explicit account is not taken of the stochastic nature of various inputs. This is especially the case for landscape phenomena whose evolution is governed by complex, interdependent processes characterized by differing time scales. Within the context of our problem, for example, it is difficult to imagine representations of either small scale (and hence short term) perturbations to the landsurface z(x,y,t) or short term variations in rainfall R that are not stochastic. Given such inputs, the flow of water over the surface z(x,y,t), as modeled by equation (9), is itself assured to be stochastic, as are any short term effects of such variations on the entrainment and transport of sediment. As we show below, however, random effects such as these that occur over relatively short time scales can neither be ignored nor averaged out at longer time scales, since their effects may accumulate in non-trivial ways.

Plausible sources for random inputs to the models presented in [Smith et al., 1997a, Smith et al., 1997b] and represented in (9), (10) include, at the least, (1) the form of the initial surface; (2) the structure of the eroding substrate; and (3) the inputs of rainfall. In the investigations described in [Smith et al., 1997b], and again in the present analysis, the initial surfaces from which channelized forms evolve are randomly perturbed. Not only is this necessary for breaking the symmetry of the surfaces, and hence allowing channel-forming instabilities to become effective, but accords with the preceding argument on the impossibility of modeling, in a deterministic manner, the processes and forms leading to small surface perturbations.

In the investigations described in [Smith et al., 1997b], and again in the present analysis, neither the eroding substrate nor the rainfall inputs are represented in terms of stochastic effects. It is, therefore, valid to ask whether this apparent lack of time-dependent stochastic variation in the model inputs leads to any significant lack of generality in our modeling results. We do not believe this to be the case. First we note compelling empirical evidence from a variety of studies (see, for example, Morisawa [Morisawa, 1964] and Schumm, Mosley, and Weaver [Schumm et al., 1987]) indicating that the configuration of emergent drainage patterns is strongly determined by initial, and essentially random, surface configurations.

Second, and perhaps of greater intrinsic interest, is the fact that the mechanisms of water flow embedded in equation (9) are capable of generating variations in flows at short time scales that may be interpreted as constituting noise. In particular, we note that equation (9) is capable of generating shock waves that provide a continuous source of variation to the system. It is, for example, well-known that flows described by the analogous kinematic wave approximation of the Saint-Venant equations are characterized by shocks [Dingman, 1984, Singh, 1996]. The kinematic wave approximation is known to hold under conditions of overland flow on steep slopes and for variations in discharge that are of small amplitude or long wavelength (see, for example, [Dingman, 1984, Singh, 1996]).

In Appendix A we show that the water flow equation (9) becomes a one-dimensional equation describing the flow of water down the gradient of the water surface

\frac{∂h}{∂t} − \frac{5}{3η^{2}} h^{2/3}|∇H |^{\frac{1}{2}} \frac{∂h}{∂s} = [ \frac{|∇H |^{2}ΔH − \frac{1}{2}(H_{x}^{2}H_{xx}+ 2H_{x}H_{y}H_{xy}+ H_{y}^{2}H_{yy})}{|∇H |^{\frac{5}{2}}} ] \frac{h^{5/3}}{η^{2}}
+ \frac{R}{η^{2}}.  (13)

in which s parametrizes this direction. We provide a proof of the existence of shocks in the flows described in this equation in Appendix A, and show below that such shocks lead to white noise that seeds the instabilities in the flow of sediment in equation (10) on a longer time-scale. Two distinct mechanisms lead to such noise. First, on convex surfaces the water flow equation (13) has solutions that develop shocks. Second, the sediment flow equation (10) develops solutions that are not smooth and that possess singularities in a derivative of H. These singularities, which we may interpret as waterfalls or rapids, produce shocks in the flow of water. While both of those mechanisms are present on the convex surfaces that emerge and exist during the early stages of the evolution of the land surface, the second shock producing mechanism dominates the concave surfaces characterizing the mature stage of evolution.

In addition to these sources of noise, small round-off errors in the numerical computations provide a continuous source of small random inputs to both the flow of water over the surface and the erosion of the surface. Since the numerical procedures that underly the current analysis (described in Section 3.1) are more accurate than the procedures used in [Smith et al., 1997b], they do not require the suppression of machine noise with the use of large amounts of artificial viscosity in order to obtain convergent computations. Instead there are very small amounts of numerical viscosity present in these more accurate methods, enough to suppress round-off errors but small enough to give an accurate presentation of the evolution of the small scales.

In summary, we have identified two mechanisms capable of generating noise in the system, on the intermediate time-scale defined in Section 2.1, in addition to the initial random perturbations to the surface and the extremely small numerical noise. While these randomly driven disturbances occur at the short time-scales defined by τ = t/η^{2} and τ = t/η, their cumulative effect is felt in the sediment flow on the long time-scale of erosion through the dependence on h and the term −η\frac{∂h}{∂t}. They act as seeds for the instabilites in the sediment flow and are amplified into noise in the sediment flow on the long time-scale in Section 2.1. These effects constitute the random variations in the model (9)-(10) that drive the dynamics of the system.

2.5 Complex Surfaces, Structure Functions, and Scaling Relations

Surfaces and interfaces modeled by nonlinear PDEs driven by noise are, in general, too complex to allow an analytic representation. In many cases the solutions that exist do not even possess derivatives. Hence various methods have been developed to represent the structure, and particularly the statistical structure, of such surfaces. An important set of characterizations involve structure functions that take the form

ρ_{k}(x,t) =< |F (x + y ,t) − F (y,t)|^{k} >  (14)

in which F (x,t) is some function characterizing a surface, x = (x,y) are lag variables, <> represents the expectation operator taken over some ensemble of (randomly selected) surfaces, and k = 1,2,... is a positive integer. If these structure functions possess a scaling, then we can conclude that the system is statistically self-similar and has the same structure and perhaps dynamics on all length scales including arbitrarily small ones.

Certain cases of the structure function are especially important in applications, including the case k = 1, in which ρ_{1}(x,t) represents the mean value of the differences of the function at different distances apart, and the case k = 2, which is often employed in terms of its square root

W_{F}(x,t) ≡ ρ_{2}(x,t)^{\frac{1}{2}} = < |F (x + y ,t) − F (y,t)|^{2} > ^{\frac{1}{2}}  (15)

to represent the root mean square of the elevation differences as a function of different distances of separation | x |. This function, which is known variously as the variogram, heightheight correlation function, roughness function, or width function [Edwards and Wilkinson, 1982, Kardar et al., 1986, Buldyrev et al., 1992], is often used to characterize the roughness of the surface in the direction of growth. We may also interpret W_{F}(x,t) as measuring the correlation between heights at points separated by a distance | x |. We note that the width function and the spatial autocorrelation function

C_{F}(x,t) =< F (y + x ,t)F (y,t) >  (16)

are related by

C_{F}(x,t) = σ_{F}^{2}(t) − \frac{1}{2}W_{F}^{2}(x,t),  (17)

in which σ_{F} is the variance ^{9}. An increase in W corresponds to an increase in correlations between two points on a surface. However, the scalings of C and W^{2} are, in general, different. Brownian motion, for example, is characterized by σ_{F}(t) ∼ t^{1/2} but by W_{F} ∼ t^{1/4}.

9In practice the variance is always finite for topographic surfaces.

Important questions concern the manner in which structure functions and their analogs vary with changes in the values of the lag variable x and time, and in particular, whether their variation may be represented in terms of simple polynomial functions. In many applications, the width function possesses particularly simple asymptotic forms representing scaling invariances in the system. It is often the case that the behavior of a surface is well described by asymptotic forms of (15) that are applicable in different space-time regimes of the system’s dynamics (see, for example, [Buldyrev et al., 1992, Paczuski, 1995, Sneppen, 1998]). Such forms have the representation

W_{F}(x,t) ∼ t^{β}f (t/|x|^{z}),  (18)

in which t is time, x is the lag variable, β is known as the temporal roughness exponent, and z is called the dynamic exponent. An interpretation of these scaling relations is that there is an equivalence between the temporal and the spatial scales given by t ∼ |x|^{z} that allows one to define various space-time regimes of the process.

The roughening of the surface, such as occurs for example when channelized forms emerge on fluvially eroding surfaces, is a transient phenomenon. In the space-time regime of transient behavior the roughness of a surface, as measured by W_{F} ∼ t^{β}f(t/x^{z}), typically increases. This regime is characterized by f(y) ∼ y^{−β} for y << 1, from which it follows that the width function may be approximated by a simple power law

W_{F}(x,t) ∼ t^{β}  (19)

in which β is the temporal roughness exponent. Since the temporal and spatial roughness exponents are related by^{10}

χ = zβ,  (20)

the width function in this case may be approximated by the relationship

W_{F}(x,t) ∼ | x |^{χ}  (21)

in which χ is the spatial roughness exponent. It is this growth in spatial roughness that is easiest to compute and is estimated numerically below. The length of the transient time interval is of the order t ∼ L^{z}, where L is the size of the system^{11}. Following this transient period, there may be a space-time regime of stationary, or saturated, behavior in which the roughness of the surface fluctuates about statistically steady values. In this regime, which is characterized by f(y) ∼ constant, y >> 1, pairs of points on the surface that are separated by a distance |x| ∼ L are characterized by differences in the expected values of the width function of order of |x|^{χ}.

10This indicates that the dynamic exponent z may be interpreted as the constant of proportionality linking the temporal and spatial scales.

11This indicates that the system takes time t = L^{z} to completely reorganize itself.

In surface growth models in which W_{F}(x,t) measures the roughness of the surface in the direction of the growth, χ may be interpreted as the fractional part of the dimension of the surface, i.e., the embedding or avalanche dimension is d + χ, where d is an integer. Roughly speaking the dimension d characterizes the spatial extent of the avalanche and χ characterizes the height of its “bulge”, see Sneppen [Sneppen, 1998].

If the roughness exponent χ does not depend linearly on k then the scaling is called multifractal, while if the scaling exponent depends on the direction of the lag vector, the scaling is said to be anisotropic (see, for example, [Klinkenberg and Goodchild, 1992], [Lavallée et al., 1993]). For an arbitrary structure function of order k, for example, this may take the form ρ_{k}((x,y)) = constant × (|x|^{kχ_{x}}+ |y|^{kχ_{y}}) in which χ_{x} ≠ χ_{y}. In the case of single-fractal and isotropic scalings, only two scaling exponents are required to determine all others (see, for example, [Sneppen, 1998]). The fractal dimension, for example, may be represented as D_{f} = n − χ, where n is the total (spatial) dimension, while the embedding or avalanche dimension may be represented by D = d +χ, where d is the “dimension” of the interface^{12}. The avalanche dimension relates the spatial extent of an avalanche to its temporal growth, |x| ∼ t^{D} ^{13}.

An important class of statistical phenomena may be characterized in terms of processes that take the form of fractional Brownian motion [Mandelbrot and Ness, 10]. The variance of such processes, which are of relevance for the analysis presented in Section 5.3, satisfies scaling laws of the form

σ ∼ t^{H}f(t/|x|^{z})  (22)

in which 0 < H < 1 is known as the Hurst exponent. This exponent is the leading exponent in the rate of increase of the correlation function and is, in general, different from the temporal roughness exponent β that gives the temporal rate of increase of the structure functions.

The usual Brownian motion is a special case of fractional Brownian motion with Hurst exponent 1/2 and with time increments that are uncorrelated in time. For fractional Brownian motion with Hurst exponent H > 1/2, increasing trends in the past lead to a further increases in the future while for fractional Brownian motion with Hurst exponent H < 1/2, decreasing trends in the past lead to further decreases in the future for arbitrarily large times. It is important to note that systems characterized by regular Brownian motion achieve a statistically stationary state with respect to their time-increments in the long run, whereas for fractional Brownian motion H ≠ 1/2, the long-time state is not stationary. However, in both cases the long-time behavior is different from the transient behavior, which explains the difference between the Hurst coefficient and the temporal roughness exponent β. It follows that for Brownian motion (H = 1/2) the (generalized) time derivative is statistically stationary white noise.

121 for a fractal line, 2 for a fractal surface, etc..

13We may define an avalanche as a burst of activity in surface growth, in which a significant number of sites move together.

2.6 Estimating the Scaling Exponents from Experimental Data

In the analysis presented below, we employ width functions to characterize the correlation structure of: (1) the water depth h in terms of the width function W_{h}; (2) the water surface elevation H in terms of the width function W_{H}; and (3) the water surface slope | ∇H | in terms of the width function W_{|∇H|}. Because of the relation (20), we only require two exponents in characterizing scaling relations involving these width functions and we employ the temporal and spatial roughness exponents β and χ. The spatial roughness exponent is computed from the scaling of the width function (15). We also make use of the fact that we may represent various scaling relations in terms of scaling of the (temporal) power spectrum, see [Maslov et al., 1994]. For present purposes, it suffices to note that

Ĉ = S(k) ∼ \frac{1}{k^{1+2β}}  (23)

in which S(k) represents the power spectrum at frequency k. Employing the power spectrum, we may compute the temporal roughness exponent β independently of the spatial roughness coefficient χ, and thus also obtain the dynamic coefficient z = χ/β using equation (20).

In validating and interpreting our results, we make use of known results concerning the scaling of solutions to equations in the same universality classes as those we examine below.

3 The Two Scalings: Numerical Results

Numerical simulations of the evolution of eroding surfaces generated by equations (9), (10) provide convincing evidence for the emergence of two fundamental landscape scalings. We briefly discuss the numerical methods underlying our analysis and then describe the scaling results.

3.1 Numerical Methods for Accurate Scaling Relations

While the numerical methods employed by Smith, Merchant, and Birnir [Smith et al., 1997b] in solving equations (6), (7) lead to solutions that represent the correct, separable landscapes, they turn out to be insufficient for deriving numerically accurate scaling results.

These solution methods, described in detail in [Smith et al., 1997b], involve alternative computations of equilibrium water flows and very small erosional changes to the land surface, using an explicit, two-step MacCormack predictor-corrector finite-difference scheme adapted for two-dimensional overland flow [Zhang and Cundy, 1989, Smith and Merchant, 1995]. This scheme involves the use of artificial viscosity to avoid numerical oscillations [Fennema and Chaudhry, 1990]. By suppressing the noise that is amplified by the instability of the initial surface, however the small scale dynamics become inaccurate. Hence these methods are inadequate for numerically accurate estimates of the scaling parameters.

For the purposes of the current investigation, we designed a significantly more accurate numerical scheme capable of representing the correct small scale dynamics and based on the existence theory of the nonlinear PDEs. By solving the hyperbolic equation (9) in two dimensions using an upwind scheme, we render the scheme stable and eliminate the need for artificial viscosity.^{14} A Crank-Nicholson scheme may then be applied in computing the nonlinear gradient flow (10) and by solving the resulting implicit matrix equation. This computation is significantly less expensive than the previous computations because in the explicit scheme significant amounts of computation are required in order to dampen instabilities in the water depth by artificial viscosity^{15}.

We note in passing that the comparison with the previous numerical scheme is an interesting study in numerical analysis and we refer the reader to [Birnir et al., 1999] for further details on the new numerical scheme. This code will also be made public.

3.2 The Numerical Simulations

Employing the new code, we ran a full suite of computational experiments simulating the erosion of a section of a linear ridge satisfying the initial and boundary conditions (12), (11), which are the same as those employed in [Smith et al., 1997a, Smith et al., 1997b]. These experiments were run using various refinements of the grids and over all relevant parameter values of the model (9, 10) to provide a basis for deriving the scaling relationships that we describe in Section 3.3.

In Figure 1-4, we show typical configurations of eroding surfaces, water depth, and slopes of the water surface for values of the sediment transport parameters γ = 2,δ = 2 at two characteristic times: when 10% of the surface is eroded and when 60% of the surface is eroded, which is well into the separable regime. In these figures, which indicate that the

14An upwind scheme still produces small amounts of numerical viscosity but of magnitude less than the artificial viscosity in [Fennema and Chaudhry, 1990].

15We implemented the scheme on a parallel supercomputer at the University of California, Santa Barbara, employing as a final step the implementation of a four color scheme to parallelize the procedure for maximum computational advantage.

Figure 5: Log-log plot of water depth at a given cross-section

x − axis runs laterally across the ridge and the y − axis runs longitudinally down the ridge, we represent the water surface H = z + h, in Figures 1 and 3; we represent the water depth h and the gradient of the water surface | ∇H |, in Figures 2 and 4, again at 10% and 60% of the surface eroded.

3.3 The Basic Scaling Results

The scaling results that we discuss in this section are based on numerical solutions obtained from experiments with sediment transport parameters γ = 2,δ = 2. We note that there is no significant variation of our results in a whole neighborhood of such values of γ and δ and that experiments starting with different initial conditions result in the same scaling results.

Scaling results were obtained by computing the width function (15) over sections of the ridge either in the x-direction transverse to the ridge or in the y-direction parallel to the longitudinal axis of the ridge. Estimates of the exponent χ of the asymptotic forms (21) were provided by the slope of log-log plots of width functions averaged over an ensemble of results from randomly seeded runs. Numerical estimates of χ were obtained by least squares fits. In Figures 5 and 6, we illustrate typical log-log plots taken at a given transverse section of the ridge for the water depth h and water surface H respectively. The R^{2}-value associated with each plot provides a measure of the appropriateness of scaling relation.

It is convenient to view the values of the scaling parameters as a function of their location

Figure 6: Log-log plot of the water surface at a given cross-section

on the ridge. For example, the values characterizing scalings in directions transverse to the ridge may be viewed as a function of the location on the longitudinal axis of the ridge χ_{x}= χ_{x}(y). We may then represent the scaling relations found over sections of the ridge by plotting both the estimated value of χ and the associated R^{2}-value as a function of location on the ridge. This is illustrated in Figures 7-10, in which we represent the scalings for h, H, and | ∇H |, at different cross-sections of the ridge at the two erosion times.

The scaling results from solutions to (9) and (10) for a time of 10% eroded are shown in Figures 7 and 8 as a function of distance upslope. These figures indicate that, initially, the structure functions of the water depth (W_{h}) and the gradient of the water surface (W_{|∇H|}) scale with the roughness exponent χ_{x} = 0.5 in the lateral (or x) direction. We note that this value characterizes a large range of y-locations on the ridge. It should be noted that, while in Figure 8 the numerically computed value of the exponent of |∇H| is χ = 0.25, the true value of this exponent is 0.5. The apparent discrepancy is caused by the initial slope H_{y}= c whose cross-term appears when one squares the components of |∇H| and then takes the square root, as explained in Appendix C.

As the solutions converge towards the separable forms defining the mature landscape of equation (8), the above scaling disappears to be replaced by another scaling. In a neighborhood of the separable landscape, at 60% of the surface eroded, the water surface scales in the x-direction as shown in Figure 9 with the scaling exponent χ_{x} = 0.7, again over a large range of y-locations on the ridge. While the corresponding numerically computed scaling of |∇H| is χ_{x} = 0.35, we note again that the true value of this exponent is 0.7 rather than the smaller value, for the reasons stated above and explained in Appendix C, the role of the

Figure 7: The scaling for the water depth

constant slope c now being played by the gradient of the separable water surface ∇H_{2}.

We note, therefore, that at 10% of the landsurface eroded the water depth h and the gradient of the water surface |∇H| scale, whereas at 60% of the landsurface eroded the water surface H itself and the gradient of the water surface |∇H| show the same scaling. In both cases |∇H| controls the scaling; in the former case it controls h; in the latter case it controls H. Our results indicate that H does not scale over any significant range of y-values in the first case, due to the presence of the initial slope, whereas h does not scale in the latter case. This clearly indicates that two distinct processes lead to the two different scalings.

In characterizing some of the details of these scalings, we may interpret the longitudinal axis of the ridge as a time scale and may think of the figures as representing increasing time values. Figures 7-10 indicate that the scaling is uniform neither in y nor in time. Initially, or at 10% of the landsurface eroded, the roughness exponent of W_{h} is approximately constant (0.5) for most of the upslope direction but deviates at the top where the boundary exerts its influence, as shown in Figure 7. As time increases this plateau disappears and one obtains different exponents as functions of y. Similarly, the width function of the water surface W_{H} initially shows different values of χ_{x} as a function of y but as the solution approaches the separable surface, a plateau is formed where the roughness exponent is constant (0.7) for a large portion of the upslope directions, as shown in Figure 9. Thus the existence of a plateau in χ_{x}, as a function of y, implies that the ensemble average possesses a scaling. These figures show the slope (or scaling exponent) as a function of the distance y up the ridge.

Figure 8: The first scaling of the gradient of the water surface

The scaling is anisotropic in the sense that there is no plateau in χ_{y} as a function of x, but this is more due to the fact that the upslope direction is really a time-scale, rather than a real anisotropy. We examined the slope-free part of the gradient H_{x}, in Appendix C, in various directions but found little evidence of significant anisotropy. We also examined both higher and lower structure functions ρ_{k}, k ≠ 2, or moments, of h,H and ∇H but found no evidence that the first scaling is multifractal. The second scaling on the other hand seems to be multifractal but this will be explored in a later publication. The spatial scalings were checked by computing the spatial power spectrum. The temporal power spectrum was very noisy due to the unevenness of the temporal evolution in the upslope direction and machine noise. It only gave reasonable results for the first scaling and then only for the water depth h.

4 The Emergence of Channels and the First Fundamental Scaling of Fluvial Landscapes

We now consider the emergence, early in the evolution of surfaces determined by the model (9)-(10), of scalings in the water depth h and in the slope of the water surface | ∇H | that are associated with the emergence of channelized flows. In particular, we focus on explaining how the fundamental scalings of 0.5 in the width functions W_{h} and W_{|∇H|} are driven by noise generated by divergences in the sediment flow that are seeded by instabilities in the water flow. Since this process is intimately related to the emergence of channels, we view

Figure 9: The scaling of the water surface

the nature and origin of this fundamental diffusive scaling as providing significant insight into the manner in which the early growth of channel-like features is dominated by random influences. In particular, these considerations provide a model of how random influences determine the location of channels and, as such, may be used in justifying the assumptions that underly the applications of various discrete stochastic modeling approaches to fluvial networks.

4.1 The Emergence of Channels and the Fundamental Channelization Instability

Results from previous analyses of models of the form (1), (2), provide insight into the earliest stages of channel evolution. Smith and Bretherton [Smith and Bretherton, 1972] analyzed versions of equations (6), (7) in which the constitutive relation for sediment transport was generalized to q_{s} = F (q_{w},| ∇H |), with \frac{∂F}{∂q_{w}} > 0 and \frac{∂F}{∂S} > 0 ^{16}. Based on a linearized analysis of these equations, they showed that if \frac{F}{q_{w}} < \frac{∂F}{∂q_{w}}, then small amplitude disturbances of small lateral scale will initially grow very rapidly, with the initial growth rate being greatest for disturbances with the smallest wavelengths. This criterion is equivalent to the condition γ > 1 for the particular sediment transport law embedded in equation (5). This latter fact leads to the particularly simple criterion that channel incision takes place when the sediment transporting capacity of water increases at a faster rate than the

16The term \frac{∂h}{∂t} in equation (7) was not considered in these analyses.

Figure 10: The second scaling of the gradient of the water surface

discharge of water, and occurs faster for shorter wavelength disturbances because the convergence of water flows is greater for disturbances with steeper lateral gradients. These results were reconfirmed by Loewenherz [Loewenherz-Lawrence, 1991], who also derived the full asymptotic solution in the small wavelength limit.

The simulations of Smith, Merchant, and Birnir [Smith et al., 1997b] confirmed numerically that the highest initial growth rates occur at the highest frequencies. They also found, however, that the growth rates of lower frequencies became dominant after a very short initial period, presumably as a result of saturation of the nonlinearities in the erosion process. The significance of these results for the present analysis is that they show that the instabilities of the initial erosion process are highly susceptible to being driven by random influences.

We now provide an alternative analysis of these early channel forming instabilities of (1),

(2) which provides additional insight into the underlying processes. In particular, we linearize the equations (9) and (10) about an initial surface, linear in y and uniform in the x direction, and a constant water depth H = cy + d,h = d ^{17}, where c and d are constants. We note that, while these initial conditions differ from those in the stability analyses of [Smith and Bretherton, 1972, Loewenherz-Lawrence, 1991], they lead to essentially the same results but can be regarded as one step earlier in the evolution from a ridge with constant slope. As such they provide an additional interpretation of the channel forming instabilities.

17Strictly speaking we should let the water depth go to zero at the top of the ridge to satisfy the boundary conditions there. This, however, makes no difference to the results of the analysis.

We let

h = d + ϵa(x,y,t), H = cy + d + ϵb(x,y,t),

and include a (very) small stochastic variation in the rainfall rate

R = R_{o}+ ϵr(x,y,t).

The linearized system of equations becomes

η^{2} \frac{∂a}{∂t} = \frac{5}{3}d^{\frac{2}{3}}c^{\frac{1}{2}} \frac{∂a}{∂y} + \frac{d^{\frac{5}{3}}}{c^{\frac{1}{2}}}(\frac{∂^{2}b}{∂x^{2}}+ \frac{1}{2} \frac{∂^{2}b}{∂y^{2}}) + r,  (24)

\frac{∂b}{∂t} = \frac{5}{3}d^{\frac{2}{3}}c^{\frac{1}{2}}(\frac{1}{η} + 2c^{\frac{5}{2}}d^{\frac{2}{3}})\frac{∂a}{∂y} + \frac{d^{\frac{5}{3}}}{c^{\frac{1}{2}}}(\frac{1}{η} + c^{\frac{5}{2}}d^{\frac{5}{3}})\frac{∂^{2}b}{∂x^{2}}

+ \frac{d^{\frac{5}{3}}}{c^{\frac{1}{2}}}(\frac{1}{2η} + \frac{5}{2}c^{\frac{5}{2}}d^{\frac{5}{3}})\frac{∂^{2}b}{∂y^{2}} + \frac{r}{η} .  (25)

The numerical results reported in Section 3.3 indicate that the variation of the water surface b does not scale, initially at least, with a roughness exponent that is constant over a range of y values. However, the gradient ∇b does scale from which it follows that the equations generating the scaling must be equation (24) and the equation^{18}

\frac{∂∇b}{∂t} = \frac{5}{3}d^{\frac{2}{3}}c^{\frac{1}{2}}(\frac{1}{η} + 2c^{\frac{5}{2}}d^{\frac{2}{3}})\frac{∂∇a}{∂y} + \frac{d^{\frac{5}{3}}}{c^{\frac{1}{2}}}(\frac{1}{η} + c^{\frac{5}{2}}d^{\frac{5}{3}})\frac{∂^{2}∇b}{∂x^{2}}

+ \frac{d^{\frac{5}{3}}}{c^{\frac{1}{2}}}(\frac{1}{2η} + \frac{5}{2}c^{\frac{5}{2}}d^{\frac{5}{3}})\frac{∂^{2}∇b}{∂y^{2}} + \frac{∇r}{η} .  (26)

Since η = d/100c ∼ h/H is small, the term ∇r /η can be large if r is not exceedingly small. However in most cases r can be taken to be very small and therefore we have dropped it in the equation below, still keeping in mind that very small fluctuations are always present in landscape evolution.

A straight-forward stability analysis shows that the initial unchanneled surface is unstable on the long time-scale of the sediment flow. We set \frac{∂h}{∂t} = 0 in the equations (9) and (10), or equivalently use the slow equations (6) and (7), with \frac{∂h}{∂t} = 0. On taking the gradient, equation (24) becomes

18We have merely applied the ∇ operator to equation 25.

\frac{5}{3}d^{\frac{2}{3}}c^{\frac{1}{2}} \frac{∂∇a}{∂y} = − \frac{d^{\frac{5}{3}}}{c^{\frac{1}{2}}}(\frac{∂^{2}∇b}{∂x^{2}} + \frac{1}{2} \frac{∂^{2}∇b}{∂y^{2}} )  (27)

and a substitution into (26) gives the equation

\frac{∂∇b}{∂t} = \frac{5}{3}d^{\frac{2}{3}}c^{\frac{1}{2}}[(−2 + d)\frac{∂^{2}∇b}{∂x^{2}} + (−\frac{1}{2} + 3d)\frac{∂^{2}∇b}{∂y^{2}} ],  (28)

which is a negative heat equation if d < 1/6. Recall that d is the initial water depth that can be taken to be very small. The dispersion relation for (28)

ω = \frac{5}{3}d^{\frac{2}{3}}c^{\frac{1}{2}}[(2 − d)k_{1}^{2}+ (\frac{1}{2} − 3d)k_{2}^{2}],

shows that all the spatial frequencies are unstable and that the highest frequencies grow the fastest. We note that this is essentially the same instability as discovered by [Smith and Bretherton, 1972, Loewenherz-Lawrence, 1991], but expressed in terms of the unstable growth of arbitrarily small (spatial) frequencies in the slope of the surface rather than in the elevation of the surface. In particular, this result is nicely in accordance with an interpretation of the stability results of Smith and Bretherton [Smith and Bretherton, 1972] that suggests that smaller perturbations grow at faster rates because, for a given amplitude of perturbation, they possess steeper slopes which lead to more rapid convergence of eroding water flows.

The significance of this result follows from our preceding conclusion that fluctuations of arbitrarily small frequencies are always present, and that these instabilities always lead to the exponential amplification of noise seeds, from the water flow, and the consequent channelization of the initial surface.

4.2 The Generation of Noise During Initial Channelization

We now provide and justify a hypothesis that instabilities in the water surface provide a significant source of white noise that seeds the instabilities of the erosion process. In particular, we show that shocks in the form of hydrolic jumps form in the flow of water over initially convex portions of the ridges between emerging channels. Although a concavity in the longitudinal profile emerges early near the base of the ridge, moves longitudinally towards the top of the ridge, and then expands laterally across the emerging valleys (for an analysis of this process, see [Smith et al., 1997a]), the lateral profiles of the ridges between emerging channels take an initial form that is typically convex in their upper parts. Two upridge sequences of lateral profiles that indicate, respectively, the existence of initially convex sections of the interfluvial ridges at 10% of the surface eroded and the dominance of fully concave lateral profiles at 60% of the surface eroded are shown in Figure 11.

We base our analysis of the mechanism that provides a significant source of white noise on these convex sections of ridge by linearizing equations (9) and (10) around a convex (interfluvial) profile. In particular, we assume that our perturbed solutions take the form

h(x,y,t) = h_{1}(x,y,t) + v(x,y,t), H(x,y,t) = H_{1}(x,y,t) + u(x,y,t)

where H_{1}(x,y,t) represents a convex portion of an interfluvial ridge, h_{1}(x,y,t) the depth of water flow over this portion, and v(x,y,t),u(x,y,t) are the small perturbations to these quantities. The linearized equations are

η^{2} \frac{∂v}{∂t} = ∇ · [ \frac{5}{3}h1^{\frac{2}{3}} \frac{∇H_{1}}{|∇H_{1}|^{1/2}} v] + ∇ · [h1^{\frac{5}{3}} \frac{∇u}{|∇H_{1}|^{1/2}} − \frac{1}{2}h1^{\frac{5}{3}}(∇H_{1}· ∇u)\frac{∇H_{1}}{|∇H_{1}|^{3/2}} ] , (29)

\frac{∂u}{∂t} = ∇ · [ \frac{h1^{\frac{5}{3}}}{|∇H_{1}|^{1/2}} (\frac{1}{η} + h1^{\frac{5}{3}}|∇H_{1}|^{5/2})∇u] + ∇ · [ \frac{5}{3}h1^{\frac{2}{3}} \frac{∇H_{1}}{|∇H_{1}|^{1/2}} (\frac{1}{η} + 2h1^{\frac{2}{3}}|∇H_{1}|^{5/2})v]

+ ∇ · [h1^{\frac{5}{3}}(\frac{−1}{2η} + 2h1^{\frac{5}{3}}|∇H_{1}|^{5/2})(∇H_{1}· ∇u)\frac{∇H_{1}}{|∇H_{1}|^{3/2}} ] ,  (30)

The first equation is a hyperbolic equation for v driven by ∇u. The second equation is a parabolic equation for u driven by v and ∇v.

Because the flow of water occurs down the gradient of the surface H_{1}(x,y,t), we may view the first equation as a hyperbolic PDE in one space dimension. Namely, if we let u denote the unit vector in the direction of the gradient of the water surface, we can write the first equation in the form

η^{2} \frac{∂v}{∂t} = \frac{5}{3}h1^{\frac{2}{3}}|∇H_{1}|^{1/2}u · ∇v + g(x,y,t)v + f(x,y,t).  (31)

If we then let the scalar s parametrize the direction of the gradient, we may rewrite this equation in the form

η^{2} \frac{∂v}{∂t} = \frac{5}{3}h1^{\frac{2}{3}}|∇H_{1}|^{1/2} \frac{∂v}{∂s} + g(x,y,t)v + f(x,y,t).  (32)

Figure 11: Transverse Sections at 10m separation on the Eroding Ridges at 10% and 60% Eroded The last equation is analyzed in Appendicies A and B, where it is shown that it develops shocks if the profiles about which we linearize are convex or have knick-point singularities.

After operating with the gradient, we may lump the terms in (30) causing the instability and involving the water depth together into a single noise term r. This term is the noise seeded by the shocks in (29) and generated by the sediment divergences in (30) that amplify the instabilities produced by (29). These instabilities do not lead to exponential growth because the water changes on a faster scale than the sediment, see Section 2, but instead give rise to white noise in the sediment flow. We obtain, as a result, a heat equation driven by white noise,

\frac{∂(∇u)}{∂t} = νΔ(∇u) + r(x,y,t),  (33)

where ν = \frac{d^{\frac{5}{3}}}{ηc^{\frac{1}{2}}} . We want to emphasize here that we have provided an intuitive explaination

of the noise-making generation in the sediment flow but no detailed mathematical model for it. The transition from equations (29) and (30) to (33) or the modeling of the amplification of the fast-time noise in the water flow by the divergences in the sediment flow to produce white noise, must be taken to be a hypothesis. We hope to be able to provide a detailed model of the noise-generating mechanism in future publications. However, as will be shown below with this hypothesis we can explain the numerically observed scalings.

The solution to last equation (33) is well-known to be Brownian motion

∇u = ∫_{0}^{t}∫_{IR_{2}} \frac{exp(−|x − y|^{2}/4ν(t − τ))}{(4π(t − τ))} dW(y,τ),  (34)

where x = (x,y), y = (x′,y′) and dy = dx′dy′is the volume element in IR^{2} and formally the white noise can be understood as dW(y,τ) = r(y,τ)dydτ. Adding this solution to ∇H_{1} we obtain a random walk

∇H = ∇H_{1}+ ϵ∇u,

of the slopes of the water surface down the channelizing surface (H_{1}), driven by white noise. The water depth follows suit, by equation (9).

4.3 The Scaling of the Channelization Process

We are now in a position to interpret the scaling exponents that emerge during the initial stages of the erosion process by using well-known results from the theory of the evolution of complex interfaces. We note that equation (33) is the Edward-Wilkinson equation of interface growth theory, which is known to model surfaces over which the deposition rate equals the removal (or evaporation) rate (see [Sneppen, 1998]). The values of the exponents of the associated width functions are

β = \frac{1}{4}, χ = \frac{1}{2}, z = 2.

We note that β is different from the Hurst exponent, which has the value 1/2 for a pure random walk, because β measures the transient towards a pure random walk. The resulting motion may be interpreted as a random walk occurring on top of another random walk. The gradients of the water surface walk randomly in the downslope direction and the nodes where ∇H(x + y) − ∇H(y) vanishes do another random walk on top of that. Hence the configuration of the emerging channel forms may be characterized in terms of random walks. We note that this important finding provides a physically-based justification for various stochastic modeling approaches to the emergence of network forms such, for example, as described in [Scheidegger, 1967a, Takayasy et al., 1991]. The surprising observation is that the stochastic and the deterministic approach lead to the same conclusions and can be considered to be equivalent methods to model the same phenomena; the time-evolution of a fluvial network.

This surface sits in a 3-dimensional space and the fractal dimension is D_{f} = 3 − χ = 5/2. The avalanche dimension is the same D = 2+χ = 5/2. The temporal roughness coefficient β may be obtained from the power spectrum. The reason is that the erosion is initially greatest towards the bottom of the slope and then the channelization migrates upstream and ends when the concavity, forced by the boundary conditions at the bottom of the slope, reaches the top of the slope, see [Smith et al., 1997a, Smith et al., 1997b]. This implies that the y direction is really a time-scale as discussed above. That is, instead of measuring the percentage of the landsurface that has been eroded, we can simply measure how far upstream the channelization process has migrated.

5 Self-Similar Landscapes and the Second Fundamental Scaling of Fluvial Landscapes

In further support of our hypothesis that the observed scalings arise from identifiable physical processes, we now show that the second scaling of 0.7 in Section 3.3 may be interpreted as characterizing the erosion processes in the mature, separable landscapes. These processes are driven by noise that is seeded by shocks in the water flow over the concave slopes of landscapes such as those shown in Figure 11, and as explained in Appendix B. These water instabilities are then fed into the divergences in the sediment flow and since the water instabilities continue to vary on a faster scale, the sediment divergences produce noise in the sediment flow. That is different water profiles with different noise are present whenever a significant amount of sediment is eroded. However, there is little water on top of the concave ridges in Figure 11, see Figures 2 and 4, and there the noise is quenched. On the lower reaches of the concave slopes and in the valleys the same noise generating mechanism as discussed in the previous section is at work, so the noise is not uniformely distributed over mature landscapes. This results in a colored noise and scaling which is very different from the previous scaling. In particular, we provide evidence indicating that the resulting behavior is characteristic of a SOC system describing the driven motion of an interface in random media [Leschorn, 1993] and [Sneppen, 1998].

5.1 The Emergence of Self-Similar Landscapes

Smith, Birnir, and Merchant [Smith et al., 1997a] showed that the model (6)-(7) possesses separable solutions H(x,y,t) = H_{0}(x,y)T_{H}(t) and q_{w}(x,y,t) = q_{0}(x,y). These are associated with an optimality principle by which mature eroding surfaces evolve to minimize a simple function of the sediment flux over the surface, subject to two constraints. In their analysis of numerical solutions of (6), (7), Smith, Merchant, and Birnir [Smith et al., 1997b] showed that after a characteristic erosion time^{19} the surfaces exhibit clear and convincing evidence of convergence towards such separable solutions, and the corresponding satisfaction of the optimality criterion.

As shown in Section 3, this stage of development is characterized by a dominant scaling of the width function with exponent 0.7. We now analyze the model (9), (10) in order to determine the origin and physical basis for this scaling.

5.2 The Origin of the Second Diffusive Scaling

The scaling of the water flow down the slope of a mature landscape and the resulting scaling of the water (and land) surface is different from the scaling of channel formation in Section 4. To find and analyze this scaling we linearize equations (9) and (10) about the separable solutions (8) representing the mature landscapes of valleys and ridges^{20}

H = H_{2}(x,y,t) + ϵu(x,y,t), h = h_{2}(x,y.t) + ϵv(x,y,t),

where H_{2} = H_{o}(x,y)T (t), h_{2} = h_{o}(x,y)T^{− \frac{3}{10}} (t), are the separable solutions of the equations (9, 10). The form of the equations that we obtain for u(x,y,t), v(x,y,t) by this linearization process is the same as that characterizing the early period of channel emergence discussed in Section 4.2, namely (29) and (30), but with the terms H_{1}(x,y,t), h_{1}(x,y,t) replaced by the separable solutions H_{2}, h_{2}.

19The characteristic time is 45% of the surface eroded for planar initial surfaces with γ = δ = 2.

20As in our previous analyses, we employ sediment transport parameter values γ = δ = 2, noting that similar results hold for parameter values in a neighborhood of these.

Again, the first equation is a hyperbolic equation for v driven by ∇u and the second is a parabolic equation for u driven by v and ∇v. Since the first equation is really a hyperbolic PDE in one space dimension, exactly as in Section 4.2, we let the scalar s parametrize the direction of the gradient to obtain

η^{2} \frac{∂v}{∂t} = \frac{5}{3}h2^{\frac{5}{3}}|∇H_{2}|^{1/2} \frac{∂v}{∂s} + g(x,y,t)v + f(u,t).  (35)

A straight-forward analysis of this equation, presented in Appendix B, shows that its solutions develop shocks, for separable surfaces with (slope) singularities, since concave slopes with knick points are the dominant feature of the mature separable landscape, as illustrated in Figure 11.

Thus the water flow develops shocks at the top as it flows down the separable ridges, which then propagate down the slopes. In addition on the lower reaches of the slopes and in the valleys the instabilities discussed in Section 4 are at work. The sediment flow is driven by the pull of gravity down the slopes and because there is no water on top of the concave ridges, see Figures 2 and 4, there the surface is pinned. Thus the devergences in the sediment flow magnify the seed in the water flow into a highly colored noise. We make this hypothesis and consequently turn the equation (29), with the separable solutions H_{2}(x,y,t), h_{2}(x,y,t) replacing H_{1}(x,y,t), h_{1}(x,y,t), into a diffusion equation driven by quenched noise or the Langevin equation

\frac{∂u}{∂t} = Δu + F + r(u,t),  (36)

where F is the driving force, r is the quenched noise, and λ is a parameter. Moreover, it is very different from equation (33) since the noise term

\frac{5}{3}h2^{\frac{2}{3}} \frac{∇H_{2}}{|∇H_{2}|^{1/2}} (\frac{1}{η} + 2h2^{\frac{2}{3}}|∇H_{2}|^{5/2}) · ∇v

is not spatially homogeneous. The water depth h_{2} becomes very small at various locations on the separable surface, and at these locations the noise is quenched. One can think of these spots as pinning the surface and creating small local minima in the sediment flow. The last term in equation (29) plays the role of the forcing. It is known, see [Sneppen, 1998], [Leschorn, 1993], that such processes produce polynomial scaling with the spatial roughness exponent close to 0.75. This implies that the water surface

H = H_{o}(x,y)T (t) + ϵu(x,y,t)

executes a random motion driven by noise produced by the divergence in the sediment flow, seeded by shocks ∇h_{o}(x,y)T^{− \frac{3}{10}} (t) + ϵ∇v(x,y,t) propagating down the separable surface on a faster time-scale, the surface being pinned at randomly distributed sites. The slope ∇H evolves in the same fashion.

5.3 The Scaling of Mature Landscapes

We compute the exponents using the formulas in Section 2.5, starting with the numerical value of the spatial and temporal roughness exponents from Section 3. The values of the exponents are

β = 0.4, χ = 0.7, z = 1.8.

The fractal dimension becomes D_{f} = 3−χ = 2.3. The avalanche dimension is D = 2+χ = 2.7. The values of these exponents are close to the corresponding values for a driven interface in random media in 1+2 dimensions, see [Leschorn, 1993] and [Sneppen, 1998], and described by equation (36),

β = 0.475, χ = 0.75, z = 1.58.

The computed spatial roughness coefficient χ is slightly smaller (0.7) than the spatial roughness coefficient (0.75) of the driven interface, however this may be due to the fact that the system crosses over from one type of a scaling (Edward-Wilkinson) to another (driven interface in random media) and if we would follow it for a longer time (than 60% eroded) the spatial roughness component would eventually make it up to the latter value.

It is, in fact, well-known that landscapes generated by discrete landscape models do not reproduce realistic landscape unless random pinning that models quenched noise is included in the model, see [Calderri et al., 1996]. The latter models also give rise to values of the scaling exponent β similar to ours. Thus two remarkable features of the continuous model (9 ,10) are apparent. First, it generates the noise from the water flowing down the separable surface (8) and shocking as it passes over surface singularities. Second, it produces the pinning of the surface at random sites where these shocks and their consequent noise in the sediment flow is absent.

We note with interest that driven interface in random media is known to be an example of an SOC system ^{21}. Hence we may interpret the landscapes described by our advective model as SOC systems, at least at stages in the landscape evolution where the presence of the separable landscape is felt.

21As noted in [Sneppen, 1998]:

A system is a SOC system if the scaling (of the structure) functions during the saturated stage

is different from the scaling during the transient stage and the activity during the saturated

stage occurs in bursts that are self-organized during the transient.

6 Discussion and Significance of the Results

The preceding results on the origins of scaling in solutions to the models (9)-(10) have significant implications both for the evolution of many of the empirical regularities observed in fluvial landscapes and for the manner in which one may model such phenomena. We now describe a few of the more important implications.

6.1 The Two Processes Forming the Landscape

A picture of two processes shaping the formation of general landscapes emerges from the previous sections, although for most landscapes important details such as weathering, vegetation and mass wasting are missing. We call the first process the Channelization Process. It takes place on unchanneled surfaces and is characterized by a Brownian motion of the slopes in the water (or land-) surface marching down to form a network of channels and ridges. This process characterizes the initial stage of the landscape evolution (youth) and is associated with the spatial roughness exponent χ = 0.5. The system then gradually makes a transition to a second process that we call the Maturation Process. This process emerges much later when a system of large valleys and concave ridges close to a landscape determined by separable solutions has emerged, characterizing a mature landscape. The maturation process evolves the landscape towards the slowly decaying but self-similar separable landscape. Thus it characterizes the maturation (middle-age) of the landscape and is associated with the spatial roughness exponent of χ = 0.75. ^{22}

One has to put those two processes together to understand the shape of river basins. Namely, the channelization process controls the shape of the river basin along the main river, whereas the maturation process controls the shape perpendicular to the main river, see Figure 12. The main river is formed by a biased random walk down a slope ( a drunk that stumbles in the downhill direction) and this makes river networks analogus to directed percolation networks (DP) although the scaling exponents are in general different. In the direction perpendicular to the main river one runs into the separable slopes for the mature landscape. Thus for the river basin which is the projection of the surface in three dimensions to the plane the (planar) avalance dimension becomes D = 1 + χ. The exponent one measures the size of basin l along the main river and the exponent χ measures the size of its bulge

22We believe this to be the true value of χ althought the transcient numerically determined value of 0.7 is somewhat lower at 60% eroded.

Figure 12: The basin of the Amazon river

l^{χ} (width), in the perpendicular direction, see Figure 12, where l is the length of the main river.

6.2 The Origin of Observed Landscape Width Scalings

We note with great interest that empirical results for the height-height correlation function from Ethiopia, Somalia, and Saudi Arabia, as well as values generated for synthetic landscapes as reported by Weissel and Pratson [Weissel and Pratson, 1994], fall in the range 0.5-0.7. The Saudi Arabian data is virtually isotropic with the value 0.5, whereas some of the Somalian data is anisotropic, with values roughly 0.5 in one direction and 0.7 in the other. Although empirical scaling results of this form tend to cluster around a relatively small set of values, these values provide some validation for our theory particularly since our models may be interpreted as representing advective, transport-limited (or badland) conditions, which presumably hold in these areas.

The theory presented in this paper suggests that the spatial roughness scaling exponents for landscapes that are approximately described by our advection models should lie in a neighborhood of 0.5 and 0.7, the latter exponent growing to the value 0.75 for older landscapes. Young channelizing surfaces should exhibit isotropic scaling with exponent close to 0.5, whereas more mature surfaces should be anisotropic. This anisotropic scaling should have the scaling exponent close to 0.7 on the slopes of the separable ridges but on channelizing surfaces along the main rivers the scaling exponent should be close to 0.5.

6.3 The Origin of Hack’s Law:

The preceding results allow us to derive some of the fundamental scaling results that are known to characterize fluvial landsurfaces. In particular, the avalanche dimension derived in Section 5.3 allows us to derive Hack’s Law relating the length of a river l to the area A of the basin that it drains. This is the area of the river network that is given by the (planar) avalanche dimensions

A ∼ l^{D}

and the (planar) avalanche dimensions is D = 1 + χ. This relation says that if the length of the main river is l then the width of the basin in the direction, perpendicular to the main river, is l^{χ}. Stable scalings for the surface emerge together with the emergence of the separable solutions. We note that in this case χ = 0.75, hence we obtain

l ∼ A^{\frac{1}{1+χ}}  (37)
≈ A^{0.57}  (38)

a number that is in excellent agreement with observed values of the exponent of Hack’s law of 0.568, see [Gray, 1961]. We note, in particular, that χ is not the roughness coefficient of a fractal Brownian motion. We also note that the area A is the same as the area of voids (areas without branches) in a DP network (see [Sneppen, 1998]) and that Hack’s Law is just the relation between the area of the voids and the length of the DP backbone. Thus one can understand the "oval shape" of most river basins, see [Rodriguez-Iturbe and Rinaldo, 1997], as caused by the (second) scaling of the mature separable landscape, that rules the surfaces in the upper reaches of the river valleys and in the direction perpendicular to the main rivers.

6.4 The Scaling of Exceedence Probabilities

The probability P[A ≥ a] that the area A drained by a channel exceeds a is known as the exceedence probability. Empirical values for the exponent of this probability when it is approximated by the power law

P[A ≥ a] ≈ a^{ρ},  (39)

which is a good approximation to observed variation, are close to −0.4 (see, for example, [Rodriguez-Iturbe and Rinaldo, 1997]).

Using the scaling relation perpendicular to the meanderings of the main river and the two scaling above we can easily deduce this scaling of the exceedence probability once we know Hack’s law. The idea to keep in mind is that a river basin is similar to a pure DP network, in that when two branches of the river network approach each other, a small subbasin of the larger river basin is created. The subbasins without a river branch are the voids of the DP network, discussed above. The time at which two branches originating at the same location join, which is called the self-interaction time β_{si}, can be measured in terms of the length l of the main river branch. We let the dimension of the network along the main river branch (backbone of the DP) be one and the dimension in the cut perpendicular to the main river branch be χ as above. Then

P[A ≥ a] = P[l ≥ t]

since by Hack’s Law A = l^{1+χ} and a = t^{1+χ} , t being the length of the main branch of the subbasin having area a. However

∫_{0}^{t}P[l ≥ s]ds = ∫_{0}^{t}s^{−β_{si}} ds ∼ t^{1−β_{si}} ,

and

1 − β_{si} = −χ ,

see [Sneppen, 1998]. Thus using Hack’s law again we obtain

P[A ≥ a] ≈ t^{−χ}  (40)
≈ a^{− \frac{χ}{χ+1}}  (41)
≈ a^{−0.4} ,  (42)

for χ = 0.75. Notice that the exponent is the negative of the dynamic roughness exponent β = 0.4 in Section 5.3.

6.5 A Bridge Between Deterministic and Stochastic Modeling

Our analysis of the origin of the fundamental scaling with exponent 0.5 in the width function indicates that it arises as a result of white noise being generated by divergences in the sediment flow seeded by shocks in the water flow. In particular, we found that the evolution of the width functions of both the slopes and the water depths is well-characterized in terms of random walk processes. Hence the behavior and configuration of the emerging patterns of channels and valleys described by the nonlinear equations (9) and (10) and observed in the studies of Smith, Merchant, and Birnir [Smith et al., 1997b] is that resulting from a random walk.

We note that these results represent an important bridge between the deterministic and stochastic lines of modeling that we discussed in Section . From the stochastic modeling literature, it is well-known that random walk models of erosion in the plane lead to relatively realistic models of key aspects of river networks (see, for example, [Scheidegger, 1967b, Takayasy et al., 1991]) and are capable of generating branching networks that obey Horton’s laws (see [Rodriguez-Iturbe and Rinaldo, 1997] for a summary). We believe, in fact, that the mechanisms that we have shown to underly the fundamental 0.5 scaling provide the first demonstration that physically-based PDE models lead to tree-like drainage structures whose configuration is determined by a random walk process. Hence we may view our results as providing a physically-based justification for the assumptions underlying the various stochastic modeling approaches. We note, however, that further investigation is required in order to lay a fully-acceptable foundation for such assumptions.

6.6 A Model of the Stages of Drainage Basin Evolution

Based on the preceding results and the findings of our two previous papers [Smith et al., 1997a, Smith et al., 1997b], we now suggest an end-to-end theory of the evolution of landscapes that are described by our advective models. In particular, this theory characterizes landscape development in terms of an evolutionary sequence involving three stages and associated time scales ^{23}. In a first stage, occurring over a relatively short time-scale, the mechanisms underlying the first scaling of 0.5 become effective and a channelized network that may be interpreted in terms of random walks emerges. Hence in this stage, the foundation of random channel network configurations are laid down. This stage may be viewed as constituting the youthful stage of W.M.Davis’ grand model of landscape evolution [Davis, 1902].

A second stage, occurring at a longer time scale, involves a self-similar cascade towards decaying separable solutions. We may view the separable solutions as a transient attractor (see, for example, [Smith et al., 1997a]) towards which the system is drawn. The mechanisms underlying the second scaling become effective during this stage, leading to a scaling of the width function of 0.75 as well as to various scaling laws, such as Hack’s Law and the Law of Exceedence Probabilities. We may view this as constituting the second, or mature, stage of the Davisian model. A third and final stage involves a short time scale over which the separable solutions collapse as the potential energy in the system approaches small values and the stable global attractor H = 0 exerts its influence. This stage may be viewed as the final, or old age, stage in Davis’ model.

We realize that may elements are missing from this model, including such processes as tectonic and isostatic uplift, gravity-driven diffusion, weathering, subsurface water movement, vegetation effects, and a host of lesser processes. Nevertheless, we believe that our advective model captures the essence of processes that lead to the evolution of dendritic, fluvial structures and their most important scaling properties.

23Clearly, the three stages do not have well-defined beginnings and endings, but they do mark significant differences in the qualitative dynamics of the system.

7 Summary

The simple advective model of fluvial erosion (9)-(10) provides us with a compelling explanation for the fundamental processes present in landscape evolution. This model is capable of representing the process of channel initiation. Instabilities in unstructured initial surfaces lead to a channelization process which we are justified in viewing as the Edward-Wilkinson transient of Brownian motion. This is a noisy process: small, initial, random perturbations amplified by the nonlinear instabilities lead to noise. More importantly, however, the model (9)-(10) incorporates mechanisms capable of generating large noise. The relatively thin film of water flowing over the ridges between channels “folds over” to forms hydrolic jumps as it moves over initially convex portions of the interfluvial ridges. The derivatives of these shock waves constitute a white noise process that seeds the divergences of sediment in equation (10). The divergences in the sediment flow seeded by the white noise in the water flow create a noise-driven gradient flow for the sediment with a noise term that is white both in space and time. The solution of this equation, linearized about the initial surface, consists of Brownian motion. Thus the projections of the channels onto the plane execute a random walk and the locations on the surface where the slope equals the average slope do a random walk on top of the first one. This explains why river networks are observed to possess the characteristics of random networks and, in particular, are characterized by the relation (55) between the bifurcation and length ratios in Horton’s laws.

This random walk is, however, not the only process influencing the shape of river basins. As the channels merge into rivers a mature landscape consisting of a pattern of concave valleys and cusp-shaped ridges emerges. These surfaces, represented by the separable solutions of the evolution equations, are characterized by another stochastic process that differs from a random walk. The reason for this difference is that the mechanism generating the noise is now different from the mechanism generating noise in a channelizing surface. First, the separable surfaces are characterized by concavities whereas the ridges between the channels initially involve upper convexities. This means that shock waves in the water flow, originating from surface singularities, propagate down the surface. Most importantly, the vanishing of the water depth, for example on the ridge tops, on the separable surface, quenches the noise. These fast-flowing water instabilities are then fed into the sediment divergences evolving on a slower scale resulting in another noise-driven gradient flow for the sediment. However the noise is not white in space anymore and the sediment flow over the separable surfaces may be characterized as diffusion driven by quenched noise. The properties of this process, typical for a driven interface in a random media [Leschorn, 1993] and [Sneppen, 1998], are very different from those of Brownian motion.

Combining these two processes, the channelization and the maturation processes, we obtain a strikingly simple model of drainage basin evolution that is analogous to directed percolation networks. The channels execute a random walk of length l in the longitudinal direction under the channelization process. In a direction perpendicular to this random walk, however, the maturation process is dominant in regions in which separable solutions have emerged, giving us the width l^{χ} of the basin, where χ = 0.75 is the spatial roughness coefficient of the maturation process in Section 5.3. This accounts for the oval shape of river basins (see [Rodriguez-Iturbe and Rinaldo, 1997]), and leads to a derivation of Hack’s law (see Section 6.3) as well as the observed distribution of the exceedence probabilities, see Section 6.4.

These investigations lead to interesting sets of geomorphological and mathematical problems that we intend to address in future publications. In particular we intend to investigate whether the maturation process is multi-fractal, see [Lavallée et al., 1993]. While our models capture two fundamental processes, they also raise the question as the role and significance of processes of landscape evolution that are not as yet included in our models. In particular, it is of interest to conjecture whether different multifractal processes characterize more complex models. We also note that the analysis of the water flow in Section 2.4 and Appendix A provides complete information only for water flows of very small depth. If the rainfall rate, and hence the water depth, is significantly greater then more complex fluid equations are probably justified, and this may change the noise-generating mechanism and the stochastic driving of the erosion process. The two processes work in tandem in the general situation discussed in this paper but it is also desirable to find special situations where each process can be studied in isolation.

Acknowledgements

The authors wish to thank Kim Sneppen, Daniel Lavallée, A. Brad Murrey and A. Bertozzi for helpful comments. This work was supported by NSF grants DMS-9704874, BCS-9819095, DMS-9628606, and PHY-9601954, the last of which funded the use of UCSB’s Crunch supercomputer, as well as grants from the UCSB Academic Senate and the Icelandic Science Foundation.

A The Generation of Shocks in the Water Flow Equation

Equation (9) describes the water flow down a water surface represented by H. We rewrite this equation in the form

η^{2} \frac{∂h}{∂t} − ∇ · [uh^{\frac{5}{3}} |∇H |^{\frac{1}{2}}] = R,  (43)

where u = \frac{∇H}{|∇H|} denotes the unit vector in the direction of the gradient of the water surface. If we expand the gradient term, then this equation can be written in the form

\frac{∂h}{∂t} − \frac{5}{3η^{2}} h^{2/3}|∇H |^{\frac{1}{2}}u · ∇h = Q\frac{h^{5/3}}{η^{2}} + \frac{R}{η^{2}},  (44)

where

Q = [ \frac{|∇H |^{2}ΔH − \frac{1}{2}(H_{x}^{2}H_{xx}+ 2H_{x}H_{y}H_{xy}+ H_{y}^{2}H_{yy})}{|∇H |^{\frac{5}{2}}} ]

and using that u·∇h = \frac{∂h}{∂s}, where s parametrizes the downstream direction given by u, we obtain a a one-dimensional PDE

\frac{∂h}{∂t} − \frac{5}{3η^{2}} h^{2/3}|∇H |^{\frac{1}{2}} \frac{∂h}{∂s} = Q\frac{h^{5/3}}{η^{2}} + \frac{R}{η^{2}}.  (45)

where Q is the function defined above. It is well-known that such one-dimensional equations can either develop shocks or their solutions may blow up in finite time. We now show that for a small rainfall rate R the solutions do indeed develop a shock, which is a singularity in the s derivative, in finite time.

We first set R = 0 and write down the characteristic equations for the equation (45),

\frac{dt}{dτ} = 1, \frac{ds}{dτ} = −\frac{5}{3η^{2}} h^{2/3}|∇H |^{\frac{1}{2}}

(46)

\frac{dh}{dτ} = Q\frac{h^{5/3}}{η^{2}} = −p(s,τ)h^{5/3}.

The function p is positive

p(s,τ) = −Q > 0,

if the quadradic form Q(H_{xx},H_{xy},H_{yy}) is negative definite and it is as long as H is convex (concave in mathematical terms) and ∇H ≠ 0. Now we solve the h equation in (46) and obtain

h(s,t) = \frac{h_{o}(s)}{(1 + \frac{2}{3}h_{o}^{2/3}(s)P(s,t))^{3/2}}  (47)

where

P(s,t) = ∫_{0}^{t}p(s(τ) + s,τ)dτ,

increases along the characteristics because p > 0. Thus h(s,t) does not blow up but decreases along the characteristics. Next we compute the characteristics making an approximation that simplifies the computation. The equation determining the characteristics is

\frac{ds}{dτ} = − \frac{5}{3η^{2}} h^{2/3}|∇H |^{\frac{1}{2}} = − \frac{5}{3g(s)η^{2}} h^{2/3},  (48)

if H evolves on a slower time-scale than h. But this was precisely the conclusion of the discussion in Section 2.1. Thus we can separate the variables and solve the above equation

G(s)|_{s_{o}}^{s} = ∫_{s_{o}}^{s}g(s)ds = −\frac{1}{η^{2}} ∫_{0}^{t}h^{2/3}(s(τ),τ) dτ.

Next we expand the solution h in (48) above in powers of h(s_{o}),

h(s,t) = h(s_{o}) − h^{5/3}(s_{o})∫_{0}^{t}p(s(τ),τ) dτ + O(h^{7/3}(s_{o})),

and substitute the result into the equation for G, to get

G(s) − G(s_{o}) = −\frac{5}{3} \frac{h^{2/3}(s_{o})}{η^{2}}t + O(h^{5/3}(s_{o})/η^{2}).

This implies that two characteristics emanating from two distinct points s_{1} > s_{o} meet at the time

t = \frac{3}{5} \frac{η^{2}(G(s_{1}) − G(s_{o}))}{h^{2/3}(s_{1}) − h^{2/3}(s_{o})} + O(h(s_{o}),h(s_{1})).  (49)

Moreover, these characteristics carry different values of h by equation (47) producing a shock at the location

s = G^{−1} (G(s_{o}) − \frac{5}{3} \frac{h^{2/3}(s_{o})}{η^{2}}t + O(h^{5/3}(s_{o})/η^{2}))

at the time t. In one-dimensional models the water depth is an increasing function of s so t is positive because G is also increasing, see the example below. This shock is a singularity in the derivative of h so we can think of it as a travelling hydrolic jump.

The characteristics that we obtain are not straight lines as for the nonlinear wave equation, due to the term |∇H|^{1/2}. For example when g(s) = |∇H|^{−1/2} = s^{1/4}, G(s) = \frac{4}{5}s^{5/4}, the characteristics are given by the formula

s^{5/4} = s_{o}^{5/4}− \frac{25}{12} \frac{h^{2/3}(s_{o})}{η^{2}}t + O(h^{5/3}(s_{o})/η^{2})

and the shock develops at the time

t = \frac{12}{25} \frac{η^{2}(s_{1}^{5/4}− s_{o}^{5/4})}{h^{2/3}(s_{1}) − h^{2/3}(s_{o})} .

The analysis is similar for small values of R, that can be treated perturbatively, but if the rainfall rate becomes large, it dominates and prevents shocks from developing as shown above. In summary, the water flow equation possesses two mechanisms producing shocks in the water surface, one caused by the nonlinearities in h the other by the singularities in the water (or land) surface H. We may interpret these latter shocks as rapids or knick points, that cause shocks in the water depth and these shocks propagate downstream, as illustrated in Appendix B.

Figure 13: The characteristics for the shock caused by surface singularities.

B Surface Singularities and Shocks

We provide an illustrative example of how shocks are generated, on concave separable surfaces with singularities in their slopes. We note that these examples may be generalized to non-separable surfaces. We consider the one-dimensional hyperbolic equation for the evolution of the variation of the water surface

η^{2} \frac{∂v}{∂t} + \frac{5}{3} \frac{∂}{∂s}(−h1^{\frac{2}{3}}|∇H_{1}|^{1/2}v) = f(x,y,t),  (50)

where

\frac{1}{g(s)} = −\frac{5}{3}h1^{\frac{2}{3}}|∇H_{1}|^{1/2} = −η^{2} \frac{∂q_{w}}{∂h}

is either a strictly decreasing or a strictly increasing function. We consider approximations to initial and separable landscapes in which we assume that the landscape is locally uniform in directions perpendicular to the gradient of the water surface.

Close to a separable landscape uniform in the direction, u_{w}^{⊥}, where u_{w} = −\frac{∇H}{|∇H|},^{24}

g(s) = −s^{−1/4} .  (51)

is a reasonable model. We make the change of variables

a = \frac{5}{3}h1^{\frac{2}{3}}|∇H_{1}|^{1/2}v = −\frac{v}{g(s)}

24The separable landscape is of course not uniform in the direction of u^{⊥} but this example still gives the right information.

then the equation can be written in the form

\frac{∂a}{∂t} + \frac{1}{g(s)η^{2}} \frac{∂a}{∂s} = f(x,y)/η^{2}.  (52)

If we now let

\frac{ds}{dτ} = \frac{1}{g(s)η^{2}} ,

define the characteristics and g = \frac{∂G}{∂s}, then equation (52) can be solved along characteristics,

a(s,t) = a_{0}(G^{−1} (G(s) − \frac{t}{η^{2}} )) + \frac{1}{η^{2}} ∫_{0}^{t}f(G^{−1} (G(s) − \frac{t − τ}{η^{2}}),τ)dτ.

The important equation is the equation for the characteristics

\frac{ds}{dτ} = \frac{1}{g(s)η^{2}} = −\frac{5}{3η^{2}} h1^{\frac{2}{3}}|∇H_{1}|^{1/2}.  (53)

It says that if the derivative of |∇H|^{1/2} has a singularity at the origin then the derivative of h will obtain a singularity there and this singularity will propagate downstream along the characteristics

s = G^{−1} (G(s_{o}) − \frac{t}{η^{2}} ),

see Figure 13. This constitutes the second type of shock formation, the shocks that are caused by singularities, rapids or knick points, in the surface, in distinction to the first type of shock formation discussed in Section 2.4 and Appendix A.

C The Scaling of the Gradient

In this appendix we explain how the presence of a profile such as the initial uniform surface with constant slope, or the profile of the separable landscape, affects the scaling of the gradient of the water surface |∇H|. Consider

|∇H|(y + x) = [H_{x}^{2}(y + x) + H_{y}^{2}(y + x)]^{1/2}
∼ [(ex^{χ} + O(x^{2χ}) + ··· )^{2} + (c + dx^{χ} + O(x^{2χ}))^{2}]^{1/2}
= (c^{2}+ 2cdx^{χ} + O(x^{2χ}))^{1/2},

in a statistical sense where c,d,e, are constants. By the inequality,

(a + b)^{1/2} ≤ a^{1/2}+ b^{1/2},

< (|∇H|(y + x) − |∇H|(y))^{2} >^{1/2}  = (c^{2}+ 2cdx^{χ} + O(x^{2χ}))^{1/2}− c
≤ (2cdx^{χ} + O(x^{2χ}))^{1/2} = x^{χ/2}(2cd + O(x^{χ}))^{1/2}.

Similarly

< (|∇H|(y + x) − |∇H|(y))^{2} >^{1/2}  = (c^{2}+ 2cdx^{χ} + O(x^{2χ}))^{1/2}− c
≥ (2cdx^{χ} + O(x^{2χ}))^{1/2} = −x^{χ/2}(2cd + O(x^{χ}))^{1/2}.

by the inequality

(a + b)^{1/2} ≥ a^{1/2}− b^{1/2}.

This shows that the statistical scaling of |∇H| has the exponent χ/2, whenever χ is the spatial roughness exponent.

D River Basins and Horton’s Laws

Consider the projection of the river network onto the x − y plane. Whereas the roughening of the surface above is controlled by the width function, the lateral location of a channel is determined by a pure random walk in the downridge directions. Since the flow is in the downstream direction, permissible motions at each step occur in an arc subtending 180^{o}. The corresponding exponents for this random walk are

β = \frac{1}{2}, χ = 1, z = 2,

where β is now just the Hurst exponent. Now the tributary channels sit in the 2-dimensional plane and the fractal dimension of the river network is D_{f} = 2 − χ = 1. The corresponding (planar) avalanche dimension is D = 1 + χ = 2. If we take a one-dimensional slice perpendicular to the direction of the channel, then its position on this slice is that of a one-dimensional random walker and the water surface above this one-dimensional slice is determined by another random walker in the vertical direction (see Section 2.5.)

This analysis implies that many of the properties of river networks may be derived from a random walk. For example, the avalanche dimension for a pure random walk also gives us a simple proof of a well-known link between the Horton relations for streams. Consider the Strahler-Horton ordering for streams ^{25} Horton [Horton, 1945] discovered that if N(a) denotes the number of streams of order a and L(a) is their mean length then the ratios

\frac{L(a + 1)}{L(a)} = R_{L}, \frac{N(a)}{N(a + 1)} = R_{B},  (54)

are approximately constant over the same river networks, where R_{L}is called the length ratio and R_{B} is called the bifurcation ratio. In [Barbera and Rosso, 1989], this was related to the avalanche (embedding) dimension of a river network

D = \frac{ln(R_{B})}{ln(R_{L})} = 2,

by the value of D computed above. Thus we obtain

R_{B}= R_{L}^{2}.  (55)

This relationship is borne out by observations where typically R_{B}∼ 4 and R_{L} ∼ 2.

251. streams without tributaries have order one; 2. when two streams of order a join, they form a stream of order a + 1; 3. when two streams of different order join, the resulting stream inherits the higher order of the two.

References

[Bak and Paczuski, 1995] Bak, P. and Paczuski, M. (1995). Complexity, contingency, and criticality. Proceedings of the National Academy of Science, 92:6689–6696.

[Bak et al., 1988] Bak, P., Tang, C., and Wiesenfeld, K. (1988). Self-organized criticality. Physical Review A, 38:364–374.

[Barbera and Rosso, 1989] Barbera, L. and Rosso, R. (1989). On the fractal dimension of river networks. Water Resour. Res, 25(4):735–741.

[Birnir et al., 1999] Birnir, B., Merchant, G., and Smith, T. (1999). Ill-posedness and stochasticity in scientific computations. University of Calfornia, Santa Barbara, preprint.

[Birnir et al., 1998] Birnir, B., Smith, T., and Merchant, G. (1998). Scaling laws and the emergence of channelized drainage patterns in a class of non-linear landscape evolution models. Proceedings, American Geophysical Union Chapman Conference on Fractal Scaling, Non-linear Dynamics, and Chaos in Hydrologic Systems.

[Buldyrev et al., 1992] Buldyrev, S. et al. (1992). Anomolous surface roughening in porous media: Experiment and model. Physical Review A, 45:8313–8316.

[Calderri et al., 1996] Calderri, G., Maritan, A., Banavar, J., Giacometti, A., Rodriguez-Iturbe, I., and Rinaldo, A. (1996). Cellular models for river networks. Massachusettes Institute of Technology, preprint.

[Davis, 1902] Davis, W. (1902). Geographical Essays. Ginn, Boston, 777 pp. (Reprinted 1954, Dover, New York).

[Dingman, 1984] Dingman, S. (1984). Fluvial Hydrology. W. H. Freeman and Company, New York.

[Edwards and Wilkinson, 1982] Edwards, S. and Wilkinson, D. (1982). Proc. Roy. Soc. London, page 17.

[Fennema and Chaudhry, 1990] Fennema, R. and Chaudhry, M. (1990). Explicit methods for 2-d transient free-surface flows. Journal of Hydraulic Engineering, 116:1013–1034.

[Fowler, 1997] Fowler, A. (1997). Mathematical Models in the Applied Sciences. Cambridge University Press, Cambridge.

[Gray, 1961] Gray, D. (1961). Interrelationships of watershed characteristics. Journal of Geophysics Research, 66(4):1215–1223.

[Horton, 1945] Horton, R. (1945). Erosional development of streams and their drainage basins: a hydrophysical approach to quantitative morphology . Geol. Soc. Am. Bull., 56:275–370.

[Howard, 1994] Howard, A. (1994). A detachment-limited model of drainage basin evolution. Water Resources Research, 30:2261–2285.

[Howard et al., 1994] Howard, A., Dietrich, W., and Seidl, M. (1994). Modeling fluvial erosion on regional to continental scales. Journal of Geophysical Research, 99:13971– 13986.

[Julien and Simons, 1985] Julien, P. and Simons, D. (1985). Sediment transport capacity of overland flow . American Society for Agricultural Engineers, 28:755–762.

[Kardar et al., 1986] Kardar, M., Parisi, G., and Zhang, Y.-C. (1986). Dynamic scaling at finite temperatures. Physical Review Letters, 56:889.

[Klinkenberg and Goodchild, 1992] Klinkenberg, B. and Goodchild, M. (1992). The fractal properties of topography. Earth Surf. Processes Landforms, 17:217–234.

[Lavallée et al., 1993] Lavallée, D., Lovejoy, S., Schertzer, D., and Ladoy, P. (1993). Nonlinear variability of landscape topography: Multifractal analysis and simulation. Fractals in Geometry, edited by N. S. Lam and L. De Cola Prentice Hall, Englewood Cliffs, NJ, pages 158–192.

[Leschorn, 1993] Leschorn, H. (1993). Interface depinning in disordered mediumnumerical results. Physica A, 195:324–335.

[Loewenherz-Lawrence, 1991] Loewenherz-Lawrence, D. (1991). Stability and the initiation of channelized surface drainage: a reassessment of the short wavelength limit. Journal of Geophysical Research, 96:8453–8464.

[Mandelbrot, 1983] Mandelbrot, B. (1983). The Fractal Geometry of Nature. Freeman, New York, New York.

[Mandelbrot and Ness, 10] Mandelbrot, B. and Ness, J. V. (10). Fractional brownian motions, fractional noises and applications. SIAM Rev., 1968:422–437.

[Maslov et al., 1994] Maslov, S., Paczuski, M., and Bak, P. (1994). Avalanches and 1/f noise in evolution and growth models. Physical Review Letters, 73:2162–2165.

[Morisawa, 1964] Morisawa, M. (1964). Development of drainage systems on an upraised lake floor. American Journal of Science, 262:340–354.

[Paczuski, 1995] Paczuski, M. (1995). Dynamic scaling: Distinguishing self-organized from generically critical systems. Physical Review E, 52:2137–2140.

[Rodriguez-Iturbe and Rinaldo, 1997] Rodriguez-Iturbe, I. and Rinaldo, A. (1997). Fractal River Basins: Chance and Self-Organization. Cambridge University Press, Cambridge.

[Sapozhnikov and FouFoula-Georgiou, 1996] Sapozhnikov, V. and FouFoula-Georgiou, E. (1996). Do current landscape evolution models show self-organized criticality. Water Resources Research, 32:1109–1112.

[Scheidegger, 1967a] Scheidegger, A. (1967a). A stochastic model for drainage patterns in an intermontane trench. Bulletin of the International Association of Hydrologists, 12:15–20.

[Scheidegger, 1967b] Scheidegger, A. (1967b). A stochastic model for drainage patterns in an intermontane trench. Bulleting of the Association of Scientific Hydrologists, 12:15– 20.

[Schroeder, 1991] Schroeder, M. (1991). Fractals, Chaos, Power Laws. Freeman, New York, New York.

[Schumm et al., 1987] Schumm, S., Mosley, M., and Weaver., W. (1987). Experimental Fluvial Geomorphology. John Wiley and Sons, New York.

[Singh, 1996] Singh, V. (1996). Kinematic Wave Modeling in Water Resources: Surface Water Hydrology. Wiley and Sons, New York.

[Smith et al., 1997a] Smith, T., Birnir, B., and Merchant, G. (1997a). Towards an elementary theory of drainage basin evolution: I. the theoretical basis. Computers and Geosciences, 23(8):811–822.

[Smith and Bretherton, 1972] Smith, T. and Bretherton, F. (1972). Stability and the conservation of mass in drainage-basin evolution . Water Resources Research, 8:1506–1529.

[Smith and Merchant, 1995] Smith, T. and Merchant, G. (1995). Conservation principles and the initiation of channelized surface flows. in Natural and Anthropomorphic Influences in Geomorphology, ed. J. Costa, A. J. Miller, K. W. Potter, and P. Wilcock, pages 1–25.

[Smith et al., 1997b] Smith, T., Merchant, G., and Birnir, B. (1997b). Towards an elementary theory of drainage basin evolution: II. a computational evaluation. Computers and Geosciences, 23(8):823–849.

[Smith et al., 2000] Smith, T., Merchant, G., and Birnir, B. (2000). Transient attractors: towards a theory of the graded stream for alluvial and bedrock channels. Computers and Geosciences, 26(5):531–541.

[Smoller, 1994] Smoller, J. (1994). Shock Waves and Reaction-Diffusion Equations. Springer Verlag, New York.

[Sneppen, 1998] Sneppen, K. (1998). Fractals and intermittancy in dynamics of large systems. Nordic Institute for Theoretical Physics, preprint.

[Takayasy et al., 1991] Takayasy, H., Takayasy, M., Provata, A., and Huber, G. (1991). Statistical models of river networks. Journal of Statistical Physics, 65:725–740.

[Weissel and Pratson, 1994] Weissel, J. and Pratson, L. (1994). The length-scaling of topography. Journal of Geophysical Research, 99:13997–14012.

[Willis, 1924] Willis, J. (1924). Fractals and power laws in biology, sociology, and physics. Nature, February:178.

[Zhang and Cundy, 1989] Zhang, W. and Cundy, T. (1989). Modeling of two-dimensional overland flow . Water Resources Research, 25:2019–2035.

[Zipf, 1949] Zipf, G. (1949). Human Behavior and the Principle of Least Effort. Addison-Wesley, Cambridge, MA.

53 Figure 2: The Water Depth h and the Slope of the Water Surface | ∇H | of a Typical Run at 10% of the Surface Eroded 55 Figure 4: The Water Depth h and the Slope of the Water Surface | ∇H | of a Typical Run at 60% of the Surface Eroded