a
b
c
For all Figures, clicking on the img will display an enlarged version
From Abrahart, R.J. (ed.) (1998). Proceedings of the Third International Conference on GeoComputation (CD), University of Bristol, 17-19 September 1998
David Favis-Mortlock
Environmental Change Unit, University of Oxford, 5 South Parks Road, Oxford
OX1 3UB, United Kingdom.
John Boardman
School of Geography and Environmental Change Unit, University of Oxford, 5 South
Parks Road, Oxford OX1 3UB, United Kingdom
Tony Parsons
Department of Geography, University of Leicester, University Road, Leicester
LE1 7RH, United Kingdom
Bruce Lascelles
Department of Geography, University of Leicester, University Road, Leicester
LE1 7RH, United Kingdom
The initial location and subsequent spatial development of hillslope rill systems cannot be forecast by current models of soil erosion by water. Modelled estimates of erosion rate can thus only represent some spatially generalised result e.g. a mean erosion rate for a field. Yet within the field, the severity of erosion will vary markedly from place to place. Such spatial variability results in part from the way in which rill systems develop. Favis-Mortlock (1996, 1998) noted that some rills are more 'successful' and grow in preference to others, with whom they compete for runoff. This suggested an evolutionary analogy. Based on work on self-organising dynamic systems, a model ('RillGrow 1') was constructed to evaluate this concept. The model applies simple rules at a millimetre scale, to govern the iterative interaction between microtopography, runoff routing and soil loss. No explicit separation is made between rill and interrill processes, which are deemed instead to be end-members of a continuum.
Initial results proved the model capable of realistically simulating both the initiation of a rill network and its subsequent development. Networks reproduce many 'emergent' features of observed systems on a scale of centimetres to metres, including a decrease of mean rill spacing and a non-linear increase of total erosion with steepening slope. In a subsequent validation study (Favis-Mortlock et al., 1998; Favis-Mortlock et al., in preparation; Guerra et al., in preparation), digital photogrammetry was used to generate DEMs of microtopography before and after simulated rainfall on a small field plot. Patterns of erosive flow generated by RillGrow compared well with observations.
Despite these encouraging results, the original RillGrow has several fundamental limitations. It does not operate in a true time domain; important process descriptions (e.g. deposition) are ignored; and the hydraulics of rill initiation are oversimplified. Thus a new version of the model ('RillGrow 2') is being developed, in conjunction with a programme of laboratory-based experimental work. The RillGrow 2 model, even at this early stage in its development, appears not only to be able to replicate the successes of the earlier model, but to exceed them.
Keywords: modelling, soil erosion, self-organisation, emergence, rill initiation
Current models for soil erosion by water such as WEPP (Lane and Nearing, 1989) and EUROSEM (Morgan et al., 1998) consider an eroding hillslope to be conceptually partitioned into rill and interrill areas. Rills are envisaged as concentrating interrill runoff into hydraulically efficient zones that transport much of the flow and sediment off the hillslope; rills are represented as channels and the laws of hydraulics are applied to flow in these channels. Soil loss in rill and interrill areas is calculated separately and summed to yield a prediction of total soil erosion.
This conceptual model derives from the work of Meyer and Wischmeier (1969). It works well in many respects (e.g. Favis-Mortlock, 1998a). However, any predictive model for soil erosion which adopt this approach will inevitably inherit a number of theoretical and practical shortcomings.
There is thus a need for a modelling approach which is capable of explicitly specifying the spatial location of the initiation and subsequent development of rills.
Favis-Mortlock (1996; 1998b) noted that some rills are more 'successful' than others, in the sense that such rills preferentially grow in size and sustain flow throughout a rainfall event. In one sense, they can be envisaged as 'competing' for runoff. Thus from an initial population of many tiny rills (i.e. micro-rills or 'traces') only a subset will subsequently develop into larger rills and eventually form part of a connected network (Table 1).
| Category of rill | Rate of growth | Effectiveness during rainfall event |
| Successful | Higher | Becomes major carrier for runoff and eroded soil for part of hillslope; may 'capture' weaker rills |
| Unsuccessful | Lower | Becomes less and less important as a carrier for runoff and sediment; may eventually be 'captured' or become completely inactive |
Table 1. 'Successful' and 'unsuccessful' rills. From Favis-Mortlock (1996)
If the 'successful' microrills prosper at the expense of their 'unsuccessful' siblings, then it might be expected that at the end of the rainfall event only the 'successful category will have left readily observable evidence of their existence. In addition, if the assumption is made that 'success' depends upon the rate of growth of the rill, then 'success' is likely to also depend largely on the volume of runoff which passes through the growing rill, since a greater flow will result in more erosion within the rill. Thus the probability of ‘success’ (and hence the development of observable rills) may be expected to increase downslope, as upslope catchment area increases.
| Position on slope | Rill morphology |
| Upper | Small and very short 'trace' rills |
| Middle | 'Discontinuous' rills, well defined but usually less than a metre in length |
| Lower middle to bottom | Continuous rills |
Table 2. Rill morphology and position on slope. From data in Evans (1995)
This might be expected to produce a downslope spatial patterning of rill morphology similar to that described by Evans (1995) from field observation of sandy and silty soils in the UK. Where slope length is sufficiently great, similar sequences have been noted in laboratory studies (e.g. Bryan and Poesen, 1989). In another long flume study, Slattery and Bryan (1992) expressed a 'competitive' view of rill initiation by noting that "Rill incision always arose through the formation of a knickpoint, although not all knickpoints became rills".
While downslope position is undoubtedly important in determining rill 'success' or 'failure' at any point on a hillslope, this is a static attribute which cannot change during a rainfall event. Thus some other factor must be responsible for the observed temporal dynamism of rill 'competition'.
Favis-Mortlock (1996; 1998b) suggested that microtopography could fit this role. It seems reasonable to expect that the vagaries of microtopography might also play an important part in determining the upslope supply of runoff and sediment at any point on the hillslope. But microtopography has an important extra attribute: it is not temporally static, since erosional processes must (by definition) modify the soil's surface (e.g. Zobeck and Onstad, 1987). As a result, runoff during the latter part of a rainfall event will flow over a soil surface which may be very different from that encountered by runoff earlier in the storm. In this sense, erosion creates its own surface.
Erosive modification of microtopography thus constitutes a feedback loop. This might be expected to operate in a positive sense, since the most 'successful' rills (i.e. those conveying the most runoff) will modify the local microtopography to the greatest extent, and so will most effectively amplify their chances of capturing and conveying future runoff (cf. Figure 2 in Lenton, 1998). The existence of a positive feedback loop is one prerequisite for the existence of self-organising dynamical systems (e.g. Casti, 1994; Cohen and Stewart, 1994; Coveney and Highfield, 1995; Bak, 1997).
Self-organising dynamical systems tend toward greater systemic orderliness over time: in other words, the system produces a local decrease in entropy (cf. Ruelle, 1993) which is the result of interactions between its components. Such system-generated orderliness is described as 'emergent': emergent properties are those that become apparent at a certain level of system complexity, but do not exist at lower levels (C.D. Broad, early 1920s: quoted by Capra, 1997). While the majority of research into these types of systems is of recent date, the underlying concept has a long history. Self-organising systems are essentially synergistic — 'greater than the sum of their parts' — a notion which dates back to Immanuel Kant in the 1790s (Capra, 1997).
Several workers have investigated the applicability of concepts of self-organisation to natural systems, including landscape evolution (e.g. Green, 1993; Perez-Trejo, 1993; Coulthard et al., 1996) and other aspects of geomorphology (e.g. Scheidegger, 1994; Thorn and Welford, 1994; Phillips, 1996, 1997; Hergarten, 1997; Murray and Paola, 1997). Favis-Mortlock (1996; 1998b) suggested that the observable features of hillslope erosion may be manifestations of a underlying self-organising system (Table 3). Note that these features cover a wide range of spatial scales.
| Component | Role in a self-organising systems view of erosion | Scale |
| Raindrop | Input | Millimetres |
| Microtopography | Input | Vertical: millimetres Horizontal: millimetres to centimetres |
| Microrill | 'Emergent' output | Millimetres to centimetres |
| Rill | 'Emergent' output | Centimetres to tens of metres |
| Hillslope rill system | 'Emergent' output | Tens to hundreds of metres |
Table 3. A self-organising systems view of hillslope erosion
From this perspective, several commonly observed responses of erosional systems are also emergent outputs of this self-organising system.
A partial test of this hypothesis is inductive: a modelling approach to hillslope based upon concepts of self-organisation should also produce these responses as emergent outputs.
A model ('RillGrow 1') was constructed (Favis-Mortlock, 1996; 1998b) which applied simple rules to the movement of individual runoff 'packets' on a grid of microtopographic heights. Initial results proved the model capable of realistically simulating both the initiation of a rill network and its subsequent development. Networks reproduce many 'emergent' features of observed systems on a scale of centimetres to metres, including a decrease of mean rill spacing and a non-linear increase of total erosion with steepening slope. In a subsequent validation study (Favis-Mortlock et al., 1998; Favis-Mortlock et al., in preparation; Guerra et al., in preparation), digital photogrammetry was used to generate DEMs of microtopography before and after simulated rainfall on a small field plot. Patterns of erosive flow generated by RillGrow 1 compared well with observations.
Despite these encouraging results however, the original RillGrow model has a number of fundamental limitations. It does not operate in a true time domain; important process descriptions (e.g. deposition) are ignored; and the hydraulics of rill initiation are oversimplified. Thus a new version of the model ('RillGrow 2') is being developed, in conjunction with a programme of laboratory-based experimental work. While work on RillGrow 2 is ongoing, the model and some early results are presented here.
As in RillGrow 1, the model's spatial domain is an x-y grid (as in so-called cellular automata models: Wolfram, 1982). This represents the bare soil area on which erosion occurs. Each element of the grid is a cell: each cell is assumed to be square in plan and vertically homogenous, and each cell has a height which represents the microtopographic elevation at the centre of the cell. Also as in the first version of the model, runoff is considered to be discretised into 'packets', which move from cell to cell. Each packet may have any depth, but is assumed to uniformly cover the entire cell. No distinction is made between interrill and rill areas.
RillGrow 2 differs from the earlier version however in that multiple packets are routed concurrently across the grid. The model thus operates in a true time domain. Table 4 gives the model's main routing algorithm.
For each iteration:
|
Table 4. The main routing algorithm used in RillGrow 2
Raindrop volumes are sampled from a normal distribution of known mean and standard deviation. Splash erosion is not yet modelled. Infiltration is also not yet implemented, thus (as in RillGrow 1) all water on the grid is assumed to be rainfall excess. Two components are summed to give a travel velocity for each cell: a 'base' component and a depth-dependent component. 'Base' velocity for each cell is sampled from a normal distribution of known mean and standard deviation. Currently there is no spatial correlation of 'base' packet velocities. At present, the depth-dependent component of packet velocity is calculated arbitrarily as a function of packet depth raised to an integer power (currently two).
As in the first version of the model, a logistic (i.e. s-curve) relationship between unit sediment load and stream power is used to estimate the soil loss resulting from packet flow. This was derived experimentally by Nearing et al. (1997) and has the form:
| log(qs) = {a + b eg +d log(w )} / {1 + eg +d log(w )} |
(1) |
where:
| qs = unit sediment load (g s-1 cm-1) |
| a , b , g , d are constants |
and:
| w = r w g S q |
(2) |
where:
| w = stream power (kg s-3) |
| r w = density of water (kg m-3) |
| g = gravitational acceleration (m s-2) |
| S = slope |
| q = unit discharge of water (m2 s-1) |
As in the first version of the model, there are several caveats concerning the use of this relationship.
One modification which is necessary to this relationship in RillGrow 2 concerns the attenuating effect of water depth. If a small volume of outflow occurs from a cell with a relatively large depth of water, it seems reasonable to expect that the resulting shear stress — and hence detachment — will be lower than in the case of an identical outflow from a shallower water depth, all else being equal. An arbitrary expression is used to estimate this depth-attenuation as follows:
| qsactual = qs × e(k× W(E-W)) |
(3) |
where:
| qsactual = unit sediment load (g s-1 cm-1) used to estimate soil loss |
| qs = unit sediment load (g s-1 cm-1) as calculated by equation (1) |
| W = water depth on the cell (cm) |
| E = a pre-specified 'effective' water depth (cm) |
| k = constant |
Deposition is not yet implemented, so that transport capacity is effectively assumed to be infinite.
RillGrow 2 is implemented in C++. Both input and output values for each cell are in a format which can be read by the Idrisi Geographical Information System (Eastman, 1997). The model also produces 'real-time' graphics.
Twenty-six simulations were carried out using four DEMs: these are summarised in Table 5. The 'smooth' DEM has no microtopographic relief. The 'inset' DEM also has no relief over 8/9 of its area, but has a fractally-rough (Favis-Mortlock, 1996) area inset into the central portion (Figure 1). The 'laser' DEM (Figure 1) was produced by laser-scanning a level soil pan (Favis-Mortlock, 1996; 1998b), and the 'photo' DEM (Figure 1) was derived by photogrammetric survey of a small field plot (Favis-Mortlock et al., 1998; Guerra et al., in preparation). Maximum slope on this area is approximately 7%.
Approximately the same total rainfall was applied for every simulation. Each simulation was replicated using two rainfall intensities ('high' and 'low') were used (Table 5). For the 'smooth' and 'inset' DEMs, each run was further replicated with differing standard deviations of 'base' packet velocity (2 mm sec-1 and zero). Replicate runs for the 'laser' DEM were carried out for a range of slope angles (Table 5). Runs for the 'photo' DEMs were replicated using differing random number seeds.
| DEM name | DEM description | DEM grid: x size by y size by cell side (mm) |
Slope angle (%) |
Rainfall rate (mm hr-1) |
Standard deviation of rainfall rate (mm hr-1) |
Number of iterations |
| SMOOTH | Smooth, 1.5 m x 4 m | 300 x 800 x 5 |
20 |
10 |
2 |
5000 |
| INSET | Smooth with fractal area inset, 1.5 m x 1.5 m | 300 x 300 x 5 |
20 |
10 |
2 |
5000 |
| LASER | Laser scanned soil pan (see Favis-Mortlock, 1998), 2.2 m x 2.2 m | 220 x 220 x 10 |
0 to 30 in 5% increments |
10 |
2 |
10 000 |
| PHOTO | From a field plot, using photogrammetry (see Guerra et al., in prep.), 0.42 m x 0.8 m | 167 x 319 x 2.5 |
0 |
10 |
2 |
10 000 |
| SMOOTH | Smooth, 1.5 m x 4 m | 300 x 800 x 5 |
20 |
2 |
0.4 |
25 000 |
| INSET | Smooth with fractal area inset, 1.5 m x 1.5 m | 300 x 300 x 5 |
20 |
2 |
0.4 |
25 000 |
| LASER | Laser scanned soil pan (see Favis-Mortlock, 1998), 2.2 m x 2.2 m | 220 x 220 x 10 |
0 to 30 in 5% increments |
2 |
0.4 |
50 000 |
| PHOTO | From a field plot, using photogrammetry (see Guerra et al., in prep.), 0.42 m x 0.8 m | 167 x 319 x 2.5 |
0 |
2 |
0.4 |
10 000 |
Table 5. Simulation details
|
a |
b |
|
c |
Figure 1. Plan views of the initial surfaces
of three DEMs used in this study: (a) 'inset' DEM (b) laser-scanned DEM
(c) photogrammetrically-derived DEM. Note that the illustrations are not
to scale: see Table 5 for sizes. The downslope end is to the bottom in each
case
For all Figures, clicking on the img will display an enlarged version |
Details of other model parameters are given in Table 6. These were used for all runs (with the exception of the standard deviation of 'base' packet velocity: this was set to zero for some runs using the 'smooth' and 'inset' DEMs).
| Parameter |
Value |
| Mean raindrop volume |
5.0 mm3 |
| SD of raindrop volume |
2.0 mm3 |
| Timestep |
0.1 sec |
| Mean 'base' packet velocity |
10.0 mm sec-1 |
| Standard deviation of 'base' packet velocity |
2.0 mm sec-1 |
| Soil bulk density |
1.3 g cm-3 |
| Angle of repose |
40.0 % |
Table 6. Other parameters
Simulated water depth is shown in Figures 2 to 6. In all cases, greater depths are indicated by a darker shade of blue and are relative rather than absolute (i.e. a similar colour may refer to different depths in different diagrams). Flow is from top to bottom.
 |
