Stability Analysis of Numerical Methods for a 1.5-Layer Shallow-Water Ocean Model

Abstract

A 1.5-layer reduced-gravity shallow-water ocean model in spherical coordinates is described and discretized in a staggered grid (standard Arakawa C-grid) with the forward-time central-space (FTCS) method and the Leap-frog finite difference scheme. The discrete Fourier analysis method combined with the Gershgorin circle theorem is used to study the stability of these two finite difference numerical models. A series of necessary conditions of selection criteria for the time-space step sizes and model parameters are obtained. It is showed that these stability conditions are more accurate than the Courant-Friedrichs-Lewy (CFL) condition and other two criterions (Blumberg and Mellor, 1987; Casulli, 1990, 1992). Numerical experiments are proposed to test our stability results, and numerical model that is designed is also used to simulate the ocean current.

1. Introduction

The shallow-water model is a set of partial differential equations (PDEs), which derived from the principles of conservation of mass and conservation of momentum (the Navier-Stokes equations). Because the horizontal length scale is much greater than the vertical length scale, under this condition, the conservation of mass implies that the vertical velocity of the fluid is very small. It can be shown from the momentum equations that horizontal pressure gradients are due to the displacement of the pressure surface (or free surface) in a fluid, and that vertical pressure gradients are nearly hydrostatic [1]. The vertical integrating allows the vertical velocity to be removed from the equations; this is a classical derivation of the shallow-water system.

The situations in fluid dynamics where the horizontal length scale is much greater than the vertical length scale are very common; that is to say, the vertical acceleration of the fluid can be negligible. The flow of water over a free surface is a ubiquitous physical phenomenon that has aroused many scientists and engineers’ interest. For instance, if we consider the Coriolis forces in shallow-water model (the Coriolis term exists because we describe flows in a reference frame fixed on earth), this set of equations is particularly well suited for the study and numerical simulations of a large class of geophysical phenomena, such as atmospheric flows, ocean circulation, coastal flows, tides, tsunamis, and river and lake flows [2–10].

The shallow-water equations are derived from the Navier-Stokes equations that are nonlinear partial differential equations, which describe the motion of fluids. The nonlinearity makes most problems difficult or impossible to solve and it is the main contributor to the turbulence. Mathematicians and physicists believe that the turbulence can be found through an understanding of solutions to the Navier-Stokes equations. However, in the mathematical field, mathematicians have not yet proven that in three dimensions solutions always exist (existence), or that if they do exist, then they do not contain any singularity (that is smoothness). These are called the Navier-Stokes existence and smoothness problems; this is one of the seven most important open problems (the Millennium Prize Problems) in mathematics [11]. Therefore, it is also a challenge to make substantial progress toward the exact solution of shallow-water equations.

Research on numerical methods for the solution of the shallow-water system has attracted much attention; numerical simulation is an effective tool in solving them and a great variety of numerical methods have been developed to solve this system [12–22]. The numerical models which are based on the shallow-water equations, especially for finite difference numerical models, have been successfully applied to study the ocean circulation; for example, the low-frequency variability and bifurcation structure of wind-driven ocean circulation [5–7], the shallow-water model for the study of the Gulf Stream and its extension region [23–25], the Kuroshio current and its extension system [26–30], and so on. Although few people discuss the stability of numerical models, only several papers give the stability conditions of the simple formulation and linearized shallow-water equations [20, 31–36]. In this paper, we use the discrete Fourier analysis method and the Gerschgorin circle theorem to study the stability of the shallow-water numerical models and give a series of necessary conditions for the selection criteria of time step size.

The remainder of the paper is organized as follows. In Section 2, the brief description of 1.5-layer reduced-gravity shallow-water ocean model has been introduced. In Section 3, we use the FTCS method and the Leap-frog finite difference scheme to solve the shallow-water equations in a staggered grid. In Section 4, the discrete Fourier method combined with the Gerschgorin circle theorem is used to analyze the stability of these two numerical methods. In Section 5, numerical examples are given to test our results. The conclusions are given in Section 6.

2. Matematical Model

The mathematical derivation of the shallow-water equations can be found in many fluid dynamics books [37] and has already been presented by many authors. In our study, 1.5-layer shallow-water equations in spherical coordinates are nondimensionalized using the length scale 𝑟₀, which is radius of the earth, the mean depth of the upper layer 𝐻, a characteristic horizontal velocity scale 𝑈, a time scale 𝑟₀/𝑈, and a wind-stress scale 𝜏₀ (see [27, 38]).

The nondimensional equations are governed by the following reduced-gravity, nonlinear partial differential equations: 𝜀𝜕𝑢+𝑢𝜕𝑡cos𝜃𝜕𝑢𝜕𝜙+𝑣𝜕𝑢=𝜕𝜃−𝑢𝑣tan𝜃−𝑣sin𝜃−𝜀𝐹cos𝜃𝜕ℎ∇𝜕𝜙+𝐸²𝑢𝑢−cos²𝜃−2sin𝜃cos²𝜃𝜕𝑣𝜕𝜙+𝐹_𝜙,𝜀𝜕𝑣+𝑢𝜕𝑡cos𝜃𝜕𝑣𝜕𝜙+𝑣𝜕𝑣𝜕𝜃+𝑢²tan𝜃+𝑢sin𝜃=−𝜀𝐹𝜕ℎ∇𝜕𝜃+𝐸²𝑣𝑣−cos²𝜃+2sin𝜃cos²𝜃𝜕𝑢𝜕𝜙+𝐹_𝜃,𝜕ℎ+1𝜕𝑡𝜕()cos𝜃ℎ𝑢+𝜕()𝜕𝜙ℎ𝑣cos𝜃𝜕𝜃=0, (1)(2)(3)

