Physically based modelling of

shallow landslide sediment yield

at a catchment scale
A. Burton 7 J. C. Bathurst

greatly increased by human activities (for example,

through land use change) often with detrimental consequences downstream, such as reservoir siltation, channel
aggradation and instability, and loss of aquatic habitat.
Consequently, as the upland areas prone to landsliding
become increasingly developed, the need increases for
methodologies which can predict, at the basin scale, the
effects of development on landslide incidence and the resulting sediment yield. Basin management strategies can
then be identified which minimize unwanted development impacts, in advance of any development being implemented.
Such aims are increasingly achieved with the support of
mathematical models which incorporate the available understanding of erosion and sediment yield processes and
the degree to which they are affected by human activities.
Because of the need to predict the impacts of change in
advance of the change taking place, only the physically
based, spatially distributed approach to modelling can be
used (for example, Abbott and others 1986a). Until very
recently, our ability to model basin erosion and sediment
yield on a physical and a spatially distributed basis has
been limited to the effects of raindrop impact and overland flow (for example, Park and others 1982; Wicks and
Bathurst 1996). Similar progress has not been achieved
for landslide erosion and sediment yield modelling. InKey words Catchment model 7 GIS analysis 7
deed it is only very recently that any kind of model of
Landslide model 7 Model resolution 7 Sediment
the relationship between landslide erosion and basin seyield
diment yield has been attempted (James 1985; Ziemer
and others 1991a, 1991b). This paper therefore presents a
physically based, spatially distributed model for simulating the basin-scale erosion and subsequent sediment
Introduction
yield arising from landsliding. The paper reviews the
processes which need to be incorporated in the model,
Landsliding, as a form of mass movement, is one of the
presents the model detail and structure, and illustrates
principal processes of hillslope erosion and it can therethe capabilities of the model through a hypothetical apfore play an important role in determining river-basin se- plication to a landslide active catchment in Scotland. A
diment yield. The incidence of landsliding, and thus the
partially completed version of the model was described in
magnitude of the erosion and sediment yield, can be
Burton and Bathurst (1994).
The model considers only rain- and snowmelt-triggered
shallow landslides, that is, those that fail in a translational rather than a rotational manner. Individually these
Received: 30 October 1996 7 Accepted: 25 June 1997
may involve areas of a few hundred square metres and
depths of 12 m; however, large numbers of such slides
A. Burton (Y) 7 J. C. Bathurst
may occur during a single major rainfall event such as a
Water Resource Systems Research Laboratory,
cyclone. In this paper, landslide refers to a shallow hillsDepartment of Civil Engineering, University of Newcastle,
Newcastle upon Tyne, NE1 7RU, UK
lope or soil mass failure at a localized site, landslide eroAbstract A shallow landslide erosion and sediment
yield component, applicable at the basin scale, has
been incorporated into the physically based, spatially distributed, hydrological and sediment transport modelling system, SHETRAN. The component
determines when and where landslides occur in a
basin in response to time-varying rainfall and
snowmelt, the volume of material eroded and released for onward transport, and the impact on basin sediment yield. Derived relationships are used
to link the SHETRAN grid resolution (up to 1 km),
at which the basin hydrology and final sediment
yield is modelled, to a subgrid resolution (typically
around 10100 m) at which landslide occurrence
and erosion is modelled. The subgrid discretization,
landslide susceptibility and potential landslide impact are determined in advance using a geographic
information system (GIS), with SHETRAN then
providing information on temporal variation in the
factors controlling landsliding. The ability to simulate landslide sediment yield is demonstrated by a
hypothetical application based on a catchment in
Scotland.

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

89

sion refers to the consequent loss or release of material

from landslide sites, and landslide sediment yield refers
to that part of the material which, through onward transport from the landslide sites, arrives at the outlet from a
specified area (for example a hillslope or a river basin).

Modelling framework
A past obstacle to the development of a landslide erosion
and sediment yield model has been the lack of a suitable
modelling framework representing the geotechnical, hydrological and hydraulic principles involved. However,
physically based, spatially distributed, hydrological modelling systems are now available which incorporate not
only the overland and channel flow routing needed for
sediment transport modelling but the simulation of soil
moisture conditions, essential in determining the potential for rain- and snowmelt-triggered landsliding. This capability has been achieved through an integrated surface
and subsurface approach to basin modelling, of which the
Systme Hydrologique Europen (SHE) (Abbott and others 1986a, 1986b; Bathurst and others 1995) is perhaps
the most advanced example intended for practical application. SHE is a grid-based, finite-difference modelling
system which represents the major processes of the land
phase of the hydrological cycle. At the Water Resource
Systems Research Laboratory, University of Newcastle
upon Tyne, it has been further developed into SHETRAN,
a water flow, sediment transport and contaminant migration modelling system (Ewen 1995), and it is as a component of SHETRAN that the landslide erosion and sediment yield model is being built.
It is emphasized that the overall aim of the new component is simulation of sediment yield at the basin scale; it
is not the intention to simulate detailed changes in channel geometry or hillslope geomorphology, except in as far
as is required for determining sediment yield. It is also
stressed that the new component is not intended to be an
advanced geotechnical research model or to simulate the
fine geotechnical detail of individual hillslope failures. Its
innovative nature lies in its integration of generally accepted geotechnical techniques in a hydrological and sediment yield modelling framework, to give a capability
for determining the sediment yield from multiple hillslope failures occurring over periods of time ranging from
single rainstorms to several years.

Component specification
The new component is required to simulate the effects of
landslide erosion on the basin sediment yield regime. It
must therefore determine:
1. When and where landslides occur, accounting for the
integrated effects of the relevant controls, especially
the changing soil saturated zone thickness as a trigger.
90

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

