An implicit method to solve Saint Venant equations

Journal o f Hydrology, 24 (1975) 171--185 Q North Holland Publishing Company, Amsterdam -- Printed in The Netherlands

AN IMPLICIT METHOD TO SOLVE SAINT VENANT

EQUATIONS

FRANCESCO GRECO and LORENZO PANATTONI IBM Scientific Center, Pisa (Italy) (Received April 29, 1973; accepted April 4, 1974)

ABSTRACT Greco, F. and Panattoni, L., 1975. An implicit method to solve Saint Venant equations. J. Hydrol., 24: 171--185. An implicit numerical method for solving Saint Venant equations has been defined for an application relating to the river Arno. This method exploits the linearity in the discharge of the mass equation, by means of which it is possible to express the discharge as a function of the water level and to use this expression in the equation of motion. Then the resulting non-linear equation for a grid element on the (x, t) plane contains only the water levels as unknown quantities. The solution of the system of equations for the entire reach, by the use of the Newton iteration method, is facilitated by the particular form of the matrix of the coefficients. The rapid rate of convergence and the limited storage allocation are characteristics of this implicit scheme. This method has been compared with other implicit methods based on the same grid of points.

INTRODUCTION The Saint Venant equations constitute a system of two non-linear partial differential equations. When applied to natural or artificial watercourses to describe the phenomenon of flood propagation, these equations cannot be s o l v e d , e x c e p t in v e r y s p e c i a l cases, b y a n a l y t i c a l m e t h o d s ( Y e v d j e v i c h , 1 9 6 4 ) . I t is t h e r e f o r e n e c e s s a r y t o a d o p t n u m e r i c a l m e t h o d s b a s e d o n t h e u s e o f a f i n i t e d i f f e r e n c e s c h e m e ; t h e r e s u l t i n g a l g e b r a i c e q u a t i o n s , in v i e w o f t h e c o m plexity of the phenomenon, must be solved by means of digital computers. V a r i o u s n u m e r i c a l m e t h o d s o f s o l u t i o n , d i f f e r i n g o n e f r o m a n o t h e r in t h e numerical scheme adopted, have already been developed (Strelkoff, 1970). This study describes an original method based on an implicit scheme and a d o p t e d in a s t u d y o n t h e r i v e r A r n o . SAINT VENANT EQUATIONS The Saint Venant equations are a system of two hyperbolic partial differe n t i a l e q u a t i o n s a n d t h e y m a y b e w r i t t e n as f o l l o w s :

172 ~Q 3z --- +B-----q=O 3x ~t

~z

(1)

_1 ( ~Q/A

-- +Sf+ 5x g

---+ Ot

Q ~Q/A ) A

Ox

=0

(2)

in which x = c u r r e n t abscissa (positive in the d i r e c t i o n o f flow); t = time; Q = discharge across a section; B = w i d t h o f w a t e r surface; z = w a t e r level as r e f e r r e d t o a h o r i z o n t a l d a t u m ; g = a c c e l e r a t i o n o f gravity; A = cross-sectional area; Sf = friction head losses; q = lateral in- or o u t f l o w {positive f o r inflow, negative for o u t f l o w ) . While it is possible to refer to t h e classical writings o n hydraulics for the d e d u c t i o n o f these equations, it should briefly be stated t h a t eq.1 expresses the c o n s e r v a t i o n o f mass and eq.2 the c o n s e r v a t i o n o f m o m e n t u m . The Manning f o r m u l a has been used t o express friction head losses: Q2

(3)

Sf = n 2 A2R4/3

