Some factors controlling gully growth in fine-grained sediments: a model applied in southeast Spain

Catena 40 Ž2000. 127–146 www.elsevier.comrlocatercatena

Some factors controlling gully growth in fine-grained sediments: a model applied in southeast Spain M.J. Kirkby ) , L.J. Bull

1

School of Geography, UniÕersity of Leeds, Leeds, LS2 9JT, UK Received 23 March 1998; received in revised form 22 July 1999; accepted 25 August 1999

Abstract Gullies in the Guadalentin catchment, SE Spain, are formed in a variety of lithologies, ranging from schists to marls. There are clear differences in valley and channel morphology between these two types. Gully erosion in schists produces large amounts of clastic material that restricts sediment transport, leading to a relatively smooth transition from gullies to hillslopes, and to broad alluvial valley bottoms with laterally migrating channels. Gullies in marls contain only limited gravels, so gully extension leads to sharply incised gully heads. A two-dimensional model has been developed to investigate and simulate the distribution of gully heads in the Guadalentin. The model runs on a smooth initial surface with an added fractal perturbation. Without some perturbation of the initial surface or material properties, symmetry prevents the formation of any channels. Fractal perturbation was used here to provide initial irregularities at all scales. This paper explores the dependence of model response on the form of the sediment transport law. Initially, sediment transport was considered to be transport-limited and so depends on the relationship between wash and creeprsplash. This was modified by Ži. altering the exponents for gradient and area to conform with experimental results, Žii. adding a critical tractive power for the initiation of rillwash, Žiii. including a more explicit treatment of mass movements, Živ. including terms to reduce the dependence of the model on the computational grid size, Žv. providing a basis for explicitly integrating from individual storms to average long-term behaviour and, most critically Žvi. introducing the ‘effective bedload fraction’ Žebf. as a simplifying concept to allow for the selective transportation of eroded sediment. This modification allows the model to represent removal of fines in marl-rich catchments, without the need to introduce a fully explicit travel distance formulation. The steep gully heads and vertical side walls developed in marl catchments in the Guadalentin basin are qualitatively well-represented by the model. Preliminary

) 1

Corresponding author. Tel.: q44-113-233-3300; fax: q44-113-233-3308; e-mail: [email protected] E-mail: [email protected].

0341-8162r00r$ - see front matter q 2000 Elsevier Science B.V. All rights reserved. PII: S 0 3 4 1 - 8 1 6 2 Ž 9 9 . 0 0 0 7 7 - 6

128

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

tests showed that steep-walled gullies are generated where the ebf falls to below about 0.5. If a tractive power threshold is also introduced, gentler headcuts are produced. Selective transportation is thus seen as the critical determinant of steep headcuts in gully systems. q 2000 Elsevier Science B.V. All rights reserved. Keywords: Gully; Erosion; Two-dimensional model; Sediment transport

1. Introduction Development of channel and gully heads is important in understanding and predicting landscape evolution, morphology and response to climatic change ŽDietrich and Dunne, 1993.. Modelling the development of rilling and gullying began in earnest in the 1970s with most models in this period concentrating on the stochastic reproduction of the badlands landscape Že.g., Scheideger, 1970; Howard, 1971; Smart and Moruzzi, 1971a,b; Freeze, 1972; Ahnert, 1976.. The 1980s saw the introduction of models based on process laws in an attempt to understand the theories behind gully initiation ŽCordova et al., 1983; Hugus and Mark, 1984; Kirkby, 1986, 1987; Howard, 1988; Howard and MacLane, 1988.. More recently, the 1990s have then seen the development of models that combine processes occurring on hillslopes and within channels to reproduce gully development ŽKirkby, 1991, 1992, 1993, 1994, 1995; Willgoose et al., 1991a,b,c, 1994; Howard, 1994, 1997.. A more complete review of these processes and models is included in the work of Bull and Kirkby Ž1997.. The models of Willgoose et al. Ž1991a. and Howard Ž1994. illustrate recent attempts at high resolution process-based simulation of the slope and channel development at the basin scale. The model developed by Willgoose et al. Ž1991a; b; c; 1994. is a large scale representation of catchment evolution involving channel network growth and elevation evolution. It assumes that overland flow, as well as flow in channels is transport-limited. It brings together a model of erosion processes that has been theoretically and experimentally verified at small scales and a physically based conceptualisation of channel growth processes ŽWillgoose et al., 1991a.. Willgoose et al. Ž1991a. determine elevation changes by means of a continuity equation for sediment transport. Network growth is controlled by the spatial distribution of a channel initiation function. The channel initiation function is non-linearly dependent on the contributing area, average runoff rate and local hillslope gradient. Mass transport processes considered are fluvial sediment transport and mass movement Žincluding creep, rainsplash and landslides.. A point is defined as being in the channel when selected flow and transport processes exceed a threshold value. Channel growth is governed by hillslope form and processes that occur upslope of the channel head by their ‘channel initiation function’. If the channel head function exceeds a threshold value at some point, the channel head advances to that point. The channel head initiation function is dependent on discharge and gradient, and the threshold is dependent on the resistance of the catchment to channelisation. Crucially, the model explicitly incorporates the interaction between hillslopes and the growing channel network on the basis of observed mechanisms. However, one major problem here is that channel heads can only advance and not retreat.

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

