Multivariate hydrologic time series analysis

Multivariate hydrologic time series analysis C. T. H S U and K. A D A M O W S K I Department of Civil Engineering, University of Ottawa, Ottawa, Ontario, K IN 9B4, Canada The objective of this study is to establish a multivariate watershed hydrologic system model involving meteorological data as the input and river flow as the output of the system. Monthly hydrological time series of runoff, temperature and precipitation were selected for analysis. A firstorder autoregressive-moving average (ARMA) transfer function model was found adequate to describe the multivariate watershed hydrologic system for the monthly runoff and meteorological time series. The results also indicated that the coordinated use of the meteorological and hydrometric data would enhance the accuracy of estimation of runoff characteristics.

INTRODUCTION In the past two decades, "joint mapping '~ or coordinated use of the meteorological and hydrometric data, from the simplest coordinated tracing of precipitation and runoff isolines 2 to the most sophisticated deterministic watershed model 3 have been proved to greatly enhance the accuracy of estimation of river flow characteristics. This branch of deterministic (or parametric) response hydrology is in full progress at present 3"4. Its counterpart - stochastic system modelling, however, still remains in the state of investigating isolated hydrologic phenomenon, in most cases, precipitation or river flow series. Therefore, it is felt that in order to explore the physical implication of a stochastic system model and to improve its efficiency in describing river flow characteristics, it would be desirable to apply the same principle of joint mapping, i.e., coordinated use of meteorological and hydrometric data to stoclaastic modelling of a watershed hydrologic system. In modelling the stochastic behaviour of a multivariate watershed hydrologic system, meteorological and runoff series are considered as the input and the output, respectively, of the system under study. Temperature and precipitation are the only meteorological series under consideration because they are the most commonly available meteorological data in Canada. It is reasonable to assume that the effect of other hydrologic variables such as soil moisture content and evaporation are also taken care of by the dynamic nature of the system, i.e., they depend on the current and previous states of temperature and precipitation conditions. A linear system can be described in the time domain by the impulse response function and in the frequency domain by the frequency response function. These system response functions for a multivariate system can be determined through multiple cross-correlation and multiple cross-spectral analyses. Multiple cross-spectral analysis bears analogy to the conventional multiple regression analysis except that for the former regression goes on separately at each frequency and that regression coefficients take complex values rather than real values, which enables one to learn more about the underlying relationship among various time series. The other advantage of the analysis in the frequency domain is that the 0 3 0 9 1708/81~020083 1352.00

© 1981 CML Publications

smoothing in the computation of spectrums reduces excessive sampling variations. Finally, the system is represented in the form of transfer function models. System transfer functions are constructed from the impulse response functions. The theory and approach adopted in this study are mainly based on the works by Jenkins and Watts 5 and Box and Jenkins 6. THEORETICAL BACKGROUND A time-variant multivariate linear system may be written in the matrix form: :x

(xq + 1(0) = f (h~ + l.~(u))(xj~t - u))du + (:o + l(t)) 0

(i= 1, 2 . . . . . r:

j = 1, 2 . . . . . q)

(1)

where (xq + 1(0) is the column vector of'r' output variables, (xj4t)) is a column vector of'q" input variables, (zq+ 1(0) is a column vector of'r" residual variables and (hq+,-j4u)) is a r x q impulse response matrix. In the case of the multivariate system analysis of runoff, temperature and precipitation series, equation (1) then becomes: 2t

X

x3(t)=f h31(u)xi(t-u)du+f h32(u)x2(t-u)du+z3(t) o o

(2)

It can be shown 5 that the minimum mse (mean square error) estimates of the impulse response functions are obtained from the following Wiener-Hopf equations:

t;~'~k(u-v)~hik(v))dv

(7~i(u))o

1i=3:

j, k = l , 2)

Adv. Water Resources, 1981, Volume 4, June

(3)

83

Muhivariate hydrologic time series analysis: C. T. Hsu and K. Adamowski where ;'j~lu)) is the lagged cross-variance between the output X~(t) and the input X)(t-u), 7~k(U) is the lagged covariance between the inputs X~It-u) and Xk(tt and hit(v) is the transpose of a matrix. By the Fourier transformations of equation (3), the estimation equation for the frequency response functions may be written:

S s,(D = (S )k(I) X H ,k(D ) (i = 3:

j, k = 1, 2)

With given impulse response functions, the transfer functions may be obtained from the relationship:

h3j~B)- d(B)

and vice versa. For systems with superimposed noise, equation (7) may be modified as: X3.n = E

K2., 20'))

j= I r=O

h3j.rxj.n-r+ Y3.n

(7a)

and the transfer function model, equation (8) then becomes 2

d(B)x3. . = ~ wj~B)xj., + z,

(8c)

j=l

where Y3,, and Z3. . are residuals in the form of either white noise or some other known processes.

{51 N U M E R I C A L ANALYSIS

where K2.12(f)=(H31(f}S31{jC)+ H32(fJS32~cj)/S33(f)

(6)

is called the squared multiple coherency spectrum. The multiple coherency spectrum measures the proportion of the output spectrum which can be predicted from the inputs. Hence equation 16) is the frequency domain analogue of the multiple correlation coefficient in the time domain. In the time domain, discrete-time modelling is more suitable for many design and planning purposes and a discrete-time approximation of the system model (eqn. (2)) can be obtained in the form:

2± X3.n = E

j=lr=O

h3j.rXj,n-r

(7)

It is, however, more desirable to design an adequate but 'parsimonious' model, i.e., to employ the smallest possible number of parameters yet adequately still representing the real system. This could be accomplished through an autoregressive-moving average (ARMA) type transfer function model: 2

d(B)x,= ~ wj~B)xs.,

(8a)

j=l

or more explicitly

( 1 - d, B -

. - dqBq)xi.,, =

~. (Ws.o - ws. lB -- ... w~.,B')xs. .

xi(t) -qi(t)-mi(t) S~(t)

i = 1, 2, 3

(10)

where q~(t) are the original observed values of the time series: mi(t) and S~(t) are the long-term monthly means and standard deviations of the time series and x~(t) are transformed stationary components with zero means and unit variances; the subscripts '1', '2' and '3' stand for temperature, precipitation and runoff respectively. The auto- and cross-covariance functions were estimated by 1 N-k

i,j=1,2,3

(ll)

(8b)

j=l

where d's and w's are constants; 'q" and "r' are the orders of the ARMA transfer function and are usually relatively small compared to the terms required for the impulse response functions.

Adv. Water Resources, 1981, Volume 4, June

The multivariate hydrologic time series analysis was carried out for eight natural regime watersheds in western Canada. Four stations are in the Pacific region: three in the Fraser River Basin (08KA007, 08KA005 and 08KB001) and one lake outlet station (08LA008). The other four are from Prairie stations: two in the mountainous Prairies (05AB021 and 05AB002) and two in the central and eastern Prairies (05GG001 and 05TD001). Both monthly runoff data and monthly meteorological data, i.e., temperature and precipitation, for the period 1952-1966 were obtained from the hydrological data bank resulting from the hydrometric network planning study for western and northern Canada by the Department of the Environment v. Second-order stationarity for all the time series were achieved by removing the seasonal components in monthly means and variances from the original time series by:

co~k) = N .El=Xi.nX. +k

