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 N2. 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 N2 is maximal), thickness, for example, sharpness (being the vertical extension of the N2 profile around the mean thermocline depth), and intensity (being the value of N2 at the mean thermocline depth). Although the mean stratification in the equatorial Pacific exhibits important zonal variations, here a mean buoyancy vertical profile N2(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, N2(z) will be referred to as the mean stratification. From N2(z), a basis of vertical functions Fn with associated speed cn can be derived assuming the separation into horizontal and vertical variables as follows:

 
formula

where (u, υ) are the anomalies of zonal and meridionnal currents, respectively; p is the pressure; and (Xz, Yz) 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 β plane as shown:

 
formula

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:

 
formula

with . The qn0 and qnm are the projection of the nth baroclinic mode on the Kelvin and mth meridional Rossby wave structures, respectively, and the ψm terms are the Hermite functions. A differential equation for the ocean dynamics can then be written for each baroclinic mode n:

 
formula
 
formula

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 cn, while the second equation (m > 0) accounts for the free westward propagation of Rossby waves with phase speed −cn/(2m + 1). Odd (even) Rossby waves have a meridional symetric (asymetric) structure. A damping is introduced as a Rayleigh-type friction with e-folding rate ɛn 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 m. We assume an infinitely extended wall in the north–south direction for the eastern boundary xE and a no mass flow condition for the western boundary xW. Therefore, the theorical values of reflection efficiency in the case of a Kelvin wave (m = 0) and the first meridional-mode Rossby wave (m = 1) can be derived as

 
formula

where we use the conditions at xE and at xW.

A constant mixed layer depth of 50 m is assumed, and for each vertical mode a projection coefficient Pn is introduced. By taking H0 the ocean bottom depth, we have

 
formula

where , τx is the zonal wind stress anomalies, and hmix 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 qnm. 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 qnm. The use of several baroclinic modes introduces slower-propagating waves and therefore lower frequencies in the model. One could alternatively increase 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 cn) 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).

Table 1.

Ocean model parameters.

Ocean model parameters.
Ocean model parameters.

b. Atmospheric component

A statistical relationship between SST anomalies along the equator TE and τx 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 TE and τx of the SODA data during 1958–2007, which leads to the modes ek(x) and fk(x, y) associated to the time series αk and βk, such that

 
formula

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

 
formula

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 ek. 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 (f1 and f2) 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 Ln 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 FK and FR are longitudinal profiles independent of the baroclinic mode order that is considered.

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 TE anomalies, thermocline depth anomalies hE, and zonal surface currents anomalies uE are averaged in the equatorial strip and therefore depend only on longitude and time. For TE, we have

 
formula

where

 
formula

Here ɛT is a thermodynamical damping coefficient of 125 days−1, is the mean upwelling, hmix 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 hE and uE from ocean dynamics (see section 2a). Following Dewitte (2000), hE and uE are estimated as

 
formula

with Fn(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 TSUB (the temperature at hmix) 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(Bh0), where h0 = 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.

Fig. 1.

Zonal profiles of FT (solid line) and FZA (dotted line). FT is adimensionalized by and FZA by , where H = 150 m, L is the basin width (m) ranging from 125°E to 80°W, DT = 7.5°C, and c1 is the speed of the first baroclinic mode (m s−1).

Fig. 1.

Zonal profiles of FT (solid line) and FZA (dotted line). FT is adimensionalized by and FZA by , where H = 150 m, L is the basin width (m) ranging from 125°E to 80°W, DT = 7.5°C, and c1 is the speed of the first baroclinic mode (m s−1).

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 3ac). 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 (qnm, TE) 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 μ, 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).

Fig. 2.

Plot of eigenvalues of the TD model. Axes are frequency (yr−1) and growth rate (yr−1). We represent μ with increasing dot size as well as with color value (right scale), ranging from 0 to 1.5. For a coupling equal to unity, eigenvalues are overcontoured in black thick-lined open circles. Five ocean basin modes from the uncoupled case are identified (marked by OBM1, OBM2, OBM3, or an asterisk) as well as the ENSO mode in the coupled case (see text). The plot is symetric along the y axis, as all complex (i.e., oscillating) eigenvalues of the system come with their corresponding complex conjugate.

