On the Stability of an Intermediate Coupled Ocean-Atmosphere Model

Abstract

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 [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 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 [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 upper-ocean flow 𝐮=(𝑢,𝑣) is governed by reduced-gravity shallow-water equations: 𝜕𝑢𝜕𝑡+𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝑦−𝑓𝑣=−𝑔^′𝜕ℎ+𝜏𝜕𝑥𝜌₀𝐻^𝑥𝜕+𝐾²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²−𝛾𝑢,𝜕𝑣𝜕𝑡+𝑢𝜕𝑣𝜕𝑥+𝑣𝜕𝑣𝜕𝑦+𝑓𝑢=−𝑔^′𝜕ℎ+𝜏𝜕𝑦𝜌₀𝐻^𝑦𝜕+𝐾²𝑣𝜕𝑥²+𝜕²𝑣𝜕𝑦²−𝛾𝑣,𝜕ℎ+𝜕𝑡𝜕𝐻𝑢+𝜕𝑥𝜕𝐻𝑣𝜕𝑦=0, (1)(2)(3) in which 𝑢,𝑣 are the horizontal velocities; 𝑓=𝑓₀+𝛽𝑦 is the Coriolis parameter of the equatorial betaplane approximation, where 𝑓₀=2Ωsin𝜃, 𝛽=(2Ωcos𝜃)/𝑟₀ 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 𝜕𝑇𝜕𝑡+𝑢_𝑠𝜕𝑇𝜕𝑥+𝑣_𝑠𝜕𝑇𝜕𝜕𝑦=𝜈²𝑇𝜕𝑥²+𝜕²𝑇𝜕𝑦²−𝜆𝑇−𝑇₀+𝑄𝜌₀𝐶_𝑝𝐻₀−1𝐻₀𝑤_𝑒𝐻𝑤_𝑒𝑇−𝑇_𝑒, (4) 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 𝑤_𝑒 (𝐻(𝑥)=1, if 𝑥>0, 𝐻(𝑥)=0, if 𝑥<0) 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 𝑟_𝑠𝑢_𝑒−𝑓𝑣_𝑒=𝜏_𝑥𝜌₀𝐻₀,𝑟_𝑠𝑣_𝑒−𝑓𝑢_𝑒=𝜏_𝑦𝜌₀𝐻₀, (5) 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: 𝑇_𝑒=0.5𝑇+𝑇₅₀+2𝜕𝑇₅₀ℎ𝜕𝑧^′, (6) where 𝑇₅₀ is the observed mean temperature, 𝜕𝑇₅₀/𝜕𝑧 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]): 𝜏𝑄=𝜏𝑄+𝜇𝒜𝑇^′𝜔ℬ𝑇^′+𝜁𝜒, (7) 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, 𝒜𝑇^′=𝜓,𝑇^′=_Ω𝜓_{(𝑥,𝑦;𝑥′,𝑦′)}𝑇^′_{𝑥′,𝑦′,𝑡}𝑑𝑥^′𝑑𝑦^′,ℬ𝑇^′=𝜙,𝑇^′=_Ω𝜙_{(𝑥,𝑦;𝑥′,𝑦′)}𝑇^′_{𝑥′,𝑦′,𝑡}𝑑𝑥^′𝑑𝑦^′, (8) 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 𝑖,𝑗=0,1,…, 𝑛=0,1,…. 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 𝛿_𝑡𝛼_𝑖,𝑗=𝛼^𝑛+1_𝑖,𝑗−𝛼^𝑛_𝑖,𝑗,𝛿Δ𝑡_𝑥𝛼_𝑖,𝑗=𝛼^𝑛_𝑖+1,𝑗−𝛼^𝑛_{𝑖−1,𝑗},𝛿2Δ𝑥^𝑦𝛼_𝑖,𝑗=𝛼^𝑛_𝑖,𝑗+1−𝛼^𝑛_{𝑖,𝑗−1},𝛿2Δ𝑦_2𝑥𝛼_𝑖,𝑗=𝛼^𝑛_𝑖+1,𝑗−2𝛼^𝑛_𝑖,𝑗+𝛼^𝑛_{𝑖−1,𝑗}Δ𝑥²,𝛿²_𝑦𝛼_𝑖,𝑗=𝛼^𝑛_𝑖,𝑗+1−2𝛼^𝑛_𝑖,𝑗+𝛼^𝑛_{𝑖,𝑗−1}Δ𝑦². (9)(10)(11)