in which 𝑢 is the zonal velocity, 𝑣 is the meridional velocity, 𝜙 is the coordinate in the zonal direction, 𝜃 is the coordinate in the meridional direction, and ℎ is the thickness of the upper ocean. The parameters in the equations are 𝜀=𝑈/2Ω𝑟₀, 𝐸=𝐴_𝐻/2Ω𝑟²₀, and 𝐹=𝑔^′𝐻/𝑈², where Ω is the rotation rate of the earth, 𝐴_𝐻 is the lateral friction coefficient, and 𝑔^′ is the reduced gravity.

The terms 𝐹_𝜙 and 𝐹_𝜃 in (1)-(2) are defined as follows: 𝐹_𝜙=𝛼𝜏_𝜙ℎ−𝜇𝑢,𝐹_𝜃=𝛼𝜏_𝜃ℎ−𝜇𝑣, (4)

where 𝛼=𝜏₀(2Ω𝜌𝐻𝑈) is wind stress coefficient (𝜏₀ is the amplitude of wind stress) and 𝜏_𝜙 and 𝜏_𝜃 are the zonal component and meridional component of wind stress, respectively. In addition, 𝜇=𝛾/(2Ω), where 𝛾 is the interfacial friction coefficient.

3. Finite Difference Schemes

If one discretizes the domain to a grid with equally spaced points with a spacing of Δ𝜙 in the 𝜙-direction, Δ𝜃 in the 𝜃-direction, and Δ𝑡 in the 𝑡-direction, we define 𝑢^𝑛_𝑖,𝑗=𝑢(𝜃_𝑖,𝜙_𝑗,𝑡_𝑛) with 𝜃_𝑖=𝜃₀+𝑖Δ𝜃, 𝜙_𝑗=𝜙₀+𝑗Δ𝜙, and 𝑡_𝑛=𝑛Δ𝑡 for 𝑖=0,1,…,𝑀, 𝑗=0,1…,𝑁, and 𝑛=1,2,…,𝑄, where 𝑀, 𝑁, and 𝑄 are positive integers. The variables 𝑢, 𝑣, and ℎ are evaluated at a staggered grid (standard Arakawa C-grid) as shown in Figure 1; then the shallow-water equations can be solved by using finite difference method.

Figure 1 

Staggered grid on space domain.

3.1. The FTCS Method

The forward difference approximation is used for the time derivative and the central difference approximation for the spatial derivatives. The difference approximation of (1)–(3) is given by 𝑢^𝑛+1_{𝑖+1/2,𝑗}−𝑢^𝑛_{𝑖+1/2,𝑗}+𝑢Δ𝑡^∗𝑢^𝑛_{𝑖+1/2,𝑗+1}−𝑢^𝑛_{𝑖+1/2,𝑗−1}+𝑣2Δ𝜙_∗𝑢^𝑛_{𝑖+3/2,𝑗}−𝑢^𝑛_{𝑖−1/2,𝑗}2Δ𝜃−𝑓^∗𝑣^𝑛_{𝑖,𝑗+1/2}𝐹=−^∗ℎ^𝑛_𝑖,𝑗+1−ℎ^𝑛_{𝑖,𝑗−1}2Δ𝜙+𝐸^∗𝑢^𝑛_{𝑖+1/2,𝑗+1}−2𝑢^𝑛_{𝑖+1/2,𝑗}+𝑢^𝑛_{𝑖+1/2,𝑗−1}Δ𝜙²+𝑢^𝑛_{𝑖+3/2,𝑗}−2𝑢^𝑛_{𝑖+1/2,𝑗}+𝑢^𝑛_{𝑖−1/2,𝑗}Δ𝜃²−𝐸^∗₁𝑣^𝑛_{𝑖,𝑗+3/2}−𝑣^𝑛_{𝑖,𝑗−1/2}2Δ𝜙+𝐻^∗𝜏_𝜙−Φ𝑢^𝑛_{𝑖+1/2,𝑗},𝑣^𝑛+1_{𝑖,𝑗+1/2}−𝑣^𝑛_{𝑖,𝑗+1/2}+𝑢Δ𝑡^∗𝑣^𝑛_{𝑖,𝑗+3/2}−𝑣^𝑛_{𝑖,𝑗−1/2}+𝑣2Δ𝜙_∗𝑣^𝑛_{𝑖+1,𝑗+1/2}−𝑣^𝑛_{𝑖−1,𝑗+1/2}2Δ𝜃+𝑓^∗𝑢^𝑛_{𝑖+1/2,𝑗}𝐹ℎ=−^{𝑛𝑖+1,𝑗}−ℎ^𝑛_{𝑖−1,𝑗}2Δ𝜃+𝐸^∗𝑣^𝑛_{𝑖,𝑗+3/2}−2𝑣^𝑛_{𝑖,𝑗+1/2}+𝑣^𝑛_{𝑖,𝑗−1/2}Δ𝜙²+𝑣^𝑛_{𝑖+1,𝑗+1/2}−2𝑣^𝑛_{𝑖,𝑗+1/2}+𝑣^𝑛_{𝑖−1,𝑗+1/2}Δ𝜃²+𝐸^∗₁𝑢^𝑛_{𝑖+1/2,𝑗+1}−𝑢^𝑛_{𝑖+1/2,𝑗−1}2Δ𝜙+𝐻^∗𝜏_𝜃−Φ𝑣^𝑛_{𝑖,𝑗+1/2},ℎ^𝑛+1_𝑖,𝑗−ℎ^𝑛_𝑖,𝑗+𝑢Δ𝑡^∗ℎ^𝑛_𝑖,𝑗+1−ℎ^𝑛_{𝑖,𝑗−1}+𝑣2Δ𝜙_∗ℎ^𝑛_𝑖+1,𝑗−ℎ^𝑛_{𝑖−1,𝑗}+ℎ2Δ𝜃^∗1𝑢^𝑛_{𝑖+1/2,𝑗+1}−𝑢^𝑛_{𝑖+1/2,𝑗−1}+ℎ2Δ𝜙_∗𝑣^𝑛_{𝑖+1,𝑗+1/2}−𝑣^𝑛_{𝑖−1,𝑗+1/2}2Δ𝜃−ℎ^∗₂𝑣^𝑛_{𝑖,𝑗+1/2}=0, (5)