Fig. 2.

Plot of eigenvalues of the TD model. Axes are frequency (yr−1) and growth rate (yr−1). We represent μ with increasing dot size as well as with color value (right scale), ranging from 0 to 1.5. For a coupling equal to unity, eigenvalues are overcontoured in black thick-lined open circles. Five ocean basin modes from the uncoupled case are identified (marked by OBM1, OBM2, OBM3, or an asterisk) as well as the ENSO mode in the coupled case (see text). The plot is symetric along the y axis, as all complex (i.e., oscillating) eigenvalues of the system come with their corresponding complex conjugate.

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 cn and cn/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 (L/cn + 3L/cn)−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 μ.

Fig. 3.

Dependence of the eigenvalues of the leading eigenmode (ENSO mode) on the relative μ. Solid curve is frequency (yr−1); dotted curve is growth rate (yr−1). The parameters are as in Fig. 2, but μ ranges from 0.3 to 1.3.

Fig. 3.

Dependence of the eigenvalues of the leading eigenmode (ENSO mode) on the relative μ. Solid curve is frequency (yr−1); dotted curve is growth rate (yr−1). The parameters are as in Fig. 2, but μ ranges from 0.3 to 1.3.

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.

Fig. 4.

Time–longitude plots of (a) SST anomalies (°C), (b) thermocline depth anomalies (m) and (c) zonal current anomalies (m s−1) along the equator for the ENSO mode at μ = 1 (see text). The contour interval is 0.5°C for the SST, 5 m for H, and 0.1 m s−1 for U. Shading is for positive values.

Fig. 4.

Time–longitude plots of (a) SST anomalies (°C), (b) thermocline depth anomalies (m) and (c) zonal current anomalies (m s−1) along the equator for the ENSO mode at μ = 1 (see text). The contour interval is 0.5°C for the SST, 5 m for H, and 0.1 m s−1 for U. Shading is for positive values.

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.

Fig. 5.

Eigenvector amplitude as a function of longitude and μ for the ENSO mode. We consider the anomalies of (a) SST (°C), (b) H (m), and (c) U (m s−1). The contour interval is 0.2°C for the SST, 2 m for H, and 0.05 m s−1 for U.

Fig. 5.

Eigenvector amplitude as a function of longitude and μ for the ENSO mode. We consider the anomalies of (a) SST (°C), (b) H (m), and (c) U (m s−1). The contour interval is 0.2°C for the SST, 2 m for H, and 0.05 m s−1 for U.

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 (P1 ~ P2, 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).

Fig. 6.

(a) As in Fig. 5b, but for the contribution of the baroclinic modes to H anomalies (m). (left to right) Contribution of the first, second, and third baroclinic modes is contoured with a 2-m interval. Shading is for values greater than 6 m. (b) As in Fig. 5c, but for the contribution of the baroclinic modes to U anomalies (cm s−1) contoured with a 4 cm s−1 interval. Shading is for values greater than 12 cm s−1.

Fig. 6.

