## Abstract

El Niño–Southern Oscillation (ENSO) is driven by large-scale ocean–atmosphere interactions in the equatorial Pacific and is sensitive to change in the mean state. Whereas conceptual models of ENSO usually consider the depth of the thermocline to be influential on the stability of ENSO, the observed changes in the depth of the 20°C isotherm are rather weak, on the order of approximately 5 m over the last decades. Conversely, change in stratification that affects both the intensity and sharpness of the thermocline can be pronounced. Here, the two-strip conceptual model of An and Jin is extended to include three parameters (i.e., the contribution of the first three baroclinic modes) that account for the main characteristics of the mean thermocline vertical structure.

A stability analysis of the model is carried out that indicates that the model sustains a lower ENSO mode when the high-order baroclinic modes (M2 and M3) are considered. The sensitivity of the model solution to the coupling efficiency further indicates that, in the weak coupling regime, the model allows for several ocean basin modes at low frequency. The latter can eventually merge into a low-frequency and unstable mode representative of ENSO as the coupling efficiency increases. Also, higher baroclinic modes project more energy onto the ocean dynamics for the same input of wind forcing. Therefore, in this study’s model, a shallower, yet more intense mean thermocline may still sustain a strong (i.e., unstable) and low-frequency ENSO mode. Sensitivity tests to the strength of the two dominant feedbacks (thermocline vs zonal advection) indicate that the presence of high-order baroclinic modes favors the bifurcation from a low-frequency regime to a higher-frequency regime when the zonal advective feedback is enhanced.

It is suggested that the proposed formalism can be used to interpret and measure the sensitivity of coupled general circulation models to climate change.

## 1. Introduction

El Niño–Southern Oscillation (ENSO), the well-known climate fluctuation, is driven by large-scale ocean–atmosphere interactions in the equatorial Pacific. It has dramatic impacts in many regions surrounding the Pacific Ocean. In recent decades, many efforts have been dedicated to the observation and modeling of ENSO. In particular, the Tropical Ocean and Global Atmosphere (TOGA) program (McPhaden et al. 1998) has led to the development of a wide variety of models and paradigms that allowed improving our understanding of this phenomenon (Neelin et al. 1998). From the early work of Bjerknes (1969), who first observed positive ocean–atmosphere feedbacks leading to El Niño (and oppositely to La Niña) growth, various ENSO theories were proposed that consider negative oceanic feedbacks associated with the decay and phase transition of such events. Those negative oceanic feedbacks were based on, among others, the delayed effect of planetary equatorial wave propagation on SST anomalies (Battisti and Hirst 1989; Schopf and Suarez 1988) and/or on the oceanic heat transfer from the equatorial region toward the off-equatorial regions and vice versa [cf. recharge–discharge oscillator of Jin (1997a,b)]. The oscillatory nature of ENSO could therefore be linked to the coupled nature of the equatorial ocean–atmosphere system, in which the ocean provides the memory of the system. In a simplified framework, this has been addressed through the theorical formalism of basinwide coupled instabilities (Neelin et al. 1994), which permits the exploration of the phase space of the ENSO mode. In particular, this approach permitted establishing that ENSO (characterized here by a growth rate and an oscillation frequency) is very sensitive to the characteristics of the oceanic mean state (An and Jin 2000; An and Wang 2000; Fedorov and Philander 2001). The latter actually controls the relative strength of the dominant feedbacks in the tropical Pacific—namely, the thermocline feedback, which relates the thermocline depth anomalies to SST anomalies [dominant in the Niño-3 region (5°S–5°N, 150°–90°W)]; and the zonal advective feedback, which designates the rate of SST changes along the equator as mostly controlled by anomalous zonal advection of mean SST [dominant in the Niño-4 region (5°S–5°N, 150°E–150°W)]. This conceptual view has been particularly fruitful for interpreting more complex systems (Mechoso et al. 2003), which include the coupled general circulation models (CGCMs) of the Coupled Model Intercomparison Project phase 3 (CMIP3) (Meehl et al. 2007), and their sensitivity to change in mean state. In particular, the diversity of the mean state characteristics in these models goes along with a wide range of ENSO dynamics (Van Oldenborgh et al. 2005; Guilyardi 2006; Belmadani et al. 2010). For instance, in these models, a cooler SST mean state tends to be associated with a larger (reduced) contribution of the zonal advective (thermocline) feedback strength (Belmadani et al. 2010). Changes in ENSO variability in recent decades (Yeh et al. 2009; Lee and McPhaden 2010) also suggest that a warmer state may influence the ENSO dynamics through the change in the balance between the two feedbacks. However, the diversity of the sensitivity of models to the increase in greenhouse gases (An et al. 2008; Yeh et al. 2010) makes it difficult to identify which parameter change of the mean state is critical. An et al. (2008) show that the change in the vertical contrast—that is, the change in difference between the subsurface temperature and SST—turns out to be more influential on ENSO than changes in the mean SST itself. This indicates that the vertical ocean stratification is a key parameter for understanding the sensivity of ENSO to global warming (Yeh et al. 2010).

In linear conceptual models of ENSO [Jin 1997b; An and Jin 2001, the two-strip Jin model (JN model); Fedorov and Philander 2001)], three parameters usually define the mean state: 1) the mean upwelling; 2) the mean zonal SST gradient across the equatorial Pacific; and 3) the mean thermocline depth, which constrains the wave dynamics. Although the mean thermocline depth (as inferred from the depth of the 20°C isotherm) determines to some extent the characteristics of the mean stratification, paradoxically, observations and oceanic reanalyses reveal a small change in the mean thermocline depth at decadal-to-interdecadal time scales (Wang and An 2001; Dewitte et al. 2009). Conversely, change in stratification in the vicinity of the thermocline can be significant, reflecting the rich vertical structure variability of the equatorial Pacific (Dewitte et al. 1999; Moon et al. 2004). In particular, a slight flattening thermocline, as experienced after the 1970s, is associated with an increased stratification in the central equatorial Pacific that corresponds to a 1.4°C decrease at 150 m, 180° (cf. Fig. 1 of Dewitte et al. 2009). In CGCMs, increased stratification may actually lead to an ENSO amplitude rise, although the thermocline flattens and thus reduces the Bjerknes feedback (Dewitte et al. 2007). Although the change in ENSO may result from a rectified effect by the change in mean state (Jin et al. 2003; Rodgers et al. 2004; Dewitte et al. 2007), it has been argued that the stability of ENSO may be altered in a warmer climate, meaning that the balance between the thermocline and zonal advective feedbacks is modified (Yeh et al. 2009; Collins et al. 2010). Before addressing the climate change problem, it seems necessary to first better understand the sensitivity of ENSO to the representation of the thermocline in a linear framework. This is the objective of this paper, which builds upon the works by Jin (1997b) and An and Jin (2001) to document the stability of ENSO in a more realistic framework. In particular, it proposes an extension of the Jin two-strip model, which takes into account the peculiarities of the equatorial Pacific thermocline—namely, its depth, thickness (e.g., sharpness), and intensity—which had not been considered explicitly in the original Jin model. The latter are introduced in the model in the form of the contribution of the first three baroclinic mode structure functions. In that sense the model can be considered as a reduced and linearized version of the model by Dewitte (2000), in the same manner that the Jin two-strip model is a reduced and linearized version of the model by Zebiak and Cane (1987). The background motivation of this study is to provide a theoretical framework within the linear approximation for the understanding of the impact of climate change on ENSO.