in which 𝑢^∗=𝑢^𝑛_{𝑖+1/2,𝑗}/cos𝜃_𝑖+1/2,  𝑣^∗=𝑣^𝑛_{𝑖,𝑗+1/2}, 𝑓^∗=𝑢^𝑛_{𝑖+1/2,𝑗}tan𝜃_𝑖+1/2+(1/𝜀)sin𝜃_𝑖, 𝐹^∗=𝐹/cos𝜃_𝑖, 𝐸^∗=𝐸/𝜀, 𝐸^∗₁=2𝐸sin𝜃_𝑖/𝜀cos²𝜃_𝑖, ℎ^∗=ℎ^𝑛_𝑖,𝑗, 𝐻^∗=𝛼/𝜀ℎ^∗, Φ=(𝜇/𝜀)−(𝐸/𝜀cos²𝜃_𝑖+1/2), ℎ^∗₁=ℎ^∗/cos𝜃_𝑖, and ℎ^∗₂=ℎ^∗tan𝜃_𝑖.

3.2. The Leap-Frog Scheme

The Leap-frog differences are used for time derivatives and centered differences for space derivatives; the diffusion terms are lagged by one time step following the previous studies [33, 35]; we obtain 𝑢^𝑛+1_{𝑖+1/2,𝑗}−𝑢^𝑛−1_{𝑖+1/2,𝑗}+𝑢2Δ𝑡^∗𝑢^𝑛_{𝑖+1/2,𝑗+1}−𝑢^𝑛_{𝑖+1/2,𝑗−1}+𝑣2Δ𝜙_∗𝑢^𝑛_{𝑖+3/2,𝑗}−𝑢^𝑛_{𝑖−1/2,𝑗}2Δ𝜃−𝑓^∗𝑣^𝑛_{𝑖,𝑗+1/2}𝐹=−^∗ℎ^𝑛_𝑖,𝑗+1−ℎ^𝑛_{𝑖,𝑗−1}2Δ𝜙+𝐸^∗𝑢^𝑛−1_{𝑖+1/2,𝑗+1}−2𝑢^𝑛−1_{𝑖+1/2,𝑗}+𝑢^𝑛−1_{𝑖+1/2,𝑗−1}Δ𝜙²+𝑢^𝑛−1_{𝑖+3/2,𝑗}−2𝑢^𝑛−1_{𝑖+1/2,𝑗}+𝑢^𝑛−1_{𝑖−1/2,𝑗}Δ𝜃²−𝐸^∗₁𝑣^𝑛_{𝑖,𝑗+3/2}−𝑣^𝑛_{𝑖,𝑗−1/2}2Δ𝜙+𝐻^∗𝜏_𝜙−Φ𝑢^𝑛_{𝑖+1/2,𝑗},𝑣^𝑛+1_{𝑖,𝑗+1/2}−𝑣^𝑛−1_{𝑖,𝑗+1/2}+𝑢2Δ𝑡^∗𝑣^𝑛_{𝑖,𝑗+3/2}−𝑣^𝑛_{𝑖,𝑗−1/2}+𝑣2Δ𝜙_∗𝑣^𝑛_{𝑖+1,𝑗+1/2}−𝑣^𝑛_{𝑖−1,𝑗+1/2}2Δ𝜃+𝑓^∗𝑢^𝑛_{𝑖+1/2,𝑗}𝐹ℎ=−^{𝑛𝑖+1,𝑗}−ℎ^𝑛_{𝑖−1,𝑗}2Δ𝜃+𝐸^∗𝑣^𝑛−1_{𝑖,𝑗+3/2}−2𝑣^𝑛−1_{𝑖,𝑗+1/2}+𝑣^𝑛−1_{𝑖,𝑗−1/2}Δ𝜙²+𝑣^𝑛−1_{𝑖+1,𝑗+1/2}−2𝑣^𝑛−1_{𝑖,𝑗+1/2}+𝑣^𝑛−1_{𝑖−1,𝑗+1/2}Δ𝜃²+𝐸^∗₁𝑢^𝑛_{𝑖+1/2,𝑗+1}−𝑢^𝑛_{𝑖+1/2,𝑗−1}2Δ𝜙+𝐻^∗𝜏_𝜃−Φ𝑣^𝑛_{𝑖,𝑗+1/2},ℎ^𝑛+1_𝑖,𝑗−ℎ^𝑛−1_𝑖,𝑗+𝑢2Δ𝑡^∗ℎ^𝑛_𝑖,𝑗+1−ℎ^𝑛_{𝑖,𝑗−1}+𝑣2Δ𝜙_∗ℎ^𝑛_𝑖+1,𝑗−ℎ^𝑛_{𝑖−1,𝑗}+ℎ2Δ𝜃^∗1𝑢^𝑛_{𝑖+1/2,𝑗+1}−𝑢^𝑛_{𝑖+1/2,𝑗−1}+ℎ2Δ𝜙_∗𝑣^𝑛_{𝑖+1,𝑗+1/2}−𝑣^𝑛_{𝑖−1,𝑗+1/2}2Δ𝜃−ℎ^∗₂𝑣^𝑛_{𝑖,𝑗+1/2}=0. (6)

