Climate model experiments are analyzed to elucidate if and how the changes in mean climate in response to doubling of atmospheric CO2 (2xCO2) influence ENSO. The processes involved the development, transition, and decay of simulated ENSO events are quantified through a multimodel heat budget analysis. The simulated changes in ENSO amplitude in response to 2xCO2 are directly related to changes in the anomalous ocean heat flux convergence during the development, transition, and decay of ENSO events. The weakening of the Walker circulation and the increased thermal stratification, both robust features of the mean climate response to 2xCO2, play opposing roles in ENSO–mean climate interactions. Weaker upwelling in response to a weaker Walker circulation drives a reduction in thermocline-driven ocean heat flux convergence (i.e., thermocline feedback) and, thus, reduces the ENSO amplitude. Conversely, a stronger zonal subsurface temperature gradient, associated with the increased thermal stratification, drives an increase in zonal-current-induced ocean heat flux convergence (i.e., zonal advection feedback) and, thus, increases the ENSO amplitude. These opposing processes explain the lack of model agreement in whether ENSO is going to weaken or strengthen in response to increasing greenhouse gases, but also why ENSO appears to be relatively insensitive to 2xCO2 in most models.
Increasing greenhouse gas (GHG) experiments coordinated by the Coupled Model Intercomparison Project phase 3 (CMIP3) do not agree whether the El Niño–Southern Oscillation (ENSO) is going to strengthen or weaken. Whether ENSO has changed due to recent observed warming is also controversial according to the observational record (e.g., Trenberth and Hoar 1997; Harrison and Larkin 1997; Rajagopalan et al. 1997). For these reasons, the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4) concluded that there is no consistent indication of discernible changes ENSO amplitude in response to increasing GHGs (Meehl et al. 2007). Given that ENSO is the dominant mode of tropical variability, the lack of agreement among models is an important source of uncertainty for projecting future regional climate change throughout the Pacific basin (IPCC AR4).
In contrast, the CMIP3 models largely agree in the response of the mean ocean climate, that is, the background ocean conditions over which ENSO variability occurs. When forced with increasing GHGs, the great majority of models simulate a shoaled, less tilted, and sharper thermocline; weaker zonal currents; and weaker upwelling (Vecchi and Soden 2007; DiNezio et al. 2009). These robust ocean responses are driven by a weakening of the Walker circulation, for which there is observational evidence (Vecchi et al. 2006), and by increased thermal stratification in the upper ocean. ENSO theory indicates that any of these changes in the mean climate can lead to changes in the strength of the ENSO feedbacks, and thus ENSO amplitude; yet their direct influence on ENSO simulations in CMIP3 climate models is not evident (Vecchi and Wittenberg 2010; Collins et al. 2010).
Theoretical, observational, and modeling studies have linked changes in the thermocline with changes in ENSO amplitude. The linear instability analysis of Fedorov and Philander (2001) showed that a sharper thermocline leads to weaker ENSO amplitude in a simple coupled ocean–atmosphere model. This result contradicted previous results from general circulation model (GCM) experiments of Münnich et al. (1991), which showed increased ENSO variability. The Fedorov and Philander model showed that the increased stratification also leads to changes in the mean climate that render ENSO less unstable. This result has not been confirmed by coupled GCM experiments. In contrast, enhanced ENSO variability in response to increase of GHGs is generally attributed to a sharper thermocline in coupled GCM experiments (e.g., Timmermann et al. 1999; Park et al. 2009).
Conversely, the results of Fedorov and Philander (2001) indicate that a shallower thermocline could lead to enhanced ENSO variability. Observations, in contrast, suggest that the strong ENSO events of the 1980s and 1990s could be a result of a deepening of the thermocline after the 1976 climate shift (Guilderson and Schrag 1998) or a sharper thermocline due to GHG related warming (Zhang et al. 2008). However, the observational evidence is not conclusive because 1) there is evidence of strong ENSO activity before the twentieth century (e.g., Grove 1988) and 2) ENSO has been relatively quiet during the first decade of the twenty-first century despite continued warming. Coupled GCMs exhibit a robust relationship between increased ENSO amplitude and reduced vertical diffusivity (i.e., a sharper thermocline) in the equatorial thermocline (Meehl et al. 2001). This relationship explains why the previous generation of ocean models, which had very diffuse thermoclines, simulated much weaker ENSO variability than observed.
All models participating in CMIP3 simulate a sharper thermocline in response to increasing GHGs, yet not all of them simulate a stronger ENSO. Other physical processes, such as shoaling of the thermocline, weaker upwelling, or warmer sea surface temperature (SST) could also have an amplifying or damping effect on ENSO. Thus, it is reasonable to hypothesize that, depending on the balance of these changes, ENSO could strengthen or weaken (Vecchi and Wittenberg 2010; Collins et al. 2010). A few studies, however, have actually attempted to isolate and quantify the contribution from each feedback (e.g., van Oldenborgh et al. 2005; Philip and van Oldenborgh 2006; Kim and Jin 2011a,b). Philip and van Oldenborgh (2006) used a simplified SST equation to show that shoaling of the thermocline enhances ENSO variability, but the warmer mean SST results in stronger atmospheric damping. Kim and Jin (2011b), used the Bjerknes (BJ) index to show how, depending on the balance among the different ENSO feedbacks, the changes in mean climate are directly related to whether ENSO strengthens or weakens in response to increasing GHGs. However, because the BJ index is computed for the Niño-3 region, the results do not indicate the spatial patterns involved in ENSO–mean climate interaction.
In this paper we also quantify the contribution from the robust changes in the mean climate on ENSO as simulated by CMIP3 models. In section 2 we present the climate model experiments. In section 3 we perform a heat budget analysis of ENSO variability directly from the output of an ensemble of preindustrial and CO2-doubling global warming climate simulations performed with 10 CMIP3 coupled GCMs. The heat budget is computed as a balance between the heat storage rate, the advective heat flux convergence, and the net atmospheric heat flux. In contrast with the studies discussed above, computing the advective terms at every grid point allows us to explore the spatial pattern of the interaction between ENSO anomalies and changes in mean ocean climate. The methodology also allows us to closely balance the heat budget in all models, increasing our confidence in the attribution of the ENSO changes. This is a key advantage over previous methodologies, which do not necessarily satisfy the requirement of a balanced heat budget. Finally, in section 4 we use the changes in the heat budget to show how robust changes in the mean ocean climate drive opposing effects resulting in the wide range of changes that the ENSO amplitude exhibited by the CMIP3 models in response to increasing GHGs. Results are discussed and conclusions are drawn in section 5.
2. Global warming experiments
In this study we analyze both changes in ENSO variability and in the mean climate of the equatorial Pacific as simulated in climate model experiments coordinated by CMIP3. A preindustrial control experiment is used as a baseline climate. An idealized experiment in which atmospheric CO2 is doubled (2xCO2) with respect to preindustrial levels is used to compute the changes in ENSO and the mean climate. For all models, the preindustrial climate experiment was forced with 280 ppm CO2 and 760 ppb CH4. This is the “picntrl” experiment in the CMIP3 database. See Table 1 for a list of climate models used in this study.
The idealized 2xCO2 experiment starts from the picntrl experiment, increasing CO2 concentrations at a rate of 1% yr−1 from 280 ppm until doubling at 560 ppm on year 71. Then the experiment is run 150 additional years with constant 2xCO2 forcing. This is the “1pctto2x” experiment in the CMIP3 database. All ENSO statistics and heat budgets for the 2xCO2 climate are computed using model output from the last 150 years of the 1pctto2x experiment. The models still exhibit warming trends of less than 0.4 K (100 yr)−1 during the last 150 years of this experiment. However, these trends are small compared with the warming of about 2K during the first 71 years when the GHG forcing is largest. The 2xCO2 changes in the mean climate are computed by differencing the annual-mean climatology from the 2xCO2 (1pctto2x) experiment minus the annual-mean climatology from the preindustrial (picntrl) experiment. The 2xCO2 changes in ENSO are computed by differencing the ENSO statistics during the 150 years of quasi-equilibrated 2xCO2 climate (1pctto2x) minus the ENSO statistics during the 500-yr preindustrial (picntrl) climate.
In the next section we analyze the ocean processes involved in the growth of ENSO events in the unperturbed preindustrial climate. We first compute ENSO anomalies with respect to the climatological seasonal cycle. Then, we regress these anomalies on the tendency of the Niño-3 index (∂N3/∂t index) so as to estimate the magnitude and spatial pattern of the physical processes involved in the development phase of ENSO events. More details on this procedure can be found in the appendix.
Throughout this study we focus on those aspects of the ENSO mechanisms and their response to 2xCO2 that appear in the multimodel mean. To provide an indication of how robust these signals are across the different models, we also indicate where models agree with the sign of the multimodel mean anomaly or change (e.g., Fig. 2, nonstippled areas). This estimate of robustness does not provide information about how close the model anomalies/changes are to the multimodel mean and, thus, is not useful to detect outliers. However, it remains useful in our study because much of the debate on the sensitivity of ENSO to increasing GHGs has been on the sign of the amplitude change (i.e., weaker or stronger). In addition, we analyzed the response by each individual model to avoid making erroneous conclusions from the multimodel mean.
3. Robust ENSO mechanisms
a. Recharge mode
All models simulate thermocline anomalies with spatial pattern and time evolution indicating their fundamental role in the development, transition, and decay of ENSO events. In all models, thermocline depth anomalies and sea surface temperature anomalies (SSTA) are in quadrature throughout the ENSO cycle (Fig. 1, red and blue lines respectively). The multimodel composite shows that the thermocline deepens (red line) about 10 months before the maximum warming of the cold tongue (blue line). The thermocline shoals about a year later after the peak of the warm ENSO event, driving the transition into the cold phase of the ENSO cycle.
The multimodel composite heat budget [Eq. (A1)], also shows that the anomalous heat storage rate (= ρ0cpH∂T′/∂t) results almost entirely from the anomalous ocean heat flux convergence (Fig. 1, gray and dashed black lines respectively). In contrast, the net air–sea heat flux damps SSTA throughout the entire ENSO cycle (green line). The anomalous ocean heat flux convergence results from anomalous temperature advection by resolved and parameterized ocean currents along with the effect of subgrid-scale processes, such as mixing and entrainment. Because only monthly-mean fields were archived by CMIP3, can only be computed as a residual between the heat storage rate , and . However, we also estimate the contribution to from anomalous temperature advection by resolved currents (Fig. 1, black line). The close correspondence of and in the multimodel and in each individual composite shows that the advection by resolved currents is a good approximation of the total effect of ocean dynamics on the heat budget on ENSO time scales. Note that does not include mixing or entrainment but it includes the nonlinear terms from monthly-mean fields. More details on how and are computed are given in the appendix.
The multimodel composite also shows in phase with (Fig. 1) indicating that ocean dynamics, and in particular the equatorial thermocline, plays a fundamental role in the generation of ENSO events in all models. Note that deepening of the thermocline prior the development of a SSTA is approximately in phase with . For this reason, throughout our analysis we regress anomalies on the tendency of the Niño-3 index (∂N3/∂t index) so as to capture the magnitude and spatial pattern of the different physical processes driving . More details on these regressions can be found in the appendix.
The spatial pattern of the deepening of the thermocline during the development of ENSO events exhibits a zonal mean deepening along the equatorial waveguide (Fig. 2a). The models also simulate increased sea level consistent with a deeper thermocline (Fig. 2b). Thus, both the phasing among , , and SSTA (Fig. 1, red, black, and gray lines respectively) and the spatial pattern of prior to the development of ENSO events are consistent with the recharge oscillator (Jin 1997) or with the delayed oscillator (Schopf and Suarez 1988; Battisti 1988; Suarez and Schopf 1988; Battisti and Hirst 1989).
The multimodel regressions of the thermocline-driven anomalous surface stratification, anomalous zonal currents, and anomalous upwelling show how ENSO interacts with the mean climate of the equatorial Pacific. The deepening of the thermocline drives anomalously weak stratification, ∂T′/∂z, in the upper 100 m of the ocean over the central Pacific (Fig. 3a, colors) where the climatological upwelling is also strongest (Fig. 3a, contours). This indicates that during the development phase of ENSO events the anomalous ocean heat flux convergence (hereafter ENSO heat flux convergence) results from the vertical advection of thermocline temperature gradient anomalies by climatological upwelling (i.e., ) (Battisti 1988; Battisti and Hirst 1989). The zonal currents during the development phase (estimated from the regressions) exhibit eastward anomalies located in the eastern Pacific (Fig. 3b, colors). In the presence of the climatological zonal SST gradient (Fig. 3b, contours), these anomalies also contribute to the ENSO heat flux convergence (i.e., ).
Wind anomalies are negligible during the recharge or development phase, thus the current anomalies (u′), estimated with the regressions cannot be driven by local winds, which only weaken when the ENSO SSTA is developed. However, the regressions are consistent with the dynamics of the recharge mode, which has associated zonal current anomalies (Kirtman 1997; Clarke 2010), since it is a packet of Kelvin waves reflected from the western boundary as a result of the wind stress curl (WSC)-forced Rossby waves. Geostrophy and the meridional gradients in the thermocline anomalies can also lead to zonal current anomalies, however not on the equator (Jin et al. 2006).
Vertical velocity during developing ENSO events (w′) exhibit anomalous downwelling located in the eastern Pacific (Fig. 3c, colors). This downwelling is not a response to local winds since the trade winds do not weaken until the ENSO SSTAs develop but is consistent with the convergence of the anomalous zonal currents at the eastern boundary. The meridional currents at this stage of the ENSO cycle are only significant on the coast (not shown), suggesting coastally trapped Kelvin waves. These anomalous currents diverge on the equator driving upwelling; thus, the anomalous downwelling suggested by w′ (Fig. 3c) can only be explained by the convergence of u′ due to the recharge mode (Fig. 3b). In the presence of the climatological stratification (Fig. 3c, contours), w′ must also contribute to the ENSO heat flux convergence.
b. Linear ENSO heat budget
The anomalous heat flux convergence associated with anomalous thermocline , zonal currents , and upwelling , are estimated as the advection of temperature anomalies (primed quantities) by climatological fields (bar quantities) as:
We use resolved monthly-mean ocean fields to compute these terms of the linear heat budget because these are the highest resolution ocean data available in the CMIP3 database. The primed quantities are anomalies with respect to the climatological annual cycle.
The multimodel regressions of these fields on the ∂N3/∂t index indicate that, during the development of ENSO events, the anomalous ocean heat flux convergence due to advection of the upper-ocean temperature anomaly by climatological upwelling , is concentrated in a narrow band in the central equatorial Pacific (Fig. 4a). Note that the largest coincides to where the climatological upwelling, , is strongest (Fig. 3a, contours). The anomalous heat flux convergence due to advection of the climatological upper-ocean temperature by anomalous zonal currents , is strongest in the eastern Pacific (Fig. 4b), coincident with the location of u′ (Fig. 3b, colors). The flux is strongest in the eastern Pacific close to the coast of South America (Fig. 3c), coincident with the location of the anomalous downwelling w′ (Fig. 3c, colors). Note that we do not include the effect of meridional currents in the heat flux convergence because these terms tend to cancel each other on ENSO time scales and their magnitude is relatively smaller to the terms of (1).
In all models is strongest over a narrow area in the equatorial waveguide coinciding approximately with the operational Niño-3.4 region. This region is where coupling among SST, wind, and thermocline anomalies is strongest due to the presence of east–west gradients in the climatological SST and thermocline depth (Schopf and Suarez 1988). We define a Niño-3.4m region located in the central equatorial Pacific (2.5°S–2.5°N, 180°–110°W) where is largest and, thus, SSTAs are more likely to drive anomalous winds and close the Bjerknes feedback loop. Note that this Niño-3.4m region is narrower and more westward than the observational definition so as to account for SST biases in the models. Also note that we use a slightly different index, the Niño-3 index, to quantify the amplitude of ENSO events and to capture the spatial pattern of the variables involved in development phase of events. The regressions are not sensitive to the index used because the two indices have tendencies that are highly correlated.
4. ENSO response to global warming
The coupled models analyzed here do not agree in the sign of the changes in ENSO amplitude in response to global warming as reported by previous studies (van Oldenborgh et al. 2005; Philip and van Oldenborgh 2006; Guilyardi 2006; Merryfield 2006). The intermodel differences in the changes in ENSO amplitude are directly linked to the intermodel differences in the 2xCO2 change in ENSO ocean heat convergence, (Fig. 5). Here we compute the change in ENSO amplitude as the difference in standard deviation of the (dimensional) N3 index between the 2xCO2 and preindustrial climates. The models with stronger ENSO amplitudes in the 2xCO2 climate (y axis) exhibit increased , (x axis), and vice-versa. For instance, GFDL CM2.1 simulates an increase in ENSO amplitude of about 0.2 K along with an increase in of about 7 W m−2 and FGOALS-g1.0 simulates a reduction in ENSO amplitude of about 0.5 K commensurate with a reduction in of about 7 W m−2.
The close relationship between the 2xCO2 changes in ENSO amplitude and in ENSO heat flux convergence is not unexpected because, as discussed in section 3, SST anomalies not only result from, but also drive the changes in, via the Bjerknes feedback. Thus a cause-and-effect link cannot be immediately established. Moreover, because the 2xCO2 climate is computed from just 150 years, the 2xCO2 changes could arise from unforced centennial changes in ENSO amplitude. A recent modeling study using the GFDL CM2.1 has suggested that multidecadal and centennial changes in ENSO amplitude are possible, even in the absence of external forcing, (Wittenberg 2009). Thus, given the shortness of the 1% to CO2-doubling (1pctto2x) experiment, the changes in ENSO amplitude computed from the last 150 years may not isolate the response to 2xCO2 forcing.
To determine whether the changes in amplitude are due to 2xCO2 forcing we compare them with estimates of centennial changes in ENSO amplitude from the preindustrial control experiments. The range of possible multidecadal and centennial unforced changes in ENSO amplitude is computed as the standard deviation between the different ENSO amplitudes during overlapping 100-yr periods taken every 50 years from the preindustrial experiments. These estimates of uncertainty are shown in Fig. 5 as vertical error bars. Most of the models exhibit changes in amplitude that are larger than the range of unforced centennial changes and therefore are attributable to 2xCO2. The large uncertainty exhibited by GFDL CM2.1 is consistent with the results of Wittenberg (2009), yet the 2xCO2 change in ENSO amplitude is very likely to be externally forced because it is larger than the unforced 1σ range of ENSO amplitudes.
The spatial patterns of the 2xCO2 changes in ENSO heat flux convergence also correspond with the spatial pattern of the 2xCO2 changes in ENSO amplitude ΔSSTA (Fig. 6). The models that simulate stronger ENSO amplitude in the 2xCO2 climate (GFDL CM2.1, MRI CGM2.3.2a) show a pattern of positive ΔSSTA (Fig. 6a) and positive (Fig. 6c) concentrated in the central Pacific. The models that simulate weaker ENSO in the 2xCO2 climate [Community Climate System Model, version 3.0 (CCSM3.0, FGOALS-g1.0, IPSL CM4), show a pattern of negative ΔSSTA (Fig. 6b) and negative (Fig. 6d) concentrated in the central Pacific. Note that the models with stronger ENSO in the mean climate have ΔSSTA and displaced.
Changes in ENSO amplitude and the associated can result from changes in the branch of the Bjerknes feedback loop involving SST and wind changes, even in the absence of changes in background ocean conditions. This involves the response of the equatorial trade winds to a given SST anomaly and depends mostly on how the Walker circulation responds to latent heat release associated with convective precipitation. The sensitivity of these processes can certainly change as the tropical atmosphere warms up in response to the 2xCO2 forcing. We quantify the strength of the wind–SST coupling by defining a coupling coefficient as the regression coefficient between the monthly anomalies of zonal surface wind stress in the Niño-4 region (5°S–5°N, 140°E–160°W) and SST in the Niño-3.4m region (Guilyardi 2006). A large coupling coefficient indicates a stronger response of the trade winds for the same magnitude of SSTA. Some of the models analyzed here exhibit large changes in coupling coefficient in the 2xCO2 climate. However, these changes are not related to the changes in the (intermodel r = −0.16, Fig. 7) nor ENSO amplitude (intermodel r = −0.08, figure not shown). For instance, IPSL CM4 and FGOALS-g1.0 exhibit increases in coupling of 25% and 9%, respectively, but they fail to translate into increased ENSO amplitude in the 2xCO2 climate. Note that, with the exception of MRI CGM2.3.2a, the majority of the models exhibit increased or unchanged coupling coefficient. The enhanced wind response to a given SSTA could result from increased latent heat release in a warmer climate due to nonlinearity of the Claussius–Clapeyron equation. A cogent explanation for this is lacking in the literature and is beyond the scope of this study.
In contrast, the changes in are related to the changes in and . In general, the models with increased ENSO amplitudes also exhibit an increase of all three terms of the linear heat budget. Note that the intermodel are well captured by the intermodel changes in advective heat flux convergence (Fig. 8a). This allows us to use the linear decomposition of the heat budget to attribute the changes in ENSO amplitude. The models show a close relationship with (Fig. 8b) and (Fig. 8d) averaged over the Niño-3.4m region (intermodel r = 0.82 and r = 0.71 respectively). We compare the changes in heating averaged over Niño-3.4m region because this is where the resulting SST changes are most effective at influencing the atmospheric circulation, closing the ENSO feedback loop. Not all models exhibit downwelling anomalies in the central Pacific (not shown): this is why not all the models show a close relationship with upwelling , (Fig. 8c). Particularly, the models with reduced ENSO amplitude in the 2xCO2 climate do not exhibit changes in (Fig. 8c, models CCSM3.0, FGOALS-g1.0, IPSL CM4).
However, cannot be used to attribute the 2xCO2 changes in ENSO amplitude without entering into a circular argument because of the Bjerknes feedback. For instance, according to (1a), can change through changes in the mean upwelling or changes in the anomalous stratification . However, only the former is directly related to the 2xCO2 changes in mean climate, while is to the change in ENSO amplitude.
The influence of the changes in the mean climate on ENSO becomes clear when the changes in each term of the linear ENSO heat flux convergence (1) are computed:
Throughout this paper the delta notation refers to 2xCO2 climate changes and primed quantities are ENSO anomalies, for example, w′ are the upwelling anomalies with respect to the monthly-mean seasonal cycle, which in the equatorial band are dominated by ENSO variability. Thus, the Δ operator applied to a primed quantity indicates a 2xCO2 change in an ENSO anomaly. Conversely, a delta applied to a bar quantity indicates a change in mean climate.
Equation (2) shows that the changes in cannot be immediately used to attribute changes in ENSO amplitude because the second term in each of the integrals on the right-hand side includes 2xCO2 changes in ENSO anomalies (Δ∂T′/∂z, Δu′, Δw′), thus leading to a circular argument. However, the first term in the integrand of (2) involves the 2xCO2 changes in the mean climate (, , ) and the ENSO anomalies in the control climate (, u′, w′). Thus, these terms can be used to quantify the effect of the changes in mean climate on ENSO heat flux convergence as
This expression can be interpreted as the heat flux convergence that results from the interaction of ENSO in the unperturbed climate (primed quantities) and the changes in the mean climate in response to 2xCO2 (deltas of bar quantities). According to (3) this anomalous heat convergence is due to changes in climatological upwelling , changes in climatological zonal temperature gradient , and changes in climatological stratification . Here we focus on the effect of the changes in the mean ocean climate on ENSO amplitude; however, ENSO amplitude can change due to other processes, such as wind–SST coupling and atmospheric damping. These changes will also lead to a change in via changes in the ENSO anomalies , Δu′, and Δw′ [second term in Eq. (2)].
The changes in ocean heat flux convergence due to the changes in the mean climate—that is, due to changes in the climatological upwelling, zonal temperature gradient, and stratification—are robust among the seven models that have a realistic thermocline feedback (Figs. 9 and 10). The first term in (3), the change ENSO heat convergence due to changes in climatological upwelling, is negative; that is, it acts to reduce and thus weaker ENSO amplitude (Fig. 8a and Fig. 10, blue bars). This response results from weaker climatological upwelling in the 2xCO2 climate (i.e., < 0), driven by weakening of the Walker circulation (Vecchi and Soden 2007; DiNezio et al. 2009). The second term in (3) is positive in the upper thermocline and negative in the lower thermocline (Fig. 8b). The resulting increase in ENSO heat flux convergence in the surface layer (Fig. 10, light blue bars) is not a result of a stronger SST gradient, but of a stronger subsurface zonal temperature gradient (Fig. 11a). Note that this zonal temperature gradient occurs because the time-mean thermocline shoals in the 2xCO2 climate, also explaining the anomalous cooling below the thermocline (Fig. 8b). The third term in (3), is positive, that is, an increase in , due to a sharper thermocline in the 2xCO2 climate (Fig. 8c). However, note that this response is restricted to the eastern boundary where anomalous downwelling occurs during the growth of ENSO events (Fig. 3c).
The models do not agree on the combined effect of the three processes represented by , despite agreeing on the sign of each individual process. However, is directly related to the changes in (r = 0.84, Fig. 12), which is the main contributor to . This relationship is evident in models with large changes in ENSO amplitude, such as CCSM3.0, FGOALS-g1.0, and GFDL CM2.1. The reduction in ENSO amplitude in response to 2xCO2 simulated by CCSM3.0 and FGOALS-g1.0 occurs because the effect of weaker mean equatorial upwelling dominates. All models simulate reduced ENSO heat flux convergence due to weaker mean equatorial upwelling (Fig. 10, dark blue bars); however it only leads to weaker ENSO in those models (CCSM3.0, FGOALS-g1.0, IPSL CM4) where this term dominates. This effect is less pronounced in GFDL CM2.1, thus allowing ENSO to strengthen via the effect of the sharper and shallower thermocline on the zonal advection and upwelling terms (Fig. 10, light blue and green bars). Unlike the majority of the models, the downwelling anomalies simulated by GFDL CM2.1 and GFDL CM2.0 during ENSO events extend into the central Pacific (not shown). For this reason, ENSO is more sensitive to changes in stratification in this model (Fig. 10, green bars).
There are two exceptions to this explanation for the diverging ENSO responses simulated by this ensemble of climate models. The changes in simulated by MRI CGM2.3.2a cannot be explained by . However, it is possible that the stronger ENSO in the 2xCO2 climate, despite the cooling effect of , is driven by the (unrealistic) positive net atmospheric heat flux (not shown). The changes in the mean ocean climate results in stronger in the CNRM-CM3 (Fig. 12, dot 9); however, this fails to translate into stronger ENSO amplitude in the 2xCO2 climate. In this model the changes in and are confined to the eastern boundary where the coupling is ineffective in amplifying the changes.
5. Discussion and conclusions
According to this heat budget analysis of the CMIP3 models, ENSO can either weaken or strengthen via changes in the equatorial Pacific Ocean in response to 2xCO2. The changes in ENSO amplitude in the 2xCO2 climate can be directly attributed to 2xCO2 forcing because they are larger than unforced centennial changes estimated from the control climate. Whether ENSO amplitude increases or decreases depends on a subtle balance between the changes in advection of the upper-ocean temperature anomaly by climatological upwelling versus advection of the climatological upper-ocean temperature by anomalous upwelling and zonal currents. The weakening of the Walker circulation and changes in the thermocline in response to 2xCO2 play opposing roles in this balance. In the 2xCO2 climate, the advection of the upper-ocean temperature anomaly by climatological upwelling decreases as the equatorial climatological upwelling weakens in response to weakening of the Walker circulation/trade winds. In contrast, advection of the climatological upper-ocean temperature by anomalous zonal currents increases as the subsurface zonal temperature gradient strengthens due to a sharper thermocline.
Previous studies also reported diverging ENSO responses but they attributed it to different mechanisms (Philip and van Oldenborgh 2006; Kim and Jin 2011b). Their results show a stronger sensitivity to the changes in stratification and in atmospheric damping, which act to increase and decrease ENSO variability, respectively. In contrast, we find that the intermodel differences in ENSO amplitude are mainly the result of a diverging balance between a weaker thermocline feedback and a stronger zonal advection and upwelling feedback. These studies fitted the model variables into a simplified SST equation (Philip and van Oldenborgh 2006) or to the recharge oscillator (Kim and Jin 2011b). Our heat budget decomposes the changes in the temperature equation directly from the models output, thus preserving the spatial correlation between the changes in the mean climate and the ENSO anomalies. This approach also allows us to quantify the different ENSO mechanisms without making any a priori assumptions about their role in ENSO variability.
The BJ index used by Kim and Jin (2011b) is very well suited to estimate the strength of the feedbacks, but fails to preserve the spatial patterns of the ENSO anomalies, which are shown here to be important in the interaction between ENSO and the background climate change. For instance, their methodology averages the model variables over the Niño-3 region, thus the spatial correlation between background climate and ENSO anomalies may be lost. This could be problematic for the upwelling feedback, which is confined to the eastern boundary in the climate models. Thus, averaging over the entire eastern Pacific may render their methodology sensitive to basinwide changes in stratification. Moreover, these studies find an important role for atmospheric damping, thus weakening the ENSO. However, unlike observations, atmospheric fluxes play a smaller role in ENSO variability simulated by the models in the preindustrial climate (Wittenberg et al. 2006). This model bias may render the models insensitive to the changes in atmospheric damping, which should lead to a weaker ENSO.
Myriad mechanisms can give rise to ENSO variability in models. It is not clear whether the real world ENSO is governed by these same mechanisms or that the balance among them is realistic. Therefore, our conclusions cannot be directly extrapolated to project how the real world ENSO will change in response in to increasing GHGs. The existence of well-known biases in the mean cimate, such as the cold tongue and the double ITCZ biases, can be responsible for altering the balance of processes, and therefore sensitivity to 2xCO2. For instance, an excessively strong zonal SST gradient due to the “cold tongue” bias could make the zonal advection feedback stronger in the models. Moreover, the cold tongue bias also results in peak ENSO SSTA that are located off the eastern boundary where the upwelling anomalies occur. This could make the upwelling feedback less sensitive to 2xCO2 changes in stratification. Thus, the real world ENSO could be more sensitive to a sharpening of the equatorial thermocline, and stronger ENSO events become stronger in response to global warming. Furthermore, it is well known that coupled climate models underestimate the role of atmospheric damping (e.g., Wittenberg et al. 2006; Lloyd et al. 2010). For instance, in the majority of the models the transition from warm to cold events is driven by ocean heat flux convergence with a very small contribution from atmospheric fluxes (see the composite ENSO heat budget for CCSM3.0; green lines in Fig. A2b).
Another common bias of ENSO simulations is the lack of asymmetry between warm and cold ENSO events (An et al. 2005; Zhang et al. 2009). Studies focusing on the nonlinear aspects of ENSO and its rectification effect into the mean climate have suggested that the approach followed here—that is, understanding the 2xCO2 response of ENSO as a result of forced changes in the mean climate—may be inherently limited. This alternative view looks at ENSO events as regulators of the stability of the mean climate—specifically the temperature contrast between the warm-pool SST and the thermocline water down below (Sun and Zhang 2006). This regulatory effect is tied to ENSO asymmetry or more generally to the nonlinearity of ENSO dynamics. The models analyzed here exhibit a wide range of asymmetry. For instance, MRI CGCM2.3.2 and GFDL CM2.1 simulate stronger warm events, CCSM3.0 simulates very symmetric events, and CNRM-CM3 simulates stronger cold events; yet the link of the asymmetry and the 2xCO2 response is not evident. Moreover, all models agree on the forced response of the mean climate to 2xCO2 despite the lack of agreement in ENSO response or ENSO asymmetry. More research is evidently needed to bridge these complementary views of ENSO–mean climate interactions.
We have not considered whether changes in high frequency variability, such as the MJO and westerly wind bursts (WWBs), or nonlinearities can result in ENSO changes. Observations suggest that random weather noise helps sustain an otherwise damped ENSO mode (e.g., Penland and Sardeshmukh 1995; McPhaden and Yu 1999; Thompson and Battisti 2000, 2001; Kessler 2002). We have not considered nonlinear terms in the heat budget, which can act as a positive or negative feedback to ENSO (Münnich et al. 1991; Jin et al. 2003; An and Jin 2004). The sensitivity of these processes to global warming and whether changes in their statistics could lead to changes in ENSO amplitude has not been studied in detail.
The heat budget analysis indicates that the 2xCO2 changes in the mean ocean climate play an important role in the changes in ENSO amplitude. The oceanic dynamical response to weakening of the Walker circulation and the increased thermal stratification associated with the surface-intensified ocean warming play opposing roles in the ENSO response. The weakening of mean equatorial upwelling in response to a weaker Walker circulation/trade winds drives a reduction in ocean heat convergence. A stronger mean zonal (subsurface) temperature gradient associated with the increased stratification drives increased ocean heat convergence.
A very tight relationship has been found between intermodel differences in the ENSO response to 2xCO2 and the meridional shape of the zonal wind anomalies in the control climate (Merryfield 2006). According to our analysis, the ENSO weakens in response to 2xCO2 in those the models where the thermocline feedback dominates over the zonal advection feedback. Moreover, these models also have zonal wind anomalies that are meridionally narrower compared with the models where ENSO strengthens. The narrow wind anomalies lead to stronger WSC anomalies and stronger recharge/discharge explaining why the thermocline feedback dominates over the advection feedback in these models. In contrast, the models with wider wind anomalies have relatively weaker WSC anomalies and thermocline deepening during the recharge phase, thus ENSO is less sensitive to the weaker climatological upwelling. As a result, ENSO strengthens in these models due to the stronger zonal advection feedback. This is the same idea put forth by Neale et al. (2008) to explain why a change in the convection scheme in CCSM3 results in wider wind anomalies shifting ENSO variability from being an oscillation to a series of events.
The roles played by weakening of the Walker circulation and the sharper thermocline presented here can be easily understood by contrast with the effect of these mechanisms on the response of the mean climate. In the mean response, the weaker Walker circulation drives a warming tendency opposed by a cooling tendency due to a sharper thermocline (DiNezio et al. 2009). Since ENSO is a perturbation of the mean climate, opposite roles should be expected from these mechanisms. This is what effectively occurs, with a weaker ENSO driven by a weaker Walker circulation and stronger ENSO due to a sharper and shallower thermocline. Note that the sharper thermocline plays a less central role in the ENSO response because its effect is restricted to the eastern boundary where the w′ is largest. An exception to this is the GFDL CM2.1, which simulates the ENSO with downwelling anomalies in the central Pacific and thus is more sensitive to changes in stratification.
The ENSO heat budget presented here has advantages compared with methodologies used by previous studies. Our methodology allows us to compute the contribution of different oceanic processes to the heat budget directly from model output without making assumptions on the origin of ENSO variability. Moreover, we consider the spatial patterns of the ENSO anomalies and the changes in mean climate when we compute their effect on the heat budget. This feature of our methodology becomes very useful to explore the impact of well-known model biases, which are very likely to influence the sensitivity of the simulated ENSO to global warming.
The thermocline feedback, which according to our results is expected to weaken, is still the basis of how El Niño events grow, regardless whether ENSO is self-sustained or noise-driven. However, changes in the statistics of the stochastic forcing and the details of the interaction between high and low frequency modes needs to be considered in order to fully characterize the sensitivity of ENSO to increasing CO2. Moreover, the CMIP3 climate models simulate a too weak atmospheric damping of ENSO anomalies compared with observations. Therefore, the real world ENSO could also weaken due to enhanced atmospheric damping in a warmer climate.
Despite the very large uncertainty associated with the model projections of ENSO changes, it is clear that the sensitivity of ENSO depends on the balance of weaker upwelling driven by the weakening of the Walker circulation and by changes in thermocline depth and sharpness. These two responses have different sensitivities to global warming because the weakening of the Walker circulation is governed by the response of the hydrological cycle. In contrast, the increase in stratification depends on how the surface warming is diffused into the deep ocean. These results indicate that both ENSO simulation and the sensitivity and patterns of tropical climate change need to be improved in order to have reliable projections of ENSO amplitude for the twenty-first century.
We acknowledge the international modeling groups participating in CMIP3 for providing their data for analysis and the Program for Climate Model Diagnosis and Intercomparison (PCMDI) and the IPCC Data Archive at IPCC Data Archive at Lawrence Livermore National Laboratory (LLNL) for archiving the data. PCMDI and the IPCC Data Archive at LLNL are supported by the Office of Science of the U.S. Department of Energy. This research was carried out in part under the auspices of the Cooperative Institute for Marine and Atmospheric Studies (CIMAS), a Joint Institute of the University of Miami and the National Oceanic and Atmospheric Administration (NOAA), Cooperative Agreement NA10OAR4320143. PDN and AC were funded by DOE Grant DESC0004897 and NOAA Grant NA10OAR4310204 during this study.
ENSO Heat Budget
To reveal the ocean processes that influence the amplitude of ENSO events and its sensitivity to GW we focus on the growing phase of ENSO events. We use the tendency of the N3 index, our ∂N3/∂t index, to study the growth of events. The N3 index is computed for each individual model using SST anomalies computed with respect to a climatological seasonal cycle averaged over a box in the equatorial cold tongue. This box spans the east–central equatorial Pacific between 5°N and 5°S, 180° and 90°W and is shifted westward with respect to the conventional Niño-3 region to account for biases in the coupled models. This same N3 region is used for all models. Before computing the time derivative, the N3 indices are bandpass filtered with cutoff frequencies between 18 months and 8 yr so as to capture interannual variability.
Consider the heat budget for a surface ocean layer with constant depth H,
where ρ0cp = 4.1 × 106 J m−3 K−1 is the ocean density times the specific heat of seawater, ∂T/∂t is the tendency of the vertically averaged temperature, Qnet is the net atmospheric heat flux, and Qocn is the convergence of ocean heat transport. Averaging the heat budget (A.1) over the so-called Niño-3 region and computing anomalies by removing the mean seasonal cycle, we obtain the tendency of the N3 index on the left-hand side. For this reason, in this study, we linearly regress the different variables involved in the heat budget on the normalized ∂N3/∂t index in order to diagnose the ocean and atmospheric processes involved in the growth of ENSO events. The dimensional ∂N3/∂t index is computed as centered differences using monthly mean time series and then normalized by the standard deviation to obtained the normalized ∂N3/∂t index.
The ∂N3/∂t index peaks during the development of warm ENSO events (El Niño), during the transition into cold events (La Niña), and during the decay of cold events. The regression of anomalies on the normalized ∂N3/∂t index contains information of these three instances during the life cycle of ENSO. Thus the regressions assume that the spatial patterns of warm and cold ENSO events are symmetric. Moreover, the asymmetry between warm and cold events results from nonlinear terms in the temperature equation, therefore our methodology only estimates the heating due to linear terms. In other words, this methodology assumes that warm and cold ENSO events result from the same physical processes. Observations exhibit warm events with larger amplitude and propagation characteristics than cold events, thus rendering this assumption inadequate; however, it is reasonable for the simulated ENSO events in most of the CMIP3 models due to the lack of skewness between warm and cold events (van Oldenborgh et al. 2005).
The multimodel regression of the heat storage rate Qt = ρ0cpH∂T/∂t on ∂N3/∂t (Fig. 1b) shows a spatial pattern in close agreement with the spatial pattern of the multimodel regression of SSTA on the N3 index (Fig. A1a). This result illustrates how the developed SSTA pattern (Fig. A1a) results from the time integration of the heat storage rate (Fig. A1b). The depth of integration, H, used to compute the anomalous heat storage rate is 100 m. The next section discusses why this value is adequate to capture the subsurface changes influencing SSTA during ENSO events. The heat budget (A.1) indicates that anomalies in could either result from anomalies net atmospheric heat fluxes, , or anomalous convergence of heat due ocean currents, . The latter can be computed as a residual between and using (A.1). The multimodel regressions of and on the ∂N3/∂t index (Fig. A1c) show close agreement in spatial pattern (spatial correlation: 0.99) and magnitude (Fig. A1b). This result is not unexpected but confirms that the heat storage rate associated with growing ENSO events, and hence the amplitude of the developed events, is entirely due to oceanic processes. In other words, in the models as in the actual tropical Pacific, atmospheric fluxes do not play a role during the growth of ENSO events.
The dominant role of the ocean dynamical processes during an ENSO cycle is clearly seen in the evolution of composites of SSTA, , and averaged over the Niño-3 region (Fig. A2). All models simulate negligible when the tendency of SSTA is largest, thus explains the growth of SSTA entirely. For this reason, leads SST by a quarter of a cycle. Moreover intermodel differences in the magnitude of averaged over the N3 region and scaled by the average duration of the growing phase are consistent with the respective ENSO amplitude as measured by the standard deviation of the dimensional N3 index (Fig. A3a). However, cannot be readily used to attribute changes in ENSO because it is computed as a residual from (A.1).
The ocean heat flux convergence computed using resolved monthly mean ocean currents, Qadv, approximates Qocn very well (Fig. A2, compare solid and dashed black lines). We compute Qadv using monthly mean fields of temperature T, horizontal currents (u, υ), and upwelling w following the methodology of DiNezio et al. (2009):
The spatial pattern of the ocean heat flux convergence during the development of ENSO events computed using (A.2) (Fig. A1b) is strikingly similar to the estimate computed as a residual from (A.1) (Fig. A1d). Note that also captures the magnitude and phasing of throughout the ENSO cycle in all models (Fig. A2).
The advective ENSO heat flux convergence , estimated using (A.2) captures the intermodel differences averaged over the Niño-3 region (Fig. A3b). Moreover, the spatial correlation between the multimodel and is 0.98 with models ranging from 0.92 (CNRM-CM3) to 0.99 [CCCma Coupled General Circulation Model, version 3.1 (CGCM3.1)]. Three models [the Model for Interdisciplinary Research on Climate 3.2 (MIROC3.2), CGCM3.1, and INM Coupled Model, version 3 (INM-CM3)] simulate averaged over the Niño-3 region of less than 20 W m−2, compared with the remaining models where it is larger than 30 W m−2. Moreover, as we show in section 4, this is due to a much weaker thermocline feedback, possibly because the zonal structure of the mean thermocline prevents the interannual anomalies from propagating to the east where the thermocline is shallow and coupling with SST and winds is more effective. The choice of the depth of integration, H, and limitations of using a constant depth layer are discussed next.
The total heat convergence due to monthly mean currents averaged over this Niño-3.4m region is closely related the sum of , , and (Fig. A4a). Note that we compute (y axis) using all three components of the monthly mean velocity field [see Eq. (A.2)], including meridional currents. In contrast the linear (x axis) is the sum of , , and as defined in (1). Moreover, does not necessarily need to balance the heat budget because it does not include the effect of mixing, parameterized eddies, and submonthly resolved currents. In contrast, includes all oceanic processes because it is computed as a residual from the heat storage rate and the atmospheric heat fluxes. The appendix shows how nearly balances the heat budget on ENSO time scales, thus can be used to study the interaction of ENSO and the changes in mean climate due to 2xCO2.
The models also exhibit differences in how the advective terms of the linear heat budget (1) contribute to the development of ENSO events. The advective heat flux convergence, , is dominated by , and model values ranging from 5 to 40 W m−2 (Fig. A4b). Anomalous zonal currents also contribute to (Fig. A4c) with values of ranging from 5 to 20 W m−2. In contrast, is negligible over Niño-3.4m in all models (not shown), with the exception of CGCM3.1 in which dominates with values of 8 W m−2. MIROC3.2, CCCma-CGCM3.1, and INM-CM3 simulate much weaker owing to a much weaker (Fig. A4b). These models simulate much smaller ENSO thermocline depth anomalies that the models with stronger ENSO events; yet, their climatological thermocline is as sharp. In contrast, the models with weak ENSO (in the control climate) exhibit localized steep east–west gradients or “thermocline jumps,” which could suppress the eastward propagation of thermocline anomalies associated with Kelvin waves and, hence, diminish ENSO variability (Spencer et al. 2007).
Sensitivity of the heat budget to the depth of integration
Estimating the ocean heat flux divergence on a constant depth layer (A.2), while being physically consistent, poses limitations to fully describe the influence of some of the oceanic processes in the heat budget of the ocean mixed layer. Using a constant depth layer could fail to capture the changes involving the thermocline because of its east–west tilt. For instance, the anomalous stratification associated with deepening of the thermocline prior to warm ENSO events does not occur on a constant depth surface, but follows the east–west tilt of the climatological thermocline instead.
The depth dependence of these processes can also be analyzed by computing the temperature tendency and advection terms at each three-dimensional grid point. An equatorial section of the temperature tendency (Fig. A5b) and the advection of temperature by zonal and vertical velocity (Fig. A5c) regressed on the ∂N3/∂t index shows anomalous convergence of heat uniformly distributed in the upper 100 m in the central and eastern Pacific. For this reason we use H = 100 to vertically integrate the heat storage rate in (A.1) and the ocean heat divergence due to resolved currents (A.2).
The vertical distribution of the temperature tendencies associated with the thermocline zonal current and downwelling anomalies shows more details of the mechanisms discussed in section 3 (Fig. 4). The temperature tendencies associated with changes in thermocline, zonal currents, and upwelling do not depend strongly on the depth of integration H. The temperature tendencies due to anomalous temperature gradients occur in the upper 100 m (Fig. A6a) and are largest in the central equatorial Pacific where the climatological equatorial upwelling is strongest.
The temperature tendencies due to the anomalous zonal currents are large below the surface (Fig. A6b) where the largest climatological zonal temperature gradients are located. The zonal current anomalies are strongest at the surface between 150° and 90°W (Fig. 3b) where, unlike observations, the zonal SST gradient is weak. This occurs because the equatorial cold tongue extends too far to the west in coupled climate models. However, is large in the subsurface owing to the zonal temperature gradient associated with the east–west tilt of the thermocline. This is a clear example of how biases in the simulation of the mean climate can result in an unrealistic balance among ENSO mechanisms. The temperature tendencies due to the anomalous upwelling are large close to the eastern boundary (Fig. A6c) where the downwelling anomalies and the climatological stratification are large (Fig. 3c).