The paper is organized as follows: section 2 provides a description of the model. Section 3 is dedicated to a stability analysis of the model and documents the sensitivity of the ENSO mode to several key parameters accounting for the mean state change. Section 4 is a discussion followed by concluding remarks.

## 2. Model description

The model used in this study is called the Thual–Dewitte (TD) model. It can be viewed as an extension of Jin’s (1997b) two-strip model. It is based on similar physics, except that it has a refined representation of the mean thermocline through the consideration of several baroclinic modes (instead of one mode in the Jin model). It is thus composed of six equatorial strips extending from 125°E to 80°W for the oceanic component. Like the Jin model, only one equatorial strip is considered for the mixed layer model, and a statistical relationship is considered for the atmosphere. Also, the TD model can be considered as a reduced and linearized version of the model by Dewitte (2000), in the same manner the Jin two-strip model is a reduced and linearized version of the model by Zebiak and Cane (1987). This implies that nonlinear advection is omitted in the formulation of the model and that entrainment temperature is symmetrical and linear with respect to thermocline depth anomalies. The TD model also does not consider a seasonally varying mean state, unlike the Dewitte (2000) model. Finally, the Dewitte (2000) model considers an infinite number of meridional Rossby modes, whereas the TD model considers only the first meridional Rossby mode. All model parameters were derived from the Simple Ocean Data Assimilation (SODA) reanalysis over the period 1958–2007 (Carton and Giese 2008), and anomalies are considered with respect to the mean over this period. The following details each component of the model.

### a. Ocean dynamics

We consider an ocean at rest in hydrostatic balance with a reference state that is described by the mean density *ρ*_{0} and Brunt–Väisälä frequency *N*^{2}. The Brunt–Väisälä frequency (or mean buoyancy) is proportional to the vertical gradient of density and therefore defines various characteristics of the mean thermocline, such as its depth (being the depth where *N*^{2} is maximal), thickness, for example, sharpness (being the vertical extension of the *N*^{2} profile around the mean thermocline depth), and intensity (being the value of *N*^{2} at the mean thermocline depth). Although the mean stratification in the equatorial Pacific exhibits important zonal variations, here a mean buoyancy vertical profile *N*^{2}(*z*) typical of the central Pacific (170°–120°W) is chosen because the linear formalism cannot easily account for the zonal changes in stratification. Note that the central Pacific (170°–120°W) is the area that mostly controls ocean dynamics in our model, as the wind stress forcing is maximal in such an area. For clarity, *N*^{2}(*z*) will be referred to as the mean stratification. From *N*^{2}(*z*), a basis of vertical functions *F _{n}* with associated speed

*c*can be derived assuming the separation into horizontal and vertical variables as follows:

_{n}where (*u*, *υ*) are the anomalies of zonal and meridionnal currents, respectively; *p* is the pressure; and (*X _{z}*,

*Y*) are the zonal and meridionnal wind stress vertical derivatives, respectively. The ocean dynamics thus reduce to an equation of evolution for each vertical mode on the equatorial

_{z}*β*plane as shown:

In this paper, the focus is on the long wavelength–low-frequency equatorial waves (Kelvin and Rossby waves), which are nondispersive. The meridional momentum equation is thus approximated by the geostrophic balance, and therefore the meridional wind stress contribution is neglected. The system (1) can be projected onto the basis of the theoretical meridionnal functions associated with an infinitely unbounded ocean (Cane and Sarachik 1976), and the contribution of the baroclinic modes to zonal currents and pressure reads as follows:

with . The *q _{n}*

_{0}and

*q*are the projection of the

_{nm}*n*th baroclinic mode on the Kelvin and

*m*th meridional Rossby wave structures, respectively, and the

*ψ*terms are the Hermite functions. A differential equation for the ocean dynamics can then be written for each baroclinic mode

_{m}*n*:

where we use and

In the absence of wind forcing, the first equation (*m* = 0) accounts for the free eastward propagation of Kelvin equatorial waves with phase speed *c _{n}*, while the second equation (

*m*> 0) accounts for the free westward propagation of Rossby waves with phase speed −

*c*/(2

_{n}*m*+ 1). Odd (even) Rossby waves have a meridional symetric (asymetric) structure. A damping is introduced as a Rayleigh-type friction with

*e*-folding rate

*ɛ*for each vertical mode following Dewitte (2000). Furthermore, a perfect reflection of the waves at the meridional boundaries is considered that leads to an explicit formulation of the reflection efficiencies (Boulanger and Menkes 1995). For simplicity, the reflection efficiency depends only on the horizontal modes

_{n}*m*. We assume an infinitely extended wall in the north–south direction for the eastern boundary

*x*and a no mass flow condition for the western boundary

_{E}*x*. Therefore, the theorical values of reflection efficiency in the case of a Kelvin wave (

_{W}*m*= 0) and the first meridional-mode Rossby wave (

*m*= 1) can be derived as

where we use the conditions at *x _{E}* and at

*x*.

_{W}A constant mixed layer depth of 50 m is assumed, and for each vertical mode a projection coefficient *P _{n}* is introduced. By taking

*H*

_{0}the ocean bottom depth, we have

where , *τ _{x}* is the zonal wind stress anomalies, and

*h*

_{mix}is the mixed layer depth.

Therefore, by considering *n* vertical modes and *m* horizontal orders, the ocean dynamics can be depicted by *n* × *m* strips of variables *q _{nm}*. In our model, we consider the contribution of the first three vertical baroclinic modes (

*n*= 1, 2, 3). These modes are the most energetic propagating modes in the equatorial Pacific that best describe the wave dynamics (Dewitte et al. 1999, 2009). As in Kang and An (1998), only the contribution of the Kelvin and first-meridional Rossby waves (

*m*= 0, 1) are taken into account, considering that they explain a large amount of variance of the baroclinic modes (Dewitte et al. 1999, 2003). Therefore, the ocean dynamics is depicted by six strips of variables

*q*. The use of several baroclinic modes introduces slower-propagating waves and therefore lower frequencies in the model. One could alternatively increase

_{nm}*m*instead of

*n*to introduce slower-propagating Rossby waves (at speed ); however, it would keep the propagation speed of the Kelvin waves (at speed

*c*) unchanged, and the introduced Rossby waves would be dissipated relatively quickly because of their slower phase speed (shorter wavelength). The parameters of the ocean model are summarized in Table 1. A more detailed description of the ocean dynamics formalism used in this article can be found in Clarke (2008).

