The earliest research of an equatorial wave system was done by Matsuno (1966). As the Coriolis parameter $\,f$ is 0 and its derivative $\beta = {{{\rm d}\,f} / {{\rm d}y}}$ reaches a maximum at the equator, fluctuations in equatorial waters become the most distinctive motion in marine dynamics. In the equatorial ocean, there are Kelvin waves travelling eastward and Rossby waves travelling westward, as well as inertial gravity waves of high frequency and mixed Rossby inertial-gravity waves. These waves constitute the Matsuno wave system, also known as the equatorial wave system. The adjustment of Kelvin and Rossby waves in the equatorial wave system will produce the oscillation of the interannual time scale and can explain El Niño-Southern Oscillation (ENSO) circulation. Thus, the importance of the fluctuations in the equatorial waters has gradually been recognized.

Waves which propagate in parallel with the shoreline but become obvious only in the vicinity of the shoreline were first studied by Lord Kelvin and are thus named after their discoverer. Obviously the constraint of shoreline is an important factor. Matsuno (1966) found that similar waves also exist even without the constraint of shoreline when the Coriolis parameter changes with space, a condition in which a wave form is only obvious close to the equator and there is only one eastward moving wave. Here, the equator plays the role of shoreline. Hirst (1989) and Yamagata and Masumoto (1989) theoretically proved that disturbances propagating to the east exist in an ocean-atmosphere coupled model, which they believe is induced by eastward moving equatorial Kelvin waves. We also found that subsurface ocean temperature anomalies (SOTA) in the tropical Pacific propagate from west to east along the thermocline of the equator in years of both higher-than-average plum rain precipitation and lower-than-average plum rain precipitation in the middle and lower reaches of the Changjiang (Yangtze) River (Lu and Zhang, 2009), which suggests that this may be related with eastward moving Kelvin waves in the equatorial Pacific.

Nevertheless, as there are coastlines in the eastern and western tropical Pacific, Kelvin waves will inevitably be reflected by the eastern boundaries of the ocean as they travel eastward, disrupting their mode of propagation, which has been used to explain the ENSO (Chao, 1993). This indicates that travelling Kelvin waves in an equatorial ocean cannot satisfy lateral boundary conditions from the perspective of mathematical physics, nor is it the solution of that bounded ocean. All of the three major tropical oceans, i.e., the Pacific, the Indian Ocean and the Atlantic Ocean, have boundaries. Thus it is of great significance to seek the solution of equatorial oceans with boundaries in the equator wave system. Does the solution exist? If it does, what are its spatial structure and vibrational frequencies? How is it connected with climate change? These are key issues of physical oceanography and climate science. Existing research in this regard, however, is lacking, and the above questions remain largely unanswered. Our earlier analyses (Zhang, 2006; Zhang and He, 2005; Lu et al., 2014) of complex EOF (empirical orthogonal function) for the flow field anomalies of the tropical Pacific and the tropical Indian Ocean show that the spatial field of anomalous flow fields, in the form of zonal flows, is trapped near the equator and decays rapidly away from it. This suggests that these equatorial flow field anomalies are very likely to be a wave packet solution satisfying the requirement of the wave system of bounded equatorial oceans. This wave packet solution may play an important role in the flow field anomalies in the equatorial ocean. An ideal ocean model is used in this paper to answer these questions and prove the above hypothesis that the mode of the complex EOF for the flow field anomalies is the wave packet solution of the bounded equatorial oceans. Since the upper ocean above the thermocline can be seen as positive, linearized shallow disturbance equations (positive) with approximation in the equatorial $\beta $ plane is used in this paper. By introducing reduced gravity (Zhang, 1995) and taking the lateral boundary conditions of the east and west coasts in an ideal ocean into account, the paper seeks the analytical solution of wave packet perturbation (anomaly) of the bounded ocean in a climatic scale, which is necessary and highly significant for understanding its shape and characteristics and for explaining the nature of the complex EOF analysis results of the flow field anomalies made earlier on the tropical Pacific and Indian Oceans.

