## 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

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 *n*th time step is assumed:

and

where *i* = , *ĥ ^{n}*, and

*û*are the amplitude of

^{n}*h*and

*u*at the

*n*th 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), *R*^{2} = *gH*Δ*t*^{2}*X*^{2} = (*CX*Δ*t*)^{2}, *S*^{2} = *ν*Δ*tX*^{2}, 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

or

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 *R*^{2} = (*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 **x**_{2} can be found from

and

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 *S ^{2}* = 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 − *S*^{2} = 1 − *ν*Δ*tX*^{2} 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 *R ^{2}* = 1, from (10b) we obtain

Since the maximum value of |*λ*| is greater than 1 for *S*^{2} > 0, the FL is also unstable for *R*^{2} = 1, (i.e., 2Δ*x* wave with Co = 0.5). The amplification factor for 2Δ*x* wave = 2 for *ν* = 0.25 (i.e., *S*^{2} = 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

*ĥ*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Δ

^{n}*t*oscillations, the magnitude of

*ĥ*increases linearly with time for 2Δ

^{n}*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*η* sin^{2}(*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.

## 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

**,**

**,**

**,**

## 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