2. The volume of material eroded and released for onward transport, and the spatial extent of landsliding
across the basin.
3. The impact on sediment yield at the basin outlet, accounting for transport of material from the landslide
site and areas of sediment deposition within the basin.
The relevant physical processes are now briefly reviewed
so that the equations and rules which form the basis of
the new component can be identified.

Shallow landslide occurrence

Factor of safety
Typically a shallow landslide consists of a slip along an
interface dividing a shallow upper soil layer from an underlying stronger and often less permeable lower soil
layer or bedrock. The soil is subject to two major opposing influences: the downslope component of soil weight,
which acts to shear the soil along a potential failure plane
parallel to the hillslope; and the resistance of the soil to
shearing (known as its shear strength). The relationship
between the two influences is expressed as a factor of safety, FS, where
FS p

Resistance of soil to shearing

Downslope component of soil weight

(1)

If the downslope weight exceeds the resistance to shearing (FS~1), the hillslope is expected to fail. Factor of safety analysis therefore forms the basis for simulating
landslide occurrence in the new component.
Controls on landslide occurrence
Various factors influence slope stability and thus landslide occurrence. Usually they act in combination and
slope stability should not therefore be considered in
terms of just one individual control.
Shear strength is controlled by such factors as the cohesive forces between soil particles, the binding action of
tree roots, and the frictional forces resulting from surface
friction and interlocking between grains of soil or blocks
of rock. In soils or rocks containing water, the frictional
forces are strongly affected by the water pressure in the
voids or pores between the grains or blocks. Typically the
higher the pore-water pressure, the lower are the frictional resistance and the shear strength; increased soil water
content also increases the bulk weight of the soil. Hillslope stability is thus strongly dependent on soil water
content and, in particular, the thickness of the saturated
soil zone relative to the soil depth. Important controls on
landslide occurrence are therefore the input of moisture
from rain and snowmelt (for example, Caine 1980), interception and transpiration of moisture by vegetation (for
example, Greenway 1987) and the concentration of
groundwater by topographic and geological features such
as topographic depressions or hollows (for example, Sidle
and others 1985; Reneau and Dietrich 1987a).

The downslope component of soil weight is a function of

hillslope gradient. Field surveys and theory (for example,
Sidle and others 1985) suggest that there is a lower slope
limit of about 257 (about 50%) for shallow landslides in
saturated soils with typical angles of friction, not reinforced by root binding. Vegetation is of particular importance in determining hillslope stability. Root binding provides the main natural mechanism by which the stability
of a given slope can be increased, while interception and
transpiration have a direct influence on soil moisture
content (Greenway 1987). Changes in vegetation associated with land use change are one of the principal
means by which human activities affect hillslope stability.

others 1998). The range of dimensions can span two or

three orders of magnitude but typically there is a clustering of sizes towards the lower end of the range (for example, Reneau and Dietrich 1987b; Wieczorek 1987).
Reneau and Dietrich (1987b), for example, note that the
scars in their study are typically 1020-m long, 710-m
wide, 0.71.1-m deep and less than 200 m 3 in volume.
The location of the failure plane may lie: within the soil,
below the depth of major rooting; at an interface between
soils of different packing densities; at the bedrock interface; or at an interface corresponding to a decrease in hydraulic conductivity from the overlying soil to the lower
layer.

Integration of geotechnical and hydrological

models
For the factor of safety stability analysis to be applied to
determine landslide occurrence across a basin, a framework is needed to provide the relevant parameters and
variables. From above, these include the basin characteristics (topographic, soil, geological and vegetation properties), rainfall and snowmelt input, and basin hydrological
response in terms of the soil moisture conditions, all spatially and temporally distributed. The framework is provided by an appropriate hydrological modelling system
(in this case SHETRAN). The geotechnical stability analysis then forms an overlay to the hydrological model, enabling the occurrence of landslides in time and space to be
determined as a function of soil moisture content, varying in response to rainfall and snowmelt.
The feasibility of modelling landslide occurrence by applying geotechnical stability analysis to a spatially discretized representation of a basin has been demonstrated by
Ward and others (1981, 1982), Okimura and Kawatani
(1987), Mizuyama (1991), Montgomery and Dietrich
(1994) and Wu and Sidle (1995). Topographic variation is
allowed between discretization elements and the factor of
safety analysis is applied to each element to give the potential for failure. It may be noted, though, that the
quoted models are limited to small basins of a few square
kilometres or less, while the SHETRAN component is intended to apply at scales ranging from less than a square
kilometre to around 500 km 2.

Additional slope failure

Subsequent upslope failure may be triggered as the
ground slope of the upslope soil mass is effectively increased by the removal of the support provided by the
initially failed material. Subsequent downslope failure
may occur as the initially failed material is mobilized as a
debris flow and rides over the soil in its path (for example, Reneau and Dietrich 1987b). Additional soil erosion
is then caused either by scour along the flow path or by
the sudden loading (and failure) of the downslope material by the debris flow.

Transport of failed material

The impact of landslide erosion on basin sediment yield
depends on whether the eroded material is transported to
the channel network. If the material deposits immediately
downslope from the landslide scar, it can enter the channel network only if the landslide is adjacent to the channel or if it can be carried there by overland flow. If the
material breaks up and, through mixing with water,
evolves into a debris flow (for example, Johnson 1984;
Takahashi 1991), it can travel a considerable distance and
so increase its potential for discharging directly into the
channel network. A method is then required for calculating the percentage of the landslide material which is delivered to the channel, that is, the percentage delivery.