Suppose the ocean is divided homogeneously in the vertical direction into upper and lower layers, which is bounded by thermocline. The densities of the two layers are constant ${\rho _1}$ and ${\rho _2}$ respectively after the introduction of reduced gravity (Zhang and He, 2005). Suppose the lower ocean is still, then the governing equation of the upper ocean may be considered as satisfying linearized shallow disturbance equations (positive) with approximation in equatorial the $\beta $ plane.

where $\bar H$ is the average depth of the upper layer and is a constant; $\eta $ represents a small disturbance near $\bar H$; $\hat g = \left( {{\rho _2} - {\rho _1}} \right)g{\rm{/}}{\rho _2}$, is reduced gravity and $g$ is the gravity. In the non-dimensionalization of the above equations, “′” represents the non-dimensional variables. $u' = u/\bar U$, $v' = v/\bar V$ and $\delta \eta ' = \Delta \eta /\Delta \bar \eta $, are introduced for dependent variables, with “-” representing the corresponding scale. $x' = x/{L_x}$, $y' = y/{L_y}$ and $t' = t/\bar t$ are introduced for independent variables. Here ${L_x}$, ${L_y}$ and $\bar t$ represent the scale of corresponding variables. Then the corresponding dimensionless equations (Xuan et al., 2014; Wu, 2002) can be obtained as follows:

where $\varepsilon = {L_y}/{L_x}$ and it is equivalent to long wave approximation when $\varepsilon = 0$. The symbol “′” representing non-dimensional variables in Eq. (2) is omitted for the convenience of writing and so is in the following text unless otherwise indicated. Suppose the wave solution is

where Φ represents the structure of η in the y (meridional) direction.

Then

The solution of Eq. (4) is as follows. First considering the simplest Kelvin wave, we set $V \equiv 0$. Thus, according to the existing literature (Xuan et al., 2014; Wu, 2002), the dispersion relation of $\omega $ and $k$ can be obtained as

It can be seen that the wave propagates eastward, and the solution of the wave is

where $\tilde U$ is a random integral constant; and $\tilde \varPhi = \tilde U\omega /k$, is the amplitude of Kelvin wave.

Next, we take general circumstances into account. According to literature (Xuan et al., 2014; Wu, 2002), similarly, by eliminating $U$ and $\varPhi $ from Eq. (4), we can get the second-order linear equation with variable coefficients as

This is a Hermite equation, which means its solution with boundaries, i.e., the solution with physical significance, contains Hermite polynomials. This solution must satisfy the following condition (Wu, 2002; Zhang and Zhang, 2006):

The above formula is also the dispersion relation of $\omega $ and $k$. The solution meeting the requirement of this dispersion relation is

were $\hat V\,$ is also a random integral constant and ${{\rm{H}}_n}(y)$ are Hermite polynomials. The first six polynomials are ${{\rm{H}}_0}(y) = 1$, ${{\rm{H}}_1}(y) = 2y$, ${{\rm{H}}_2}(y) = 4{y^2} - 2$, ${{\rm{H}}_3}(y) = {2^3}{y^3} - 12y$, ${{\rm{H}}_4}(y) = {2^4}{y^4} - $$ 48{y^2} + 12$, ${{\rm{H}}_5}(y) = {2^5}{y^5} - 160{y^3} + 120y$, respectively. By eliminating $\varPhi $ from Eq. (4), we can get$\left({{\omega ^2} - {k^2}} \right)U + \omega yV - k{\rm{d}}V/{\rm{d}}y = 0$. By substituting Eq. (9a) in this equation and using derivative formula of Hermite polynomials ${\rm{d}}{{\rm{H}}_n}/{\rm{d}}y = 2n{{\rm{H}}_{n - 1}}$ and recursive formula ${{\rm{H}}_{n + 1}} = 2y{{\rm{H}}_n} - 2n{{\rm{H}}_{n - 1}}$, we can get the following equation:

