Department of CMl Engineering
MORPHOLOGICAL COMPUTATIONS
Lecturenotes f l O a
prof.dr.ir. M. d e Vries
The s p e c i a l l e c t u r e s on sedimenttransport both f o r the C i v i l Engineering Department
and f o r the I n t e r n a t i o n a l Course f o r Hydraulic Engineering concentrate on mor-p h o l o g i c a l commor-putations f o r r i v e r mor-p r o b l e m s . Therefore the o l d l e c t u r e n o t e s (de V r i e s , 1971) are now replaced by new ones c o n t a i -ning a l s o c o n s i d e r a t i o n s based on research c a r r i e d out i n recent years. These l e c t u r e notes serve f o r an advanced course and
mutatis mutandis are based on the assumption
t h a t the reader i s f a m i l i a r w i t h some basic knowledge on sedimenttransport i n a l l u v i a l channels.
1. Basic considerations 1 1.1. I n t r o d u c t i o n 1 1.2. Propagation o f disturbances 1 1.3. Discussion of c e l e r i t i e s 5 m . R e s t r i c t i o n s 10 2. A n a l y t i c a l models 12 2.1. I n t r o d u c t i o n 12 2.2. Simple wave model 12 2.3. Parabolic model 14 2.4. Morphological time-scale 15 2.5. Hyperbolic model 21 3. Numerical models 22 3.1. I n t r o d u c t i o n 22 3.2. Boundary conditions 22 3.3. Method o f c h a r a c t e r i s t i c s 24 3.4. F i n i t e d i f f e r e n c e s 26 3.5. A p p l i c a t i o n s • 30 3.5. Comments 37 4. Suspended load t r a n s p o r t 4.1. I n t r o d u c t i o n 38 4.2. Theory 38 4.3. Morphological computations 14.5 5. On t r a n s p o r t measurements 5.1. I n t r o d u c t i o n 46 5.2. Measurement o f bedload i|6
5.3. Measijrement o f suspended load 50 5.4. Washload and bedmaterial load 51
5.5. Tracer measurements. 54
- references - l i s t of symbols
66 67
1. Basic c o n s i d e r a t i o n s 1.1. Introduc t i o n .
In r i v e r e n g i n e e r i n g morphological computations g r a d u a l l y become
a t o o l t o f o r e c a s t degradation and aggradation due t o n a t u r a l causes or human i n t e r f e r e n c e . T h e . f r u i t f u l a p p l i c a t i o n o f the
computa-t i o n a l computa-techniques i s only guarancomputa-teed i f
( i ) the morphological phenomena are understood
( i i ) some i n s i g h t i s present i n t o the r e l i a b i l i t y o f the computa-t i o n a l r e s u l computa-t s .
The p h y s i c a l understanding o f the phenomena i s i n the f i r s t approach served i n the f i r s t chapter. The r e s u l t s are used t o describe
analytiaal models (Chapter 2) and nimeriaal models (Chapter 3) both
having t h e i r ovm m e r i t s and shortcomings.
The exlrension o f the p r e s e n t l y a v a i l a b l e models t o cases i n which suspended load predominates i s t r e a t e d i n Chapter 4.
F i n a l l y Chapter 5 discusses some basic aspects o f the measurements of sedimenttransport. These measurements are v i t a l w i t h respect t o the a p p l i c a b i l i t y of the morphological computations f o r p r a c t i c a l cases.
I t has t o be remarked t h a t these l e c t u r e n o t e s due t o r e s t r i c t i o n
i n s i z e do n o t discuss e x p l i c i t l y the roughness o f a l l u v i a l channels, although t h i s t o p i c c e r t a i n l y has t o be considered when morphological computations have t o be c a r r i e d out.
1.2. Propagation o f disturbances.
To understand the morphological behaviour o f a r i v e r b e d , a schematized case i s considered. A s t r a i g h t wide a l l u v i a l channel i s considered. A l l parameters are f u n c t i o n s o f t i m e - ( t ) and place (x ). The dependent v a r i a b l e s are supposed t o be:
• f l o w v e l o c i t y u ( x , t ) • sedimenttransport s ( x , t )
• waterdepth a ( x , t ) ® b e d l e v e l z ( x , t )
g r a i n s i z e (D) are supposed t o be constant. The f o l l o w i n g ana-l y s i s can be given (de V r i e s , 1959, 1965, 1969).
I n t h i s case f o u r equations are necessary t o describe the bed-v a r i a t i o n s :
® The equations of motion and c o n t i n u i t y f o r the f l u i d ® The equations o f motion and c o n t i n u i t y f o r the s e d i
-ment . Hence
a t ^ 8x
s = f ( u )
The f i r s t two equations represent the equation o f motion and c o n t i n u i t y f o r water. Equation ( 3 ) represents the equation o f c o n t i -n u i t y , whereas Eq. ( 4 ) describes the equatio-ns o f motio-n f o r the sediment i n i t s most simple form t o make the f o l l o w i n g deduction as t r a n s p a r e n t as p o s s i b l e .
Equations ( 3 ) and ( 4 ) can e a s i l y be combined:
a_z d f ( u ) 9u _ n at du ax ~
The equations ( 1 ) , ( 2 ) and ( 5 ) are t h r e e l i n e a r equations i n the. p a r t i a l d e r i v a t i v e s o f the remaining t h r e e dependent v a r i a b l e s u, z and a.
Besides a l s o the equations f o r the t o t a l d i f f e r e n t i a l s du, da and dz are a v a i l a b l e d t ^ + dx ^ - du aa , , aa at + 87 = ,^ az _^ , az a t aïï
-The equations ( 1 ) , ( 2 ) , ( 5 ) , ( 6 ) , ( 7 ) and ( 8 ) form a system of s i x l i n e a r equations i n the s i x p a r t i a l d e r i v a t i v e s
This system reads i n m a t r i x f o r m , i f d f ( u ) / d u = f ^ :
1 u 0 g 0 g au/9t 1 R 0 a 1 u 0 0 3u/9x 0 0 f u 0 0 1 0 9a/9t = 0 d t dx 0 Ó 0 , 0 9a/9x du 0 0 d t dx 0 0 9z/9t da 0 0 0 0 dt dx 9z/9x dz ( 9 )
This system can e a s i l y be used t o analyse the p o s s i b i l i t i e s o f the propagation o f disturbances. Consider a d i s c o n t i n u i t y i n t h e s i x
d e r i v a t i v e s . This w i l l be i n d i c a t e d by ind^gterminancy o f t h e s o l u -t i o n of -these d e r i v a -t i v e s . Accordan-t -t o Cramer's r u l e o f de-terminan-ts t h i s indeterminancy r e s u l t s from v a n i s h i n g o f both numerator and denominator.
Hence disocntinuities in the partial derivatives can.only be present i f 1 0 0 u a f u d t dx 0 1 0 0 g u 0 0 d t dx 0 0 def 0 0 1 0 0 0 0 0 ! 0 = 0 d t dx (10) For d t ?i 0 and c = dx/dt t h i s y i e l d s - c^ + 2 UG^ + (ga - u^ + g f ^ ) c - ug f ^ = 0
I t i s convenient t o d e f i n e three dimensiónless parameters (j) = c/u = r e l a t i v e c e l e r i t y
F = u//ga = Froude number
\l) = f^/a^ = dimensionless transportparameter
Combination o f Eqs. (11) and (12) gives
())^ - 2<1)^ t ( 1 - F" ^pF ^)(}) + \l)F~'^ = 0
(11)
(12)
Before considering Eq. (13) i n d e t a i l , the f o l l o w i n g remarks have t o be made.
( i ) For the time being i n . u s i n g Eq. (10) only the c o n d i t i o n f o l l o w i n g from the denominator of Cramer's r u l e i s used. The c o n d i t i o n s f o l l o w i n g f o r the numerators get f u r t h e r a t t e n t i o n i n par. 3.3.
( i i ) The parameter i|; becomes c l e a r i f f o r f ( u ) a simple power law i s used
s = f ( u ) = m u " i i . n ^ . - - , , i . • ^^^^
- ^ n — " ^ I t can be shown e a s i l y t h a t i n t h i s case
1^ = n
J
( 1 5 )Hence i|; i s p r o p o r t i o n a l t o the u s u a l l y small r a t i o o f sedimenttransport and w a t e r t r a n s p o r t . The c o e f f i c i e n t of p r o p o r t i o n a l i t y i s equal t o the exponent o f the powerlaw of Eq. ( 1 4 ) .
I t appears t h a t f o r r e a l i s t i c values o f F and ij; as present i n n a t u r e , the cubic equation (Eq. 13) has always thvee veal roots. Figure 1 gives a g r a p h i c a l r e p r e s e n t a t i o n of the three r o o t s r, ^
1,2,3 as a f u n c t i o n o f Froude number and f o r l i n e s o f equal values of the parameter ^.
1.3. Discussion o f c e l e r i t i e s .
Equation (13) and i t s g r a p h i c a l r e p r e s e n t a t i o n i n F i g . 1 r e q u i r e a s u b s t a n t i a l d i s c u s s i o n .
( i ) I n t h e f i r s t place the case w i t h a fixed bed (no t r a n s p o r t ) has t o be considered. This case i s reached by i n s e r t i n g f^=0 i n Eq. ( 1 3 ) , vgl v.,(3.-i^l>' ip*.-Sf|*'(( i - p"'^
According t o Eq. (13) t h i s leads t o
Or i n a dimensional form (16) =1,2 = ^ ± * ^ F i g . 2. C e l e r i t i e s f o r f i x e d bed! • ( 1 7 ) Equation (16) i s g r a p h i c a l l y represented i n F i g . 2. Hence the s o l u t i o n s o f (j)^(or c ) f o r f i x e d bed are s p e c i a l cases f o r the general case o f a mobile bed.
Apparently here (|)„ ( o r c„) can be i d e n t i f i e d w i t h the c e l e r i t y o f a small d i s t u r -bance a t the bed. For a f i x e d bed n a t u r a l l y = 0.
Note t h a t the use o f a loga-r i t h m i c scale i n F i g s . 1 and 2 forces t o p l o t (J)^ .
( i i ) For moderate Froude numbers according t o F i g . 1 only one r o o t (c() ) i s i n f l u e n c e d by the parameter . This i s due t o the s m a l l values o f ijj present i n n a t u r e . For, the Rhine i n the Netherlands
^ = 10~^ t o lO"^. Thus i n t h i s case a l s o
*1„2 = 1 ± (18)
Or w i t h Eq. (18) and « 1
* 3 = — ^ - - r V ^ i ^ (19)
^ - ( 1 - F ^) ^ ^
(f)^ 2 ^ ® "the normal c e l e r i t i e s a t t h e w a t e r l e v e l as f o r f i x e d bed. Apparently i n t h i s case:
^^12 ^ ® normal a
® (j)^ i s a p o s i t i v e c e l e r i t y o f a ( s m a l l .') disturbance a t t h e bed.
( i i i ) For aritioal flow (F = 1) the equation (13) becomes
(i>^ - 2<^^ - + ^ = 0 ' - ' ' ' ( 2 0 )
I t can be seen from F i g . 1 t h a t (^^ i s not a f f e c t e d by t h e m o b i l i t y o f the bed, hence (J)^ = 2. Moreover (^^ ^^"^ ^^l almost equal.
Thus
-((> 1*2*3 = 'f' (21)
and
<i>2^3 = + ( 2 2 )
( i v ) For superoriHoal flow (F > 1) and not t o o close t o c r i t i c a l f l o w F i g . 1 shows t h a t two r o o t s are again not a f f e c t e d by the m o b i l i t y o f the bed. The two c e l e r i t i e s f o r the w a t e r l e v e l are now both p o s i t i v e . However, as again
-({)^(j)2(|)3 = (23)
the t h i r d r o o t , t h a t can be i d e n t i f i e d w i t h the c e l e r i t y o f a small disturbance a t the bed i s negative. This i s i n accordance the behaviour o f mtidun.es propagating opposite t o t h e flow d i r e c t i o n .
( v ) Considering now again F i g . 1 f o r i n c r e a s i n g Froude number, t h e f o l l o w i n g can be s t a t e d . The p o s i t i v e c e l e r i t y o f the watersurface remains f o r • r e a l i s t i c values o f i|; u n l f f a c t e d by the m o b i l i -t y o f -the bed. The o-ther -two c e l e r i -t i e s change r o l e s . The nega-tive surface c e l e r i t y f o r F « 1 becomes the negative b e d c e l e r i t y f o r
F » 1. On t h e other hand the p o s i t i v e b e d c e l e r i t y f o r F « 1 becomes the p o s i t i v e surface c e l e r i t y f o r F » 1.
Near c r i t i c a l f l o w t h e c e l e r i t i e s <j) and (f) together d e t e r
-2. 3
mine propagation a t surface and bed and f o r F = 1 they become equal i n absolute sence.
The above given discussion o f the c e l e r i t i e s has l e d t o the n o t i o n t h a t the system o f equations ( 9 ) can be reduced sub-s t a n t i a l l y f o r moderate Froude numbersub-s (de V r i e sub-s , 1959).
For the case F^ << 1 the m o b i l i t y of the bed does not e f f e c t the watermovement.
X ^X ^X F i g . 3. Possible schematizations ( f o r F < 1 ) .
Figure 3 shows the x - t plane f o r the general case o f the pro-pagation o f disturbances i n an open channel w i t h a mobile bed. Two extreem schematizations are p o s s i b l e .
Tidal aomputations
Usually w i t h o u t s t a t i n g i t e x p l i c i t l y , f o r t i d a l computations i t i s assumed t h a t c^ ^ leading t o the assumption c^ = 0. I n other words f o r t i d a l computations the m o b i l i t y o f t h e bed i s neglected. Figure 1 shows t h a t f o r moderate Froude numbers t h i s i s j u s t i f i e d .
Movphologiaal aomputations
I f time depending changes o f a r i v e r b e d are s t u d i e d a q u i t e
d i f f e r e n t schematization i s p o s s i b l e . Now from Ic I » c_ i t has 1,2 3 t o be concluded |c^^^2l °° ov dt 0. Thus the f l o w can be c o n s i -dered quasi-steady or,, i n other words the discharge i s not dependent on X.
This case can also be reached by n e g l e c t i n g 9u/9t from the equation o f motion and 9a/9t from the equation o f c o n t i n u i t y of t h e f l u i d .
Hence, f o r moderate Froude numbers 4 Oi u 0 g 0 g 9u/9t R 0 a 0 1 u 0 0 9u/9x 0 0 f u 0 0 1 • 0 9a/9t 0 d t dx 0 0 0 0 aa/9x du 0 0 d t dx 0 0 9z/9t da 0 0 0 0 d t dx 9z/9x dz (24)
As p r e v i o u s l y the c h a r a c t e r i s t i c c e l e r i t i e s can again be found by supposing t h e main determinant o f Eq. (24) being equal zero. This y i e l d s indeed
d t = 0 thus <j)^ 2 * °° (25) and
I n the case o f moderate Froude numbers, the system o f equations reduces t o u 9u/9x + g 9a/9x + g 9z/9x = R ua = q ( t ) 9z/9t t 9s/9x = 0 s = s ( u ) (27) ( 2 8 ) (29) (30)
By e l i m i n a t i n g a and s t h i s leads t o two d i f f e r e n t i a l equations f o r t h e remaining depending variables- u and z vis
and
U'^
9u 9z
Remarks
( i ) Equation (31) i s an o r d i n a r y d i f f e r e n t i a l equation, d e s c r i b i n g a backwater curve.
( i i ) Note t h a t n e g l e c t i n g o f 9a/9t from Eq. 2 leads t o
a du/öx + u 6a/öx = aq/9x = 0 or q = q ( t )
Thus q may s t i l l vary i n t i m e ; t h i s i s important f o r p r a c t i c a l a p p l i c a t i o n s .
1.4. R e s t r i c t i o n s .
The a p p l i c a t i o n o f the above given basic considerations t o a n a t u r a l r i v e r can only be c a r r i e d out w i t h great care. Consider a f r e e l y meandering a l l u v i a l r i v e r , The f o l l o w i n g remarks can then be made.
( i ) E:^<odable banks. For a n a t u r a l channel the e r o d a b i l i t y o f t h e banks may play a r o l e . This r e q u i r e s an a d d i t i o n a l " equation
as the width B ( x , t ) becomes an e x t r a dependent v a r i a b l e . As y e t t h i s equation i s not known and t h e r e f o r e the above given theory i s r e s t r i c t e d t o cases i n which the banks are f i x e d or r e s i s t a n t t o s i g n i f i c a n t erosion.
( i i ) Meanders. The c o n s i d e r a t i o n s are given f o r a s t r a i g h t channel whereas i n nature curved channels are present. This has a number of consequences. I t has been shown (de V r i e s , 1961) t h a t there e x i s t phaselags between depth ( a ) , f l o w v e l o c i t y (u) and mean g r a i n s i z e ( 5 ) along a l i n e p a r a l l e l t o t h e banks ( F i g . 4 ) .
I n the f l o w d i r e c t i o n the depth r e a c t s i n the f i r s t place on the curvature o f the banks. The f l o w v e l o c i t y r e a c t s i n the second place whereas the g r a i n s i z e comes i n t h e l a s t place.
J I i • 1 i
-1 2
3 4 5
XThe presence o f these phaselags cannot be seen u s u a l l y from the d i r e c t observations due t o a l a r g e amount o f s c a t t e r . However, s t a t i s t i c a l methods can be used t o show the t r e n d . As a consequence i t can be stated t h a t t h e mean values a, ü and D are n o t present a t the same place. Hence the a p p l i c a -b i l i t y of t r a n s p o r t f o r m u l a e ( u s u a l l y derived from experiments i n s t r a i g h t l a b o r a t o r y flumes) i s r e s t r i c t e d .
Moreover, as the presence o f a curved channel i s not p o s t u l a t e d i n t h e basic c o n s i d e r a t i o n s , phenomena l i k e bedlevel changes at r i v e r c r o s s i n g s due t o a v a r y i n g discharge are not reproduced a u t o m a t i c a l l y .
( i i i ) Gradation of sediment. I n applying a t r a n s p o r t e q u a t i o n l i k e s = f ( u ) i t i s assumed i m p l i c i t l y t h a t the bedmaterial i s uniform or n e a r l y uniform. I f s i g n i f i c a n t gradation i s present then also D ( x , t ) . The equation o f c o n t i n u i t y f o r the sediment has t o be w r i t t e n f o r each s i z e f r a c t i o n seperately. L i t t l e research i s y e t c a r r i e d out i n t h i s r e s p e c t .
( i v ) Type of transport. The basic c o n s i d e r a t i o n s are based on the assumption t h a t the local sedimenttransport i s a f u n c t i o n o f the local h y d r a u l i c c o n d i t i o n s . This seems c o r r e c t f o r a case i n which bedload t r a n s p o r t dominates. The assumption i s no longer v a l i d i f s u b s t a n t i a l suspended load i s present. Chapter 4 c o n s i -ders t h i s case.
( v ) Alluvial roughness. I n f a c t the h y d r a u l i c roughness as e.g. expressed i n t h e Chezy c o e f f i c i e n t also i s a dependent v a r i a b l e . An accurate p r e d i c t i o n o f the roughness f o r steady f l o w i s s t i l l not possible,(ASCE 1971). For non-steady f l o w p r e d i c t i o n methods do n o t even e x i s t . Therefore here the roughness i s t r e a t e d
as a constant or as a f u n c t i o n o f x o n l y . The value C(x) i s taken from the known s i t u a t i o n a t t = 0 .
2. A n a l y t i c a l models. 2.1. I n t r o d u c t i o n .
A f t e r s i m p l i f y i n g t h e basic .equations f o r t h e non-steady bed-movement f o r t h e case i n which the m o b i l i t y of the bed does not
i n f l u e n c e the c e l e r i t i e s ( c . o f the w a t e r l e v e l , present f o r moderate Froudenumbers a number of a n a l y t i c a l models can be de-duced each a p p l i c a b l e f o r a c e r t a i n number o f cases.
The a p p l i c a b i l i t y depends on the f u r t h e r schematizations made. As the basic equations i n u and z are h i g h l y non l i n e a r (see Eqs. (31) and ( 3 2 ) ) i t i s c l e a r t h a t a n a l y t i c a l s o l u t i o n s can only be expected when l i n e a r i z a t i o n i s c a r r i e d o u t .
The f o l l o w i n g paragraphs give some examples o f a n a l y t i c a l models and the a p p l i c a b i l i t y i s discussed.
2.2. Simple-wave model.
From Eq. ( 3 1 ) and (32) a f i r s t order d i f f e r e n t i a l equation i n z can be derived e a s i l y by e l i m i n a t i n g 3u/9x. The r e s u l t i s
As 3t g df/du gq/u"- - u 9z 9x g df/du _ u df/du . a ^ - u 1 -ga R df/du gq/u^- u -1 = u 1 - F
J
l
(33) (34) Equation (33) can be w r i t t e n as 9z , 9z _ G ^ + c ^ = R . -9t 9x g (35) As R = u |u (36) i t y i e l d s 9z 9 t + c 92 9x u_ C^a a (3©)This i s a simple wave equation i n which c determines t h e propagation of a bedondulation z ( x , t ) and a takes care o f the damping.
For a = O a bedwave w i t h o u t damping i s described. This can be shown as f o l l o w s . Suppose z = z ( y ) w i t h y = x - c t i s a s o l u t i o n o f Eq. (37) w i t h a = 0.
9z/9t = (dz/dy) ( 9 y / 9 t ) = - c dz/dy ' '^>'
and 2^^^ 9z/9#= (dz/dy) (9y/9x) = dz/dy ' ,
Hence , ' ^, ' Oty^ az/9t + c 9z/9x = { - c dz/dy} + c {dz/dy} = 0
Thus z ( y ) w i t h y = x - c t i s a s o l u t i o n o f Eq. (37) f o r a = 0.
The simple wave equation can be used t o study the v a r i a t i o n o f bedloadtransport along a r i v e r b e d .
Suppose a bedform i s propagated w i t h c e l e r i t y Cj^^the h y d r a u l i c r e s i s t a n c e i s neglected, then 9z/9t + 9z/9x = 0 But a l s o 9z/9t + 9s/9x = 0 Combination leads t o 9s/9x - 9z/9x = 0 I n t e g r a t i o n y i e l d s
/
ll
^dx'
= 0^ l^^dx^ (38) Hence s ( x ) = z ( x ) + s^ (39)i f the i n t e g r a t i o n s t a r t s at a through (z = 0) where a t r a n s p o r t i n suspension ( s ^ ) may be present t h a t does not p a r t i c i p a t e i n the propaga-t i o n o f propaga-the bedform. I f propaga-t h i s suspended load can be neglecpropaga-ted, propaga-than propaga-the b e d l e v e l v a r i a t i o n r e f l e c t s d i r e c t l y according t o Eq. (39) the v a r i a t i o n i n t r a n s p o r t .
Or, f o r a t r i a n g u l a r dune: '^''^^X '""^\
'crest = 2 ^ and s = 1 cj^ H (i+o) i n which H denotes the height o f the t r i a n g u l a r dune.
I n p r a c t i c e
^ c r e s t ^ 1«8 s and s 0.6 c^ H (41)
I n par. 5.2. these considerations w i l l be elaborated f u r t h e r t o discuss some g u i d e l i n e s f o r the measurement o f bedload t r a n s p o r t .
The a p p l i c a b i l i t y o f the simple-wave model i s r e s t r i c t e d . I t was used already by Exner (1925) t p describe the propagation Of bedforms.
2.3. Parabolic model.
The pax-abolio model d e r i v e d from the Eqs. (31) and (32) has a wider a p p l i c a b i l i t y than the simple-wave model. The d e r i v a t i o n as r e p o r t e d by Vreugdenhil and de V r i e s (1973) w i l l be given here.
Suppose t h e watermovement i s steady and uniform d u r i n g t r a n s i e n t stages of t h e bed. The equation o f motion then reduces t o
u _u_ C*a 9z 3x Or Thus 9 z u" 3^z - o u2 3u 9ïï^ " 3x Together w i t h Eg. ( 3 2 ) 3 z d f ( u ) 3u _ 3 t du 3x " " e l i m i n a t i o n o f 3u/3x can be c a r r i e d o u t .
This leads t o the p a r a b o l i c d i f f e r e n t i a l equation
(42) (43) (44) f^-O ( 4 5 ) w i t h 3 z ^ 9^z _ n 3 t 3P- - 0 ^ _ 1 C^q d f ( u ) / d u 1 u ds/du K - _ _ _ _ _ _ (46) (47) i n which t h e s u b s c r i p t o r e f e r s t o the o r i g i n a l ( u n i f o r m ) s i t u a t i o n .
A f t e r l i n e a r i z a t i o n , p o s s i b l e f o r u !^ u
K ftj - i ^ ds/du 3 i
o r , more p a r t i c u l a r l y f o r s = m u" K % | n f
The p a r a b o l i c model was already derived q u a l i t a t i v e l y by C u l l i n g (1960). However, he obtained no expression f o r t h e " d i f f u s i o n " c o e f f i c i e n t K. The expressions found by Ashida and Michug (1971) are s i m i l a r t o the ones i n Eqs. (46) and ( 4 8 ) . However, they seem t o l i n e a r i z e t o e a r l y d w i n g the d e r i v a t i o n and t h e r e f o r e obtained i n Eq. (48) a c o e f f i c i e n t 2 i n s t e a d o f ~.
I t has been shown (Vreugdenhil and de V r i e s , 1973) t h a t the p a r a b o l i c model i s o n l y v a l i d f o r large values o f x.
As a r u l e o f thumb
X > 3 a / i ( 5 0 ) Note t h a t a / i i s equal t o t h e l e n g t h L o f a r i v e r r e a c h f o r which the
d i f f e r e n c e i n w a t e r l e v e l i s j u s t equal t o t h e depth a. '
Although t h i s r e s t r i c t s the a p p l i c a b i l i t y o f the p a r a b o l i c model f o r p r a c t i c a l problems, t h i s model has a s u b s t a n t i a l m e r i t .
This becomes c l e a r by r e a l i z i n g t h a t the d e r i v a t i o n o f the p a r a b o l i c model from Eqs. ( 4 2 ) and (45) c o n s i s t s i n d i f f e r e n t i a t i o n , w i t h respect t o X only. This i m p l i e s t h a t the time t i s t r e a t e d as a parameter.
Therefore the p a r a b o l i c model reads.Lin general
9 t - ^ ( * ) 3 3 ^ = 0 ( 5 1 )
Hence t h i s model can be used t o study a n a l y t i c a l l y the behaviour o f
z ( x , t ) f o r a vivevregime, a l b e i t t h a t the r e s t r i c t i o n o f Eq. (50) makes a p p l i c a t i o n o f Eq. (51) only p o s s i b l e f o r l a r g e distances and l a r g e times. An example i s given i n par. 2.4.
2.4. Morphological time-scale.
An i n t e r e s t i n g a p p l i c a t i o n o f the p a r a b o l i c model i s the d e f i n i t i o n o f a morphologiaal time-scale f o r aggradation or degradation i n r i v e r s (de V r i e s , 1973, 1975). The complete d e r i v a t i o n w i l l be given here.
(48)
Consider a r i v e r and suppose t h i s r i v e r i s d i s c h a r g i n g i n t o a lake ( F i g . 5 ) .
At t = 0 the lake l e v e l i s supposed t o drop over a d i s -tance Az. This leads t o de-gradation o f t h e r i v e r b e d u n t i l f o r t -> <» the e n t i r e bed i s lowered w i t h Az.
I t i s convenient t o take x i n t h i s case p o s i t i v e i n
upstream d i r e c t i o n and t o take z = 0 along the o r i g i n a l bed-l e v e bed-l .
This has no i n f l u e n c e on t h e d i f f e r e n t i a l equation. Thus
(52) F i g . 5, D e f i n i t i o n sketch.
z ( x , t ) can be derived from
at
w i t h the boundary c o n d i t i o n s : z(x,o) = 0 z ( 0 , t ) = -Az (53) (54)The s o l u t i o n can be found by a p p l y i n g Laplace transforms, Define z 0 Then = ƒ e"P^ z ( x , t ) d t 9_z 9x -pt 3 z _ 9 - p t _ 9 z f „ d t = 7 r - } ^ Z d t = - r — 9x 9 x / 9x (55) (56) and also 9^2/9x2 9 2 i / 9 x 2 (57) Moreover 9z I t 00 = ) e"P* ( 9 z / 9 t ) d t = cp 00 -pt e ^ z + p] z e~P* d t 0 0 (58)
a p p l y i n g the r u l e f o r p a r t i a l i n t e g r a t i o n / u dv = uv - / v du I t f o l l o w s from Eq. (58) t h a t 9z/3t = pz as l i m e z = 0 t 00 and l i m e z = 0 t 0
i n which t h e boundary c o n d i t i o n o f Eq. (53) i s used. Hence the transformed d i f f e r e n t i a l equation reads
pz - K 3 ^ = 0
(59)
(50)
This i s an o r d i n a r y d i f f e r e n t i a l equation w i t h s o l u t i o n s
z(x,p) = A^(p) exp (X^x) + A2(p) exp (XjX)
Combination of Eqs.. (60) and (.61) y i e l d s
(61)
or
p - K 0 (62)
(63)
The p o s i t i v e value o f X (say X^) w i l l lead t o i n f i n i t e l y l a r g e values o f z. This i s p h y s i c a l l y impossible. Therefore A^(p) = 0 has t o be selected t o exelude t h i s p h y s i c a l i m p o s s i b i l i t y .
The i n t e g r a t i o n constant A2(p) can be derived from the t r a n s f o r m a t i o n of the boundary c o n d i t i o n o f Eq. (54)
z ( 0 , t ) ƒ e-P* (-Az)dt = -Az ƒ e-P* d t = - ^ (64)
Thus w i t h Eq. (61)
z(0,p) = - ^ = A2(p)
The s o l u t i o n o f Eq. (60) t h e r e f o r e becomes
I t can be shown (see e.g. Carslaw and Jaeger, 1963) t h a t the Laplace transform o f a f u n c t i o n e r f c = 4-. / exp(-u2)du (66) reads exp (- av^) P thus by assuming a =
x/
Z
k"
the s o l u t i o n reads f i n a l l y z ( x , t ) = - Az e r f c (67) (68) +x 2/Kt (69)The oomplementary errorfunotion as d e f i n e d by Eq. (67) has been com-puted f o r v a r i o u s values o f the argument (see e.g. Crank, 1957). Some data are given i n Table.1.
y -1.0 -0,5 -0.2 -0.1 0 0.1 0.2 _ 0.5 1..0 2.0
e r f c y -0.16 -0.48 -0.78 -0.89 1.00 0.89 0.78 0-.48 0.16 0.005 Table 1. Complementary e r r o r f u n c t i o n (two decimals )b.
As has been shown i n par. 2.3. the p a r a b o l i c model i s a l s o v a l i d f o r a v a r y i n g discharge leading t o K = K ( t ) .
In t h i s case instead of Eq. (69) the s o l u t i o n becomes z(x,T) = - Az e r f c
ƒ K ( t ) d t
For the d e f i n i t i o n o f a morphological time-scale the f o l l o w i n g reasoning can be given. Consider a standard l e n g t h L . The question can be r a i s e d : "How long does i t take (T_.) before a t the s t a t i o n x = L the r i v e r b e d i s
m m lowered by 50% o f the f i n a l value (thus i/E;)?".
According t o Table 1 the argument of the complementary e r r o r f u n c t i o n o f Eq. (70) has t o be 0.48. Thus
1
(71)
I t i s convenient t o replace T by the number o f years N .
In d e f i n i n g w i t h Eq. (49)
1 year 1 year
Y = / K ( t ) d t i / s ( t ) d t ( 7 2 ) 0 0
i n which B denotes the w i d t h of the r i v e r . Note t h a t the l a s t i n t e g r a l of Eq. (72) denotes the y e a r l y sedimenttransport.
Hence w i t h Eqs. (71) and (72)
In Table 2 the morphological timescale f o r some r i v e r s i s given. The c h a r a c t e r i s t i c l e n g t h has been taken 200 km.'.This i s necessary t o make the p a r a b o l i c model a p p l i c a b l e . Note from Table 2 t h a t i n some cases i n f a c t has s t i l l been chosen too s m a l l .
Remarks
( i ) The r i v e r s mentioned i n Table 2 show a considerable d i f f e r e n c e i n timescale. The Danube example shows t h a t f o r one r i v e r the morpho-l o g i c a morpho-l time-scamorpho-le can d i f f e r from pmorpho-lace t o p morpho-l a c e ; notabmorpho-ly because o f a change i n slope and g r a i n s i z e .
( i i ) I t i s c l e a r t h a t the accuracy o f the N^ values i s d i r e c t l y governed by the accuracy by which the y e a r l y sedimenttransport i s given.
RIVER STATION (approx. d i s t a n c e t o sea) D mm i KlO~^ 3 a / i km N m c e n t u r i e s Rhine (Netherlands) Zaltbommel 2 1.2 100 20 Magdalena (Columbia) Puerto B e r r i o (730 km) 0.33 5 30 2 Danube (Hungary) DunaiPemete (1826 km) 2 3.5 40 10 Danube (Hungary) Nagymaros (1695 km) 0.35 0.8 180. 2.6 Danube (Hungary) Kunaujvaros (1581 km) 0.35 0.8 180 .1.5 Danube (Hungary) Baja (1480 km) 0.26 0.7 210 0.6 Tana (Kenya) Bura (230 km) 0.32 3.5 50 2 Apure (Venezuela) San Fernando 0.35 0.7 200 4.4 Mekong ( T h a i l a n d ) Pa Mong 0.32 1.1 270 1.3 para Serang (Indonesia) Godong 0.25 0.25 50 2.0
Table 2. Morphological time-scales ( a f t e r de V r i e s , 1975) f o r L = 200 km. m
This accuracy seems r e l a t i v e l y l a r g e f o r the Rhine and the Danube, because the values are based on l a r g e amounts o f f i e l d data. For the Magdalena River the y e a r l y t r a n s p o r t has been e s t a b l i s h e d v i a the Engelund-Hansen formula. . This formula was t e s t e d f o r t h i s
r i v e r e x t e n s i v e l y (MITCH, 1973). For Tana R, Apure R. and Mekong R. a l s o the Engelund-Hansen formula was applied because they have almost the same g r a i n s i z e as the Magdalena River,
( i i i ) Two r i v e r s viz the Rhine-near Zaltbommel and the Danube River near Dunarem-ete seem t o be very slow. I t i s a question whether i n t h i s case the r i v e r s are even i n e q u i l i b r i u m ; r e l a t i v e l y frequent human i n t e r f e r e n c e forces the r i v e r b e d t o change.
2.5. Hyperbolic model.
The disadv aitage o f the p a r a b o l i c model viz only v a l i d a t r e l a t i v e l y l a r g e d i s t a n c e s , has l e d t o a deduction w i t h o u t assuming uniform f l o w from the beginning (Vreugdenhil and de V r i e s , 1973).
A h y p e r b o l i c model has been obtained.
The general expression f o r K i s r a t h e r elaborate here but a f t e r l i n e a r i -z a t i o n again K i s given by Eq. (48) or ( 4 9 ) .
U n f o r t u n a t e l y Eq. (74) can only be derived f o r a constant discharge. Hence, c o n t r a r y t o t h e p a r a b o l i c model the h y p e r b o l i c model cannot be used t o
study the e f f e c t of the whole r i v e r regime.
Moreover, Eq. (74) can only be solved a n a l y t i c a l l y f o r a r e s t r i c t e d number of boundary c o n d i t i o n s . Even i n the case o f a success, the a n a l y t i c a l s o l u -t i o n i s So elabora-te -t h a -t a compu-ter i s r e q u i r e d -t o ge-t i n s i g h -t i n -t o -the s o l u t i o n . I t i s then o f course more a t t r a c t i v e t o solve the d i f f e r e n t i a l equations d i r e c t l y n u m e r i c a l l y w i t h o u t the r e s t r i c t i o n s present i n the d e r i v a t i o n o f Eq. ( 7 4 ) .
c 3x9t K 8 ^
3. Numerical models. 3.1. I n t r o d u c t i o n .
S o l u t i o n o f the basic equations f o r degradation and aggradation n u m e r i c a l l y can lead t o a t t r a c t i v e i n f o r m a t i o n f o r r i v e r e n g i n e e r s i n f o r e c a s t i n g the behaviour o f the b e d l e v e l z(.x,t).
Here the numerical s o l u t i o n w i l l be discussed f o r the case o f mode-r a t e Fmode-roude numbemode-rs which mode-r e q u i mode-r e s accomode-rding t o pamode-r. .1.3. the
solu-t i o n of u - ^ + g ^ = R d f 8u 8z du 8x 8t (75) (76)
However, before d i s c u s s i n g numerical models, i t i s necessary t o dis-cuss a t some l e n g t h the boundary c o n d i t i o n s t h a t may apply.
3.2. Boundary c o n d i t i o n s .
The equations (75) and (76) have t o be solved over an i n t e r v a l 0 < X < L. to
region
of
influence
F i g . 6. x - t diagram. A number o f boundary c o n d i t i o n s i s r e q u i r e d . ( i ) Initial condition. At t = 0 the s i t u a t i o n , notably z ( x , 0 ) . has t o be known. ( i i ) Downstream condition. At x = L f o r any t the w a t e r l e v e l has t o be known ( r a t i n g c u r v e ) . Obvious-l y when the b e d Obvious-l e v e Obvious-l a t x = L changes also the r a t i n g curve changes. I n p r a c t i c e t h i s i m p l i e s t h a t the downstream boundary has t o be placed so f a r downstream t h a tchanges a t z ( L , t ) do n o t take place w i t h i n t h e time o f i n t e r e s t .
( i i i ) Upstream boundary. F i r s t of a l l a t x = 0 t h e value of q ( t ) has t o be known. Together w i t h the downstream boundary c o n d i t i o n t h i s i s s u f f i c i e n t t o solve f o r s u b c r i t i c a l f l o w t h e d i f f e r e n t i a l equation f o r the backwatercurve ( i n
f a c t Eq. (75) when the b e d l e v e l z ( x , t ) i s known).
Moreover at x = O i n f o r m a t i o n on the sedimenttransport has t o be a v a i l a b l e : s ( O j t ) has t o be known. I t can e a s i l y be seen t h a t i n -stead o f s ( 0 , t ) i t i s a l s o s u f f i c i e n t t o know z ( 0 , t ) . I n the l e t t e r case a ( 0 , t ) and thus u ( 0 , t ) i s known. Via s = s ( u ) again s ( 0 , t ) i s known.
The r e q u i r e d boundarycondition s ( 0 , t ) or z ( 0 , t ) , however, creates great d i f f i c u l t i e s . Only seldom an accurate boundary c o n d i t i o n can be given. From F i g . 6 i t can be seen t h a t any e r r o r i n the upstream boundary i s propagating i n s i d e the x t diagram. The r e g i o n o f i n f l u -ence i s bordered by the c h a r a c t e r i s t i c f o l l o w i n g from c^ and through ( 0 , 0 ) . For instance f o r t > t ^ the r e s u l t s a t XQ become a f f e c t e d . I t i s now p o s s i b l e t o make use o f the f a c t t h a t c i n many cases i s
o very small.
I f Cg = Ikm/a then t h e b e d l e v e l a t x^ can be forecasted an e x t r a
year i f the upstream boundary i s moved upstream by only.one k i l o m e t e r . I t i s c l e a r t h a t t h i s i s more e f f e c t i v e f o r slow r i v e r s than f o r
fast rivers.'
Special a t t e n t i o n has t o be paid t o places where d i s c o n t i n u i t i e s
For reasons t o be explained l a t e r i t i s convenient t o place an
inter-nal boundary a t the d i s c o n t i n u i t y
i f the place o f t h i s c o n t i n u i t y i s known and f i x e d .
This can be explained by the case of an i n t a k e a t x = x^ ( F i g . 7 . ) .
For the sake of s i m p l i c i t y the r i v e r d i s c h a r g e Q and the discharge through the i n t a k e are supposed t o be constant. An i n t e r n a l boundary i s placed at X = XQ and the basic equations are solved over the i n t e r v a l s 0 < x < x and X < X < L using the c o n d i t i o n s present a t x = x .•
i n the b e d l e v e l are present.
a
Q
Q ( t ) Q ( t ) - A Q
0 + > Q _ S+= S.
I t r e q u i r e s some consequent reasoning t o f i n d the c o n d i t i o n s at XQ. A l o g i c a l sequence has t o be f o l l o w e d . Marking parameters j u s t upstream of XQ w i t h + and the ones j u s t downstream w i t h -, the f o l l o w i n g can be stated ( i ) Dïsohavge Q ( i i ) Transport S ( i i i ) Velooity u ( i v ) Depth ( v ) Waterlevel ( v i ) Bedlevel z_, Aq -> discontinuous S_ ; no withdrawal o f sediment u_ ; due t o s = s ( u ) and ( i i ) a_ ; due t o ( i ) and ( i i i ) and
a_ + z_ ; no h y d r a u l i c jump -> continuous
z_ ; due t o ( i v ) and ( v ) discontinuous
continuous continuous
g ^ discontinuous
3.3.1 Method o f c h a r a c t e r i s t i c s .
The s o l u t i o n o f the basic equations f o r moderate Froude-numbers (Eqs. 31 and 32) can i n p r i n c i p l e be c a r r i e d out by means o f the method of c h a r a c t e r i s t i c s D e f i n i n g
2 G = u - gq/u
Eqs. (31) and (32) together w i t h the expressions f o r the t o t a l d i f f e r e n t i a l s du and dz from the system.
0 G 0 g » 3u/3t R 0 f u 1 0 3u/3x _ 0 d t dx 0 0 3z/3t du 0 0 dx 3z/3x dz
The main determinant of the system (78) set equal zero y i e l d s d t = 0 and Gdx + f ^ g d t = 0 °^ f - ^ - u^ \{JU °(3) " d t " u - gq/u2 " 1 - f2 as found e a r l i e r . (78) (79) (80)
According t o Cramer's r u l e a l l other determinants have t o be zero t o a r r i v e a t indetermancy.
This leads t o the r e l a t i o n s v a l i d along the c h a r a c t e r i s t i c s i n the x - t plane. These r e l a t i o n s are ordinary d i f f e r e n t i a l equations. One c h a r a c t e r i s t i c r e l a t i o n can be found d i r e c t l y because along dt = 0 the ( o r d i n a r y ) d i f f e r e n t i a l equation o f the backwatercurve i s v a l i d du dx dz dx u j u C^a f o r d t =- 0
The c h a r a c t e r i s t i c r e l a t i o n along the other system o f characteris-t i c s ( characteris-t h e ones f o l l o w i n g from c^) can f o r inscharacteris-tance be found from
d t (Rdx - gdz) = 0 0 R 0 0 0 1 0' d t du 0 0 0 dz d t dx (81) (82) Thus f o r dt ^ 0: dz d f u l u dx ^
c V -
dt d t dx u a 1 (83) For F « 1 t h i s reduces t o : f o r d t dx q dz dx t (84) dx d t u ^ f u q dé d t u ^ f • u " C^q (85)Early computations have indeed been c a r r i e d out by means o f the method of c h a r a c t e r i s t i c s . Although the method i s a t t r a c t i v e w i t h respect t o accuracy, there are large disadvantages.
This i s r e l a t e d t o the presence o f " s h o c k l i k e " f e a t u r e s .
At a d i s c o n t i n u i t y f o r instance i n the b e d l e v e l z a t x = x^ the c e l e r i t y i s d i s c o n t i n e o u s . The o r i g i n a l d i f f e r e n t i a l equations are not v a l i d . Two s i t u a t i o n s are p o s s i b l e
( i ) The d i s c o n t i n u i t y i s present a t a f i x e d place. An example i s given i n Fig. 7. At the i n t a k e z_^ < z_ the i n t e r n a l boundary a t x =
-c,>c
t i e s are met.
( i i ) The d i s c o n t i n u i t y moves as a f u n c t i o n of time. I n t h i s case d u r i n g the computations the l o c a t i o n o f the shock has t o be sought and . then s h o c k f i t t i n g can take place. This r e q u i r e s a good d e a l o f
bookkeeping and i n t e r p o l a t i o n i n the computerprogram. An automat i c procedure by means of a d i f f e r e n automat meautomathod seems more a p p r o p r i -ate (see par. 3.4.).
An example o f a moving shock i s given i n Fig. 8.
A dredged trench perpendicular t o the main stream w i l l g r a d u a l l y be f i l l e d
w h i l e moving downstream.
Note t h a t under the assumptions made i n par. 1.4. only the case o f predominant bedload can be considered here.
At the upstream boundary c_^ > c , t h i s i m p l i e s t h a t the slope w i l l be here as steep as possible ( n a t u r a l s l o p e ) . Fig. 8. Dredged t r e n c h .
At t h e downstream boundary c_^ < c^. Here the slope w i l l g r a d u a l l y become more g e n t l e .
This i s s i m i l a r t o what i n gasdynamics i s c a l l e d an expansion Wave.
3.4. F i n i t e d i f f e r e n c e s . c,<c.
The wish t o have automatic s h o c k f i t t i n g has l e d t o the use of f i n i t e difference-methods (Vreugdenhil and de V r i e s , 1967).
The method used a t present (1976) i s s l i g h t l y d i f f e r e n t from the method published o r i g i n a l l y f o r reasons which become c l e a r l a t e r .
I t has t o be r e c a l l e d t h a t the f o l l o w i n g system o f equations has t o be solved: du dx dz _ u^ dx ~ ^ (86) and 9z _ _ d f 8u 9t " du ' 9x (87)
I f t h e boundary c o n d i t i o n s are s t a t e d p r o p e r l y , then the computations can be c a r r i e d out w i t h the f o l l o w i n g a l t e r n a t i n g steps
Step I Compute u with Eq. (86) for known z (and hence dz/dx)
Step II Compute z with Eq. (87) for known u (and hence du/dx)
The f i r s t step does not r e q u i r e f u r t h e r e x p l a n a t i o n . Any s u i t a b l e com-puterprogram f o r backwaterciarves w i l l do.
The second step needs some a t t e n t i o n . F i r s t o f a l l i t has t o be r e c a l l e d t h a t i t regards here a conservation law (vis the equation o f c o n t i n u i t y of t h e sediment).
Equation (87) can be w r i t t e n as
(88)
A m o d i f i e d Lax scheme can be used t o solve t h i s equation n u m e r i c a l l y . F o l l o w i n g the n o t a t i o n o f F i g . 9 both 9z/9t and 9f/9x are replaced by d i f f e r e n c e s . The f o l l o w i n g
difference equation i s used t o r e
-present the differentialequation of Eq. ( 8 8 ) .
X
Fig. 9. D i f f e r e n c e scheme.
\ - Q Oiz^ + ( 1 - a)z^ + 1 azg} f ^ - f
(89)
The r e l a t i o n between Eqs. (88) and (89) can be seen from Taylor expan-sions arround the p o i n t ( x , t ) .
For instance
(90)
Together w i t h the expansions o f z^ and Zg t h i s leads t o the f o l l o w i n , approximation f o r t h e f i r s t term o f Eq. ( 8 8 ) .
<V', 9,2
ff ^ If + i ( A t ) ^ 1^9^2 + O { ( A t ) 2 } - i a ( A x ) 2 P9^2 -1- O {(Ax)3(At)--'} (91) >-l.
I n a s i m i l a r way
^3 " ^ 1 9 f
2Ax 0 { ( A x ) 2 } (•92)
Combination o f Eqs. (91) and C92) y i e l d s
• l A t 9 z 9 f
F t 9x
£ z ^ (Ax)2 9^2 9 t ^
-i f the h-igher order terms are neglected.
The terms of the r i g h t h a n d side can be combined as follows-According t o Eq. (35)
(93)
| | + c | ^ = R ^
9 t 9x g (94)
I f i t i s now assumed t h a t locally c i s constant and R can be neglected then 9^2 9 ^ - c 9^2 9 x 9 t - c 9x 92 9t 9x - c 92 9x = c2 )2: 9"P (95) Combination of Eqs. (93) and (95) y i e l d s
M + M = ( A x ) ' 9t 9x " 2At 9^2 93F" (96) Or, w i t h y = c A t / A x 9 2 9 f _ ( A x ) ' 9 t 9x " 2At a - .\x' 9^2 91? (97)
The f o l l o w i n g remarks can be made.
( i ) Apparently the d i f f e r e n c e approach i s consistent^ because f o r Ax -> 0, l e a v i n g Ax/At = c/y constant the differential equation i s reached from the difference equation.
( i i ) With regard t o the stability t h e f o l l o w i n g remarks can be made. Following the procedure o f Von Neumann (see Richtmyer and Morton, 1967) the c o n d i t i o n f o r s t a b i l i t y i s
y ' < a < l (98) This c o n d i t i o n i s reached i n assuming t h a t a t the l e v e l t an
e r r o r i n z and u i s present. These e r r o r s are developed i n F o u r i e r s e r i e s o f which one term i s
z ( x , t ) = z ( t ) e x p { i k x } u ( x , t ) = G ( t ) e x p { i k x }
At the l e v e l t + A t t h e e r r o r s become z ( x , t + A t ) = e ^ z ( t ) exp { i k x } u ( x , t + A t ) = e^_^ü(t) exp { i k x }
The e r r o r i s not magnified i f e < 1 and
ty e
u
< 1 f o r every k, This leads t o Eq. (98) as a c o n d i t i o n o f s t a b i l i t y .
( i i i ) Apparently Eq. (98) s t a t e s t h a t p o i n t 4 i n F i g . 9 has tó be
s i t u a t e d w e l l below the c h a r a c t e r i s t i c through p o i n t 1. This could be expected as t h e d i f f e r e n c e scheme i s explicit.
( i v ) For a = 1 the d i f f e r e n c e scheme i s s i m i l a r t o t h e one Lax (1954, 1960) a p p l i e d f o r t h e equation o f motion f o r an i d e a l f l u i d . I n t h a t case the r i g h t - h a n d member o f the d i f f e r e n c e equation i s l i k e a v i s c o s i t y term. This i s why t h i s d i f f e r e n c e scheme i s sometimes r e f e r r e d as a pseudo-visaosity approach. Note t h a t the e x t r a term i s o f a numerical nature and n o t o f a p h y s i c a l one!
(v) The r i g h t h a n d side o f Eq. (97) introduces some numevioal damping In order t o a r r i v e a t accurate r e s u l t s , i t i s important t h a t t h i s term i s small compared t o the other terms o f the equations. This damping can be reduced i f Ax i s selected s m a l l and/or i f 3 = a -i s s m a l l . Th-is -i s e s p e c -i a l l y the case -i f 3^z/9x^ -i s l a r g e .
The above given d i f f e r e n c e scheme does take care o f automatic s h o c k f i t t i n g , t h e r e f o r e i n p r i n c i p l e t h e r e i s no'need t o use i n t e r n a l boundaries f o r cases l i k e i n F i g . 7. However, nevertheless these i n t e r n a l boundaries are used. The reason can be made c l e a r from Eq. ( 9 7 ) . At places where s h o c k l i k e f e a t u r e
are present automatic s h o c k f i t t i n g leads t o l o c a l l y l a r g e values o f 9^z/9x^. Hence the r i g h t h a n d side value o f Eq. (97) can be l a r g e . Thus a t these shocks t h e numerical damping i s r e l a t i v e l y l a r g e . Therefore i f the l o c a t i o n of the shock i s knovm, l i k e i n the case of the i n t a k e , i t i s s t i l l a t t r a c t i v e t o use an i n t e r n a l boundary.
A disadvantage i n p r a c t i c e i s the high degree o f n o n - l i n e a r i t y o f the t r a n s p o r t f u n c t i o n . Small d i f f e r e n c e s i n b e d l e v e l (and thus i n u) cause a r e l a t i v e l a r g e change i n s (and hence i n c ) . In order t o a r r i v e a t s t a b i l i t y the c o n d i t i o n o f Eq. (98) has t o be fulfilléd i n the e n t i r e i n t e r v a l 0 < x < L. Thus the value o f At, once Ax has been s e l e c t e d , has t o be based on the maximum value o f c present.
-2 -3
For 3 a value o f 10 or 10 i s selected and a i s selected a c c o r d i n g l y . I t i s not a t t r a c t i v e t o vary a w i t h x because then the equation o f c o n t i n u i t y i s not f u l f i l l e d anymore.
In p r a c t i c e i t i s not easy t o avoid numerical damping and nevertheless have an a t t r a c t i v e speedy computation.
Other d i f f e r e n c e schemes have been suggested i n l i e u o f the one d i s -cussed here; For instance Perdreau and Cunge (1971) proposed an upstream-d i f f e r e n c e scheme; t h i s one appears, however, not t o be more accurate.
'3.5. A p p l i c a t i o n s .
Figure 10 shows a case o f a r i v e r l o c a l l y c o n s i s t i n g o f two channels, surrounding an i s l a n d . The depth i s r e s t r i c t e d ; t h e r e f o r e plans have been made t o close one channel i n order t o increase the depth i n the other channel. The question i s now r a i s e d a t which speed the r i v e r w i l l a d j u s t i t s e l f t o the new s i t u a t i o n or whether i t i s necessary t o c a r r y out some dredging i n order t o o b t a i n s u f f i c i e n t depth f o r n a v i g a t i o n soon.
The computations show t h a t indeed dredging i s r e q u i r e d . Otherwise the investments f o r the c l o s i n g o f f w i l l only be e f f e c t i v e a f t e r some years.
The f o l l o w i n g remarks can be made:
( i ) A s e l e c t i o n Ax = 250 m i s combined w i t h a time step At = 5 days. This means t h a t the use o f an expHoit scheme h a r d l y hampers the s e l e c t i o n o f a l a r g e time step. A l a r g e r t i m e s t e p , p o s s i b l e
O- TIME IN YEARS
t
12000-E
DISTANCE IN km
the s e l e c t i o n o f a large time step. A l a r g e r time s t e p , possible f o r an i m p l i c i t scheme would not be favourable. This i s due t o the f a c t t h a t then the boundary c o n d i t i o n Q ( t ) cannot be i n t r o -duced c o r r e c t l y .
• ( i i ) The r e s u l t s are obtained by i n t r o d u c i n g a regime and not a constant discharge,
( i i i ) The computations are based on the formula o f Meyer-Peter and Mueller, supposed t o be a p p l i c a b l e f o r the r e l a t i v e l y coarse bedmaterial. Sofar the t r a n s p o r t f u n c t i o n s = f ( u ) has been used t o make the deductions transparent. However, f o r r e a l computations a r e a l i s t i c t r a n s p o r t f o r m u l a has t o be used, a p p l i c a b l e t o the r i v e r considered.
The next case considers morphological computations c a r r i e d out f o r the B i f u r c a t i o n Pannerden ( F i g . ! ! ) •
The question was r a i s e d t o what extent a h i g h f l o o d on the Rhine might lead t o serious erosion on the upper p a r t o f the Pannerden Channel and hence might lead t o a t e m p o r a r i l y l a r g e c a p a c i t y o f t h i s r i v e r . This might lead t o the s i t u a t i o n t h a t d u r i n g a second high f l o o d soon a f t e r the f i r s t one the discharge through the Pannerden Channel might become so l a r g e t h a t the downstream embankments might bê too low. The p o s s i b i l i t y o f erosion becomes clear from the f l o w p a t t e r n . At two places water coming from the high waterbed (and c a r r y i n g h a r d l y any sediment) i s e n t e r i n g the low waterbed i . e . i n c r e a s i n g l o c a l l y the sedimenttransport by erosion of the bed.
The mathematical model b u i l t f o r t h i s problem represented the (mobile) low waterbed whereas the discharges coming from the high waterbed were introduced as boundary c o n d i t i o n s . The model was c a l i b r a t e d by means o f scarce observations made d u r i n g and a f t e r the high f l o o d of 1926. The design f l o o d ( F i g . 11) was used t o compute b e d l e v e l v a r i a t i o n s along the Pannerden Channel.
The r e s u l t s are represented i n F i g . 11 as f l u c t u a t i o n s arround the i n i t i -a l b e d l e v e l . Indeed serious scour occurs ne-ar s t -a t i o n s 869 -and 870 where
water without.sediment i s d i s c h a r g i n g intö the low water bed. However, the r e l a t i v e l y low speed a t which morphological changes i n the Rhine branches take place i s c l e a r l y demonstrated. The eroded sediment s e t t l e s t e m p o r a r i l y j u s t downstream of the places w i t h erosion. This means t h a t the e f f e c t o f the inovease o f the discharge c a p a c i t y i s
/
only a l o c a l e f f e c t . I t i s balanced by a local deorease and hence dovmstream o f the f i r s t few k i l o m e t e r s o f the Pannerden Channel no r e s u l t i n g increase o f the discharge capacity has t o be f e a r e d .
F i n a l l y a case w i l l be discussed a t some l e n g t h i n which the use o f the morphological computations i n the approach o f r i v e r problems i s demonstrated. I t regards here the Tana River (Kenya). I n the framework o f a new i r r i g a t i o n p r o j e c t a f i x e d weir has been designed. Part o f the r i v e r d i s c h a r g e Q ( t ) i s used f o r i r r i g a t i o n .
The sediment i s supposed n o t t o enter the i r r i g a t i o n s y s t e m and i f so be caught by means o f a sandtrap and from time t o time f l u s h e d i n t o the r i v e r again. Hence the r i v e r downstream o f the w e i r f i n a l l y has t o c a r r y the same sediment w i t h a smaller discharge.
Moreover, a f t e r c o n s t r u c t i o n o f the weir sedimentation w i l l take place upstream o f the w e i r . This may i n the beginning
lead t o erosion downstream o f the weir due t o a smaller sediment supply t o t h i s reach and i n s p i t e o f the r e d u c t i o n o f the d i s -charge.
Before s t a r t i n g a time and money consuming numerical model, some basic questions can be answered r a t h e r simply.
( i ) W i l l the f i n a l changes o f the r i v e r be so l a r g e t h a t d e t a i l e d computations are necessary?
( i i ) W i l l the p e r i o d i n which these changes take place be l a r g e or s m a l l compared t o the l i f e t i m e o f the p r o j e c t ?
Thus as a f i r s t approach some rough approximations w i l l be made. N a t u r a l l y from the Tana River no long s e r i e s o f observations on water- and sediment movement are present.
Jan. m s May 455 Sept. 80
Febr. 94 June 223 Oct. 118
March 102 J u l y 125 Nov. 379
A p r i l 262 August- 95 Dec. 342 Table 3. Average monthly discharges a t Garissa (m / s ) .
At Garissa, upstream o f the dam s i t e the average monthly d i s -charges are known. As the discharge f l u c t u a t e s n o t very much d u r i n g a month (lower r i v e r ) the data o f Table 3 can be con-sidered a f a i r estimate o f the p r o b a b i l i t y d i s t r i b u t i o n o f the discharge.
Some c h a r a c t e r i s t i c s o f the r i v e r are given: -4
w i d t h B = 80 m; slope i = 3.5 x 10 ; g r a i n s i z e D^^ = 0.32 mm. From data q f other r i v e r s t h e r e i s good reason t o assume t h a t
the t r a n s p o r t f o r m u l a o f Engelund-Hansen (1967) can be a p p l i e d . For the f i n a l erosion and sedimentation i t i s assumed t h a t C-value and g r a i n s i z e remain the same.
This i m p l i e s t h a t the Engelund-Hansen formula can be w r i t t e n as a f u n c t i o n o f discharge and slope only
S ^ . 5/3
I n d i c a t i n g the present slope as i ^ and the f i n a l f u t u r e slope w i t h i ^ , the c o n d i t i o n t h a t the y e a r l y sediment t r a n s p o r t w i l l f i n a l l y be equal t o the present one w i l l lead t o :
Q//3 = (Q,-AQ)^/^ i , ^ / ^ (91) i = l ^ i = l
i n which denote the monthly discharges.
As the o v e r a l l slope o f a r i v e r does not change w i t h the d i s c h a r -ge i t f o l l o w s
12 Z Qi 1=1 5/3 12 I (Q^-AQ) 5/3
From the a c t u a l data f o r the proposed AQ f o r t h i s i r r i g a t i o n scheme i t appeared t h a t i-^/iQ = 1.045 t o 1.055.
This small value has t o be combined w i t h the i n f o r m a t i o n t h a t the damsite i s s i t u a t e d about 225 t o 230 km from the ocean. Hence the f i n a l sedimentation downstream o f the weir can amount t o 0(4m), t h a t i s o f the same order o f magnitude as the h e i g h t o f the proposed w e i r .
An almost s i m i l a r approach i s possible f o r e s t i m a t i n g the f i n a l sedimentation upstream o f the w e i r . The u l t i m a t e depth ( a ) o f the r i v e r upstream o f the weir i s l i n k e d t o the o r i g i n a l h e i g h t o f the weir ( a ^ ) the waterdepth above the c r e s t o f the weir (AH) and the f i n a l sedimentation (Aa) according t o
a. = a. + AH. - Aa i l l
The Engelund-Hansen formula has now t o be expressed i n terms o f Q and a according t o
S 'v. Q^a ^
F i n a l e q u i l i b r i u m means now
ZQ.^ a.:^ = EQ.^a, + AH. - Aa)"^ (
1 l O 1 1 1
Using Eq. (93) and expressing AH as a f u n c t i o n o f Q by means o f the c h a r a c t e r i s t i c s o f the weir t h i s leads t o an i m p l i c i t equation f o r Aa. This equation can e a s i l y be solved e.g. by means o f the
regula falsi. For the Tana case a^ = 4m led t o Aa 3m.
These two answers: downstream sedimentation 0(4in) and upstream sedimentation 0(3m) lead t o t h e conclusion t h a t d e t a i l e d computa-t i o n s are necessary. This i s computa-the more computa-the case as according computa-t o Table 2 t h e Tana River r e a c t s r e l a t i v e l y quick.
In t h i s case numerical computations had t o be c a r r i e d out w i t h scarce data. Again t h e v a r i a t i o n s o f the b e d l e v e l have been com-puted w i t h respect t o the present s i t u a t i o n .
The f o l l o w i n g r e s u l t s were obtained f o r a weir w i t h a h e i g h t o f 4 m. ( i ) Upstream o f the weir the b e d l e v e l w i l l r i s e f i n a l l y 2 m i n
a p e r i o d o f about 5 years.
( i i ) Just downstream o f the weir w i t h i n t h e f i r s t year a f t e r the c o n s t r u c t i o n t h e l a r g e s t , erosion (4 t o 5m) can be present. ( i i i ) This erosion w i l l l a t e r be f o l l o w e d by sedimentation o f some
meters w i t h respect t o the present b e d l e v e l . However, i t w i l l take almost a centuary before t h i s i s the case.
The degree o f erosion and sedimentation i s s t r o n g l y i n f l u e n c e d by the h e i g h t o f the w e i r .
3.6. Comments.
The above given methods on morphological computations s t i l l have some serious draw-backs. However, i f a p p l i e d c a r e f u l l y they can give some important i n f o r m a t i o n i n the design o f riverimprovements.
As i s demonstrated by the three examples o f par. 3.5. the methods have been used as t o o l s t o answer s p e c i f i c questions. I t cannot be ex-pected t h a t the computational methods, based on r a t h e r crude sche-m a t i z a t i o n s o f a cosche-mplex p h y s i c a l phenosche-menon, can give a d e t a i l e d p r e d i c t i o n o f the bedtopography o f a r i v e r .
Much research w i l l have t o be c a r r i e d out before t h i s u l t i m a t e g o a l i s achieved. I n the mean time the mathematical models a v a i l a b l e a t present serve t o the p r e d i c t i o n o f some tendencies o f morphological behaviour o f r i v e r beds.
Combined w i t h the use o f aoale models, capable o f p r e d i c t i n g much more d e t a i l a l b e i t f o r r e s t r i c t e d reaches, i t i s p o s s i b l e t o f o r e c a s t mor-p h o l o g i c a l changes i n r i v e r b e d s due t o human i n t e r f e r e n c e .
4. Suspended load t r a n s p o r t . 4.1. I n t r o d u c t i o n .
The morphological computations discussed i n the previous Chapters d e a l w i t h the case i n which the i n f l u e n c e o f suspended load can be neglected. I n t h i s Chapter t h e e x t e n t i -on o f the p o s s i b i l i t i e s f o r these computati-ons are discussed. Therefore the theory o f suspended load both f o r uniform and non-uniform f l o w are discussed'in par. 4.2., f o l l o w e d by remarks on morphological computations i n par. 4,3.
4.2. Theory.
The c l a s s i c a l approach f o r suspended load i s based on the work o f H. Rouse f o l l o w i n g a suggestion o f Von Karman.
A (dynamic) e q u i l i b r i u m i s supposed t o be present: p a r t i c l e s tend t o s e t t l e due t o g r a v i t y ( f a l l v e l o c i t y W) and they are t r a n s p o r t e d upwards due t o t u r b u l e n t d i f f u s i o n .
The equation o f c o n t i n u i t y can be set f o r the c o n c e n t r a t i o n ( c ) .
. c W f s ^ | | = 0 ( 9 6 )
Here z i s used, d i f f e r e n t from the previous Chapters, f o r the v e r t i c a l o r d i n a t e , z^ being the b e d l e v e l now.
The d i f f u s i o n c o e f f i c i e n t f o r the sediment w i l l i n p r i n c i p l e be d i f f e r e n t from the one f o r momentum ( e ) according t o
e = KU z ( l - .^//z/^) ( 9 7 )
In p r a c t i c e e J^i^ e i s s e l e c t e d , which means t h a t Eqs. ( 96) ^ d ( 97 ) can be combined. I n t e g r a t i o n y i e l d s
Z . _^ fa - z"l
c = constant (98)
w i t h Z = W/KU . For e = ge the ( c o n s t a n t ) f a c t o r 3 i s a l s o i n t r o
-X s
duced i n Z.
Hunt (1954) has shown t h a t i n f a c t Eq. ( 96 ) should read
( 1 - c)cW + e ^ = 0 ( 9 9 ) dz
The concentration d i s t r i b u t i o n o f Eq. (98 ) i s a r e l a t i v e one. An i n t e g r a t i o n c o n s t a n t has t o be known. Since E i n s t e i n (1950) i t has been common p r a c t i c e t o s e l e c t t h i s i n t e g r a -t i o n cons-tan-t a -t a l e v e l close -t o -the bed. E i n s -t e i n de-r i v e d the c o n c e n t de-r a t i o n neade-r the bed fde-rom h i s bedload equa-t i o n . However, equa-there are some imporequa-tanequa-t argumenequa-ts equa-t o be made against t h i s procedure:
( i ) Near the bed the presence o f the bedform ( r i p p l e s or dunes) becomes important.
( i i ) Near the bed t h e concentration i s n o t s m a l l ; hence Eq. ( 99 ) should be used.
( i i i ) For r e l a t i v e l y l a r g e concentrations the f a l l v e l o c i -t y i s i n f l u e n c e d by -the c o n c e n -t r a -t i o n .
A good r e l a t i o n i s
^ = ( 1 - c ) ^ - ^ ^ f o r 1 < Re = 1^ < 500 (100)
( i v ) According t o Eq. ( 9 7 ) apparently e(0) ~ 0 l e a d i n g t o C(0) =00,
This cannot be t r u e because there i s exchange o f s e d i -ment a t the bed, moreover c i s equal the value f o r
max ^ l o o s e l y packed sand.
T r i a l s have been made t o d e r i v e e from observations. The obser-s
v a t i o n s o f Coleman (1970) have been reproduced i n Figs. 12 and 13.
In F i g . 12 also the values o f e according t o Eq. (96) have been p l o t t e d . The data show a l a r g e s c a t t e r . This i s due t o the f a c t t h a t e has been derived from Eq. (96) by differentiating measu¬ red values o f c. However, some trends are c l e a r
( i ) The parameter W/u seems t o i n f l u e n c e .
X s
( i i ) Near the watersurface e / e . s
I t has t o be n o t i c e d t h a t Coleman d i d h i s observations i n a l a b o r a t o r y flume w i t h fixed bed. Hence .his data on e near t h e
s
bed are not r e a l i s t i c . This i s not the case f o r t h e data from the Enoree River ( F i g . 13). Here, however, as i t regards f i e l d data the accuracy cannot be l a r g e e i t h e r .
1.0 0.4 a i 0.06 0.04 0.01 a o o e 0.004 VON KARMAN POWER LAW LOGARITHMIC A / u . a i o i o 0 4 6 2 S a i 7 9 (D a s 3 e @ o j o e e a s s s 9 0291 ® 0 6 9 6 @ a 3 4 2 O o e i o ® a 3 7 e s OS07 ©
AFTER COLEMAN 1970 (MEASUREMENTS ANDERSON 1942)
t i o n r e s p e c t i v e l y ; are the momentary f l o w v e l o c i t i e s . This leads t o
f -^lïï ^ " l ^ > ^ |r ^"2^> + = 0 (101)
D e f i n i n g = u^ + u^' and C = c + c' then averaging i n time leads t o
|£ + |_ {u^e} +
1^
{u^c} + ^ { U 3 C - Wc) +I t i s assumed t h a t the averaging takes place over period 6 l a r g e enough t o get r
I t can now be defined
l a r g e enough t o get r i d o f the t u r b u l e n t f l u c t u a t i o n s .
- u. ' 'c' = 1 'c' = 2 - 'c' = 3 l l ' ^ x + ^12^y + ^13^z 21'^x + ^22^y + ^23^z (103) Sl'^x + ^32^y + ^33°z
i n which the s u b s c r i p t s x, y and z denote d i f f e r e n t i a t i o n . For the two-dimensional case i n the x,z-plane a l l d e r i v a t i v e s w i t h respect t o y disappear.
Hence, combination o f Eqs. (102) and (103) thus gives
- I? <^31=x * =33\ï = °
The c o n t i n u i t y equation f o r t h e f l u i d reads here
3u 3u
3 ^ + 3 ^ = ° (105)
I t i s common t o assume t h a t the axes o f the coordinate system are taken along the p r i n c i p a l axes o f the d i f f u s i o n tensor.
Hence the d i f f u s i o n tensor ^ reduces t o the s c a l e r (see e.g. Graf, 1971, p. 164).
Thus
9c ^ 9c . „ 9 c „ 9c 9 r 9 c , _ 9 ^ r 9c-,_ .-„„v 9 t ^ ^ 9ÏÏ ^3 9 i • ^ 9¥ 9ÏÏ ^"^1 9 l ^ 9^ ^^3 9 2 ^ " ° ^^^^^
Further s i m p l i f i c a t i o n can be obtained as f o l l o w s . ( i ) For n e a r l y uniform f l o w u^ !v 0.
( i i ) For O(e^) = O(e^) i t can be assumed
fe(^ii>« fc fc>
-because the l e n g t h L considered i n x d i r e c t i o n i s much l a r g e r than the depth a. •
Hence
If fc-fc =3 fc' = °
For steady uniform f l o w the concentration does not vary i n x d i r e c -t i o n .
For t h a t case Eq. (109) reduces a u t o m a t i c a l l y t o
9 {We + e_
If}
= 0 (110)9 z 3 dz or
Wc + e_ |2. = constant (111)
As Eq. ( I l l ) expresses the v e r t i c a l f l u x o f the sediment i t can be stated t h a t the i n t e g r a t i o n constant equals zero (no sediment passes the w a t e r s u r f a c e ) .
Thus Eq. ( I l l ) i s s i m i l a r t o Eq. ( 9 6 ) .
In order t o get i n s i g h t i n t o , the adaption o f the concentration d i s t r i b u t i o n f o r the case o f 3c/9x = 0, Kerssens (1974) c a r r i e d out computations f o r Eq. (109) f o r 9c/9t ^ 0 .
Thus
(112)
The s o l u t i o n o f Eq. (112) r e q u i r e s boundary c o n d i t i o n s ,
Hc(oi2) C ( U ) The c o n d i t i o n s are r e q u i r e d f o r c ( 0 , z ) ; c ( x , z ^ ) ; c ( x , a ) and c ( L , z ) . As t h e i n t e n t i o n was t o o b t a i n i n f o r m a t i o n about F i g . 14. Adaption o f c ( x , z ) .
the l e n g t h L over which t h e c o n c e n t r a t i o n d i s t r i b u t i o n adapts i t s e l f t o the e q u i l i b r i u m case i t i s l o g i c a l t h a t c(0,|i) i s chosen as simple as p o s s i b l e : c(0,z) = c^ = constant.
For z = a n a t u r a l l y the c o n d i t i o n reads t h a t t h e v e r t i c a l f l u x i s zero. For X = L the e q u i l i b r i u m - c o n c e n t r a t i o n d i s t r i b u t i o n i s taken.
The boundary c o n d i t i o n near the bed r e q u i r e s some remarks. For bedload
transport i t can be assumed t h a t the transport r e a c t s d i r e c t l y on a
change i n the f l o w c o n d i t i o n s . For suspended load i t i s now assumed
t h a t the aonoentration near the bed reacts d i r e c t l y on a change i n the f l o w c o n d i t i o n s . This concentration c ( x , z ^ ) i s f o r the case o f uniform f l o w ( F i g . 14) ,a constant c^. The value o f c^ can be derived from c ( L , z )
as
J u ( z ) c ( z ) = t o t a l t r a n s p o r t (113)
As f o r X <» both u ( z ) and c ( z ) are known and the t o t a l t r a n s p o r t can be estimated by means o f any adequate t r a n s p o r t f o r m u l a , Eq. (113) contains only one unknown viz c ( z ^ ) .
