Abstract

The conceptual recharge oscillator model intuitively established in Part I is derived from a dynamical framework of a Cane–Zebiak type model for tropical ocean–atmosphere interaction. A two-strip approximation to the equatorial ocean dynamics and one-strip approximation to the SST dynamics are employed to obtain a stripped-down coupled model that captures the main physics of the Cane–Zebiak type model. It is shown that the conceptual recharge oscillator model can be obtained from the stripped-down coupled model with a two-box approximation in the zonal direction or a low-frequency approximation to filter out high-frequency modes. Linear solutions of the stripped-down model are analytically solved and the dependence of coupled modes on various model parameters is delineated. In different parameter regimes, the stripped-down coupled model describes a coupled-wave mode and a mixed SST–ocean-dynamics mode that results from the merger of a nonoscillatory ocean adjustment mode with an SST mode. These two coupled oscillatory modes undergo a further merger. In the neighborhood of this merger, the leading mode of the system becomes a generalized mixed mode. It is suggested that the slow ENSO regime can be best characterized by this generalized mixed mode whose essential physics are described by the conceptual recharge oscillator model proposed in Part I.

1. Introduction

The El Niño–Southern Oscillation (ENSO) has been recognized as a natural coupled basinwide oscillation of the tropical Pacific ocean–atmosphere system. The positive feedback processes of tropical ocean–atmosphere interaction and the recharge/discharge process of the equatorial heat content through ocean adjustment dynamics are hypothesized as the growth and the phase-transition mechanisms of ENSO, respectively (Bjerknes 1969; Wyrtki 1975, 1986; Cane and Zebiak 1985, BWCZ hereafter). ENSO-like variability has been simulated in a variety of coupled models from the intermediate coupled models (e.g., Zebiak and Cane 1987, hereafter ZC) and the hybrid coupled models (Neelin 1990; Latif et al. 1991) to the state of the art coupled GCMs (e.g., Philander et al. 1992).

Theoretical studies have focused on the coupled instability of the basic climatological state of the tropical ocean–atmosphere system to provide the understanding of the origin of the ENSO-like variability in various coupled models and the ENSO phenomenon in nature (e.g., Philander et al. 1984; ZC; Hirst 1986, 1988; Battisti 1988; Battisti and Hirst 1989, BH hereafter; Suarez and Schopf 1988 and Schopf and Suarez 1988, SS hereafter; Wakata and Sarachik 1991; Neelin 1991; Jin and Neelin 1993a,b and Neelin and Jin 1993, hereafter JN). A great deal of progress has been made through these studies. It has been recognized that the subsurface ocean dynamics, oceanic Ekman layer dynamics, and SST thermodynamics are major contributors to the coupled interannual variation modes with different characteristics (such as westward and eastward propagating modes and standing oscillatory modes) in coupled tropical ocean–atmosphere models.

Along with the development of instability theories, a number of prototype models have also been proposed to describe the gross features of the low-frequency oscillation in the coupled ocean–atmospheric system (e.g., McWilliams and Gent 1978; Graham and White 1988; SS; BH; Cane et al. 1990 and Münnich et al. 1991, CMZ hereafter). Among them, the concise delayed oscillator model (SS; BH) best presents the basic features of the ENSO-like oscillation and describes the coupled positive feedback for the growth of SST anomalies over the central to eastern Pacific and the phase transition due to subsurface ocean wave dynamics. Further attempts at improving the conceptual delayed oscillator have focused on detailed analyses of oceanic wave dynamics (Schopf and Suarez 1990; CMZ; Schneider et al. 1995). For example, at the so-called fast-SST limit, which assumes the equatorial SST anomaly to be in quasi-equilibrium with subsurface thermocline fluctuation, CMZ discovered that the coupled process modifies the ocean basin modes (Cane and Moore 1981) to form a coupled wave oscillator described by mapping equations. In fact, CMZ showed that the BH delayed oscillator is an approximation to the coupled wave oscillator at the fast-SST limit, and thus the physics of the coupled wave oscillator and the BH delayed oscillator are the same. However, using a simplified version of the ZC model, JN showed that the coupled wave mode of CMZ and the coupled SST modes of Neelin (1991) are combined to form a mixed SSTocean dynamics mode in the most relevant parameter regimes. Although analytical solutions have been obtained by JN and CMZ in different parameter regimes surrounding the most relevant part of parameter space for the mixed mode, the merger of the two kinds of modes was only demonstrated numerically by JN. It was further speculated in JN that the delayed oscillator is an ad hoc representation of the mixed mode, which was apparently in contradiction to the findings of CMZ.

The conceptual recharge oscillator model proposed in Part I of this paper provides not only an analytical approach to the mixed SST–ocean dynamics mode but also a direct and concise representation of the BWCZ hypotheses. In the recharge oscillator model, a nonoscillatory ocean dynamics mode describing the adjustment process of the subsurface ocean to wind forcing and an SST mode for the SST thermodynamics merge into the mixed oscillatory mode through the coupled interaction. This mixed mode gives rise to the ENSO-like recharge oscillator whose phase transition is clearly depicted as the result of the recharge/discharge process of the equatorial heat content. Different from the delayed oscillator, the conceptual recharge oscillator model does not explicitly require a prescribed time delay. Yet the recharge oscillator model embodies the BH delayed oscillator, which can also be interpreted physically as a recharge oscillator.

In Part II, the recharge oscillator model is derived from a ZC-type coupled model through a number of approximations to provide a solid ground of the conceptual recharge oscillator established in Part I. Section 2 gives the derivation of the stripped-down version of the ZC-type model. Section 3 presents a derivation of the conceptual model from the two-strip coupled model and some numerical solutions of the latter to demonstrate that the basic characteristics of the recharge oscillator are captured by the conceptual model. Section 4 describes regime behavior of the analytical eigen solutions of the two-strip model. Further discussion is given in section 5 on the different approximations leading to the prototype models for the recharge oscillator, the coupled wave oscillator, and the BH delayed oscillator, and a summary is given in section 6.

2. A two-strip ZC-type model

The subsurface ocean dynamics in the ZC-type coupled model are governed by a linear shallow water model. Advective nonlinearity in the shallow water model is neglected because the applicability of shallow water dynamics to the three-dimensional ocean would fail before such terms become important. The small effects of SST on the dynamics through density changes are also neglected. This kind of shallow water model for ocean dynamics under the long-wave approximation on the equatorial β-plane can be written as (e.g., Cane and Sarachik 1981)

 
formula

where u and υ are the zonal and meridional upper layer ocean currents, h is the departure of the upper-layer thickness from the reference depth H and represents a thermocline depth anomaly field, and g′ is the reduced gravity parameter that gives rise to an internal gravity speed as c0 = g′H. The meridional component of the wind stress is neglected because the wind perturbations are predominantly zonal.

System (2.1) is then nondimensionalized by the oceanic Rossby deformation radius Lo = c0 /β for y, the ocean basin width L for x, the Kelvin wave crossing time L/c0 for t, the reference depth H for h, the Kelvin wave speed c0 for u and υ, and ρHc20/L for the wind stress variables τx, respectively. By further eliminating velocity variables in (2.1), a single equation for the thermocline depth is obtained:

 
formula

Throughout this paper, H is chosen as 150m; L /c0 and the weak damping rate εm are taken as 2 months and (2.5 yr)−1, respectively. The boundary conditions (e.g., Cane and Sarachik 1981) of no normal motion at the east boundary and zero integrated mass flux at the west boundary can be expressed as

 
formula

by using the second equation of (2.1). The shallow water system (2.2) and (2.3) retains all equatorial long Rossby waves and the Kelvin wave.

For simplification of the dynamics, hemispheric symmetry of the system is assumed and only the changes of thermocline depth in an equatorial strip and an off-equatorial strip will be taken into consideration. Such a simplification will be referred to hereafter as a two-strip approximation. Near the equator, for instance, within one oceanic Rossby radius of deformation, the thermocline depth approximately has the form of h(x,y,t ) = he(x,t )+y2Δh(x,t ) when the system is symmetric about the equator. The meridional variation of thermocline depth near the equator is consistent with a zonal current in the equatorial strip. Here he is the value of h at the equator and Δh is related to the zonal current along the equator. Letting h|y = 1 = [he(x,t ) + hn(x,t )]/2, where hn is the thermocline depth in the off-equator strip and is evaluated separately, one gets Δh = (hnhe)/2. The zonal current along the equator is u = hehn. Substituting this form of h into (2.2), the ocean dynamics equation along the equator (y = 0) can be then approximately written as

 
(∂t + εm)(hehn) + ∂xhe = τxe,
(2.4)

