\documentclass[12pt]{article}
\usepackage{times}
\usepackage{mathptm}
%\usepackage{modeqnos}
\usepackage{epsfig,t1enc,rotating}
\topmargin 0.5 in
\headheight 0in
%\footheight .25in
\oddsidemargin .0 in
\evensidemargin .0 in
\textheight 8.00 in
\textwidth 6.00 in
\headsep=0in
\pagestyle{empty}
\pagenumbering{roman}
\setlength{\parindent}{0cm}
\parskip=.2cm

%\bibliographystyle{apalike}

\setlength{\baselineskip}{5mm}

\pagenumbering{arabic}
\pagestyle{plain}
\setcounter{page}{1}

\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{prop}[thm]{Proposition}
\newtheorem{lem}[thm]{Lemma}

\date{}

\def\ra{\rightarrow}
\def\gD{\Delta}
\def\gd{\delta}
\def\gb{\beta}
\def\gs{\sigma}
\def\ga{\alpha}
\def\go{\omega}
\def\gg{\gamma}
\def\gl{\lambda}
\def\nn{\nonumber}
\def\nin{\noindent}
\def\ul{\underline}
\def\Bbb R{{\rm I\!R}}


\begin{document}

\begin{center}

{\bf {\Large THE SCALING OF FLUVIAL LANDSCAPES}}

\vspace{.4in}

{\bf Bj\"{o}rn Birnir~\,$\!^\ddagger$~~~
Terence R. Smith~\,$\!^\dagger$\,$\!^*$~~~
George E. Merchant~\,$\!^\dagger$}
\end{center}


