A nonlinear hydrologic cascade

Journal of Hydrology, 51 (1981) 283--293

283

Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands

A NONLINEAR HYDROLOGIC CASCADE

VIJAY P. SINGH* and SOMKID BUAPENG

School of Engineering and Applied Science, The George Washington University, Washington, DC 20052 (U.S.A.) Groundwater Division, Department of Mineral Resources, Bangkok (Thailand) (Accepted for publication October 10, 1980)

ABSTRACT Singh, V.P. and Buapeng, S., 1981. A nonlinear hydrologic cascade. In: L.R. Beard (GuestEditor), Water for Survival. J. Hydrol., 51: 283--293. A nonuniform, nonlinear, hydrologic cascade is developed for watershed runoff modeling, and its state--space representation is provided. By applying the cascade model to natural watersheds, an appraisal of its potential in runoff modeling is made. The model parameters are studied regarding their sensitivity. Some common features between the proposed model and several existing runoff models are indicated.

INTRODUCTION

The growing interest in the development, calibration, testing and application of nonlinear runoff models stems primarily from recognition of the nonlinearity of the rainfall--runoff process. The past two decades have witnessed a wide variety of nonlinear approaches to r u n o f f modeling. Of them all storage models appear to be most prominent, in a large measure for their simplicity. Models of Diskin (1964), Kulandaiswamy (1964), K.P. Singh (1964), Dooge (1967), Prasad (1967), Mein et al. (1974), Reed et al. (1975) and V.P. Singh {1976) are some of the examples of nonlinear r u n o f f models. Like all storage models, these models are based principally on: (1) the spatially lumped form of the continuity e q u a t i o n ; a n d (2) the storage--discharge relationship. The differences amongst the various storage models arise from the difference in consideration of this storage relationship. A striking feature of the storage models is that the watershed geometry is n o t considered explicitly. No consideration is given to the spatial variability of rainfall either. In this study our objective is to develop a nonlinear runoff model of storage type and test it on natural watersheds.

* Now with Department of Civil Engineering, Mississippi State University, Mississippi State, MS 39762, U.S.A.

284 NONLINEAR HYDROLOGIC CASCADE MODEL FORMULATION The nonlinear cascade represents a watershed by a series of unequal nonlinear storage elements or reservoirs as shown in Fig. 1. The governing equations for a nonlinear element are a spatially lumped continuity equation and a nonlinear storage--discharge relationship, which can be written respectively as: and

p = q+ds/dt

q = ks x

(1),(2)

where p is inflow to the storage element in cm/hr., q outflow from the storage element in cm/hr., s storage in the storage element in cm, t time in hours, d s / d t rate of change of storage in the element, k storage characteristic parameter, and x and index of nonlinearity. Eqs. 1--2 can be combined to yield a single differential equation relating inflow and storage: ds/dt

Pl (t)

kl'x sl(t)

(3)

p--ks ~

=

p2(t)

l

l k2, x S~(t)

p3(t)

Q2(t } -

k3, x s3(t)

pn(t )

qn(t ) kn,x Sn(t)