129

The drainage basin model of Howard Ž1994. combines diffusive Žmass wasting and rainsplash. plus advective Žfluvial erosion. processes. This model differs from previous advection–dispersion models ŽKirkby, 1971, 1986; Ahnert, 1976, 1987; Willgoose et al., 1991a,b. by assuming that headwater channels are detachment-limited. Potential erosion or deposition due to diffusive processes is given by the spatial divergence of the vector flux of regolith movement. The rate of movement is expressed as two terms, one for creep andror splash and one for mass movement. Fluvial erosion is assumed to consist of two processes. For the first, erosion is considered to be detachment-limited in steep channels flowing on bedrock or regolith in which the bedload sediment flux is less than a capacity load. The second is applied to low gradient alluvial channels where erosion is transport-limited, and the erosion rate due to detachment is assumed to be proportional to the shear stress exerted by the dominant discharge. Processes are considered to be continuous in time, representing an implicit integration over the frequency distribution. One advantage of this model is thus that both fluvial and slope processes can occur in each cell and that the location of channel heads is defined by a morphometric criterion. This paper presents a two-dimensional simulation model of landform evolution by wash and creeprsplash. The same process law is applied to all points in the model space to allow the development of a valley network with minimal constraints in order to represent erosion on a surface of constant composition. This model follows the approach taken by Armstrong Ž1976., Dunne and Aubry Ž1986. and Kirkby Ž1986. and concentrates on gully incision into previously unchanneled hillslopes. This model is used to investigate the way in which variations in the governing sediment transport laws influence landscape evolution. A two-dimensional simulation model is used because it has the advantage of simplicity, yet still provides useful information regarding the way in which hillslope and channel processes interact to form the landscape and set its characteristic scale.

2. Formulation of the gully simulation model The model focuses on the conditions required to initiate gully heads. In real landscapes, these occur at two distinct scales, as rills or ephemeral gullies, and as the heads of major valleys. However, the two types grade into one another in some locations, for example, where ephemeral gullies converge. Conceptual approaches in this model are based on work by Smith and Bretherton Ž1972., Willgoose et al. Ž1991a; b; c., Dietrich and Dunne Ž1993., and Kirkby Ž1994., in which they explore the balance between enlargement and infilling in the neighbourhood of the channel. For single storms, sediment transport is assumed to be related to the storm runoff intensity. The formation of both rills and ephemeral gullies is seen as a response to individual storms in which the runoff intensity is able to exceed the critical stress for grain entertainment or vegetation mat incision. Over longer periods, many of these small channel features are assumed to be refilled by tillage and splash in smaller storms, so that they do not become a part of the permanent gully network. The governing equations can be written in the following forms ŽKirkby, 1992.. For a flow strip of variable width, w the following statements below apply.

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

130

Ž1. The continuity equation allows for the conservation of mass. The expression is in volume terms, so that elevation should be corrected to rock-equivalent units. Ez 1 E Ž Sw . q s0 Ž 1. Et w Ex where z is elevation, w is flow strip width, S is actual sediment transport Žgenerally at less than the capacity, C ., and x is horizontal distance measured along the flow strip. Ž2. The sedimentation balance equation gives the total rate of sediment pick-up from its rate of detachment and redeposition. If the material is consolidated, then the detachment rate is reduced by a fraction a - 1. The difference between total and partial differentials on the left-hand side is only important in transient flows, and is not considered to be relevant in this application. 1 d Ž Sw . S sa Dy Ž 2. w dx h where D is the rate of sediment detachment and h is the sediment travel distance. Ž3. For wash and splash, this equation has been simplified using the concept of the Effective Bedload Fraction Žebf.. This is formally defined as:

ž

ÝSrh ebj s

d

ÝCrh d

/

ÝSrh s

d

Ž 3.

ÝD d

where the summations are taken over each grain-size class and C is the transporting capacity. For fines, S < C and the contribution to the ebf is small, while for coarser material for which movement is transport-limited, the contribution to the ebf is almost total. Ž4. The sediment transport law for rainsplash, rainflow and rill erosion is taken from the family of empirical expressions for sediment transport in a storm: c