in which n is the Manning c o e f f i c i e n t and R t h e h y d r a u l i c radius. In a d d i t i o n to the a p p r o x i m a t i o n s implicitly c o n t a i n e d in eqs.1 and 2, see (Evangelisti, 1 9 6 6 ) f o r e x a m p l e , the t e r m expressing lateral in- and o u t f l o w s has been o m i t t e d in t h e e q u a t i o n o f m o t i o n . FINITE DIFFERENCE SCHEME Fig.1 illustrates the grid in t h e (x, t) plane used for the n u m e r i c a l integration o f eqs.1 and 2; the vertical lines r e p r e s e n t the sections along t h e river at abscisse Xl, x2 . . . . . XN, selected in o r d e r t o o b t a i n a geometrical d e s c r i p t i o n o f the s t r e a m (the distance b e t w e e n t w o c o n s e c u t i v e sections is generally variable); x l is t h e abscissa end section u p s t r e a m and x N the end section d o w n s t r e a m . T h e h o r i z o n t a l lines r e f e r to times. T h o u g h n o t essential, t h e / x t intervals are here regarded as c o n s t a n t .

k~xl-

Xr

X~

X3 ......

Xt

Fig.1. Computational grid in (x, t) plane.

X~.~

....

XM

- X

173 Assuming that the values of the quantities z and Q are known at the time (k -- 1)At, the solution of eqs.1 and 2 consists in the calculation of the same quantities in the next time period kA t. At the initial time, the water levels and discharges (Initial Conditions) are supplied by a backwater computation. Taking as y any of the quantities: Q, A, B, z, Sf, or Q/A, its average value and that of its derivatives in the specific grid element i as defined by points (i, k -- 1), (i + 1, k -- 1), (i, k) and (i + 1, k) may be written as follows: Y:

n

2

+(l--n)

(4)

2

5y ~ ~ (y~÷, _ y ~ ) + .1. .--n . (y,~-,'--y~ _,) ~)x Ax Ax

1

Oy Ot

--

2at

,,_,

(yk

- - Yi÷, + y l [ _ y~[-

(5)

)

(6)

in which ~ is the weight coefficient with which the different variables and their derivatives are averaged in relation to time. In agreement with the results obtained b y other writers (Cunge and Wegner, 1964) it has been seen that the value of V must lie between 0.5 and 1 to ensure the stability of the implicit scheme. The value of V does n o t have any significant effect on the results, provided the above condition is respected. For the sake of simplicity, therefore, the value ~ = 1 has been chosen. By replacing eqs.4, 5 and 6 in eqs.1 and 2, and with reference to the generic grid element i, the following algebraic system of two equations containing the four unknowns z~, Q~, z~÷,, Q~÷, is obtained:

Fi-B~+B~÷'2

Axi (z~--zT-' +z~÷,--z~-')+Q~.,,., 2At Ax,

+

,

Axi

( Q~

'

1

Q~÷,

Q~

[Q~. !

-Q~-

qi=O

(7)

,_

Qik, '

Q~

Q~-' ) +

,

A,+,

A7

A?-'

Q~

in which qi is the lateral inflow in the reach xi+ 1 - - x ; and AX i = xi+ 1. METHOD OF SOLUTION The method of solution adopted is based on the linearity in the discharge e,f the mass equation. It is possible then, by means of this equation, to express

174 at a given time discharges as a f unct i on of the water levels along the reach and of the discharge at a given point of the reach. The shape o f the resulting system of the algebraic non-linear equations slightly depends on the b o u n d a r y conditions chosen. Since the variations from one case to a n o t h e r are easily obtainable, only the case in which the water levels z,(t) at the upstream section and the rating curve f(zN, QN) = 0 at the downstream section are given as b o u n d a r y conditions has been dealt with here. In this case, from eq.7 it is possible to express explicitly the discharge Q7 in the section i at time k A t as a function of the level z~ at the same point and o f values Q/k÷, and z~÷, at the same instant in section i + 1, thus obtaining:

Q7 = QIf+' + BIf +2B~+, Ax i2At (z~ -- z7 -' + z/k+, -- z ai÷,-' ) -- Cli

(9)

F r o m this stage onwards, for the sake of simplicity, the variables w i t h o u t upper indexes are referred to time kA t. By deducing the value of QN as a function of the depth z N in section N f r om the relationship f{ZN, QN) = O, and then applying eq.9 f r om this section onwards to all subsequent grid elements, it is possible to express the discharge in a generic section as a function o f the water levels at the same instant in all sections between this generic section and section N. Thus by replacing the discharge value indicated in eq.8 with the discharge expressed as a f unc t i on o f the water levels, calculated in the manner described, a system of N -- 1 equations for the N -- 1 unknow ns z2, z3 . . . . . z N is obtained: F1 (Z2, Z3 . . . . . .

ZN) = 0

F2 (Z2, Z3 . . . . . .

ZN) = 0

F3 (Z3, Z4,.

ZN ) = 0

(10) F N - 2 ( Z N - 2' Z N - , , ZN ) = O F N , (ZN_ ,, Z N ) = 0

The Newto n iterative m e t h o d has been a d o p t e d to solve system 10 of nonlinear equations. This m e t h o d consists in assigning trial values for u n k n o w n water levels and in calculating the corrections from the following system of linear equations, in which the derivatives are naturally calculated in correspondence with the trim values:

175

5F~ ~F~ OF1 ~F~ dz2 + - -- dz3 + . . . + - dZN_ , + dz N = - - F I Oz2 ~za ~ZN_, ~Z N

~F2 ~F2 _ ~F2 ~F2dz: + dz3 + . . . + - dZN_ 1 + - -- d z N = - - F2 ~ZN_~ ~Z g Oz2 Oz3

__

OF 3

~P3

+~F3

-- dz 3 + . . . + d z N_, ~)z3 ~}z N ,

•

~FN-ldzN_ ' +SFN-ldz N ~ZN_,

(11)

- - dz N =--F3 OZN

.

|

FN

OZN

The coefficients are c o m p u t e d from eqs.8 and 9:

~JFi 5z~+,

-

OGi ~zi÷ ,

~ F i _ ~)Gi ~zi

~}zi

~Gi ~Qi ~)Qi ~zi÷,

+

÷

~Gi 5Qi,,

OQi÷~ (12)

~zi+,

OGi ~)Qi

