The matched diffusivity technique applied to kinematic cascades. II. Analysis of model performance

Journal of Hydrology, 121 (1990) 363-377

363

Elsevier Science Publishers B.V., Amsterdam

[31 THE M A T C H E D D I F F U S I V I T Y T E C H N I Q U E A P P L I E D TO KINEMATIC CASCADES. II. A N A L Y S I S OF MODEL PERFORMANCE

B.H. SCHMID

Institut fi~r Hydraulik, Gew~sserkunde und Wasserwirtschaft, Technische Universitiit Wien, Karlsplatz 13, A-1040 Vienna (Austria) (Received November 11, 1989; accepted after revision April 7, 1990)

ABSTRACT Schmid, B.H., 1990. The matched diffusivity technique applied to kinematic cascades. II. Analysis of model performance. J. Hydrol., 121: 363-377. Performance of the Muskingum-Cunge (MC) method is investigated by comparing numerical results with semianalytical solutions. Special attention is paid to cases involving discontinuous solutions. The adequacy of the MC approach is demonstrated for turbulent kinematic flows with and without the influence of reasonably strong shocks. The matched diffusivity technique also yields satisfactory results for laminar flows without shocks, whereas application of this method to shock-affected laminar flow cannot be advised. INTRODUCTION

Comparison of recorded and computed quantities (runoff, water depths) carried out in Part I of this paper showed Muskingum-Cunge (MC) routing to be the most suitable matched diffusivity method for application to kinematic cascades. As model validation could cover only a relatively small number of cases, a more detailed analysis of MC model performance is needed to complete the picture. This paper, therefore, deals with such an analysis, which was conducted with the aid of a semianalytical reference model to be described in the next section. Special attention has been given to cases characterized by the occurrence of shocks. As pointed out by Croley and Hunt (1981), errors produced by many currently available kinematic cascade catchment models may be because of the presence of discontinuities. Although shocks are mostly weak in this context, they can entail numerical problems with regard to stability and convergence. Errors from this source frequently remain undetected within the framework of a complex catchment model. Unfortunately, shock formation on kinematic cascades is the rule rather than the exception. Recalling, for instance, the criterion derived by Kibler and Woolhiser (1970, 1972), p+ _ ak ~'wh , a~" wk 0022-1694/90/$03.50

(1)

© 1990 Elsevier Science Publishers B.V.

364

B.H. SCHMID

with w d e n o t i n g the p l a n e width, a the k i n e m a t i c p a r a m e t e r (see P a r t I or eqn. (2) below) and index k the flow plane number. A c c o r d i n g to Kibler and W o o l h i s e r (1970, 1972), shocks are g e n e r a t e d for P~ > 1. Since wk and ak are i n d e p e n d e n t of wk_ 1 and a~ 1, it is u n l i k e l y t h a t all values of P~ for k = 2 . . n will be lower t h a n or equal to u n i t y for cascades comprising several (n) planes. One v a l u e of Ps(k - 1,k) > 1 is e n o u g h to s t a r t a first-order shock, in which case the model used must be able to cope with the discontinuity. It is for this r e a s o n t h a t the discussion of shocks figures p r o m i n e n t l y in this paper. Before the simulation, it should be recognized t h a t model input is no longer given as in P a r t I, but must be chosen in a realistic and m e a n i n g f u l way. To do this, some c o n s i d e r a t i o n must be given to the f r a m e w o r k w i t h i n which the model is to be operated, i.e. o v e r l a n d flow h y d r o l o g y and hydraulics. Reports of experiments usually c o n t a i n i n f o r m a t i o n on the flow regime ( t u r b u l e n t or laminar). For s y n t h e t i c simulations, w h e t h e r a n a l y t i c or n u m e r i c a l , a suitably based decision on what the flow regime will be like in the c i r c u m s t a n c e s defined must be taken. In this context, the Reynolds n u m b e r Re

v'h

-

-

q

(2)

where v is the m e a n velocity, h is the w a t e r depth, v is the k i n e m a t i c viscosity (here 1.15 × 10 6m2s 1) and q is the discharge per unit width, may serve as a guideline. R e p o r t e d values of the t r a n s i t i o n Reynolds n u m b e r s c a t t e r widely (mostly b e t w e e n 100 and 1000). Considering t h a t r a i n d r o p impact tends to i n t r o d u c e a significant a m o u n t of t u r b u l e n c e and t h a t the h i g h e r t r a n s i t i o n Reynolds n u m b e r s are usually valid for smooth surfaces, a choice of Retr = 150 is justified ( P e t r a s c h e c k , 1978). MC r o u t i n g as described in P a r t I uses a single-valued r e l a t i o n s h i p q = q ( h ) of the form q -- a . h m

(3)

P a r a m e t e r s a and m are different for t u r b u l e n t and l a m i n a r flows, respectively. If M a n n i n g ' s law is applied to the t u r b u l e n t case, the k i n e m a t i c p a r a m e t e r s are o b t a i n e d from 1 a

=

-"

x//So

n

(4)

where n is M a n n i n g ' s coefficient in SI units, so is the bed slope and m =

5/3

(5)

In the l a m i n a r case, use of the D a r c y - W e i s b a c h r e s i s t a n c e law finally leads to 8 "g'so

a

-

K'r

(6)

where g is 9 . 8 1 m s -2, K is a surface-related p a r a m e t e r and m =

3.0

(7)

THE MATCHEDD1FFUSIVITYTECHNIQUE

I1

365

A SEMIANALYTICALREFERENCE MODEL Before the description of this model, a short remark on its property of being semianalytical is called for. As will be seen from the equations given below, the prefix 'semi' indicates only the use of iteration techniques for the solution of certain nonlinear equations within an otherwise entirely analytical framework. As in other contexts the word 'semianalytical' has been used for numerical models containing some analytical components, it is pointed out that this is not the case here. The purpose of the model is the solution of the kinematic wave equation on two successive planes, assuming constant effective rainfall intensity and including the effect of shocks, if any. As, under the conditions defined, shocks will be relevant only in connection with the rising limb of the hydrograph (Croley and Hunt, 1981) and shocks constitute an important aspect of this analysis, it was decided to restrict the investigation to the equilibrium case. The kinematic wave equation to be solved for both planes should be recalled: ~q ?h + ~x ~t

--

i~

(8)

Substitution of eqn. (3) into eqn. (8) yields ~h --+

?h

a.m'h"

~t

1 _ _

~x

=

ie

(9)

and the corresponding system of characteristics is given by dh dt dx _

dt

ie

(10)

a.m.h

m 1

(11)

Initial and boundary conditions for the first (upper) plane are trivial (index denotes plane number): h,(x,,t-0)

= 0

(12)

-

(13)

and h , ( x l = O,t)

0

As shocks will not be encountered on the uppermost plane of a cascade with i~ constant (Croley and Hunt, 1981), the solution of eqns. (10) and (11) can be found immediately and yields the following result: l i e ' t f ° r t _< te, (14) L l ,t) hl(Xl ) lie tel for t _> tel with tel = ( L 1 / a l %"m 1,V,~ ~ as the time of equilibrium for Plane 1. Kinematic flow on Plane 2 is more complex. Whereas the initial condition is still the trivial one, h2(x2,t - 0) = 0, the boundary condition is not:

B.H. SCHMID

366

0,t) = ( wl

h2(x 2 =

' a l ~ 1/m

"hl(X 1

Ll,t )

(15)

\W2 "a2/ with w, and w2 the respective widths. Following the example of Kibler and W o o l h i s e r (1970, 1972), a shock p a r a m e t e r is defined in a c c o r d a n c e with eqn. (1): po

wl "al

(16)

w 2 "a 2

F u r t h e r m o r e , it was found c o n v e n i e n t to d e n o t e p~,m by "P~m". As explained above, Ps > 1 indicates shock f o r m a t i o n on the lower plane, and Ps ~< 1 allows shocks to be excluded from f u r t h e r analysis. Accordingly, s u b s e q u e n t model description will deal with these two cases separately.

(a) Ps <~ 1, case without the occurrence of shocks F i g u r e 1 yields a s u r v e y of the c h a r a c t e r i s t i c s in the (Xz,t) plane. It is a p p a r e n t t h a t up to t h r e e different solution domains must be t a k e n into a c c o u n t for water-level c o m p u t a t i o n at some a r b i t r a r i l y chosen time tw and location x2j. These domains are s e p a r a t e d from each o t h e r by the limiting c h a r a c t e r i s t i c and by a c h a r a c t e r i s t i c o r g i n a t i n g from the point (x2 = 0, t - te,), respectively. The values of Xe2 and xo,~2 can be o b t a i n e d from xe2(tw)

a2 .'rote 1 .{[tw _

=

tel "(1

Psm)]m - P~ 'C~}

