Abstract

This paper shows that in the linearized shallow-water equations, the numerical schemes can become weakly unstable for the 2Δx wave in the C grid when the Courant number is 1 in the forward–backward scheme and 0.5 in the leapfrog scheme because of the repeated eigenvalues in the matrices. The instability can be amplified and spread to other waves and smaller Courant number if the diffusion term is included. However, Shuman smoothing can control the instability.

1. Introduction

The forward–backward in time and the leapfrog schemes have been applied to the shallow-water equations by Mesinger and Arakawa (1976), Haltiner and Williams (1980), and the internal gravity waves by Sun (1980, 1984) and many others. However, they can be unstable for 2Δx waves because of the repeated eigenvalues when Courant number (Co) is 0.5 in the leapfrog scheme (LF) or Courant number is 1 in the forward–backward schemes (FB) in C grid because of the existence of repeated eigenvalues. The weak instability of the LF for a simplified wave equation at 2Δx wave and Co = 0.5 has also been discussed by Durran (1999). For an unstaggered grid, the critical Courant number is 1 (Durran 1999). Here, we provide a detailed discussion on instability of the FB scheme at Co = 1. It is also found that the instability will be amplified and spread to the longer waves if the diffusion terms are added in both schemes. On the other hand, Shuman smoothing can be applied to control the instability for both schemes in the shallow-water equations.

2. Numerical schemes and eigenvalues

The 1D linearized shallow-water equations are

 
formula
 
formula

where H is the mean depth, g is gravity, h and u are depth perturbation and velocity, and ν is viscosity. The control parameter is δ = 0 or 1, because viscosity is frequently applied in the momentum equations, but may not be applied in the equation of mass field h (or density in the compress fluid; Versteeg and Malalasekera 1995).

In the staggered C grids, the finite-difference equations for the FB scheme become

 
formula

where q = p ± ½, Δx is the spatial interval between the p and p + 1 grid points. If a wave-type solution at the nth time step is assumed:

 
formula

and

 
formula

where i = , ĥn, and ûn are the amplitude of h and u at the nth time step; k is the wavenumber. The difference equations become

 
formula
 
formula

Which can be written as

 
formula

The eigenvalue of (6) is

 
formula

where X = sin(kΔx/2)/(Δx/2), R2 = gHΔt2X2 = (CXΔt)2, S2 = νΔtX2, and . For simplicity, we set g = 1, H = 1, and Δx = 1.

Following the same procedure, the difference equations for the LF scheme are

 
formula

Using (4a) and (4b), we obtain

 
formula
 
formula

The eigenvalue is

 
formula

or

 
formula

if δ = 1.

3. Results and discussion

a. Results without viscosity

Without viscosity (i.e., ν = 0), the FB is neutrally stable, |λ| = 1, as long as the Courant number Co = CΔtx < 1 according to (7); and the LF is also neutrally stable when Co = CΔtx < 0.5 according to (10a) (Mesinger and Arakawa 1976; Haltiner and Williams 1980).

However, the FB has repeated eigenvalues λ1,2 = /2 = −1 for 2Δx waves at Co = 1, because R2 = (CΔtX)2 = [CΔt sin(kΔx/2)/Δx/2]2 = 4. Equation (6) becomes

 
formula

An eigenvector corresponding to λ1 = −1 is . The generated eigenvector x2 can be found from

 
formula

and

 
formula

If we define a matrix , then, we can have

 
formula

which is in Jordan block form (Nobel 1969), and since , we obtain

 
formula
 
formula

Hence, the FB becomes weakly unstable, the magnitude linearly increases with time step and with 2Δt oscillation.

When Co = 0.5, (i.e., R = 1), and S2 = 0, (i.e., ν = 0), the LF scheme also has the repeated eigenvalues: λ1 = λ2 = −i and λ3 = λ4 = i. Hence, it becomes weakly unstable, as does the FB scheme. The weakly instability of the LF in a simplified wave equation for 2Δx wave at Co = 0.5 has been discussed by Durran (1999) and will not be repeated here.