The Leap-frog method allows the direct calculation of 𝑢^𝑛+1_{𝑖+1/2,𝑗},  𝑣^𝑛+1_{𝑖,𝑗+1/2},  ℎ^𝑛+1_𝑖,𝑗 from the known values at time levels 𝑛 and 𝑛−1, which are the explicit difference equations.

4. Stability Analysis

In this section we present the stability conditions of the finite difference numerical models by using the discrete Fourier analysis method and the Gerschgorin circle theorem. The specific analysis of procedure is given in the following parts.

First of all, we set 𝐗^𝑛_𝑖,𝑗=(𝑢^𝑛_{𝑖+1/2,𝑗},𝑣^𝑛_{𝑖,𝑗+1/2},ℎ^𝑛_𝑖,𝑗)^𝑇, and after some rearrangement, (5) can be written as 𝐗^𝑛+1_𝑖,𝑗=𝐴₁𝐗^𝑛_𝑖+1,𝑗+𝐶₁𝐗^𝑛_{𝑖−1,𝑗}+𝐵₁𝐗^𝑛_𝑖,𝑗+𝐷₁𝐗^𝑛_𝑖,𝑗+1+𝐸₁𝐗^𝑛_{𝑖,𝑗−1}+𝐞_𝑛, (7)

in which 𝐞_𝑛=(Δ𝑡𝐻^∗𝜏_𝜙,Δ𝑡𝐻^∗𝜏_𝜃,0)^𝑇, 𝐴₁=−Δ𝑡𝑣^∗+2Δ𝜃Δ𝑡𝐸^∗Δ𝜃²000−Δ𝑡𝑣^∗+2Δ𝜃Δ𝑡𝐸^∗Δ𝜃²−Δ𝑡𝐹2Δ𝜃0−Δ𝑡ℎ^∗−2Δ𝜃Δ𝑡𝑣^∗,𝐶2Δ𝜃₁=Δ𝑡𝑣^∗+2Δ𝜃Δ𝑡𝐸^∗Δ𝜃²000Δ𝑡𝑣^∗+2Δ𝜃Δ𝑡𝐸^∗Δ𝜃²Δ𝑡𝐹02Δ𝜃Δ𝑡ℎ^∗2Δ𝜃Δ𝑡𝑣^∗,𝐵2Δ𝜃₁=1−Δ𝑡Φ−2Δ𝑡𝐸^∗1Δ𝜙²+1Δ𝜃²Δ𝑡𝑓^∗0−Δ𝑡𝑓^∗1−Δ𝑡Φ−2Δ𝑡𝐸^∗1Δ𝜙²+1Δ𝜃²00Δ𝑡ℎ^∗₂1,𝐷₁=−Δ𝑡𝑢^∗+2Δ𝜙Δ𝑡𝐸^∗Δ𝜙²−Δ𝑡𝐸^∗₁−2Δ𝜙Δ𝑡𝐹^∗2Δ𝜙Δ𝑡𝐸^∗₁−2Δ𝜙Δ𝑡𝑢^∗+2Δ𝜙Δ𝑡𝐸^∗Δ𝜙²0−Δ𝑡ℎ^∗₁2Δ𝜙0−Δ𝑡𝑢^∗,𝐸2Δ𝜙₁=Δ𝑡𝑢^∗+2Δ𝜙Δ𝑡𝐸^∗Δ𝜙²Δ𝑡𝐸^∗₁2Δ𝜙Δ𝑡𝐹^∗−2Δ𝜙Δ𝑡𝐸^∗₁2Δ𝜙Δ𝑡𝑢^∗+2Δ𝜙Δ𝑡𝐸^∗Δ𝜙²0Δ𝑡ℎ^∗₁02Δ𝜙Δ𝑡𝑢^∗.2Δ𝜙 (8)

In a similar fashion, (6) can be written as the following matrix equations: 𝐗^𝑛+1_𝑖,𝑗=𝐴₂𝐗^𝑛_{𝑖−1,𝑗}+𝐵₂𝐗^𝑛_𝑖+1,𝑗+𝐶₂𝐗^𝑛_{𝑖,𝑗−1}+𝐷₂𝐗^𝑛_𝑖,𝑗+1+𝐸₂𝐗^𝑛_𝑖,𝑗+𝐹₂𝐗^𝑛−1_𝑖+1,𝑗+𝐻₂𝐗^𝑛−1_{𝑖−1,𝑗}+𝑀₂𝐗^𝑛−1_𝑖,𝑗+1+𝑁₂𝐗^𝑛−1_{𝑖,𝑗−1}+𝑃₂𝐗^𝑛−1_𝑖,𝑗+𝐞_𝑛, (9)

