# On the Stability of an Intermediate Coupled Ocean-Atmosphere Model

Authors: Tianxu Zhao; Guang-an Zou
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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 Section3, 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 flowu=(u,v) is governed by reduced-gravity shallow-water equations:
(1)∂u∂t+u∂u∂x+v∂u∂y-fv=-g′∂h∂x+τxρ0H+K(∂2u∂x2+∂2u∂y2)-γu,(2)∂v∂t+u∂v∂x+v∂v∂y+fu=-g′∂h∂y+τyρ0H+K(∂2v∂x2+∂2v∂y2)-γv,(3)∂h∂t+∂Hu∂x+∂Hv∂y=0,
in which u,v are the horizontal velocities; f=f0+βy is the Coriolis parameter of the equatorial betaplane approximation, where f0=2Ωsinθ, β=(2Ωcosθ)/r0 and r0 is radius of the earth; A density difference between the upper and lower layers results in a reduced gravity constant g′=g(ρ2-ρ1)/ρ0 (where g is the acceleration of gravity, ρ1 and ρ2 are the constant upper and lower layer densities, respectively, and ρ0 is a mean density); τx, τy are the wind stress over the sea surface; H=D+h is the upper layer thickness, where D is the undisturbed layer thickness, h is the interface displacement; K 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(4)∂T∂t+us∂T∂x+vs∂T∂y=ν(∂2T∂x2+∂2T∂y2)-λ(T-T0)+Qρ0CpH0-1H0weH(we)(T-Te),
where T is sea surface temperature (SST), us, vs are the horizontal velocities in the mixed layer, T0 is average SST, ν is horizontal heat diffusion coefficient, λ is adjustable coefficient, Q is downward heat flux, H0 is mean depth of mixed layer, we is the entrainment velocity, H(we) is a Heaviside step function of we (H(x)=1, if x>0, H(x)=0, if x<0) and Te is entrainment water temperature.In (4), the horizontal velocity us=(us,vs) in the surface mixed layer is separated into two components: us=u+ue(H-H0)/H, where u=(u,v) is upper-ocean flow velocity, surface Ekman flow ue=(ue,ve) is determined by Coriolis force and wind-stress force
(5)rsue-fve=τxρ0H0,rsve-fue=τyρ0H0,
where rs is the Rayleigh friction coefficient.The entrainment velocitywe is determined as divergence of surface flow we=H0∇·us and the temperature of entrained water Te is taked as a linear form:
(6)Te=0.5(T+T-50)+2∂T50∂zh′,
where T-50 is the observed mean temperature, ∂T50/∂z is the vertical gradient of T- at 50 m depth, and h′ 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)(τQ)=(τ-Q-)+(μ𝒜(T′)ωℬ(T′))+(ζχ),
where τ- and Q- 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)𝒜(T′)=(ψ,T′)=∫‍∫Ω‍ψ(x,y;x′,y′)T(x′,y′,t)′dx′dy′,ℬ(T′)=(ϕ,T′)=∫‍∫Ω‍ϕ(x,y;x′,y′)T(x′,y′,t)′dx′dy′,
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Δx in the x-direction, Δy in the y-direction, and Δt in the t-direction, and we define ui,jn=u(xi,yj,tn) with xi=x0+iΔx, yj=y0+jΔy, tn=nΔt for i,j=0,1,…, n=0,1,…. The governing equations will be solved on a staggered (Arakawa C-) grid. That is, the u, v are evaluated at the cell boundaries and the h and T 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(9)(δtα)i,j=αi,jn+1-αi,jnΔt,(10)(δxα)i,j=αi+1,jn-αi-1,jn2Δx,(δyα)i,j=αi,j+1n-αi,j-1n2Δy,(11)(δx2α)i,j=αi+1,jn-2αi,jn+αi-1,jnΔx2,(δy2α)i,j=αi,j+1n-2αi,jn+αi,j-1nΔy2.Therefore, the difference approximation of the partial differential equations (1)–(4) is given by
(12)(δtu)i+1/2,j+u¯*(δxu)i+1/2,j+v¯*(δyu)i+1/2,j-f*vi,j+1/2n=-g′(δxη)i,j+τxρ0H*+K((δx2u)i+1/2,j+(δy2u)i+1/2,j)-γui+1/2,jn,(δtv)i,j+1/2+u¯*(δxv)i,j+1/2+v¯*(δyv)i,j+1/2+f*ui+1/2,jn=-g′(δyη)i,j+τyρ0H*+K((δx2v)i,j+1/2+(δy2v)i,j+1/2)-γvi,j+1/2n,(δth)i,j+u¯*(δxh)i,j+v¯*(δyh)i,j+H*((δxu)i+1/2,j+(δyv)i,j+1/2)=0,(δtT)i,j+u¯s*(δxT)i,j+v¯s*(δyT)i,j=ν((δx2T)i,j+(δy2T)i,j)-λTi,jn-T*((δxu)i+1/2,j+(δyv)i,j+1/2)+C(T0,Q)
in which u¯*=(ui+1,jn+ui-1,jn)/2, v¯*=(vi,j+1n+vi,j-1n)/2, f*=f0+jβΔy, H*=D+hi,jn, u¯s*=u¯*+ue(H-H0)/H, v¯s*=v¯*+ve(H-H0)/H, T*=H(we)(Ti,jn-Te*), C(T0,Q)=λT0+(Q/ρ0CpH0).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 thatXi,jn=(ui+1/2,jn,vi,j+1/2n,hi,jn,Ti,jn)T, after some rearrangement, (12) can be written as
(13)Xi,jn+1=AXi+1,jn+CXi-1,jn+BXi,jn+DXi,j+1n+EXi,j-1n+en
in which en=(τx/ρH*,τy/ρH*,0,C(T0,Q))T,(14)A=(-Δtu¯*2Δx+ΔtKΔx20-Δtg′2Δx00-ru¯*2+rKΔx200-ΔtH*2Δx0-Δtu¯*2Δx0-ΔtT*2Δx00-Δtu¯s*2Δx+ΔtνΔx2),C=(Δtu¯*2Δx+ΔtKΔx20Δtg′2Δx00ru¯*2+rKΔx200ΔtH*2Δx0Δtu¯*2Δx0ΔtT*2Δx00Δtu¯s*2Δx+ΔtνΔx2),B=(b1Δtf*00-Δtf*b1000010000b2),(15)b1=1-Δtγ-2ΔtK(1Δx2+1Δy2),b2=1-Δtγ-2Δtν(1Δx2+1Δy2),(16)D=(-Δtv*2Δy+ΔtKΔy20000-Δtv*2Δy+ΔtKΔy2-Δtg′2Δy00-ΔtH*2Δy-Δtv¯*2Δy00-ΔtT*2Δy0-Δtvs*2Δy+ΔtνΔy2),E=(Δtv*2Δy+ΔtKΔy20000Δtv*2Δy+ΔtKΔy2Δtg′2Δy00ΔtH*2ΔyΔtv¯*2Δy00ΔtT*2Δy0Δtvs*2Δy+ΔtνΔy2).Definition 1.
The two-dimensional discrete Fourier transform ofu∈ℓ2 is the function u†∈L2([-π,π]×[-π,π]) defined by (see[30]):(17)u†(ξ,θ)=12π∑m,k=-∞∞‍exp(-i(mξ+kθ))um,k
for ξ,θ∈[-π,π].By using the discrete Fourier analysis to transform both sides of (9), we can obtain
(18)X†n+1(ξ,θ)=G(ξ,θ)X†n(ξ,θ)+en′,
where en′=(1/2π)∑m,k=-∞∞exp(-i(mξ+kθ))en, and the growth matrix
(19)G(ξ,θ)=Aexp(iξ)+Cexp(-iξ)+B+Dexp(iθ)+Eexp(-iθ)=(a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a44),
in which
(20)a11=a22=1-Δtγ-2ΔtK(1Δx2+1Δy2)+2ΔtK(cosξΔx2+cosθΔy2)-iΔt(u¯*sinξΔx+v¯*sinθΔy),a12=Δtf*,a13=-iΔtg′sinξΔx,a21=-Δtf*,a23=-iΔtg′sinθΔy,a14=a24=a34=a43=0,a31=-iΔtH*sinξΔx,a32=-iΔtH*sinθΔy,a33=1-iΔt(u¯*sinξΔx+v¯*sinθΔy),a41=-iΔtT*sinξΔx,a42=-iΔtT*sinθΔy,a44=1-Δtλ-2Δtν(1Δx2+1Δy2)+2Δtν(cosξΔx2+cosθΔy2)-iΔt(u¯s*sinξΔx+v¯s*sinθΔy).Theorem 2.
Letρ(G) be the spectral radius of the matrix G. If there exist constants Δx0, Δy0, Δt0, and C, independent of Δx, Δy, Δt, and ξ, θ, such that
(21)|ρ(G)|≤1+CΔt
for 0<Δt≤Δt0, 0<Δx≤Δx0, and 0<Δy≤Δy0, and for all ξ,θ∈[-π,π], then the scheme (18) is stable in the ℓ2 norm.Theorem 3 (Gerschgorin Circle Theorem, see [31]).
LetA be a complex n×n matrix, aij are elements of matrix A. For each eigenvalue λ of A, there exists an s such that
(22)|λ-ass|≤ρs,s=1,2,…,n,
where ρs=∑j=1,j≠sn|asj| are the sum of the absolute values of the elements in the sth row except for the diagonal element.By Gerschgorin Circle Theorem, ifλ is an eigenvalue of G, we have
(23)|λ-a11|≤|a12|+|a13|=|Δtf*|+|-iΔtg′sinξΔx|,|λ-a22|≤|a21|+|a23|=|-Δtf*|+|-iΔtg′sinθΔy|,|λ-a33|≤|a31|+|a32|=|-iΔtH*sinξΔx|+|-iΔtH*sinθΔy|,|λ-a44|≤|a41|+|a42|=|-iΔtT*sinξΔx|+|-iΔtT*sinθΔy|.Using the backwards triangular inequality, we takeC<1/Δt, if the time-space steps and model parameters yield the following conditions:
(24)ΔtΔx(|u¯*|+g′)+Δt|v¯*|Δy+Δt(|f*|+|γ|)+4ΔtK(1Δx2+1Δy2)<1,ΔtΔy(|v¯*|+g′)+Δt|u¯*|Δx+Δt(|f*|+|γ|)+4ΔtK(1Δx2+1Δy2)<1,ΔtΔx(|u¯*|+|H*|)+ΔtΔy(|v¯*|+|H*|)<1,ΔtΔx(|u¯s*|+|T*|)+ΔtΔy(|v¯s*|+|T*|)+Δt|λ|+4Δtν(1Δx2+1Δy2)<1.Then we have|λ|≤1+CΔt, so the explicit finite difference scheme (13) is conditionally stable. Consequently, the time-step size satisfies the following inequality:
(25)Δt<min{1a,1b,1c,1d}.
In which a=(1/Δx)(max|u¯*|+g′)+(max|v¯*|/Δy)+(|f*|+|γ|)+4K((1/Δx2)+(1/Δy2)), b=(1/Δy)(max|v¯*|+g′)+(max|u¯*|/Δx)+(|f*|+|γ|)+4K((1/Δx2)+(1/Δy2)), c=(1/Δx)(max|u¯*|+|H*|)+(1/Δy)(max|v¯*|+|H*|), d=(1/Δx)(max|u¯s*|+max|T*|)+(1/Δy)(max|v¯s*|+max|T*|)+|λ|+4ν((1/Δx2)+(1/Δy2)).When the rotation, dissipation, and friction terms are neglected, we have(26)Δt<min{(Δxmax|u¯*|+Δymax|v¯*|),(Δxmax|u¯s*|+Δymax|v¯s*|)}.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 f, the interfacial friction coefficient γ, the lateral eddy viscosity coefficient K, and the adjustable coefficient λ. Parameter values are shown in Table 1.Table 1
The values of parameters in the model.
Parameter
Value
Ω
7.29
×
1
0
-
5 s−1
D
100 m
K
2.5
×
1
0
-
4 m2 s−1
ν
2.5
×
1
0
-
4 m2 s−1
H
0
50 m
g
′
4.17 cm s−2
ρ
0
1000 kg m−3
γ
2.5
×
1
0
-
3 day−1
λ
5
×
1
0
-
3 day−1
C
p
4.2
×
1
0
3 J Kg−1 K−1The model domain extends form30° S to 30° N in latitude and from 120° E to 80° W in longitude. The model resolution is 2° in x direction and 1° in y 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 Δt<14306s; from the CFL condition, we have Δt<189031s.Running the model for one year, Figure1 gives the results of zonal velocity u and meridional velocity v (at the latitude 0° N) by changing with time-step Δt(dt).Figure 1
The figures show the values of the zonal velocityu and the meridional velocity v at the latitude 0° N, with dt1=600s; dt2=1200s; dt3=3600s; dt4=7200s; dt5=10800s; dt6=12000s.When we takeΔt=14306s, Figure 2 shows the results of zonal velocity u 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 Δt>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) shows the values of the zonal velocityu at the latitude 0° N, and (b) at the longitude 180° E.
(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 time-step sizeΔt=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.Figure3 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 u with the variation of lateral eddy viscosity coefficient K. Figure 6 exhibits the meridional velocity v (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 velocityu with the variation of K, in which K1=2×10-3m2s-1; K2=3×10-3m2s-1; K3=4×10-3m2s-1; K4=5×10-3m2s-1; K5=6×10-3m2s-1; K6=7×10-3m2s-1.Figure 6
The figure shows the results of meridional velocityv (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.


## 5.1. Example 1

According to the stability criterion (25), we obtain Δt<14306s; from the CFL condition, we have Δt<189031s.Running the model for one year, Figure1 gives the results of zonal velocity u and meridional velocity v (at the latitude 0° N) by changing with time-step Δt(dt).Figure 1
The figures show the values of the zonal velocityu and the meridional velocity v at the latitude 0° N, with dt1=600s; dt2=1200s; dt3=3600s; dt4=7200s; dt5=10800s; dt6=12000s.When we takeΔt=14306s, Figure 2 shows the results of zonal velocity u 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 Δt>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) shows the values of the zonal velocityu at the latitude 0° N, and (b) at the longitude 180° E.
(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 time-step sizeΔt=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.Figure3 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 u with the variation of lateral eddy viscosity coefficient K. Figure 6 exhibits the meridional velocity v (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 velocityu with the variation of K, in which K1=2×10-3m2s-1; K2=3×10-3m2s-1; K3=4×10-3m2s-1; K4=5×10-3m2s-1; K5=6×10-3m2s-1; K6=7×10-3m2s-1.Figure 6
The figure shows the results of meridional velocityv (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.



---

*Source: 658345-2013-12-03.xml*