Error Dynamics and Its Use in Developing Optimal

Space-Time Discretization Schemes for Linearized

Rotating Shallow Water Equations
Swagata Bhaumik, Tapan K. Sengupta, Manoj K. Rajpoot*
Department of Aerospace Engineering, Indian Institute of Technology Kanpur, Kanpur
*Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur
208016, India

Abstract
In the present study, we have derived the numerical error propagation equation for
the dispersive two-dimensional linearized rotating shallow water equations (LR-
SWE). In deriving this error equation for LRSWE, we have followed an analysis
similar to that in Sengupta et al. [JCP, 226, 1211-1218 (2007)] for non-dispersive
1D convection equation. Here, also we show that for LRSWE, error and signal
do not follow the same dynamics − contrary to von Neumann analysis for linear
dynamical equations. Thus, the present study is an exercise in showing the error
dynamics for a dispersive multi-dimensional system, which has not been reported
before. We specifically show that the numerical forcing to the error equation is
due to numerical dissipation, phase and dispersion errors, as was also the case for
1D convection equation. Present error equation quantifies the error for various
collocated/staggered grids (A- to E-grid) given in Mesinger & Arakawa [GARP
publication, 17, 43-64 (1976)]. Such an error equation also enables one to opti-
mize space and time discretization schemes together for high accuracy comput-
ing, as is done here for four-stage, fourth order Runge-Kutta scheme coupled with
staggered compact scheme on C-grid.
Keywords: Linearized rotating shallow water equations (LRSWE), Compact
schemes, Runge-Kutta time integration method, Error dynamics, Dispersion
relation preserving (DRP) scheme

Preprint submitted to Applied Mathematics and Computation June 25, 2012