where τxe is the wind stress anomaly evaluated in the equatorial strip. Incidentally, with the approximate form of h near the equator, Eq. (2.4) is the same as the zonal momentum equation evaluated along the equator.

Equation (2.4) describes the Kelvin wave signal along the equator with an inclusion of the effect of Rossby waves because hn is the thermocline depth associated with Rossby wave signals in the off-equatorial strip. For the off-equatorial strip centered at yn, the first term in the left-hand side of equation (2.2) dominates the second term to balance the third term when yn is large (yn ≥ 2, for instance). Thus, the second term can be neglected, the ocean adjustment dynamics in this off-equatorial strip is dominated by long Rossby waves, and the equation for hn can be approximately written as

 
(∂t + εm)hn− ∂xhn/y2n = ∂y(τx/y)|y=yn.
(2.5)

This long Rossby wave equation has been used, for example, by Kessler and McPhaden (1995) to diagnose the thermocline variability off-equator, for instance at 5°N. When yn = 2 (about 5°–6°N), the Rossby wave speed is one-quarter of the Kelvin wave speed. It is close to one-third for the first free Rossby wave mode in the shallow water system (2.2). To some extent, the long Rossby wave in (2.5) can be viewed as a collective representation of the equatorial Rossby waves in the strip. Equations (2.3)–(2.5) are the two-strip approximation to (2.1). They qualitatively capture the basic dynamics of equatorial ocean adjustment under the Kelvin and long Rossby waves. They are somewhat similar to equations for the oceanic Kelvin wave and the first equatorial Rossby wave obtained by applying two-mode truncation to system (2.1). However, the two-strip approximation is better than the two-mode truncation because the severely truncated model will suffer great error in representing the wind stress structure and the boundary conditions, due to the well-known slow convergence of the Hermite function expansion.

The boundary conditions for the two-strip ocean dynamics model can be written as

 
hn(1,t) = rEhe(1,t), he(0,t) = rWhn(0,t ),
(2.6)

where xE = 1 and xW = 0 are used and rW and rE are two parameters introduced for considering various boundary conditions. For instance, when rE = 1, the eastern boundary condition in (2.6) is equivalent to that in (2.3), which requires a complete reflection of the Kelvin wave signal at the eastern boundary into Rossby waves and thus a zero normal velocity at all latitudes. When rE < 1, the eastern boundary in (2.6) becomes less restricted than that in (2.3) and represents a more general condition of an imperfect reflection. It allows an equatorial zonal mass flux at the eastern boundary, proportional to uE = (1 − rE)he(1, t ), and implies the mass exchange between the equatorial strip and off-equatorial regions through eastern boundary layer currents. The western boundary condition in (2.3) requires the integrated mass flux to vanish. In the equatorial band, the zonal fluxes associated with the incoming Rossby waves and the reflected Kelvin wave normally yield a net zonal flux proportional to uW = −(1 − rW)hn(0, t), which is compensated by western boundary currents. When rW = 3/5, the western boundary condition in (2.6) is equivalent to the integral in (2.3) under the two-strip approximation with h = he, hn, and 0 at y = 0, yn, and 2yn and beyond, respectively. More generally, rW < 1 is valid because of energy leakage through the western boundary under conditions in (2.3) (see a proof in appendix B) and rw can also be viewed as a measure of the reflection efficiency of Rossby waves at the western boundary. Thus, both rW and rE are chosen as model parameters that conveniently measure the contributions to the recharging/discharging of the equatorial heat content through the boundary currents.

Under a slowly varying wind forcing, the Kelvin and Rossby waves eventually establish a zonal pressure gradient and set up the quasi-Sverdrup balance in the equatorial strip. The Rossby waves generated in the off-equatorial strip will continuously be reflected to Kelvin waves at the western boundary. The associated total meridional transport by both the interior ocean current across the basin and the mass fluxes at the eastern and western boundaries will determine the change of the zonal mean thermocline depth in the equatorial strip, which resides in the memory of the ocean dynamics as it is emphasized in the conceptual model of Part I. Therefore, the two-strip model describes one of the key elements in the BWCZ hypothesis, a phase transition mechanism through the recharge/discharge of the equatorial heat content.

The unforced shallow water system under a long-wave approximation gives slightly damped free oscillations, known as ocean basin modes (Cane and Moore 1981) and scattering modes (JN), due to the reflections of Kelvin and Rossby waves. Collectively they will be referred to as ocean adjustment modes. The two-strip model captures the gross feature of the free ocean adjustment modes. One can obtain these free modes by assuming a solution in the form of B1eσt+(σ+εm)y2nx for hn and B2eσt−(σ+εm)x + hn/(1+y2n) for he. Substituting this solution form into (2.4)–(2.5) and using boundary condition (2.6), one finds that for nontrivial solutions, the eigenvalue σ must satisfy the following dispersion relationship:

 
formula

The final eigen solution can be written as

 
formula

where B is arbitrary amplitude and c.c. is the complex conjugate of the term in the front of it. These free ocean adjustment modes of the two-strip model are much like the damped standing oscillations of a string with two slip boundary conditions. The second term in the right-hand side of (2.7a) always has a negative real part. Thus these ocean adjustment modes are decaying due to the weak damping εm and the energy loss accompanying the mass exchanges through the western and eastern boundary reflections represented by rE and rW. There is a purely decaying basinwide adjustment mode (j = 0), which results from the fact that basinwide mean of h is a free-mode solution of system (2.1). In this two-strip approximation, this mode can be generally viewed as an approximation to a number of low-frequency gravest scattering modes found in JN. The decaying uncoupled ocean adjustment mode in the conceptual model is an explicit indication of its potential role in the coupled dynamics. The oscillatory modes ( j = ± 1,2,3, ···) are similar to the ocean basin modes first found by Cane and Moore (1981). The one with the largest spatial scale ( j = ±1) has a period of the basin crossing time for the equatorial Kelvin wave and the Rossby wave. At yn = 2, the period is about 10 months. Clearly, the two-strip model successfully captures the equatorial ocean adjustment dynamics in the ocean basin under a wind stress forcing.

The change of SST may be described by the thermodynamics of a constant depth mixed layer embedded in the upper-layer ocean. This simple approach was proven to be successful, for instance, in the ZC model. Motivated by the fact that the strongest SST anomaly in response to upwelling, advection, and thermocline depth changes is confined in a narrow band along the equator, Neelin (1991) and JN used an equatorial-strip approximation to the simple thermodynamics equation of the SST. Adopting this equatorial-strip approximation to the SST equation, and considering only the dominant processes, one gets the following equatorial SST anomaly equation linearized about an upwelling climate state:

 
formula

where

 
formula

The equatorial SST anomaly Te in (2.8) is nondimensionalized by ΔT, which is set as 7.5°C, and the time variable and thermocline depth are nondimensionalized in the same way as in the ocean dynamics equation. The terms on the right-hand side of (2.8) represent a Newtonian-type damping process with the damping rate εT and the vertical advection of the equatorial SST anomaly Te and subsurface temperature anomaly Tsub by the mean upwelling w1 of the equatorial strip. Here, H1.5 is a constant depth at which w1 is best characterized. The subsurface temperature anomaly, which follows the thermocline depth change, is essential in controlling the SST of the central to eastern equatorial Pacific where upwelling is strong and the thermocline is shallow. As in ZC, Tsub is parameterized by a nonlinear function of the thermocline depth. Linearizing this function with respect to the observed climate state, Tsub becomes linearly proportional to he with a factor γ0 = ∂hTsub. Values of c(x) and γ(x) depend on longitude because of nonuniform longitudinal distributions of upwelling and γ0. Other anomalous vertical and horizontal advection terms abbreviated in (2.8), although quantitatively important in determining the equatorial SST anomaly (e.g., Battisti 1988), will be ignored throughout this paper for the sake of simplicity. This simplified SST equation is adequate for understanding the basic features of interannual coupled modes dominated by the thermocline feedback process, whereas the coupled modes with the full SST equation of the equatorial strip have been numerically explored in JN and Jin et al. (1996).

In the dynamics of the tropical lower troposphere, mechanical balances are between the pressure gradient, the Coriolis force, and boundary layer friction, while the thermodynamic balance is between diabatic heating and divergent flow acting in both moisture convergence and adiabatic cooling to a zero-order approximation. Although simple models based on these balances differ in physical interpretation (Lindzen and Nigam 1987; Neelin and Held 1987; Zebiak 1986), all approximately give an equivalent Gill-type model (Gill 1980) through proper scaling (Neelin 1989). When heating induced by an equatorial SST anomaly is assumed to be proportional to the SST perturbation with a Gaussian meridional structure depending on the atmospheric Rossby radius, the Gill-type atmospheric dynamics results in a largely symmetric zonal wind stress anomaly. It can be approximately written as

 
τx(x,y,t ) = μA(Te,x) exp(−(yLo/La)2/2),
(2.9)