| Figure 2. Simulated water depth at the end of a 'smooth' DEM simulation. High-intensity rainfall was used, with the standard deviation of 'base' packet velocity set to 2 mm sec-1 |
For all runs with the 'smooth' DEM, average water depth increased gradually and linearly downslope. However, instantaneous depth during some simulations presented patterns as seen in Figure 2. Wave-like patterns formed on the lower part of the slope early on in the run. Individual 'waves' migrated downslope during the course of the simulation, moving at a rather slower rate than the packets making them up (i.e. the packets moved 'through' the waves). During the simulation with high-intensity rainfall and zero 'base' packet velocity similar wave-like patterns were formed, but here they were rather less distinct and began further down slope. No patterns formed during the simulations with low-intensity rainfall.
 |
| Figure 3. Simulated water depth at the end of an 'inset' DEM simulation. Rainfall was high-intensity rainfall, and standard deviation of 'base' packet velocity set to zero |
Wave-like patterns also formed in some runs with the 'inset' DEM (Figure 3). Here they occurred only on the lower part of the slope in the smooth areas. They did not form in the fractally-rough area: here flow instead tended to concentrate into narrow zones. Similar wave-like patterns also formed during the simulation with high-intensity rainfall and 2 mm sec-1 'base' packet velocity, where they were rather more distinct and began further up-slope; flow in the fractally-rough area was similar to that shown in Figure 3. No wave-like patterns formed during the simulations with low-intensity rainfall.
Figures 4 and 5 show average (not instantaneous) water depths on the laser-scanned DEM for the high-intensity and low-intensity simulations respectively. At slope angles of zero and 5%, runoff tended to pond; rainfall intensity made little difference here. As slope angle was increased, there was a greater tendency for water to concentrate into channelised flow. While channelisation was very clear at the highest slope angles, the form of the channel networks differed for the high-intensity and low-intensity simulations. For the high-intensity simulations, there was a tendency to form braided networks of flow paths. Networks were more dendritic in form for the low-intensity simulations. In both cases, flow paths tended to be straighter on the steeper slopes.
a |
b |
c |
d |
e |
f |
g |
Figure 4. Average water depth at the end of simulations on the laser-scanned DEM; high-intensity rainfall. Slope angle is (a) 0% (b) 5% (c) 10% (d) 15% (e) 20% (f) 25% (g) 30% |
a |
b |
c |
d |
e |
f |
g |
Figure 5. Average water depth at the end of simulations on the laser-scanned DEM; low-intensity rainfall. Slope angle is (a) 0% (b) 5% (c) 10% (d) 15% (e) 20% (f) 25% (g) 30% |
Extensive ponding occurred on the photogrammetrically-derived DEM (Figure 6). This was broadly similar for both low-intensity and high-intensity rainfalls.
 |
