cibotetal2005cd

The ocean component is the global configuration ORCA2 of the OPA8.2 ocean general circulation model (OGCM) (Madec et al. 1998), a hydrostatic primitive equation model with a free surface formulation (Roullet and Madec 2000). The horizontal resolution is 2° in longitude and varies from 0.5° at the equator to 2° cos  towards the pole from the tropics, where denotes latitude. The meridional grid spacing is increased near the equator to improve equatorial dynamics. To overcome the singularity at the North Pole, two poles, located on Asia and North America, have replaced the northern point of convergence. There are 31 levels along the vertical direction, with the highest resolution (10 m) in the upper 150 m. The bottom topography and the coastlines are derived from a global atlas (Smith and Sandwell 1997).

The Jackett and McDougall (1995) equation of state is used. It provides density as a function of potential temperature, salinity and pressure. Vertical eddy diffusivity and viscosity coefficients are computed from a 1.5 turbulent closure scheme that allows an explicit formulation of the mixed layer as well as minimum diffusion in the thermocline (Blanke and Delecluse 1993). The nature of horizontal mixing of momentum is of a Laplacian type with an eddy viscosity coefficient of 40,000 m² s^–1, which reduces in the tropics to reach 2,000 m² s^–1 at the equator, and increases in depth in the deep ocean. The lateral mixing of tracers (temperature and salinity) is isopycnal and a parameterization of eddy-induced velocity with spatially varying coefficient is implemented. Zero fluxes of heat and salt and no-slip conditions are applied at solid boundaries.

The model includes a full thermodynamical sea-ice component, the Hibler-type dynamic–thermodynamic Louvain-La-Neuve sea-ice model, developed at the Université Catholique de Louvain (UCL) by Fichefet and Morales Maqueda (1997).

The atmospheric component is the third version of the ARPEGE-Climat Atmospheric General Circulation Model (AGCM) (ARPEGE is the acronym for Action de Recherche Petite Echelle Grande Echelle) developed at Météo-France (Déqué et al. 1994). The standard configuration of the climate version uses a T63 triangular horizontal truncation. Diabatic fluxes and nonlinear terms are calculated on a Gaussian grid of about 2.8° longitude × 2.8° latitude. The vertical is discretized over 31 levels (20 levels in the troposphere) using a progressive vertical hybrid coordinate extending from the ground up to about 34 km (7.35 hPa). The model time-step is 30 min for this resolution and time discretization scheme.

ARPEGE-Climat includes all basic atmospheric physical parameterizations as well as comprehensive land surface processes via the ISBA model (Douville 1998). The convection is represented by a mass flux scheme with detrainment as proposed by Bougeault (1985). The radiative package follows the Fouquart and Morcrettes scheme (Morcrette 1990). The long-wave radiation parameterization includes the effect of trace gases (CH₄, N₂O, CFC_11-12) as well as CO₂, O₃ and H₂O. Parameterizations of boundary layer and turbulence are based on vertical exchange coefficients that are computed as functions of the local Richardson number according to the method adopted by Louis et al. (1982).

The ARPEGE and ORCALIM (ORCA/Louvain Ice Model) models are coupled through the OASIS 2.5 coupler developed at CERFACS (Valcke et al. 2000), which ensures the time synchronization between the GCMs and performs the spatial interpolation from one grid to another. Daily mean SST, sea-ice extent and surface albedo are provided to the atmosphere GCM, and daily mean surface fluxes of heat, momentum and fresh water (including run-off and calving) are provided to the ocean GCM. Global flux conservation is performed when passing the heat and water fluxes. All concentrations of greenhouse gases are fixed at pre-industrial values and no flux corrections are applied. The atmospheric initial state is in January and is taken from an uncoupled integration of ARPEGE. The initial ocean state is taken from Levitus for temperature and salinity, and is zero for the velocity field.

From the initial state, the integration is performed over 70 years to allow the system to adjust and to reach a quasi-equilibrium state in the upper layers of the ocean. We then integrate the coupled model for 200 years. We only consider the last 200-year period for the subsequent analyses, as they appear free of any systematic drift in the Pacific (the coupled simulation is labelled PR1 thereafter).

