Research Article  Open Access
Stability Analysis of Numerical Methods for a 1.5Layer ShallowWater Ocean Model
Abstract
A 1.5layer reducedgravity shallowwater ocean model in spherical coordinates is described and discretized in a staggered grid (standard Arakawa Cgrid) with the forwardtime centralspace (FTCS) method and the Leapfrog 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 timespace step sizes and model parameters are obtained. It is showed that these stability conditions are more accurate than the CourantFriedrichsLewy (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 shallowwater model is a set of partial differential equations (PDEs), which derived from the principles of conservation of mass and conservation of momentum (the NavierStokes 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 shallowwater 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 shallowwater 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 shallowwater equations are derived from the NavierStokes 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 NavierStokes 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 NavierStokes 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 shallowwater equations.
Research on numerical methods for the solution of the shallowwater 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 shallowwater equations, especially for finite difference numerical models, have been successfully applied to study the ocean circulation; for example, the lowfrequency variability and bifurcation structure of winddriven ocean circulation [5–7], the shallowwater 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 shallowwater 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 shallowwater 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.5layer reducedgravity shallowwater ocean model has been introduced. In Section 3, we use the FTCS method and the Leapfrog finite difference scheme to solve the shallowwater 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 shallowwater equations can be found in many fluid dynamics books [37] and has already been presented by many authors. In our study, 1.5layer shallowwater 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 windstress scale (see [27, 38]).
The nondimensional equations are governed by the following reducedgravity, nonlinear partial differential equations:
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 , , 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:
where 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, , 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 , , and , where , , and are positive integers. The variables , , and are evaluated at a staggered grid (standard Arakawa Cgrid) as shown in Figure 1; then the shallowwater equations can be solved by using finite difference method.
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
in which , , , , , , , , , , and .
3.2. The LeapFrog Scheme
The Leapfrog 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
The Leapfrog method allows the direct calculation of , , from the known values at time levels and , 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 , and after some rearrangement, (5) can be written as
in which ,
In a similar fashion, (6) can be written as the following matrix equations:
where ,
Definition 1. The twodimensional discrete Fourier transform of is the function defined by (see [39])
for .
We begin by taking the discrete Fourier transform of both sides of (7), and we can obtain
By making the change of variables we get
Similarly, we have
Thus using the expressions (13)(14) in (12) leads to
where , and the growth matrix
in which
Remark 2. The discrete Fourier transform is used to deal with the vector . 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 threedimensional discrete Fourier transform of is the function defined by
for .
Taking the discrete Fourier transform (18) to the both sides of (9), similar to the derivation of (15), we have
where the growth matrix
in which
Theorem 4. The difference schemes (7) and (9) are stable in the norm; then there exist constants , , , and , independent of , , and , so that
for , and , , 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 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
If is an eigenvalue of , by using Gerschgorin circle theorem and the triangular inequality, one has
Then the timespace steps size and model parameters yield the following conditions (one takes , where ):
Therefore, one has , and the FTCS scheme is conditionally stable. The condition is as follows:
in which , , and .
Similarly, one obtains the stable condition of the explicit Leapfrog finitedifference scheme (9):
where , , , and are the same as in expression (28).
Remark 6. The derivations of stability conclusions in this study are still valid for both Agrid and Bgrid; the results depend mainly on the choice of vector ; for example, in the Agrid, we take . In addition, it is easy to prove that the stability conditions derived from Cgrid are the same for both Agrid and Bgrid.
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 ()
This is the CourantFriedrichsLewy condition (CFL condition; see [13, 34, 42]) in twodimensional case. Assuming the terms , and in the expression (26), then we have
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 , the first and second expressions in (27) are given as
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 timestep 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 Leapfrog finitedifference 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 spacestep size . Standard parameter values in the shallowwater model are shown in Table 1.

5.1. Example 1
We take , the zonal velocity , and the meridional velocity . According to the expression (28), it is easy to obtain ; multiplying by the time scale, we have s. In the light of the CFL condition (30), we have s. With the expression (31), we obtain s. From the stability criterion (32), we get s.
Case 1. Setting a timestep size 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)
(a)
(b)
Case 2. When we choose the step size 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)
Case 3. The model is run with a time step size of 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)
(a)
(b)
(c)
It is obvious that the stability condition (28) is reasonable, because the numerical model is unstable when we take the time step size 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.5layer shallowwater 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 timespace 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 timedependent solution that evolves in one day. We will make an attempt to use this explicit shallowwater numerical model to simulate the Kuroshio current and its extension system in further studies.
6. Conclusions
The FTCS and the Leapfrog finite difference scheme for solving 1.5layer shallowwater 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 atmosphereocean 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 (111495jch).
References
 C. B. Vreugdenhil, Numerical Methods for Shallow Water Flow, Kluwer Academic Publishers, Boston, Mass, USA, 1994. View at: Publisher Site
 S.I. Iga and Y. Matsuda, “Shear instability in a shallow water model with implications for the Venus atmosphere,” Journal of the Atmospheric Sciences, vol. 62, no. 7, part 2, pp. 2514–2527, 2005. View at: Publisher Site  Google Scholar  MathSciNet
 J. Lambaerts, G. Lapeyre, and V. Zeitlin, “Moist versus dry barotropic instability in a shallowwater model of the atmosphere with moist convection,” Journal of the Atmospheric Sciences, vol. 68, no. 6, pp. 1234–1252, 2011. View at: Publisher Site  Google Scholar
 S. Jiang, F. F. Jin, and M. Ghil, “Multiple equilibria, periodic and aperiodic solutions in a winddriven doublegyre, shallowwater model,” Journal of Physical Oceanography, vol. 25, pp. 764–786, 1995. View at: Google Scholar
 F. Primeau, “Multiple equilibria and lowfrequency variability of the winddriven ocean circulation,” Journal of Physical Oceanography, vol. 32, no. 8, pp. 2236–2256, 2002. View at: Publisher Site  Google Scholar  MathSciNet
 E. Simonnet, M. Ghil, K. Ide, R. Temam, and S. Wang, “Lowfrequency variability in shallowwater models of the winddriven ocean circulation. I. Steadystate solution,” Journal of Physical Oceanography, vol. 33, no. 4, pp. 712–728, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 F. W. Primeau and D. Newman, “Bifurcation structure of a winddriven shallow water model with layeroutcropping,” Ocean Modelling, vol. 16, no. 34, pp. 250–263, 2007. View at: Publisher Site  Google Scholar
 O. B. Andersen, “Shallow water tides in the northwest European shelf region from TOPEX/POSEIDON altimetry,” Journal of Geophysical Research C: Oceans, vol. 104, no. 4, pp. 7729–7741, 1999. View at: Publisher Site  Google Scholar
 S. B. Yoon, C. H. Lim, and J. Choi, “Dispersioncorrection finite difference model for simulation of transoceanic tsunamis,” Terrestrial, Atmospheric and Oceanic Sciences, vol. 18, no. 1, pp. 31–53, 2007. View at: Google Scholar
 N. V. Chemetov, F. Cipriano, and S. Gavrilyuk, “Shallow water model for lakes with friction and penetration,” Mathematical Methods in the Applied Sciences, vol. 33, no. 6, pp. 687–703, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 J. Carlson, A. Jaffe, and A. Wiles, The Millennium Prize Problems, American Mathematical Society, 2006.
 S. Chippada, C. N. Dawson, M. L. Martinez, and M. F. Wheeler, “A Godunovtype finite volume method for the system of shallow water equations,” Computer Methods in Applied Mechanics and Engineering, vol. 151, no. 12, pp. 105–129, 1998. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 F. Alcrudo and P. GarciaNavarro, “Highresolution Godunovtype scheme in finite volumes for the 2D shallowwater equations,” International Journal for Numerical Methods in Fluids, vol. 16, no. 6, pp. 489–505, 1993. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 R. A. Walters, “Numerically induced oscillations in finite element approximations to shallow water equations,” International Journal for Numerical Methods in Fluids, vol. 3, no. 6, pp. 591–604, 1983. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 V. Casulli, “Semiimplicit finite difference methods for the twodimensional shallow water equations,” Journal of Computational Physics, vol. 86, no. 1, pp. 56–74, 1990. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 G.F. Lin, J.S. Lai, and W.D. Guo, “Finitevolume componentwise TVD schemes for 2D shallow water equations,” Advances in Water Resources, vol. 26, no. 8, pp. 861–873, 2003. View at: Publisher Site  Google Scholar
 M. L. Batteen and Y. J. Han, “On the computational noise of finitedifference schemes used in ocean models,” Svenska Geofysiska Föreningen, vol. 33, no. 4, pp. 387–396, 1981. View at: Google Scholar  MathSciNet
 V. Casullia and R. A. Waltersb, “An unstructured grid, threedimensional model based on the shallow water equations,” International Journal for Numerical Methods in Fluids, vol. 32, pp. 331–348, 2000. View at: Publisher Site  Google Scholar
 W. C. Thacker, “Irregular grid finitedifference techniques: Simulations of oscillations in shallow circular basins,” Journal of Physical Oceanography, vol. 7, no. 2, pp. 284–292, 1977. View at: Publisher Site  Google Scholar
 V. Casulli and R. T. Cheng, “Semiimplicit finite difference methods for threedimensional shallow water flow,” International Journal for Numerical Methods in Fluids, vol. 15, no. 6, pp. 629–648, 1992. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 Y. Xing and C.W. Shu, “High order finite difference WENO schemes with the exact conservation property for the shallow water equations,” Journal of Computational Physics, vol. 208, no. 1, pp. 206–227, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. D. Nair, S. J. Thomas, and R. D. Loft, “A discontinuous Galerkin global shallow water model,” Monthly Weather Review, vol. 133, no. 4, pp. 876–888, 2005. View at: Publisher Site  Google Scholar
 G. Quattrocchi, S. Pierini, and H. A. Dijkstra, “Intrinsic lowfrequency variability of the Gulf Stream,” Nonlinear Processes in Geophysics, vol. 19, no. 2, pp. 155–164, 2012. View at: Publisher Site  Google Scholar
 M. J. Schmeits and H. A. Dijkstra, “Bimodal behavior of the Kuroshio and the Gulf Stream,” Journal of Physical Oceanography, vol. 31, no. 12, pp. 3435–3456, 2001. View at: Publisher Site  Google Scholar
 R. J. Greatbatch, X. Zhai, M. Claus, L. Czeschel, and W. Rath, “Transport driven by eddy momentum fluxes in the Gulf Stream Extension region,” Geophysical Research Letters, vol. 37, no. 24, Article ID L24401, 2010. View at: Publisher Site  Google Scholar
 B. Qiu and W. Miao, “Kuroshio path variations south of Japan: bimodality as a selfsustained internal oscillation,” Journal of Physical Oceanography, vol. 30, no. 8, pp. 2124–2137, 2000. View at: Publisher Site  Google Scholar
 Q. Wang, M. Mu, and H. A. Dijkstra, “Application of the conditional nonlinear optimal perturbation method to the predictability study of the Kuroshio large meander,” Advances in Atmospheric Sciences, vol. 29, no. 1, pp. 118–134, 2012. View at: Publisher Site  Google Scholar
 S. Pierini, “A Kuroshio extension system model study: decadal chaotic selfSustained oscillations,” Journal of Physical Oceanography, vol. 36, no. 8, pp. 1605–1625, 2006. View at: Publisher Site  Google Scholar
 F. Primeau and D. Newman, “Elongation and contraction of the western boundary current extension in a shallowwater model: a bifurcation analysis,” Journal of Physical Oceanography, vol. 38, no. 7, pp. 1469–1485, 2008. View at: Publisher Site  Google Scholar
 W. Kramer, H. A. Dijkstra, S. Pierini, and P. J. Van Leeuwen, “Measuring the impact of observations on the predictability of the kuroshio extension in a shallowwater model,” Journal of Physical Oceanography, vol. 42, no. 1, pp. 3–17, 2012. View at: Publisher Site  Google Scholar
 G. M. Almeida, P. Bates, J. E. Freer, and M. Souvignet, “Improving the stability of a simple formulation of the shallow water equations for 2D flood modeling,” Water Resources Research, vol. 48, p. W05528, 2012. View at: Publisher Site  Google Scholar
 S. Sankaranarayanan and M. L. Spaulding, “Dispersion and stability analyses of the linearized twodimensional shallow water equations in boundaryfitted coordinates,” International Journal for Numerical Methods in Fluids, vol. 42, no. 7, pp. 741–763, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Israeli, N. H. Naik, and M. A. Cane, “An unconditionally stable scheme for the shallow water equations,” Monthly Weather Review, vol. 128, no. 3, pp. 810–823, 2000. View at: Publisher Site  Google Scholar
 R. Heihes and D. A. Randall, “Numerical integration of the shallowwater equations on a twisted icosahedral grid—partI: basic design and results of tests,” Monthly Weather Review, vol. 123, pp. 1862–1880, 1995. View at: Publisher Site  Google Scholar
 A. Grammeltvedt, “A survey of finitedifference schemes for the primitive equations for a barotropic fluid,” Monthly Weather Review, vol. 97, pp. 384–404, 1969. View at: Google Scholar
 T. C. Rebollo, A. D. Delgado, and E. D. Nieto, “A family of stable numerical solvers for the shallow water equations with source terms,” Computer Methods in Applied Mechanics and Engineering, vol. 192, no. 12, pp. 203–225, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 J. D. Anderson, Computational Fluid Dynamics: The Basics with Application, McGrawHill Science, 1995.
 H. A. Dijkstra, Nonlinear Physical Oceanography: A Dynamical Systems Approach to Thelarge Scale Ocean Circulation and El Nino, Springer, Dordrecht, The Netherlands, 2nd edition, 2005.
 J. W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, vol. 22 of Texts in Applied Mathematics, Springer, New York, NY, USA, 1998. View at: MathSciNet
 J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, The Wadsworth & Brooks/Cole Mathematics Series, Wadsworth & Brooks, Pacific Grove, Calif, USA, 1989. View at: MathSciNet
 B. Engquist and A. Majda, “Radiation boundary conditions for acoustic and elastic wave calculations,” Communications on Pure and Applied Mathematics, vol. 32, no. 3, pp. 314–358, 1979. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 D. R. Durran, Numerical Methods for Fluid Dynamics, vol. 32 of Texts in Applied Mathematics, Springer, New York, NY, USA, 2nd edition, 2010. View at: Publisher Site  MathSciNet
 A. F. Blumberg and G. L. Mellor, “A description of a threedimensional coastal ocean model,” in ThreeDimensional Coastal Ocean Models, Coastal and Estuarine Sciences, vol. 32, pp. 1–16, 1987. View at: Google Scholar
Copyright
Copyright © 2013 Guangan Zou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.