Two-dimensional shallow water flow identification

Two-dimensional shallow water flow identification A. W. Heemink Data Processing Division, Rijkswaterstaat, N~jverheidsstraat I, 2288 BB Rijswijk, The Netherlands (Received June 1987: revised October 1987)

A discrete time-invariant Kalman filter for the identification and prediction of two-dimensional shallow water flow using observations of the water level registered at some locations, has been developed. The filter is based on a set of difference equations derived from the linear two-dimensional shallow water equations using the finite difference scheme proposed by Sielecki. By introducing a system noise process, we can embed the difference equations into a stochastic environment. This enables us to take into account the uncertainties of these equations. A Chandrasekhar-type algorithm is employed to obtain the steady-state filter. In this way the fact that the noise is less spatially variable than the underlying process can be exploited to reduce the computational burden. The capabilities of the filter are illustrated by applying it to the six-hours-ahead prediction of storm surges in the North Sea. The results show excellent filter performance, and, with respect to the results of the underlying deterministic model which were achieved without using the water-level measurements available, the improvement obtained by filtering the measurements is substantial.

Keywords: mathematical models, shallow water equations, Kalman filtering, storm surges, identification, data assimilation Introduction The delta area in the southwestern part of the Netherlands lies below or just above sea level. Many dikes and a storm surge barrier, which has many large gates that will be closed only in case of extreme high water, must guarantee the safety of the people living there. Accurate water-level predictions, approximately six hours ahead, are necessary to decide what precautionary actions have to be taken to protect the dikes. Furthermore, since the storm surge barrier will have to be closed a few hours before the expected occurrence of extreme high water~ accurate predictions are also required to ensure that the barrier will be closed in time. Another reason for the need of accurate predictions of the water level a few hours ahead is that with the increase in size of ships using the Euro Channel and Meuse Channel to Rotterdam harbour, this port has become inaccessible to many ships, except during a short period at high water. To allow safe passage there is a demand for accurate tidal predictions. Previously, tidal prediction methods have been either statistical or deterministic in nature. Applying staffs-

© 1988 Butterworth Publishers

tical methods--i.e., simple empirical or black-box models derived from series of observations--it is possible to use on-line measurements of the water level. This feature is very important, considering the fact that the number of on-line measurement stations in the North Sea is increasing rapidly. Furthermore, statistical methods are very easy to employ and can be implemented on very small computers. Simple linear regression models have been derived, for instance, by Rossiter ~ to predict the meteorological effect at various locations along the eastern coast of England. To predict the water height in inland waters such as Lake St. Clair, Budgell and EI-Shaarawi 2 applied a Box and Jenkins 3 transfer function model. Deterministic methods, including the analytical or numerical solutions to the hydrodynamic equations, do not have the possibility of incorporating on-line measurements. However, they have a more physical basis and provide a more realistic description of the water movement. Linear North Sea models, neglecting the tide and only describing the meteorological effect, were developed by Heaps 4 and Timmerman ~ among

Appl. Math. Modelling, 1988, Vol. 12, April

109

Two-dimensional shallow water flow identification: A. W. Heemink

others; nonlinear tidal models, for instance, were developed by Flather 6 and Voogt. 7 Since the statistical and deterministic approaches have their specific advantages, the use of a method that combines both approaches in order to obtain their features is preferable. Such a method is the Kalman filter. This stochastic filter is based on a deterministic model and also has the capability of correcting model predictions using on-line information. In order to use a Kalman filter for the identification of shallow water flow, the deterministic model describing this water movement is embedded into a stochastic environment by introducing a system noise process. In this way it is possible to take into account the inaccuracies of the deterministic system. By using a Kalman filter, the information provided by the stochastic dynamic model and the noisy measurements taken from the actual system are combined to obtain an optimal (least squares) estimate of the state of the system. Having identified the shallow water flow using the measurements available, one can compute predictions by the underlying deterministic model, where the estimated state of the system is used as the initial condition for the prediction. Since this initial condition has been filtered by using the measurements available, the predictions are improved with respect to the results of the underlying deterministic model, which were achieved without using these measurements. A more detailed discussion on the predictive capability of Kaiman filters for hyperbolic-type systems can be found in Heemink. 8 Kalman filtering can also be considered to be a data assimilation technique. These techniques have been the subject of considerable interest in meteorology (Bengtsson, Ghil, and K/illrng). The most common data assimilation technique used in numerical weather prediction is optimal interpolation. Here some estimates of the error statistics of the numerical model are used to correct the results of the model using the measurements available. Kalman filtering is similar to optimal interpolation. However, with a Kalman filter the error statistics of the numerical model are determined by using the stochastic extension of the model. In this way the correction produced by the Kalman filter is guaranteed to be consistent with this stochastic model, even in case of very irregular flow patterns caused by complicated geometries. The use of optimal interpolation does not guarantee this consistency and, in such cases, will probably yield unrealistic corrections or instabilities. Since the original work of Kalman, ~° Kalman filters have been successfully used in numerous applications. Most of these filters were developed for the determination of satellite orbits and for the navigation of submarines, aircraft, and spaceships. In the last decennium Kalman filter techniques have also gained acceptance in meteorology (Ghil et al.J~), oceanography (MillerJ2), and in several areas of hydraulics and water resources (Chao-lin Chiul3). However, these identification techniques have seldom been applied to tidal prediction problems. Budgell and Unny t4 have

