Observational, Theoretical and Computational Research on the Climate System
|© Springer-Verlag 2005|
Mechanisms of tropical Pacific interannual-to-decadal variability in the ARPEGE/ORCA global coupled model
Carole Cibot1, Eric Maisonnave1, Laurent Terray1 and Boris Dewitte2
|(1)||Sciences de lUnivers au CERFACS, CERFACS/CNRS URA 1875, 31057 Toulouse Cedex 1, France|
|(2)||Laboratoire dEtude en Géophysique et Océanographie Spatiale, IRD - Centre de Nouméa, 101 promenade Roger Laroque- Anse Vata, BP A5, 98848 Nouméa Cedex, New Caledonia, France|
Received: 25 February 2004 Accepted: 25 November 2004 Published online: 4 May 2005
Abstract Spatial and temporal structures of interannual-to-decadal variability in the tropical Pacific Ocean are investigated using results from a global atmosphere–ocean coupled general circulation model. The model produces quite realistic mean state characteristics, despite a sea surface temperature cold bias and a thermocline that is shallower than observations in the western Pacific. The periodicity and spatial patterns of the modelled El Niño Southern Oscillations (ENSO) compare well with those observed over the last 100 years, although the quasi-biennial timescale is dominant. Lag-regression analysis between the mean zonal wind stress and the 20°C isotherm depth suggests that the recently proposed recharge-oscillator paradigm is operating in the model. Decadal thermocline variability is characterized by enhanced variance over the western tropical South Pacific (~7°S). The associated subsurface temperature variability is primarily due to adiabatic displacements of the thermocline as a whole, arising from Ekman pumping anomalies located in the central Pacific, south of the equator. Related wind anomalies appear to be caused by SST anomalies in the eastern equatorial Pacific. This quasi-decadal variability has a timescale between 8 years and 20 years. The relationship between this decadal tropical mode and the low-frequency modulation of ENSO variance is also discussed. Results question the commonly accepted hypothesis that the low-frequency modulation of ENSO is due to decadal changes of the mean state characteristics.
The El Niño Southern Oscillation (ENSO) phenomenon is the dominant mode of climate variability at interannual timescale. Active in the tropical Pacific, it is a coupled ocean–atmosphere mode, which influences nearly all regions of the globe. Although an increasing number of studies have focused on ENSO for the last 20 years, the detailed balance between several mechanisms interplaying in the ENSO is far from being completely understood and the phenomenon remains difficult to forecast. Two factors may explain the difficulties in predicting ENSO: many climate variability models consider the tropical Pacific as a basin isolated from the extra-tropical regions and prescribe a fixed mean state. However, the tropical Pacific is linked to extra-tropical latitudes via ocean–atmosphere interactions, on both interannual and decadal timescales. Furthermore, several studies have found evidence of decadal climate variability modes in the tropical Pacific (Wang 1995; Goddard and Graham 1997). For example, the unusually warm conditions prevailing during the period 1980–2000 over the tropical Pacific were associated with the strongest El Niño ever observed (1982/1983 and 1997/1998). Latif et al. (1997) have suggested that the decadal mode dominating the 1990s is a source of ENSO irregularity and is therefore an important factor in predicting the ENSO. Furthermore, most ENSO prediction models failed to predict the conditions in the tropical Pacific after 1991 (Ji et al. 1996) and do not take into account intrinsic characteristics of ENSO that evolved during the last century. For example, the sea surface temperature (SST) related to ENSO has changed since the 1976–1977 climate regime shift (Trenberth and Hurrell 1994). Before 1977, positive sea surface temperature anomalies (SSTA) were initiated in the eastern equatorial Pacific and were moving westward (Rasmusson and Carpenter 1982). After 1977, the warm events appeared first in the central Pacific and then spread eastward (Wang 1995).
Despite the progress made over the past few decades and the recent successful international programmes like Tropical Ocean Global Atmospheres/Coupled Ocean Atmosphere Response Experiment (TOGA/COARE) programme, the mechanisms responsible for the dynamics of ENSO as a coupled phenomenon are not yet well understood. In particular, owing to the scarcity of observations covering the tropical Pacific, over both space and time, it is quite difficult to study the links between ENSO and decadal variability from observations alone.
Theoretical and observational studies have proposed various mechanisms for the origin of decadal variability in the tropical Pacific. They may be classified into four categories. The first one suggests that ENSO-like decadal variability is due to extra-tropical–tropical interactions, in the ocean and atmosphere, involving North Pacific subtropical cells (STC, Gu and Philander 1997; Kleeman et al. 1999). However, the hypothesis of Gu and Philander (1997) on the subduction of temperature anomalies has been criticized by observational (Schneider et al. 1999) and modelling (Nonaka and Xie 2000) studies, which argue that subducted temperature anomalies do not actually reach the equator. The second hypothesis proposes that Pacific decadal variability is generated primarily within the North Pacific (Latif and Barnett 1994; Barnett et al. 1999a) and may influence the tropics via atmospheric (Barnett et al. 1999b; Pierce et al. 2000) and/or oceanic (McCreary and Lu 1994; Liu et al. 1994; Lysne et al. 1997) teleconnections. The third suggests that tropical decadal variability is only due to intrinsic tropical processes (Knutson and Manabe 1998). The fourth states that tropical low-frequency variability is linked to the atmospheric stochastic forcing (Kirtman and Schopf 1998), with the ocean having a reddening effect.
Among the tropical processes proposed, the study of Luo and Yamagata (2001) has suggested that positive SST anomalies in the eastern equatorial Pacific are linked to a negative wind stress curl in the southwest of the tropical Pacific. This negative wind stress curl causes the thermocline to shoal, leading to negative subsurface temperature anomalies. The temperature anomalies then move in a northwesterly direction to the western central equatorial Pacific, turn eastward and propagate along the subsurface region to reach the eastern Pacific to reverse the process.
Apart from the long list of preceding mechanisms, some recent studies have proposed causal relationships between decadal-scale climate anomalies and low-frequency modulation of interannual ENSO-like variability (Knutson et al. 1997; Timmermann et al. 1999), mainly through nonlinear processes (Timmermann and Jin 2002; Rodgers et al. 2004).
The main purpose of this paper is to investigate the mechanisms responsible for decadal variability in the tropical Pacific, as well as its links with the low-frequency modulation of the ENSO variability. To perform this study, a 200-year integration of a global coupled ocean–atmosphere model as well as available surface observations are considered.
The paper is organized as follows. Section 2 is devoted to a brief description of the coupled model and simulation, the observed data-sets and the statistical techniques used. In Sect. 3, we present a validation of the model in terms of mean seasonal cycle and interannual variability within the tropical band. In Sect. 4 we examine in detail the mechanisms responsible for the tropical decadal variability that is simulated by the model and in Sect. 5 we study its links with the low-frequency modulation of ENSO. A summary and concluding remarks are given in Sect. 6.
2 Model description, observed data-sets and statistical tools
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 m2 s–1, which reduces in the tropics to reach 2,000 m2 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 (CH4, N2O, CFC11-12) as well as CO2, O3 and H2O. 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 Model validation: coupled model climatology
3.1 Mean and seasonal cycle of SST and wind stress
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.
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.
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).
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.
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.
4 Decadal thermocline variability
Thermocline processes are likely to play an important role in tropical variability at decadal timescales. The subduction hypothesis suggests that temperature anomalies subduct at mid-latitudes and propagate in the thermocline towards the tropics at a decadal timescale (Gu and Philander 1997). Baroclinic Rossby waves associated with thermocline variability may also contribute to decadal timescale processes via ocean circulation adjustments (Capotondi and Alexander 2001). This section analyses the decadal thermocline variability in the coupled model integration. We first determine the main centres of thermocline decadal variability and then investigate the mechanisms responsible for these maxima.
We have first used the STD of D20 (annually averaged) as a representative of the thermocline variability (Fig. 8a). In the Pacific, the simulated thermocline exhibits two main centres of variability. A first pattern is equatorially trapped, between 5°S and 5°N, and represents a global thermocline displacement all along the equator, with a maximum around 130°W of more than 20 m. The second pattern is confined to the western South Pacific region, along 8°S between 160°E and 170°W, with slightly smaller amplitude (18 m). We also note weaker variability along 10°N. The central core along the equator is dominant at ENSO timescale, with variability concentrated between 2 years and 5 years (not shown). The western South Pacific lobe is representative of lower frequency variability. Although this simulated thermocline variability is still dominated by the interannual timescale, its wavelet spectrum is more energetic at a lower frequency, with energy around the 10-year timescale (Fig. 8b, c). To remove the variability associated with the simulated ENSO timescale, the model data are low-pass filtered with a 7-year cut-off. Considering the wavelet analysis of the southwest subtropical thermocline variability (Fig. 8c), the 7-year period coincides with a relative minimum of the energy spectrum. Furthermore, this cut-off value has been previously used in several studies (Knutson and Manabe 1998; Capotondi and Alexander 2001; Luo and Yamagata 2001; Montecinos et al. 2003).
The resulting low-frequency subsurface temperature variability (Fig. 9a) shows maximum amplitude in the southwest tropical Pacific as anticipated, with values reaching 0.8°C around 10°S–5°S/170°W–140°W. It represents 40–60% of the total interannual variability in the southwest part of the tropical Pacific region (Fig. 9b). The vertical structure of the simulated low-frequency variability exhibits a marked maximum between 100 m and 150 m, located around the mean thermocline position (Fig. 9c).
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.
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
4.3.2 Diabatic forcing
Such temperature anomalies on mean isopycnal surfaces may also result from diabatic forcing through anomalous surface heat and freshwater fluxes. In that case, temperature and salinity surface anomalies can lead to local variations of the density gradient in the thermocline due to vertical mixing.
The low-frequency variability of heat and fresh water fluxes appears to be very weak in the region of interest (not shown). Significant anomalies of heat and water fluxes would lead to uncompensated temperature changes and be associated with mass exchange between adjacent layers of fluid, yielding variations in layer thickness. To confirm that the fresh water and heat flux do not play a major role, we have calculated the layer thickness between the 24.6 and 25.4 surfaces (Fig. 11c). In the model, the standard deviation of the layer thickness shows weak variability (lower than 2 m) in the southwest tropical Pacific, representing a small fraction of the total =25 depth STD.
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 (WE) is calculated from the wind stress =(x,y,0) and given by:
where is the unit vector in the vertical direction, is the gradient operator, f is the Coriolis parameter, varying with the latitude, and 0 is the mean density of seawater (0=1,025 kg m–3). To avoid unrealistic values near the equator, where f vanishes, the estimation of WE is restricted to the regions outside the 4°S–4°N latitudinal band.
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).
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).
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.
5 Links with the low-frequency fluctuations of ENSO variability
Following a commonly accepted hypothesis, decadal modes (or decadal modulations of the mean state) are expected to modulate the characteristics of ENSO variability on decadal timescale. However, the recent results of Timmermann (2003) and Rodgers et al. (2004), based on CGCM simulations, are inconsistent with this hypothesis. Rodgers et al. (2004) have suggested that, on the contrary, the decadal ENSO amplitude modulations (DEAMs) could influence the tropical decadal variability. They have identified ENSO nonlinearities (spatial asymmetries between El Niño and La Niña events), which can account for the dominant structure of the decadal variability in their model. Timmermann (2003) has shown that DEAMs and the existence of a decadal tropical mode are more likely two sides of one coin (nonlinear tropical dynamics) than generating one another. So the interaction between the decadal mode of the model and low frequency modulations of ENSO is now investigated.
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).
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.
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).
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.
A 200-year simulation of the ARPEGE-ORCA coupled model has been used to investigate the characteristics of interannual-to-decadal variability in the tropical Pacific. Despite a cold bias and a too shallow thermocline in the western Pacific, the simulated SST seasonal cycle is close to the observed one, although not contrasted enough spatially and in time. The mean equatorial subsurface thermal structure is realistic given that the model is not flux adjusted. The thermocline is however too diffuse. The model simulates ENSO events fairly close to those observed and well represented compared to others coupled models. However, the interannual variability is characterized by too strong and regular events, with a dominant biennial signal. A joint D20 and zonal wind stress analysis suggests that the dominant mechanism exhibits properties similar to those characterizing the recharge–discharge paradigm proposed by Jin (1997). This mechanism also seems to operate at the decadal timescale in the model, although further investigations are needed to firmly demonstrate it. Moreover, the model presents realistic low frequency modulations of ENSO with the same decadal timescale as the observations.
Decadal variability is investigated via tropical subsurface temperature (Tsub25) analysis. The main feature is a maximum pattern of variability in the southwestern tropical Pacific, associated with strong wind stress curl anomaly south of the equator in the central Pacific. The wind stress appeared to be induced by SST anomalies in the eastern equatorial Pacific. This result is interesting, considering that the coupled model seems to reproduce the decadal variability characteristics highlighted by Luo and Yamagata (2001) from re-analysis products. Like in their study, the Tsub25 variability is shown to result from an adiabatic process. The propagation of spiciness anomalies from the southern subduction region is, however, not detected in the simulation, which compromises full inter-comparison with the mechanism proposed by Luo and Yamagata (2001). In addition, Rossby-wave propagation from the southwestern region to the equator is too fast to explain the timescale of the decadal variability depicted in our study. It is likely that different or additional mechanisms than the one proposed by Luo and Yamagata (2001) are operating in the simulation, which deserves further investigation. Decadal interactions with the mid-latitudes could be at work in the model and further analysis should determine if this tropical decadal variability interacts with the north Pacific Decadal Oscillation or the Antarctic Oscillation for example. Nonlinear interactions between the seasonal cycle and the ENSO mode within the tropics are also a possible source of the decadal variability, which may help to resolve the complete cycle of the decadal variability described here. The model decadal variability patterns are shown to be highly correlated to the low frequency modulations of ENSO, in such a way that periods of high ENSO activity are characterized by a flattening of the equatorial thermocline, as recently found in other coupled models (Timmermann 2003; Rodgers et al. 2004). This is a counterintuitive result, as this particular mean state characteristic is expected to accompany reduced ENSO activity. It then questions the commonly accepted hypothesis that low-frequency modulations of ENSO are due to decadal changes of the mean state. We have shown that the spatial patterns associated with periods of high ENSO activity are more likely the result of ENSO nonlinearities (asymmetric spatial structure between El Niño and La Niña events). These nonlinearities (also found in observations) can play a greater role in the model, because they are favoured by stronger and more regular ENSO events than observed. So other mechanisms, which operate in observations and may influence the decadal variability, are expected to be less visible in the model. The model study of Rodgers et al. (2004) has suggested that, in contrast to the common hypothesis, the dominant decadal variability could be induced by the low frequency fluctuations of ENSO. The lag-correlation analysis did not allow us to conclude on the advance or the delay between the decadal variability phases and the low-frequency modulation of ENSO, and at that stage great care must be taken in suggesting a phase relationship. It seems more appropriate to reason in terms of timescale interactions rather than in terms of direct causal relationships. As a matter of fact, recent analyses on these model outputs (B. Dewitte et al., to be submitted) suggest that the ENSO mode may interact with a fast-coupled mode (intra-seasonal timescale), mostly during the cold phases and in the western Pacific. This generates additional nonlinearity and can produce aperiodicity in the system, and consequently ENSO modulation. This timescale interaction could thus be a source of the decadal variability that may come into play with the processes described in this study and allow the description a of self-termination mechanism of a hypothetical decadal coupled mode. These aspects are currently under further investigation.
Acknowledgements We would like to thank the two anonymous reviewers for useful comments as well as K. Rodgers, J. P. Boulanger and J. Picaut for helpful discussions. The authors are also grateful to A. Weaver, C. Cassou and C. Périgaud for their comments on the manuscript. This study was supported by the French Programme National dÉtude Du Climat (PNEDC).