+ ---

~ Qi

(13)

~zi

and, f o r ] > i + l : ()Fi

--

~Gi

~Qi

+

~Gi

~Qi.~

(14)

From eq.9 it follows:

(15)

~Qi/~z~ = ~Q~.,/~z~

from which it may be deduced that: ~Qi/~zj = ~Q,1/~z i = ...

= ~Qj_~/~zj

(16)

Bearing in mind the last relationship, eq.14 may be written as follows:

~zi

~zj

\~i

+ ~Qi+,

i.e. as a product of two terms, one dependent on the index j in relation to the water level z on the basis of which function F is derived, and the other on the index i in relation to the grid element in which function F is calculated. From this it is deduced t h a t the elements of the coefficients matrix may be grouped together in only four vectors:

176 =- a i

(18)

Ol'i/OZ i

= bi

(19)

OFi/Oz j

=- c i d y

(20)

Ori/Ozi+

'

The coefficients matrix for system 11 is a Hessemberg matrix and, using these symbols, may be written: al

d3cl

b2

d,acl

...

d N lcl

dNCl

a2

d4c2

. • •

tiN_ ~c2

dN¢2

b3

a3

•• •

d N_~03

dNC3

b N -,.

a N -2

d N C N -:

bN

aN

M =

- I

(21)

-

In c o m p a c t terms system 11 may be written as follows: M. Az = F

(22)

in which Az is the vector of the unknow ns dz, and F the vect or of the known terms. This expression is solved by means of the direct factorization m e t h o d (Isaacson and Keller, 1966), with which, by expressing the matrix M as a p r o d u c t of the o t h e r two matrices L and U: M=L.U

(23)

the solution o f system 11 may be achieved in t w o successive steps: L. Y = F

(24)

and U . Az = Y

(25)

where Y is an intermediate solution vector. In the case in question matrices L and U have the following form:

