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.
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
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
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:
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
Which can be written as
The eigenvalue of (6) is
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
The eigenvalue is
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Δt/Δx < 1 according to (7); and the LF is also neutrally stable when Co = CΔt/Δx < 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
An eigenvector corresponding to λ1 = −1 is . The generated eigenvector x2 can be found from
If we define a matrix , then, we can have
which is in Jordan block form (Nobel 1969), and since , we obtain
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).
In the LF, if R2 = 1, from (10b) we obtain
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.
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.
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.
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.
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.
Corresponding author address: Wen-Yih Sun, Department of Earth and Atmospheric Sciences, Purdue University, 550 Stadium Mall Dr. West, Lafayette, IN 47907-2051. Email: email@example.com