## Abstract

The equatorial tropical Pacific climate system is a delicate coupled system in which winds driven by gradients of sea surface temperature (SST) within the basin interact with the ocean circulation to maintain SST gradients. This results in a time mean state having a strong zonal temperature contrast along the equator with an eastern cold tongue and a western warm pool. By the same coupled processes, interannual variability, known as the El Niño–Southern Oscillation (ENSO), is present in the Pacific. This variability can be attributed to an oscillatory coupled mode, the ENSO mode, in the equatorial ocean–atmosphere system. Using a Zebiak–Cane-type intermediate coupled model, the coexistence of an eastern cold tongue in the annual mean state and ENSO in the Pacific climate system is investigated. The ENSO mode arises as a robust oscillatory mode on a coupled mean state and becomes unstable if the cold tongue of the mean state is sufficiently strong. The origin of this mode, its propagation mechanism, its sensitivity to parameters, and its relation to the spatial structure of the annual mean state are considered.

## 1. Introduction

One of the first models by which the interannual variability in the tropical Pacific, known as El Niño–Southern Oscillation (ENSO), was reasonably simulated was that of Zebiak and Cane (1987). An annual mean state or seasonal cycle of both ocean and atmosphere is obtained from data, and within the model the evolution of anomalies with respect to this reference state is computed. This coupled model produces recurring ENSO events, as a basinwide variation in the coupled ocean–atmosphere system, that are irregular in both amplitude and spacing but favor a 3–4-yr period. The Zebiak–Cane model and its variants (Battisti 1988; Neelin 1991) turn out to be very well suited to determine the details of the evolution of ENSO events and have been widely used to study their origin and physics. Over the past decade, ENSO theory has advanced to a relatively mature stage (Neelin et al. 1998).

For a spatially constant annual mean state, Hirst (1986, 1988) found that oceanic equatorial modes can be destabilized by coupled processes. Further work into the nonlinear equilibration of these instabilities (Battisti and Hirst 1989; Schopf and Suarez 1988; Suarez and Schopf 1988) led to the BHSS delayed action oscillator paradigm (Schopf and Suarez 1990). In this metaphor for ENSO dynamics, only one oscillatory mode is associated with the interannual variations on the annual mean state. This mode becomes unstable at sufficiently large coupling strength—that is, the amount of wind stress per SST anomaly—and equilibrates due to nonlinear effects. Ocean adjustment processes cause a delayed negative feedback, with central importance for a reflection of Rossby waves along the western boundary. These processes can be modeled within one nonlinear differential-delay equation. Alternative simple models, for example, the coupled wave oscillator model, were derived from point-coupling models that combine the coupled feedbacks with linear wave dynamics (Cane et al. 1990; Münnich et al. 1991).

In this view of the ENSO mode, adjustment processes of sea surface temperature are relatively unimportant (the fast SST limit). However, it appeared that interannual oscillations can also occur (Neelin 1991) even when ocean wave timescales are very fast compared to those determining the SST field (the fast-wave limit). These SST (or thermal) modes (which decay when there is no ocean dynamics) are modified through ocean–atmosphere dynamics and can be destabilized through coupled processes. Hence, coupled ocean–atmosphere models may exhibit a variety of different oscillatory modes depending on the conditions set by parameters in the model.

In the uncoupled case, distinct sets of modes occur that are primarily related to the timescales of SST change (SST modes) and to the timescales of ocean adjustment (ocean dynamics modes). The connection between both classes of modes in parameter space was explored by Jin and Neelin (1993a,b) and Neelin and Jin (1993) in a Zebiak–Cane-type model using an equatorial strip approximation for the SST equation. At realistic coupling strength, the coupled modes are best described by mixed SST–ocean dynamics modes. The spatial structure of the ENSO mode is set by properties of a stationary SST mode that is unstable over a large range in parameter space at strong coupling. With decreasing coupling strength, this mode merges with oscillatory modes originating from ocean dynamics. In this way, ocean subsurface dynamics controls the period of oscillation. This standing–oscillatory regime as a standing SST mode perturbed by ocean wave dynamics provides a framework of understanding ENSO as an unstable coupled mode.

However, the picture of the mergers of the different modes is very complicated in the realistic regime (Jin and Neelin 1993a). The reason for this is the presence of resolution-dependent modes (called scatter modes) in the shallow-water model on the equatorial *β*-plane arising from discretization of an essentially continuous spectrum (Moore 1968). The classic ocean basin modes (Cane and Moore 1981) arise as peaks in this spectrum. Hence, there has been a desire to develop more simple models capturing the essentials of the ENSO mode. In a two-strip approximation of the shallow-water dynamics, Jin (1997b) shows that all adjustment due to the scatter modes is described within one stationary mode. As coupling is increased, this stationary mode merges with a stationary SST mode to give an oscillatory unstable mode. The resulting mode has the property that the zonally averaged thermocline anomaly and the eastern Pacific SST anomaly are not in phase. At the transition phase of zero SST anomaly, the average thermocline depth is shallower or deeper than its climatological value. Hence, within one cycle of the oscillation the equatorial heat content is discharged and recharged once. These features have been observed in data (Wyrtki 1975) and in models (Cane 1992) and have led to a complementary view of the ENSO mode as a recharge oscillator. Jin (1997b) shows that the ocean wave oscillator, BHSS delayed oscillator, and recharge oscillator can all be derived from the two-strip approximation using reasonable assumptions. Furthermore, each basic oscillation can be captured by a simple two-dimensional dynamical system (Jin 1997a). This unifies the different paradigms of the oscillatory nature of the ENSO mode, and it is simultaneously consistent with the Jin and Neelin (1993a) picture of the unstable modes, although the spectrum of scatter modes has been considerably reduced.

In all studies cited above, the annual mean state is constructed—for example, by using flux correction (Neelin and Dijkstra 1995)—and is independent of coupled processes governing the ENSO mode. In reality, coupled processes also determine the annual mean state (and the seasonal cycle). A relatively small zonally constant wind stress, which is thought to arise from factors outside the basin, causes a weak zonal SST gradient, which is amplified by coupling to give the zonal equatorial warm pool–cold tongue structure. This was recently referred to as the “climatological version of the Bjerknes hypothesis” (Dijkstra and Neelin 1995a). Again, an intermediate model was successfully used to demonstrate this and to identify areas in parameter space, where a fully coupled climatology can be found (Dijkstra and Neelin 1995a) with the correct east–west structure. It was also shown that this state is unique and that previously found regimes of multiple steady states are model artifacts (Neelin and Dijkstra 1995).