Various data-sets are used to assess the model performance. For SST maps and long-term SST-based ENSO indexes, the Hadley Centres SST data-set (HadISST, 1871–1999, version 1.1, Rayner et al. 2003) is used, which includes data up to year 2000. For evaluating the model wind stress climatology, monthly wind stress data were taken from the ERS-1 and ERS-2 data-sets, from January 1992 to December 2000, obtained from CERSAT (http://www.ifremer.fr/cersat/fr/data/download/download.htm). For evaluating the model thermocline structure, we use the upper ocean temperature data from the TAO data-set (20°C isotherm depth), from January 1992 to December 2002. These were obtained from the TAO website (http://www.pmel.noaa.gov/tao) as well as the Levitus climatology for temperature and salinity (Antonov et al. 1998; Boyer et al. 1998).

Spatial and temporal characteristics of the simulated leading modes of Pacific variability are examined using empirical orthogonal functions (EOFs), maximum covariance analysis (MCA) using singular value decomposition (SVD) techniques, lead-lag regressions and cross-correlations, for various atmospheric and oceanic fields.

The statistical significance of the SVD covariance and correlation is tested with a Monte Carlo method. The atmospheric time-series are randomly scrambled and the SVD is re-calculated. This operation is repeated a hundred times, and the original SVD results are compared with the hundred scrambled ones, giving the significance. The significance of cross-correlations and regressions is based on a hypothesis testing, using a Fisher transformation. To account for serial correlation, the effective sample-size is estimated following the method proposed by Chelton (1983) (see also Emery and Thomson 2001).

The wavelet transform is used to determine the dominant modes of variability and how they vary in time. The wavelet power spectrum gives a measure of the time-series variance as a function of time and frequency. The Morlet wavelet is used, and the transformation is performed in Fourier space using the method described in Torrence and Compo (1998). To reduce wrap-around effects, each time-series is padded with zeros. A Fourier red noise background is used to derive the significance levels and confidence intervals of the wavelet power spectrum. To investigate decadal variability, the data are yearly averaged using a 7-year fourth-order Butterworth low pass filter (Christoph et al. 1995; Murakami 1979).

3.1 Mean and seasonal cycle of SST and wind stress

Figure 1 shows the annual mean SST and wind stress for both model and observations. Overall, the simulated patterns compare reasonably well with the observed ones; however, close inspection reveals some important differences. The model cold bias in the western equatorial Pacific is associated with a cold bias in the regions of Inter-tropical Convergence Zone (ITCZ) and South Pacific Convergence Zone (SPCZ), 1–2°C owing to stronger off-equatorial branches of the trade winds. Moreover, these intensified trade winds may also explain a reduced extension of the equatorial warm pool in boreal spring (not shown).

Fig. 1 Left column mean SST (°C) for a PR1 simulation and c HadISST. Contours intervals are 1°C. Right column mean wind stress (10^–2 N m^–2) for b PR1 simulation and d ERS 1-2. Contour lines indicate the wind stress module. Contour interval is 2.5×10^–2 N m^–2

---

The warm bias present in the southeastern Pacific basin is mainly owing to the lack of stratocumulus clouds in the model (Terray 1998). This warming, greater than 2°C in the upwelling region off South America, could also be partially due to the underestimation of the meridional wind component (Fig. 1b, d). The model presents a too zonal SPCZ, which extends too far to the east. Easterly winds are reversed in the central Pacific around 10°S, between 170°W and 120°W, leading to weak northeasterly winds and to the eastward extension of the convergence zone.

The coupled simulation presents a seasonal cycle of SST in the equatorial Pacific reasonably close to the observed one (Fig. 2), although less pronounced both in amplitude and spatial extension. The cold tongue appears prematurely and extends too far west. The warming in April is not well reproduced, as the model simulates stronger-than-observed equatorial easterlies at that time. The simulated maximum of SST occurs with a delay of 1 month, and is 1.2°C too weak.

Fig. 2 Seasonal cycle of SST anomalies (°C) averaged between 5°S and 5°N for a PR1 simulation and b HadISST data-set. Contour interval is 0.25°C

---

Previous studies have suggested the importance of a realistic simulation of the equatorial thermocline structure in ENSO modelling. The thermocline plays a critical role in the so-called Bjerknes mechanism, which links the equatorial SST zonal gradient to the strength of the zonal winds, through equatorial upwelling. The strength of the Bjerknes feedback depends on oceanic characteristics such as thermocline depth, thickness and zonal slope. For instance, it has been suggested that a shallow thermocline depth associated with an intense (small thickness) thermocline, is favourable for generating ENSO events with large amplitude.

The simulated annual mean upper ocean temperature along the equator is reasonably realistic, given that the model is not flux-adjusted. However the model simulates a thermocline shallower than the observations in the western and central equatorial Pacific (Fig. 3), leading to cooler-than-observed SSTs in the western Pacific and a weaker-than-observed thermocline zonal slope in the central Pacific. Furthermore, the simulated mean thermocline is too diffuse, a characteristic bias of current ocean models While the observed seasonal cycle of the 20°C isotherm depth (hereafter named D20) is weak, the simulated one exhibits a noticeable cycle (Fig. 4). The model D20 is shallower than the observed one along the longitude band of 170–110°W in spring, leading to colder-than-observed SSTs. These differences are associated with stronger easterlies in the eastern equatorial Pacific at that time. The model eastern shallowing is associated with a deepening of the D20 in the western Pacific a month later, leading to a strong thermocline slope along 170°W in April–June. A similar feature is present also in the observations, but it occurs earlier and with smaller amplitude.

Fig. 3 Mean temperature (°C) averaged between 2°S and 2°N for a PR1 simulation and b Levitus climatology. Contour interval is 2°C

---

Fig. 4 Seasonal cycle of D20 (m) averaged between 2°S and 2°N for a PR1 simulation and b TAO observations. Contour interval is 10 m

---

 

The leading EOFs of monthly mean SST anomalies for both model and observations are shown in Fig. 5. The dominant interannual simulated signal in the tropical Pacific appears to be close to that observed. The spatial representation of ENSO by the model is realistic and compares well with results from other coupled models simulations (Achuta Rao and Sperber 2002). Nevertheless, several differences between the modelled and observed pattern emerge. First, the simulated ENSO variability is larger than the observed one, with a SST anomaly standard deviation (STD) of 1.3°C compared to 0.8°C for the observations (not shown). Second, the modelled variability is more equatorially confined and extends too far westward in the Pacific (Fig. 5a, b). The model also exhibits a very regular ENSO behaviour with a timescale of 2–3 years, instead of 2–7 years for observations (Fig. 5c, d). These characteristics represent typical systematic errors in global coupled models. The high ENSO amplitude simulated in the quasi-biennial band can be related to the mean state of the model, as proposed by Meehl et al. (2001). A shallow mean thermocline could explain the large amplitude of the model ENSO. This could also be due to the biases of the simulated seasonal cycle (upwelling, zonal current and SST), which lead to the enhancement of a fast coupled mode that interacts with the ENSO mode (B. Dewitte et al., to be submitted), and to the bias of the longitudinal extent of the peak in the interannual wind stress variance near the equator (An and Wang 2000).

Fig. 5 Upper panel spatial pattern of the dominant EOF mode of SST for a PR1 simulation and b HadISST. Contour interval is 0.2°C. Explained variance is indicated on the top. Lower panel wavelet power spectrum (°C²) of the associated time series for c PR1 simulation and d HadISST. The contour levels are chosen so that 75, 50, 25 and 5% of the wavelet power is above each level, respectively. The cross-hatched region is the cone of influence, where zero padding has reduced the variance. Black contour is the 5% significance level, using a red noise (autoregressive lag1) background spectrum

---

 

To further explore the relationship between SST and thermocline behaviour during ENSO, an EOF analysis is performed on the D20 (Fig. 6). The first two EOF modes account for 32% and 24% of the total variance, respectively. They vary on the same interannual timescale as the SST leading mode, with a preferred period between 2 years and 4 years. The second D20 mode is first described. It represents the classical thermocline slope mode showing a deepening/shoaling of the thermocline in the equatorial eastern Pacific during El Niño/La Niña events (Fig. 6b). It leads the SST mode by 2–3 months (with a correlation of 0.88, Fig. 6d) in agreement with a Kelvin-type response to westerly winds in the west central Pacific. Conversely, the first D20 mode is characterized by a large shoaling/deepening in the western and central Pacific, spreading between 7°N and 15°S, during El Niño/La Niña events (Fig. 6a), more representative of the forced Rossby wave response to interannual winds. An opposite sign anomaly is confined along the Peruvian coast, east of 110°W. It lags the second D20 mode by 6 months, and the first SST mode by 3–4 months (with a correlation of 0.92, Fig. 6c). As suggested by Jin (1997), this mode represents the slow adjustment process associated with the discharge–recharge of the zonal mean equatorial heat content. Figure 7 shows lag-regression maps between the equatorial zonal mean wind stress anomalies and the D20 field along the equatorial Pacific. Positive lags indicate that the wind stress anomalies are leading the thermocline depth anomalies by the appropriate time lag. Assuming westerly wind stress anomalies (positive sign), the oceanic response near zero-lag shows a deepening of the thermocline in the eastern Pacific and a shallowing in the western Pacific 2–3 months later, leading to the flattening of the thermocline slope around 150°W. These phase relationships suggest that the Kelvin and Rossby waves play a role in the thermocline adjustment process in the eastern and western Pacific, respectively. As time proceeds, the shallowing of the thermocline that was initially located in the western Pacific, propagates to the east, where the thermocline initially deepened. About 10 months after the anomalous wind stress maximum, the thermocline is seen to shallow in all longitudes. The shoaling of the thermocline then leads to a surface cooling trend, as the mean upwelling is more efficient to draw cold sub-thermocline water towards the ocean surface. The SST anomalies will eventually change sign and reverse the ENSO cycle.

Fig. 6 Upper panel spatial pattern of the two first EOF modes of D20: a mode 1, contour interval is 3 m, b mode 2, contour interval is 5 m. Explained variance is indicated on the top. Lower panel lag-correlation between the SST leading PC and the c first D20 PC and the d second D20 PC. Dashed line indicates the 95% significant level

---

Fig. 7 Lag-regression of the D20 (m) on the zonal wind stress averaged between 140°E and 90°W along the equator. Contour interval is 2 m. Thick contour indicates the 99% significant level

---

 

Observations show a very similar behaviour during the 1997–1998 El Niño event (http://www.pmel.noaa. gov/tao/jsdisplay/plots/assortgif/uwnd_sst_iso20_anom_mon_nocap. gif), with a large shoaling of the thermocline along the whole equatorial Pacific during the year following the eastern surface warming. This chain of mechanisms can be related to the recharge–discharge oscillator model for ENSO proposed by Jin (1997). He proposed a mechanism based on the build-up of equatorial heat content, and its discharge phase during El Niño events. During the warm phase, westerly wind anomalies and equatorial eastern Pacific warm SST anomalies are associated with a divergence of the zonally integrated Sverdrup transport, which leads to the discharge of equatorial heat content and to the associated shallowing of the thermocline.

Further analysis of the model interannual variability can be found in Dewitte et al. (to be submitted). In the following, we concentrate on the decadal mode simulated by the model, assuming a priori no timescale interactions.

The spatial and temporal structure of the main decadal variability mode is now investigated for the subsurface temperature data. In order to describe the different mechanisms responsible for the simulated decadal variability, the analysis is focused on the temperature along the mean isopycnal surface associated with the main thermocline variability.

The temperature field has thus been projected onto the three-dimensional density field, and, in particular, upon the isopycnal surface representing the main thermocline. In the region where the simulated decadal thermocline variability is maximum, the mean D20 is located between the =24.5 and =25 isopycnal surfaces (Fig. 9c). The =24.5 isopycnal outcrops in the central equatorial Pacific (around 110°W) several years during the coupled simulation, so we have chosen to use the temperature projected onto the =25 surface. However, repeating the analysis for the temperature projected onto the =24.5 surface shows the same results. The temperature is then projected along the time-averaged isopycnal depth (and denoted by the acronym Tsub25). Note that the temperature anomalies used in the following analysis are thus not necessarily salinity compensated.

The leading decadal mode of tropical temperature variability on the =25 isopycnal is isolated via an EOF analysis (Fig. 10a). This mode accounts for 47% of the total decadal variance. The main centre of action is located in the southwest tropical Pacific, with a maximum of variability between the equator and 10°S and from 160°E to 140°W. The eastern Pacific is dominated by values of opposite sign with maxima located off-equator. The wavelet spectrum (Fig. 10b) shows enhanced power at periods between 8 years and 20 years as well as a very low-frequency modulation of the decadal variance.

Fig. 10 a Spatial pattern of the dominant EOF mode of Tsub25. Contour interval is 0.15°C. Dashed contours represent negative values. b Wavelet power spectrum of the time-series associated (°C²). The contour levels are chosen so that 75, 50, 25 and 5% of the wavelet power is above each level, respectively. The cross-hatched region is the cone of influence, where zero padding has reduced the variance. Black contour is the 5% significance level, using a red noise (autoregressive lag1) background spectrum. Thick dashed lines indicate periods of 7 and 20 years, respectively

---

 

The subsurface temperature anomalies characterising the decadal mode may have various origins. For instance, Schneider (2000) has suggested a mechanism for decadal variability based upon the anomalous advection of compensated temperature and salinity anomalies (spiciness anomalies) along the isopycnals located in the main thermocline. Others have suggested the influence of adiabatic displacements of the thermocline associated with a purely mechanical forcing from the winds, or the role of diabatic forcing from anomalous surface heat and fresh water fluxes.

4.3.1 Spiciness anomalies origin

 

In order to quantify the role of the spiciness anomalies, we compare the STD of the temperature anomalies on the mean, =25 isopycnal (not shown but with a spatial structure very close to the leading EOF mode) to the STD of the temperature along the time-varying, =25 isopycnal (which is a measure of the amplitude of the spiciness anomalies) (Fig. 11a). The decadal variability due to spiciness is negligible in the southwest region (lower than 0.1°C, less than 15% of the Tsub25 variability) in contrast with large amplitudes associated with the subduction regions (located in the eastern subtropical Pacific). The spatial structure of the =25 isopycnal depth (Fig. 11b) indicates that the Tsub25 anomalies in the southwest tropical Pacific are associated with relatively large vertical thermocline displacements. This result suggests that the model temperature variability in the southwest tropical Pacific is primarily explained by variations of the thermocline depth rather than anomalous advection of spiciness anomalies. Note that this result may be specific to the present simulation. Indeed, such a limitation is not restricted to the ORCALIM-ARPEGE model. Most of the coupled model simulations performed so far fall into one preferred mode of variability, which is overestimated compared to the many modes of variability that matter for the observed evolution of climate.

Fig. 11 a Standard deviation of the temperature on the time-varying =25 surface (°C). Contour interval is 0.05°C. b STD of the =25 surface depth (m). Contours interval is 1.5 m. c STD of the 24.6–25.4 layer thickness (m). Contour interval is 1.5 m

---

4.3.3 Ekman forcing

So the low-frequency simulated thermocline variability has been shown to be primarily due to adiabatic processes associated with wind forcing. A MCA between the Tsub25 and Ekman pumping is used to identify clearly the relation between the two fields (Fig. 12). The Ekman pumping (W_E) is calculated from the wind stress =(_x,_y,0) and given by:

(1)

where is the unit vector in the vertical direction, is the gradient operator, f is the Coriolis parameter, varying with the latitude, and ₀ is the mean density of seawater (₀=1,025 kg m^–3). To avoid unrealistic values near the equator, where f vanishes, the estimation of W_E is restricted to the regions outside the 4°S–4°N latitudinal band.

Fig. 12 a Spatial pattern of the leading MCA mode for Tsub25 (shading) (°C) and Ekman pumping (thick contours) (10^–4 cm s^–1). Contour interval for W_ek is 0.25×10^–4 cm s^–1. The explained covariance (scf) and significance level are indicated on the top. b Time-series associated. The dashed line corresponds to Tsub25 and the full line to W_ek. Correlation and significance level are indicated on the top

---

 

The spatial representation of the leading covarying mode of the subsurface temperature (Fig. 12) is very similar to the first Tsub25 EOF mode (Fig. 10a). The associated Ekman pumping variability is mainly located between the Equator and 10°S and centred around 150°W. The pattern is slightly shifted to the east of the subsurface temperature variability. This mode of variability accounts for 47% of the total low-frequency variability of the Ekman pumping, and also resembles its first EOF mode (not shown). A positive Ekman pumping induces surface divergence, resulting, consequently, in a cooling of subsurface temperature. With an explained covariance of 90% and a correlation of 0.95 between the temperature on the mean =25 isopycnal and the Ekman pumping scores, it appears that the two fields are strongly related.

These patterns are consistent with the mechanism proposed by Luo and Yamagata (2001): the low-frequency subsurface temperature anomalies are strongly related to local wind stress anomalies. From their study of atmospheric and oceanic re-analysis data, they have suggested that positive SST anomalies in the eastern Pacific induce negative wind stress curl associated with the atmospheric teleconnection in the south tropical Pacific. This negative wind stress curl causes the thermocline to shoal, leading to negative subsurface temperature anomalies, with a pattern of variability similar to the one presented here. These authors have then proposed that the temperature anomalies could then move in a northwesterly direction relative to the western and/or northward relative to the central equatorial Pacific, turn eastward, move with anomalous easterlies, and reach the eastern Pacific to reverse the process. They found a period of 14 years for the complete cycle. However, our results show a slight difference with their work: they have found a pattern of wind stress curl anomalies located more westward, and oriented in the southeast–northwest direction from 20°S–150°W to 5°S–160°W.

The positive Ekman pumping anomaly originates in the model from a positive surface temperature anomaly in the eastern equatorial Pacific, similarly to Luo and Yamagata (2001). Lag-regression analysis between the SVD score of the Tsub25 field and the (7–25 years) band-pass filtered SST is used to illustrate the decadal evolution of the surface temperature signature (Fig. 13). Near zero-lag (Fig. 13b, c), the resulting spatial structure is a zonal dipole, with positive SST anomalies in the eastern Pacific and weaker negative SST anomalies in the western equatorial and South Pacific. When the SST field leads or lags the decadal Tsub25 mode by 5 years, the spatial SST pattern is similar but with opposite sign. This suggests a timescale of about 10–12 years for the decadal cycle, which is slightly shorter than the period of 14 years found by Luo and Yamagata (2001).

Fig. 13 Lag-regression of SST (°C) on the Tsub25 MCA time-series, for lags a –6 years, b –1 year, c 0 year and d 5 years. Contour interval is 0.05°C. Dashed contours represent negative values. The thick contour indicates the 95% significance level. Negative lags mean SST leads Tsub25

---

 

In order to examine the behaviour of the equatorial heat content (HC) associated with this decadal process, a further lag-regression analysis between the SVD score of the Tsub25 field and the (7–25 years) band-pass filtered HC in the first 300 m is performed. The result indicates that the decadal mode is associated with an east–west seesaw of equatorial HC anomalies (Fig. 14a). At zero-lag, a large pattern of positive HC anomalies develops in the western and central equatorial Pacific (from 140°E to 130°W). East of 130°W, weaker anomalies of the opposite sign lead the western anomalies by 1 year. The positive western anomalies are connected 5 years later to eastern positive anomalies, suggesting an eventual eastward propagation leading to the opposite phase. This result suggests that the ENSO recharge–discharge mechanism proposed by Jin (1997), which operates in the model at the interannual timescale (Fig. 7), also operates at decadal timescale, and could participate to the phase reversal of the decadal mode. Note that like Tsub25, the maximum variability of HC at a decadal timescale is located in the southwestern tropical Pacific around 7°S (Fig. 14b). Furthermore, the HC along the equator lags (by less than 1 year) HC at 7°S and the amplitude of HC along the equator is lower (by ~40%) than the amplitude at 7°S (Fig. 14c).

Fig. 14 Lag-regression of HC (10⁸ W) on the Tsub25 MCA time-series during the period from –10 years to +10 years. Negative lags mean HC leads Tsub25. a HC is averaged between 5°S and 5°N. Contour interval is 0.3×10⁸ W. Shaded region indicates the 95% significance level. b HC (10¹⁵ W) is integrated between 120°E and 80°W. Contour interval is 0.5×10¹⁵ W. c From the previous zonal integration, HC (10¹⁵ W) is averaged in two latitudinal bands: 12°S–4°S (thick line) and 4°S–4°N (dashed line)

---

 

4.3.4 Propagation of the subsurface temperature anomalies

 

Luo and Yamagata (2001) have argued for a propagation of subsurface temperature anomalies from the southwest tropical Pacific to the eastern equatorial Pacific. They have speculated that in the western Pacific, the north-westward and northward motion could be due to Rossby wave propagation and/or due to mean flow advection. In a recent study, White et al. (2003) propose, on the basis of re-analysis data, the existence of a tropical decadal mode involving westward propagation of Rossby waves in the southern hemisphere, along 15°S. They have found a propagation of ~5 years from west of Tahiti to the western boundary near Australia. As a simple test of coupled Rossby wave propagation, we have followed the low-pass filtered anomalies of the =25 isopycnal depth, in the southwest and equatorial regions, from the region of maximum variability to the equator (Fig. 15a). A 100 years of the model are represented. Propagations are quite fast (Fig. 15b), less than 2 years between the southwest region and the equatorial Pacific, which cannot explain the timescale of the decadal mode in the present study. Along the equator, the anomalies do not reach the eastern boundary of equatorial Pacific and become evanescent around 120–130°W. At last, it was not possible to identify, as in Luo and Yamagata (2001), the propagation of subsurface temperature anomalies by means of coupled Rossby waves or by mean flow advection, which suggests that a different mechanism is at work for decadal variability in the present simulation.

Fig. 15 a Position of sections used in the subsequent Hovmuller diagram: y=6°S (from 80°W to 160°E), x=160°E (from 6°S to the equator), and equator (from 160°E to 80°W). b Low-pass filtered anomalies of the =25 isopycnal depth (m), along the lines shown in (a). Contour interval is 3 m; the zero-line is not represented

---

 

First the decadal ENSO modulations are characterized. Annual means of SST anomalies averaged over the Niño3 region (5°S–5°N, 150°W–90°W) exhibit clear low frequency modulations in addition to the bi-annual signal (Fig. 16a). To represent these modulations, the interannual (2–7 years) wavelet variance of the Niño3 SSTA (referred to as N3Var index) is computed (Fig. 16a), following the Torrence and Webster (1998) method. Also used in Timmermann (2003), this method provides a robust representation of the ENSO variance, by the possibility of selecting a preferred timescale (i.e. 2–7 years for ENSO) to calculate the variance. As expected, periods of high N3Var index amplitude correspond to periods with strong El Niño and La Niña events. A wavelet analysis indicates the N3Var index is modulated on a decadal time-scale, in the 10–20 years spectral band (Fig. 16b).

Fig. 16 Niño3 SST anomalies (°C) (thin line) and N3Var index (°C²) (thick line) for a PR1 simulation and b HadISST. The wavelet power spectrum of the N3Var index is represented for c PR1 simulation and d HadISST. The contour levels are chosen so that 75, 50, 25 and 5% of the wavelet power is above each level, respectively. The cross-hatched region is the cone of influence, where zero padding has reduced the variance. Black contour is the 5% significance level, using a red-noise (autoregressive lag1) background spectrum

---

 

This low-frequency modulation of the simulated ENSO is realistic, as the N3var index calculated for SST observations varies on the same decadal timescale. This is especially true for the last 50 years, characterized by a high ENSO activity (Fig. 16d). We note the stronger amplitude of the simulated N3Var index (Fig. 16a, c), due to the intrinsic bias of this coupled model, which overestimates ENSO amplitude.

The hypothesis that the decadal mode influences ENSO variability is then tested. To determine the decadal background state changes associated with periods of high ENSO variance, regression of low-pass filtered =25 isopycnal depth and SST on the N3Var index is made (Fig. 17a, b). We observe a pattern of =25 isopycnal depth variability in the southwestern Pacific, associated with a shoaling of the thermocline, very similar to the decadal mode (dominant EOF mode of the low-pass filtered D20, not shown). SST presents a west–east equatorial gradient, with an eastern warming and a western cooling. Like =25 isopycnal depth, this pattern is closed to the decadal mode, in particular the second low-pass filtered EOF mode (not shown). Another interesting behaviour is the reduction of the thermocline slope along the equator, which is associated with a thermocline deepening and subsurface warming in the eastern equatorial Pacific. According to Zebiak and Cane (1987) and Meehl et al. (2001), these particular mean state characteristics (thermocline flatness, reduced zonal gradient of SST) are expected to accompany reduced ENSO activity. Therefore, the simulated thermocline slope reduction cannot be a cause for the associated high ENSO variance, in opposition with the hypothesis that the decadal background state influences ENSO variability. Two possibilities remain: the decadal mode is caused by low-frequency modulations of ENSO, or both are the result of a third cause. As a nonlinear analysis, composites have been made to differentiate the thermocline and SST behaviour during the periods of high and weak ENSO activity. Results (not shown) indicate that periods of high ENSO activity are associated with the same decadal background characteristics as shown by the regression analysis. However, it is interesting to note that these characteristics are not symmetric: periods of weak ENSO activity do not exhibit significant patterns.

Fig. 17 Regression of a low-pass filtered =25 surface depth and b SST on the N3Var index. Dashed contours represent negative values. The thick contour indicates the 95% significance level. Contour intervals are a 1 m and b 0.05°C. c Lag-correlation between the N3Var index and the SVD score of Tsub25. Full line indicates the 95% significance level, dashed line the 99% level. Negative lags mean the N3Var index leads Tsub25

---

 

To determine their phase relationship, the lag-correlation between the SVD score of Tsub25 and the N3Var index is represented in Fig. 17c. The maximum of correlation (0.65) is significant at the 99% level. It occurs at zero-lag and at 1-year lag (N3var ahead Tsub25). Since low-pass filtered annual means are used, this lag of 1 year is not significant enough to conclude on a causal relation.

Following the hypothesis that low frequency modulations of ENSO generate decadal variability, Rodgers et al. (2004) suggest the decadal mode is due to an asymmetry between the anomaly patterns associated with El Niño and La Niña states, this asymmetry reflecting nonlinearities in ENSO variability. Since the results presented previously are similar to the findings of Rodgers et al. (2004), the existence of an asymmetry in our coupled model is tested from composites of El Niño and La Niña events. As the simulated ENSO events are phase-locked on the seasonal cycle, composites have been calculated with annual means centred on December, to better fit the ENSO variability and characteristics at interannual timescale. The composite for warm events and the composite for cold events are first calculated, then summed. The pattern obtained reflects the spatial asymmetries between the El Niño and La Niña events. This residual is calculated for the =25 isopycnal depth and the SST (Fig. 18).

Fig. 18 Residual patterns of the asymmetry between El Niño and La Niña events for a =25 isopycnal depth and b SST. Contour intervals are a 3 m and b 0.1°C. Dashed contours represent negative values

---

 

The SST residual exhibits a zonal dipole along the equator, with positive anomalies in the east, negative anomalies in the west and a zero crossing near 170°W. It results from stronger anomalies during warm events rather than during cold ones, and a more westward extension of the negative anomalies during cold events. The =25 isopycnal depth pattern shows a flattening of the thermocline along equator, and a shoaling in the southwestern Pacific. Changes in the thermocline slope are a characteristic of warm and cold events, but appear less pronounced during cold events. In contrast, the southwestern anomaly is only present during warm events (not shown). So these residual patterns result from the stronger amplitude of El Niño events compared to La Niña events. But the main characteristic is that they are very similar to the structures associated with periods of high ENSO activity. So the counterintuitive decadal background characteristics associated with periods of high ENSO activity can be interpreted as the residual of the spatial asymmetry between El Niño and La Niña events. Keeping in mind that these spatial structures are also quasi-identical to the decadal mode ones, it suggests, as described in Rodgers et al. (2004), that the dominant decadal mode of the model could be induced by the low frequency fluctuations of ENSO.

These results are consistent with observations. During the last two decades, the equatorial mean state exhibited thermocline flatness and a reduced zonal gradient of SST (Wang and An 2001) whereas the amplitude of ENSO became stronger, and its nonlinearity as well (Jin et al. 2003). The SST asymmetry between El Niño and La Niña events is also present in the HadISST observations, but with a less pronounced negative western anomaly. However, as the model variability is too strong and too regular, the residual pattern signature can act stronger in the model than in observations, and could be sufficient to induce the dominant decadal pattern of the model.