(17)

and ,'m ] . m

x0.e2(tw)

a2 te

tw

(18)

F o r any x2j ~< xe2(G) the c o r r e s p o n d i n g w a t e r level is ll]m

[

(19)

(x2j,tw) = ,i~ "x2j + p~ %.m .tern1 L a2

t

, t~(L~) .....-~

. ~ . I'~''~ ~ w

7"'~-

t~,~(Lt)

- - - ' - ~ =-

•/ .11"~ _ t°, •/ / i'/l /A X2

h~(x~ = O, t)

0

xe2(tw)

xo~2(t~,)

L2

Fig. 1. Schematic survey of the characteristics related to Plane 2 (case without shock formation).

THE MATCHED D1FFUSIVITY TECHNIQUE

lI

367

whereas for xe2(tw) < x2j < x0,e2(tw) a n o n l i n e a r e q u a t i o n in t2,0 must be solved first (e.g. by means of the N e w t o n - R a p h s o n algorithm):

x2j :

a2

i~ " {[Psm .ie "t2~o + ie "(tw

t2.0)]m - P~ . (i~ .t,~0)m}

(20)

with (x2 = 0, t = t,~0) as the s t a r t i n g point of the c h a r a c t e r i s t i c considered. h2 (x2.j, tw) can t h e n be computed from