{\small
$^\ddagger$ Department of Mathematics,
University of California at Santa Barbara,
California 93106, USA
and University of Iceland, Science
Institute, 3 Dunhaga, 107 
\ Reykjav{\'{\i}}k\\
$^\dagger$ 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}

\vspace{.4in}

%\maketitle

\vspace{.4in}

\begin{center}
ABSTRACT
\end{center}


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 {\em 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 {\em 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. 

\newpage

{\section{Introduction}}
\label{sec:intro}

In three preceding papers \cite{SBM97,SMB97, SMB99},
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 \cite{RR97}, 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 {\em 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 \cite{SBM97,SMB97}
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 \cite{S91,Z49, W24,M83}.
It is now known that the evolution of surfaces
whose dynamics are driven by various forms of noise
are frequently characterized by {\em 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 {\em universality classes} of,
behaviors (see, for example, Sneppen \cite{Sn98}.)
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 \cite{SBM97,SMB97} as systems driven
by noise or stochastic processes.
First, we characterize the statistical structure
of eroding surfaces and flows in terms
of various {\em 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. \cite{BTW88,BP95,P95,Sn98}.
Although this concept has yet to find an appropriate
mathematical formulation, it is based on
the underlying idea of an {\em 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 \cite{SFG96}
and by specific asymptotic forms of the structure functions
\footnote{We believe, in fact,
that our theory may provide some
of the foundation for a {\em mathematical theory}
of self-organized criticality.}.

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 {\em 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 {\em 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.

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 \cite{WP94},
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
\footnote{In 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.}.

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 
{\em early channel initiation} involves the formation of {\em shocks}
or hydrolic jumps
\footnote{A 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 \cite{Smo94}).}
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 {\em 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.

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 \cite{SBM97},
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 {\em 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 \ref{sec:model} 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 \ref{sec:numer},
provide a basis for the empirical
and analytical investigations presented
in the remaining sections of the paper.
In Section \ref{sec:channels},
we provide a theoretical explanation for the first scaling 
relation that is associated with the emergence
of channelization on unstructured surfaces. 
In Section \ref{sec:selfsim},
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 \ref{sec:disc},
while Section \ref{sec:summary} 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 
\ref{sec:shockequation}.
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  \ref{sec:sr} and \ref{sec:grad} respectively and 
a derivation of a relationship between Horton's laws is contained in
Appendix \ref{sec:Horton}.

\section{Modeling Eroding Surfaces and their Scaling Laws}
\label{sec:model}

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.

\subsection{A Family of Landscape Evolution Models}
\label{sec:family}

The scaling results presented in this paper are derived
from the same family of models studied analytically and numerically
by Smith, Birnir, and Merchant \cite{SBM97},
Smith, Merchant, and Birnir \cite{SMB97},
and Birnir, Smith, and Merchant \cite{BSM98}.
These models represent the {\em advective entrainment 
and transport of sediment in transport limited} 
conditions \cite{Ho94}.
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 \cite{HDS94})
\footnote{Whether they capture essential aspects
at large scales of spatial resolution is an open issue.}.
While it is straightforward to extend 
these models to represent other processes 
that are significant in the evolution of fluvial landscapes,
such as gravity-driven 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
\footnote{We plan to explore the effects of such processes
on scaling in later papers.}.

We now provide a summary description
of the family of models.
A more detailed description of their derivation
is provided in \cite{SBM97}.
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)$
\begin{eqnarray}
\label{eq:m1}
\vspace{.1in}
~-R ~+~\frac{\partial h}{\partial t}
~=~\!\nabla \cdot \left( {\bf u}_w ~q_w \right),~~~~~~
\vspace{.1in}
\end{eqnarray}
\begin{eqnarray}
\label{eq:m2}
\vspace{.1in}
\frac{\partial z}{\partial t}
~=~\!\nabla \cdot \left( {\bf u}_w ~ q_s \right)
\vspace{.1in}
\end{eqnarray}
in which $h=h(x,y,t)$ is the depth of water
varying continuously over the landsurface,
$-{\bf u}_w = -\frac{\!\nabla H }{\!\mid \!\nabla H \!\mid}$
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 (\ref{eq:m1}), (\ref{eq:m2}).
These time scales
may be derived by transforming the variables 
of equations (\ref{eq:m1}), (\ref{eq:m2}) 
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
(\ref{eq:m1}), (\ref{eq:m2}),
\footnote{A useful set of such values
characterize the representative scales
when the basic model (\ref{eq:m1}, \ref{eq:m2})
is used to represent the erosion of a small,
linear ridge.}
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.

A {\em 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 (\ref{eq:m1}), (\ref{eq:m2})
and dropping the asterisks denoting dimensionless variables,
we obtain
\begin{eqnarray}
\vspace{.1in}
\frac{\partial h}{\partial t}
&=& \!\nabla \cdot \left( {\bf u}_w ~q_w \right)
+ \frac{[x][R]}{[q_w]} R \nonumber \\
\vspace{.1in}
\frac{\partial H}{\partial t} -\frac {[h]}{[H]} \frac{\partial h}{\partial t} 
&=& \frac{[h][q_s]}{[H][q_w]} \!\nabla \cdot \left( {\bf u}_w ~ q_s \right).
\nonumber
\vspace{.1in}
\end{eqnarray}
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, \cite{Fo97}).
The significance of this fact, as we show in
Section \ref{sec:noise} and in Appendix \ref{sec:shockequation},
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{\partial H}{\partial t} \approx
\frac{[h]}{[H]} \frac{\partial h}{\partial t} \]
which indicates that variations in
the water surface $H$ are {\em 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 {\em 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 (\ref{eq:m1}), (\ref{eq:m2}) we
obtain
\begin{eqnarray}
\vspace{.1in}
\frac{[q_s]}{[q_w]} \frac{\partial h}{\partial t}
&=& \!\nabla \cdot \left( {\bf u}_w ~q_w \right)
+ R \nonumber \\
\vspace{.1in}
\frac{\partial H}{\partial t} 
&=& \frac{[h]}{[H]}  \!\nabla \cdot \left( {\bf u}_w ~ q_s \right)
+  \frac{[h]}{[H]} \frac{\partial h}{\partial t}
\nonumber
\vspace{.1in}
\end{eqnarray}
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{\partial h}{\partial 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
\cite{Fo97}, 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 {\em 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 (\ref{eq:m1}), (\ref{eq:m2}) 
and again using the scaling relation $[x][R]/[q_w] = 1$,
we obtain
\begin{eqnarray}
\vspace{.1in}
\frac{[h][q_s]}{[H][q_w]} \frac{\partial h}{\partial t}
&=& \!\nabla \cdot \left( {\bf u}_w ~q_w \right)
+ \frac{[x][R]}{[q_w]} R \nonumber \\
\vspace{.1in}
\frac{\partial H}{\partial t} 
&=& \!\nabla \cdot \left( {\bf u}_w ~ q_s \right)
+  \frac {[h]}{[H]} \frac{\partial h}{\partial t}.
\nonumber
\vspace{.1in}
\end{eqnarray}
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
\begin{eqnarray}
\label{eq:long}
\vspace{.1in}
0 &\approx& \!\nabla \cdot \left( {\bf u}_w ~q_w \right)
+ \frac{[x][R]}{[q_w]} R 
\vspace{.1in}
\end{eqnarray}
as was done, for example, in the analyses presented
in \cite{SB72,Lo91,SBM97,SMB97}.
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] {\partial h}/{\partial t}$
as representing random variations that drive
the sediment flow.

\subsection{Previous Results for Long Time Scale Approximations}
\label{sec:prevanal}

Smith, Birnir, and Merchant \cite{SBM97,SMB97} analyzed
a subfamily of the models (\ref{eq:m1}), (\ref{eq:m2})
that was obtained with the use of the following assumption.
First, the long time scale approximation (\ref{eq:long})
to the basic model was assumed, together with the additional
scaling relationship $[h]/[H] = [q_s]/[q_w] \equiv \eta$.
Second, a Manning-type constitutive relation 
\begin{eqnarray}
\label{eq:manning}
q_w = n h^\frac{5}{3}{\!\mid \!\nabla H \!\mid }^\frac{1}{2}
\end{eqnarray}
where n, the inverse of the channel roughness, is a constant,
was adopted to represent the flux of water. 
Third, a power law relation 
\begin{eqnarray}
\label{eq:power}
q_s = k h^{\frac{5}{3}\gamma} {\!\mid \!\nabla H \!\mid }^{\frac{\beta}{2} + \delta
}\end{eqnarray}
was employed in representing the flux of sediment $q_s$ \cite{JS85}.
These assumptions lead to the following dimensionless equations:
\begin{eqnarray}
\vspace{.1in}
\label{eq:m3}
-R ~=~\!\nabla \cdot \left[ \frac{\!\nabla H }{\!\mid \!\nabla H \!\mid}
h^\frac{5}{3}~{\!\mid \!\nabla H \!\mid }^\frac{1}{2} \right]
\vspace{.1in}
\end{eqnarray}

\begin{eqnarray}
\vspace{.1in}
\label{eq:m4}
\frac{\partial H}{\partial t}-\eta\frac{\partial h}{\partial t}
~=~\!\nabla \cdot \left[ \frac{\!\nabla H }{\!\mid \!\nabla H \!\mid }
h^{\frac{5}{3} \gamma} ~{\!\mid \!\nabla H
\!\mid}^{\frac{\gamma}{2} + \delta} \right].
\vspace{.1in}
\end{eqnarray}
The initial and boundary conditions employed in these
analyses modeled a linear ridge 
undergoing erosion as a result of a uniform rainfall
(see equations (\ref{eq:m6}), (\ref{eq:m7}) below.)

Applying the numerical methods described in \cite{SMB97}
to equations (\ref{eq:m3})-(\ref{eq:m4}), it was shown
\cite{SMB97} that initially unstructured but 
randomly perturbed two-dimensional 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
\footnote{A variety of numerical experiments
have shown that channelization occurs
on surfaces with a large array of initial configurations.}.
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
\footnote{We note, in particular, the sharpness of the 
longitudinal ridges, or divides, separating the valleys.}.\\




The convergence to these forms occurs from initial surfaces that are
both planar and non-planar.
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 $\gamma$ and $\delta$.
These mature landscapes, consisting of stable patterns 
of concavities, valleys, ridges, and associated flows,
are well-described by a class of solutions to the nonlinear 
PDEs (\ref{eq:m3})-(\ref{eq:m4}) for cases in which
$ \eta \frac{\partial h}{~~\partial t}$
may be ignored in equation (\ref{eq:m4})
relative to $\frac{\partial H}{\partial t}$ \cite{SBM97}.
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
\begin{eqnarray}
\label{eq:sep}
\vspace{.1in}
H(x, y, t) ~=~T_{H}(t) H_o(x, y),~~~~~ h(x, y, t) ~=~T_{h}(t) h_o(x, y);
\vspace{.1in}
\end{eqnarray}
and they are stable \cite{SBM97}.
They are also characterized in terms of a
variational principle \cite{SBM97}, 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.

\subsection{The Models Underlying the Present Analysis}

The analysis presented in the current paper
is similarly based on the use of equations (\ref{eq:m1}), (\ref{eq:m2}),
the Manning-type constitutive relation
(\ref{eq:manning}) for the flux of water,
a power law relation (\ref{eq:power}) for the flux of sediment,
and the initial and boundary conditions prescribed in equations
(\ref{eq:initial}), (\ref{eq:boundary}).
For reasons that we discuss below, however,
{\em we no  longer assume that the term
$\frac{\partial h}{\partial t}$ in equation (\ref{eq:m1})
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 \ref{sec:family}.
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
\begin{eqnarray}
\vspace{.1in}
\label{eq:m6}
 {\eta^2} \frac{\partial h }{\partial t}~=~\!\nabla \cdot \left[ \frac{\!\nabla H }{\!\mid \!\nabla H \!\mid}
h^\frac{5}{3}~{\!\mid \!\nabla H \!\mid }^\frac{1}{2} \right]~+~R,
\vspace{.1in}
\end{eqnarray}

\begin{eqnarray}
\vspace{.1in}
\label{eq:m7}
\frac{\partial H}{\partial t}
~-~ {\eta} \frac{\partial h }{\partial t}
~=~\!\nabla \cdot \left[ \frac{\!\nabla H }{\!\mid \!\nabla H \!\mid }
h^{\frac{5}{3} \gamma} ~{\!\mid \!\nabla H
\!\mid}^{\frac{\gamma}{2} + \delta} \right]
\vspace{.1in}
\end{eqnarray}
in which we have adopted the scaling relation
$[h]/[H]=[q_w]/[q_s] \equiv \eta$.
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 $\tau = t/\eta^2$
and the intermediate time scale version
by applying the the transformation $\tau = t/\eta$.

We employ the same boundary and initial conditions
as \cite{SBM97,SMB97} 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)\in {\Bbb R}^2| 0 \leq x \leq L, \; 0 \leq y \leq W \},
$$
with initial conditions
\begin{eqnarray}
\label{eq:initial}
h(x,y,0) &=& d(y),~~ d(0)=h_o,~~ d(W)=0,\nonumber\\
H(x,y,0) &=& cy+h_o,~~ \; 0\le y \le W
\end{eqnarray}
and boundary conditions
\begin{eqnarray}
\label{eq:boundary}
h(x,W,t) &=& 0,\nonumber\\
H(x,0,t) = &h_0& = h(x,0,t)
\end{eqnarray}
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 $\nabla 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\cdot \nabla H(x,w,t) = \infty.
$$ 
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.

\subsection{The Role of Noise and Its Generation}
\label{sec:noise}

The reasons for retaining the term 
${\eta^2} \frac{\partial h }{\partial t}$
in the water flow equation (\ref{eq:m6}) 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 \cite{Sn98} 
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 (\ref{eq:m6}), 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 \cite{SBM97,SMB97}
and represented in  (\ref{eq:m6}), (\ref{eq:m7}) 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 \cite{SMB97},
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 \cite{SMB97},
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 \cite{Mo64} and Schumm, Mosley, and Weaver \cite{SMW87})
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 (\ref{eq:m6}) are capable of
generating variations in flows at short time scales 
that may be interpreted as constituting noise.
In particular, we note that equation (\ref{eq:m6})
is capable of generating {\em 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 \cite{Di84,Si96}. 
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, \cite{Di84,Si96}).

In Appendix \ref{sec:shockequation} we show that the water flow 
equation (\ref{eq:m6}) becomes a one-dimensional equation describing the
flow of water down the gradient of the water surface
\begin{eqnarray}
\vspace{.1in}
\label{eq:s0}
  \frac{\partial h }{\partial t}~-~
\frac{5}{3\eta^2}h^{2/3}\ {\!\mid \!\nabla H \!\mid }^\frac{1}{2}\frac{\partial
h 
}
{\partial s}  ~&=&~\left[{{{\!\mid \!\nabla H \!\mid }^{2}\Delta
H-\frac{1}{2}(H_x^2H_{xx}+2H_x H_y H_{xy}+H_y^2H_{yy})}
\over{{\!\mid \!\nabla H \!\mid }^\frac{5}{2}}}
\right]{{h^{5/3}}\over{\eta^2}}\nonumber\\
~&+&~\frac{R}{\eta^2}.
\vspace{.1in}
\end{eqnarray}
in which $s$ parametrizes this direction.
We provide a proof of the existence of shocks in the flows
described in this equation in Appendix \ref{sec:shockequation},
and show below that such shocks 
lead to white noise that seeds the instabilities in the flow of 
sediment in equation (\ref{eq:m7}) on a longer time-scale.
Two distinct mechanisms lead to such noise.
First, on convex surfaces the water flow equation (\ref{eq:s0})
has solutions that develop shocks.
Second, the sediment flow equation (\ref{eq:m7}) 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 \ref{sec:nummeth})
are more accurate than the procedures used in \cite{SMB97}, 
they do not require the suppression of machine noise with the use of
large amounts of {\em 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 \ref{sec:family},
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 $\tau = t/\eta^2$ and
$\tau = t/\eta$, their cumulative 
effect is felt in the sediment flow
on the long time-scale of erosion through the 
dependence on $h$ and the term 
$-\eta {{\partial h}\over{\partial 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 \ref{sec:family}.
{\em These effects constitute the random variations
in the model (\ref{eq:m6})-(\ref{eq:m7})
that drive the dynamics of the system.}

\subsection{Complex Surfaces, Structure Functions, and Scaling Relations}
\label{sec:struct}

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
{\em structure functions} that take the form
\begin{eqnarray}
\rho_k ({\bf x}, t) \,  =  \, <|F({\bf x+ y},t)-F({\bf y},t)|^k>
\end{eqnarray}
in which $F({\bf x},t)$ is some function characterizing a surface,
${\bf 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 {\em 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 $\rho_1({\bf 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 
\begin{eqnarray}
\label{eq:e1}
W_F ({\bf x},t) \equiv  \rho_2 ({\bf x}, t)^\frac{1}{2}
\,  =  \, {<|F({\bf x+ y},t)-F({\bf y},t)|^2>}^\frac{1}{2}
\end{eqnarray}
to represent the root mean square of the elevation differences
as a function of different distances of separation $\mid{\bf x}\mid $.
This function, which is known variously
as the {\em variogram, height-height correlation function},
{\em roughness function}, or {\em width function}
\cite{EW82, KPZ86,BBC92},
is often used to characterize the roughness of the surface
in the direction of growth.
We may also interpret $W_F ({\bf x},t)$ as measuring
the correlation between heights at points
separated by a distance $\mid{\bf x}\mid$. 
We note that
the width function and the spatial autocorrelation function 
\begin{eqnarray}
\label{eq:e2}
C_F({\bf x}, t) \, = \, <F({\bf y+x},t) F({\bf y},t)>
\end{eqnarray}
are related by 
\begin{eqnarray}
\label{eq:e3}
C_F({\bf x}, t) = \sigma_F^2(t) -{1\over 2} W_F^2({\bf x},t),
\end{eqnarray}
in which $\sigma_F$ is the variance 
\footnote{In practice the variance
is always finite for topographic surfaces.}.
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
$\sigma_F(t) \sim t^{1/2}$ but by $W_F \sim t^{1/4}$.

Important questions concern the manner in which
structure functions and their analogs
vary with changes in the values
of the lag variable ${\bf 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 {\em scaling invariances}
in the system.
It is often the case that the behavior of a surface
is well described by asymptotic forms of (\ref{eq:e1}) 
that are applicable in different space-time regimes
of the system's dynamics
(see, for example, \cite{BBC92,P95,Sn98}).
Such forms have the representation
\begin{eqnarray}
\label{eq:e4}
W_F({\bf x},t)  \sim t^\beta f\left({t/{|{\bf x}|^z}}\right),
\end{eqnarray}
in which $t$ is time, ${\bf x}$ is the lag variable,
$\beta$ is known as {\em the temporal roughness exponent},
and $z$ is called {\em the dynamic exponent}.
An interpretation of these
scaling relations is that there is an 
{\em equivalence between the temporal
and the spatial scales given by $t \sim |{\bf 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 {\em roughness} of a surface, as measured by 
$W_F \sim  t^\beta  f(t/{\bf x}^z)$, 
typically increases. 
This regime is characterized by
$f(y) \sim  y^{-\beta}$ for $y << 1$,
from which it follows that the width function
may be approximated by a simple power law
\begin{eqnarray}
\label{eq:e5}
W_F({\bf x},t)  \sim t^\beta
\end{eqnarray}
in which $\beta$ is the temporal roughness exponent.
Since the temporal and 
 spatial roughness exponents are related by
\footnote{This indicates that the dynamic exponent $z$
may be interpreted as the constant of proportionality
linking the temporal and spatial scales.}
\begin{equation}
\label{eq:e7}
{\chi} = z \gb,
\end{equation}
the width function in this case may be approximated 
by the relationship
\begin{eqnarray}
\label{eq:e6}
W_F({\bf x},t)  \sim {\mid {\bf x} \mid }^\chi
\end{eqnarray}
in which $\chi$ is {\em 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 \sim L^z$, where $L$ is the size of the system
\footnote{This indicates that
the system takes time $t=L^z$ to completely reorganize itself.}.
Following this transient period, 
there may be a space-time regime of stationary,
or {\em saturated}, behavior
in which the roughness of the surface fluctuates
about statistically steady values.
In this regime, which is characterized by
$f(y) \sim$ {\em constant},$~y>>1 $,
pairs of points on the surface that are separated by
a distance $|{\bf x}| \sim L$ are characterized by
differences in the expected values of the width function 
of order of $|{\bf x}|^\chi$.

In surface growth models in which
$W_F({\bf x},t)$ measures the roughness of the surface 
in the direction of the growth,
$\chi$ may be interpreted as the fractional part 
of the dimension of the surface,
i.e., the embedding or avalanche dimension is
$d + \chi$, where $d$ is an integer.
Roughly speaking the dimension $d$ 
characterizes the spatial extent of the avalanche and
$\chi$ characterizes the height of its ``bulge'', see Sneppen \cite{Sn98}.

If the roughness exponent $\chi$ does not depend linearly on $k$ then
the scaling is called
{\em multi-fractal}, while if
the scaling exponent depends
on the direction of the lag vector,
the scaling is said to be
{\em anisotropic} (see, for example, \cite{KG92}, \cite{LD93}).
For an arbitrary structure function of order
$k$, for example, this may take the form
$ \rho_k({ (x,y)}) = constant \times (|{ x}|^{k \chi_x}+ |{ y}|^{k \chi_y})$
in which $\chi_x \neq \chi_y$. 
In the case of single-fractal and isotropic scalings,
only two scaling exponents are required to determine 
all others (see, for example, \cite{Sn98}).
The {\em fractal dimension}, for example,
may be represented as $D_f = n-\chi$,
where $n$ is the total (spatial) dimension, 
while the {\em embedding} or {\em avalanche} dimension 
may be represented by $D = d + \chi$,
where $d$ is the ``dimension'' of the interface 
\footnote{1 for a fractal line, 2 for
a fractal surface, etc..}.
The avalanche dimension relates the spatial extent of 
an avalanche to its temporal growth, 
$ |{\bf x}| \sim t^D$
\footnote{We may define an avalanche as a 
burst of activity in surface growth, 
in which a significant number of sites move together.}.

An important class of statistical phenomena may be characterized
in terms of processes that take the form
of fractional Brownian motion \cite{MaVanNess68}.
The variance of such processes,
which are of relevance for the analysis presented in Section \ref{sec:mature},
satisfies scaling laws of the form
\begin{eqnarray}
\label{eq:e8}
\sigma  \sim t^H f(t/|{\bf x}|^{z})
\end{eqnarray}
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 $\gb$ 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 \neq 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 $\gb$.
It follows that for Brownian motion ($H=1/2$) the (generalized) time
derivative is statistically stationary {\em white noise}. 

\subsection{Estimating the Scaling Exponents from Experimental Data}
\label{sec:exp}

In the analysis presented below,
we employ width functions to characterize
the correlation structure of:
(1) the {\em water depth} $h$ in terms of the
width function $W_{h}$;
(2) the {\em water surface elevation} $H$ in terms of 
the width function $W_{H}$; and
(3) the {\em water surface slope} $\mid \nabla H \mid$ 
in terms of the width function $W_{\mid \nabla H \mid }$.
Because of the relation (\ref{eq:e7}),
we only require two exponents 
in characterizing scaling relations involving
these width functions and
we employ the temporal and spatial roughness
exponents $\gb$ and $\chi$. The spatial 
roughness exponent is computed from the scaling
of the width function (\ref{eq:e1}).
We also make use of the fact that we may
represent various scaling relations
in terms of scaling of the (temporal)
power spectrum, see \cite{MPB94}.
For present purposes, it suffices to note that
\begin{equation}
\label{eq:e9}
\hat C ~=~ S(k) \sim \frac{1}{k^{1+2\gb}}
\end{equation}
in which $S(k)$ represents the power spectrum
at frequency $k$.
Employing the power spectrum, we may compute
the temporal roughness exponent $\gb$ independently of the spatial
roughness coefficient $\chi$, and thus also obtain the dynamic
coefficient $z = \chi/\gb$ using equation (\ref{eq:e7}).

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.


\section{The Two Scalings: Numerical Results}
\label{sec:numer}

Numerical simulations of the evolution
of eroding surfaces generated by equations
(\ref{eq:m6}), (\ref{eq:m7})
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.

\subsection{Numerical Methods for Accurate Scaling Relations}
\label{sec:nummeth}

While the numerical methods employed
by Smith, Merchant, and Birnir \cite{SMB97}
in solving equations (\ref{eq:m3}), (\ref{eq:m4})
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 \cite{SMB97},
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 \cite{ZC89,SM95}.
This scheme involves the use of artificial viscosity
to avoid numerical oscillations \cite{FC90}.
By suppressing the noise that is amplified 
by the instability of the initial surface, however
the {\em 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 (\ref{eq:m6}) in two dimensions
using an {\em upwind scheme}, we render 
the scheme stable and eliminate the
need for artificial viscosity.\footnote{An upwind scheme still produces
small amounts of numerical viscosity but of magnitude less than the 
artificial viscosity in \cite{FC90}.} 
A Crank-Nicholson scheme may then be applied in computing 
the nonlinear gradient flow (\ref{eq:m7})
 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
\footnote{We implemented the scheme on a parallel
supercomputer at the University of California, Santa Barbara,
employing as a final step the implementation of
a {\em four color scheme} to parallelize the procedure for maximum 
computational advantage.}.

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 \cite{BMS99}
for further details on the new numerical scheme. 
This code will also be made public. 

\subsection{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 
(\ref{eq:boundary}), (\ref{eq:initial}),
which are the same as those employed in \cite{SBM97,SMB97}. 
These experiments were run using various refinements of the grids
and over all relevant parameter values of the model
(\ref{eq:m6}, \ref{eq:m7}) to provide a basis
for deriving the scaling relationships that we describe
in Section \ref{sec:results}.

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 $\gamma=2,\delta=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 
$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
$\mid \nabla H \mid $, in Figures 2 and 4, again
at 10\% and 60\% of the surface eroded.

\subsection{The Basic Scaling Results}
\label{sec:results}

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

Scaling results were obtained by computing
the width function (\ref{eq:e1}) 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 $\chi$ of
the asymptotic forms (\ref{eq:e6}) 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 $\chi$ were obtained by least squares fits.
In Figures \ref{fig:logdepth} and \ref{fig:logsurf},
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.

\setcounter{figure}{4}
\begin{figure}[tb]
\centerline{\psfig{figure=Figures/figure53.ps,height=3in}}
\vspace{-7pt}
\caption{Log-log plot of water depth at a given cross-section}
\vspace{-7pt}
\label{fig:logdepth}
\end{figure}
\vspace{0.1in}

\begin{figure}[tb]
\centerline{\psfig{figure=Figures/figure7.ps,height=3in}}
\vspace{-7pt}
\caption{Log-log plot of the water surface at a given cross-section}
\vspace{-7pt}
\label{fig:logsurf}
\end{figure}
\vspace{0.1in}

It is convenient to view the values 
of the scaling parameters as a function of their
location 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 $\chi_x = \chi_x(y)$.
We may then represent the scaling relations
found over sections of the ridge by plotting
both the estimated value of $\chi$
and the associated $R^2$-value as a function
of location on the ridge.
This is illustrated in 
Figures \ref{fig:depth}-\ref{fig:grad2}, 
in which we represent  the scalings
for $h$, $H$, and  $\mid\nabla H\mid$,
at different cross-sections of the ridge at the two
erosion times.

\begin{figure}[tb]
\centerline{\psfig{figure=Figures/figure1.10.ps,height=3in}}
\vspace{-7pt}
\caption{The scaling for the water depth}
\vspace{-7pt}
\label{fig:depth}
\end{figure}  

\begin{figure}[tb]
\centerline{\psfig{figure=Figures/figure2.10.ps,height=3in}}
\vspace{-7pt}
\caption{The first scaling of the gradient of the water surface}
\vspace{-7pt}
\label{fig:grad1}
\end{figure}
  
The scaling results from solutions to
(\ref{eq:m6}) and (\ref{eq:m7}) for a time
of 10\% eroded are shown in Figures 
\ref{fig:depth} and \ref{fig:grad1} 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_{\mid\nabla H\mid}$)
scale with the roughness exponent  
$\chi_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 \ref{fig:grad1}
the {\em numerically computed} value of the exponent of 
$|\nabla H|$ is $\chi = 0.25$, the {\em 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 $|\nabla H|$ 
and then takes the square root, as explained in Appendix \ref{sec:grad}. 

As the solutions converge towards the separable forms 
defining the mature landscape of equation (\ref{eq:sep}), 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 \ref{fig:surface} with the scaling exponent $\chi_x = 0.7$,
again over a large range of $y$-locations on the ridge.
While the corresponding {\em numerically computed} scaling 
of $|\nabla H|$ is $\chi_x = 0.35$, we note again that
the {\em true} value of this exponent is $0.7$
rather than the smaller value, for the reasons stated above
and explained in Appendix \ref{sec:grad}, the role of the constant slope
$c$ now being played by the gradient of the separable water 
surface $\nabla H_2$.  

We note, therefore, that at 10\% of the landsurface eroded 
the water depth $h$ and the gradient of the 
water surface $|\nabla H|$ scale, 
whereas at 60\% of the landsurface eroded the water surface $H$
itself and the gradient of the water surface $|\nabla H|$ show the same 
scaling. In both cases $|\nabla 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.

\begin{figure}[tb]
\centerline{\psfig{figure=Figures/figure3.60.ps,height=3in}}
\vspace{-7pt}
\caption{The scaling of the water surface}
\vspace{-7pt}
\label{fig:surface}
\end{figure}
\vspace{0.1in}

\begin{figure}[tb]
\centerline{\psfig{figure=Figures/figure4.60.ps,height=3in}}
\vspace{-7pt}
\caption{The second scaling of the gradient of the water surface}
\vspace{-7pt}
\label{fig:grad2}
\end{figure}

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 \ref{fig:depth}.
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 $\chi_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 \ref{fig:surface}.
Thus the existence of a plateau in $\chi_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.  

The scaling is anisotropic in the sense that there is no plateau in $\chi_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
\ref{sec:grad}, in various directions but found 
little evidence of significant anisotropy. 
We also examined both higher and lower 
structure functions $\rho_k,~k\neq 2$, 
or moments, of $h,H$ and $\nabla H$  
but found no evidence that the first scaling is multifractal. 
The second scaling on the other hand seems to be multi-fractal
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$.



\section{The Emergence of Channels and the First Fundamental Scaling of Fluvial
Landscapes}
\label{sec:channels}

We now consider the emergence, 
early in the evolution of surfaces determined 
by the model (\ref{eq:m6})-(\ref{eq:m7}), 
of scalings in the water depth $h$ and
in the slope of the water surface $\mid \nabla H \mid$  
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_{\mid \nabla H \mid}$
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 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.

\subsection{The Emergence of Channels and the Fundamental
Channelization Instability}
\label{sec:prev}

Results from previous analyses of models 
of the form (\ref{eq:m1}), (\ref{eq:m2}),
provide insight into the earliest stages of channel evolution.
Smith and Bretherton \cite{SB72} analyzed versions
of equations (\ref{eq:m3}), (\ref{eq:m4})
in which the constitutive relation for sediment transport
was generalized to $q_s=F(q_w,\mid \nabla H \mid)$,
with $\frac{\partial F}{\partial q_w} > 0$
and $\frac{\partial F}{\partial S} > 0$
\footnote{The term $\frac{\partial h}{\partial t}$ in equation
(\ref{eq:m4}) was not considered in these analyses.}.
Based on a linearized analysis of these equations,
they showed that if $\frac{F}{q_w} < \frac{\partial F}{\partial 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 $\gamma > 1$
for the particular sediment transport law embedded
in equation (\ref{eq:power}).
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 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 \cite{Lo91},
who also derived the full asymptotic solution 
in the small wavelength limit.

The simulations of Smith, Merchant, and Birnir \cite{SMB97}
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 
(\ref{eq:m1}), (\ref{eq:m2}) which provides additional
insight into the underlying processes.
In particular, we linearize the equations (\ref{eq:m6}) and 
(\ref{eq:m7}) about an initial  surface, linear in $y$ and uniform
in the $x$ direction,  and a constant water depth 
$H= cy+d, h = d$
\footnote{Strictly 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.}, 
where $c$ and $d$ are constants.
We note that, while these initial conditions differ from
those in the stability analyses of \cite{SB72,Lo91}, 
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.

We let 
$$
h = d  + \epsilon a(x,y,t), ~~~H = cy + d + \epsilon b(x,y,t),
$$ 
and include a (very) small stochastic variation in the rainfall rate
$$
R = R_o + \epsilon r(x,y,t).
$$
The linearized system of equations becomes
\begin{equation}
\label{eq:flat}
\vspace{.1in}
 {{\eta}^2} \frac{\partial a}{\partial t} ~=~\frac{5}{3}d^\frac{2}{3}
 c^\frac{1}{2}\frac{\partial a}{\partial y}
~+~{{d^\frac{5}{3}}\over{c^\frac{1}{2}}}(\frac{\partial^2 b}{\partial x^2}
+\frac{1}{2}\frac{\partial^2 b}{\partial y^2})
 ~+~ r,
\end{equation}

\begin{eqnarray}
\vspace{.1in}
\label{eq:notflat}
\frac{\partial b}{\partial t} 
~=~\frac{5}{3}d^\frac{2}{3}
c^\frac{1}{2}({1\over\eta} + 2 c^\frac{5}{2}d^\frac{2}{3})\frac{\partial
a}{\partial y}~&+&~
{{d^\frac{5}{3}}\over{c^\frac{1}{2}}}
({1\over\eta} + c^\frac{5}{2}d^\frac{5}{3})\frac{\partial^2 b}{\partial x^2}\nonumber\\
~&+&~{{d^\frac{5}{3}}\over{c^\frac{1}{2}}}
({1\over{2\eta}} + \frac{5}{2} c^\frac{5}{2}d^\frac{5}{3})\frac{\partial^2 b}{\partial
y^2}
~+~{r\over\eta}.
\vspace{.1in}
\end{eqnarray}

The numerical results reported in Section \ref{sec:results}
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 $\nabla b$ does scale
from which it follows that the equations generating the 
scaling must be equation (\ref{eq:flat}) and the equation
\footnote{We have merely applied the $\nabla$ operator
to equation \ref{eq:notflat}.}
\begin{eqnarray}
\label{eq:grad}
\vspace{.1in}
\frac{\partial \nabla b}{\partial t} 
~=~\frac{5}{3}d^\frac{2}{3}
c^\frac{1}{2}({1\over\eta} + 2 c^\frac{5}{2}d^\frac{2}{3})\frac{\partial
\nabla a}{\partial y}~&+&~
{{d^\frac{5}{3}}\over{c^\frac{1}{2}}}
({1\over\eta} + c^\frac{5}{2}d^\frac{5}{3})\frac{\partial^2 \nabla b}{\partial
x^2}\nonumber\\
~&+&~{{d^\frac{5}{3}}\over{c^\frac{1}{2}}}
({1\over{2\eta}} + \frac{5}{2} c^\frac{5}{2}d^\frac{5}{3})\frac{\partial^2
\nabla b}{\partial
y^2}
~+~{{\nabla r}\over\eta}.
\vspace{.1in}
\end{eqnarray}
Since $\eta = d/100c \sim h/H$ is small, the term ${\nabla r}/\eta$ 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 ${{\partial h}\over{\partial t}} = 0$ 
in the equations (\ref{eq:m6}) and (\ref{eq:m7}), 
or equivalently use the slow equations (\ref{eq:m3}) and
(\ref{eq:m4}), with ${{\partial h}\over{\partial t}} = 0$. 
On taking the gradient, equation  (\ref{eq:flat}) becomes
\begin{equation}
\label{eq:blow}
\vspace{.1in}
  \frac{5}{3}d^\frac{2}{3}
 c^\frac{1}{2}\frac{\partial \nabla a}{\partial y}~=~-
{{d^\frac{5}{3}}\over{c^\frac{1}{2}}}(\frac{\partial^2 \nabla b}{\partial x^2}
+\frac{1}{2}\frac{\partial^2 \nabla b}{\partial y^2})
\end{equation}
and a substitution into (\ref{eq:grad}) gives the equation
\begin{equation}
\label{eq:backheat}
\vspace{.1in}
\frac{\partial \nabla b}{\partial t} 
~=~\frac{5}{3}d^\frac{2}{3}
c^\frac{1}{2}[(-2+d)\frac{\partial^2 \nabla b}{\partial
x^2}+(-\frac{1}{2}+3d)\frac{\partial^2 \nabla b}{\partial y^2}],
\end{equation}
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 (\ref{eq:backheat})
$$
\omega ~=~\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 \cite{SB72,Lo91},
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
\cite{SB72} 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. 

\subsection{The Generation of Noise During Initial Channelization}
\label{sec:shocks}

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 {\em 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 \cite{SBM97}),
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 \ref{fig:sections}.

\begin{figure}[tbp]
\centerline{\psfig{figure=Figures/sections.ps,height=6.5in}}
\vspace{-7pt}
\caption{Transverse Sections at 10m separation on the Eroding Ridges
         at 10\% and 60\% Eroded}
\vspace{-7pt}
\label{fig:sections}
\end{figure}
\vspace{0.1in}

We base our analysis of the mechanism that provides a
significant source of white noise
on these convex sections of ridge by
linearizing equations (\ref{eq:m6}) and (\ref{eq:m7}) 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
\begin{equation}
\vspace{.1in}
\label{eq:39}
 {\eta}^2 \frac{\partial v }{\partial t} ~=~
 \nabla \cdot \left[\frac{5}{3}h_1^\frac{2}{3}{{\nabla H_1}\over{|\nabla 
 H_1|^{1/2}}}v\right] ~+~ \nabla \cdot \left[h_1^\frac{5}{3}{{\nabla
u}\over{|\nabla 
 H_1|^{1/2}}}-\frac{1}{2}h_1^\frac{5}{3}(\nabla H_1 \cdot \nabla u)
{{\nabla H_1}\over{|\nabla H_1|^{3/2}}}\right],
\end{equation}

\begin{eqnarray}
\vspace{.1in}
\label{eq:40}
\frac{\partial u}{\partial t}
&=&\nabla \cdot \left[{{h_1^\frac{5}{3}}\over{|\nabla
H_1|^{1/2}}}({1\over{\eta}}
+h_1^\frac{5}{3}|\nabla H_1|^{5/2})\nabla u\right]
~+~\nabla \cdot \left[\frac{5}{3}h_1^\frac{2}{3}{{\nabla H_1}\over{|\nabla 
 H_1|^{1/2}}}({1\over{\eta}}
+2h_1^\frac{2}{3}|\nabla H_1|^{5/2})v\right]\nonumber\\[.1in]
&+&\nabla \cdot \left[h_1^\frac{5}{3}({-1\over{2\eta}}
+2h_1^\frac{5}{3}|\nabla H_1|^{5/2})(\nabla H_1 \cdot \nabla u)
{{\nabla H_1}\over{|\nabla H_1|^{3/2}}}\right],
\vspace{.1in}
\end{eqnarray}
The first equation is a hyperbolic equation for $v$ driven by $\nabla u$.
The second equation is a parabolic equation for $u$ driven by $v$ and 
$\nabla 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 $\bf u$ denote the unit vector in the direction
of the gradient of the water surface,
we can write the first equation in the form
\begin{equation}
\vspace{.1in}
\label{eq:41}
 {\eta}^2 \frac{\partial v }{\partial t} ~=~\frac{5}{3}h_1^\frac{2}{3}|\nabla 
 H_1|^{1/2} {\bf u} \cdot \nabla v ~+~ g(x,y,t) v ~+~ f(x,y,t).
\vspace{.1in}
\end{equation}
If we then let the scalar $s$ parametrize the direction of the gradient,
we may rewrite this equation in the form
\begin{equation}
\vspace{.1in}
\label{eq:42}
 {\eta}^2 \frac{\partial v }{\partial t} ~=~\frac{5}{3}h_1^\frac{2}{3}|\nabla 
 H_1|^{1/2} \frac{\partial v }{\partial s}  ~+~ g(x,y,t) v ~+~ f(x,y,t).
\vspace{.1in}
\end{equation}
The last equation is analyzed in Appendicies \ref{sec:shockequation} and 
\ref{sec:sr}, 
where it is shown that 
{\em 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 (\ref{eq:40}) 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 (\ref{eq:39}) 
and generated by the sediment divergences in (\ref{eq:40}) 
that amplify the instabilities produced by (\ref{eq:39}).
These instabilities do not lead to exponential growth because
the water changes on a faster scale than the sediment, see Section
\ref{sec:model}, but instead give rise to white noise
in the sediment flow. 
We obtain, as a result, a heat equation driven by white noise,
\begin{equation}
\vspace{.1in}
\label{eq:17}
\frac{\partial (\nabla u)}{\partial t}
~=~ \nu \gD (\nabla u) ~+~ {\bf r}(x,y,t),
\end{equation}
where $\nu ={{d^\frac{5}{3}}\over{\eta 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 (\ref{eq:39})
and (\ref{eq:40}) to (\ref{eq:17}) 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 {\em 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 (\ref{eq:17}) is well-known to be Brownian 
motion
\begin{equation}
\vspace{.1in}
\label{eq:18}
\nabla u =  \int_0^t \int_{\Bbb R^2} {exp({-|{\bf x}-{\bf 
y}|^2/4\nu(t-\tau))}\over {(4\pi (t-\tau))}} \, d{\bf W}({\bf y},\tau),
\end{equation}
where ${\bf x} = (x,y)$, ${\bf y} = (x',y')$ and $dy = dx'dy'$ is the 
volume element in $\Bbb R^2$ and formally the white noise can be 
understood as $d{\bf W}({\bf y},\tau)= {\bf r}({\bf y},\tau) dy d\tau$. 
Adding this solution to $\nabla H_1$ we obtain {\em a random walk
$$
\nabla H = \nabla H_1 + \epsilon \nabla 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 (\ref{eq:m6}).

\subsection{The Scaling of the Channelization Process}
\label{sec:chan}

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 (\ref{eq:17})
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 \cite{Sn98}). 
The values of the exponents of the associated width functions are
$$
\gb = \frac{1}{4},~~\chi =\frac{1}{2},~~z = 2.
$$
We note that $\gb$ is different from the Hurst exponent,
which has the value $1/2$ for a pure random walk,
because $\gb$ 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
$\nabla H(x+y) - \nabla H(y)$ vanishes
do another random walk on top of that.
Hence {\em 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 \cite{S67,TTPH91}. 
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-\chi = 5/2$. The avalanche dimension is the same $D=2+\chi= 5/2$.
The temporal roughness coefficient $\beta$
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 \cite{SBM97, SMB97}. 
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. 

\section{Self-Similar Landscapes and the Second Fundamental Scaling of Fluvial
Landscapes}
\label{sec:selfsim}

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 \ref{sec:results}
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 \ref{fig:sections}, and as explained in 
Appendix \ref{sec:sr}. 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 \ref{fig:sections}, 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 \cite{Le93}
and \cite{Sn98}. 

\subsection{The Emergence of Self-Similar Landscapes}
\label{sec:prevss}

Smith, Birnir, and Merchant \cite{SBM97}
showed that the model (\ref{eq:m3})-(\ref{eq:m4})
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
{\em 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 (\ref{eq:m3}), (\ref{eq:m4}),
Smith, Merchant, and Birnir \cite{SMB97}
showed that after a characteristic erosion time
\footnote{The characteristic time is 45\% of the surface eroded
for planar initial surfaces with $\gamma=\delta=2$.}
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 \ref{sec:numer},
this stage of development is characterized by 
a dominant scaling of the width function with exponent $0.7$.
We now analyze the model (\ref{eq:m6}), (\ref{eq:m7})
in order to determine the origin and physical
basis for this scaling.

\subsection{The Origin of the Second Diffusive Scaling}
\label{sec:second}

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 \ref{sec:channels}.
To find and analyze this scaling we linearize equations (\ref{eq:m6})
and (\ref{eq:m7}) about the separable solutions (\ref{eq:sep}) 
representing the mature landscapes of valleys and ridges
\footnote{As in our previous analyses, we employ sediment transport
parameter values $\gg = \gd = 2$, noting
that similar results hold for
parameter values in a neighborhood of these.}
$$
H = H_2(x,y,t) + \epsilon u(x,y,t),~~h = h_2(x,y.t) + \epsilon 
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 (\ref{eq:m6}, \ref{eq:m7}).
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 (\ref{eq:39}) and (\ref{eq:40}),
but with the terms $H_1(x,y,t),~~h_1(x,y,t)$
replaced by the separable solutions
$H_2,~~h_2$.

Again, the first equation is a hyperbolic equation for 
$v$ driven by $\nabla u$ and
the second is a parabolic equation for $u$ driven by $v$ and 
$\nabla v$. Since the first equation 
is really a hyperbolic PDE in one space 
dimension, exactly as in Section \ref{sec:shocks}, 
we let the scalar $s$ parametrize the direction of the gradient
to obtain
\begin{equation}
\vspace{.1in}
\label{eq:22}
 {\eta}^2 \frac{\partial v }{\partial t} ~=~\frac{5}{3}h_2^\frac{5}{3}|\nabla 
 H_2|^{1/2} \frac{\partial v }{\partial s}  ~+~ g(x,y,t) v ~+~ f(u,t).
\vspace{.1in}
\end{equation}
A straight-forward analysis of this equation,
presented in Appendix \ref{sec:sr}, 
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 \ref{sec:channels}
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 (\ref{eq:39}), 
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 {\em quenched noise} or 
the Langevin 
equation 
\begin{equation}
\label{eq:ren}
\frac{\partial u }{\partial t} ~=~  \Delta u + F + r(u,t),
\end{equation}
where $F$ is the driving force, $r$ is the quenched noise, and $\gl$ is a
parameter.
Moreover, it is very different from equation (\ref{eq:17}) 
since the noise term
$$
\frac{5}{3}h_2^\frac{2}{3}{{\nabla H_2}\over{|\nabla 
 H_2|^{1/2}}}({1\over{\eta}}
+2h_2^\frac{2}{3}|\nabla H_2|^{5/2})\cdot \nabla 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 (\ref{eq:39}) plays the role of the forcing.
It is known, see \cite{Sn98}, \cite{Le93}, 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)+ \epsilon u(x,y,t)
$$
executes a random motion driven by noise produced by the divergence in the 
sediment flow, seeded by shocks
$\nabla h_o(x,y) T^{-\frac{3}{10}}(t) + \epsilon
\nabla v(x,y,t)$ propagating down the separable surface on a faster time-scale,
the surface being pinned 
at randomly distributed sites.  The slope $\nabla H$ evolves in the same
fashion.

\subsection{The Scaling of Mature Landscapes}
\label{sec:mature}

We compute the exponents using the formulas in Section \ref{sec:struct}, 
starting with the numerical value of the spatial and temporal roughness 
exponents from Section \ref{sec:numer}.  The values of the exponents are 
$$
\gb = 0.4,~~\chi = 0.7,~~z = 1.8.
$$
The fractal dimension becomes $D_f = 3-\chi = 2.3$.
The avalanche dimension is $D = 2+\chi = 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  \cite{Le93} and \cite{Sn98}, and described by equation (\ref{eq:ren}),
$$
\gb = 0.475,~~\chi = 0.75,~~z = 1.58.
$$
The computed spatial roughness coefficient $\chi$ 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 \cite{Ca96}. 
The latter models also give rise to values of the scaling exponent
$\gb$ similar to ours. Thus two remarkable 
features of the continuous model (\ref{eq:m6} ,\ref{eq:m7}) are apparent.
First, it generates the noise from the water 
flowing down the separable surface (\ref{eq:sep}) 
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
\footnote{As noted in \cite{Sn98}:
\begin{quote} 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.  \end{quote}
}.
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. 

 

\section{Discussion and Significance of the Results}
\label{sec:disc}

The preceding results on the origins of scaling
in solutions to the models (\ref{eq:m6})-(\ref{eq:m7})
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.

\subsection{The Two Processes Forming the Landscape}
\label{sec:process}
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 
{\em 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 $\chi = 0.5$.
The system then gradually makes a transition to a second process that 
we call the {\em 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 $\chi = 0.75$.
\footnote{We believe this to be the true value of $\chi$ althought the 
transcient numerically determined value of $0.7$ is somewhat lower at
$60$\% eroded.}


\begin{figure}[tb]
\centerline{\psfig{figure=Figures/basin.ps,height=3in}}
\vspace{-7pt}
\caption{The basin of the Amazon river}
\vspace{-7pt}
\label{fig:basin}
\end{figure}


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 \ref{fig:basin}.
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~+~\chi$. The exponent 
one measures 
the size of basin $l$ along the main river and the exponent $\chi$ measures 
the size of its
bulge $l^\chi$ (width), in the perpendicular direction, 
see Figure \ref{fig:basin}, where $l$ is the length of the main river.


\subsection{The Origin of Observed Landscape Width Scalings}
\label{sec:emp}

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 \cite{WP94},
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 {\em 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.  

\subsection{ The Origin of Hack's Law:~}
\label{sec:Hack}

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 \ref{sec:mature} 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 \sim l^D
$$
and the (planar) avalanche dimensions is $D = 1+\chi$. 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^\chi$.
Stable scalings for the surface emerge together with
the emergence of the separable solutions.
We note that in this case $\chi=0.75$, hence we obtain
\begin{eqnarray}
\label{hack}
l &\sim& A^{1\over{1+\chi}}\\
  & \approx& A^{0.57} 
\end{eqnarray}
a number that is in excellent agreement with observed
values of the exponent of Hack's law of $0.568$, see \cite{Gr61}. 
We note, in particular, that $\chi$ 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 \cite{Sn98}) 
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 \cite{RR97},
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. 


\subsection{The Scaling of Exceedence Probabilities}
\label{sec:emp1}

The probability $P[A \geq 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
\begin{eqnarray}
P[A \geq a] \approx a^\rho,
\end{eqnarray}
which is a good approximation to observed variation,
are close to $-0.4$ (see, for example, \cite{RR97}).

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 $\beta_{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 $\chi$ as above. Then
$$
P[A \geq a] = P[l \geq t]
$$
since by Hack's Law 
$A = l^{1+\chi}$ and $a = t^{1+\chi}$, $t$ being the length
of the main branch of the subbasin having area $a$.
However
$$
\int_0^t P[l \geq s] ds = \int_0^t s^{-\beta_{si}}ds \sim t^{1-\beta_{si}},
$$
and 
$$
1-\beta_{si} = {-{\chi}
},
$$
see \cite{Sn98}.
Thus using Hack's law again
we obtain
\begin{eqnarray}
P[A \geq a] &\approx& t^{-{\chi}}\\
            &\approx& a^{-\frac{\chi}{\chi+1}}\\
            &\approx& a^{-0.4},
\end{eqnarray}
for $\chi = 0.75$. Notice that the exponent is the negative of the 
dynamic roughness exponent $\beta = 0.4$ in Section \ref{sec:mature}.

\subsection{A Bridge Between Deterministic and Stochastic Modeling}
\label{sec:hor}

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
(\ref{eq:m6}) and (\ref{eq:m7}) and 
observed in the studies of Smith, Merchant, and Birnir \cite{SMB97}
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 \ref{sec:intro}.
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, \cite{Sc67,TTPH91}) and
are capable of generating branching networks that obey
Horton's laws (see \cite{RR97} 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.

\subsection{A Model of the Stages of Drainage Basin Evolution}
\label{sec:stages}

Based on the preceding results and the
findings of our two previous papers \cite{SBM97,SMB97},
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
\footnote{Clearly, the three stages
do not have well-defined beginnings and endings,
but they do mark significant differences
in the qualitative dynamics of the system.}.
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 {\em youthful} stage of W.M.Davis' grand model
of landscape evolution \cite{Da02}.

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 {\em transient attractor} (see, for example, \cite{SBM97})
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
{\em 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 {\em global attractor} $H~=~0$ exerts its influence.
This stage may be viewed as the final, or {\em 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.

\section{Summary}
\label{sec:summary}

The simple advective model of fluvial erosion 
({\ref{eq:m6})-({\ref{eq:m7}) 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 ({\ref{eq:m6})-({\ref{eq:m7}) 
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 ({\ref{eq:m7}). 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 (\ref{eq:ratios}) 
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 \cite{Le93} and
\cite{Sn98},
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^\chi$ of the basin, 
where $\chi = 0.75$ is the spatial roughness coefficient
of the maturation process in Section \ref{sec:mature}. 
This accounts for the oval shape of river basins
(see \cite{RR97}), 
and leads to a derivation of Hack's law (see Section \ref{sec:Hack})
as well as the observed distribution of the 
exceedence  probabilities, see Section \ref{sec:emp1}.

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 \cite{LD93}. 
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 \ref{sec:noise} and Appendix \ref{sec:shockequation} 
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.
 
\section*{Acknowledgements}
The authors wish to thank Kim Sneppen, 
Daniel Lavall\'ee, 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.

\newpage

\appendix
\section{The Generation of Shocks in the Water Flow Equation}
\label{sec:shockequation}

Equation (\ref{eq:m6}) describes the water flow down a water surface
represented by H. 
We rewrite this equation in the form
\begin{eqnarray}
\vspace{.1in}
\label{eq:s1}
 {\eta^2} \frac{\partial h }{\partial t}~-~
\!\nabla \cdot \left[ {\bf u}
h^\frac{5}{3}~{\!\mid \!\nabla H \!\mid }^\frac{1}{2} \right]~=~R,
\vspace{.1in}
\end{eqnarray}
where ${\bf u} = \frac{\!\nabla H }{\!\mid \!\nabla H \!\mid}$ 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
\begin{eqnarray}
\vspace{.1in}
\label{eq:s2}
  \frac{\partial h }{\partial t}~-~
\frac{5}{3\eta^2}h^{2/3}\ {\!\mid \!\nabla H \!\mid }^\frac{1}{2} {\bf 
u}\cdot\nabla 
h ~=~Q{{h^{5/3}}\over{\eta^2}}~+~\frac{R}{\eta^2},
\vspace{.1in}
\end{eqnarray}
where
$$
Q ~=~ \left[{{{\!\mid \!\nabla H \!\mid }^{2}\Delta H-\frac{1}{2}
(H_x^2H_{xx}+2H_x H_y H_{xy}+H_y^2H_{yy})}
\over{{\!\mid \!\nabla H \!\mid }^\frac{5}{2}}}
\right]
$$
and using that ${\bf u}\cdot\nabla h ~=~  \frac{\partial h }{\partial s}$, where
$s$ parametrizes the downstream direction given by ${\bf u}$, we obtain a
a one-dimensional PDE
\begin{eqnarray}
\vspace{.1in}
\label{eq:s3}
  \frac{\partial h }{\partial t}~-~
\frac{5}{3\eta^2}h^{2/3}\ {\!\mid \!\nabla H \!\mid }^\frac{1}{2}\frac{\partial h}
{\partial s}  ~=~Q{{h^{5/3}}\over{\eta^2}}~+~\frac{R}{\eta^2}.
\vspace{.1in}
\end{eqnarray}
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 (\ref{eq:s3}),
\begin{eqnarray}
\label{eq:s5}
\vspace{.1in} 
\frac{dt}{d\tau} = 1, ~~~\frac{ds}{d\tau} =-\frac{5}{3\eta^2}h^{2/3}\
{\!\mid \!\nabla H \!\mid }^\frac{1}{2}\nonumber\\
\\
\frac{dh}{d\tau} = Q{{h^{5/3}}\over{\eta^2}} = -p(s,\tau) h^{5/3}.\nonumber
\vspace{.1in}
\end{eqnarray}
The function $p$ is positive
$$
p(s,\tau) = -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
$\nabla H \neq 0$.  
Now we solve the $h$ equation
in (\ref{eq:s5}) and obtain
\begin{equation}
\label{eq:s6}
h(s,t) = {h_o(s)\over{(1+\frac{2}{3}h_o^{2/3}(s) P(s,t))^{3/2}}}
\end{equation}
where
$$
P(s,t) = \int_0^t p(s(\tau)+s,\tau) d\tau,
$$
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
\begin{equation}
\label{eq:s7}
\frac{ds}{d\tau} ~=~-\frac{5}{3\eta^2}h^{2/3}\
{\!\mid \!\nabla H \!\mid }^\frac{1}{2}~=~ -\frac{5}{3g(s)\eta^2}h^{2/3},
\end{equation}
if $H$ evolves on a slower time-scale than $h$. But this was precisely the 
conclusion of the discussion in Section \ref{sec:family}. Thus we can
separate the variables and solve the above equation
$$
G(s) |_{s_o}^s = \int_{s_o}^s g(s) ds =
 -\frac{1}{\eta^2}\int_0^t h^{2/3}(s(\tau),\tau)\  d\tau. 
$$

Next we expand the solution $h$ in (\ref{eq:s7}) above in powers of $h(s_o)$,
$$
h(s,t) = h(s_o) - h^{5/3}(s_o)\int^t_0 p(s(\tau),\tau)\ d\tau + 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)}{\eta^2} t + O(h^{5/3}(s_o)/\eta^2).
$$
This implies that two characteristics emanating from two distinct points
$s_1 > s_o$ meet at the time
\begin{equation}
\label{eq:s8}
t = \frac{3}{5}{{\eta^2(G(s_1)-G(s_o))}\over{h^{2/3}(s_1)-h^{2/3}(s_o)}} + 
O(h(s_o),h(s_1)).
\end{equation}
Moreover, these characteristics carry different values of $h$ by
equation (\ref{eq:s6}) producing a shock at the location
$$
s = G^{-1}(G(s_o)-\frac{5}{3}\frac{h^{2/3}(s_o)}{\eta^2} t + O(h^{5/3}(s_o)/\eta^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 $|\nabla H|^{1/2}$. For example when $g(s) = |\nabla 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)}{\eta^2} t 
+ O(h^{5/3}(s_o)/\eta^2)
$$
and the shock develops at the time
$$
t = \frac{12}{25}{{\eta^2(s_1^{5/4}-s_o^{5/4})}
\over{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
{\em 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 \ref{sec:sr}.

\newpage


\begin{figure}[tb]
\centerline{\psfig{figure=Figures/char.ps,height=2in}}
\vspace{-7pt}
\caption{The characteristics for the shock caused by surface singularities.}
\vspace{-7pt}
\label{fig:sr}
\end{figure}
\vspace{0.1in}

\section{Surface Singularities and Shocks}
\label{sec:sr}

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 {\em one-dimensional} hyperbolic equation for the evolution
of the variation of the water surface

\begin{equation}
\vspace{.1in}
\label{eq:55}
 {\eta}^2 \frac{\partial v }{\partial t} ~+~\frac{5}{3}\frac{\partial
}{\partial s}(-h_1^\frac{2}{3}|\nabla 
 H_1|^{1/2} v) ~=~ f(x,y,t),
\vspace{.1in}
\end{equation}
where 
$$
{1\over{g(s)}} = -\frac{5}{3}h_1^\frac{2}{3}|\nabla 
 H_1|^{1/2} = -{\eta}^2\frac{\partial q_w}{\partial 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, ${\bf u}_w^\bot$, where
${\bf u}_w= -\frac{\nabla H}{|\nabla H |}$,\footnote{The separable landscape
is of course not uniform in the direction of ${\bf u}^\bot$ but this example 
still gives the right information.}
\begin{equation}
\label{eq:convex1}
g(s) = -s^{-1/4}.
\end{equation}
is a reasonable model. We make the change of variables
$$
a = \frac{5}{3}h_1^\frac{2}{3}|\nabla  H_1|^{1/2}v = -\frac{v}{g(s)}
$$
then the equation can be written in the form
\begin{equation}
\label{eq:sr}
\frac{\partial a }{\partial t} ~+~\frac{1}{g(s)\eta^2}\frac{\partial a
}{\partial s} ~=~ f(x,y)/\eta^2.
\end{equation}	
If we now let
$$
\frac{ds}{d\tau} = \frac{1}{g(s)\eta^2},
$$
define the characteristics and $g = \frac{\partial G
}{\partial s}$, then equation (\ref{eq:sr}) can be solved 
along characteristics,
$$
a(s,t) = a_0(G^{-1}(G(s)-\frac{t}{\eta^2})) + \frac{1}{\eta^2}\int_0^t 
f(G^{-1}(G(s)-\frac{t-\tau}{\eta^2}),\tau)d\tau.
$$
The important equation is the equation for the characteristics 
\begin{equation}
\label{eq:char}
\frac{ds}{d\tau} = \frac{1}{g(s)\eta^2} = -\frac{5}{3 \eta^2}h_1^\frac{2}{3}|\nabla 
H_1|^{1/2}.
\end{equation}

It says that if the derivative of $|\nabla 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}{\eta^2}),
$$
see Figure \ref{fig:sr}. 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 \ref{sec:noise} and Appendix 
\ref{sec:shockequation}.

\newpage

\section{The Scaling of the Gradient}
\label{sec:grad}

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 $|\nabla H| $. 
Consider
\begin{eqnarray*}
|\nabla H | (y+x) &=& [H_x^2(y+x)+H_y^2(y+x)]^{1/2}\\
\sim [(ex^\chi + O(x^{2\chi})+\cdots)^2&+&(c+dx^\chi+O(x^{2\chi}))^2]^{1/2}\\
=(c^2+2cdx^{\chi}&+&O(x^{2\chi}))^{1/2},
\end{eqnarray*}
in a statistical sense
where $c, d, e,$ are constants. 
By the inequality,
$$
(a+b)^{1/2} \le a^{1/2}+b^{1/2},
$$
\begin{eqnarray*}
<(|\nabla H | (y+x)-|\nabla H | (y))^2>^{1/2} 
&=& (c^2+2cdx^\chi+O(x^{2\chi}))^{1/2}-c \\
\le (2cdx^\chi+O(x^{2\chi}))^{1/2} &=&  x^{\chi/2}(2cd+O(x^{\chi}))^{1/2}.
\end{eqnarray*}
Similarly
\begin{eqnarray*}
<(|\nabla H | (y+x)-|\nabla H | (y))^2>^{1/2} 
&=& (c^2+2cdx^\chi+O(x^{2\chi}))^{1/2}-c \\
\ge (2cdx^\chi+O(x^{2\chi}))^{1/2} &=&  -x^{\chi/2}(2cd+O(x^{\chi}))^{1/2}.
\end{eqnarray*}
by the inequality
$$
(a+b)^{1/2} \ge a^{1/2}-b^{1/2}.
$$
This shows that the statistical scaling of $|\nabla H|$ has the exponent
$\chi/2$, whenever $\chi$ is the spatial roughness exponent. 
	
\newpage

\section{River Basins and Horton's Laws}
\label{sec:Horton}

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
$$
\gb = \frac{1}{2},~~\chi =1,~~z = 2,
$$
where $\gb$ 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-\chi =1 $. The corresponding (planar) avalanche dimension is  
$D=1+\chi= 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 \ref{sec:struct}.)

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
\footnote{
1. 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.} 
Horton \cite{ho45}
discovered that if $N(a)$ denotes the number of streams of order $a$ and
$L(a)$ is their mean length then the ratios
\begin{equation}
\label{eq:Horton}
{{L(a+1)}\over{L(a)}}= R_L,~~{{N(a)}\over{N(a+1)}}= R_B,
\end{equation}
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 \cite{BaRo89}, this was related to the avalanche (embedding) 
dimension of a river network
$$
D = {{ln(R_B)}\over{ln(R_L)}}~=~2,
$$
by the value of $D$ computed above. Thus we obtain
\begin{equation}
\label{eq:ratios}
R_B = R_L^2.
\end{equation}
This relationship is borne out by observations
where typically $ R_B \sim 4$ and $R_L \sim 2$. 
%\bibliographystyle{apalike}
%\bibliography{ref4}
\input{sca17.bbl}
%\bibliography{ref}
\setcounter{figure}{0}
\begin{figure}[tb]
\centerline{\psfig{figure=Figures/tenpercentH1.ps,height=9in}}
\vspace{-7pt}
\caption{The Water surface $H$ of a Typical Run at 10\% of the Surface Eroded}
\vspace{-7pt}
\label{fig:tenperH}
\end{figure}


\begin{figure}[tb]
\centerline{\psfig{figure=Figures/tenpercent1.ps,height=9in}}
\vspace{-7pt}
\caption{The Water Depth $h$ and the Slope of the Water Surface
$\mid \nabla H \mid$ of a Typical Run at 10\% of the Surface Eroded}
\vspace{-7pt}
\label{fig:tenper}
\end{figure}

\begin{figure}[tb]
\centerline{\psfig{figure=Figures/sixtypercentH1.ps,height=9in}}
\vspace{-7pt}
\caption{The Water surface $H$ of a Typical Run at 60\% of the Surface Eroded}
\vspace{-7pt}
\label{fig:sixtyperH}
\end{figure}


\begin{figure}[tb]
\centerline{\psfig{figure=Figures/sixtypercent1.ps,height=9in}}
\vspace{-7pt}
\caption{The Water Depth $h$ and the Slope of the Water Surface
$\mid \nabla H \mid$ of a Typical Run at 60\% of the Surface Eroded}
\vspace{-7pt}
\label{fig:sixtyper}
\end{figure}
\end{document}