qn(t L

Fig. 1. A nonlinear hydrologic cascade with distributed input. The nonlinear cascade of Fig. I can then be represented by a system of equations: ds l / d t

= Pl --kls~

ds:/dt

= P2 + klSXl --k2S~

ds3/dt

= Ps + k : s ~ - k 3 s

dsi/dt

= Pi + kj-1 s~-i -- k j s f

dsn-1/dt ds./dt

~

(4)

P , - 1 + kn-2 Sn-2 x x -- kn-18n-1 = p . + k . _ l S ~ _ l -- k . s ~

where n denotes the number of elements in the cascade. Note t h a t the parameter x does n o t change from one element to another whereas k does. Eq. 4 can be written more conveniently and concisely in the matrix form as: =

P + KBS

(5)

285 where

u

dsl /dt

Pl

kl

Sl

ds2/dt

P2

k2

82

P3

k3

83

dsJdt

i,

_-

d s . -1 ~dr ds./dt

k

P~

=

;

kj

Pn -I

k . -1

Pn

k.

s1

8 . -

Sn

and S~ -1

0

0

...

0

0

84 -1

-- s~ -1

0

...

0

0

0

sX~-1

-- s~ -1

...

0

0

0

0

0

x-I 8n-1

- - 8n~-1

B=

O u r main p u r p o s e is t o o b t a i n the relationship b e t w e e n the final o u t f l o w q . and t h e lateral inflows P l , P2, P3, • • • , Ps, • • •, P . - 1 , P . . T h e storage--discharge relationship f o r the cascade can t h e n be expressed in the m a t r i x f o r m as:

(6)

Q = KCS where

Q=

m

ql

s~ -1

0

q2

0

s~ -~

qj

C =

0

0

...

0

0

0

0

x-1 Sn-1

0

0

s~ -1

q, a n d K a n d S are as d e f i n e d previously• Eqs. 5 a n d 6 p r o v i d e a state variable r e p r e s e n t a t i o n o f the cascade• F o r s i m u l a t i o n o f r u n o f f , h o w e v e r , the cascade with l u m p e d i n p u t was c o n s i d e r e d ; t h a t is, the inflows P2, P3, • • •, P , in P are zero a n d P l is positive.

286 Cascade o p e r a t i o n

The cascade operation can be summarized as follows: (1) Specify parameters h i , h2 . . . . . hn ; x and n. (2) Select a time interval A t = ti+s -- ti, i = 0, 1 . . . . (3) At the beginning of time, t = 0, B and S are zero. The inputs P2, P3, • • • , Pn in th e matrix form are all zero; only p 1 is positive. Given t he input p atter n (rainfall-excess) use eq. 5 to c o m p u t e S. So the value of S at t = 0 is obtained. This is the initial condition. (4) I n p u t into the j t h storage element at the beginning of a particular time interval At is an impulse of input, equal to the sum of input p~., the surface inflow qj._~, and the surface o u t f l o w of qj. The impulse input causes an instantaneous change in storage of the j t h element. Responses of the elements to the impulse inputs determine the out put s and states of the elements at the end o f the time interval At; the new states are the initial conditions for the following time increment. (5) F o r given initial states of the system, B, S and P at the beginning of a particular time interval At, S can be obtained from eq. 5. The state of the system, at the end of this time interval At, is the initial state at the beginning of the n e x t time interval At and is given by (S + S). Outflows from the n elements at the end of the time interval can be obtained from eq. 6. Derivation of other models

Several linear as well as nonlinear models can be deduced from this proposed cascade by restricting its parameters and i nput distribution. If h 1 = h2 = • . . = hn, the result is a u n i f o r m l y nonlinear model (Dooge, 1967; V.P. Singh, 1976). If x = 1, hi =/=k i + l , i = 1 , . . . , (n -- 1), the resulting model is a linear distributed parameter model. If x = 1, h~ = k2 = • . . = kn, t he resulting mo d el is a cascade o f linear reservoirs, know n popularly as t he Nash (1957) model. B ot h linear and non-linear models can be made bot h lumped i n pu t models and distributed i nput models (Dooge, 1959).

APPLICATION TO NATURAL WATERSHEDS For the mo del application thirty-eight natural agricultural watersheds were selected. These watersheds represent several regions of the U.S.A. Of t h e m sixteen watersheds are near Riesel (Waco), Texas; five near Hastings, Nebraska; five near Coshocton, Ohio; eight near Oxford, Mississippi; and one each in Watkinsville, Georgia, McCredie, Missouri; Fenni m ore, Wisconsin; and Ralston Creek, Iowa. For c o m p l e t e i nform at i on on watershed characteristics see Shelburne and Singh (1976) or U.S.D.A. publications entitled H y d r o l o g i c D a t a f o r E x p e r i m e n t a l A g r i c u l t u r a l W a t e r s h e d s in the U n i t e d States.

287

Determination o f mean areal rainfall Rainfall--runoff data of all 38 watersheds were obtained from the aforementioned U.S.D.A. publications. These publications are released almost every year and contain for each watershed the largest yearly event. Although there may be more than one raingage on a watershed, data are normally available in these publications for a centrally located raingage indicating that this represents the mean areal rainfall. For consistency this practice was followed for each watershed.

Determination o f rainfall-excess Rainfall-excess is inflow to the model and is determined by subtracting infiltration from rainfall. Philip's (1957) equation was used to estimate infiltration, which can be written as: f = a + ~fit -°'s

(7)

where f is infiltration rate in cm/hr., fi a parameter, dependent on soil characteristics and initial soil moisture content, a a parameter, dependent on soil characteristics, and t time in hours. The parameter a was determined from watershed soil characteristics. To account for antecedent soil moisture conditions, the parameter ~ was allowed to vary with each rainfall event and was determined by preserving the mass continuity of water.

Choice o f objective function The following objective function was used in this study: M

F = ~ [Qp0 ( / ) - - Qpe(J)]2 ~ min

(S)

i=1

where F is objective function, Qp0 (J) observed hydrograph peak in cm/hr. for the jth element, Qpe (J) estimated hydrograph peak in cm/hr, for the jth element and M the number of rainfall--runoff events in the optimization set. The choice of this objective function is based on the findings of V.P. Singh (1975a, b, c, 1 9 7 6 ) a n d Shelburne and Singh (1976). Because it only requires the hydrograph peak from each rainfall--runoff event, computationally it is more efficient. When F is divided by the number of events in the optimization set, the mean square error can be obtained which reflects on the average error occurring in t h e performing of optimization.

Determination o f computation time interval To study the effect of c o m p u t a t i o n time interval on model parameters and its predictive performance a watershed, 4-H, near Hastings, Nebraska, was selected. This watershed has eighteen events. These events were divided

288

I E X

X IX

X 2

X

Z Ld

X

1 x 3 0

20

40

60

80

100

ERROR FUNCTION (%)

Fig 2. T i m e interval vs. e r r o r f u n c t i o n o f p r e d i c t e d r u n o f f o n w a t e r s h e d 4-H, Hastings, Nebraska.

into two sets; one set, called the optimization set, consisted of ten events; the other set, called the prediction set, consisted of eight events. The number of elements n was restricted to three and x was fixed at 1.5 (V.P. Singh, 1976). Then parameters k~, k2 and k3 were optimized for the optimization set o f events by the modified Rosenbrock--Palmer algorithm (Rosenbrock, 1960; Palmer, 1969; Himmelblau, 1972) using the objective function of eq. 8. The optimization was performed for different time intervals of 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 and 4.5 min. It was found that the parameters k 1, k 2 and k 3 changed with the time interval and the value of F increased with the increase in At as shown in Fig. 2. To see the effect of time interval on model performance runoff hydrographs were predicted for the events in the prediction set, using optimized values of k~, k2 and k 3 corresponding to each computation time interval. The relative error in both hydrograph peak and its time was found to increase with increasing At. These results led to the choice o f 1 min. as c o m p u t a t i o n time interval. For more details see Buapeng and Singh (1977).

Parameter optimization The nonlinear cascade has (n + 2) parameters: n, x, k l , k 2 , . . . , kn. The parameter x reflects the degree of nonlinearity in surface runoff. Based on the findings of V.P. Singh (1976), x was fixed at 1.5. The parameter n, denoting the number of storage elements in the cascade, must be greater than one to obtain proper hydrograph shape and will depend on the topographic complexity. Again, based on the studies by V.P. Singh (1976), and $helburne and Singh (1976), n was fixed at 3. Thus parameters for optimization are k l, k2 and k 3 . As before, rainfall--runoff events of each of the 38 watersheds were divided into an optimization set and a prediction set. Parameters k 1, k 2 and k 3 were then optimized for each optimized set of events b y the modified Rosenbrock--Palmer algorithm (Rosenbrock, 1960; Palmer, 1969; Himmelblau, 1972) in accordance with the objective function of eq. 8.

289

Hydrographic prediction Hydrographs were predicted for the events in the prediction sets of several watersheds, utilizing the optimized parameter values. For a few sample watersheds hydrograph peak characteristics are given in Table I. It is clear that the predicted hydrograph peak characteristics are in reasonable agreement with observed peak characteristics. A comparison of observed and predicted hydrographs for two sample events is shown in Figs. 3 and 4. These figures indicate that the cascade can simulate the entire hydrograph well.

T 100C

10.00 T

" ~ 7-5c LIJ

~

500

-J

250

w~ i,-z ~ <~

750 h b0 z 5.00 rY

XI.Xo~i0 XX_XXiOo xxXO o

< Z

×~

ooo

2 50

o

×

zc~

_ _ ~ . ° o o p ,o o x 25

50

75

100

i

L

125

x

150

0.00

175

TIME ( m l n ) ~

Fig. 3. H y d r o g r a p h p r e d i c t i o n for rainfall e v e n t o f A p r i l 2 4 , 1 9 5 7 o n w a t e r s h e d Y-IO, Riesel (Waco), Texas (k I = 0.1045; k2 = 0.1874; ks = 0.2; n = 3). -= rainfall i n t e n s i t y ; - - - = i n f i l t r a t i o n rate; x x x = o b s e r v e d r u n o f f ; o o o = p r e d i c t e d r u n o f f .

10"001

7"50

~

-

>od-fl

.d~, ~ Z

x

X

o o

b_ h 0 50O Z

o

Oo ××o ~o

250

XxO 0

--

<

oo

xo oLd:o: . . . . Z . . . . . .o 2-E~-~-2 0

20

40

60

80

x _ 100 120 TIME, rain

_

:

00o

140

Fig. 4. H y d r o g r a p h p r e d i c t i o n for rainfall e v e n t o f M a y 1 3 , 1 9 5 7 o n w a t e r s h e d Y - 8 , R i e s e l ( W a c o ) , T e x a s ( k l = 0 . 1 1 2 7 ; k2 = 0 . 2 0 1 2 ; k s = 0 . 3 ; n = 3). - - ---- rainfall i n t e n s i t y ; - - - = i n f i l t r a t i o n rate; x x x = o b s e r v e d r u n o f f ; o o o = p r e d i c t e d r u n o f f .

290

TABLE

I

Comparison

of observed and predicted

hydrograph

peak and its time

Watershed name

Date

Observed hydrograph peak (cm/hr.)

Predicted hydrograph peak (cm/hr.)

Relative error (%)

Observed hydrograph peak time (min.)

Predicted hydrograph peak time (min.)

Relative error (%)

H a s t i n g s 2-H

8-11-39 5-20-49 6-12-65 6-12-65 6-29-65 8-7-46 8-7-42 6-16-50 9-7-42 9-5-46

2.8194 0.7036 8.8138 2.1565 2,0676 3,7592 2.5197 0.1735 3.5306 2.2885

0.9774 0.1351 7.9924 2.1880 1.5737 1,5174 1,15120 0.0040 1.8000 1.3593

65.33 80.79 9,32 1.46 23.89 59.64 39.99 97.69 49.02 40.60

10.0 3.0 9.0 13.0 9.0 7.0 10 2.0 15.0 9.0

10.0 20.0 24.0 15.0 9.0 10.0 23.0 43.0 12.0 12.0

0 566.0 166.0 23.0 0 42.0 - 130.0 2,050.0 20.0 33.0

H a s t i n g s 4-H

8-11-39 6-20-42 9-5-46 6-1-51 7-13-52 6-12-58 6-12-65 6-12-65

4.5212 5.8166 3.8862 6.7564 9.1948 0.9195 9.7028 6.1468

3.4254 5.0194 2.1212 6.3084 9.9094 0.5000 8.2274 4.9645

24.24 13.71 45.42 6.63 -- 7 . 7 7 45.62 14.21 19.23

5.0 8.0 12.0 123.0 21.0 14.0 19.0 7.0

6.0 11.0 8.0 131,0 19.0 11,0 15~0 12,0

20.0 37.0 33.0 6.0 9.0 21.0 21.0 71.0

Riesel (Waco) W-1

3-12-53 6-23-59 3-29-65

2.8702 4.8006 5.8750

3.1407 4.6877 7.9949

9.42 2.35 36.08

59.0 87.0 101.0

26,0 111.0 71,0

55.0 27.0 29.0

Y-8

5-13-57 3-23-65 6-18-61

5.6642 5.7122 0.1986

5.5050 6.9353 0.4147

2.81 21,41 - 108.80

19.0 61.0 100

28,0 73,0 72,0

47.0 19.0 28.0

Y-10

4-24-57 5-13-57 6-9-62 3-23-65

6.8580 4.8514 1.0008 6.9248

6.7709 4.4727 1.3655 7.8860

1.27 7,81 36,45 13.88

35.0 26.0 38.0 69,0

36.0 33.0 38.0 81,0

2.0 26.0 0 17.0

H a s t i n g s W-8

6-7-42 5-11-44 7-6-45 6-1-51 6-26-52 7-6-52 6-16-57

0.2718 0.4521 0.4267 1.1328 0.4491 0.2195 0.6579

0.1693 0.0819 0.2626 1.2051 0.4790 0.0720 0.9788

37.69 81.89 38,46 --6.38 --6.65 67.19 ---48.78

64.0 75.0 71.0 154,0 106.0 62.0 212.0

48.0 112.0 55.0 168.0 88.0 85.0 149.0

25.00 49.33 22.54 --9.09 16.98 -37.10 29.72

Coshocton 97

4-25-61 6-18-40 3-911-64 6-12-57 6-28/29-40 6-4-41 9-23-45 7-11-46 6-28-57

1.3919 0.2946 0.4729 0.6604 0.5182 0.9144 0.8204 0,5359 1.8390

1.6C09 0,3684 0.6060 0.5747 0.6941 0.6623 0.8292 0.3902 1.5310

--15.01 25.05 --28.12 12,97 --33.96 27.57 --1.07 27.19 16.75

130.0 122.0 394.0 150.0 150.0 156.0 120.0 212.0 416.0

108.0 57.0 393.0 80.0 113.0 62.0 139.0 100.0 379.0

16.92 53.28 0,25 46.67 24.67 60.26 --15.83 52.83 8.89

291

Discussion of prediction results Prediction errors in hydrograph peak and its time were high in several cases. Browsing through the rainfall--runoff record and hydrograph computation indicated that it was due to inaccurate determination of rainfall-excess and poor synchronization between rainfall and runoff observations. Poor hydrograph prediction was particularly evident for the events having very complex rainfall distribution indicating the model sensitivity to rainfallexcess. Accurate determination of rainfall-excess is a major problem in rainfall--runoff modeling. In this study Philip's (1957) equation was used to determine infiltration. There are two parameters in this equation. These parameters depend on soil characteristics and initial soil moisture content. The latter is n o t usually known or difficult to quantify. In most cases these parameters cannot be estimated accurately.

Parameter correlation An a t t e m p t was made to correlate k l , k 2 and k3 with physically measurable watershed characteristics. These characteristics included average watershed slope, slope of the main stream, watershed area, width of watershed, length of main stream, drainage density, shape factor and stream order. These watershed characteristics were considered to affect the runoff hydrograph (Sherman, 1932; Edson, 1951; Taylor and Schwarz, 1952; Strahler, 1957; Gray, 1961, 1970; Black and Cronn, 1975). Regrettably, the correlations between the parameters and watershed characteristics were n o t encouraging. This led to the belief that parameters depended on other factors also, including perhaps rainfall characteristics, antecedent soft moisture conditions and the interactions between them. CONCLUDING REMARKS

The following remarks can be made a b o u t this study: (1) The nonlinear hydrologic cascade can predict watershed runoff satisfactorily. (2) Computation time interval is important for model performanc e and can be taken as I min. (3) Reliable correlations could n o t be established between the model parameters and watershed characteristics. This needs further investigation. (4) The model parameters have no correlation between themselves. ACKNOWLEDGMENT

The work u p o n which this study is based was supported in parts by funds provided through the New Mexico Water Resources Research Institute by

292

U.S. Department of Interior, Office of Water Resources and Technology, as authorized under the Water Resources Research Act of 1964, Public Law 88-379 as amended, under project number 3109-206. REFERENCES Black, P.E. and Cronn, J.W., 1975. Hydrograph responses to watershed model size and similitude relationships. J. Hydrol., 26: 225--565. Buapeng, S. and Singh, V.P., 1977. Studies on rainfall--runoff modeling, 7. A nonlinear hydrologic cascade. N.M. Water Resour. Res. Inst., N.M. State Univ., Las Cruces, N.M., Water Resour. Res. Inst., Rep. No. 87, 69 pp. Diskin, M.H., 1964. A basic study of the linearity of the rainfall--runoff process in watersheds. Ph.D. Thesis, University of Illinois, Urbana. Ill. Dooge, J.C.I., 1959. A general theory of the unit hydrograph. J. Geophys. Res., 64: 241-256. Dooge, J.C.I., 1967. A new approach to nonlinear problems in surface water hydrology: hydrologic systems with uniform nonlinearity. Int. Assoc. Sci. Hydrol., Publ. 76: 409--413. Edson, C.G., 1951. Parameters for relating unit hydrographs to watershed characteristics. Trans. Am. Geophys. Union, 32: 591--596. Gray, D.M., 1961. Synthetic unit hydrographs for small watersheds. J. Hydraul. Div., Proc. Am. Soc. Cir. Eng., 87(HY4}: 33--53. Gray, D.M., 1970. Handbook on the Principles of Hydrology. Water Information Center, Ottawa, Ont., pp. 5.1--5.55. Himmelblau, D.M., 1972. Applied Nonlinear Programming. McGraw-Hill, New York, N.Y., pp. 158--167. Kulandaiswamy, V.C., 1964. A basic study of the rainfall excess--surface runoff relationship in basin system. Ph.D. Thesis, University of Illinois, Urbana, Ill. Mein, R.G., Laurenson, E.M. and McMahon, T.A., 1974. Simple nonlinear model for flood estimation. J. Hydraul. Div., Proc. Am. Soc. Civ. Eng., 1 0 0 ( H Y l l ) : 1507--1518. Nash, J.E., 1957. The form of the instantaneous unit hydrograph. Int. Assoc. Sci. Hydrol., Publ. No. 45, 3: 1114--1121. Palmer, J.R., 1969. An improved procedure for orthogonalizing the search vectors in Rosenbrock's and Swann's direct search optimization methods. Comput. J., 12: 69-71. Philip, J.R., 1957. The theory of infiltration, 1. The infiltration equation and its solution. Soil Sci., 83: 345--357. Prasad, R., 1967. A no.nlinear hydrologic system response model. J. Hydraul. Div., Proc. Am. Soc. Civ. Eng., 93(HY4): 201--221. Reed, D.W., Johnson, P. and Firth, J.M., 1975. A nonlinear rainfall--runoff model, providing for variable lag time. J. Hydrol., 25: 295--305. Rosenbrock, H.H., 1960. An automatic method for finding the greatest or least value of a function. Comput. J., 3(3): 175--184. Shelburne, K.L. and Singh, V.P., 1976. Studies on rainfall--runoff modeling, 4. Estimation o f parameters of two mathematical models of surface runoff. N.M. Water Resour. Res. Inst., N.M. State Univ., Las Cruces, N.M., Water Resour. Res. Inst., Rep. No. 76, 98 pp. Sherman, L.K., 1932. The relation of runoff to size and character of drainage basins. Trans. Am. Geophys. Union, 13: 332--339. Singh, K.P., 1964. Nonlinear instantaneous unit-hydrograph theory. J. Hydraul. Div., Proc. Am. Soc. Civ. Eng., 90(HY2): 313--347. Singh, V.P., 1975a. A laboratory investigation of surface runoff. J. Hydrol., 25(2): 187-200.

293 Singh, V.P., 1975b. Hybrid formulation of kinematic wave models of watershed runoff. J. Hydrol., 27: 33--40. Singh, V.P., 1975c. Estimation and optimization of kinematic wave parameters. Water Resour. Bull., 11(6): 1091--1102. Singh, V.P., 1976. Studies on rainfall--runoff modeling, 5.A uniformly nonlinear hydrologic cascade model. 78, N.M. Water Resour. Res. Inst., N.M. State Univ., Las Cruces, N.M., Water Resour. Res. Inst., Rep. No. 78, 47 pp. Strahler, A.N., 1957. Quantitative analysis of watershed geomorphology. Trans. Am. Geophys. Union, 38: 913--920. Taylor, A.B. and Schwarz, H.E., 1952. Unit hydrograph lag and peak flow related to basin characteristics. Trans. Am. Geophys. Union, 33: 235--246.