**Department of Mathematics**

## High-order compact finite difference schemes

## for the spherical shallow water equations

### Sarmad Ghader and Jan Nordstrm

### LiTH-MAT-R–2013/09–SE

### High-order compact finite difference schemes for the spherical

### shallow water equations

### Sarmad Ghader

*a∗*

### and Jan Nordstr¨

### om

*b*

*a _{Institute of Geophysics, University of Tehran, Iran}*

*b _{Division of Computational Mathematics, Department of Mathematics, Link¨}_{oping University, Sweden}*

SUMMARY

This work is devoted to the application of the super compact finite difference (SCFDM) and the
combined compact finite difference (CCFDM) methods for spatial differencing of the spherical shallow
water equations in terms of vorticity, divergence and height. Five high-order schemes including the
fourth-order compact, the sixth-order and eighth-order SCFDM and the sixth-order and eighth-order
CCFDM schemes are used for spatial differencing of the spherical shallow water equations. To advance
the solution in time, a semi-implicit Runge-Kutta method is used. In addition, to control the
non-linear instability and avoiding the polar problem a high-order spatial filter is proposed. An unstable
barotropic mid-latitude zonal jet is employed as an initial condition. For the numerical solution of the
elliptic equations in the problem, a direct hybrid method which consists of using a high-order compact
scheme for spatial differencing in the latitude coordinate and a fast Fourier transform in longitude
coordinate is utilized. The convergence rate for all methods is studied and verified. Qualitative
and quantitative assessment of the results, such as measures of maximum vorticity gradient, power
spectrum of total energy, relative change in potential enstrophy and potential palinstrophy, reveal
that the sixth-order and eighth-order CCFDM and the sixth-order and eighth-order SCFDM methods
lead to a remarkable improvement of the solution over the fourth-order compact method. It is also
shown that the performance of the sixth-order and eighth-order CCFDM methods are superior to the
sixth-order and eighth-order SCFDM methods. Copyright c*⃝ 2013 John Wiley & Sons, Ltd.*
key words: Compact finite difference; High-order methods; Spherical shallow water equations;
Numerical accuracy

1. Introduction

The shallow water equations are widely used in various oceanic and atmospheric problems. This model is applied to a fluid layer of constant density in which the horizontal scale of the flow is much greater than the layer depth (e.g., [1, 2]). Even though, the dynamics of two dimensional shallow water model is less general than the three dimensional general circulation models, it is often preferred because of its mathematical and computational simplicity. The

*∗*_{Correspondence to: Sarmad Ghader, Institute of Geophysics, University of Tehran, North Kargar Ave., Tehran,}

spherical shallow water equations are usually used as the first step to examine new numerical algorithms to be used in more complex global climate and weather prediction models.

Compact finite-difference schemes are simple and powerful ways to reach the objectives of high accuracy and low computational cost (e.g., [3, 4, 5]). Compared with the traditional explicit finite-difference schemes of the same-order, compact schemes are significantly more accurate along with the benefit of using smaller stencil sizes. The compact finite difference methods can not be formulated on so called summation-by-parts (SBP) form [6, 7, 8, 9] which is a serious drawback for finite domain problems since the stability can not be proved. However, in this work we only consider periodic domains and this drawback is not important.

One of the main problems with the spherical harmonic method [10], which is the dominant global modelling method, is its rapid increase of computational cost with resolution. Therefore, examination of other high-order schemes, such as compact finite difference methods, as alternatives is important. Previous applications of compact schemes to idealized models of the atmosphere and oceans have shown that they are a promising alternative for the numerical simulation of geophysical fluid dynamics problems (e.g., [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]). Motivated by the remarkable performance of the sixth-order and eighth-order super compact finite difference methods compared with the fourth-order compact method in Ref. [18], this paper is devoted to the application of two families of high-order compact finite difference methods to the spatial differencing of the spherical shallow water equations. The super compact finite difference (SCFDM) and the combined compact finite difference (CCFDM) methods [26, 27, 21, 22, 23, 24, 25] are considered in this study.

The remainder of this paper is organized as follows. The spherical shallow water equations are presented in section 2. Details of the semi-implicit Runge-Kutta time advancing method and formulation of the high-order compact methods used for spatial differencing of the spherical shallow water equations are given in sections 3 and 4, respectively. The procedure used for the numerical solution of the elliptic equations is described in section 5. The spatial filter formulation is given in section 6. The order of convergence for the high-order compact methods used in this study is presented in section 7. Details of the test case used for the numerical solution, numerical results and discussions are given in section 8. Finally, section 9 presents the conclusions.

2. The shallow water equations

The momentum and mass continuity representation of the spherical shallow water equations in vector form can be written as (e.g., [28])

**Dv**

*Dt* *+ f ˆ***z****× v + c**

2_{∇˜h = 0}_{(1)}

*∂˜h*

*∂t* +* ∇ · [(1 + ˜h)v] = 0.* (2)

**In (1) and (2), v = uˆi + vˆj is the horizontal velocity vector with u and v being the velocity****components in longitudinal (λ) and latitudinal (ϕ) directions, respectively. ˆi and ˆj are the unit****vectors in longitudinal and latitudinal directions, respectively. D()/Dt = ∂()/∂t + (v**· ∇)() is*the substantial derivative and f = 2Ω sin ϕ is the Coriolis parameter in which Ω is the rotation*
rate of the earth. ˜*h = (h−H)/H is the dimensionless perturbation depth where h is the actual*

*depth where H is the domain area average depth. c =* *√gH is the gravity wave speed and g*

denotes the acceleration due to gravity. The unit vector in vertical direction is denoted by ˆ**z.**

Equations (1) and (2) can be rewritten in terms of relative vorticity, horizontal divergence
and height as
*∂ζ*
*∂t* = *Sζ* (3)
*∂δ*
*∂t* +*H˜h = Sδ* (4)
*∂˜h*
*∂t* *+ δ* = *Sh.* (5)
In equations (3)-(5),
*Sζ* =* −∇ · [(ζ + f)v], Sδ= f (ζ− f˜h) − ∇*2(

*2 2 )*

**|v|***=*

**− ∇ · (−ζvˆi + ζuˆj) − u ˆβ, S**h*and*

**−∇ · (˜hv)***H = c*2

*∇*2

*− f*2is the modified Helmholtz operator, ˆ

*β = (df /dϕ)/a = (2Ω/a) cos ϕ and a*

is the radius of the earth. The relative vorticity and divergence are defined as

*ζ =* 1
*a cos ϕ*
(
*∂v*
*∂λ* *−*
*∂u cos ϕ*
*∂ϕ*
)
*,* *δ =* 1
*a cos ϕ*
(
*∂u*
*∂λ* +
*∂v cos ϕ*
*∂ϕ*
)
*.* (6)

In addition, to find the velocity components and to close the system of equations the Helmholtz
**decomposition v = ˆz***×∇ψ+∇χ is employed where ψ and χ denote streamfunction and velocity*

potential, respectively. To find the streamfunction and velocity potential, the following Poisson equations must be solved:

*∇*2_{ψ = ζ,}* _{∇}*2

_{χ = δ.}_{(7)}

3. The semi-implicit Runge-Kutta method

In the present work the semi-implicit Runge-Kutta (hereafter, SIRK) time advancing method developed in [29] for the flux convergence form of the shallow water equations in Cartesian coordinates, is extended to the vorticity, divergence and height representation of the spherical shallow water equations. For the vorticity-divergence form of the shallow water equations, the SIRK formulation is applied to equations (4) and (5) [30]. The Runge-Kutta formulation used in the present work is a third-order explicit three step method [31].

Equations (4) and (5) can be rewritten in the following vector form:

**∂q**

*∂t* * = L(q) + N (q)* (8)

**where q = (δ, ˜**h)T_{, L(q) = (}_{−H˜h, −δ)}T_{and N (q) = (S}

*δ, Sh*)*T. The superscript T denotes*

transpose. The semi-implicit approach which employs the third-order Runge-Kutta method for time advancement of equation (8) leads to

**q***n+αk − qn*

*αk∆t*=

**L(q***n+αk*

**n**_{) + L(q}_{)}2

**+ N (q***n+αk−1*

_{)}

_{(9)}

*where the superscript n denotes the time step. The three steps of the Runge-Kutta method*
*are written in a compact form where subscript k (k = 1, 2, 3) denotes the three steps of the*
*Runge-Kutta method and the values of αk* *are α*0 *= 0, α*1*= 1/3, α*2*= 1/2 and α*3= 1. The
following auxiliary variables are defined

¯
*δn+αk*_{=} *δ*
*n+αk _{+ δ}n*
2

*,*¯

*hn+αk*

_{=}˜

*hn+αk*

_{+ ˜}

*2*

_{h}n*.*(10)

