Research Article  Open Access
On the Stability of an Intermediate Coupled OceanAtmosphere Model
Abstract
The explicit finite difference scheme for solving an intermediate coupled oceanatmosphere equations has been proposed and discussed. The discrete Fourier analysis within Gerschgorin circle theorem is applied to the stability analysis of this numerical model. The stability criterion that we obtained includes advection, rotation, dissipation, and friction terms, without any assumptions, which is also including the CourantFriedrichsLewy (CFL) condition as a special case. Numerical sensitivity experiments are also carried out by varying the model parameters.
1. Introduction
The geophysical motions in both the atmosphere and ocean can be described by partial differential equations (PDEs), in recent years, which have become the very popular and important tools in the study of climate change, weather forecasting, and climate prediction. However, both the coupling process of the atmosphere to the ocean and the corresponding PDEs are very complicated; it is almost impossible to find the exact solution of coupled oceanatmosphere equations. Research on numerical simulation of the oceanatmosphere system has aroused many scientists and engineers’ interest; a great variety of numerical methods (especially for finite difference method) have been developed to solve this PDEs system [1–23]. Nevertheless, as we know, there are few people who gives the stability analysis of these complicated numerical models from the mathematical point of view.
The main purpose of this study is to introduce and solve an intermediate coupled oceanatmosphere PDEs. The stability analysis of numerical method has been taken into consideration by the discrete Fourier analysis combined with Gerschgorin circle theorem. Compared with the time step allowed by the CFL stability criterion, our stability bounds are more accurate and effective. Numerical examples are also presented to test the sensitivity of model. This paper is organized as follows. In the next section, the brief description of an intermediate coupled oceanatmosphere mathematical model has been introduced. In Section 3, the explicit finite difference scheme is used to solve the PDEs. The stability conditions are given by the numerical analysis in Section 4. Then the results of experiments are presented and discussed in Section 5. Finally, conclusions are drawn in Section 6.
2. The Mathematical Model
The intermediate coupled oceanatmosphere model used here includes the model of ocean fluid dynamics, the mixedlayer thermodynamics, and the empirical atmospheric model. This coupled oceanatmosphere model is a modified version of the intermediate coupled model (ICM) developed by Chang [24] and Wang et al. [25], which is an extension of 1.5layer reduced gravity system that includes the physics of the surface mixed layer and allows the prediction of the sea surface temperature. ICM had been successfully used to the study of El NiñoSouthern Oscillation (ENSO), to simulate the seasonal cycle of sea surface temperature in the tropical Pacific Ocean [24–28]. The main purpose of present study is to investigate the stability analysis of numerical model from the mathematical analysis point of view.
The upperocean flow is governed by reducedgravity shallowwater equations: in which are the horizontal velocities; is the Coriolis parameter of the equatorial betaplane approximation, where , and is radius of the earth; A density difference between the upper and lower layers results in a reduced gravity constant (where is the acceleration of gravity, and are the constant upper and lower layer densities, respectively, and is a mean density); , are the wind stress over the sea surface; is the upper layer thickness, where is the undisturbed layer thickness, is the interface displacement; is the lateral eddy viscosity coefficient; is the interfacial friction coefficient.
Assuming that the effects of compressibility and salinity are of secondary importance in the mixedlayer, the thermodynamics equation can be expressed by where is sea surface temperature (SST), , are the horizontal velocities in the mixed layer, is average SST, is horizontal heat diffusion coefficient, is adjustable coefficient, is downward heat flux, is mean depth of mixed layer, is the entrainment velocity, is a Heaviside step function of (, if , , if ) and is entrainment water temperature.
In (4), the horizontal velocity in the surface mixed layer is separated into two components: , where is upperocean flow velocity, surface Ekman flow is determined by Coriolis force and windstress force where is the Rayleigh friction coefficient.
The entrainment velocity is determined as divergence of surface flow and the temperature of entrained water is taked as a linear form: where is the observed mean temperature, is the vertical gradient of at m depth, and is the fluctuation of thermocline.
The atmospheric part of the model is written as the sum of the annual mean wind field (or climatology of the monthly mean when the seasonal cycle is involved), coupled feedback, and atmosphere internal variability (see [29]): where and are usually prescribed from observations; and represent coupling between the atmosphere and ocean, which are empirical functions of SST anomaly; and are the coupling strength. The SSTforced surface wind stress and heat flux are determined by and , which can be written as linear integral operators over the entire domain , that is, where and are given by a simple dynamical model of atmosphere.
3. The Explicit Finite Difference Scheme
The domain is discretized with a spacing of in the direction, in the direction, and in the direction, and we define with , , for , . The governing equations will be solved on a staggered (Arakawa C) grid. That is, the , are evaluated at the cell boundaries and the and points are in the center grids.
The forward difference approximation is used for the time derivative, and the central difference approximation for the spatial derivatives. The finite difference operators are defined as
Therefore, the difference approximation of the partial differential equations (1)–(4) is given by in which , , , , , , , .
Equations (12) are finite difference equations which represent the original partial differential equations expressed in (1)–(4).
4. Stability Criterion
In this section, the stability analysis of numerical schemes will be investigated by the discrete Fourier analysis within Gerschgorin circle theorem. A series of necessary conditions are also obtained. The details are given as follows.
Assuming that , after some rearrangement, (12) can be written as in which ,
Definition 1. The twodimensional discrete Fourier transform of is the function defined by (see[30]): for .
By using the discrete Fourier analysis to transform both sides of (9), we can obtain where , and the growth matrix in which
Theorem 2. Let be the spectral radius of the matrix . If there exist constants , , , and C, independent of , , , and , , such that for , , and , and for all , then the scheme (18) is stable in the norm.
Theorem 3 (Gerschgorin Circle Theorem, see [31]). Let be a complex matrix, are elements of matrix . For each eigenvalue of , there exists an such that where are the sum of the absolute values of the elements in the th row except for the diagonal element.
By Gerschgorin Circle Theorem, if is an eigenvalue of , we have
Using the backwards triangular inequality, we take , if the timespace steps and model parameters yield the following conditions:
Then we have , so the explicit finite difference scheme (13) is conditionally stable. Consequently, the timestep size satisfies the following inequality: In which , , , .
When the rotation, dissipation, and friction terms are neglected, we have
This is CourantFriedrichsLewy (CFL) condition [32–34]; obviously, it is only a special case in our results.
5. Numerical Results
Sensitivity studies are a necessary and important part of the developments of mathematical models of geophysical fluid systems. They can help to reveal aspects of the model that will most profitably benefit from further refinement, as well as providing insights into the fundamental dynamics of these complicated fluid systems [18]. In this study, we consider the sensitivity of an intermediate coupled oceanatmosphere model to change in the Coriolis parameters , the interfacial friction coefficient , the lateral eddy viscosity coefficient , and the adjustable coefficient . Parameter values are shown in Table 1.