b. Results with viscosity

Diffusion has been frequently applied to control the noise of the short waves in equations, because the amplification factor is λ = 1 − S2 = 1 − νΔtX2 when the forward-in-time and centered-in-space scheme is applied to the diffusion equation (Sun 1982). In the FB, viscous terms create unstable 2Δx and 3Δx waves at Co = 1, and unstable 2Δx wave at Co = 0.8. However, the growth rate becomes less than 1 for longer waves and the maximum damping occurs at 4Δx or 3Δx waves (Fig. 1).

Fig. 1.

Growth rate for forward–backward scheme as function of wavelength (Δx) and viscosity (ν) at (a) Co = 1 and (b) Co = 0.8.

Fig. 1.

Growth rate for forward–backward scheme as function of wavelength (Δx) and viscosity (ν) at (a) Co = 1 and (b) Co = 0.8.

In the LF, if R2 = 1, from (10b) we obtain

 
formula

Since the maximum value of |λ| is greater than 1 for S2 > 0, the FL is also unstable for R2 = 1, (i.e., 2Δx wave with Co = 0.5). The amplification factor for 2Δx wave = 2 for ν = 0.25 (i.e., S2 = 0.5), as shown in Fig. 2, because the viscous term enhances the 2Δt oscillation for the short waves (wavelength ≤ 3Δx). Instability also shows up for 2Δx wave when ν > 0.13 for Co = 0.4. Viscosity has the largest damping at 4Δx waves when Co = 1 in FB and Co = 0.5 in LF and at 3Δx waves when Co = 0.8 in FB and Co = 0.4 in LF. The effect of viscosity on FB and LF are similar. Direct integrating (5) and (9) confirms the viscous term amplifying instability of 2Δx or 3Δx waves. The instability becomes weaker with δ = 0.

Fig. 2.

Growth rate for leapfrog scheme as function of wavelength (Δx) and viscosity (ν) at (a) Co = 0.5 and (b) Co = 0.4.

Fig. 2.

Growth rate for leapfrog scheme as function of wavelength (Δx) and viscosity (ν) at (a) Co = 0.5 and (b) Co = 0.4.

c. Results with Shuman smoothing, but without viscosity

If we integrate (5) and (9) directly with g = 1, H = 1, (i.e., C = 1), and Δx = 1. The initial values of u and h are ĥ0 = 1.0 and û0 = 0.0 for different wavelengths, which represents the wave at the peak potential energy, but without kinetic energy initially and is a popular initial condition.

With ν = 0 and Co = 1, the time sequences of ĥn integrated from FB as function of wavelength are shown in Table 1; the values of ĥn integrated using LF with Co = 0.5 are shown in Table 2. In the leapfrog scheme, the first time step is integrated with a centered difference in space and forward in time. In addition to 2Δt oscillations, the magnitude of ĥn increases linearly with time for 2Δx wave in both schemes.

Table 1.

Time sequence of height perturbation as function of wavelength for forward–backward scheme (ν = 0, η = 0, and Co = 1).

Time sequence of height perturbation as function of wavelength for forward–backward scheme (ν = 0, η = 0, and Co = 1).
Time sequence of height perturbation as function of wavelength for forward–backward scheme (ν = 0, η = 0, and Co = 1).
Table 2.

Time sequence of height perturbation as function of wavelength for leapfrog scheme (ν = 0, η = 0, and Co = 0.5).

Time sequence of height perturbation as function of wavelength for leapfrog scheme (ν = 0, η = 0, and Co = 0.5).
Time sequence of height perturbation as function of wavelength for leapfrog scheme (ν = 0, η = 0, and Co = 0.5).