C s k 1 r 2L a q k 2 L b Ž q L y Q . Ž 4. where C s Dh is the sediment transporting capacity, r is storm rainfall, q is overland flow storm discharge per unit width Žm2 sy1 ., L is the local gradient, Q is the flow power threshold for entrainment Žm2 sy1 . and k 1 , k 2 , a, b, and c are constants. The first term represents rainsplash, and the second channelled wash Žtermed rillwash below.. Ž5. The sediment transport law for mass movements. These are treated as a continuous process ŽKirkby, 1984., following Eqs. Ž1. and Ž2. above, but with a separate sediment budget, and with different rates of detachment and travel distance, as follows: D M s D 0 L Ž L y L0 . where S ) 0 Ž 5. h M s h 0r Ž L T y L . where L - L T . 3. The ebf In this formulation, the second equation is a sedimentation balance, and the transporting capacity, C s Dh, is not generally equal to the sediment transport rate, S. Over the

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

131

range of grain sizes, the sediment travel distance increases as size decreases, so that transport is effectively flux-limited for coarse debris, and supply-limited for fine material. The response and evolution of the coarse- and fine-grained rambla systems studied in SE Spain is thus very different. For coarse-grained systems, such as Rambla Nogalte, sediment transport is largely flux-limited, and headwaters remain well-connected with their source areas. For fine-grained systems, gully enlargement produces sediment that travels far from its source, and shows a delicate dependence on hydraulic and hydrological conditions. These cases are usefully distinguished by the ebf, without requiring a full analysis of the selective transportation. For the simplest case, with a uniform sediment transport environment Žconstant D and h above. downstream, the actual sediment transport downstream follows the approximate form: S s Dh 1 y exp Ž yxrh .

Ž 6.

This expression shows sediment transport rising to its capacity over a distance scaled to the travel distance, h. Using this form, ebf may be estimated in the form: ebf s

Ý Di 1 y exp Ž yxrh i . Ý Di

.

Ž 7.

This provides an explicit, though approximate, expression for the change in ebf downstream in relation to the travel distances of the various grain-size classes, and their proportions in the source material Žthe Di ., assuming equal mobility for particle entrainment. Fig. 1 shows the estimated ebf weighings for different grain sizes for the relationship in Eq. Ž7.. It may be seen that the ebf is increasingly dominated by the coarsest fractions in the downstream direction.

Fig. 1. The ebf weighings downstream for each grain-size class.

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

132

Fig. 2. Estimated ebf changes downstream for different grain-size mixtures.

Combining these estimates for log-normally distributed grain-size mixtures, Fig. 2 shows the total effect on the estimated ebf. The distance units in these figures are arbitrary, but are estimated to be of the order of 1 m, so that the rate of increase in ebf declines steadily downstream in the zone near to the channel head. It is for this reason that the use of an approximately constant value of the ebf is considered to provide an effective basis for a simplified operational model, and contains most of the essential elements of a model with full grain-size disaggregation. In field use, the definitions of ebf in Eq. Ž3. or Eq. Ž7. are not easy to apply, but a good estimate of the ebf may be obtained by comparing the grain-size distributions in the channel bed with that of the source deposits from which the bed material is derived. An enrichment factor ŽEF. may be calculated for each grain-size class as: EFi s

Frequency of ith size class in channel bed material Frequency of ith size class in source material

Ž 8.

The ebf may then be estimated as the reciprocal of the maximum EF.

4. Conditions for stream head formation due to wash processes The general form of wash transport for a single storm of r generating discharge per unit width q, has been taken in Eq. Ž4. above as: C s Srebf s k 1 r 2L a q k 2 L b Ž q L y Q .

c

Ž 9.

where q s aŽ r y t . for a really uniform runoff above threshold t. The exponents a, b, and c are poorly constrained by data on direct process measurements, but may be better

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

133

constrained by experimental data on stream head position. Limiting stability considerations are also relevant, in that material can be moved by very low hydraulic forces as the threshold angle is approached. Thus, the gradient terms should properly be replaced by LrŽ1 y LrL T ., with a cut-off to prevent infinite values above gradient L T . Experimenting with different forms, the preferred final form is with exponents a s c s 1; b s 0 above. With this form, the critical condition for unstable hollow enlargement is obtained from the Smith and Bretherton Ž1972. criterion, aŽŽES .rŽEa.. ) S, as follows. The critical threshold gradient is given by the condition of equality, aŽŽES .rŽEa.. ) S, so that: 2

S s k 1 r 2L q k 2 Ž q L y Q . s k 1 r 2L q k 2 a Ž r y t . L y Q

2

Ž 10 .

ES a

Ea

s 2 k 2 aŽ r y t . L aŽ r y t . L y Q

At equality,

