System Identification and Simulation of Chesapeake Bay and Delaware Bay Canal Hydraulic Behavior

System Identification and Simulation of Chesapeake Bay and Delaware Bay Canal Hydraulic Behavior B.B. Hsieh Maryland Department of Natural Resources, Annapolis, M D 21.401, USA

ABSTRACT The Kalman f i l t e r combined w i t h one-dimensional s h a l l o w w a t e r wave e q u a t i o n i s q p p l i e d t o d e s c r i b e t h e t i d a l h y d r a u l i c behavi o r of t h e Chesapeake Bay and Delaware Bay c a n a l (C&D c a n a l ) . T h i s model j o i n s two s e p a r a t e e s t u a r i n e systems and one branched i n t h e Elk River a r e a . The f i e l d measurements from s u r f a c e w a t e r e l e v a t i o n s and t i d a l c u r r e n t s are used t o c a l i b r a t e t h i s s t o c h a s t i c - d e t e r m i n i s t i c model. The agreement between f i e l d o b s e r v a t i o n s and f i l t e r e d e s t i m a t e i s b e t t e r t h a n f i e l d and t h e p u r e l y numeric a l s o l u t i o n . The model performance i s e v a l u a t e d by t h e s t a t i s t i c a l methods. The c o n s i d e r a t i o n s f o r e x t e n d i n g t h e modeling proc e s s t o t h e s y s t e m n o i s e and t h e measurement n o i s e p r o v i d e much b e t t e r p r e d i c t i o n s and c a n be used t o s i m u l a t e t h e s p e c i a l f o r cing events. INTRODUCTION

The i n c r e a s e of s h i p p i n g a c t i v i t i e s through t h e waterway and t h e c o n c e r n of baywide r e s o u r c e s , r e v e a l e d t h e i m p o r t a n t i s s u e of t h e dynamic of t h e C&D c a n a l . P a r t i c u l a r l y , t h e c o r r e c t e s t i m a t i o n of t i d a l f l u x from t h i s i n t e r s y s t e m c o u l d a d d r e s s much more t r a n s p o r t b e h a v i o r of t h e d i s t r i b u t i o n of a contaminant o r a f i s h egg and l a r v a e , Hydrodynamic e q u a t i o n s have been s u c c e s s f u l l y used by many r e s e a r c h e r s ( G a r d n e r and P r i t c h a r d ( 1 9 7 4 ) ; H u n t e r ( 1975); R i v e s and P r i t c h a r d ( 1 9 7 8 ) , among o t h e r s ) f o r c o n s t r u c t i n g a d e t e r m i n i s t i c n u m e r i c a l s o l u t i o n of t h e C&D c a n a l . However, a s u b s t a n t i a l e r r o r c o u l d be produced by t h e p u r e l y s t o c h a s t i c met e o r o l o g i c a l f o r c i n g . This c o n t r a s t t o t h e p u r e l y s t o c h a s t i c model i s i n c a p a b l e of p r o v i d i n g p r e d i c t i o n s o v e r t h e s p a t i a l l y v a r y i n g d i s t r i b u t i o n of l o n g wave p r o p a g a t i o n . B u d g e l l ( 198 1 ) developed a s t o c h a s t i c - d e t e r m i n i s t i c model f o r a p p l y i n g t o t h e G r e a t Bay e s t u a r y . T h i s t y p e of approach c o n s i d e r s t h e s y s t e m n o i s e , such a s t h e s p e c i f i e d c r o s s - s e c t i o n a l a r e a , flow f i e l d , boundary c o n d i t i o n s , and t h e measurement n o i s e i s m a i n l y caused by i n s t r u m e n t s e n s i t i v i t y . B u d g e l l ( 1982) i n d i c a t e d w i t h o u t C O K -

313

314 r e c t i o n s , t h i s modeling e r r o r c o u l d be propagated through time and s p a c e by t h e d e t e r m i n i s t i c model. A means of computing t h e o p t i m a l c o r r e c t i o n t o be a p p l i e d t o t h e computed y a l u e s a t e a c h time s t e p i s t h r o u g h t h e use of t h e Kalman F i l t e r . I n t h i s p a p e r , t h e s t o c h a s t i c f i l t e r d e r i v e d from t h e Kalman f i l t e r i s c o n s t r u c t e d around t h e one-dimensional e q u a t i o n s gov e r n i n g motion i n An open c h a n n e l . The boundary c o n d i t i o n s a r e u s i n g t h e s u r f a c e w a t e r e l e v a t i o n a t Reedy P o i n t a t i t s e a s t e r n end and a t Old Town P o i n t a s i t s w e s t e r n e n d . One t i d a l measurement s t a t i o n a t Chesapeake C i t y and t h e c u r r e n t m e t e r s a t Summit Bridge are used f o r model c a l i b r a t i o n and f i l t e r e d estimate c a l c u l a t i o n . The a p p l i c a t i o n of t h i s model c a n be extended t o d e s i g n f u t u r e m o n i t o r i n g programs and t o s i m u l a t e t h e s t o r m e v e n t which may c a u s e abnormal t i d a l f l u c t u a t i o n s . MODEL DESCRIPTION