_{n}### b. Atmospheric component

A statistical relationship between SST anomalies along the equator *T _{E}* and

*τ*between 11°S and 11°N is used to build the atmospheric model. Here SST anomalies along the equator are calculated from a meridional average between 1°S and 1°N and therefore depend only on longitude and time. The statistical relationship is derived from the singular value decomposition (SVD) (Bretherton et al. 1992) between

_{x}*T*and

_{E}*τ*of the SODA data during 1958–2007, which leads to the modes

_{x}*e*(

_{k}*x*) and

*f*(

_{k}*x*,

*y*) associated to the time series

*α*and

_{k}*β*, such that

_{k}We retain the first two SVD modes (*k* = 1, 2) that account for 66% and 15% of the covariance, respectively. The first (second) SVD mode accounts for 70% (28%) of the total variance for the SST and 31% (30%) of the wind stress. The first SVD mode captures a remote interaction between eastern positive SST anomalies and central positive wind stress anomalies, while the second SVD mode captures a more local and delayed interaction between central positive SST anomalies and positive (negative) wind stress anomalies slightly west (east) of the SST anomalies (not shown). Therefore, in our model the atmospheric relationship reduces to

where is a regression coefficient defined as , where the angle brackets here denote the time average over the entire length of the used sample, 1958–2007; and at a given time is derived by projecting SST anomalies onto the SST eigenvector *e _{k}*. The atmospheric relationship is somewhat similar to the one of An and Jin (2000), but it is based on a SVD rather than an empirical orthogonal function (EOF) decomposition. As in An and Jin (2000) and other comparable studies, a coupling coefficient (or efficiency)

*μ*is introduced to explore the sensitivity of the coupled instabilities to the strength of the air–sea interaction;

*μ*= 1 corresponds to realistic ocean–atmosphere coupling at interannual time scale. Here, the projection of wind stress forcing patterns (

*f*

_{1}and

*f*

_{2}) on the equatorial waves meridional profiles is only slightly dependent on the baroclinic mode order (not shown). This is consistent with the fact that the atmospheric Rossby radius is much larger than the oceanic Rossby radius

*L*of the gravest baroclinic modes (assuming that the zonal wind stress anomalies along the equator can be fitted to a Gaussian profile in the meridional direction). This implies that the forcing term of the right-hand side of (2) in section 2a can be approximated by and for the Kelvin and first-Rossby equation, respectively, where

_{n}*F*and

_{K}*F*are longitudinal profiles independent of the baroclinic mode order that is considered.

_{R}### c. Mixed layer model

An equatorial strip is used for the thermodynamic equation controlling the rate of SST anomaly changes. The formulation is identical to An and Jin (2001), but here the equatorial strip extends from 1°S to 1°N. The *T _{E}* anomalies, thermocline depth anomalies

*h*, and zonal surface currents anomalies

_{E}*u*are averaged in the equatorial strip and therefore depend only on longitude and time. For

_{E}*T*, we have

_{E}where

Here *ɛ _{T}* is a thermodynamical damping coefficient of 125 days

^{−1}, is the mean upwelling,

*h*

_{mix}is constant (=50 m),

*γ*is an efficiency factor associated with the vertical thermal advection, and is the mean zonal gradient of SST along the equatorial Pacific. We derive

*h*and

_{E}*u*from ocean dynamics (see section 2a). Following Dewitte (2000),

_{E}*h*and

_{E}*u*are estimated as

_{E}with *F _{n}*(0) = 1. The angle brackets denote a meridional average in the equatorial strip and is the mean depth of the equatorial thermocline. The first, second, and third terms on the right-hand side of (3) correspond to the damping feedback, thermocline feedback (FT), and zonal advective feedback (FZA), respectively, in the equatorial strip. The thermocline feedback accounts for the anomalous advection of subsurface temperature anomalies by the mean upwelling, while the zonal advective feedback represents the anomalous advection of mean SST by zonal currents. Thus, FT and FZA depend on the characteristics of the mean state of the equatorial Pacific. The values of those parameters are derived from SODA. They are plotted in Fig. 1 in a nondimensionnal form to ease the comparison with previous studies. To derive (mean vertical velocity), the linear model of Dewitte (2000) is forced by the mean climatological wind stress from SODA with wave parameters derived from the vertical model decomposition of the mean stratification of SODA (see Table 1). From this model a mean upwelling of = 0.99 m day

^{−1}is estimated by zonal averaging over the domain 140°–90°W. We use an empirical function connecting

*h*(defined as the depth of the 20°C isotherm) and the subsurface temperature

*T*

_{SUB}(the temperature at

*h*

_{mix}) as (Zebiak and Cane 1987; Battisti and Hirst 1989), where the coefficients

*A*and

*B*are estimated by the least squares method from the SODA data. After the linearization of the derivative , the efficiency factor for the thermocline feedback is thus

*γ*(

*x*) =

*A*(

*x*)

*B*(

*x*). Such a formulation is different from An and Jin (2000), who further introduce a constant coefficient as

*γ*=

*AB*cosh

^{−2}(

*Bh*

_{0}), where

*h*

_{0}= 70 m. However

*γ*(

*x*) is around 0.15°C m

^{−1}in the eastern Pacific, which is in good agreement with the obtained values in An and Jin (2000). A Savitsky–Golay filter (Savitzky and Golay 1964) is used to smooth the thermocline feedback profile and to calculate the derivative in the zonal advective feedback term.

## 3. Stability analysis of the simulated ENSO

A stability analysis of the TD model is performed and the results are analyzed in this section. Because the mean stratification influences both the mean circulation and the variability, its impact on the ENSO mode cannot be inferred in a straightforward way. Here, we focus on interpretating the model’s behavior in light of previous comparable theoretical studies, as well as depict the model’s sensitivity to the mean stratification. First, for a fixed mean stratification, we document the dependency of the ENSO mode to the strength of the ocean–atmosphere coupling, as in previous comparable studies (sections 3a–c). Then we document the dependency of the ENSO mode to the projection of wind stress forcing on the different baroclinic modes (section 3d). Finally, we document the dependency of the ENSO mode to the relative strength of the thermodynamical feedbacks (section 3e).

### a. Frequency and growth rate

The solutions of the TD model are calculated numerically. It consists of seeking the eigenvectors and eigenvalues of the Jacobian matrix built from the linear differential equations described above (see section 2), where the vector (*q _{nm}*,

*T*) is the unknown of our system. The eigenvalues provide the growth rate and frequencies of the leading coupled modes. They are plotted in Fig. 2 for different values of

_{E}*μ*, ranging from 0 to 1.5. Only the eigenmodes with a growth rate greater than −1.8 yr

^{−1}and a frequency less than 1.5 yr

^{−1}are plotted. The results indicate that the TD model is able to sustain a variety of oscillatory and stationary modes that can be stable (negative growth rate) or unstable (positive growth rate).