2

84

(9)

(4)

where So.U')is the cross-spectrum between the output Xj4t) and the input XM), SSk(]) is the cross-spectrum between the inputs Xj~t) and Xk(t); and Hik(f) is the frequency response function corresponding to the impulse response function h~k(U). In addition to estimating the frequency response functions, it is necessary to characterize the properties of the residual variable Z3(t}. This is done by determining the residual spectrum $33 in the following form: $ 3 3 ( f ) = $33{if"~ 1 --

j = 1, 2

The auto- and cross-correlation functions were then obtained by

Cij(k) r°Ik)=/C,,(O)Cj,I O ) x /

i , j = l, 2, 3

{12)

Multivariate hydroloaic time series analysis." C. T. Hsu and K. Adamowski It has been shown 5 that variances of sample auto- and cross-correlations may both be approximated by:

1

Var(ro(k))~ ~

(13)

for a white noise or for two uncorrelated discrete processes with either or both being white noises. The sample estimate of autospectrum was obtained by

2,(f):z(1-a'~2,}X/v-2q / !

where v is defined by equation (16) and q = 1.. The above statistics were also computed between input series, i.e., temperature and precipitation, denoted by subscript '12'. The frequency response function of the system, equation (2) was obtained by solving equation (4):

M Sii= cii(O)+ 2 ~ w(k)c.(k)cos 2rcfk

/'~3 100 =

k=l

i=l, 2, 2 (1½~
(14)

where w(k) in the Parson lag window and M = N/4 ~ 48 (months). The 100~1" a)~o confidence interval for Silo0 is

v~,,00

~
x~(a/2)

x~(1 - (a/2))

.6)

(17)

where

k=1

To obtain expressions for the gains and phases, it is necessary to take the modulus and argument of equation (24). Then

and

~310q = tan - l -/~31 "431

where

k=l

S11S22 - [$1212

~31(f)

- 0j300

Cj3(f~

(19)

-2 -2 Cj300 + 0j300

K j3~O= S'=jfljqS--~3(~)

(20)

(21)

By taking analogy to the multiple correlation coefficient, z-transformation was carried out for the cross coherency spectrum s, 200 =½1n l +IKI

1 -tKi

Q13S22-Q23C12 - C23Q12 s11522 --i~x 212

(26b)

The gain and phase for Ba200 may also be obtained from equations (25) and (26) by exchanging '1' and '2'. Approximate confidence limits for gain and phase spectrum are given by: (1 - K 2

(27)

and

~(j)+sin- l it~2qF2q., - { 1-]~2

The squared coherency was computed by: -2

(26a)

(18)

The phase spectrum was then obtained by: .00 = t a n - '

(25b)

a31(J)= C13S33 +023012-C23C12

m

M

Qj300 = ~ w(kXc~a(k)-c3j&))sin2xfk

(24b)

(25a)

(-Sj3(f) = Cj300 -- iQj3(f)

M Cj3(f')=cj3 (0)+ E w(k}{cj3(k)+c3flk))cos2xfk

(24a)

/.'~32{f} = $2 3(J')S 1 l(f) - S 13(f)$2 ~00 S1100S22(f ) - S1200

The smoothed cross-spectrum may be expressed by:

}=1, 2 (-1/2~f<~1/2)

200 -- $2 300S1200 )-,~11(f)L~22(.f) - S1200 2

S1300S2

(15)

where x,,(1-(a/2)) and x,,(a/2) are Z,2, values at the 100(1 -(a/2))% and 100(a/2)~o levels respectively with the degree of freedom for Parson's lag window 5

)

(23)

(22)

and the approximate 100(1 -a)~o confidence interval for 200 is:

)

(28)

where F2q.v_ 2q is an F-statistic with degrees of freedom (2q, v - 2q) and q = 25. The impulse response functions in the time domain were then obtained by taking the inverse Fourier transformation of the corresponding frequency response functions. Therefore 1/2

1/2

h3j.,=2 f ,43fifJcos 2nfrdf+2 f /~3flJ)sin 2hi?dr 0 o

.j=l, 2

(29)