where $\hat U = \hat V$, is an arbitrary constant. After $V$ and $U$ are solved, the following equation can be obtained according to Eq. (4a):

where $\hat \varPhi = \hat U = \hat V$, is an arbitrary constant.

It can be seen that factor ${{\rm{e}}^{ - \scriptsize\displaystyle\frac{1}{2}{y^2}}}$ exists in the expression of both Eq. (6) and Eq. (10). This shows that it decays faster than exponential decay with increasing distance from the equator, which means the disturbance or anomalies may be trapped near the equator. The following is a brief discussion of the dispersion Eq. (8), which is a cubic equation for frequency $\omega $. Three real $\omega $ can be obtained with a given wavenumber $k$ (or wavelength) and an integer $n \text{≥} 1$. They correspond to a couple of inertial gravitational waves and westward Rossby waves, respectively. The former two are fast waves and the latter is a slow wave. The frequencies of the former two and the latter are absolutely distinguishable. When $n = 0$, one solution of Eq. (8) is $\omega = - k$, which corresponds to Rossby-inertial gravitational waves. Besides, there are two inertial gravitational waves propagating eastward and westward, respectively. $n$ must be a nonnegative integer in Eq. (8). However, if $n = - 1$, the solution of Eq. (8) is $\omega = k$, which is Eq. (5) above, namely, the dispersion relation of Kelvin wave. Thus, the dispersion relation of Kelvin wave is also included in Eq. (8).

$U, \;V\;$ and $\varPhi $ in Eqs (6) and (9) are substituted into Eq. (3) and the real part is taken to obtain the solution of the disturbances of waves propagating along the direction of x in the unbounded ocean, i.e., the solution of anomalies. In practice, however, the Pacific, Indian and Atlantic Oceans are invariably bounded at the equator by eastern and western coasts, while none of the above wave solutions satisfy the boundary conditions. Suppose the ocean is bounded by the eastern and western coasts in the north-south direction and its length is $L = {L_x}$, then the boundary conditions can be represented as

Then the wave solution satisfying the boundary Eq. (10) can be obtained. As the purpose of this paper is to discuss the interannual and decadal variabilities of abnormal fluctuations, the gravity inertia wave which propagates faster than the Kelvin wave is ruled out. Thus, the Kelvin wave is the only wave propagating eastward. It can be seen from Eq. (6a) that $U$ reaches the maximum at the equator. After ruling out the gravity inertia wave propagating westward, the mixed Rossby inertial-gravity wave and the Rossby wave are the only westward travelling waves. For the former, as can be seen from Eq. (9b), ${{\rm{H}}_1}(y) = 2y$ when $n = 0$, and $U$ is 0 at the equator where $y = 0$, obviously making it impossible be superposed on Kelvin waves to satisfy the boundary Eq. (10). Therefore, Kelvin waves superposed on Rossby waves must be considered to see if it satisfies the boundary Eq. (10). To explore the interannual and decadal variabilities of the abnormal fluctuations as mentioned above, the lowest frequency modes should be considered first. Now examine the structure of Rossby waves with the lowest frequency in the direction of $y$. Here, we might take $n = 1$. In addition, as fluctuations are trapped near the equator, it is only necessary to study the situations near the equator. If $n = 1$, noticeably ${{\rm{H}}_1} = 2y$ and ${{\rm{H}}_2} = 4{y^2} - 2$, it can thus be derived that:

The structure of Eq. (11a) is the same with that of Eq. (6a) in the direction of $y$. When $y \approx 0$, it can be written as

As Eq. (11b) shows, $V \approx 0$. It can also be found that the structure of $\varPhi $ is the same with Eq. (6b) in the direction of $y$. Thus we can consider the superposition of eastward travelling Kelvin waves and westward travelling Rossby waves when $n = 1$. Here, take the frequency ${\omega _0}$ of Kelvin wave, wave number ${k_1}$ and integral constant $\tilde U$. Then the solution of eastward travelling waves in an unbounded ocean is (substitute Eq. (6a) into Eq. (3) and take the real part of it)