where 𝐞_𝑛=(2Δ𝑡𝐻^∗𝜏_𝜙,2Δ𝑡𝐻^∗𝜏_𝜃,0)^𝑇, 𝐴₂=−𝐵₂=Δ𝑡𝑣^∗0Δ𝜃00Δ𝑡𝑣^∗Δ𝜃Δ𝑡𝐹0Δ𝜃Δ𝑡ℎ^∗Δ𝜃Δ𝑡𝑣^∗,𝐶Δ𝜃₂=−𝐷₂=Δ𝑡𝑢^∗Δ𝜙Δ𝑡𝐸^∗₁Δ𝜙Δ𝑡𝐹^∗−Δ𝜙Δ𝑡𝐸^∗₁Δ𝜙Δ𝑡𝑢^∗0Δ𝜙Δ𝑡ℎ^∗₁0Δ𝜙Δ𝑡𝑢^∗,𝐸Δ𝜙₂=−2Δ𝑡Φ2Δ𝑡𝑓^∗0−2Δ𝑡𝑓^∗−2Δ𝑡Φ002Δ𝑡ℎ^∗₂0,𝐹₂=𝐻₂=2Δ𝑡𝐸^∗Δ𝜃²0002Δ𝑡𝐸^∗Δ𝜃²0,𝑀000₂=𝑁₂=2Δ𝑡𝐸^∗Δ𝜙²0002Δ𝑡𝐸^∗Δ𝜙²0,𝑃000₂=1−4Δ𝑡𝐸^∗1Δ𝜙²+1Δ𝜃²0001−4Δ𝑡𝐸^∗1Δ𝜙²+1Δ𝜃²0.001 (10)

Definition 1. The two-dimensional discrete Fourier transform of 𝑢∈ℓ₂ is the function 𝑢^†∈𝐿₂([−𝜋,𝜋]×[−𝜋,𝜋]) defined by (see [39]) 𝑢^†=1𝜉,𝛽2𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝑢_𝑚,𝑘 (11)
for 𝜉,𝛽∈[−𝜋,𝜋].

We begin by taking the discrete Fourier transform of both sides of (7), and we can obtain 𝐗^†𝑛+1=1𝜉,𝛽2𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐗^𝑛+1_𝑚,𝑘=12𝜋^∞_{𝑚,𝑘=−∞}×𝐴exp−𝑖𝑚𝜉+𝑘𝛽₁𝐗^𝑛_𝑚+1,𝑘+𝐶₁𝐗^𝑛_{𝑚−1,𝑘}+𝐵₁𝐗^𝑛_𝑚,𝑘+𝐷₁𝐗^𝑛_𝑚,𝑘+1+𝐸₁𝐗^𝑛_{𝑚,𝑘−1}+𝐞_𝑛=𝐴₁2𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐗^𝑛_𝑚+1,𝑘+𝐶₁2𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐗^𝑛_{𝑚−1,𝑘}+𝐵₁2𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐗^𝑛_𝑚,𝑘+𝐷₁2𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐗^𝑛_𝑚,𝑘+1+𝐸₁2𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐗^𝑛_{𝑚,𝑘−1}+12𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐞_𝑛. (12)

By making the change of variables 𝑙=𝑚±1 we get 12𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐗^𝑛_𝑚±1,𝑘=12𝜋^∞_{𝑙,𝑘=−∞}exp−𝑖(𝑙∓1)𝜉+𝑘𝛽𝐗^𝑛_𝑙,𝑘()1=exp±𝑖𝜉2𝜋^∞_{𝑙,𝑘=−∞}exp−𝑖𝑙𝜉+𝑘𝛽𝐗^𝑛_𝑙,𝑘=exp(±𝑖𝜉)𝐗^†𝑛𝜉,𝛽. (13)

Similarly, we have 12𝜋^∞_{𝑚,𝑘=−∞}exp−𝑖𝑚𝜉+𝑘𝛽𝐗^𝑛_𝑚,𝑘±1=exp±𝑖𝛽𝐗^†𝑛𝜉,𝛽. (14)

Thus using the expressions (13)-(14) in (12) leads to 𝐗^†𝑛+1𝜉,𝛽=𝐺₁𝜉,𝛽𝐗^†𝑛𝜉,𝛽+𝐞^′_𝑛, (15)

where 𝐞^′_𝑛=(1/2𝜋)∑^∞_{𝑚,𝑘=−∞}exp(−𝑖(𝑚𝜉+𝑘𝛽))𝐞_𝑛, and the growth matrix 𝐺₁𝜉,𝛽=𝐴₁exp(𝑖𝜉)+𝐶₁exp(−𝑖𝜉)+𝐵₁+𝐷₁exp𝑖𝛽+𝐸₁=𝑎exp−𝑖𝛽₁₁𝑎₁₂𝑎₁₃𝑎₂₁𝑎₂₂𝑎₂₃𝑎₃₁𝑎₃₂𝑎₃₃, (16)