Therefore, the difference approximation of the partial differential equations (1)–(4) is given by 𝛿_𝑡𝑢_{𝑖+1/2,𝑗}+𝑢^∗𝛿_𝑥𝑢_{𝑖+1/2,𝑗}+𝑣^∗𝛿_𝑦𝑢_{𝑖+1/2,𝑗}−𝑓^∗𝑣^𝑛_{𝑖,𝑗+1/2}=−𝑔^′𝛿_𝑥𝜂_𝑖,𝑗+𝜏_𝑥𝜌₀𝐻^∗+𝐾𝛿²_𝑥𝑢_{𝑖+1/2,𝑗}+𝛿²_𝑦𝑢_{𝑖+1/2,𝑗}−𝛾𝑢^𝑛_{𝑖+1/2,𝑗},𝛿_𝑡𝑣_{𝑖,𝑗+1/2}+𝑢^∗𝛿_𝑥𝑣_{𝑖,𝑗+1/2}+𝑣^∗𝛿_𝑦𝑣_{𝑖,𝑗+1/2}+𝑓^∗𝑢^𝑛_{𝑖+1/2,𝑗}=−𝑔^′𝛿_𝑦𝜂_𝑖,𝑗+𝜏_𝑦𝜌₀𝐻^∗+𝐾𝛿²_𝑥𝑣_{𝑖,𝑗+1/2}+𝛿²_𝑦𝑣_{𝑖,𝑗+1/2}−𝛾𝑣^𝑛_{𝑖,𝑗+1/2},𝛿_𝑡ℎ_𝑖,𝑗+𝑢^∗𝛿_𝑥ℎ_𝑖,𝑗+𝑣^∗𝛿_𝑦ℎ_𝑖,𝑗+𝐻^∗𝛿_𝑥𝑢_{𝑖+1/2,𝑗}+𝛿_𝑦𝑣_{𝑖,𝑗+1/2}𝛿=0,_𝑡𝑇_𝑖,𝑗+𝑢^∗_𝑠𝛿_𝑥𝑇_𝑖,𝑗+𝑣^∗_𝑠𝛿_𝑦𝑇_𝑖,𝑗=𝜈𝛿²_𝑥𝑇_𝑖,𝑗+𝛿²_𝑦𝑇_𝑖,𝑗−𝜆𝑇^𝑛_𝑖,𝑗−𝑇^∗𝛿_𝑥𝑢_{𝑖+1/2,𝑗}+𝛿_𝑦𝑣_{𝑖,𝑗+1/2}+𝐶𝑇₀,𝑄 (12) in which 𝑢^∗=(𝑢^𝑛_𝑖+1,𝑗+𝑢^𝑛_{𝑖−1,𝑗})/2, 𝑣^∗=(𝑣^𝑛_𝑖,𝑗+1+𝑣^𝑛_{𝑖,𝑗−1})/2, 𝑓^∗=𝑓₀+𝑗𝛽Δ𝑦, 𝐻^∗=𝐷+ℎ^𝑛_𝑖,𝑗, 𝑢^∗_𝑠=𝑢^∗+𝑢_𝑒(𝐻−𝐻₀)/𝐻, 𝑣^∗_𝑠=𝑣^∗+𝑣_𝑒(𝐻−𝐻₀)/𝐻, 𝑇^∗=𝐻(𝑤_𝑒)(𝑇^𝑛_𝑖,𝑗−𝑇^∗_𝑒), 𝐶(𝑇₀,𝑄)=𝜆𝑇₀+(𝑄/𝜌₀𝐶_𝑝𝐻₀).

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 𝐗^𝑛_𝑖,𝑗=(𝑢^𝑛_{𝑖+1/2,𝑗},𝑣^𝑛_{𝑖,𝑗+1/2},ℎ^𝑛_𝑖,𝑗,𝑇^𝑛_𝑖,𝑗)^𝑇, after some rearrangement, (12) can be written as 𝐗^𝑛+1_𝑖,𝑗=𝐴𝐗^𝑛_𝑖+1,𝑗+𝐶𝐗^𝑛_{𝑖−1,𝑗}+𝐵𝐗^𝑛_𝑖,𝑗+𝐷𝐗^𝑛_𝑖,𝑗+1+𝐸𝐗^𝑛_{𝑖,𝑗−1}+𝐞_𝑛 (13) in which 𝐞_𝑛=(𝜏_𝑥/𝜌𝐻^∗,𝜏_𝑦/𝜌𝐻^∗,0,𝐶(𝑇₀,𝑄))^𝑇,−𝐴=Δ𝑡𝑢^∗+2Δ𝑥Δ𝑡𝐾Δ𝑥²0−Δ𝑡𝑔^′0𝑟2Δ𝑥0−𝑢^∗2+𝑟𝐾Δ𝑥²−00Δ𝑡𝐻^∗2Δ𝑥0−Δ𝑡𝑢^∗0−2Δ𝑥Δ𝑡𝑇^∗2Δ𝑥00−Δ𝑡𝑢^∗_𝑠+2Δ𝑥Δ𝑡𝜈Δ𝑥²,𝐶=Δ𝑡𝑢^∗+2Δ𝑥Δ𝑡𝐾Δ𝑥²0Δ𝑡𝑔^′00𝑟2Δ𝑥𝑢^∗2+𝑟𝐾Δ𝑥²00Δ𝑡𝐻^∗02Δ𝑥Δ𝑡𝑢^∗02Δ𝑥Δ𝑡𝑇^∗2Δ𝑥00Δ𝑡𝑢^∗_𝑠+2Δ𝑥Δ𝑡𝜈Δ𝑥²𝑏,𝐵=¹Δ𝑡𝑓^∗00−Δ𝑡𝑓^∗𝑏₁000010000𝑏₂,𝑏₁1=1−Δ𝑡𝛾−2Δ𝑡𝐾Δ𝑥²+1Δ𝑦²,𝑏₂1=1−Δ𝑡𝛾−2Δ𝑡𝜈Δ𝑥²+1Δ𝑦²,−𝐷=Δ𝑡𝑣^∗+2Δ𝑦Δ𝑡𝐾Δ𝑦²0000−Δ𝑡𝑣^∗+2Δ𝑦Δ𝑡𝐾Δ𝑦²−Δ𝑡𝑔^′02Δ𝑦0−Δ𝑡𝐻^∗−2Δ𝑦Δ𝑡𝑣^∗02Δ𝑦0−Δ𝑡𝑇^∗2Δ𝑦0−Δ𝑡𝑣^∗_𝑠+2Δ𝑦Δ𝑡𝜈Δ𝑦²,𝐸=Δ𝑡𝑣^∗+2Δ𝑦Δ𝑡𝐾Δ𝑦²0000Δ𝑡𝑣^∗+2Δ𝑦Δ𝑡𝐾Δ𝑦²Δ𝑡𝑔^′002Δ𝑦Δ𝑡𝐻^∗2Δ𝑦Δ𝑡𝑣^∗002Δ𝑦Δ𝑡𝑇^∗02Δ𝑦Δ𝑡𝑣^∗_𝑠+2Δ𝑦Δ𝑡𝜈Δ𝑦². (14)(15)(16)

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