By substituting the auxiliary variables (10) into equation (9), the following equations for the auxiliary variables can be found

(*H −* 4
*α*2
*k∆t*2
)¯*hn+αk* _{=} 2
*αk∆t*
*(δn− S _{h}n+αk−1) + S_{δ}n+αk−1−* 4

*α*2

*k∆t*2 ˜

*hn*(11) ¯

*δn+αk*

_{=}2

*αk∆t*(˜

*hn− ¯hn+αk*

_{) + S}n+αk−1*h*

*.*(12)

*Equation (11) is a modified Helmholtz equation. The unknown variables δn+αk* _{and ˜}_{h}n+αk_{are}

found by substitution of the auxiliary variables into equation (10). The time discretization of the vorticity equation (3) by the third-order Runge-Kutta method leads to:

*ζn+αk _{= ζ}n_{+ α}*

*k∆tS*
*n+α _{k−1}*

*ζ* (13)

*where the definitions of n and αk* are given above.

4. Spatial differencing

For spatial differencing of the first and second derivatives in equations (11), (12), (13) and the diagnostic equations (7), two families of compact finite difference schemes, i.e., the super compact finite difference method (SCFDM) and the combined compact finite difference method (CCFDM) are used.

In the SCFDM, approximations of the first (second) derivatives depends on all the nodal values and their odd (even) derivatives, but the CCFDM approximates the first and second derivatives simultaneously. See references [26, 27, 16] and [21, 22, 23, 24, 25] for more information on details of the derivation and alternative forms of the SCFDM and CCFDM, respectively.

*4.1. The super compact method (SCFDM)*

The SCFDM relations for approximation of the first derivative of an arbitrary function Φ in a uniform grid are

where
**F***j* =
Φ*<1> _{j}*
Φ

*<3>*.. . Φ

_{j}*<2M*

_{j}*−1>*

*F*

**, E***j*=

*D*0Φ

*j*0 .. . 0

*, Q*f = 1

*d* 1 1! 1 3! 1 5!

*. . .*1

*(2M− 1)!*

*−d*2 2

*D*2 1 2! 1 4!

*. . .*1

*[2(M− 1)]!*0

*−d*

_{2}2

*D*2 1 2!

*. . .*1

*[2(M− 2)]!*

*. . .*

*. . .*

*. . .*

*. . .*

*. . .*0 0

*. . .*

*−d*

_{2}2

*D*2 1 2! In (14), the difference operators are defined as

*D*0Φ*j*= (Φ*j+1− Φj _{−1})/2d,*

*D*2Φ

*j*= (Φ

*j+1− 2Φj*+ Φ

*j*2

_{−1})/d*in which d denotes the grid spacing and Φ<l> _{j}*

*= dl*(

*∂l*

_{Φ}

*∂xl*

)

*j*

. Qf *is a M × M matrix, F and*

**E**F

*are M dimensional vectors. The expression Φ<2l*

_{j}*−1>/d2l−1*

*(l = 1, 2, . . . , M ) approximates*

*∂2l−1*_{Φ}

*∂x2l _{−1}*

*with an accuracy of order 2(M− l + 1) [26, 27].*

In a similar way, the SCFDM relations for approximation of the second derivative on a uniform grid are

Qs**S***j***= E**S*j* (15)
in which
**S***j*=
Φ*<2> _{j}*
Φ

*<4>*.. . Φ

_{j}*<2M >*

_{j}*S*

**, E***j*= 1 2

*D*2

_{Φ}

*j*0 .. . 0

*, Q*s= 1

*d*2 1 2! 1 4! 1 6!

*. . .*1

*(2M )!*

*−d*2 2

*D*2 1 2! 1 4!

*. . .*1

*[2(M− 1)]!*0

*−d*2 2

*D*2 1 2!

*. . .*1

*[2(M− 2)]!*

*. . .*

*. . .*

*. . .*

*. . .*

*. . .*0 0

*. . .*

*−d*

_{2}2

*D*2 1 2! Here, Qs

*is a M*S

**× M matrix, S and E***are M dimensional vector. The expression Φ<2l>j*

*/d2l*

approximates *∂*
*2l _{ϕ}*

*∂x2l* *with an accuracy of order 2(M− l + 1) [26, 27].*

*Using different values of M (M = 1, 2, 3,· · · ) in equations (14) and (15) leads to different*
orders of the SCFDM relations for approximation of the first and second derivatives. For
*example, M = 1 leads to the conventional second-order central finite difference formulation.*
*The fourth-order Pade type compact finite difference are found by using M = 2. Higher values*
*of M leads to the higher order SCFDM relations.*

*4.1.1. The SCFDM relations* *Using M = 2, 3, 4 in equation (14) leads to the fourth-order,*
sixth-order and eighth-order SCFDM relations for the approximation of the first derivative.
The fourth-order relation reads

1
6Φ
(1)
*j _{−1}*+
2
3Φ
(1)

*j*+ 1 6Φ (1)

*j+1−*1

*2d*(Φ

*j+1− Φj−1) = 0,*(16)