In the uncoupled case (*μ* = 0), the oscillating eigenmodes consist of ocean basin modes (OBMs). In our model with no ocean–atmosphere interaction (*μ* = 0), each baroclinic mode is independent of the others (there is no mode interference) and has associated Kelvin and Rossby-1 waves, which propagate freely at speeds *c _{n}* and

*c*/3, respectively. Therefore, each baroclinic mode produces an ocean basin mode of its own, whose frequency is consistent with the time for a Kelvin wave and a Rossby-1 wave to reflect and travel across the basin [the frequency is (

_{n}*L*/

*c*+ 3

_{n}*L*/

*c*)

_{n}^{−1}, where

*L*is the length of the basin]. The ocean basin mode frequency is 1.2, 0.73, and 0.5 yr

^{−1}(see Table 1) for the first, second, and third baroclinic modes (M1–M3), respectively. In Fig. 2, we identify OBM1, OBM2, and OBM3 (at

*μ*= 0), which are associated with the baroclinic modes M1, M2, and M3, respectively. Two other ocean basin modes can be observed in Fig. 2 in the uncoupled case (marked by asterisks), whose frequencies (1.46 and 1 yr

^{−1}) are the double of the frequency of OBM2 and OBM3 (0.73 and 0.5 yr

^{−1}, respectively). These solutions are generally associated with large zonal wavenumbers, which are difficult to relate to some aspects of the observed variability in the equatorial Pacific. They also remain secondary for the stability of the system in term of growth rate.

We now focus on the system’s behavior in the coupled case (*μ* > 0). Here, the baroclinic modes are no longer independent because of ocean–atmosphere interaction, and therefore the eigenmodes take into account the contribution of every baroclinic mode. For each *μ*, the leading eigenmode of the system is defined as the eigenmode with maximal growth rate. With an increasing coupling coefficient, the leading eigenmode switches from a nonoscillatory damped mode (*μ* ≤ 0.3) to an oscillatory eigenmode (0.3 ≤ *μ* ≤ 1.3) before breaking into two nonoscillatory unstable modes (*μ* > 1.3, not shown). Such behavior is similar to the one observed in the Jin two-strip model (An and Jin 2001), although for a different range of *μ*. Our focus is on the range of values for coupling efficiency where the leading eigenmode is oscillatory (0.3 ≤ *μ* ≤ 1.3), which we define as the ENSO mode of our model. In Fig. 2, the ENSO mode is approximately in the range of frequency and growth rate of 0.05 to 0.5 and −1 to 0.8 yr^{−1}, respectively. The frequency and growth rate of the ENSO mode are further displayed in Fig. 3 as a function of *μ* (0.3 ≤ *μ* ≤ 1.3). For *μ* = 1, the ENSO mode is unstable with a period of approximately 3.7 yr, which is consistent with the contemporary ENSO characteristics (Trenberth and Stepaniak 2001) and the dominant frequency of the ENSO mode in the Dewitte’s (2000) model (namely, 3.3 yr). In the following sections, we further analyze the spatial characteristics (eigenfunctions) of the ENSO mode of our model and its sensitivity to *μ*.

### b. Structure of the ENSO mode

In this section, we analyze the eigenfunctions of the ENSO mode of our model. They are presented for the SST, thermocline depth *H*, and zonal current *U* anomalies of the model’s equatorial strip in Fig. 4 for a value of coupling efficiency *μ* = 1. In such representation, all eigenfunctions are sinusoidal. All fields are dominated by a standing ENSO oscillation for a period of 3.7 yr. The warm phase is associated with positive SST that departs from the central Pacific (160°W) and propagates both eastward and westward. West of 160°E, the SST is of opposite sign because of the cooling effect of the zonal advective feedback in this region (cf. Fig. 1). Along the equator, analytical Kelvin (Rossby-1) wave solutions have a contribution that is larger on *H* (*U*) (Cane and Sarachik 1976). Figure 4 indeed suggests that the transition from warm to cold phase is ensured by Rossby-1 wave reflections at the western boundary (with positive *U*) into Kelvin waves (with negative *H*), consistent with the delayed oscillator theory (Schopf and Suarez 1988). Features relevant to the recharge–discharge theory (Jin 1997a,b) are also observed, which result from oceanic wave adjustment in the model: the warm water volume (WWV) (or zonal average of *H*) is maximal at the beginning of the warm phase (when *H* is positive all along the equator), thus leading the El Niño signal in the east by approximately a quarter of a period.

The characteristics of the ENSO mode in the TD model can be interpreted in light of the coupled instability theory, in which the *H* mode (*U* mode) refers to the instabilities arising from a dominant thermocline feedback (zonal advective feedback) process. The *H* mode and the *U* mode impact on SST can be observed in Fig. 4, where they are responsible for the eastward SST anomaly propagation from 160° to 90°W and the westward SST anomaly propagation from 160°W to 125°E, respectively. Here, the *H* mode is dominant, which produces maximal SST in the eastern Pacific. The feedback mechanisms of the two instabilities (the *H* mode and the *U* mode) are different in many ways. The *H* mode, associated with the thermocline feedback, arises from the merging of a damped SST mode and ocean adjustment mode (An and Jin 2001). Such a mode was originally discussed as the delayed oscillator mode by Schopf and Suarez (1988). The zonal advective feedback produces the *U* mode by destabilizing the gravest ocean basin mode (An and Jin 2001). It is related to the SST mode introduced by Neelin et al. (1994). Jin and An (1999) further demonstrated that both instabilities are dynamically linked via the geostrophic balance between the upper-ocean zonal current and the meridional gradient of the thermocline. In our model, because the geostrophic balance is not considered explicitly, the dynamical coexistence of the two instabilities has to emerge through a different process.

Our results indicate that the TD model is able to simulate a realistic ENSO mode at *μ* = 1 that is consistent with former theoretical studies (Jin 1997a,b; An and Jin 2001). Figure 4 can be compared to Fig. 8 of An and Jin (2001), which shows the results of the stability analysis of their model [model parameters are derived from the National Centers for Environmental Prediction (NCEP) oceanic reanalysis covering the period from January 1980 to December 1995]. Although the TD model and the Jin model are based on a different formalism for ocean dynamics, the eigenmodes characteristics are comparable. The TD model, however, produces more intense SST anomalies in the east because of a stronger thermocline feedback. The sensitivity of the TD model to *μ* is depicted in the next section.

### c. Sensitivity to coupling coefficient