By using the discrete Fourier analysis to transform both sides of (9), we can obtain 𝐗^†𝑛+1(𝜉,𝜃)=𝐺(𝜉,𝜃)𝐗^†𝑛(𝜉,𝜃)+𝐞^′_𝑛, (18) where 𝐞^′_𝑛=(1/2𝜋)∑^∞_{𝑚,𝑘=−∞}exp(−𝑖(𝑚𝜉+𝑘𝜃))𝐞_𝑛, and the growth matrix 𝐺()()()()()=𝑎𝜉,𝜃=𝐴exp𝑖𝜉+𝐶exp−𝑖𝜉+𝐵+𝐷exp𝑖𝜃+𝐸exp−𝑖𝜃₁₁𝑎₁₂𝑎₁₃𝑎₁₄𝑎₂₁𝑎₂₂𝑎₂₃𝑎₂₄𝑎₃₁𝑎₃₂𝑎₃₃𝑎₃₄𝑎₄₁𝑎₄₂𝑎₄₃𝑎₄₄, (19) in which 𝑎₁₁=𝑎₂₂1=1−Δ𝑡𝛾−2Δ𝑡𝐾Δ𝑥²+1Δ𝑦²+2Δ𝑡𝐾cos𝜉Δ𝑥²+cos𝜃Δ𝑦²−𝑖Δ𝑡𝑢^∗sin𝜉+Δ𝑥𝑣^∗sin𝜃,𝑎Δ𝑦₁₂=Δ𝑡𝑓^∗,𝑎₁₃=−𝑖Δ𝑡𝑔^′sin𝜉,𝑎Δ𝑥₂₁=−Δ𝑡𝑓^∗,𝑎₂₃=−𝑖Δ𝑡𝑔^′sin𝜃,𝑎Δ𝑦₁₄=𝑎₂₄=𝑎₃₄=𝑎₄₃𝑎=0,₃₁=−𝑖Δ𝑡𝐻^∗sin𝜉Δ𝑥,𝑎₃₂=−𝑖Δ𝑡𝐻^∗sin𝜃,𝑎Δ𝑦₃₃=1−𝑖Δ𝑡𝑢^∗sin𝜉+Δ𝑥𝑣^∗sin𝜃,𝑎Δ𝑦₄₁=−𝑖Δ𝑡𝑇^∗sin𝜉Δ𝑥,𝑎₄₂=−𝑖Δ𝑡𝑇^∗sin𝜃,𝑎Δ𝑦₄₄1=1−Δ𝑡𝜆−2Δ𝑡𝜈Δ𝑥²+1Δ𝑦²+2Δ𝑡𝜈cos𝜉Δ𝑥²+cos𝜃Δ𝑦²−𝑖Δ𝑡𝑢^∗_𝑠sin𝜉+Δ𝑥𝑣^∗_𝑠sin𝜃.Δ𝑦 (20)