in which 𝑎₁₁=𝑎₂₂=1−Δ𝑡Φ−2Δ𝑡𝐸^∗1Δ𝜙²+1Δ𝜃²+2Δ𝑡𝐸^∗cos𝛽Δ𝜙²+cos𝜉Δ𝜃²𝑢−𝑖Δ𝑡^∗sin𝛽+𝑣Δ𝜙^∗sin𝜉,𝑎Δ𝜃₁₂=Δ𝑡𝑓^∗−𝑖Δ𝑡𝐸^∗₁Δ𝜙sin𝛽,𝑎₁₃=−𝑖Δ𝑡𝐹^∗𝑎Δ𝜙sin𝛽,₂₁=−Δ𝑡𝑓^∗+𝑖Δ𝑡𝐸^∗₁Δ𝜙sin𝛽,𝑎₂₃=−𝑖Δ𝑡𝐹𝑎Δ𝜃sin𝜉,₃₁=−𝑖Δ𝑡ℎ^∗₁Δ𝜙sin𝛽,𝑎₃₂=Δ𝑡ℎ^∗₂−𝑖Δ𝑡ℎ^∗𝑎Δ𝜃sin𝜉,₃₃𝑢=1−𝑖Δ𝑡^∗𝑣Δ𝜙sin𝛽+^∗.Δ𝜃sin𝜉 (17)

Remark 2. The discrete Fourier transform is used to deal with the vector 𝐗^𝑛_𝑖,𝑗=(𝑢^𝑛_{𝑖+1/2,𝑗},𝑣^𝑛_{𝑖,𝑗+1/2},ℎ^𝑛_𝑖,𝑗)^𝑇. Therefore, the individual variables 𝑢^∗, 𝑣^∗, and ℎ^∗ in the coefficient matrix of expression (7) can be treated as constants. Consequently, the coefficient matrixes can be extracted from the Fourier transform in the expression (12). These are similar to the frozen coefficients approach for discussing the stability of numerical solution of variable coefficient partial differential equation (see [40–42]).

Definition 3. The three-dimensional discrete Fourier transform of 𝑢∈ℓ₂ is the function 𝑢^†∈𝐿₂([−𝜋,𝜋]×[−𝜋,𝜋]×[−𝜋,𝜋]) defined by 𝑢^†=1𝜉,𝛽,𝜑2𝜋^∞_{𝑚,𝑘,𝑛=−∞}exp−𝑖𝑚𝜉+𝑘𝛽+𝑛𝜑𝑢^𝑛_𝑚,𝑘 (18)
for 𝜉,𝛽,𝜑∈[−𝜋,𝜋].

Taking the discrete Fourier transform (18) to the both sides of (9), similar to the derivation of (15), we have 𝐗^†𝑛+1𝜉,𝛽,𝜑=𝐺₂𝜉,𝛽,𝜑𝐗^†𝑛𝜉,𝛽,𝜑+𝐞^′_𝑛, (19)

where the growth matrix 𝐺₂𝜉,𝜃,𝛽=𝐴₂exp(−𝑖𝜉)+𝐵₂exp(𝑖𝜉)+𝐶₂exp−𝑖𝛽+𝐷₂exp𝑖𝛽+𝐸₂+𝐹₂exp𝑖𝜉−𝜑+𝐻₂exp−𝑖𝜉+𝜑+𝑀₂exp𝑖𝛽−𝜑+𝑁₂exp−𝑖𝛽+𝜑+𝑃₂=𝑏exp−𝑖𝜑₁₁𝑏₁₂𝑏₁₃𝑏₂₁𝑏₂₂𝑏₂₃𝑏₃₁𝑏₃₂𝑏₃₃, (20)

in which 𝑏₁₁=𝑏₂₂𝑢=−2Δ𝑡Φ−2𝑖Δ𝑡^∗𝑣Δ𝜙sin𝛽+^∗+Δ𝜃sin𝜉1−4Δ𝑡𝐸^∗1Δ𝜙²+1Δ𝜃²+exp−𝑖𝜑2Δ𝑡𝐸^∗Δ𝜙²𝑖+exp𝛽−𝜑+exp−𝑖𝛽+𝜑2Δ𝑡𝐸^∗Δ𝜃²𝑖𝑏exp𝜉−𝜑+exp−𝑖𝜉+𝜑₁₂=2Δ𝑡𝑓^∗−2𝑖Δ𝑡𝐸^∗₁𝑏Δ𝜙sin𝛽,₁₃=−2𝑖Δ𝑡𝐹^∗Δ𝜙sin𝛽,𝑏₂₁=−2Δ𝑡𝑓^∗+2𝑖Δ𝑡𝐸^∗₁𝑏Δ𝜙sin𝛽,₂₃=−2𝑖Δ𝑡𝐹Δ𝜃sin𝜉,𝑏₃₁=−2𝑖Δ𝑡ℎ^∗₁𝑏Δ𝜙sin𝛽,₃₂=2Δ𝑡ℎ^∗₂−2𝑖Δ𝑡ℎ^∗𝑏Δ𝜃sin𝜉,₃₃𝑢=exp−𝑖𝜑−2𝑖Δ𝑡^∗𝑣Δ𝜙sin𝛽+^∗.Δ𝜃sin𝜉 (21)

Theorem 4. The difference schemes (7) and (9) are stable in the ℓ₂ norm; then there exist constants Δ𝑥₀, Δ𝑦₀, Δ𝑡₀, and 𝐶, independent of Δ𝑥, Δ𝑦, and Δ𝑡, so that 𝜌(𝐺)≤1+𝐶Δ𝑡 (22)
for 0<Δ𝑡≤Δ𝑡₀, and 0<Δ𝑥≤Δ𝑥₀, 0<Δ𝑦≤Δ𝑦₀, and for all 𝜉,𝜃,𝜑∈[−𝜋,𝜋]. 𝜌(𝐺) is the spectral radius of the matrix 𝐺 (𝐺₁ and 𝐺₂) [39].