In the following, as a first step, we document the sensitivity of the dynamical adjustment of the model to *μ*. Figure 5 presents the eigenvectors amplitude of SST, *H*, and *U* for the ENSO mode as a function of longitude and *μ* (0.3 ≤ *μ* ≤ 1.3). The results indicate that the TD model sustains distinct regimes depending on the coupling coefficient. The amplitude of SST is maximal in the eastern Pacific and shows little sensitivity to the coupling coefficient. The amplitude of *H* increases with the coupling coefficient in both the eastern and western Pacific. The amplitude of *U*, however, shows a different sensitivity, with strong zonal currents in the central Pacific for both a small coupling coefficient (*μ* ~ 0.3) and a strong coupling coefficient (*μ* ~ 1.3). First, for a small coupling coefficient (*μ* ~ 0.3), the ENSO mode is very close to OBM3 from the uncoupled case (see Fig. 2). As a result the model produces relatively strong zonal currents (*U* ≥ 0.3 m s^{−1}) in the central Pacific, which is identical to the feature of an ocean basin mode (An 2005). Second, for a large coupling coefficient (*μ* ~ 1.3), the ENSO mode has a period that is very large compared to the time delays induced by the propagation of equatorial Kelvin and Rossby-1 waves (see section 3a). Therefore, there is no apparent zonal propagation on the dynamical fields *U* and *H* (not shown). This absence of zonal propagations leads to an increase of *U* in the central Pacific through the contribution of Rossby-1 waves, which are locally forced by the atmosphere. This peculiar behavior of the system with a strong coupling coefficient is in contrast to the behavior of the Jin model which produces strong zonal currents only for a small coupling coefficient [see Fig. 5 of An (2005) for a comparison]. It is attributed to different ocean dynamics adjustments in the Jin model and the TD model. The Jin model assumes a geostrophic balance between the equatorial and the off-equatorial strips as well as a locally wind-driven component to derive zonal current anomalies along the equator (Jin and An 1999), whereas in the TD model zonal current anomalies are derived explicitly from the Kelvin and Rossby-1 wave contributions. Therefore, differences between the two models originate mostly from their sensitivity to *μ*. In particular, for a similar value of the coupling coefficient, the TD model sustains a lower-frequency ENSO mode than the Jin model, which is mainly attributable to the contribution of higher-order baroclinic modes. Details of the differences in sensitivity to *μ* of the TD model and the Jin model are provided in the appendix.

We now investigate the role of vertical structure variability in the model, documenting the dynamical adjustment with regard to the contribution of the baroclinic modes.

### d. Role of baroclinic modes

Considering the importance of high-order baroclinic modes (*n* = 2, 3) in setting the sensitivity of eigenmodes to the coupling efficiency described above, in this section we detail their specific role in the transition phase between leading instabilities and in setting the characteristic time scales of variability of the ENSO mode. The eigenvectors’ amplitude of *H* and *U* for M1–M3 is displayed in Fig. 6 as a function of longitude and *μ* (0.3 ≤ *μ* ≤ 1.3). Each baroclinic mode exhibits maximum amplitude for *H* in the eastern and western Pacific and maximum amplitude for *U* in the central Pacific. Interestingly, differences arise regarding the baroclinic modes’ sensitivity to coupling. For a small coupling coefficient (0.3 ≤ *μ* ≤ 0.8), the amplitude of M3 (for *U* and *H*) is dominant because the ENSO mode is close to OBM3 (see Fig. 2), which has a dominant amplitude on M3 (not shown). However, for a large coupling coefficient (0.8 ≤ *μ* ≤ 1.3), M2 is dominant both in terms of *U* and *H*. Interestingly, the amplitude of M1 remains overall secondary compared to the other baroclinic modes. This illustrates the importance of higher baroclinic modes in setting the characteristics of the ENSO mode in our model. First, although wind stress forcing projects similarly onto the first two baroclinic modes (*P*_{1} ~ *P*_{2}, see Table 1), M2 is overall the main contributor to the ENSO mode. Second, when varying *μ*, the baroclinic mode contribution is modified through change in the nature of the ENSO mode from an ocean basin mode regime to a strong coupling regime (see section 3c).

It is also worth investigating the inverse problem, which consists of determining how changes in the baroclinic mode contribution only may alter the nature of the leading eigenmode. To document such sensitivity, the solution of the TD model is sought for a wide range of values of the projection of wind stress forcing on each baroclinic mode *P _{n}* while the other parameters derived from the mean stratification remain constant. We introduce a new varying sensitivity parameter

*δ*while

*μ*is kept constant (

*μ*= 1). In such a test, the wind stress projection coefficients (see section 2a) vary as

*P*

_{1}(1 −

*δ*),

*P*

_{2}(1 +

*δ*/2), and

*P*

_{3}(1 +

*δ*/2) for the baroclinic modes M1, M2, and M3, respectively. At

*δ*= −2 (

*δ*= 1), the wind stress projects solely on M1 (M2 and M3) while we roughly conserve the wind stress energy input, defined here as the sum of

*P*

_{1}+

*P*

_{2}+

*P*

_{3}. It is worth noting here that M1 (

*n*= 1) accounts for the density difference between the upper layer and the bottom ocean, and therefore it mostly represents the mean thermocline depth, whereas higher-order baroclinic modes M2 and M3 (

*n*= 2, 3) stand for the mean thermocline sharpness and intensity (Dewitte et al. 2007, 2009). Simply put, a deeper thermocline tends to increase the first baroclinic mode contribution, while a sharper and more intense thermocline leads to an increased contribution of the second and third baroclinic modes. A stability analysis is then performed for each value of

*δ*in the range (−1.5 ≤

*δ*≤ 1), and the frequency and growth rate of the solutions are presented in Fig. 7. For an increasing parameter

*δ*(

*P*

_{2}and

*P*

_{3}favored over

*P*

_{1}), the ENSO mode is displaced toward a more unstable and low-frequency regime. For instance, a larger

*δ*corresponds to a sharper and shallower thermocline, which will tend to trap more wind stress forcing in the surface layers because of the finer meridional scale of variability of the high-order baroclinic modes. Note that in our model, a larger

*δ*is roughly equivalent to a larger

*μ*in terms of the growth rate and frequency of the simulated ENSO mode (by comparing Figs. 3 and 7). However,

*μ*is arbitrary in our model, while

*δ*can be related to the characteristics of the thermocline, and in that sense has a physical meaning.

This emphasizes the role of mean stratification in setting the stability characteristics of the ENSO mode.

### e. Sensitivity to thermodynamical feedbacks

Changes in mean state associated with natural variability or climate change have the potential to impact the ENSO dynamics through changes in the main ENSO feedbacks (An and Jin 2001). In the TD model, the mean stratification *N*^{2}(*z*) controls the wind stress forcing of ocean dynamics through the relative contribution of the baroclinic modes (see section 3d), but it also constrains some aspects of the mixed layer thermodynamics because it influences the strength of the thermocline feedback and the zonal advective feedback (see section 2c). For instance, mean upwelling conditions (which contribute to the thermocline feedback) depend on the characteristics of the mean thermocline. The more intense (i.e., thinner) or shallower the thermocline, the larger the for the same input of wind stress forcing. Also, the dynamical relationship between changes in the thermocline depth and changes in the subsurface temperature at the base of the mixed layer, described by *γ* in the thermocline feedback (see section 2c), is sensitive to the mean stratification. A shallower thermocline leads to a stronger *γ*, which is observed, for example, in the eastern equatorial Pacific in comparison to the central equatorial Pacific (see Fig. 1). Furthermore, a more intense thermocline leads to a stronger *γ*, as subsurface temperature anomalies become more sensitive to the vertical displacements of the isopycnes. Conversely, the sensitivity of the zonal advective feedback to the mean stratification is not straightforward. However, the mean zonal gradient of SST (see section 2c) is related to ocean–atmosphere heat flux adjustments, which are also constrained by the mean stratification characteristics.