Theorem 2. Let 𝜌(𝐺) be the spectral radius of the matrix 𝐺. If there exist constants Δ𝑥₀, Δ𝑦₀, Δ𝑡₀, and C, independent of Δ𝑥, Δ𝑦, Δ𝑡, and 𝜉, 𝜃, such that 𝜌(𝐺)≤1+𝐶Δ𝑡 (21) for 0<Δ𝑡≤Δ𝑡₀, 0<Δ𝑥≤Δ𝑥₀, and 0<Δ𝑦≤Δ𝑦₀, 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 𝜆−𝑎_𝑠𝑠≤𝜌_𝑠,𝑠=1,2,…,𝑛, (22) where 𝜌_𝑠=∑^𝑛_{𝑗=1,𝑗≠𝑠}|𝑎_𝑠𝑗| 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 𝜆−𝑎₁₁≤𝑎₁₂+𝑎₁₃=Δ𝑡𝑓^∗+−𝑖Δ𝑡𝑔^′sin𝜉,Δ𝑥𝜆−𝑎₂₂≤𝑎₂₁+𝑎₂₃=−Δ𝑡𝑓^∗+−𝑖Δ𝑡𝑔^′sin𝜃,Δ𝑦𝜆−𝑎₃₃≤𝑎₃₁+𝑎₃₂=−𝑖Δ𝑡𝐻^∗sin𝜉+Δ𝑥−𝑖Δ𝑡𝐻^∗sin𝜃,Δ𝑦𝜆−𝑎₄₄≤𝑎₄₁+𝑎₄₂=−𝑖Δ𝑡𝑇^∗sin𝜉+Δ𝑥−𝑖Δ𝑡𝑇^∗sin𝜃.Δ𝑦 (23)

Using the backwards triangular inequality, we take 𝐶<1/Δ𝑡, if the time-space steps and model parameters yield the following conditions: Δ𝑡Δ𝑥𝑢^∗+𝑔^′+Δ𝑡𝑣^∗𝑓Δ𝑦+Δ𝑡^∗+𝛾1+4Δ𝑡𝐾Δ𝑥²+1Δ𝑦²<1,Δ𝑡Δ𝑦𝑣^∗+𝑔^′+Δ𝑡𝑢^∗𝑓Δ𝑥+Δ𝑡^∗+𝛾1+4Δ𝑡𝐾Δ𝑥²+1Δ𝑦²<1,Δ𝑡Δ𝑥𝑢^∗+𝐻^∗+Δ𝑡Δ𝑦𝑣^∗+𝐻^∗<1,Δ𝑡Δ𝑥𝑢^∗_𝑠+𝑇^∗+Δ𝑡Δ𝑦𝑣^∗_𝑠+𝑇^∗|𝜆|1+Δ𝑡+4Δ𝑡𝜈Δ𝑥²+1Δ𝑦²<1. (24)