h,2(x2y,t~)

=

Psm

ie "Q0 + ie' (tw - Q0)

(21)

x2j >/ X0.e2,finally, is the trivial case of

h,2(x~j,tw)

= ie'tw

(22)

H y d r o g r a p h c o m p u t a t i o n at x2 = L2 follows the same lines, with x2j = L2 and the i n t e r v a l limits now c o n s t i t u t e d by

to.~,2(L~)

(LJa2 "i~ ')~''

(23)

and to~(L~)

=

a~ : io

1 + P~ " tom~

+ (1 -

P~o) " to~

(24)

(b) P~ > 1, case influenced by shocks on P l a n e 2 The domains described in subsection (a) r e m a i n formally u n c h a n g e d . The limiting c h a r a c t e r i s t i c is, however, replaced by the shock path, across which a d i s c o n t i n u i t y of the solution exists. A m e t h o d of shock p a t h c o m p u t a t i o n was published by Kibler and Woolhiser (1970, 1972), a modified version of which was i n c o r p o r a t e d into the model discussed here. Once the l o c a t i o n of the shock is known, the solution can proceed in the same way as described above. C O M P A R I S O N OF S E M I A N A L Y T I C A L AND M C - B A S E D R E S U L T S

F o u r different series of s i m u l a t i o n runs were carried out in the course of this investigation. Series A t r e a t s the case of t u r b u l e n t flow on a two-plane cascade of uniform slope 0.05 and rainfall i n t e n s i t y 1.5 mm min 1. The upper plane is c h a r a c t e r i z e d by a l e n g t h of 12 m, M a n n i n g s ' s n of 0.05 and v a r y i n g width Wl. The lower plane has the same length, slope and r o u g h n e s s as the u p p e r one, but the width is fixed at 6.0 m. F r o m eqn. (16) it is clear t h a t shocks will form only if wl > 6.0 m. F i g u r e 2a shows the h y d r o g r a p h s for t h r e e different values of w 1 ~< 6m, including wl = w2 = 6 m, which is e q u i v a l e n t to the case of a single plane of 24-m total length. C o r r e s p o n d i n g values of c a l c u l a t e d w a t e r depth are given in Fig. 2b. In all t h r e e cases, the space i n c r e m e n t Ax was chosen as 2.0 m and the time step At as 2 s. T h r o u g h o u t the rest of this paper, full lines d e n o t e s e m i a n a l y t i c a l results, and dashed lines indicate MC-based solutions.

368

B.H. SCHMID 4,0ul=bm

-"'.-

~2.

wl wl

-

i

I m 2 m

1,0-

0.0 (a)

I

i

I

i

I

I

51

111

151

211

251

311

T i m (s) 0.50-

wl = 5 m

6. 411-

wl

= 3m

w1

2 m

D.30E

D, 2 n -

0,10-

6, DO

(b)

I 5I

I lII

I

I

I

151

211

25I

Time

"l Ill

(s)

Fig. 2. Simulation runs of Series A, without shocks. Upper plane width, w~, is varied. (a) Hydrographs. (b) Water depths corresponding to cascade outflow.

Figures 2a and 2b show clearly t h a t s e m i a n a l y t i c a l and n u m e r i c a l results m a t c h v e r y well, which is confirmed by a graph of w a t e r depth vs. distance x (x = 0 at the u p p e r m o s t point of the cascade) displayed in Fig. 2c (wl - 3 m). Results derived obviously indicate t h a t the case of k i n e m a t i c flow on a two-plane cascade w i t h o u t the o c c u r r e n c e of shocks can easily and reliably be modelled on the basis of the MC method. This s t a t e m e n t will be s u b s t a n t i a t e d by o t h e r examples (Series B D).

THE MATCHED

[]

D[FFUSIVITY TECHNIQUE

369

0.40-

~1

!

5

o o~

-

_

"~ O.20.c:

b0 s

•

T

l

I

I

I

I

l

2

4

5

8

11

I

12 5.)

I

I

l

I

I

l

14

lb

1D

21

22

24

(c) Fig. 2. (cont.)(c) L o n g i t u d i n a l sections r e l a t e d to w~ - 3 m.

Series A c o m p r i s e s t h r e e m o r e s i m u l a t i o n runs, e a c h w i t h P~ > 1 a n d t h e r e f o r e affected by s h o c k p r o p a g a t i o n . R e s p e c t i v e widths wl w e r e c h o s e n as 9 m (Ps = 1.5), 12 m (P~ = 2.0) and 15m (Ps = 2.5). C o m p u t e d s h o c k p a t h s a r e s h o w n in Fig. 3a, a n d c o r r e s p o n d i n g h y d r o g r a p h s c a n be seen in Fig. 3b. L o n g i t u d i n a l sections for wl = 15m are g i v e n in Fig. 3c. S p a c e a n d t i m e i n c r e m e n t s were c h o s e n as follows: wl = 9.0m: Ax = 2 . 0 0 m a n d A t = 2 s for b o t h planes; wl = 12.0m: Ax = 2.00m a n d At = l s for the u p p e r plane, Ax = 0.25m and At = l s for the lower one; w~ = 15.0m: Ax a n d At i d e n t i c a l to t h o s e g i v e n for wl = 12 m. Again, t h e a g r e e m e n t b e t w e e n n u m e r i c a l a n d s e m i a n a l y t i c a l r e s u l t s is v e r y good. As c a n be seen from Fig. 3c, the MC s o l u t i o n shows m i n o r o s c i l l a t i o n s i m m e d i a t e l y d o w n s t r e a m of the shock. T h e s e are too small to influence the h y d r o g r a p h s (Fig. 3b), w i t h the e x c e p t i o n of the case w~ = 9m, w h e r e the c o a r s e s t grid h a s b e e n used. In t h a t example, the d a s h e d line displays a slight b u l g e a r o u n d t = 100 s, w h i c h b e c o m e s m o r e p r o n o u n c e d as Ps i n c r e a s e s (space a n d t i m e i n c r e m e n t s b e i n g equal). It can, h o w e v e r , be r e d u c e d a n d t h u s k e p t w i t h i n r e a s o n a b l e limits by the choice of a finer grid, as h a s b e e n done for Wl = 12 m a n d wl = 15 m. At first glance, this u n s t a b l e b e h a v i o u r is surprising, as t h e w e i g h t i n g p a r a m e t e r 0 did n o t exceed 0.5, as r e q u i r e d f r o m the v o n N e u m a n n a n a l y s i s (see P a r t I, appendix). H o w e v e r , o n l y l i n e a r s t a b i l i t y h a s b e e n proved, w h e r e a s a s h o c k is a n e s s e n t i a l l y n o n l i n e a r p h e n o m e n o n . M o r e o v e r , e v e n L a x ' s e q u i v a l e n c e t h e o r e m , which, c o n s c i o u s l y or not, forms

370

B.H.SCHMID

141112D1DO-

bB~-4m2._l

;=

1

2

I---.1--t~, 3

4

5

5

7

9

9

16

11

12

x2 (m)

(a) 7-

wl ,, L S l 5 ,,"

~1. = 12 m

~i

le'l " 9 an

II

I

I

I

I

I

I

51

111

151

211

I

IN

I

311

Timg (s) (b) Fig. 3. Shock-affected turbulent flow cases for Series A, (a) Shock paths. (b) Hydrographs.

the basis of most finite difference computations, need not apply here (Abbott, 1979). To put it differently, and perhaps more clearly, consistency and stability no longer guarantee convergence. As, however, the parasitic oscillations encountered vanish, if Ax and At are chosen small enough, this is a case of nonlinear instability rather than one of a stable but non-converging solution. Although the tendency towards unstable behaviour is not strongly pronounced in turbulent overland flow cases, an adequate choice of grid must take this

THE MATCHED DIFFUSIVITY TECHNIQUE

|.

lI

371

71-

3|l

s

|.bO-

B. 5 1 ~ D. 4 n -

D. 2 0 -

ill" OO

51 s

! 2

| 4

I b

I O

I 1O

I 12 x (m)

I 14

I 15

I 1.9

I 20

! 22

I 24

(c) Fig. 3. ( c o n t . ) ( c ) L o n g i t u d i n a l s e c t i o n s r e l a t e d to w 1 -

15 m.

property into account whenever shocks are expected. As the linear stability analysis, unfortunately, does not offer any guidelines except 0 ~< 0.5, the proper grid size has to be determined from numerical experiments. However, no problems of this kind are encountered in cases not affected by shocks. The simulation runs of Series B describe kinematic flow (turbulent regime) on two successive planes, the upper one measuring 30 m in length and the lower one 20m. Width amounts to 10.0m for both, and roughness is given by Manning's n = 0.4. The slope of the lower plane is kept at 0.10 in all cases except the last one discussed (0.05), whereas the upper plane slope is the input quantity varied ceteris paribus. Effective rainfall intensity amounts to 1.0 mm min Under these conditions the shock formation criterion, eqn. (16), can be simplified to yield sl > s2

(25)

where sl is the slope of the upper plane and s2 is that of the lower plane. Figure 4a, therefore, shows the hydrographs related to two cases without the occurrence of shocks. The longitudinal sections displayed in Fig. 4b belong to sl = 0.05. MC computation was carried out with At = 0.5min and Ax = 2 m. The plots clearly confirm the adequacy of the model. Proceeding from these results, sl was increased to study shock-affected kinematic flow on Plane 2. No change in the space and time increments was necessary to deal with the case of sl = 0.20. F u r t h e r increase in s~(sl = 0.30) and therefore in shock strength (Ps = 1.73 instead of 1.41), however, leads to a

372

B.H. SCHMID

9-

7-

sl = l, ll/

/ /

5-

/

~5-

~4§ "~" 3-

2~.I

,

,

,

I

,

[

5

II

15 Time (min)

21

25

31

(a) 1.8D-

211 m i n

-

I. b l -

~

1,4l1,28"~ 1. l l -

/

11 mi~

u

1,4D0,20m. i~lj

I

I

I

I

I

I

I

I

1

I

5

16

15

2D

25 x (m)

30

35

40

45

56

[b) Fig. 4. Simulation runs or Series B, cases without shocks. Upper plane slope, s~, is varied. (a) Hydrographs. (b) Longitudinal sections related to s~ = 0.05.

reappearance of the bulge (see Fig. 5), the underlying cause of which has already been discussed. To control this undesirable property of the MC algorithm, Ax and At were varied both jointly and separately. As expected, hx may still be chosen freely with respect to the upper plane, as no shock is encountered there. Performance of the MC algorithm for Plane 2, however, strongly depends on Ax2 as well as

THE MATCHED D1FFUS1VITY TECHNIQUE

373

II

9-

8?bU, 5-

~_'I-

21il

J 5

i 11

i 20

!

15

Tim

i 25

i ]l

(-,in)

Fig. 5. Series B, hydrograph related to upper plane slope s~ = 0.30 and increments At Ax - 2.0m.

0.5 min,

9-

87bin ,%

U, 5~o - 4 6

3_

2-

!

i 5

~ 10

i 15

TiN

i 20

J 25

i :]e