110 Appl. Math. Modelling, 1988, Vol. 12, April

developed a Kalman filter to estimate tides in branched estuaries. This filter is based on the one-dimensional shallow water equations and is used to estimate water levels and velocities at spatially distributed measurement locations as well as between these locations. Heemink and de Jong a5and ten Brummelhuis, de Jong, and Heemink j6 have devised a Kalman filter procedure for the on-line identification and prediction of water levels in an estuary during stormy periods. Their approach is based on a nonlinear one-dimensional model, while the physical parameters in the model, such as bottom friction and wind stress, are conceived as stochastic variables that are estimated together with the water levels and velocities. In this way the model is being adapted continuously to changing conditions. The extension of the one-dimensional time-varying filters just described to two space dimensions does not give rise to conceptual problems, but requires a very good insight into the complex two-dimensional filtering problem. (Here we agree with Polis, ~7 who said that those who develop methods which supposedly work for two-dimensional problems and then treat onedimensional examples to illustrate their approaches, are just kidding themselves.) Furthermore, extending the one-dimensional filters to two dimensions would impose an unacceptably greater computational burden. Therefore, in order to obtain a computationally efficient Kalman filter for two-dimensional shallow water flow identification, in this study a time-invariant filter based on a linear deterministic model is developed. In this case the time-consuming filter equations do not have to be computed over again as new measurements become available, but need only be solved once. As a consequence, the computations can be carried out offline on a large computer. Furthermore, for the computation of a time-invariant filter special algorithms have been developed to reduce drastically the amount of computations. The capabilities of the filter are illustrated in this paper by applying it to the six-hoursahead prediction of storm surges in the North Sea.

Deterministic model The deterministic model is based on the linear shallow water equations, describing the motion of long waves in water where the depth is very little compared to the length of the wave. They are given by the momentum equations Ou --+g Ot

Oh

u

-~x-

fv+h~

- y

V 2 cos - D

+ Ov

+

Oh

+ :"

+

v

-

10pa pw Ox

-0

(1)

V 2 sin q,

1 apa +---= 0 pw Oy

(2)

Two-dimensional shallow water flow identification: A. W. Heemink and the continuity equation Oh --+ dt

O(Du) a(Dv) + =0 Ox Oy

(3)

- (Dij + Dij+ OUik] ~] At 4Ay [ ( D i j + ! + h i + ,

(4)

At an open boundary a radiation-type condition is applied: v± +_ ~ / g / D h = 0

(5)

By means of equation (5) we prescribe the in-going Riemann invariant. For the linear equations (I)-(3), only one condition is required at an open boundary. To obtain a discrete system representation of the model. The equations are discretised. Defining a spacestaggered grid Gt and using the scheme proposed by Sielecki ivan der Houwent9), we obtain the difference equations At , , ~ ; ' = ,~.j - gTA-~xih,.j - h L ,.A At

T V 2 c o s ~b

Di d + Di,j+t

At Vk*l = V ~ j _ - * -- h~j-i) ": ' g 2Av (h,j

_At - : ~

..k+l

*+t

*+~

(ui.j_ ~ + u~ + ~,j_,

k+l)

+ l l i + l , j + lli. j

At < 2 A x A y X / g D ( A x 2 + Ay2)

At dp. pw

(9)

the scheme is stable. Without bottom friction the scheme is not dissipative. If the difference scheme (6)-(8) is embedded in a Kalman filter, it is important that the very short waves that may be introduced by filtering the noisy measurements be eliminated (HeeminkS). Therefore, to improve stability we introduce a smoothing operator. At each time step the velocities are smoothed according to the equations u,.~-=(l-(2+

Ax 2

c ((Ax)

+

*

+

+ \kAy/ {Axe'- k i~ + ~-~y ) u,j+ I -

2+2

(10)

)

~

c vi.j +

\\Ay,I

, (ax'x'- , \ + vi-,.j + v~+,.j + Ik-~Y) vi.J+,l/

,

-

(8)

Here Ax and Ay are the distance steps, At is the time step, and i,j,k are integer indices such that u~j,*v~i,, and h~.j a r e the finite difference approximations of u((2i - 1)bx,2jAy,kAt), v(2iAx,(2j - l)Ay,kAt), and h(2iAx,2jAy,kAt), respectively, and Di.j = D((2i - l)Ax, (2j - 1)Ay). The pressure, wind speed and direction as functions of space and time are computed by a numerical model of the atmosphere. The grid of this model is much coarser than the grid Gj; therefore, interpolation is necessary to obtain the wind speed and direction in the grid points of Gt. The finite difference approximation possesses second-order accuracy in space and first-order accuracy in time. It can be shown that for A ,~ D/At, f ~ l/At, and

v~'j=

+ f - ~ i v , _ , j + v L , ~ + , +v~,.~+ Oi.j+ * I) _ 2Athtt~j

I , j + IlVi,j+ ~, k+ t I

- (Dij + Di+ i.j)V~k.] I]

