Li et al. Progress in Earth and Planetary Science
https://doi.org/10.1186/s40645-021-00432-9(2021) 8:43

Progress in Earth andPlanetary Science

RESEARCH ARTICLE

Open Access

The vertical structure of annual wave
energy flux in the tropical Indian Ocean
Zimeng Li1* , Hidenori Aiki2,3*, Motoki Nagura4 and Tomomichi Ogata3

Abstract
A recently developed energy flux diagnosis scheme, which incorporates a smooth connection between the tropicaland subtropical zones, is used in the present study to investigate vertically propagating waves in the tropical Indian
Ocean (IO) based on the result of a linear, continuously stratified ocean model driven by climatological wind
forcing. This extended diagnosis reveals deep-reaching eastward energy fluxes at the equator which develop fourtimes per year and are associated with equatorial Kelvin waves (KWs) generated by semiannual winds. The authors
find that the downward transfer of wave energy is particularly deep in the southern Bay of Bengal (SBoB). This
downward flux is attributed to off-equatorial Rossby waves and appears four times per year, maximizing its
amplitude during November–December. Southwesterly winds in the Arabian Sea intensify eastward energy flux of
KWs at mid-depth, which maximizes in amplitude in August. This is contrastive to KW energy flux at the surface
which peaks in May. These mid-depth equatorial KW packets subsequently arrive at the eastern boundary of the IO
and are diffracted poleward to produce downward energy flux in November and December detected in the SBoB.
Keywords: Indian Ocean, Kelvin wave, Rossby wave, Seasonal variation, Vertical structure, Wave energy

1 Introduction
In the tropical regions, the dynamics of the atmosphere–
ocean interactions dictating the growth and morphologyof large-scale surface currents strongly depend on the subsurface oceanic structure. In turn, variability in the surface
currents themselves can induce deep and persistent flows
which remain trapped along the equator (Luyten andSwallow 1976; Hayes and Milburn 1979; McCreary 1984;
Rothstein et al. 1985; Kessler and McCreary 1993; Dewitte
and Reverdin 2000). In the central equatorial Indian
Ocean (IO), wind forcing develops on a semiannual timescale. Further to the west and in the regions off the equator, the monsoonal phase transition produces a seasonal
cycle in wind forcing (Reppin et al. 1999; Schott andMcCreary 2001). The former semiannual wind forcing induces equatorial Kelvin waves (KWs) with a downwelling* Correspondence: lizimeng1995@gmail.com; aiki@nagoya-u.jp
1
Graduate School of Environmental Studies, Nagoya University, Furo-cho,
Chikusa-ku, Nagoya 464-8601, Japan
2
Institute for Space-Earth Environmental Research, Nagoya University,
Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan
Full list of author information is available at the end of the article

phase in the surface layer in both boreal spring and autumn and a corresponding upwelling phase in both boreal
winter and summer, leading to an eastward transfer ofwave energy four times per year (Li and Aiki 2020, hereafter LA20). Striking the eastern boundary of the tropical
IO, KWs are first retro-reflected, then diffracted poleward
and subsequently radiate westward-propagating equatorial
and midlatitude Rossby waves (RWs) (Mark and Huang
2004; Han 2005; Fu 2007).
Previous studies have investigated the vertical propagationof equatorial KWs and RWs in the Pacific and Atlantic
Oceans on intraseasonal to interannual time scales (McPhaden 1982; Lucas and Firing 1985; McPhaden and Ripa 1990;
Kessler and McCreary 1993; Brandt and Eden 2005; Thierry
et al. 2006; Ishizaki et al. 2014). Vertically propagating waves
are also observed in the tropical IO. Owing to the dominanceof semiannual wind forcing, variability in mid-depth of theequatorial IO below the main pycnocline is characterized bysemiannual variability (Luyten and Roemmich 1982; Naguraand McPhaden 2016; Huang et al. 2018a, b; Zanowski andJohnson 2019; Chen et al. 2020). They consist of equatorial
KWs and long RWs, which carry energy to the east and to© The Author(s). 2021 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,
which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate ifchanges were made. The images or other third party material in this article are included in the article's Creative Commonslicence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commonslicence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain
permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

(2021) 8:43

Li et al. Progress in Earth and Planetary Science

Page 2 of 19

the depth, and to the west and to the depth, respectively.
Waves are reflected at the eastern and western boundaries of
the IO, and boundary-generated waves further propagate
downward (Nagura and McPhaden 2016; Huang et al. 2019).
In particular, mooring observations in the southern Bay of
Bengal (SBoB) at 5° N, 90.5° E show a robust half-year cycle
between 200 and 1200 m depth (Huang et al. 2019). The
subsurface mooring at the same location used in Chen et al.
(2017) reveals remarkable intraseasonal meridional currents
across SBoB that are attributed to westward-propagating
RWs. In the off-equatorial region, variability below thepycnocline has annual periodicity, owing to the predominance of annual wind forcing related to monsoon (Johnson
2011; Nagura 2018; Huang et al. 2019).
These previous studies examined vertical propagation ofwaves using ray tracing technique and/or vertical mode
decomposition. Ray tracing technique requires the assumption that the WKB approximation is valid (i.e., thevertical scale of variation in buoyancy frequency is larger
than vertical wavelength of waves). This assumption is
generally satisfied below the pycnocline, but not in andabove the pycnocline. Owing to this, energy propagation
cannot be traced from the surface, where winds excite
wave energy, to the depth below the pycnocline using the
ray tracing technique. This raises an essential difficulty ininterpreting the generation mechanism of mid-depth
waves. Vertical mode decomposition does not require
such an assumption, but vertically propagating waves consist of multiple baroclinic modes (Table 1), and examination of each baroclinic mode is not informative about
vertical propagation. An alternative is sensitivity experiments using a numerical model, but this method is relatively expensive computationally.
These problems can be solved by the use of energy
flux newly defined by Aiki et al. (2017; hereafter
AGC17) and Aiki et al. (under revision in Journal of the
Atmospheric Sciences, hereafter JAS2021). Their scheme
does not use the WKB approximation and is applicableseamlessly to depths above, in, and below the pycnocline. It does not require quasi-geostrophy and covers
variability both in the equatorial and subtropical regions.

It does not require Fourier transform either and thus
can be used in coastal regions. The scheme is diagnostic
and provides a computationally cheaper way comparedwith sensitivity experiments using an elaborated numerical model. The AGC17 scheme allowed Ogata and Aiki
(2019), LA20, Song and Aiki (2020), and Toyoda et al.
(2021) to obtain wholly new perspectives on the horizontal transfer routes of waves in the tropical oceans.
We extend their work to a three-dimensional systemand describe vertical wave propagation.
The depth-dependent version of the seamless AGC17
scheme has not been applied to the analysis of oceanic
waves before the present study. Our initial focus will
mainly be on seasonal variability in the IO, while interannual variations will be investigated in the future work.
The manuscript is laid out in the following fashion: Section 2 describes the formulation of an adopted model,
method for energy flux computation, and temporal andspatial characteristics of simulated fields. Section 3 provides an analysis of the vertical structure of energy fluxes
associated with KWs and RWs, while Section 4 gives asummary of the main results.

2 Materials and methods
AGC17 employed a single-layer system based on theshallow-water equations to derive depth-independent expressions of the energy flux. As it has the advantage of asmooth tropical–subtropical transition, the AGC17
scheme can effectively determine the horizontal distribution of energy flux pointing in the direction of group velocity of waves. JAS2021 has extended the formulation of
AGC17 to cover the case of a baroclinic system in the atmosphere. Their essential contribution is the derivation of
the exact expression of the vertical component of energy
flux (this is not described in the present manuscript). The
expressions given below for the horizontal component of
depth-dependent energy flux represent a straightforward
extension of the content of AGC17 and are summarized
below by the adoption of an oceanic framework that thegeneral reader should find convenient.