| Figure 6. Average water depth on the photogrammetrically-derived DEM at the end of the simulation with low-intensity rainfall |
Simulated soil loss is shown in Figures 7 to 10. Again, depths are relative rather than absolute. Flow is from top to bottom.
Total soil loss on the 'smooth' DEM increased linearly downslope (not shown).
On the 'inset' DEM (Figure 7), soil loss also increased linearly downslope on the smooth portion. On the fractally rough area it was concentrated into narrow zones.
 |
| Figure 7. Total soil loss for the 'inset' DEM simulation: high-intensity rainfall and zero 'base' packet velocity |
a |
b |
c |
d |
e |
f |
g |
Figure 8. Patterns of soil loss at the end of simulations on the laser-scanned DEM; high-intensity rainfall. Slope angle is (a) 0% (b) 5% (c) 10% (d) 15% (e) 20% (f) 25% (g) 30% |
Soil loss on the laser-scanned DEM to some extent mirrors flow patterns. Little soil loss occurred at low slopes, irrespective of rainfall intensity. As slope angle increased there was an increased tendency for erosion to be concentrated into narrow zones. Erosion in these zones began at a series of isolated points, frequently in a broad along-contour zone along the areas of slightly higher microtopography (cf. Figures 1b and 8c). These points subsequently tended to link up. There are however clear differences in pattern for the steeper slope angles for differing rainfall intensities. As with the flow patterns, erosional patterns are more braided in form for the high-intensity runs, and more dendritic for the low-intensity runs. Similarly, these patterns tend to become straighter as slope angle increases. On the more dendritic networks, there is some tendency for increased depth downstream from a junction.
a |
b |
c |
d |
e |
f |
g |
Figure 9. Patterns of soil loss at the end of simulations on the laser-scanned DEM; low-intensity rainfall. Slope angle is (a) 0% (b) 5% (c) 10% (d) 15% (e) 20% (f) 25% (g) 30% |
Figure 10 shows the relationship between erosion patterns and microtopography. Note how the eroded channels (darker) avoid the higher (light-coloured) spots (aggregates).
 |