Since the North Sea is rather deep, the nonlinear effects are small and have been neglected. A numerical model of the North Sea and adjacent waters, based on the equations just described, is used to predict the water levels along the Dutch coast (TimmermanS). Because of the linearity of the model, in practice the complex astronomical tide is predicted separately by harmonic analysis (Godin~a). The numerical model is only used to describe the meteorological effects which are superimposed on the atmospherical tide. The wave motion is completely described by equations (1)-(3), provided that initial values and boundary conditions are given. At a closed boundary the velocity normal to the boundary is zero: v± = 0

(7)

At 4Ax [(Di+ t ,j + D i + l,j+ IJtli+ ,. k+I.j i

hk+!i,j = h~.j -

where h = water level //,/3 = water velocities in the x- and ydirections, respectively D = depth of the water f = Coriolis parameter A = linear bottom friction coefficient 7 = wind friction coefficient V = wind speed q,= wind direction Pw = density of water P a = atmospheric pressure g = acceleration of gravity

At C]pa p,~ dy

- 2At Av~j - TV 2 sin $ Di.j + Di+ t.~

(6)

v~.j_.

ill)

where u~. and v~. are the smoothed values of u/*.j and v~j, respectively, and c < ½Ay2/(Ax: + Ay2) is a positive constant that is chosen as small as possible to guarantee that the finite difference scheme is still an accurate approximation of the equations (1)-(3). By defining equations il0) and (I I), we introduce numerical viscosity. It is easy to show that the smoothing operators (10) and i l l ) are the second-order approxi-

Appl. Math. Modelling, 1988, Vol. 12, April

111

T w o - d i m e n s i o n a l s h a l l o w w a t e r f l o w identification: A. IN. H e e m i n k mations of the diffusion terms

At =

4cAx2fO2u 02u'~ At \ Ox2 + Oy2]

v .j -

(12) -

.f

g

(h

j -

At (u/.j_ k+t l +

,,k +,,.j-, + -i+ . k + ,.j J + u~.f ') -i+

and _ 2At A v ~ j - yV2 sin ~b + Wv k+l

4cAx 2 (02v 02"~ At \Ox 2 + Oy2]

respectively. The very short waves which cannot be represented by the finite difference approximations (l)(3)Jn any case are now dissipated. At a closed boundary where the velocity normal to the boundary is zero, no special boundary scheme is required. The discretised open boundary condition is /-~/+li.j + ~ g / D h~.f ~ = 0 -

AX,~÷, = BX,~ + u,k+ ,

(15)

where A and B are coefficient matrices, the state X,, is an n-vector consisting of the water levels and velocities in all the grid points at time kAt, and u, k represents the meteorological input at time kAt.

Stochastic d y n a m i c system representation The deterministic model (15) is not perfect. Errors are introduced by fluctuations in the meteorological input (or, which amounts to the same thing, by the poorly known influence of the wind on the tidal movement described by the wind stress coefficient) and by the open boundary conditions which have to be specified. Furthermore, errors associated with the neglected nonlinearities, the specification of the uncertain bottom friction coefficient, and the numerical approximation also contribute to considerable uncertainty. Therefore, we embed the finite difference equations (6)-(8) into a stochastic environment by adding white Gaussian system noise. Assuming that the most important errors of the finite difference equation are caused by an erroneous wind input, the system noise processes Wu ~ and Wvk are introduced in, respectively, equations (6) and (7) as follows:

.t:j

-

g

At

k

(h,.j

-

At vk

+ fT(

,-,.j + vL,.j+, + v ij + vtj+,)

- 2At Au~j - yV 2 cos 0 + Wu k + Dij + D,-j+ i

At Op, pw Ox

The processes W t : and Wv k are assumed to be location invariant with statistics

E{Wuk(x,y)} = 0

(16)

112 Appl. Math. Modelling, 1988, Vol. 12, April

(18)

E{ Wuk(x l ,y l) Wuk(x2,Y2) } = 0"2exp [--aX/(xl -- x2)2 + (Yl -- Y2)2]

(14)

-

Note that (14) is a zero-order approximation. In addition it is assumed that along the boundary the velocity parallel with the boundary is zero. This assumption is necessary only because of the discretised Coriolis terms. The finite difference equations (1)-(3), the smoothing operators (10) and (11), and the boundary condition (14) can be rewritten as a discrete system

=

Di.j + Di+lu

(13)

AtOp, pw Ox (17)

E{W.k(x,y)} = 0

(19) (20)

E{ Wvk(x l,y O W,k(x2,y_,)} = cr2exp[--aX/(xl --x2) 2 + (Yl --Y2) 2]

E{Wuk(x~,yOWvk(x2,y2)} = 0

(21)' (22)

The parameters o~ and a have to be specified. Here we note that the covariances are parameterized rather simply. H o w e v e r , as long as our knowledge of the system noise processes remains very poor, it is not useful to improve these parametrizations. The noise processes Wu k and Wv ~ are introduced to model the uncertainty associated with the momentum equations. The continuity equation is assumed to be perfect. The parameter a in equations (19) and (21) can be interpreted as a measure of the spatial variability o f the errors concerned with the momentum equations. In choosing a suitable value for a, it has to be taken into account that finite difference schemes are not able to represent short waves accurately (i.e., waves with a wavelength of the order of 2Ax). Therefore, in order to obtain meaningful solutions and to avoid instabilities, the energy of these noise waves should be limited. As a consequence, a has to be chosen sufficiently small. Here we note that for small values of a, the smoothing factor c, introduced to eliminate instabilities, can be chosen small too. The uncertainty associated with the open boundary condition (14) can be introduced according to h ~.j

-

•

•

(23)

where the white Gaussian system noise W~j represents the uncertain part of the boundary. W~ij.is assumed to be independent of the system noise processes W t : and Wv k. By means of equation (23) we consider the ingoing Riemann invariant to be corrupted by noise. Of course, there are other ways of introducing uncertainty to an open boundary. A more general stochastic open-boundary treatment can be found in Heemink. 8 Using measurements taken along the boundary, we can estimate h~j by the filter. If there are no observations available in the neighborhood of the boundary, the filter is not able to estimate these components of the state vector accurately, and the erroneous bound-

Two-dimensional shallow water flow identification: A. W. Heemink ary condition propagates into the model. However, as soon as measurements are available, the errors caused by the erroneous boundary condition that has been prescribed can be corrected by filtering these measurements. Note that although the ingoing Riemann invariant cannot be measured in nature, it can be identified by the filter using water-level observations. If the system noise is defined on a grid G2, which need not be space-staggered, the stochastic dynamic model can be described by AX,,+, = BX,, + u,,+, + AW,,~,

X,o = Xo

(24)

where W,, is a p-vector consisting of the noise components at the boundary and at the grid points of G2. The covariance matrix Q of W,, can be derived easily from equations (18)-(23). The initial condition Xo is also assumed to be Gaussian with statistics:

E{Xo}=

ko

E{[Xo - Xo] [Xo - X0]T} = Po

(25)

The matrix A represents a sequence of linear operations required to interpolate the system noise at the grid points of G,. Since the most important model errors are caused by fluctuations in the meteorological input, G2 is chosen to coincide approximately with the grid of the atmospheric model, yielding p <~ n. An important advantage of the fact that the grid Gz of the system noise is coarser than the grid GE is that the energy of the short waves introduced on Gj by the filter is limited. If, in addition, the parameter a is small and the bottom friction is sufficiently strong, it may not be necessary to introduce numerical viscosity to eliminate the short waves. Assuming that measurements of the water level are available at m grid points of G~, the observation equation can be derived easily: Z,~ = MX,, + V,,

(26)

where the m-vector Z,, contains the measurements taken at time kAt, and V,, is the white Gaussian measurement noise at time k A t with covariance matrix R. In practice, we may assume that the measurement errors at different locations are mutually independent and have equal variance, so that R is diagonal with elements r'-. If necessary, equation (25) can easily be modified to account for the fact that not all measurements are available exactly at grid points. It is further assumed that the initial state Xo, the system noise W,,, and the measurement noise V,, are mutually independent. The approach described above has been used to develop a filter based on a numerical model describing the water movement in the entire North Sea and the Channel. In F i g u r e 1 the grid G~ covering this area is shown. Here Ax = 18.5 km and Ay = 19 kin, and as a consequence of condition (9) At = 10 min. Furthermore, the parameters in the model are chosen to be the linear bottom friction coefficient h = 2.4 x 10 -3, the Coriolis parameter f = 1.25 x 10 -4, and the wind friction coefficient 3' = 3.2 x 10 -6. The grid G_, and

Figure I

The model of the North Sea

the measurement locations in the North Sea that have been used in this study are also shown in F i g u r e 1.

Kalman filtering It is desired to combine the measurements taken from the actual system and modelled by relation (26) with the information provided by the system model (24) in order to obtain an optimal estimate of the system state X,,. If X(kll) is defined as the least squares estimate of X,~ based on the measurements Z,,, Z, 2. . . . . Z,,, and P(k[I) represents the covariance matrix of the estimation error, recursive filter equations to obtain these quantities can be summarized as follows. The optimal state estimate is propagated from measurement time tk-~ to measurement time tk by the equations A X ( k ] k - 1) = B X ( k - l [ k AP(k[k-

I)a v =BP(k-

- 1) + u,~

l [ k - I)B T + A Q A T

(27) (28)

At measurement time tk, the measurement Z,, becomes available. The estimate is updated by the equations "X(klk) = ~((kJk - I) + K(k)[Z,, - A,ffX(klk - 1)] P(k[k) = P(k[k - 1) - K ( k ) m P ( k l k

- 1)

(29) (30)

where K(k) = P(k[k - I ) M T [ M P ( k l k

- 1)M T + R ] -

is the filter gain. The initial condition for the recursion

Appl. Math. Modelling, 1988, Vol. 12, April

113

Two-dimensional shallow water f l o w identification: A. W. Heemink

is given by X:(010) = Xo

e(010) = Po

(31)

The filter just described is the celebrated Kalman filter for a discrete problem formulation. The derivations of these filter equations can be found, for instance, in Maybeck's work. 2° The filter has a predictor-corrector structure. Based on all previous information, a prediction of the state vector at time tk is made by means of equations (27) and (28). Once this prediction is known, it is possible to predict the next measurement by means of equation (26). When this measurement has become available, the difference between this measurement and its predicted value is used to update the prediction of the state vector by means of equations (29)-(30). The model (24)-(26) is not the most general one. Using the Kalman filter algorithm (27)-(31),we can allow all the coefficient matrices A, B, A, and M as well as the covariances Q and R to be time varying. This is a very favourable property of this algorithm. However, using it for the identification of two-dimensional shallow water flow would impose an unacceptable computational burden. In order to obtain a computationally efficient filter, simplifications have to be introduced. Unfortunately, there are serious problems with the more obvious simplifications one may consider. One possibility is to divide the domain of interest into subregions and to apply the filter algorithm (27)(31) to each subregion individually. The main difficulty associated with this approach is that only the measurement stations located inside a subregion are used by the filter to correct the estimates in this subregion. The information available in other subregions cannot be incorporated to improve these estimates. As a consequence, strong variations may occur between the components corresponding to subregion boundaries and those corresponding to grid points immediately outside. Hence, at the boundaries of the subregions, artificial, high gradients are likely to be introduced. Although this Kalman filter approach has been successfully applied to predict diffusion processes (Fronza, Spirito, and Tonie]li2~), we may expect serious difficulties if it is employed to predict the tidal movement. Since the system of the shallow water equations describing this movement is hyperbolic, the distortions introduced by the partial filtering would propagate through the entire domain of the problem, probably causing unsatisfactory filter performance or instabilities. Another approach would be to calculate the covariance matrix on a coarser grid than that used for the computation of the estimates. As a consequence, the correction achieved by the filter is also available on the coarser grid and has to be extended to the finer grid to obtain the estimates. In case of a complicated boundary, this procedure is very complex and prone to numerical difficulties. Moreover, a covariance matrix calculated on the basis of a coarser mesh would be the covariance matrix for a model with different geometry. But as Miller ~2pointed out, it would be the right solution to the wrong problem.

114

Appl. Math. Modelling, 1988, Vol. 12, April

Another possibility for simplification would be to make an assumption about the structure of the covariance matrix. If it is assumed that errors at large distant points are not correlated, the covariance matrix is banded down its diagonal. The sparseness of the covariance matrix can easily be exploited to reduce the computation time and the storage requirements of the algorithm (27)-(31). In our case, however, the assumption just mentioned is not realistic because the errors are highly correlated in space. Due to the difficulties described above, the filter that is developed in this study has been based on the time-invariant deterministic model (15). The main advantage of using this linear model is that if the noise statistics are time invariant the filter gain becomes time invariant as well. Here we note that this assumption concerning the noise statistics is not very restrictive because very little is known (yet) about the time-varying behaviour of the noise processes, and therefore these processes are usually assumed to be stationary. In case the filter gain is time invariant the timeconsuming second moment calculations need only be solved once. These computations can be performed off-line on a large computer. Furthermore, since the system noise has been defined on a coarser grid than that used for the computation of the estimates, high dimensionality of the filter equations can be avoided by using a discrete form of the Chandrasekhar-type algorithm. This algorithm was first proposed by Morf, Sidhu, and Kailath 22 and has been used in numerous applications. It uses the fact that for certain initial conditions the incremental covariance has a rank p .(p is the dimension of the system noise process) and can be factorized as follows: P(klk) - P(k - Ilk - 1) = S(k)L(k)S(k) T

(32)

where P(klk) is the n x n covariance matrix of the state estimate at time k based on the observations available up to and including time k (n is the dimension of the system state), S(k) is an n x p matrix, and L(k) is a p x p matrix. For the model (24)-(26) the recursive equations to obtain the steady-state filter gain are A Y ( k + 1) = BS(k)

(33)

G(k + 1) = G(k) + Y(k + 1) x L(k) Y(k + 1)TMT

(34)

R'(k + l) = R'(k) + M Y ( k + 1) × L(k)Y(k + i)TM T

(35)

K(k + 1) = G(k + 1)R'(k + l ) - t

(36)

S(k + 1)= Y(k + l ) -

(37)

K(k + I)MY(k + l)

L(k + 1) = L(k) + L(k)Y(k + 1)T x M T R ' ( k ) - t M Y ( k + 1)L(k)

(38)

where K(k) is the filter gain at time kAt and the initial condition for the recursion is given by AY(I)=A R~(0)=R

G(0)=0 L(0)= Q

(39)

Two-dimensional shallow water flow identification: A. W. Heemink Equations (33)-(38) are iterated until IIg(k + 1) - K(k)ll < ellK(k)ll

(40)

where e is prespecified. Since the underlying deterministic model (15) is of the hyperbolic type, the number of iterations depends on the travelling times of the waves in the model and, therefore, on the size of the domain of the problem. In some applications the Chandrasekhar-type algorithms were shown to be prone to numerical difficulties (Van Dooren and Verhaegen23). A detailed discussion on the numerical properties of the algorithm (33)-(38) can be found in Heemink. 8 It is shown that numerical difficulties can be avoided by properly choosing the smoothing factor and the spatial variability of the system noise. If K is the desired steady-state filter gain, it is easy to verify that X(/l/) can be calculated recursively by

X(I + lit + I) = [ 1 - KM]A-'B'JI(III) + gz,,

(41)

Here A - t is neither computed not stored; it serves to represent a sequence of linear operations. The filter algorithm is implemented on a CDC Cyber 205 vector processor. The CPU time required for the computation of the steady-state filter based on the model described in this paper consisting of more than 900 active calculation points is approximately 12 s. From these results it is to be expected that when using a vector processor it is possible to combine the steadystate filter approach with very large deterministic models.

Filtering results Before the filter can be safely applied to prototype situations, it is necessary to demonstrate that it can perform adequately under known conditions. In this case by means of a random generator, data sets are created using a " t r u t h m o d e l . " By employing a Kalman filter, based on a "filter m o d e l " which differs from the truth model, it is possible to investigate whether the filter is able to reconstruct the water levels and velocities using the simulated data, despite the differences between the filter model and the truth model. The importance of these experiments is that the true water levels and velocities are known, and, consequently, the performance of the filter can be evaluated quantitatively under a variety of circumstances. It is obvious that when using such a limited number of measurements as shown in Figure 1, the entire wave motion cannot be reconstructed uniquely (i.e., the system is not observable). The spatial variability of the water movement that can be identified depends on the spatial distribution of the measurement locations. In the southern part of the North Sea the number of tide gauges available is relatively large. As a consequence, the filter is able to produce a detailed picture of the water movement. H o w e v e r , in the remote parts of the model this is not the case. Fortunately, owing to the Coriolis forces, external surges created in the northern part of the model propagate approximately like Kelvin

waves along the English coast toward the south (Dronkers24). We may expect that by using the measurements available as shown in Figure 1, this wave motion, which is in essence one-dimensional, can be identified accurately. Some experiments using simulated data can be found in Heemink 25 and, in more detail, in Heemink. 8 The results have shown excellent filter performance and that the filter is able to correct an erroneous meteorological input or errors introduced by the open boundary condition. In this paper we describe an application of the Kalman filter to field data gathered during the stormy period January 26, 1983 to February 2, 1983. The wind velocities and pressure fields during this period have been computed by an atmospheric model of the K.N.M.I. (Royal Netherlands ~Meteorological Institute). The measurements were available each hour and are assumed to have been taken at the nearest grid point of the model. Because of the linearity of the model, the astronomical tide has been eliminated by harmonic analysis. The standard deviation of the measurement noise is chosen to be r = 0.03 m. The system noise is introduced at each measurement. The values of the system noise parameters tr 2 = 100 m2/s 2 and a = 0.90 are determined from the data available from January 26 to January 29, 1983 by a combination of physical intuition and trial and error; i.e., one employs the filter for various values of these parameters until one gets satisfactory filter performance. The performance of the filter can be judged by monitoring the residuals R, k, defined as the difference between the measurements and the prediction of these measurements based on all previous information: R,, = Z,, - MX(k[k - I)

(42)

It is easy to verify that

E{R,,} = 0

(43)

E{R,,R T} = MP(kIk - I)M r + R Since the theoretical statistics of the residuals R,, are known, the actual residuals can be monitored and compared with this description. By checking whether the residuals indeed possess their theoretical statistical properties we are able to judge whether the mathematical model satisfactorily describes the real system behaviour. In fact, this is the calibration of the filter. Fortunately, the filter results turned out not to be very sensitive to the choice of the noise parameters. The poorly known noise statistics is one of the major problems of the application of Kalman filters, and therefore adaptive filters often have been employed (Jazwinski26). Using these filters, we examine the residuals to establishwhether they actually possess their theoretical statistical properties, and, if not, we adopt the statistics of the noise processes. However, in our case these filters are not successful. Since the filter results are not very sensitive to the noise statistics, many residuals have to be analysed before a statistically significant conclusion can be drawn. Unfortunately, during stormy periods the conditions change

Appl. Math. Modelling, 1988, Vol. 12, April

115

Two-dimensional shallow water f l o w identification: A. W. Heemink

rapidly, so the filter is not able to adapt the model in time. In fact, the storm is over before the filter is aware of the fact that there is something wrong with the covariance matrices of the noise processes. The results of the three-hours-ahead predictions of the meteorological effect at a number of measurement stations during the period January 26, 00:00 to January 29, 23:00 are shown in Figure 2. Here the predictions of the deterministic model without using the data and the filter predictions are compared with the observations. During this period two relatively small surges occurred. The six-hours-ahead predictions at IJmuiden can be found in Figure 3. Here, since we are mainly interested in the predictions during high and low water, we also present the predictions with the astronomical tide included. Similar results obtained during January 30, 00:00 to February 2, 23:00, are shown in Figures 4-5. During this period one small surge and one extremely large surge occurred. The results show that, although the deterministic model produces erroneous predictions, the filter is able to predict relatively accurately. Here we note that by

..... OETEI~INISTIC P ~ : O I C T I ~ --FILTF~PRf:OICTI~

HOUR

-2. nEASUL~'f:NT DET~INISTIC P~01CT 10N

..... - -

FILTER E~EOICTION

-3

,.........

i~ . . . . . . . .

~ ........

~ . . . . . . . . ;9 . . . . . . . . ~ . . . . . . . . ~ . . . . . . . . ~e . . . . . . . . ~i . . . . . . . . ~ . . . . . . . . ;89

HOUR

Figure 3 January

"'---. IO

29

o,," 38

46

50

G9

70

89

90

189

~ ......

189

HOUR

rEkSLIREtIENT OETE~IJNISTIC PRf:OICTION FILTF~ PREDICTION

..... - -

" ' - . . . . . o .o" ........

ie. . . . . . . .

i~ . . . . . . . . ~i . . . . . . . . ~g . . . . . . . . ~i . . . . . . . . ~i . . . . . . . . ~fi . . . . . . . . ~ . . . . . .

HOtlt

• Z

.....

I~TF3allNIST IC PREDICTION

- -

FILTER PREDICTION

1

-I

. . . . . . . . . ; i . . . . . . . . ~ . . . . . . . . ~ . . . . . . . . ; i . . . . . . . . ~ . . . . . . . . ~ . . . . . . . . ~i . . . . . . . . ~ . . . . . . . . ~ . . . . . . . . 189 HOUR

Figure 2 26,

00:001:

T h r e e - h o u r s - a h e a d predictions (0 h o u r s = J a n u a r y (al.Lowestoft;

(b).

Southend;

(c).

IJmuiden

116 Appl. Math. Modelling, 1988, Vol. 12, April

Six-hours-ahead 26,

predictions

at

IJmuiden

(0

hours

=

00:00)

increasing the domain of the North Sea model or by using a more detailed, possibly nonlinear model, the deterministic predictions can be improved. However, the largest errors of the deterministic model are caused by the erroneous meteorological input; therefore no such model, linear or nonlinear, would produce accurate results. (This is no criticism of atmospheric models. Predicting meteorological conditions is extremely complicated owing to the chaotic behaviour of the atmosphere). Of course, the filter results also depend on the meteorological input. However, the influence of an erroneous input on the predictions produced by the filter is limited. Since the model is of the hyperbolic type, only input errors inside the domain of dependence of the predictions degrade the performance of the filter (see Figure 6). The size of this domain depends on the prediction interval. By increasing this interval the influence of an erroneous meteorological input becomes more important. Besides the erroneous meteorological input, prediction errors are caused by the fact that we have used a linear model and have neglected the nonlinear interaction between tide and surge. However, as in the case of the erroneous meteorological input, only the nonlinearities inside the domain of dependence of the predictions degrade the filter performance. Unfortunately, the most important nonlinear effects are introduced in the southern part of the North Sea and, consequently, for the six-hours predictions, are inside the domain of dependence of these predictions. Therefore, we may

Two-dimensional shallow water flow identification: A. W. Heemink L

• ..... - -

nE~SU~dlBa Pff:DICTI ON FILTE~P~DICTI~ ~'tlNISTIC

2

• ." ",..

......

0ET~IINISTIC I ~ 0 I C T I ~

.

~'~

•

i~.

i! •

""

I

°"

......... ii ........ ~i ....... ~ ........ ;i ........ ~ ........ ~ ........ ~i. . . . . ~ ........ ~ ........ ;~

•

"-t

U

t

-L

.....

0EI'E.~IHISTIC P~EI)I L'T10N

.

.-.

........

, .........

I0

~

8

Figure

5

.

.........

2tl HOUR

. ..........

30

Six-hours-ahead

. ........

~

. ........

, ...................

50

G8

predictions

70

at I J m u i d e n

, .........

, .........

88

9O

(0 h o u r s

,

=

January 30, 00:00) -I

. . . . . . . . . {e . . . . . . . . ~ . . . . . . . . ~ . . . . . . . . ;g . . . . . . . . ~ . . . . . . . . ~ . . . . . . . . ~ . . . . . . . . ~a. . . . . . . . ~

........ ioe

~

,our 4 Three-hours-ahead p r e d i c t i o n s (0 h o u r s = J a n u a r y 30, 0 0 : 0 0 ) (a). L o w e s t o f t ; (b). S o u t h e n d ; (c). I J m u i d e n

CHARACTERISTIC LINES

Figure

PREOICTION

t~

~

expect that the filter predictions can be improved by the use of a nonlinear model.

PRESENT

T

j

IPREDICTIO N

~

__so__.T

Conclusions and recommendations In this paper a steady-state Kalman filter based on a linear two-dimensional model of the North Sea and the Channel has been described. The filter has been tested extensively using field data. The results show excellent filter performance, especially if we take into account that the number of measurements available (as yet) has been limited. With respect to the results of the underlying deterministic model achieved without using the water-level measurements available, the improvement obtained by filtering these measurements is substantial. Since the deterministic North Sea model that has been developed in this paper is linear, the nonlinear interaction between tide and surge has not been taken into account. Therefore, to improve the predictions it is intended to extend the filter to the weakly nonlinear case and to combine it with the nonlinear model of the North Sea and adjoining waters that has recently been developed by Voogt. 7

STATION

D o m a i n o f d e p e n d e n c e o f t h e p r e d i c t i o n s in c a s e o f o n e s p a c e d i m e n s i o n (in c a s e o f t w o s p a c e d i m e n s i o n s w e h a v e , i n s t e a d o f t w o c h a r a c t e r i s t i c lines, a c h a r a c t e r i s t i c c o n e ) . Figure

6

In the derivation of the filter equations it has been assumed that the system noise is white (i.e., uncorrelated in time). In practice, this condition is not likely to be satisfied. Jazwinski 26 describes a procedure to take into account an exponential autocorrelation of the system noise by augmenting the state vector with the noise components. In our case the increase of the computational effort would be very limited, and we may expect that by implementing this feature the accuracy of the predictions can be further improved.

Acknowledgment The author wishes to acknowledge stimulating discussions with Dr. B. de Jong and Professor P. J. Zandbergen of the University of Twente.

Appl. Math. Modelling, 1988, Vol. 12, April

117

Two-dimensional shallow water f l o w identification: A. W. Heemink

References 1 2 3 4 5 6 7 8 9 10 11

12 13 14

Rossiter, J. R. Research on methods of forecasting storm surges on the east and south coasts of Great Britain. Quart. J. Roy. Met. Soc. 1959, 85, 262 Budgell, W. P. and EI-Shaarawi, A. Time series modelling of storm surges on a medium-sized lake. In Marine Forecasting, ed. J. C. J. Nihoul. Elsevier, Amsterdam, 1978, p. 197 Box, G. E. P. and Jenkins, J. M. Time Series Analysis: Forecasting and Control. Holden-Day, San Francisco, 1970 Heaps, N. S. A two-dimensional numerical sea model. Philos. Trans. A 1969, 265, 93 Timmerman, H. On the importance of atmospheric pressure gradients for the generation of external surges in the North Sea. Deutsche Hydro. Z. 1975, 28, 62 Flather, R. A. Tidal model of the north-west European continental shelf. M#m. Soc. R. Sci. Li#ge 1976, 10, 141 Voogt, L. Een Getijmodel van de Noordzee gebaseerd op de Jonsdap-1976 Meeting, rapport WWKZ-84G.006, 1985 (in Dutch) Heemink, A. W. Storm surge prediction using Kalman filtering, Vol. 46. Rijkswaterstaat Communications, The Hague, 1986 Bengtsson, L., Ghil, M. and K~illrn, E. eds. Dynamic Meteorology: Data Assimilation Methods. Springer-Verlag, New York, 1981 Kalman, R. E. A new approach to linear filter and prediction theory, J. Bas. Engng. 1960, 82, 17 Ghil, M., Cohn, S. E., Tavantzis, J., Bube, K., and lsaacson, E. Application of estimation theory to numerical weather prediction. In Dynamic Meteorology: Data Assimilation Methods, eds. L. Bengtsson, M. Ghil, and E. K~ill~n. Springer-Verlag, New York, 1981, p. 139 Miller, R. N. Toward the application of the Kalman-Bucy filter to regional open ocean modelling. J. Physical Oceanography, 1986, 16, 72 Chao-lin Chiu (ed.) Application of Kalman filter theory to hydrology, hydraulics and water resources. University of Pittsburg, Pittsburg, 1978 Budgell, W. P. and Unny, T. E. A stochastic-deterministic model for predicting tides in branched estuaries. Proceedings

118 Appl. Math. Modelling, 1988, Vol. 12, April

15

16

17

18 19 20 21 22 23 24 25

26

of the Third International Symposium on Stochastic Hydraulics, Tokyo, 1980 Heemink, A. W. and de Jong, B. The use of Kalman-Bucy filters in forecasting the water-levels in the Dutch coastal area. Proceedings of the Fourth International Conference on Finite Elements in Water Resources. Springer-Verlag, Berlin, 1982 ten Brummelhuis, P. G. J., de Jong, B. and Heemink, A. W. On-line prediction of water-levels in an estuary using Kalman filters. Proceedings of the Fourth International Symposium on Stochastic Hydraulics, Urbana-Champaign, 1984 Polis, M. P. The distributed system parameter identification problem: A survey of recent results. Proceeding of the Third IFAC Symposium on the Control of Distributed Parameter Systems, Toulouse, 1982 Godin, G. The Analysis of Tides. University of Toronto Press, Toronto, 1972 van der Houwen, P. J. Finite Difference Methods for Soh,ing Partial Differential Equations. Mat. Cen. Tracts, vol. 20, Mat. Cen., Amsterdam, 1968 Maybeck, P. S. Stochastic Models, Estimation and Control, Vol. I. Academic Press, New York, 1979 Fronza, G., Siprito, A. and Tonielli, A. Real-time forecast of air pollution of episodes in the Venetian region. Part 2. The Kalman predictor. Appl. Math. Modelling 1979, 3, 409 Morf, M., Sidhu, S. S., and Kailath, T. Some new algorithms for recursive estimation in constant, linear, discrete time systems. IEEE Trans. Autom. Control 1974, AC-19, 315 van Dooren, P. and Verheagen, M. On the reliability and efficiency of different Kalman filter algorithms. Part I. Theoretical results. IEEE Trans. Autom. Control 1986, AC-31,907 Dronkers, J. J. Tidal Computations in Rivers and Coastal Waters. North-Holland, Amsterdam, 1964 Heemink, A. W. Application of Kalman filtering to tidal flow prediction. Proceedings of the IFAC/IFORS Symposium on System Analysis Applied to Water and Related Land Resources, Lisbon, 1985 Jazwinski, A. H. Stochastic Processes and Filtering Theory. Academic Press, New York, 1970