where La is the atmospheric Rossby radius of deformation and A(Te,x) is a nonlocal operator that relates the SST anomaly in the equatorial strip to the wind stress anomaly. As in Part I, a relative coupling coefficient μ is introduced as a free parameter for exploring the regime behavior of the system. The stress terms in the equations (2.4) and (2.5) can thus be expressed as

 
formula

where θ is about 1.0∼1.25, depending on yn and the ratio of the oceanic and atmospheric Rossby radii of deformation.

The variable A(Te,x) can be determined based on the Gill-type model dynamics (e.g., JN) or through an empirical relation of the equatorial wind stress and SST anomaly. Using linear regression, Deser and Wallace (1990) showed that, corresponding to the SST anomaly of the central to eastern Pacific, the atmospheric wind stress anomaly patterns exhibit an extended patch of the wind anomaly from the central to western Pacific with little anomaly in the other part of the equatorial Pacific band. For simplicity, a simple nonlocal relationship between the SST anomaly and the wind stress anomaly is assumed here:

 
A(Te,x) = b0TeEf(x),
(2.11)

where TeE is the SST anomaly averaged over the eastern half of the equatorial strip and f (x) is a normalized function whose zonal integration is unity. The nondimensional factor b0 is 2.5 as estimated in Part I.

Equations (2.4)–(2.6), (2.8), (2.10), and (2.11) form a simple linearized two-strip coupled model that captures the main physical processes incorporated in the ZC-type model. It allows analytical solutions and insight into the physics of coupled interannual modes arising from the ocean–atmospheric interaction around the climatological basic state.

3. Recharge oscillation in the two-strip model

a. Two-box approximation

SST anomalies in reality and in the ZC-type models are mainly trapped in the central to eastern Pacific where the feedback processes for the SST change are strongest due to the shallow thermocline depth and strong upwelling in the climatological basic state, whereas little SST anomaly occurs in the western Pacific. The equatorial belt thus can be roughly divided into an eastern and a western half, and the SST anomaly in the western Pacific is negligible. Then, the SST equation (2.8) averaged over the eastern half of the basin can be written as

 
formula

where subscript E denotes variables over the eastern half of the basin. Equation (3.1) is the same as the SST equation (2.4) of Part I.

Using upstream differences for (2.4) and (2.5) [evaluating the time-derivative terms in (2.4) at the middle point of the eastern basin, the time-derivative term in (2.5) at the middle point of the western basin, and the x-derivative terms in both equations at the center of these two points by finite difference], one obtains a crude two-box approximation to the ocean dynamics:

 
formula

Relations of hew = rwhnw and hnE = rEheE from the boundary condition (2.6) and Eq. (2.10) have been used, and AM(TeE) is the wind stress averaged over the middle half of the basin and Δx = 1/2 is the half-basin width. An additional term, (1 − rE)(∂/∂t + εm)heE, on the right-hand side of the first equation of (3.2) is neglected because the timescale of the mode retained in the approximation is much longer than the basin crossing time of the Kelvin wave. The above approximation reduces the ocean adjustment process to the simplest version as proposed in the conceptual model in Part I. By one-step differencing along the Rossby and Kelvin wave characteristics, we have neglected the detailed wave propagation process. The two-box truncation thus filters out all oscillatory ocean adjustment modes, but still roughly describes the basinwide mass adjustment. In particular, the slowly and purely decaying free mode of the ocean adjustment is retained. When the system is forced by an anomalous wind stress, it describes the recharge/discharge process of the zonal mean equatorial thermocline depth. The energy loss due to the mass exchange through the boundaries (rE ≤ 1, rw < 1) is clearly a major contributor to the decay timescale of ocean adjustment, and the condition of 0 < θrE ≪1 is another factor for characterizing the slow basinwide ocean adjustment.

From (2.11), AM(TeE) is proportional to the SST anomaly. By assuming that the wind stress anomalies largely occur over the middle part of the basin and the zonal average of the wind stress is dominated by the domain average in the middle half of the basin, one gets

 
AM(TeE) = μb0TeE /Δx.
(3.3)

Slightly reorganizing the Eqs. (3.1)–(3.3), one reduces the two-strip model to exactly the same form of the conceptual recharge oscillator model as proposed in Part I,

 
formula

by denoting rwhnw as hw, TeE as TE, b = b0μ, R = γb − c, and replacing the parameters r, c, α, and γ as

 
formula

Thus, through systematic and reasonable simplifications, the conceptual model in Part I is obtained from the dynamic framework of a ZC-type coupled model. This derivation also reveals the relevant ranges of the parameters in the conceptual model. For instance, with rE = 0.9, rW = 0.75 (these two numbers are roughly estimated from a standard ZC coupled simulation), θ = 1.2, cE = 1, γE = 1, and Δx = 0.5, one finds that the parameters (r,c,b0,α,γ) have similar values to those in the conceptual model of Part I.

b. Numerical solutions

The relevance of the box model approximations can be demonstrated by comparing the characteristics of the solutions between the two-strip model and the conceptual model in Part I. The following simple form of the coefficient in the SST equation and a simple function f(x) describing the patch of wind stress,

 
formula

are used to mimic the basic state conditions applied to the box model; x1 = 1/4 and x2 = 3/4 are chosen. The linear eigen solutions of the two-strip model with the specifications in (3.6) are calculated here numerically with a high resolution in x (50 grid points), whereas analytical solutions will be further discussed in the next section. It should be pointed out that the coupled solution is insensitive to the shape of f(x) (e.g., CMZ; Schopf and Suarez 1990). The solutions of the model are virtually unchanged when other functions—for instance, a half-sine function sin [(xx1)π/(x2x1)] π/2(x2x1)—are used for describing the wind patch.

An example of the dependence of eigenmodes on the coupling coefficient is shown in Fig. 1. At zero coupling, the eigenvalues in the plot correspond to one SST mode and two oceanic adjustment modes, including the nonoscillatory and gravest ocean wave modes. When the coupling coefficient increases, the decay rate of the SST mode is effectively reduced due to the positive feedback process. This damped SST mode eventually merges with the damped nonoscillatory ocean adjustment mode to produce a mixed oscillatory mode, which is the only mode becoming unstable at sufficiently strong coupling. This mode breaks down into two real modes when the coupling coefficient is further increased. This behavior is very similar to what has been seen in Fig. 2 of the conceptual model in Part I. The gravest coupled wave mode, on the other hand, is a damped mode throughout the range of the coupling coefficient in this case. The instability of the mixed mode of the system gives rise to a low-frequency Hopf bifurcation with a period of about 4 yr near the bifurcation point. A finite amplitude solution is achieved by adding a term of −enh3e, the same form of simple thermocline nonlinearity, into the SST equation (2.8). A typical solution shown in Figs. 2a and 2b exhibits a standing oscillation in both the thermocline depth and the SST. The time series of the thermocline depth at the western Pacific and the SST in the eastern Pacific and their trajectory in the phase space (Figs. 2c and 2d) are also in good agreement with the conceptual model. In fact, the time lag between the maximum negative thermocline depth in the western Pacific and the maximum positive SST in the eastern Pacific is shortened to about 3 months from 6 months for the solutions of the conceptual model in Part I. Thus, compared with the result from the conceptual model (Fig. 3 of Part I), the shape of the elliptical limit cycle here resembles somewhat better the result of the ZC model and the observed feature shown in the Part I.

Fig. 1.

Dependence of the eigenvalues of the first two modes on the relative coupling coefficient. The heavy (thin) solid curve is for growth rate, and the heavy (thin) dashed curve is for frequency of the leading (second) mode. Dividing π/3 by the frequencies yields period in years. The parameters are set as rW = 0.75, rE = 0.9, εm = 1/15, θ = 1.2, and yn = 2.

Fig. 1.

Dependence of the eigenvalues of the first two modes on the relative coupling coefficient. The heavy (thin) solid curve is for growth rate, and the heavy (thin) dashed curve is for frequency of the leading (second) mode. Dividing π/3 by the frequencies yields period in years. The parameters are set as rW = 0.75, rE = 0.9, εm = 1/15, θ = 1.2, and yn = 2.

Fig. 2.

The time–longitude plots of the (a) SST anomaly (in °C) and (b) thermocline depth (in m) along the equator, and the (c) time series plot and (d) trajectory plot for the thermocline depth and SST averaged over the western and eastern half of the basin for the nonlinear solution with μ = 0.9 and en = 3. The other parameters are the same as in Fig. 1.