{min)

Fig. 6. Series B, hydrograph for s 1 - 0.30 obtained with At = 0.25min and Ax2 = 0.5m. At. I t w a s f o u n d t h a t a d e c r e a s e i n Ax2 a l o n e o r i n A t a l o n e d o e s n o t i m p r o v e the numerical outcome notably. Simultaneous choice of smaller space and time increments, however, could be shown to produce satisfactory results, as can be s e e n f r o m F i g . 6.

374

B.H. SCHMID

9-

7-

g4-

21-

D

I

[

I

I

I

I

5

1D

15

21

25

31

Ti~

(min)

(a)

2.2D2. IN-

21 ani.

t.,a2

t. 41. ~ 1. 211-

'.811:O. bO: ~

~

"L

0.402

~

5 enin

4 min

__

D.20O, DD

(b)

i

i

I

a

I

I

I

!

!

I

5

ID

15

20

25 x Cm)

30

35

40

45

50

Fig. 7. Series B, results for s z = 0.40 and s2 = 0.05. (a) Hydrograph. (b) Longitudinal sections. F o r Series B, it m a y be n o t e d t h a t s h o c k s t r e n g t h as reflected by the p a r a m e t e r Ps depends on the s q u a r e - r o o t of the slope (see eqns. (4) and (16), respectively) so t h a t the effect of a n y increase in sl will be damped considerably. The last case t r e a t e d w i t h i n Series B was aimed at g e n e r a t i n g a fairly s t r o n g shock, w h i c h was a c c o m p l i s h e d by the choice of sl = 0.40 and s2 = 0.05, yielding Ps = 2.83. Ax2 and At were kept at t h e i r previous values of 0.5 m and