the sixth-order relations are
{
Φ*j+1− Φj−1− 2dΦ*
(1)
*j* *−*
*d*3
60Φ
(3)
*j−1−*
*3d*3
10Φ
(3)
*j* *−*
*d*3
60Φ
(3)
*j+1*= 0
Φ(1)* _{j+1}− 2Φ*(1)

*+ Φ*

_{j}*(1)*

_{j}

_{−1}−d_{12}2Φ(3)

_{j}_{−1}−5d_{6}2Φ(3)

_{j}*−d*

_{12}2Φ(3)

*(17) and the eighth-order relations are found as*

_{j+1}= 0,
Φ*j+1− Φj _{−1}− 2dΦ*
(1)

*j*

*−*

*d*3 3Φ (3)

*j*

*−*

*2d*5 5040Φ (5)

*j−1−*

*2d*5 126Φ (5)

*j*

*−*

*2d*5 5040Φ (5)

*j+1*= 0 Φ(1)

*(1)*

_{j+1}− 2Φ*+ Φ(1)*

_{j}*2Φ(3)*

_{j}_{−1}− d

_{j}*−*

_{360}

*d*4 Φ(5)

_{j}_{−1}−7d_{90}4Φ

*(5)*

_{j}*−*

_{360}

*d*4 Φ(5)

*= 0 Φ(3)*

_{j+1}*(3)*

_{j+1}− 2Φ*+ Φ(3)*

_{j}*2 12Φ (5)*

_{j}_{−1}−d*j−1−*

*5d*2 6 Φ (5)

*j*

*−*

*d*2 12Φ (5)

*j+1= 0.*(18)

*Similarly, using M = 2, 3, 4 in equation (15) leads to the fourth-order, sixth-order and*
eighth-order SCFDM relations for the approximation of the second derivative. The
fourth-order relation is
1
12Φ
(2)
*j−1*+
5
6Φ
(2)
*j* +
1
12Φ
(2)
*j+1−*
1
*d*2(Φ*j+1− 2Φj*+ Φ*j−1) = 0,* (19)
the sixth-order relations are

{
Φ*j+1− 2Φj*+ Φ*j−1− d*2Φ
(2)
*j* *−*
*d*4
360Φ
(4)
*j−1−*
*7d*4
90Φ
(4)
*j* *−*
*d*4
360Φ
(4)
*j+1*= 0
Φ(2)* _{j+1}− 2Φ*(2)

*+ Φ(2)*

_{j}*2 12Φ (4)*

_{j}_{−1}−d*j−1−*

*5d*2 6 Φ (4)

*j*

*−*

*d*2 12Φ (4)

*j+1= 0,*(20)

and the eighth-order relations are found as
Φ*j+1− 2Φj+1*+ Φ*j _{−1}− d*2Φ
(2)

*j*

*−*

*d*4 12Φ (4)

*j*

*−*

*d*6 20160Φ (6)

*j*

_{−1}−*3d*6 1120Φ (6)

*j*

*−*

*d*6 20160Φ (6)

*j+1*= 0 Φ(2)

*(2)*

_{j+1}− 2Φ*+ Φ(2)*

_{j}*2*

_{j}_{−1}− d_{Φ}(4)

*j*

*−*

*d*4 360Φ (6)

*j−1−*

*7d*4 90Φ (6)

*j*

*−*

*d*4 360Φ (6)

*j+1*= 0 Φ(4)

*(4)*

_{j+1}− 2Φ*+ Φ(4)*

_{j}

_{j}_{−1}−d_{12}2Φ

*(6)*

_{j}

_{−1}−5d_{6}2Φ

*(6)*

_{j}*−d*

_{12}2Φ(6)

*(21)*

_{j+1}= 0.*In these equations, the superscript (k) denotes kth spatial derivative.*

*4.2. The combined compact method (CCFDM)*

The CCFDM relations for the first and second derivatives on a uniform grid can be written as

Qc**T***j***= G***j* (22)
where
**T***j* =
Φ*<1> _{j}*
Φ

*<2>*.. . Φ

_{j}*<2N*

_{j}*−1>*Φ

*<2N >*

_{j}*= *

**, G**j*D*0Φ

*j*

*D*2Φ

*j*0 .. . 0

Qc=
*a*
1! 0
*a*
3! 0 *. . .*
*a*
*(2N− 1)!* 0
0 *2b*
2! 0
*2b*
4! *. . .* 0
*2b*
*(2N )!*
*−D*0 *a*
1! 0
*a*
3! *. . .* 0
*a*
*(2N* *− 1)!*
*−D*2 _{0} *2b*
2! 0 *. . .*
*2b*
*(2N )!* 0
0 *−D*0 *a*
1! 0 *. . .* 0
*a*
*(2N* *− 3)!*
0 *−D*2 _{0} *2b*
2! *. . .*
*2b*
*(2N− 2)!* 0
0 0 *−D*0 *a*
1! *. . .* 0
*a*
*(2N* *− 5)!*
0 0 *−D*2 _{0} _{. . .}*2b*
*(2N− 4)!* 0
..
. ... ... ... ...
0 0 0 0 *. . .* *. . .* *. . .*
*and a = 1/d, b = 1/d*2_{. For the CCFDM method, Q}

c *is a 2N × 2N matrix and T and G are*

*2N dimensional vectors. Equation (22) approximates both the first and second derivatives of*

*Φ with an accuracy of order 2N .*

*Substitution of different values of N (N = 1, 2, 3,· · · ) in equation (22) leads to various orders*
of the CCFDM relations for approximation of the first and second derivatives simultaneously.
*Similar to the SCFDM, N = 1 and N = 2 lead to the second-order central finite difference*
*and fourth-order Pade type compact finite difference relations. Higher values of N leads to the*
higher order CCFDM relations.

*4.2.1. The CCFDM relations* *For the CCFDM method by using N = 3, 4 in equation (22),*
the sixth-order and eighth-order relations for the approximation of the first and second
derivatives are obtained. The sixth-order relations are

{
7
16(Φ
(1)
*j+1*+ Φ
(1)
*j−1*) + Φ
(1)
*j* *−*
*d*
16(Φ
(2)
*j+1− Φ*
(2)
*j−1*)*−*
15
*16d*(Φ*j+1− Φj−1*) = 0
9
*8d*(Φ
(1)
*j+1− Φ*
(1)
*j−1*)*−*
1
8(Φ
(2)
*j+1*+ Φ
(2)
*j−1*) + Φ
(2)
*j* *−*
3
*d*2(Φ*j+1− 2Φj*+ Φ*j−1) = 0,*
(23)

and the eighth-order relations are found as
19
32(Φ
(1)
*j+1*+ Φ
(1)
*j−1*) + Φ
(1)
*j* *−*
*d*
8(Φ
(2)
*j+1− Φ*
(2)
*j−1*) +
*d*2
96(Φ
(3)
*j+1*+ Φ
(3)
*j−1*)*−*
35
*32d*(Φ*j+1− Φj−1*) = 0
29
*16d*(Φ
(1)
*j+1− Φ*
(1)
*j−1*)*−*
5
16(Φ
(2)
*j+1*+ Φ
(2)
*j−1*) + Φ
(2)
*j* +
*d*
48(Φ
(3)
*j+1− Φ*
(3)
*j−1*)*−*
4
*d*2(Φ*j+1− 2Φj*+ Φ*j−1*) = 0
*−105*
*16d*2(Φ
(1)
*j+1*+ Φ
(1)
*j _{−1}*) +
15

*8d*(Φ (2)

*j+1− Φ*(2)

*j*)

_{−1}*−*3 16(Φ (3)

*j+1*+ Φ (3)

*j*) + Φ (3)

_{−1}*j*+ 105

*16d*3(Φ

*j+1− Φj*(24)

_{−1}) = 0.*4.3. The numerical grid*

The shallow water equations are solved on the doubly periodic longitude-latitude grid system
(0*≤ λ ≤ 2π, −π/2 ≤ ϕ ≤ π/2) of reference [32]. The grid points in this system are arranged*
such that no grid points are located at the poles. This grid system has the advantage of
avoiding singularity at the poles and enabling the periodic boundary conditions to be imposed

in both longitude and latitude directions. It should be noted that the sign of a quantity must
be adjusted properly when using the periodic boundary condition in the latitude direction,
depending on whether the quantity is a scalar or a vector. The grid points are defined as
*(λi, ϕj) = (i∆λ,−π/2 + (j − 1/2)∆ϕ) where ∆λ = 2π/Nλ* *and ∆ϕ = π/Nϕ. Nλ* *and Nϕ* are

the number of grid points in latitude and longitude directions, respectively.

5. The elliptic equations

The SIRK formulation used to advance the shallow water equations in time leads to a set of modified Helmholtz equations (11) at each successive step of the Runge-Kutta method. In addition, to find the velocity components, the two Poisson equations given in (7) must be solved at each step.

The elliptic equations (11) and (7) are solved using a direct method. The procedure is similar to one developed in [33, 34]. However, the following changes are applied to make the procedure easier. The spatial discretization of the derivatives in longitude are performed by discrete Fourier series and the FFT (Fast Fourier Transform) technique. In the latitudinal direction we discretize using different orders of the SCFDM and CCFDM schemes.

*5.1. The Helmholtz equation*

The modified Helmholtz equation (11) can be written as

*∇*2_{¯}_{h}_{− β¯h = R}_{(25)}

where the Laplacian operator is

*∇*2_{=} 1
*a*2_{cos}2_{ϕ}*∂*2
*∂λ*2+
1
*a*2* _{cos ϕ}*
(

*∂*

*∂ϕcos ϕ*

*∂*

*∂ϕ*)

*.*

*In (25), β is a function of ϕ and R denotes the known right hand side of the elliptic equation.*
The longitudinal parts of ¯*h and R in (25) are approximated by truncated Fourier series. By*

substitution of the Fourier components ¯*h = ˆhmeimλ* *and R = ˆRmeimλ* *with wave number m*

into equation (25) we have
(
*−m*2
*a*2_{cos}2* _{ϕ}− β*
)
ˆ

*h +*1

*a*2

*d*2ˆ

_{h}*dϕ*2

*−*

*tan ϕ*

*a*2

*dˆh*

*dϕ*= ˆ

*R.*(26)

Now, the fourth-order compact, the sixth- and eighth-order SCFDM and CCFDM methods are used to discretize the first and second latitudinal derivatives in equation (26). We arrive at the following block tridiagonal system of equations

AS_{j}**U**S* _{j}_{−1}*+ BS

_{j}**U**S

*+ CS*

_{j}

_{j}**U**S

_{j+1}**= D**S

*(27)*

_{j}for the different methods are
**U**C4* _{j}* =
ˆ

*hj*ˆ

*h*(1)

*ˆ*

_{j}*h*(2)

* *

_{j}*SC6*

** , U***j*= ˆ

*hj*ˆ

*h*(1)

*ˆ*

_{j}*h*(2)

*ˆ*

_{j}*h*(3)

*ˆ*

_{j}*h*(4)

* *

_{j}*SC8*

**, U***= ˆ*

_{j}*hj*ˆ

*h*(1)

*ˆ*

_{j}*h*(2)

*ˆ*

_{j}*h*(3)

*ˆ*

_{j}*h*(4)

*ˆ*

_{j}*h*(5)

*ˆ*

_{j}*h*(6)

* *

_{j}*CC6*

**, U***= ˆ*

_{j}*hj*ˆ

*h*(1)

*ˆ*

_{j}*h*(2)

* *

_{j}*CC8*

** , U***j*= ˆ

*hj*ˆ

*h*(1)

*ˆ*

_{j}*h*(2)

*ˆ*

_{j}*h*(3)

* *

_{j}*.*(28) The superscripts C4, SC6, SC8, CC6, CC8 denote the fourth-order compact, the sixth-order SCFDM, the eighth-order SCFDM, the sixth-order CCFDM and the eighth-order CCFDM methods, respectively. Equation (27) is a three-point implicit formulation which solution leads to a linear block tridiagonal system of equations. The block size for the fourth-order compact is 3

*× 3, for the sixth-order SCFDM is 5 × 5, for the eighth-order SCFDM is 7 × 7, for the*sixth-order CCFDM is 3

*× 3 and for the eighth-order CCFDM is 4 × 4.*

*To close the system at south pole (j = 1) and north pole (j = Nϕ*) boundary conditions are

needed. The boundary values are obtained by using the symmetry constraint of the Fourier
coefficients (e.g., [33]). For the grid system used in the present work, the symmetry constraints
*imply that at j = 1 we have ˆh(k)*_{0} = (*−1)k*(*−1)m*ˆ*h(k)*_{1} *and for j = Nϕ*ˆ*h*

*(k)*

*Nϕ*+1= (*−1)*

*k*_{(}_{−1)}m_{ˆ}_{h}(k)

*Nϕ*.

*The superscript (k) denotes as before the kth derivative (k = 0, 1, 3,· · · ) and m denotes the*
wave number. There is one sign change due to odd and even wave numbers and another sign
change related to the spatial derivative.

By inserting these values into equation (27) at boundaries and solving the resulting block
tridiagonal system of equations, the values of ˆ*h are obtained and then by using FFT, the values*

of ¯*hi,j* are found. The details of matrix coefficients AS*j*, BS*j*, CS*j* **and right hand side vector D**S*j*

are given in appendix I.

*5.2. The Poisson equation*

The Poisson equations in (7) have the following form

*∇*2_{ψ = ζ.}_{(29)}

The procedure for the numerical solution of the Poisson equations is similar to the procedure for
*the modified Helmholtz equation above and can be reached by letting β = 0 in equation (26).*
*For nonzero wave numbers (m* *̸= 0) equation (29) can be solved in a similar manner as the*
modified Helmholtz equation.

*The procedure for the zeroth Fourier mode (m = 0) is different. It can be seen that by using*

*β = 0 and m = 0 in (26), ˆh vanishes and the Poisson equation for wave number m = 0 is*

singular, since the solution is unique within a constant. In fact, a Poisson equation on a sphere must satisfy the following compatibility condition (e.g.,[35, 33])

∫ *2π*
0

∫ *π/2*
*−π/2*

*ζ cos ϕdϕdλ = 0* (30)

*By substitution of the Fourier components in equation (29), and considering m = 0 we have*
1
*a*2
*d*2*ψ*¯
*dϕ*2 *−*
*tan ϕ*
*a*2
*d ¯ψ*
*dϕ* = ¯*ζ* (31)

in which ¯*ψ = ˆψ(m = 0) and ¯ζ = ˆζ(m = 0). After application of the Fourier method in longitude*

direction the equivalent form of the compatibility condition (30) for the latitudinal direction

is _{∫}

*π/2*
*−π/2*

¯

*ζ cos ϕdϕ = 0.* (32)

To impose the compatibility condition, first the fourth-order compact, the sixth- and
eighth-order SCFDM and the sixth- and eighth-eighth-order CCFDM methods are used to discretize the
first and second latitudinal derivatives in equation (31) in the same way as for (26), and then
the right hand side of equation (31) in each discrete version is replaced by ¯*ζj− P where P is a*

*constant. P can be found by using the discrete version of the compatibility condition (32) as*

*P =*

*j=N*∑*ϕ*

*j=1*

¯

*ζjcos ϕj.* (33)

Now, by finding the values of ¯*ψj* for the zeroth Fourier mode and having the values of ˆ*ψj*

for other wave numbers we employ the FFT to find the solution of the Poisson equation. The matrix coefficients and right hand side vector for different compact methods used to discretize equation (31) are also described in appendix I.

6. The spatial filter

The generation of fine-scale vortical structures by the nonlinear advection terms is a dominant feature of the shallow water equations. A problem facing any grid-based numerical method developed for the shallow water equations is misinterpretation of the fine-scale structures due to insufficient resolution and the potential for non-linear numerical instability. In this work, a simple spatial filter is developed to control the non-linear instability. In addition, on a regular latitude-longitude system the lines of constant longitude converge at the poles which leads to a high concentration of grid points near poles. The cluster of grid points forces a short time step to satisfy the CFL stability condition and this is often referred to as the polar problem [36]. The spatial filter can act as a smoother to avoid the polar problem.

*The general form of a one-dimensional spatial filter (for example in latitude, ϕ) to filter an*
arbitrary field, Υ, can be written as

˜
Υ*j* = Υ*j+ σp(∆ϕ*2)*p*
(
*∂*2_{Υ}
*∂ϕ*2
)*p*
*j*
*,* *p = 1, 2, 3, . . .* (34)

where ˜Υ*j* denotes the filtered values. This formulation is used to construct a low-pass filter

*which removes wave-lengths shorter than 2∆ϕ. Using higher values of p moves the filtering*
effect towards the shorter waves. The order of accuracy of the filter is determined by the finite
difference method used to discretize the even derivative in equation (34).

*The order of even derivatives (p) and the discretization method in equation (34) can be*
selected by the requirements that (i) the dissipation is kept minimal for stable solutions and
(ii) the desired order of accuracy is not degraded by the action of the filter.

*The unknown coefficient σp* is determined by setting the location of cutoff in the transfer

*function of the filter (e.g., [37]). The transfer function, Tf*, of the filter is found using Fourier

transformation of equation (34) and is defined as ˜*Υ(eiκϕ) = Tf(κ)eiκϕ*, where i =

*√*
*−1 and*
*κ denotes wave number. For a low-pass filter which removes the wave-lengths shorter than*

*2∆ϕ it is necessary to have Tf(κ*max*) = 0 where κ*max *= π/∆ϕ. Therefore, the spatial filter*
coefficient can be found as

*σp*=
(
*−1*
*T*2*(κ*max)
)*p*
(35)
*where T*2*is the transfer function of the second derivative. Table I presents values of σp* for the

various schemes used in this study.

*Table I. values of the spatial filter coefficient, σp*, for different methods.

Method *σp*
Fourth-order compact (*−1*_{6} )*p*
Sixth-order SCFDM (* _{7.05882353}−1* )

*p*Sixth-order CCFDM (

*−1*)

_{9.6}*p*Eighth-order SCFDM (

*)*

_{7.67741935}−1*p*Eighth-order CCFDM (

*)*

_{9.84615385}−1*p*

The same spatial differencing method that is applied to discretize the governing equations is
used to estimate the even derivative in the spatial filter (34). In addition, the spatial filter
is sequentially applied in latitude and longitude directions at each time step of the time
*integration. We have used p = 3 for the filter power in both latitude and longitude directions.*
*However, for the longitude direction we have used a slightly more diffusive filter with p = 2*
near the poles where *0.7π*_{2} *≤ ϕ ≤* *π*_{2} and*−π*_{2} *≤ ϕ ≤ −0.7π*_{2} , to better remove the restriction of
the polar problem on the time step.

7. The order of convergence

The spatial order of convergence of high-order compact schemes used in this work is assessed by considering a Helmholtz equation with a known analytical solution. In fact, at a fixed time the main part of the total algorithm consists of the solution of the elliptic equations and it is therefore a meaningful test.

The following elliptic equation is solved

*∇*2_{Ψ}_{− aΨ = R}

Ψ (36)

*where R*Ψ*= 2ν sin ϕ− 30K sin ϕ cos*4*ϕ cos 4λ + a*3*ν sin ϕ− a*3*K cos*4*ϕ sin ϕ cos 4λ and a is the*
*radius of the Earth. The constants ν and K are defined as ν = K = 7.848× 10−6*s*−1*.

The Rossby-Haurwitz wave [28], is the analytical solution of equation (36) given by
Ψ =*−a*2*ν sin ϕ + a*2*K cos*4*ϕ sin ϕ cos 4λ.* (37)
To investigate the order of convergence, the numerical solution of equation (36) generated
by the 4th-order compact, 6th-order SCFDM, 6th-order CCFDM, 8th-order SCFDM and
8th-order CCFDM methods is compared with the analytical solution (37). The following
normalized global error is used for the error measurement [28]:

*l*2(Ψ) =
{
*I*(*|Ψ*n*(λ, ϕ)− Ψ*E*(λ, ϕ)|*
2)}*1/2*
*{I (|Ψ*E*(λ, ϕ)|*
2_{)}*}1/2* (38)

where Ψnand Ψ*E* *are the numerical and exact solutions, respectively. The function I is defined*

as
*I(Ψ) =* 1
*4π*
*Nλ*
∑
*k=1*
*Nϕ*
∑
*l=1*
Ψ* _{k,l}cos ϕ∆ϕ∆λ* (39)

*in which Nλ* *and Nϕ* are the number of grid points in longitude and latitude directions and

Ψ*k,l* *is the argument at grid point (λk, ϕl*).

For the numerical solution of equation (36), the compact finite difference methods are
only applied for discretization in the latitudinal direction (see section 5). Therefore, to find
*the order of convergence we use ∂Ψ/∂ϕ instead of Ψ. Figure 1 presents the normalized*
*global error l*2 *for ∂Ψ/∂ϕ at different uniform horizontal resolutions in longitude and latitude*
directions calculated by the 4th-order compact, 6th-order SCFDM, 6th-order CCFDM,
8th-order SCFDM and 8th-8th-order CCFDM schemes. Figure 1 clearly shows that the convergence
rates of the different solutions are in agreement with theoretical order of convergence. In
addition, it is seen that the 6th-order CCFDM method produces less error than the 6th-order
SCFDM method and the 8th-order CCFDM generates less error than the 8th-order SCFDM
method, respectively.

8. Numerical results

*8.1. Initial condition*

The test case proposed in reference [38] which is an unstable barotropic mid-latitude zonal jet is used as the starting point for the numerical solution and evaluation of the spherical shallow water equations. The initial vorticity filed of the unstable barotropic zonal jet is

*ζ =* 1

*a cos ϕ*

(

*u sin ϕ− cos ϕ∂u*
*∂ϕ*

)

*.* (40)

*For this test case v = 0 and δ = 0 and*

*u =*
0 for *ϕ≤ ϕ*0
*u*max
*en*
( 1
*(ϕ− ϕ*0*)(ϕ− ϕ*1)
) for *ϕ*0*< ϕ < ϕ*1
0 for *ϕ≥ ϕ*1
(41)

**N**_{λ}** ( = 2 N**_{φ}** )**
**L2**
** (**
**E**
**rr**
**o**
**r)**
100 200 300 400 500
10-16
10-14
10-12
10-10
10-8
10-6
10-4
10-2
100 **4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**P=4**
**P=6**
**P=8**

*Figure 1. Normalized l*2 *error of ∂Ψ/∂ϕ field versus number of grid points obtained using 4th-order*
compact method, 6th-order SCFDM, 6th-order CCFDM, 8th-order SCFDM and 8th-order CCFDM.
*Lines with slops 4, 6 and 8 are also shown. A uniform grid resolution is used in both λ and ϕ directions.*

is used in which the constants are

*u*max= 80 ms*−1,* *ϕ*0=
*π*
7*,* *ϕ*1=
*π*
2 *− ϕ*0*,* *en*= exp[
*−4*
*(ϕ*1*− ϕ*0)2
*].*
The following relation is used to find the initial height field

*h =* 1
*g*
{
*gh*0*−*
∫ *ϕ*
*au(ϕ′*)
[
*f +tan(ϕ*
*′*_{)}
*a* *u(ϕ*
*′*_{)}]* _{dϕ}′*
}
(42)

*where the value of constant h*0needs to be found such that the domain area average of depth is

*H = 10000 m. Here, we have used the trapezoid rule to evaluate the integral in equation (42)*

*numerically and finding h to machine precision.*

Next, the following perturbation is added to the basic height in equation (42) to trigger the barotropic instability

*h′(λ, ϕ) = hpcos(ϕ)e−(λ/γ)*

2

*e−[(ϕ*2*−ϕ)/ ¯β]*2_{,}_{for} *− π < λ < π* _{(43)}

*where λ is longitude and the constants are hp= 120 m, γ = 1/3, ϕ*2*= π/4 and ¯β = 1/15.*

*8.2. Grid resolution and time step*

*The grid resolutions used in the numerical experiments reported here are Nλ* *× Nϕ* =

128*×64, 256×128, 512×256. To find the time steps a fixed Courant number of less than unity*
based on gravity-wave speed is used. Using small time steps reduces the frequency distortion

*of the semi-implicit time-stepping scheme on gravity waves. The time steps ∆t = 270, 60, 15*
seconds have been used for the successive grids.

*8.3. Results and discussions*

We begin by providing some qualitative results for the vorticity field. Figure 2 presents the time
*evolution of the vorticity field for Nλ× Nϕ*= 512*× 256 resolution at t = 0, 2, 4, 6 days (1 day*

= 86400 s) calculated by the eighth-order CCFDM. A qualitative comparison of these results
with those presented in [38, 17] indicates the validity of computations. In addition, figures 3
and 4 provide a qualitative comparison of the vorticity field calculated by the compact schemes
*used in this study at time t = 6 for Nλ×Nϕ*= 256*×128 and Nλ×Nϕ*= 512*×256 resolutions.*

It can also be seen that the sixth-order and eighth-order CCFDM methods produce better
results than the sixth-order and eighth-order SCFDM methods. Furthermore, it is evident
that the fourth-order compact method can not reproduce the structure of flow in particular
*the high gradients of vorticity field at Nλ× Nϕ* = 256*× 128 resolution. It can be seen that*

the solutions generated by the sixth-order and eighth-order SCFDM and CCFDM methods are very similar and a quantitative assessment of the solutions is needed to better inspect and understand the properties of different methods.

**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(t = 2)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(t = 0)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(t = 6)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(t = 4)**

Figure 2. The time evolution of the vorticity field at 2-day time intervals for unstable barotropic jet
*generated by the eighth-order CCFDM. The grid resolution is Nλ× Nϕ* = 512*× 256 and contour*

interval is 1*× 10−5* s*−1*. Solid and dashed lines are for positive and negative contours, respectively.
The zero contour is not shown and only the Northern Hemisphere is shown in figure.
Figure 5 presents time evolution of magnitude of the maximum vorticity gradient for different

**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(a)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(c)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(b)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(d)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(e)**

*Figure 3. The vorticity field at time t = 6 days, for (a) 4th-order compact method, (b) 6th-order*
SCFDM, (c) 6th-order CCFDM, (d) 8th-order SCFDM and (e) 8th-order CCFDM. The grid resolution
*is Nλ×Nϕ*= 256*×128 and the contour interval is 1×10−5*s*−1*. Solid and dashed lines are for positive

and negative contours, respectively.

methods. Clearly, the sixth-order and eighth-order CCFDM exhibit more levels of the vorticity gradients than the other methods. As sharp gradients of vorticity field are a feature of the unstable barotropic mid-latitude zonal jet test case used in this study [38, 17], it can be seen that sixth-order and eighth-order CCFDM represent the vorticity field in a better way than the other methods.

Figure 6 shows the power spectrum of total energy, which is the spectrum of 1_{2}*h(u*2* _{+ v}*2

_{) +}1

2*gh*

2* _{, at time t = 6 days for different methods and different resolutions. It is seen that the}*
spectrum of total energy for different methods can be differentiated only at the lower end of the
scales. The theory in [1] suggests that at inertial range, the spectrum of the energy decays as

**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(a)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(c)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(b)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(d)**
**Longitude**
**L**
**a**
**ti**
**tu**
**d**
**e**
-150 -100 -50 0 50 100 150
0
20
40
60
80
**(e)**

*Figure 4. The vorticity field at time t = 6 days, for (a) 4th-order compact method, (b) 6th-order*
SCFDM, (c) 6th-order CCFDM, (d) 8th-order SCFDM and (e) 8th-order CCFDM. The grid resolution
*is Nλ×Nϕ*= 512*×256 and the contour interval is 1×10−5*s*−1*. Solid and dashed lines are for positive

and negative contours, respectively.

*K−3, K being the total wave number. Figure 6 indicates that at the inertial range all methods*
used in this study exhibit this feature of the solution.

*By inspecting the time evolution of potential enstrophy C*2 *and potential palinstrophy P*
defined by
*C*2=
1
2
⟨
(1 + ˜*h)Q*2
⟩
*,* *P =* 1
2
⟨
(1 + ˜*h)|∇Q|*
⟩
(44)
*we can compare the global spatial distribution of potential vorticity (Q = (ζ +f )/(1+˜h)) in the*

solutions. In equation (44) the domain area average is denoted by*⟨ ⟩. The time variations of the*
*percentage changes in C*_{2}*′* *and P , i.e. (C*_{2}*′(t)− C*_{2}*′(0)) /C*_{2}*′(0) and (P (t)− P (0)) /P (0) multiplied*

**Time (day)**
**M**
**a**
**x**
**im**
**u**
**m**
** v**
**o**
**rt**
**ic**
**it**
**y**
** g**
**ra**
**d**
**ie**
**n**
**t**
0 1 2 3 4 5 6
0.001
0.002
0.003
0.004
0.005
0.006
0.007
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(a)**
**Time (day)**
**M**
**a**
**x**
**im**
**u**
**m**
** v**
**o**
**rt**
**ic**
**it**
**y**
** g**
**ra**
**d**
**ie**
**n**
**t**
0 1 2 3 4 5 6
0.001
0.002
0.003
0.004
0.005
0.006
0.007
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(b)**
**Time (day)**
**M**
**a**
**x**
**im**
**u**
**m**
** v**
**o**
**rt**
**ic**
**it**
**y**
** g**
**ra**
**d**
**ie**
**n**
**t**
0 1 2 3 4 5 6
0.002
0.003
0.004
0.005
0.006
0.007
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(c)**

Figure 5. Time evolution of maximum vorticity gradient for different methods at grid resolutions (a)
*Nλ× Nϕ*= 128*× 64, (b) Nλ× Nϕ*= 256*× 128 and (c) Nλ× Nϕ*= 512*× 256*

*by 100, where C*_{2}*′* *= C*2*−*
⟨

*f*2⟩* _{/2 are shown in figures 7 and 8 for different methods and different}*
grid resolutions. We have subtracted the time-independent contribution of⟨

*f*2⟩

_{/2 from C}2 to have a better quantification of the effects of damping/discretization error in destroying the global conservation of potential enstrophy.

Figures 7 and 8 show that the sixth-order and eighth-order CCFDM methods have better
global conservation properties of the potential enstrophy than the other methods. For the
*sixth-order and eighth-order CCFDM methods at Nλ× Nϕ*= 512*× 256 resolution, a less than*

*1% reduction in C*_{2}*′* during the first 3 days is followed by a much stronger loss of potential
*enstrophy to near 5% at time t = 6. For the sixth-order and eighth-order SCFDM methods at*

*Nλ× Nϕ*= 512*× 256 resolution, it is observed that a less than 2% reduction in C*2*′* during the
*first 3 days is followed by a stronger loss of the potential enstrophy to near 7% at t = 6. For*

**Log (wavenumber)**
**E**
**n**
**e**
**rg**
**y**
** p**
**o**
**w**
**e**
**r **
**s**
**p**
**e**
**c**
**tr**
**u**
**m**
0 0.5 1 1.5 2 2.5
0
5
10
15
20
25
30
35
40
45
50
**Theory**
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(a)**
**Log (wavenumber)**
**E**
**n**
**e**
**rg**
**y**
** p**
**o**
**w**
**e**
**r **
**s**
**p**
**e**
**c**
**tr**
**u**
**m**
0 0.5 1 1.5 2 2.5
0
5
10
15
20
25
30
35
40
45
50
**Theory**
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(b)**
**Log (wavenumber)**
**E**
**n**
**e**
**rg**
**y**
** p**
**o**
**w**
**e**
**r **
**s**
**p**
**e**
**c**
**tr**
**u**
**m**
0 0.5 1 1.5 2 2.5
0
5
10
15
20
25
30
35
40
45
50
**Theory**
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(c)**

*Figure 6. Power spectrum of total energy at time t = 6 day for (a) Nλ× Nϕ* = 128*× 64, (b)*

*Nλ× Nϕ*= 256*× 128 and (c) Nλ× Nϕ*= 512*× 256 resolutions. In each panel the straight line plots*

*a K−3*spectrum, for comparison.

*the fourth-order compact method at Nλ×Nϕ*= 512*×256 resolution, a less than 3% reduction*

*in C*_{2}*′* during the first 3 days is followed by a much stronger loss of the potential enstrophy
*to near 9% in at time t = 6. It can be seen that a slightly better conservation of potential*
enstrophy is obtained when the order of spatial differencing of each scheme is increased.

For potential palinstrophy, as the grid resolution increases and the solution converges, a
*rapid growth of disturbances in the unstable jet which reaches a maximum at time t = 6 days*
can be observed. This is in agreement with the results presented in figure 5. Furthermore, it
is observed that a slightly higher values of potential palinstrophy are obtained when the order
of spatial differencing of each scheme is increased.

**V1**
**P**
**o**
**te**
**n**
**ti**
**a**
**l **
**e**
**n**
**s**
**tr**
**o**
**p**
**h**
**y**
** (**
**re**
**la**
**ti**
**v**
**e**
**)**
0 1 2 3 4 5 6
-50
-45
-40
-35
-30
-25
-20
-15
-10
-5
0
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(a)**
**V1**
**P**
**o**
**te**
**n**
**ti**
**a**
**l **
**e**
**n**
**s**
**tr**
**o**
**p**
**h**
**y**
** (**
**re**
**la**
**ti**
**v**
**e**
**)**
0 1 2 3 4 5 6
-50
-45
-40
-35
-30
-25
-20
-15
-10
-5
0
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(b)**
**V1**
**P**
**o**
**te**
**n**
**ti**
**a**
**l **
**e**
**n**
**s**
**tr**
**o**
**p**
**h**
**y**
** (**
**re**
**la**
**ti**
**v**
**e**
**)**
0 1 2 3 4 5 6
-50
-45
-40
-35
-30
-25
-20
-15
-10
-5
0
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(c)**

Figure 7. Time evolution of the percentage of relative change in potential enstrophy for different
*methods at (a) Nλ× Nϕ* = 128*× 64, (b) Nλ× Nϕ* = 256*× 128 and (c) Nλ× Nϕ* = 512*× 256*

resolutions.
defined by
*E =*
⟨
1
2*h(u*
2* _{+ v}*2

_{) +}1 2

*gh*2 ⟩ (45)

*where h = H(1 + ˜h) and*

*⟨⟩ denotes the domain area average. Figure 9 shows a 20 days time*

*evolution of the percentage changes in E′defined by (E′(t)− E′(0)) /E′*(0) multiplied by 100,
*where E′* *= E−*⟨*c*2* _{H}*⟩

*⟨*

_{/2. The time-independent part}*2*

_{c}*⟩*

_{H}*to have a better representation of changes in total energy. It is seen that the sixth-order and*

_{/2 has been subtracted from E}*eighth-order CCFDM methods at time t = 20 days shows approximately 7 percent reduction*in total energy, the eighth-order SCFDM method exhibits approximately 9 percent reduction in total energy, the sixth-order SCFDM method shows approximately 10 percent reduction in

**V1**
**P**
**a**
**li**
**n**
**s**
**tr**
**o**
**p**
**h**
**y**
** (**
**re**
**la**
**ti**
**v**
**e**
**)**
0 1 2 3 4 5 6
-80
-60
-40
-20
0
20
40
60
80
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(a)**
**V1**
**P**
**a**
**li**
**n**
**s**
**tr**
**o**
**p**
**h**
**y**
** (**
**re**
**la**
**ti**
**v**
**e**
**)**
0 1 2 3 4 5 6
-80
-60
-40
-20
0
20
40
60
80
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(b)**
**V1**
**P**
**a**
**li**
**n**
**s**
**tr**
**o**
**p**
**h**
**y**
** (**
**re**
**la**
**ti**
**v**
**e**
**)**
0 1 2 3 4 5 6
-80
-60
-40
-20
0
20
40
60
80
**4th-order compact**
**6th-order CCFDM**
**6th-order SCFDM**
**8th-order CCFDM**
**8th-order SCFDM**
**(c)**

Figure 8. Time evolution of the percentage of relative change in potential palinstrophy for different
*methods at (a) Nλ×Nϕ*= 128*×64, (b) Nλ×Nϕ*= 256*×128 and (c) Nλ×Nϕ*= 512*×256 resolutions.*

total energy and the fourth-order compact method shows approximately 13 percent reduction in total energy. We conclude that the sixth-order and eighth-order CCFDM methods better conserve total energy compared to the other methods.

To summarize, the qualitative and quantitative assessment of the results indicate noticeable improvement of the results by using the sixth-order and the eighth-order SCFDM and CCFDM schemes over the fourth-order compact method in terms of global conservation and solution quality.

**Time (day)**
**C**
**h**
**a**
**n**
**g**
**e**
** i**
**n**
** t**
**o**
**ta**
**l **
**e**
**n**
**e**
**rg**
**y**
** (**
**re**
**la**
**ti**
**v**
**e**
**)**
0 5 10 15 20 25
-15
-10
-5
0
5
**4th-order Compact**
**6th-order SCFDM**
**6th-order CCFDM**
**8th-order SCFDM**
**8th-order CCFDM**

Figure 9. Time evolution of the percentage of relative change in total energy for different methods at
*Nλ× Nϕ*= 256*× 128 resolution.*

9. Conclusions

The focus of this work was on the application and qualitative and quantitative comparison of two families of high-order compact finite difference methods, namely, the super compact finite difference and the combined compact finite difference methods, to spatial differencing of the spherical shallow water equations in terms of vorticity, divergence and height on a regular latitude-longitude grid. They were compared with the standard fourth-order compact scheme. In addition, a semi-implicit Runge-Kutta method was developed for time advancing of the vorticity, divergence and height representation of the spherical shallow water equations.

Furthermore, to control the non-linear instability and avoiding the polar problem a spatial filter was developed. The assessment of the five high-order compact schemes was carried out for the test case developed in reference [38]. For numerical solution of the elliptic equations appeared in formulation of the problem the direct procedure proposed in references [33, 34] with minor changes was extended to the high-order compact finite difference schemes.

A study of the convergence rate showed that all spatial schemes convergence in agreement with the formal order of accuracy of the method. For qualitative assessment of the solutions the measurements of maximum vorticity gradient, power spectrum of total energy, relative change in potential enstrophy and potential palinstrophy were used.

Qualitative and quantitative measurements for the test case employed in this work, shows that using the the sixth-order and eighth-order CCFDM and the sixth-order and eighth-order SCFDM for spatial differencing of the spherical shallow water equations lead to a noticeable improvement of the solution over the fourth-order compact method. However, the performance of sixth-order and order CCFDM methods is better than the sixth-order and

eighth-order SCFDM methods.

Furthermore, the computational cost of the spatial differencing for finite difference methods used in this study is linearly proportional to the number of grid points, therefore, the methods in terms of computational cost from the lowest to highest are: the fourth-order compact, the sixth-order CCFDM, the order CCFDM, the sixth-order SCFDM and the eighth-order SCFDM. Therefore, in terms of qualitative and quantitative measurements including computational cost we conclude that the sixth-order and eighth-order CCFDM methods are superior to the sixth-order and eighth-order SCFDM methods for spatial differencing of the spherical shallow water equations.

ACKNOWLEDGEMENTS

Authors would like to thank university of Tehran, Link¨oping university and National Supercomputer Center in Sweden (NSC) for supporting this research work.

APPENDIX

I. Details of block tridiagonal system for numerical solution of elliptic equations
The matrix coefficients AS*j*, BS*j*, CS*j* **and right hand side vector D**S*j* of the block tridiagonal system

obtained from different compact method used in this study are presented here.
*I.1. Helmholtz equation*

*The matrix coefficients for the fourth-order compact method for 2 < j < Nϕ−1* are

BC4*j* =
*α*0*j* *− tan ϕ*2 *j* 1
3 0
2
*d*2 0
5
6
* , A*C4
*j* =
01 0 0
*2d*
1
6 0
*−1*
*d*2 0
1
12
* , C*C4
*j* =
*−1*0 0 0
*2d*
1
6 0
*−1*
*d*2 0
1
12
*where d is grid space and*

*αj*= *−m*
2
cos2_{ϕ}*j* *− a*
2
*βj* (46)

*and coefficients for j = 1 and j = Nϕ*are

BC41 =
*α*1 *− tan ϕ*1 1
(*−1)m*
*2d*
2
3+
(*−1)m+1*
6 0
*−(−1)m*
*d*2 +* _{d}*22 0
(

*−1)m*12 + 5 6

* , C*C4 1 =

*−1*0 0 0

*2d*1 6 0

*−1*

*d*2 0

_{12}1 BC4

*Nϕ*=

*αNϕ*

*− tan ϕNϕ*1

*−(−1)m*

*2d*2 3 + (

*−1)m+1*6 0

*−(−1)m*

*d*2 + 2

*d*2 0

*−(−1)*

*m*12 + 5 6

* , A*C4

*Nϕ*= 01 0 0

*2d*1 6 0

*−1*

*d*2 0 1 12

**and the right hand side vector D**

*j*

*for 1 < j < Nϕ*is

**D**C4*j* =
*a*
2_{R}_{ˆ}
*j*
0
0

*The matrix coefficients for the sixth-order SCFDM method for 2 < j < Nϕ−1*are
BSC6*j* =
*αj* *− tan ϕj* 1 0 0
0 *−2d* 0 *−*_{10}3*d*3 0
0 *−2* 0 *−*5_{6}*d*2 0
*−2* 0 *−d*2 0 *−*_{90}7*d*4
0 0 *−2* 0 *−*5_{6}*d*2
*, A*
SC6
*j* =
0 0 0 0 0
*−1 0 0 −d*3
60 0
0 1 0 *−d*_{12}2 0
1 0 0 0 *−d*4
360
0 0 1 0 *−d*2
12
CSC6*j* =
0 0 0 0 0
1 0 0 *−d*3
60 0
0 1 0 *−d*2
12 0
1 0 0 0 *−*_{360}*d*4
0 0 1 0 *−d*2
12
*and coefficients for j = 1 and j = Nϕ*are

BSC61 =
*α*1 *− tan ϕ*1 1 0 0
*−(−1)m* * _{−2d}*
0

*−3d*

_{10}3

*−*(

*−1)m+1*

_{60}

*d*3 0 0

*−2 + (−1)m+1*0

*−5d*

_{6}2

*−*(

*−1)m+1*

_{12}

*d*2 0

*−2 + (−1)m*

_{0}

*2*

_{−d}_{0}

*4 90*

_{−}7d*−*(

*−1)m*4 360 0 0

_{d}*−2 + (−1)m*0

*−5d*

_{6}2

*−*(

*−1)*

_{12}

*md*2 CSC61 = 0 0 0 0 0 1 0 0

*−d*3 60 0 0 1 0

*−d*

_{12}2 0 1 0 0 0

*−d*4 360 0 0 1 0

*−d*2 12 BSC6

*Nϕ*=

*αNϕ*

*− tan ϕNϕ*1 0 0 (

*−1)m*

*−2d*0

*−3d*

_{10}3

*−*(

*−1)m+1*

_{60}

*d*3 0 0

*−2 + (−1)m+1*0

*−5d*

_{6}2

*−*(

*−1)m+1*

_{12}

*d*2 0

*−2 + (−1)m*

_{0}

*2*

_{−d}_{0}

*4 90*

_{−}7d*−*(

*−1)m*4 360 0 0

_{d}*−2 + (−1)m*0

*−5d*

_{6}2

*−*(

*−1)*

_{12}

*md*2 ASC6

*N*= 0 0 0 0 0

_{ϕ}*−1 0 0 −d*3 60 0 0 1 0

*−d*2 12 0 1 0 0 0

*−d*4 360 0 0 1 0

*−d*

_{12}2

**and the right hand side vector D**

*j*

*for 1 < j < Nϕ*is

**D**SC6*j* =
*a*2*R*ˆ*j*
0
0
0
0

*The matrix coefficients for the eighth-order SCFDM method for 2 < j < Nϕ−1*are
BSC8*j* =
*αj* *− tan ϕj* 1 0 0 0 0
0 *−2d* 0 *−d*_{3}3 0 *−2d*_{126}5 0
0 *d* 0 *d*3
2 0
*7d*5
180 0
0 0 0 *d*3 0 *5d*_{12}5 0
*−2* 0 *−d*2 _{0} *−d*4
12 0 *−6d*
6
2240
0 0 *d*2 0 *d*_{2}4 0 *7d*_{180}6
0 0 0 0 *d*4 0 *5d*_{12}6
ASC8*j* =
0 0 0 0 0 0 0
*−1* 0 0 0 0 *−2d*_{5040}5 0
0 *−d*_{2} 0 0 0 _{720}*d*5 0
0 0 0 *−d*_{2}3 0 *d*5
24 0
1 0 0 0 0 0 _{40320}*−2d*6
0 0 *−d*_{2}2 0 0 0 _{720}*d*6
0 0 0 0 *−d*_{2}4 0 *d*6
24
CSC8*j* =
0 0 0 0 0 0 0
1 0 0 0 0 *−2d*_{5040}5 0
0 *−d*_{2} 0 0 0 _{720}*d*5 0
0 0 0 *−d*_{2}3 0 *d*5
24 0
1 0 0 0 0 0 _{40320}*−2d*6
0 0 *−d*_{2}2 0 0 0 _{720}*d*6
0 0 0 0 *−d*_{2}4 0 *d*6
24
*and coefficients for j = 1 and j = Nϕ*are

BSC81 =
*α*1 *− tan ϕ*1 1 0 0 0 0
*−(−1)m* _{−2d}_{0} *−d3*
3 0 *−2d5*126 *−*
(*−1)m+12d5*
5040 0
0 *d−*(*−1)m+1d*_{2} 0 *d3*_{2} 0 *7d5*_{180} +(*−1)m+1d5*_{720} 0
0 0 0 *d*3* _{−}*(

*−1)m+1d3*2 0

*5d5*12 + (

*−1)m+1*24 0

*−2 + (−1)m*

_{0}

*2*

_{−d}_{0}

*−d4*12 0

*−6d6*2240

*−*(

*−1)m2d6*40320 0 0

*d*2

*−*(

*−1)md2*2 0

*d4*2 0

*7d6*180+ (

*−1)md6*720 0 0 0 0

*d*4

*−*(

*−1)md4*

_{2}0

*5d6*

_{12}+(

*−1)m*

_{24} CSC81 = 0 0 0 0 0 0 0 1 0 0 0 0

*−2d*

_{5040}5 0 0

*−d*

_{2}0 0 0

*d*5 720 0 0 0 0

*−d*

_{2}3 0

*d*

_{24}5 0 1 0 0 0 0 0

_{40320}

*−2d*6 0 0

*−d*

_{2}2 0 0 0

*d*6 720 0 0 0 0

*−d*

_{2}4 0

*d*

_{24}6 BSC8

*Nϕ*=

*α*

_{Nϕ}*− tan ϕNϕ*1 0 0 0 0 (

*−1)m*

*−2d*0

*−d3*

_{3}0

*−2d5*

_{126}

*−*(

*−1)m+12d5*

_{5040}0 0

*d−*(

*−1)m+1d*2 0

*d3*2 0

*7d5*180 + (

*−1)m+1d5*720 0 0 0 0

*d*3

*−*(

*−1)m+1d3*2 0

*5d5*12 + (

*−1)m+1*24 0

*−2 + (−1)m*

_{0}

*2*

_{−d}_{0}

*−d4*12 0

*−6d6*2240

*−*(

*−1)m2d6*40320 0 0

*d*2

*−*(

*−1)md2*2 0

*d4*2 0

*7d6*180+ (

*−1)md6*720 0 0 0 0

*d*4

*(*

_{−}*−1)md4*2 0

*5d6*12 + (

*−1)m*24

ASC8*Nϕ* =
0 0 0 0 0 0 0
*−1* 0 0 0 0 *−2d*_{5040}5 0
0 *−d*_{2} 0 0 0 _{720}*d*5 0
0 0 0 *−d*_{2}3 0 *d*5
24 0
1 0 0 0 0 0 _{40320}*−2d*6
0 0 *−d*_{2}2 0 0 0 _{720}*d*6
0 0 0 0 *−d*_{2}4 0 *d*6
24
**and the right hand side vector D***j* *for 1 < j < Nϕ*is

**D**SC8*j* =
*a*2*R*ˆ*j*
0
0
0
0
0
0

*The matrix coefficients for the sixth-order CCFDM method for 2 < j < Nϕ−1* are

BCC6*j* =
*α*0*j* *− tan ϕ*1 *j* 10
6
*d*2 0 1
* , A*CC6
*j* =
150 0 0
*16d*
7
16
*d*
16
*−3*
*d*2 *−98d* *−1*8
* , C*CC6
*j* =
*−15*0 0 0
*16d*
7
16
*−d*
16
*−3*
*d*2
9
*8d* *−1*8
*and coefficients for j = 1 and j = Nϕ*are

BCC61 =
*α*1 *− tan ϕ*1 1
(*−1)m*15
*16d* 1 +
(*−1)m+1*7
16
(*−1)md*
16
6
*d*2 *−*
(*−1)m*3
*d*2 *−(−1)*
*m+1*_{9}
*8d* 1*−*
(*−1)m*
8
* , C*CC6
1 =
*−15*0 0 0
*16d*
7
16 *−d*16
*−3*
*d*2
9
*8d*
*−1*
8
BCC6*Nϕ* =
*αNϕ* *− tan ϕNϕ* 1
*−(−1)m*_{15}
*16d* 1 +
(*−1)m+1*_{7}
16
*−(−1)m _{d}*
16
6

*d*2

*−*(

*−1)m*

_{3}

*d*2 (

*−1)m+1*

_{9}

*8d*1

*−*(

*−1)m*8

* , A*CC6

*Nϕ*= 150 0 0

*16d*7 16

*d*16

*−3*

*d*2

*−98d*

*−1*8

*and the right hand side vector for 1 < j < Nϕ*is

**D**CC6*j* =
*a*
2_{R}_{ˆ}
*j*
0
0

*The matrix coefficients for the eighth-order CCFDM method for 2 < j < Nϕ−1*are

BCC8*j* =
*αj* *− tan ϕj* 1 0
0 1 0 0
8
*d*2 0 1 0
0 0 0 1
* , A*CC8
*j* =
0 0 0 0
35
*32d*
19
32
*d*
8
*d*2
96
*−4*
*d*2 *−29 _{16d}*

*−5*

_{16}

*−d*

_{48}

*−105*

*16d*3

*−10516d*2

*−158d*

*−3*16 CCC8

*j*= 0 0 0 0

*−35*

*32d*19 32

*−d*8

*d*2 96

*−4*

*d*2 29

*16d*

*−5*16

*d*48 105

*16d*3

*−10516d*2 15

*8d*

*−3*16

*and coefficients for j = 1 and j = Nϕ*are
BCC81 =
*α*1 *− tan ϕ*1 1 0
(*−1)m*35
*32d* 1 +
(*−1)m+1*19
32
(*−1)md*
8
(*−1)m+1d*2
96
8
*d*2 *−*
(*−1)m*_{4}
*d*2 *−(−1)*
*m+1*_{29}
*16d* 1*−*
(*−1)m*_{5}
16
*−(−1)m+1 _{d}*
48

*−(−1)m*

_{105}

*16d*3

*−(−1)*

*m+1*

_{105}

*16d*2

*−(−1)*

*m*

_{15}

*8d*1

*−*(

*−1)m+1*

_{3}16

*, C*CC8 1 = 0 0 0 0

*−35*

*32d*19 32

*−d*8

*d*2 96

*−4*

*d*2 29

*16d*

*−5*16

*d*48 105

*16d*3

*−10516d*2 15

*8d*

*−3*16 BCC8

*N*=

_{ϕ}*αNϕ*

*− tan ϕNϕ*1 0

*−(−1)m*

_{35}

*32d*1 + (

*−1)m+1*

_{19}32

*−(−1)m*8 (

_{d}*−1)m+1*2 96 8

_{d}*d*2

*−*(

*−1)m*

_{4}

*d*2 (

*−1)m+1*

_{29}

*16d*1

*−*(

*−1)m*

_{5}16 (

*−1)m+1*48 (

_{d}*−1)m*

_{105}

*16d*3

*−(−1)*

*m+1*

_{105}

*16d*2 (

*−1)m*

_{15}

*8d*1

*−*(

*−1)m+1*

_{3}16

*, A*CC8

*N*= 0 0 0 0 35

_{ϕ}*32d*19 32

*d*8

*d*2 96

*−4*

*d*2

*−2916d*

*−5*16

*−d*48

*−105*

*16d*3

*−105*2

_{16d}*−15*

_{8d}*−3*

_{16}

*and the right hand side vector for 1 < j < Nϕ*is

**D**CC8*j* =
*a*2*R*ˆ*j*
0
0
0

*I.2. Poisson equation*

*The matrix coefficients for the Poisson equation for nonzero wave numbers (m* *̸= 0) are found by*
*setting βj* *= 0 in (46). For wave number zero (m = 0) equation (31) is used and we need to modify*

the first row of BS*j* **and D**S*j*.

REFERENCES
*1. Pedlosky J. Geophysical fluid dynamics. Springer-Verlag, 1987.*

*2. Vallis G.K. Atmospheric and oceanic fluid dynamics: Fundamentals and large-scale circulation. Cambridge*
University Press, 2006.

*3. Kreiss H. O., Oliger J. Comparison of accurate methods for the integration of hyperbolic equations. Tellus*
**1972; 24:199–215.**

4. Hirsch S. R. Higher order accurate difference solutions of fluid mechanics problems by a compact differencing
**technique. J. Comput. Phys. 1975; 19:90–109.**

* 5. Lele S.K. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 1992; 103:16–42.*
6. Kreiss H.O., Scherer G. Finite element and finite difference methods for hyperbolic partial differential
equations, in: C. De Boor (Ed.), Mathematical Aspects of Finite Elements in Partial Differential Equation,
Academic Press, New York, 1974.

7. Carpenter M.H., Nordstr¨om J., Gottlieb D. A stable and conservative interface treatment of arbitrary
**spatial accuracy, J. Comput. Phys. 1999; 148:341-365.**

8. Mattsson K., Nordstr¨om J. Summation by parts operators for finite difference approximations of second
**derivatives, J. Comput. Phys. 2004; 199:503-540.**

* 9. Strand B. Summation by parts for finite difference approximation for d/dx, J. Comput. Phys. 1994; *
110:47-67.

*10. Durran, D.R. Numerical methods for fluid dynamics with applications to geophysics. Second Edition,*
Springer, 2010.

11. Navon I.M., Riphagen H.A. An implicit compact fourth-order algorithm for solving the shallow water
**equations in conservative-law form. Mon. Wea. Rev. 1979;107:1107–1127.**

*12. Chang H.H., Shirer N., Compact spatial differencing techniques in numerical modeling. Mon. Wea. Reu.*
**1985; 113:409–423.**

13. Spotz W.F., Taylor M.A., Swarztrauber P.N. Fast shallow-water solvers in latitude-longitude coordinates.
**Journal of Computational Physics, 1998; 145:432–444.**

*14. Blayo E. Compact finite difference schemes for ocean models: 1.ocean waves. Journal of Computational*
**Physics 2000; 164:241–257.**