The stability properties of the fully coupled mean states may differ considerably from their constructed counterparts. When coupling is varied, the basic state changes simultaneously with the most unstable perturbations. The absence of multiple equilibria, that is, transcritical bifurcations, forbids the existence of a neutral stationary SST mode, which appears central in the Jin and Neelin (1993a) picture. Hence, one may question whether the ENSO mode is also strongly modified in the fully coupled case. This question was first addressed by Jin (1996) using the simple (box) recharge oscillator model. In this fully coupled model, oscillatory instabilities still persist and have ENSO characteristics. The stability of coupled annual mean states has been investigated in the Jin and Neelin (1993a) model by Dijkstra and Neelin (1999). In this case, the uncoupled ocean dynamics spectrum is filtered such that only ocean basin modes remain. It turns out that a mixed SST–ocean dynamics mode is able to destabilize, but both its structure and period appear very sensitive to parameters. Apparently, one can no longer easily destabilize this mode, since in the regime where this is possible, the climatology has not the correct spatial structure. This is in contrast with the results in models with a constructed mean state, where this regime can always be found, because the mean state is a solution for all values of the coupling strength.

In this paper, the stability of fully coupled climatologies is studied within a model that closely approximates the Zebiak–Cane model. Simultaneous sensitivity of the annual mean state and ENSO mode is determined by using techniques of numerical bifurcation theory (Dijkstra and Neelin 1995b). The model is briefly described in section 2. In section 3, it is shown that ENSO-like unstable modes exist on a fully coupled climatology that has a reasonable warm pool–cold tongue structure. The stability of this climatology is investigated, and modes are traced in parameter space to clarify their origin. In section 4, the robustness of these results is examined, and a topological view of the ENSO mode is presented in section 5, followed by the last section, where the results are summarized.

## 2. The fully coupled Zebiak–Cane model

The formulation of the Zebiak–Cane model is well described elsewhere (Zebiak and Cane 1987). The ocean model consists of a shallow-water layer of mean depth *H* with an embedded mixed layer of fixed depth *H*_{1} in an ocean basin that has zonal length *L* and is meridionally unbounded. The flow in the shallow-water layer with velocity field (*u, υ, w*) and thermocline displacement *h* is determined by a reduced-gravity model in the long-wave limit and forced by a wind stress vector ** τ** = (

*τ*

^{x},

*τ*

^{y}). Friction in the shallow-water component is linear with damping coefficient

*r**. The zonal velocity vanishes at the eastern boundary, and a zero mass flux condition is applied at the western boundary; all quantities are bounded far from the equator. An Ekman layer model with linear damping coefficient

*ε*

^{*}

_{s}, is driven by the wind stress

*τ*and determines the Ekman velocities (

*u*

_{s},

*υ*

_{s}).

Ocean processes affect the sea surface temperature in a way modeled by the SST equation:

where *u*_{1} = *u*_{s} + *u* and *υ*_{1} = *υ*_{s} + *υ* are the horizontal velocities in the mixed layer, *w*_{1} = *w*_{s} + *w* the vertical velocity just below the mixed layer, and ℋ is the Heaviside function. The quantity *T*_{0} is the radiation equilibrium temperature, which is the temperature that is attained in absence of ocean dynamics; Newtonian cooling is represented by the coefficient *ε*_{w}. The subsurface temperature *T*_{s}(*h*) depends on *h* according to

where *h*_{0} and *H*∗ control the steepness and the offset of the *T*_{s} profile, and *T*_{so} is the characteristic temperature being upwelled into the surface layer.

The ocean model is coupled to a Gill (1980) atmosphere model with horizontal velocities (*u*_{a}, *υ*_{a}), geopotential height *ϕ,* and with linear damping coefficient *A**. The atmosphere is driven by heat fluxes from the ocean that depend linearly on the anomalies of sea surface temperature *T* with respect to the radiation equilibrium temperature *T*_{0}, with proportionality constant *α*_{T}.

Under the assumptions that (i) there exists a horizontally uniform radiative–convective equilibrium state and that (ii) in the Tropics the radiative and convective processes tend to bring the atmosphere to a thermodynamical equilibrium—whereas the dynamics is responsible for the departure of this equilibrium—the simple heating parameterization and the Gill-type dynamical framework provide the zeroth-order description of dynamics and thermodynamics of the tropical atmosphere. In that sense, all components of heating are condensed into the parameter measuring the restoration of the radiative–convective equilibrium. This is a very crude parameterization. Nevertheless, considering other simplifications in the coupled model, it is not very productive to employ very detailed parameterization schemes of the heating. It should be pointed out that the heat flux into the ocean and the heating in the atmosphere are consistent, in the sense that without dynamical processes in the atmosphere and ocean, the SST is also in the radiative–convective equilibrium state. The interpretation of heating rate is different here from the conventional models because we have assumed a uniform radiative convective equilibrium. In our formulation of the Gill model, there is little departure, in terms of heating, from the heating rate of the equilibrium state in the convective regions (e.g., western Pacific), but there is a net heating deficit in the nonconvective regions (e.g., cold tongue region). If one takes this into consideration, there is no inconsistency, at least qualitatively, with the observations.

Wind stresses are related linearly to surface velocities with proportionality constant *γ*_{τ}. As in Neelin and Dijkstra (1995), the dimensionless zonal wind stress forcing *τ*^{x} is decomposed into

and the meridional wind stress component *τ*^{y} is neglected. As discussed in detail in Dijkstra and Neelin (1995a), the easterly wind stress, *τ*^{x}_{ext}, is interpreted as that part of the wind stress that does not depend on the coupled feedbacks within the basin. Even in the absence of longitudinal SST gradients, there would still be a Hadley and a Walker circulation, driven by zonally symmetric latitudinal gradients and zonally asymmetric land–sea contrast, respectively. Its inclusion is necessary, because the Gill model has deficiencies in simulating the zonal mean part of the atmospheric flow. The part of the wind stress due to coupling *τ*^{x}_{c} is proportional to the coupling strength.

For the forcing of the Ekman layer, an extra parameter *δ*_{s} is introduced,

The parameter *δ*_{s} measures the effect of the Ekman velocities on the SST evolution through coupled feedbacks, as in Jin and Neelin (1993a).

