## Abstract

A new conceptual model for ENSO has been constructed based upon the positive feedback of tropical ocean–atmosphere interaction proposed by Bjerknes as the growth mechanism and the recharge–discharge of the equatorial heat content as the phase-transition mechanism suggested by Cane and Zebiak and by Wyrtki. This model combines SST dynamics and ocean adjustment dynamics into a coupled basinwide recharge oscillator that relies on the nonequilibrium between the zonal mean equatorial thermocline depth and wind stress. Over a wide range of the relative coupling coefficient, this recharge oscillator can be either self-excited or stochastically sustained. Its period is robust in the range of 3–5 years. This recharge oscillator model clearly depicts the slow physics of ENSO and also embodies the delayed oscillator (Schopf and Suarez; Battisti and Hirst) without requiring an explicit wave delay. It can also be viewed as a mixed SST–ocean dynamics oscillator due to the fact that it arises from the merging of two uncoupled modes, a decaying SST mode and a basinwide ocean adjustment mode, through the tropical ocean–atmosphere coupling. The basic characteristics of this recharge oscillator, including the relationship between the equatorial western Pacific thermocline depth and the eastern Pacific SST anomalies, are in agreement with those of ENSO variability in the observations and simulations with the Zebiak–Cane model.

## 1. Introduction

Interaction between the tropical ocean and atmosphere produces interannual climate variability dominated by the El Niño–Southern Oscillation (ENSO) phenomenon, which has been the object of much research (e.g., Bjerknes 1969; Wyrtki 1975; Rasmusson and Carpenter 1982; Cane and Zebiak 1985; Cane et al. 1986; Graham and White 1988; Philander 1990; Barnett et al. 1991; Philander et al. 1992). Theory for the ENSO phenomenon has approached a mature stage during the past decade (Neelin et al. 1994). ENSO modeling has advanced to the point where predictions are being made on a regular basis (Cane et al. 1986; Zebiak and Cane 1987, hereafter ZC; Barnett et al. 1988; Keppenne and Ghil 1992; Latif et al. 1994; Ji et al. 1994).

Bjerknes (1969) first hypothesized that ENSO is a result of ocean–atmosphere interaction in the tropical Pacific from his analysis of the empirical relations of El Niño and the Southern Oscillation. He recognized that the equatorial SST zonal gradient drives the easterlies over the tropical Pacific. These easterlies in turn create the cold SST over the eastern Pacific and thus strengthen the SST gradient. It is this Bjerknes positive feedback process of tropical ocean–atmosphere interaction that sustains either a warm or a cold SST anomaly. Wyrtki (1975) realized that during El Niño the oceanic anomalies in sea-level data are basinwide and dynamical. He suggested that the buildup of sea level (an indicator of heat content) over the western Pacific is related to the strengthening of the trades, and the accumulated warm water flows eastward in the form of Kelvin waves to give birth to an El Niño event. However, this Bjerknes–Wyrtki hypothesis did not explain the salient cyclic nature of ENSO. It took about another decade of research before Cane and Zebiak (1985) and Wyrtki (1986) envisioned a new hypothesis (Cane 1992a):

An earlier version (Cane and Zebiak 1985; Cane et al. 1986) emphasized the recharging of the equatorial “reservoir of warm water” as a necessary precondition for the initiation of a warm event. On the basis of his analysis of sea level data, Wyrtki (1986) developed a similar hypothesis. The aftermath of a warm event leaves the thermocline along the equator shallower than normal (i.e., equatorial heat content is low and SST is cold; this is the La Niña phase). Over the next few years the equatorial Kelvin waves allowed by linear equatorial ocean dynamics can move enough of the warm water to the eastern end of the equatorial Pacific to initiate the next event.

The combined BWCZ (short for Bjerknes–Wyrtki–Cane–Zebiak) hypothesis indicated that ENSO is a natural basinwide oscillation of the tropical Pacific ocean–atmosphere system. Both the positive feedback of the coupled interaction of the ocean–atmosphere system of the tropical Pacific and the memory of the system in the subsurface ocean dynamical adjustment are essential for ENSO. The ZC model that incorporates these two elements is the first and perhaps still among the most successful coupled models in simulating, understanding, and predicting ENSO.

The theoretical understanding gained over the past decade is largely through mechanistic studies based on simple coupled models, which often include shallow-water or two-layer ocean models coupled to steady shallow-water-like atmospheric models with heavily parameterized physics (Lau 1981; Philander et al. 1984; Anderson and McCreary 1985; Cane and Zebiak 1985; Zebiak and Cane 1987; Hirst 1986, 1988; Battisti 1988; Battisti and Hirst 1989, BH hereafter; Suarez and Schopf 1988; Schopf and Suarez 1988, SS hereafter; Yamagata and Masumoto 1989; Wakata and Sarachik 1991; Neelin 1991; Jin and Neelin 1993a,b; Neelin and Jin 1993, hereafter JN). The basic assumption underlying these studies is that ENSO-like variability can be considered as resulting from a coupled oscillatory instability of the basic climatological state of the tropical ocean–atmosphere system. Different conceptual models of coupled unstable modes have been put forward to explain the diverse results obtained in various coupled models. For example, under the assumption that the SST anomaly is in quasi-equilibrium with subsurface thermocline fluctuation, Cane et al. (1990, CMZ hereafter) discovered that the coupled process modifies the ocean basin modes (Cane and Moore 1981) to form a standing unstable oscillatory mode, a coupled wave oscillator that is strongly related to oceanic wave dynamics. In contrast, a coupled “SST mode” (Neelin 1991; JN; Hao et al. 1993) that depends on the feedback of surface-layer processes with quasi-equilibrium ocean dynamics is found to be a good candidate to explain the propagating interannual oscillations simulated in certain coupled models (Neelin 1990; Meehl et al. 1990).

The detailed connections between these two seemingly contrary idealizations were quantitatively mapped out by JN using a stripped-down version of the ZC model. It was found that the regimes of the unstable ocean dynamics modes of CMZ and the coupled SST modes of Neelin (1991) represent two extremes: the fast SST limit and the fast wave limit. In the former case, the timescale of SST change is much shorter than that of the oceanic dynamics, and SST change is thus essentially controlled by subsurface feedback processes. In the latter case, the timescale of oceanic adjustment is much shorter than that of SST, so that the subsurface ocean dynamics is in quasiequilibrium with wind stress forcing and SST change is dominated by surface layer processes. It is made clear by JN that these two types of coupled modes merge continuously across parameter space. In particular, in the physically most relevant part of parameter space, they are strongly mixed to form a “mixed” SST–subsurface ocean dynamics mode. From this unifying framework, the diversity of interannual behavior found in different models can be viewed in terms of slight “distortions” in the mixed SST–subsurface ocean dynamics modes, arising from different approximations of model physics.

The concise and toylike delayed oscillator (SS; BH) model incorporates two elements: 1) the coupled positive feedback for SST anomalies over the central to eastern Pacific through changes of the thermocline and upwelling and 2) memory of subsurface ocean dynamics. Because of the latter, the thermocline anomalies consist of the contributions both from local wind stress forcing and from a delayed signal from the western boundary produced by the reflection of remotely generated Rossby waves. The delayed oscillator requires a time delay due to the equatorial wave propagation and a strong wave reflection at the western boundary. In fact, it was shown by CMZ that the coupled wave oscillator at the fast SST limit can be approximated into a BH delayed oscillator when the wave reflection is ignored at the east boundary. This has led to an ENSO paradigm with a focus on the equatorial wave propagation instead of the recharge–discharge of equatorial heat content suggested in the earlier BWCZ hypotheses. The simple and popular delayed oscillator, although challenged (Li and Clarke 1994) about the role of propagation of the equatorial ocean waves, was argued to be consistent with observation (Mantua and Battisti 1994) and the ENSO-like variability in the ZC model. On the other hand, it is speculated by JN that the delayed oscillator is an ad hoc approximation to the mixed SST–ocean dynamics mode, which is the most relevant mode to ENSO. However, a clear and simplified representation of this mixed mode and demonstration of its close linkage with the BWCZ hypothesis are yet to be established.