QI

-]

b2 ~2

b3~3 (26)

L =

. bN-i

j ~N-I

177

1 UI2 /-/13 1

U=

UI,N_ t

U23

U?,N_ 1

1

U3, N_

(27)

UN - 2,N -1 1

in which the elements ~ and u may be obtained from the following formulae:

U:j

=Cldj/~.l

~1

=al

uij

= (cid i

Qi

j=2,3 ..... N--1

=ai--biui-,,i

i=2,3

-- biui_,,j)/~ i

.....

j = i + 1,...

(28) t

N--1 , N -- 1

The solution of steps 24 and 25 is immediate in view of the triangular shape of matrices L and U.

l

7,

7~

e.

e-

:

-

/

/

~2

~o

24

o

-

........

3

......

2

~2

'~

1

~e

24

I

6

~2

T I M [ - HOURS

Fig.2. Flood of 27 January 1 9 4 8 : 1 = observed hydrograph at Callone; 2 = observed hydrograph at S. Giovanni alla Vena; 3 = computed hydrograph at S. Giovanni alla Vena.

178

APPLICATIONS

The method described has been applied to the study of flood routing in a reach of the river Arno (Gallati et al., 1971) between the hydrometers at Callone (54.03 km from the mouth of the river) and S. Giovanni alla Vena (35.81 km from the mouth), by means of which, in the events of flooding, hourly recordings are made of the levels of the river. The tributary Era flows into the reach at about 41 km from the river mouth. The river bed is described by cross-sections about 200 m one from the other surveyed by the Pisa section of the Italian Hydrography Office in the period 1953--55. In view of the considerable variability of the sections each one is defined geometrically by a set of points of its boundary. By means of the coordinates of these points it is possible to determine the various geometrical quantities which appear in eqs.7 and 8. The computing system was programmed on the IBM 3 6 0 / 6 7 computer at the National University Electronic Computing Center (CNUCE) at Pisa University. To test its validity the water levels at S. Giovanni alla Vena, calculated starting from the hourly known levels at Callone, were compared with the recorded ones.

co

:E

•

f

4"Y ,'7 "•

i

I

6

-

u+-4~l

i

I

12

'

i

n

-

.

n

,

.

.

I

.

u

2

.

~

I

I

~

18

TIM[

I

24

:

*

:

:

~-

~

6

J

,I

i

12

Hou,s

Fig.3. F l o o d o f 21 D e c e m b e r 1 9 6 0 : 1 = observed hydrograph at Callone; 2 = observed hydrograph at S. Giovanni alla Vena; 3 = c o m p u t e d hydrograph at S. Giovanni alia Vena.

179

r

~

=/

=

6

12

18

24

6

12

18

24 TIM[

6

12

18

24

6

12

18

24

HOUriS

Fig.4. F l o o d of 27 J a n u a r y 1 9 4 8 : 1 = o b s e r v e d h y d r o g r a p h at L e o n c i n i ; 2 = o b s e r v e d h y d r o g r a p h at S. G i o v a n n i alla V e n a ; 3 = c o m p u t e d h y d r o g r a p h at S. G i o v a n n i alla Vena.

o=

-

/

'

\,,,\/

.

o

........

I : ~ ',,, I, ~ ', : ~i~:~ : ~ i ;; ',;; ',;; I: :;:; ', : : ',: : ',: '. I ::I'.: 12 24 12 24 12 24 12 14 T IME

23

,.:

o

', : : ',: :', '.;I: :', : : ', : : ', ;; ', ; :;; ,~,, 12 24 12 14 12

Houe~

Fig.5. F l o o d o f 2 4 N o v e m b e r 1 9 4 9 : 1 = o b s e r v e d h y d r o g r a p h at L e o n c i n i ; 2 = o b s e r v e d h y d r o g r a p h at S. G i o v a n n i alia V e n a ; 3 = c o m p u t e d h y d r o g r a p h at S. G i o v a n n i alla Vena.

180 The results obtained for two flood events are given in Figs.2 and 3. The same m e t h o d was also used for a shorter reach of the river: from the h y d r o m e t e r at Leoncini {40.60 km from the river mouth), immediately downstream from the confluence with the Era, to the h y d r o m e t e r at S. Giovanni alla Vena. The results are shown in Figs.4 and 5. The value 0.037, which was deduced directly from experimental data (Gallati et al., 1971), was attributed to the Manning coefficient. COMPARISON WITH OTHER METHODS AND RESULTS In order to have an objective idea of the efficiency of the m e t h o d two other algorithms were used on the same implicit scheme, one prepared by R. Amaftiesi and F. Jonescu for a study on the river Danube (Amaftiesi and Jonescu, 1968), and the other by M. Amein and C.S. Fang for an application on the river Neuse (Amein and Fang, 1970). The two methods are briefly described below. From now on, for the sake of brevity, the three methods will be identified by the names of the rivers studied. Danube method

By assigning trial values to the water levels in the N sections (in the original method the trial values were assigned to the discharges), it is possible, for each grid element, to get the values of the discharge in the two sections that define the element from eqs.7 and 8. In this way to the generic section i are assigned two discharge values, resulting from the solution of eqs.7 and 8 for grid elements i -- 1 and i, respectively. For this section it is possible to define the quantity A Qi, i.e. the difference between the two discharge values thus calculated, which generally has a value other than zero, because the trial values of levels z i_,, z i and z i+, used for the calculation of the discharges do not generally coincide with the real ones. For obvious reasons of continuity, for each section it is necessary to set A Q; = O. The system of N -- 1 equations with N -- 1 unknowns is therefore obtained: AQ2

= f2(z2, z3) = 0

AQ3

= f3(z2, z3, za) = 0

(29)

AQN_,

= fN-,(ZN-~,ZN-,,ZN)=O

AQN

= f N ( Z ~ - , " ZN ) = 0

For the last section one of the two discharge values is obtained from the solution of grid element N -- 1 and the other from the rating curve.

181 System 29 is non-linear, and can be solved using the Newton method. The solution of the resulting linear system is very simple because, in view of the special form of the equations, the matrix of the coefficients that is obtained is tridiagonal, and can therefore be solved explicitly in accordance with the method of direct factorization already mentioned. The reason why the original m e t h o d of Amaftiesi and Jonescu has been modified is that with the above-described approach, for two adjacent sections that define a grid element, it is possible to write the discharges as a function of the levels in an explicit manner, as eqs.7 and 8 constitute an algebraic system of second order in the u n k n o w n s Qi and Qi÷," The two values that are obtained for each one of these have opposite signs and are about equal in absolute value; and therefore for physical reasons it remains certain that there is only one solution. It should be noted that during the iterations necessary for the solution of system 29, it may happen that for two adjacent sections the level of the downstream section is higher than that of the upstream section. When this situation arises the values calculated for the discharges relating to the grid element in question are not real. In order to avoid this, it is necessary to check continuously, and if necessary to correct the values of the water levels, increasing then the number of iterations and consequently the computing time.

Neuse method The algorithm used by M. Amein and C.S. Fang (1970) considers the variables z and Q at each vertex of the grid in Fig.1 as independent variables. In this way, intending usually all the variables related to time kA t and applying eqs.7 and 8, the following system of 2 ( N - 1) equations with 2N u n k n o w n quantities is obtained: Fl(Zl, Ql, z2, Q2) = 0 Gl(zl, Q1, z2, Q:) = 0

Fi(zi, Qi, zi+l, Qi+l) = 0 Gi(zi, Qi, Zi÷l, Qi+, )= 0

FN-, (ZN-,, QN-,' zN' QN ) = 0 GN-, (ZN-~, QN-~, zN, QN ) = 0

(30)

182

By adding the two boundary conditions:

zt(t) =5(t) I f(Zg, QN ) = 0

(31)

l

where 2{t) = known water levels at the upstream section, it becomes possible to solve the system using the Newton method. In view of the special form of the equations, the resulting coefficients matrix may be considered as block tridiagonal, where each block is a square matrix 2 × 2 (O'Loughlin and Short, 1971), and may therefore be easily solved by using the m e t h o d of direct factorization.

Adequacy o f results In order to get an objective standard of comparison between the various methods, for each flood event considered, the mean square error, the m a x i m u m positive and negative deviation, between the levels calculated and those measured, were computed. From this standpoint the three methods proved to be more or less equivalent, as, for example, may be seen from Table I, which contains the results of flood events simulated with the three methods; for each event the mean square error, the m a x i m u m positive deviation and the m a x i m u m negative deviation are shown. TABLE I C o m p a r i s o n a m o n g t h e t h r e e m e t h o d s : deviations o f the c o m p u t e d waves f r o m t h e measured ones Flood

Method

Mean square error (m)

Reach L e o n c i n i - - S . Giovanni 27 Jan. 1948 Arno 0.1138 Danube 0.1142 Neuse 0.1135 24 Nov. 1949

Arno Danube Neuse

Reach Callone--S. Giovanni 27 Jan. 1948 Arno Danube Neuse 21 Dec. 1960

Arno Danube Neuse

Max. deviations (m) -positive negative

0.2581 0.2581 0.2559

--0.2039 --0.2037 --0.2063

0.1467 0.1467 0.1468

0.3352 0.3351 0.3348

--0.2430 --0.2428 --0.2428

0.1059 0.1084 0.1073

0.0743 0.0958 0.0701

--0.1914 --0.1953 --0.1933

0.2747 0.2730 0.2741

0.4463 0.4386 0.4424

--0.3388 --0.3412 --0.3407

183 TABLE II Comparison among the three methods: number of iterations and computing time Flood

Method

Iterations

Reach Leoncini--S. Giovanni 27 Jan. 1948 Arno 1.73 Danube 2.37 Neuse 1.86

Time (sec • 10 -2) 0.77 1.24 0.79

3 Jan. 1949

Arno Danube Neuse

1.78 2.35 1.87

0.85 1.27 0.85

24 Nov. 1949

Arno Danube Neuse

1.77 2.34 1.84

0.81 1.26 0.81

2.19 3.27 2.35

1.13 1.69 0.90

Reach Callone--S. Giovanni 27 Jan. 1948 Arno Danube Neuse 17 Feb. 1960

Arno Danube Neuse

2.31 3.28 2.41

1.21 1.63 0.93

21 Dec. 1960

Arno Danube Neuse

2.23 3.06 2.26

1.18 1.54 0.88

Speed o f convergence The mean number of iterations necessary to calculate the condition of the river a t a given t i m e f o r a n u m b e r o f f l o o d e v e n t s is s h o w n in T a b l e II. F r o m t h i s t a b l e it m a y be seen t h a t t h e A r n o a n d N e u s e m e t h o d s r e q u i r e p r a c t i c a l l y t h e s a m e n u m b e r o f i t e r a t i o n s , a l t h o u g h t h e first m e t h o d is s l i g h t l y m o r e c o n v e r g e n t ; t h e D a n u b e m e t h o d on t h e o t h e r h a n d d e m a n d s a m a r k e d l y higher number of iterations. In a d d i t i o n t o t h e c a u s e p r e v i o u s l y seen in t h e d e s c r i p t i o n o f t h e D a n u b e m e t h o d , a n o t h e r r e a s o n f o r t h i s r e s u l t is t h e f a c t t h a t t h e f u n c t i o n s f o r w h i c h t h e z e r o s are s o u g h t are d i f f e r e n t in t h e t w o cases, since t h e S a i n t V e n a n t e q u a t i o n s a r e n o t u s e d d i r e c t l y in t h e D a n u b e m e t h o d . T h e s l i g h t d i f f e r e n c e b e t w e e n t h e A r n o a n d N e u s e m e t h o d s can m o r e o v e r be j u s t i f i e d b y t h e f a c t t h a t , u n l i k e t h e N e u s e m e t h o d , in w h i c h t h e t r i a l values o f t h e levels a n d o f t h e d i s c h a r g e s are i n d e p e n d e n t , in t h e A r n o m e t h o d t h e values o f t h e d i s c h a r g e are c o n g r u e n t w i t h t h e t r i a l values o f t h e levels.

184

Computing times All three methods lead back to the solution of a system of non-linear equations, specifically N -- 1 equations for the Arno and Danube methods and 2(N -- 1) equations for the other method. From Table II it may be seen that the Danube method, as was to be expected, requires longer computing times. As regards the other methods, although the speed of convergence is slightly higher in the Arno method, the calculation times are more or less the same for the shorter reach, while they become shorter with the Neuse method as the length of the reach increases, that is to say as the number of sections increases.

Storage req uiremen t Another factor to be considered is the a m o u n t of storage required by the various methods. In view, however, of the particular form of the equations from this aspect also, while the Arno and Danube methods need less storage, there are no appreciable differences as only a few vectors, and no matrices, have to be stored. CONCLUSIONS An implicit m e t h o d has been developed for the solution of the Saint Venant equations in a study relating to flood routing in open channels. The method is based on the observation that the u n k n o w n Q appears linearly in the mass equation which may therefore be easily integrated. As it is thus possible to express the discharge values as functions of the water levels, the only remaining unknowns in the problem are the values of the water levels along the entire reach studied. The resulting system of equations may be easily solved with the Newton method, in view of the particular form of the coefficients matrix. The particular advantages offered by the method have proved to be the rapid convergence and limited storage requirement. In an application relating to the river Arno, the method has also been proved efficient in the simulation of flood conditions. ACKNOWLEDGEMENTS This study forms part of a research conducted by the IBM Scientific Center in Pisa in connection with the setting up of a mathematical model of the river Arno. Specific thanks go to Prof. U. Maione, Head of the Institute of Hydraulics of Pavia University and also to Mario Gallati, Assistant Professor at the same Institute for their valuable advice and discussions.

185

REFERENCES Amaftiesi, R. and Jonescu, F., 1968. Un module math(imatique pour le calcul de la propagation des ondes de crue. Proc. Journ. Hydraul., Xes, Paris, Rep. 9-II. Amein, M. and Fang, C.S., 1970. Implicit flood routing in natural channels. J. Hydraul. Div. A.S.C.E., 96(HY12): 2481--2500. Cunge, J.A. and Wegner, M., 1964. Int(!gration num(irique des (!quations d'(icoulement de barr~ de Saint Venant par un sch(ima implicite de diffeirences finies. La Houille Blanche, 1 : 33--38. Evangelisti, G., 1966. Mathematical models of open channel flow. Proc. NATO Adv. Stud. Inst., Bressanone. Gallati, M., Greco, F., Malone, U., Panattoni, L. and Todini, E., 1971. Flood routing in open channels: a survey on the Arno river. IBM Scientific Center of Pisa, Tech. Rep. 17 : 40 pp. Isaacson, E. and Keller, H.B., 1966. Analysis of Numerical Methods. Wiley, New York, N.Y., 541 Pl O'Loughlin, E.M. and Short, D.L., 1971. Implicit flood routing in natural channels. J. Hydraul. Div., A.S.C.E., 97(HY10): 1771--1774. Strelkoff, T., 1970. Numerical solution of Saint-Venant equations. J. Hydraul. Div., A.S.C.E., 96(HY1): 223--252. Yevdjevich, V.M., 1964. Bibliography and Discussion of Flood Routing Methods and Unsteady Flow in Channels. U.S. Govt. Print. Off., Washington, D.C., 235 pp.