Effect of vegetation
Landslides triggered on forested slopes (planar or gully)
Landslide dimensions and volume release such energy and mass that a debris flow nearly always develops. This usually both scours all the material
of eroded material
in its path and continues to move downslope until the
gradient falls below that needed to maintain flow. Some
Two contributions to the total amount of material reexamples are shown in Eschner and Patric (1982) and
leased as a result of landsliding can be identified: materi- DeGraff and others (see Fig. 20 1989). For grass cover the
al from the landslide itself at the initial point of failure,
root-binding effect is weak and failure occurs at condiand material derived from subsequently triggered upslope tions only slightly exceeding soil strength. If the failure
and downslope failures.
occurs in a gully or ephemeral channel (which acts to
concentrate a supply of water), there is still a high likeliCharacteristic landslide dimensions
hood that the material will move downslope as a debris
As indicated by their scars, landslides tend to have typflow. If the failure occurs on a planar slope, though, there
ical dimensions in any region (see Table 2 Burton and
is a much lower likelihood of debris flow generation; the

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

91

material is more likely to drain and stiffen, forming a tail

to the landslide scar.
Debris flow runout distance
If a debris flow moves onto a slope too gentle to support
continued transport, its ability to reach the channel network depends on the distance over which it comes to a
halt the runout distance. No generally satisfactory formula for runout distance has yet been devised. Theoretically based relationships (for example, Takahashi 1991)
are generally inappropriate for the relatively simple level
of modelling envisaged here but several empirical approaches can be identified in the literature. Four such
formulae were tested using hypothetical and measured
hillslope data by Bathurst and others (1997). The formulae are simple relationships which enable runout distance
to be estimated from hillslope geometry. The fraction of
the landslide material which is delivered to the river system is then calculated as:

component for routing to the basin outlet. Because of

limitations in either current knowledge or computing
power, simplifications are necessary to ensure that all the
processes are incorporated into the component and that
the component satisfies its specification requirements.
These simplifications provide fertile ground for further
research.

Dual resolution
The simulation time and memory requirements of SHETRAN are proportional to the number of grid squares (or
rectangles) used in representing spatial distribution in a
basin. For river basins of interest for landslide modelling,
which may typically vary in size up to 500 km 2, the grid
resolution is limited to 0.51 km. However, the preceding
review indicates that landslides and the topographic features which relate to their occurrence typically have plan
dimensions of 1020 m. If SHETRAN were to be used
with a grid resolution of this size, simulations would be
limited to basins smaller than 1 km 2. To overcome this
(WPDz)
(2) problem, a dual resolution approach has been adopted in
FDEL p
W
which the basin hydrology is modelled at the SHETRAN
grid (coarse) resolution while landslide erosion is modwhere W is the runout distance and Dz is the distance
between the point at which debris flow deposition begins elled at a subgrid (fine) resolution. This involves splitting
and the nearest reach of channel in the downslope direc- the process of landslide modelling into two parts, repretion. Bathurst and others (1997) also tested a fifth, statis- sented by the subcomponents GISLIP and SHESLIP
(Fig. 1).
tically based, methodology developed by Ward (1994) to
GISLIP is a geographical information system (GIS) analyestimate percentage delivery directly. None of the apsis which is applied separately from SHETRAN and prior
proaches was accurate over the full range of deliveries
and parallel use of two models was recommended to pro- to the time-varying simulation. It identifies regions of the
vide uncertainty bounds on the estimated percentage de- basin that are at risk from landslides and the soil saturalivery. However, the best compromise between simplicity tion conditions critical for triggering a landslide at a given point. The failure criteria are determined using factor
and reliability was a model based on a study by Vandre
of safety analysis. GISLIP also determines the potential
(1985) in which runout distance is given as:
sediment yield consequences of a landslide at each point,
(3) including: the initial volume of failed material, the volWpa Dy
where Dy is the elevation difference between the head of ume of any material subsequently scoured downslope,
the slide and the point at which deposition begins; and a deposition of eroded material on the ground surface, the
is an empirically derived fraction (set at 0.4 in this case). trajectory of any debris flow and, should the trajectory
intercept the channel network, the volume of sediment
In general, deposition from debris flows tends to begin
delivered for onward channel routing. The analysis is
once the slope falls to around 6107, at least in steep
based on terrain data (in the form of digital terrain modchannel networks. The lowest gradient at which debris
els) that are assumed to be time invariant during the siflow occurs is around 357 (for example, Ikeya 1981;
mulation. GISLIP is designed for use with a fine-resoluJohnson 1984; Benda 1985; Vandre 1985).

Model development
Background to the methodology
The preceding review provides the equations, rules and
process descriptions which form the basis of the SHETRAN landslide erosion and sediment yield component.
In the component, the occurrence of shallow landslides is
determined as a function of the time- and space-varying
soil saturation conditions simulated by SHETRAN. The
volume of eroded material is determined and the proportion of this material reaching the channel network is then
Fig. 1
The dual resolution approach to landslide modelling
calculated and fed to the SHETRAN sediment transport
92

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

tion grid containing on the order of one million squares;

it can therefore represent basins of up to 500 km 2 with a
spatial resolution appropriate for landslide simulation.
SHESLIP carries out the time-varying simulation of landslide occurrence in conjunction with the SHETRAN simulation of basin hydrology and at the SHETRAN grid
(coarse) resolution. At each time step the coarse-grid soil
saturation conditions are compared with the critical conditions required for failure of subgrid (that is, GISLIP
grid) elements (already determined by GISLIP). For this
comparison to be made it is necessary to determine the
subgrid saturated zone thickness as a function of the
SHETRAN grid saturated zone thickness and to this end
a disaggregation technique involving a wetness index is
used. If a subgrid element is then simulated as failing,
the relevant debris flow trajectory and sediment delivery
already defined by GISLIP are applied to the SHETRAN
simulation to provide the sediment yield at the basin
scale. In this way GISLIP provides predetermined information on a look-up basis for the time-varying SHESLIP/SHETRAN simulation.

where
Lp

q0
gsat
gm
c(1Pm)
cm
gw d
gw
gw

(5)

and FS is the factor of safety (FS~1 unsafe, FS61 safe);

Cs is the effective soil cohesion; Cr is the root cohesion; f
is the effective angle of internal friction of soil on impermeable layer; d is the soil depth above failure plane; b
is the slope angle; q0 is the vegetative surcharge per unit
plan area; gsat is the weight density of saturated soil; gm
is the weight density of soil at field moisture content; gw
is the weight density of water; m is the relative saturated
depth (thickness of saturated zone divided by soil depth,
above the failure plane). Van Westen and Terlien (1996)
note that this is the only model which calculates slope instability on a pixel basis and that it is therefore very suitable for use in a raster-based GIS.
Of these terms, most are spatially variable but it is assumed that only m is time-varying, that is, the factor of
safety is a function of m, FS{FS(m). Providing that at
least one of Cs, Cr or tanf is positive, then it can be
shown that FS 1 is well defined.
Assuming that the value of every term, except for m, is
known or can be estimated for each fine-grid square, a
GIS analysis (GISLIP)
critical relative saturated depth, mi c, can be determined
for each fine-grid square, i, where mi cpFSi 1(1). The critObtaining a hydrologically sound DEM
ical relative saturated depth, for square i, is the value of
The first requirement of GISLIP is a fine-resolution digi- m that puts the square on the borderline between safe
tal elevation model (DEM) that is hydrologically sound.
and unsafe. Consequently, the failure condition for each
That is, a clearly defined drainage direction exists for ev- grid square can be written in terms of the time-varying
ery grid square in the DEM. Commonly DEMs suffer
relative saturated depth: for mi^mi c, the slope is safe;
from two faults which prevent this requirement from be- and for mi 1 mi c, the slope is unsafe. In addition, if mi c 1 1
ing satisfied without treatment: 1) spurious pits and
then a square is assumed unconditionally safe; and if
dams; and 2) flat regions. Pits are local elevation minima mi c^0 then a square is assumed unconditionally unsafe.
consisting of a number of connected squares, all adjacent To obtain these conditions, GISLIP determines mi c over
squares to which are at a higher elevation. Dams are spu- the whole catchment.
rious obstructions that block an otherwise well-drained
valley. Flat regions are connected areas of two or more
squares in which all the squares have the same elevation
Wetness index
and at least one square has no clearly defined drainage
direction. The first step in GISLIP is therefore to remove In the time-varying simulation, soil saturation conditions
these faults. The dams problem is overcome by artificially are modelled at the coarse-grid resolution of SHETRAN;
raising the elevation of the surface behind the dam, thus in particular, each SHETRAN grid element is characterised by a single value of saturated zone thickness at each
restoring the original downhill form of the region. Flat
regions are treated by defining an artificial drainage net- time step. Within the area represented by each coarsegrid element, however, there is likely, in reality, to be
work for the region.
considerable spatial variability as a function of topography, soil characteristics and vegetation effects. A wetness
Fine-grid failure conditions
The critical soil saturation conditions for landslide occur- index is therefore being used to link the coarse-scale saturated zone thickness with a subgrid distribution, prorence across the basin are determined through an inverviding appropriate variability in saturated zone thickness
sion of the standard factor of safety analysis (Eq. 1). For
at the fine-grid resolution. At present the effect of only
an infinite slope (the assumption of which is generally
topography is considered, using the index of Beven and
accepted as the basis for modelling shallow landslides)
Kirkby (1979) which takes into account slope and accuthe factor of safety is calculated as (for example, Ward
mulated upslope area. This is defined for each fine-grid
and others 1981):
square as:
2 (CscCr)
(LPm) tan f
c
ai
gw d sin (2 b)
tan b
(6)
(4) Iipln tan b
FS p
i
L

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