Table 1 List of experiment parameters. The basin mode period is defined as 4L/c(n) where L = 6110 km is the zonal distance of thepﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Indian Ocean at the equator (from 45 to 100°E). The wind-coupling thickness is defined as hmix Hbottom =αðnÞ , where hmix= 35 m and
Hbottom= 5500 m
Baroclinic mode

1st

2nd

3rd

4th

5th

6th

Gravity wave speed (m/s)

2.63

1.57

0.956

0.706

0.557

0.467

Equatorial deformation radius (km)

339

262

205

176

156

143

Equatorial inertial period (day)

9.4

12.1

15.6

18.1

20.4

22.3

Basin mode period (year)

0.29

0.49

0.81

1.1

1.39

1.66

α(n) (non-dimensional)

0.39

0.41

0.27

0.29

0.23

0.22

Wind-coupling thickness (m)

1123

1062

1644

1505

1933

1955

Li et al. Progress in Earth and Planetary Science(2021) 8:43

Page 3 of 19

2.1 Formulation of wave energy flux

In order to investigate the vertical structure of energy
transfers, we aim to derive depth-dependent expressionsof energy flux. We use the primitive equations for a Boussinesq fluid on an equatorial β plane (Gill 1982) appropriate for linear waves propagating in a rotating frame ofreference and with zero mean flow. Let x, y, z, and t define
independent variables in a Cartesian coordinate system
and increase eastward, northward, and vertically upward
respectively. The quantities u, v, and w represent the corresponding three-dimensional components of velocity (afull list of variables is given in Table 2). The equations
may then be written as
∂u
∂Φ
−fv þ
¼ 0;
∂t
∂x

ð1aÞ

∂v
∂Φ
þ fu þ
¼ 0;
∂t
∂y

ð1bÞ

∂Φ
ρ
¼ −g ;∂zρ0

ð1cÞ

−1 ∂2 Φ;N 2 ∂z∂t

ð1dÞ

g dρ
;
ρ0 dz

ð1eÞ

w¼

N2 ¼ −∂u ∂v ∂w
þ þ
¼ 0;
∂x ∂y ∂z

ð1fÞTable 2 List of symbols used in Section 2.1, where A is arbitrary
quantities written dimensionally. The coordinate axes x, y, z
increase eastward, northward, and vertically upward, respectively
f = f0 + βyCoriolis parameter
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
N ¼ −ðg=ρ0 Þðdρ=dzÞ Buoyancy frequency
x, y, z, t

Cartesian coordinates〈〈a, b, c〉〉

Vector with x, y and z-axis components a, b,
and c

U ≡ 〈〈u, v, w〉〉

Velocity vector∇ ≡ 〈〈∂x, ∂y, ∂z〉〉

Gradient operatorΦ

Geopotentialq

Potential vorticity q ≡ vx − uy + [(f/N2)Φz]z

φ

Solution of (6a)φapp

Solution of (7a)ðu2 þ v 2 þ Φ2z =N2 Þ=2Wave energy density

μ = kx + mz − ωtWave phasek

Zonal wavenumbermVertical wavenumberω

Wave frequency

A

Phase average of A

where f = f0 + βy is the Coriolis parameter, Φ = Φ(x, y, z,t) corresponds to geopotential, and ρ = ρ(x, y, z, t) is theperturbation of potential density. The symbols N = N(z)and g are the buoyancy frequency and the acceleration
owing to gravity, respectively. The symbols ρ0 and ρðzÞrepresent the reference density and the backgrounddensity, respectively. Manipulation of (1a)–(1f) yields aprognostic equation for potential vorticity (PV: symbolized as q) to read∂ ∂v ∂u ∂f ∂Φ
þvβ ¼ 0:− þ∂t ∂x ∂y ∂z N 2 ∂z
|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}

ð2Þ

≡q

This relation is generally applicable to waves at all latitudes, such as midlatitude RWs, midlatitude inertia–
gravity waves (IGWs), and equatorial waves (i.e., equatorial RWs and IGWs, equatorial Rossby-gravity waves(RGWs, i.e., Yanai waves) and equatorial KWs (Matsuno
1966; Yanai and Maruyama 1966), if we applied it to an
equatorial β plane (f0=0)). Equatorial KWs (i.e., v = 0) arecharacterized by q = 0. From (1a)–(1c), the wave energy
density may be represented in prognostic fashion by" 2 #∂1 2
1∂Φu þ v2 þ 2
∂t 2
∂z
N

þ ∇∙ uΦ; vΦ; wΦ
¼ 0;

ð3aÞ

∂ ∂ ∂
; ∂y ; ∂zii and the overbar signifies phase
where ∇ ≡ hh∂x
averaging (i.e., for a sinusoidal wave, A ¼ 0 for A = u, v,
w and Φ), or low-pass time filtering (considering thiscase, we retain the local time derivative in (3a), thereby
permitting slow time variations in the general case). Theenergy flux vector in (3a)

UΦ ¼



uΦ; vΦ; wΦ ;ð3bÞis codirectional with the group velocity of midlatitude
IGWs in the atmosphere and ocean (Gill 1982). However, the energy flux of (3b) is not aligned with thegroup velocity of midlatitude RWs (Longuet-Higgins
1964).
In order to obtain an expression of the energy flux thatis codirectional with group velocity of all wave types,
AGC17 (and also JAS2021) has revisited the equations
to define energy flux of both equatorial and midlatitude
waves. The derivation is briefly described below using
linear waves on an equatorial β plane and zero mean
flow. Let k, m, and ω be zonal and vertical wavenumbers
and wave frequency, respectively. The zonal componentof group velocity of equatorial waves may be written as

Li et al. Progress in Earth and Planetary Science∂ω
2ω2 þ ω=k;
¼ 3 2
∂k 2ω m = kN 2 þ 1(2021) 8:43

Page 4 of 19ð4Þ

which has been written in a dimensional form and m denotes vertical wavenumber.
The group velocity times wave energy density is rewritten in a depth-dependent framework as" 2 #∂ω 1 21∂Φu þ v2 þ 2
∂k 2
∂z
N

!
∂ Φφ 1 ∂2 u
¼ uΦ þ
φ ;þ
∂y 2β ∂t 2φ≡−vμ;k þ 2ω3 m2 = βN 2

ð5aÞ
ð5bÞ

where the scalar quantity φ may be referred to as
pseudo-streamfunction and the subscript μ signifies partial derivative with respect to wave phase (details are
given in Supplementary Material of this manuscript).
We note that, for any type of equatorial waves, there is a
precise identity between the meridional profile of zonal
energy flux on the right-hand side of (5a) and the quantity ð∂ω=∂kÞðu2 þ v2 þ Φ2z =N 2 Þ=2, which is illustrated inJAS2021. The definition of φ in (5b) makes use of a Fourier expansion. In a depth-dependent form, (5b) can bereformulated such that the parameters k, m, and ω areabsent, which yields an inversion equation for PV to
read

∂2
∂2
φ
þ
∂x2 ∂y2  

∂
1 ∂
∂2
2
þ
f φþ3 2φ∂t∂z N 2 ∂z
¼ q;



UΦ þ




∂2
∂2
∂ f 2 ∂ app
app
þ
þ
φ
¼ q;
φ
∂x2 ∂y2∂z N 2 ∂z

ð6aÞ

!
! ++
∂ Φφ 1 ∂2 u
∂ Φφ 1 ∂2 u
þ
þ
φ ;−
φ ;0
;
∂y 2
β ∂t 2∂x 2
β ∂t 2

ð6bÞ
whose horizontal component is referred to as the level 0(i.e. exact) expression for waves at all latitudes. The vertical component of (6b) represents a rough approximation of the vertical group velocity times wave energy forequatorial waves and its exact expression for midlatitude
waves. See JAS2021 for the details of (6a)–(6b).
AGC17 considered eliminating the term representing
the second-order time derivative from their exact inversion equation for PV. This corresponds to (6a) in the

ð7aÞ

which may be solved using the set of a Dirichlet boundary condition (φapp = 0) at coastlines and bottom topography and a Neumann boundary condition (∂φapp/∂z = 0)
at the sea surface. What AGC17 has called the level-2(i.e., approximate) expression for energy flux is extended
here to a depth-dependent form as
UΦ þ