In this section, we introduce the new sensitivity parameter *α*, which controls the ratio of the FZA over the FT in the SST equation

The mean stratification influences the relative strength of the zonal advective and thermocline feedbacks, thus the *α* parameter. However, we assume that *α* has a different effect on the ENSO mode than the other parameters (*c _{n}*,

*P*,

_{n}*F*) characterizing the mean stratification because it influences the ENSO mode through the modification of the mixed layer thermodynamics instead of the ocean dynamics. Therefore, we assume that

_{n}*α*can be modified to explore the phase space of the ENSO mode in the model while the other parameters (

*c*,

_{n}*P*,

_{n}*F*) remain fixed. Note also that this parameter has already been used to explore various regimes of the Jin model in previous studies (An and Jin 2000, 2001; An 2005). For

_{n}*α*= −1, only the thermocline feedback is at work, while for

*α*= 1, the rate of SST changes is mainly controlled by the zonal advective feedback. The case

*α*= 0 corresponds to the reference run described previously. We choose to work with the two parameters (

*α*,

*μ*) to explore the phase space of the ENSO mode as a function of the strength of the ENSO feedbacks. The analysis is similar to the one by An (2005) but in a multibaroclinic mode context.

The results of the stability analysis of the model in terms of growth rate and frequency are presented in Fig. 8 for different values of *α* and *μ*. High-frequency eigenmodes (>0.75 yr^{−1}) are not contoured, as they remain secondary for the stability of the system in term of growth rate. In Fig. 8 we identify two types of eigenmodes called EG2 and EG3. By definition, the eigenmodes of EG2 (EG3) depart in the uncoupled case from the OBM2 (OBM3) (see section 3a). With an increasing *μ*, the eigenmodes of EG3 are eventually displaced toward a low-frequency regime, while the eigenmodes of EG2 keep a somewhat similar frequency in the range of 0.55–0.73 yr^{−1}. This has implications for the stability of the ENSO mode of our model, defined here as the leading oscillating eigenmode of the system (see section 3a). For a negative *α* (FT favored), the ENSO mode belongs to EG3 and the stability is very similar to the previous case (Fig. 2). However, for strong positive values of *α* (*α* ≥ 0.7), the ENSO mode eventually switches (bifurcates) from EG3 to EG2 with significant changes of its growth rate and frequency. To better illustrate such a feature, the frequency and growth rate of the leading eigenmode (ENSO mode) are contoured in Fig. 9 as a function of *α* and *μ*. Here, the ENSO mode is alternatively in EG3 or in EG2, or is nonoscillating (NO, not contoured). The transition between EG3 and EG2 is marked by a slight jump in frequency and a local minimum of growth rate. Near such a transition, the second leading eigenmode has a similar growth rate compared to the first one (not shown); therefore, the system has two favored frequencies (around 0.5 and 0.6 yr^{−1}). It is worth noting that the ENSO mode can even switch to EG1 for very strong (hypothetical) coupling values and *α* ~ 1 (not shown). The representation in Fig. 9 gives insight into the system’s sensitivity to both *α* and *μ*. When *α* decreases (FT favored), the leading eigenmode is displaced toward a more unstable and low-frequency regime. For *α* near unity (FZA favored), the system’s frequency remains close to the frequency of OBM2 (0.73 yr^{−1}), independent of the value of the coupling efficiency parameter. Such a result is in agreement with An (2005), who suggests that the *H* mode (see section 3b) induced by the thermocline feedback is much more sensitive to *μ* than the *U* mode induced by the zonal advective feedback.

In the following, we focus on the spatial structure of the ENSO mode in the phase space (*α*, *μ*). The localization of SST anomalies depends mostly on *α* and very little on *μ*, with a favored thermocline (zonal advective) feedback producing maximal SST anomalies in the eastern (central) Pacific (not shown). Also, when *α* increases so does the amplitude of *H* and *U* while the thermocline pivot (defined here as the longitude of the minimum amplitude of *H* along the equator) is progressively displaced to the west, from 145°W at *α* = −1 to 175°E at *α* = 1 (not shown). These changes reflect changes in the WWV anomalies, which characterize the delayed response of the ocean with respect to the atmosphere in the coupled system (Meinen and McPhaden 2000). Here, the WWV anomaliess are defined as the zonal and meridional average of *H* in the equatorial strip of the model (1°N–1°S, 125°E–80°W). The amplitude of the WWV anomalies is displayed in Fig. 10a as a function of *α* and *μ*. When the ENSO mode is displaced toward a low-frequency regime (see also Fig. 9), the amplitude of the WWV anomalies decreases and the *H* anomalies reduce to a zonal seasaw mode or “tilt” mode, consistent with Clarke (2010). This is illustrated by the phase difference of *H* between the eastern Pacific (90°W) and the western Pacific (180°) (gray dashed lines on Fig. 10a). At low frequencies, the eastern and western Pacific are almost in phase opposition (180°), which leads to weak WWV anomalies. Also, the WWV anomalies are related to the recharge–discharge mechanism proposed by Jin (1997a). Jin’s (1997a) simple model for ENSO is based on two linear first-order ordinary differential equations for eastern SST and heat content (~WWV) anomalies. Therefore, in this model the WWV anomalies lead the eastern SST anomalies by a quarter of a period (90°) [see also Burgers et al. (2005) for a similar approach]. The phase difference between the eastern SST and WWV anomalies is contoured in Fig. 10b as a function of *α* and *μ* for the TD model. When the thermocline feedback is favored (*α* ~ −1), the TD model reproduces the recharge–discharge mechanism of Jin (1997a) with a phase difference between SST and WWV anomalies close to 90°. However, when the zonal advective feedback is favored (*α* ~ 1), the WWV and eastern SST anomalies are almost in phase, indicating that the recharge–discharge mechanism of Jin (1997a) is not at work. It is also interesting to note from Fig. 10a that when the recharge–discharge mechanism is operating in the TD model (*α* ~ −1), the amplitude of the WWV anomalies is smaller than when the recharge–discharge mechanism is not operating (*α* ~ 1).

Considering the sensivity of the recharge–discharge mechanism in the TD model to *α*, it is worth investigating how M1–M3 contribute to the WWV anomalies. For each baroclinic mode, the WWV anomalies (WWV* _{n}* for

*n*= 1, 2, 3) were estimated and analyzed (not shown), where the total WWV anomalies result from WWV = WWV

_{1}+ WWV

_{2}+ WWV

_{3}. The results indicate that first, WWV