The one-dimensional hydrodynamic e q u a t i o n s ( t h e c o n t i n u i t y equat i o n and momentum e q u a t i o n ) combined w i t h a s e t of e q u a t i o n s def i n i n g flow c o n d i t i o n s a t j u n c t i o n s and b o u n d a r i e s c o n s t i t u t e the d e t e r m i n i s t i c component of t h e model. The m a t r i x form may be e x p r e s s e d a s A(n,n+l) X(n+l)

=

B(n,n+l)X(n) + G(n,n+l) U(n+l)

(1)

where A ( n , n + l ) and B ( n , n + l ) a r e 2Nx2N c o e f f i c i e n t s m a t r i c e s which depend on t h e schemes and N g r i d p o i n t s . T X (n+l) =(ul,Vl, u ,qn) 2NxI v e l o c i t y and s u r f a c e e l e v a t i o n v e c t g r s a t time l e v e l ( n + l ) .

...,

X(n)

=2NxI v e l o c i t y and s u r f a c e e l e v a t i o n v e c t o r a t t i m e l e v e l n.

U(n+ I )

= v e c t o r c o n t a i n i n g boundary c o n d i t i o n s a t t i m e level (n+l).

G ( n , n + l ) = Nxr m a t r i x which s p e c i f i e s which e q u a t i o n s a r e r boundary c o n d i t i o n s . E q u a t i o n ( I ) can be e x p r e s s e d i n an a l t e r n a t i v e form u s i n g t h e i n v e r s e of A(n , n + 1 ) . X ( n + l ) = g ( n , n + l i X(n) + G * ( n , n + l ) U ( n + l )

(2)

where & & n , n + l ) = A I : ( n , n + l ) B ( n , n + l ) G (n,n+l)= A (n,n+l) G(n,n+l) J a w i n s k i ( 1970) developed t h e a l g o r i t h m of t h e d i s c r e t e system f o r computing an o p t i m a l f i l t e r based on o b s e r v a t i o n s . The p r e d i c t i o n s t e p c o n s i s t s of X ( n + l / n ) = s ( n , i l + l ) X(n/n) P ( n + l / n ) = p ( n , n + ~ )P ( n / n )

+ G*(n,n+l) U ( n + l ) T

5

(n,n+l) +Tq(n+l)fkT

(3)

(4)

The measurement u p d a t e s t e p c o n s i s t of X(n+l/n+l) = X(n+l/n)

+ K(n+l)(Z(n+l)-

H X(n+l/n))

(5)

315 P(n+l/n+l) = P(n+l/n)

-

(6)

K(n+l) H P ( n + l / n )

where E : ( n + l ) Q(n+ I )

= t h e Kalman Gain =system n o i s e c o v a r i a n c e = t r a n s i t i o n m a t r i x of t h e system n o i s e R(n+ I ) =measurement n o i s e c o v a r i a n c e Z(n+ 1) =observations H =measurement t r a n s i t i o n m a t r i x X ( n + l / n ) =one-step ahead p r e d i c t i o n v e c t o r P ( n + l / n ) = t h e e r r o r c o v a r i a n c e of X ( n + l / n ) X ( n + l / n + l ) = t h e f i l t e r e d e s t i m a t e of X ( n + l ) P ( n + l / n + l ) = t h e e r r o r c o v a r i a n c e of X ( n + l / n + l )

T

EQUATIONS

FIELD HEASURE.

_INITIAL _ _ _C.lloUNDARYc.

-

P

ONE-STEP AHEAD

NLMERICAL

ONE-STEP AHEAD PREDICTION ERROR

APPRROXIPY\TlON . L PREDICTION

T

1

NOH-LIN TERY ITERATION

C-

FILTERED __*

ESTIIKhTE

4

HEASURE.

-INITIAL - - -'.-

'

'"

REDICTION C

ERROR COVARIASCE

1

1 , KAUIAN GAIN

c

nEASUREHENT NOISE ESTI~MATION

I

SYSTEX

XOISE

COY.

F i g u r e I . Block diagram f o r a s t o c h a s t i c - d e t e r m i n i s t i c model. The C&D CANAL STUDY I n o r d e r t o c a l c u l a t e t h e t i d a l e l e v a t i o n and c u r r e n t v e l o c i t y from t h e hydrodynamic e q u a t i o n s , t h i s s y s t e m i s d i v i d e d i n t o a number of r i v e r r e a c h over a branched e s t u a r i n e waterway. ( F i g u r e 2 ) The upper p o r t i o n of t h e E l k River i s c o n s i d e r e d as the extending t i d a l behavior a r e a . A f t e r t h e cubic s p l i n e i n t e r p o l a t i o n p r o c e s s and t h e e s t i m a t i o n of c u r r e n t v e l o c i t y from 6 mooring