Theorem 5 (Gerschgorin circle theorem, see [39]). Suppose 𝐴=(𝑎_𝑖𝑗) is a general 𝑛×𝑛 matrix, and 𝜌_𝑠=∑^𝑛_{𝑗=1,𝑗≠𝑠}|𝑎_𝑠𝑗| is the sum of the absolute values of the elements in the 𝑠th row except for the diagonal element. For each eigenvalue 𝜆 of 𝐴, there exists an 𝑠 such that 𝜆−𝑎_𝑠𝑠≤𝜌_𝑠,𝑠=1,2,…𝑛. (23)
If 𝜆 is an eigenvalue of 𝐺₁, by using Gerschgorin circle theorem and the triangular inequality, one has |𝜆|−𝑎₁₁≤𝑎₁₂+𝑎₁₃,|𝜆|−𝑎₂₂≤𝑎₂₁+𝑎₂₃,|𝜆|−𝑎₃₃≤𝑎₃₁+𝑎₃₂. (24)
Then the time-space steps size and model parameters yield the following conditions (one takes 𝐶≤𝜖/Δ𝑡, where 0<𝜖≤1): 𝜖Δ𝑡≤𝑓^∗+|Φ|𝜖,Δ𝑡≤ℎ^∗₂,𝑓𝜉=𝛽=0Δ𝑡^∗+|Φ|+4𝐸^∗Δ𝜙²+4𝐸^∗Δ𝜃²≤𝜖𝜉=𝛽=𝜋,Δ𝑡𝑢Δ𝜙^∗+𝐸^∗₁+𝐹^∗+Δ𝑡𝑣^∗|Φ|+𝑓Δ𝜃+Δ𝑡^∗+2Δ𝑡𝐸^∗Δ𝜙²+2Δ𝑡𝐸^∗Δ𝜃²≤𝜖,Δ𝑡𝑢Δ𝜙^∗+𝐸^∗₁+Δ𝑡𝑣Δ𝜃^∗+|𝐹|+Δ𝑡|Φ|+𝑓^∗+2Δ𝑡𝐸^∗Δ𝜙²+2Δ𝑡𝐸^∗Δ𝜃²≤𝜖,Δ𝑡𝑢Δ𝜙^∗+ℎ^∗₁+Δ𝑡𝑣Δ𝜃^∗+ℎ^∗+Δ𝑡ℎ^∗₂𝜋≤𝜖𝜉=𝛽=2. (25)(26)(27)
Therefore, one has |𝜆|≤1+𝐶Δ𝑡, and the FTCS scheme is conditionally stable. The condition is as follows: 𝜖Δ𝑡≤min𝑓^∗+|Φ|,𝜖ℎ^∗₂,𝜖𝛿₁,𝜖𝛿₂,𝜀𝛿₃,𝜖𝛿₄, (28)
in which 𝛿₁=|𝑓^∗|+|Φ|+(4|𝐸^∗|/Δ𝜙²)+(4|𝐸^∗|/Δ𝜃²), 𝛿₂=(1/Δ𝜙)(|𝑢^∗|+|𝐸^∗₁|+|𝐹^∗|)+(|𝑣^∗|/Δ𝜃)+|Φ|+|𝑓^∗|+(2|𝐸^∗|/Δ𝜙²)+(2|𝐸^∗|/Δ𝜃²),𝛿₃=(1/Δ𝜙)(|𝑢^∗|+|𝐸^∗₁|)+(1/Δ𝜃)(|𝑣^∗|+|𝐹|)+|Φ|+|𝑓^∗|+(2|𝐸^∗|/Δ𝜙²)+(2|𝐸^∗|/Δ𝜃²), and 𝛿₄=(1/Δ𝜙)(|𝑢^∗|+|ℎ^∗₁|)+(1/Δ𝜃)(|𝑣^∗|+|ℎ^∗|)+|ℎ^∗₂|.
Similarly, one obtains the stable condition of the explicit Leap-frog finite-difference scheme (9): 𝜖Δ𝑡≤min2𝑓^∗+|Φ|,𝜖2ℎ^∗₂,𝜖2𝛿₁,𝜖2𝛿₂,𝜖2𝛿₃,𝜖2𝛿₄, (29)
where 𝛿₁, 𝛿₂, 𝛿₃, and 𝛿₄ are the same as in expression (28).

Remark 6. The derivations of stability conclusions in this study are still valid for both A-grid and B-grid; the results depend mainly on the choice of vector 𝐗^𝑛_𝑖,𝑗; for example, in the A-grid, we take 𝐗^𝑛_𝑖,𝑗=(𝑢^𝑛_𝑖,𝑗,𝑣^𝑛_𝑖,𝑗,ℎ^𝑛_𝑖,𝑗)^𝑇. In addition, it is easy to prove that the stability conditions derived from C-grid are the same for both A-grid and B-grid.
As a matter of fact, when the rotation, eddy viscosity, wind stress, and interfacial friction are neglected, the second expression in (27) can be written as (𝜖=1) Δ𝑡𝑢^∗+𝑣Δ𝜙Δ𝑡^∗+|𝐹|Δ𝜃Δ𝑡Δ𝜃≤1. (30)
This is the Courant-Friedrichs-Lewy condition (CFL condition; see [13, 34, 42]) in two-dimensional case. Assuming the terms |𝑓^∗|=0, |Φ|=0 and 𝜖=1 in the expression (26), then we have 4Δ𝑡𝐸^∗Δ𝜙²+4Δ𝑡𝐸^∗Δ𝜃²≤1. (31)
In fact, this is the same as the conditions identified by Blumberg and Mellor in 1987 [43]. When the rotation, wind stress and interfacial friction terms are neglected and set 𝜖=1, the first and second expressions in (27) are given as Δ𝑡𝑢^∗+𝑣Δ𝜙Δ𝑡^∗+1Δ𝜃2Δ𝑡𝐸^∗Δ𝜙²+Δ𝑡𝐸^∗Δ𝜃²≤1. (32)
This condition is the same as [15, 20] that given by Casulli and Cheng. The stability criteria (30)–(32) have been widely applied to the numerical model for the selection of the time-step size. However, these three conditions are only special cases in our results.

