The explicit finite difference scheme for solving an intermediate coupled ocean-atmosphere 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 Courant-Friedrichs-Lewy (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 ocean-atmosphere equations. Research on numerical simulation of the ocean-atmosphere 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 [123]. 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 ocean-atmosphere 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 ocean-atmosphere 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 ocean-atmosphere model used here includes the model of ocean fluid dynamics, the mixed-layer thermodynamics, and the empirical atmospheric model. This coupled ocean-atmosphere 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.5-layer 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ño-Southern Oscillation (ENSO), to simulate the seasonal cycle of sea surface temperature in the tropical Pacific Ocean [2428]. The main purpose of present study is to investigate the stability analysis of numerical model from the mathematical analysis point of view.

The upper-ocean flow is governed by reduced-gravity shallow-water 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 mixed-layer, 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 upper-ocean flow velocity, surface Ekman flow is determined by Coriolis force and wind-stress 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 SST-forced 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 two-dimensional 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 time-space steps and model parameters yield the following conditions:

Then we have , so the explicit finite difference scheme (13) is conditionally stable. Consequently, the time-step size satisfies the following inequality: In which , , , .

When the rotation, dissipation, and friction terms are neglected, we have

This is Courant-Friedrichs-Lewy (CFL) condition [3234]; 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 ocean-atmosphere 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. Free-slip boundary conditions in velocity and no-flux 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 time-step (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.

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 time-step 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 ocean-atmosphere 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 ocean-atmosphere 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 ocean-atmosphere system more clearly. This intermediate coupled ocean-atmosphere 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.


The authors thank the editor and the reviewers for their helpful and constructive suggestions.