_{1}leads WWV

_{2}, consistent with the faster oceanic adjustment time scale associated with M1 compared to M2. For the same reason, WWV

_{2}leads WWW

_{3}. Second, on average in the (

*α*,

*μ*) phase space, the second baroclinic mode mostly contributes to the total WWV anomalies (by almost 40%). As a result, the total WWV anomalies tend to be in phase with WWV

_{2}.

## 4. Discussion and conclusions

In this study, we propose an extension of the Jin two-strip model (Jin 1997b; An and Jin 2001) that includes the dominant characteristics of the vertical structure variability of the ocean stratification, in order to explore the sensitivity of the ENSO mode to change in mean stratification. While conceptual ENSO models usually rely on shallow-water dynamics in a two-layer framework, here we test to what extent the characteristics of mean stratification can impact the ENSO mode considering the contribution of the first three baroclinic modes. As with Jin (1997a,b), the model solution consists of a mode that shares many characteristics with the observed ENSO in terms of frequency and stability. Both elements of the recharge–discharge theory and the delayed oscillator theory are observed in the model, which result from oceanic wave adjustment. As in An and Jin (2001), the stability of the ENSO mode depends on *μ*, which controls the coupling efficiency between the ocean and the atmosphere. In the standard case (*μ* = 1), the dominant mode is unstable and oscillates at a period of 3.7 yr, which is consistent with the contemporary ENSO characteristics (Trenberth and Stepaniak 2001). Interestingly, for a similar value of the coupling coefficient, the TD model sustains a lower-frequency ENSO mode than in the Jin model that is mainly attributable to the longer oceanic adjustment time scales associated with the contribution of the higher-order baroclinic modes (*n* = 2, 3) in the TD model.

In the model the contribution of each baroclinic mode to the thermocline depth and zonal current anomalies is sensitive to the coupling coefficient. For the intermediary range of values of the coupling coefficient, the system is dominated by the contribution of the second baroclinic mode, while the contribution of the first baroclinic mode remains secondary. As with Dewitte (2000), this stresses the role of the second baroclinic mode in setting the variability time scales and stability characteristics of the ENSO mode. At weak coupling, however, the system is dominated by the contribution of the third baroclinic mode. Also, in a hypothetical case where the projection of wind stress forcing is accounted for essentially through higher-order baroclinic modes, the ENSO mode is displaced toward a more unstable and low-frequency regime. The latter suggests that a shallower, yet more intense mean thermocline may still sustain a strong (e.g., unstable) and low-frequency ENSO mode. We further introduce the sensitivity parameter *α*, which controls the ratio of the thermodynamical feedbacks in our model, namely, the thermocline and zonal advective feedbacks. Depending on the value of *α*, the system can produce different El Niño events with either the eastern or the central SST anomalies in the case of a favored thermocline or zonal advective feedback, respectively. For a favored zonal advective feedback, the frequency exhibits little sensivity to the coupling coefficient, as mentioned in An (2005). For a favored thermocline feedback, the recharge–discharge mechanism is at work with the second baroclinic mode controlling (in terms of phase and magnitude) the WWV anomalies.

Within certain ranges of parameters in our model, the consideration of high-order vertical modes allows for the existence of several ocean basin modes, which permits the bifurcation behaviors of the solution within the interannual frequency band. Thus, several leading modes with similar growth rate but a different frequency can coexist. This is consistent with the wide spectrum of ENSO variability. It suggests that the spreading of ENSO variability in the frequency domain may be explained to some extent within the linear formalism through the consideration of the peculiarities of the ocean stratification, without invoking chaotic behavior or nonlinear processes (Timmermann et al. 2003).

We now review a few limitations of our approach. Considering the critical role of high-order baroclinic modes in setting the ENSO mode characteristics and its sensitivity to *μ*, one may wonder about the implication of truncating to the first three baroclinic modes. Figure 11a shows the ENSO mode as simulated by the TD model in terms of frequency and growth rate, depending on the number of baroclinic modes considered, from 1 to 5. For *n* = 3 the results are identical to Fig. 2. As more baroclinic modes are considered, more energy is added to the system; therefore, the ENSO mode is slightly displaced toward a more unstable and low-frequency regime. Interestingly, the stability is quite singular in the case of *n* = 1, in which the ENSO mode remains in a high-frequency and strongly damped regime. The latter stresses the importance of higher baroclinic modes (*n* ≥ 1) for energy balance and dynamical adjustement of the tropical Pacific system (see also the appendix). However, for *n* ≥ 4, the higher-order modes’ contribution is essentially trapped in the upper layer and corresponds to the locally forced variability (Dewitte et al. 1999); therefore their representation may be inaccurate within the simple formalism of the TD model. Alternatively, the variability for *n* ≥ 4 can be taken into account considering a simple frictional equation (Blumenthal and Cane 1989) forced by (where *P _{i}* is the wind projection coefficient for the baroclinic mode

*i*and

*h*

_{mix}is the mixed layer depth), which represents the part of the wind stress forcing that does not project onto the gravest baroclinic modes (here, the first three baroclinic modes). This leads to an estimation of zonal Ekman currents as shown:

where *ρ*_{0} is the mean density, *e*_{ECK} is a damping rate coefficient of 1.5 days^{−1}, and *τ _{x}* is the zonal wind stress anomalies. For the sake of simplicity, Ekman currents were not considered in this study. Considering them in the model does not lead to qualitative changes of the results (not shown), which indicates that only the baroclinic modes having clear propagating characteristics (

*n*= 1, 2, 3) are relevant here. Similarily, one may wonder what is the implication of truncating to the first Rossby wave (

*m*= 1) in the TD model. An experiment is performed in which the contribution from the third and fifth Rossby waves (

*m*= 3, 5) is added to the model, and the results are displayed in Fig. 11b. Note that the Rossby waves for

*m*= 2, 4 are meridionally asymmetrical and are therefore omitted. The contribution from Rossby-3 and Rossby-5 waves increases the strength of zonal advective processes in the TD model, and as a consequence it leads to a slight increase of the frequency of the ENSO mode (see section 3e). The results, however, suggest that the contribution from higher-order Rossby waves is not critical for the stability of the TD model.

Another limitation is associated with the representation of SST in the model that omits the meridional variability scales, because, like in the Jin model, only one equatorial strip for SST is considered. Considering that the first baroclinic mode extends more meridionally than the higher-order baroclinic modes (in terms of the Kelvin wave), this may somehow lead to a greater contribution of the first baroclinic mode to SST change in the equatorial band (5°N–5°S) than considered here. The consideration of two or more strips for SST off the equator would, however, require considering the contribution of meridional advection, which may burden the formalism proposed here. Also, specific coupling between baroclinic modes and the atmospheric (SVD) modes may be at work in the real word (Yeh et al. 2001), which stresses the sensivity of ENSO to the meridional scale of variability of SST. Further study is required to document this process in a more realistic framework.