Fig. 2.

The time–longitude plots of the (a) SST anomaly (in °C) and (b) thermocline depth (in m) along the equator, and the (c) time series plot and (d) trajectory plot for the thermocline depth and SST averaged over the western and eastern half of the basin for the nonlinear solution with μ = 0.9 and en = 3. The other parameters are the same as in Fig. 1.

Fig. 3.

Dependence of the growth rate (solid curve) and frequency (dashed curve) of the leading mode from Eq. (4.4) on the relative coupling coefficient for δ = 1 (a), δ = 0.01 (b), and δ = 100 (c). In plot (b), growth rates are rescaled by δ. The other parameters in these plots are the same as in Fig. 1.

Fig. 3.

Dependence of the growth rate (solid curve) and frequency (dashed curve) of the leading mode from Eq. (4.4) on the relative coupling coefficient for δ = 1 (a), δ = 0.01 (b), and δ = 100 (c). In plot (b), growth rates are rescaled by δ. The other parameters in these plots are the same as in Fig. 1.

4. Analytical eigen solutions

The two-strip coupled model allows analytical solutions for certain x dependence of the coefficients in the SST equation (2.8), such as specified in (3.6). In this case, one can actually approximate the thermocline depth in the eastern half of the equatorial strip by its value at the eastern boundary. This is a good approximation for the leading coupled oscillator mode because the thermocline depth of the leading mode is almost uniform over the eastern Pacific as shown in Fig. 2b. With this simplification, one can use Eq. (3.1) to replace (2.8), and the two-strip model can be rewritten as

 
formula

where heE equals he(0,t ) and the thermocline depths of the equatorial and off-equatorial strip satisfy the boundary condition (2.5). Assuming an eigen solution of (4.1) has the form

 
formula

and substituting (4.2) into (4.1), one obtains

 
formula

With the boundary conditions e(0) = rw n(1) and n(1) = rE e(1), a dispersion relationship is obtained for nontrivial solutions of (4.1) as follows:

 
formula

where x* is the center position of the wind patch and

 
formula

For the convenience of discussion, a free parameter, δ, the ratio between the timescales of SST change and ocean wave dynamics, is again introduced following the approach in JN. Here, δ−1 is attached to the time derivative of SST, δ−1 = 0 represents the fast-SST limit, and δ−1 ≫1 represents the fast-wave limit.

When the system is uncoupled, or at the low coupling limit, there exists an uncoupled SST mode with a damping rate σT = −δcE and a set of uncoupled ocean adjustment modes whose eigenvalues are the same as expressed in (2.7). In the coupled case, a new mode is formed as shown in Fig. 3a; that is, the decaying SST mode merges with the decaying nonoscillatory ocean adjustment mode ( j = 0) to form an oscillatory mixed mode that becomes unstable at a moderate coupling coefficient. Under the same choice of the model parameters, the eigenvalues of the leading mode from (4.4) are very similar to the results shown in Fig. 1. The ranges of the coupling coefficient for the leading mixed mode to exist are slightly different in Fig. 3a from Fig. 1 because of the approximation of average thermocline depth over the eastern half of the equatorial Pacific by the thermocline depth at the eastern boundary. The difference is not essential because it can be eliminated by rescaling of the relative coupling coefficient.

It has been noted that the recharge oscillator, or the mixed mode, does not exist at the fast-wave limit and the fast-SST limit in the conceptual model as discussed in Part I. The leading eigen mode of the two-strip model at these two limits is shown in Figs. 3b and 3c. At the fast-wave limit, all the ocean adjustment modes are over damped as coupling increases, and only the SST mode is destabilized and eventually becomes a growing mode at a strong coupling coefficient. At the fast-SST limit, the leading unstable mode comes from the destabilization of the gravest oscillatory ocean basin mode. The coupling process not only leads the mode to its instability but also reduces its frequency. The system thus gives rise to an unstable oscillatory mode that is the coupled wave oscillator found in CMZ. The unstable coupled wave oscillator again breaks down and leads to a purely growing mode at a moderate coupling. This coupled wave mode at fast-SST limit has a relatively higher frequency near the Hopf bifurcation (around 1 yr).

Between these two extreme limits, there exist two mergers of different kinds of modes. The first merger (a multiplicity-two degeneracy of real modes, see appendix A) gives birth to the mixed SST–ocean dynamics mode. It occurs in the parameter space away from the fast-wave limit. Farther away from the fast-wave limit, there is a second kind of merger (a multiplicity-two degeneracy of complex modes, see appendix A) between the mixed SST–ocean dynamics mode and the coupled wave mode originated from the gravest ocean basin mode to form a generalized mixed SST–ocean dynamics mode. It is general in the sense that both basinwide ocean adjustment and the detailed wave dynamics shape the characteristics of the mode in terms of its instability and frequency. However, this generalization is not essential because the conceptual recharge oscillator model involves one merger, whereas this generalization involves two mergers to characterize a series of mergers that could occur between the two extreme limits, as was shown in JN in a more complex model.

Because these degeneracies (mergers of modes of different origins) mark the boundaries of qualitative changes in the characteristics of the coupled modes in the two-strip model, tracing these degeneracies can provide a sketch of the regime diagram in the parameter space. This can be done with the technique described in appendix A. For the first kind of degeneracy, it occurs on a (k−1)-dimensional surface in k-dimensional parameter space (see appendix A). For example, on the k = 2 parameter plane (μ, δ), this surface becomes a curve (or a few curves); the second kind of degeneracy occurring on a k−2 surface appears as a point (or a few points) on the parameter plane (μ, δ). A typical regime diagram is shown in Fig. 4. The small circle with a cross inside labels the point where two complex modes merge and the thin line passing it divides the regime into two parts. In the regime where δ is smaller than that indicated by the thin line, the leading mode that becomes unstable as coupling increases comes from the mixed mode resulting from the first kind of merger. In this part of parameter space, the dash–dot and solid curves indicate the occurrence of the first kind of merger and the breakdown of the mixed mode. It is clear that the mixed mode regime does not extend to the fast-wave limit. Similar to those shown in Fig. 1, all coupled wave modes are decaying and never break down into real modes in this part of parameter space. For the standard value of δ = 1, the mixed mode is the leading mode undergoing the Hopf bifurcation. In the regime where δ is larger than that indicated by the thin line, the leading mode has its root from the gravest ocean basin mode. Although the first kind of merger still exists, the mixed mode from the first kind of merger is always a decaying mode and never breaks down into a real mode again. The leading mode undergoing the Hopf bifurcation is the gravest coupled wave mode, which also breaks down, as indicated by the solid curve that extends continuously from the fast-SST limit toward the fast-wave limit. Thus the second kind of merger indicates the switch between the two kinds of oscillatory modes as a leading mode. In the neighborhood of this merger, the leading mode undergoing the Hopf bifurcation of the system has the characteristics of both the mixed SST–ocean dynamics mode from the first merger and the coupled wave mode destabilized from the gravest ocean basin mode.

Fig. 4.

The regime diagram obtained from Eq. (4.4) in the parameter space of (δ, μ) with cE = 1, γE = 1.0, rW = 0.75, rE = 0.9, εm = 1/15, θ = 1.2, Δx = 0.5, x* = 0.5, and yn = 2 (δ = 0, ∞ correspond to the fast-wave limit and fast-SST limit, and δ = 1 is the standard case). The solid curve marks the breaking down of the leading oscillatory mode into two real modes. The dot–dash curve marks formation of the mixed oscillatory mode that exists in the regime to its right. In the regime between these two curves, there are two leading oscillatory modes, the mixed oscillatory mode and the gravest coupled wave mode. The dashed curve indicates the stability boundary, and one oscillatory mode is unstable in the stripped regime. The circle with a cross indicates the degeneracy point where the two leading oscillatory modes merge to share the same eigenvalue. The thin line passing this point distinguishes the unstable oscillatory mode as either the mixed oscillatory mode (in the regime stripped by heavy solid lines) or the gravest coupled wave mode (in the regime stripped by heavy dotted lines).

Fig. 4.

The regime diagram obtained from Eq. (4.4) in the parameter space of (δ, μ) with cE = 1, γE = 1.0, rW = 0.75, rE = 0.9, εm = 1/15, θ = 1.2, Δx = 0.5, x* = 0.5, and yn = 2 (δ = 0, ∞ correspond to the fast-wave limit and fast-SST limit, and δ = 1 is the standard case). The solid curve marks the breaking down of the leading oscillatory mode into two real modes. The dot–dash curve marks formation of the mixed oscillatory mode that exists in the regime to its right. In the regime between these two curves, there are two leading oscillatory modes, the mixed oscillatory mode and the gravest coupled wave mode. The dashed curve indicates the stability boundary, and one oscillatory mode is unstable in the stripped regime. The circle with a cross indicates the degeneracy point where the two leading oscillatory modes merge to share the same eigenvalue. The thin line passing this point distinguishes the unstable oscillatory mode as either the mixed oscillatory mode (in the regime stripped by heavy solid lines) or the gravest coupled wave mode (in the regime stripped by heavy dotted lines).