In Part I of this two-part paper, a simple prototype model explicitly consistent with the BWCZ hypothesis is shown to give a harmonic oscillator that is at the heart of the mixed SST–ocean dynamics mode suggested by JN. This oscillator can be systematically reduced from a ZC-type model (Part II). In section 2, the conceptual harmonic oscillator based on the BWCZ hypothesis for ENSO is formulated and discussed. In section 3, the excitation of the recharge oscillator is examined in linearly unstable and stable regimes with considerations of the nonlinearity and stochastic forcing. Further discussion of the recharge oscillator and a summary are given in sections 4 and 5, respectively.

## 2. A conceptual model for the recharge oscillator

The thermocline depth of the equatorial ocean is one of the main variables that needs to be considered for understanding ENSO dynamics. On the ENSO timescale, the leading ocean dynamical balance is between the pressure gradient force and wind stress over the equatorial band, for instance, within a couple of oceanic Rossby radii of deformation. In other words, the zonal pressure gradient force accompanying the thermocline depth tilt along the equator is largely in a Sverdrup balance with the equatorial wind stress force (e.g., Cane and Sarachik 1981; Philander 1990). A numerical ocean model forced with a slowly changing wind stress anomaly can be used to verify this quasi-equilibrium relation (e.g., Schneider et al. 1995). For simplicity, one can take the reduced gravity shallow-water model as a base for the upper ocean dynamics. This quasibalance thus leads to a simple relationship,

where *h*_{W} denotes the thermocline depth anomaly in the western Pacific, for instance, within one oceanic Rossby radius of deformation from the equator; *h*_{E} is the thermocline depth anomaly in the equatorial eastern Pacific; and *τ̂* is proportional to the zonally integrated wind stress in this band. This relation approximately holds for the anomalies associated with ENSO, since ENSO related *τ̂* changes slowly in time relative to the timescale of the equatorial oceanic Kelvin waves and is roughly centered at the equator with latitudinal variations on the scale of the atmospheric Rossby radius. However, this leading balance of forces only constrains the east–west contrast of the thermocline depth. The absolute depth at the western Pacific or the mean thermocline depth over the equatorial band is not constrained by this balance. The mean thermocline depth depends on the mass adjustment of the entire tropical Pacific ocean and may not be in equilibrium with the slowly varying wind forcing. In the nonequilibrium between the zonal mean thermocline depth and the wind stress forcing resides the memory of the subsurface ocean dynamics and makes the fast wave limit approximation inappropriate for an ENSO-like mode that relies on the subsurface ocean memory (e.g., Cane 1992b).

There are two equivalent qualitative approaches to this adjustment process. One is from the wave propagation point of view. The ocean mass adjustment is completed through the oceanic Kelvin and Rossby waves, which propagate in opposite directions and are forced by the basin boundaries to change into opposite wave characteristics through reflections. The damping effect due to the leakage of energy via the boundaries and other damping processes allow this adjustment to finally settle into a quasi-equilibrium in the equatorial region. In this view, one of the two unknowns in Eq. (2.1) can be related to a boundary condition that determines wave reflections as shown in a number of theoretical studies (e.g., Cane and Sarachik 1981). The other approach is to consider that the equatorial wave propagation process is relatively fast for establishing this thermocline slope that extends to the off-equatorial region as a result of the broadness of the atmospheric wind system. The Coriolis force becomes important off the equatorial band, and therefore there will be Sverdrup transport either pumping the mass in or out of the equatorial region depending on wind forcing, as hypothesized by Cane and Zebiak (1985) and Wyrtki (1986). Under linear shallow-water dynamics, this Sverdrup transport is accomplished by the Rossby waves. The zonally integrated effect of this Sverdrup transport of mass or, equivalently, heat content (in the sense of shallow-water dynamics) will result in the deepening or shoaling of the western Pacific thermocline depth (Wyrtki 1986). Thus, although the thermocline tilt along the equator is set up quickly to balance the equatorial wind stress as expressed by (2.1), the thermocline depth of the warm pool takes time to adjust to the zonal integrated meridional transport, which is related to both the wind stress and its curl off the equatorial band.

Using the second approach and assuming that the adjustment timescale is much longer than that for a Kelvin wave crossing the basin, one can symbolically describe this adjustment process by the following equation:

This equation focuses on the thermocline depth changes averaged over the western equatorial Pacific during the basinwide adjustment because the tropical wind anomaly associated with ENSO is largely over the western to central Pacific (Deser and Wallace 1990). The first term on the right-hand side represents the ocean adjustment. It is simply characterized by a damping process with a rate *r* that collectively represents the damping of the upper ocean system through mixing and the equatorial energy loss to the boundary layer currents at the east and west sides of the ocean basin. This will become more clear in Part II of this paper where this relation is derived from basic dynamical principles through a number of simplifications. The wind forcing *F*_{τ} is related to the zonally integrated wind stress and its curl. It represents the Sverdrup transport across the basin. To a large extent *F*_{τ} is also proportional to *τ̂*, that is, *F̂*_{τ} = *α**τ̂*, which approximately holds for the ENSO-related wind stress anomaly with a broad meridional scale. One also expects that *F*_{τ} is a weak forcing because only a part of the wind stress forcing is involved in the slow adjustment process, whereas the other part is in the quasi-Sverdrup balance. A small *r* and a weak forcing *F*_{τ} are consistent in describing the slow basinwide ocean adjustment process. The minus sign for the wind forcing term comes from the fact that a westerly wind stress anomaly will lead to a shallower thermocline over the western Pacific, whereas a strengthened trade will result in a buildup of the warm pool as suggested by Wyrtki (1975, 1986). The Eq. (2.2) can be rewritten as