THE MATCHED DIFFUSIVITY TECHNIQUE

II

375

TABLE 1 Laminar cases: survey of input data Input

Series C

Series D

Rainfall intensity (mmmin ]) Length of upper plane (m) Width of upper plane (m) Length of lower plane (m) Width of lower plane (m) Slope of upper plane Slope of lower plane Roughness parameter K (both planes)

1.0 5.0 3.0 5.0 3.0 Varied 0.10

1.0 6.0 3.0 4.0 3.0 Varied 0.10

100

1000

0.25m, respectively. Results are s h o w n in Figs. 7a and 7b, and once more d e m o n s t r a t e the validity of the MC c o n c e p t in the c o n t e x t of t u r b u l e n t k i n e m a t i c c a s c a d e flow c o m p u t a t i o n . The r e m a i n i n g s i m u l a t i o n runs, d e n o t e d as Series C and D, are dedicated to the a n a l y s i s of MC model p e r f o r m a n c e u n d e r l a m i n a r flow conditions. I n p u t d a t a h a v e been s u m m a r i z e d in Table 1. Two examples were selected to i l l u s t r a t e the modelling results for P~ ~< 1, one t a k e n from Series C, the o t h e r from Series D. F i g u r e s 8 and 9 show excellent a g r e e m e n t b e t w e e n s e m i a n a l y t i c a l and n u m e r i c a l solutions, the l a t t e r being o b t a i n e d with Ax2 = 0.50m and At = 2 s in both cases. F o r shock-affected l a m i n a r flow on a two-plane cascade, however, it soon became clear t h a t MC r o u t i n g does not p r o d u c e s a t i s f a c t o r y results even for 1.51-