1. Introduction
The von Neumann analysis is often used for error analysis of linear dynamical
systems with constant coefficients [1, 2]. However, it has been shown in [3], that
this is incorrect for the one-dimensional convection equation. For high accuracy
computing, it is imperative that a correct error equation to be framed, as shown in
[3] for a non-dispersive system. It has been noted that for a homogeneous gov-
erning equation an equivalent numerical forcing is created due to dispersion and
dissipation errors as indicated later in Eq. (9). In this equation, the last term on
right hand side is due to dissipation, while the first two terms are due to numeri-
cal dispersion. We note that this is for a physically non-dispersive system. This
clearly indicates the necessity for developing appropriate error equation for indi-
vidual systems. In the present work, we report the error dynamics of a physically
dispersive system, namely linearized rotating shallow water equations (LRSWE).
Dispersion error depends solely on the numerical dispersion relation vis-à-
vis physical dispersion relation. Additionally, a major motivation of the present
work is in identifying the sources of error for this dispersive system as applied to
various grid arrangements. The latter aspect is important as the use of collocated
and staggered grid can create different error sources.
Role of grid in deciding error has been investigated in [4] via a semi-discrete
analysis. The authors proposed five different grid configurations to be used for
shallow water equations in primitive variable and recommended a staggered grid
(C-grid) using second order central difference (CD2 ) spatial discretization scheme.
Here, we will refer to these grids as Arakawa’s grids. In [5], shallow water equa-
tions were studied using a vorticity and divergence of velocity formulation on
a collocated grid and reported superior results. Thus, it appears that staggering
does not always lead to a better performance irrespective of formulation and dis-
cretization methods. Present study incorporates all the aspects, simultaneously for
LRSWE.
It is equally important to note that the numerical dispersion is correctly ob-
tained by considering space and time discretization together [3, 6, 7] and not by
semi-discrete analysis [4, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Such
semi-discrete analysis has lead to observations [21] that “through Fourier analysis,
one can calculate the phase and amplitude error of a given method as a function
of the wavenumber. However, this information can be difficult to interpret”. One
of the reason for such observations is due to the inability to identify the correct
numerical dispersion relation. In the literature, only in [3, 7, 22] correct numerical
dispersion relation has been used for 1D convection equation. In [22], only the

2
numerical phase speed has been identified as the parameter for error propagation.
It is now known [19, 23, 24] that the group velocity is the more relevant propa-
gation property for error evolution. A correct numerical dispersion relation and
associated error propagation equation has only been presented in [3].
It is possible to analyse numerical methods, by recasting the discrete equation
into an equivalent differential equation. This has been termed as modified equation
approach. In the modified equation approach [25, 26, 27, 28, 29, 30, 31], one
looks at the equivalent differential equation in the physical plane via semi-discrete
and full discretization analysis. In [30], this has been classified as Π and Γ-form,
respectively. The analysis in [3] is a complementary Γ-form approach developed
in the spectral plane and not based on the leading truncation error.
In section 2, we present a brief introduction related to numerical dispersion
relation, its consequences for dispersion relation preserving (DRP) analysis and
error dynamics for 1D convection equation, as reported in [3]. In section 3, we
introduce the numerical dispersion relation for LRSWE and define various nu-
merical parameters, as in [32]. We have derived the error propagation equation
in section 4, for dispersive LRSWE for different Arakawa’s grid shown in Fig. 1.
It is established that similar to non-dispersive case, here also, error and signal do
not follow the same dynamics. Finally, we used this error propagation equation
to design a DRP scheme on one of the Arakawa’s grid (C-grid), using four-stage,
fourth order Runge-Kutta (RK4 ) time integration method with staggered compact
scheme given in [33] and its optimized version reported in [32]. Present opti-
mized Runge-Kutta scheme is given in section 5. Moreover, to demonstrate the
efficiency of optimized scheme, we have computed the numerical error and its 2D
FFT using present optimized RK method and compare the results with classical
RK4 time integration method. Furthermore, we have also computed the numerical
solution of LRSWE by using optimized and classical schemes.

2. Dispersion Relation and Its Consequences for 1D Convection Equation

Here, we provide a brief introduction on numerical dispersion relation and its
consequences related to dispersion relation preservation (DRP) analysis [3]. For
this we consider first the model equation as a linear non-dispersive 1D convection
equation given by
∂u ∂u
+c = 0, c > 0 (1)
∂t ∂x
If the unknown u is represented by its bi-lateral Fourier-Laplace R transform at
the jth node of a discrete uniform mesh with width ∆x as, u(x j , t) = U(k, t) eikx j dk,

3
 0
then one can represent the spatial derivative u j by
h 0i Z
uj = ikeq Ueikx j dk (2)
numerical

where, keq denote the equivalent wavenumber for the first order spatial derivative
[23].
The first order spatial derivative can also be represented as [23], u j = ∆x
0 1 PN
l=1 C jl ul ,
where N is the total number of grid points and [C]-matrix can be obtained for dif-
ferent methods, as given in [34, 35] for some representative discretization schemes
in finite difference, finite volume and finite element methods. Using spectral no-
tations, this derivative can be written as [34],
Z N
0 1 X
uj = C jl Ueik(xl −x j ) eikx j dk (3)
∆x l=1

From Eqs. (1) and (3), denoting the quantity ( c∆t ∆x

) as the CFL (Courant-Friedrichs-
Lewy) number (Nc ), one can define a nodal numerical amplification factor as,
G j = U j (k, t(n+1) )/U j (k, t(n) ). For the classical four-stage, fourth order Runge-
Kutta (RK4 ) scheme nodal numerical amplification factor (G j ) is given as [7, 24],

A2j A3j A4j

Gj = 1 − Aj + − + (4)
2 6 24
where, A j = Nc l=1
PN
C jl eik(xl −x j ) .
If the initial condition for Eq. (1) is given as
Z
u(x j , t = 0) = A0 (k) eikx j dk , (5)

then one can write the general solution at any arbitrary time by
Z
uj =
n
A0 (k) (|G j |)n ei(kx j −nβ j ) dk (6)

where, A0 (k) denotes the Fourier transform of the input spectrum and |G j |2 =
[(G j )2real + (G j )2imag ] and tan β j = −[(G j )imag /(G j )real ], with (G j )real and (G j )imag are
the real and imaginary parts of G j , respectively. The numerical phase speed (cN )
can be obtained from nβ j = kcN t. Although, c is same for all k, this shows that
cN is k-dependent. Hence, the numerical solution is dispersive, despite the non-
dispersive nature of Eq. (1). Implications of this difference are very important for
error, as shown in [3].

4
 The physical dispersion relation for Eq. (1) is given by ω = kc, where, ω is
the physical circular frequency. DRP by definition implies consideration of space-
time discretization together for unsteady problems. Despite this, in many earlier
attempts DRP schemes were designed by minimizing leading truncation error [36,
37] for spatial discretization only. Therefore, if one consider the space and time
discretization independently to obtain keq and numerical circular frequency (ωN )
separately, as in [23, 34, 38], then the numerical dispersion relation is erroneously
taken as
ωN = c keq (7)
with, ωN solely depends on spatial discretization (keq ). However, If one consider
space and time discretization together as in [3, 7], then the correct numerical dis-
persion relation is given by
ωN = k cN (8)
where cN is the numerical phase speed and k is the independent variable.
These two different approaches to obtain the numerical dispersion relation
have profound consequences as explained in [3]. First, if one considers the space
and time discretization independent to each other, the numerical dispersion rela-
tion Eq. (7) leads to numerical group velocity given by
dωN dkeq
=c
dk dk
Instead, if one considers space and time discretization together, then the nu-
merical group velocity using correct numerical dispersion relation Eq. (8) is given
as
dωN dcN
= cN + k
dk dk
This is the correct expression for the numerical group velocity. In this case, the
non-dimensional numerical phase speed and group velocity at the jth node for
Eq. (1) was obtained in [3, 7] as
βj
" #
c 
N VgN 1 dβ j
= ; =
c j ω∆t c j Nc d(k∆x)
Using the correct numerical dispersion relation, the propagation equation for nu-
merical error [e = u(x, t) − uN ] for Eq. (1) was obtained in [3] as
∂e ∂e cN  ∂uN
! Z
VgN − cN 
 Z 
0
+ c = −c 1 − − ik0 A0 [|G|]t/∆t eik (x−cN t) dk0 dk
∂t ∂x c ∂x k

5
 Z
ln |G|
− A0 [|G|]t/∆t eik(x−cN t) dk (9)
∆t
This expression is independent of any specific numerical discretization scheme
and is generic to all discretization methods, representing error evolution of nu-
merical processes in contributing to overall error. The right hand side of Eq. (9)
consists of three terms: The first term is due to the phase error, the second term is
due to the dispersion error and the last term is due to numerical dissipation error.
For neutrally stable schemes (|G| = 1), there is no contribution from the last term
and if there are no phase and dispersion errors, i.e., the numerical group velocity
(VgN ) and the numerical phase speed cN must be identical to physical values, then
the first two terms on the right hand side of Eq. (9) will be identically zero. The
von Neumann analysis [1, 2] for linear partial differential equations suggests that
error and signal follow same dynamics and is therefore valid only for neutrally
stable algorithm with no phase and dispersion errors in the limit ∆x, ∆t → 0.
However, such conditions can only prevail for very limited ranges of vanishingly
small space and time steps. Thus, von Neumann analysis is not correct for solving
linear PDEs with finite ∆x and ∆t. It is necessary that error propagation equation
has to be developed for every equation types and this is done here for LRSWE.

3. Numerical Properties of Linearized Rotating Shallow Water Equations

In the Cartesian coordinates, LRSWE is given as the set of following equations
[39, 40, 41],
∂u ∂η
− f v + g =0 (10)
∂t ∂x
∂v ∂η
+ f u + g =0 (11)
∂t ∂y
∂η  ∂u ∂v 
+ H + =0 (12)
∂t ∂x ∂y
where, u and v are the velocity components in x- and y-directions, respectively;
η denotes the displacement of the free surface from constant mean depth H in
z-direction; f is the coriolis parameter; g is the acceleration due to gravity. First
order equations, Eqs. (10) - (12), represents a dispersive model system which can
be written as a single equation for any one of the dependent variables, η, u and v.
Here, a single equation in η is considered, which is given by
∂ ∂
" #
+ f − gH∇ η = 0
2 2
(13)
∂t ∂t2

6
where ∇2 denotes the Laplacian operator in the (x, y)-plane.
Here, we provide a brief introduction of dispersion analysis of LRSWE as
reported in [32]. The dispersion relation for Eqs. (10) - (12) or Eq. (13) [32, 42]
is given by
ω ω2 − c2 |K|
 ~ 2 + f 2
 = 0 (14)
Roots of this equation are given as
q
ω1 = 0 ; ω2,3 = ± f 2 + c2 |K| ~2
√
where, K~ = îk x + ĵky and c = gH, with k x and ky as wavenumber components in
x- and y-directions, respectively. Here, ω1 represents the geostrophic mode, while
ω2,3 represent the inertia-gravity modes [41].
Physical group velocity components in the x- and y-directions for inertia-
gravity modes can be obtained from ω2,3 as [40, 42],
∂ω2,3 k x c2 ∂ω2,3 ky c2
(Vgx )2,3 = =±q ; (Vgy )2,3 = =±q (15)
∂k x ~ ∂ky ~
f + c |K|
2 2 2 f + c |K|
2 2 2

From Eq. (15), exact absolute group velocity (Vg ) for inertia-gravity modes
are given as q
q
Vg = Vgx ~
2 + V 2 = (c | K|)/
gy
2 ~2
f 2 + c2 |K|

which makes an angle, θex = tan−1 (Vgy /Vgx ) with the x-axis.
Similarly, the physical phase speed components in x- and y-directions are ob-
tained as q q
ω2,3
2 2 ~
f + c |K| 2
ω2,3 ~2
f 2 + c2 |K|
(cex )2,3 = =± ; (cey )2,3 = =± (16)
kx kx ky ky
Next, numerical equivalent of these quantities for any general numerical spa-
tial discretization schemes are obtained for any two level explicit time integration
methods. Here, the numerical properties of LRSWE are obtained by considering
space-time discretization together, as given in [32]. One can write Eqs. (10) - (12)
in matrix notation by denoting Z = [u, v, η]T as

∂Z ∂Z ∂Z
+ [A]Z + [B] + [C] =0 (17)
∂t ∂x ∂y
where,

7
      
 0 − f 0   0 0 g   0 0 0 
[A] =  f 0 0  , [B] =  0 0 0  and [C] =  0 0 g 
     
0 0 0 H00 0H0
     

Let, Ẑ = [Û, V̂, Ê]T , where, Û, V̂ and Ê denote the Fourier-Laplace trans-
R Rη, respectively. Using Fourier-Laplace transform one can repre-
forms of u, v and
sent Z as, Z = Ẑ ei (kx x+ky y) dk x dky . The first order spatial derivative using any
discretization scheme can be represented as

"
∂Z
= i(k x )eq Ẑei (kx x+ky y) dk x dky
∂x
"
∂Z
= i(ky )eq Ẑei (kx x+ky y) dk x dky (18)
∂y
where (k x )eq and (ky )eq are numerical equivalent wavenumbers for obtaining first
order derivatives in x- and y-directions, respectively. Hence, upon using numerical
spatial discretization schemes in Eq. (17) and using Eq. (18), we get the following

∂Ẑ
+ [D] Ẑ = 0 (19)
∂t
where,  
 0 − feq i(k x )eq g 
[D] =  feq 0 i(ky )eq g 
 
i(k x )eq H i(ky )eq H 0
 

with feq = f T x (k x ∆x) T y (ky ∆y), where T x and T y are the spectral transfer func-
tions associated with the chosen midpoint interpolation scheme in staggered grids
for all variables in x- and y-directions, respectively. For collocated grids T x =
T y = 1, as it does not require any midpoint interpolation. Here, ∆x and ∆y are
uniform grid spacings in x- and y-directions, respectively. The numerical prop-
erties of any space-time discretization method can be evaluated for the system in
Eq. (17). For any two time-level temporal discretization method, one can relate
the values of Ẑ at (n + 1)th time-level to nth time-level values as

Ẑ n+1 = [PE ] Ẑ n (20)

where [PE ] is the evolution matrix. For classical RK4 time integration method, the
2 (∆t)3
evolution matrix [PE ] is obtained as [PE ] = I − (∆t) [D] + (∆t)
2
[D] 2
− 6
[D]3 +

8
(∆t)4
24
[D]4 . It can be shown that the modal amplification factors, G = [G1 , G2 , G3 ]T
are the eigenvalues of [PE ]. For any general spatial discretization scheme appli-
cable to uniform grids, it can be shown [32] that
1 1
(k x )eq = ζ(k x ∆x) and (ky )eq = ψ(ky ∆y)
∆x ∆y
where, ζ and ψ are functions of k x ∆x and ky ∆y, respectively.
Let us denote a ≡ Ncx ζ(k x ∆x) and b ≡ Ncy ψ(ky ∆y), where Ncx = c∆t ∆x
and
Ncy = c∆t
∆y
are the CFL numbers based on grid spacing in x- and y-directions,
respectively. Using symbolic toolbox operations, eigenvalues of [PE ] for classical
RK4 time integration scheme are calculated as

G1 = 1, G2, 3 = δ ∓ i  (21)

where, δ = 1 − 12 γ2 + 241 γ4 ;  = γ − 16 γ3 with γ = (a2 + b2 + p2 ) and p = feq ∆t.

p

Here, G1 corresponds to the geostrophic mode and G2, 3 correspond to the inertia-
gravity modes.
Numerical circular frequency, ωN , for a particular mode can be calculated as
1
ωN = − tan−1 Gimag /Greal


∆t
where, Gimag and Greal are the imaginary and real parts of the modal amplification
factor, respectively. From Eq. (21), these can be obtained as

1 −1 
 
(ωN )1 = 0, (ωN )2,3 = ± tan (22)
∆t δ
Numerical group velocity components in x- and y-directions, i.e., VgN x and VgNy
for inertia-gravity modes are obtained as
∂ωN ∂ωN
VgN x = ; VgNy = (23)
∂k x ∂ky
and the numerical phase speeds in x- and y-directions, i.e., cN x and cNy , can be
obtained from
ωN ωN
cN x = ; cNy = (24)
kx ky

9
 Using Eqs. (22) and (23), one can obtain numerical group velocity components
in x- and y-directions for inertia-gravity modes as
  1 dξ  dζ ∆x dT x 
VgN x =±c a + p Ty
2,3 γ dγ d(k x ∆x) Lr d(k x ∆x)
  1 dξ  dψ ∆y dT y 
VgNy =±c b + p Tx (25)
2,3 γ dγ d(ky ∆y) Lr d(ky ∆y)
  √
where, ξ = tan−1 Gimag /Greal and Lr (= gH/ f ) is the Rossby radius [41].
From Eq. (25), one calculates absolute numerical group velocity (VgN ) and its
angle of propagation (θN ) for inertia-gravity modes as
" #
−1 (VgNy )2,3
q
VgN 2,3 = (VgN x )2,3 + (VgNy )2,3 ; (θN )2,3 = tan
2 2


(26)
(VgN x )2,3
Similarly, from Eqs. (22) and (24), numerical phase speed components in x- and
y-directions, can be obtained as
ξ ξ
(cN x )2,3 = ± c ; (cNy )2,3 = ± c (27)
Ncx k x ∆x Ncy ky ∆y
Next, we obtain the governing equation for the error propagation of LRSWE
on different grid arrangements.

4. Error Dynamics for LRSWE

Here, we obtain the error propagation equation for dispersive LRSWE (Eqs. (10)
- (12)). Let {ZNn } = {unN , vnN , ηnN }T and {ẐNn } = {Û Nn , V̂Nn , Ê Nn }T denote the numerical
solution of LRSWE at the nth time level for any two time-level temporal discretiza-
tion scheme and its corresponding bilateral Fourier-Laplace transforms, respec-
tively. Then using Eq (20), one can write

{ẐNn } = [PE ]n {Ẑ 0 } (28)

where, {Ẑ0 } = {Û0 , V̂0 , Ê0 }T is the bilateral Fourier-Laplace transform of the pre-
scribed initial condition {Z0 } = {u0 , v0 , η0 }T . Let, el ’s be the right eigenvectors of
PE with corresponding eigenvalues Gl for l = 1, 2, 3. Here, the subscript l denotes
three different modes. Then one can express Ẑ0 , as
3
X
{Ẑ0 } = κl {el } (29)
l=1

10
where, κl ’s are constants those can be evaluated. Using Eq. (28) - (29) one can
write the general expression for {ẐNn } as
3
X
{ẐNn } = κl (Gl )n {el } (30)
l=1

From Eq. (29) one can also define spectral weights al , bl and cl ’s as al =
κl e1l /Û0 ; bl = κl e2l /V̂0 and cl = κl e3l /Ê0 with l = 1, 2, 3, where eil is the ith
element of lth eigenvector {el }. One can also notice that 3l=1 al = 3l=1 bl =
P P

l=1 cl = 1. Hence, one obtains

P3

3
X 3
X 3
X
Û Nn = al (Gl )n Û0 , V̂Nn = bl (Gl )n V̂0 and Ê Nn = cl (Gl )n Ê0
l=1 l=1 l=1

These relations enable one to write the general expression for {ZNn } = {unN , vnN , ηnN }T ,
as
" X 3
uN =
n
al (Gl )n Û0 ei(kx x+ky y) dk x dky
l=1
" X
3
vnN = bl (Gl )n V̂0 ei(kx x+ky y) dk x dky
l=1
" X
3
ηnN = cl (Gl )n Ê0 ei(kx x+ky y) dk x dky (31)
l=1

Since Gl ’s are complex quantities for inertia-gravity modes, one can write Gnl
as
Gnl = |Gl |n e−i ωNl (n∆t) = |Gl |t/∆t e−i ωNl t
Using this expression one can construct the following relation

∂Gnl
= Gnl Rl (32)
∂t
 (G ) 
where, Rl = (θl − i ωNl ), with θl = ln∆t|Gl | , ωNl = − ∆t1 tan−1 (Gll )imag
real
.
Combining Eqs. (10) - (12) with Eqs. (31) and (32), one obtains
" X 3
∂ηN  ∂u
N ∂vN  h
+ H + = (Gl )n cl Rl Ê0 + H ik x al Û0
∂t ∂x ∂y l=1

i
+iky bl V̂0 ei(kx x+ky y) dk x dky (33)

11
 " X
3
∂uN ∂ηN h
− f vN + g = (Gl )n al Rl Û0 − f bl V̂0
∂t ∂x l=1
i
+ik x gcl Ê0 ei(kx x+ky y) dk x dky (34)
" X
3
∂vN ∂ηN h
+ f uN + g = (Gl )n bl Rl V̂0 + f al Û0
∂t ∂y l=1
i
+iky gcl Ê0 ei(kx x+ky y) dk x dky (35)
By using Eqs. (33) - (35) and performing the following operations
∂ ∂ ∂ 
(33) − H (34) + (35)
∂t ∂x ∂y
and after simplifying, one gets
" X
3
∂2 ηN  ∂v
N ∂uN  h
~ 2

− gH ∇ ηN + H f
2
− = (Gl )n cl Ê0 R2l + gH|K|
∂t2 ∂x ∂y l=1

i
+ H f ik x bl V̂0 − iky al Û0 ei(kx x+ky y) dk x dky (36)

Once again using Eqs. (34) - (36) and performing following operations
∂ ∂ ∂ 
(36) − H (35) + (34)
∂t ∂x ∂y
after simplifying, one finally obtains
" X
3
∂  ∂2 ηN 
2 ∂uN
 ∂vN  h
~ 2

− gH ∇ ηN − H f
2
+ = (Gl )n cl Ê0 R2l + gH|K|
∂t ∂t2 ∂x ∂y l=1

i
+ H f 2 ik x al Û0 − iky bl V̂0 ei(kx x+ky y) dk x dky (37)

Using Eq. (33), one can simplify Eq. (37) as

∂  ∂2 ηN 
− gH ∇2
η N + f 2
η N
∂t ∂t2
" X 3
~ 2 + f 2 dk x dky
 
= (Gl )n cl Rl Ê0 R2l + gH |K| (38)
l=1

12
 Let us define an error as, eη = ηexact − ηN , and write the error propagation
equation for the variable η using Eq. (38) as

3 "
∂ ∂2 eη
" # X
+ f eη − gH∇ eη = −
2 2
cl (Gl )n T l Eˆ0 ei(kx x+ky y) dk x dky (39)
∂t ∂t2 l=1

where, T l = Rl R2l + gH |K| ~ 2 + f 2
for l = 1, 2, 3.

Equation (39) represents the error propagation dynamics for dispersive LR-
SWE, for any general numerical schemes and grid employed for the solution
method. One obtains identical error propagation equation for the other two phys-
ical variables u and v. Similar to 1D convection equation [3], here also, one notes
that the error and the signal do not follow the same dynamics, despite the fact
that the governing equations are homogeneous and linear. One notes that for the
geostrophic mode, G1 = 1 and hence θ1 = ωN1 = 0. Therefore, error will only be
forced due to inertia-gravity modes (G2 and G3 ).
Also, for inertia-gravity modes (l=2, 3), T l ’s represent complex conjugate pair
and therefore, one needs to evaluate the modulus of T l for l = 2, while |T 2 | = |T 3 |.
For the first inertia-gravity mode, T 2 is given by

T 2 = R2 R22 + ω22 = R2 R22 + c2 (k2x + ky2 ) + f 2

   
(40)
h (G ) i
where, R2 = ∆t1 (ln |G2 | − iβ2 ), with β2 = − tan−1 (G22 )imag
real
= ωN2 ∆t and c2 = gH.
One can represent T 2 in terms of |G2 | and ωN2 as

1
T2 = T∗ (41)
(∆t)3 2
with
h
3
2 i
T 2∗ = ln |G2 | − (∆t)2 3ω2N2 − ω22 ln |G2 |−3i ∆t ωN2 ln |G2 | + i(∆t)3 ωN2 (ω2N2 −ω22 )

In the error equation (Eq. (39)), T 2 forces the error due to dissipation, phase mis-
match and dispersion error. Here, errors are due to terms involving ln |G2 | and
(ωN2 − ω2 ). In T 2 , the terms involving ln |G2 | represent forcing by numerical am-
plification factor for |G2 | , 1, i.e., for schemes which are not neutrally stable.
Once again, this is in conformity with the observation in [3] for 1D convection
equation. Additionally, the fourth term in T 2∗ involving (ωN2 − ω2 ) contributes to
error via mismatch of physical and numerical dispersion relation. Consequently,

13
if the numerical scheme is neutrally stable with no phase/dispersion error, then
T 2 ≡ 0. Therefore, one can conclude that for DRP schemes error will be minimal,
when it is neutrally stable with negligible phase and dispersion errors.
To quantify error on different grids, one rewrites T 2∗ in terms of non-dimensional
wavenumbers (k x ∆ and ky ∆) and CFL numbers (Nc x = Ncy = Nc) with uniform
grid spacing (∆x = ∆y = ∆) as

T 2∗ = X + iY (42)

where,
h
2 n o i
X = ln |G2 | ln |G2 | − 3(β2 )2 + (Ncx )2 (k x ∆)2 + (ky ∆)2 + ( f ∆t)2
and h n o i

2
Y = −β2 3 ln |G2 | − (β2 )2 + (Ncx )2 (k x ∆)2 + (ky ∆)2 + ( f ∆t)2
In the present case we report the
√ result for the case for which grid spacing is
less than the Rossby radius (Lr = (gh)/ f ). One can note that for this case C-grid
performs well as compared to other grids [5]. Thus, in the present case we use the
developed error equation to design optimal DRP schemes on C-grid, which can be
extended to any other grid in a similar fashion. Next, we report the minimization
of the numerical error on Arakawa’s C-grid for first inertia-gravity mode.

5. Optimized Space-Time Discretization Schemes for Least Error for LR-

SWE
To minimize error for LRSWE, we consider an initial conservative spectrum
of η0 given by, Eˆ0 = 1, which corresponds to seeking an impulse response for
the error dynamics. Thus, the right hand side of Eq. (39), which defines the error
forcing is strictly dependent upon ∆t and T 2∗ . Thus, seeking minimum error for
LRSWE is equivalent to minimizing the objective function

F = F1 + F2 (43)

with, Z σπ Z σπ
F1 = |X| dk x dky
0 0
and Z σπ Z σπ
F2 = |Y| dk x dky
0 0

14
where, X and Y are as given in Eq. (42) and 0 ≤ σ ≤ 1. In the present case we
have used σ = 1/2.
In the present case, we use optimized staggered compact scheme (OSCS)
for spatial discretization with optimized interpolation scheme [32] and optimize
fourth order, four-stage Runge-Kutta (RK4 ) time integration scheme as the base-
line time integration method. Thus, to minimize the numerical error, one can
construct a perturbed RK4 scheme for which the evolution matrix in Eq. (20) is
given by
(∆t)2
PE = I − (∆t)[D] + [D]2 − α2 (∆t)3 [D]3 + α1 (∆t)4 [D]4 ,
2
where, α1 and α2 are the free parameters leaving the resultant Runge-Kutta scheme
as second order accurate. Thus, the search space for optima is (α1 , α2 )-plane.
Once we have the optimal values for free parameters α1 and α2 , corresponding
algorithm for optimized Runge-Kutta (ORK) scheme in a low-storage form for
∂u
∂t
= L(u) is given as
α 
1
u(1) = un + ∆t L(un )
α2
u(2) = un + (2α2 ) ∆t L(u(1) )
1
u =u +
(3) n
∆t L(u(2) )
2
un+1 = un + ∆t L(u(3) ) (44)
In this process, first we have considered the variation of |F| in the (α1 , α2 )-
plane as shown in Fig. 2(i). One notes from Fig. 2(i) that the minimum value
of error is almost independent of the variation of parameter α2 in the admissible
range. Therefore, a single parameter optimization is performed by fixing α2 = 1/6
(its classical value), and locate an optimal value for α1 , as shown in Fig. 2(ii).
Notionally, this is a third order accurate method. In Fig. 2, results are shown
plotted for a value of Nc = 0.2. Variation of the optimal value for parameter
α1 with CFL number is shown in Fig. 3. One can notice from Fig. 3, that the
variation of parameter α1 with higher value of Nc is almost negligible.
In Fig. 4, numerical properties of the present optimized scheme along with its
classical counterpart are shown for LRSWE when using the C-grid. It is noted that
after optimization for Nc = 0.2, the area under |G2 | = 1 contour increases dras-
tically, thereby improving its neutral stability property tremendously. It should
be also noted that, such optimization procedure does not significantly alter the
numerical phase speed and group velocity properties.

15
6. Dynamics of Numerical Error and Numerical Solution of LRSWE using
Classical/Optimized Schemes
Here, we demonstrate the quantitative behavior of numerical error with respect
to the particular initial condition on A- and C-grids. For this purpose, we have
chosen the following initial condition

u(x, y, t = 0) = 0
v(x, y, t = 0) = 0
 (x−x0 )2 
+(y−y0 )2
η(x, y, t = 0) = e
−α 15 sin(k x x + ky y) (45)

with α = 0.05.
The numerical error e (=ηexact − ηnum ) evolution for the initial condition given
by Eq. (45) which is centered at (ko ∆x, ko ∆y) = (1.2, 1.2) is studied. Numerical
error is computed by considering the difference between the exact and numerical
solution at particular time instants. For the numerical solution of LRSWE, we
have considered the product of Coriolis parameter f and a time step ∆t as, f ∆t =
1.2 × 10−2 with grid spacing in both the directions ∆x = ∆y = 0.3 in the domain
of size [0, 154] × [0, 154]. The corresponding exact solution of the LRSWE with
respect to the same initial condition (Eq. (45)) is computed using 2D FFT.
Numerical error at different time instants for Nc = 0.2 on C-grid using clas-
sical RK4 method along with 6th order staggered compact scheme (SCS) [33] are
shown in the left column of Fig. 5 at indicated time instants. Moreover, numeri-
cal error on C-grid using present optimized RK method with optimized staggered
compact scheme (OSCS) [32] at different time instants are shown in the right col-
umn of Fig.5. Corresponding 2D FFT of numerical error shown in Fig. 6 for the
C-grid using optimized and classical schemes. One can notice form the maximum
and minimum values marked in each sub-frames in Figs. 5-6, that numerical er-
ror on C-grid using the optimized Runge-Kutta method with optimized staggered
compact scheme (OSCS) scheme [32] is lesser as compared to that incurred by
classical RK4 method with staggered compact scheme given in [33].
To further assess the enhanced efficiency of the optimized scheme, next we
consider the propagation of 2D wave-packet following LRSWE. In this process, to
demonstrate the superiority of optimized scheme on C-grid in the physical plane,
we have solved the LRSWE by considering propagation of a 2D wave-packet
given by given by Eq. (45) with α = 0.05 in a domain of size [0, 360] × [0, 360]
using uniform grid spacing with ∆x = ∆y = 0.3. Wave-packet propagation angle,

16
θ = 45o , is used corresponding to the aspect ratio of λ = ∆y/∆x = 1. Moreover, we
have chosen product of Coriolis parameter f and a time step ∆t as, f ∆t = 1.2×10−2
and Nc x = Ncy = 0.2. This packet is centered at (k0 ∆x, k0 ∆y) = (1.4, 1.4).
We have computed the numerical solutions of LRSWE at different time in-
stants following the 2D wave-packet given by Eq. (45) using non-optimal Runge-
Kutta time integration method along with optimized staggered compact scheme
(OSCS) [32] for spatial discretization, which correspond to the α1 = 0.03 and
α2 = 1/6. We have also computed the numerical solutions using classical Runge-
Kutta scheme along with staggered compact scheme with staggered interpolation
technique reported in [33] and present optimized Runge-Kutta method with OSCS
spatial discretization scheme [32], as shown in Fig. 7. Top frame of Fig. 7 shows
the initial condition. Results corresponding to non-optimal Runge-Kutta method
are shown in frames (i) and (ii) of Fig. 7. Frames (iii) and (iv) are due to classical
Runge-Kutta method and results with optimized Runge-Kutta method are shown
in frames (v) and (vi) of Fig. 7. One can conclude from Fig. 7 that due to lesser
dissipation error for the optimized scheme (as noted from Fig. 4), attenuation in
the amplitude of packet is lesser than the other two cases, as noted from the max-
imum values. Numerical results in Figs. 5-7, clearly indicate the efficiency of
optimized schemes.

7. Conclusions
In the present work, we have derived the error propagation equation for the
case of 2D dispersive LRSWE. Similar to 1D convection equation [3], here also
it is shown that error and signal does not follow the same dynamics and once
again this is contrary to assumption made by von Neumann for the case of linear
PDEs. Developed error propagation equation for LRSWE is used to design true
DRP low-storage RK schemes. This is performed here on Arakawa’s C-grid. The
present optimization strategy can be extended to the other Arakawa’s grids as well
as other prevailing schemes used to solve LRSWE [43, 44] also. In Figs. 5-6, we
have also plotted the numerical error with respect to the initial condition given by
Eq. (45) and its 2D FFT at different time instants for C-grid with Nc = 0.2. From
these results we have noted that C-grid using present optimized RK method with
OSCS scheme provides least error as compared to the other cases reported here.
Furthermore, to show the effectiveness of the optimized scheme, we have con-
sidered the propagation of a 2D wave-packet following LRSWE. Numerical re-
sults are shown at t = 10.8 and t = 14 in the Fig. 7 for non-optimal, classical
and optimized schemes. By looking at the maximum and minimum values of the

17
amplitude of wave-packet, as marked in the respective frames, one can notice that
present scheme performs better than the other schemes considered in Fig. 7. This
is in accordance with the numerical properties of optimized scheme, which has
significantly larger neutrally stable region as compared to the classical Runge-
Kutta scheme.
The final form of the error equation for numerical error, eη = ηexact − ηN , for
free surface disturbance η is given as

3 "
∂ ∂2 eη
" # X
+ f eη − gH∇ eη = −
2 2
cl (Gl )n T l Eˆ0 ei(kx x+ky y) dk x dky
∂t ∂t2 l=1

where, T l = Rl R2l + gH |K| ~ 2 + f 2
for l = 1, 2, 3 with Rl = (θl − i ωNl ), θl = ln |Gl |

 (G )  ∆t
and ωNl = − ∆t1 tan−1 (Gll )imag
real
.
The final scheme for least numerical error for solving LRSWE (Eqs. (10) -
(12)) on C-grid for temporal and spatial discretizations is given next. The algo-
rithm for optimized Runge-Kutta method is given as

α 
1
u (1)
=u +
n
∆t L(un )
α2
u = u + (2α2 ) ∆t L(u(1) )
(2) n
1
u(3) = un + ∆t L(u(2) )
2
un+1 = un + ∆t L(u(3) )

with α1 = 0.0428799995 and α2 = 1/6.

For the spatial discretization of LRSWE on C-grid, the optimized staggered
compact scheme (OSCS) [32] is given as

0 0 0
u
j+3/2− u j−3/2  u
j+1/2 − u j−1/2

ρ1 u j−1 + u j + ρ1 u j+1 = b1 +a1
3∆ ∆

where a1 = 38 (3 − 2ρ1 ), b1 = 18 (−1 + 22ρ1 ) with ρ1 = 0.18 and ∆ denotes the grid
spacing.
As for staggered girds one also requires mid-point interpolation of the vari-
ables, the optimized interpolation scheme [32] is given as

18
 b2   a 
2

ρ2 û j−1 + û j + ρ2 û j+1 = u j+3/2 + u j−3/2 + u j+1/2 + u j−1/2
2 2

where û j denotes interpolated value of u at jth node with a2 = 18 (9 + 10ρ2 ), b2 =

1
8
(−1 + 6ρ2 ) where ρ2 = 0.32.

References
[1] J. G. Charney, R. Fjørtoft, J. von Neumann, Numerical integration of the
barotropic vorticity equation, Tellus 2 (1950) 237–254.

[2] J. Crank, P. Nicolson, A practical method for numerical evaluation of so-

lutions of partial differential equations of the heat-conduction type, Proc.
Camb. Phil. Soc. 43 (1947) 50–67.

[3] T. K. Sengupta, A. Dipankar, P. Sagaut, Error dynamics: Beyond von Neu-

mann analysis, J. Comput. Phys. 226 (2007) 1211–1218.

[4] F. Mesinger, A. Arakawa, Numerical methods used in atmospheric models,

Vol. 1, GARP Publ. Ser. No. 17, 1976. WMO, Geneva.

[5] D. A. Randall, Geostrophic adjustment and the finite-difference shallow-

water equations, Mon. Wea. Rev. 122 (1994) 1371–1377.

[6] M. K. Rajpoot, T. K. Sengupta, P. K. Dutt, Optimal time advancing disper-

sion relation preserving schemes, J. Comput. Phys. 229 (2010) 3623–3631.

[7] T. K. Sengupta, A. Dipankar, A comparative study of time advancement

methods for solving Navier-Stokes equations, J. Sci. Comput. 21 (2004)
225–250.

[8] J.-M. Beckers, E. Deleersnijder, Stability of a FBTCS scheme applied to the

propagation of shallow-water inertia-gravity waves on various space grids,
J. Comput. Phys. 108 (1993) 95–104.

[9] J. K. Dukowicz, Mesh effects for Rossby waves, J. Comput. Phys. 119
(1995) 188–194.

19
[10] M. G. G. Foreman, A two-dimensional dispersion analysis of a selected
methods for solving the linearized shallow water equations, J. Comput. Phys.
56 (1984) 287–323.

[11] W. G. Gray, D. R. Lynch, On the control of noise in finite element tidal

computations: A semi-implicit approach, Comput. Fluids. 7 (1979) 47–67.

[12] R. F. Henry, Richardson-Sielecki schemes for the shallow-water equations

with application to Kelvin waves, J. Comput. Phys. 41 (1981) 389–406.

[13] D. K. Lilly, A proposed staggered grid system for numerical integration of

dynamic equations, Mon. Wea. Rev. 89 (1961) 59–65.

[14] M. S. F.-Ravinovitz, Computational dispersion properties of horizontal stag-

gered grids for atmospheric and ocean models, Mon. Wea. Rev. 119 (1991)
1624–1639.

[15] A. Sielecki, An energy-conserving difference scheme for the storm surge

equations, Mon. Wea. Rev. 96 (1968) 150–156.

[16] Y. Song, T. Tang, Dispersion and group velocity in numerical schemes for
three-dimensional hydrodynamic equations, J. Comput. Phys. 105 (1993)
72–82.

[17] J. Thuburn, Rossby wave dispersion on the C-grid, Atmos. Sci. Let. 8 (2007)
37–42.

[18] J. Thuburn, Numerical wave propagation on the hexagonal C-grid, J. Com-

put. Phys. 227 (2008) 5836–5858.

[19] R. Vichnevetsky, Energy and group velocity in semi discretizations of hy-

perbolic equations, Math. Comp. in Simulation (1981) 333–343.

[20] R. C. Wajsowicz, Free planetary waves in finite-difference numerical meth-

ods, J. Comput. Phys. 16 (1986) 773–789.

[21] D. W. Zingg, Comparison of high accuracy finite-difference methods for

linear wave-propagation, SIAM J. Sci. Stat. Comput. 22 (2000) 476–502.

[22] G. H. Haltiner, R. T. Williams, Numerical Prediction and Dynamic Meteo-

rology, John Wiley & Sons 2nd Ed., 1980.

20
[23] T. K. Sengupta, G. Ganeriwal, S. De, Analysis of central and upwind com-
pact schemes, J. Comput. Phys. 192 (2003) 677–694.

[24] R. Vichnevetsky, J. B. Bowles, Fourier Analysis of Numerical Approxima-

tions of Hyperbolic Equations, SIAM Stud. Appl. Math., SIAM Philadel-
phia, 1982.

[25] S.-C. Chang, A critical analysis of the modified equation technique of Warm-
ing and Hyett, J. Comput. Phys. 86 (1990) 107–126.

[26] D. F. Griffiths, J. M. Sanz-Serna, On the scope of the method of modified

equations, SIAM J. Sci. Stat. Comput. 7 (1986) 994–1008.

[27] F. F. Grinstein, L. G. Margolin, W. J. Rider, Implicit Large Eddy Simulation:

Computing Turbulent Fluid Dynamics, Cambridge Univ. Press, 2007.

[28] C. W. Hirt, Heuristic stability theory for finite-difference equations, J. Com-

put. Phys. 2 (1968) 339–355.

[29] P. J. Roache, Computational Fluid Dynamics, Hermosa publishers, Albu-

querque, New Mexico, 1972.

[30] Y. I. Shokin, The Method of Differential Approximation, Springer-Verlag,

1983.

[31] R. F. Warming, B. J. Hyett, The modified equation approach to the stability

and accuracy analysis of finite difference methods, J. Comput. Phys. 14
(1974) 159–179.

[32] M. K. Rajpoot, S. Bhaumik, T. K. Sengupta, Solution of linearized rotating

shallow water equations by compact schemes with different grid staggering
strategies, J. Comput. Phys. 231 (2012) 2300–2327.

[33] S. Nagarajan, S. K. Lele, J. H. Ferziger, A robust high-order compact method

for large eddy simulation, J. Comput. Phys. 19 (2003) 392–419.

[34] T. K. Sengupta, Fundamentals of Computational Fluid Dynamics, Univ.

Press, Hyderabad, 2004.

[35] T. K. Sengupta, Y. G. Bhumkar, M. K. Rajpoot, V. K. Suman, S. Saurabh,

Spurious waves in discrete computation of wave phenomena and flow prob-
lems, Appl. Math. Comput. 218 (2012) 9035–9065.

21
[36] C. Bogey, C. Bailly, On the application of explicit spatial filtering to the
variables or fluxes of linear equations, J. Comput. Phys. 225(2) (2007) 1211–
1217.

[37] C. K. W. Tam, J. C. Webb, Dispersion-relation-preserving finite difference

schemes for computational acoustics, J. Comput. Phys. 107 (1993) 262–281.

[38] L. N. Trefethen, Group velocity in finite difference schemes, SIAM Rev. 24

(1982) 113–136.

[39] A. E. Gill, Atmosphere-Ocean Dynamics, International Geophysics Series,

Academic Press, New York, 1982.

[40] P. K. Kundu, I. M. Cohen, Fluid Mechanics, Elsevier, 2008.

[41] J. Pedlosky, Geophysical Fluid Dynamics, Springer Verlag, New York, 1979.

[42] T. K. Sengupta, M. K. Rajpoot, S. Saurabh, V. V. S. N. Vijay, Analysis of

anisotropy of numerical wave solutions by high accuracy finite difference
methods, J. Comput. Phys. 230 (2011) 27–60.

[43] A. F. Shchepetkin, J. C. McWilliams, The regional oceanic modeling sys-

tem (ROMS): a split-explicit, free-surface, topography-following-coordinate
oceanic model, Ocean Modelling 9 (2005) 347–404.

[44] A. F. Shchepetkin, J. C. McWilliams, Computational kernel algorithms for

fine-scale, multi-process, long-term oceanic simulations, Handbook of Nu-
merical Analysis, Vol. XIV: Computational Methods for the Ocean and At-
mosphere Elsevier Science, doi:10.1016/S1570-8659(08)01202-0, 2008. P.
G. Ciarlet, editor, R. Temam and J. Tribbia, guest eds.

22
 u, v, η u, v, η ( η η
u, v, η η
j+1
B-grid u, v
A-grid
u, v
u, v, η u, v, η u, v, η η η η
j
u, v u, v

j-1
u, v, η u, v, η u, v, η η η η
i+1
i-1
∆
i
∆

η u η u η v η v
η η
C-grid v u
v v D-grid u
u
η u η u η η v η v η
v v u u u
v
η u η u η η v η v η
∆ ∆
u, v j+1/2 ϕ, µ, η ϕ, µ, η (
η u, v ϕ, µ, η
j+1/2 j+1
E-grid Z-grid
η u, v
η ϕ, µ, η ϕ, µ, η
j ϕ, µ, η
j
∆

u, v η u, v ϕ, µ, η ϕ, µ, η ϕ, µ, η
j-1
i-1/2 j-1/2
i i+1/2
i-1 ∆ i i+1

Fig. 1 Different types of Arakawa’s grids (A- to E-grid) considered for the finite
difference solution of the LRSWE with ∆ representing the grid spacing of the
variables for which finite difference is taken. Here, Z-grid is due to Randall [5]
with φ denoting the divergence of velocity and µ is the vorticity.

23
 (i)

0.05 0.0002

Error
α2

0.1

0.15 0.00015

0.2

0.08 0.06 0.04

α1
0.02 0

(ii) 0.003

0.002
Error

Classical value
for α1
0.001 Present optimum

0
α 1= 0.0428799995

-0.001
0 0.02 0.04 α1 0.06 0.08

Fig. 2 (i) Two parameter optimization in the (α1 , α2 )-plane for Nc = 0.2; (ii)
Single parameter optimization for LRSWE for Nc = 0.2 with dash-dotted line
represent F1 and the solid line represent F2 , as given by Eq. (43). The Vertical
dashed line indicates the classical value of parameter α1 .

24
 0.0432

0.0428
α1

0.0424

0.042

0.0416

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nc

Fig. 3. Variation of optimal value of parameter α1 with different values of Nc for

the case of single parameter optimization for LRSWE.

25
 Optimized Classical
|G| |G| Classical
3 3

0.999865
2 2 0.999865
ky∆y

ky∆y
1
Min = 0.9986
1 1 Max = 1.0
Min = 0.9990
Max = 1.0 1

1 kx∆x 2 3 1 kx∆x 2 3

(VgN/Vg) (VgN/Vg)
3 3
Min = 1.73×10-16 Min = 1.73×10-16
Max = 1.0063 Max = 1.0063

2 2
ky∆y

ky∆y

1 1

1 kx∆x 2 3 1 kx∆x 2 3

(cNx/cex) (cNx/cex)
3 Min = 0.8597 3 Min = 0.8597
Max = 1.0024 Max = 1.0024

2 2
ky∆y

ky∆y

1 1

1 kx∆x 2 3 1 kx∆x 2 3

Fig. 4 Contour plots for C-grid showing |G|, (VgN /Vg ) & (cN x /cex ) using the op-
timized Runge-Kutta method with optimized staggered compact scheme (OSCS)
[32] (shown in the left frames) and classical Runge-Kutta scheme with staggered
compact scheme reported in [33] (shown in the right frames). Here, the numerical
properties are shown plotted for Nc = 0.2 and f ∆t = 1.2 × 10−2 .

26
 Num. error for RK4-SCS Num. error for optimized
scheme RK-OSCS scheme
min: -0.0023 max: 0.0023 min: -0.0017 max: 0.0017
150 150
t = 0.05 t = 0.05
100 100
Y

Y
50 50

0 0
0 50 100 150 0 50 100 150
X X
min: -0.0438 max: 0.0430 min: -0.0316 max: 0.0322
150 150
t = 1.0 t = 1.0
100 100
Y

50 50

0 0
0 50 100 150 0 50 100 150
X X
min: -0.0813 max: 0.0824 min: -0.0561 max: 0.0550
150 150
t = 2.5 t = 2.5
100 100
Y
Y

50 50

0 0
0 50 100 150 0 50 100 150
X X
min: -0.1074 max: 0.1084 min: -0.0786 max: 0.0781
150 150
t = 4.5 t = 4.5
100 100
Y
Y

50 50

0 0
0 50 100 150 0 50 100 150
X X

Fig. 5 Numerical error on C-grid at indicated time instants using: (i) RK4 time
integration scheme along with 6th order staggered compact scheme (SCS) [33]
(shown in the left column); (ii) present optimized RK method with optimized
staggered compact scheme (OSCS) [32] (shown in the right column) for f ∆t =
1.2 × 10−2 , λ = ∆y/∆x = 1 and Nc = 0.2.
27
 FFT of numerical error for FFT of numerical error for
RK4 - SCS scheme optimized RK - OSCS scheme
min: 0 max: 0.0080 min: 0 max: 0.0056
t = 0.05 t = 0.05
1.5 1.5
ky

ky
1 1

1 kx 1.5 1 kx 1.5

min: 0 max: 0.1564 min: 0 max: 0.1148

t = 1.0 t = 1.0
1.5 1.5
ky

ky

1 1

1 kx 1.5 1 kx 1.5

min: 0 max: 0.4000 min: 0 max: 0.2792

t = 2.5 t = 2.5
1.5 1.5
ky

ky

1 1

1 kx 1.5 1 kx 1.5

min: 0 max: 0.7023 min: 0 max: 0.5211

t = 4.5 t = 4.5
1.5 1.5
ky

ky

1 1

1 kx 1.5 1 kx 1.5

Fig. 6 Zoomed view of the 2D FFT of the corresponding numerical error on C-grid
shown in Fig. 5 at indicated time instants.

28
 t=0
250 min = -0.969331
max = 0.987081
200

Y
150

100
100 150 X 200 250

Non-optimal RK method with OSCS scheme [32]

(i) t = 10.8 (ii) t = 14
min = -0.360671 min = -0.324764
250 250
max = 0.362508 max = 0.322009

200 200
Y

150 150

100 100
100 150 X 200 250 100 150 X 200 250

Classical RK method with staggered compact scheme [33]

(iii) t = 10.8 (iv) t = 14
250
min = -0.469865 250
min = -0.457875
max = 0.470552 max = 0.457011

200 200
Y

150 150

100 100
100 150 X 200 250 100 150 X 200 250

Optimized RK method with OSCS scheme [32]

(v) t = 10.8 (vi) t = 14
250
min = -0.470054 250
min = -0.462006
max = 0.476650 max = 0.453330

200 200
Y

150 150

100 100
100 150 X 200 250 100 150 X 200 250

Fig. 7 Numerical solution of LRSWE on C-grid for 2D wave-packet, given in

Eq. (45) using indicated numerical methods with grid aspect ratio λ = ∆y/∆x = 1,
Nc = 0.2 and f ∆t = 1.2 × 10−2 .
29