Equations (2.1) and (2.3) give a gross description of the basinwide equatorial oceanic adjustment under anomalous wind stress forcing of low frequencies (with timescales larger than the basin crossing time of oceanic Kelvin waves). Equation (2.3) differs from the delayed oscillator model (BH) in description of slow ocean dynamics. In the latter, *h*_{W} is simply proportional to −*τ̂*(*t* − *η*), which results from the reflection of forced equatorial oceanic Rossby waves with a delay *η* (the time for Rossby wave propagating from the center of the forced region to the western boundary). This wave-delay description introduces an infinite number of degrees of freedom (representing different spatial scales of Rossby waves), yet suffers a degeneracy at *τ̂* ≡ 0 and fails to properly describe the ocean adjustment dynamics. In (2.3), the explicit wave-propagating process is omitted, and the collective role of equatorial waves in achieving the quasi-equilibrium adjustment dynamics is accounted in the recovery timescale 1/*r.* The basinwide slow adjustment dynamics is simply characterized by one degree of freedom and the degeneracy in the wave-delay description is avoided. Moreover, Eqs. (2.1) and (2.3) generally represent the quasi-balance between the thermocline tilt and wind stress and the slow dynamic renewal of the warm pool heat content. The whole basinwide distribution of heat content, residing in the memory of subsurface ocean dynamics, is thus explicitly taken into consideration. The slow buildup of the western Pacific warm pool in terms of thermocline depth, or the recharging of the entire equatorial heat content, takes place during the phase of strengthened trade wind (*τ̂* < 0), whereas a weakened trade results in the gradual discharging of the western Pacific warm pool and the reducing of the entire equatorial thermocline depth. Thus, (2.3) gives a simple and clear description of the equatorial heat content recharge–discharge process, which is the essential phase-transition mechanism of ENSO as suggested by Cane and Zebiak (1985) and Wyrtki (1986).

The variation of SST during ENSO is largely confined within the central to eastern equatorial Pacific. The SST anomaly in this region strongly depends on the local thermocline depth that determines the temperature of the subsurface water, because this water is pumped up into the surface layer to control the SST by the climatological upwelling associated with the climatological trade wind along the equator. Changes in the trade wind intensity in response to the SST anomaly may also further reinforce the SST anomaly by altering upwelling and horizontal advection. The mean climatological upwelling and heat exchange between the atmosphere and ocean tend to damp out the SST anomaly. Although the details of all these processes can be complicated (e.g., ZC, BH, JN), they can be roughly depicted in a simple equation for the SST anomaly *T*_{E}, averaged over the central to eastern equatorial Pacific:

The first term on the right-hand side is the relaxation of SST anomaly toward climatology (or zero anomaly) caused by the above-mentioned damping processes with a collective damping rate *c.* The second and third terms are the thermocline upwelling and advective feedback processes, respectively; *τ*_{E} is wind stress averaged over the domain where the SST anomaly occurs; *γ* and *δ _{s}* are thermocline and Ekman pumping feedback coefficients.

Atmospheric response to a warm SST anomaly of the central to eastern Pacific is a westerly wind over the central to western equatorial Pacific and an easterly anomaly to the east of the SST anomaly. There is an overall westerly (easterly) anomaly for a positive (negative) SST anomaly averaged over the entire basin of the equatorial band, but a much weaker westerly (easterly) anomaly averaged over the eastern half of the basin (e.g., Deser and Wallace 1990). This allows the simple approximate relations of the wind stress and SST anomalies:

where *b* and *b*′ are coupling coefficients.

Combining the Eqs. (2.1)–(2.5), one obtains a simple linear coupled system with both the subsurface ocean-adjustment dynamics and the surface-layer SST dynamics:

where *R* = *γ**b* + *δ*_{s}*b*′ − *c.* Clearly, *R* collectively describes the Bjerknes positive feedback hypothesis of the tropical ocean–atmospheric interaction and leads to instability when growth rate (*R* − *r*)/2 >0, while *α**b**γ*, representing the recharging mechanism, allows oscillation when frequency *ω* =*α**b**γ* − (*r* + *R*)^{2}/4 is real. Thus, this system combines the mechanisms of the BWCZ hypotheses to describe a coupled recharge oscillator. The growth rate of the oscillation depends mainly on the term *γ**b,* which represents the thermocline feedback, because the Ekman upwelling feedback parameter *δ*_{s}*b*′ is relatively small due to the weak local wind stress averaged over the eastern basin. Thus this term will be ignored for simplicity throughout the paper by setting *δ*_{s}*b*′ = 0 or *R* = *γ**b* − *c.*

The physics of this simple recharge oscillator of (2.6) is clear as schematically illustrated in Fig. 1. An initial positive SST anomaly induces a westerly wind forcing over the central to western Pacific. The anomalous slope of the equatorial thermocline is promptly set up to be proportional to the wind stress and thus to the SST anomaly. The deepening of the thermocline in the eastern Pacific results in a positive feedback process that amplifies this anomaly and brings the oscillation to a mature phase as shown in Fig. 1a (the mechanisms of sustaining a finite amplitude of the recharge oscillation will be discussed in section 3). At the same time, the wind stress also gradually reduces thermocline depth in the western Pacific and leads to a negative zonal mean thermocline depth across the Pacific as the result of the divergence of zonal integrated Sverdrup transport. This process can be viewed as the discharge of the zonal mean equatorial heat content because thermocline depth is highly related to the dynamical part of the upper ocean heat content (e.g., Zebiak 1989). This discharge also gradually reduces the thermocline depth in the eastern Pacific and eventually leads to a cooling trend for the SST anomaly when the warming due to the still positive thermocline depth anomaly in the eastern Pacific is balanced off by the damping process of the SST. Thus the warm phase evolves to the transition phase as shown in Fig. 1b. At this time, the SST anomaly cools to zero. The east–west thermocline tilt diminishes because it is always in a quasi-equilibrium balance with the zonal wind stress, which has disappeared following the SST anomaly. However, the entire equatorial Pacific thermocline depth and thus the eastern Pacific thermocline depth is anomalously shallow because of the discharge of the equatorial heat content during the warm phase. This zonal mean thermocline is not in quasi-equilibrium with the wind stress and has already started to increase as the result of the ocean adjustment after the diminished wind stress. It is this anomalous shallow thermocline depth at the transition that allows anomalous cold water to be pumped into the surface layer by climatological upwelling; the SST anomaly then slides into a negative phase. Once the SST anomaly becomes negative, the cooling trend proceeds because the negative SST anomaly will be further amplified through the positive thermocline feedback. That is, the enhanced trades, in response to the cold SST anomaly, deepen the thermocline depth in the western Pacific and lift the thermocline depth up in the east. Thus the oscillation develops into its mature cold phase as shown in Fig. 1c. At the same time, the zonal mean thermocline depth over the equatorial Pacific is deepening, as a result of the recharging of the equatorial heat content due to the strengthened trades. This reverses the cooling trend after the mature cold phase and brings it to another transition phase as shown in Fig. 1d. When the cold SST anomaly reduces to zero, the positive zonal mean thermocline depth generated by the recharging process will lead the SST anomaly evolving back to another warm phase.

During the cycle, a positive western Pacific thermocline depth anomaly leads to a warm SST anomaly in the eastern Pacific. This can also be viewed as the thermocline depth anomaly of the western Pacific being negatively proportional to the SST anomaly of the eastern Pacific with a time lag. Thus, the nonequilibrium between zonal mean thermocline depth and the wind stress forcing, owing to the slow basinwide ocean adjustment, serves as the memory of the coupled system. The discharge of the equatorial heat content in the warm phase and recharge at the cold phase serve as the phase-transition mechanism from a warm event to a cold event and vice versa. This conceptual basinwide coupled recharge oscillator model is, perhaps, the simplest representation of the entire BWCZ hypothesis.

To facilitate some quantitative analysis, one needs to estimate the values of the parameters in the coupled system (2.6). The collective damping rate *c* is dominated by the mean climatological upwelling, which yields a typical damping timescale of about 2 months (e.g., JN); *γ* is related to both the mean climatological upwelling and the sensitivity of subsurface ocean temperature to the thermocline depth. It is chosen to give an SST change rate of 1.5°C over 2 months (which is the upwelling timescale) per 10 m of thermocline depth anomaly over the eastern Pacific. This value of *γ* is close to that parameterized in the ZC model. The collective damping rate *r* in the ocean adjustment includes not only the weak linear damping [about (2.5 years)^{−1} in ZC model], but more importantly, also the damping effect due to the loss of energy to the boundary currents of the west and east boundary layers (see Part II for details). The latter could lead to a much shorter damping timescale. Parameter *r* is set as (8 months)^{−1}, which still gives a damping timescale much longer than the typical 2-month-crossing timescale of the first baroclinic Kelvin wave. If one assumes that for a given steady wind stress forcing, the zonal mean thermocline anomaly of this linear system is about zero at the equilibrium state, that is, *h*_{E} + *h*_{W} = 0, then from (2.1) and (2.3), one finds that *α* shall be about half of *r*. Finally, parameter *b* is a measure of the thermocline slope, which is in balance with the zonal wind stress produced by the SST anomaly. It is chosen to give 50 m of east–west thermocline depth difference per 1°C of the SST anomaly. This is a high-end estimation of parameter *b*, whereas for a typical ENSO event the value may be somewhat smaller. For the convenience of discussion, a so-called relative coupling coefficient *μ* is introduced as

where *b*_{0} is the high-end estimation of *b,* and *μ* is referred to as a relative coupling coefficient that will be varied in the range 0 to 1.5 embracing the uncoupled and strongly coupled cases. System (2.6) is hereafter nondimensionalized by scales of [*h*] = 150 m, [*T*] =7.5 K, and [*t*] = 2 months for anomalous thermocline depth, SST, and the time variable, respectively. Accordingly, parameters *c, r,* and *α* are scaled by [*t*]^{−1}, and parameters *γ* and *b*_{0} by [*h*][*t*]/[*T*] and [*T*]/[*h*]. Their nondimensional values are *c* = 1, *γ* = 0.75, *r* = 0.25, *α* = 0.125, and *b*_{0} = 2.5.

The detailed analyses of the model sensitivity to this set of parameters will be elaborated in Part II of the paper. With the values of the parameters given above, the dependence of the eigen modes of the coupled system (2.6) on the relative coupling coefficient can be analytically solved as

and the results are shown in Fig. 2. When the relative coupling coefficient *μ* is weak (*μ* < *μ*_{1}), the system has two decaying modes that can be identified at *μ* = 0 as the uncoupled SST model and the ocean-adjustment mode respectively. These two modes eventually merge into an oscillatory mode as the coupling coefficient increases to *μ* > *μ*_{1}. When the coupling is further increased to *μ* > *μ*_{2}, the oscillator breaks down to give two modes. One is a purely growing mode because the strong coupling through Bjerknes feedback results in a growth rate being too fast to be linearly checked by the slow ocean adjustment process. The other is a real mode whose growth rate decreases rapidly and becomes negative when *μ* is larger than *μ*_{2}. As is shown in Fig. 2, for a wide range of moderate coupling between (*μ*_{1}, *μ*_{2}), the system does support an oscillatory mode with a period mostly in the range of 3–5 years. This oscillatory mode is supercritical if *μ* > *μ*_{c} = 2/3 and subcritical if the coupling coefficient is smaller than the critical value *μ*_{c}. As will be shown in section 3, in the supercritical oscillatory regime, adding nonlinearity to the coupled system will limit the linear growth to a self-excited coupled oscillation. In the subcritical oscillatory regime, a coupled oscillatory solution can also be sustained when sources of stochastic forcing are taken into consideration.

Before going into the nonlinear solutions or stochastically forced solutions, one can look at the linear neutral oscillation at the critical coupling. When *μ* = 2/3, the growth rate in (2.8) is zero and the critical frequency is *ω*_{c} = 3/32, and (2.6) gives a perfect harmonic oscillation:

The amplitude of the oscillation depends only on the initial conditions of the eastern Pacific SST *T* ^{0}_{E} and the western Pacific thermocline depth *h*^{0}_{W}. The western Pacific thermocline depth anomaly in the solution is also written to be negatively proportional to the SST with a time lag *η*.

An example of the solution (2.9) is shown in Fig. 3 with a period around 3.4 years and a roughly 6 month time lag between the positive maximum of the SST in the eastern Pacific and the negative maximum of the thermocline depth at the western Pacific. The trajectory of the solution on the phase plane of (*h*_{W}, *T*_{E}) exhibits an elliptical-shaped orbit centered at the origin and rotating clockwise. The ellipticity and rotating direction of the solution along the orbit heavily depend on the lag. When this lag is a quarter of the period, the shape of the limit cycle is a perfect circle. A smaller lag makes the orbit become more elliptical. If the lag is negative and smaller than a quarter of the period, the rotation direction reverses. From the BWCZ hypothesis, only a positive lag *η* smaller than a quarter of the oscillation period is possible for the harmonic oscillation of (2.6), simply because both *r* and *α**b* are positive to characterize the damping effect in the ocean adjustment and the buildup of the western Pacific thermocline for a strengthening trade. In other words, the BWCZ hypothesis predicts this feature of an oscillation for ENSO well as shown in Fig. 3. The solution of recharge oscillator qualitatively agrees with the features seen in Fig. 4 in terms of the orientation and rotation direction of the trajectory on the phase plane of the observed equatorial sea level anomaly and the observed anomalous NINO3 index (SST averaged over an equatorial band of the eastern Pacific) during the period 1979–1991. Similar features can also be found in a typical solution from the standard ZC model in an active ENSO episode as shown in Fig. 5. Unlike the results in Fig. 3a, the observed trajectory and that of the ZC model solution do not precisely follow closed cycles centered at the origin of the phase plane. These discrepancies come from the nonlinearity and irregularity of ENSO variability in the observation and the ZC model, which are not captured by the harmonic oscillation of the conceptual model. Nevertheless, the fundamental features of ENSO variability, the ellipticity of the phase space trajectory, the rotation direction, and the periodicity are well depicted by the simple linear dynamics. The qualitative agreement in all plots indicates the relevance of the conceptual model.

## 3. Excitation of the recharge oscillator

### a. Self-excitation

When the coupling intensity is sufficiently strong (*μ* > 2/3), the linear growth rate in (2.8) becomes positive, and system (2.6) gives an unstable oscillatory mode and thus a Hopf bifurcation from the climate state about which the model has been linearized. Nonlinear processes will limit the linear growth and lead to a finite amplitude oscillation. As shown by Battisti and Hirst (1989), the dominant nonlinearity in the ZC model is in the thermocline feedback, that is, a very shallow or very deep thermocline will not linearly increase the cooling or warming because of the nonlinear dependence of the SST on the thermocline depth in subsurface temperature parameterization. Such a nonlinear dependence reflects the nonlinear vertical distribution of the temperature in the tropical upper ocean. One can consider such a nonlinear dependence of subsurface temperature on the thermocline depth by adding a term −*e*_{n}*h*^{3}_{E} into the SST equation of (2.6) (e.g., BH). Using (2.1), one gets a simple nonlinear model that is valid in the neighborhood of the climatological state about which the model has been linearized:

Near the Hopf bifurcation point, one can assume *μ* = *μ*_{c} + Δ with 0 < Δ ≪ 1. This Δ represents a small departure from the critical point *μ*_{c}. Then the finite amplitude oscillation of system (3.1) can be approximately solved using a perturbation method such as the one outlined in the appendix or some other methods (e.g., Iooss and Joseph 1990). The final solution can be expressed as

where

and *η* is the same as defined in (2.9). Clearly, when the system is supercritical (Δ > 0), the solution is also a perfect harmonic oscillation and very similar to the linear neutral solution (2.9). However, the amplitude does not depend on initial conditions but on the supercriticality because the solution is now a limit cycle. In other words, the recharge oscillator of (3.1) is now self-excited. The period is weakly controlled by the supercriticality and largely determined by *ω*_{c}, the frequency at the bifurcation point. The frequency correction factor Ω and the amplitude factor B are functions of *e*_{n} and other model parameters, and the details can be found in the appendix. Because the correction on the frequency is reduced by the nonlinearity as discussed in the appendix, one can virtually ignore the small frequency correction. The time lag *η*_{ω} between the SST and thermocline depth thus approximately equals *η*. This limit cycle solution then becomes the same as the linear solution (2.9) in terms of characteristics of the orbit in the phase space, time evolution of the SST, and thermocline depth as shown in Fig. 3 after a simple rescaling of the amplitude. Thus, in the neighborhood of the Hopf bifurcation, nonlinearity does not alter basic dynamics of the recharge oscillator as described by (2.6).

If one pushes the coupling parameter to an extremely strong range, for instance *μ* > *μ*_{p} = 16/15, system (3.1) will have two steady-state solutions in addition to the limit cycle solution of the nonlinear recharge oscillator, as shown in the appendix. These two steady-state solutions, however, are unstable until *μ* > 19/15. When *μ* > 19/15, they become stable states and the nonlinear recharge oscillator disappears. Thus, in this extreme regime, the coupled model settles to either a warm state or a cold state. These new states are far away from the climatology state about which the anomaly coupled model is constructed. They often represent artifacts of anomaly coupled models and thus the danger of the so-called flux corrections implied in these models (Neelin and Dijkstra 1995). However, the simple model (3.1) favors the Hopf bifurcation as the first bifurcation. The simple nonlinearity actually further extends the oscillatory regime far beyond the range for linear oscillations to avoid the danger of setting the model to a permanent cold or warm state in the relevant range of the relative coupling coefficient. The detailed form of nonlinearity crucially determines the bifurcation structure beyond the neighborhood of the first Hopf bifurcation. Thus the bifurcation diagram shown in the appendix for (3.1) may be significantly altered with the inclusion of advective nonlinearity in the SST and ocean dynamics and other nonlinear processes, for instance, in the atmospheric responses to SST anomalies. Nevertheless, the stable steady-state solutions tend to appear far away from the realistic parameter regime. This is consistent with the fact that anomaly coupled models in ENSO modeling are successful in simulating ENSO-like interannual variability in the realistic parameter range without facing the great danger of being trapped into a spurious permanent warm or cold state.

### b. Stochastic excitation

When the coupling intensity is not sufficiently strong (i.e., *μ* < 2/3), the linear solution of (2.6) decays because the growth rate in (2.8) is negative. In this case extending (2.6) to the nonlinear system (3.1) does not help to sustain the oscillation of the system. However, the recharge oscillator still can be excited by stochastic forcing. The stochastic excitation of a decaying coupled oscillator has been suggested as a scenario for ENSO by a number of authors (McWilliams and Gent 1978; Lau 1985). Viewing the ENSO as episodic, Wyrtki (1986) also proposed this scenario in conjunction with the recharge mechanism for ENSO. More recently, Penland and Sardeshmukh (1995) further argued in favor of this scenario based on the inverse modeling.

There are abundant transient atmospheric disturbances of much shorter decorrelation timescales than the ENSO timescale. These transient disturbances can exert forcing on the coupled system through wind stress and heat flux. They can be conceptually considered as sources of noise agitation to the slow ENSO dynamical systems (2.6) and (3.1). Thus the linear conceptual model (2.6) can be extended as follows:

where *ξ*_{1} and *ξ*_{2} represent the random wind stress forcing added to *τ̂* in equations (2.1) and (2.3) and random heating added into equation (2.4), respectively. For the validity of the quasi-Sverdrup balance and slow adjustment dynamics described in (2.1) and (2.3), *ξ*_{1} should not include high-frequency variability (with timescale shorter than the basin crossing time of oceanic Kelvin waves). For mathematical simplicity, *ξ*_{1} and *ξ*_{2} are assumed as two independent sources of white noise forcing with 〈*ξ*_{k}(*t* + *s*), *ξ*_{j}(*t*)〉 = *σ*^{2}_{j}*δ*(*s*)*δ*_{k,j} (〈 〉 denotes the ensemble mean) and *σ*^{2}_{1}, *σ*^{2}_{2} as their variances for mathematical simplicity. The strong dynamical filtering effect of coupled system (3.4) to high-frequency forcing makes the high-frequency part in *ξ*_{1} have little impact on the response spectra of the system (3.4).

With the system in the subcritical regime *μ* < *μ*_{c}, (3.4) describes the damped recharge oscillator with the stochastic excitation, and the solutions are bounded with finite variance. Using the Fourier transform pair

one can obtain the ensemble-averaged power spectra and the cross spectrum of the SST and thermocline depth of the system as follows:

where the asterisk denotes the complex conjugate.

The spectra of the SST and thermocline depth exhibit strong peaks near the eigen frequency (*ω* = 27/16) of (2.6) for the slightly damped case under the stochastic forcing in heating (Fig. 6a). In the more damped case, the peaks are much less evident and of much less amplitude (Fig. 6b). In the nonoscillatory regime these peaks disappear. The results under the stochastic wind forcing are very similar to those under the stochastic heating forcing (not shown). Integrating Eq. (3.4) following a standard numerical method (e.g., Penland and Matrosova 1994) with a small time step (6 h) to ensure a sufficient sample interval on the noise, one can obtain simulated realizations of SST and thermocline depth as shown in Fig. 7. Clearly, under the same magnitude of forcing, the recharge oscillator can be vigorously excited in the slightly damped case, whereas in the more heavily damped case, although still being excited, the oscillation is less coherent with a much weaker amplitude. These simulated realizations are typical because the spectra shown in Fig. 6 can be also approached by taking an ensemble average of the power spectra calculated from the simulated finite time realizations with a standard Fourier transform method if the number in the ensemble average is sufficiently large (over 100 simulations, not shown). These simulated evolutions of SST and thermocline depth are again similar to the result shown in Fig. 3a and the observation as shown in Fig. 4, particularly when the time series are smoothed and plotted in the same manner (not shown).

Inverse Fourier transforms of the spectra in (3.6) give the lagged autocorrelations of SST and thermocline depth, 〈*T*_{E}(*t* + *s*), *T*_{E}(*t*)〉 and 〈*h*_{W}(*t*), *h*_{W} (*t* + *s*)〉 (*s* denotes the lag), and their lagged cross correlation, 〈*T*_{E}(*t*), *h*_{W}(*t* + *s*)〉. They can be analytically obtained from residue calculus and the results under stochastic heating are shown in Fig. 8. The autocorrelations decay with the lag time on the timescale of the negative growth rate of the linear oscillator. The lag interval between the first and second zeros of autocorrelation, which can be most clearly identified in the slightly damped case in Fig. 8a, corresponds to a half of the period inferred from the eigen frequency of the linear damped oscillator. The cross correlation indicates the phase relation between the thermocline depth and SST. For the cases shown in Figs. 8a and 8b, the thermocline depth at the western Pacific is negatively correlated to the SST with a time delay of about 4–5 months, slightly shorter than 6 months of deterministic oscillations of (2.9) and (3.2). Under stochastic wind stress forcing, this lag is shorter, only about 3 months in the case of *μ* = 0.5 as shown in Fig. 9, and about 4 and 2 months in the cases of*μ* = 0.4 and 0.6 (not shown), respectively. Thus, stochastic forcing may play some role in shaping the ensemble average phase relation between the thermocline depth anomaly of the western Pacific and the SST anomaly in the eastern Pacific.

Although an ensemble of large size is needed to achieve the analytical spectra (3.6) through numerical simulations, nearly identical results of autocorrelations and cross correlation shown in Fig. 8 can be accurately calculated using the time series as shown in Fig. 7. This is because the autocorrelations and cross correlation are integrated quantities of the spectra and thus are robust. For the comparison with the observation, the autocorrelations and cross correlation of the NINO3 index and the sea level averaged at two islands, Truk and Pago Pago, where there are long records, are shown in Fig. 10. There are some differences between the results in Fig. 10 and in Figs. 8 and 9; for instance, the period exhibited in the autocorrelations of Fig. 10 is about 4 years, only about a half year longer than that in Figs. 8 and 9. The time delay indicated in the cross correlation in Fig. 10 is about 2 months, which is slightly shorter than that in Figs. 8 and 9. Although the results in Fig. 10 are subject to some error due to the short record of the dataset, they agree well with those in Figs. 8 and 9.

The analytical solutions of the spectra, autocorrelations, and cross correlations cannot be obtained in the supercritical regime where nonlinearity has to be taken into consideration for a bounded solution with the stochastic forcing. However, the same numerical integration method can be applied to obtain the numerical simulations and a typical solution is shown in Fig. 11a. This nonlinear limit cycle solution perturbed by the stochastic forcing is quite similar to the weakly damped linear oscillation shown in Fig. 7a under the same stochastic agitation. In fact, for this particular case, the autocorrelations and cross correlation (Fig. 11b) calculated from the time series in Fig. 11a are almost identical to those in Fig. 8a. It should be pointed out that the decay rate of the autocorrelations of the perturbed limit cycle solution by a stochastic forcing depends on the nonlinearity factor *e _{n}*. The smaller

*e*is, the slower the autocorrelations decay with the lag time, and the more regular the oscillatory solutions of the system become. Thus stochastic forcing works as a source of dissipation in the nonlinear model, and the dissipation rate depends upon both the amplitude of variance of the stochastic forcing and the degree of the nonlinearity of the system. The nonlinear factor used in Fig. 11 is close to the degree of the nonlinearity in a ZC-type model (e.g., BH). The similarity of the autocorrelations from the result in Fig. 11b to that in Fig. 8a illustrates the difficulty of distinguishing two views: the view of ENSO as a limit cycle perturbed by noise and the view of ENSO as a decaying oscillator sustained by noise. However, with the consideration of both the noise excitation and self-excitation, the basic BWCZ mechanism gives rise to the recharge oscillator in a wide parameter range that embraces the possible realistic values of the coupling coefficient. Although whether ENSO is actually self-excited or sustained by noise will be difficult to determine from the short records, the salient low-frequency cyclic nature of the ENSO can be explained by the coupled basinwide recharge oscillator of the tropical ocean–atmosphere system.

_{n}## 4. Discussion

### a. Zonal mean thermocline depth and recharge oscillator

The zonal mean equatorial thermocline depth is of essential importance in forming the coupled recharge oscillator. The strong dependence of the ENSO-like oscillation on the zonal mean equatorial heat content was first discovered by Zebiak and Cane (1987) in their numerical experiments with the ZC model. If the zonal mean thermocline depth in their coupled simulation is artificially removed in the calculation of SST change, the interannual variability in the model will completely disappear. Artificially reducing it by half and doubling it will result in significant changes in the characteristics of the interannual oscillations of the model. They found that the periods of the ENSOlike oscillation changed from about 4 years to about 5–6 years and about 2 years, respectively, for these two cases.

These sensitivity experiments can be analyzed with the conceptual coupled model (2.6) to elucidate the role of the zonal mean thermocline depth in the formation of the recharge oscillator. With (2.1), (2.5), and (2.7), the zonal mean thermocline depth can be expressed as

Following the experiments conducted by ZC, one replaces *h*_{E} by *h*_{E} − *λ**ĥ* in the SST equation, and then the modified coupled system (2.6) becomes

If *λ* = 1, the zonal mean thermocline depth is removed in the calculation of thermocline feedback to SST, and the SST equation becomes decoupled from the ocean dynamics equation. The recharge oscillator is completely destroyed. Instead, the system has a decaying ocean adjustment mode and a coupled nonoscillatory SST mode. If the coupling is sufficiently strong, the latter will become an unstable mode. When the nonlinearity is included, this unstable mode leads to new steady-state solutions. If *λ* = 1/2, half of the zonal mean thermocline depth is artificially removed. In this case, ocean adjustment dynamics and SST dynamics are still coupled. The critical coupling coefficient is at *μ*_{c} = 8/9, and the frequency at this coupling coefficient is *ω*_{c} = 2/33/32 for the otherwise same parameter choice as for the normal case (*λ* = 0). Comparing this with *ω*_{c} = 3/32 at *μ*_{c} = 2/3 for the normal case, there is indeed an increase of period of about 50% from about 3.4 years to about 5 years, corresponding to both frequencies at different critical coupling coefficients. Similarly, if one artificially doubles the zonal mean thermocline depth by setting *λ* = −1, the critical coupling coefficient will be *μ*_{c} = 4/9 and corresponding frequency will be *ω*_{c} = 3/32 14/9, which gives a period shorter than that of the normal case. As is shown in the appendix, the frequencies of nonlinear oscillations in the supercritical regime tend to be held close to the value at critical coupling. Thus this analysis is grossly in agreement with the results of the ZC experiments with a fixed coupling coefficient. Moreover, it clearly demonstrates the importance of the equatorial zonal mean heat content or zonal mean thermocline depth in the formation of the recharge oscillator.

### b. Recharge oscillator and delayed oscillator

It will be shown in Part II that the BH delayed oscillator and the conceptual recharge oscillator are different approximations to a slow coupled mode of the tropical ocean–atmosphere system. The main difference, however, is in the emphases on the phase transition mechanisms. The delayed oscillator requires an explicit time delay associated with the equatorial wave propagation and a strong western boundary wave reflection. Because the forced Rossby and Kelvin waves have opposite signs in their associated thermocline depth anomaly fields and the boundary reflection does not alter the signs in these fields, the locally forced Kelvin wave and the Kelvin wave remotely generated from the western boundary reflection of the forced Rossby wave have different signs in their associated thermocline depth anomaly fields. The contributions from these two parts to the thermocline depth in the equatorial eastern Pacific are both related to the local wind forcing. The reflected part has a time delay owing to the time taken for the Rossby wave to propagate from its forcing region to the western boundary and the reflected Kelvin wave to propagate from the western boundary to the central and eastern Pacific. The reflected Kelvin waves overtake at some later time the forced Kelvin waves in the eastern Pacific, which turns the thermocline depth associated with the forced Kelvin wave around. It is the western boundary reflection that accomplishes the phase transition in the delayed oscillator.

In the recharge oscillator model, the detailed wave propagation process is not emphasized. Equatorial wave dynamics are collectively viewed as an ocean adjustment process to redistribute mass and heat under changing wind stress forcing. The emphasis here is the quasi-Sverdrup balance near the equator and nonequilibrium of the zonal mean pressure field (or thermocline depth) under a slowly changing wind forcing. The equatorial wave dynamics quickly establishes the pressure gradient for maintaining the quasi-Sverdrup balance. For instance, a westerly wind anomaly is almost in quasiequilibrium with a positive tilt of thermocline anomaly. This provides the pressure gradient for a meridional geostrophic current away from the equator throughout the basin so as to drain the mass and heat content slowly out of the equatorial band. Thus, during a warm phase with westerly wind stress anomaly, the equatorial ocean zonal mean heat content is discharging. Similarly, a cold phase corresponds to the recharging of the equatorial ocean zonal mean heat content. It is this recharge–discharge process associated with basinwide ocean adjustment that leaves an anomalous shallow (deep) equatorial thermocline after a warm (cold) event that serves as the phase transition mechanism proposed by Cane and Zebiak (1985) and Wyrtki (1986).

In fact, the recharge oscillator embodies the BH delayed oscillator, despite the differences in the views on the transition mechanisms. The conceptual model (2.6) and its nonlinear version (3.1) can be rewritten into the same mathematical form as the BH delayed oscillator under some special conditions. As shown in the neutral solution (2.9) and nonlinear solution (3.2), anomalies of the thermocline depth of the western Pacific are negatively proportional to the SST of the eastern Pacific with a time delay. Substituting the thermocline depth into the SST variable with the time delay, the neutral solution (2.9) approximately satisfies the following delay equation:

Here, *b̂*= *α**b*_{0}*μ*_{c}*γ **μ*/*μ*_{c}; *R* and *η* are the same as in (2.6) and (2.9). Equation (4.3) takes precisely the same form of the linear delayed oscillator suggested by BH based on the numerical fitting to the ENSO-like oscillation in the ZC model. Although (4.3) and (2.6) share the same neutral solution (2.9) only when *μ* = *μ*_{c}, the leading eigen mode of (4.3) is a good approximation to the oscillatory mode in (2.6) for a wide range of coupling coefficient as shown in Fig. 12, even if the time delay *η* is artificially fixed in (4.3). The significant differences, comparing the results in Fig. 2 and Fig. 12, only occur for low values of the coupling coefficient.