Then we have |𝜆|≤1+𝐶Δ𝑡, so the explicit finite difference scheme (13) is conditionally stable. Consequently, the time-step size satisfies the following inequality: 1Δ𝑡<min𝑎,1𝑏,1𝑐,1𝑑. (25) In which 𝑎=(1/Δ𝑥)(max|𝑢^∗|+𝑔^′)+(max|𝑣^∗|/Δ𝑦)+(|𝑓^∗|+|𝛾|)+4𝐾((1/Δ𝑥²)+(1/Δ𝑦²)), 𝑏=(1/Δ𝑦)(max|𝑣^∗|+𝑔^′)+(max|𝑢^∗|/Δ𝑥)+(|𝑓^∗|+|𝛾|)+4𝐾((1/Δ𝑥²)+(1/Δ𝑦²)), 𝑐=(1/Δ𝑥)(max|𝑢^∗|+|𝐻^∗|)+(1/Δ𝑦)(max|𝑣^∗|+|𝐻^∗|), 𝑑=(1/Δ𝑥)(max|𝑢^∗_𝑠|+max|𝑇^∗|)+(1/Δ𝑦)(max|𝑣^∗_𝑠|+max|𝑇^∗|)+|𝜆|+4𝜈((1/Δ𝑥²)+(1/Δ𝑦²)).

When the rotation, dissipation, and friction terms are neglected, we have Δ𝑡<minΔ𝑥max𝑢^∗+Δ𝑦max𝑣^∗,Δ𝑥max𝑢^∗_𝑠+Δ𝑦max𝑣^∗_𝑠. (26)

This is Courant-Friedrichs-Lewy (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 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.

Table 1 

The values of parameters in the model.

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.

5.1. Example 1

According to the stability criterion (25), we obtain Δ𝑡<14306s; from the CFL condition, we have Δ𝑡<189031s.

Running the model for one year, Figure 1 gives the results of zonal velocity 𝑢 and meridional velocity 𝑣 (at the latitude 0^∘ N) by changing with time-step Δ𝑡(dt).

Figure 1 

The figures show the values of the zonal velocity 𝑢 and the meridional velocity 𝑣 at the latitude 0∘ N, with dt1=600s; dt2=1200s; dt3=3600s; dt4=7200s; dt5=10800s; dt6=12000s.

When we take Δ𝑡=14306s, 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 Δ𝑡>14306s, 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)

Figure 2 

(a) shows the values of the zonal velocity 𝑢 at the latitude 0∘ N, and (b) at the longitude 180∘ E.

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 Δ𝑡=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 𝑣 (cm/s) with the variation of adjustable coefficient 𝜆 (lambda).

Figure 3 

The figure gives the results of SST (°C) with the variation of angle 𝜃 (theta), in which 𝜃1=0∘; 𝜃2=30∘; 𝜃3=45∘; 𝜃4=60∘; 𝜃5=90∘; 𝜃6=135∘.

Figure 4 

The figure gives the results of SSH with the variation of 𝛾 (gamma), in which 𝛾1=1.5×10−3day−1; 𝛾2=2.5×10−3day−1; 𝛾3=3.5×10−3day−1; 𝛾4=4.5×10−3day−1; 𝛾5=5.5×10−3day−1; 𝛾6=6.5×10−3day−1.

Figure 5 

The figure gives the results of meridional velocity 𝑢 with the variation of 𝐾, in which 𝐾1=2×10−3m2s−1; 𝐾2=3×10−3m2s−1; 𝐾3=4×10−3m2s−1; 𝐾4=5×10−3m2s−1; 𝐾5=6×10−3m2s−1; 𝐾6=7×10−3m2s−1.

Figure 6 

The figure shows the results of meridional velocity 𝑣 (cm/s) with the variation of 𝜆 (lambda), where 𝜆1=3×10−3day−1; 𝜆2=4×10−3day−1; 𝜆3=5×10−3day−1; 𝜆4=6×10−3day−1; 𝜆5=7×10−3day−1; 𝜆6=8×10−3day−1.

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.

Acknowledgments

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