316 meters a t t h r e e !.evels, a d a t a s e t w i t h 3 t i d a l h e i g h t s , one f r e s h w a t e r i n f l o w and c r o s s - s e c t i o n a l mean v e l o c i t y i s c r e a t e d . T h i s d a t a f i l e c o v e r s t h e p e r i o d of June 14, 1984 t o J u l y 5 , 1984 with half-hourly b a s i s . Four p o i n t s i m p l i c i t f i n i t e - d i f f e r e n c e scheme w i t h w e i g h t i n g f a c t o r e q u a l t o 0 . 5 5 i s used t o s o l v e t h e m a t r i x system. Twenty p a i r s of c o n t i n u i t y e q u a t i o n s a n d momentum e q u a t i o n s p l u s t h r e e boundary c o n d i t i o n s and t h r e e j u n c t i o n e q u a t i o n s form a 46x46 m a t r i x s y s t e m . The i n s t a b i l i t y f o r i n v e r s e t h i s m a t r i x is u s u a l l y o b t a i n e d due t o band w i d t h of m a t r i x i s produced. The t e c h n i q u e f o r decomposing a m a t r i x i n t o upper t r i a n g l e and lower t r i a n g l e c a n avoid d i r e c t i n v e r s e c o m p u t a t i o n . The n o n - l i n e a r i t y of t h e momentum e q u a t i o n i s s o l v e d by t h e Newton-Raphson method. (Mickle and S z e , 1972)

J

F i g u r e 2 . The C&D Canal E s t u a r i n e System

i

SYSTEM PARAMETERS The s p e c i f i c a t i o n s of t h e model s t r u c t u r e a r e i m p o r t a n t paramet e r s t o be i d e n t i f i e d b e f o r e t h e model i s r u n . system noise covariance I n o r d e r t o o b t a i n t h e independent s y s t e m n o i s e c o v a r i a n c e

t h e d a t a s e t i s d i v i d e d i n t o two d a t a f i l e s . F i r s t s e t ( l 0 t i d a l c y c l e s ) i s used t o g e n e r a t e t h e s y s t e m n o i s e c o v a r i ance w h i l e t h e second s e t i s used t o c a l i b r a t e t h e model. The ARIMA p r o c e s s i s used t o g e n e r a t e t h e white n o i s e from t h e f i e l d d a t a f o r both v a r i a b l e s . I t i s assumed t h a t t h e system n o i s e c o v a r i a n c e f o r c u r r e n t v e l o c i t y remains t h e s a m i n a l l g r i d systems The system n o i s e c o v a r i a n c e f o r t i d a l h e i g h t s w i l l remain t h e c o n s t a n t f o r some c e r t a i n g r i d p o i n t s only. 2 2 = 0.0038 m =O.OOll m Q7- 14 1-6

317 '15-23

=0.0036 m

2

2

= 0.0028 m / s e c

Q

( 2 ) measurement n ise c o v a r i a n c e 9 R R= 0.000225 m

U

=

2

2 2 0.0001 m / s e c