1.4l-

sl = 1 . 1 . 0 / /

in %.

_

/ /

I 15

,,_ 1.3OLI o I:: -i

o~ ! . 2 1 -

O. 11-

i. li I

i 11

i 21

i 31

t 41

i 51

i bl

i 71 Tim

I I IXI I I (s)

i III

t t 211 I n

i s i I ] i i 141 151

Fig. 8. Simulated hydrographs of Series C, cases w i t h o u t shock formation.

376

B.H. SCHMID

O. bO-

0.50-

1.411-

51 =

l. l O #

/'-

]

f

~_B. 36o r -1

" g. 2D-

6. l O -

ll'lid

'

'

I

'

3D

'

I

'

'

bO

I

'

90

'

I

126

'

'

I

150 Time

'

(s)

'

I

100

'

'

I

21D

'

'

I

240

:

'

I

'

2TO

'

I

lid

Fig. 9. S i m u l a t e d h y d r o g r a p h s of Series D, c a s e s w i t h o u t s h o c k formation.

O. b O -

0,50-

~_ | . 3 J -

~: 6.2DL' O. l D -

6.0!

I 10

I 26

I 36

I 46

I 56

i 5D

J TO

I OD

Time (s)

, 90

J i , J i i 100 1 1 6 12D 130 146 156