Take the frequency $ - {\omega _0}$ of Rossby wave, when $n = 1$, wave numbers ${k_2}$ and integral constant $\hat U$. Then the solution of westward travelling waves in an unbounded ocean is (substitute Eq. (11a) into Eq. (3) and take the real part of it)

Equations (12a) and (12b) are both solutions of Eq. (2). The equations are homogeneous linearity equations and satisfy a superposition principle. Therefore, the solution ${u_0} = {u_1} - {u_2}$ after superposing two waves is also the solution of Eq. (2). Now, take integral constants ${\tilde U_1} = {U_0}$ and ${\hat U_2} = {U_0}\left({\omega _0^2 - k_2^2} \right)/\left({2{k_2}} \right)$. Then

Equation (5) shows that ${k_1} = {\omega _0}$$\cdot \varepsilon \left[ {{{\left({ - {\omega _0}} \right)}^2} - k_2^2} \right] - {k_2}/\left({ - {\omega _0}} \right) = 3$ can be gotten from Eq. (8). For convenience’s sake while not sacrificing generality, we take a long wave approximation, i.e., take $\varepsilon = 0$. Here, the dispersion relation of Rossby wave is ${k_2} = 3{\omega _0}$ and ${k_2} = 3{k_1}$ can be derived. Substituting this into Eq. (13), we can get

In the above equation, the integral constant is still denoted by ${U_0}$. The equation is a wave packet and satisfies Eq. (2). Here, we have

Equation (14a) can now be examined to see whether it satisfies the boundary Eq. (10). Obviously, it satisfies the conditions that $x = 0, \quad {u_0} = 0$. In addition, when $x = 1$, if it satisfies the relation ${k_1} = l{\rm{\text{π} }}/2$, here $l$ being an integer and a mode number, then it also satisfies the boundary condition $x = 1, \; {u_0} = 0$. We can therefore see that wave number ${k_1}$ of Kelvin wave cannot be arbitrary, and can only take discrete values satisfying the relationship ${k_1} = l{\rm{\text{π} }}/2$ (for $l$ is a positive integer). Once mode number $l$ is determined, then wave number ${k_1}$ is fixed, and according to Eq. (5), the frequency is also determined. Therefore, the absolute value of the frequency of the corresponding Rossby wave equals the frequency of the Kelvin wave and the wave number of the Rossby wave is three times of that of the Kelvin wave.

Finally, in order to discuss and compare with actual ocean situations, the above dimensionless solution is converted to a dimensional one. The following dimensional quantity is

where $L$ is the length of the bounded ocean. The wave number of the Kelvin wave can only be${k_1} = l{\rm{\text{π} }}/L$. Thus the solution satisfying Eq. (1) and boundary Eq. (15) at the same time near the equator is

Note that all variables and parameters in the above Eq. (15) are dimensionless. It can be seen from Eq. (16) that the solution is a wave package and its envelope is a sine wave unchanged with time, of which the wave number is $l{\rm{\text{π} }}/L$ and the wave length is $2L/l$. The carrier frequency is a sine wave propagating westward and its wave number is $l{\rm{\text{π} }}/(2L)$, the wave length is $4L/l$, the wave speed is ${c_0} = \sqrt {\hat g\bar H} $, the oscillation frequency is $\sigma = l{\rm{\text{π} }}{c_0}/(2L) = $$l{\rm{\text{π} }}\sqrt {\hat g\bar H} /(2L) $ and the oscillation period is $T = 4L/(l{c_0}) = 4L/$$\left( l\sqrt {\hat g\bar H} \right)$. The frequency $\sigma $ and the period $T$ is also the respective frequency and period of the bounded equatorial ocean wave packet.