( 3 ) f r i c t i o n f a c t o r ( i n i t i a l ) = O . 020 ( 4 ) i n i t i a l e r r o r covariance P(n/n

(5) i n i t i a l conditions MODEL CALIBRATION AND EVALUATION

The s t o c h a s t i c component of t h e model c a n n o t be i n t r o d u c e d u n t i l t h e d e t e r m i n i s t i c component a l o n g w i t h f r i c t i o n f a c t o r a r e d e t e r m i n e d . It is more p r a c t i c a l t o use s e v e r a l t i d a l c y c l e s (35 ) d a t a f o r c a l i b r a t i n g t h e f r i c t i o n f a c t o r i n s t e a d of u s i n g t h e e n t i r e f i l e t o d o t h i s a d j u s t m e n t . The b e s t s o l u t i o n of f r i c t i o n f a c t o r i s 0.021. The c r i t e r i o n o f t h i s c a l i b r a t i o n i s t o s e l e c t t h e minimum sum of s q u a r e e r r o r between t h e o b s e r v a t i o n s and model computed v a l u e s . The f i l t e r e d e s t i m a t e of s u r f a c e e l e v a t i o n s and c u r r e n t v e l o c i t i e s a r e c a l c u l a t e d by t h e f i l t e r i n g p r o c e s s . The comparison from hour 200 t o 245 between t h e observed v a l u e s and f i l t e r e d p r e d i c t i o n a r e shown by F i g u r e 3 ( a ) and 3 ( b ) . The improvementsare made by t h i s c a l c u l a t i o n . The more s i g n i f i c a n t a d v a n t a g e s of t h e f i l tered estimate f o r current v e l o c i t i e s is attributed t o the noisy s y s t e m of t h e v e l o c i t y v a r i a b l e s . An even b e t t e r performance c o u l d be made i f t h e d a t a s e t i s s e l e c t e d by t h e w i n t e r month.

- _ _f i l t e r e d

1.lj

::jLo ------

-1.1

F i g u r e 3 ( a ) Half-Hourly T i d a l Heights (meters)

F i g u r e 3 ( b ) Half-Hourly C u r r e n t V e l o c i t y (m/sec)

The performance of t h e model is e v a l u a t e d by t h e s t a t i s t i c a l methods. The harmonic a n a l y s i s i s used t o examine t h e c o n s i s t e n c y of t i d a l c o n s t i t u t e s between f i e l d measurements and t h e numeric a l s o l u t i o n . The amplitude and phase a n g l e f o r major components of d i u r n a l and s e m i - d i u r n a l f r e q u e n c y bands are c a l c u l a t e d . It i s found t h a t t h e dominant components, such a s M2, S 2 , and N 2 , d o n o t v a r y t o o much (5-10 p e r c e n t ) . The F - t e s t is used t o s e a r c h how t h e s t o c h a s t i c - d e t e r m i n i s t i c approach i s s u p e r i o r t o p u r e l y numeri-

318 c a l s o l u t i o n . Both t i d a l e l e v a t i o n s ( a t Chesapeake C i t y ) and t i d a l c u r r e n t s c a t Summit B r i d g e ) show t h e i r s i g n i f i c a n c e a t 0.01 l e v e l . (Table 1) STATION CHESAPEAKE CITY

SUMMIT BRIDGE

DEGREE OF FREEDOM

58 5 58

STO-DETER MODEL

NUMERICAL MODEL

F-VALUES

0.01589

0.06938

4.36

**

0.004 19

0.000938

2.24

**

Table 1. E s t i m a t e d model p r e d i c t i o n e r r o r c o v a r i a n c e . CONCLUSIONS The one-dimensional hydrodynamic e q u a t i o n s combined w i t h t h e Kalman f i l t e r i s used t o mode1iAg t h e t i d a l h y d r a u l i c b e h a v i o r of t h e C&D c a n a l . The f i l t e r e d e s t i m a t e from t h e computation a l g o r i t h m s f o r t h e s o l u t i o n i s b e t t e r t h a n t h e n u m e r i c a l approx i m a t i o n o n l y . The e x t e n s i o n work involved more sampling s t a t i o n s f o r b o t h s u r f a c e w a t e r e l e v a t i o n s and c u r r e n t v e l o c i t i e s an? c a n i d e n t i f y a b e t t e r p i c t u r e of t h e C&D c a n a l waterway. REFERENCES 1. Budgell W.P.

( 198 I ) , A S t o c h a s t i c - D e t e r m i n i s t i c Model f o r E s t i m a t i n g T i d e s i n Branched E s t u a r i e s , Manuscript Report S e r i e s No. 10. Ocean Science and Surveys C e n t r a l Region, Canada Centre f o r I n l a n d W a t e r s , B u r l i n g t o n , O n t a r i o .

2 . Budgell W . P. ( 1 9 8 2 ) , A Dynamic-Stochastic Approach f o r Modeling Advection-Dispersion P r o c e s s e s i n Open Channel. I n : Time S e r i e s Methods i n H y d r o s c i e n c e s , (Ed. El-Shaarwi and E s b e r t y ) , E l s e v i e r S c i e n t i f i c P u b l i s h i n g Company, Amsterdam, N e t h e r l a n d s , pp. 244-263.

3 . Gardner G . B . and P r i t c h a r d (19741, V e r i f i c a t i o n and Use of a Numerical Model. of t h e C&D c a n a l , C B I , T e c h n i c a l Report 87., Ref 74-7. 4 . H u n t e r , J . R . ( 1 9 7 5 1 , A One-Dimnsional Dynamic and Kinematic Numerical Model S u i t a b l e f o r C a n a l s and E s t u a r i e s , C B I , S p e c i a l Report 4 7 , Ref. 75-10. 5 . Mickle M. H . and Sze T . W . (19721, O p t i m i z a t i o n i n Systems Engineering, I n t e x t Educational Publishers, Scranton, Pennsylvania. 6 . R i v e s , S . K . and P r i t c h a r d D . W . (19781, A d a p t a t i o n of J . R . H u n t e r s ' One-Dimensional Model t o t h e C&D Canal System, C B I , S p e c i a l Report 6 6 , Ref 78-6.