Fig. 10. Series C, e x a m p l e of shock-affected h y d r o g r a p h for l a m i n a r flow case, s 1 = 0.18.

fairly w e a k shocks. The e x p o n e n t m is 3.0 h e r e instead of 5/3 for t u r b u l e n t flows so t h a t the i n h e r e n t n o n l i n e a r i t y is now m a r k e d l y s t r o n g e r t h a n before. A l t h o u g h a considerable n u m b e r of s i m u l a t i o n r u n s was carried out, the c o n c l u s i o n r e m a i n e d essentially the same, as must be d r a w n from Fig. 10 (Ax2 = 0.33m, At = 0.5s). Clearly, the case of shock-affected l a m i n a r flow c o n s t i t u t e s a limit to the applicability of MC routing.

THE MATCHEDDIFFUSIVITYTECHNIQUE

II

377

SUMMARY AND CONCLUSIONS

A large number of simulations was carried out to make a full comparison between semianalytical and MC-based computations of kinematic flow on a two-plane cascade. In the course of this investigation, MC routing was shown to produce perfectly satisfactory results for turbulent flow cases, both with and without reasonably strong shocks. It was found necessary, however, to decrease space and time increments as shock strength increases. Favourable results were also obtained for laminar flow without shocks. Shock-affected laminar kinematic flow, however~ proved to be outside the range of applicability of the MC model. REFERENCES Abbott, M.B., 1979. Computational Hydraulics. Pitman, London, 143. Croley II, T.E. and Hunt, B., 1981. Multiple-valued and non-convergent solutions in kinematic cascade models. J. Hydrol., 49:121 138. Kibler, D.F. and Woolhiser, D.A., 1970. The kinematic cascade as a hydrologic model. Colorado State University, Fort Collins, CO, Hydrol. Paper 39, 27 pp. Kibler, D.F. and Woolhiser, D.A., 1972. Mathematical properties of the kinematic cascade. J. Hydrol., 15:131 147. Petrascheck, A., 1978. Die Berechnung des Oberfl/ichenabflusses von Fl~ichenelementen (in German). 0sterr. Wasserwirtschaft, 30 (3/4): 65 72.