The multiple squared coherency between the runoff and the meteorological series was obtained from equation (6) K23.12(f}: $22[$3112+ SI 11S3212- 2Re(S12S23S31) S33(511 $22 - [ S l 212) (30)

Adv. Water Resources, 1981, Volume 4, June

85

Muhi~'ariate hydrolo,qic time series analysis: C. T. Hsu and K. Adamowski where

Re(S12523S31)~~- (~l 2 ~'23C" 13 "[- ~'12Q2.,,Q

13

--

(31)

Q,2Q2_,c,3 +Q_,2Q,3c23 The residual spectrum was computed using equation (5) (32)

~ 3 3 ( / ) = ~33(/'}[ 1 - - / ~ 2 . 1 2 [ f ) )

As in the multiple regression analysis, it is useful to measure the partial coherency spectrum between the output and one of the input processes after allowance is made for the effect of the other input processes. It can be shown 5 that the squared partial coherency spectrum can be obtained from the relationship:

1

-2 1 - / ~ 2 a 12(J) _K13.2ff ) - I - K 2z3[y) 1 --/~2 12(1)

1- / ( ~ 3 . , ( J ) - ] = ~

(1 - d l B - d 2 B 2 ) = 2

Z (Wj.o-Wj. I B - w j 2 B ' 2 ) x i . - h + z ,

(36)

j=l

where z,, is a residual term to be modelled by an ARMA (2,0) process: (1-01B-02B2)z,=a,

(33b)

The determination of the model parameters was carried out in the same manner as for the previous case except that the parameter optimization was performed on the following expression:

1,2

j = 1, 2

Second-order A R M A transjer Junction The second-order ARMA transfer function is defined by:

(33ai

The approximate variance for system impulse response functions has been shown 8 to be:

Varlh3j'r)=N1 f Sjj(I"~3 -K22)dJ

The above function minimization was carried out using Rosenbrock's algorithm of rotating coordinates ~1'~2 which allows the principal axis of optimum seeking always in the line of the steepest descendent and hence expedites the rate of convergence to some required accuracy.

(34)

(36a)

E ( x 3 . n - d l *x 3 . n - 1 -d'~x3.n_ 2 2 Z (,,'*0 - , , 7 , 8 j=l

- ,,.Z, 82).

,,._ h; -

-12

( 0'~ - O*B)z, _ 1)2--+min The approximate confidence limits for multiple and partial coherency spectrums were obtained in the same manner using equations (22) and (23) except that in this case q = 2 in equation (23). It has already been shown 9~° that most monthly runoff series may be successfully modelled by a first-order autoregressive process or AR MA(1, 0) process. Therefore only two ARMA type transfer function schemes of firstorder and second-order respectively were attempted for the identification of the underlying multivariate hydrologic system between monthly runoff and monthly meteorological time series, i.e.. temperature and precipitation. First-order A R M A tran,~ler Junction The first-order ARMA transfer function is defined by: 2

1 1 - d l B ) x 3 . , = ~ (Wi.o--%.lB)xi.,_h,+ z,,

(35)

j=l

where b i is an integer denoting lag time in months, and z, is a residual term to be represented by an ARMA[1, 0) process ( 1 - ~o a B).-.,,

= a,

(35a)

The initial estimates of the model parameters were obtained using the relationship (9) and the criterion given by equation (34). The final values of the parameters were those which minimized the variance of the noise A,, i.e.,

The efficiency of the models was then evaluated by the following procedures: (i) Comparison between the estimated impulse response function from equation (29) and the theoretical impulse response function using equation (9). Approximate standard deviations of the respective estimated system impulse response functions were computed by taking square root of the variances defined by equation (34). (ii) Comparison between estimated and theoretical frequency response functions, i.e., estimated and theoretical gain and phase spectra. Confidence limits of gain and phase spectrums were computed using equations (27) and (28) respectively. (iii) Examination of autocorrelations of noise series A, and its cross-correlations with input series X 5,, and X2. .. One basic assumption concerning the noise term of the above schemes is that it is of the nature of a purely random process, or white noises and it is uncorrelated with the input series. (iv) Examination of the unexplained portion of the variance or the variance of noise series A, due to the proposed models. Along with the presentation and discussion of the results in the following section, the above procedures are illustrated in graphical forms by a typical example at 05AB002 - - Willow Creek near Nolan, Alberta. RESULTS AND DISCUSSION

2

E(x3,. - d*x3.n- 1 - Z t w T . o j=l

_,,_,

--.rain

86

Adv. Water Resources, 1981, Volume 4, June

The mean and variance of the raw time series were computed for runoff, temperature and precipitation re-

Multivariate hydrologic time series analysis: C. T. Hsu and K. Adamowski The seasonal components ( F i g . 1) w e r e t h e n r e m o v e d f r o m t h e o r i g i n a l t i m e s e r i e s u s i n g e q u a t i o n i10). T h e resulting stochastic component would yield zero mean and unit variance throughout the study period and between months. Sample estimates of the autocorrelations for the transformed runoff, temperature and precipitation series were c o m p u t e d u s i n g e q u a t i o n s (1 1) a n d ( 121 a n d a r e i l l u s t r a t e d in F i g . 2. A u t o c o r r e l a t i o n s u p t o t h e first t w o l a g s a r e

s p e c t i v e l y f o r e a c h w a t e r s h e d a n d l i s t e d in T a b l e 1. Harmonic analysis was carried out for runoff, temperat u r e a n d p r e c i p i t a t i o n s e r i e s a n d is p r e s e n t e d in T a b l e 1. It can be observed that the annual cycle and its subharmonics, i.e., 12~ 6, 4~ 3, 2 . 4 a n d 2 ( m o n t h s ) h a v e a c c o u n t e d f o r 92.0~'o t o 97.3~o o f t h e t o t a l v a r i a n c e o f t h e t e m p e r a t u r e ; w h e r e a s o n l y 25.7~°~o t o 66.2~o f o r t h e p r e c i p i t a t i o n s e r i e s a n d in t h i s c a s e 33.5~o t o 9 3 . 2 0 0 f o r t h e runoff series.

Table 1. Summary t?fyeneral statistics of h)'droloyic time seriex inrohed in multit'ariate water.shed system ident(/ication 05AB002

Station identification

05AB021

900 12

Drainage area [sq. milesl Met. Stations available

446 3 Hydrological parameters

Statistical parameters Mean Variance 12 m o n t h s 6 4 3 2.4 2 Total

Explained seasonal variance ~V,,,)

Residual variance Auto and cross-corrc coeffs

Runoff

Temperature

Precipitation

Zero order 1st 2nd Zero order I st 2nd Zero order I st 2nd

Runoff

Temperature

Precipitation

Runoff

Temperature

0.20 in 0.14 in 2

38.6 F 255.2 ( F ) 2

1.65 in 1.60 in 2

0.30 in 0.29 in-'

37.3 F 229.5 ( F)-"

1.92 in 2.10 in-"

25.68 6.31 1.03 0.30 0.11 0.06 33.48

91.15 0.14 0.06 0.11 0.31 0.24 92.01

27.44 3.37 1.16 1.19 2.92 1.34 37.42

27.16 8.11 1.34 0.25 0.09 0.03 36.98

90.97 0.34 0.31 0.17 0.36 0.39 92.53

29.93 2.57 0.22 0.14 3.69 0.78 37.34

0.09 in 2

20.41 F):

1.00 in 2

0.18 in-'

17.1 I F) 2

1.31 in"

1.00 0.51 0.29 -0.06 0.09 0.08 0.14 - 0.15 - 0.17

-0.06 - 0.32 - 0.19 1.00 0.05 0.03 - 0.55 0.05 0.08

0.14 0.32 0.20 -0.55 0.01 - 0.05 1.00 - 0.01 0.03

1.00 0.57 0.39 +0.01 0.06 0.05 0.12 - 0.14 - 0.10

-0.10 - 0.26 - 0.23 1.00 0.03 0.04 - 0.50 0.04 0.10

0.12 0.27 0.22 -0.50 0.02 - 0.0 I 1.00 - 0.01 - 0.07

Station identification Drainage area Met. Stations available

05GG001

05TD001

46100 182

6030 7

Precipitation

Hydrological parameters Statistical parameters

Runoff

Temperature

Precipitation

Runoff

Temperature

Precipitation

Mean Variance

0.22 in 0.04 in 2

35.8 F 413.0 ( F ) 2

1.38 in 1.11 in"

0.45 in 0.11 in-'

30.5 F 620.5 I F)-'

1.43 in 0.82 in-'

66.75 4.26 0.14 1.02 0.66 0.16 73.01

94.90 0.28 0.04 0~02 0.11 0.15 95.68

51.37 13.24 0.03 0.46 0.97 0.10 66.17

31.18 6.28 1.62 0.87 0.30 0.06 40.31

96.44 0.46 0.09 0.05 ().0X 0.17 97.29

45.69 3.59 0.35 0.83 1.27 2.25 53.98

0.01 in-'

17.9 ( F):

0.38 in 2

0.07 in 2

16.8 ( F)-'

0.38 in-'

1.00 0.64 0.42 0.07 0.04 - 0.01 0.01 - 0.11 0.04

0.07 - 0.20 - 0.24 1.00 0.09 0.06 - 0.44 - 0.12 - 0.14

0.01 0.30 0.25 -0A4 -0.13 0.03 1.00 0.06 - 0.05

1.00 0.79 0.69 -0.10 -0.07 - 0.02 0.13 0.11 0.07

- 0.10 - 0. I g - 0.13 1.00 0. tl - 0.01 - 0.16 - 0.09 - 0.09

0.16 0.09 -0.16 -0.12 - 0.01 1.00 0.01 - 0.03

12 months 6 4 3 2.4 2 Total

Explained seasonal variance

Residual variance Auto and cross-corr. coeffs.

Run off"

Temperature

Precipitat ion

Zero order I st 2nd Zero order 1st 2nd Zero order I st 2nd

0.13

Adv. Water Resources, 1981, Volume 4, J u n e

87

Multivariate hydroloqic time series analysis: C. T. Hsu and K. Adamowski Table l (continued). Summary o] general statistics of hydrologic time series involved in multirariate watershed system identification Station identification

08KA005

Drainage area Met. Stations available

08KA007

2690 3

615 1 Hydrological parameters

Precipitate

Runoff

38.0 ~F 254.4 [ F) 2

2.04 in 0.93 in 2

3.05 in 11.6 in 2

35.5 F 262.9 [ F):

2.58 in 2.08 in 2

73.65 16.06 3.25 0.23 0.06 0.01 93.26

92.10 0.68 0,09 0.05 0.34 0.14 93.41

19.47 4.21 0.47 0.25 0.42 0.88 25.70

67.15 18.65 5.31 0.97 0.13 0.05 92.25

92.54 0.90 0.15 0.10 0,34 0.21 94.23

24.82 4.25 0.03 0.03 0.72 0.60 30.46

0.56 in 2

16.8 ~F) 2

0.69 in 2

0.90 in 2

15.1 ( F ) 2

1.45 in 2

1.00 0.33 0.19 0.31 -0.001 -0.02 0.02 -0.02 0.05

0.31 - 0.01 -0.05 1.00 0.27 0.12 -0.33 -0.13 - 0.001

0.02 0.15 0.10 -0.33 -0.16 -0.11 1.00 -0.03 - 0.05

1.00 0.38 0.16 0.27 0.001 -0.001 0.01 0.05 0.08

0.27 0.001 -0.10 1.00 0.29 0.12 -0.33 -0.10 - 0.05

0.01 0.14 0.08 -0.33 -0.20 -0.11 1.00 -0.001 0.02

Statistical parameters

Runoff

Mean Variance

2.90 in 8.32 in 2

Explained seasonal variance ("o!

12 months 6 4 3 2.4 2 Total

Residual variance Auto and cross-corr, coeff.

Runoff

Temperature Precipitation

Zero order I st 2nd Zero order

1st 2nd Zero order 1st 2nd

Temperature

Station identification Drainage area Met. Stations available

Temperature

08KB001

08LA008

12500 20

1780 4

Precipitation

Hydrologic parameters Statistical parameters

Runoff

Temperature

Precipitation

Runoff

Temperature

Precipitation

Mean Variance

2.68 in 4.92 in 2

37.8 ~F 252.6 ( F ) :

2.24 in 8.01 in 2

0.75 in 0.75 in 2

40.2 ' F 22.1 ( F ) 2

1.84 in 0.95 in 2

67.19 15.03 5.09 0.37 0.37 0.01 88.05

92.68 0. 75 0.08 0.03 0.24 0.17 93.95

16.46 I 1.85 2.38 0.08 4.05 1.31 36.13

46.03 24.28 8.78 2.21 0.63 0.08 82.01

92.42 0.28 0,02 0.09 0.37 0.34 93.52

12.86 12.77 3.07 2.07 7.02 1.31 29.10

0.59 in 2

15.3 ( F ) 2

0.51 in 2

0.13 in 2

14.4 ( F ) 2

0.58 in 2

1.00 0.43 0,17 0.20 0.03 -0.02 O. 19 0.07 - 0.06

0.20 -0.09 -0.06 1.00 0.29 0.13 - 0.35 -0.17 - 0.09

0.19 0.29 0.10 - 0.35 - 0.14 -0.12 1.00 -0.03 - 0.04

1.00 0.66 0.46 0.04 0.04 0.07 0.003 -0.10 - O. 15

0.04 -0.09 -0.14 1.00 0.34 0.16 - 0.35 -0.10 - 0.09

0.003 0.25 0.15 - 0.35 - 0.28 -0.09 1.00 0,04 0,003

Explained seasonal variance I",,)

12 months 6 4 3 2.4 2 Total

Residual variance

Auto and cross-corr. coeff.

Runoff

Temperature Precipitat ion

Zero order 1st 2nd Zero order I st 2nd Zero order 1st 2nd

listed in Table 1. It is observed that all the runoff series fall into the autoregressive-moving average scheme. However, temperature and precipitation showed no obvious deviation from a purely random process except for temperatures in the Pacific mountainous region indicating persistence in the time series. Smoothed autospectrums for the transformed runoff and meteorological series were computed using equation (14) and are illustrated in Fig. 3, whose patterns further confirm the above observation.

88

Adv. Water Resources, 1981, Volume 4, June

Cross-correlations among the transformed runoff, temperature and precipitation series were estimated using equations (11) and (12) and are illustrated in Fig. 4. Crosscorrelations up to the first two lags are given in Table 1. It is observed that the cross-correlations between the runoff and the meteorological series are not stable except at the three stations 05AB002, 05AB021 and 05GG001. These are the three watersheds with denser meteorological network (Table 1), which provides more representative pictures of the meteorological condition in the watershed.

Multivariate hydrologic time series analysis: C. Z Hsu and K. Adamowski

Co-spectra, quadrature spectra, phase functions and coherencies were computed using equations (18), (19), (20) and (21) respectively and are illustrated in Figs. 5 and 6. They generally support the findings from the above discussion on cross-correlation functions, i.e. meteorological series generally lead the runoff series (negative phase spectrum over most frequencies), significant coherencies exist at low frequencies between runoff and meteorological series and among meteorological series. As in multiple regression analysis, frequency response functions serve as the regression coefficients in the frequency domain. They measure the relative contri-

4.0~

RUNOFF ~ PI~CIP - -

~0©

/

A

vz~

\

.//

.

.

-.

.

"F

~- . . . .

i¢~

.i

w

0,~0

,

F

i

A

M

,

,

.

J

J

A

o.

S

0

N

D

MONTH

_

_:__

o,___ ~e

io.

---;2Y-z'-\

/

" ,.f:

F -~c,~

- - - -

E o~o,/ X

o,i

o.ool

i

~

RUNOFF PR1ECIP I

o.~

I

TE/~' .... o.o~

F

M

A

M

J

b

J

A

I

O

N

•

IO.0

a"

-ore.

Figure 1. (a) Monthly means; (b) monthly variances i.

/N

u

RUNOFF ,(

o.,1~o o.~

,

~

,

,

,

,

-o.l~

LAGIkK~ITHS) l I

I I

I

l

ii

80

FREQUENCY (CYCLES PER 4 e M O N T H S )

Figure 3.

Variance spectrum

TEMPERATURE

i~ooo o.~oo o.~oo

TEMP-- PR(CIP o.eoo oAoo

o,oo~

o.o~: (MONTHS) -O,II~

V

:,sJ ~ . t o J

-~

i.o~o

~-5-.~"~

LAG (MONTHS)

.O.~I

PRECIPITATION TEMP-- RUNOFF

o.~oo

oJoo

o,o N

~

°

v

to

~

fs ~ ' ~ ' "

eo

"~/

Ves

L A G (MONTHS)

Figure 2.

Autocorrelations

It is also observed from Table 1 and Fig. 4, however, that the meteorological series generally lead the runoff series, i.e. peaks occur at positive lags. Cross-correlations between the temperature and precipitation series are stable. This is expected because both temperature and precipitation are measured at the same station. All the meteorological time series have high cross-correlations at the zero lag.

f

LAG (M(~THS)

*.ooo o.eoo

5v

~

,P"

.o

-'~.

PRECIP-- RUNOFF

" x ..,-~,X'x_. .+"-'xM .go

.~s

-*o

.s " ~ J

LAG |MONTHS) o ~

.o,lmo

Figure 4.

Cross-correlations

Adv. Water Resources, 1981, Volume 4, June

89

Multivariate hydrologic time series analysis: C. 77 Hsu and K. Adamowski Is.oo

,

,

,

, ,,,,,

,

i

i

i

damped sine wave, which is the characteristic of an ARMA transfer function. Phase functions, however, were hardly resembling any of the known systems, let alone to detect absolute delay in physical time. This was also confirmed by Dowling ~3 and a time domain system response function is therefore required for the determination of physical time lags and for the construction of a time domain transfer function model which is to be discussed in the following sub-section. Partial coherency spectrums were computed using equation (33) and are illustrated in Fig. 9. As expected, they were lower than the corresponding cross-coherency

i , * l

,111"

TEMP-- PRECIP

iO~ 31r

5.00 Z

•

0.~

IlL

-Z'~"

-I0.~

O.OI

O. I

O.S

16.00

FREQUENCY (CYCLES PER MONTH)

4it

I~.00

io.c¢ s-it

TF_EE~ - - RUNOFF . . . .

t,r

Io .00 P R E Q P

I

RUNOFF

.

I

Z

,/'/'~J~,,:,,

- - .

t~r

Z

"r

e~

~

\

/

\

°~

~

.~-

G.

- ! ,rr"

<

RUI~)FF-TEMP OEK~EIq'VED IST

-31r

0,.

2ND

,llr 95%

-41i,"

-Io.1~

ORDERARMA

~

ORDER

. . . .

CONFIDENCE

ARMA LIMITS

~

.

.lll.~

4~

O, Oi

.l&l~

i

*

i

i ~ l i l *

,

,

i

O, I

O.!

0.11

FI~EQU~NCY(CYCLESPER MONTH)

Figure 5.

03

FREQUENCY (CYCLES PER MONTH)

i i , i

l&O0

Phase junctions

i

i

i

illlll

i

i

t

iii

4~r io.o0

211"

I~O T E MP--

OIK>O

PRECIP

¢.1

OlOO

0700 o.10o

oz

OlO0

"

.

~,p--d

o.loo

-IO, O0 o,~o

I

I

I

i

I

5

Io

18

ItO

i4

2NO ORDER A R ~

-41r

95% CO~'I~E

ollLO FREQUENCY

(CYCLES

PER

48

MONTHS)

i

i

I

....

LiMiTS ~ _ ~

I i iiii

O~OI

i

I

O.I

I I i O,5

(CYCLES PER MONTH)

F R ~ C Y

Figure 7.

I

Phase of./i'equeney response Junctions

J.ooo TEMP

-- RUNOFF

. . . .

PRECIP

-- RUNOFF

~

oloo G~o

,

io °°"-=

"°",o= I " .... , r ....

R U N O F F- T I I ~ P 0 B ~ R ~ 4 ~

r• ,

\ 4 \~4"

,>7 v

,' ,'

A I\

{CYCLES

2ND ~

14

AR uA ....

,

,,,,,,

,

,

,,,,,i

- -

Y PER

\

/

,

~.-

95% CONFI0~NCE LII~TS __.__

Squared eohere,cie.s

Adv. Water Resources, 1981, Volume 4, June

\

\

48 MONTHS)

--

2 NO O R O [ R A R M A . . . .

9 5 % CONFIDIrNCIr L I M T S - - . -

/, (-W-,, \.-,, / ,,~

"-,,"

but(on of variance to the output process from one of the input processes allowing for the effect of other input processes at various frequencies. Gain and phase for the frequency response functions were estimated from equations (25) and (26) and the results are illustrated in Figs. 7 and 8. Most of the gain functions revealed a system of exponential decay or

90

,

RUNO~r-pREClP ~

I ST 0R0(R ARMA --

x, FREOUENCY

Figure 6.

\

- -

i

!

~(

. . . . . . . . .

0OI

I

Oi

I

I

I

I ;]1

O5

FREQUENCY ( ~ r CLES PiER M O N T H )

Figure 8.

Gains of jrequence

O~l o Ol

/\ A o,

ot

FRE~JENCYICYCI.ES PER klONTH)

response junctions

Muhirariate hydrologic time series analysis: C. Z Hsu and K. Adamowski UPO'ER1'5% CONFI~NC[t . i r

iooo oloo i

o.moo

~,= .... a.z

v,, " ? . - : - ,

o.i~

gao,

m

oloo 5

io

is

zo

14

FREQUENCY (CYCLES PER 48 MONTHS}

o.)oo ~z ~eoo

TEMP -RUNOFF . . . . PRECIP-RUNOFF - -

.Jw

po

,./ oloo

o.~

/%,',,',..,A.

:.:-.... ,

,,,

,ll \

, ',.-'i,',',,./'?,._'_:.4\ s

io

is

Io

i4

FREQUENCY (CYCLES PER 48 MONTHS )

Figure 9.

Multiple and partial squared coherencies

spectrums (Fig. 6) but not significantly different. Multiple coherency spectrums between the runoff and the meteorological series were obtained from equations (30) and (31) and are also illustrated in Fig. 9. It can be seen from Fig. 9 that the coherency was generally improved through the multiple correlation regression of the runoff on the meteorological time series, i.e. temperature and precipitation. It can also be observed that gain spectrums of the frequency response functions and multiple coherencies all indicate high values at low frequencies around the annual cycle (Figs. 7, 8 and 9). Recalling the fact 9 that monthly runoff series fall into a first-order autoregressive scheme which has rich variance at low frequencies, the above observation further confirms the finding in that previous paper that a natural watershed is a 'damped" system which damps out the high frequency variations of the meteorological inputs through its storage effect. Impulse response functions were obtained by taking the Fourier transformation of the corresponding frequency response functions using equation (29) and the results are illustrated in Fig. 10. Their patterns did follow the form of exponential decay or damped sine wave. This justifies the postulation of an ARMA type transfer function model. Physical time lags were determined based on the criterion (9). Two types of ARMA transfer functions were attempted, namely, first-order and second-order. Initial values of the transfer function parameters were estimated, using the relationship (34), respectively for the two models. The values for the noise model parameters were obtained by fitting ARMA(1, 0) and ARMA(2, 0) processes respectively to the two noise series. The initial model parameters are listed in Table 2. The final values of the model parameters were then obtained through model parameter optimization based on Rosenbrock's hill-climbing techniques. The final model parameters and their approximate standard deviations are also given in Table 2. It can be seen that most of the watershed systems could be represented by both

first-order and second-order ARMA transfer function models (Table 2 and Figs. 7, 8, 10 and 11). Furthermore, from the results in Table 2 it can be observed that for the Prairie region the initial parameters generally agree remarkably well with the final parameters after parameter optimization. The estimates of time lags are also consistent in both cases. Therefore, system impulse response function seems to be an efficient estimator for the estimation of transfer function model parameters. The estimation of model parameters was also attempted using least square fitting directly on the data, which equivalently is to solve the simultaneous WienerHopf equations using auto and cross-correlations of the involved time series. The results, however, were not as good and consistent and therefore are not presented here. It is also observed from Table 2 that the improvement from fitting a second-order ARMA transfer function model instead of a first-order model to the watershed system of monthly runoff, temperature and precipitation series is not significant. By examining the corresponding standard deviations of estimates, higher order weights in the second-order model are relatively small and hardly to be of any statistical significance. However, for the Pacific mountain stations, (he initial and final estimates of transfer function parameters as well as the two sets of model parameters are not consistent. This may be attributed to the scarcity of the available meteorological data in these mountainous basins, with the meteorological network density from 445 to 897 (sq. miles per station) and most of meteorological data sampled at lower to medium band of the elevation (1465 to 4180 It). It is further observed that the weight for the previous state of runoff (Table 2) is one of the same order of magnitude as the autoregressive coefficient for the univariate runoff model l~'. This leads to the postulation that the disturbance term in the univariate model may be modelled by the terms at the right-hand side of the ARMA transfer function model. By comparing the explained RUNOFF-TEMP OBSERVED IST ORDER ARMA ~ ORDER ARMA . . . . 9 5 % CONFIDENCE LIMITS ~ - ~

O3O0

2ST

oloo i

~Joo

d'K " ' / , ""JA

~ o.o~ -o.ioo -o.z~

d LAG (MONTHS)

.o,!u~

RUNOFF-TEMP OBSERVED PST ORDER A R M A - - ~

/~ osoo

--'~-J ~

o.loo

•

o.loo •

~

/,,

/%

"s

-%.,-

2 ST ORDER ARMA . . . . ",~'\

v \

9 5 % CONFIDENCE LIMITS - - - - -

'\ .\

-'~.1 ',, ~ , - - " v-../\

i

/

,,-._..

,_,-.1

/-,.-

,

-

-

1

,,v - " . \

,\

.o

_

,

i -..

-oaoo ~o.zoo •

LA~;H/S

)

Figure 10. Impulse response JUnctions

Adv. Water Resources, 1981, Volume 4, June

91

M u l t i v a r i a t e h y d r o l o g i c t i m e series a n a l y s i s : C. T. H s u a n d K . A d a m o w s k i v a r i a n c e s b y t h e u n i v a r i a t e r u n o f f m o d e P o t o t h o s e listed in T a b l e 2, t h e i m p r o v e m e n t f r o m c o o r d i n a t e d u s e o f t h e m e t e o r o l o g i c a l a n d h y d r o m e t r i c d a t a is o b v i o u s . F o r w a t e r s h e d w i t h l a r g e s t o r a g e , i.e. 0 5 T D 0 0 1 , h o w e v e r , t h e i m p r o v e m e n t is n o t p r o n o u n c e d . I n t h i s c a s e , t h e l a r g e s t o r a g e effect p l a y s a d o m i n a n t r o l e in t h e g e n e r a t i n g process of s y s t e m r e s p o n s e a n d the c o n t r i b u t i o n f r o m e x t e r n a l d i s t u r b a n c e s will b e d a m p e d o u t b e f o r e it b e c o m e s significant. Table 2.

Summary on multivariate watershed system identification

Station identification .

05AB002

xa: runoff; x2: precipitation; x,: temperature (1 - d , B)x3, = 1 " 1 o - wl IBlxll-h, +

(W2o- v"21B)x2t - h: +

dl wlO wl I bI W2o 'iA'21 b2

O~z, _ ~+ a,

01

( 1 - d I B - d: B2)x3, = (Wto - v,1t B - wt2B2)Xlt_~,, + (W2o - w , l B - - w22B2)x2t_h2 +

(01B + 02B2)z, + a,

Final estimate

Std. dev. estimate

Initial estimate

Final estimate

Std. dev. estimate

0.58 0.16 0.24 0 O.19 --0.08

0.65 0.06 0.23 0 0.30 --0.11

0.05 0.06 0.06

0.73 0.10 0.24 0 0.27 -0.09

0.05 0.06 0.06

0.06 0.06

0.80 0.20 0.28 0 0.20 -0.06

0.08

- 0.36

0

dI d2 'a'lo wl I wl 2 bI "'20 w21 "'22 b2 01 02

x3: runoff; x2: precipitation; x~: temperature

01ZI- I "+"t~t

0.43 0.21 0.16 0.22 0,002 0 0,19 -0,11 - 0.04 0 -0.10 -0.15

(Wlo -- Wl I B -- w I 2 B 2 ) x I t - b , +

(w2o - w2~B - w22B2)x2t-h: + (OiB+O2B2z, +a,

0

dI w¢o wll bI W2o w21 b2 (tI

0

- 0 32

0.08

yes 50 0.43 0.16 0.06 0.20 0.08 0 0.31 -0.19 - 0.005 0 -0.10 -0.08

0.23 0.17 0.06 0.07 0.08 0.06 0.09 0.02 0.24 0.10

0.75 0.06 0.20 0.26 0.007 0 0.20 -0.07 - 0.04 0 -0.37 -0.15

0.69 0.02 0.10 0.21 0.06 0 0.30 -0.09 - 0.03 0 -0.32 -0.10

yes 52

yes 52

05GG001

05TD001

0.10 0.07 0.06 0.06 0.06 0.06 0.07 0.06 0.12 0.09

Final estimate

Std. dev. estimate

Initial estimate

Final estimate

Std. dev. estimate

0.65 0.13 0.17 0 0.24 0.02 1 -0.11

0.72 0.10 0.16 0 0.25 0.02 1 -0.18

0.05 0.05 0.06

0.74 -0.11 -0.06 1 0.16 - 0.08 0 -0.09

0.85 -0.09 -0.02 1 0.05 - 0.03 0 -0.22

0.04 0.04 0.04

0.06 0.03 0.09

yes 52 dI d2 "'lo Wl I wl 2 bI w2o w2~ w22 b2 0, 02

0.06 0.06

Initial estimate

Model validity at 95% level Explained residual variance (')ol (1 -- d 1B -d2B2)x3, =

- 0.30 yes 50

Station identification

(W2o- w2~B)x2t - l,: +

0

- 0.24

Model validity at 95% level Explained residual variance (%;)

(1 -dlB)x3, =(Wio- w1l)xll_t,, +

05AB021

Initial estimate

Model validity at 95% level Explained residual variance (%)

0.69 - 0.02 0.13 0.18 0.08 0 0.24 0.03 0.10 1

A d v . W a t e r R e s o u r c e s , 1981, V o l u m e 4, J u n e

0.09 0.09 0.05 0.06 0.06 0.06 0.05 0.05

0.69 0.09 -0.11 --0.05 0.01 1 0.16 - 0.08 0.04 0

0

0.05 0.09

-0.06 -0.09

-0.16 -0.19

1

0.01 -0.16 yes 55

0.04 0.04 0.08

yes 67 0.57 0.20 0.08 0.11 O.13 0 0.23 - 0.02 0.11

-0.13 0.004

Model validity at 95% level Explained residual variance (%)

92

Under the assumption that the multivariate ARMA t r a n s f e r f u n c t i o n mod61 is j u s t t h e e x t e n s i o n o f a u n i variate ARMA process by introducing the contribution from meteorological inputs, some interesting findings can b e m a d e . C o n s i d e r i n g t h e c a s e o f A R M A ( 1 , 0) p r o c e s s fitted d i r e c t l y f r o m t h e t r a n s f o r m e d r u n o f f series, s t a n d a r d d e v i a t i o n s o f e s t i m a t e s o f a u t o r e g r e s s i v e c o e f f i c i e n t s for t h e e i g h t s t a t i o n s i n c l u d e d in t h e m u l t i v a r i a t e s y s t e m a n a l y s i s w e r e f o u n d r e s p e c t i v e l y t o b e 0.06, 0.06, 0.06, 0.04,

0.78 0.10 -0.09 --0.02 0.06 1 0.06 - 0.05 0.06

yes 68

0.20 0.20 0.04 0.03 0.04 0.04 0.04 0.04 0.21 0.08

M u l t i v a r i a t e h y d r o l o g i c t i m e series a n a l y s i s : C. T. H s u a n d K . A d a m o w s k i Table 2 [continued).

Summary on multivariate watershed system identification

Station identification

08KA005

x3: runoff; x2: precipitation; xl: temperature

( 1 - d I B)X3,

=

(W1o - - Wl 1 n ) x l t - h L+

(w2 o - w21B)x2t - h: +

01z,- 1+a,

dl Wlo wl i bI w20 w21 b2 01

Initial estimate

Final estimate

Std. dev. estimate

Initial estimate

Final estimate

Std. dev. estimate

0.35 0.49 0.12 0 0.20 -0.21 0 -0.02

0.54 0.44 0.22 0 0.15 -0.15 0 -0.21

0.08 0.07 0.08

0.27 0.40 0.03 0 0.13 -0.19 0 0.14

0.48 0.38 0.19 0 0.09 -0.15 0 -0.08

0.10 0.07 0.07

Model validity at 95~o level Explained residual variance (~o)

(Wlo-wllB-w12B2)xlt_~,+

(W2o-w21B-w22B2)x2t_~,2+ (01B + 02B2)z, + a,

0.17 0.10 0.49 0.03 - 0.02 0 0.20 -0.24 - 0.15

b2

0

0

01 02

0.14 --0.06

-0.02 -0.41

Station identification x3: runoff; x2: precipitation; xl: temperature wlIBjXlt_ht

(W2o -- w21B)x2t_h2 +

01z,- 1 +a,

+

dI Wlo wll bI W2o w21 b2 01

(Wlo - w I ~B -

w 12BZ)Xlt

-h, +

(W2o - w21B - w22B2)x2t _b2 + (01B + 02B2)z, + a,

Model validity at 95~o level Explained residual variance (~o)

0.29 0.40 0.44 0.12 O.15 0 0.14 -0.18 - 0.05

0.08 0.09 0.07 0.07 0.08 0.07 0.07 0.06

0.29 -0.08 0.40 0.04 0.03 0 0.13 -0.19 - 0.03 0

0.06 0.09

O.lO 0.11

0.44 0.15 0. 37 0.14 O.13 0 0.09 -0.17 0.0005

0.15 0.13 0.07 0.08 0.08 0.06 0.07 0.004

0

0.06 --0.21

yes 35

yes 32

08KB001

08LA008

0.15 0.10

Final estimate

Std. dev. estimate

Initial estimate

Final estimate

Std. dev. estimate

0.52 0.41 0.14 0 0.40 - 0.23 0 -0.18

0.41 0.36 0.15 0 0.28 - 0.22 0 -0.05

0.09 0.07 0.07

0.74 0.17 0.10 0 0.26 0.02 1 -0.12

0.72 0.10 0.09 0 0.25 0.05 1 -0.10

0.07 0.06 0.06

0.06 0.07 0.11

yes 38

dI d2 wlo Wl 1 Wl 2 bI W2o w, t w, 2 b, 01 02

0.12

Initial estimate

Model validity at 95~0 level Explained residual variance (~o) (1 - d l B - d 2 B 2 ) x 3 , =

0.10

0.07 0.07

? 29

dI d2 wl0 wll wl 2 bI W2o w21 w22

Model validity at 95~o level Explained residual variance (~o)

-

0.07 0.07

yes 30

(1 - d , B - d 2 B 2 ) x 3 t =

(1 -dlB)x3, =(Wlo

08KA007

0.34 0.13 0.41 0.06 -- 0.05 0 0.40 - 0.30 - 0.07 0 - 0.02 -0.25

0.09

yes 52 0.03 0.42 0.39 --0.02 0.08 0 0.31 - 0.39 - 0.13 0 0.23 -0.50

yes 44

0.06 0.05

0.04 0.07 0.06 0.02 0.06 0.06 0.06 0.07 0.07 0.08

0.44 0.35 0.18 0.51 0.04 0 0.26 - 0.05 0.03 1 0.20 -0.26

0.39 0.33 0.10 0.05 0.05 0 0.27 - 0.06 - 0.03

0.13 0.12 0.05 0.05 0.05 0.06 0.07 0.05

I

0.23 -0.26

O.13 0.11

yes 53

Model validity at 95~o level was based on various statistical tests at 95~o confidence level on comparison between observed and theoretical impulse response functions and frequency response functions, on cross-correlations between the noise and input series and on randomness of the noise series. "T indicates failure in any of the above tests to validate the proposed model. Standard deviations of estimates of time lags "b" are not computed because they can take only integer values. 0.07, 0.07, 0.07 and 0.06. Comparing these values with the corresponding values given in Table 2 for the first-order A R M A t r a n s f e r f u n c t i o n m o d e l , it c a n b e o b s e r v e d t h a t for Prairie stations the autoregressive coefficients obtained from the multivariate system analysis have smaller standard errors than those obtained from the univariate

r u n o f f a n a l y s i s e x c e p t f o r t h e l a r g e s t o r a g e w a t e r s h e d , i.e. 0 5 T D 0 0 1 , t h e i m p r o v e m e n t is n o t s i g n i f i c a n t . I n c o n t r a s t , t h e P a c i f i c m o u n t a i n s t a t i o n s all i n d i c a t e h i g h e r s t a n d a r d errors in the multivariate case due to the noise introduced by insufficient sampling of meteorological data in these basins. Thus, with reasonable coverage of meteorological

A d v . W a t e r R e s o u r c e s , 1981, V o l u m e 4, J u n e

93

Multivariate hydrologic time series analysis: C. T. Hsu and K. Adamowski • o.soo ~

NOISE - - NOISE

9 5 % CONFIDENCE L I M I T S ~ - ~ FOR WHITE NOISES (W.N.)

[4)

0.41~

-osoo

(5) i.ooo TEMP - - NOISE

o.ooo O

O.SOO

Oo

..o.zoo

.

o~

.

.

.

.

9 5 % CONFIDENCE LIMITS ~ FOR TWO UNCORRELATED PROCESSES WITH EITHER OR BOTH BEING W.N.

,,,.

__,7"~

(6)

..,.._,^7,

~u) -o4oo Oo~ °o~Loo

(7)

-o.1~o °1.~o

~aoo

PRECIP m NOISE

o.~ B

klJ

9 5 % CONFIDENCE LIMITS ~ - ~ FOR TWO UNCORRELATED PROCESSF..S WITH EITHER OR BOTH BEING W,N.

0~00

(8)

Figure 11.

Auto- and cross-correlations

network in the area, the coordinated use of meteorological and hydrometric data also improves the estimation of parameters in the univariate runoff model. In conclusion, a first-order ARMA transfer function model is considered adequate to model a watershed system of monthly runoff, temperature and precipitation time series. This coordinated use of the meteorological and hydrometric data has been shown to enhance the accuracy of estimation of runoff characteristics. CONCLUSIONS (11

(2)

(3)

94

The transformed stationary runoff time series fall into the autoregressive-moving average regime. However, the transformed stationary temperature and precipitation time series showed no significant difference from a purely random process except for temperatures in the Pacific mountainous region indicating persistence in the time series. Significant cross-correlation exists between temperature and precipitation at the zero lag. Crosscorrelations between the runoff and the meteorological time series are in some cases not stable because of the scarcity of the available meteorological data. In general, meteorological series were found to lead the runoff series. Most of the gain functions revealed a system of exponential decay or damped sine wave, which is the characteristic of an ARMA transfer function. Phase functions, however, did not resemble any of the known systems. Thus, a time domain system response function is required for the determination

Adv. Water Resources. 1981, Volume 4. June

(9)

of physical time lags and for the construction of a time domain transfer function model. Most of the significant coherencies and gains of the frequency response functions occur at low frequencies around the annual cycle. This suggests a "damped" watershed system which transforms the meteorological inputs into the runoff with rich variance at low frequencies. Impulse response function was found to be an efficient estimator for the estimation of transfer function parameters and for the determination of physical time lags. The multivariate watershed system of monthly runoff, temperature and precipitation time series may be adequately represented by a first-order ARMA transfer function model. However, model instability may be introduced due to the insufficient spatial sampling of meteorological data. The weight for the previous state of runoff in the multivariate system transfer function is of the same order of magnitude as the autoregressive coefficient in the univariate runoff model. This leads to the postulation that the disturbance term in the univariate model may be modelled by the terms at the right-hand side of the ARMA transfer function model. Coordinated use of the meteorological and hydrometric data would enhance the accuracy of estimation of runoff characteristics. For watershed with large storage, the improvement was not pronounced. In this case, the large storage effect plays a dominant role in the generating process of system response and the contribution from external disturbances will be damped out before it becomes significant. For multivariate watershed hydrologic system modelling, a combined deterministic-statistical approach could be developed. For example, instead of considering temperature and precipitation as separate input variables, a transformed variable of net precipitation obtained by subtracting from the precipitation the amount of evapotranspiration computed from some known physical or empirical equations such as Penman's equation 14, Turc's formula ~5 etc. may be used as the input variable to the watershed system model. This will establish more meaningful link between the univariate and the multivariate system models and provide more insight in the physical significance of the proposed stochastic system model. The identification of the watershed hydrologic system should not be considered as the ultimate goal of the hydrologic system analysis. The practical application of the proposed models is at least two-fold: to estimate the runoff characteristics or synthesize the runoff time series at the ungaged site through the correlation regression between model parameters and watershed physiographic characteristics 16, and to make short-term runoff forecast from the proposed system models following the procedures demonstrated by Box and Jenkins 6.

ACKNOWLEDGEMENTS This project was carried out with financial assistance

Muhit'ariate hydrolo#ic time series analysis: C. Z Hsu and K. Adamowski

provided by the National Research Council of Canada. grant A7467.

7

8 9

REFERENCES 1 2

3 4 5 6

Solomon, S. I. Joint mapping, WMO Casebook on Hydrological Network Design Practics, Sect. III-2.1, 1972 Coulson, A. Estimating runoff in southern Ontario. Tech. Bull. No. 7, Inland.Waters Branch, Department of Energy, Mines and Resources, Ottawa, 1967 Crawford, N. H. and Linsley, R. K. A conceptual model of the hydrologic cycle, IASH Publ., No. 63, 1964, pp. 573-587 Hsu, C. T. A mathematical model for continuous basin runoff simulation, MSc Thesis, Queen's University, (1970) Jenkins, G. M. and Watts, D. G. Spectral Analysis and its Application, Holden Day, San Francisco, 1968 Box, G. E. P. and Jenkins, G. M. Time Series Analysis Forecastin9 and Control, Holden Day, San Francisco, 1970

10 11 1"2 13

14 15 16

Shawinigan Eng. Co. Hydrometric network plannmy stud)' lor western and northern Canada, Report for the Canada Department of Energy, Mines and Resources, 1970 Hannon, E. J. Multiple Time Series, John Wiley, New York, 1970 Hsu, C. T. and Adamowski, K. Physically based stochastic model of a watershed system. Int. Syrup. Rait!lalI-Runq[L Mississippi State Unit'. 1981 Hsu, C. T. A systems approach to hydrologic time series analysis, PhD Thesis, University of Ottawa 1197% Rosenbrock, H. H, An automatic method for finding the greatest or least value of a function, Computer J. 1960, 3, 173 Wilde. D. J. Optin, um Seekim! Method. Prentice-Halt, Englewood Cliffs, 1964 Dowling. J. M. A note on the use of spectral analysis to detect leads and lags in annual cycles of water quality, 14ater Resourc. Res. 1974. 10. 343 Penman, H. L. Natural evaporation from open waters, bare soft and grass, Proc. R. Soc. 1948. A193. 120 Turc, k. Le bilan d'eau des sols. Relations entre les precipitation, I'evaporation et l'ecoulement. Ann. A qron. 1954, 5, 491 Pentland, R. L. and Cuthbert, D. R. Operational hydrology for ungauged streamflows by the grid squares technique, Water Resourc. Res. 1971, 7. 283

Adv. W a t e r R e s o u r c e s , 1981, Volume 4, J u n e

95