- Email: [email protected]
Dispersion and stability analyses of the linearized two-dimensional shallow water equations in Cartesian coordinates
Ocean Engineering 30 (2003) 2343–2362 www.elsevier.com/locate/oceaneng
Dispersion and stability analyses of the linearized two-dimensional shallow water equations in Cartesian coordinates S. Sankaranarayanan a,, Malcolm L. Spaulding b a
Applied Science Associates, 70 Dean Knauss Drive, Narragansett, RI 02882, USA Ocean Engineering, University of Rhode Island, Narragansett, RI 02882, USA
b
Received 18 November 2002; accepted 16 April 2003
Abstract In the present study, a Fourier analysis is used to develop expressions for phase and group speeds for both continuous and discretized, linearized two-dimensional shallow water equations, in Cartesian coordinates. The phase and group speeds of the discrete equations, discretized using a three-point scheme of second order, five-point scheme of fourth order and a three-point compact scheme of fourth order in an Arakawa C grid, are calculated and compared with the corresponding values obtained for the continuous system. The threepoint second-order scheme is found to be non-dispersive with grid resolutions greater than 30 grids per wavelength, while both the fourth-order schemes are non-dispersive with grid resolutions greater than six grids per wavelength. A von Neumann stability analysis of the two- and three-time-level temporal schemes showed that both schemes are stable. A wave deformation analysis of the two-time-level Crank–Nicolson scheme for one-dimensional and two-dimensional systems of shallow water equations shows that the scheme is nondispersive, independent of the Courant number and grid resolution used. The phase error or the dispersion of the scheme decreases with a decrease in the time step or an increase in grid resolution. # 2003 Elsevier Ltd. All rights reserved. Keywords: Stability; Dispersion; Cartesian coordinates; Compact difference
Corresponding author. Formerly at Ocean Engineering, University of Rhode Island, Narragansett, RI 02882, USA. Tel.: +1-401-7896224; fax: +1-401-7891932. E-mail address: [email protected] (S. Sankaranarayanan).
0029-8018/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0029-8018(03)00103-3
2344
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
1. Introduction The dispersion and stability analyses for various schemes used in solving shallow water equations have been carried out in the past using a Fourier analysis of the differencing techniques, while solving the governing equations in a Cartesian coordinate system (Leendertse, 1967; Messinger and Arakawa, 1976). Foreman (1984) estimated the accuracy of some finite difference and finite element techniques, by comparing the numerical and analytical plane wave solutions using a Fourier analysis. Song and Tang (1993) analyzed the dispersion and group velocity for linearized three-dimensional hydrodynamic equations discretized using a hybrid scheme on Arakawa’s B and C spatial grid types. It is noted that Arakawa and Lamb (1977) and Song and Tang (1993) evaluated the dispersion from a spatially discretized system, while retaining the time derivative in its differential form; Leenderstse (1967) and Foreman (1984), however, evaluated the numerical dispersion for fully discretized equations, from their eigenvalues obtained through a Fourier analysis. Messinger and Arakawa (1976) found that fourth-order numerical schemes have lower numerical dispersion when compared with traditional second-order schemes. Hirsch (1975) used a three-point compact scheme of fourth order, whereas traditional fourth-order methods require a five-point stencil. In compact difference methods, the function and all of its derivatives are considered as unknowns for obtaining higher-order accuracy with the minimum number of nodes, whereas traditional finite difference methods require additional nodes to achieve higher-order accuracy. Compact difference methods thus avoid the difficulties of having to consider additional fictitious nodes, near the boundaries, to maintain higher-order accuracy. Hirsch (1975) also noted that coefficients of the truncation error term are smaller in the compact scheme, for the same differencing order. Later, Lele (1992), while solving the one-dimensional advection equation, found that the three-point compact scheme has marginally better dispersion properties than the traditional five-point scheme. The three-point second-order scheme is employed in most coastal and ocean models. The five-point fourth-order scheme (Messinger and Arakawa, 1976) and the three-point fourth-order compact scheme (Lele, 1992) were shown to possess less dispersion when solving the advection equation. The three-time-level (Muin and Spaulding, 1996) and two-time-level Crank–Nicolson scheme (Casulli and Cattani, 1994) are second-order accurate in time and space and have been successfully used in coastal modeling. In the present study, the three spatial discretization schemes and two temporal discretization schemes discussed above are considered for stability and dispersion analysis. One-dimensional and two-dimensional wave deformation analyses, based on a complex propagation factor developed by Leendertse (1967), are used to estimate the amplitude and phase errors of the fully discretized equations.
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2345
2. Governing equations The linearized, depth-averaged, two-dimensional shallow water equations, in Cartesian coordinates are given as: Continuity equation: @f @UD @VD þ þ ¼0 @t @x @y Momentum equation in the x-direction: @UD @f ¼ gD @t @x Momentum equation in the y-direction: @VD @f ¼ gD @t @y
ð1Þ
ð2Þ
ð3Þ
where U(x, y, t) and V(x, y, t) are the vertically averaged velocities in the x and y directions, respectively, fðx; y; tÞ is the water surface elevation, g is the acceleration due to gravity, and D is the total water depth. 3. Dispersion analysis of the continuous system Taking the case of constant depth, Eqs. (1)–(3) become: Continuity equation: @f @U @V þD þD ¼0 @t @x @y Momentum equation in the x-direction: @U @f þg ¼0 @t @x Momentum equation in the y-direction: @V @f þg ¼0 @t @y
ð4Þ
ð5Þ
ð6Þ
A continuous Fourier mode is introduced in each field variable, assuming that the governing equations are defined on an infinite spatial domain or with periodic boundary conditions on a finite domain, and represented by U ¼ U0 ei½kx xþky yxt V ¼ V0 ei½kx xþky yxt f ¼ f0 ei½kx xþky yxt
ð7Þ
2346
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
Eqs. (4)–(6) reduce to 0 10 1 ix 0 igkx U0 @0 ix igky A@ V0 A ¼ 0 iDkx iDky ix f0 where kx and ky are the wave numbers in the x and y direction, respectively. A non-trivial solution of Eq. (8) requires that x3 xgD kx2 þ ky2 ¼ 0 Solution of Eq. (9) gives three frequencies: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 x ¼ 0; gD ky þ kx
ð8Þ
ð9Þ
ð10Þ
The x-directed group velocity CGx can be obtained by taking a partial derivative of Eq. (10) with respect to the wave number in the x-direction (kx) as 9 8 > > > > = < dx gDkx ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ð11Þ CGx ¼ ¼ 0; > dkx > > ; : gD ky2 þ kx2 > The y-directed group velocity CGy can be obtained by taking a partial derivative of Eq. (10) with respect to the wave number in the y-direction (ky) as 9 8 > > > > = gDky dx < ð12Þ CGy ¼ ¼ 0; rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > dky > > ; : gD ky2 þ kx2 > 3.1. Dispersion analysis for a second-order finite difference scheme In this section, expressions for phase and group speeds are developed for linearized two-dimensional shallow water equations, discretized using a second-order centered finite difference scheme on an Arakawa C grid. A second-order finite difference approximation for qf =qx can be written as 3 fiþ1=2 fi1=2 @f 1 3 @ f ðdxÞ ¼ ð13Þ @x i 24 @x3 dx The spatially discretized equations on an Arakawa C grid are given by: Continuity equation: nþ1 nþ1 nþ1 nþ1 Ujþ1=2;m Uj1=2;m Vj;mþ1=2 Vj;m1=2 df þD ¼0 þD dt ð j;mÞ dx dy
ð14Þ
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
x-Momentum equation: nþ1 g fjþ1;m fj;m dU ¼0 þ dt ð jþ1=2;mÞ dx y-Momentum equation: nþ1 g fj;mþ1 fj;m dV ¼0 þ dt ð j;mþ1=2Þ dy
2347
ð15Þ
ð16Þ
A discrete Fourier mode is introduced in each field variable such that Uj;m ¼ U0 ei½ jkx dxþmky dyxt V ¼ V ei½ jkx dxþmky dyxt j;m
ð17Þ
0
fj;m ¼ f0 e
i½ jkx dxþmky dyxt
Eqs. (14)–(16) reduce to 0 B ix B B B0 B B @ 2iD sinðkx dx=2Þ dx
0 ix 2iD sinðky dy=2Þ dy
1 2ig sinðkx dx=2Þ C0 1 dx C 2ig sinðky dy=2Þ C U0 C@ V0 A ¼ 0 C dy C f0 A ix
ð18Þ
A non-trivial solution of Eq. (18), assuming dy ¼ dx, requires 2 kx dx 2 ky dx x ðdxÞ 4xgD sin þ sin ¼0 2 2 3
2
ð19Þ
Solution of Eq. (19) gives three values of frequencies: ( x¼
2 0; dx
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi) ky dx kx dx 2 2 gD sin þ sin 2 2
ð20Þ
The phase speed is obtained from Eq. (20) by using the relation C ¼ x=keq , qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi keq ¼ kx2 þ ky2 . The x-directed group velocity CGx can be obtained by taking a partial derivative of Eq. (20) with respect to the wave number in the x-direction (kx) as 9 8 > > > > > > > > = < gD sinðkx dxÞ ffi ð21Þ CGx ¼ 0; sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > > k dx k dx y x > > > > 2 gD sin2 þ sin2 ; : 2 2
2348
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
The y-directed group velocity CGy can be obtained by taking a partial derivative of Eq. (20) with respect to the wave number in the y-direction (ky) as 9 8 > > > > > > > > = < gD sinðky dxÞ ffi ð22Þ CGy ¼ 0; sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > > k dx k dx y x > > > > 2 gD sin2 þ sin2 ; : 2 2 3.2. Dispersion analysis for a fourth-order finite difference scheme, using a five-point stencil Higher-order finite difference schemes can be derived by using the method of undefined coefficients (Kowalik and Murty, 1993). The following form of stencil can be used to define a fourth-order finite difference operator for the first derivative: fiþ1=2 fi1=2 fiþ3=2 fi3=2 @f þ b2 ð23Þ ¼ b1 @x i dx 3dx Expanding both sides of Eq. (23) in a Taylor series and collecting the terms of firstand third-order derivatives of f yields two equations, b1 þ b2 ¼ 1;
b1 þ 9 b2 ¼ 0
ð24Þ
The system of Eqs. (24) is solved to obtain a fourth-order finite difference approximation for qf =qx as 5 @f 9 fiþ1=2 fi1=2 1 fiþ3=2 fi3=2 3 @ f dx4 ¼ ð25Þ @x i 8 8 640 @x5 dx 3dx The finite difference approximation given by Eq. (25), when applied to Eqs. (4)–(6), assuming dy ¼ dx, yields: Continuity equation: 3 2 nþ1 nþ1 nþ1 nþ1 9 U U U U jþ1=2;m j1=2;m jþ3=2;m j3=2;m df D 5 þ 4 dt ð j;mÞ 8 dx 3dx 3 2 nþ1 nþ1 nþ1 nþ1 9 V V V V j;mþ1=2 j;m1=2 j;mþ3=2 j;m3=2 D 5 þ 4 8 dx 3dx ¼0
ð26Þ
x-Momentum equation:
3 2 nþ1 nþ1 nþ1 nþ1 9 f f f f jþ1;m j;m jþ2;m j1;m dU D 5¼0 þ 4 dt ð jþ1=2;mÞ 8 dx 3dx
ð27Þ
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2349
y-Momentum equation:
3 2 nþ1 nþ1 nþ1 nþ1 9 f f f f j;mþ1 j;m j;mþ2 j;m1 dV D 5¼0 þ 4 dt ð j;mþ1=2Þ 8 dx 3dx
ð28Þ
The discrete Fourier modes given by Eq. (17) are introduced in each field variable such that Eqs. (26)–(28) reduce to 1 0 10 ix 0 p1 g U0 @0 ix p2 g A@ V0 A ¼ 0 ð29Þ p1 D p2 D ix f0 where i q3 9q1 ; 4 3 kx dx ; q1 ¼ sin 2 3kx dx ; q3 ¼ sin 2
p1 ¼
i q4 9q2 ; 4 3 ky dx q2 ¼ sin ; 2 3ky dx q4 ¼ sin 2
p2 ¼
ð30Þ
A non-trivial solution of Eq. (29) requires that x3 þ
# gD " 2 x q4 q23 729 q21 þ q22 þ 54ðq1 q3 þ q2 q4 Þ ¼ 0 144
Solution of Eq. (31) gives three frequencies, n p3 o x ¼ 0; 3
ð31Þ
ð32Þ
where p3 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gDð36q22 þ q61 þ 12q42 þ 12q41 þ q62 þ 36q21 Þ
ð33Þ
The x-directed group velocity CGx can be obtained by taking a partial derivative of Eq. (32) with respect to the wave number in the x-direction (kx) as gDq5 q1 q21 þ 2 q21 þ 6 dx ¼ 0; CGx ¼ ð34Þ dkx 2p3 where kx dx q5 ¼ cos 2
ð35Þ
and q1 to q4 remain as defined in Eq. (30). The y-directed group velocity CGy can be obtained by taking a partial derivative of Eq. (32) with respect to the wave number in the y-direction (ky) as gDq6 q2 q22 þ 2 q22 þ 6 dx ¼ 0; CGy ¼ ð36Þ 2p3 dky
2350
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
where q2 and p3 are defined as above and ky dx q6 ¼ cos 2
ð37Þ
3.3. Dispersion analysis for a fourth-order compact difference scheme, using a threepoint stencil A three-point compact difference approximation of fourth order for the first derivative can be derived by using a relationship of the following form (Lele, 1992): fiþ1=2 fi1=2 @f @f @f ð38Þ þ þ a1 ¼ b1 a1 @x i1 @x i @x iþ1 dx The relationships between the coefficients a1 and b1 can be obtained by matching the Taylor series coefficients for the first- and third-order terms. 24a1 ¼ b1 ð39Þ 2a1 þ 1 ¼ b1 ; The system of Eqs. (39) is solved to obtain a compact finite difference approximation of fourth order for qf =qx: 5 1 @f @f 1 @f 12 fiþ1=2 fi1=2 17 4 @ f dx þ þ ¼ 22 @x i1 @x i 22 @x iþ1 11 5280 @x5 dx ð40Þ The fifth-order unmatched coefficient determines the formal truncation error of this approximation. The finite difference approximations given by Eq. (40), when applied to Eqs. (4)–(6), yield: Continuity equation: @f @U @V þD þD ¼0 ð41Þ @t ðj;mÞ @x ðj;mÞ @y ðj;mÞ Momentum equation in the x-direction: @U @f þg ¼0 @t ðjþ1=2;mÞ @x ðjþ1=2;mÞ Momentum equation in the y-direction: @V @f þg ¼0 @t ðj;mþ1=2Þ @y ðj;mþ1=2Þ
ð42Þ
ð43Þ
Auxiliary equations are needed to define the relationship between the variables and their derivatives appearing in Eqs. (41)–(43). For example, the relationship between ð@U=@xÞðj;mÞ and Ujþ1=2;m , can be written as ( ) Uðjþ1=2;mÞ Uðj1=2;mÞ 1 @U @U @U þ 22 þ ¼ 24 @x ðj1;mÞ @x ðj;mÞ @x ðjþ1;mÞ dx ð44Þ
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2351
Similarly, a relationship between ð@f=@xÞðjþ1=2;mÞ and fj;m can be written in the form 1 24
(
@f @x
@f þ 22 @x ðj1=2;mÞ fðjþ1;mÞ fðj;mÞ ¼ dx
þ
ðjþ1=2;mÞ
@f @x
)
ðjþ3=2;mÞ
ð45Þ
If a discrete Fourier mode is introduced in U, Eq. (44) reduces to
@U @x
¼ ðj;mÞ
24i sin kx 2dx Uðj;mÞ 11 þ cosðkx dxÞ
ð46Þ
Similarly, a discrete Fourier mode introduced for f reduces Eq. (45) to
@f @x
24i sin kx 2dx f ¼ 11 þ cosðkx dxÞ ðjþ1=2;mÞ ðjþ1=2;mÞ
ð47Þ
If discrete Fourier modes are introduced for all spatial derivatives appearing in Eqs. (41)–(43), using Eqs. (46) and (47) and assuming dy ¼ dx, Eqs. (41)–(43) reduce to 0
ix @0 q7 D
0 ix q8 D
1 10 q7 g U0 q 8 g A@ V 0 A ¼ 0 ix f0
ð48Þ
where
kx dx 24i sin 2 ; q7 ¼ 11 þ cosðkx dxÞ
ky dx 24i sin 2 q8 ¼ 11 þ cosðky dxÞ
ð49Þ
Solution of Eq. (48) gives three frequencies: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9 gD 36q21 þ q22 q41 þ 36q22 þ q42 q21 24q22 q21 = 2 x ¼ 0; 12 : ; q2 6 q21 6 8 <
ð50Þ
where q1 and q2 are defined as in Eq. (30). The x-directed group velocity CGx can be obtained by taking a partial derivative of Eq. (50) with respect to the wave number in the x-direction (kx) as ( !) dx 6gD dx q10 12q9 dx q1 q5 þ ¼ 0; 2 CGx ¼ ð51Þ 2 dkx q9 q2 6 q21 6 q2 6 q2 6 2
1
2352
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gD 36q21 þ q22 q41 þ 36q22 þ q42 q21 24q22 q21 ¼ 2q22 q31 q5 24q22 q1 q5 þ 36q1 q5 þ q42 q1 q5
q9 ¼ q10
ð52Þ
The y-directed group velocity CGy can be obtained by taking a partial derivative of Eq. (50) with respect to the wave number in the y-direction (ky) as ( !) dx 6gD dx q11 12q9 dx q2 q6 þ CGy ¼ ¼ 0; 2 ð53Þ 2 dky q9 q2 6 q21 6 q2 6 q2 6 2
where q1, q2, q6 and q9 remain as defined above and q11 ¼ q2 q41 q6 24q2 q21 q6 þ 36q2 q6 þ 2q32 q21 q6
1
ð54Þ
4. Comparison of phase and group speeds for second-order and fourth-order schemes The phase speed gives the speed and direction of individual waves and the group speed the speed of energy propagation. The normalized phase speed, defined as the ratio of the phase speed of the numerical wave to that of the analytical wave, is computed, and it is used as a measure to estimate the accuracy of the solution. The contours of normalized phase speed, as a function of the non-dimensional wave numbers in the x and y directions (kx dx and ky dx) for the second-order scheme are shown in Fig. 1. It is seen that the three-point second-order scheme (3p-FD) compares well with the analytical solution for non-dimensional wave numbers (kx dx) ranging from zero to 0.2, which corresponds to a grids per wavelength ratio ranging from infinity to 30 (Fig. 1). Fig. 2 shows the contours of the normalized phase speed for the five-point fourth-order scheme (5p-FD), and it is seen that the scheme is non-dispersive for non-dimensional wave numbers ranging from 0 to 1, which corresponds to a grids per wavelength ratio ranging from infinite to six grids per wavelength, with the three-point compact scheme (3p-CD) showing (Fig. 3) a marginally better comparison with the analytical solution than the five-point scheme. Figs. 4, 5 and 6 show, respectively, the normalized group speeds for the threepoint second-order, five-point fourth-order and three-point compact fourth-order schemes. It is seen that group velocities obtained using the three-point secondorder (3p-FD) scheme compare well with the analytical solution for non-dimensional wave numbers ranging from 0 to 0.2, which corresponds to a grids per wavelength ratio ranging from infinity to 30, while the group velocities obtained using five-point and compact fourth-order schemes compare well with the analytical solution for non-dimensional wave numbers ranging from 0 to 1, which corresponds to a grids per wavelength ratio ranging from infinity to 6.
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2353
Fig. 1. Contours of normalized phase speeds for the three-point second-order scheme.
5. Relative merits of various numerical schemes in hydrodynamic modeling A Fourier analysis of the spatially discretized shallow water equations, using the three spatial discretization schemes, namely, the three-point second-order scheme, five-point fourth-order scheme and three-point fourth-order compact scheme, shows that both fourth-order schemes are non-dispersive for grid resolutions greater than six grids per wavelength, whereas the second-order scheme is nondispersive for grid resolutions greater than 30 grids per wavelength.
Fig. 2. Contours of normalized phase speeds for the five-point fourth-order scheme.
2354
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
Fig. 3. Contours of normalized phase speeds for the fourth-order compact scheme.
In order to obtain unconditional numerical stability, the velocity derivatives in the continuity equation and the pressure gradient terms in the momentum equations are formulated implicitly in time. Such a technique results in a linear system of 3NxNy equations, while using second-order three-point and fourth-order fivepoint schemes. Even though the five-point fourth-order scheme has better dispersive properties, retaining the fourth-order accuracy becomes difficult near the boundaries, which are especially important in modeling narrow coves in typical
Fig. 4. Contours of normalized group speeds for the three-point second-order scheme.
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2355
Fig. 5. Contours of normalized group speeds for the five-point fourth-order scheme.
estuarine systems. The use of three-point compact differences results in a much larger linear system of 7NxNy equations, with the seven unknowns being u, v, f, qf=qx, qf=qy, qu=qx and qv=qy, which is computationally expensive. In order to reduce the system of equations to NxNy, Madala and Piacsek (1977) substituted the time-discretized momentum equations into the continuity equations to obtain a Helmholtz equation in terms of the surface elevation. The momentum
Fig. 6. Contours of normalized group speeds for the fourth-order compact scheme.
2356
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
equations are then solved explicitly to obtain the velocities. This technique has proved to be very efficient and has been used by Swanson (1986), Muin and Spaulding (1996) and others while solving two-dimensional systems of shallow water equations in Cartesian and boundary-fitted coordinates. The use of a Madala and Piacsek type implicit scheme with the three-point compact difference scheme results in a larger system of equations due to the auxiliary equations for the second-order derivatives in the Helmholtz equation. Thus, it is seen that the three-point fourth-order compact scheme is very accurate and has the advantages of using only three grid points. However, the large linear system of equations required for its implementation makes it at least eight times more expensive in terms of computational time than the second-order scheme. 6. Two-dimensional stability analysis A generalized second-order (in space and time) scheme, when applied to Eqs. (4)–(6), gives Continuity equation: n n1 ð1 þ bÞfnþ1 j;m ð2b þ 1Þfj;m þ bfj;m dt nþ1 n Ujþ1=2;m Uj1=2;m Ujþ1=2;m Uj1=2;m þ Dð1 hÞ þ Dh dx dx nþ1 n Vj;mþ1=2 Vj;m1=2 Vj;mþ1=2 Vj;m1=2 þ Dh þ Dð1 hÞ dy dy ¼0
ð55Þ
x-Momentum equation: nþ1 n n1 ð1 þ bÞUjþ1=2;m ð2b þ 1ÞUjþ1=2;m þ bUjþ1=2;m
dt nþ1 nþ1 fjþ1;m fj;m fjþ1;m fj;m þ gð1 hÞ þ gh dx dx ¼0
ð56Þ
y-Momentum equation: nþ1 n n1 ð1 þ bÞVj;mþ1=2 ð2b þ 1ÞVj;mþ1=2 þ bVj;mþ1=2
dt nþ1 nþ1 fj;mþ1 fj;m fj;mþ1 fj;m þ gh þ gð1 hÞ dy dy ¼0 where b and h are the weighting factors for the second-order scheme.
ð57Þ
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2357
After introducing the spatially and temporally discrete Fourier modes for all the variables as ^ n ¼ U0 ei½ jkx dxþmkx dxnx dt U j;m ^ n ¼ V0 ei½ jkx dxþmkx dxnx dt V j;m ^fn ¼ f0 ei½ jkx dxþmkx dxnx dt
ð58Þ
j;m
the fully discretized equations reduce to 1 " # 2ia1 g ð1 hÞk2 þ hk ð1 þ bÞk2 ð2b þ 1Þk þ b 0 C B dt dx C B " # C B 2 2 C B 2ia g ð1 hÞk þ hk ð1 þ bÞk ð2b þ 1Þk þ b 2 C B0 C B dy dt C B " # " # C B @ 2ia1 D ð1 hÞk2 þ hk 2ia2 D ð1 hÞk2 þ hk ð1 þ bÞk2 ð2b þ 1Þk þ b A dx dy dt 0 1 U0 B C
@ V0 A ¼ 0 ð59Þ 0
f0 where a1 ¼ sinðkx dx=2Þ, a2 ¼ sin ky dy=2 . A non-trivial solution of Eq. (59), assuming dy ¼ dx, requires that 4 k ð8hp12 þ 4h2 p12 þ b2 þ 4p12 þ 1 þ 2bÞ þ k3 ð8hp12 8h2 p12 4b2 2 6bÞ þ k2 ð6b2 þ 4h2 p12 þ 1 þ 6bÞþkð2b 4b2 Þ þ b2 ðk2 þ k2 b 2kb k þ bÞ ¼ 0
ð60Þ
where p12 ¼ sin
2
kx dx 2
þ sin
2
ky dx 2
c2o
ð61Þ
The three-time-level temporal scheme given in Muin and Spaulding (1996) for two dimensions can be obtained by substituting b ¼ 1=2 and h ¼ 0 in Eqs. (55)–(57). Thus, for the three-time-level scheme, Eq. (60) becomes ðk 1Þð3k 1Þ a11 k4 24k3 þ 22k2 8k þ 1 ¼ 0
ð62Þ
where ky dx kx dx a11 ¼ sin2 þ sin2 16c2o þ 9 2 2 pffiffiffiffiffiffiffi co ¼ gD dt=dx and co is the Courant number.
ð63Þ
2358
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
Solution of Eq. (62) gives six eigenvalues: 9 8 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 a12 a14 a15 2a13 > > > > > > > > > > a11 > > p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > = < 6 a12 a14 þ a15 þ 2a13 > ð64Þ fkg ¼ a11 > > > > > > > > > > 11 > > > > ; : 3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where a12 ¼ 4 1=ð9 a11 Þ, a13 ¼ 9 a11 , a14 ¼ a13 ð72 7a11 Þ;, a15 ¼ 216 33a11 þ a211 and a11 is defined by Eq. (63). The modulus of the eigenvalues, as a function of non-dimensional wave numbers in the x and y directions for all the eigenmodes, was seen to be less than one, thus indicating that the three-level scheme is unconditionally stable for all nondimensional wave and Courant numbers. A two-time-level Crank–Nicolson scheme can be obtained by putting b ¼ 0 and h ¼ 1=2 in Eq. (55)–(57). Thus, Eq. (62) becomes ðp11 þ 1Þk3 þ ðp11 3Þk2 ðp11 3Þk ðp11 þ 1Þ ¼ 0
ð65Þ
where 2 kx dx 2 ky dx p11 ¼ sin þ sin c2o 2 2 pffiffiffiffiffiffiffi where co ¼ gD dt=dx and co is the Courant number. Solution of Eq. (66) gives three eigenvalues: 8 pffiffiffiffiffiffi 9 < 1 p11 2i p11 = fkg ¼ 1 þ p11 : ; 1
ð66Þ
ð67Þ
The moduli of the eigenvalues given by Eq. (67) are less than one for all wave and Courant numbers, thus indicating that the scheme is unconditionally stable.
7. Two-dimensional wave deformation analysis The amplitude and phase components of the Fourier series representing the computed wave are compared with the real solution using the complex propagation factor defined by Leendertse (1967). The complex propagation factor gives the relation between the computed and real waves, in terms of their amplitudes and phases after a certain time interval, thus providing an estimate of the accuracy of the numerical scheme. The change in amplitude and phase lag of the wave from its original value, as time advances from t to t þ dt, is given by its eigenvalue (k). The eigenvalue of the
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2359
analytical wave is given by k ¼ eð ix dtÞ
ð68Þ
Taking the analytical frequency of the one-dimensional system from Eq. (10), the eigenvalue of the analytical solution can be written as pffiffiffiffiffi k ¼ eð i dt kx gDÞ ð69Þ The complex propagation factor, T(kx dx), is defined as the complex ratio of the computed wave in amplitude and phase to the physical wave after a time interval in which the physical wave propagates over its wavelength. The modulus of the propagation factor is a measure of the dissipation of the computed wave. The argument of the propagation factor gives a measure of the dispersion or phase lag of the computed wave when compared to the real wave. The propagation factor can thus be expressed as 0
Tðkx dxÞ ¼
eiðx tþkx xÞ eiðxtþkx xÞ
ð70Þ
where x ¼ 2p=kx , t ¼ 2p=x and x and x0 are the frequencies of the analytical and the numerical waves, respectively. Eq. (70) can also be written as 0
Tðkx dxÞ ¼ ei2pððx =xÞ1Þ The phase lag of the propagation factor can be calculated as 1 tan ½ImðkÞ=ReðkÞ argðTðkx dxÞÞ ¼ 2p 1 co kx dx
ð71Þ
ð72Þ
The propagation of a physical wave having a frequency of x over its wavelength L requires m ¼ 2p=x dt ¼ 2p=co kx dx time steps. The modulus of the propagation factor is given by m jknum j ¼ jknum jm ð73Þ jTðkx dxÞj ¼ jkana j pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 2 k ¼ eð i dt kx gDðkx þky ÞÞ ð74Þ Knowing the eigenvalue of the numerical scheme from Eq. (67) for the two-timelevel scheme, the modulus and the phase lag of the propagation factor can be calculated. The modulus of the propagation factor for the Crank–Nicolson scheme is found be unity for all non-dimensional wave and Courant numbers. The phase lag of the propagation factor is found to increase with increase in Courant number and an increase in grid resolution, as shown in Figs. 7–10. The numerical dispersion can be reduced by decreasing the time step and increasing the grid resolution in both directions.
2360
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
Fig. 7. Contours of phase lag of propagation factor for k1 in the Crank–Nicolson scheme for a Courant number of 0.5.
8. Conclusions A Fourier analysis of the spatially discretized shallow water equations shows that the five-point fourth-order and three-point fourth-order compact schemes are less dispersive than the conventional three-point second-order scheme, for a given grid resolution and Courant number. The use of a three-point compact scheme is
Fig. 8. Contours of phase lag of propagation factor for k1 in the Crank–Nicolson scheme for a Courant number of 1.
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2361
Fig. 9. Contours of phase lag of propagation factor for k1 in the Crank–Nicolson scheme for a Courant number of 2.
computationally expensive, while it will be difficult to maintain fourth-order accuracy near the boundary while using the five-point scheme. A von Neumann stability analysis of the three-time-level scheme for one-dimensional and twodimensional systems of shallow water equations shows that the scheme is unconditionally stable.
Fig. 10. Contours of phase lag of propagation factor for k1 in the Crank–Nicolson scheme for a Courant number of 5.
2362
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
A wave deformation analysis of the two-time-level Crank–Nicolson scheme for one-dimensional and two-dimensional systems of shallow water equations shows that the scheme is non-dissipative, independent of the Courant number and grid resolution used. The phase error or the dispersion of the scheme is found to be small for grid resolutions greater than 20 grids per wavelength and Courant numbers less than 1. Acknowledgements The first author wishes to thank the University of Rhode Island (URI) for providing him with a graduate fellowship for a year. References Arakawa, A., Lamb, V.R., 1977. Computational design of the basic dynamical processes of the UCLA General Circulation Model. Methods in Computational Physics 17, 173–265. Casulli, V., Cattani, E., 1994. Stability, accuracy and efficiency of a semi-implicit method for threedimensional shallow water flow. Computers and Mathematics with Applications 27 (4), 99–112. Foreman, M.G.G., 1984. A two-dimensional dispersion analysis of selected methods for solving linearized shallow water equations. Journal of Computational Physics 56, 287–323. Hirsh, R.S., 1975. Higher order accurate difference solution of fluid mechanics problems by a compact differencing technique. Journal of Computational Physics 9 (1), 90–109. Kowalik, Z., Murty, T.S., 1993. Numerical Modeling of Ocean Dynamics. World Scientific, Singapore. Lele, S.K., 1992. Compact finite difference schemes with spectral-like resolution. Journal of Computational Physics 103, 16–42. Leendertse, J.J., 1967. Aspects of computational model for long period water wave propagation, Report RM-5294-PR, The Rand Corp., Santa Monica, CA. Madala, R.V., Piacsek, S.A., 1977. A semi implicit numerical model for barocline oceans. Journal of Computational Physics 23, 167–178. Messinger, F., Arakawa, A., 1976. Numerical Methods Used in Atmospheric Models, vol. I. GARP Publication Series, No. 17. WMO–ICSU Joint Organizing Committee. Muin, M., Spaulding, M.L., 1996. Two-dimensional boundary-fitted circulation model in spherical coordinates. Journal of Hydraulic Engineering, ASCE 122 (9), 512–521. Song, Y., Tang, T., 1993. Dispersion and group velocity in numerical schemes for three-dimensional hydrodynamic equations. Journal of Computational Physics 105, 72–82. Swanson, J.C., 1986. A three-dimensional numerical model system of coastal circulation and water quality. Ph.D. Thesis, University of Rhode Island, Kingston, RI.
Dispersion and stability analyses of the linearized two-dimensional shallow water equations in Cartesian coordinates S. Sankaranarayanan a,, Malcolm L. Spaulding b a
Applied Science Associates, 70 Dean Knauss Drive, Narragansett, RI 02882, USA Ocean Engineering, University of Rhode Island, Narragansett, RI 02882, USA
b
Received 18 November 2002; accepted 16 April 2003
Abstract In the present study, a Fourier analysis is used to develop expressions for phase and group speeds for both continuous and discretized, linearized two-dimensional shallow water equations, in Cartesian coordinates. The phase and group speeds of the discrete equations, discretized using a three-point scheme of second order, five-point scheme of fourth order and a three-point compact scheme of fourth order in an Arakawa C grid, are calculated and compared with the corresponding values obtained for the continuous system. The threepoint second-order scheme is found to be non-dispersive with grid resolutions greater than 30 grids per wavelength, while both the fourth-order schemes are non-dispersive with grid resolutions greater than six grids per wavelength. A von Neumann stability analysis of the two- and three-time-level temporal schemes showed that both schemes are stable. A wave deformation analysis of the two-time-level Crank–Nicolson scheme for one-dimensional and two-dimensional systems of shallow water equations shows that the scheme is nondispersive, independent of the Courant number and grid resolution used. The phase error or the dispersion of the scheme decreases with a decrease in the time step or an increase in grid resolution. # 2003 Elsevier Ltd. All rights reserved. Keywords: Stability; Dispersion; Cartesian coordinates; Compact difference
Corresponding author. Formerly at Ocean Engineering, University of Rhode Island, Narragansett, RI 02882, USA. Tel.: +1-401-7896224; fax: +1-401-7891932. E-mail address: [email protected] (S. Sankaranarayanan).
0029-8018/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0029-8018(03)00103-3
2344
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
1. Introduction The dispersion and stability analyses for various schemes used in solving shallow water equations have been carried out in the past using a Fourier analysis of the differencing techniques, while solving the governing equations in a Cartesian coordinate system (Leendertse, 1967; Messinger and Arakawa, 1976). Foreman (1984) estimated the accuracy of some finite difference and finite element techniques, by comparing the numerical and analytical plane wave solutions using a Fourier analysis. Song and Tang (1993) analyzed the dispersion and group velocity for linearized three-dimensional hydrodynamic equations discretized using a hybrid scheme on Arakawa’s B and C spatial grid types. It is noted that Arakawa and Lamb (1977) and Song and Tang (1993) evaluated the dispersion from a spatially discretized system, while retaining the time derivative in its differential form; Leenderstse (1967) and Foreman (1984), however, evaluated the numerical dispersion for fully discretized equations, from their eigenvalues obtained through a Fourier analysis. Messinger and Arakawa (1976) found that fourth-order numerical schemes have lower numerical dispersion when compared with traditional second-order schemes. Hirsch (1975) used a three-point compact scheme of fourth order, whereas traditional fourth-order methods require a five-point stencil. In compact difference methods, the function and all of its derivatives are considered as unknowns for obtaining higher-order accuracy with the minimum number of nodes, whereas traditional finite difference methods require additional nodes to achieve higher-order accuracy. Compact difference methods thus avoid the difficulties of having to consider additional fictitious nodes, near the boundaries, to maintain higher-order accuracy. Hirsch (1975) also noted that coefficients of the truncation error term are smaller in the compact scheme, for the same differencing order. Later, Lele (1992), while solving the one-dimensional advection equation, found that the three-point compact scheme has marginally better dispersion properties than the traditional five-point scheme. The three-point second-order scheme is employed in most coastal and ocean models. The five-point fourth-order scheme (Messinger and Arakawa, 1976) and the three-point fourth-order compact scheme (Lele, 1992) were shown to possess less dispersion when solving the advection equation. The three-time-level (Muin and Spaulding, 1996) and two-time-level Crank–Nicolson scheme (Casulli and Cattani, 1994) are second-order accurate in time and space and have been successfully used in coastal modeling. In the present study, the three spatial discretization schemes and two temporal discretization schemes discussed above are considered for stability and dispersion analysis. One-dimensional and two-dimensional wave deformation analyses, based on a complex propagation factor developed by Leendertse (1967), are used to estimate the amplitude and phase errors of the fully discretized equations.
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2345
2. Governing equations The linearized, depth-averaged, two-dimensional shallow water equations, in Cartesian coordinates are given as: Continuity equation: @f @UD @VD þ þ ¼0 @t @x @y Momentum equation in the x-direction: @UD @f ¼ gD @t @x Momentum equation in the y-direction: @VD @f ¼ gD @t @y
ð1Þ
ð2Þ
ð3Þ
where U(x, y, t) and V(x, y, t) are the vertically averaged velocities in the x and y directions, respectively, fðx; y; tÞ is the water surface elevation, g is the acceleration due to gravity, and D is the total water depth. 3. Dispersion analysis of the continuous system Taking the case of constant depth, Eqs. (1)–(3) become: Continuity equation: @f @U @V þD þD ¼0 @t @x @y Momentum equation in the x-direction: @U @f þg ¼0 @t @x Momentum equation in the y-direction: @V @f þg ¼0 @t @y
ð4Þ
ð5Þ
ð6Þ
A continuous Fourier mode is introduced in each field variable, assuming that the governing equations are defined on an infinite spatial domain or with periodic boundary conditions on a finite domain, and represented by U ¼ U0 ei½kx xþky yxt V ¼ V0 ei½kx xþky yxt f ¼ f0 ei½kx xþky yxt
ð7Þ
2346
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
Eqs. (4)–(6) reduce to 0 10 1 ix 0 igkx U0 @0 ix igky A@ V0 A ¼ 0 iDkx iDky ix f0 where kx and ky are the wave numbers in the x and y direction, respectively. A non-trivial solution of Eq. (8) requires that x3 xgD kx2 þ ky2 ¼ 0 Solution of Eq. (9) gives three frequencies: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 x ¼ 0; gD ky þ kx
ð8Þ
ð9Þ
ð10Þ
The x-directed group velocity CGx can be obtained by taking a partial derivative of Eq. (10) with respect to the wave number in the x-direction (kx) as 9 8 > > > > = < dx gDkx ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ð11Þ CGx ¼ ¼ 0; > dkx > > ; : gD ky2 þ kx2 > The y-directed group velocity CGy can be obtained by taking a partial derivative of Eq. (10) with respect to the wave number in the y-direction (ky) as 9 8 > > > > = gDky dx < ð12Þ CGy ¼ ¼ 0; rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > dky > > ; : gD ky2 þ kx2 > 3.1. Dispersion analysis for a second-order finite difference scheme In this section, expressions for phase and group speeds are developed for linearized two-dimensional shallow water equations, discretized using a second-order centered finite difference scheme on an Arakawa C grid. A second-order finite difference approximation for qf =qx can be written as 3 fiþ1=2 fi1=2 @f 1 3 @ f ðdxÞ ¼ ð13Þ @x i 24 @x3 dx The spatially discretized equations on an Arakawa C grid are given by: Continuity equation: nþ1 nþ1 nþ1 nþ1 Ujþ1=2;m Uj1=2;m Vj;mþ1=2 Vj;m1=2 df þD ¼0 þD dt ð j;mÞ dx dy
ð14Þ
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
x-Momentum equation: nþ1 g fjþ1;m fj;m dU ¼0 þ dt ð jþ1=2;mÞ dx y-Momentum equation: nþ1 g fj;mþ1 fj;m dV ¼0 þ dt ð j;mþ1=2Þ dy
2347
ð15Þ
ð16Þ
A discrete Fourier mode is introduced in each field variable such that Uj;m ¼ U0 ei½ jkx dxþmky dyxt V ¼ V ei½ jkx dxþmky dyxt j;m
ð17Þ
0
fj;m ¼ f0 e
i½ jkx dxþmky dyxt
Eqs. (14)–(16) reduce to 0 B ix B B B0 B B @ 2iD sinðkx dx=2Þ dx
0 ix 2iD sinðky dy=2Þ dy
1 2ig sinðkx dx=2Þ C0 1 dx C 2ig sinðky dy=2Þ C U0 C@ V0 A ¼ 0 C dy C f0 A ix
ð18Þ
A non-trivial solution of Eq. (18), assuming dy ¼ dx, requires 2 kx dx 2 ky dx x ðdxÞ 4xgD sin þ sin ¼0 2 2 3
2
ð19Þ
Solution of Eq. (19) gives three values of frequencies: ( x¼
2 0; dx
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi) ky dx kx dx 2 2 gD sin þ sin 2 2
ð20Þ
The phase speed is obtained from Eq. (20) by using the relation C ¼ x=keq , qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi keq ¼ kx2 þ ky2 . The x-directed group velocity CGx can be obtained by taking a partial derivative of Eq. (20) with respect to the wave number in the x-direction (kx) as 9 8 > > > > > > > > = < gD sinðkx dxÞ ffi ð21Þ CGx ¼ 0; sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > > k dx k dx y x > > > > 2 gD sin2 þ sin2 ; : 2 2
2348
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
The y-directed group velocity CGy can be obtained by taking a partial derivative of Eq. (20) with respect to the wave number in the y-direction (ky) as 9 8 > > > > > > > > = < gD sinðky dxÞ ffi ð22Þ CGy ¼ 0; sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > > k dx k dx y x > > > > 2 gD sin2 þ sin2 ; : 2 2 3.2. Dispersion analysis for a fourth-order finite difference scheme, using a five-point stencil Higher-order finite difference schemes can be derived by using the method of undefined coefficients (Kowalik and Murty, 1993). The following form of stencil can be used to define a fourth-order finite difference operator for the first derivative: fiþ1=2 fi1=2 fiþ3=2 fi3=2 @f þ b2 ð23Þ ¼ b1 @x i dx 3dx Expanding both sides of Eq. (23) in a Taylor series and collecting the terms of firstand third-order derivatives of f yields two equations, b1 þ b2 ¼ 1;
b1 þ 9 b2 ¼ 0
ð24Þ
The system of Eqs. (24) is solved to obtain a fourth-order finite difference approximation for qf =qx as 5 @f 9 fiþ1=2 fi1=2 1 fiþ3=2 fi3=2 3 @ f dx4 ¼ ð25Þ @x i 8 8 640 @x5 dx 3dx The finite difference approximation given by Eq. (25), when applied to Eqs. (4)–(6), assuming dy ¼ dx, yields: Continuity equation: 3 2 nþ1 nþ1 nþ1 nþ1 9 U U U U jþ1=2;m j1=2;m jþ3=2;m j3=2;m df D 5 þ 4 dt ð j;mÞ 8 dx 3dx 3 2 nþ1 nþ1 nþ1 nþ1 9 V V V V j;mþ1=2 j;m1=2 j;mþ3=2 j;m3=2 D 5 þ 4 8 dx 3dx ¼0
ð26Þ
x-Momentum equation:
3 2 nþ1 nþ1 nþ1 nþ1 9 f f f f jþ1;m j;m jþ2;m j1;m dU D 5¼0 þ 4 dt ð jþ1=2;mÞ 8 dx 3dx
ð27Þ
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2349
y-Momentum equation:
3 2 nþ1 nþ1 nþ1 nþ1 9 f f f f j;mþ1 j;m j;mþ2 j;m1 dV D 5¼0 þ 4 dt ð j;mþ1=2Þ 8 dx 3dx
ð28Þ
The discrete Fourier modes given by Eq. (17) are introduced in each field variable such that Eqs. (26)–(28) reduce to 1 0 10 ix 0 p1 g U0 @0 ix p2 g A@ V0 A ¼ 0 ð29Þ p1 D p2 D ix f0 where i q3 9q1 ; 4 3 kx dx ; q1 ¼ sin 2 3kx dx ; q3 ¼ sin 2
p1 ¼
i q4 9q2 ; 4 3 ky dx q2 ¼ sin ; 2 3ky dx q4 ¼ sin 2
p2 ¼
ð30Þ
A non-trivial solution of Eq. (29) requires that x3 þ
# gD " 2 x q4 q23 729 q21 þ q22 þ 54ðq1 q3 þ q2 q4 Þ ¼ 0 144
Solution of Eq. (31) gives three frequencies, n p3 o x ¼ 0; 3
ð31Þ
ð32Þ
where p3 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gDð36q22 þ q61 þ 12q42 þ 12q41 þ q62 þ 36q21 Þ
ð33Þ
The x-directed group velocity CGx can be obtained by taking a partial derivative of Eq. (32) with respect to the wave number in the x-direction (kx) as gDq5 q1 q21 þ 2 q21 þ 6 dx ¼ 0; CGx ¼ ð34Þ dkx 2p3 where kx dx q5 ¼ cos 2
ð35Þ
and q1 to q4 remain as defined in Eq. (30). The y-directed group velocity CGy can be obtained by taking a partial derivative of Eq. (32) with respect to the wave number in the y-direction (ky) as gDq6 q2 q22 þ 2 q22 þ 6 dx ¼ 0; CGy ¼ ð36Þ 2p3 dky
2350
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
where q2 and p3 are defined as above and ky dx q6 ¼ cos 2
ð37Þ
3.3. Dispersion analysis for a fourth-order compact difference scheme, using a threepoint stencil A three-point compact difference approximation of fourth order for the first derivative can be derived by using a relationship of the following form (Lele, 1992): fiþ1=2 fi1=2 @f @f @f ð38Þ þ þ a1 ¼ b1 a1 @x i1 @x i @x iþ1 dx The relationships between the coefficients a1 and b1 can be obtained by matching the Taylor series coefficients for the first- and third-order terms. 24a1 ¼ b1 ð39Þ 2a1 þ 1 ¼ b1 ; The system of Eqs. (39) is solved to obtain a compact finite difference approximation of fourth order for qf =qx: 5 1 @f @f 1 @f 12 fiþ1=2 fi1=2 17 4 @ f dx þ þ ¼ 22 @x i1 @x i 22 @x iþ1 11 5280 @x5 dx ð40Þ The fifth-order unmatched coefficient determines the formal truncation error of this approximation. The finite difference approximations given by Eq. (40), when applied to Eqs. (4)–(6), yield: Continuity equation: @f @U @V þD þD ¼0 ð41Þ @t ðj;mÞ @x ðj;mÞ @y ðj;mÞ Momentum equation in the x-direction: @U @f þg ¼0 @t ðjþ1=2;mÞ @x ðjþ1=2;mÞ Momentum equation in the y-direction: @V @f þg ¼0 @t ðj;mþ1=2Þ @y ðj;mþ1=2Þ
ð42Þ
ð43Þ
Auxiliary equations are needed to define the relationship between the variables and their derivatives appearing in Eqs. (41)–(43). For example, the relationship between ð@U=@xÞðj;mÞ and Ujþ1=2;m , can be written as ( ) Uðjþ1=2;mÞ Uðj1=2;mÞ 1 @U @U @U þ 22 þ ¼ 24 @x ðj1;mÞ @x ðj;mÞ @x ðjþ1;mÞ dx ð44Þ
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2351
Similarly, a relationship between ð@f=@xÞðjþ1=2;mÞ and fj;m can be written in the form 1 24
(
@f @x
@f þ 22 @x ðj1=2;mÞ fðjþ1;mÞ fðj;mÞ ¼ dx
þ
ðjþ1=2;mÞ
@f @x
)
ðjþ3=2;mÞ
ð45Þ
If a discrete Fourier mode is introduced in U, Eq. (44) reduces to
@U @x
¼ ðj;mÞ
24i sin kx 2dx Uðj;mÞ 11 þ cosðkx dxÞ
ð46Þ
Similarly, a discrete Fourier mode introduced for f reduces Eq. (45) to
@f @x
24i sin kx 2dx f ¼ 11 þ cosðkx dxÞ ðjþ1=2;mÞ ðjþ1=2;mÞ
ð47Þ
If discrete Fourier modes are introduced for all spatial derivatives appearing in Eqs. (41)–(43), using Eqs. (46) and (47) and assuming dy ¼ dx, Eqs. (41)–(43) reduce to 0
ix @0 q7 D
0 ix q8 D
1 10 q7 g U0 q 8 g A@ V 0 A ¼ 0 ix f0
ð48Þ
where
kx dx 24i sin 2 ; q7 ¼ 11 þ cosðkx dxÞ
ky dx 24i sin 2 q8 ¼ 11 þ cosðky dxÞ
ð49Þ
Solution of Eq. (48) gives three frequencies: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9 gD 36q21 þ q22 q41 þ 36q22 þ q42 q21 24q22 q21 = 2 x ¼ 0; 12 : ; q2 6 q21 6 8 <
ð50Þ
where q1 and q2 are defined as in Eq. (30). The x-directed group velocity CGx can be obtained by taking a partial derivative of Eq. (50) with respect to the wave number in the x-direction (kx) as ( !) dx 6gD dx q10 12q9 dx q1 q5 þ ¼ 0; 2 CGx ¼ ð51Þ 2 dkx q9 q2 6 q21 6 q2 6 q2 6 2
1
2352
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gD 36q21 þ q22 q41 þ 36q22 þ q42 q21 24q22 q21 ¼ 2q22 q31 q5 24q22 q1 q5 þ 36q1 q5 þ q42 q1 q5
q9 ¼ q10
ð52Þ
The y-directed group velocity CGy can be obtained by taking a partial derivative of Eq. (50) with respect to the wave number in the y-direction (ky) as ( !) dx 6gD dx q11 12q9 dx q2 q6 þ CGy ¼ ¼ 0; 2 ð53Þ 2 dky q9 q2 6 q21 6 q2 6 q2 6 2
where q1, q2, q6 and q9 remain as defined above and q11 ¼ q2 q41 q6 24q2 q21 q6 þ 36q2 q6 þ 2q32 q21 q6
1
ð54Þ
4. Comparison of phase and group speeds for second-order and fourth-order schemes The phase speed gives the speed and direction of individual waves and the group speed the speed of energy propagation. The normalized phase speed, defined as the ratio of the phase speed of the numerical wave to that of the analytical wave, is computed, and it is used as a measure to estimate the accuracy of the solution. The contours of normalized phase speed, as a function of the non-dimensional wave numbers in the x and y directions (kx dx and ky dx) for the second-order scheme are shown in Fig. 1. It is seen that the three-point second-order scheme (3p-FD) compares well with the analytical solution for non-dimensional wave numbers (kx dx) ranging from zero to 0.2, which corresponds to a grids per wavelength ratio ranging from infinity to 30 (Fig. 1). Fig. 2 shows the contours of the normalized phase speed for the five-point fourth-order scheme (5p-FD), and it is seen that the scheme is non-dispersive for non-dimensional wave numbers ranging from 0 to 1, which corresponds to a grids per wavelength ratio ranging from infinite to six grids per wavelength, with the three-point compact scheme (3p-CD) showing (Fig. 3) a marginally better comparison with the analytical solution than the five-point scheme. Figs. 4, 5 and 6 show, respectively, the normalized group speeds for the threepoint second-order, five-point fourth-order and three-point compact fourth-order schemes. It is seen that group velocities obtained using the three-point secondorder (3p-FD) scheme compare well with the analytical solution for non-dimensional wave numbers ranging from 0 to 0.2, which corresponds to a grids per wavelength ratio ranging from infinity to 30, while the group velocities obtained using five-point and compact fourth-order schemes compare well with the analytical solution for non-dimensional wave numbers ranging from 0 to 1, which corresponds to a grids per wavelength ratio ranging from infinity to 6.
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2353
Fig. 1. Contours of normalized phase speeds for the three-point second-order scheme.
5. Relative merits of various numerical schemes in hydrodynamic modeling A Fourier analysis of the spatially discretized shallow water equations, using the three spatial discretization schemes, namely, the three-point second-order scheme, five-point fourth-order scheme and three-point fourth-order compact scheme, shows that both fourth-order schemes are non-dispersive for grid resolutions greater than six grids per wavelength, whereas the second-order scheme is nondispersive for grid resolutions greater than 30 grids per wavelength.
Fig. 2. Contours of normalized phase speeds for the five-point fourth-order scheme.
2354
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
Fig. 3. Contours of normalized phase speeds for the fourth-order compact scheme.
In order to obtain unconditional numerical stability, the velocity derivatives in the continuity equation and the pressure gradient terms in the momentum equations are formulated implicitly in time. Such a technique results in a linear system of 3NxNy equations, while using second-order three-point and fourth-order fivepoint schemes. Even though the five-point fourth-order scheme has better dispersive properties, retaining the fourth-order accuracy becomes difficult near the boundaries, which are especially important in modeling narrow coves in typical
Fig. 4. Contours of normalized group speeds for the three-point second-order scheme.
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2355
Fig. 5. Contours of normalized group speeds for the five-point fourth-order scheme.
estuarine systems. The use of three-point compact differences results in a much larger linear system of 7NxNy equations, with the seven unknowns being u, v, f, qf=qx, qf=qy, qu=qx and qv=qy, which is computationally expensive. In order to reduce the system of equations to NxNy, Madala and Piacsek (1977) substituted the time-discretized momentum equations into the continuity equations to obtain a Helmholtz equation in terms of the surface elevation. The momentum
Fig. 6. Contours of normalized group speeds for the fourth-order compact scheme.
2356
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
equations are then solved explicitly to obtain the velocities. This technique has proved to be very efficient and has been used by Swanson (1986), Muin and Spaulding (1996) and others while solving two-dimensional systems of shallow water equations in Cartesian and boundary-fitted coordinates. The use of a Madala and Piacsek type implicit scheme with the three-point compact difference scheme results in a larger system of equations due to the auxiliary equations for the second-order derivatives in the Helmholtz equation. Thus, it is seen that the three-point fourth-order compact scheme is very accurate and has the advantages of using only three grid points. However, the large linear system of equations required for its implementation makes it at least eight times more expensive in terms of computational time than the second-order scheme. 6. Two-dimensional stability analysis A generalized second-order (in space and time) scheme, when applied to Eqs. (4)–(6), gives Continuity equation: n n1 ð1 þ bÞfnþ1 j;m ð2b þ 1Þfj;m þ bfj;m dt nþ1 n Ujþ1=2;m Uj1=2;m Ujþ1=2;m Uj1=2;m þ Dð1 hÞ þ Dh dx dx nþ1 n Vj;mþ1=2 Vj;m1=2 Vj;mþ1=2 Vj;m1=2 þ Dh þ Dð1 hÞ dy dy ¼0
ð55Þ
x-Momentum equation: nþ1 n n1 ð1 þ bÞUjþ1=2;m ð2b þ 1ÞUjþ1=2;m þ bUjþ1=2;m
dt nþ1 nþ1 fjþ1;m fj;m fjþ1;m fj;m þ gð1 hÞ þ gh dx dx ¼0
ð56Þ
y-Momentum equation: nþ1 n n1 ð1 þ bÞVj;mþ1=2 ð2b þ 1ÞVj;mþ1=2 þ bVj;mþ1=2
dt nþ1 nþ1 fj;mþ1 fj;m fj;mþ1 fj;m þ gh þ gð1 hÞ dy dy ¼0 where b and h are the weighting factors for the second-order scheme.
ð57Þ
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2357
After introducing the spatially and temporally discrete Fourier modes for all the variables as ^ n ¼ U0 ei½ jkx dxþmkx dxnx dt U j;m ^ n ¼ V0 ei½ jkx dxþmkx dxnx dt V j;m ^fn ¼ f0 ei½ jkx dxþmkx dxnx dt
ð58Þ
j;m
the fully discretized equations reduce to 1 " # 2ia1 g ð1 hÞk2 þ hk ð1 þ bÞk2 ð2b þ 1Þk þ b 0 C B dt dx C B " # C B 2 2 C B 2ia g ð1 hÞk þ hk ð1 þ bÞk ð2b þ 1Þk þ b 2 C B0 C B dy dt C B " # " # C B @ 2ia1 D ð1 hÞk2 þ hk 2ia2 D ð1 hÞk2 þ hk ð1 þ bÞk2 ð2b þ 1Þk þ b A dx dy dt 0 1 U0 B C
@ V0 A ¼ 0 ð59Þ 0
f0 where a1 ¼ sinðkx dx=2Þ, a2 ¼ sin ky dy=2 . A non-trivial solution of Eq. (59), assuming dy ¼ dx, requires that 4 k ð8hp12 þ 4h2 p12 þ b2 þ 4p12 þ 1 þ 2bÞ þ k3 ð8hp12 8h2 p12 4b2 2 6bÞ þ k2 ð6b2 þ 4h2 p12 þ 1 þ 6bÞþkð2b 4b2 Þ þ b2 ðk2 þ k2 b 2kb k þ bÞ ¼ 0
ð60Þ
where p12 ¼ sin
2
kx dx 2
þ sin
2
ky dx 2
c2o
ð61Þ
The three-time-level temporal scheme given in Muin and Spaulding (1996) for two dimensions can be obtained by substituting b ¼ 1=2 and h ¼ 0 in Eqs. (55)–(57). Thus, for the three-time-level scheme, Eq. (60) becomes ðk 1Þð3k 1Þ a11 k4 24k3 þ 22k2 8k þ 1 ¼ 0
ð62Þ
where ky dx kx dx a11 ¼ sin2 þ sin2 16c2o þ 9 2 2 pffiffiffiffiffiffiffi co ¼ gD dt=dx and co is the Courant number.
ð63Þ
2358
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
Solution of Eq. (62) gives six eigenvalues: 9 8 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 a12 a14 a15 2a13 > > > > > > > > > > a11 > > p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > = < 6 a12 a14 þ a15 þ 2a13 > ð64Þ fkg ¼ a11 > > > > > > > > > > 11 > > > > ; : 3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where a12 ¼ 4 1=ð9 a11 Þ, a13 ¼ 9 a11 , a14 ¼ a13 ð72 7a11 Þ;, a15 ¼ 216 33a11 þ a211 and a11 is defined by Eq. (63). The modulus of the eigenvalues, as a function of non-dimensional wave numbers in the x and y directions for all the eigenmodes, was seen to be less than one, thus indicating that the three-level scheme is unconditionally stable for all nondimensional wave and Courant numbers. A two-time-level Crank–Nicolson scheme can be obtained by putting b ¼ 0 and h ¼ 1=2 in Eq. (55)–(57). Thus, Eq. (62) becomes ðp11 þ 1Þk3 þ ðp11 3Þk2 ðp11 3Þk ðp11 þ 1Þ ¼ 0
ð65Þ
where 2 kx dx 2 ky dx p11 ¼ sin þ sin c2o 2 2 pffiffiffiffiffiffiffi where co ¼ gD dt=dx and co is the Courant number. Solution of Eq. (66) gives three eigenvalues: 8 pffiffiffiffiffiffi 9 < 1 p11 2i p11 = fkg ¼ 1 þ p11 : ; 1
ð66Þ
ð67Þ
The moduli of the eigenvalues given by Eq. (67) are less than one for all wave and Courant numbers, thus indicating that the scheme is unconditionally stable.
7. Two-dimensional wave deformation analysis The amplitude and phase components of the Fourier series representing the computed wave are compared with the real solution using the complex propagation factor defined by Leendertse (1967). The complex propagation factor gives the relation between the computed and real waves, in terms of their amplitudes and phases after a certain time interval, thus providing an estimate of the accuracy of the numerical scheme. The change in amplitude and phase lag of the wave from its original value, as time advances from t to t þ dt, is given by its eigenvalue (k). The eigenvalue of the
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2359
analytical wave is given by k ¼ eð ix dtÞ
ð68Þ
Taking the analytical frequency of the one-dimensional system from Eq. (10), the eigenvalue of the analytical solution can be written as pffiffiffiffiffi k ¼ eð i dt kx gDÞ ð69Þ The complex propagation factor, T(kx dx), is defined as the complex ratio of the computed wave in amplitude and phase to the physical wave after a time interval in which the physical wave propagates over its wavelength. The modulus of the propagation factor is a measure of the dissipation of the computed wave. The argument of the propagation factor gives a measure of the dispersion or phase lag of the computed wave when compared to the real wave. The propagation factor can thus be expressed as 0
Tðkx dxÞ ¼
eiðx tþkx xÞ eiðxtþkx xÞ
ð70Þ
where x ¼ 2p=kx , t ¼ 2p=x and x and x0 are the frequencies of the analytical and the numerical waves, respectively. Eq. (70) can also be written as 0
Tðkx dxÞ ¼ ei2pððx =xÞ1Þ The phase lag of the propagation factor can be calculated as 1 tan ½ImðkÞ=ReðkÞ argðTðkx dxÞÞ ¼ 2p 1 co kx dx
ð71Þ
ð72Þ
The propagation of a physical wave having a frequency of x over its wavelength L requires m ¼ 2p=x dt ¼ 2p=co kx dx time steps. The modulus of the propagation factor is given by m jknum j ¼ jknum jm ð73Þ jTðkx dxÞj ¼ jkana j pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 2 k ¼ eð i dt kx gDðkx þky ÞÞ ð74Þ Knowing the eigenvalue of the numerical scheme from Eq. (67) for the two-timelevel scheme, the modulus and the phase lag of the propagation factor can be calculated. The modulus of the propagation factor for the Crank–Nicolson scheme is found be unity for all non-dimensional wave and Courant numbers. The phase lag of the propagation factor is found to increase with increase in Courant number and an increase in grid resolution, as shown in Figs. 7–10. The numerical dispersion can be reduced by decreasing the time step and increasing the grid resolution in both directions.
2360
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
Fig. 7. Contours of phase lag of propagation factor for k1 in the Crank–Nicolson scheme for a Courant number of 0.5.
8. Conclusions A Fourier analysis of the spatially discretized shallow water equations shows that the five-point fourth-order and three-point fourth-order compact schemes are less dispersive than the conventional three-point second-order scheme, for a given grid resolution and Courant number. The use of a three-point compact scheme is
Fig. 8. Contours of phase lag of propagation factor for k1 in the Crank–Nicolson scheme for a Courant number of 1.
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
2361
Fig. 9. Contours of phase lag of propagation factor for k1 in the Crank–Nicolson scheme for a Courant number of 2.
computationally expensive, while it will be difficult to maintain fourth-order accuracy near the boundary while using the five-point scheme. A von Neumann stability analysis of the three-time-level scheme for one-dimensional and twodimensional systems of shallow water equations shows that the scheme is unconditionally stable.
Fig. 10. Contours of phase lag of propagation factor for k1 in the Crank–Nicolson scheme for a Courant number of 5.
2362
S. Sankaranarayanan, M.L. Spaulding / Ocean Engineering 30 (2003) 2343–2362
A wave deformation analysis of the two-time-level Crank–Nicolson scheme for one-dimensional and two-dimensional systems of shallow water equations shows that the scheme is non-dissipative, independent of the Courant number and grid resolution used. The phase error or the dispersion of the scheme is found to be small for grid resolutions greater than 20 grids per wavelength and Courant numbers less than 1. Acknowledgements The first author wishes to thank the University of Rhode Island (URI) for providing him with a graduate fellowship for a year. References Arakawa, A., Lamb, V.R., 1977. Computational design of the basic dynamical processes of the UCLA General Circulation Model. Methods in Computational Physics 17, 173–265. Casulli, V., Cattani, E., 1994. Stability, accuracy and efficiency of a semi-implicit method for threedimensional shallow water flow. Computers and Mathematics with Applications 27 (4), 99–112. Foreman, M.G.G., 1984. A two-dimensional dispersion analysis of selected methods for solving linearized shallow water equations. Journal of Computational Physics 56, 287–323. Hirsh, R.S., 1975. Higher order accurate difference solution of fluid mechanics problems by a compact differencing technique. Journal of Computational Physics 9 (1), 90–109. Kowalik, Z., Murty, T.S., 1993. Numerical Modeling of Ocean Dynamics. World Scientific, Singapore. Lele, S.K., 1992. Compact finite difference schemes with spectral-like resolution. Journal of Computational Physics 103, 16–42. Leendertse, J.J., 1967. Aspects of computational model for long period water wave propagation, Report RM-5294-PR, The Rand Corp., Santa Monica, CA. Madala, R.V., Piacsek, S.A., 1977. A semi implicit numerical model for barocline oceans. Journal of Computational Physics 23, 167–178. Messinger, F., Arakawa, A., 1976. Numerical Methods Used in Atmospheric Models, vol. I. GARP Publication Series, No. 17. WMO–ICSU Joint Organizing Committee. Muin, M., Spaulding, M.L., 1996. Two-dimensional boundary-fitted circulation model in spherical coordinates. Journal of Hydraulic Engineering, ASCE 122 (9), 512–521. Song, Y., Tang, T., 1993. Dispersion and group velocity in numerical schemes for three-dimensional hydrodynamic equations. Journal of Computational Physics 105, 72–82. Swanson, J.C., 1986. A three-dimensional numerical model system of coastal circulation and water quality. Ph.D. Thesis, University of Rhode Island, Kingston, RI.