In the following, δ is always set to 1. In this case, it can be easily shown using (4.4) that if the degeneracies in Fig. 4 are traced and plotted on the parameter plane (cEμ, cE) with cE as the vertical axis and cEμ as the horizontal axis, the regime diagram plot will be identical to Fig. 4. Therefore, the collective damping rate cE in the SST thermodynamics, replacing the role played by δ, becomes a crucial parameter in determining the nature of the leading mode of the system. To explore the sensitivity of the coupled modes to other parameters of the two-strip model, one can trace the Hopf bifurcation and the second kind of degeneracy in a parameter space of high dimensions. A three-dimensional (3D) subspace is convenient for display of a 2D surface on which the Hopf bifurcation takes place and a 1D curve on which the mergers of the second kind occur. For instance, in the 3D subspace (μ, rE, rW), the Hopf bifurcation occurs at the critical coupling μ = μc(rE, rW), which is a 2D surface in this parameter space. The contour plot of μc is the projection of this surface on the parameter plane (rE, rW) as shown in Fig. 5a. At every point of the surface, there is a critical frequency σi = ωc(rE, rW), which is also shown in the contour plot Fig. 5b. The frequency and the critical coupling coefficient depend on the reflection parameters that measure the contributions to equatorial mass recharge/discharge by the western and eastern boundary currents. There is a low limit of about 0.4∼0.7 for the reflection parameter rW in order to have the leading mode of the system reach the Hopf bifurcation. When rW is below this value, the oscillation mode either disappears or remains stable because the mass flux through the western boundary tends to cancel the interior meridional transport that normally dominates the recharge/discharge process. When rW is close to 1, the mass flux through the western boundary is small, and the recharge/discharge process is more effeciently carried out by the interior meridional transport. Thus the frequency of the mode is relatively high and the value of the critical coupling coefficient is also relatively low. When rW is small, the mass flux through the western boundary is large and may even overwhelm the interior meridional transport. In this situation, a westerly equatorial wind stress anomaly associated with a positive SST anomaly in the eastern Pacific, for instance, leads to the recharging of the equatorial heat content. This recharging will further amplify the SST anomaly and a phase transition thus does not take place. This is also clearly shown in the expression (3.5), because when rW = 0, the mixed oscillatory mode does not exist in the two-box model. Thus the western boundary plays an important role in the existence of the recharge oscillator, which is consistent with the findings of BH based on their version of the ZC model and the BH delayed oscillator model.

Fig. 5.

Contour plots of critical values of the (a) coupling coefficient and (b) corresponding frequency of the leading mode of Eq. (4.4) at the Hopf bifurcations on the parameter plane (rE, γW) with εm = 1/15, θ = 1.2, yn = 2, γE = 1.0, and cE = 1. The heavy curve distinguishes two regimes where the mixed oscillatory mode (below the curve) or the coupled wave mode (above the curve) is the leading mode reaching the neutral stability.

Fig. 5.

Contour plots of critical values of the (a) coupling coefficient and (b) corresponding frequency of the leading mode of Eq. (4.4) at the Hopf bifurcations on the parameter plane (rE, γW) with εm = 1/15, θ = 1.2, yn = 2, γE = 1.0, and cE = 1. The heavy curve distinguishes two regimes where the mixed oscillatory mode (below the curve) or the coupled wave mode (above the curve) is the leading mode reaching the neutral stability.

The recharge oscillator continuously changes its frequency and instability in the full range of the reflection parameter of the eastern boundary. As illustrated in the next section, the BH delayed oscillator can be treated as the special case of the recharge oscillator with rE = 0. Reflection parameter rE measures the eastern boundary mass flux, which enhances the total recharge/discharge of equatorial heat content. A large rE implies a small mass flux through the eastern boundary and thus a decreased rate of total recharge/discharge of equatorial heat content. This is well illustrated by the expression of damping rate r and forcing factor α in (3.5) for the conceptual model. Thus the frequency of the recharge oscillator is reduced and a stronger coupling coefficient is also required for bringing the mode to the Hopf bifurcation as rE increases.

There are two regimes divided by the heavy curve as shown in Figs. 5a and 5b. This heavy curve is the projection of the curve on which the merger of the mixed oscillatory mode and the coupled wave mode occurs in the 3D parameter space (μ, rE, rW) onto the parameter plane (rE, rW). The coupling coefficients corresponding to the heavy curve are different from the values indicated by the contour plot in Fig. 5b because this merger and the Hopf bifurcation always happen at different coupling coefficients. In the regime above the curve, the leading mode undergoing the Hopf bifurcation comes from the coupled wave mode; below the curve, the leading mode is the mixed mode. The regime of coupled wave mode disappears in the plot if the damping rate cE of the SST anomaly is decreased, but it expands if cE increases. When cE is larger than 2, the whole regime is dominated by the coupled wave mode as shown in Fig. 6. Thus a strong damping rate cE pushes the system toward the fast-SST limit where the leading mode tends to have a higher frequency at the Hopf bifurcation.

Fig. 6.

Same as Fig. 5 except cE = 2.

Fig. 6.

Same as Fig. 5 except cE = 2.

The dependence on the other parameters can be explored by dropping out one of the parameters in the above 3D subspace of the parameters and adding a new one to form another 3D parameter space. For example, the wind patch as the atmospheric response to the SST anomalies is highly idealized for simplicity. The dependence of the eigenmodes on the width and central location of the wind patch provides a measure of the sensitivity of the modes in the two-strip model to the wind patch specifications. This sensitivity was discussed by Schopf and Suarez (1990) and CMZ for coupled wave modes. As shown by Schopf and Suarez (1990), the leading mode of the system is insensitive to the width of the wind patch, particularly when the frequency of the leading mode is low because sK ≈ 1 and sR ≈ 1 become good approximations and the dependence on the width disappears. As found by CMZ, the leading eigenmode is more sensitive to the location of the wind patch x*. Thus it is chosen to replace rE to form a new 3D parameter space (μ, rW, x*). As shown in Fig. 7 for a small rE, the gravest coupled wave mode is the leading mode undergoing the Hopf bifurcation when the wind patch is in the western Pacific, whereas the mixed mode is the leading mode when the wind patch is located in the central to eastern Pacific. This is also true for a large rE, except the regime undergoing the Hopf bifurcation is slightly shrunken compared with the result in Fig. 7 because a large rE implies less mass flux by the eastern boundary currents, and the recharge–discharge process can be more easily overwhelmed by the mass transport through the western boundary currents.

Fig. 7.

Same as Fig. 5 except on the parameter plane (x*, γW) with rE = 0.1.

Fig. 7.

Same as Fig. 5 except on the parameter plane (x*, γW) with rE = 0.1.

5. Relations among recharge oscillator, delayed oscillator, and coupled wave oscillator

If the weak damping εm is ignored, the second equation in (4.1) can be integrated along the Rossby wave characteristics. Similarly, the first equation in (4.1) then can be integrated along the Kelvin wave characteristics. Thus, the partial differential equations in (4.1) can be written in simple delay integral equations that can be further approximated by delay equations as follows:

 
formula

and

 
formula

Due to the fact that y2n + 1 ≫ 1, the terms divided by this factor have been ignored in the approximation equation (5.1). The first integral (5.1) and the integral in (5.2) have also been simplified by the mean value approximations. The following two delay equations can then be obtained by using the boundary conditions in (2.6):

 
formula

and

 
formula

where heW and heE denote he(0,t ) and he(1,t), respectively.

Delay equation (5.4) or Eqs. (5.1) and (5.3), together with the SST equation in (4.1) form a coupled delay-differential set whose dispersion equation can be written as

 
formula

This dispersion relation can be directly obtained from (4.4) by employing the same approximations used in (5.1) and (5.2). The mean value approximation is equivalent to sK ≈ 1 and sR ≈ 1, which are valid for modes with small |σ| or for a narrow wind patch. The terms divided by y2n + 1 in (4.4) can be ignored by the same consideration for the approximation in (5.1). The two dispersion equations (4.4) and (5.5) give qualitatively similar results for the dependence of the leading eigen modes on the model parameters as illustrated by an example shown in Fig. 8. In the following, it will be illustrated that this coupled delay differential model can be further reduced to various prototype models proposed for ENSO.