93

that would be attained if possible, and that it actually

corresponds to a soil moisture deficit zi which is truncated to either 0 or di as appropriate. The value of zb is
(IPIi)
(7) considered to be the effective mean soil moisture deficit
(zbi Pzb) p
f
within the SHETRAN rectangle, thus relating the potential soil moisture deficits between the contained fine
where Ii is the topographic index for square i; ai is the
squares by Eq. 7. The following procedure is then used to
accumulated runoff area per unit contour length for
account for this problem in order to relate the finesquare i; bi is the slope of square i; zbi is the soil moissquare failure criterion to the average relative saturated
ture deficit of square i; zb and I are the the mean values depth of the SHETRAN rectangle.
of zb and I respectively over fine-resolution squares in
For a fine square, i, that is conditionally unsafe,
the specified region; and f is a basin constant that relates 0^m ci^1, the soil moisture deficit, zb ci will lie in the
transmissivity to depth. The prime on the z indicates that valid range. Therefore:
no account is taken of the possibility that unrealistic val(8)
zb cipz cipdi (1Pm ci)
ues of z may result from the use of Eq. 7. For this stage
of the analysis, GISLIP uses the multiflow direction algorithm, as described in Quinn and others (1991), to deter- From Eq. 7,
mine the topographic index for each fine-grid square in
(IPIi)
(9)
zb cipdi (1Pm ci) P
the catchment.
f
and is related to soil moisture deficit in a specified region by:

Coarse-grid failure conditions

