- Email: [email protected]
A distribution function approach to modelling basin sediment yield
Journal of Hydrology, 65 (1983) 239--257
239
Elsevier Science Publishers B.V., Amsterdam -- Printed in The Netherlands
A DISTRIBUTION FUNCTION APPROACH TO MODELLING BASIN SEDIMENT YIELD
R.J. MOORE and R.T. CLARKE
Institute of Hydrology, Wallingford, Oxon, OX10 8BB (Great Britain) (Accepted for publication July 13, 1982)
ABSTRACT Moore, R.J. and Clarke, R.T., 1983. A distribution function approach to modelling basin sediment yield. In: I. Rodriguez-Iturbe and V.K. Gupta (Guest-Editors), Scale Problems in Hydrology. J. Hydrol., 65: 239--257. A new approach to basin sediment yield modelling is proposed based on distribution function theory which provides both a plausible description of sediment removal and translation as supply- and transport-limited processes, and a model structure suitable for automatic parameter optimisation by gradient-based procedures. The mathematical description of accumulation of sediment, and its removal by direct runoff, is capable of representing the exhaustion of sediment supplies during a storm event, and the increasing availability of sediment as the inter-storm period lengthens. Translation of water and sediment to the basin outlet is accomplished using a linear convolution integral operation, and an inverse Gaussian probability density function is proposed as a suitable twoparameter form for the instantaneous unit hydrograph. The relation of this distribution to the convection--diffusion equation, and its extension to represent settling out of sediment from suspension, are discussed. An extension of the distribution function model to incorporate a drainage (or baseflow) component is introduced as one means of representing hydrographs with protracted recessions.
INTRODUCTION
Many sediment transport models proposed in the literature attempt to represent the processes of sediment degradation and aggradation in main river channels (Raudkivi, 1967; Graf, 1971; Bettes and White, 1981; Soni, 1981). Few researchers, however, have addressed the more challenging problem of modelling the processes operating at a basin scale which: (1) make sediment available on the land surface; and (2) remove sediment from the land surface; only then is sediment available for transport by channel flow. Wolman (1977) succinctly draws attention to this limited state of knowledge at the basin scale by stating that "Little is k n o w n about sequential processes involved in the systems of erosion and sedimentation, and practice and theory require attention to unsteady or discontinuous erosion and transportation as sediments move from source through channel systems with intermittent periods of storage. While climatic and hydrologic variation markedly affect yield, transport and deposition thresholds of erosion of cohesive materials and sequences of such effects remain unclear."
0022-1694/83/$03.00
© 1983 Elsevier Science Publishers B.V.
240 The present paper takes a holistic basin-scale view of the sediment system, proposing a new approach to sediment modelling which incorporates the processes of sediment accumulation and detachment within a river basin, and sediment transport to the basin outlet. The precise modes of detachment, the source areas of sediment removal, and the exact paths taken by the translated sediment are not considered but rather a lumped basin model is developed using mean removal rates, and mean travel paths and times. The model is intended for use at a time scale of a day or less and attempts to represent the transitory nature of sediment supply, detachment, and transport processes as a function of meteorological variables, and thereby obtain predictions of basin sediment yield at short time intervals based on meteorological measurements. At any instant in time the sediment yield of a basin will be governed by the supply of sediment available for transport, and the capacity of the basin for transporting sediment to the basin outlet. Clearly a storm following soon after a large storm will be restricted in terms of its ability to generate sediment at the basin outlet by the limited sediment volume available for transport. Sediment will be made available by erosion processes in the inter-storm period, and if this period is long, then the transport capacity will be dictated by the storm magnitude which will constitute the limiting factor in determining the basin sediment yield. A model of basin sediment yield should therefore be capable of representing the increasing availability of sediment as the inter-storm period lengthens, and also reflect the power of the storm to transport the available sediment. No explicit distinction will be made between transport of sediments of different grade; rather the parameters defining the functional relationships describing sediment supply and transport will change in value depending on whether the model is used to predict, for example, suspended sediment or bedload yields. Calibration of the model using real data is used to determine these parameter values, although it is hoped that, ultimately, values may be inferred from basin characteristics. Negev (1967) has developed a basin sediment yield model which incorporates some of the basic hydrological principles discussed above, essentially as an extension to the Stanford Watershed Model (Crawford and Linsley, 1966). Two major shortcomings of this model, however, are: (1) its use of many parameters; and (2) the non-differentiable functional relationships employed. The model is not robust in the sense that model parameters display interdependence and no one set of parameter values minimises the objective function used for model calibration. P.R. Johnston and Pilgrim (1976) and Moore and Clarke (1981) provide critiques of these problems in the c o n t e x t of conceptual rainfall--runoff models, and similar points of criticism may be made with regard to conceptual sediment yield models. Indeed the principle of parsimony of parameters becomes even more important in developing sediment yield models as extensions of existing rainfall--runoff models, because this of necessity leads to even more parameters being
241 employed. The Negev model, in c o m m o n with other models, regards sediment as a threshold limiting p h e n o m e n o n and incorporates threshold elements which lead to difficulties in parameter optimisation, and preclude the use of fast and reliable gradient-based procedures. The new approach presented here avoids the use of non-differentiable threshold elements and employs a gradient-based algorithm for parameter optimisation.
OBJECTIVES The sediment model represents an extension of the distribution function approach to rainfall--runoff modelling developed by Moore and Clarke {1981, 1982). The model extension attempts to present a plausible description of basin sediment removal b y taking into account: (1) The varying availability of sediment over time. (2) The transport of available sediment by a runoff limiting process. The model should be capable of providing adequate predictions of basin sediment yield at time intervals of a day or less, given only rainfall and potential evaporation as data inputs; runoff and suspended sediment data should only be required for model calibration and validation purposes. By satisfying these requirements the model should prove to be a useful tool for predicting the basin sediment yield from storms for which neither flow nor sediment records exist, and in this context will be valuable in extending streamflow and sediment records for periods when only meteorological measurements are available. As few model parameters as possible should be employed, and these should be capable of physical interpretation, in the hope that regionalisation of the model to ungauged basins will be possible. Further, the model should be able to use flow records only -- where they exist -- in the absence of sediment records, to improve the precision of the model's prediction of basin sediment yield.
A REVIEW OF THE RAINFALL--RUNOFF MODEL Since the new approach to rainfall--runoff modelling developed by Moore and Clarke {1981) plays an intrinsic part in the basin sediment yield model, a review of its main features will be given first. A key feature of the approach is the substitution of the single store, c o m m o n l y employed in rainfall--runoff models to represent the processes of interception and soil moisture storage, by an infinite number (statistical population) of stores. This results in a model which can be expressed by a simple algebraic expression which is differentiable at all points in the parameter space, thereby allowing fast and reliable gradient procedures to be used for parameter estimation. The continuous nature of the model structure together with its use of few parameters (two in its basic form)leads to an objective function with a unique minimum,
242
weighting function, 1/G s
'~
flsl~oslexp[s/o's~}
-\
fisl
stores contributing to direct runoff
A th
W
/storage
element
B
si
de[~~s~ ' A W ,
Fig. 1. Definition diagram of the storage distribution function model of direct runoff. and whose value is sensitive to the values taken by all the parameters (the parameters show little interdependence}. A single store (termed a storage element} may be visualized as a narrow tube of depth, s, closed at the b o t t o m and open at the top, and which is filled by rainfall, Pi, generating direct runoff, q, when the tube overflows, and which is depleted by evaporation, Ei, until empty. A weighting, dF = f(s)ds, is given to tubes of a particular depth to reflect their frequency of occurrence in the basin, f(s) being a probability density function (p.d.f.) of an appropriately chosen form; thus f(s)ds is the probability of a store having a depth between s and (s + ds). If stores of different depth are arranged in order of depth with the shallowest on the left and the deepest on the right, with their open ends positioned at the top and all at the same level, then a horizontal line through their tops (AB) and a sloping line through their bottoms (AA ') will form a wedge-shaped diagram as in Fig.l; the weighting attached to stores in a particular depth range is illustrated on the left of this diagram for the case when f(s) is an exponential p.d.f. If each tube is subject to the same net rainfall (rainfall minus potential evaporation, ~i = P~ --El), 7ri ~ 0, then all stores of depth s ~ lri will fill completely, the remainder being filled by an a m o u n t ~i, and having a deficit ( s - ~i} giving rise to a sloping water level line (WW') in the wedge diagram, intersecting the top horizontal line at the point where the vertical depth s is 7ri. If in a dry interval the net rainfall is negative (~i+t < 0), due to evaporation, and if --~i+1 < 7ri, then all stores with depths less than or equal to --~÷~ will be emptied of water, the remainder being depleted by an a m o u n t -- 7ri+t giving rise to a horizontal water level surface in the wedge diagram, cutting the sloping b o t t o m line of the wedge at a point where the store depth is -- 7ri+1 . Fig. 2 shows how an alternating sequence of wet and dry periods gives rise to a water surface profile across the wedge comprised of sloping and horizontal segments. Rainfall which spills from full stores is assumed to form direct runoff, and thus tubes with depths less than that at the point at which the sloping line intersects the horizontal line will generate direct runoff; the store
243 I. T T : - 2
\
i02=2 I* C 1 = ~
3. f'r:2
2. ~=-2
\\ ~ = 4
I°,:'
!C1=0°
4. l~= - 1
5. 1 1 = $
Ic:. '
Ic,=
I
oo
Fig. 2. Updating of contents segments, Ck, and deficit segments, Dk, for an initially saturated basin after five successive periods during which the net rainfalls, hi, are (-- 2, --2,2,--1,5).
depth at this point will be denoted by C*(t). This sloping line intersecting the top horizontal line will move from left to right as a function of the net rainfall, 7ri, and in this way defines the contributing zone of direct runoff. The proportion of stores with depths less than C*(t) will define the proportion of the basin generating direct runoff, and is given by the distribution function: C*(t)
F{C* (t)} = Prob{s•C*(t)}
=
I-
f(z)dz
(1)
0
If A denotes the catchment area, then the contributing area is:
Ac(t ) = AF{C*(t)}
(2)
and if 7ri is the net rainfall at time t, then the direct runoff is:
q(t) = max (O,zri)A¢(t)
(3)
For example, if the distribution of store depths is exponential, then f(s) = as 1 exp (--s/as) and:
q(t) = max (03ri)A[1 -- exp (--C*(t)/os)]
(4)
where the parameter os may be interpreted as the mean store depth. The translation of direct r u n o f f to the basin outlet may be accomplished
244
by considering the bivariate distribution, d F = f(s,t)dsdt, where t is the time taken for direct r u n o f f from storages with depths between (s, s + ds) to reach the basin outlet. Then the translated r u n o f f at the basin outlet at time t is given by: t C*(T)
Q(t) = A [- t- max {O,n(r)}f(s, t--r)dsdr o
(5)
o
For the degenerate case when s and t are independent, d F = f(s)f(t)dsdt, and: t (6)
Q(t) = [ q ( r ) f ( t - - r ) d T 0
where q(T) is given by eq. 3. Moore and Clarke (1981) present a simple closed-form solution to eq. 6 when s and t are assumed to be exponential. The parameter a t of the p.d.f, of translation times, f(t) = a t 1 exp (-- t/ot), is interpreted as the mean translation time.
THE SEDIMENT MODEL
Having reviewed the distribution function approach to rainfall--runoff modelling we can now begin the extension to model basin sediment yield. The time-varying availability of sediment within the basin and its detachm e n t by direct runoff is considered first, the translation of the detached sediment to the basin outlet being dealt with later.
Accumulation and removal of sediment The model extension takes, as a starting point, the following proposition about the availability of sediment: after an element of basin area has ceased to contribute to r u n o f f (say at time to ), sediment that is available for future removal begins to be generated at a m a x i m u m rate R 0, and thereafter at a rate defined by a two-parameter exponentially-decreasing curve (Fig. 3):
R(t) = R o e x p [ - - K ( t - - t 0 ) ]
(7)
where K is a decay rate parameter. At time t the depth of sediment will be given by: t
t
d(t) = | R ( T - - t o ) d T = Ro [ e x p [ - - K ( T - - t 0 ) ] d T to
to
= Ro K-1 [1 -- e x p { - - K ( t - - to)}]
(8)
In the absence o f rainfall, sediment will continue to accumulate at an exponentiaUy-decreasing rate, asymptotically approaching a m a x i m u m depth
245 Ro
Accumulation rate R [t]
Time since last rainfall, t t o
Fig. 3. Sediment accumulation rate curve. R0 / K . . . . : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Sediment depth I o r massl, d It]
/ Time since last rainfall, t~ t o
Fig. 4. Sediment depth (or mass) curve. d(oo) = R o / K (Fig. 4). This depth would then remain on an element of basin
area, until such time as that element contributed to runoff, when the sediment would be completely removed; thereafter sediment would be regenerated at the instantaneous rate R (t). The definition diagram of store depths (Fig. 1) may now be supplemented by another showing the depth of sediment associated with stores of a particular depth (Fig. 5). Here, sediment has built up over progressively longer periods, Tr, depending on when past contributing areas of differing magnitude, Ac(t), last encroached over storage elements of a particular depth to remove sediment. In Fig. 5 the runoff contributing zone expands over stores with depths s ~ C*, and previous rainfall events have caused the contributing zone to extend over stores with depths s ~ Cr. The depth of sediment over any element of basin area will be a function of the time Tr that has elapsed since stores with depths s ~ Cr were last encroached by the contributing zone of runoff, which is determined by past C*-values; clearly T~ 1 ~ T r. In general, the ith step in the sediment profile at time t will be characterised by: (1) C~, the contents of the last contents segment, C*, that extended up to that step; and (2) Tt, the time over which the sediment has accumulated.
246 cIs, T1S = o o ~ c2s, T2s ~ ] c3s,T3 s r-- I c4S,T4 s
t. ~1 st lores cont r~bluting directrunoff
_ _ [ R°/K
! 1 " ~
~ ~ r e c t ,u~o~ T1 s time units previously Fig. 5. Definition diagram of the accumulation-removal component of the sediment distribution function model. These two quantities give the areal extent (using eq. 2) and the depth of each stepped layer (using eq. 8 with t -- to = Tr). Thus during a wet period the contents segment of depth C*(t), such that C~ < C*(t) ~ Ci~_l, will identify the depth of sediment:
d ( t ) = Ro g - I [ 1 - - e x p (--KT~)]
(9)
to be removed. The rate at which sediment is removed in the infinitesimal interval (t, t + d t ) will therefore be: C* ( t )+max( O,TQdt ) 1 r ( t , t + d t ) = ~-/ f d(t+dt)f(s)ds
c*(t) and in the limit as dt --> O, then:
r(t) = max (O,ni)d(t)
f{c*(t)}
(10)
For an exponential distribution of store depths, f(s) = oi I exp (--s/os), and:
r(t) = m a x (Ojri)d(t)a~ 1 exp [ - - C * ( t ) / o s ]
(11)
An important and realistic feature of this representation of basin sediment accumulation and detachment is that a large storm, occurring after a long dry period, will remove a great deal of sediment; a similar storm, following only shortly afterwards, will remove very little.
Translation o f s e d i m e n t to the basin outlet In developing the translation c o m p o n e n t of the model it will be initially assumed that sediment is intimately mixed with the direct runoff, q(t), so that sediment is translated to the basin outlet in the same manner as runoff. Also the translation time t will be taken to be independent of the store depth s, and s and t assumed to have exponential p.d.f.'s. Then the rate at
247
which sediment is removed from the basin, after translation to the basin outlet, will be: t
S(t) = j f(t--T)r(T)dT
(12) o and for an exponential distribution of translation times, this may be expressed in simple recursive form: t+At
S(t + At) = exp (-- A t / o t ) S ( t ) + I
f(t--
T)r(T)dT
(13)
t
The time interval At is chosen such that C*(~) in the expression (10) or expression (11) for r(t) varies according to:
C*(T) = C*(t) + max (0,Tri)(r--t),
t < r < ~ t + At
(14)
and also satisfies the condition imposed by the sediment profile that:
Cr < C*(T) <~C2-1,
t < T <~t + At
(15)
A supplementary expression to eq. 12 or eq. 13 is required when the contents segment C*(t) jumps instantaneously from C*(t-) to C*(t ÷) when a deficit segment is fully replenished by rain, and starts to control r u n o f f and sediment generation (see Moore and Clarke, 1981, fig. 2); the superscripts -- and + are used to denote the contents immediately before and after the jump occur, respectively. In this case sediment from an area identified by the contents range [C* (t-), C* (t÷)] is instantaneously removed, so the depth of sediment removed is given by:
c*(t ÷) d*(t) = [- d(t)f(s)ds C*(t-) = d(t) [exp{-- C* (t-)/a, } -- exp{-- C* (t+)/% }]
(16)
The above integral must in general be evaluated in parts if more than one stepped sediment profile is removed as the contents increases from C*(t-) to C*(t+), such that d(t) is n o t constant. The sediment removal rate at time t ÷ after translation is:
S(t +) = S(t-) + otl d*(t)
(17)
since for f(t) = ot I exp (-- t/at), the initial response to a unit input is ot 1 Finally, it can be shown that when f(s) and f(t) are exponential densities with parameters, %, the mean store depth, and, or, the mean translation time, then eq. 13 (supplemented if necessary by eqs. 16 and 17) may be written as:
S(t) = exp(-- A t / o t ) S ( t - - At) + I(t) (+ o~-ld*(t)) where
(18)
248 I ( t ) = max (O,lri)d(t)a~ 1 a t 1 exp [--{C*(t -- At) + 7r~At}/us] Y
and Y = At{1--exp(--x)}/x;
x
=
COAt;
(~9 =
(Os--TfiOt)/OsO
t
Although basin sediment availability, detachment and transport have been discussed in terms of sediment depth over the basin, in practice the basin sediment yield is required in units of mass. Depth may be equated to the mass equivalent, however, since the expressions derived remain valid. All quantities have been assumed to be in basin-scale units so that the basin area A does not appear, for example, in eq. 18. When both rainfall and runoff are in basin scale units (e.g., mm hr_ ] ), then A in eq. 5 should be replaced by unity.
MODEL EXTENSIONS
We have seen that the distribution function theory of storage depths together with an accumulation profile of sediment mass (or an equivalent. depth) can be used to determine the sediment load and direct r u n o f f to be transported to the basin outlet. For the purposes of illustration, it has been assumed that: (1) sediment and r u n o f f are intimately mixed during transport to the basin outlet; and (2) the translation process can be characterised by an exponential distribution of travel times. The choice of an appropriate travel time distribution is, however, not immediately obvious. Moore and Clarke (1981) selected the exponential distribution of translation times as a simple one-parameter translation model for preliminary analysis. One shortcoming of this model observed when applied to rainfall--runoff data from a number of small upland basins was its inability to represent recessions which decayed rapidly initially, but which later decayed at a much slower rate. Gamma and Weibull distributions of travel times were considered as more suitable candidates due to their positively skewed and unimodal shapes, but little improvement was obtained; since the tails of these two distributions are not too different from exponential they are also unable to represent the "heavy-tailed" nature of observed recessions. Two alternative approaches to overcome the model's shortcomings will be considered here. The first looks to hydrodynamics and the convection-diffusion equation for motivation, and the second to an enhancement of the distribution function model of storage depths to incorporate a drainage (or baseflow) component. The manner whereby each approach allows the intimate mix of sediment and flow assumption to be relaxed is discussed. The convection--diffusion equation
The wide use of the convection--diffusion equation in transport theory suggests that we should look here for motivation in developing a suitable
249
travel time distribution for flow and/or sediment. The convection--diffusion equation may be expressed as the parabolic partial differential equation:
½02
(O2p/Ox2) -- v(Op/Ox) = Op/bt
(19)
where p - - p (x,t) may represent the translated flow (p ~ Q (t)) or sediment load (p ~ S ( t ) ) at time t and at a distance x from its point of origin; 02 is the variance parameter ([L 2 T - l ] ) , and u the drift parameter ([L T - I ] ) . In hydrology these parameters are usually described, using terminology from hydrodynamics, as the convection or wave celerity parameter, u, and D = o 2/2 as the diffusivity. When used to represent translation of flow in an open channel ( p = - Q ( t ) ) , eq. 19 with o z = A o C 2 H ~ / Q o , u = 3 Q o / Z A o , also results from neglecting the inertia terms in the linearised Saint-V+nant equations (Hayami, 1951; Eagleson, 1970); here Q0, H0 and A 0 are the reference flow, depth, and cross-sectional area, and C is the Ch6zy eoefficient (other friction laws could have been assumed). The solution of eq. 19 for a Dirac delta function input at x = 0, t = 0 gives the impulse response function, equivalent to a probability density function:
f(t,x;u,o)
-
x
o(27rt3)1/2 e x p [ - - ( x - - u t ) 2 / 2 o 2 t ] ,
t>O,
u>O
(20)
This solution is well known in first passage time problems occurring in statistics (e.g., Cox and Miller, 1965, p.221), where it is the p.d.f, of the first passage time T for a Wiener process starting at 0 to reach an absorbing barrier at the point x, where u is the positive drift and 02 the variance of the Wiener process. In the present context, T may be interpreted as the random time taken for a water drop (or sediment particle) to reach the basin outlet (the absorbing barrier), x is the distance crossed, and x/~, is the mean translation time of the water drop (or sediment particle). Reparameterisation of the density (eq. 20), such that # = x/u, k = x 2/o 2 , leads to the density:
f(t; la,~) =
~-~) 0,
exp[-- ~(t -- p)2/2U2 t],
t :> 0 otherwise
(21)
where the parameters p and ~ are positive and of dimension [T]. The transformation has essentially led from a Lagrangian description of the transport phenomenon, in which both space and time coordinates are considered, to a Eulerian one in which x is a specified fixed distance. The density in this form is more appropriate for basin-scale problems, when the fixed distance x may be considered as a characteristics length of the basin. Tweedie (1957) termed the density (eq. 21) an inverse Gaussian p.d.f.; Folks and Chhikara (1978) provide a review of its development, and N.L. Johnston and Kotz (1970) summarise its properties. The mean and the variance of this distribution are p, and p3/~, so the parameter p may be interpreted physically as the mean translation time. The distribution is positively skewed and unimodal, with the mode (or time-to-peak) given by:
250
tm ---- 2h/[3 + {9 + 4(X/p)2} '/2]
(22)
An important special case is obtained for zero drift (no convection) when p -+ oo (v = 0) and: f ( t ; h ) = (h/27rt 3)1/2 exp(-- X/2t)
(23)
which is the impulse response function corresponding to the diffusion equation: !2 02 ( 3 2 P / 3 x 2 ) -=
3p/Ot
(24)
.25
f It;p,hl
r
I 1
i 2
3
t
Fig. 6. The inverse Gaussian probability density function for various values of the drift parameter, p, with k = 1. This is a density function having a single parameter, X, and which is unimodal (tin = X/3) and positively skewed. The more general two-parameter inverse Gaussian density (eq. 21) is depicted in Fig. 6 with X = 1, for several values of p. Both densities offer promise as suitable functions for representing the transport of water and/or sediment to the basin outlet. The distribution functions corresponding to these densities will be needed and are given by:
F(t;p,h) = (P [(~/t) 1/2 (-- 1 + t/p)] + exp(2h/p)@[-- (h/t) 1/2 (1 + t/p)] F(t;X) = 2(1)[-- (•/t) 1/2 ]
(25a) (255)
where (P(') is the normal distribution function: Z
(P(z) = (27r)-I/2 j
exp(--u:/2)du
(26)
Alternatively these distribution functions may be expressed in terms of the complement of the error function:
251
Z
erfc(z) = 1 -- err(z) = 1 -- 2(2n) -1/2 ] exp(-- u 2 ) d u oo
0
2(2n) -1/2 f exp(-- u 2 ) d u
(27)
Z
using the i d e n t i t y q)(z) = ½erfc(-- z / v ~ ), to give:
F(t; p, X) = ½[erfc {(X/2t)l/2(1 -- t/p)} + exp(2X/p) erfc {(~./2t)1/2(1 + t/u)}] (28a) F ( t ; X ) = erfc {(X/2t) 1/2}
(28b)
The form of the distribution f u n c t i o n (28a) is illustrated in Fig. 7 for a range of p a r a m e t e r values. D
.25
f
.5
~:[t;p,A) b3 O
~'
'
~t
Fig. 7. The inverse Gaussian distribution function for various values of the drift parameter, p, with k = 1.
Inverse Gaussian distribution function models It was shown earlier t h a t the translation of flow and sediment to the basin o u t l e t led to c o n v o l u t i o n integral expressions: t
Q(t) = ~ q ( r ) f ( t - r ) d r and
(6)
o t
S(t) = ~ r(T)f(t--T)dT
(12)
0
respectively. We will n o w consider h o w these convolutions can be evaluated for the case w h e n the density of translation times, f(t), is given by the inverse Gaussian of eq. 21. Since no simple recursive solution exists, except
252
for the case of exponential densities, the integral is evaluated as the sum of all t h e M t past At (strictly A t i ) increments during which net rainfall ~(T) ~ 0 . Thus for the case of sediment translation: Mt S(t)
It, t~
=
(29)
i=o
where It, t*
t* =
f
r(r)[(t
(30)
r)dr
t~- A t i
and t* is the time at the end of the ith wet interval of duration Ati when lr(r) ~ 0. A simple and computationally efficient approximation to the integral in eq. 30 is the mid-point approximation: (31a)
It, ¢ = r(ti, m ) f ( t -- ti.r, ) A t i
where ti, m = t * - ½Ati. Alternatively if r ( ' ) is assumed to be constant in the interval of duration Ati, and equal to its mid-point value r(ti,m ), then the solution to the integral (30) is: It, ¢ = r(ti, m )[F(t~) -- F ( t * -- At/)]
(31b)
where F ( ' ) is the inverse Gaussian distribution function (28a). S e d i m e n t transport as a c o n v e c t i o n - - d i f f u s i o n process w i t h settling
Although the convection--diffusion equation and its solution in the form of the inverse Gaussian distribution provide an attractive approach to representing both flow and sediment translation processes, in practice the response of the t w o processes can differ considerably. Sediment contrasts with a tracer in solution in not being intimately mixed with the transporting medium, possessing a significant gravity c o m p o n e n t which causes sediment to settle out of suspension at lower flows at a rate dependent on the grade and density of the sediment particles and on the flow r~gime itself. We will consider the affect on the solution of the convection--diffusion equation when an additional term to account for settling of sediment out of suspension is incorporated. The extended equation may be written as: ~pl~t = D (()2pit'X2)
-- l/(()p/c~X) --
s ( p , t)
(32)
where s - s ( p , t ) = OG/~t is the settling rate term; G is the mass of sediment dropped; and p is the translated sediment load (in cases where p is the concentration, c [M L - a ] , then s = V-~o1 ~ G / 3 t where V0 is a characteristic volume). Two simple settling models will be considered. The first, called linear settling, simply assumes that the load settled out, G, is a linear function of the translated load, p, such that:
253 G = klp+k
2
(33)
and eq. 32 becomes: O(~p/Ot) = D(~2pfi)x 2 ) -- u(~p/~x)
(34)
where 0 = 1 + k 1 (or 0 = 1 + k 1/V o i f p is concentration rather than load). It is clear from eq. 34 that the solution obtained for a Dirac delta function input leads to the same inverse Gaussian density (21), except with # and replaced by p* = pO and ~* = ~0. Thus if it is assumed that flow and sedim e n t are both governed by the convection--diffusion equation with celerity u, and diffusivity D, but with sediment subject to additional linear settling, then the above analysis shows how p and ~ are modified to account for this phenomenon and how the underlying inverse Gaussian distribution is retained. The second settling model to be considered assumes that the settling load rate is a linear function of load (or concentration), such that: ~G/~t = k i p
(35)
It can be shown (using results obtained, e.g., by Ogata, 1958; Gupta and Greenkorn, 1973) that the unit-step function for the resulting parabolic partial differential equation: ~p/Ot = D ( O 2 p / O x 2
) - - v(Op/Ox) - - k * p
(36)
where k* = kl i f p is load, and k* = k l / V o i f p is concentration, is: F(t;p,h,~) = ~[exp{~(1
a)}erfc{[~]
+ e x p { ~ (1 + ~)} erfc ( t~-~]
(1 + p t ) } l
(37)
where a = (1 + 2k*p 2/h) 1/2 . This is a defective distribution function (Feller, 1966, p.127) since the settling p h e n o m e n o n introduced into the convection-diffusion equation results in a loss of sediment so that F(oo) < 1. The corresponding p.d.f, can be easily shown to be: f(t;p,h,a)
= [~-~]
exp
~pp2t
t--
(38)
this is the inverse Gaussian density (21) with # replaced by p* = p/a, and unchanged. Consequently assuming that either of the two settling models (33) and (35) hold, then the inverse Gaussian p.d.f, remains appropriate, although with modified parameter vlaues. The constant factor exp[hp -~ (1 -~)] which makes the distribution function defective and which was used as the integrating factor to obtain the p.d.f. (38), can be ignored in the distribution function model since it appears simply as a linear constant. An extension of the second settling model to include the re-suspension of sediment by scour, for example, may be expressed by the equation:
254 OG/~t = k l p - k 2 G
(39)
Lapidus and Amundson (1952) obtain a step-function solution to the resulting partial differential equation, which involves a modified Bessel function of the first kind, is not of closed form, and does not reduce to an inverse Gaussian density function. An approximate closed-form solution was obtained by Ogata (1958), which is equivalent to eq. 37 except for a factor exp(-- k2 t), and therefore also departs from an inverse Gaussian form. Dimensionless analysis and regionalisation Distribution function models based on the inverse Gaussian distribution offer particular promise for regionalisation to ungauged basins using basin characteristics because of: (1) the use of only two parameters; and (2) the physical interpretation given to these parameters through the distribution's relation to the convection--diffusion equation and its extension to include settling o u t of sediment in suspension. To allow for easier comparison of parameter values obtained for basins of differing scale, it is helpful to develop the model in terms of dimensionless variables. The convection--diffusion equation: ~p/~t = D ( ~ 2 p / ~ x 2 ) - v(~p/~x) can be transformed to dimensionless form by introducing the dimensionless variables P = P /Po, X = x /L, and T = vt /L, where the normalising parameters are, P0, the "normal" sediment load (or concentration, or flow), and, L, the characteristic length. This transformation results in the dimensionless convection--diffusion equation to be solved at X = 1 (x = L): 3P/~T
=
l~e 1
(32P/3X 2 ) -- 3P/3X
(40)
where P~ = vL/D is the P~elet number, a dimensionless measure of the relative importance of diffusion and convection. The corresponding dimensionless form of the associated inverse Gaussian density and distribution funetions (21) and (28a) are:
fl (T;P~) = (Pc / 4rrT3 )1/2 e x p [ - - P e (1 -- T) 2/4T]
(41)
and 1[ F I ( T ; P ~ ) = ~ erfe
2(T/pe)l/2
+ exp(Pe) erfc
2(T/Pe)l/21
(42)
also Pe = 2X/p. Drainage from storage elements Although the inverse Gaussian distribution of travel times provides an attractive means of representing the heavy-tailed nature of observed
255
hydrograph recessions, a more conventional solution might a t t e m p t to incorporate a baseflow c o m p o n e n t operating in parallel with direct runoff. The introduction of this p h e n o m e n o n into the distribution function approach to rainfall--runoff modelling will be treated next. It will be assumed that the baseflow c o m p o n e n t does not contribute to basin sediment yield: invoking this assumption accords with the abrupt recession characteristic of the sediment load. In developing the distribution function approach to direct r u n o f f generation it was assumed that the basin was made up of a statistical population of storage elements. Each element was envisaged as a narrow tube of depth s, having a closed bottom and an open top, spillage of water from a full tube giving rise to direct runoff. If the tube is now considered to be open at the b o t t o m allowing drainage to occur at a constant rate 7 until the tube is empty, then the total drainage, b(t) from the population of storage elements at time t can be calculated as follows. The line WW' in Fig. 8 indicates the position of the water level surface across stores of different depth at some time t. Drainage of an a m o u n t 7 occurs in a unit time interval from all stores with contents Ci >~ 7, and by an amount, Ci or (s -- Di), equal to the storage element's water c o n t e n t otherwise. By considering drainage from each of the partitioned segments indicated by vertical dashed lines in Fig. 8, the following general expression for the basin drainage rate at time t can be obtained: kc
b(t) =
Di+l +Ci
kc
~
f
i= n
J Di+l +Ci+I
Ci+Di
(s--D/)f(s)ds + ~
Cif(s)ds
f J
i= n
Ci+Di+ 1
Dn+7
+
(s--D.)f(s)ds+ ; , f(s)ds D n +Cn
A~
I I I
(43)
Dn + T I
i
--
I
Dkd-1
Okd-2
~Z i
I
----B
l
I
l
"\I
eko-,, t ~Z ~ I
I I
I ~ l l l
Al
Fig. 8. D e f i n i t i o n d i a g r a m for t h e storage d i s t r i b u t i o n drainage rate, 7.
function
model
with constant
256
where n is such that Cn < 7 ~< Cn-l. The analytical solutions for the three forms of integral occurring in eq. 43 when the density of store depths is exponential (f(s) = o[ 1 exp(-- s/as) ) are: Ul
f
(s--D)f(s)ds
u0
(u o + o s ) e x p ( - - U o / O s ) - - ( u
1 0+s )
exp(--ul/os)
-- D [exp{-- u0/o~) -- exp(-- u 1/o s )]
Ul
f Cf(s)ds = C[exp(--
u0
/o~)-- exp [-- (t -- T)/O b)]dT
(44)
u0
Tf(s)ds = 7 e x p ( - - u 0 / O s ) u0
Basin drainage consequently occurs at a constant rate 7 when the basin is maintained in a saturated state by rainfall, but decreases as a function of the water c o n t e n t of each storage element and the depth density function assumed for the population of storage elements. Translation of the drainage c o m p o n e n t to the basin outlet may be considered in the same manner as for direct runoff; for an exponential distribution of translation times, f b ( t ) = o~ 1 e x p ( - - t / o b ) , then the translated drainage, or baseflow, is given by: t
QD (t) = exp(-- A t / O b ) Q b { t - - At) + ~ b(~r)og 1 e x p ( - - ( t - - T)/OD)] dT (45) and because o f the principle of superposition for linear systems is given by : t--At
Q { f ) = Qs(t) + Qb(t)
(46)
where Qs(t) is the translated direct runoff. The closed-form solution of the integral in eq. 45 involving b(T) is complicated because C and D are functions of time, i.e. :
C(T) ---- C ( t ) + • i ( r - - t ) ;
D(T) = D ( t ) - - r i ( T - - t
)
where ~i = Pi - - E i - - 7 . In practice the trapezoidal approximation:
[b(t -- At)o~' exp(-- At~oh ) + b(t)] At~2 may be used. For the distribution function model of sediment load the equations developed previously remain unchanged except that the definition of net rainfall, ~i, is modified from Pi -- Ei to Pi -- Ei -- 7. CONCLUSION
Although the extension runoff modelling to model ing, it remains to be seen here perform satisfactorily
of the distribution function theory of rainfall-basin sediment yield may be theoretically appealwhether the various alternative forms presented in practical applications. Consequently, a forth-
257
coming paper will deal with the application of the distribution function model using hourly suspended sediment and flow data from the Creedy Basin in southwest England (Walling and Webb, 1981).
REFERENCES Bettes, R. and White, W.R., 1981. Mathematical simulation of sediment movement in streams. Proc. Inst. Civ. Eng., 71 : 879--892. Cox, D.R. and Miller, H.D., 1965. The Theory of Stochastic Processes. Chapman and Hall, London, pp. 220--222. Crawford, N.H. and Linsley, R.K., 1966. Digital simulation in hydrology: Stanford Watershed Model IV. Stanford Univ., Stanford, Calif., Dep. Civ. Eng., Tech. Rep. No. 39, 210 pp. Eagleson, P.S., 1970. Dynamic Hydrology. McGraw-Hill, New York, N.Y., pp. 363--364. Feller, W., 1966. An Introduction to Probability Theory and its Applications, Vol. II. Wiley, New York, N.Y., p. 127. Folks, J.L. and Chhikara, R.S., 1978. The inverse Gaussian distribution and its statistical application -- a review. J.R. Stat. Soc. Set. B., 40(3): 263--289. Graf, W.H., 1971. Hydraulics of Sediment Transport. McGraw-Hill, New York, N.Y., 513 pp. Gupta, S.P. and Greenkorn, R.A., 1973. Dispersion during flow in porous media with bilinear adsorption. Water Resour. Res., 9(5): 1357--1368. Hayami, S., 1951. On the propagation of flood waves. Disaster Prevent. Res. Inst., K y o t o Univ., Kyoto, Bull. No. 1, 16 pp. Johnston, N.L. and Kotz, S., 1970. Distributions in Statistics: Continuous Univariate Distributions 1. Houghton--Mifflin, Boston, Mass. Johnston, P.R. and Pilgrim, D.H., 1976. Parameter optimisation for watershed models. Water Resour. Res., 12(3): 477--486. Lapidus, L. and Amundson, N.R., 1952. Mathematics of adsorption in beds, VI. The effect of longitudinal diffusion in ion exchange and chromatographic columns. J. Phys. Chem., 56: 984--988. Moore, R.J. and Clarke, R.T., 1981. A distribution function approach to rainfall--runoff modelling. Water Resour. Res., 17(5): 1367--1382. Moore, R.J. and Clarke, R.T., 1982. A distribution function approach to modelling basin soil moisture deficit and streamflow. In: V.P. Singh (Editor), Statistical Analysis of Rainfall and Runoff. Water Resources Publications, F o r t Collins, Colo., pp. 173--190. Negev, M., 1967. A sediment model on a digital computer. Stanford Univ., Stanford, Calif., Dep. Cir. Eng., Tech. Rep. No. 76, 89 pp. Ogata, A., 1958. Dispersion in porous media. Ph.D. Dissertation, Northwestern Univ., Evanston, Ill., 121 pp. Raudkivi, A.J., 1967. Loose Boundary Hydraulics. Pergamon, Oxford, 331 pp. Soni, J.P., 1981. Unsteady sediment transport law and prediction of aggradation parameters. Water Resour. Res., 17(1 ): 33--40. Tweedie, M.C.K., 1957. Statistical properties of inverse Gaussian distributions, I. Ann. Math. Stat., 28: 362--377. Walling, D.E. and Webb, B.W., 1981. The reliability of suspended sediment load data. Proc. Symp. on Erosion and Sediment Transport Measurement, Florence, June 1981. Int. Assoc. Hydrol. Sci., Publ. No. 133,177--194. Wolman, M.G., 1977. Changing needs and opportunities in the sediment field. Water Resour. Res., 13(1): 50--54.
239
Elsevier Science Publishers B.V., Amsterdam -- Printed in The Netherlands
A DISTRIBUTION FUNCTION APPROACH TO MODELLING BASIN SEDIMENT YIELD
R.J. MOORE and R.T. CLARKE
Institute of Hydrology, Wallingford, Oxon, OX10 8BB (Great Britain) (Accepted for publication July 13, 1982)
ABSTRACT Moore, R.J. and Clarke, R.T., 1983. A distribution function approach to modelling basin sediment yield. In: I. Rodriguez-Iturbe and V.K. Gupta (Guest-Editors), Scale Problems in Hydrology. J. Hydrol., 65: 239--257. A new approach to basin sediment yield modelling is proposed based on distribution function theory which provides both a plausible description of sediment removal and translation as supply- and transport-limited processes, and a model structure suitable for automatic parameter optimisation by gradient-based procedures. The mathematical description of accumulation of sediment, and its removal by direct runoff, is capable of representing the exhaustion of sediment supplies during a storm event, and the increasing availability of sediment as the inter-storm period lengthens. Translation of water and sediment to the basin outlet is accomplished using a linear convolution integral operation, and an inverse Gaussian probability density function is proposed as a suitable twoparameter form for the instantaneous unit hydrograph. The relation of this distribution to the convection--diffusion equation, and its extension to represent settling out of sediment from suspension, are discussed. An extension of the distribution function model to incorporate a drainage (or baseflow) component is introduced as one means of representing hydrographs with protracted recessions.
INTRODUCTION
Many sediment transport models proposed in the literature attempt to represent the processes of sediment degradation and aggradation in main river channels (Raudkivi, 1967; Graf, 1971; Bettes and White, 1981; Soni, 1981). Few researchers, however, have addressed the more challenging problem of modelling the processes operating at a basin scale which: (1) make sediment available on the land surface; and (2) remove sediment from the land surface; only then is sediment available for transport by channel flow. Wolman (1977) succinctly draws attention to this limited state of knowledge at the basin scale by stating that "Little is k n o w n about sequential processes involved in the systems of erosion and sedimentation, and practice and theory require attention to unsteady or discontinuous erosion and transportation as sediments move from source through channel systems with intermittent periods of storage. While climatic and hydrologic variation markedly affect yield, transport and deposition thresholds of erosion of cohesive materials and sequences of such effects remain unclear."
0022-1694/83/$03.00
© 1983 Elsevier Science Publishers B.V.
240 The present paper takes a holistic basin-scale view of the sediment system, proposing a new approach to sediment modelling which incorporates the processes of sediment accumulation and detachment within a river basin, and sediment transport to the basin outlet. The precise modes of detachment, the source areas of sediment removal, and the exact paths taken by the translated sediment are not considered but rather a lumped basin model is developed using mean removal rates, and mean travel paths and times. The model is intended for use at a time scale of a day or less and attempts to represent the transitory nature of sediment supply, detachment, and transport processes as a function of meteorological variables, and thereby obtain predictions of basin sediment yield at short time intervals based on meteorological measurements. At any instant in time the sediment yield of a basin will be governed by the supply of sediment available for transport, and the capacity of the basin for transporting sediment to the basin outlet. Clearly a storm following soon after a large storm will be restricted in terms of its ability to generate sediment at the basin outlet by the limited sediment volume available for transport. Sediment will be made available by erosion processes in the inter-storm period, and if this period is long, then the transport capacity will be dictated by the storm magnitude which will constitute the limiting factor in determining the basin sediment yield. A model of basin sediment yield should therefore be capable of representing the increasing availability of sediment as the inter-storm period lengthens, and also reflect the power of the storm to transport the available sediment. No explicit distinction will be made between transport of sediments of different grade; rather the parameters defining the functional relationships describing sediment supply and transport will change in value depending on whether the model is used to predict, for example, suspended sediment or bedload yields. Calibration of the model using real data is used to determine these parameter values, although it is hoped that, ultimately, values may be inferred from basin characteristics. Negev (1967) has developed a basin sediment yield model which incorporates some of the basic hydrological principles discussed above, essentially as an extension to the Stanford Watershed Model (Crawford and Linsley, 1966). Two major shortcomings of this model, however, are: (1) its use of many parameters; and (2) the non-differentiable functional relationships employed. The model is not robust in the sense that model parameters display interdependence and no one set of parameter values minimises the objective function used for model calibration. P.R. Johnston and Pilgrim (1976) and Moore and Clarke (1981) provide critiques of these problems in the c o n t e x t of conceptual rainfall--runoff models, and similar points of criticism may be made with regard to conceptual sediment yield models. Indeed the principle of parsimony of parameters becomes even more important in developing sediment yield models as extensions of existing rainfall--runoff models, because this of necessity leads to even more parameters being
241 employed. The Negev model, in c o m m o n with other models, regards sediment as a threshold limiting p h e n o m e n o n and incorporates threshold elements which lead to difficulties in parameter optimisation, and preclude the use of fast and reliable gradient-based procedures. The new approach presented here avoids the use of non-differentiable threshold elements and employs a gradient-based algorithm for parameter optimisation.
OBJECTIVES The sediment model represents an extension of the distribution function approach to rainfall--runoff modelling developed by Moore and Clarke {1981, 1982). The model extension attempts to present a plausible description of basin sediment removal b y taking into account: (1) The varying availability of sediment over time. (2) The transport of available sediment by a runoff limiting process. The model should be capable of providing adequate predictions of basin sediment yield at time intervals of a day or less, given only rainfall and potential evaporation as data inputs; runoff and suspended sediment data should only be required for model calibration and validation purposes. By satisfying these requirements the model should prove to be a useful tool for predicting the basin sediment yield from storms for which neither flow nor sediment records exist, and in this context will be valuable in extending streamflow and sediment records for periods when only meteorological measurements are available. As few model parameters as possible should be employed, and these should be capable of physical interpretation, in the hope that regionalisation of the model to ungauged basins will be possible. Further, the model should be able to use flow records only -- where they exist -- in the absence of sediment records, to improve the precision of the model's prediction of basin sediment yield.
A REVIEW OF THE RAINFALL--RUNOFF MODEL Since the new approach to rainfall--runoff modelling developed by Moore and Clarke {1981) plays an intrinsic part in the basin sediment yield model, a review of its main features will be given first. A key feature of the approach is the substitution of the single store, c o m m o n l y employed in rainfall--runoff models to represent the processes of interception and soil moisture storage, by an infinite number (statistical population) of stores. This results in a model which can be expressed by a simple algebraic expression which is differentiable at all points in the parameter space, thereby allowing fast and reliable gradient procedures to be used for parameter estimation. The continuous nature of the model structure together with its use of few parameters (two in its basic form)leads to an objective function with a unique minimum,
242
weighting function, 1/G s
'~
flsl~oslexp[s/o's~}
-\
fisl
stores contributing to direct runoff
A th
W
/storage
element
B
si
de[~~s~ ' A W ,
Fig. 1. Definition diagram of the storage distribution function model of direct runoff. and whose value is sensitive to the values taken by all the parameters (the parameters show little interdependence}. A single store (termed a storage element} may be visualized as a narrow tube of depth, s, closed at the b o t t o m and open at the top, and which is filled by rainfall, Pi, generating direct runoff, q, when the tube overflows, and which is depleted by evaporation, Ei, until empty. A weighting, dF = f(s)ds, is given to tubes of a particular depth to reflect their frequency of occurrence in the basin, f(s) being a probability density function (p.d.f.) of an appropriately chosen form; thus f(s)ds is the probability of a store having a depth between s and (s + ds). If stores of different depth are arranged in order of depth with the shallowest on the left and the deepest on the right, with their open ends positioned at the top and all at the same level, then a horizontal line through their tops (AB) and a sloping line through their bottoms (AA ') will form a wedge-shaped diagram as in Fig.l; the weighting attached to stores in a particular depth range is illustrated on the left of this diagram for the case when f(s) is an exponential p.d.f. If each tube is subject to the same net rainfall (rainfall minus potential evaporation, ~i = P~ --El), 7ri ~ 0, then all stores of depth s ~ lri will fill completely, the remainder being filled by an a m o u n t ~i, and having a deficit ( s - ~i} giving rise to a sloping water level line (WW') in the wedge diagram, intersecting the top horizontal line at the point where the vertical depth s is 7ri. If in a dry interval the net rainfall is negative (~i+t < 0), due to evaporation, and if --~i+1 < 7ri, then all stores with depths less than or equal to --~÷~ will be emptied of water, the remainder being depleted by an a m o u n t -- 7ri+t giving rise to a horizontal water level surface in the wedge diagram, cutting the sloping b o t t o m line of the wedge at a point where the store depth is -- 7ri+1 . Fig. 2 shows how an alternating sequence of wet and dry periods gives rise to a water surface profile across the wedge comprised of sloping and horizontal segments. Rainfall which spills from full stores is assumed to form direct runoff, and thus tubes with depths less than that at the point at which the sloping line intersects the horizontal line will generate direct runoff; the store
243 I. T T : - 2
\
i02=2 I* C 1 = ~
3. f'r:2
2. ~=-2
\\ ~ = 4
I°,:'
!C1=0°
4. l~= - 1
5. 1 1 = $
Ic:. '
Ic,=
I
oo
Fig. 2. Updating of contents segments, Ck, and deficit segments, Dk, for an initially saturated basin after five successive periods during which the net rainfalls, hi, are (-- 2, --2,2,--1,5).
depth at this point will be denoted by C*(t). This sloping line intersecting the top horizontal line will move from left to right as a function of the net rainfall, 7ri, and in this way defines the contributing zone of direct runoff. The proportion of stores with depths less than C*(t) will define the proportion of the basin generating direct runoff, and is given by the distribution function: C*(t)
F{C* (t)} = Prob{s•C*(t)}
=
I-
f(z)dz
(1)
0
If A denotes the catchment area, then the contributing area is:
Ac(t ) = AF{C*(t)}
(2)
and if 7ri is the net rainfall at time t, then the direct runoff is:
q(t) = max (O,zri)A¢(t)
(3)
For example, if the distribution of store depths is exponential, then f(s) = as 1 exp (--s/as) and:
q(t) = max (03ri)A[1 -- exp (--C*(t)/os)]
(4)
where the parameter os may be interpreted as the mean store depth. The translation of direct r u n o f f to the basin outlet may be accomplished
244
by considering the bivariate distribution, d F = f(s,t)dsdt, where t is the time taken for direct r u n o f f from storages with depths between (s, s + ds) to reach the basin outlet. Then the translated r u n o f f at the basin outlet at time t is given by: t C*(T)
Q(t) = A [- t- max {O,n(r)}f(s, t--r)dsdr o
(5)
o
For the degenerate case when s and t are independent, d F = f(s)f(t)dsdt, and: t (6)
Q(t) = [ q ( r ) f ( t - - r ) d T 0
where q(T) is given by eq. 3. Moore and Clarke (1981) present a simple closed-form solution to eq. 6 when s and t are assumed to be exponential. The parameter a t of the p.d.f, of translation times, f(t) = a t 1 exp (-- t/ot), is interpreted as the mean translation time.
THE SEDIMENT MODEL
Having reviewed the distribution function approach to rainfall--runoff modelling we can now begin the extension to model basin sediment yield. The time-varying availability of sediment within the basin and its detachm e n t by direct runoff is considered first, the translation of the detached sediment to the basin outlet being dealt with later.
Accumulation and removal of sediment The model extension takes, as a starting point, the following proposition about the availability of sediment: after an element of basin area has ceased to contribute to r u n o f f (say at time to ), sediment that is available for future removal begins to be generated at a m a x i m u m rate R 0, and thereafter at a rate defined by a two-parameter exponentially-decreasing curve (Fig. 3):
R(t) = R o e x p [ - - K ( t - - t 0 ) ]
(7)
where K is a decay rate parameter. At time t the depth of sediment will be given by: t
t
d(t) = | R ( T - - t o ) d T = Ro [ e x p [ - - K ( T - - t 0 ) ] d T to
to
= Ro K-1 [1 -- e x p { - - K ( t - - to)}]
(8)
In the absence o f rainfall, sediment will continue to accumulate at an exponentiaUy-decreasing rate, asymptotically approaching a m a x i m u m depth
245 Ro
Accumulation rate R [t]
Time since last rainfall, t t o
Fig. 3. Sediment accumulation rate curve. R0 / K . . . . : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Sediment depth I o r massl, d It]
/ Time since last rainfall, t~ t o
Fig. 4. Sediment depth (or mass) curve. d(oo) = R o / K (Fig. 4). This depth would then remain on an element of basin
area, until such time as that element contributed to runoff, when the sediment would be completely removed; thereafter sediment would be regenerated at the instantaneous rate R (t). The definition diagram of store depths (Fig. 1) may now be supplemented by another showing the depth of sediment associated with stores of a particular depth (Fig. 5). Here, sediment has built up over progressively longer periods, Tr, depending on when past contributing areas of differing magnitude, Ac(t), last encroached over storage elements of a particular depth to remove sediment. In Fig. 5 the runoff contributing zone expands over stores with depths s ~ C*, and previous rainfall events have caused the contributing zone to extend over stores with depths s ~ Cr. The depth of sediment over any element of basin area will be a function of the time Tr that has elapsed since stores with depths s ~ Cr were last encroached by the contributing zone of runoff, which is determined by past C*-values; clearly T~ 1 ~ T r. In general, the ith step in the sediment profile at time t will be characterised by: (1) C~, the contents of the last contents segment, C*, that extended up to that step; and (2) Tt, the time over which the sediment has accumulated.
246 cIs, T1S = o o ~ c2s, T2s ~ ] c3s,T3 s r-- I c4S,T4 s
t. ~1 st lores cont r~bluting directrunoff
_ _ [ R°/K
! 1 " ~
~ ~ r e c t ,u~o~ T1 s time units previously Fig. 5. Definition diagram of the accumulation-removal component of the sediment distribution function model. These two quantities give the areal extent (using eq. 2) and the depth of each stepped layer (using eq. 8 with t -- to = Tr). Thus during a wet period the contents segment of depth C*(t), such that C~ < C*(t) ~ Ci~_l, will identify the depth of sediment:
d ( t ) = Ro g - I [ 1 - - e x p (--KT~)]
(9)
to be removed. The rate at which sediment is removed in the infinitesimal interval (t, t + d t ) will therefore be: C* ( t )+max( O,TQdt ) 1 r ( t , t + d t ) = ~-/ f d(t+dt)f(s)ds
c*(t) and in the limit as dt --> O, then:
r(t) = max (O,ni)d(t)
f{c*(t)}
(10)
For an exponential distribution of store depths, f(s) = oi I exp (--s/os), and:
r(t) = m a x (Ojri)d(t)a~ 1 exp [ - - C * ( t ) / o s ]
(11)
An important and realistic feature of this representation of basin sediment accumulation and detachment is that a large storm, occurring after a long dry period, will remove a great deal of sediment; a similar storm, following only shortly afterwards, will remove very little.
Translation o f s e d i m e n t to the basin outlet In developing the translation c o m p o n e n t of the model it will be initially assumed that sediment is intimately mixed with the direct runoff, q(t), so that sediment is translated to the basin outlet in the same manner as runoff. Also the translation time t will be taken to be independent of the store depth s, and s and t assumed to have exponential p.d.f.'s. Then the rate at
247
which sediment is removed from the basin, after translation to the basin outlet, will be: t
S(t) = j f(t--T)r(T)dT
(12) o and for an exponential distribution of translation times, this may be expressed in simple recursive form: t+At
S(t + At) = exp (-- A t / o t ) S ( t ) + I
f(t--
T)r(T)dT
(13)
t
The time interval At is chosen such that C*(~) in the expression (10) or expression (11) for r(t) varies according to:
C*(T) = C*(t) + max (0,Tri)(r--t),
t < r < ~ t + At
(14)
and also satisfies the condition imposed by the sediment profile that:
Cr < C*(T) <~C2-1,
t < T <~t + At
(15)
A supplementary expression to eq. 12 or eq. 13 is required when the contents segment C*(t) jumps instantaneously from C*(t-) to C*(t ÷) when a deficit segment is fully replenished by rain, and starts to control r u n o f f and sediment generation (see Moore and Clarke, 1981, fig. 2); the superscripts -- and + are used to denote the contents immediately before and after the jump occur, respectively. In this case sediment from an area identified by the contents range [C* (t-), C* (t÷)] is instantaneously removed, so the depth of sediment removed is given by:
c*(t ÷) d*(t) = [- d(t)f(s)ds C*(t-) = d(t) [exp{-- C* (t-)/a, } -- exp{-- C* (t+)/% }]
(16)
The above integral must in general be evaluated in parts if more than one stepped sediment profile is removed as the contents increases from C*(t-) to C*(t+), such that d(t) is n o t constant. The sediment removal rate at time t ÷ after translation is:
S(t +) = S(t-) + otl d*(t)
(17)
since for f(t) = ot I exp (-- t/at), the initial response to a unit input is ot 1 Finally, it can be shown that when f(s) and f(t) are exponential densities with parameters, %, the mean store depth, and, or, the mean translation time, then eq. 13 (supplemented if necessary by eqs. 16 and 17) may be written as:
S(t) = exp(-- A t / o t ) S ( t - - At) + I(t) (+ o~-ld*(t)) where
(18)
248 I ( t ) = max (O,lri)d(t)a~ 1 a t 1 exp [--{C*(t -- At) + 7r~At}/us] Y
and Y = At{1--exp(--x)}/x;
x
=
COAt;
(~9 =
(Os--TfiOt)/OsO
t
Although basin sediment availability, detachment and transport have been discussed in terms of sediment depth over the basin, in practice the basin sediment yield is required in units of mass. Depth may be equated to the mass equivalent, however, since the expressions derived remain valid. All quantities have been assumed to be in basin-scale units so that the basin area A does not appear, for example, in eq. 18. When both rainfall and runoff are in basin scale units (e.g., mm hr_ ] ), then A in eq. 5 should be replaced by unity.
MODEL EXTENSIONS
We have seen that the distribution function theory of storage depths together with an accumulation profile of sediment mass (or an equivalent. depth) can be used to determine the sediment load and direct r u n o f f to be transported to the basin outlet. For the purposes of illustration, it has been assumed that: (1) sediment and r u n o f f are intimately mixed during transport to the basin outlet; and (2) the translation process can be characterised by an exponential distribution of travel times. The choice of an appropriate travel time distribution is, however, not immediately obvious. Moore and Clarke (1981) selected the exponential distribution of translation times as a simple one-parameter translation model for preliminary analysis. One shortcoming of this model observed when applied to rainfall--runoff data from a number of small upland basins was its inability to represent recessions which decayed rapidly initially, but which later decayed at a much slower rate. Gamma and Weibull distributions of travel times were considered as more suitable candidates due to their positively skewed and unimodal shapes, but little improvement was obtained; since the tails of these two distributions are not too different from exponential they are also unable to represent the "heavy-tailed" nature of observed recessions. Two alternative approaches to overcome the model's shortcomings will be considered here. The first looks to hydrodynamics and the convection-diffusion equation for motivation, and the second to an enhancement of the distribution function model of storage depths to incorporate a drainage (or baseflow) component. The manner whereby each approach allows the intimate mix of sediment and flow assumption to be relaxed is discussed. The convection--diffusion equation
The wide use of the convection--diffusion equation in transport theory suggests that we should look here for motivation in developing a suitable
249
travel time distribution for flow and/or sediment. The convection--diffusion equation may be expressed as the parabolic partial differential equation:
½02
(O2p/Ox2) -- v(Op/Ox) = Op/bt
(19)
where p - - p (x,t) may represent the translated flow (p ~ Q (t)) or sediment load (p ~ S ( t ) ) at time t and at a distance x from its point of origin; 02 is the variance parameter ([L 2 T - l ] ) , and u the drift parameter ([L T - I ] ) . In hydrology these parameters are usually described, using terminology from hydrodynamics, as the convection or wave celerity parameter, u, and D = o 2/2 as the diffusivity. When used to represent translation of flow in an open channel ( p = - Q ( t ) ) , eq. 19 with o z = A o C 2 H ~ / Q o , u = 3 Q o / Z A o , also results from neglecting the inertia terms in the linearised Saint-V+nant equations (Hayami, 1951; Eagleson, 1970); here Q0, H0 and A 0 are the reference flow, depth, and cross-sectional area, and C is the Ch6zy eoefficient (other friction laws could have been assumed). The solution of eq. 19 for a Dirac delta function input at x = 0, t = 0 gives the impulse response function, equivalent to a probability density function:
f(t,x;u,o)
-
x
o(27rt3)1/2 e x p [ - - ( x - - u t ) 2 / 2 o 2 t ] ,
t>O,
u>O
(20)
This solution is well known in first passage time problems occurring in statistics (e.g., Cox and Miller, 1965, p.221), where it is the p.d.f, of the first passage time T for a Wiener process starting at 0 to reach an absorbing barrier at the point x, where u is the positive drift and 02 the variance of the Wiener process. In the present context, T may be interpreted as the random time taken for a water drop (or sediment particle) to reach the basin outlet (the absorbing barrier), x is the distance crossed, and x/~, is the mean translation time of the water drop (or sediment particle). Reparameterisation of the density (eq. 20), such that # = x/u, k = x 2/o 2 , leads to the density:
f(t; la,~) =
~-~) 0,
exp[-- ~(t -- p)2/2U2 t],
t :> 0 otherwise
(21)
where the parameters p and ~ are positive and of dimension [T]. The transformation has essentially led from a Lagrangian description of the transport phenomenon, in which both space and time coordinates are considered, to a Eulerian one in which x is a specified fixed distance. The density in this form is more appropriate for basin-scale problems, when the fixed distance x may be considered as a characteristics length of the basin. Tweedie (1957) termed the density (eq. 21) an inverse Gaussian p.d.f.; Folks and Chhikara (1978) provide a review of its development, and N.L. Johnston and Kotz (1970) summarise its properties. The mean and the variance of this distribution are p, and p3/~, so the parameter p may be interpreted physically as the mean translation time. The distribution is positively skewed and unimodal, with the mode (or time-to-peak) given by:
250
tm ---- 2h/[3 + {9 + 4(X/p)2} '/2]
(22)
An important special case is obtained for zero drift (no convection) when p -+ oo (v = 0) and: f ( t ; h ) = (h/27rt 3)1/2 exp(-- X/2t)
(23)
which is the impulse response function corresponding to the diffusion equation: !2 02 ( 3 2 P / 3 x 2 ) -=
3p/Ot
(24)
.25
f It;p,hl
r
I 1
i 2
3
t
Fig. 6. The inverse Gaussian probability density function for various values of the drift parameter, p, with k = 1. This is a density function having a single parameter, X, and which is unimodal (tin = X/3) and positively skewed. The more general two-parameter inverse Gaussian density (eq. 21) is depicted in Fig. 6 with X = 1, for several values of p. Both densities offer promise as suitable functions for representing the transport of water and/or sediment to the basin outlet. The distribution functions corresponding to these densities will be needed and are given by:
F(t;p,h) = (P [(~/t) 1/2 (-- 1 + t/p)] + exp(2h/p)@[-- (h/t) 1/2 (1 + t/p)] F(t;X) = 2(1)[-- (•/t) 1/2 ]
(25a) (255)
where (P(') is the normal distribution function: Z
(P(z) = (27r)-I/2 j
exp(--u:/2)du
(26)
Alternatively these distribution functions may be expressed in terms of the complement of the error function:
251
Z
erfc(z) = 1 -- err(z) = 1 -- 2(2n) -1/2 ] exp(-- u 2 ) d u oo
0
2(2n) -1/2 f exp(-- u 2 ) d u
(27)
Z
using the i d e n t i t y q)(z) = ½erfc(-- z / v ~ ), to give:
F(t; p, X) = ½[erfc {(X/2t)l/2(1 -- t/p)} + exp(2X/p) erfc {(~./2t)1/2(1 + t/u)}] (28a) F ( t ; X ) = erfc {(X/2t) 1/2}
(28b)
The form of the distribution f u n c t i o n (28a) is illustrated in Fig. 7 for a range of p a r a m e t e r values. D
.25
f
.5
~:[t;p,A) b3 O
~'
'
~t
Fig. 7. The inverse Gaussian distribution function for various values of the drift parameter, p, with k = 1.
Inverse Gaussian distribution function models It was shown earlier t h a t the translation of flow and sediment to the basin o u t l e t led to c o n v o l u t i o n integral expressions: t
Q(t) = ~ q ( r ) f ( t - r ) d r and
(6)
o t
S(t) = ~ r(T)f(t--T)dT
(12)
0
respectively. We will n o w consider h o w these convolutions can be evaluated for the case w h e n the density of translation times, f(t), is given by the inverse Gaussian of eq. 21. Since no simple recursive solution exists, except
252
for the case of exponential densities, the integral is evaluated as the sum of all t h e M t past At (strictly A t i ) increments during which net rainfall ~(T) ~ 0 . Thus for the case of sediment translation: Mt S(t)
It, t~
=
(29)
i=o
where It, t*
t* =
f
r(r)[(t
(30)
r)dr
t~- A t i
and t* is the time at the end of the ith wet interval of duration Ati when lr(r) ~ 0. A simple and computationally efficient approximation to the integral in eq. 30 is the mid-point approximation: (31a)
It, ¢ = r(ti, m ) f ( t -- ti.r, ) A t i
where ti, m = t * - ½Ati. Alternatively if r ( ' ) is assumed to be constant in the interval of duration Ati, and equal to its mid-point value r(ti,m ), then the solution to the integral (30) is: It, ¢ = r(ti, m )[F(t~) -- F ( t * -- At/)]
(31b)
where F ( ' ) is the inverse Gaussian distribution function (28a). S e d i m e n t transport as a c o n v e c t i o n - - d i f f u s i o n process w i t h settling
Although the convection--diffusion equation and its solution in the form of the inverse Gaussian distribution provide an attractive approach to representing both flow and sediment translation processes, in practice the response of the t w o processes can differ considerably. Sediment contrasts with a tracer in solution in not being intimately mixed with the transporting medium, possessing a significant gravity c o m p o n e n t which causes sediment to settle out of suspension at lower flows at a rate dependent on the grade and density of the sediment particles and on the flow r~gime itself. We will consider the affect on the solution of the convection--diffusion equation when an additional term to account for settling of sediment out of suspension is incorporated. The extended equation may be written as: ~pl~t = D (()2pit'X2)
-- l/(()p/c~X) --
s ( p , t)
(32)
where s - s ( p , t ) = OG/~t is the settling rate term; G is the mass of sediment dropped; and p is the translated sediment load (in cases where p is the concentration, c [M L - a ] , then s = V-~o1 ~ G / 3 t where V0 is a characteristic volume). Two simple settling models will be considered. The first, called linear settling, simply assumes that the load settled out, G, is a linear function of the translated load, p, such that:
253 G = klp+k
2
(33)
and eq. 32 becomes: O(~p/Ot) = D(~2pfi)x 2 ) -- u(~p/~x)
(34)
where 0 = 1 + k 1 (or 0 = 1 + k 1/V o i f p is concentration rather than load). It is clear from eq. 34 that the solution obtained for a Dirac delta function input leads to the same inverse Gaussian density (21), except with # and replaced by p* = pO and ~* = ~0. Thus if it is assumed that flow and sedim e n t are both governed by the convection--diffusion equation with celerity u, and diffusivity D, but with sediment subject to additional linear settling, then the above analysis shows how p and ~ are modified to account for this phenomenon and how the underlying inverse Gaussian distribution is retained. The second settling model to be considered assumes that the settling load rate is a linear function of load (or concentration), such that: ~G/~t = k i p
(35)
It can be shown (using results obtained, e.g., by Ogata, 1958; Gupta and Greenkorn, 1973) that the unit-step function for the resulting parabolic partial differential equation: ~p/Ot = D ( O 2 p / O x 2
) - - v(Op/Ox) - - k * p
(36)
where k* = kl i f p is load, and k* = k l / V o i f p is concentration, is: F(t;p,h,~) = ~[exp{~(1
a)}erfc{[~]
+ e x p { ~ (1 + ~)} erfc ( t~-~]
(1 + p t ) } l
(37)
where a = (1 + 2k*p 2/h) 1/2 . This is a defective distribution function (Feller, 1966, p.127) since the settling p h e n o m e n o n introduced into the convection-diffusion equation results in a loss of sediment so that F(oo) < 1. The corresponding p.d.f, can be easily shown to be: f(t;p,h,a)
= [~-~]
exp
~pp2t
t--
(38)
this is the inverse Gaussian density (21) with # replaced by p* = p/a, and unchanged. Consequently assuming that either of the two settling models (33) and (35) hold, then the inverse Gaussian p.d.f, remains appropriate, although with modified parameter vlaues. The constant factor exp[hp -~ (1 -~)] which makes the distribution function defective and which was used as the integrating factor to obtain the p.d.f. (38), can be ignored in the distribution function model since it appears simply as a linear constant. An extension of the second settling model to include the re-suspension of sediment by scour, for example, may be expressed by the equation:
254 OG/~t = k l p - k 2 G
(39)
Lapidus and Amundson (1952) obtain a step-function solution to the resulting partial differential equation, which involves a modified Bessel function of the first kind, is not of closed form, and does not reduce to an inverse Gaussian density function. An approximate closed-form solution was obtained by Ogata (1958), which is equivalent to eq. 37 except for a factor exp(-- k2 t), and therefore also departs from an inverse Gaussian form. Dimensionless analysis and regionalisation Distribution function models based on the inverse Gaussian distribution offer particular promise for regionalisation to ungauged basins using basin characteristics because of: (1) the use of only two parameters; and (2) the physical interpretation given to these parameters through the distribution's relation to the convection--diffusion equation and its extension to include settling o u t of sediment in suspension. To allow for easier comparison of parameter values obtained for basins of differing scale, it is helpful to develop the model in terms of dimensionless variables. The convection--diffusion equation: ~p/~t = D ( ~ 2 p / ~ x 2 ) - v(~p/~x) can be transformed to dimensionless form by introducing the dimensionless variables P = P /Po, X = x /L, and T = vt /L, where the normalising parameters are, P0, the "normal" sediment load (or concentration, or flow), and, L, the characteristic length. This transformation results in the dimensionless convection--diffusion equation to be solved at X = 1 (x = L): 3P/~T
=
l~e 1
(32P/3X 2 ) -- 3P/3X
(40)
where P~ = vL/D is the P~elet number, a dimensionless measure of the relative importance of diffusion and convection. The corresponding dimensionless form of the associated inverse Gaussian density and distribution funetions (21) and (28a) are:
fl (T;P~) = (Pc / 4rrT3 )1/2 e x p [ - - P e (1 -- T) 2/4T]
(41)
and 1[ F I ( T ; P ~ ) = ~ erfe
2(T/pe)l/2
+ exp(Pe) erfc
2(T/Pe)l/21
(42)
also Pe = 2X/p. Drainage from storage elements Although the inverse Gaussian distribution of travel times provides an attractive means of representing the heavy-tailed nature of observed
255
hydrograph recessions, a more conventional solution might a t t e m p t to incorporate a baseflow c o m p o n e n t operating in parallel with direct runoff. The introduction of this p h e n o m e n o n into the distribution function approach to rainfall--runoff modelling will be treated next. It will be assumed that the baseflow c o m p o n e n t does not contribute to basin sediment yield: invoking this assumption accords with the abrupt recession characteristic of the sediment load. In developing the distribution function approach to direct r u n o f f generation it was assumed that the basin was made up of a statistical population of storage elements. Each element was envisaged as a narrow tube of depth s, having a closed bottom and an open top, spillage of water from a full tube giving rise to direct runoff. If the tube is now considered to be open at the b o t t o m allowing drainage to occur at a constant rate 7 until the tube is empty, then the total drainage, b(t) from the population of storage elements at time t can be calculated as follows. The line WW' in Fig. 8 indicates the position of the water level surface across stores of different depth at some time t. Drainage of an a m o u n t 7 occurs in a unit time interval from all stores with contents Ci >~ 7, and by an amount, Ci or (s -- Di), equal to the storage element's water c o n t e n t otherwise. By considering drainage from each of the partitioned segments indicated by vertical dashed lines in Fig. 8, the following general expression for the basin drainage rate at time t can be obtained: kc
b(t) =
Di+l +Ci
kc
~
f
i= n
J Di+l +Ci+I
Ci+Di
(s--D/)f(s)ds + ~
Cif(s)ds
f J
i= n
Ci+Di+ 1
Dn+7
+
(s--D.)f(s)ds+ ; , f(s)ds D n +Cn
A~
I I I
(43)
Dn + T I
i
--
I
Dkd-1
Okd-2
~Z i
I
----B
l
I
l
"\I
eko-,, t ~Z ~ I
I I
I ~ l l l
Al
Fig. 8. D e f i n i t i o n d i a g r a m for t h e storage d i s t r i b u t i o n drainage rate, 7.
function
model
with constant
256
where n is such that Cn < 7 ~< Cn-l. The analytical solutions for the three forms of integral occurring in eq. 43 when the density of store depths is exponential (f(s) = o[ 1 exp(-- s/as) ) are: Ul
f
(s--D)f(s)ds
u0
(u o + o s ) e x p ( - - U o / O s ) - - ( u
1 0+s )
exp(--ul/os)
-- D [exp{-- u0/o~) -- exp(-- u 1/o s )]
Ul
f Cf(s)ds = C[exp(--
u0
/o~)-- exp [-- (t -- T)/O b)]dT
(44)
u0
Tf(s)ds = 7 e x p ( - - u 0 / O s ) u0
Basin drainage consequently occurs at a constant rate 7 when the basin is maintained in a saturated state by rainfall, but decreases as a function of the water c o n t e n t of each storage element and the depth density function assumed for the population of storage elements. Translation of the drainage c o m p o n e n t to the basin outlet may be considered in the same manner as for direct runoff; for an exponential distribution of translation times, f b ( t ) = o~ 1 e x p ( - - t / o b ) , then the translated drainage, or baseflow, is given by: t
QD (t) = exp(-- A t / O b ) Q b { t - - At) + ~ b(~r)og 1 e x p ( - - ( t - - T)/OD)] dT (45) and because o f the principle of superposition for linear systems is given by : t--At
Q { f ) = Qs(t) + Qb(t)
(46)
where Qs(t) is the translated direct runoff. The closed-form solution of the integral in eq. 45 involving b(T) is complicated because C and D are functions of time, i.e. :
C(T) ---- C ( t ) + • i ( r - - t ) ;
D(T) = D ( t ) - - r i ( T - - t
)
where ~i = Pi - - E i - - 7 . In practice the trapezoidal approximation:
[b(t -- At)o~' exp(-- At~oh ) + b(t)] At~2 may be used. For the distribution function model of sediment load the equations developed previously remain unchanged except that the definition of net rainfall, ~i, is modified from Pi -- Ei to Pi -- Ei -- 7. CONCLUSION
Although the extension runoff modelling to model ing, it remains to be seen here perform satisfactorily
of the distribution function theory of rainfall-basin sediment yield may be theoretically appealwhether the various alternative forms presented in practical applications. Consequently, a forth-
257
coming paper will deal with the application of the distribution function model using hourly suspended sediment and flow data from the Creedy Basin in southwest England (Walling and Webb, 1981).
REFERENCES Bettes, R. and White, W.R., 1981. Mathematical simulation of sediment movement in streams. Proc. Inst. Civ. Eng., 71 : 879--892. Cox, D.R. and Miller, H.D., 1965. The Theory of Stochastic Processes. Chapman and Hall, London, pp. 220--222. Crawford, N.H. and Linsley, R.K., 1966. Digital simulation in hydrology: Stanford Watershed Model IV. Stanford Univ., Stanford, Calif., Dep. Civ. Eng., Tech. Rep. No. 39, 210 pp. Eagleson, P.S., 1970. Dynamic Hydrology. McGraw-Hill, New York, N.Y., pp. 363--364. Feller, W., 1966. An Introduction to Probability Theory and its Applications, Vol. II. Wiley, New York, N.Y., p. 127. Folks, J.L. and Chhikara, R.S., 1978. The inverse Gaussian distribution and its statistical application -- a review. J.R. Stat. Soc. Set. B., 40(3): 263--289. Graf, W.H., 1971. Hydraulics of Sediment Transport. McGraw-Hill, New York, N.Y., 513 pp. Gupta, S.P. and Greenkorn, R.A., 1973. Dispersion during flow in porous media with bilinear adsorption. Water Resour. Res., 9(5): 1357--1368. Hayami, S., 1951. On the propagation of flood waves. Disaster Prevent. Res. Inst., K y o t o Univ., Kyoto, Bull. No. 1, 16 pp. Johnston, N.L. and Kotz, S., 1970. Distributions in Statistics: Continuous Univariate Distributions 1. Houghton--Mifflin, Boston, Mass. Johnston, P.R. and Pilgrim, D.H., 1976. Parameter optimisation for watershed models. Water Resour. Res., 12(3): 477--486. Lapidus, L. and Amundson, N.R., 1952. Mathematics of adsorption in beds, VI. The effect of longitudinal diffusion in ion exchange and chromatographic columns. J. Phys. Chem., 56: 984--988. Moore, R.J. and Clarke, R.T., 1981. A distribution function approach to rainfall--runoff modelling. Water Resour. Res., 17(5): 1367--1382. Moore, R.J. and Clarke, R.T., 1982. A distribution function approach to modelling basin soil moisture deficit and streamflow. In: V.P. Singh (Editor), Statistical Analysis of Rainfall and Runoff. Water Resources Publications, F o r t Collins, Colo., pp. 173--190. Negev, M., 1967. A sediment model on a digital computer. Stanford Univ., Stanford, Calif., Dep. Cir. Eng., Tech. Rep. No. 76, 89 pp. Ogata, A., 1958. Dispersion in porous media. Ph.D. Dissertation, Northwestern Univ., Evanston, Ill., 121 pp. Raudkivi, A.J., 1967. Loose Boundary Hydraulics. Pergamon, Oxford, 331 pp. Soni, J.P., 1981. Unsteady sediment transport law and prediction of aggradation parameters. Water Resour. Res., 17(1 ): 33--40. Tweedie, M.C.K., 1957. Statistical properties of inverse Gaussian distributions, I. Ann. Math. Stat., 28: 362--377. Walling, D.E. and Webb, B.W., 1981. The reliability of suspended sediment load data. Proc. Symp. on Erosion and Sediment Transport Measurement, Florence, June 1981. Int. Assoc. Hydrol. Sci., Publ. No. 133,177--194. Wolman, M.G., 1977. Changing needs and opportunities in the sediment field. Water Resour. Res., 13(1): 50--54.