\pretolerance=1400
\magnification=\magstephalf
\baselineskip=0.50cm
\nopagenumbers
\font\bigbrm=cmbx16
\font\bigrm=cmr14
\font\medrm=cmr12
\overfullrule=0pt
%\font\bigrm=cmr10 scaled\magstep2
%\font\medrm=cmr10 scaled\magstep1
%\font\rmnor=cmr12
%\font\twbf=cmbx12
%\font\bigbf=cmbx12 scaled\magstephalf
\def\ss{\smallskip\noindent}
\def\ms{\medskip}
\def\bs{\bigskip}
\def\bsn{\bigskip\noindent}
\def\hah{\hfill&}
\def\hac{\hfill &\cr}
\def\ni{\noindent}
\def\itb{\item{$\bullet$}}
\def\itc{\item{$\circ$}}
\def\bull{$\bullet\;\;$}

\ni
{\bigbrm  Multiscale Green's function Monte Carlo approach
to erosion modeling and its application to land use optimization.}

\bsn
 Lubos Mitas and Helena Mitasova

\ni
University of Illinois at Urbana-Champaign,
Urbana, Illinois, USA

%National Center for Supercomputing Applications,
%405 N. Mathews Ave,
%
% Geographic Modeling Systems Laboratory, Department of Geography,
%220 Davenport Hall, University of Illinois at Urbana-Champaign,
%Urbana, Illinois, USA

