On the Stability of an Intermediate Coupled Ocean-Atmosphere Model

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.


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

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 [24][25][26][27][28].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 u = (, V) is governed by reducedgravity shallow-water equations: in which , V are the horizontal velocities;  =  0 +  is the Coriolis parameter of the equatorial betaplane approximation, where  0 = 2Ω sin ,  = (2Ω cos )/ 0 and  0 is radius of the earth; A density difference between the upper and lower layers results in a reduced gravity constant   = ( 2 −  1 )/ 0 (where  is the acceleration of gravity,  1 and  2 are the constant upper and lower layer densities, respectively, and  0 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),   , V  are the horizontal velocities in the mixed layer,  0 is average SST, ] is horizontal heat diffusion coefficient,  is adjustable coefficient,  is downward heat flux,  0 is mean depth of mixed layer,   is the entrainment velocity, (  ) is a Heaviside step function of   (() = 1, if  > 0, () = 0, if  < 0) and   is entrainment water temperature.In (4), the horizontal velocity u  = (  , V  ) in the surface mixed layer is separated into two components: u  = u + u  ( −  0 )/, where u = (, V) is upper-ocean flow velocity, surface Ekman flow u  = (  , V  ) 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   =  0 ∇⋅u  and the temperature of entrained water   is taked as a linear form: where  50 is the observed mean temperature,  50 / is the vertical gradient of  at 50 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; A and B 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 A and B, 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.

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   =  0 + Δ,   =  0 + Δ,   = Δ for ,  = 0, 1, . ..,  = 0, 1, . ... The governing equations will be solved on a staggered (Arakawa C-) grid.That is, the , V 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).

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.
By using the discrete Fourier analysis to transform both sides of (9), we can obtain where Then we have || ≤ 1+Δ, so the explicit finite difference scheme ( 13) is conditionally stable.Consequently, the timestep size satisfies the following inequality: In which  = (1/Δ)(max When the rotation, dissipation, and friction terms are neglected, we have This is Courant-Friedrichs-Lewy (CFL) condition [32][33][34]; obviously, it is only a special case in our results.

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 30 ∘ S to 30 ∘ N in latitude and from 120 ∘ E to 80 ∘ W in longitude.The model resolution is 2 ∘ in  direction and 1 ∘ in  direction.Free-slip boundary conditions in velocity and no-flux boundary conditions in temperature are used in the ocean model.
Running the model for one year, Figure 1 gives the results of zonal velocity  and meridional velocity V (at the latitude 0 ∘ N) by changing with time-step Δ(dt).
When we take Δ = 14306 s, Figure 2 shows the results of zonal velocity  at the latitude 0 ∘ N and the longitude 180 ∘ E, respectively, and the maximum speed is 3.5 m/s, which is not in accordance with the actual condition.When we take Δ > 14306 s, 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.

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 Δ = 10800 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 0 ∘ N and the longitude 180 ∘ 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 V (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.

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

Table 1 :
The values of parameters in the model.