∂ φapp Φ∂ φapp Φ
;−
;0
∂y
2
∂x
2

;

ð7bÞ

which has been used in the overall model analyses of thepresent study.
Equation (7b) may be applied also to equatorial KWs
for which UΦ becomes identical to the product of groupvelocity and wave energy density. Namely, no term for
additional rotational flux is present. A characteristic of
KWs is v = 0, which yields q = 0 following (2). Since theright-hand side of (7a) vanishes, the expression of energy
flux for KWs given in (7b) reduces to UΦ. This result is
consistent with the nature of gravity waves.
We consider the Helmholtz decomposition of thelevel-2 wave energy flux to read,
−

∂
∂
∂ app
φ Φ=2 ;
R− D ¼ ρ0 uΦ þ ρ0
∂y
∂x
∂y

ð8aÞ

∂
∂
∂ app
φ Φ=2 ;
R− D ¼ ρ0 vΦ−ρ0
∂x
∂y
∂x

ð8bÞ

þ

where q has been defined in (2). The three-dimensional
components of wave energy flux are written in a vector
form as,**present study. Let φapp be the solution of the approximated PV inversion equation to read,where the scalar quantities R corresponds to the rotational component of energy flux and is referred to as energy flux streamfunction. The quantity D in (8a)–(8b) is
associated with the divergent components of energy fluxand is referred to as energy flux potential. The energy
flux potential becomes positive in regions where winds
supply wave energy (see Appendix for the expression ofenergy input by winds). The energy flux potential becomes negative in regions where wave energy is dissipated. As shown by LA20, the distribution of potential issmooth and nicely shows large-scale characteristics, incontrast to the distributions of the wind input and dissipation terms that manifest small-scale and locally signindefinite signals. This is the reason why the present
study recommends to use the distribution of the energy
flux potential, instead of its pre-inversion terms. We
shall show the distribution of the energy flux potential
and streamfunction, as given by (8a)–(8b), in the IO
later in the manuscript.

Li et al. Progress in Earth and Planetary Science(2021) 8:43

2.2 Model setup

We simulated climatological variability in the IO using aclassical linear ocean model (LOM, see Appendix for details). This model yields both the three-dimensional andbaroclinic-mode profiles of velocity vector and geopotential, that is aimed at facilitating our analysis of simulated wave dynamics. This kind of model has been usedby previous studies to examine vertical propagation ofwaves (Huang et al. 2018a, 2019, 2020; Chen et al. 2020)
and is dynamically consistent with our formulation ofenergy flux that is derived from linear equations (Eq.
(1a)–(1f)). The wind forcing for the model is derived
from the 10 m wind data of the European Centre for
Medium-Range Weather Forecasts Reanalysis (ECMWF
ERA-40) by applying the Large and Pond (1981) bulk
formula to wind data from September 1957 to August
2002. While the experiments described in LA20 assumed
no motion near the ocean bottom and examined the first
three baroclinic modes, here, we assume a free-slip bottom boundary condition (Chelton and Schlax 1996) asin traditional studies and focus on the first six baroclinic
modes. This change of model setup is required in order
to investigate the vertical propagation of waves in thepresent study. Each of our six numerical experiments involved integration over a period of 20 years (inclusive of
the spin-up phase from a state of no motion). Investigation of the fundamental baroclinic modes begins when
the system reaches a confirmed climatological equilibrium. The subsequent analysis is then based on the final
100 “snapshots” of model output. Since these are separated by 3.65 days, we make use of data from the lastyear of each run. Below, all results are shown in thedepth coordinate, which are reconstructed from themodel outputs of the first six baroclinic modes. The result of the LOM experiment is justified in the next subsection by comparing with the output from a highresolution climatological Ocean General Circulation
Model (OGCM) simulation for the semi-global ocean.
We have used monthly current and wind stress data derived from the OGCM For the Earth Simulator (OFES)
forced by National Centers for Environmental Prediction
reanalysis winds (Masumoto et al. 2004). OFES is on thebasis of Modular Ocean Model version 3. The output ismonthly mean during the last 6 years (from 45 to 50
years) of 50 year climatological integration with a horizontal resolution of 0.1°× 0.1° and 54 vertical levels.
2.3 Simulated currents and waves

We first briefly outline some key features from the LOM
experiment. The model reproduces the Somali jet, whichflows along the East African coast and is driven by largescale seasonal monsoonal winds, and the Wyrtki (1973)
jet (e.g., Nagura and McPhaden 2010a), which represents
an along-equator response to semiannual westerly winds

Page 5 of 19

during the monsoon transition period (Fig. 1). To investigate how the monsoon circulation affects velocity variability, we have constructed the longitude–depth
sections of zonal velocity along the equator (left panels
in Fig. 2) and the latitude–depth sections of zonal velocity at 70° E (right panels in Fig. 2). Semiannual winds
induce east-upward tilt in the phase line of zonal velocity (yellow line in Fig. 2a, e), which can explain theformation of the Wyrtki jet in the following monsoon
transition period. A semiannual variation of baroclinic
flow is dominant below about 200 m depth along theequator in response to the surface-layer Wyrtki jets. An
eastward undercurrent near 150 m depth along theequator in February (Fig. 2a), referred to as Equatorial
Undercurrent, is generated by equatorial KWs and long
RWs, as investigated by Chen et al. (2015a, 2019). Simulated zonal velocity along 70° E is trapped at the equatorwith a vertical phase propagation, which coincides withthe latitude-depth plots of climatological zonal currents
along 80° E using ocean reanalysis data in Huang et al.
(2018b). They pointed out that vertically propagating
RWs generated by the reflection of wind-forced KWs atthe eastern boundary of IO are responsible for the vertical phase propagation of mid-depth zonal currents on
the basis of mooring and reanalysis data, and this conclusion was further confirmed by Zanowski and Johnson(2019) using velocity and isotherm displacement obtained from Argo floats.
We compare the seasonal anomalies of current and windstress derived from OFES (Fig. 3a, b) and our climatologicalmodel experiments (Fig. 3c, d). The Hovmöller longitudetime diagrams along the equator at the sea surface indicatethat the eastward-propagating KW signals are less pronounced in surface zonal velocity and more pronounced ingeopotential anomaly (Fig. 3a, c, e). Semiannual winds cause
velocity reversals four times per year. This is consistent with
Yu and McPhaden (1999), Ding et al. (2009) and Nagura and
McPhaden (2010a, b), who pointed out that variability in surface zonal velocity along the equator is dominated by RWs,
and that in SSH by equatorial KWs. Geopotential anomalies
are positive near the surface during boreal spring and fall,
which are associated with wind-forced downwelling KWs.
Both the horizontal section at the surface (Fig. 3a, c) and thevertical section in the upper 2000 m (Fig. 3b, d) of zonal velocity derived from OGCM (Fig. 3a, b) and simulated by
LOM (Fig. 3c, d) show upward phase signals and reversals of
zonal currents four times per year with a noticeable semiannual cycle. These results agree with Huang et al. (2018a) and
Chen et al. (2020), who demonstrated that the semiannual
cycle is the predominant feature of mid-depth zonal currents
in the equatorial IO. Geopotential anomalies (colors in Fig.
3f) are also characterized by the semiannual cycle with an
upward phase propagation, which is attributable to windgenerated upwelling (during the monsoon seasons) andLi et al. Progress in Earth and Planetary Science(2021) 8:43

Page 6 of 19