\bsn
\bsn
\bsn
\bsn
\bsn
\bsn
%{\bf Abstract }
%New generation of simulation tools for modeling soil erosion, sediment
%transport and deposition by overland water flow in complex landscapes 
%is presented. The simulations are based on the solution
%of bivariate continuity equations describing
%water and sediment transport over 3D terrain with variable climatic, 
%soil and land cover conditions. Underlying equations are solved by 
%advanced computational approaches based on stochastic projection techniques 
%(Green's function Monte Carlo), providing the robustness and flexibility
%necessary for complex conditions and multiscale implementation. 
%%with spatially variable accuracy/resolution.
%The processing, analysis and visualization of data and results are performed 
%using an emerging multidimensional dynamic GIS technology.
%The possibilities of using the simulations as a tool for finding optimal
%land use patterns with minimized erosion risk are evaluated by comparing
%computer optimized land use scenario with the traditional
%landuse design.

\bsn
\bsn
{\bigrm
1 Introduction
}
\bsn
Efforts to balance economic development with environmental protection
have increased the demand for simulation tools which enable predictions
of the human
impact on the landscape. In order to prevent irreversible changes and avoid
costly, ineffective solutions, the simulation tools should provide
detailed spatial and temporal distributions of modeled phenomena.
Statistical averages for entire study areas or predictions only for a
certain point, such as a watershed outlet, are often insufficient.
Effectiveness of land management decisions aimed at preventing negative
impacts of soil erosion in complex landscapes can be significantly
improved by detailed predictions of erosion and deposition patterns
for proposed land use alternatives.
Recent advances in Geographic Information Systems (GIS)
technology, especially support for modeling
with multivariate functions (Mitasova et al., 1995; Mitas et al., 1997),
along with the exponential growth in computational power,
stimulate the shift from empirical, lumped models
to physically-based, distributed ones
(Moore {\sl et al.} 1993, Maidment 1996, Saghafian {\sl et al.} 1995).

In spite of a significant progress in this area of research,
applications of distributed models
are still rather laborious and often results do not have sufficient detail,
accuracy and realism for land management purposes. 
In our research we have tried to address some of these
problems by focusing on: 

\item{(a) } description of processes by 'first principles' relations in a 
 bivariate form allowing us to incorporate impact of spatial variability in
 rainfall excess, terrain, soil and cover conditions;

\item{(b) } use of robust solvers (Green's function Monte Carlo)  
 with multiscale implementation which 
minimize the preprocessing of input data and support modeling
at spatially variable resolution;

\item{(c) }investigations of computer simulated land use designs for 
finding optimal land use patterns with minimized erosion risk;

\item{(d) } increase in efficiency by extensive use of GIS for processing, 
analysis and visualization of the data and results, 
as described in detail by
Mitasova {\sl et al.} (1995); Mitas {\sl et al.} (1997).

First, we briefly describe water and sediment flow 
 models based on the solution of bivariate first principles 
equations by multiscale Monte Carlo method.
Then, we illustrate the capabilities of the proposed approach
to simulate erosion/deposition patterns for a study area with 
spatially variable terrain, 
soil and cover conditions for different land use designs
aimed at improving
the effectiveness of erosion prevention.
Finally, we present an example of water flow simulation 
with spatially variable resolution using heterogeneous
elevation data.

 
\bsn
\bsn
{\bigrm
2 Methods
}
\bsn
The methodological framework for the simulation of human
impact on erosion/deposition is based upon the
description of water flow and sediment transport processes
by first principles equations, a concept outlined previously,
most often for a one dimensional case, for example,
by Bennet, (1974).
Within our approach, inputs and outputs of simulations
are represented by continuous multivariate functions 
discretized as grids, 
 as opposed to homogeneous hillslope segments or subwatersheds
used in more traditional approaches.
To fully incorporate the impact of spatial variability in
terrain, soils and cover, we describe the water and sediment
flow as bivariate vector fields rather than using
a 1D flow routing in the steepest slope direction
from a planar hillslope segment common in
many distributed models (e.g., Moore et al., 1993; 
Flacke et al., 1990).
Advanced GIS technology
is used to support the processing, analysis and visualization
of the multiscale data and simulation results (Mitasova et al., 1995; 
Mitas et al., 1997).

\bsn
\bsn
{\bf
2.1 Overland water flow
}
\bsn
The model used in this paper is described by Mitas and Mitasova (1998)
therefore here we briefly present only its principles.
A bivariate shallow water flow continuity equation for a rainfall event 
(e.g., Saghafian {\sl et al.} 1995) is given by

$$
{\partial h({\bf r},t) \over \partial t} =
 i_e({\bf r},t) - \nabla \cdot [h({\bf r},t){\bf v}({\bf r},t)]
\eqno (1)
$$

\ni
where $h({\bf r},t)\;[m]$ is the depth of overland flow,
$i_e({\bf r},t)\;[m/s]$ is the rainfall excess, 
{\bf v}({\bf r},t) is velocity and
${\bf r}=(x,y)\; [m]$ is the position. The Eq. (1)
 is coupled with the momentum conservation equation
 (in the diffusive wave approximation) which together with
 the Manning's relation between the depth and velocity 
 create a closed system of equations.
 The steady state form of Eq. (1) is given by

$$
-{\varepsilon\over 2 }\nabla^2 [h^{5/3}({\bf r})]
+\nabla \cdot {\bf q}({\bf r}) = i_e({\bf r})
\eqno (2)
$$

\ni
where ${\bf q}({\bf r},t)\;[m^2/s]$ is the unit flow discharge
(water flow per unit width) and $\varepsilon$ is a 
diffusion coefficient.
The diffusive wave effects are incorporated approximately by
the term $ \propto \nabla^2 [h^{5/3}({\bf r})]$.
Such an
 incorporation of diffusion in the water flow simulation is not new
 and a similar term
 has been obtained in derivations of diffusion-advection equations for overland
 flow, e.g.,  by Dingman (1984).
 In our reformulation we use
 a modified diffusion term, which depends on $h^{5/3}({\bf r})$ instead of
 $h({\bf r})$. Eq. (2) has the advantage of being {\sl linear}
 in the function $h^{5/3}({\bf r})$ 
 which enables us to solve it by means of the Green's function method
 using stochastic (Monte Carlo) techniques as described later.

 The diffusion constant which we have used is rather small 
 (approximately one order of
 magnitude smaller than the reciprocal Manning's coefficient) and
 therefore the resulting flow is close to the kinematic regime. However,
 the  diffusion term 
improves the kinematic solution by filling and overcoming shallow depressions
common in digital elevation models (DEM),
reducing thus the need for manual modification of a DEM.
It also produces smooth flows around slope discontinuities
or abrupt changes in land cover which are
typical for anthropogenic landscapes. 


\bsn
\bsn
{\bf
2.2 Sediment transport by overland flow
}
\bsn
Overland water flow is the driving force for hillslope erosion
 which includes
sediment entrainment, transport and deposition.
 The continuity of sediment mass is given by

$$
{\partial [\rho_sc({\bf r},t)h({\bf r},t)] \over \partial t} 
-{\omega \over 2} \nabla^2 [\rho_sc({\bf r},t)h({\bf r},t)]+
\nabla\cdot {\bf q}_s({\bf r},t) = {\rm sources - sinks}=
D({\bf r},t)
\eqno(3)
$$

\ni
where
${\bf q}_s({\bf r},t) = \rho_sc({\bf r},t){\bf q}({\bf r},t)
\; [kg/(ms)]$ is the sediment flow rate per unit width,
$c({\bf r},t)\;[{\rm particle}/m^3]$ is sediment concentration,
$\rho_s\;[kg/{\rm particle}]$ is mass per sediment particle,
and $\omega [m^2/s] $ characterizes local dissipation (diffusion) processes. 
For a steady state case the equation is 

$$
-{\omega \over 2} \nabla^2 [\rho_sc({\bf r})h({\bf r})]+
\nabla\cdot {\bf q}_s({\bf r}) = D({\bf r}) .
 \eqno(4)
$$

\ni
The sources/sinks term is derived from the assumption that
the detachment and deposition rates are proportional to the difference between
the sediment transport capacity and the actual
sediment flow rate (Foster and Meyer, 1972):

$$
D({\bf r})
=\sigma({\bf r}) \bigl[ T({\bf r})- |{\bf q_s}({\bf r})|\bigr]
\eqno(5)
$$

\noindent
where
$T({\bf r})\;[kg/(ms)]$ is the sediment transport capacity,
and
$\sigma({\bf r})\;[m^{-1}]$ is the first order reaction term
 dependent on soil and cover properties.
 The expression for $\sigma({\bf r})=D_c({\bf r})/T({\bf r})$
is obtained from the following relationship (Foster and Meyer, 1972):

$$
D({\bf r})/D_c({\bf r}) + |{\bf q_s}({\bf r})|/T({\bf r}) = 1
\eqno(6)
$$

\ni
 which states that the ratio of erosion rate to the detachment
 capacity $D_c ({\bf r})\;[kg/(m^2s)]$
 plus the ratio of the sediment flow to the sediment transport capacity is
 a conserved quantity (unity).
Eq. (6), proposed by Foster and Meyer (1972),
is based on observed relationship between soil detachment and transport
described e.g., by Meyer and Wischmeier
(1969). This concept is used in several erosion models including WEPP
(Haan et al., 1994; Flanagan and Nearing, 1995).
The qualitative arguments, experimental observations and values for
 $\sigma({\bf r})$   are discussed by Foster and Meyer, (1972).

The sediment transport capacity $T({\bf r})$ and detachment capacity
$D_c({\bf r})$ represent maximum potential
sediment flow rate and maximum potential detachment rate, respectively,
%by overland flow
and are functions of a shear stress (Foster and Meyer, 1972):

$$
T({\bf r})=K_t({\bf r}) \bigl[\tau({\bf r})\bigr]^p
=K_t({\bf r}) \bigl[\rho_w\, g h({\bf r}) \sin \beta({\bf r})
\bigr]^p
\eqno(7)
$$

$$
D_c({\bf r})=K_d({\bf r}) \bigl[\tau({\bf r}) - \tau_{cr}({\bf r})\bigr]^q
=K_d({\bf r}) \bigl[\rho_w\, g h({\bf r}) \sin \beta({\bf r})
-\tau_{cr}({\bf r})\bigr]^q
\eqno (8)
$$

\noindent
where $\tau({\bf r})=\rho_w\, g h({\bf r}) \sin \beta({\bf r})
\; [Pa]$ is the shear stress,
$ \beta \; [deg] $ is the slope angle,
$p$ and $q$ are exponents,
$ K_t({\bf r}) \; [s]$ is the effective transport capacity coefficient,
%($[s]$ for $p=1$),
$K_d({\bf r}) \; [s/m]$ is the effective erodibility
(detachment capacity coefficient),
%($[s/m]$ for $p=1$),
$\rho_w\, g$  is the hydrostatic pressure
of water with the unit height,
$ g=9.81 \;[m/s^2]$ is the gravitational acceleration,
$\rho_w=10^3 \;[kg/m^3] $ is the  mass density of water,
and
$\tau_{cr}({\bf r})\;[Pa]$ is the critical shear stress.
The parameters and adjustment factors for the estimation of
$D_c ({\bf r}), \,  T({\bf r})$
are functions of soil and cover properties, and their values
for a wide range of soils, cover, agricultural and erosion prevention practices
are being developed within the WEPP model.

Recently, Nearing et al. (1997) presented an improved fit to several sets
of experimental data by relating sediment load
$q_s [g/ms]$ to stream power $\omega [g/s^3]$:
$$
\log q_s = A +{Be^{(C+D \log\omega)}
            \over {1+e^{(C+D \log\omega)}}}
\eqno(9)
$$
 with the constants $A=-34.47, \, B=38.61, \, C=0.845, \, D=0.412$.
Based on the conditions of experiments
it was suggested that Eq. (9) could be a reasonable estimation
of the sediment transport capacity. The Eq. (9)
 can be rewritten to the following form:
$$
T({\bf r}) \approx |{\bf q_s}({\bf r})|=
 a_0 \exp
 \left[-{b \over {1 + [\omega({\bf r})/ \omega_0]^{d}}}\right]
\eqno(10)
$$
where $a_0=1380\, [kg/s m],\,
b=88.90, \, d=0.179, \,
\omega_0=8.89.10^{-6} \, [kg/s^3], \,$
$ \omega=\rho_w\, g |{\bf q}({\bf r})| \sin \beta({\bf r})$
is the stream power $[kg/s^3]$, (note the
 change of units from $gramm$ to $kg$ for  $q_s$ and $\omega$ 
when compared with Eq. (9)).
This form of the Eq. (9) allows us to define a physical interpretation
of the constants, as $a_0$ represents a saturated sediment load for
infinitely large stream power,
$\omega_0$ is a "reference" stream power,  $b=88.90$ and $d=0.179$
are dimensionless
 exponents. Strictly speaking, the choice of the constants corresponds
 to the experimental results used in the fit and could be different
in other cases, e.g., an effective transport capacity coefficient
 analogous to the one in Eq. (7)
 has to be included for different covers, etc.
 The important difference between the Eqs. (7) and (9, 10) is that
through the stream power
$\omega ({\bf r}) =\tau({\bf r}) |{\bf v}({\bf r})|$
 the effect of flow velocity is directly incorporated
 into the  transport capacity.
For complex terrain and cover conditions, the flow velocity varies
and can change dramatically with varying location, so we can
 expect differences in predicted patterns of erosion/deposition
 when using Eq. (10) for the sediment transport capacity.

\bsn
\bsn
{\bf
2.3 Green's function Monte Carlo method 
}
\bsn
The continuity Eqs. (1) and (3) are traditionally solved by finite 
element or finite difference methods. 
These methods often require special data structures (e.g. meshes for
finite element methods) or have problems with numerical stability for complex,
spatially variable conditions.
As a robust and flexible alternative to these methods we have proposed to use
a stochastic approach to the solution, based on Green's function Monte Carlo method.
 Because of space constraints, we consider only the steady state 
 cases given by Eqs. (2) and (4), however, 
methodology is similar also for the time dependent cases. 
 The Eqs. (2), (4) have a general form in which
 a linear differential operator $\cal O$  acts on a nonnegative function
 $\gamma({\bf r})$ (either $h({\bf r})$ or 
 $\rho_sc({\bf r})h({\bf r})$),
  while on the right hand side, there is a source term ${\cal S}({\bf r})$
$$
{\cal O} \gamma({\bf r})= {\cal S}({\bf r})
\eqno(11)
$$
Using the Green's function $G({\bf r},{\bf r'},p)$ the solution can be
 expressed as
$$
\gamma({\bf r})= \int_0^{\infty}\int G({\bf r},{\bf r'},p)
 {\cal S}({\bf r'}) d{\bf r'}dp
\eqno(12)
$$
while 
$G({\bf r},{\bf r'},p)$
is given by the following equation and an initial condition
$$
{\partial G({\bf r},{\bf r'},p)
\over \partial p} =-{\cal O} G({\bf r},{\bf r'},p)\, ;
\;\; G({\bf r},{\bf r'},0) =\delta ({\bf r}-{\bf r'})
\eqno(13)
$$
 where $\delta$ is the Dirac function. In addition,
 we assume that the spatial region
 is a delineated drainage basin with zero boundary condition which is
 fulfilled by $G({\bf r},{\bf r'},p)$. 

Eqs. (1-4)  can be 
 interpreted as Fokker-Planck stochastic processes (Gardiner, 1985) 
 with diffusion and drift components. Such an interpretation
 opens new possibilities to solve these equations through a 
 {\it simulation} of the underlying process
  utilizing stochastic methods (Gardiner, 1985). This type of
 Monte Carlo approach is one of the modern alternatives to finite element
 or finite difference approaches and is being explored
  in computational fluid dynamics or
 in quantum Monte Carlo methods for solving the Schr\"odinger
 equation (see
Mitas, 1996, and references therein).
 Very briefly, the solution is obtained as follows. 
 A number of sampling points distributed according 
 to the source $ {\cal S}({\bf r'})$ is generated.
 The sampling points are then propagated
 according to the function $G({\bf r},{\bf r'},p)$ and averaging 
 of path samples provides 
 an estimation of the actual solution $ \gamma ({\bf r})$ with
 a statistical accuracy  
 proportional to $1/\sqrt{M}$ where $M$ is the number of samples. 
 Fig 1. illustrates the solution with increasing 
number of walkers (animated illustration of  this method
  can be found in Mitas {\sl et al.} (1997)).

 The Monte Carlo technique has several unique advantages when compared
 with more traditional methods.
 It is very robust and enables studies 
 for spatially complex cases  with minimum of 
 manual preprocessing of input data.
Moreover, rough solutions, which identify the major sediment
concentrations and erosion/deposition patterns can be estimated quickly,
allowing us to carry out preliminary quantitative studies or
 to rapidly extract qualitative trends by parameter scans.
 In addition,  Monte Carlo methods are tailored to
the new generation of computers as they provide scalability from
a single workstation to large parallel
 machines due to the independence of sampling points. Therefore, the methods
 are useful both for everyday exploratory work using a desktop computer
 and for large, cutting-edge applications using high performance
computing. 

\bsn
\bsn
{\bf
2.4 Multiscale implementation for data with spatially variable resolution
}
\bsn
With the growing capabilities to collect geospatial data 
from various sources using different technologies 
the data set representing the studied landscape can be very heterogeneous
with different coverage, resolution, detail and accuracy. To fully
exploit the best data available for a given area 
algorithms should be able to make effective use
of these heterogeneous data sets, maximize the accuracy and detail of
predictions while performing computations efficiently.

To address this need 
we have reformulated the solution through the Green's function given 
by Eq. (12)  for accommodation of {\it spatially variable accuracy}
 and {\it resolution}.  The integral (12) can be multiplied by a {\it reweighting } 
 function  $W({\bf r})$ 
$$
W({\bf r})\gamma({\bf r})= \int_0^{\infty}\int W({\bf r})G({\bf r},{\bf r'},p)
 {\cal S}({\bf r'})  d{\bf r'}dp =
\int_0^{\infty}\int G^*({\bf r},{\bf r'},p)
 {\cal S}({\bf r'})  d{\bf r'}dp 
\eqno(14)
$$
 which is equal to the appropriate increase in accuracy ($W({\bf r})>1$)
 in the regions of interest while it is 
  unity elsewhere. The function $W({\bf r})$ can change (abruptly or smoothly)
  between regions with unequal  
 resolutions  and in fact, can be optimally adapted to the quality of
 input data (terrain, soils, etc) so that the accurate solution 
 is calculated only in the regions with correspondingly accurate inputs.
 The reweighted Green's function $G^*({\bf r},{\bf r'},p)$, in effect,
  introduces higher density of sampling points in the region with
 large $W({\bf r})$. The statistical noise will be spatially
 variable as
    $\approx 1/[W({\bf r})\sqrt{M}]$, where $M$ is the average number
 of samples  resulting in the accuracy increase for the areas
 with $W({\bf r}) > 1$.

The Eqs. (1-4) describe the water and sediment flow
at a spatial scale equal or larger than an average distance
between rills (i.e., grid cell size $\ge$ 1m) and therefore
the presented approach allows us to perform landscape scale
simulations at variable spatial resolutions from one to hundreds of meters,
depending on the complexity and importance of
studied subregions. 

\bsn
\bsn
{\bigrm
3 Results
}

\bsn
\bsn
{\bf
3.1 High resolution erosion and deposition patterns for spatially
complex conditions
}

\bsn
We have evaluated the capabilities of the presented model 
using an $\approx 1 \, km^2$ subarea of the Scheyern experimental farm 
(Auerswald {\sl et al.} 1996).
The measured elevation data were interpolated to a $2\, m$ resolution DEM
by regularized spline with tension (Mitasova \& Mitas, 1993)
and land cover and soil data were used to estimate erodibility
and transportability parameters based on the literature
(Foster \& Meyer 1972, Flanagan \& Nearing 1995).

First, we have applied the model for traditional land use (Fig. 2a)
for two different situations:
\item{(a) } dense grass in the meadow area, bare soil in the arable area 
and an extreme storm event; 
\item{(b) } vegetation cover everywhere and a lower intensity rainfall event.

\ni
We compared the results of simulations  
with spatial distributions of observed  colluvial deposits (Fig. 3a) and
linear erosion features digitized from aerial photographs (Fig. 2a).

For case (a), prevailing 
{\sl detachment limited} erosion ($\sigma \to 0$)
is predicted for the bare soil area,
due to the high transporting capacity of fast moving water.
The net erosion $D({\bf r}) \approx D_c({\bf r})$
and almost all detached sediment is transported to stream
while deposition is restricted to small concave areas and channels 
(Fig. 2 b,c).
This represents a situation close to the one observed
after an extreme storm event in 1993 when 
extensive rilling occurred, with only 7 \% of eroded sediment
deposited within the area (Auerswald {\sl et al.} 1996).

For the case (b), with smaller transport capacities
and ($\sigma > 1$) the erosion process is close to 
{\sl sediment transport capacity limited} case %($\sigma \to \infty$),
when $|{\bf q}_s({\bf r})| \approx T({\bf r})$. 
The erosion rates are lower and the model predicts large extent 
of areas with deposition (Fig. 3 b,c).
Such a behavior is close to the observed distribution of colluvial 
deposits (Fig. 3 a). For this case, the net erosion
and deposition can be approximated as 

$$ 
D({\bf r}) = \nabla\cdot |{\bf q}_s({\bf r})| \approx
\nabla\cdot \bigl[T({\bf r}){\bf s_0}({\bf r})\bigr]
\eqno(15)
$$

\ni
where ${\bf s_0}({\bf r})$ is the
unit vector in the steepest slope direction and the erosion and 
deposition pattern is significantly inlfuenced by terrain shape.
The impact of topography can be demonstrated
by substituing for  $T({\bf r})= K_t \rho_w\, g h({\bf r})
\sin \beta({\bf r})$ and deriving the net erosion and deposition
as a function of water depth and terrain curvatures:

$$
D({\bf r}) = K_t\rho_w\, g \Bigl\{ [\nabla h({\bf r})]
\cdot {\bf s_0 }({\bf r}) \sin\beta({\bf r}) -
h({\bf r})\bigl[\kappa_p({\bf r})+\kappa_t({\bf r})\bigr] \Bigr\}
\eqno (16)
$$

where $\kappa_p({\bf r})$ is the profile curvature ( curvature
in the direction of the steepest slope),
$\kappa_t({\bf r})$ is the tangential curvature (curvature in the
direction perpendicular to the gradient,
 i.e., tangential to the terrain
surface isoline; for curvature definitions
see Mitasova and Hofierka, 1993).
Therefore, according to (16)
the spatial distribution of erosion/deposition
is given by the directional derivative of the overland flow depth
(first term)  {\sl and} by the local geometry of terrain (second term),
including {\sl both} profile and tangential curvatures (Fig. 4).
The commonly used univariate formulations of erosion models
 neglect the influence of tangential curvature (Fig. 4),
underestimating thus the deposition in areas with $\kappa_t<0$
(tangential concavity) and underestimating erosion in areas with $\kappa_t>0$
(tangential convexity) as demonstrated by Mitas and Mitasova (1998).
 Although the above analysis strictly applies to the case when
 $p=1$ in Eq. (7), it is possible to derive similar expression
 with a general exponent $p\not= 1$ and the qualitative conclusions
 remain the same.

From the point of view of land use management, it is important to
note that for both simulations (a), (b)
the highest rates of net erosion as well as net deposition
were predicted in hollows with high concentrated
sediment flow (Figs. 2 b,c; 3 b,c).  Field measurements confirm that this
area has the thickest layers of colluvial deposits
but also large linear erosion features were observed here
after a strong storm (Auerswald {\sl et al.} 1996).
The second highest erosion is predicted on upper convex parts
of hillslopes where the highest loss of radio-tracers and
the lowest yields were observed (Auerswald, personal communication).
Increased erosion is predicted also
for bare narrow stripes below the grass areas,
where water accelerates after depositing the sediment.
The major difference in spatial patterns between the two modeled 
 cases is the spatial extent of erosion/deposition.
In the simulation (a) 93\% of the area experienced erosion
while simulation (b) predicted erosion only for 70\% of the area,
with extent of deposition close to the observed
spatial distributions of colluvial deposits.
Deposition was also predicted
at upper edges of meadows, where the 
rills follow  the borderline between
the grass and bare soil areas and where deposits were found
in the upper convex part of the grass hillslope (Figs. 2, 3).

\bsn
\bsn
{\bf
3.2  Land use design for improved erosion prevention
}

\bsn
The presented bivariate erosion model can be used for analyzing and designing
the placement of selected erosion protection measures based on land cover,
as illustrated by
the following simple example. First, we used the model to identify
locations with the highest erosion risk, assuming a uniform land use. Then,
the protective grass cover was distributed to the high risk areas
while preserving the extent of grass cover at the original 30\% of the area
(Fig. 5). We performed a simulation with the new land use
to evaluate its effectiveness. The results demonstrate that
the new design has a potential to dramatically reduce soil loss and
sediment loads in the ephemeral streams (Fig. 5).
The crest in sediment flow in
the valley disappears and is replaced by light deposition within
the grassways, while the maximum and total rates of erosion are
significantly reduced.  We have found, that the effectiveness
of this design depends on differences in roughness, combination of
 very smooth bare soil and very
dense grassway resulted in predictions of higher
erosion along the borders of the grass waterway.
%, from 1.7 to 0.3 $kg/m^2s$, and from
%635 to 306 $kg/s$, respectively.
%total excess? erosion 24 to 14 $kg/s$) (Table 2).
%deposition -4.6 to -0.86, -611 to 292,
%max sed.flow.rate 34 to 3.5$kg/ms$
%reduced to half.
It is interesting to note, that the land use design obtained
by this rather simple computational procedure, using only the
elevation data, had several common features with the
sustainable land use design proposed and implemented
in 1993 at the farm, based on extensive experimental work
(Auerswald et al., 1996). 
The design uses significantly
higher proportion of permanent grass cover and fallow and increases
cover/roughness in the hop field in the center of the valley. 
The results show that this design keeps increased amounts of soil moisture
 and, at higher cost, reduces the net erosion even further.

\bsn
\bsn
{\bf
3.3 Simulations with spatially variable resolution and accuracy.
}
\bsn
The necessity for spatially heterogeneous modeling
of fluxes can arise due to various reasons, for example, when the
modeled region is represented by various data sets at different
resolutions or when the region is very large
with spatially variable complexity and 
only a subarea of high spatial variability needs to be simulated
at high resolution while for the rest of the reagion lower resolution
is sufficient.

We illustrate the concept of multiscale modeling using a 16 km$^2$
region represented by a 10m resolution DEM (400c x 400r) with
a smaller 1 km$^2$ subregion mapped with more detail and represented by
a 2m resolution DEM (680c x 395r). Resampling the DEM for entire 
region to 2m resolution would lead to 4 million grid cell DEM 
requiring substantial increase in computational requirements 
with little gain in detail for low resolution area. With spatially 
variable resolution we simulate the water flow at 10m resolution 
for the entire area and save the walkers entering the high resolution
subregion. Water flow for the subregion is then simulated at 2m
resolution using the saved resampled walkers as inputs, as illustrated
by Fig. 6. Probably the most significant research and application 
potential of this approach is in multiscale/multiprocess simulations
which would allow us to run simulations using the models which are the
most appropriate for the given scale. For example, by transfering
the water and sediment from average overland flow model for low resolution
area to rill-interril-channel model for subregion represented at
high, submeter resolution.

\bsn
\bsn
{\bigrm
4 Conclusions 
}

\bsn
 The stochastic method of solving the first principles
 equations using Green's function Monte Carlo technique
 provided us with a valuable research tool of a much needed
 robustness and flexibility. It enabled us to investigate
 several important issues such as impact of complex terrain
on spatial patterns of erosion/deposition under different
regimes of sediment transport, impact of land use design
on effective erosion prevention and the possibility to
simulate flow at spatially variable resolution.
 
Comparisons of simulation results with field data indicate
that the study area at Scheyern experimental farm exhibits
 colluvial deposition patterns consistent with prevailing 
transport capacity limited regime, modeled for the conditions
with vegetative cover and small event.
 On the other hand, as we conjecture from our preliminary
investigations, the detachment limited case seems to dominate
 the creation of short term
erosion features, such as rilling  and gully formation,
which are related to large rainstorm events, bare soil 
conditions and high transport capacities concentrated 
to short time intervals.
The bivariate model was also capable to predict the location of
thalweg erosion by simulating high rates of both erosion and
deposition in valleys.

 A suggestion of Nearing at al (1997) that
the sediment loads from rills may be more strongly influenced by
a sediment transport limit rather than by the soil detachment, in general,
 agrees with our results. This seems to be true  especially
in complex terrain conditions where transport capacity changes
significantly due to variations in terrain shape and
cover affecting significantly
the distribution and amplitudes of the water flow field.
 Our calculations and analysis also suggest that the
 sediment transport capacity plays a more important role than
 anticipated by the previous research which focused on erodibility
 as the key control quantity. Obviously, a subtle and spatially
 variable interplay between erodibility and transport capacity can
 influence the processes in a profound way.
 We believe that this complexity
 clearly points out towards the importance of
 the high resolution 2D simulations.

The presented model enabled us to study the impact 
of different land cover designs on the erosion processes
and estimate their effectivenes for erosion protection,
demonstrating the possibilities to create
powerful tools for computer aided sustainable land
use design.

The presented model was implemented as an independent
computational module named
SIMWE (SIMulation of Water Erosion) and linked to  the Geographic
Resources Analysis System (GRASS).


\bigskip
\ni
{\bf Acknowledgments}
We would like to express our gratitude to K. Auerswald
from the Technical University
Muenchen, Germany,
and to S. Warren from the USA Construction
Engineering Research Laboratory (USA CERL) for providing the
experimental data used in this paper and for valuable
discussions. We acknowledge the support by
 the Strategic Environmental Research
and Development Program (SERDP) and the USA CERL.
 We thank Douglas Johnston from the GMS Laboratory at the University
of Illinois for his support and discussions, and to William Brown
for his contribution to the new GRASS tools used in this paper.

\vfill
\eject
\bsn
{\bigbrm Bibliography}

\bsn
\bsn
\bsn
\bsn
\bsn
\bsn
\bsn
\bsn
\ss
\hangindent=2pc \hangafter=1
Auerswald, K., A. Eicher, J. Filser, A. Kammerer, M. Kainz, R. Rackwitz, J. Schulein,
H. Wommer, S. Weigland, \& K. Weinfurtner (1996)  The Scheyern project of the FAM.
In: {\sl Development and Implementation of Soil Conservation Strategies for
Sustainable Land Use} (ed. by H. Stanjek) (Int. Congress of ESSC, Tour Guide II),
25-68, Technische Universitaet Muenchen, Freising-Weihenstephan, Germany.

\ni
\hangindent=2pc \hangafter=1
Bennet, J. P., Concepts of Mathematical Modeling of Sediment Yield,
{\sl Water Resources Research}, 10, 485-496, 1974.

\ni
\hangindent=2pc \hangafter=1
Dingman, S. L., {\sl Fluvial hydrology}, Freeman, New York, 1984

\ni
\hangindent=2pc \hangafter=1
Flacke, W., K. Auerswald, and L. Neufang, Combining a modified USLE
with a digital terrain model for computing high resolution maps of
soil loss resulting from rain wash, {\sl Catena}, 17, 383-397, 1990.

\ni
\hangindent=2pc \hangafter=1
Flanagan, D.C. \& M.A. Nearing (eds.)  (1995)
{\sl Water Erosion Prediction Project}.
Report no. 10, National Soil Erosion Lab., USDA ARS, West Laffayette, 
Indiana, USA.

\ni
\hangindent=2pc \hangafter=1
Foster, G. R. \& L. D. Meyer (1972) A closed-form erosion equation for upland areas.
In: {\sl Sedimentation: Symposium to Honor Prof. H.A. Einstein}
(ed. by H. W. Shen), 12.1-12.19, Colorado State University, Ft. Collins, USA.

\ni
\hangindent=2pc \hangafter=1
Gardiner, C. W. (1985) {\sl Handbook of Stochastic Methods for Physics, Chemistry, and the
Natural Sciences}. Springer, Berlin.

\ni
\hangindent=2pc \hangafter=1
Haan, C. T., B. J. Barfield, and J. C. Hayes, {\sl Design Hydrology and
Sedimentology for Small Catchments}, pp. 242-243, Academic Press,  1994

\ni
\hangindent=2pc \hangafter=1
Maidment, D. R. (1996) Environmental modeling within GIS.
In {\sl GIS and Environmental Modeling: Progress and Research Issues}
(ed. by M. F Goodchild, L. T. Steyaert, \& B. O. Parks), 315-324, GIS World, Inc.,
Ft. Collins, USA.

\ni
\hangindent=2pc \hangafter=1
Meyer, L.D., and W. H. Wischmeier, W.H., Mathematical simulation
of the process of soil erosion by water. {\sl Transactions of the ASAE},
12, 754-758, 1969.

\ni
\hangindent=2pc \hangafter=1
Mitas, L. (1996) Electronic structure by Quantum Monte Carlo:
atoms, molecules and solids. {\sl Computer Physics Communications} {\bf 97}, 107-117.

\ni
\hangindent=2pc \hangafter=1
Mitas L., \&  Mitasova, H. (1998) Distributed soil erosion simulations
for effective landuse managment. {\sl Water Resources Research} {\bf 23}, 505-516.

\ni
\hangindent=2pc \hangafter=1
Mitas, L., W. M. Brown \& H. Mitasova  (1997) Role of dynamic cartography in simulations
of landscape processes based on multi-variate fields.
{\sl Computers \& Geosciences} {\bf 23}, 437-446.

\ni
\hangindent=2pc \hangafter=1
Mitasova,  H. \& L. Mitas (1993) Interpolation by Regularized Spline with Tension:
I. Theory and  implementation. {\sl Mathematical Geology} {\bf 25}, 641-55.

\ni
\hangindent=2pc \hangafter=1
Mitasova, H. \& J. Hofierka (1993) Interpolation by Regularized Spline with Tension:
 II. Application to terrain modeling and surface geometry  analysis,
{\sl Mathematical  Geology} {\bf 25}, 657-69.

\ni
\hangindent=2pc \hangafter=1
Mitasova, H., L. Mitas, W. M. Brown, D. P. Gerdes, I. Kosinovsky \& T. Baker (1995)
Modelling spatially and temporally distributed phenomena: new methods and tools for GRASS
GIS. {\sl Int. Journal of Geographical Information Systems} {\bf 9}, 433-446.

\ni
\hangindent=2pc \hangafter=1
Moore, I. D.,  A. K. Turner, J. P. Wilson, S. K. Jensen \& L. E. Band (1993)
GIS and land surface-subsurface process modeling. In
{\sl Geographic Information Systems and Environmental Modeling} (ed. by M. F Goodchild,
L. T. Steyaert, \& B. O. Parks), 196-230, Oxford University Press,  New York.

\ni
\hangindent=2pc \hangafter=1
Nearing M.A., Norton L.D., Bulgakov D.A., Larionov G.A., West L.T. \& Dontsova K.M.
(1997) Hydraulics and erosion in eroding rills. {\sl Water Resources Research} {\bf 33}, 865-876.

\ni
\hangindent=2pc \hangafter=1
Saghafian, B., P. Y. Julien \& F. L. Ogden (1995) Similarity in catchment response
1. Stationary rainstorms, {\sl Water Resources Research} {\bf 31}(6), 1533-1541.


\vfill
\eject

\ss
{\bf Figure captions}

\ni
{\sl Note:} Color versions of figures are available at:

\ni
{\sl http://www2.gis.uiuc.edu/papers/unesco99.html}


\bsn
{\sl Fig. 1 Monte Carlo solution of sediment transport equation
using increasing number of walkers: (a) sediment flow rate, (b)
net erosion and deposition

\bsn
Fig. 2 Terrain model with (a) traditional land use and
observed rills and gullies, (b) simulated sediment flow,
(c) erosion/deposition pattern for bare soil in the
agricultural field and a 70 mm/hr rainfall excess.

\bsn
Fig. 3 Terrain model with (a) observed depths of colluvial deposits,
(b) simulated sediment flow,
(c) erosion/deposition pattern for vegetative cover in the
agricultural field and a 20 mm/hr rainfall excess.

\bsn
Fig. 4 Univariate and bivariate computation of net erosion
and deposition for transport capacity limiting case:
: (a) erosion/deposition as a change
 of 1D sediment flow along flow line:$ d |{\bf q_s}| / ds$,
(b) erosion/deposition as a divergence of 2D vector field
representing sediment flow:$ div {\bf q_s}$.

\bsn
Fig. 5 Terrain model with (a) computer designed land use,
(b) simulated sediment flow,
(c) erosion/deposition pattern for bare soil in the
agricultural field and a 70 mm/hr rainfall excess.

\bsn
Fig. 6 Multi-scale simulation of steady state water depth:
(a)water depth at combined 10m and 2m resolution for the entire
study area, (b) selected subarea at 10m resolution, (c) the
 same selected
subarea at 2m resolution.}


\vfill
\eject
{\bf Corresponding author:}

Lubos Mitas, PhD

Research Scientist

National Center for Supercomputing Applications

405 N. Mathews Ave

University of Illinois at Urbana-Champaign

Urbana, IL 61801, USA

Office ph.: + 217 244 1971

Fax: + 217 244 2909

E-mail: lmitas@ncsa.uiuc.edu

\bye