| Figure 10. Patterns of soil loss simulated on the laser-scanned DEM: low-intensity rainfall and 20% slope angle |
Simulated soil loss on the photogrammetrically-derived DEM was very low (not shown).
 |
| Figure 11. Total runoff and soil loss for the simulations with the laser-scanned DEM |
Total runoff volume (Figure 11) increases with slope angle up to about 10% slope; it is approximately constant thereafter (i.e. is independent of slope angle for the steeper slopes). Total runoff is higher for the low-intensity (i.e. longer) runs.
By contrast, total soil loss increases non-linearly with increased slope angle, more so for the simulations with high-intensity rainfall so that soil loss is markedly (about 25%) higher for the higher rainfall intensities on the steepest slopes.
For the following analysis, an arbitrary depth of soil loss (> 0.5 mm) was assumed to be 'rill' soil loss, with smaller depths representing 'interrill' soil loss (cf. Favis-Mortlock, 1996; 1998b).
a |
b |
| Figure 12. Average rill spacings for the simulations which used the laser-scanned DEM, (a) high-intensity rainfall (b) low intensity rainfall |
In general, average rill spacing decreases with increasing slope angle (Table 7), and decreases downslope (Figure 12). Absolute spacings are greater for the runs with low-intensity rainfall. For both intensities, tote that there are two zones (at about 1100 and 1600 mm downslope distance) where rill spacing temporarily increase. These correspond to area of relatively higher relief (cf. Figure 1b).
| Slope angle (%) |
Rainfall intensity |
|
|
High |
Low |
|
| 0 |
- |
- |
| 5 |
538 |
543 |
| 10 |
379 |
342 |
| 15 |
128 |
179 |
| 20 |
61 |
111 |
| 25 |
38 |
85 |
| 30 |
31 |
76 |
Table 7. Average rill spacings (mm) on the lowest third of the slope for the simulations which used the laser-scanned DEM
a |
b |
c |
d |
e |
f |
g |
Figure 13. Total rill-interrill soil loss for the simulations which used the laser-scanned DEM, high-intensity rainfall. Slope angle is (a) 0% (b) 5% (c) 10% (d) 15% (e) 20% (f) 25% (g) 30%. Results for low-intensity rainfall (not shown) are substantially similar |
Figure 13 shows that rill soil loss increases downslope strongly when slope angles are steep, whereas interrill soil loss increases only up to a slope angle of about 15%, and is virtually constant on steeper slope. Results for the high-intensity and low-intensity rainfall differed little. Note that temporary peaks in rill erosion are again visible, but this time slightly downslope of the peaks in rill spacing (about 1300 and 2000 mm downslope distance: cf. Figure 12).
a |
b |
| Figure 14. Time-series of runoff volume and wetted area for the simulations which used the laser-scanned DEM (a) high-intensity rainfall (b) low-intensity rainfall. Each point is itself an average of values from three iterations. Note the different x-axis scales |
Some time-series results are shown in Figure 14, for high-intensity and low-intensity rainfalls. In both cases runoff volume increases approximately linearly at first, then reaches a steady state value (but with wide fluctuations about this). The wetted area (i.e. the number of runoff packets) increases in a more non-linear way at first, reaches a steady-state value at about the same time as runoff volume, and then begins a gradual decline.
a |
b |
| Figure 15. Cumulative runoff and soil loss for the simulations which used the laser-scanned DEM (a) high-intensity rainfall (b) low-intensity rainfall. Note the different x-axis scales |
By contrast, cumulative runoff and soil loss increase linearly throughout. Despite the deliberate similarity between total runoff volumes for high-intensity and low-intensity simulations, total soil loss differed markedly (about 25% lower for the low-intensity runs).
 |