Fig. 1 Left panels: the zonal component of simulated velocity at the sea surface (color shading with a unit of m/s), the zonal component of wind stress in themodel (solid and dotted contours for eastward and westward anomalies, respectively, with an interval of 0.02 N/m2), and the velocity vector (arrows in units ofm/s) in the middle of a February, c May, e August, and g November from “Year 20” of our climatological model experiments. When under 0.2 m/s, themagnitude of simulated velocity vectors in left panels is indicated either by arrow length or, for greater values, by heavy arrows of uniform length. Right panels(b, d, f, and h): same as left panels except for simulated sea surface geopotential anomaly (color shading with a unit of m2/s2) and PV-inverted function at thesea surface (solid and dotted contours represent positive and negative values, respectively, with an interval of 20,000 m2/s). The PV-inverted function may beinterpreted as geostrophic streamfunction in off-equatorial regions. The sea surface quantities have been reconstructed from the model results of fundamental
six baroclinic modes (Appendix)

downwelling (during the monsoon transition period) KWs.
The propagation of waves has a deterministic role on thevertical profile of geopotential anomalies (see Figure S1 in
Supplementary Material). The semiannual signals of

geopotential anomalies are weaker than that of zonal velocity,which are attributed to the interaction between KWs and
RWs in the central Indian Ocean (Yuan and Han 2006). The
upward phase signals of both zonal velocity and geopotential

Li et al. Progress in Earth and Planetary Science(2021) 8:43

Page 7 of 19

Fig. 2 The zonal component of simulated velocity (color shading with a unit of m/s) and the vertical component of simulated velocity (solid- anddotted-contours indicate positive and negative values, respectively, with an interval of 4 × 10−6 m/s) shown by snapshots in the middle of a, b
February, c, d May, e, f August, and g, h November from “Year 20” of our climatological model experiments. Zonal-vertical section at the equator(left panels) and meridional-vertical section at 70° E (right panels)

anomaly shift gradually in depth, and variability in the deeper
layer leads to that in the upper layer. This is regarded as theindication of downward energy transfers by oceanic waves(McCreary 1984; Miyama et al. 2006; Pujiana and McPhaden2020).

3 Results and discussions
This paper presents results from the first oceanic application of the depth-dependent version of the AGC17
scheme. In the following, all results associated with the
AGC17 scheme were obtained using this extendedLi et al. Progress in Earth and Planetary ScienceFig. 3 (See legend on next page.)(2021) 8:43

Page 8 of 19

Li et al. Progress in Earth and Planetary Science(2021) 8:43

Page 9 of 19(See figure on previous page.)
Fig. 3 Hovmöller zonal–time diagrams at the sea surface of the equator (left panels) and depth–time diagrams at 70°E of the equator (right panels) from
OGCM OFES (a, b) and “Year 20” of our climatological model experiments (c–f), respectively. Color shading in the panels (a–d) shows the output of zonal
component of velocity based on OGCM (a, b) and LOM (c, d) with a unit of m/s. Color shading in the panels (e–f) shows simulated geopotential anomaly witha unit of m2/s2. Panels (a, c) include the zonal component of wind stress (solid and dotted contours for positive and negative values, respectively, with aninterval of 0.1 N/m2) from OGCM and LOM, respectively. Panels (b, d) include the meridional component of velocity based on OGCM and LOM (solid- anddotted-line contours for positive and negative values, respectively, with an interval of 0.02 m/s). Panel (f) includes the vertical component of simulated velocity(solid and dotted contours for positive and negative values, respectively, with an interval of 4 × 10−6 m/s). Solid blue lines in the panels (c, e) indicate thetheoretical phase speeds of both Kelvin waves and equatorial long Rossby waves in the second baroclinic modeversion. We use the pressure flux (hereafter PF) scheme,
as given by ρ0 UΦ , for comparison with the AGC17
scheme, which is identical to flux defined by Eq. (3b) except for ρ0. While the PF scheme yields meridionally alternating energy flux vectors (Fig. 4a), the AGC17
scheme is able to show basin-wide cyclonic patterns,
which represent the depth-integrated transfer routes of
KW and RW energy (Fig. 4b). These waves are driven bywind energy input in the vicinity of the Maldives (from60°E to 80°E), the Arabian Sea and the southern Bay of
Bengal (SBoB), as shown by the positive values ofenergy-flux potential (indicated by the contours of Fig.
4b). Near the African coasts, energy flux streamfunction
derived from the AGC17 scheme (color shadings in Fig.
4b) reveals cyclonic circulation that is localized in nature, whereas on the east side of the basin, the cyclonic
circulation develops on significantly larger spatial scales
in each hemisphere, as identified by LA20.
In terms of model performance at the key depths of 550
m and 1100 m, the PF scheme fails to reveal important energy fluxes with distinct eastward KW signatures andsimilarly misses a series of westward signals associatedwith midlatitude RWs (Fig. 4c, e). In contrast, the AGC17
scheme is able to show the eastward transfer of wave energy along the equator associated with equatorial KWs
and also the westward transfer of wave energy in offequatorial regions associated with midlatitude RWs (Fig.
4d, f). Note that pressure variability along the equator represents KWs and zonal velocity variability along the equator are caused by RWs (Yu and McPhaden 1999; Ding
et al. 2009; Nagura and McPhaden 2010a, b). The AGC17
scheme is appropriate when tracking the group velocity ofbaroclinic waves at all latitudes, particularly as it incorporates a realistic smooth transitioning between the tropicaland subtropical regions. The contours in Fig. 4d and fshow the vertical component of energy flux. Note that
downward energy propagation is large in magnitude in
the SBoB, the south Arabian Sea and the region between
the equator and 5° S in 55°–75° E. In particular, the downward transfer of wave energy reaches its greatest depth in
SBoB, as is described below.
Along the latitude-depth plane at 70° E, the westward
flux of RW energy is symmetrically distributed acrossthe equator (color shading in Fig. 5a). The vertical fluxof wave energy at 90° E (contours in Fig. 5b) reaches a
deep level at around 3°–5° N, corresponding to thedownward flux of RW energy in the SBoB (Fig. 4d, f).
The characteristics of this downward flux1 are investigated by comparing time–depth diagrams at 70° E and90° E along 3° N (Fig. 5c, d). The zonal component ofenergy flux based on the AGC17 scheme shows negative
peaks four times per year, which is more distinct at 90°
E than at 70° E (color shadings in Fig. 5c, d). The verticalcomponent of energy flux at 90° E (contours in Fig. 5d)
also shows negative peaks four times per year that become deepest in November–December. These features
can be generated by either wind forcing in the SBoB orby the reflection and diffraction of equatorial KWs as
they encounter the eastern boundary, as will be discussed in the next paragraph. The time–depth diagramsof the energy flux at 70° E and 90° E on the equator (Fig.
5e, f) show peaks of eastward energy flux when zonal
velocity exhibits upward phase propagation (Fig. 3d).
The Hovmöller longitude–time plots at the equator(panels on the left in Fig. 6) and at 3° N (panels on theright) shows zonal energy flux derived from the AGC17
scheme (color shadings). Zonal energy flux along theequator shows four beams per year at both the surface
and at 550 m depth, which suggests a group velocity
consistent with the second baroclinic mode KWs (bluesolid lines). In contrast, the westward transfers of wave
energy at the equator and at 3° N indicate group velocities consistent with both reflected equatorial RWs anddiffracted midlatitude RWs. The surface eastward flux
along the equator, which is intense in both February andMay, is associated with upwelling and downwelling
KWs, respectively (Fig. 6a). Energy input rate by wind
forcing (Fig. 7) is predominant along the African coast
during the northeast monsoon and southwest monsoon
than that of during boreal spring and fall. Monsoonal
winds excite the upwelling KWs along the African coastand bring the divergence towards the equator. It is1

The vertical energy flux at 103° E, 6° S shows a similar deepest
reaching pattern in November–December (not shown), whichindicates that the reflection and diffraction of equatorial KWs affectboth Northern and Southern Hemispheres.

Li et al. Progress in Earth and Planetary Science(2021) 8:43

Page 10 of 19