Fig. 8.

Same as Fig. 5 except for the dispersion equation (5.5).

Fig. 8.

Same as Fig. 5 except for the dispersion equation (5.5).

a. Recharge oscillator

The stripped-down model facilitates the derivation of the conceptual recharge oscillator model as shown in section 3. The delay equations (5.1) and (5.3) and the SST equation in (4.1) can also be approximated by a similar system of two first-order ordinary differential equations. This can be achieved by assuming that the timescale of the leading coupled mode is much longer than 1 + y2n, the period of the gravest free ocean basin mode. Thus one can use Taylor expansions to first order to express every term in Eq. (5.3) at time t and obtain

 
formula

where high orders of the Taylor expansion are ignored. Similarly, (5.1) can be approximated as

 
formula

The last term in Eq. (5.6) and the last two terms in Eq. (5.7) are all relatively small and can be ignored under the long timescale or low-frequency approximation. Then two first-order ordinary differential equations for the SST anomaly in the eastern Pacific and thermocline anomaly in the western Pacific are obtained as follows:

 
formula

This approximate form of the stripped-down coupled model is similar to that obtained under the two-box approximation. It should be noted that keeping the derivative terms in the right-hand side of (5.6)–(5.7) will also give a similar equation to (5.8) with slightly different coefficients. When the nonlinear thermocline feedback term needs to be considered, it can be simply added into (5.8) and (3.4) as well. The two-box truncation and the Taylor expansion under the low-frequency approximation both filter out all the fast and oscillatory ocean basin modes to allow a simple description of the low-frequency mixed SST–ocean dynamics mode that depicts the recharge oscillator mechanism.

The difference between the coefficients in (3.4) and (5.8) is due to the finite differencing in the two-box model approximation. As discussed in section 4 and shown in (3.5), the reflection parameters that measure the mass exchanges of the equatorial strip and off-equatorial regions through boundary currents are of great impact on the recharge oscillator. The critical coupling coefficient and corresponding frequency of (5.8) now can be expressed as

 
formula

where r = (1 − rWrE)/(1 + y2n) is the decay timescale of the zero-frequency ocean adjustment mode. Thus the frequency of the recharge oscillator at the critical coupling depends on the decay rate of the uncoupled ocean adjustment mode, the damping rate cE in the SST thermodynamics, and the crossing time (1 + y2n) of the equatorial Kelvin and Rossby. However, the condition 0 < rW(θrE) ≪ 1 is necessary for the recharge oscillator to exist and to have a period much longer than all these three intrinsic timescales of the system. A large rW is not a necessary but a sufficient condition for the existence of the recharge oscillator and it tends to shorten its period. However, rW ≠ 0 is a necessary condition, because otherwise the recharge oscillator will be destroyed. On the other hand, a large rE increases the period of the recharge oscillator that still exists even for rE = 0. Thus qualitatively, (5.9) agrees with the analyses in section 4 of the sensitivity of the leading mode to model parameters.

b. BH delayed oscillator

There have been several attempts to justify the BH delayed oscillator (CMZ 1990; Neelin et al. 1994). One assumption for establishing the BH delayed oscillator is to ignore the eastern boundary reflection (e.g., CMZ). By setting rE = 0 in (5.4), the thermocline depth at the eastern boundary only includes two parts, the locally forced Kelvin wave signal and the Kelvin wave from the western boundary due to the reflection of the forced Rossby wave. Using the SST equation of (4.1), one gets the familiar BH delayed oscillator

 
formula

where a small time delay, 1 − x*, related to the Kelvin wave propagation in the second term of the right-hand side of (5.4), has been ignored. The delay 1 + y2nx* in the last term of (5.10) is the time taken for the Rossby wave propagating from the center of the wind patch to the western boundary and the reflected Kelvin wave crossing the basin.

The physical mechanism of the phase transition of this delayed oscillator model can be interpreted by the recharge mechanism. Assuming that the zonal mean thermocline depth approximately equals the average of the thermocline depths at the eastern and western boundaries, it can be expressed as

 
formula

if rE = 0. This zonal mean thermocline depth results from the recharge/discharge due to meridional mass transport associated with the forced Rossby wave and mass fluxes through the western and eastern boundary currents. After the peak warm phase, ĥ becomes negative, implying that the discharge of the equatorial heat content has taken place. When the SST anomaly following a warm phase approaches zero, the east–west tilt of the thermocline depth is also near zero (if the small time lag 1 − x* is ignored, it becomes precisely zero); thus the thermocline depth at the eastern Pacific is the same as the zonal mean that is negative due to the slow discharge process during the warm phase. This slowness of the recharge/discharge process allows the negative zonal mean thermocline depth at the time of zero SST anomaly, which implies that the time delay 1 + y2nx* is smaller than a quarter period of the oscillation. This negative zonal mean thermocline depth brings the SST to the cold phase due to anomalous cold water being pumped into the surface layer by the climatological Ekman pumping, which is clearly expressed by (5.10) with a negative change rate of the SST at the time of the zero SST anomaly. Similarly one can understand the cold to warm phase transition. Thus the BH delayed oscillator can be viewed as a special case of the recharge oscillator. Figure 8 also shows clearly that the leading mode at rE = 0 is indeed embodied in the regime for the recharge oscillator.

c. Coupled wave oscillator at fast-SST limit

The coupled wave oscillator derived at the fast-SST limit, as shown by CMZ, can be formulated in the form of mapping equations. A similar and simple mapping equation is obtained:

 
formula

by using (5.4) and the relation of the thermocline depth anomaly and SST anomaly in the eastern Pacific in this limit. This simple mapping model describes the essential physics of the coupled-wave oscillator formulated by CMZ who more carefully considered the shallow water ocean dynamics in the coupled model. Similar mapping equations were vigorously derived and discussed in Schopf and Suarez (1990). The dispersion relation of (5.12) can be obtained by replacing σ + cE with cE in (5.5), and the behavior of the leading mode is very similar to that shown in Fig. 3c. When the reflections at both the western and eastern boundaries are strong or rWrE is close to 1, the oscillatory ocean wave modes are less damped and can be easily destabilized through the coupling. The oscillating mechanism is inherited from the uncoupled free mode which depends on the two-way propagation of the equatorial waves in the ocean basin. In this case, Sverdrup balance is not valid even in the equatorial strip, which is consistent with the fact that the coupled wave mode at the fast wave limit tends to be a high-frequency mode.

In the system (5.12), the wind stress is directly related to the thermocline depth at the eastern Pacific because the SST anomalies of the eastern Pacific are slaved by the thermocline depth. When the eastern boundary reflection is ignored (rE = 0), this wind stress, on one hand, continuously generates the forced Kelvin wave signals [the second term in the right-hand side of (5.12)] that enhance the thermocline depth anomaly and thus act as a positive feedback. On the other hand, the wind stress also generates the forced Rossby wave signals that are reflected as Kelvin wave signals [the last term in the right-hand side of (5.12)] at the western boundary. The reflected Kelvin wave signal tends to reverse the thermocline depth anomaly in the eastern Pacific, but this negative feedback is retarded due to wave propagation processes. Thus in this case, the coupled wave oscillator shares the original interpretation proposed by BH for the BH delayed oscillator.

As shown in CMZ, the mapping model (5.12) can also be approximately reduced to the form of the BH delay oscillator by setting rE = 0 and taking the Taylor expansion for the second term on the right-hand side of (5.12):

 
formula

Now, by considering the zonal mean thermocline depth in the same manner as in the previous section for the other version of the BH delayed oscillator, one can explain the oscillation mechanism of (5.13) with the recharge oscillator mechanism as well. However, at the fast-SST limit, the period of the coupled wave mode is about a year at Hopf bifurcation. It was shown by CMZ that the inclusion of higher-order Rossby wave modes results in a much longer period through a nonlinear period doubling route. As shown in JN, the inclusion of high-order Rossby waves allows the lower frequency ocean adjustment mode to be merged with the SST mode to produce a low-frequency mixed SST–ocean-dynamics mode. In the two-strip model, the zero-frequency ocean adjustment mode can be viewed as the representation of the low-frequency ocean adjustment modes, and its merging with the SST mode gives rise to the slow recharge oscillation. Thus the scenario of producing a slow coupled wave mode at the fast-SST limit through period doubling is an unnecessary conceptual complication for the understanding of the slow ENSO variation.

The coupled-wave mode and the mixed SST–ocean-dynamics mode of the two strip-model exist in different parameter regimes and they merge to form a generalized mixed mode that leads to the recharge oscillation. The coupled wave oscillator and the BH delayed oscillator are special approximations to the generalized mixed mode. They all can be physically interpreted by the recharge oscillation mechanism. This discussion may resolve the apparent conceptual differences among the different prototype models.