As a last limitation, it is worth mentioning the simplicity of the formalism with regard to the realism of the thermocline in terms of its east–west contrast. The thermocline slope (deep in the central Pacific and shallow in the eastern Pacific) is not accounted for in our model, which considers a zonally uniform mean stratification. Because of modal dispersion, a greater contribution of higher baroclinic modes is expected in the eastern basin (Dewitte et al. 1999). Although such an effect could be taken into account within the linear formalism used here (cf. Busalacchi and Cane 1988), it is not clear to what extent it may modify the sensitivity of the model solution to the mean stratification characteristics documented here. This would deserve further study.

Despite these limitations, the model presented here may contain the appropriate physics to document the interaction between mean state changes and ENSO variability. Fedorov and Philander (2001) already showed that the depth of the mean thermocline can modify the ENSO characteristics (growth and frequency) and that with a realistic value for the thermocline depth, the ENSO mode can be significantly modified. Here, two other ingredients that depict the thermocline structure—namely, its sharpness and intensity—are considered in a similar theoretical framework. This has implications for the study of the ENSO modulation in the contemporary climate. In particular, the vertical stratification is expected to increase during the future greenhouse warming (An et al. 2008; Yeh et al. 2009), which tends to increase the higher-order baroclinic mode contributions (Dewitte et al. 2009). To what extent change in stratification associated with global warming can impact the balance between feedback processes is an issue of debate in the community. At this stage, it is interesting to note that most Intergovernmental Panel on Climate Change (IPCC) models exhibit a high sensitivity to global warming with drastic change in stratification (Timmerman et al. 1999; An et al. 2008; Yeh et al. 2009, 2010). The extended version of the Jin model proposed here can provide a theoretical framework for interpreting such diversity in model behavior. This is the topic of future research.

## Acknowledgments

S. Thual was supported by Centre National de la Recherche Scientifique (CNRS) and Conseil Régional Midi-Pyrénées, Contract 022009. S.-I. An was supported by Grant RACS 2010-2601 from the Korea Meteorological Administration Research and Development program. The two anonymous reviewers are thanked for their constructive comments.

### APPENDIX

#### Sensitivity to the Coupling Coefficient: Comparison with the Jin Model

In this appendix we describe the differences in sensitivity to the coupling coefficient *μ* of the TD model and the Jin two-strip (JN) model (Jin 1997b; An and Jin 2001). Both models use similar ocean thermodynamics and statistics for the atmosphere; although, in the TD model, the sea surface temperature and wind stress anomalies are statistically related through singular value decomposition (SVD), rather than empirical orthogonal function (EOF) formalism as in An and Jin (2000). The two models mostly differ in their ocean dynamics. The Jin model assumes a geostrophic balance between the equatorial and the off-equatorial strips to derive zonal current anomalies along the equator, whereas in the TD model zonal current anomalies are derived explicitly from the Kelvin and Rossby-1 wave contributions. Also, the JN model proposes a collective view of the equatorial Rossby waves propagating at one-quarter (instead of one-third) of the Kelvin wave speed for the first baroclinic mode, as well as a reduced domain (140°E–90°W) compared to the TD model (125°E–80°W).

To ease the comparison between the two models, we introduce a reduced version of our model that considers only the contribution of the first vertical baroclinic mode (*n* = 1). Such a reduced model (TD-M1 model) serves as a benchmark for the comparison to the Jin model, because the “one-mode” approximation in our model is equivalent to consider shallow-water dynamics in a two-layered ocean. The same atmospheric and mixed layer models from sections 2b and 2c are used in the JN, TD, and TD-M1 models. For clarity, it is convenient to use the same fraction of the total wind stress forcing that is projected onto the baroclinic ocean in all three models. In the case of the TD model, the *P _{n}* coefficients give the fraction of wind stress forcing that is projected onto the corresponding vertical baroclinic modes, where the total fraction is expressed as . Therefore, the first, second, and third baroclinic modes project 19%, 18%, and 7%, respectively, of the total wind stress forcing and 44% in total (see Table 1). In the case of the JN model, which comprehends a single baroclinic mode and considers a mean thermocline depth at 150 m with a 50-m mixed layer depth, the momentum forcing that is projected onto the baroclinic ocean accounts for 33% of the total wind stress forcing. To compare the TD model and the TD-M1 model with the Jin model, the coupling efficiency is chosen so that wind stress forcing projecting onto the baroclinic ocean is adjusted to reach 33% of the total wind stress forcing, as in the JN model. The latter is achieved by adjusting/scaling the ocean–atmosphere coupling coefficient.

A stability analysis of the two models is performed for a set of values of *μ*. The eigenvalues of the ENSO mode (defined here as the dominant oscillating eigenmode of the system) are presented in Fig. A1 for the three models as a function of *μ*. To ease the comparison, eigenvalues are omitted for the range of values of coupling efficiency where the ENSO mode is not the leading eigenmode of the system. In all three models, the frequency (growth rate) of the ENSO mode continuously decreases (increases) with *μ*. Both the JN model and the TD-M1 sustain a low-frequency ENSO mode (with a frequency between 0.1 and 0.5 yr^{−1}) only when values of the coupling efficiency are large (2.5 ≤ *μ* ≤ 3.5). This suggests that the stability of the JN and TD-M1 models are somewhat similar, although the JN model also exhibits solutions at a lower value of the coupling efficiency (*μ* ≤ 2.5). For the intermediary range of value of the coupling efficiency (*μ* ~ 1), the JN model sustains a high-frequency ENSO mode (~1 yr^{−1}) because it has a fast oceanic adjustment time related to the phase speed of the first baroclinic mode (*c*_{1} = 2.7 m s^{−1}). By comparison, the TD model sustains a lower-frequency ENSO mode (~0.4 yr^{−1}) around *μ* ~ 1 because higher baroclinic modes (*c*_{2} = 1.6 m s^{−1} and *c*_{3} = 1.1 m s^{−1}) introduce longer oceanic adjustment time scales. A similar comparison is made in Dewitte (2000) but with an extended version of the Cane–Zebiak model. Nevertheless, the reason why the JN model and the Cane–Zebiak model simulate a proper time scale of ENSO is that the ENSO time scale can also be modified by the other model’s parameters, such as the sensivity of entrainment temperature to the thermocline depth, for example (Dewitte and Perigaud 1996). Intermediate coupled models, such as the Cane–Zebiak model, and simple models, such as the JN model, and the TD model all relate to some point on an empirical formulation that best fits their needs to reproduce realistic ENSO characteristics. Most of the assumptions are usually made on the mixed layer thermodynamics, especially when nonlinearities are taken into account. Here, it is suggested that higher baroclinic mode contribution (*n* = 2, 3) introduces longer oceanic adjustment time scales, thus leading to larger and more realistic ENSO time scales (2–7 yr). This does not exclude the possibility to reduce the ocean dynamics to a single baroclinic mode for simplicity; however, if so, its characteristics should be representative of the whole baroclinic motion.