Fig. 4 a, b Depth-integrated energy flux streamfunctions (color shaded in units of 109 W = GW), depth-integrated energy flux potentials (solid
contours for positive values and dotted contours for negative values; contour interval = 0.05 GW), and depth-integrated energy flux vectors(arrows in units of W/m). c, d The zonal component of the energy flux (color shaded in units of W/m2) and the energy flux vector (arrows inunits of W/m2) at 550 m depth. e, f Same as c and d except that depth is now 1100 m. All quantities represent annual mean values taken from“Year 20” of the model experiments, and estimated by the pressure flux scheme (left panels) and the depth-dependent version of the Aiki,
Greatbatch, and Claus Level-2 scheme (right panels). Panel (b) includes two groups of cyclonic circulations of wave energy shown by green
arrows that are demonstrated in LA20. Panels (d and f) include the vertical component of the energy flux, shown by contours (solid for positiveand dotted for negative values) at an interval of 10−4 W/m2. When under 1000 W/m, the magnitude of the depth-integrated energy flux vectorsin (a, b) is indicated either by arrow length or, for greater values, by heavy arrows of uniform length. Arrows in the series c, d, e, and f aredisplayed in a similar fashion to (a) and (b), with the magnitude of the energy flux indicated either by arrow length (when under 2 W/m2) or byheavy arrows of uniform length for greater valuesLi et al. Progress in Earth and Planetary ScienceFig. 5 (See legend on next page.)(2021) 8:43

Page 11 of 19

Li et al. Progress in Earth and Planetary Science(2021) 8:43

Page 12 of 19(See figure on previous page.)
Fig. 5 The zonal component of the energy flux (color shaded in units of W/m2) and the vertical component of the energy flux (solid contours for
positive values and dotted contours for negative values (a, b) with an interval of 10−4 W/m2, (c–f) at magnitude levels of 2, 5, 10, 20, 30, 50
×10−4 W/m2) as estimated from the depth-dependent version of the Aiki, Greatbatch, and Claus Level-2 scheme, as given by (7b). Annual meansin the meridional-vertical sections are shown in a at 70° E and b at 90° E. Hovmöller time–depth diagrams at 3° N for c 70° E and d 90° E and on
the equator for e 70° E and f 90° E. All quantities are derived from “Year 20” of the climatological model experiments associated with first six
fundamental baroclinic modes

consistent with the remarkable energy input by wind
forcing at the African coast during boreal winter andsummer at the equator and 3° N at the surface. In contrast, the eastward flux at 550 m depth along the equatoris intense in both May and August (Fig. 6c) and reaches
1100 m depth in November (Fig. 6e), as shown by thetime–depth diagram of Fig. 5e.
The poleward diffraction of equatorial KWs at theeastern boundary of the IO, followed by the radiation ofoff-equatorial RWs four times per year, is clearer at 550
m depth (Fig. 6d) than at the surface (Fig. 6b). This isbecause energy fluxes at the surface are forced by winds(contours in Fig. 6a, b) and represent the superpositionof forced and freely propagating waves. Energy input by
winds is more significant at 3° N than at the equator,
which may be confirmed by the distribution of energy
flux potential (contours in Fig. 4b). Returning to theissue of the SBoB downward energy flux that peaks in
November–December, we note that energy input by
winds peaks in February and June between 80 and 90° E(Fig. 6b). It is less likely that wind forcing is the primaryfactor to generate SBoB downward flux in November–
December. We note that the SBoB downward flux is induced by the boundary diffraction of mid-depth equatorial KW beam which peaks in August (Fig. 6c). As
explained in the previous paragraph, eastward energy
fluxes along the equator at 550 m depth lag behind thoseat the surface by 3 months.
This is further investigated using snapshots of energy
fluxes in the longitude–depth sections along the equator(Fig. 8, left) and along 3° N (Fig. 8, right). Zonal energy
flux at the equator propagates eastward at depths between 300 and 1000 m in August (colours in Fig. 8e).
This eastward peak in the western basin is followed by
the development of an eastward energy flux in the central basin at around 900 m depth in November (Fig. 8g).
These results indicate that the southwest monsoon over
the African coast directly causes the August maximumin eastward energy flux at 550 m depth in Fig. 6c. Southwesterly winds over the Arabian Sea drive a strong
northward East African Coastal Current along the continental coasts. The current crosses the equator and isaccelerated by winds in July, while the divergence of thesurface flow excites upwelling KWs along the equator.
The eastward transfer of KW energy is amplified by thesouthwest monsoon and appears not only at the sea

surface but also at mid-depth. The boundary-diffracted
RWs associated with these equatorial KWs dominate thedownward penetration of wave energy in the SBoB, withsignals reaching the greatest depths in November–December. Contours in Fig. 8 show the vertical componentof energy flux. At 3° N (right panels), energy flux is
downward all the year and represents an intense westward energy flux associated with off-equatorial RWs
around 75–90° E. Both westward and downward energy
fluxes reach greatest depths in November (Fig. 8h). The
above results highlight propagation and reflection ofoceanic waves and clearly determine their origins. Thus
the depth-dependent version of the AGC17 scheme isable to reveal the properties of waves at subsurface andprovides a better understanding of the processes governing the transfer of wave energy.

4 Conclusions
In order to study the effect of annual and semiannual
wind forcing on wave propagations in mid-depth in thetropical IO, we have estimated wave energy flux using
output from a linear model driven by climatologicalwind forcing. We have compared the performance of the
PF scheme to that of the depth-dependent version of the
AGC17 scheme summarized herein. The PF scheme represents meridionally alternating signals through thebasin at the surface and shows neither eastward KW energy flux along the equator nor the westward transfer ofenergy flux related to off-equatorial RWs. The AGC17
scheme provides a better treatment of energy flux ofboth equatorial KWs and RWs and smoothly links thetropical and subtropical regimes.
Eastward energy flux shows distinct semiannual variations along the equator at subsurface, peaking four times
per year, which indicate that energy given by surface
monsoonal winds reaches below the pycnocline. The
eastward flux at the surface on the equator is intense inboth February and May associated with upwelling anddownwelling KWs forced by semiannual winds, respectively. Eastward flux at 550 m depth at the equator is intense in both May and August and reaches down to1100 m depth in November. In off-equatorial regions,westward energy flux associated with RWs exhibits four
clearly defined peaks.
The SBoB is the region where wave energy penetrates
deepest. This downward energy flux is clearest at aroundLi et al. Progress in Earth and Planetary ScienceFig. 6 (See legend on next page.)(2021) 8:43

Page 13 of 19

Li et al. Progress in Earth and Planetary Science(2021) 8:43

Page 14 of 19(See figure on previous page.)
Fig. 6 Hovmöller zonal–time diagrams for the zonal component of energy flux (color shaded in units of W/m2) at a, b sea surface, c, d 550 m
depth, and e, f 1100 m depth. Left panels show the equatorial behavior, while right panels are for 3° N. Energy flux has been estimated using thedepth-dependent version of the Aiki, Greatbatch, and Claus Level-2 scheme, as given by (7b). Panels (a, b) include wind-forced wave energy
input (solid contours for positive values, dotted contours for negative values; contour interval = 0.002 W/m2) following (12). All quantities arefrom “Year 20” of the model experiments associated with first six fundamental baroclinic modes. In the left-hand panels, theoretical phase speeds
calculated for both Kelvin waves and equatorial long Rossby waves are shown as solid blue lines

Fig. 7 Energy input rate by wind forcing ( W/m2) near the African coast from “Year 20” of our climatological model experiments associated withthe sum of 6 baroclinic modes: a mid-February, b mid-May, c mid-August, and d mid-November. Spatially integrated values (35–60° E, 10° S–10°
N) of the wind input rate are shown in the title (109 W = GW)Li et al. Progress in Earth and Planetary Science(2021) 8:43

Page 15 of 19

Fig. 8 Zonal-vertical sections at the equator (left) and at 3° N (right) from “Year 20” of our climatological model experiments: a, b mid-February, c,d mid-May, e, f mid-August, and g, h mid-November. The zonal component of the energy flux (color shaded in units of W/m2), the verticalcomponent of the energy flux (solid contours for positive values, dotted contours for negative values; with an interval of 4 ×10−4 W/m2) and theenergy flux vector (arrows in units of W/m2) estimated from the depth-dependent version of the Aiki, Greatbatch, and Claus Level-2 scheme, as
given by (7b). When under 3 W/m2, the magnitude of energy flux is indicated either by arrow length or, for greater values, by heavy arrows of
uniform length