L

k1 r 2

k 1 r 2rk 2 q

k2

2

L s a Ž r y t . L y Q 2 or

(Ž k r rk . q 4Ž Q aŽ r y t . . 2

2

1

2

2

2 Ž aŽ r y t . .

2

The final expression of Eq. Ž10. is illustrated in Fig. 3, with and without the correction for threshold gradient L T . It may be seen that with the threshold gradient, the form of this relationship has the same general form of those found in field work by Poesen et al. Ž1997. and Vandekerckhove et al. Žsubmitted. Ž L ( Ay0 .2 ., while also

Fig. 3. Critical gradient as a function of downslope distance for single storms.

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

134

being consistent with the work of Montgomery and Dietrich Ž1989. Ž L ( Ay1 . for angles well below the threshold. The effect of different sized storms andror different runoff thresholds is to shift the whole diagram to the left or right, but not to change its overall form. Thus, for large storms channel enlargement occurs far upstream, whereas in smaller storms Ž r smaller. or under drier conditions Ž t larger., the condition is met farther downstream while wash processes produce channel infilling upstream. There is a second condition to be met, namely that the runoff is great enough to exceed the flow power threshold, Q . This may be written as: aŽ r y t . L ) Q

Ž 11 .

It may be seen immediately, from the third expression in Eq. Ž10., that this condition is met for all positive solutions for the critical gradient, so that this condition is never limiting. Significant change in channel heads takes place episodically, with major change in infrequent large storms interspersed with long periods of quiescence or more passive modification of the channel forms. For these long-term changes, two approaches are possible. The process may be modelled as a continuous process, integrating implicitly or explicitly over the frequency distribution of storms; or by treating each complete storm as a single iteration. The former approach emphasises the effects of progressive climate or land use change; the latter emphasises the structure of the frequency distribution and the effect of particular histories of storm sequences. It should be noted that Eq. Ž10. above is for single storms, and should be aggregated over a series of storms, weighed by their frequencies, to give a long-term evolution model. This approach follows the method by Kirkby Ž1993. Žpp. 267–289., although the equations in each case are different due to the different transport law exponents used here. For exponentially distributed storms, the number density of storms with a storm Žor daily. rainfall of r is: N

nŽ r . s

r0

r

ž /

exp y

Ž 12 .

r0

where N is the number of rain days and r 0 is the mean rain per rain day. Summing the above expression over the spectrum of storms, the total sediment transport is: `

S s k1 L

H0 r

2

`

nŽ r . d r q k2

H1

X

aŽ r y t . L y Q

2

nŽ r . d r

Ž 13 .

where hX s h q ŽŽQ .rŽ a L... The range of integration is set by the lower limit of sediment producing storms. Storms in excess of t produce runoff, but no sediment until tX is reached. Completing the integration: Ssk Lq

aL

ž / u

2

exp Ž y2F .

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

135

where k1

t

( ž /

u

k2

exp

2 r0

and

Q Fs

Ž 14 .

2 a L r0

In this expression, u is the constant distance at which rillwash begins to dominate, for unit gradient and zero threshold. F is a dimensionless variable which contains both area, a, and gradient, L. The average position of the stream head instability is again given by the point at which aŽŽES .rŽEa.. s S. For the transport law above, the critical gradient and area are linked by the implicit relationship:

L

a

ž / u

2

s Ž 1 q 2F .

y1

exp Ž 2F .

and with more general exponents a, b, and c,

LŽ bqcya.

a

ž / u

c

s Ž c y 1 q 2F .

y1

exp Ž 2F .

Ž 15 .

These expressions define a unique relationship between critical gradient and area. If gradient is corrected for a threshold stable angle, the long-term relationship between them is qualitatively very similar to that shown above in Fig. 3. As in this figure, the effect of a larger flow power threshold is to reduce the exponent of distance or area at low gradients. Thus, both the behaviour in single storms ŽFig. 3. and the cumulative effect integrated over the frequency distribution produce morphological responses which are very similar. Comparison of the form of the relationship in Fig. 3 with the results of Vandekerckhove et al. Žibid. also provides an important and sensitive confirmation of the general form of the sediment transport equation adopted. Comparing average stream head position with the position in individual storms, a critical storm size may be defined for which the position is similar to the average position. Where the mean rain per rainday, r 0 , is much less than the storage capacity, t, the critical storm, rc is only slightly above the storage capacity. Where r 0 ) 1r2 t a better estimate is rc s 1r2 t q 2 r 0 . In general, trrc s 1 y expwytrŽ2 r 0 .x.

5. Dependence on grid scale The problem addressed in this section is that changes in grid size can make a substantial difference to the way in which the model runs. This is attributed to the fact that unit area a Žarea drained per unit contour length. is calculated in the model with reference to the cell width, L. Where flow is unconcentrated, this has little impact, but

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

136

where there is an implicit channel Žremembering that there are no defined channels as such in the model formulation. of width w - L, then the flow will concentrate into this. Thus, for an area A Žmade up of cells drained by whatever appropriate algorithm., the simplest model computes a s ArL, giving the fluvial sediment transport and net denudation scaled to Ž aru. 2 s Ž ArLu. 2 . A neighbouring cell, with a small area, will add little to the sediment transport, so that the average, over the two neighbouring cells, is 1r2Ž ArLu. 2 . For a recomputation with doubled cell size, the same total is calculated as Ž ArLu. 2 , which is only half the first value. Thus, the coarser grid produces less intense erosion along channelways, although the same on unchannelled slopes. There are two possible approaches to dealing with this, which should preferably produce similar results. In the first, a correction is made for estimated channel width, from the stream head down to the point at which the channel width is greater than the cell width. In the second, the correction is made via an estimate of the random roughness, which determines the degree of flow concentration. In both cases, the model must be given a dynamic for changing width or roughness over time. Taking a channel width approach based on hydraulic geometry, channel width behaves roughly as kA0.5 for k ; 0.001, although this is only valid once channels are clearly established beyond a critical threshold, say A c , associated with critical width wc , and scaled in the model to the distance u at which rillwash becomes dominant. An algebraic form for the effective width that is dimensionally consistent and avoids discontinuities is: wc s

L L q kA1r2

ž

LA1r2 c A1r2

/

q kA1r2 .

Ž 16 .

Below the critical value of A s A c , this expression is truncated at the cell width L, corresponding to unchannelled hillslopes. Above this value, the effective width falls to a minimum, and then rises as A0.5, before levelling out again at the cell width, corresponding to channels which are wide relative to the cell width, for which the second term in the final bracket is dominant. It is assumed that within the cell, the channel is randomly located relative to the cell boundaries. Empirically, it has been found that the response for different cell sizes is most uniform if the critical area, A c , takes the form A c s 10y4 u 3 Ly2 where u, the distance at which rillwash becomes the dominant process, L is in metres, and A c in m2 . The effect of using this expression is shown in Fig. 4. The reduction in effective width is responsible for a substantial increase in sediment transport along channelways through the greater concentration of the flow. It may be seen that the effective width is largely independent of the number of cells across the 1000 = 1000 m2 area. A roughness approach can also be based on hydraulic geometry. RMS roughness can be calculated for a rectangular channel of width d and width w within an otherwise undissected grid cell of side L. Calculation of RMS roughness j also gives zero above the stream head, rises to a maximum and then falls progressively, in a way that is qualitatively similar to that shown in Fig. 1. If the channelling is reinterpreted as roughness across the whole cell, the roughness can be used to create a critical unit area, A c s LjÕrI, where Õ is the overland flow velocity and I is the runoff intensity, so that

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

137

Fig. 4. Effective channel width from Eq. Ž16. in terms of drainage area and grid resolution.

the sediment transport behaves, over the rough surface as AŽ A q A c .rL2 , in place of the a 2 term above. Equating these two expressions, the roughness is estimated from: 2

Ac

Lqw s wc A L qw w



0

y 1 for w G wc , and zero otherwise.

Ž 17 .

For the growth of roughness, the equations used in the MEDRUSH model ŽKirkby et al., 1998. form a starting point. dh dt

s V y m j Ž 1 q zrz c . q l z

Ž 18 .

where j is the roughness and the final term is associated with channel erosion, with l being proportional to the local rate of denudation Žyve for deposition.. In the current form of the model, no explicit expression for roughness growth is included, since the three-dimensional dynamics of the model already provide for unstable hollow growth and infilling. The channel width function has, however, been explicitly included.

6. Computational procedure The computational procedure is outlined in Fig. 5. The initial surface is a pseudofractal surface obtained by successive mid-point interpolated values. This surface is

138

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

Fig. 5. Outline flow diagram for the two-dimensional simulation model.

generated within a square frame at elevation zero; beginning with a centre-point height, an initial standard deviation, and an excess fractal dimension. Successive sets of four points are obtained by linear interpolation from the four corners of the square in which they lie. A normally distributed random variate, with a mean of zero and a standard deviation which reduces with the size of the square and the excess dimension, is added to the interpolated value. This initial surface is then linearly combined with a rectilinear surface that consists of a uniform slope, falling to zero at the base. Each point in the grid is considered to be the

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

139

mid-point of a square cell. The upper edge of the grid is taken to be the divide, with symmetrical removal of sediment on the opposing slope. The slope base is considered to be a base level of fixed position and zero elevation, at which all material is removed. The two sides of the slope are allowed to connect to one another to create a periodic boundary condition. This initial surface is scanned to infill any closed depressions before beginning the simulation. At each model time step, the elevation of each point is compared to that of all its immediate neighbours, and the direction of flow is assigned to the lowest of these points. All flow is then assumed to follow the path to the lowest of these points. Multiple flow path algorithms have also been used, but not in the runs presented here. Drainage area at each point in the grid is evaluated starting at each element, and adding unity to the area of this element and to all elements to which it successively drains down the length of the slope. Sediment transport is then evaluated from each grid cell using Eqs. Ž2. – Ž5.. Each element provides a rate of erosion for the element itself and a rate of deposition for the element to which it drains. The solution procedure is stabilised by choosing a time step that allows the gradient to change by a maximum of 20% at the worst point. This condition, which ensures computational stability, is not applied where points differ by less than 4 mm, so that reversals of slope through divide migration can occur. Elevations are finally updated and the next time step begun. At chosen intervals, the data matrices are displayed andror stored for analysis of erosion rates, drainage areas and network patterns.

7. Model results This type of simulation model has much in common with other two-dimensional simulations models ŽKirkby, 1994; Willgoose et al., 1994; Howard, 1997, etc.., and is able to generate simulated landscapes which share many qualitative features of real landscapes. Such landscapes develop valley networks with plausible drainage densities, long profiles and hillslope profiles, even though the present model formulation does not explicitly make a demarcation between ‘hillslopes’ and ‘channels’. The qualitative factors which are explored here are the cross profiles of mature gullies and their evolution; and the morphology and extension of valley heads. Preliminary tests show that steep-walled gullies are generated as the ebf falls below about 0.5, and that the introduction of a non-zero tractive power threshold appears to create less steep headcuts. Fig. 6 gives two examples, both with a zero tractive power threshold, but with values for ebf of 1.0 and 0.1, showing the effect of this term on the morphology of gully growth. It may be seen that the reduction of ebf allows much more sediment to be evacuated, but also dramatically changes the valley form. In Fig. 7, the long profile forms are shown for the two cases shown in Fig. 6. With a lower ebf, the valley floor is flatter, whereas the valley head remains steep, which helps to maintain its continued back-cutting. The profile form for ebf s 0.1 is also shown after 100 years, when the valley head is close to that for ebf s 1.0 after 1000 years. It may be seen that flat gully floors and steep heads are developed throughout the length of the gully. In

140 M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

Fig. 6. Comparison of simulated gullies. ŽLeft. ebf s 0.1 after 1000 years; ŽRight. ebf s1.0 after 10 000 years. Other parameters and initial surface identical.

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

141

Fig. 7. Profiles of the simulated gullies shown in Fig. 6. Gullies with low ebf have markedly lower valley floor slopes and sharp valley heads. Cross profile is shown with expanded horizontal scale. Note shallow fill along valley axis.

transverse profile, the valley sides are held close to the limiting angle for mass movements Ž20% in this example., and there is a shallow Ž0.75 m. fill along the valley axis. Gullies with low ebf have markedly lower valley floor slopes and sharp valley heads.

8. Comparison of simulated and real landscapes Here, the model is compared with badland catchments in southeast Spain. Validation of models of drainage basin development is rendered more difficult by limited opportunity to compare predicted with observed landform evolution, however, surveying and digital flights within the study catchments allow comparison of gully morphology and drainage densities. Three catchments are being used to validate the model: the Rambla Nogalte, Rambla Salada and Rambla de Torrealvilla ŽFig. 8.. Rambla Nogalte has an area of 190 km2 , and drains an area of mainly metamorphic rocks and conglomerates derived from them. Its headwaters are well-connected with the highland areas, sediment transport is largely transport-limited because there are significant amounts of gravel throughout the system. The ebf is estimated to be about 0.5 from comparison of gully bed and side wall composition. Dominant processes operating are creeprwash and rilling on the hillslopes of the sub-catchments. These processes move sediment into the main channel to await transport during large flood events. Sediment

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

Fig. 8. Catchment locations in southern Spain.

142

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

143

awaiting transport from the tributary catchments tends to clog up tributaries and disconnects the sub-catchments from the main channel. The development of steep-sided arroyo type channels is limited. Rambla de Torrealvilla drains 113 km2 of varied lithologies, including large areas of marls and other sedimentary rocks. It is characterised by large areas of undissected pediments, within which there are well-developed box-shaped channels Žarroyos. with steep sides and active lateral and headward retreat. Gully beds have little gravel. It is clear that gully expansion and refilling have occurred episodically on more than one occasion, and it is surmised that the entire pediment surface reflects periods of mainly depositional activity. Small upland areas of resistant rocks persist, but are generally not well-connected to the arroyo systems. Dominant processes operating at the channel heads are rainsplash and mass failure. Splash operates on near vertical slopes and flat areas between gully heads to move fines downslope. In some instances, gravel sheets cover the surface area between gully heads, protecting the marls and limiting erosion from rainsplash. Mass failures at the gully heads are driven by the development of tension cracks that are gradually enlarged to produce slab failures of the vertical faces. This sediment then awaits transport at the base of steep gully walls. The ebf in this system is estimated as very low Ž- 0.01., and the system is dominated by detachment processes. Rambla Salada drains 132 km2 of mainly marls, although with some inter-beds of conglomerate. Arroyo forms are well-developed, with some gravel where coarse beds outcrop. The ebf is estimated at 0.05–0.2. This catchment differs from the Torrealvilla in that relatively little undissected pediment survives between the arroyo systems. The dominant processes operating include mass failure and wash processes. Mass failures are predominantly driven by the development of pipes and their subsequent collapse, although the development of tension cracks is also evident at some channel heads. Piping may be partly responsible for producing the more heavily dissected landscape. Model outputs simulate the landscape morphology within the study catchments well. Fig. 9 shows a digital image with 40 cm pixels of box like arroyos in the Rambla de Torrealvilla. The image shows steep-sided ephemeral channels with sharp heads, the valley floors of which have been made into agricultural terraces. The plan form of these channels is similar to the output from the model illustrated in Fig. 6. The Rambla de Torrealvilla has a low ebf for which the model predicts high sediment transport, low gradient valley floors and steep gully heads and side walls ŽFig. 7., supported by field observations. In contrast, simulations for Rambla Nogalte have been assigned a high ebf, resulting in lower rates of sediment transport and less abrupt distinctions between the ephemeral channels and contributing hillslopes. This is again supported by field observations. The morphological impact of a storm of 70–150 mm Žaccording to site. rainfall in September 1997 has been surveyed in the three catchments. It appears to have generated substantial flows in some of the small gullies Žup to 100 mmrh peak runoff rate in one case. but only to have caused minor and localised gully extension Ž0–2 m.. The duration of intense flow in major streams also appears not to have provided a connected throughput of water and sediment as an integrated system ŽShannon and Richardson, personal communication.. This storm, which is provisionally estimated to have a

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

144

Fig. 9. Digital image of arroyos in the Rambla De Torrealvilla, SE Spain.

recurrence interval of 10–20 years, is therefore considered to lie close to the threshold for major gully extension.

9. Conclusions The model gives qualitatively plausible forms. The next stage of development is to use it to explore the sensitivity of the three catchment systems to changes over time in two distinct senses. First, the frequency distribution of the existing rainfall regime has been analysed to estimate the magnitude, frequency and cumulative impact of formative floods. Second climatic scenarios for the future or recent past can be similarly analysed to indicate expected changes in relation to past land use history and future policy. Finally, and taking a longer time perspective of thousands rather than tens of years, the next direction for modelling may be based around the debate concerning the age of badlands and rates of gully development. Although badlands are generally assumed to be dynamic with high rates of erosion, this may not be the norm, and models should also be able to incorporate periods of relative stability. Modelling may also provide an intermediate approach to understanding channel head location by using simple relationships that can be predicted theoretically. However, in the long run transport laws used in more complicated models need testing with field data to validate results. This also highlights the need for a comprehensive review of research to enable fieldwork and models to be developed in conjunction with each other to improve our understanding.

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

145

References Ahnert, F., 1976. Brief description of a comprehensive three dimensional process-response model of landform development. Zeitschrift fur ¨ Geomorphologie Supplement Band 25, 29–49. Ahnert, F., 1987. Process-response models of denudation at different spatial scales. Catena Supplement 10, 31–50. Armstrong, A.C., 1976. A three-dimensional simulation of slope forms. Zeitschrift fur ¨ Geomorphologie Supplement Band 25, 20–28. Bull, L.J., Kirkby, M.J., 1997. Gully processes and modelling. Progress in Physical Geography 21 Ž3., 354–374. Cordova, J.R., Rodriguez-Iturbe, I., Vaca, P., 1983. On the development of drainage networks. International Association of Hydrological Sciences 137, 35–42. Dietrich, W.E., Dunne, T., 1993. The channel head. In: Beven, K., Kirkby, M.J. ŽEds.., Channel Network Hydrology. Wiley, Chichester, pp. 175–219. Dunne, T., Aubry, B.F., 1986. Evaluation of Horton’s theory of sheetwash and rill erosion on the basis of field experiments. In: Abrahams, A.D. ŽEd.., Hillslope Processes. Allen and Unwin, London, pp. 31–53. Freeze, R.A., 1972. Role of subsurface flow in generating runoff 1. Base flow contributions to channel flow. Water Resources Research 8, 1272–1283. Howard, A.D., 1971. Simulation stream networks by headward growth and branching. Geografiska Annaler 3, 29–50. Howard, A.D., 1988. Groundwater sapping experiments and modelling. In: Howard, A.D., Kochel, R.C., Holt, H.E. ŽEds.., Sapping Features of the Colorado Plateau, a Comparative Planetary Field Guide. NASA Scientific and Technical Information Division, p. 108. Howard, A.D., 1994. A detachment-limited model of drainage basin evolution. Water Resources Research 30, 2261–2285. Howard, A.D., 1997. Badland morphology and evolution: interpretation using a simulation model. Earth Surface Processes and Landforms 22, 211–227. Howard, A.D., MacLane, C.G., 1988. Erosion of cohesionless sediment by groundwater seepage. Water Resources Research 24, 1659–1674. Hugus, M.K., Mark, D.M., 1984. Spatial data processing for digital simulation of erosion. Technical Papers of the Fall Convention of the American Society of PhotogrammetryrAmerican Congress on Surveying and Mapping, pp. 683–693. Kirkby, M.J., 1971. Hillslope process-response models based on the continuity equation. Institute of British Geographers Special Publication 3, 15–20. Kirkby, M.J., 1986. A two dimensional simulation model for slope and stream evolution. In: Abrahams, A.D. ŽEd.., Hillslope Processes. Allen and Unwin, London, pp. 53–73. Kirkby, M.J., 1987. Modelling some influences of soil erosion, landslides and valley gradient on drainage density and hollow development. Catena Supplement 10, 1–11. Kirkby, M.J., 1991. Sediment travel distance as an experimental and model variable in particulate movement. Catena Supplement 19, 111–128. Kirkby, M.J., 1992. An erosion-limited hillslope evolution model. Catena Supplement 23, 157–187. Kirkby, M.J., 1993. Longterm interactions between networks and hillslopes. In: Beven, K., Kirkby, M.J. ŽEds.., Channel Network Hydrology. Wiley, Chichester, pp. 255–293. Kirkby, M.J., 1994. Thresholds and instability in stream head hollows: a model of magnitude and frequency for wash processes. In: Kirkby, M.J. ŽEd.., Process Models and Theoretical Geomorphology. Wiley, Chichester, pp. 295–352. Kirkby, M.J., 1995. Modelling the links between vegetation and landforms. Geomorphology 13, 319–335. Kirkby, M.J., Abrahart, R.J., McMahon, M.L., Shao, J., Thornes, J.B., 1998. MEDALUS soil erosion models for global change. Geomorphology. Montgomery, D.R., Dietrich, D.R., 1989. Source areas, drainage density and channel initiation. Water Resources Research 25, 1907–1918. Poesen, J., Oostwoud Wijdenes, D., Vandekerckhove, L., 1997. MEDALUS III Project 4: Ephemeral Channels and Rivers. Fourth Interim Report, 3–9.

146

M.J. Kirkby, L.J. Bull r Catena 40 (2000) 127–146

Scheideger, A.E., 1970. Theoretical Geomorphology. Prentice-Hall, London. Smart, I.S., Moruzzi, V.L., 1971a. Computer simulation of Clinch Mountain drainage networks. Journal of Geology 79, 572–584. Smart, I.S., Moruzzi, V.L., 1971b. Random-walk model of stream development. Journal of Research and Development 15, 197–203. Smith, T.R., Bretherton, F.P., 1972. Stability and the conservation of mass in drainage basin evolution. Water Resources Research 8, 1506–1529. Vandekerckhove, L., Poesen, J., Oostwoud Wijdenes, D., Figurido, T., submitted. Topographic thresholds for ephemeral gully initiation in intensively cultivated areas of the Mediterranean. Catena. Willgoose, G., Bras, I., Rodriguez-Iturbe, I., 1991a. A coupled network growth and hillslope evolution model 2: nondimensionalization and applications. Water Resources Research 27, 1685–1696. Willgoose, G., Bras, I., Rodriguez-Iturbe, I., 1991b. A coupled channel network growth and hillslope evolution model 1: theory. Water Resources Research 27, 1671–1684. Willgoose, G., Bras, I., Rodriguez-Iturbe, I., 1991c. Results from a new model of river basin evolution. Earth Surface Processes and Landforms 16, 237–254. Willgoose, G., Bras, I., Rodriguez-Iturbe, I., 1994. Hydrogeomorphology modelling with a physically based river basin evolution model. In: Kirkby, M.J. ŽEd.., Process Models and Theoretical Geomorphology. Wiley, Chichester, pp. 3–22.