A failure criterion for each fine-grid square is required in
terms of the soil saturation conditions over the parent
SHETRAN grid rectangle, for comparison with the timevarying simulated conditions. This is achieved by combining the soil saturation failure criteria for each finegrid square, with the distribution of soil moisture within
the SHETRAN grid rectangle provided by the topographic
index, Eqs. 6 and 7. Equation 7 describes the subgrid distribution of moisture conditions relative to the mean
condition for a region, in this case a coarse-grid rectangle. It therefore provides a method by which, given a saturated zone thickness at one fine-grid square, it is possible for the saturated zone thicknesses at all the other
fine-grid squares in the parent coarse-grid rectangle to be
calculated. Consequently, knowing the critical saturated
zone thickness for a fine-grid square, the corresponding
values of saturated zone thickness can be found for all
the other fine squares. Aggregated together these give the
critical saturated zone thickness at the coarse-grid scale
which corresponds to the critical condition for the particular fine-grid square. This is a fixed time-invariant value. GISLIP carries out this calculation for each fine-grid
square before the time-varying simulation and stores the
results in the look-up table. During the time-varying
SHESLIP/SHETRAN simulation, all that is then necessary
to determine whether a fine-grid square is at its critical
condition at any time step, is to compare the corresponding critical saturated zone thickness for the parent
coarse-grid rectangle with the actual (time-varying) thickness determined by SHETRAN at that time step.
Although Eq. 7 provides a means to distribute moisture
within a SHETRAN rectangle, it takes no account of the
possibility that, given a set of the other three variables,
both zbi and Zb can be calculated to have an unrealistic
value in the sense that they might exceed the soil depth
or be negative. In order to use this approach within the
model, zbi is considered to be a potential value, such that
zbi outside the valid range is considered to be the value
94

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

is the effective mean soil moisture deficit in the SHETRAN (coarse) rectangle to be critical for the fine square
i. However, in this formulation, the deficit is calculated
without considering whether the corresponding moisture
contents implied for the other fine-grid squares are realistic. A further adjustment is therefore required. Given
the effective mean soil moisture deficit at the coarse-grid
scale, the potential soil moisture deficit can be obtained
from Eq. 7 for the other fine squares j for the conditions
at which fine square i has a critical moisture deficit:
zb cijpdi (1Pm ci) P

(IiPIj)
f

(10)

This value is adjusted in the model as follows:

zijb c is negative, corresponds to exfiltration, so zij c is set
to zero.
zijb c is greater than the soil depth, corresponds to negative moisture, so zij c is set to dj.
Otherwise zijb c is considered to be valid and zij c is set to
the same value.
That is:
0,
z p dj,
zb cij,
c
ij