3–5° N, 90° E that we attributed to the downward flux of
RW energy. The energy input due to wind forcing at 3°
N peaks between 80 and 90° E in February and June andcannot account for downward energy flux inNovember–December. An alternative explanation is thatthe boundary diffraction of mid-depth equatorial KWs is
responsible for the downward penetration of energy in
the SBoB. The AGC17 scheme has shown that eastward

Li et al. Progress in Earth and Planetary Science(2021) 8:43

energy flux at 550 m depth at the equator peaks in August, which induces boundary-diffracted RWs at thisdepth at 3° N in November–December. The peak in August may in part be due to the maximum of eastward
energy flux in May at the surface on the equator. The
maximum mid-depth eastward flux lags behind that of
the surface eastward flux by 3 months, suggesting adownward propagation of wave energy.
A further important contribution is made by winds ofthe southwest monsoon in the Arabian Sea, which generates the divergence of energy at the surface near theequator along the African coastline and excites upwelling equatorial KWs in July. This process can have a direct impact on the formation of the peak in eastward
energy flux associated with KWs in August below 300 m
depth. These equatorial KW packets are diffracted poleward along the eastern boundary and become the sourceof off-equatorial RWs. The latter is responsible for the
downward flux of wave energy which reaches the greatest depth in November–December in the SBoB.
The advanced analytical and numerical techniques developed in this study can be extended to interannual
timescales. Horii et al. (2008) found that the subsurface
signals preceded the surface signals by 3 months during
the initial stage of the positive IO dipole event (Saji et al.
1999) in 2006, which hints a vertical propagation. The
subsurface upwelling phase of wind-forced KWs potentially affects the development of the IO dipole by alteringsea surface temperature (Chen et al. 2016). The propagation and reflection process of these wave energy is also
responsible for the significant intraseasonal oscillation in
the equatorial eastern IO (Chen et al. 2015b) and within
the Bay of Bengal (Chen et al. 2018). Future work in thisdirection will enhance our knowledge of the highly influential interactions and exchanges between the tropicaland extratropical regions (particularly during active
phases of the IO dipole) and will therefore contribute toa better understanding of climate variation.

5 Appendix
5.1 Linear Ocean Model (LOM)

We simulated climatological variability in the IO using alinear ocean model (LOM). The model has been discretized in a spherical coordinate system with a grid spacing of (1/4)° in both zonal and meridional directions.
The model domain extends from 35 to 120° E and from
25° S to 20° N and adopts realistic coastlines except for
the lack of connection to the Southern Ocean and thePacific Ocean. The LOM is written for the developmentof zonal velocity u(n)(x, y, t), meridional velocity v(n)(x, y,t), and geopotential Φ(n)(x, y, t) associated with the generic n-th baroclinic mode. For each baroclinic modethese are related throughPage 16 of 19

ðn Þ

τx
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ρ0 hmix H bottom
þ SGSx ;ð9aÞτypﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ρ0 hmix H bottom
þ SGSy ;ð9bÞ

ut −f vðnÞ þ ΦðxnÞ ¼ αðnÞðn Þ

vt þ f uðnÞ þ ΦðynÞ ¼ αðnÞ 2 ðnÞΦt þ cðnÞ
uðxnÞ þ vðynÞ ¼ 0;

ð9cÞ

where horizontal and temporal partial derivatives are denoted by subscripts x, y, and t respectively. f denotes the
Coriolis coefficient. ρ0 = 1027 kg/m3 is the reference
value of water density. The vector (τx, τy) represents
wind stress. The mode-dependent parameters α(n) andc(n) have been calculated in the present study using theclimatological annual mean stratification of the tropicalIO (explained later in this section). The definition of α(n)is described in detail in the Appendix of LA20. The
values of α(n) and c(n) used in the present study are given
in Table 1. The present study defines the mixed layer
depth as the depth where temperature is lower than thesurface by 0.2°C and calculated it using Argo float data.
From the resulting annual mean of mixed layer depth in
the tropical IO, we set hmix = 35 m. Hbottom= 5500 m
represents the reference depth of the bottom of theocean in the tropical IO. The subgrid-scale terms, which
are labelled as SGS in (9a) and (9b), indicate the effect oflateral eddy viscosity applying the Smagorinsky et al.
(1965) scheme with a nondimensional coefficient of 0.1
after squaring.
Reconstruction of vertical profile based on the outputof LOM is done by using the eigen vectors, as in previous studies. While LA20 executed a series of numericalexperiments focusing on the first three baroclinic modes
which assumed no motion at the bottom, the present
study has investigated the first six baroclinic modes
under the assumption of a free-slip bottom boundary(Chelton and Schlax 1996; see below). The first six
modes explain about 80% of baroclinic kinetic energy
fraction in the equatorial regions of a global ocean reanalysis, as shown by Fig. 1 in Toyoda et al. (2021). The
Sturm-Liouville equation for the baroclinic normalmode decomposition associated with a stratified ocean
may be written asF ðznÞ =N 2


z

¼

−1
ðcðnÞ Þ

2 F

ðnÞ

;

ð10Þ

where F(n) = F(n)(z) and c(n)are the eigenfunction andnon-rotating gravity wave speed, respectively, of the n-th
baroclinic mode. The symbol N = N(z) in (10) is thebuoyancy frequency which we have given based on thevertical profiles of climatological annual mean salinity

Li et al. Progress in Earth and Planetary Science(2021) 8:43

and temperature in the tropical IO (20° S–20° N, 40–
110° E) from the World Ocean Atlas. We have solved(10) using a Neumann boundary condition, F ðnÞ
z ¼ 0, atboth the surface and bottom of the ocean and obtainedc(n)=2.63 m/s, 1.57 m/s, 0.956 m/s, 0.706 m/s, 0.557 m/s,
and 0.467 m/s for the first six baroclinic modes. The
vertical profiles of the eigenfunctions are shown in FigR0
2
ure 9 in Appendix where −H bottom ð F ðnÞ Þ dz ¼ H bottom andbottom depth Hbottom= 5500 m. The parameter α(n) ofwind forcing partition in (9a)–(9b) has been determinedin the same manner as that in LA20.
McCreary (1984) demonstrated that the vertical propagation related to equatorial waves can be described by thesum of baroclinic modes. The present study has reconstructed the vertical structure of the tropical IO by using
the results of model experiment associated with fundamental six baroclinic modes. Quantities in the depth

Page 17 of 19

coordinate, such as geopotential Φ(x, y, z, t), zonal velocityu(x, y, z, t), meridional velocity v(x, y, z, t), and vertical velocity w(x, y, z, t), are obtained asXΦðx; y; z; t Þ ¼
ΦðnÞ ðx; y; t Þ F ðnÞ ðzÞ;ð11aÞnuðx; y; z; t Þ ¼

XuðnÞ ðx; y; t Þ F ðnÞ ðzÞ;ð11bÞvðnÞ ðx; y; t Þ F ðnÞ ðzÞ;

ð11cÞ

n

vðx; y; z; t Þ ¼

Xnwðx; y; z; t Þ ¼

X −1  ∂ðnÞ
2 ∂t Φ ðx; y; t Þ
N
n

∂ ðnÞ
F ðzÞ ;

∂z



ð11dÞ

and have been applied to each of the 100 snapshots obtained from the final year of model output.

Fig. 9 The vertical profiles of eigenfunctions F(n) associated with fundamental six baroclinic modes in the tropical Indian Ocean, calculated bysolving the Sturm-Liouville problem. Details are in the Appendix

Li et al. Progress in Earth and Planetary Science(2021) 8:43