In the following calculations, take the standard depth of the upper ocean $\bar H = {200}\;{\rm{m}} $, take $\left( {{\rho _2} - {\rho _1}} \right){\rm{/}}{\rho _2}$ coefficient $1.704 \times {10^{ - 5}}$ and acceleration of gravity $g = {9.8}\;{{\rm{m/}}{{\rm{s}}^{\rm{2}}}} $. In this case $\hat g = {{1}}{{.67}} \times {{1}}{{{0}}^{ - 4}}$ and ${c_0} = \sqrt {\hat g\bar H} = {{{0}}{{.183}}}\;{{\rm{m/}}{{\rm{s}}}} $. Then take $\beta = {2.289 \times {{10}^{ - 11}}}\;{{{\rm{m}}^{{{ - 1}}}} \cdot {{\rm{s}}^{{{ - 1}}}}} $ and ${U_0} = {0.02}\;{{\rm{m/}}{{\rm{s}}}} $. Thus $\beta /(2{c_0}) = {{ 6}}{{.26}} \times {10^{ - 11}}$. Note that these environmental parameters are not related to the length of the bounded ocean $L$ or mode number $l$.

This paper uses linearized shallow water perturbation equations with approximation in equatorial β plane, and introduces a reduced gravity to obtain the analytical solution of the wave packet anomalies on the upper bounded equatorial ocean. The calculation results are given and compared with flow field anomalies in the tropical Pacific Ocean and the tropical Indian Ocean, thus answering the questions raised in the introduction of this paper. The main conclusions are as follows.

(1) The wave packet is the superposition of eastward travelling Kelvin waves and westward travelling Rossby waves with the slowest speed, and satisfies the boundary conditions of eastern and western coasts respectively. The configuration of the flow field and thickness anomalies correspond to positive and negative thickness field perturbations in east-west zonal flows, which is consistent with the classic equatorial Kelvin wave.

(2) The coefficient determining the degree of attenuation of the wave packet to the north and south sides of the equator is inversely proportional only to the square root of the product of the reduced gravity and the mean depth of the upper water, which means it is inversely proportional to the phase velocity of Kelvin waves in the upper waters. The wave packet has the same speed of decay when the phase velocity is fixed.

(3) The oscillation frequency of the wave packet is proportional to its mode number $l$ and the above phase velocity of Kelvin wave and is inversely proportional to the length of the equatorial ocean in the east-west direction. This oscillation frequency is the natural frequency of the ocean. Smaller mode number of the solution, lower phase velocity and longer length correspond to lower frequency and longer oscillation period; the oscillation period is the longest in Mode 1 of the solution.

(4) The flow anomalies of the wave packet of Mode 1 appear as zonal flows with the same direction most of the time, reach the maximum at the center of the equatorial ocean and decay rapidly away from the equator in the form of equatorially trapped waves.

(5) The flow anomalies of the wave packet of Mode 2 appear as zonal flows with the same direction most of the time in half of the ocean, and are always 0 at the center of the entire ocean, indicating stagnation. The flow anomalies decay to the north and south sides of the equator with the same speed as that of Mode 1.

(6) The spatial structure and oscillation period of the wave packet solution of Modes 1 and 2 are consistent with the changing periods of the surface spatial field and time coefficient of the first and second modes of the complex EOF analysis of the flow anomalies in the actual equatorial ocean. This indicates that the solution does exist in the actual ocean, and that the ENSO and the IOD are both related to Mode 2.

(7) The length of the bounded equatorial ocean is the sum of that of the tropical Indian Ocean and the tropical Pacific Ocean, so that this wave packet can explain the decadal variability (about 20 a) of the equatorial Pacific and the equatorial Indian Ocean after considering the Indonesian throughflow.

It is worth noting that, as this is a theoretical investigation, the model is simplified where necessary, (considering the east and west coasts in the North-South direction, etc.) leading to certain differences from the actual situations. It is therefore natural and understandable that our results differ slightly from those obtained through the complex EOF analysis of the actual ocean despite overall consistency.