6. Summary

A two-strip version of a ZC-type coupled model is constructed as a framework to describe the slow coupled instability of the ocean–atmosphere system of the tropical Pacific. The equatorial SST anomaly is governed by a linear thermodynamical equation dominated by thermocline feedback processes. Under the assumptions that wind stress forcing is hemispherically symmetric and is of broad meridional scale relative to the oceanic Rossby radius, a two-strip approximation to the equatorial shallow water dynamics leads to a simple model for the equatorial thermocline variation near the equatorial wave guide. The two-strip ocean model describes the role of the equatorial Kelvin and Rossby waves in the establishment of the quasi-equilibrium equatorial Sverdrup balance and the slow recharge/discharge of the equatorial heat content. The latter results from both the interior meridional transport and boundary layer transports by western and eastern boundary currents. The interior meridional transport is carried out by the forced and reflected Rossby waves. The boundary layer transports are inferred from the normal mass fluxes that can be related to inexact cancellation of normal velocities of incoming and reflected waves at the west and east boundaries and parameterized by reflection parameters. The two-strip ocean model and the one-strip SST equation are coupled to form a stripped-down coupled model through the relation between the equatorial SST anomaly and wind stress anomaly.

The two-strip model based on the equatorial wave dynamics incorporates the mechanisms of BWCZ and can be reduced to the conceptual recharge oscillator model with either a two-box spatial truncation or a low-frequency approximation. Both of these two approximations yield the conceptual model because both filter all the relative high-frequency ocean adjustment modes and retain only the zero-frequency ocean adjustment mode. Moreover, the results from the two-strip model and its two-box approximation, including the regime behavior, period, phase trajectory, evolution of the SST, and thermocline depth anomalies of the oscillatory solution, are all in good agreement. Therefore a solid dynamical ground is provided to the conceptual recharge oscillator model intuitively constructed based on BWCZ mechanisms in Part I.

The two-strip coupled model also allows an analytical approach to its eigen solutions. There are two kinds of oscillatory modes in this system. One is the mixed SST–ocean-dynamics mode that comes from the merger of a decaying SST mode and a decaying ocean-adjustment mode. Another is a coupled wave mode, which is a destabilized gravest oscillatory ocean wave mode. The coupled-wave mode exists at the fast-wave limit although its period tends to be short at this limit. These two modes are further mixed in the regime away from both the fast-wave and the fast-SST limits. The regime behavior of the two-strip model analyzed analytically agrees with that obtained numerically by JN in a more complex model.

The reflection parameters in the two-strip model indicate the role of mass fluxes through eastern and western boundary layer currents in the recharging/discharging of the equatorial heat content. The dependence of the critical coupling coefficient and frequency at the critical coupling on these parameters reveals that both of them have significant impacts on the recharge oscillator. In particular, a strong eastern boundary reflection tends to reduce the frequency of the recharge oscillator and a strong western boundary reflection tends to produce a leading mode with a relative high frequency and strong instability. The leading mode for the recharge oscillator does exist even if the eastern boundary reflection parameter is set to zero, but it disappears if the reflection parameter of the western boundary is small. All these dependencies can be understood from the point of the equatorial mass recharge/discharge mechanism.

The collective damping rate cE in the linear SST equation determines whether a recharge oscillator exhibits its behavior close to the coupled wave oscillator at the fast-wave limit or close to the simple mixed mode far away from both the fast-SST limit and the fast-wave limit. This damping rate can be much reduced, for instance, by the feedback from horizontal and vertical advective feedback through Ekman flow. The weaker the collective damping rate, the less accurate is the fast-SST limit approximation. The location of the wind patch in response to the SST anomaly also has some quantitative impact on the stability of the leading mode. When the wind patch is located in the western Pacific, the system needs a relatively strong coupling to reach the Hopf bifurcation for a self-sustained oscillation. The coupled mode leading to the recharge oscillation is also more likely to be originally from the coupled wave mode. When the wind patch is located at the central to eastern Pacific, the system needs only a moderate coupling to reach instability and the slow mixed mode is the leading mode. Nevertheless, the coupled system is much more sensitive to the damping rate cE than to the location of the wind patch.

The BH delayed oscillator can be viewed as the special case of the generalized mixed mode with no eastern boundary reflection. It can be physically interpreted by the recharge oscillation mechanism proposed by BWCZ. The form of the BH delayed oscillator, the mapping equation for the CMZ coupled wave oscillator, and the conceptual recharge oscillator model can all be reduced from the two-strip model. It becomes clear that all these prototype models share the same basic BWCZ mechanism. In fact, the BH delayed oscillator can virtually be fitted from either the coupled wave oscillator or the mixed SST–ocean dynamics oscillator near the Hopf bifurcation. Thus it can be viewed either as a particular version of the recharge oscillator model or as a fitting model. The conceptual recharge oscillator model and the analytical two-strip model provide the simplest and most constructive way to approach the mixed mode. In particular, the picture of the mixed mode resulting from a series of mergers of oceanic modes with an SST mode in JN is simplified with the conceptual model to the single merger of a decaying ocean adjustment mode with an SST mode. This simplification reveals the recharge-oscillation physics of the mixed mode. Together with its generalization to the delayed oscillator and the coupled wave oscillator at the fast-wave limit, the recharge oscillator has the advantage of being a simple paradigm of ENSO.

It should be pointed out that in more complicated coupled models, different oscillatory modes that exist at different parameter regimes of this model may coexist under the same parameter setting. For instance, in the ZC model, there are a high-frequency coupled mode with period of less than a year and a low ENSO mode of 3–4 yr period (Mantua and Battisti 1995). The basic features of the slow mode of ZC can be described by the recharge oscillator. The coexistence of these modes is due to the consideration of full SST thermodynamics in the ZC model, which allows SST anomalies to develop through anomalous zonal advection in the central Pacific and through the thermocline and Ekman feedback processes in the eastern Pacific. Further understanding of the rich behavior of ENSO lies beyond the current conceptual model. Different characteristics of ENSO events in nature and in complex models may well be related to the fact that different coupled modes can be excited through the nonlinear interaction of the coupled tropical ocean–atmosphere system or by uncoupled random disturbances.

Acknowledgments

This work was supported by National Science Foundation Grant ATM-9312888 and by National Oceanographic and Atmospheric Administration Grant GC95773. The author is grateful to Mark Cane for his detailed and constructive comments on the manuscript to bring this paper into shape. He enjoyed the stimulating discussions with Bin Wang and Roger Lukas. The author thanks Eli Tziperman and Edward Schneider for their critiques that led to the improvement of the paper, Thomas Schroeder and Diane Henderson for their careful reading and editing of the manuscript, and Xiaowei Sun for producing some of the figures. This is SOEST Contribution Number 4171.

REFERENCES