A simple fourth-order Shuman smoothing (Haltiner and Williams 1980) is applied to the prognostic variables u and h immediately after the new values are calculated at each time step. The damping factor is Sh = 1 − [2η sin2(kΔx/2)]2, where η is the smoothing coefficient. In the FB scheme, the eigenvalue λ in (14a) is replaced by λ′ = Sh × λ. The magnitude of n(λ′)n−1 starts decreasing when n becomes large. On the other hand, the magnitude of n(λ)n−1 linearly increases with n, as shown in Fig. 3. Hence, smoothing can remove the weak instability. The time sequences obtained from the direct integration with smoothing shown in Tables 3 and 4 confirm that smoothing (with η = 0.15) can control the instability of 2Δx wave in both FB and LF.

Fig. 3.

Magnitude of n(−1)n−1 (dashed line) and n(−0.91)n−1 (solid line) as a function of the nth time step, where λ′ = Sh × λ = 0.91 for η = 0.15.

Fig. 3.

Magnitude of n(−1)n−1 (dashed line) and n(−0.91)n−1 (solid line) as a function of the nth time step, where λ′ = Sh × λ = 0.91 for η = 0.15.

Table 3.

As in Table 1, but with smoothing (ν = 0, η = 0.15, and Co = 1).

As in Table 1, but with smoothing (ν = 0, η = 0.15, and Co = 1).
As in Table 1, but with smoothing (ν = 0, η = 0.15, and Co = 1).
Table 4.

As in Table 2, but with smoothing (ν = 0, η = 0.15, and Co = 0.5).

As in Table 2, but with smoothing (ν = 0, η = 0.15, and Co = 0.5).
As in Table 2, but with smoothing (ν = 0, η = 0.15, and Co = 0.5).

4. Summary

In shallow-water equations, the forward–backward and leapfrog schemes can become unstable for 2Δx wave at Co = 1 in the FB and Co = 0.5 in the LF because of the multiplicity of eigenvalue for 2Δx waves in both schemes. Instability will be amplified and spread to 3Δx waves by adding viscous terms, because it enhances the 2Δt oscillation, which creates unstable 2Δx and/or 3Δx waves for both schemes. However, Shuman smoothing can dampen the short waves instability in both schemes.

Acknowledgments

The author would like to thank Dr. M. Shieh, Dr. T. Oh, Dr. J. Lee, Ms. Y. C. Wang, and Ms. Sun for useful discussions and proofreading. The comments from the reviewers are also greatly appreciated.

REFERENCES

REFERENCES
Durran
,
D. R.
,
1999
:
Numerical Methods for Wave Equations in Geophysical Fluid Dynamics.
Springer, 465 pp
.
Haltiner
,
G. J.
, and
R. T.
Williams
,
1980
:
Numerical Prediction and Dynamic Meteorology.
2nd ed. Wiley, 477 pp
.
Mesinger
,
F.
, and
A.
Arakawa
,
1976
:
Numerical Methods Used in Atmospheric Models.
Vol. 1, GARP Publ. Series, Vol. 17, 64 pp
.
Nobel
,
B.
,
1969
:
Applied Linear Algebra.
Prentice-Hall, 523 pp
.
Sun
,
W. Y.
,
1980
:
A forward-backward time integration scheme to treat internal gravity waves.
Mon. Wea. Rev.
,
108
,
402
407
.
Sun
,
W. Y.
,
1982
:
A comparison of two explicit time integration schemes applied to the transient heat equation.
Mon. Wea. Rev.
,
110
,
1645
1652
.
Sun
,
W. Y.
,
1984
:
Numerical analysis for hydrostatic and nonhydrostatic equations of inertial-internal gravity waves.
Mon. Wea. Rev.
,
112
,
259
268
.
Versteeg
,
H. K.
, and
W.
Malalasekera
,
1995
:
An Introduction to Computational Fluid Dynamics: The Finite Volume Method.
Longman Scientific & Technical, 257 pp
.

Footnotes

Corresponding author address: Wen-Yih Sun, Department of Earth and Atmospheric Sciences, Purdue University, 550 Stadium Mall Dr. West, Lafayette, IN 47907-2051. Email: wysun@purdue.edu