| Figure 16. Measured change in soil surface elevation on the photogrammetrically-derived DEM (from Guerra et al., in preparation). Green represents areas of soil loss, yellow shows areas on increased elevation |
The extent and pattern of simulated ponding on the 'photo' DEM (Figure 6) visually resembled that which was videoed during the field experiment with a portable rainfall simulator (Guerra et al., in preparation). While no quantitative estimate of similarity was possible, this observation encouraged the following analysis, carried out as a preliminary assessment of RillGrow 2's ability to predict spatial patterns of flow.
Patterns of measured change in surface elevation (both increases and decreases) were available from the field experiment (Figure 16). These were, however the result of splash erosion, erosive flow and deposition. Since the model does not yet consider splash erosion or deposition, these patterns could not currently be compared directly against modelled estimates of soil loss. Instead, the assumption was made that where a decrease in elevation had occurred, it was the results of erosive flow rather than splash. This assumption is not unreasonable, since a relatively low intensity was used for the simulated rainfall during the field experiment.
Figure 17 shows the results of plotting simulated water depth against measured decrease in elevation for all cells. Correlation coefficients are given in Table 8.
a |
b |
| Figure 17. Scatterplots of measured decrease in elevation and simulated water depth for the photogrammetrically-derived DEM with (a) high-intensity (b) low-intensity simulated rainfall |
| Simulated rainfall intensity |
R2 |
| High |
0.38 |
| Low |
0.47 |
Table 8. Correlation between measured decrease in elevation and simulated water depth for the photogrammetrically derived DEM
Note that, while all cell values for decrease in elevation are assumed to be independent for this analysis, in reality some spatial correlation will negate this assumption to some extent; the same is true for simulated water depth.
R2 values are not particularly high, however this is a rather difficult test given that each cell is only 2.5 mm in size. Figure 17 shows that in general, cells with high measured decreases in elevation — which are most likely to have experienced erosive flow — do tend to correlate with area of high simulated water depth.
A final analysis investigates the extent to which model results can be considered to be deterministic.
 |