zb cij~0
zb cij`dj
otherwise

(11)

The actual mean soil moisture deficit across the SHETRAN rectangle is then:
z cipmeanj (z cij)

(12)

Finally, to account for the possibility that the soil depth

representation in a SHETRAN rectangle may not correspond exactly to the mean soil depth in the fine-grid
squares which make up the GISLIP representation of the
rectangle, the total volume of saturated soil in the finegrid squares is set equal to the volume of saturated soil

in the parent coarse rectangle, that is for every SHETRAN rectangle:

1. On a planar grass-covered slope, there is no onward

transport. Eroded material forms a tail to the landslide
scar, with a length typical of the region.
(13) 2. On a planar forested slope, and for landslides occurDZpdz
ring in gullies, the landslide evolves into a debris flow.
where Z is the soil moisture deficit in a SHETRAN recThe
rules applied to govern debris flow transport and setangle; D is the soil depth in a SHETRAN rectangle; and
diment
deposition are as follows:
d and z are the average values of d and z over the fine
1. For slopes greater than 107, the debris flow continues
squares contained within the rectangle.
unconditionally; all soil along the track is scoured and
If M is the relative saturated depth in a SHETRAN recadded to the eroded material from the initial landslide
tangle, then
site.
Z
2.
For slopes between 47 and 107 the debris flow comes
(14)
Mp1 P
D
to a halt either over the runout distance or on reaching the 47 slope, whichever condition is first satisfied.
Substituting for critical values in Eqs. 13 and 14 the folDeposition occurs at a rate such that all the debris
lowing expression is obtained:
flow material would be spread uniformly over the full
runout distance (even if the debris flow reaches the 47
dPz ci
c
(15)
Mi p
slope
first).
D
3. For slopes less than 47 the debris flow halts uncondiKnowing the critical soil moisture condition in fine-grid
tionally and deposits all remaining material.
square i, the corresponding critical relative saturated
The runout criterion is based empirically on Eq. 3. A dedepth for the coarse-grid rectangle containing square i,
bris flow halts once
Mi c, can be obtained. This is the value of M that correDistance travelled on slopes
Elevation lost
sponds to square i being at its critical relative saturated
(16)
`0.4
depth mipmi c. Consequently, Mi c can be determined for
between 47 and 107
on slopes`107
each fine-grid square in the basin and the failure criteriwhere the travel distances are measured along the slope.
on for each fine square can be written entirely in terms
The potential trajectory of the debris flow is simulated by
of the coarse-grid soil saturation conditions. That is, for
starting at the originating fine square and progressing
M^Mi c, the slope is safe; and for M 1 Mi c, the slope is
from square to square down the maximum gradient. The
unsafe.
A hazard map is constructed of only those fine squares in slope and runout distance are calculated using the slope
the catchment that have a potentially unsafe critical rela- distance between the centres of squares in the trajectory.
Deposition and erosion occurring between two squares
tive saturated depth, 0^mi c~1 and whose Mi c permits
are distributed equally between them.
the possibility of failure. Also, as indicated previously,
Within GISLIP, a representation of the channel network
shallow landslides do not generally occur on slopes less
than 257. Consequently, an option exists within the model within the basin relates the SHETRAN channel links to
locations within the fine-grid. If the debris flow is interto exclude all fine squares with slopes less than 257.
cepted by one of these locations then it stops immediately and all remaining material is deposited in the correTransport trajectories, erosion and deposition
For this step of the analysis, GISLIP considers a landslide sponding channel, according to Eq. 2. Equation 16 can be
occurring independently at each fine square in the hazard calibrated for a particular region, in terms of the slope
angles and the decimal fraction. It can also easily be remap. Through simulation of debris flow transport,
placed by alternative formulae (for example, from the
eroded material is partitioned between direct delivery to
survey of Bathurst and others 1997) if these are considthe stream system and deposition along the hillslope.
ered more appropriate for a given application.
As previously noted, landslide scars tend to have typical
dimensions in any region and the failure plane can lie at
a variety of locations below the surface. As a first estiCombination for SHESLIP
mate GISLIP bases the width of the initial failure upon
For each fine square in the hazard map, the details relatthe width of typical landslides in the area or in regions of ing to its potential landslide are expressed in terms of
similar geographical properties; the depth of the initial
their relevance to the time-varying simulation of the
failure is set to the soil depth. These dimensions enable
catchment by SHETRAN at the coarse resolution. It is
the volume of material released by the landslide to be de- not necessary for SHETRAN to know details about the
termined.
fine square that is failing, only details of the coarse recAs described above, there is the potential for onward
tangle containing the potential failure, the critical coarsetransport of eroded material by debris flow to be strongly grid saturation necessary for failure and the supply and
dependent on vegetation cover and whether the initial
removal of sediment to each coarse-grid square and
failure occurs on a planar slope or in a gully. In accorchannel element. The failure criterion for each fine
dance with this, GISLIP applies a rule-based transport
square in the hazard map is already defined in terms of
model:
the relative saturated depth of the parent SHETRAN rec-

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

95

tangle. Consequently, the final step in the GISLIP analysis

is to express for each fine square in the hazard map, the
fine-resolution erosion, deposition and supply to channels of sediment, in terms of sediment removal and
supply to affected coarse-resolution rectangles and channel links.

Demonstration application
A hypothetical application to the topography of a 20-km 2
area containing the Kirkton research catchment at Balquhidder in Scotland (managed for scientific purposes by
the UK Institute of Hydrology) is used to illustrate the
steps in the modelling. The catchment itself is approximately 7 km 2 in area and lies in a steep-sided glaciated
valley with elevations ranging from about 240 m to
850 m. A number of landslides have been documented in
the catchment.
The GIS analysis is based on a fine-grid DEM with a resolution of 20 m within a SHETRAN grid composed of
200-m squares. Figure 2 shows the distribution of the topographic index (Eq. 6) across the catchment; darker
shades indicate greater potential for soil saturation. Figure 3 is a map of critical relative saturated depth obtained from the factor of safety analysis using estimated
terrain data. The scale on this map is such that values in
excess of 1000 indicate squares which are unconditionally
safe and values less than 1000 indicate squares that have
the potential for failure. The map of critical relative saturated depth is then combined with the topographic index
distribution to obtain for each fine square the corresponding critical saturated depth for the SHETRAN par-

Fig. 2
Distribution of topographic index for the test area. See text for
explanation

96

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

Fig. 3
Map of critical relative saturated depth for the test area. Scale is
explained in the text

ent square. This enables a map of potential landslide sites

to be determined, for each of which the debris flow trajectories are calculated and regions of erosion and deposition are determined. The potential trajectories are
shown in Fig. 4, in which black areas indicate deposition
while the light grey shading indicates erosion.
These results were processed to determine their consequences for the coarse-resolution simulation of hydrology
and sediment transport for the 7-km 2 Kirkton catchment.
Figure 5 shows the grid and channel network representation of the catchment for the SHETRAN time-varying simulations. To test the combined SHETRAN and landslide
component model, a single landslide was selected from
the list of potential landslides and its effect upon the sediment regime of the catchment was investigated in isolation. To reduce the number of factors in the simulation,
the main SHETRAN soil erosion component was parameterised so that erosion by raindrop impact, overland flow
and bank erosion was eliminated and only a single sediment type was used in the catchment. Once soil is simulated as eroded by a landslide failure, then any sediment
deposited on the hillslope is available for transport by
overland flow; thus in the simulation the only supply of
sediment in the catchment was from the landslide failure
and the model could demonstrate the subsequent transport of the sediment by debris flow, overland flow and
channel flow.
For test conditions an initial soil saturation of 50% was
set across the whole catchment. Rainfall was then applied
at the rate of 2 mm h 1 for a period in excess of 150 h.
Figure 6a, the hydrograph at the catchment outlet, shows
a flow of less than 0.5 m 3 s 1 for about the first 25 h, after which the lower regions of the catchment are saturated
and the outflow rises in response to the overland flow.

Fig. 4
Debris flow trajectories for the test area. Shading pattern is
explained in the text
Fig. 6
Simulated a water discharge and b sediment discharge at the
channel outlet for the test landslide case

this is followed by a gently increasing sediment flux, although of a lower order of magnitude than the initial sediment discharge peak. This flux is derived from the material which is deposited by the debris flow on the hillslope and which is transported to the channel more slowly
by saturation overland flow along the lower part of the
hillslope. A very similar sediment transport response has
been measured in the field by Ries (1994) for a site in the
Himalayas, where landslide material was first deposited
near a stream channel and then transported rapidly into
the channel by heavy overland flow.

Fig. 5
SHETRAN grid representation of the Kirkton catchment
topography and channel network. Scale shows elevation in
metres

Conclusions and future

developments

The application of the SHETRAN landslide component to

the Kirkton topography is intended to demonstrate the
ability of the component to perform according to its
After 31 h the soil in the square containing the selected
specifications. It does not involve the reproduction of an
landslide site is sufficiently saturated for the landslide to observed response to a measured input and the results
occur; the resulting sediment load is deposited down the should not therefore be interpreted as a typical Kirkton
hillslope and in the channel at the bottom of the hill. Fig- catchment response. However, the test shows that the
ure 6b shows the resulting sediment transport time series overall concept of the model is sound and confirms that
in the main channel. As the debris flow deposits sedithe major elements of the landslide component are in
ment in the channel, there is an abrupt increase in sediplace and can be run in an integrated fashion to produce
ment discharge. Once this sediment has been removed,
results which have a physically realistic interpretation. In
however, there is a sharp decline in transport. In turn
particular, the component has the capability to model:

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

97

1. Landslide incidence as a function of precipitation and

catchment properties, on a spatially distributed basis;
2. Debris flow transport and sediment delivery to the
channel system, allowing also for deposition of sediment along the debris flow track and the subsequent
transport of that sediment by overland flow;
3. Transport of sediment along the river system from the
point of injection.
Overall the component constitutes an innovative approach to modelling shallow landslide erosion and sediment yield at the basin scale. Its dual resolution basis,
whereby basin hydrology is modelled at the SHETRAN
grid scale while landslide erosion is modelled at a subgrid scale, enables the impact of landslide erosion to be
studied for basins of up to 500 km 2 or so. Its use of a
GIS analysis to determine the consequences of potential
landslides prior to carrying out the time-varying simulations (storing the information in look up tables) supports a relatively fast computation procedure. Nevertheless, for this first development the simplest realistic approach to modelling has been adopted and improvements
may be envisaged in the future, especially once practical
experience in applying the model is accumulated. The
modular approach of the GIS part of the model provides
exceptional flexibility which can readily incorporate future modifications.
Validation of the model requires data sets linking multiple landslide occurrences in a basin with the rainfall or
snowmelt conditions triggering the occurrences and with
the resulting sediment yield out of the basin. This is a
demanding requirement, likely to be satisfied only in
small research basins. Even then, and much more so for
larger basins, it is unlikely that satisfactory data sets will
be available for all the physical parameters required by
the component and a degree of parameter estimation
may be required (for example, using data in the literature). There may therefore be considerable uncertainty
attached to parameter evaluation. Further, the use of a
single value for a parameter for a fine-grid square may
not correctly represent the spatial variability within that
square. However, quantification of the uncertainty and a
more accurate representation of the within-grid spatial
variability may be possible if the statistical properties of
the parameters are known. A programme of field studies
at the plot scale in a landslide prone area is currently attempting to address this particular issue (Burton and
others 1998). Similarly, representation of the initial landslide dimensions might be improved by using distributions instead of single values. These considerations form
part of the long-term research programme for the new
component.
Acknowledgements The research for this paper was carried out
as part of the MEDALUS II (Mediterranean Desertification and
Land Use) collaborative research project. MEDALUS II was
funded by the European Commission under its Environment
Programme, contract numbers EV5V-CT92-0128/0164/0165 and
0166, and this support is gratefully acknowledged. The authors
thank the reviewers for their helpful comments.

98

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

References
Abbott MB, Bathurst JC, Cunge JA, OConnell PE, Rasmussen J (1986a) An introduction to the European Hydrological System Systme Hydrologique Europen, SHE, 1 : history and philosophy of a physically-based, distributed modelling system. J Hydrol 87 : 4559
Abbott MB, Bathurst JC, Cunge JA, OConnell PE, Rasmussen J (1986b) An introduction to the European Hydrological System Systme Hydrologique Europen, SHE,
2 : structure of a physically-based, distributed modelling system. J Hydrol 87 : 6177
Bathurst JC, Wicks JM, OConnell PE (1995) The SHE/
SHESED basin scale water flow and sediment transport modelling system. In: Singh VP (ed) Computer models of watershed hydrology. Water Resources Publications, Highlands
Ranch, Colo., pp 563594
Bathurst JC, Burton A, Ward TJ (1997) Debris flow run-out
and landslide sediment delivery model tests. Proc Am Soc Civ
Eng, J Hydraul Eng 123 : 410419
Benda LE (1985) Behavior and effect of debris flows on streams
in the Oregon Coast Range. In: Bowles DS (ed) Delineation of
landslide, flash flood, and debris flow hazards in Utah. General Series Rep. UWRL/G-85/03, Utah Water Research Laboratory, Utah State University, Logan, Utah, pp 153162
Beven KJ, Kirkby MJ (1979) A physically based, variable contributing area model of basin hydrology. Hydrol Sci Bull
24 : 4369
Burton A, Bathurst JC (1994) Modelling shallow landslide
erosion and sediment yield at the basin scale. In: Armanini
A, Di Silvio G (eds) Proceedings International Association
Hydraulic Research International Workshop on Floods and
inundations related to large earth movements, University of
Padua, Padua, Italy, pp B7.1B7.13
Burton A, Arkell TJ, Bathurst JC (1998) Field variability of
landslide model parameters. Environ Geol 35 : 100114
Caine N (1980) The rainfall intensity-duration control of shallow landslides and debris flows. Geogr Ann 62A:2337
DeGraff JV, Bryce R, Jibson RW, Mora S, Rogers CT
(1989) Landslides: their extent and significance in the Caribbean. In: Brabb EE, Harrod BL (eds) Landslides: extent and
economic significance. Balkema, Rotterdam, p 68
Eschner AR, Patric JH (1982) Debris avalanches in eastern
upland forests. J For 80 : 343347
Ewen J (1995) Contaminant transport component of the catchment modelling system (SHETRAN). In: Trudgill ST (ed) Solute modelling in catchment systems. Wiley, Chichester, UK,
pp 417441
Greenway DR (1987) Vegetation and slope stability. In: Anderson MG, Richards KS (eds) Slope stability. Wiley, Chichester,
UK, pp 187230
Ikeya H (1981) A method of designation for area in danger of
debris flow. In: Davies TRH, Pearce AJ (eds) Erosion and sediment transport in Pacific Rim steeplands. International Association Hydrological Science Publ No 132, Institute of Hydrology, Wallingford, Oxon, UK, pp 576588
James LD (1985) Flood hazard measurement who has a ruler?
In: Bowles DS (ed) Delineation of landslide, flash flood, and
debris flow hazard in Utah. General Series Rep UWRL/G-85/
03, Utah Water Research Laboratory, Utah State University,
Logan, Utah, pp 313335
Johnson AM (with Rodine JR) (1984) Debris flow. In: Brunsden D, Prior DB (eds) Slope instability. Wiley, Chichester, UK,
pp 257362

Mizuyama T (1991) Sediment yield and river bed change in

mountain rivers. In: Armanini A, Di Silvio G (eds) Fluvial hydraulics of mountain regions. Lecture Notes in Earth Sciences
No. 37. Springer, Berlin Heidelberg New York, pp 147161
Montgomery DR, Dietrich WE (1994) A physically based
model for the topographic control on shallow landsliding.
Water Resour Res 30 : 11531171
Okimura T, Kawatani T (1987) Mapping of the potential surface-failure sites on granite mountain slopes. In: Gardiner V
(ed) International geomorphology 1986, part I. Wiley, Chichester, UK, pp 121138
Park SW, Mitchell JK, Scarborough JN (1982) Soil erosion
simulation on small watersheds: a modified ANSWERS model. Trans Am Soc Agric Engrs 25 : 15811588
Quinn P, Beven K, Chevallier P, Planchon O (1991) The
prediction of hillslope flow paths for distributed hydrological
modelling using digital terrain models. Hydrol Proc 5 : 5980
Reneau SL, Dietrich WE (1987a) The importance of hollows
in debris flow studies; examples from Marin County, California. In: Costa JE, Wieczorek GF (eds) Debris flows/avalanches: process, recognition, and mitigation. Geological Society of America Reviews in Engineering Geology vol 7, Geological Society of America, Boulder, Colo., pp 165180
Reneau SL, Dietrich WE (1987b) Size and location of colluvial landslides in a steep forested landscape. In: Beschta RL,
Blinn T, Grant GE, Ice GG, Swanson FJ (eds) Erosion and sedimentation in the Pacific Rim. International Association Hydrological Science Publ No 165, Institute of Hydrology, Wallingford, Oxon, UK, pp 3948
Ries JB (1994) Soil erosion in the high mountain region, Eastern Central Himalaya: a case study in the Bamti/Bhandar/Surma area, Nepal. PhD thesis, Heft 42, Institut fr Physische
Geographie der Albert-Ludwigs-Universitt, Freiburg im
Breisgau, Germany
Sidle RC, Pearce AJ, OLoughlin CL (1985) Hillslope stability
and land use. Water Resources Monograph Series 11, American Geophysical Union, Washington, DC
Takahashi T (1991) Debris flow. Balkema, Rotterdam
Vandre BC (1985) Rudd Creek debris flow. In: Bowles DS (ed)
Delineation of landslide, flash flood, and debris flow hazards
in Utah. General Series Rep UWRL/G-85/03, Utah Water Research Laboratory, Utah State University, Logan, Utah,
pp 117131

Van Westen CJ, Terlien MTJ (1996) An approach towards

deterministic landslide hazard analysis in GIS: a case study
from Manizales (Colombia). Earth Surf Processes Landforms
21 : 853868
Ward TJ, Li R-M, Simons DB (1981) Use of a mathematical
model for estimating potential landslide sites in steep forested basins. In: Davies TRH, Pearce AJ (eds) Erosion and sediment transport in Pacific Rim steeplands. International Association Hydrological Science Publ No 132, Institute of Hydrology, Wallingford, Oxon, UK, pp 2141
Ward TJ, Li R-M, Simons DB (1982) Mapping landslide hazards in forest watersheds. Proc Am Soc Civ Eng, J Geotech
Eng Div 108(GT2):319324
Ward TJ (1994) Modeling delivery of landslide materials to
streams. New Mexico Water Resources Research Institute Rep
No 288, New Mexico State University, Las Cruces, N.M.
Wicks JM, Bathurst JC (1996) SHESED: a physically-based,
distributed erosion and sediment yield component for the
SHE hydrological modelling system. J Hydrol 175 : 213238
Wieczorek GF (1987) Landslide erosion in central Santa Cruz
Mountains, California, USA. In: Beschta RL, Blinn T, Grant
GE, Ice GG, Swanson FJ (eds) Erosion and sedimentation in
the Pacific Rim. International Association Hydrological
Science Publ No 165, Institute of Hydrology, Wallingford,
Oxon, UK, pp 489498
Wu W, Sidle RC (1995) A distributed slope stability model for
steep forested basins. Water Resour Res 31 : 20972110
Ziemer RR, Lewis J, Rice RM, Lisle TE (1991a) Modeling the
cumulative watershed effects of forest management strategies.
J Environ Qual 20 : 3642
Ziemer RR, Lewis J, Lisle TE, Rice RM (1991b) Long-term sedimentation effects of different patterns of timber harvesting.
In: Peters NE, Walling DE (eds) Sediment and stream water
quality in a changing environment. International Association
Hydrological Science Publ No 203, Institute of Hydrology,
Wallingford, Oxon, UK, pp 143150

Environmental Geology 35 (23) August 1998 7 Q Springer-Verlag

99