The model domain extends form S to N in latitude and from E to W in longitude. The model resolution is in direction and in direction. Freeslip boundary conditions in velocity and noflux boundary conditions in temperature are used in the ocean model.
5.1. Example 1
According to the stability criterion (25), we obtain ; from the CFL condition, we have .
Running the model for one year, Figure 1 gives the results of zonal velocity and meridional velocity (at the latitude N) by changing with timestep (dt).
When we take , Figure 2 shows the results of zonal velocity at the latitude N and the longitude E, respectively, and the maximum speed is 3.5 m/s, which is not in accordance with the actual condition. When we take , the results overflow. The range of the CFL condition is too big; that is to say, the CFL condition is less strict and accurate than our results.
(a)
(b)
5.2. Example 2
In this numerical tests, we will discuss the sensitivity of model parameters. The model has been run for one year, and we take the timestep size s. The variable changes that result from the variation of coefficients are only given one (as well as the others) in each experiments on the model simulation.
Figure 3 gives the results of sea surface temperature (SST) with the variation of Coriolis parameters (the variation of the angle (theta)) at the latitude N and the longitude E, respectively. Figure 4 displays the sea surface height (SSH) with the variation of the interfacial friction coefficient (gamma). Figure 5 shows the meridional velocity with the variation of lateral eddy viscosity coefficient . Figure 6 exhibits the meridional velocity (cm/s) with the variation of adjustable coefficient (lambda).
We find that Coriolis parameters, lateral eddy viscosity coefficient, interfacial friction coefficient, and adjustable coefficient have an impact on the model evolutions. Therefore, taking the CFL condition as the stability criterion alone is not reasonable. On the other hand, it is also important to derive future model developments and test the numerical model.
6. Conclusion
The present study provides an intermediate coupled oceanatmosphere model, and the stability criterion are also obtained for the choice of time and space steps size, in order to ensure the errors made at one time step of the calculation do not cause the errors to increase as the computations are continued. Our results are effective and reasonable, which are including CFL condition. Numerical sensitivity experiments show that the Coriolis parameters, lateral eddy viscosity coefficient, interfacial friction coefficient, and adjustable coefficient are important to the stability of model. However, external forcing also plays a dominant role in oceanatmosphere system, for example, wind stress curl, heat flux, and so forth [18, 35, 36]; are not reported in this study. Based on the parameters sensitivity studies, the present study allows understanding the physical mechanism of the oceanatmosphere system more clearly. This intermediate coupled oceanatmosphere model will be used for further study, and the method of stability analysis in present work can also be used to the stability study of earth system numerical model.
Acknowledgments
The authors thank the editor and the reviewers for their helpful and constructive suggestions.
References
 K. Bryan, “A numerical method for the study of the circulation of the world ocean,” Journal of Computational Physics, vol. 4, no. 3, pp. 347–376, 1969. View at: Google Scholar
 R. N. Miller, Numerical Modeling of Ocean Circulation, Cambridge University Press, New York, NY, USA, 2007.
 C. Chen, H. Liu, and R. C. Beardsley, “An unstructured grid, finitevolume, threedimensional, primitive equations ocean model: application to coastal ocean and estuaries,” Journal of Atmospheric and Oceanic Technology, vol. 20, no. 1, pp. 159–186, 2003. View at: Google Scholar
 A. F. Blumberg and G. L. Mellor, “A description of a threedimensional coastal ocean circulation model,” in ThreeDimensional Coastal Ocean Models, N. S. Heaps, Ed., vol. 4 of Coastal and Estuarine Series, pp. 1–16, American Geophysical Union, 1987. View at: Google Scholar
 S. Danilov, “Two finitevolume unstructured mesh models for largescale ocean modeling,” Ocean Modelling, vol. 47, pp. 14–25, 2012. View at: Publisher Site  Google Scholar
 A. M. Davies, J. Xing, and J. E. Jones, “A model study of tidal distributions in the Celtic and Irish Sea regions determined with finite volume and finite element models,” Ocean Dynamics, vol. 61, no. 10, pp. 1645–1667, 2011. View at: Publisher Site  Google Scholar
 J. Eric Jones and A. M. Davies, “Application of a finite element model to the computation of tides in the Mersey Estuary and Eastern Irish Sea,” Continental Shelf Research, vol. 30, no. 5, pp. 491–514, 2010. View at: Publisher Site  Google Scholar
 Z. Q. Yang and T. R. Khangaonkar, “Modeling tidal circulation and stratification in Skagit River estuary using an unstructured grid ocean model,” Ocean Modelling, vol. 28, no. 1–3, pp. 34–49, 2009. View at: Publisher Site  Google Scholar
 Y. L. Zhang and A. M. Baptista, “SELFE: a semiimplicit EulerianLagrangian finiteelement model for crossscale ocean circulation,” Ocean Modelling, vol. 21, no. 34, pp. 71–96, 2008. View at: Publisher Site  Google Scholar
 Q. Wang, S. Danilov, and J. Schröter, “Finite element ocean circulation model based on triangular prismatic elements, with application in studying the effect of topography representation,” Journal of Geophysical Research, vol. 113, no. C5, 2008. View at: Publisher Site  Google Scholar
 S. Danilov, G. Kivman, and J. Schröter, “A finiteelement ocean model: principles and evaluation,” Ocean Modelling, vol. 6, no. 2, pp. 125–150, 2004. View at: Publisher Site  Google Scholar
 R. Bleck, “An oceanic general circulation model framed in hybrid isopycnicCartesian coordinates,” Ocean Modelling, vol. 4, no. 1, pp. 55–88, 2002. View at: Publisher Site  Google Scholar
 G. J. Fix, “Finite element models for ocean circulation problems,” SIAM Journal on Applied Mathematics, vol. 29, no. 3, pp. 371–387, 1975. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Gerdes, “A primitive equation ocean circulation model using a general vertical coordinate transformation 1. Description and testing of the model,” Journal of Geophysical Research, vol. 98, no. 8, pp. 14683–14701, 1993. View at: Google Scholar
 S. M. Griffies, C. Böning, F. O. Bryan et al., “Developments in ocean climate modelling,” Ocean Modelling, vol. 2, no. 34, pp. 123–192, 2000. View at: Google Scholar
 D. Haidvogel, E. Curchitser, M. Iskandarani, R. Hughes, and M. Taylor, “Global modeling of the ocean and atmosphere using the spectral element method,” AtmosphereOcean, vol. 35, pp. 505–531, 1997. View at: Google Scholar
 R. A. Walters, “Design considerations for a finite element coastal ocean model,” Ocean Modelling, vol. 15, no. 12, pp. 90–100, 2006. View at: Publisher Site  Google Scholar
 F. Bryan, “Parameter sensitivity of primitive equation ocean general circulation models,” Journal of Physical Oceanography, vol. 17, no. 7, pp. 970–985, 1987. View at: Google Scholar
 A. Rosati and K. Miyakoda, “A general circulation model for upper ocean simulation,” Journal of Physical Oceanography, vol. 11, pp. 1601–1626, 1988. View at: Google Scholar
 K. Bryan, S. Manabe, and R. C. Pacanowski, “A global oceanatmosphere climate model. Part II. The oceanic circulation,” Journal of Physical Oceanography, vol. 5, pp. 30–46, 1975. View at: Google Scholar
 Y.J. Han, “A numerical world ocean general circulation model. Part II. A baroclinic experiment,” Dynamics of Atmospheres and Oceans, vol. 8, no. 2, pp. 141–172, 1984. View at: Google Scholar
 M. J. Roberts and R. A. Wood, “Topographic sensitivity studies with a BryanCoxtype ocean model,” Journal of Physical Oceanography, vol. 27, no. 5, pp. 823–836, 1997. View at: Google Scholar
 J. A. T. Bye, “Coupling oceanatmosphere models,” EarthScience Reviews, vol. 40, no. 12, pp. 149–162, 1996. View at: Google Scholar
 P. Chang, “A study of the seasonal cycle of sea surface temperature in the tropical Pacific Ocean using reduced gravity models,” Journal of Geophysical Research, vol. 99, no. 4, pp. 7725–7741, 1994. View at: Publisher Site  Google Scholar
 B. Wang, T. Li, and P. Chang, “An intermediate model of the tropical Pacific Ocean,” Journal of Physical Oceanography, vol. 25, pp. 1599–1616, 1995. View at: Google Scholar
 P. Chang, B. Wang, T. Li, and L. Ji, “Interactions between the seasonal cycle and the Southern Oscillationfrequency entrainment and chaos in a coupled oceanatmosphere model,” Geophysical Research Letters, vol. 21, no. 25, pp. 2817–2820, 1994. View at: Google Scholar
 P. Chang, L. Ji, B. Wang, and T. Li, “Interactions between the seasonal cycle and El NinoSouthern Oscillation in an intermediate coupled oceanatmosphere model,” Journal of the Atmospheric Sciences, vol. 52, no. 13, pp. 2353–2372, 1995. View at: Google Scholar
 P. Chang, L. Ji, H. Li, and M. Flügel, “Chaotic dynamics versus stochastic processes in El NiñoSouthern Oscillation in coupled oceanatmosphere models,” Physica D, vol. 98, no. 2–4, pp. 301–320, 1996. View at: Google Scholar
 F. M. Wang and P. Chang, “A linear stability analysis of couple tropical Atlantic variability,” Journal of Climate, vol. 21, no. 11, pp. 2421–2436, 2008. View at: Publisher Site  Google Scholar
 J. W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, Springer, New York, NY, USA, 1997.
 G. H. Golub and C. F. van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, Md, USA, 3rd edition, 1996. View at: MathSciNet
 J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Brooks Cole, Pacific Grove, Calif, USA, 1989. View at: MathSciNet
 D. R. Durran, Numerical Methods for Fluid Dynamics: With Applications to Geophysics, vol. 32, Springer, New York, NY, USA, 2nd edition, 2010. View at: Publisher Site  MathSciNet
 R. Heihes and D. A. Randall, “Numerical integration of the shallowwater equations on a twisted icosahedral grid. Part I: basic design and results of tests,” Monthly Weather Review, vol. 123, no. 6, pp. 1862–1880, 1995. View at: Google Scholar
 W. G. Large, G. Danabasoglu, S. C. Doney, and J. C. Mcwilliams, “Sensitivity to surface forcing and boundary layer mixing in a global ocean model: annualmean climatology,” Journal of Physical Oceanography, vol. 27, no. 11, pp. 2418–2447, 1997. View at: Google Scholar
 D. G. Wright and T. F. Stocker, “Sensitivities of a zonally averaged global ocean circulation model,” Journal of Geophysical Research, vol. 97, no. 8, pp. 12707–12730, 1992. View at: Google Scholar
Copyright
Copyright © 2013 Tianxu Zhao and Guangan Zou. 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.