| Figure 18. Correlation between simulated average water depths for each grid cell for pairs of runs with different initial random number seeds. The photogrammetrically-derived DEM was used |
Simulations of varying duration were carried out using the 'photo' DEM: each pair of runs were identical apart from the random number seeds. The spatial patterns of average water depth produced by each pair of runs were then correlated. Figure 18 plots the correlation coefficients against run length. Correlations were zero after 10 iterations (suggesting a wholly stochastic source to the patterns of water depth produced in each run). Correlation however increased quickly with run length, with both pairs of water depth DEMs converging in similarity (suggesting an increasingly deterministic component to the results). After 1000 iterations, the pair are almost identical.
The RillGrow 2 model, even at this early stage in its development, appears not only to be able to replicate the successes of the earlier model, but to exceed them. Planform flow networks appear realistic, with their morphology and density related to slope and roughness in an apparently sensible way. Simulated rill spacing conforms with observed relationships to slope angle; the spatial balance between rill and interrill erosion amounts echoes observational evidence, as does total erosion, variation in rill depth, and time series of runoff and soil loss. Since no explicit separation is made here between rill and interrill zones, it appears that there is an alternative to adoption of the Meyer-Wischmeier (1969) rill-interrill paradigm.
The wave-like patterns on the smooth slopes represent a particularly interesting — and unexpected — model response. The greater definition and stability of patterns produced when 'base' flow velocity differs for each cell (i.e. has a non-zero standard deviation) suggests that some randomness is advantageous in producing this emergent response. The reason for the lack of such patterns at low rainfall intensities is unclear.
Also of interest is the production of both braided and dendritic flow patterns on the laser-scanned DEM. In a cellular automata-based model for braided channel flow, Murray and Paola (1997) found that dendritic flow patterns were produced when their equivalent of 'packets' were routed from cell to cell in their entirety (i.e. with no divergent flow); when 'packets' were split up between adjacent cells in the course of flow routing (i.e. divergent flow occurred), braided patterns formed. It is interesting that no such external control is necessary here to make the switch between braided and dendritic flow patterns: a shift in rainfall intensity is all that is necessary.
Results from the rather crude attempt at a spatial validation undertaken here are also encouraging. Finally, the analysis of stochasticism in the model's output implies that, despite the almost total dominance of random effects in the early stages of a simulation (i.e. the spatial sequence of raindrops upon the grid), determinism 'shines through' and eventually dominates. Since the only source of deterministic spatial information here is the DEM of microtopographic heights, this result implies that inevitable errors (e.g. lack of spatial uniformity for rainfall) may not have too serious an impact on the model's ability to estimate erosion patterns.
Perhaps the most interesting result is that such a simple model is capable of producing such a rich set of responses. This might suggest that, given an appropriate level of systemic description, much of the apparently essential (but intractable!) complexity of present-day operational erosion models can be discarded in constructing a future generation of models. Is there, therefore, an optimum level of simplicity in modelling any natural process — the level at which Nature itself operates? Perhaps this is what Albert Einstein meant by: "A theory should be as simple as possible — but no simpler".
Application of a self-organising approach to modelling rill initiation and development has produced promising results.
We would like to thank Chi-Hua Huang (NSERL) for laser scanning; Tony Guerra (Rio de Janeiro) for use of the results from the field experiment; and Jim Chandler (Loughborough) for carrying out the photogrammetry. This work was funded by the UK Natural Environment Research Council (NERC GR3/11417: Rill initiation by overland flow and its role in modelling soil erosion).
Bak, P. (1997). How Nature Works: the science of self-organized criticality, Oxford University Press, Oxford. 212 pp.
Bryan, R.B. and Poesen, J. (1989). Laboratory experiments on the influence of slope length on runoff, percolation and rill development. Earth Surface Processes and Landforms 14, 211-231.
Capra, F. (1996). The Web of Life, HarperCollins, London, UK. 320 pp.
Casti, J.L. (1994). Complexification: explaining a paradoxical world through the science of surprise, Harper Collins, New York, USA. 320 pp.
Cohen, J. and Stewart, I. (1994). The Collapse of Chaos, Penguin, Harmondsworth, Middlesex, UK. 495 pp.
Coulthard, T.J., Kirkby, M.J. and Macklin, M. (1996). A cellular automaton fluvial and slope model of landscape evolution. In, Abrahart, R.J. (ed.), Proceedings of the First International Conference on GeoComputation (Volume 1), School of Geography, University of Leeds. pp. 168-185.
Coveney, P. and Highfield, R. (1995). Frontiers of Complexity, Faber and Faber, London, UK. 462 pp.
De Ploey, J. (1989). A model for headcut retreat in rills and gullies. Catena Supplement 14, 81-86.
Eastman, J.R. (1997). IDRISI for Windows User’s Manual v2.0. Clark University, Worcester MA, USA
Evans, R. (1995). Some methods of directly assessing water erosion of cultivated land – a comparison of measurements made on plots and in fields. Progress in Physical Geography 19(1), 115-129.
Favis-Mortlock, D.T. (1996). An evolutionary approach to the simulation of rill initiation and development. In, Abrahart, R.J. (ed.), Proceedings of the First International Conference on GeoComputation (Volume 1), School of Geography, University of Leeds, UK. pp. 248-281.
Favis-Mortlock, D.T. (1998a). Validation of field-scale soil erosion models using common datasets. In, Boardman, J. and Favis-Mortlock, D.T. (eds), Modelling Soil Erosion by Water, Springer-Verlag NATO-ASI Series I-55, Berlin. pp. 89-128.
Favis-Mortlock, D.T. (1998b). A self-organising dynamic systems approach to the simulation of rill initiation and development on hillslopes. Computers and Geosciences 24(4), 353-372.
Favis-Mortlock, D.T., Guerra, A.J.T. and Boardman, J. (1998). A self-organising dynamic systems approach to hillslope rill initiation and growth: model development and validation. In, Summer, W., Klaghofer, E. and Zhang, W. (eds), Modelling Soil Erosion, Sediment Transport and Closely Related Hydrological Processes, IAHS Publication No. 249, Wallingford, UK. pp. 53-61
Favis-Mortlock, D.T., Guerra, A.J.T., Chandler, J.H., Boardman, J. and Todd, M.C. (in preparation). Initial validation of a self-organising dynamic systems approach to erosion modelling, using a rainfall simulator and digital photogrammetry.
Gilley, J.E., Finckner, S.C., Nearing, M.A. and Lane, L.J. (1989). Hydraulics of overland flow. In, Lane, L.J. and Nearing, M.A. (eds), USDA - Water Erosion Prediction Project: Hillslope Model Documentation, USDA-ARS National Soil Erosion Research Laboratory Report No. 2, NSERL, West Lafayette, Indiana, USA. pp. 9.1-9.4.
Govers, G. (1991). Rill erosion on arable land in central Belgium: rates, controls and predictability. Catena 18, 133-155.
Govers, G. and Poesen, J. (1988). Assessment of the interrill and rill contributions to total soil loss from an upland field plot. Geomorphology 1, 343-354.
Green, D.G. (1993). Emergent behaviour in biological systems. In, Green, D.G. and Bossonmaier, T.J. (eds), Complex Systems - From Biology to Computation, IOS Press, Amsterdam. pp. 24-35.
Guerra, A.J.T., Favis-Mortlock, D.T. Chandler, J.H., Gallart, F., Boardman, J. and Todd, M.C. (in preparation). An investigation of rill initiation, runoff and soil loss, using a field rainfall simulator and digital photogrammetry
Hergarten, S. (1997). A Physical Stochastical Model for Investigating Pattern Formation and Self-Organized Criticality in Erosion Processes, PhD Thesis, Rheinische Friedrich-Wilhelms-Universität Bonn, Germany. 82 pp.
Horton, R.E. (1945). Erosional development of streams and their drainage basins: a hydrophysical approach to quantitative morphology. Bulletin of the Geological Society of America 56, 275-370.
Huang, C.-H. and Bradford, J.M. (1993). Analyses of slope and runoff factors based on the WEPP erosion model. Soil Science Society of America Journal 57(5), 1176-1183.
Lane, L.J. and Nearing, M.A. (eds) (1989). USDA - Water Erosion Prediction Project: Hillslope Model, USDA-ARS National Soil Erosion Research Laboratory Report No. 2, NSERL, West Lafayette, Indiana, USA.
Lenton, T.M. (1998). Gaia and natural selection. Nature 394, 439-447.
Meyer, L.D. and Wischmeier, W.H. (1969). Mathematical simulation of the process of soil erosion by water. Transactions of the American Society of Agricultural Engineers 12, 754-758.
Morgan, R.P.C., Quinton, J.N., Smith, R.E., Govers, G., Poesen, J.W.A., Auerswald, K., Chisci, G., Torri, D. and Styczen, M.E. (1998). The European Soil Erosion Model (EUROSEM): a dynamic approach for predicting sediment transport from fields and small catchments. Earth Surface Processes and Landforms 23(6), 527-544.
Murray, A.B. and Paola, C. (1997). Properties of a cellular braided-stream model. Earth Surface Processes and Landforms 22(11), 1001-1025.
Nearing, M.A., Norton, L.D., Bulgakov, D.A., Larionov, G.A., West, L.T. and Dontsova, K. (1997). Hydraulics and erosion in eroding rills. Water Resources Research 33(4), 865-876.
Parsons, A.J. (1987). The role of slope and sediment characteristics in the initiation and development of rills. In, CNRS (eds), Processus et Mesure de l'Érosion, 211-220.
Perez-Trejo, F. (1993). Landscape response units: process-based self-organising systems. In, Haines-Young, R., Green, D.R. and Cousins, S. (eds), Landscape Ecology and Geographic Information Systems, Taylor and Francis, London, UK. pp. 87-98.
Phillips, J.D. (1996). Deterministic complexity, explanation, and predictability in geomorphic systems. In, Rhoads, B.L. and Thorn, C.F. (eds), The Scientific Nature of Geomorphology, Wiley, New York. pp. 315-335.
Phillips, J.D. (1997). Simplexity and the reinvention of equifinality. Geographical Analysis 29(1), 1-15.
Ruelle, D. (1993). Chance and Chaos, Penguin, London, UK. 195 pp.
Scheidegger, A.E. (1994). Hazards: singularities in geomorphic systems. Geomorphology 10, 19-25.
Schumm, S.A., Harvey, M.D. and Watson, C.C. (1984). Incised Channels: Morphology, Dynamics and Control, Water Resources Publications, Littleton, Colorado, USA. 200 pp.
Slattery, M.C. and Bryan, R.B. (1992). Hydraulic conditions for rill incision under simulated rainfall: a laboratory experiment. Earth Surface Processes and Landforms 17, 127-146.
Thorn, C.E. and Welford, M.R. (1994). The equilibrium concept in Geomorphology. Annals of the Association of American Geographers 84(4), 666-696.
Wischmeier, W.H. and Smith, D.D. (1978). Predicting Rainfall Erosion Losses, US Department of Agriculture, Agricultural Research Service Handbook 537, 58 pp.
Wolfram, S. (1982). Cellular Automata as Simple Self-Organising Systems, Caltech preprint CALT-68-938.
Zobeck, T.M. and Onstad, C.A. (1987). Tillage and rainfall effects on random roughness: a review. Soil and Tillage Research 9, 1-20.