(a) As in Fig. 5b, but for the contribution of the baroclinic modes to H anomalies (m). (left to right) Contribution of the first, second, and third baroclinic modes is contoured with a 2-m interval. Shading is for values greater than 6 m. (b) As in Fig. 5c, but for the contribution of the baroclinic modes to U anomalies (cm s−1) contoured with a 4 cm s−1 interval. Shading is for values greater than 12 cm s−1.

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 Pn 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 P1(1 − δ), P2(1 + δ/2), and P3(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 P1 + P2 + P3. 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 δ (P2 and P3 favored over P1), 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.

Fig. 7.

Dependence of the eigenvalues of the ENSO mode on δ, where δ ranges from −1.5 to 1. The solid curve is frequency (yr−1); the dotted curve is growth rate (yr−1).

Fig. 7.

Dependence of the eigenvalues of the ENSO mode on δ, where δ ranges from −1.5 to 1. The solid curve is frequency (yr−1); the dotted curve is growth rate (yr−1).

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 N2(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

 
.
formula

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 (cn, Pn, Fn) 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 α can be modified to explore the phase space of the ENSO mode in the model while the other parameters (cn, Pn, Fn) 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 α = −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.

Fig. 8.

Plot of eigenvalues of the TD model with sensitivity to μ and α. Axes are frequency (yr−1) and growth rate (yr−1). We represent μ (0 ≤ μ ≤ 1.5) with increasing dot size and α value (−1 ≤ α ≤ 1) with color (right scale). For a coupling equal to unity, eigenvalues are overcontoured by black thick-lined open circles. OBM2 and OBM3 from the uncoupled case are identified as well as EG2 and EG3 (see text). The plot is symetric along the y axis, as all complex (i.e., oscillating) eigenvalues of the system come with their corresponding complex conjugate.

Fig. 8.

Plot of eigenvalues of the TD model with sensitivity to μ and α. Axes are frequency (yr−1) and growth rate (yr−1). We represent μ (0 ≤ μ ≤ 1.5) with increasing dot size and α value (−1 ≤ α ≤ 1) with color (right scale). For a coupling equal to unity, eigenvalues are overcontoured by black thick-lined open circles. OBM2 and OBM3 from the uncoupled case are identified as well as EG2 and EG3 (see text). The plot is symetric along the y axis, as all complex (i.e., oscillating) eigenvalues of the system come with their corresponding complex conjugate.

Fig. 9.

(a) Frequency (yr−1) of the first leading eigenmode (ENSO mode) as a function of α and μ. Contour interval is 0.1 yr−1. Thick black lines separate areas NO, EG3, and EG2 (see text). NO eigenmodes are not represented. Shading is for eigenmodes in the area EG3. (b) Growth rate (yr−1) of the first leading eigenmode. Contour interval is 0.25 yr−1. Thick black lines separate areas NO, EG3, and EG2. Shading is for positive values.

Fig. 9.

(a) Frequency (yr−1) of the first leading eigenmode (ENSO mode) as a function of α and μ. Contour interval is 0.1 yr−1. Thick black lines separate areas NO, EG3, and EG2 (see text). NO eigenmodes are not represented. Shading is for eigenmodes in the area EG3. (b) Growth rate (yr−1) of the first leading eigenmode. Contour interval is 0.25 yr−1. Thick black lines separate areas NO, EG3, and EG2. Shading is for positive values.

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).

Fig. 10.

(a) Amplitude of WWV anomalies (m) as a function of α and μ, defined here as the zonal average of H anomalies in the equatorial strip (1°N–1°S, 125°E–80°W). Contour interval is 2 m. Thick black lines separate areas NO, EG3, and EG2 as defined in Fig. 9. Shading is for values greater than 6 m. Phase difference of H anomalies between the eastern Pacific (90°W) and the central Pacific (180°) is overplotted in thick gray dashed lines (°). (b) Phase difference (°) between Niño-3 SST anomalies and WWV anomalies as a function of α and μ. Niño-3 SST is defined here as the average of SST anomalies in the region 1°N–1°S, 150°–90°W. A positive phase difference is for WWV anomalies leading the Niño-3 SST anomalies. Shading is for values greater than 50°.

Fig. 10.

(a) Amplitude of WWV anomalies (m) as a function of α and μ, defined here as the zonal average of H anomalies in the equatorial strip (1°N–1°S, 125°E–80°W). Contour interval is 2 m. Thick black lines separate areas NO, EG3, and EG2 as defined in Fig. 9. Shading is for values greater than 6 m. Phase difference of H anomalies between the eastern Pacific (90°W) and the central Pacific (180°) is overplotted in thick gray dashed lines (°). (b) Phase difference (°) between Niño-3 SST anomalies and WWV anomalies as a function of α and μ. Niño-3 SST is defined here as the average of SST anomalies in the region 1°N–1°S, 150°–90°W. A positive phase difference is for WWV anomalies leading the Niño-3 SST anomalies. Shading is for values greater than 50°.

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 (WWVn for n = 1, 2, 3) were estimated and analyzed (not shown), where the total WWV anomalies result from WWV = WWV1 + WWV2 + WWV3. The results indicate that first, WWV1 leads WWV2, consistent with the faster oceanic adjustment time scale associated with M1 compared to M2. For the same reason, WWV2 leads WWW3. 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 WWV2.

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 Pi is the wind projection coefficient for the baroclinic mode i and hmix 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:

 
formula

where ρ0 is the mean density, eECK 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.

Fig. 11.

(a) Plot of eigenvalues for the ENSO mode with sensitivity to n considered in the model, namely, n = 1 (purple), n = 2 (blue), n = 3 (green), n = 4 (yellow), and n = 5 (red). Axes are frequency (yr−1) and growth rate (yr−1). ENSO mode is defined as the dominant oscillating eigenmode of the system. Here μ increases from 0 to 2 with an interval of 0.1. For μ = 0.8, eigenvalues are overcontoured in black thick-lined open circles. (b) As in (a), but adding the contribution of the Rossby-3 and Rossby-5 equatorial waves (m = 3 and m = 5, see text) to the model.

Fig. 11.

(a) Plot of eigenvalues for the ENSO mode with sensitivity to n considered in the model, namely, n = 1 (purple), n = 2 (blue), n = 3 (green), n = 4 (yellow), and n = 5 (red). Axes are frequency (yr−1) and growth rate (yr−1). ENSO mode is defined as the dominant oscillating eigenmode of the system. Here μ increases from 0 to 2 with an interval of 0.1. For μ = 0.8, eigenvalues are overcontoured in black thick-lined open circles. (b) As in (a), but adding the contribution of the Rossby-3 and Rossby-5 equatorial waves (m = 3 and m = 5, see text) to the 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 Pn 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 (c1 = 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 (c2 = 1.6 m s−1 and c3 = 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.

Fig. A1.

Range of dependence of the eigenvalues of the ENSO mode to μ. Solid curve is frequency (yr−1); dotted curve is growth rate (yr−1). Results are presented for the TD model (0.3 ≤ μ ≤ 1.8), the JN model (0.3 ≤ μ ≤ 3.5), and the TD-M1 model (2.3 ≤ μ ≤ 4.3). Eigenvalues are omitted for the coupling range where the ENSO mode is not oscillating. TD model and TD-M1 models are artificially modified to project the same amount of wind stress forcing on baroclinic motion as for the JN model (see text).

Fig. A1.

Range of dependence of the eigenvalues of the ENSO mode to μ. Solid curve is frequency (yr−1); dotted curve is growth rate (yr−1). Results are presented for the TD model (0.3 ≤ μ ≤ 1.8), the JN model (0.3 ≤ μ ≤ 3.5), and the TD-M1 model (2.3 ≤ μ ≤ 4.3). Eigenvalues are omitted for the coupling range where the ENSO mode is not oscillating. TD model and TD-M1 models are artificially modified to project the same amount of wind stress forcing on baroclinic motion as for the JN model (see text).

REFERENCES

REFERENCES
An
,
S.-I.
,
2005
:
Relative roles of the equatorial upper ocean zonal current and thermocline in determining the timescale of the tropical climate system
.
Theor. Appl. Climatol.
,
81
,
121
132
.
An
,
S.-I.
, and
F.-F.
Jin
,
2000
:
An eigen analysis of the interdecadal changes in the structure and frequency of ENSO mode
.
Geophys. Res. Lett.
,
27
,
2573
2576
.
An
,
S.-I.
, and
B.
Wang
,
2000
:
Interdecadal change of the structure of the ENSO mode and its impact on the ENSO frequency
.
J. Climate
,
13
,
2044
2055
.
An
,
S.-I.
, and
F.-F.
Jin
,
2001
:
Collective role of thermocline and zonal advective feedbacks in the ENSO mode
.
J. Climate
,
14
,
3421
3432
.
An
,
S.-I.
,
J.-S.
Kug
,
Y.-G.
Ham
, and
I.-S.
Kang
,
2008
:
Successive modulation of ENSO to the future greenhouse warming
.
J. Climate
,
21
,
3
21
.
Battisti
,
D. S.
, and
A. C.
Hirst
,
1989
:
Interannual variability in a tropical atmosphere–ocean model: Influence of the basic state, ocean geometry and nonlinearity
.
J. Atmos. Sci.
,
46
,
1687
1712
.
Belmadani
,
A.
,
B.
Dewitte
, and
S.-I.
An
,
2010
:
ENSO feedbacks and associated time scales of variability in a multimodel ensemble
.
J. Climate
,
23
,
3181
3204
.
Bjerknes
,
J.
,
1969
:
Atmospheric teleconnections from the equatorial Pacific
.
Mon. Wea. Rev.
,
97
,
163
172
.
Blumenthal
,
B.
, and
M. A.
Cane
,
1989
:
Accounting for parameter uncertainties in model verification: An illustration with tropical sea surface temperature
.
J. Phys. Oceanogr.
,
19
,
815
830
.
Boulanger
,
J.-P.
, and
C.
Menkes
,
1995
:
Propagation and reflection of long equatorial waves in the Pacific Ocean during the 1992–1993 El Niño
.
J. Geophys. Res.
,
100
(
C12
),
25 041
25 059
.
Bretherton
,
C. S.
,
C.
Smith
, and
M.
Wallace
,
1992
:
An intercomparison of methods for finding coupled patterns in climate data
.
J. Climate
,
5
,
541
560
.
Burgers
,
G.
,
F.-F.
Jin
, and
G. J.
Van Oldenborgh
,
2005
:
The simplest ENSO recharge oscillator
.
Geophys. Res. Lett.
,
32
,
L13706
,
doi:10.1029/2005GL022951
.
Busalacchi
,
A. J.
, and
M. A.
Cane
,
1988
:
The effect of varying stratification on low-frequency equatorial motions
.
J. Phys. Oceanogr.
,
18
,
801
812
.
Cane
,
M. A.
, and
E. S.
Sarachik
,
1976
:
Forced baroclinic ocean motions. 1. Linear equatorial unbounded case
.
J. Mar. Res.
,
34
,
629
665
.
Carton
,
J. A.
, and
B. S.
Giese
,
2008
:
A reanalysis of ocean climate using Simple Ocean Data Assimilation (SODA)
.
Mon. Wea. Rev.
,
136
,
2999
3017
.
Clarke
,
A. J.
,
2008
:
An Introduction to the Dynamics of El Nino & the Southern Oscillation
.
Academic Press, 324 pp
.
Clarke
,
A. J.
,
2010
:
Analytical theory for the quasi-steady and low-frequency equatorial ocean response to wind forcing: The “tilt” and “warm water volume” modes
.
J. Phys. Oceanogr.
,
40
,
121
137
.
Collins
,
M.
, and
Coauthors
,
2010
:
The impact of global warming on the tropical Pacific Ocean and El Niño
.
Nat. Geosci.
,
3
,
391
397
.
Dewitte
,
B.
,
2000
:
Sensitivity of an intermediate ocean–atmosphere coupled model of the tropical Pacific to its oceanic vertical structure
.
J. Climate
,
13
,
2363
2388
.
Dewitte
,
B.
, and
C.
Perigaud
,
1996
:
El Niño–La Niña events simulated with Cane and Zebiak’s model and observed with satellite or in situ data. Part II: Model forced with observations
.
J. Climate
,
9
,
1188
1207
.
Dewitte
,
B.
,
G.
Reverdin
, and
C.
Maes
,
1999
:
Vertical structure of an OCGM simulation of the equatorial Pacific Ocean in 1985–94
.
J. Phys. Oceanogr.
,
29
,
1542
1570
.
Dewitte
,
B.
,
S.
Illig
,
L.
Parent
,
Y.
DuPenhoat
,
L.
Gourdeau
, and
J.
Verron
,
2003
:
Tropical Pacific baroclinic mode contribution and associated long waves for the 1994–1999 period from an assimilation experiment with altimetric data
.
J. Geophys. Res.
,
108
,
3121
,
doi:10.1029/2002JC001362
.
Dewitte
,
B.
,
S.-W.
Yeh
,
B.-K.
Moon
,
C.
Cibot
, and
L.
Terray
,
2007
:
Rectification of ENSO variability in interdecadal changes in the equatorial background mean state in a CGCM simulation
.
J. Climate
,
20
,
2002
2021
.
Dewitte
,
B.
,
S.
Thual
,
S.-W.
Yeh
,
S.-I.
An
,
B.-K.
Moon
, and
B. S.
Giese
,
2009
:
Low-frequency variability of temperature in the vicinity of the equatorial Pacific thermocline in SODA: Role of equatorial wave dynamics and ENSO asymmetry
.
J. Climate
,
22
,
5783
5795
.
Fedorov
,
A. V.
, and
S. G.
Philander
,
2001
:
A stability analysis of tropical ocean–atmosphere interactions: Bridging measurements and theory for El Niño
.
J. Climate
,
14
,
3086
3101
.
Guilyardi
,
E.
,
2006
:
El Niño-mean state-seasonal cycle interactions in a multi-model ensemble
.
Climate Dyn.
,
26
,
329
348
.
Jin
,
F.-F.
,
1997a
:
An equatorial ocean recharge paradigm for ENSO. Part I: Conceptual model
.
J. Atmos. Sci.
,
54
,
811
829
.
Jin
,
F.-F.
,
1997b
:
An equatorial ocean recharge paradigm for ENSO. Part II: A stripped-down coupled model
.
J. Atmos. Sci.
,
54
,
830
847
.
Jin
,
F.-F.
, and
S.-I.
An
,
1999
:
Thermocline and zonal advective feedbacks within the equatorial ocean recharge oscillator model for ENSO
.
Geophys. Res. Lett.
,
26
,
2989
2992
.
Jin
,
F.-F.
,
S.-I.
An
,
A.
Timmermann
, and
J.
Zhao
,
2003
:
Strong El Niño events and nonlinear dynamical heating
.
Geophys. Res. Lett.
,
30
,
1120
,
doi:10.1029/2002GL016356
.
Kang
,
I.-S.
, and
S.-I.
An
,
1998
:
Kelvin and Rossby wave contributions to the SST oscillation of ENSO
.
J. Climate
,
11
,
2461
2469
.
Lee
,
T.
, and
M. J.
McPhaden
,
2010
:
Increasing intensity of El Niño in the central-equatorial Pacific
.
Geophys. Res. Lett.
,
37
,
L14603
,
doi:10.1029/2010GL044007
.
McPhaden
,
M. J.
, and
Coauthors
,
1998
:
The Tropical Ocean-Global Atmosphere observing system: A decade of progress
.
J. Geophys. Res.
,
103
,
14 169
14 240
.
Mechoso
,
C. R.
,
J. D.
Neelin
, and
J.-Y.
Yu
,
2003
:
Testing simple models of ENSO
.
J. Atmos. Sci.
,
60
,
305
318
.
Meehl
,
G. A.
,
C.
Covey
,
T.
Delworth
,
M.
Latif
,
B.
McAvaney
,
J. F. B.
Mitchell
,
R. J.
Stouffer
, and
K. E.
Taylor
,
2007
:
The WCRP CMIP3 multimodel dataset: A new era in climate change research
.
Bull. Amer. Meteor. Soc.
,
88
,
1383
1394
.
Meinen
,
C. S.
, and
M. J.
McPhaden
,
2000
:
Observations of warm water volume changes in the equatorial Pacific and their relationship to El Niño and La Niña
.
J. Climate
,
13
,
3551
3559
.
Moon
,
B.-K.
,
S.-W.
Yeh
,
B.
Dewitte
,
J.-G.
Jhun
,
I.-S.
Kang
, and
B. P.
Kirtman
,
2004
:
Vertical structure variability in the equatorial Pacific before and after the Pacific climate shift of the 1970s
.
Geophys. Res. Lett.
,
31
,
L03203
,
doi:10.1029/2003GL018829
.
Neelin
,
J. D.
,
M.
Latif
, and
F.-F.
Jin
,
1994
:
Dynamics of coupled ocean-atmosphere models: The tropical problem
.
Annu. Rev. Fluid Mech.
,
26
,
617
659
.
Neelin
,
J. D.
,
D. S.
Battisti
,
A. C.
Hirst
,
F.-F.
Jin
,
Y.
Wakata
,
T.
Yamagata
, and
S.
Zebiak
,
1998
:
ENSO theory
.
J. Geophys. Res.
,
103
(
C7
),
14 261
14 290
.
Rodgers
,
K. B.
,
P.
Friederichs
, and
M.
Latif
,
2004
:
Tropical Pacific decadal variability and its relation to decadal modulations of ENSO
.
J. Climate
,
17
,
3761
3774
.
Savitzky
,
A.
, and
M. J. E.
Golay
,
1964
:
Smoothing and differentiation of data by simplified least squares procedures
.
Anal. Chem.
,
36
,
1627
1639
.
Schopf
,
P. S.
, and
M. J.
Suarez
,
1988
:
Vacillations in a coupled ocean–atmosphere model
.
J. Atmos. Sci.
,
45
,
549
566
.
Timmerman
,
A.
,
J.
Oberhuber
,
A.
Bacher
,
M.
Esch
,
M.
Latif
, and
E.
Roeckner
,
1999
:
Increased El Niño frequency in a climate model forced by future greenhouse warming
.
Nature
,
398
,
694
696
.
Timmerman
,
A.
,
F.-F.
Jin
, and
J.
Abshagen
,
2003
:
A nonlinear theory for El Niño bursting
.
J. Atmos. Sci.
,
60
,
152
165
.
Trenberth
,
K. E.
, and
D. P.
Stepaniak
,
2001
:
Indices of El Niño evolution
.
J. Climate
,
14
,
1697
1701
.
Van Oldenborgh
,
G. J.
,
S. Y.
Philip
, and
M.
Collins
,
2005
:
El Niño in a changing climate: A multi-model study
.
Ocean Sci.
,
1
,
81
95
.
Wang
,
B.
, and
S.-I.
An
,
2001
:
Why the properties of El Niño changed during the late 1970s
.
Geophys. Res. Lett.
,
28
,
3709
3712
.
Yeh
,
S.-W.
,
B.
Dewitte
,
J.-G.
Jhun
, and
I.-S.
Kang
,
2001
:
The characteristic oscillation induced by coupled processes between oceanic vertical modes and atmospheric modes in the tropical Pacific
.
Geophys. Res. Lett.
,
28
,
2847
2850
.
Yeh
,
S.-W.
,
J.-S.
Kug
,
B.
Dewitte
,
M.-H.
Kwon
,
B. P.
Kirtman
, and
F.-F.
Jin
,
2009
:
El Niño in a changing climate
.
Nature
,
461
,
511
514
.
Yeh
,
S.-W.
,
B.
Dewitte
,
B.
Young Yim
, and
Y.
Noh
,
2010
:
Role of the upper ocean structure in the response of ENSO-like SST variability to global warming
.
Climate Dyn.
,
35
,
355
369
,
doi:10.1007/s00382-010-0849-4
.
Zebiak
,
E.
, and
M. A.
Cane
,
1987
:
A model El Niño–Southern Oscillation
.
Mon. Wea. Rev.
,
115
,
2262
2278
.