5. Numerical Experiments

In this section, numerical examples are given to test our results. In the present study, we take the FTCS scheme, for example (because the stability criterions of the Leap-frog finite-difference scheme are similar to the FTCS scheme). The domain of integration is set as a part of the North Pacific basin (25°–35°N, 132°–140°E). We use a realistic coastline and the 200 m depth contour as the continental boundary [27]. The horizontal resolution is 0.2° × 0.2°; that is, the space-step size Δ𝜃=Δ𝜙=0.2. Standard parameter values in the shallow-water model are shown in Table 1.

Table 1 

The standard values of parameters in the model.

5.1. Example  1

We take 𝜖=0.3, the zonal velocity 𝑢=0.3, and the meridional velocity 𝑣=0.3. According to the expression (28), it is easy to obtain Δ𝑡<3/233020; multiplying by the time scale, we have Δ𝑡<821 s. In the light of the CFL condition (30), we have Δ𝑡<5790 s. With the expression (31), we obtain Δ𝑡<4.6×10⁸ s. From the stability criterion (32), we get Δ𝑡<2.2×10⁷ s.

Case 1. Setting a time-step size Δ𝑡=825 s, after running the model with 30 steps (55/8 hours), Figure 2 gives the values (the dimensionless quantity, as well as the following results) of the zonal velocity 𝑢 and the meridional velocity 𝑣 at the latitude 30°N; Figure 3 gives the results at the longitude 135°E. It is not difficult to find that the current velocities 𝑢 and 𝑣 are not in accord with the actual condition of ocean.

(a)

(b)

Figure 2 

The figures show the values of the zonal velocity 𝑢 and the meridional velocity 𝑣 at the latitude 30°N after running the model with 30 steps.

(a)

(b)

Figure 3 

The figures show the values of the zonal velocity 𝑢 and the meridional velocity 𝑣 at the longitude 135°E after running the model with 30 steps.

Case 2. When we choose the step size Δ𝑡=1200 s, the results in Figure 4 give the values of the zonal velocity 𝑢 and the meridional velocity 𝑣 at the latitude 29°N after a time integration of 5 hours. Obviously, the results are also unreasonable. Moreover, after continuing the calculation of model, we find that the results start to overflow after running the model with 17 steps (17/3 hours).

(a)

(b)

Figure 4 

The figures show the values of the zonal velocity 𝑢 and the meridional velocity 𝑣 at the latitude 29°N after a time integration of 5 hours.

Case 3. The model is run with a time step size of Δ𝑡=300 s, Figures 5 and 6 show the values of the zonal velocity 𝑢, the meridional velocity 𝑣, and the layer thickness ℎ after a time integration of 8 hours and 8 years, respectively. These results illustrate that the model is integrated for long periods of time, and the results are still reasonable.

(a)

(b)

(c)

Figure 5 

The figures show the values of the zonal velocity 𝑢, the meridional velocity 𝑣, and the layer thickness ℎ after a time integration of 8 hours.

(a)

(b)

(c)

Figure 6 

The figures show the values of the zonal velocity 𝑢, the meridional velocity 𝑣, and the layer thickness ℎ after a time integration of 8 years.

It is obvious that the stability condition (28) is reasonable, because the numerical model is unstable when we take the time step size Δ𝑡>821 s (as shown in Cases 1 and 2). On the other hand, it is easy to see that our results are more strict and accurate than the CFL condition (30) and other two criterions (31) and (32).

5.2. Example  2

In this example, the 1.5-layer shallow-water numerical model that is designed by us is used to simulate the ocean current. Based on Example  1, the ocean basin is also adopted with the part of the North Pacific basin (25°–35°N, 132°–140°E). The time-space step sizes and the standard values of parameters in the model are the same as Case 3 in Example  1. Figure 7 gives the meridional velocity of the ocean current throughout a period of the time-dependent solution that evolves in one day. We will make an attempt to use this explicit shallow-water numerical model to simulate the Kuroshio current and its extension system in further studies.

Figure 7 

The figures show the meridional velocity of the ocean current throughout a period of the time-dependent solution that evolves in one day.

6. Conclusions

The FTCS and the Leap-frog finite difference scheme for solving 1.5-layer shallow-water equations in spherical coordinates have been presented. The stability conditions of these two types of difference schemes are given, which include the CFL condition and other two criterions [15, 20, 43]. The numerical experiments are proposed for testing the stability of the FTCS scheme; the numerical results illustrate that our stability conditions are effective and reasonable. Moreover, the present stability criterion is shown to be more accurate than other criterions that this research mentioned. The theory of stability analysis in this paper can also be used to study the complex coupled atmosphere-ocean models.

Acknowledgments

This study was provided by the National Basic Research Program of China (Grant no. 2012CB417404), the National Nature Scientific Foundation of China (41230420), and the Basic Research Program of Qingdao Science and Technology Plan (11-1-4-95-jch).