Wave energy Eq. (3a) should include wind forcing anddissipation terms. These terms yield the energy flux potential defined in (8a)–(8b) in a climatological equilibrium state of wave energy. The depth-integral of thewind forcing term is written as
uðx; y; 0; t Þτ x þ vðx; y; 0; t Þτ yrﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
i
H bottom X ðnÞ h ðnÞ
x
ðn Þ
y
¼
α
u
ð
x;
y;
t
Þτ
þ
v
ð
x;
y;
t
Þτ
n
hmix
ð12Þ
which is referred to as wind input (see contours in Fig.
6a, b and color shading in Fig. 7).
Abbreviations
AGC17: Aiki et al. 2017; ECMWF ERA-40: European Centre for Medium-Range
Weather Forecasts Reanalysis; IGWs: Inertia-gravity waves; IO: Indian Ocean;
KWs: Kelvin waves; LA20: Li and Aiki 2020; LOM: Linear ocean model;
OGCM: Ocean General Circulation Model; OFES: Ocean General Circulation
Model For the Earth Simulator; PF: Pressure flux; RGWs: Mixed Rossby-gravity
waves; RWs: Rossby waves; SGS: Subgrid scale; SBoB: Southern Bay of Bengal

6 Supplementary Information
The online version contains supplementary material available at https://doi.
org/10.1186/s40645-021-00432-9.
Additional file 1. The extension of the formulation of Aiki et al. (2017).
Authors’ contributions
ZL performed the numerical experiments and carried out analysis. HA
proposed the topic and conceived the study. MN and TO helped improve
the integrity and construction of the manuscript. The authors read andapproved the final manuscript.
Funding
This study was supported by Japan Society for the Promotion of Science(JSPS) KAKENHI Grant Numbers 18H03738 and 15H02129.
Availability of data and materials
The set of model outputs obtained during the last year of each of the six
experiments performed in the present study will be archived at both https://www.diasjp.net/en/ and http://cidas.isee.nagoya-u.ac.jp/databases/index.
shtml.en upon publication of this manuscript. A number of datasets wereretrieved from the Asia-Pacific Data Research Center at the University of Hawai’i (http://apdrc.soest.hawaii.edu). These include temperature and salinity
data from the World Ocean Atlas 2013 for the calculation of gravity wave
phase speed, mixed-layer depth data from the Argo products, wind velocity
data from the European Centre for Medium-Range Weather Forecasts Reanalysis (ECMWF ERA-40) for estimation of climatological wind forcing, andoutput of velocity and zonal stress from the Ocean General Circulation
Model For the Earth Simulator (OGCM OFES).

Declarations
Competing interests
The authors declare that they have no competing interest.
Author details
1
Graduate School of Environmental Studies, Nagoya University, Furo-cho,
Chikusa-ku, Nagoya 464-8601, Japan. 2Institute for Space-Earth Environmental
Research, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan.
3
Application Laboratory, Japan Agency for Marine-Earth Science andTechnology, Kanazawa-ku, Yokohama 236-0001, Japan. 4Research Institute for
Global Change, Japan Agency for Marine-Earth Science and Technology,
Yokosuka 237-0061, Japan.

Page 18 of 19

Received: 24 December 2020 Accepted: 1 June 2021

References
Aiki H, Greatbatch RJ, Claus M (2017) Towards a seamlessly diagnosable
expression for the energy flux associated with both equatorial and midlatitude waves. Prog Earth Planet Sci 4(1):11. https://doi.org/10.1186/s40645017-0121-1
Brandt P, Eden C (2005) Annual cycle and interannual variability of the middepth tropical Atlantic Ocean. Deep Sea Res 52(2):199–219. https://doi.org/1
0.1016/j.dsr.2004.03.011
Chelton DB, Schlax MG (1996) Global observations of oceanic Rossby waves.
Science 272(5259):234–238. https://doi.org/10.1126/science.272.5259.234
https://www.jstor.org/stable/2889630
Chen G, Han W, Li Y, McPhaden MJ, Chen J, Wang W, Wang D (2017) Strong
intraseasonal variability of meridional currents near 5° N in the eastern Indian
Ocean: characteristics and causes. J Phys Oceanogr 47(5):979–998. https://doi.
org/10.1175/JPO-D-16-0250.1
Chen G, Han W, Li Y, Wang D (2016) Interannual variability of equatorial Eastern
Indian Ocean upwelling: local versus remote forcing. J Phys Oceanogr 46(3):789–807. https://doi.org/10.1175/JPO-D-15-0117.1
Chen G, Han W, Li Y, Wang D, McPhaden MJ (2015a) Seasonal-to-interannual
time-scale dynamics of the equatorial undercurrent in the Indian Ocean. J
Phys Oceanogr 45(6):1532–1553. https://doi.org/10.1175/JPO-D-14-0225.1
Chen G, Han W, Li Y, Wang D, Shinoda T (2015b) Intraseasonal variability of
upwelling in the equatorial Eastern Indian Ocean. J Geophys Res 120(11):7598–7615. https://doi.org/10.1002/2015JC011223
Chen G, Han W, Li Y, Yao J, Wang D (2019) Intraseasonal variability of theequatorial undercurrent in the Indian Ocean. J Phys Oceanogr 49(1):85–101.
https://doi.org/10.1175/JPO-D-18-0151.1
Chen G, Han W, Zhang X, Liang L, Xue H, Huang K, He Y, Li, Wang D (2020)
Determination of spatiotemporal variability of the Indian equatorialintermediate current. J Phys Oceanogr 50:3095–3108. https://doi.org/10.1175/
JPO-D-20-0042.1
Chen G, Li Y, Xie Q, Wang D (2018) Origins of Eddy kinetic energy in the Bay of
Bengal. J Geophys Res 123(3):2097–2115. https://doi.org/10.1002/2017JC013455
Dewitte B, Reverdin G (2000) Vertically propagating annual and interannual
variability in an OGCM simulation of the tropical Pacific Ocean in 1985–94. J
Phys Oceanogr 30:1562–1158. https://doi.org/10.1175/1520-0485(2000)030<
1562:VPAAIV>2.0.CO;2
Ding H, Keenlyside NS, Latif M (2009) Seasonal cycle in the upper equatorial
Atlantic Ocean. J Geophys Res 114(C9):C09016. https://doi.org/10.1029/2
009JC005418
Fu LL (2007) Intraseasonal variability of the equatorial Indian Ocean observed
from seasurface height, wind, and temperature data. J Phys Oceanogr 37(2):
188–202. https://doi.org/10.1175/JPO3006.1
Gill AE (1982) Atmosphere-ocean dynamics. Academic, London, p 662
Han W (2005) Origins and dynamics of the 90-day and 30–60-day variations in
the equatorial Indian Ocean. J Phys Oceanogr 35(5):708–728. https://doi.
org/10.1175/JPO2725.1
Hayes SP, Milburn HB (1979) On the vertical structure of velocity in the easternequatorial Pacific. J Phys Oceanogr 10:633–635. https://doi.org/10.1175/15200485(1980)010<0633:OTVSOV>2.0.CO;2
Horii T, Hase H, Ueki I, Masumoto Y (2008) Oceanic precondition and evolutionof the 2006 Indian Ocean dipole. Geophys Res Lett 35(3):L03607. https://doi.
org/10.1029/2007GL032464
Huang K, Han W, Wang D, Wang W, Xie Q, Chen J, Chen G (2018a) Features of
the Equatorial Intermediate Current associated with basin resonance in the
Indian Ocean. J Phys Oceanogr 48(6):1333–1347. https://doi.org/10.1175/JPOD-17-0238.1
Huang K, McPhaden MJ, Wang D, Wang W, Xie Q, Chen J, Shu Y, Wang Q, Li J,
Yao J (2018b) Vertical propagation of mid-depth zonal currents associatedwith surface wind forcing in the equatorial Indian Ocean. J Geophys Res
123(10):7290–7307. https://doi.org/10.1029/2018JC013977
Huang K, Wang D, Han W, Feng M, Chen G, Wang W, Chen J, Li J (2019)
Semiannual variability of middepth zonal currents along 5° N in the eastern
Indian Ocean: characteristics and causes. J Phys Oceanogr 49(10):2715–2729.
https://doi.org/10.1175/JPO-D-19-0089.1
Ishizaki H, Nakano H, Nakano T, Shikama N (2014) Evidence of equatorial Rossby
wave propagation obtained by deep mooring observations in the westernLi et al. Progress in Earth and Planetary Science(2021) 8:43

Pacific Ocean. J Oceanogr 70(6):463–488. https://doi.org/10.1007/s10872-0140247-3
Johnson GC (2011) Deep signatures of southern tropical Indian Ocean annual
Rossby waves. J Phys Oceanogr 41(10):1958–1964. https://doi.org/10.1175/
JPO-D-11-029.1
Kessler WS, McCreary JP (1993) The annual wind-driven Rossby-wave in thesubthermocline equatorial Pacific. J Phys Oceanogr 23:1192–1207. https://doi.
org/10.1175/1520-0485(1993)0231192:TAWDRW2.0.CO;2
Large WG, Pond S (1981) Open ocean momentum flux measurements inmoderate to strong winds. J Phys Oceanogr 11:324–336. https://doi.org/10.11
75/1520-0485(1981)011<0324:OOMFMI>2.0.CO;2
Li Z, Aiki H (2020) The life cycle of annual waves in the Indian Ocean as
identified by seamless diagnosis of the energy flux. Geophys Res Lett 47(2):e2019GL085670. https://doi.org/10.1029/2019GL085670
Longuet-Higgins HC (1964) On group velocity and energy flux in planetary wave
motion. Deep Sea Res 11:35–42. https://doi.org/10.1016/0011-7471(64)91080-0
Lucas R, Firing E (1985) The annual Rossby wave in the central equatorial Pacific
Ocean. J Phys Oceanogr 15(1):55–67. https://doi.org/10.1175/1520-0485(1985)0150055:TARWIT2.0.CO;2
Luyten JR, Roemmich DH (1982) Equatorial currents at semiannual period in the
Indian-Ocean. J Phys Oceanogr 12:406–413. https://doi.org/10.1175/1520-04
85(1982)012<0406:ECASAP>2.0.CO;2
Luyten JR, Swallow JC (1976) Equatorial undercurrents. Deep Sea Res 23:999–
1001. https://doi.org/10.1016/0011-7471(76)90830-5
Mark RJ, Huang B (2004) The Rossby wave as a key mechanism of Indian Ocean
climate variability. Deep Sea Res Part I 51(12):2123–2136. https://doi.org/10.1
016/j.dsr.2004.06.005
Masumoto Y et al (2004) A fifty-year eddy-resolving simulation of the world
ocean: preliminary outcomes of OFES (OGCM for the Earth Simulator). J Earth
Simulator 1:35–56
Matsuno T (1966) Quasi-geostrophic motions in the equatorial area. J Meteo Soc
Japan 44:25–43. https://doi.org/10.2151/jmsj1965.44.1_25
McCreary JP (1984) Equatorial beams. J Mar Res 42(2):395–430. https://doi.org/1
0.1357/002224084788502792
McPhaden MJ (1982) Variability in the central equatorial Indian-Ocean. I: ocean
dynamics. J Mar Res 40:157–176
McPhaden MJ, Ripa P (1990) Wave-mean flow interactions in the equatorial
ocean. Annu Rev Fluid Mech 20:167–205. https://doi.org/10.1146/annurev.
fl.22.010190.001123
Miyama T, McCreary JP, Sengupta D, Senan R (2006) Dynamics of biweekly
oscillations in the equatorial Indian Ocean. J Phys Oceanogr 36(5):827–846.
https://doi.org/10.1175/JPO2897.1
Nagura M (2018) Annual Rossby waves below the pycnocline in the Indian
Ocean. J Geophys Res 123:9405–9415. https://doi.org/10.1029/2018jc014362
Nagura M, McPhaden MJ (2010a) Wyrtki jet dynamics: seasonal variability. J
Geophys Res 115:C07009. https://doi.org/10.1029/2009JC005922
Nagura M, McPhaden MJ (2010b) Dynamics of zonal current variations associatedwith the Indian Ocean dipole. J Geophys Res 115(C11):C11026. https://doi.
org/10.1029/2010JC006423
Nagura M, McPhaden MJ (2016) Zonal propagation of near-surface zonal currents
in relation to surface wind forcing in the equatorial Indian Ocean. J Phys
Oceanogr 46(12):3623–3638. https://doi.org/10.1175/jpo-d-16-0157.1
Ogata T, Aiki H (2019) The pathway of intraseasonal wave energy in the tropicalIndian Ocean as identified by a seamless diagnostic scheme. SOLA 15:262–
267. https://doi.org/10.2151/sola.2019-047
Pujiana K, McPhaden MJ (2020) Intraseasonal Kelvin waves in the equatorial
Indian Ocean and their propagation into the Indonesian seas. J Geophys Res
125(5):e2019JC015839. https://doi.org/10.1029/2019JC015839
Reppin J, Schott FA, Quadfasel D (1999) Equatorial currents and transports in theupper central Indian Ocean: annual cycle and interannual variability. J
Geophys Res Oceans 104(C7):15495–15514. https://doi.org/10.1029/1
999JC900093
Rothstein LM, Moore DW, McCreary JP (1985) Interior reflections of a periodicallyforced equatorial Kelvin wave. J Phys Oceanogr 15(7):985–996. https://doi.
org/10.1175/1520-0485(1985)015<0985:IROAPF>2.0.CO;2
Saji NH, Goswami BN, Vinayachandran PN, Yamagata T (1999) A dipole mode in
the tropical Indian Ocean. Nature 401:360–363. https://doi.org/10.1038/43854
Schott FA, McCreary JP (2001) The monsoon circulation of the Indian Ocean.
Prog Oceanogr 51(1):1–123

Page 19 of 19

Smagorinsky J, Manabe S, Holloway JL (1965) Numerical results from a nine-level
general circulation model of the atmosphere. Mon Weather Rev 93:727–768.
https://doi.org/10.1175/1520-0493(1965)0930727:NRFANL2.3.CO;2
Song Q, Aiki H (2020) The climatological horizontal pattern of energy flux in the tropicalAtlantic as identified by a unified diagnosis for Rossby and Kelvin waves. J Geophys
Res 125:e2019JC015407. https://doi.org/10.1029/2019JC015407
Thierry V, Mercier H, Treguier AM (2006) Seasonal fluctuations in the deep central
equatorial Atlantic Ocean: Aa data-model comparison. Ocean Dyn 56(5–6):
581–593. https://doi.org/10.1007/s10236-005-0045-y
Toyoda T, Nakano H, Aiki H, Ogata T, Fukutomi Y, Kanno Y, Urakawa LS,
Sakamoto K, Yamanaka G, Nagura M (2021) Energy flow diagnosis of ENSO
from an ocean reanalysis. J Clim:1–86. https://doi.org/10.1175/JCLI-D-200704.1
Wyrtki K (1973) An equatorial jet in the Indian Ocean. Science 181(4096):262–264.
https://doi.org/10.1126/science.181.4096.262
Yanai M, Maruyama T (1966) Stratospheric wave disturbances propagating overthe equatorial Pacific. J Met Soc Japan 44(5):291–294. https://doi.org/10.2151/
jmsj1965.44.5_291
Yu X, McPhaden MJ (1999) Seasonal variability in the equatorial Pacific. J Phys
Oceanogr 29:925–947. https://doi.org/10.1175/1520-0485(1999)0290925:
SVITEP2.0.CO;2
Yuan DL, Han WQ (2006) Roles of equatorial waves and western boundaryreflection in the seasonal circulation of the equatorial Indian Ocean. J Phys
Oceanogr 36(5):930–944. https://doi.org/10.1175/JPO2905.1
Zanowski H, Johnson GC (2019) Semiannual variations in 1,000-dbar equatorial
Indian Ocean velocity and isotherm displacements from Argo data. J
Geophys Res Oceans 124(12):9507–9516. https://doi.org/10.1029/2019JC01
5342

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.