The nonlinear solution (3.2) of the system (3.1) is also an approximate solution of the nonlinear form of the BH delayed oscillator,

in the supercritical regime near the Hopf bifurcation even if *η* is fixed. Here *e*^{′}_{n} = *e*_{n}/(*μ*_{c}*b*_{0})^{3} and *â* = *α /**ω*^{2}_{c} + *r*^{2}. In fact, one can always fit a harmonic oscillator with a delayed oscillator without much error in the solutions within some range of model parameters.

As with the recharge oscillator, the equatorial zonal mean thermocline depth also plays a crucial role in the delayed oscillator. One can perform the same sort of sensitivity experiments first suggested by ZC and further analyzed in the context of the recharge oscillator in section 4a and find a similar conclusion. It will be shown in Part II that in the framework of linear shallow-water dynamics, the recharge–discharge process of equatorial heat content associated with the ocean adjustment depends on the meridional mass transport by the forced Rossby wave and the boundary mass fluxes as a result of wave reflections at the east and west sides. The BH delayed oscillator can be interpreted using the recharge oscillator by considering the recharge–discharge process resulting from the forced Rossby waves and their western boundary reflections. This is what physically underlies the mathematical similarity between the delayed oscillator and recharge oscillator. Thus the physically relevant delay argument based on wave propagation makes the delayed oscillator a different and sound approximation in the neighborhood of the Hopf bifurcation of the recharge oscillator. The delayed oscillator model can be viewed both mathematically and physically as a special version of the recharge oscillator.

The recharge oscillator is conceptually more general and yet as simple. With the recharge oscillator as a conceptual model, one need not chase the details of the ocean wave propagation and their reflections in the ocean dynamics adjustment, but rather one can focus on the slow, large-scale buildup of the western Pacific warm pool or the equatorial mean thermocline depth. This may help us to diagnose and understand the slow ENSO-like oscillation in comprehensive coupled GCMs (e.g., Philander et al. 1992) and the slow ENSO signals in the data bypassing the detailed wave timescale features.

### c. Mixed SST–ocean dynamics mode and recharge oscillator

The conceptual recharge oscillation model only captures two nonoscillatory coupled modes, an SST mode and an ocean dynamics mode, at the fast wave limit [where time derivative of thermocline depth in (2.6) is ignored] and the fast SST limit (where the time derivative of the SST equation is ignored). The fast-wave limit approximation eliminates the recharge–discharge mechanism due to the slow ocean dynamical adjustment (Cane 1992b) and thus eliminates the mechanism for the recharge oscillator. The fast-SST limit simplification is not a good approximation either, because the coupling process greatly alters the uncoupled fast damping rate of SST so that |*R*| ≪ |*c*| and the time derivative in the SST equation cannot be ignored. Although there is no oscillatory coupled mode in these limits for the coupled system (2.6) owing to the simplifications involved, it gives birth to the recharge oscillation that relies on the merger of the SST mode and the basinwide ocean adjustment mode to form an intrinsically coupled and mixed SST–ocean dynamics mode away from the two extreme limits.

The details of the mixed SST–ocean dynamics oscillator may be much more complicated, as elucidated in JN with a picture of the mixed mode as coming from a series of mergers of oceanic modes with a nonoscillatory SST mode or a propagating SST mode (Jin et al. 1996). These complications are due to the fact that, in addition to the discrete ocean basin modes (Cane and Moore 1981), there exist numerous lower-frequency leaky scattering modes as a part of the spectrum of ocean adjustment modes of shallow-water dynamics (JN). The conceptual recharge oscillator model (2.6) provides the simplest and most constructive way to conceptually approach the complicated mixed oscillator. It reduces the complicated merging picture into a single merger of a decaying ocean adjustment mode with an SST mode. This simplification allows insight into the physics that were actually envisioned by the BWCZ hypotheses. Moreover, the recharge oscillator is at the heart of the mixed SST–ocean dynamics oscillator and allows us to understand the richer behavior of coupled ENSO modes because the recharge oscillator can be more heavily weighted by a propagating SST mode or by an ocean adjustment mode in a more complicated context where the processes for these modes are included. In other words, with the analytical recharge oscillator model, one can now understand the unified picture portrayed in JN from the most relevant central regime toward the extreme limits instead of the other way around. The simplicity of the recharge oscillator model for the mixed SST–ocean dynamics mode together with its generalization to the delayed oscillator model and the fact that the coupled wave oscillator at the fast wave limit is continuously related to the mixed mode favor it being a simple paradigm for describing the main physics of ENSO as hypothesized by BWCZ.

## 5. Summary

A new conceptual model is developed to describe the nature of the ENSO variability observed in the ocean–atmosphere system of the tropical Pacific based upon the BWCZ hypothesis. The model is cast in terms of the SST anomaly over the central to eastern part of the equatorial Pacific and the thermocline depth over the warm pool region of the basin. The SST anomaly is regulated by upwelling and thermocline feedback processes, which amplify the SST anomaly as suggested by Bjerknes (1969). The recharge–discharge process of the entire equatorial zonal mean heat content is controlled by the strengthening and weakening of the trade wind system through ocean dynamic adjustment. This zonal mean heat content is out of phase with the SST anomaly, which results in further change in the trade winds. Equatorial waves actively keep the equatorial Sverdrup balance in check, whereas both the equatorial wind stress and the zonal mean thermocline depth determine the thermocline depth in the eastern equatorial Pacific, the essential regulator of the eastern Pacific SST. The out of phase relation between the zonal mean thermocline depth and the SST anomaly leads to a basinwide harmonic recharge oscillator. The basic characteristics of the recharge oscillator are in agreement with observations and the results of ZC model.

The recharge oscillator model is the simplest physical model for ENSO based on the BWCZ hypothesis in representing the amplification processes for SST and the memory of subsurface ocean in the basinwide dynamical adjustment. The latter resides in the zonal mean thermocline depth that is not in quasi-equilibrium on the ENSO timescale. The ocean–atmospheric coupling ties these processes together to form the recharge oscillator. The growth rate of this oscillator largely depends on the positive feedback processes for the SST anomaly, and its period depends on the timescale of ocean adjustment, the timescale of SST dynamics, and the coupling parameters. Over a wide range of the relative coupling coefficient, this recharge oscillator can be either self-excited or stochastically sustained with a robust period of about 3–5 years. The nonlinearity in (3.1) does not alter the characteristics of the oscillation near the Hopf bifurcation, although it extends the parameter range of the recharge oscillation. On the other hand, the ensemble average phase lag between the SST of the eastern Pacific and the thermocline depth in the western Pacific is shorter than the phase lag (about 6 months) determined by the linear eigen mode of the system, depending on the sources of stochastic forcing. Thus, for the stochastically sustained recharge oscillator, the time delay between the thermocline depth and SST depends on the dynamics of the coupled model and the noise forcing as well.

Unlike the delayed oscillator model, the recharge oscillator model ignores the explicit role of wave propagation. It is more general in treating the subsurface ocean memory in the adjustment dynamics. This recharge oscillator can be virtually fitted by the delayed oscillator proposed by BH near the bifurcation point. The latter can be interpreted physically as a version of the recharge oscillator. Furthermore, the recharge oscillator model also gives the maximum simplification to the more complicated framework for the tropical interannual variability in JN based on the detailed analysis of connection of various regimes in parameter space in relation to the mixed SST–ocean dynamics mode. This simple recharge oscillator model, with its physical clarity and conceptual generality in describing BWCZ mechanisms, serves as a better paradigm for ENSO.

Recent studies (Jin et al. 1993, 1994, 1996; Tziperman et al. 1994; Tziperman et al. 1995; Chang et al. 1995; Wang and Fang 1996) have suggested two possible scenarios for the ENSO irregularity and thus also the intrinsic limit of ENSO predictability; namely, the deterministic chaos that results from the interaction of the annual cycle in the climatological basic state with the interannual ENSO oscillator and the stochastic forcing of uncoupled “weather noise.” The new conceptual model not only provides a simple paradigm for the understanding of the regular salient cycling feature of ENSO, it can also be further used as a prototype model for gaining insight into these more challenging issues of ENSO. The role of noise forcing in shaping some features of ENSO is worth further study in gaining better understanding of, for instance, the interaction of monsoonal variability and ENSO and its predictability.

## 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. The author enjoyed stimulating discussions with Bin Wang and Roger Lukas. Thanks go to Eli Tziperman and Edward Schneider for their critiques that led to the improvement of this paper, Thomas Schroeder and Diane Henderson for their careful reading and editing of the manuscript, Stephen Zebiak for his assistance in running the ZC model, Gary Mitchum for providing the sea level data, and Xiaowei Sun for producing some of the figures. This is SOEST Contribution Number 4177.

## REFERENCES

### APPENDIX

#### Nonlinear Solutions of the Coupled Model

##### Weakly nonlinear solutions

When the relative coupling coefficient is slightly supercritical,

system (3.1) has a weak growth rate on the order of Δ. As in the general case (e.g., Iooss and Joseph 1990; Jin and Ghil 1990), near the Hopf bifurcation this slow growth will be bounded by the nonlinearity to allow an amplitude of its solution to be on the order of Δ. The weakly nonlinear solutions of the nonlinear system (3.1) near the Hopf bifurcation can be solved using a standard perturbation method, for instance, the multiple-timescale method (e.g., Pedlosky 1986). Introducing a new independent variable, *ζ* = *t*Δ, to characterize the slow growth timescale, and expanding the solution in a power series ordered by Δ, one can express the time derivative and variable as follows:

Substituting (A.1) and (A.2) into (3.1), one can find that the balance of order Δ gives the following linear equations:

Here, the fact that *R* = *r* at the critical coupling *μ* = *μ*_{c} is used. This set of equations is the same as linear system of (2.6) with the coupling coefficient *μ* = *μ*_{c}. Its solution can be written as

where *ω*_{c} is the frequency at the critical coupling, and c.c. denotes the complex conjugate of the term in the front of it. In this leading order solution, the amplitude factor *B* yet needs to be determined.

The next order (Δ^{3/2}) balance leads to the following equations:

where

The inhomogenous terms *G*_{T} and *G*_{ h} will result in following secular terms in (A5):

and the solvability condition requires that

This solvability condition yields the following equation for the amplitude factor of the leading order solution (A4):

The real and imaginary parts of the first term on the right-hand side are simply linear corrections to the growth rate and frequency related to the supercriticality in the coupling coefficient, and the real and imaginary parts in the second term are the nonlinear reduction to the growth rate and nonlinear correction to the frequency. Equation (A.9) has a simple limit cycle solution as follows:

Combining (A.4) and (A.10), one can rewrite the solution in the form as expressed in (3.2). It should be noticed that although the relative amplitude of the leading order solution (A.1) depends on Δ, its absolute value depends on the sensitivity of subsurface temperature to the thermocline depth in the parameterization. Furthermore, the nonlinear contribution significantly cancels the linear reduction of frequency in the supercritical regime so that the periods of the nonlinear solutions are close to the value at the Hopf bifurcation. A comparison of the results of the analytical solution with the numerical integration of Eqs. (3.1) is given in Fig. A1. Clearly, in the neighborhood of the Hopf bifurcation, the approximate analytical solution agrees well with the numerical solutions.

##### Fully nonlinear solutions

Nonlinear system (3.1) has two kinds of limit solutions: limit cycles and steady-state solutions (or limit points). In fact, all steady-state solutions of (3.1) can be obtained analytically. When *μ* < 16/15, the zero solution (or the climate state) is the only steady-state solution. This state becomes unstable when *μ* > 2/3, and it gives rise to the self-excitation of the recharge oscillator through the supercritical Hopf bifurcation, as analyzed in section Aa. When *μ* > 16/15, two additional steady-state solutions occur through a Pitchfork bifurcation and solutions can be expressed as

They represent a warm and a cold state respectively. Although these steady states depend on the coupling, one can analyze the stability of these new steady states analytically. The dependence of the eigenvalues of the linearized system with respect to these states on the coupling is of similar character to that in Fig. 2 but at a much stronger coupling, as shown in Fig. A2. As a matter of fact, the physics is also the same as the recharge oscillator. With the new steady states as the basic states, the thermocline feedback factor *γ* is reduced to *γ̃* = *γ* − 3*e*_{n}(*h*^{(2,3)}_{E})^{2}. The linearized equations of (3.1) with respect to these two steady states are identical. They are the same as Eqs. (2.6), except that the thermocline feedback factor *γ* becomes *γ̃*. A stronger coupling is necessary to effectively combine the SST mode and ocean adjustment mode into an oscillatory mode because *γ̃* is smaller than *γ*. These steady states are stable when *μ* > 19/15 and unstable when *μ* < 19/15. The frequency at the critical point *μ*_{pc} = 19/15 is also 3/32, the same as that at the Hopf bifurcation from the climate state at *μ*_{c} = 2/3. However, the instability near *μ*_{pc} = 19/15 does not generate limit cycles around the steady states through supercritical Hopf bifurcations; instead, it leads to the same limit cycle solution originated from the Hopf bifurcation of the climate state through a global bifurcation because the orbit of this limit cycle passes nearby the steady states. The limit cycle from the first Hopf bifurcation is the unique stable limit solution of the system between 2/3 < *μ* < 19/15 because all the three steady-state solutions are unstable and no other limit cycles exist. When *μ* > 19/15, the limit cycle collapses and the system settles to either a warm or a cold steady state. The whole bifurcation tree is illustrated in Fig. A3. The recharge oscillator exists far beyond the linear oscillatory regime. This robustness is of significance because it effectively avoids the drift to a spurious warm or cold state accompanied by the flux correction built into anomalous coupled models such as (3.1).

## Footnotes

*Corresponding author address:* Dr. Fei-Fei Jin, Dept. of Meteorology, University of Hawaii at Manoa, 2525 Correa Road, HIG 331, Honolulu, HI 96822.