REFERENCES
Battisti, D. S., 1988: The dynamics and thermodynamics of a warming event in a coupled tropical atmosphere–ocean model. J. Atmos. Sci., 45, 2889–2919.
———, and A. C. Hirst, 1989: Interannual variability in a tropical atmosphere–ocean model: Influence of the basic state, ocean geometry, and nonlinearity. J. Atmos. Sci., 46, 1687–1712.
Bjerknes, J., 1969: Atmospheric teleconnections from the equatorial Pacific. Mon. Wea. Rev., 97, 163–172.
Cane, M. A., and D. W. Moore, 1981: A note on low-frequency equatorial basin modes. J. Phys. Oceanogr., 11, 1578–1584.
———, and E. S. Sarachik, 1981: The periodic response of a linear baroclinic equatorial ocean. J. Mar. Res., 39, 651–693.
———, and S. E. Zebiak, 1985: A theory for El Niño and the Southern Oscillation. Science, 228, 1084–1087.
———, M. Münnich, and S. E. Zebiak, 1990: A study of self-excited oscillations of the tropical ocean–atmosphere system. Part I: Linear analysis. J. Atmos. Sci., 47, 1562–1577.
Deser, C., and J. M. Wallace, 1990: Large-scale atmospheric circulation features of warm and cold episodes in the tropical Pacific. J. Climate, 3, 1254–1281.
Gill, A. E., 1980: Some simple solutions for heat induced tropical circulation. Quart. J. Roy. Meteor. Soc., 106, 447–462.
Graham, N. E., and W. B. White, 1988: The El Niño cycle: Pacific ocean–atmosphere system. Science, 240, 1293–1302.
Hirst, A. C., 1986: Unstable and damped equatorial modes in simple coupled ocean–atmosphere models. J. Atmos. Sci., 43, 606–630.
———, 1988: Slow instabilities in tropical ocean basin–global atmosphere models. J. Atmos. Sci., 45, 830–852.
Jin, F.-F., and J. D. Neelin, 1993a: Modes of interannual tropical ocean–atmosphere interaction—A unified view. Part I: Numerical results. J. Atmos. Sci., 50, 3477–3502.
———, and ———, 1993b: Modes of interannual tropical ocean–atmosphere interaction—A unified view. Part III: Analytical results in fully coupled cases. J. Atmos. Sci., 50, 3523–3540.
———, ———, and M. Ghil, 1996: The interaction of ENSO and the annual cycle: Subharmonic frequency-locking and ENSO aperiodicity. Physica D, in press.
Kessler, W. S., and M. J. McPhaden, 1995: Oceanic equatorial waves and the 1991–93 El Niño. J. Climate, 8, 1757–1774.
Kubicek, M., and M. Marek, 1983: Computational Methods in Bifurcation Theory and Dissipative Structures. Springer-Verlag, 243 pp.
Latif, M., and A. Villwock, 1990: Internannual variability as simulated in coupled ocean–atmosphere models. J. Mar. Syst., 1, 51–60.
Lindzen, R. S., and S. Nigam, 1987: On the role of sea surface temperature gradients in forcing low-level winds and convergence in the Tropics. J. Atmos. Sci., 44, 2418–2436.
Mantua, N. J., and D. S. Battisti, 1994: Evidence for the delay oscillator mechanism for ENSO: The “observed” oceanic Kelvin mode in the far western Pacific. J. Phys. Oceanogr., 24, 691–699.
McWilliams, J., and P. Gent, 1978: A coupled air–sea model for the tropical Pacific. J. Atmos. Sci., 35, 962–989.
Münnich, M., M. Cane, and S. Zebiak, 1991: A study of self-excited oscillations of the tropical ocean–atmosphere system. Part II: Nonlinear cases. J. Atmos. Sci., 48, 1238–1248.
Neelin, J. D., 1989: On the interpretation of the Gill model. J. Atmos. Sci., 46, 2466–2468.
———, 1990: A hybrid coupled general circulation model for El Niño studies. J. Atmos. Sci., 47, 674–693.
———, 1991: The slow sea surface temperature mode and the fast-wave limit: Analytic theory for tropical interannual oscillations and experiments in a hybrid coupled model. J. Atmos. Sci., 48, 584–606.
———, and I. M. Held, 1987: Modeling tropical convergence based on the moist static energy budget. Mon. Wea. Rev., 115, 3–12.
———, and F.-F. Jin, 1993: Modes of interannual tropical ocean–atmosphere interaction—A unified view. Part II: Analytical results in the weak-coupling limit. J. Atmos. Sci., 50, 3504–3522.
———, M. Latif, and F.-F. Jin, 1994: Dynamics of coupled ocean-atmosphere models: The tropical problem. Annu. Rev. Fluid Mech. 26, 617–659.
Philander, S. G. H., T. Yamagata, and R. C. Pacanowski, 1984: Unstable air–sea interactions in the Tropics. J. Atmos. Sci., 41, 604–613.
———, R. C. Pacanowski, N. C. Lau, and M. J. Nath, 1992: A simulation of the Southern Oscillation with a global atmospheric GCM coupled to a high-resolution, tropical Pacific ocean GCM. J. Climate, 5, 308–329.
Schneider, E. K., B. Huang, and J. Shukla, 1995: Ocean wave dynamics and El Niño. J. Climate, 8, 2415–2439.
Schopf, P. S., and M. J. Suarez, 1988: Vacillations in a coupled ocean–atmosphere model. J. Atmos. Sci., 45, 549–566.
———, and ———, 1990: Ocean wave dynamics and the time scale of ENSO. J. Phys. Oceanogr., 20, 629–645.
Suarez, M. J., and P. S. Schopf, 1988: A delayed action oscillator for ENSO. J. Atmos. Sci., 45, 3283–3287.
Wakata, Y., and E. S. Sarachik, 1991: Unstable coupled atmosphere–ocean basin modes in the presence of a spatially varying basic state. J. Atmos. Sci., 48, 2060–2077.
Wyrtki, K., 1975: El Niño—The dynamic response of the equatorial Pacific Ocean to atmospheric forcing. J. Phys. Oceanogr., 5, 572–584.
———, 1986: Water displacements in the Pacific and the genesis of El Niño cycles. J. Geophys. Res., 91, 7129–7132.
Zebiak, S. E., 1986: Atmospheric convergence feedback in a simple model for El Niño. Mon. Wea. Rev., 114, 1263–1271.
———, and M. A. Cane, 1987: A model El Niño–Southern Oscillation. Mon. Wea. Rev., 115, 2262–2278.

APPENDIX A

The Method of Tracing Eigenvalues, Their Degeneracies, and the Hopf Bifurcation

The dispersion relation (4.4) can be more generally written as

 
F(σ, μ) = 0,
(A.1)

where σ is the eigenvalue, μ is a vector of parameters with dimension k, and F is a smooth scalar function. The roots of (A.1) give the dependence of the eigenvalues on the parameters σ(μ). Since eigenvalues can be complex, one can separate (A1) into real and imaginary parts with real functions of real variables:

 
formula

Solving eigenvalues of (A.1) becomes a problem of finding zeros in a two-dimensional variable space (σr,σi) and their dependence in the k-dimension parameter space. It can be numerically approached with a standard continuation method (e.g., Kubicek and Marek 1983).

Similarly, one can trace an algebraic multiplicity-two degeneracy occurring at the merging of two real eigenvalues so that at one side in parameter space there are two real modes and at the other side they become a pair of complex modes. This situation happens when ∂F /∂σ = 0. However, for a real mode with σi = 0, F(σr, μ) is a real function and the second equation of (A.2) is a zero identity. Thus solutions at this multiplicity-two degeneracy satisfy the following equations:

 
formula

The second equation poses one additional constraint so that solutions exist on a k−1 dimensional surface in the k-dimensional parameter space (e.g., JN). For the case of multiplicity-three degeneracy of real modes, one adds another constraint, ∂2F/∂2σr = 0, and solutions exist on a (k−2)-dimensional surface. On the other hand, to trace the complex multiplicity-two degeneracy, in addition to the set of (A2), one has two additional constraints,

 
formula

Thus this degeneracy occurs on a (k − 2)-dimensional surface in the parameter space.

One can also trace the stability boundary of the leading complex mode or the boundary in the parameter space for the first Hopf bifurcation. Let σr = 0 in (A.2), so that the bifurcation point satisfies the following equation:

 
Fr(0,σi,μ) = Re(F ), Fi(0,σi, μ) = Im(F ).  
(A.5)

Thus it occurs on a k –1 dimensional surface of the parameter space because in this case the critical frequency is constrained by two equations. All these different equation sets, (A.2), (A.3), (A.4), and (A.5), can be solved with a standard continuation method.

Tracing the degeneracies and the Hopf bifurcation is particularly useful for understanding the topology of the modes of the system in the parameter space. It provides a succinct presentation of the behavior regimes in the parameter space to indicate the dynamical origins of the modes and to reveal the physical processes that lead to their qualitative changes.

APPENDIX B

Energy Leaking at the Western Boundary

The zero integrated mass flux condition at the western boundary in (2.3) allows energy leaking in the system (2.1) and thus results in nonconservation of total energy even without any forcing and explicit damping. Following the nondimensionalization for (2.1) in section 2, and omitting the wind forcing and the damping, one can get the total energy (kinetic energy and available potential energy) equation as

 
formula

where h0 = h(0,y,t ) and h00 = h(0,0,t ). The boundary conditions in (2.3) and υh = 0 and (h0h00)2/y = 0 at y = ±∞ have been used. The total energy change rate in (B.1) is negative definite for any distribution of h, except for the case with h0 = h00 at which normal flow is zero at all latitudes of the western boundary. Therefore, the western boundary under condition (2.3) is energy leaky, which results from the filtering out of short dispersive Rossby waves and gravity waves under the long wave approximation. The case of h0 = h00 corresponds to the neutral ocean-basin modes discovered by Cane and Moore (1981). The energy leaking at the western boundary is responsible for the decay of scattering modes in the remainder spectrum (Neelin and Jin 1993) of the adjustment dynamics of (2.1).

Footnotes

Corresponding author address: Dr. Fei-Fei Jin, Dept. of Meteorology, School of Ocean and Earth Science and Technology, University of Hawaii at Manoa, 2525 Correa Road, HIG 331, Honolulu, HI 96822.