Finally, we only consider the meridionally symmetric coupled system in this paper. The processes causing equatorial asymmetry are complicated, and we consider their inclusion as a next step in a theory to understand mean state, seasonal cycle, and ENSO in one framework, which is outside the scope of this paper. With these limitations, we still think that the coupled model used here captures the zeroth-order physics of the interaction between the mean state and ENSO mode insofar that coupled processes have impact on both.

The governing equations are nondimensionalized using scales for zonal length *L,* meridional length *L*_{y}, vertical length *H,* wind stress *τ*_{0}, time *L*/*c*_{0}, horizontal velocity *c*_{0}, meridional velocity *c*_{0}*L*_{y}/*L,* vertical velocity *c*_{0}*H*_{1}/*L,* thermocline depth *H*_{1}, and temperature Δ*T.* Here, *c*_{0} = *g*′*H* which is the phase speed of the first oceanic baroclinic Kelvin mode, and *L*_{y} is the oceanic internal Rossby deformation radius. For the atmosphere, the meridional length scale is taken as *L*_{a}, and both the zonal and meridional velocities are scaled with Kelvin wave velocity *c*_{a}. This introduces several dimensionless parameters, the most important of which are the coupling strength *μ* and the atmospheric damping *ε*_{a}.

The dimensionless equations are given in appendix A, and expressions for all dimensionless parameters in the equations in terms of dimensional parameters and their standard values are given in Table 1. In comparison to the original parameter values in Zebiak and Cane (1987), the most important differences lie in the *T*_{sub} parameterization and the reference phase speeds of both the oceanic (*c*_{o}) and atmospheric (*c*_{a}) Kelvin wave, with *c*_{o} = 2 m s^{−1} and *c*_{a} = 30 m s^{−1} versus *c*_{o} = 2.9 m s^{−1} and *c*_{a} = 60 m s^{−1}. The standard value of *ε*_{a} used in the results is 1.25 and corresponds to an atmospheric damping timescale of about 5 days, which gives a smaller damping rate than that used in most models. Likely due to the deficiencies of the atmosphere model, we found best-looking results at this particular value of *ε*_{a}. The sensitivity of the model results to variations in the damping rate will be investigated in section 4c.

## 3. The warm pool–cold tongue and its bifurcation to ENSO

### a. The fully coupled steady state

At zero coupling (*μ* = 0), the ocean circulation and consequently SST are determined by the external zonal wind stress *τ*^{x}_{ext}. The amplitude and shape of *τ*^{x}_{ext} are hard to determine from observations. Although GCM runs to get an impression of these quantities have been proposed (Dijkstra and Neelin 1995a), there is no a priori justification for any assumed shape and amplitude.

This wind stress is assumed to have the form

where *α* is the ratio of Rossby deformation radii between ocean and atmosphere, which controls the meridional extension of the external wind. The amplitude *F*_{0} has a standard value of 0.2, corresponding to a dimensional value of 0.01 Pa. This rather large value of *F*_{0} is necessary to generate sufficient background upwelling in the eastern Pacific basin at the equator. This is due to the equatorially symmetric model setup and deficiencies in the atmospheric model. At each latitude, the external wind is zonally constant. All results below are for *F*_{2} = 0 corresponding to maximum trade winds at the equator. Experiments with (*F*_{0}, *F*_{2}) = (0.05, 0.1), which mimic the off-equatorial maxima of the trades near 10°*N,* showed no qualitatively different results.

The resulting temperature field at zero coupling, *T*_{ext}, and thermocline depth *h*_{ext} are shown in Fig. 1. In response to the easterlies, the temperature increases monotonically from about 25.5°C in the east to about 28.5°C in the west along the equator. The equatorial thermocline response is approximately linear, with slight off-equatorial maxima near the western boundary.

As the coupling *μ* is increased, the temperature of the Pacific cold tongue decreases due to coupled feedbacks. At low coupling, the additional wind stress due to coupling is approximately the atmosphere response to the original cooling *T*_{ext} − *T*_{0}, which enhances the external easterlies over most of the basin, leading to larger upwelling and a stronger thermocline slope (Dijkstra and Neelin 1995a). The equatorial temperature of the cold tongue (*T* − *T*_{0})_{EC} and the vertical velocity, *w*_{E}, just below the mixed layer, both at (*x* = 0.8, *y* = 0) (shown in Fig. 2), demonstrate that there is a unique steady solution as a function of *μ,* with a colder eastern Pacific and more upwelling as coupling gets stronger.

At *μ* = 0.5, the spatial structure of the coupled climatology is shown in Fig. 3. The zonal scale of the cold tongue (Fig. 3a) is set by a delicate balance of thermocline and surface layer feedbacks (Dijkstra and Neelin 1995a). The meridional extent of the cold tongue is determined both by the Ekman spreading length (*ε*^{*}_{s}/*β*) and by meridional advection. The thermocline field (Fig. 3b) is in Sverdrup balance with the wind stress giving the off-equatorial maxima and a deeper (shallower) equatorial thermocline in the west (east). This indicates that the reservoir of heat content lies off-equatorial in the central and western part of the basin. The wind stress response *u*_{a} (Fig. 3c) shows the intensification of the easterlies, with its maximum west of the cold tongue. This zonal spatial phase difference between SST and *u*_{a} is controlled by *ε*_{a}; its meridional scale is due to the atmospheric internal Rossby deformation radius, controlled by the parameter *α.* The vertical velocity structure (Fig. 3d) is clearly controlled by Ekman pumping, which gives the positive values restricted to an equatorial zone with maximum amplitude in the eastern part of the basin.

Although this mean state is a steady state of the coupled system, it still has an ad hoc nature, because of the simplicity of the model, the assumption of equatorial symmetry, and the prescribed shape of the external wind stress.

### b. Linear stability of the coupled climatology

Along a branch of steady states in Fig. 2 the linear stability is determined simultaneously, by writing the solution vector **z** (cf. appendix A) as

The vector **z** represents the climatology, **z̃** perturbations with respect to this climatology, and *σ* is the corresponding eigenvalue. In Fig. 4, the path of six modes, which become leading eigenmodes at high coupling, is plotted as a function of the coupling strength *μ* in the complex plane. In Fig. 4a, the modes are plotted in the same way as in Jin and Neelin (1993b), with a larger dot size indicating a larger value of *μ.* Both frequency and growth rate of the modes are given in yr^{−1}. In Fig. 4b only the growth rate is plotted against *μ.*

There is one stationary mode, which is the gravest mode at *μ* = 0, but its decay rate is nearly constant in *μ.* An oscillatory mode becomes unstable near *μ* = 0.503. This mode can be traced back to *μ* = 0 and appears to deform into one of the scatter modes of the shallow-water model. Although the ocean scatter modes are resolution dependent, once coupling is nonzero the character of the mode has changed significantly and is well resolved. Other modes appear hardly affected by coupling. The gravest ocean basin mode (marked with triangles) is heavily damped as *μ* is increased (Fig. 4b). The SST modes remain damped, and none of these modes was found to be interacting with the ocean dynamics modes as in the case considered in Jin and Neelin (1993a). The SST modes reside in the spectrum to the left of the modes shown in Fig. 4a. More details on the complete spectrum can be found in appendix B.

Near *μ* = 0.503, for which the mean state was shown in Fig. 3, time–longitude diagrams of the equatorial thermocline, temperature, and zonal wind anomalies of the most unstable mode are shown in Fig. 5. Note that the amplitude of the oscillation is not determined by linear stability analysis. The maximum amplitudes in the captions are relative magnitudes of the different fields; that is, the mode displays a thermocline deviation of about 10 m per degree SST anomaly.

The SST pattern shows a nearly standing oscillation for which the spatial scale is confined to the cold tongue position of the mean state. There is slight eastward propagation of the SST anomaly in the central equatorial Pacific. The thermocline anomaly shows western anomalies in heat content leading the eastern ones with the separation between positive and negative anomalies occurring within the front defined by the cold tongue; it is out of phase with the SST anomaly with a lag of about 5 months. The wind response is much broader zonally and is in phase with the SST anomaly. Although the wind maximum in the central Pacific is more to the east with respect to observations, all fields do correspond qualitatively to those observed (Latif et al. 1993;Neelin et al. 1994).

The phase relationships between wind–SST anomalies and thermocline anomalies are investigated in more detail (Fig. 6) by plotting the thermocline depth in the western Pacific *h*_{W} as well as the zonal mean equatorial thermocline displacement *h*_{ZM} against the SST anomaly in the eastern Pacific (at *x* = 0.92) over one cycle of the oscillation. At this particular location, the temperature anomaly (Fig. 5) has largest amplitude. Figure 6a shows the characteristic ENSO phase relationship between SST and thermocline anomalies with a relatively shallow (deep) thermocline in the west in the case of a warm (cold) event.

As the closed curve is traversed clockwise over one cycle of the oscillation, one clearly observes that a cold event is followed by an extremum in western thermocline anomaly. As the SST anomaly becomes zero, the western thermocline anomaly is positive in the west, and Fig. 6b shows that this also holds for the zonally averaged thermocline anomaly. Hence, the equatorial heat content is slowly built up after the cold event by the increase of the trade winds. This sets the stage for the following warm event in which the equatorial heat content is discharged. After the warm event, the zonally averaged thermocline anomaly is negative as the SST anomaly goes through zero again and the equatorial heat content is low, which causes the next cold event. Clearly, the phase relationships between the different fields of the ENSO mode are in qualitative agreement with those found in ENSO modes in Zebiak–Cane-type models (Zebiak and Cane 1987; Battisti 1988; Kirtman 1997;Schneider et al. 1995).

### c. Recharge oscillation mechanism of the ENSO mode

In this model the meridional extent of the eigenmode of the coupled (intermediate) model is well resolved. This meridional structure of the ENSO mode appears important for its period, and therefore the propagation mechanism of the oscillation is described in detail. This is done along with plots of the different fields at several phases of the oscillation relative to the period; that is, phase *t* = ^{1}/_{2} indicates the fields after half a period. In the Figs. 7–9, these plots are given for the temperature, thermocline, and zonal wind field, respectively.

The starting point of the description of the self-sustained oscillation is a positive SST anomaly (SSTA) in the eastern Pacific (early El Niño phase), as shown in Fig. 7 at *t* = 0. Westerly zonal wind anomalies to the west of the maximum in SSTA (Fig. 7) are present as can be seen in Fig. 9 at *t* = 0. The wind response amplifies the positive SSTA through Bjerknes’s feedback (*t* = ^{1}/_{16} to ^{1}/_{8}). The extent of the SST anomaly is controlled by the climatological shape of the cold tongue in both the zonal and meridional direction (cf. Fig. 3).

The equatorial thermocline response to the weaker trade winds up to *t* = ^{1}/_{8} results in a negative anomaly (i.e., negative heat content) in the western Pacific (Fig. 8, *t* = ^{1}/_{8}), which is at its minimum a few months later than the maximum of equatorial SST. As long as the positive thermocline–SST anomaly in the eastern part of the basin does not weaken, this negative anomaly cannot be discharged. However, due to ocean wave reflections at the eastern boundary, the mass fed along the equator to the eastern basin is transformed into a collective of long Rossby waves. The position of the off-equatorial maximum in the thermocline anomaly is predominantly controlled by the ratio of Rossby deformation radii between ocean and atmosphere (as will be shown below).

At the equator, the positive thermocline anomaly and consequently the SST anomalies are weakened. This reduces the east–west SST gradient, causing the anomalous westerlies to slacken (Fig. 9, *t* = ^{1}/_{8} to ^{5}/_{16}). Termination of the El Niño phase sets in, as the western warm pool discharges its previously built-up negative heat content (*t* = ^{3}/_{8} through ^{7}/_{16} in Fig. 8). This is characteristic of the recharge oscillator view of ENSO (Jin 1997a), showing a negative zonally mean thermocline anomaly at the equator during the transition from warm to cold SST anomalies. When the trades are relaxed during each transition phase between warm and cold events, anomalous cold or warm water can be pumped into the surface layer, to establish the transition from cold to warm anomalous SST. As the thermocline rises in the east, the SST anomaly becomes negative, and through Bjerknes’ feedback its amplitude increases. The trade winds recover (Fig. 9, *t* = ^{7}/_{16}) and the positive off-equatorial thermocline anomalies propagate westward (Fig. 8, *t* = ^{3}/_{8}–^{7}/_{16}).

The evolution of the meridional structure of the anomalies in the eastern Pacific (at *x* = 0.85) indicates (Fig. 10) that the largest SST anomalies are confined to the equator and the off-equatorial anomalies are in phase and have the same sign as those at the equator but are substantially weaker. The wind anomalies have a much broader structure, and off-equatorial wind anomalies are out of phase with the equatorial ones. The thermocline anomalies display a nice parabolic shape [as, for instance, assumed in the two-strip model of Jin (1997a)]. The apparent away-equator propagation is due to the fact that Rossby wave propagation speed decreases with latitude.

Time–longitude diagrams of the thermocline anomalies and meridional geostrophic velocities at 8°N (*y* = 3.8) and the anomalous wind stress curl at that particular latitude are shown in Fig. 11. The reflected Rossby waves originating at the eastern boundary are overtaken by the forced Rossby waves (Fig. 11a). Both reflected and forced Rossby waves contribute to the recharge and discharge of the equatorial thermocline. When the wind stress curl reaches zero from a positive (negative) anomaly (Fig. 11c), the Rossby waves continue to propagate, and thus the equatorial heat content continues to discharge (recharge) (Fig. 11b). It is this interval in the ocean dynamical adjustment to a changing wind stress, that is responsible for the ENSO phase transition. The physics of the ENSO mode found in the fully coupled model is similar to that in the intermediate coupled ocean–atmosphere models (Zebiak and Cane 1987; Battisti 1988) and the didactic (toy) models, derived from these intermediate models (Battisti and Hirst 1989; Schopf and Suarez 1988; Jin 1997a). The differences are in the interpretations. Within the recharge oscillator model for ENSO, one does not need to chase the detailed wave characteristics by examining the basinwide mass redistribution associated with ocean dynamical adjustment, which is completed collectively by the waves.

## 4. Robustness

Of central importance to the results above is the simultaneous existence of a sufficiently strong east–west temperature contrast, that is, a correct warm pool–cold tongue structure and a realistic ENSO mode. In the filtered model of Dijkstra and Neelin (1999) where only ocean basin modes were allowed in the ocean dynamics spectrum, this simultaneous existence turned out to be very sensitive to parameters. For example, if the coupling strength was slightly increased from a correct situation, either the spatial pattern of the ENSO mode or the climatology very soon became unrealistic.

It will be shown below that in the fully coupled model used here, sensitivity is small. Both the spatial structure of the climatology and the ENSO mode at criticality, that is, at the primary Hopf bifurcation, appear to be a very robust feature. In the next sections, the model parameters *ε*_{a}, *α,* and *F*_{0} will be changed one at a time. The path of the Hopf bifurcation will be followed simultaneously with the climatology, which is efficiently performed by continuation methods.

### a. The meridional extent of the wind response

The path of the Hopf bifurcation point in the (*α, μ*) plane and the period of oscillation at criticality are shown in Figs. 12a,b. A larger value of *α* indicates a smaller atmospheric meridional response, and for *α* → 1.0, the atmospheric meridional scale more and more resembles that of the ocean. With increasing *α,* the spatial structure of the climatology changes to become more stable; that is, one needs a larger coupling strength to excite the ENSO mode. The ENSO period decreases from about 4 yr at *α* = 0.2 to 2 yr at *α* = 0.7. This result is in agreement with those of Kirtman (1997), who found that a broader meridional wind response leads to longer ENSO periods. For *α* = 0.7, the SST climatology at the Hopf bifurcation *μ* = 2.5 (shown in Fig. 13a) shows that the spatial structure of SST is much more confined to the equator. Compared to the case *α* = 0.2, the wind response is weaker and the thermocline in the east is deeper, leading to a warmer cold tongue.

The same equatorial confinement is seen in the SST structure of the ENSO mode at the equator (Fig. 13b), and the discharge of equatorial heat content is more gradual (Fig. 13c). The spatial pattern of temperature, thermocline and zonal wind anomaly are shown in Figs. 13d–f, at phase *t* = 0 of the oscillation. A more localized wind response prevents the higher-order Rossby modes from being excited. The consequence is that only low-order Rossby modes provide the westward-propagating signal in the ocean, in this case predominantly the *n* = 1 mode. The resulting period is the 2 yr shown in Fig. 12. At smaller *α,* higher-order modes contribute to the propagation of the thermocline anomalies westward, and the period of the oscillation increases. These results are in agreement with those of Wakata and Sarachik (1992), who varied the meridional extent of the prescribed mean equatorial upwelling in a zonally periodic ocean basin and studied the effect on the dominant modes in the coupled system.

In the entire regime, the oscillation mechanism is the same as the recharge oscillator discussed in section 3c. The sensitivity of the period to *α* can be understood in terms of the ocean adjustment timescales to different meridional structures of the wind stress. This timescale decreases with increasing *α,* because fewer higher-order Rossby modes become involved. In fact, the dependence of the period and stability on *α* as shown in Fig. 12 can be understood qualitatively using the conceptual model of Jin (1997a) by examining the dependence of the ocean adjustment timescale of the model on *α.*

### b. The strength of the external wind stress

The path of the Hopf bifurcation in the (*μ, F*_{0}) plane in Fig. 14a indicates that the period and stability do not depend much on the value of *F*_{0} if its value is larger than 0.2. For small *F*_{0}, there is not sufficient upwelling to generate a correct cold tongue (Dijkstra and Neelin 1995a), stabilizing the climatology. This incorrect spatial structure of the climatology has, also, severe consequences for the spatial pattern of the most unstable mode, although the period of the oscillation (Fig. 14b) is not affected much down to *F*_{0} = 0.1.

Too large values of *F*_{0} make the eastern basin too cold because of coupled feedbacks on the externally induced zonal SST gradient. The amplitude of the cold tongue temperature anomaly for *F*_{0} = 0.4 has increased to about 11° (Fig. 15a); it is confined more to the eastern boundary and has a broader meridional extent. This is immediately reflected in the structure of the SST anomaly of the ENSO mode (Figs. 15b,d). However, there is little impact on the thermocline anomalies (Figs. 15c,e) and, hence, on the period of the oscillation.

Fortunately, the period of the ENSO mode is only weakly dependent on the particular details of *τ*^{x}_{ext}. One needs a sufficiently large externally induced SST gradient to get a correct cold tongue, and hence *F*_{0} should not be too small. Once this has been established, the spatial pattern of the ENSO mode depends on the pattern of the cold tongue, but its period is hardly affected.

### c. Atmospheric damping length

When *ε*_{a} is varied over the interval from 0.5 to 2.5, the period of oscillation is not sensitive, but the climatology stabilizes (Fig. 16). The reason for the stabilization with increasing *ε*_{a} is again the confinement of the equatorial tongue to the east (Fig. 17a), leaving a smaller area where coupling can be effective. The mechanisms of the displacement of the cold tongue with *ε*_{a} were discussed in Dijkstra and Neelin (1995a). Basically, the wind response becomes more localized to the SST field, with maximum easterly wind closer to maximum SST anomalies, leading to a shift of negative thermocline anomalies eastward.

Sea surface temperature anomalies of the ENSO mode are therefore confined to the eastern boundary for large *ε*_{a}, as can be seen in Figs. 17b,d. The wind response of the mode is more localized (Fig. 17f), but the structure of the thermocline anomalies (Figs. 17c,e) is not affected much and hence the period is not very sensitive to *ε*_{a}.

## 5. Topological view of the coupled mode

The results in the previous section indicate that the stability of the climatology with respect to the ENSO mode is controlled by the shape of the cold tongue. On the other hand, the period of the mode seems to be determined by ocean adjustment timescales that depend on the contribution of the higher-order Rossby modes to the transport of information westward. Physically, the ENSO mode can be understood as a recharge oscillator. Topologically, the mode appears to connect to one of the ocean dynamics scatter modes in the small coupling limit, as shown in section 3b.

Using the Jin and Neelin (1993a) framework to understand the origin of the ENSO modes, albeit in the constructed climatology case and using an equatorial strip approximation for the SST, it is found that a stationary SST mode is the dominant unstable mode at large coupling. This stationary mode was argued to be connected to an SST mode in the fast-wave limit (i.e., *δ* → 0 in their Fig. 22). Hence, by connection of the eigensurfaces of both type of modes, the resulting ENSO mode was characterized as a mixed SST–ocean dynamics mode.

The appearance of the unstable stationary mode was shown to be related to a transcritical bifurcation in Hao et al. (1993) and Dijkstra and Neelin (1995b). However, these transcritical bifurcations are not robust when coupled processes determine the mean state. Consequently, the stationary SST mode cannot become neutral in the fully coupled case. The question, therefore, remains whether the ENSO mode in the fully coupled case can again be characterized as a mixed SST–ocean dynamics mode.

To get as close as possible to the Jin and Neelin (1993a) results, their parameters in the *T*_{sub} parameterization have been chosen. In comparison to the standard case as in Table 1, both qualitative and quantitative changes in structure and period of the climatology and ENSO mode are small. In Fig. 18, the path of the leading eigenvalues as a function of coupling is shown for the fully coupled case. The basic structure of the modes corresponds well to that in the standard case, with one mode arising out of a bundle of ocean dynamics modes as coupling is increased and becoming unstable at about *μ*_{H} = 0.4. At high coupling (*μ* = 0.7) the oscillatory ENSO mode breaks up into two stationary modes, like in the Jin and Neelin (1993a) case. For both modes, one with an increasing and one with a decreasing growth rate, the growth rates remain positive. For even larger values of *μ,* the two stationary modes merge again to form an oscillatory pair, before becoming stable again.

The effect of coupled processes on the mean state is investigated by taking the reverse path as the one taken in Neelin and Dijkstra (1995), which went from constructed to fully coupled climatology. The climatology as in Fig. 3, obtained at *μ* = *μ*_{H}, indicated by *T* for the SST field, is taken as a reference state. The prescribed wind stress generating this climatology, say *τ*^{x}_{obs}, is obtained from the fully coupled model as

Next, the atmosphere model operator is changed such that it is only forced by SST anomalies with respect to this reference state. The path between both model formulations is then defined by a homotopy parameter *α*_{h}, and the total zonal wind stress is given by

In the constructed climatology case (*α*_{h} = 1), the paths of the leading eigenvalues with *μ* are shown in Fig. 19. In this case, the location of the primary Hopf bifurcation is (by construction) identical to that in the fully coupled case (*α*_{h} = 0). At higher coupling, the same bifurcation to stationary modes appears, but the growth rate of the lower stationary mode becomes negative. This corresponds to a transcritical bifurcation as a secondary bifurcation of the constructed climatology, which leads generically to multiple equilibriums. This is comparable to the earlier results found in Neelin and Dijkstra (1995) in a simpler model, demonstrating the occurrence of artificial stationary multiple equilibria. For this particular parameter regime the multiple equilibria are expected to be unstable, and consequently the actual impact of the transcritical bifurcation on the limit cycle arising from the Hopf bifurcation may be small.

The oscillatory mode seems to have the same origin for both types of climatologies. To investigate this, the parameter *δ* is introduced in the shallow-water component of the model. This parameter, like in Neelin and Jin (1993), measures the relative importance of ocean adjustment and thermal processes. The relation of the modes with those in the fast-wave limit (*δ* = 0) is shown in Fig. 20 and Fig. 21 for both constructed and fully coupled cases, respectively. The coupling strength has been increased to *μ* = 0.75, and the oscillatory mode has already bifurcated into two stationary modes. Note that for both cases, the climatologies are different at this value of *μ,* because in the fully coupled case this state depends on coupling, whereas in the latter case it remains unchanged, for varying values of the coupling.

As *δ* is decreased from its standard value *δ* = 1 to the fast-wave limit *δ* = 0, first the two real eigenvalues connect again and merge into an oscillatory mode near *δ* = 0.7. As this mode stabilizes with decreasing *δ,* its frequency increases. At *δ* = 0.3, the oscillatory mode becomes neutral and then rapidly decays with decreasing *δ.* In the fast-wave limit the modes remaining in the spectrum are all related to the SST equation. However, there is no merger between the SST modes and the former oscillatory coupled mode as *δ* decreases to zero.

This view of the connection between fast-wave limit and regimes dominated by ocean basin adjustment processes is the same for both cases. The parameter *δ* does not affect the climatology, and hence the differences between Fig. 20 and Fig. 21 arise from the different climatologies. In both cases an exchange of modes occurs near *δ* = 0.2 (indicated by A in both figures), but what might at first sight be interpreted from the growth rates (Fig. 21a) as a merger of the coupled mode and SST modes, is clearly discarded as a merger when viewing the frequencies (Fig. 21b). More support for this result is found by plotting the dominant modes as a function of *δ*_{sst} (Fig. 22). The mode continues without merging to the fast SST limit (*δ*_{sst} = 0), and its period is hardly affected, except near *δ*_{sst} = 0.

Our results support the idea (Jin and Neelin 1993a) that one mode, originating from an uncoupled ocean dynamics scatter mode, which gets destabilized through coupled processes, dominates the dynamics. Its timescale is set by ocean adjustment processes, and its spatial scale is heavily controlled by the climatology, instead of by the influence of SST modes. However, mixed layer processes contribute indirectly to the spatial structure of the ENSO mode, because no realistic climatology could have been obtained without them.

Within this view of the origin of the mode, the results of Jin (1997b) and of Dijkstra and Neelin (1999) can be considered as particular limits. In the two-strip model of Jin (1997b), the ocean dynamics spectrum consists of ocean basin modes and one stationary adjustment mode representing the scattering spectrum in a heavily meridionally truncated model. To get a recharge oscillator, a merger between this adjustment mode and a stationary SST mode is necessary to guarantee the right phase relationships between zonally averaged thermocline anomalies and wind stress anomalies. The former still plays an important role to give the mode its oscillatory features, and the period of oscillation is therefore dependent on *δ.* In the analysis of Dijkstra and Neelin (1999), using a model in which the oceanic scatter modes are filtered, the final unstable mode arises through a connection of an oscillatory SST mode and an ocean dynamics mode, related to the gravest ocean basin mode at small coupling. However, it appears that all characteristics of this mode are very sensitive to parameter variations and this mode is not likely to be robust, although the parameters can be tuned such that a correct climatology and a reasonable pattern of the ENSO mode occur simultaneously.

No mergers between SST modes and ocean dynamics modes were found in the present study, and moreover, the coupled mode arising from the gravest ocean basin mode does not seem to be important to induce an ENSO mode since it is heavily damped at large coupling. A rather simple picture of the origin of the ENSO mode emerges if sufficiently spatial (meridional) resolution is provided to be able to represent the correct climatology and coupled mode. When the unstable coupled mode is traced from large to zero coupling (Fig. 4) it connects to an ocean dynamics scatter mode. The sensitivity of this coupled mode to other parameters is small once the mean state and its instabilities are determined by the same processes.

## 6. Summary and discussion

Tropical coupled general circulation models have reached the stage in which the ENSO cycle is realistically simulated within the multitude of other timescale phenomena that occur in such models (Philander et al. 1992; Terray et al. 1995). Most models are also able to simulate a mean state that has correct amplitude and spatial structure. In these detailed models, both ENSO and mean state are determined by the coupled processes in the Pacific basin. An understanding of the result in terms of basic physical processes is necessary to determine a priori changes of this variability under different conditions.

Analysis of the effects of coupled processes in intermediate complexity models, such as the Zebiak–Cane model, have provided a framework to understand the origin of oscillatory coupled modes. When the mean state is constructed, and hence a solution of the model for all values of the coupling strength, the structure of the mixed SST–ocean dynamics modes provides this framework (Jin and Neelin 1993a). When the cold tongue–warm pool structure itself is determined by coupled processes, ENSO–timescale variability is connected to one robust coupled mode in the ocean–atmosphere system. The period of this mode is, for quite a large volume in parameter space, about 4 yr, when the climatology has the correct spatial structure. Both the cold tongue steady state and the characteristics of the associated ENSO mode are very robust to variations of model parameters.

The spatial pattern of the climatology is important to set the location where amplification of disturbances can occur. Only SST anomalies in the cold tongue can be amplified due to Bjerknes’s feedback. As the cold tongue warms, the trade winds relax, the equatorial thermocline tilt reduces, and the deepening of the east Pacific thermocline further amplifies the SST warming within the cold tongue. At the same time, the position of the off-equatorial waves is just outside the cold tongue. In this way, mass is exchanged from the equatorial region toward off-equatorial latitudes. This causes the discharge of equatorial heat content, which weakens the SST anomaly and, consequently, the wind anomaly. When the SST and wind anomalies diminish, this discharge process of equatorial ocean heat content continues, leading to an anomalously shallow thermocline over the entire equator. Upwelling is able to induce negative SST anomalies and therefore induces a cold phase of ENSO. This recharge oscillation mechanism, first envisioned by Wyrtki (1985) and Cane and Zebiak (1985), and later illustrated by Jin (1996, 1997a), is clearly consistent with the ENSO mode of the fully coupled model.

The ENSO mode we find has, also, features of the coupled wave oscillator model of Cane et al. (1990), since the coupled mode exists in the fast SST limit. On the other hand, it is not the gravest ocean basin mode that gets unstable at large enough coupling, but a mode that arises from the scatter spectrum of the shallow-water model. Because of this, the coupled mode obtains the phase relations as in the recharge oscillator model as described by Jin (1997a). The simplification in spatial structure of Jin (1996, 1997a,b) requires the merger between a stationary ocean dynamics mode and an SST mode to form an oscillating mode as the recharge oscillator for ENSO. When detailed spatial structure is taken into account, it is the spatial structure set by the SST processes combined with the oscillatory oceanic scatter mode to form the recharge oscillation for ENSO. Therefore the picture portrayed in Jin and Neelin (1993a) and Jin (1996, 1997a) is reconciled.

Of course, one has to be cautious with conclusions from these simple models, because, in particular, the differences between modeled mean state and the real mean state (and their different physics) may be quite large. However, the results for this simplified case provide more support that oscillatory coupled modes are relevant in the dynamics of realistic models. When conditions are supercritical, a limit cycle (the ENSO oscillation) will be the dominant mode of internal variability. But even when this oscillatory mode is not unstable—that is, the coupling strength *μ* is below the critical value for instability—it may be excited by stochastic forcing from the atmosphere (Blanke et al. 1997). The degree of sub- or supercriticality is essential whether irregularity of the ENSO cycle can be attributed to deterministic processes, that is, interaction with the seasonal cycle (Jin et al. 1994; Tziperman et al. 1994), or whether it is stochastically determined.

## Acknowledgments

This work was supported by the Dutch National Research Programme on Global Air Pollution and Climate Change (NRP) within Project 951235. All computations were performed on the CRAY C916 at the Academic Computing Centre (SARA), Amsterdam, Netherlands, within Project SC283. Use of these computing facilities was sponsored by the National Computing Facilities Foundation (NCF) with financial support from the Netherlands Organization for Scientific Research (NWO). Collaboration was initiated during visits of HD to UCLA (Los Angeles) and UH (Honolulu) in 1997 and 1998 that were sponsored by NSF Grant ATM-9521389 (J. D. Neelin) and an NWO PIONIER Grant (HD). FFJ is supported by NSF Grant ATM-9312888. Discussions with David Neelin (UCLA) are much appreciated.

## REFERENCES

### APPENDIX A

#### The Nondimensional Coupled Model and Its Discretization

The final nondimensional equations for the shallow-water component become

with the appropriate boundary conditions:

For the Ekman layer, the equations are

The nondimensional atmosphere model can be written as

and the dimensionless SST equation is given by

Variables are expanded into spectral basis functions, with Chebychev polynomials in zonal direction and Hermite functions in meridional direction. Using collocation techniques, a set of nonlinear algebraic equations is obtained for the steady states of the model. The analysis of the stability of these steady states leads to a generalized eigenvalue problem. Both steady states and their linear stability are traced through parameter space using continuation techniques as described in Dijkstra and Neelin (1995b).

To solve the equations (A1)–(A4) on the computational domain (*x, y*) ∈ [0, 1] × (−∞, ∞), a pseudospectral approach is followed. The dependent quantity *T* − *T*_{0} is thereto expanded as

where the basis functions *f*_{ij} are given by

where *y* is the oceanic meridional coordinate, the *ψ*_{j} are the Hermite functions, and the *C*_{i} are the Chebychev polynomials defined by

The atmospheric response to an arbitrary heat flux forcing *Q* is determined by writing *u*_{a} = *R* + *S* and *ϕ* = *R* − *S* in (A3). The variables *R, S,* and *V* and the forcing term *Q* = *T* − *T*_{0} are expanded into Hermite functions with the atmospheric meridional coordinate *y*_{a}. On using the recurrence relations for the Hermite functions, the system is reduced to the following set of equations for the unknowns *R*_{n}(*x*):

The forcing *Q* of the atmospheric component is obtained from the spectral representation of the SST in the ocean, and the *Q*_{n}(*x*) in (A7a) are written as

where the coefficients *c*_{nj} result from the projection of oceanic onto atmospheric Hermite functions. Solving (A7a) in terms of the coefficients *q*_{ij}, the zonal wind stress component becomes

with 𝒜_{ij} given by

The zonal wind stress forcing *τ*^{x}, (3), has atmospheric meridional scales. Hence, to obtain the correct ocean forcing, *τ*^{x} has to be projected onto the oceanic Hermite functions, that is,

with

Substitution of the expression for *u*_{a} yields

where *a*^{n}_{ij}(*x*) is given by

with 𝒜_{ij} as in (A9). The surface layer velocities in the ocean are solved directly in terms of the applied wind stress forcing, that is,

To solve the equations for the shallow-water component (A1), new variables *r* and *s* are introduced, that is, *u* = *r* + *s* and *h* = *r* − *s.* The unknowns *r* and *s* are expanded into *ψ*_{j}(*y*), and using the recurrence relations for the Hermite functions, the system (A1) is reduced to a set of equations for the unknowns *r*_{j}(*x*), which becomes

The coefficients ℰ_{j} result from the discretization of the boundary conditions and are given by

The unknowns *r*_{j} of the oceanic shallow-water component are again expanded into Chebychev polynomials, that is,

The SST equation (A4) is rewritten in terms of *T̂* = *T* − *T*_{0} and is, in addition to (A11), solved by collocation. The solution vector **z**, as presented in (6), contains the oceanic wave components *r*_{ij} and spectral coefficients of *T̂, q*_{ij} as defined in (A5). The collocation points in the meridional direction are the roots of the (*N*_{y} + 1)th Hermite function, when the sum in (A5) is truncated at *j* = *N*_{y}. The collocation points in the *x*-direction are the exterior grid points

The discretized set of equations thus contains 2*N*_{x}*N*_{y} unknowns, where *N*_{x} is the number of Chebychev polynomials and *N*_{y} the number of Hermite polynomials. The truncation of the number of basis functions in the spectral representation was determined from resolution studies of the first bifurcation point. While the results in section 3 were computed using *N*_{x} = 17 and *N*_{y} = 41, the sensitivity results presented in sections 4 and 5 were done for *N*_{x} = 11 and *N*_{y} = 21. Despite the infamous convergence properties of the Hermite functions, the patterns and periods found for *N*_{x} = 17 and *N*_{y} = 41 differed only slightly from those at lower resolution.

### APPENDIX B

#### The Complete Spectrum of the Fully Coupled Model

The complete spectrum of the uncoupled ocean–atmosphere model consists of eigenvalues originating from the SST equation and from uncoupled ocean dynamics. The latter were discussed by Moore (1968), Cane and Sarachik (1981), Cane and Moore (1981), and Neelin and Jin (1993) and are shown in Fig. B1a. The ocean dynamics spectrum consists of ocean basin modes, indicated by triangles in Fig. B1, and scatter modes. A characteristic feature of the scattermodes is their appearance as rays in the spectrum, centered around the stationary scatter mode with the smallest damping rate. The number of rays is equal to the degrees of freedom in the meridional direction, while their length is determined by the zonal resolution of the model discretization (cf. Neelin and Jin 1993). The uncoupled SST spectrum is located near Re(*σ*) = −2.8 yr^{−1} in Fig. B1a. In Fig. B1b the fully developed spectrum at large coupling (*μ* = 0.74) is plotted. The initially oscillatory instability has split into two stationary modes, which can be discerned to the right. The SST related part of the spectrum shows the development of oscillatory modes, which either grow or decay. This is more clearly seen in Fig. B1d, where, in the fast-wave limit, the ocean dynamics modes have been filtered from the spectrum. The leading modes are high-frequency asymmetric modes, related to the meridional advection terms in the SST equation. In the fast SST limit, Fig. B1c, the SST modes have been filtered, showing that, apart from the selection of the ENSO mode, the spectrum is hardly affected by coupling. The ocean basin modes, relevant in several previous studies (Cane et al. 1990; Dijkstra and Neelin 1999), are indicated by triangles.

## Footnotes

*Corresponding author address:* Dr. Paul C. F. van der Vaart, Department of Physics and Astronomy, Institute for Marine and Atmospheric Research Utrecht, Utrecht University, Princetonplein 5, 3584 CC Utrecht, Netherlands.

Email: pvaart@phys.uu.nl