Previous studies using coupled general circulation models (GCMs) suggest that the atmosphere model plays a dominant role in the modeled El Niño–Southern Oscillation (ENSO), and that intermodel differences in the thermodynamical damping of sea surface temperatures (SSTs) are a dominant contributor to the ENSO amplitude diversity. This study presents a detailed analysis of the shortwave flux feedback (αSW) in 12 Coupled Model Intercomparison Project phase 3 (CMIP3) simulations, motivated by findings that αSW is the primary contributor to model thermodynamical damping errors.
A “feedback decomposition method,” developed to elucidate the αSW biases, shows that all models underestimate the dynamical atmospheric response to SSTs in the eastern equatorial Pacific, leading to underestimated αSW values. Biases in the cloud response to dynamics and the shortwave interception by clouds also contribute to errors in αSW. Changes in the αSW feedback between the coupled and corresponding atmosphere-only simulations are related to changes in the mean dynamics.
A large nonlinearity is found in the observed and modeled SW flux feedback, hidden when linearly calculating αSW. In the observations, two physical mechanisms are proposed to explain this nonlinearity: 1) a weaker subsidence response to cold SST anomalies than the ascent response to warm SST anomalies and 2) a nonlinear high-level cloud cover response to SST. The shortwave flux feedback nonlinearity tends to be underestimated by the models, linked to an underestimated nonlinearity in the dynamical response to SST. The process-based methodology presented in this study may help to correct model ENSO atmospheric biases, ultimately leading to an improved simulation of ENSO in GCMs.
The past decade has seen steady progress in the simulation of the tropical Pacific and its interannual variability in coupled general circulation models (GCMs) (Delecluse et al. 1998; Latif et al. 2001; Davey et al. 2001; Guilyardi et al. 2009b). Advances in model formulation, such as improved parameterizations and increased horizontal and vertical resolutions, have given rise to an improved simulation of the El Niño–Southern Oscillation (ENSO) (AchutaRao and Sperber 2006). This ability to represent interannual variations in the tropical Pacific is an impressive achievement, especially given that the majority of present-day models do not use flux correction (a technique used in earlier GCMs to correct model biases) in simulating the complex ocean–atmosphere interactions involved in ENSO variability.
However, despite this progress, present-day GCMs still show deficiencies in simulating ENSO, associated with model errors in the mean climate and annual cycle of the tropical Pacific (van Oldenborgh et al. 2005; AchutaRao and Sperber 2006; Guilyardi 2006; Guilyardi et al. 2009b). Typical ENSO problems include errors in the ENSO amplitude and frequency, spatial structures that extend too far to the west (Leloup et al. 2008), and biases in the seasonal phase locking. Indeed, the ENSO properties exhibited by these models are too diverse to allow any clear projection of ENSO evolution in a warmer climate (Meehl et al. 2007; Collins et al. 2010). Furthermore, the complexity of the tropical ocean–atmosphere interactions makes it difficult to elucidate the precise origin(s) of modeled ENSO errors.
Nevertheless, many recent studies suggest that the atmosphere model plays a dominant role in determining the ENSO properties in GCMs. Guilyardi et al. (2004) used different combinations of atmosphere and ocean GCMs to show that the atmosphere model is instrumental in determining both El Niño periodicity and amplitude. Toniazzo et al. (2008), using a suite of integrations in which the atmospheric parameters of a coupled GCM are simultaneously perturbed, showed that the resulting ensemble has a range of ENSO characteristics comparable to that seen in multimodel ensembles (e.g., van Oldenborgh et al. 2005; Guilyardi 2006). Other studies have demonstrated that altering the model’s convection scheme, either by adjusting individual processes (Wittenberg et al. 2003; Wu et al. 2007; Kim et al. 2008; Neale et al. 2008) or by implementing an entirely new scheme (Guilyardi et al. 2009a) can have large impacts on the modeled ENSO properties.
Current theory describes the atmospheric processes involved in ENSO using two relatively simple linear feedbacks (Zebiak and Cane 1987; Battisti and Hirst 1989; Jin et al. 2006; Philip and van Oldenborgh 2006). First, the dynamical feedback (μ), also known as the Bjerknes feedback (Bjerknes 1969; Lin 2007), represents the remote zonal wind stress response to a given central-eastern equatorial Pacific SST anomaly; it is a positive feedback that maintains an east–west asymmetry across the equatorial Pacific. Second, a central-eastern equatorial Pacific SST anomaly will also trigger a thermodynamical response, that is, a change in the ocean–atmosphere heat flux budget. In this region, changes in the net heat flux act to dampen SST anomalies, giving rise to a negative feedback, α (Zebiak and Cane 1987; Jin et al. 2006).
Lloyd et al. [2009 (Part I)] took a step toward understanding the role of the atmosphere in the modeled ENSO by calculating these two atmospheric feedbacks in 12 of the coupled GCMs provided by the Coupled Model Intercomparison Project phase 3 (CMIP3) multimodel dataset. It was found that both μ and α are too weak compared to reanalysis datasets, indicating an error compensation in most present-day coupled GCMs. The net heat flux feedback was also shown to be inversely related to the ENSO amplitude in the models, suggesting that α is a prime candidate for explaining ENSO errors in these GCMs. This result is further confirmed by Kim and Jin (2011), who applied the Bjerknes (BJ) stability index [a linear stability analysis based on the recharge oscillator framework (Jin et al. 2006)] to show that intermodel differences in the thermodynamical damping (α) are one of the dominant contributors to the ENSO amplitude diversity in the CMIP3 models.
To shed light on the α biases, the net α feedback was split into its four individual components: the shortwave (SW), longwave (LW), latent heat (LH), and sensible heat flux feedbacks (Lloyd et al. 2009). The SW and LH flux feedbacks were found to dominate in the central-eastern equatorial Pacific region. However, whereas all models successfully simulate a negative LH flux feedback, the majority of models fail to capture the overall negative SW flux feedback observed in this region. The SW flux feedback, αSW, is thus the main source of model uncertainty in α, with 11 out of the 12 GCMs exhibiting a larger error in αSW than in any of the other heat flux feedback components. This result is supported by other studies that investigate the SW flux feedback in present-day GCMs (Sun et al. 2006, 2009; Guilyardi et al. 2009a).
Lloyd et al. (2011) followed up on Lloyd et al. (2009) by analyzing the ENSO atmospheric feedbacks in the corresponding Atmospheric Model Intercomparison Project (AMIP) simulations of the CMIP3 models, a method that can help to pin down the source of feedback biases (e.g., Sun et al. 2006; Guilyardi et al. 2009a). It was found that both μ and α exhibit improvements over their coupled counterparts and that the error compensation between the feedbacks present in the coupled GCMs is no longer evident. Although αSW is also improved in the AMIP simulations, this feedback continues to be the main source of model uncertainty in the overall α. Intermodel differences in the AMIP αSW values are dominated by biases in the eastern equatorial Pacific cloud properties, and their response to changes in the large-scale circulation during ENSO events. This important role for cloud biases is supported by Bony and Dufresne (2005), who show that the tropical cloud response to SST in present-day coupled GCMs disagrees the most with observations in regimes of large-scale subsidence, such as in this eastern equatorial Pacific region.
An open question is whether the atmospheric model biases that account for the AMIP αSW errors also account for the large coupled αSW errors revealed by Lloyd et al. (2009). Analysis of αSW is complicated by the fact that two SW flux feedback regimes coexist in the central-eastern equatorial Pacific region. First, an increased SST in a region of subsidence reduces the static stability of the atmospheric boundary layer. This acts to break up the dominant marine stratiform clouds, increase the surface SW flux, and provide a positive feedback on the SST (Klein and Hartmann 1993; Philander et al. 1996; Park and Leovy 2004; Xie 2005; Lloyd et al. 2011). Second, a SST increase in a region of ascent tends to increase the amount of convective cloud cover, reduce the surface SW flux, and hence provide a negative feedback on the SST (Ramanathan and Collins 1991; Bony et al. 1997). In this study, a “feedback decomposition method” based on these mechanisms is introduced as a means of elucidating the complex αSW model biases.
Until now, αSW has been defined as a linear feedback. This is an extension of the fact that the net heat flux feedback, α, is treated as a linear response in simple ENSO models and frameworks (e.g., Zebiak and Cane 1987; Battisti and Hirst 1989; Gordon and Corry 1991; Jin et al. 2006; Guilyardi et al. 2009a; Lloyd et al. 2009; Kim and Jin 2011; Lloyd et al. 2011). However, Zebiak and Cane (1987) admit that this linear heat flux parameterization is “clearly oversimplified and is probably incorrect in some local regions,” and Barnett et al. (1991) similarly remark that the linear approximation masks important processes such as cloud effects. There is therefore no reason to assume that the SW flux feedback is linear. Here, the SW flux feedback nonlinearity is considered for the first time.
This study of the SW flux feedback over the central-eastern equatorial Pacific is thus motivated by the following questions.
Do the atmospheric model biases that account for the AMIP αSW errors (i.e., biases in the mean cloud properties and their response to dynamical changes) also apply to the coupled αSW errors?
Does the SW flux feedback exhibit a nonlinearity that is masked by the linear αSW definition?
What explains the coupled versus AMIP differences in the αSW feedback values?
In section 2, the CMIP3 simulations analyzed in this study are presented, followed by a description of the observational and reanalysis datasets against which the models are assessed. Section 3 recaps the methodology used to calculate the SW flux feedback, and presents the coupled and AMIP αSW values. In section 4, analysis of the coupled αSW feedback begins with a diagnosis of the SW flux behavior during modeled El Niños and the introduction of a “feedback decomposition method.” In section 5, the SW flux feedback nonlinearity is diagnosed, followed by an investigation into the coupled versus AMIP αSW differences (section 6) and an analysis of the cloud properties over the eastern equatorial Pacific (section 7). Finally, the paper concludes with a summary and discussion of the results (section 8).
We use the preindustrial simulations of 12 CMIP3 GCMs, as previously analyzed in Lloyd et al. (2009) and listed in Table 1. These are long coupled simulations (at least 100 years) in which radiative forcings are held fixed at preindustrial values, allowing the models’ ENSO properties to be studied in an idealized climate scenario with no anthropogenic radiative forcing changes since the Industrial Revolution. The ENSO feedbacks are expected to change in a future warmer climate (Philip and van Oldenborgh 2006; Collins et al. 2010; Kim and Jin 2011), but the purpose of this study is to identify and explain feedback biases in the models’ intrinsic ENSO cycles. The 12 models in this CMIP3 subset represent the full model diversity in ENSO amplitude (Guilyardi et al. 2009b) and supply all the variables required to diagnose the atmospheric feedbacks. We use the first 150 years for 10 of these simulations (not including the spinup period), the full 100 years supplied for MIHR, and years 71–200 for HadGEM1, due to a data problem at the beginning of this simulation.
Nine out of these 12 models also have AMIP simulations [previously analyzed in Lloyd et al. (2011)] in which the atmosphere models are forced by observed SSTs. Although the length of the AMIP simulations varies between models, only the 19 years of data common to all models are analyzed, that is, 1980–98, which includes the two large El Niños of 1982–83 and 1997–98.
b. Reanalyses and observations
When performing intermodel comparisons, it is important to have “benchmark” datasets against which the models can be evaluated. To this end, we use the following three atmospheric reanalyses.
ERA-40: supplied by the ECMWF, covering 1958–2001 (Uppala et al. 2005). The two input SST datasets are the Met Office Hadley Centre Sea Ice and Sea Surface Temperature version 1 (HadISST1) (up to November 1981) and the NOAA/National Centers for Environmental Prediction (NOAA/NCEP) weekly SST analysis (Reynolds et al. 2002), which is used thereafter.
NCEP2: an update of the NCEP–National Center for Atmospheric Research (NCAR) reanalysis, covering 1979–2009 (Kanamitsu et al. 2002). The input SST dataset is the NOAA/NCEP weekly SST analysis (Reynolds et al. 2002).
Objectively Analyzed Air-Sea Fluxes (OAFlux): a heat flux dataset supplied by the Woods Hole Oceanographic Institution, covering 1984–2004. To calculate the sensible and latent heat fluxes, Yu and Weller (2007) combine SST, near-surface wind speed, near-surface air temperature, and near-surface specific humidity data from a number of sources, including satellite retrievals and the ERA-40 and NCEP reanalyses. This “optimal blending” of data leads to improved global estimates of the turbulent fluxes (Yu and Weller 2007). Shortwave and longwave radiative fluxes are supplied by the International Satellite Cloud Climatology Project (ISCCP), and the input SST dataset is the NOAA optimum interpolation (OI) 0.25° daily SST analysis (Reynolds et al. 2007).
In addition to the reanalyses, two observational ISCCP products are used to evaluate the models:
ISCCP D2: Provides global total, high-, mid-, and low-level cloud cover amounts for 1984–2000 (Rossow et al. 1996).
ISCCP Flux Data (FD-TOA): provides global full- and clear-sky top-of-the-atmosphere (TOA) LW and SW fluxes (1984–2000), constructed using ISCCP cloud datasets and the NASA Goddard Institute for Space Studies (GISS) radiative transfer model (Zhang et al. 2004).
It is important to remember that these ISCCP datasets have associated uncertainties due to 1) instrumental biases, 2) errors introduced when converting the raw measurements into useful data, and 3) algorithms used to fill any missing data points. Similarly, atmospheric reanalyses are far from perfect owing to biases in both numerical models and assimilated data. Nonetheless, these are the best references we have with which to compare the models.
3. Diagnosing the SW flux feedback
The net heat flux feedback (α) is defined as
where F′ is the net heat flux anomaly into the ocean and T′ is the SST anomaly, both averaged over the Niño-3 region (5°N–5°S, 150°–90°W) where the SST interannual variability is largest.
The subject of this study, αSW, is the SW flux (FSW) component of this feedback, represented by
We calculate this feedback by linearly regressing the FSW anomalies against the SST anomalies at each grid point in the tropical Pacific (a pointwise calculation) and then averaging the regression values over Niño-3 (Lloyd et al. 2009; Guilyardi et al. 2009a; Lloyd et al. 2011). This method is represented by , where regress(⋯) denotes a linear regression and the overbar denotes a spatial average over the indicated region.
This calculation method is considered more appropriate than directly regressing box-averaged quantities because it allows us to plot maps of the linear regression values and hence visualize the αSW spatial pattern (e.g., Fig. 1). The choice of calculation method does not alter our conclusions. To account for any nonlinearities arising from the annual cycle of the SW flux feedback, the αSW value is calculated as the average of the four seasonal components [December–February (DJF), March–May (MAM), June–August (JJA), and September–November (SON)].
Figure 1 shows the SW flux feedback maps for ERA-40, OAFlux, and the 12 CMIP3 coupled simulations. Blanked out points correspond to correlations of less than 0.1. ERA-40 and OAFlux are characterized by a negative SW flux feedback across the equatorial Pacific region, ranging from over −20 W m−2 C−1 in the west to weaker values as low as −5 W m−2C−1 in the east (Figs. 1a and 1b).
We note that ERA-40 and OAFlux differ in their estimates of this feedback in the eastern equatorial Pacific: the OAFlux feedback is around half the strength of the ERA-40 feedback. Mapping the ERA-40–OAFlux climatological SW flux difference (not shown) reveals that ERA-40 overestimates the OAFlux values by as much as 50 W m−2 in the far eastern tropical Pacific, off the west coast of South America. Cronin et al. (2006) attribute the excessive ERA-40 SW flux in this stratus cloud region to an overly strong SW surface cloud forcing compared to in situ buoy measurements. The OAFlux (ISCCP) SW flux field (and αSW value) is henceforth considered to be more accurate than that of ERA-40.
The spatial maps of the modeled SW flux feedbacks (Figs. 1c–n) reveal large biases in the tropical Pacific region. While most models simulate a negative SW flux feedback in the western equatorial Pacific, the feedback in the eastern equatorial Pacific is not well reproduced. In most models this is due to a region of positive feedback values extending over the eastern Pacific cold tongue, with values of over 10 W m−2 C−1 in HadCM3, HadGEM1, CNRM-CM3, ECHAM5/MPI-OM, and MRI CGCM2.3.2.
Table 1 and Fig. 2 (upper panel) present the αSW values calculated by averaging these pointwise regressions over the Niño-3 region. The values for the corresponding AMIP simulations are also shown. When averaging the regression values, all points are included regardless of their correlation so as to avoid calculating the feedback value from a small proportion of points. As a measure of the uncertainty in the feedback values Table 1 also presents the mean correlations in Niño-3 (bracketed values).
As revealed by Lloyd et al. (2009) and Fig. 1, no coupled model successfully simulates the strong negative αSW feedback found in the reanalyses, and 7 out of the 12 models actually have a positive αSW value (Table 1, column 3). Figure 2 (lower panel) presents the coupled and AMIP αSW biases (with respect to OAFlux) for the CMIP3 models. Every AMIP simulation has an improved αSW feedback value compared to the corresponding coupled simulation, as shown by Lloyd et al. (2011). We investigate these coupled-AMIP αSW differences in section 6, but first seek to understand the large coupled αSW biases (blue bars in Fig. 2, lower panel).
4. The shortwave flux response during model El Niños
In Lloyd et al. (2011), the 1997–98 El Niño was used as a case study to help understand the AMIP αSW biases. However, investigating a single El Niño event in the coupled models is made difficult by the fact that 1) no model event can be directly compared to an observed event and 2) the models themselves exhibit diverse events (both intra- and intermodel differences). This is emphasized by Fig. 3, which shows the Niño-3 SST anomalies during individual El Niños in the reanalyses and coupled simulations. Each of the faint colored lines is a single event, defined as SSTA ≥ 1.5 × SSTSD (SSTSD is the standard deviation of the Niño-3 SST anomaly time series). This criterion must be satisfied for at least three consecutive months for the anomaly to be classified as an event. Furthermore, indices that correspond to April or May anomalies are ignored so as to avoid picking out events that unrealistically peak in boreal spring. The bold red line is the event composite.
The corresponding SW flux anomalies are presented in Fig. 4. The composite ERA-40 and OAFlux events exhibit negative SW flux anomalies during El Niño, in agreement with the large negative αSW values in Table 1. The models tend to exhibit smaller SW flux anomalies than the reanalyses: only GFDL CM2.1, CNRM-CM3, and one event in ECHAM5/MPI-OM have anomalies as large as the maximum OAFlux anomaly (−50 W m−2). Furthermore, half of the models do not show consistency between the Niño-3 SW flux anomalies during modeled El Niños (Fig. 4) and the overall αSW values (Table 1, column 3). The SW flux feedback for these models must therefore be affected by the SW flux evolution during La Niña and neutral conditions. This nonlinearity is studied further in section 5.
Breaking down the El Niño shortwave flux response
To shed light on these diverse SW flux responses and better separate local and remote effects, we introduce a feedback decomposition method. This method, based on the SW flux feedback mechanism, breaks down the local SW flux response (dSW/dSST) into the following three steps.
The dynamical response to SSTs, dω500/dSST: The large-scale circulation, represented by ω500 [the vertical velocity relative to pressure levels at 500 hPa (Bony et al. 1997)] responds to changes in the SST (via changes in the atmospheric stability). For example, positive SSTAs during the 1997/98 El Niño were associated with anomalous ascent (Figs. 5 and 6b in Lloyd et al. 2009).
The cloud response to dynamics, dTCC/dω500: The total cloud cover (TCC) responds to changes in the large-scale circulation. For example, increased ascent during the 1997/98 El Niño gave rise to positive TCC anomalies over most of the equatorial Pacific (not shown).
The SW flux response to clouds, dSW/dTCC: The downward SW flux at the ocean surface responds to changes in the total cloud cover. For example, increased TCC during the 1997–98 El Niño gave rise to a decreased surface SW flux (Figs. 5 and 6c in Lloyd et al. 2009).
Using the chain rule, the product of these three responses is the SW flux response to SST anomalies:
Because each of these individual responses is assumed to be completely local, this simple framework is somewhat idealized. In particular, considering dω500/dSST as a local response ignores the fact that the eastern equatorial Pacific dynamics can also be affected by remote SST changes. Nevertheless, this method is a good starting point for understanding the source of the complex coupled αSW biases.
Table 2 presents the values of the three responses, calculated by performing linear anomaly regressions in Niño-3 for the entire time series. Multiplying together the three responses gives dSW/dSST values that are well correlated with the αSW values in Table 1 (linear correlation coefficient of 0.72 for ERA-40 and the 12 models, significant at the 0.01 level).
The response nonlinearity is considered by calculating the regression values for SSTA > 0 and SSTA < 0 separately, henceforth denoted by plus and minus superscripts respectively. In this section, the El Niño mechanisms are discussed with respect to the SSTA > 0 regression values (using the overall values does not alter the results: the correlation coefficients between the SSTA > 0 and overall regression values are 0.89, 0.97 and 0.97 for dω500/dSST, dTCC/dω500 and dSW/dTCC, respectively). The SSTA < 0 values are discussed in relation to the αSW nonlinearity in section 5.
The values in ERA-40 and NCEP2 are −15.0 and −18.4 hPa day−1 C−1, respectively, indicating increased ascent in response to warm SST anomalies. All coupled models successfully simulate a negative ω500 response to positive SST anomalies but underestimate the ERA-40 and NCEP2 values. This underestimation of the dynamical response is further demonstrated by plotting the ω500 anomalies during modeled El Niños (Fig. 5).
The other two responses, dTCC+/dω500 and dSW+/dTCC, are both negative in ERA-40 and ISCCP (Table 2, columns 3 and 4), indicating increased total cloud cover in response to increased ascent, which acts to reduce the surface SW flux. Many models exhibit biases in these responses: only six (two) models have dTTC+/dω500 (dSW+/dTCC) values that lie between the “observed” values, and nine models underestimate dSW+/dTCC with respect to ISCCP.
In Table 3, explanations are proposed for the modeled SW flux behavior during El Niños (Fig. 4), taking into account the biases in the three responses (Table 2). The final column indicates whether the model El Niño SW flux behavior agrees with its αSW value (as mentioned above). When quantifying the model biases, , dTTC+/dω500, and dSW+/dTCC (referred to as the “dynamical response,” “cloud response,” and “SW flux response to clouds,” respectively) are compared to ERA-40, ERA-40/ISCCP, and ISCCP, respectively. The model biases in each of the three responses (with respect to the observations) are plotted in Fig. 6.
The analysis presented in Table 3 and Fig. 6 suggests that an overly weak ascent response to positive SST anomalies plays a role in the SW flux feedback biases in 10 models. Furthermore, errors in the cloud cover response to the increased ascent may contribute to the SW flux feedback biases in GFDL CM2.1, HadGEM1, CNRM-CM3, and MRI CGCM2.3.2. For instance, the positive cloud response in MRI CGCM2.3.2 (i.e., reduced cloud cover in response to increased ascent) explains why this model has a positive SW flux feedback, whereas the overly strong negative cloud response in CNRM-CM3 may compensate for this model’s weak dynamical response.
Attribution of the SW flux feedback biases is complicated by the model errors in dSW+/dTCC (green circles in Fig. 6). An overly weak dSW+/dTCC value (i.e., underestimated reduction in the surface SW flux for increased total cloud cover) may contribute to the SW flux feedback biases in eight models (Table 3). Physically, model differences in the SW flux response to cloud cover are likely to be related to cloud properties such as optical depth and cloud height. For example, if the increased cloud cover during an El Niño is too optically thick, the dSW+/dTCC value will be too strong. The impact of cloud properties on dSW+/dTCC could be investigated using sensitivity experiments with different cloud parameterizations. Although such experiments are outside the scope of this study, we analyze the eastern equatorial Pacific cloud properties in section 7.
These results indicate that the model SW flux feedback biases (Fig. 4) stem from errors in the dynamical, cloud, and SW flux responses during El Niños. To analyze the relative importance of these response biases, a quantitative measure of the SW flux response to positive SST anomalies is required. This provides motivation for studying the αSW nonlinearity. Furthermore, investigating this nonlinearity will help shed light on the discrepancy between the El Niño SW flux behaviors and the overall αSW values (Table 3, column 3).
5. Nonlinearity in αSW
To diagnose the αSW nonlinearity, we introduce two new variables: and , defined as the Niño-3 linear anomaly regressions of SW flux against SST (dSW/dSST) for SSTA > 0 and SSTA < 0, respectively. The and values are plotted in Fig. 7 (red and blue bars, respectively), and a simple measure of the SW flux response nonlinearity is calculated as − (green bars in Fig. 7). This nonlinearity is strikingly large in the reanalyses and many models due to differences between negative values and positive values.
The values correspond to the El Niño SW flux responses plotted in Fig. 4. Nine out of 12 models simulate an underestimated negative , whereas the other three models (HadGEM1, MIMR, and MRI CGCM2.3.2) simulate positive feedbacks. The SW flux feedback during La Niña situations () is positive in OAFlux and all models except MIHR and CCCMA (blue bars in Fig. 7). This indicates a reduced surface SW flux in response to negative SST anomalies, enabling the La Niña to grow. OAFlux and nine models thus have and feedbacks of opposite sign.
The models tend to underestimate the OAFlux − difference of 13.2 W m−2 C−1: only GFDL CM2.1, CNRM-CM3, and ECHAM5/MPI-OM have nonlinearities within 25% of OAFlux. The strong nonlinearities in GFDL CM2.1 and CNRM-CM3 explain why these models have weak αSW values despite their strong negative SW flux response during El Niños (Fig. 4): the positive feedback during cold situations offsets the negative El Niño feedback. Similarly, the positive values in HadCM3 and IPSL CM4 can explain why these models have an overall positive αSW feedback (Table 1) despite their negative El Niño feedbacks (Fig. 4).
To understand the model differences in the SW flux feedback nonlinearity, we return to the feedback decomposition results, presented in Table 2. Plotting the nonlinearities in dω500/dSST, dTCC/dω500, and dSW/dTCC (calculated as the difference between the SSTA < 0 and SSTA > 0 regression values) against the αSW nonlinearities (not shown) gives linear correlation coefficients for the models of 0.65, −0.33, and −0.19, respectively, of which only the first is significant at the 0.05 level. This suggests that the nonlinear response of ω500 to SST is the dominant contributor to the αSW nonlinearity in the models. Those models with the largest dω500/dSST nonlinearity (e.g., GFDL CM2.1 and CNRM-CM3) also tend to have a larger αSW nonlinearity, and vice versa.
In the reanalyses and models, the dω500/dSST nonlinearities are due to a weak subsidence response for SSTA < 0 in comparison to the stronger ascent response for SSTA > 0. This dynamical nonlinearity is underestimated by all models compared to the reanalyses, mainly due to the underestimated dynamical response for SSTA > 0 (as seen in Fig. 5). The relationship between the dω500/dSST and αSW nonlinearities suggests that an improved model representation of the dynamical nonlinearity would be key to improving the simulated SW flux feedback nonlinearity.
Although the dSW/dTCC nonlinearity is small (less than 0.1 W m−2) in most models and exhibits no link with the αSW nonlinearity, there is a large dSW/dTCC nonlinearity of 0.28 (0.82) W m−2 in ISCCP (ERA-40) (Table 2, column 4). More specifically, the ISCCP (ERA-40) SW flux is 24% (92%) more sensitive to changes in cloud cover for SSTA > 0 than for SSTA < 0. This could be explained by the nonlinearity in the observed high-level cloud cover response to SST (Fig. 8). For SSTA > 0, the mean Niño-3 ISCCP high-level cloud cover increases (maximum positive anomalies of over 20%), efficiently reflecting incoming SW flux [high clouds tend to have a large optical thickness (Ramanathan and Collins 1991)] and giving rise to a strong dSW+/dTCC value. On the other hand, only small decreases in the high-level cloud cover are observed for SSTA < 0 (maximum negative anomalies of around −6%). This is because the mean annual Niño-3 high-level cloud cover in ISCCP is just 4.8%, so only small negative anomalies are possible before there are no high clouds remaining. This nonlinearity in dSW/dTCC may partly account for the αSW nonlinearity in the reanalyses.
Having shown that the dynamical response plays an important role in the modeled αSW nonlinearities, can biases in the dynamical, cloud, or SW flux responses account for the modeled and errors? Table 4 presents the linear correlations between the feedbacks and individual responses for the 12 CMIP3 models. For SSTA > 0, only dTCC+/dω500 (i.e., the total cloud cover response to dynamics for positive SSTAs) exhibits a significant relationship with . For instance, the positive feedback in MRI CGCM2.3.2 can be explained by the positive dTCC+/dω500 value, and the strongest negative dTCC+/dω500 value in CCCMA could explain why this model has a negative value despite its weak SST variability and dynamical response (Table 2). However, we note that this intermodel correlation cannot explain why the models underestimate with respect to OAFlux; earlier analysis suggests that the other two responses also play a role (see section 4 and Table 3).
On the other hand, there is no significant relationship between and any of the responses for SSTA < 0 (Table 4). It is therefore likely that biases in all three responses contribute to the intermodel differences (Fig. 7). A full analysis of the model errors would require a detailed investigation of the dynamical and cloud changes during La Niñas—a possible area of future study.
We have therefore found no simple answer to question (i) posed at the beginning of this study, “Do the atmospheric model biases that account for the AMIP αSW errors also apply to the coupled αSW errors?” In the AMIP simulations, cloud-related biases were found to play a dominant role in the SW flux response differences (Lloyd et al. 2011). Although the cloud response to dynamics does exhibit a relationship with in the coupled simulations (Table 4), analysis suggests that an underestimated dynamical response to SST changes is key to explaining the underestimated coupled αSW feedbacks (section 4 and Table 3). Question (ii), “Does the SW flux feedback exhibit a nonlinearity that is masked by the linear αSW definition?”, can be answered more confidently: large SW flux nonlinearities have been revealed, characterized by a negative (positive) SW flux feedback during El Niños (La Niñas) in OAFlux and most models. This nonlinearity is governed by a nonlinear dynamical response to SST anomalies (Table 2).
6. The effect of model coupling on αSW
In this section we investigate the third question posed at the beginning of the study: “What explains the coupled versus AMIP differences in the αSW feedback values?” Figure 9 shows the percentage changes in the overall coupled dω500/dSST, dTCC/dω500, and dSW/dTCC values compared to the corresponding AMIP values for the eight models that have an AMIP simulation and supply all required fields. Percentage changes, rather than simple differences, are calculated so as to allow a direct comparison of the three response changes.
In all models except MIHR, dω500/dSST exhibits a negative percentage change, that is, the coupled simulations have a weaker dynamical response to SST than the AMIP simulations. Percentage changes in the cloud response to dynamics, dTCC/dω500, range from 69.8% (MRI CGCM2.3.2: stronger positive coupled cloud response) to −135.3% (HadGEM1: negative AMIP cloud response changing to a positive coupled response), whereas percentage changes in the SW flux response to cloud cover, dSW/dTCC, are generally smaller but have values of 42.7%, −85.1%, and −25.2% in MIMR, IAP, and IPSL CM4, respectively.
Figure 9 thus highlights coupled–AMIP changes in all three responses. The reason for the dω500/dSST changes is discussed below. The changes in dTCC/dω500 and dSW/dTCC remain to be understood, although they are likely to be related to shifts in the mean cloud properties in the coupled simulations, as the atmosphere models (and hence parameterizations) remain the same. Furthermore, shifts in the location of the response patterns may also alter the regression values, as the averaging region (Niño-3) is kept the same for the coupled and AMIP simulations.
Figure 10 presents the relationship between the response percentage changes (Fig. 9) and the coupled–AMIP αSW differences. Linear correlation coefficients are −0.70, −0.65, and 0.57 for the dω500/dSST, dTCC/dω500, and dSW/dTCC changes, respectively, all significant only at the 0.1 level. The dSW/dTCC relationship comes mainly from IAP (the correlation ignoring this model is only 0.33). Similarly, the dTCC/dω500 relationship is set by the two outliers, HadGEM1 and MRI CGCM2.3.2 (the correlation ignoring these models is only −0.27). On the other hand, the dω500/dSST relationship is more robust: there are no clear governing models. The coupled simulations that exhibit the largest dynamical response change compared to their AMIP counterparts also tend to have a larger coupled − AMIP αSW difference. An outlier is the (flux corrected) MRI CGCM2.3.2, which exhibits a 33% weakening of the dynamical response in the coupled simulation, but the smallest change in αSW. The correlation between the dω500/dSST and αSW changes is −0.85 if MRI CGCM2.3.2 is ignored, significant at the 0.01 level. This relationship still holds if the dω500/dSST coupled–AMIP differences are used instead of the percentage changes (correlation of 0.77 ignoring MRI CGCM2.3.2, significant at the 0.05 level).
It therefore appears that differences in the dynamical response to SST are the main cause of the αSW differences between the AMIP and coupled simulations. This is not a surprising result because the atmospheric dynamics are constrained by prescribed SST forcing in the AMIP simulations but can interact with the ocean in the coupled simulations.
The dynamical response to SST: Relationship with the mean state
Having shown that coupled–AMIP changes in the dynamical response to SST are related to the coupled − AMIP αSW differences (Fig. 10), can these dynamical response changes be linked to the model mean states? Table 5 presents the mean Niño-3 ω500 values in ERA-40 and the CMIP3 models (coupled and AMIP simulations). In ERA-40, the Niño-3 region is characterized by mean subsidence (average ω500 > 0). All model simulations (except the coupled MIHR simulation) also exhibit mean subsidence, and, as expected, the coupled simulations have a larger spread in the mean ω500 values than the AMIP simulations [intermodel rms errors (with respect to ERA-40) of 7.6 and 2.6 hPa day−1 in the coupled and AMIP simulations, respectively].
Figure 11 is a scatterplot of the coupled–AMIP percentage changes in dω500/dSST (blue bars in Fig. 9) against the coupled–AMIP percentage changes in the mean Niño-3 ω500. Those models with positive (negative) ω500 coupled–AMIP percentage changes have stronger (weaker) mean subsidence in the coupled simulations. Therefore, the strong negative relationship in Fig. 11 (correlation coefficient of −0.82, significant at the 0.02 level) indicates that models with a larger mean subsidence increase in the coupled simulations (e.g., HadGEM1 and CNRM-CM3) exhibit a larger weakening in the dynamical response. This relationship still holds if the differences (rather than percentage changes) in the dω500/dSST and mean ω500 values are considered (correlation coefficient of 0.76, significant at the 0.05 level).
This link between a stronger mean subsidence and a weaker negative dynamical response is also found among the 12 CMIP3 coupled simulations: linearly regressing the mean Niño-3 ω500 coupled values (Table 5) against the dω500+/dSST and dω500-/dSST coupled values (Table 2, column 2) gives correlation coefficients of 0.50 and 0.90, respectively (significant at the 0.1 and 0.001 levels). The fact that the strongest relationship is found for the dynamical response to negative SST anomalies suggests a “saturation” in the modeled subsidence (i.e., models with a stronger mean subsidence exhibit a weaker anomalous subsidence during La Niñas). An in-depth analysis of La Niña in the coupled and AMIP simulations would be needed to find out if a similar mechanism explains the relationship in Fig. 11.
7. Cloud properties over the eastern equatorial Pacific
It was shown by Lloyd et al. (2011) that studying the modeled cloud radiative forcing (CRF) can highlight cloud biases over the eastern equatorial Pacific that impact αSW (unfortunately, the CMIP3 database does not provide high- and low-level cloud cover data for a direct comparison with ISCCP). Figure 12 shows scatterplots of shortwave (CRFSW) versus longwave (CRFLW) cloud radiative forcing (Niño-3) for the observations and coupled simulations [cf. Lloyd et al.’s (2011) Fig. 8 for the AMIP simulations]. Low clouds in these plots are positioned near the y axis, that is, where CRFLW is low (small greenhouse effect due to warm cloud-top temperatures) and CRFSW is between −20 and −50 W m−2 in ISCCP. These points correspond to the months July–December, the time of the year when low-level marine boundary layer clouds are most prevalent over the eastern tropical Pacific (Klein and Hartmann 1993).
There is a large variety of CRF behavior in the coupled models, though some of the AMIP biases highlighted by Lloyd et al. (2011) are carried through to the coupled simulations. First, the lowest clouds (points closest to the y axis) have an overly negative CRFSW in all coupled simulations except CCCMA. This bias suggests that the lowest clouds are too extensive and/or too optically thick, as also noted in the HadGEM1, IPSL CM4, and MRI CGCM2.3.2 AMIP simulations (Lloyd et al. 2011, Fig. 8). Furthermore, most coupled models also exhibit a shift of these points away from the y axis, suggesting errors in the relative amounts of high and low clouds. Only HadGEM1, ECHAM5/MPI-OM, and MRI CGCM2.3.2 underestimate the mean Niño-3 CRFLW (Table 6, column 3), as also indicated by their shift of points toward CRFLW = 0 in Fig. 12. This suggests that these three models simulate too many low clouds (or not enough high clouds), a bias that could explain the mean positive cloud response in MRI CGCM2.3.2 and HadGEM1 but not the strong negative cloud response in ECHAM5/MPI-OM (Table 2, column 3).
Another AMIP result present in the coupled simulations is the largest overestimation of the average Niño-3 CRFLW in CNRM-CM3, MIMR, and MIHR [Table 6 (column 3), Fig. 12]. These models also overestimate the mean Niño-3 total cloud cover and have the strongest average CRFSW (Table 6, columns 2 and 4), indicating the presence of too many high clouds. However, despite this, these models underestimate the negative cloud response to dynamical changes (Table 2, column 3). Furthermore, CCCMA, the model with the strongest dTCC/dω500 value, has one of the better CRF simulations compared to ISCCP (Table 6; Fig. 12i).
The lack of a clear relationship in the coupled simulations between CRFLW or CRFSW and dTCC/dω500 is supported by the weak correlation coefficients between the average Niño-3 values (−0.16 and −0.24 for CRFLW and CRFSW, respectively; not significant). Scaling CRFLW and CRFSW by the mean total cloud cover, to correct for models that simulate too many—or not enough—clouds in this region, does not alter this result (new correlations of −0.30 and −0.21, respectively; not significant). The absence of a link between the mean cloud properties and dTCC/dω500 in this analysis suggests that the cloud cover response biases may have their source in the complex model cloud schemes and/or convective parameterizations.
8. Summary and discussion
Motivated by previous studies revealing that the shortwave flux feedback over the eastern equatorial Pacific (αSW) is the primary contributor to model errors in the overall ENSO heat flux feedback, α (Lloyd et al. 2009, 2011), this study presents a detailed analysis of αSW in 12 coupled CMIP3 simulations.
To understand the source of the αSW errors, a new feedback decomposition method is introduced, breaking down the SW flux feedback into three individual, local responses: 1) the dynamical response to SST dω500/dSST, 2) the total cloud cover response to dynamics dTCC/dω500, and 3) the SW flux response to clouds dSW/dTCC.
It is shown that all coupled models underestimate the dynamical response during El Niños, a behavior that is likely to contribute to the underestimated αSW values (Fig. 5, Table 2 (column 2), Table 3). Dynamical biases play a more important role in the coupled simulations than the AMIP simulations, in which cloud-related biases were found to be the main source of αSW errors (Lloyd et al. 2011). This is because the large-scale circulation is no longer constrained by a prescribed SST forcing in the coupled simulations.
Changes in the dynamical response, dω500/dSST, between the coupled and AMIP simulations exhibit a robust statistical relationship with the coupled–AMIP αSW differences (Fig. 10). Furthermore, it is shown that the coupled–AMIP differences in dω500/dSST can be directly related to the dynamical mean state coupled–AMIP changes (Fig. 11). The coupled versus AMIP differences in the αSW values are therefore linked to a shift in the dynamical mean state when the ocean and atmosphere models are coupled, a result which further underlines the important role of the dynamics in the coupled αSW feedback.
Nevertheless, there are also large coupled model biases in the mean cloud properties over the equatorial Pacific, as supported by previous studies (Bony and Dufresne 2005; Zhang and Sun 2006; Sun et al. 2006, 2009) and found in the AMIP simulations (Lloyd et al. 2011). Biases in both the cloud response to dynamics, dTCC/dω500, and the SW flux response to clouds, dSW/dTCC, are likely to contribute to the SW flux response biases during model El Niños (Table 3). An improvement in the coupled αSW feedback will therefore only be possible with an improved simulation of both dynamical and cloud responses to SST variability in the eastern equatorial Pacific. Analysis of the model cloud properties will be facilitated in the next-generation CMIP5 simulations by the availability of cloud variables that can be directly compared to satellite observations such as ISCCP.
The observed and modeled SW flux feedback is shown to exhibit a nonlinearity that is hidden when linearly calculating the overall αSW feedback (section 5). Most models exhibit a negative (positive) SW flux feedback during El Niños (La Niñas), as found in OAFlux (Fig. 7). This nonlinearity explains why the SW flux behavior during model El Niños is not always consistent with the αSW values (Fig. 4, Table 3). In the reanalyses/observations, two physical mechanisms are proposed to explain this nonlinearity: 1) a nonlinear dynamical response to SST (Table 2), with strong anomalous ascent for positive SST anomalies and weak anomalous subsidence for negative SST anomalies, and 2) a nonlinear high-level cloud cover response to SST (Fig. 8). The SW flux feedback nonlinearity tends to be underestimated by the models (Fig. 7), linked to an underestimated nonlinearity in the dynamical response to SST (Table 2).
In the models, the main source of the SW flux feedback nonlinearity is a nonlinear dynamical response (section 5). The strong relationship between the dynamical response and the mean dynamics in the CMIP3 models (Fig. 11) therefore suggests that the model mean state may play a role in the SW flux feedback nonlinearity. Studying the evolution of the dynamics during modeled La Niñas would help to elucidate this hypothesis. An analysis of the modeled cloud response to negative SSTAs may also reveal a cloud nonlinearity such as that in the observed high clouds (Fig. 8).
A further understanding of the nonlinear SW flux response to SST thus requires a detailed analysis of the αSW mechanisms during La Niñas. This may enable us to link the SW flux feedback nonlinearity to the modeled ENSO SST skewness, that is, the relative strength of El Niño and La Niña situations. Preliminary analysis also suggests the existence of a strong relationship between the modeled ENSO amplitude and the SW flux feedback nonlinearity, the reason for which is not yet known.
These results demonstrate that the linear α parameterization used in simple models and frameworks is highly idealized, masking not only the individual heat flux feedback components but also a large SW flux feedback nonlinearity. Although the use of a linear α term in the simplest ENSO models is somewhat justified, as the observed net heat flux response is (fortuitously) close to a linear relationship [e.g., the Kim and Jin (2011) Fig. 6], it should always be borne in mind that the linear parameterization hides a number of important processes. Care should also be taken when analyzing α in ENSO frameworks such as the BJ index, as an accurate α value may result from compensating errors in the heat flux feedback components.
Another caveat associated with α is that it is defined as a local feedback. While this is likely to be a reasonable assumption for the LH flux component driven by humidity processes in the eastern equatorial Pacific region (Lloyd et al. 2011), it may not be true for the SW flux feedback. An important aspect of αSW is the dynamical response to SST anomalies, dω500/dSST. The dynamics over the equatorial Pacific are associated with two large-scale atmospheric circulations: the meridional Hadley circulation and the zonal Walker circulation. Therefore, considering dω500/dSST as a local response is an oversimplification: the dynamics over the eastern equatorial Pacific will also be affected by nonlocal SST changes.
The feedback decomposition method introduced in this study represents a first step to understanding the model αSW biases. Future work will involve developing a more complex framework that takes into account the nonlocal component of dω500/dSST, as well as any cross-correlations between the SW flux, cloud cover, and dynamical fields. This will equip us with an even more effective tool for analyzing the SW flux feedback associated with ENSO.
We propose that using the diagnostics developed in this study in conjunction with the BJ index (Jin et al. 2006; Kim and Jin 2011) and the previously developed atmospheric feedback diagnostics (Lloyd et al. 2009; Guilyardi et al. 2009a; Lloyd et al. 2011) will provide a powerful method for understanding the source of ENSO amplitude biases in the next generation of GCMs. The CMIP5 simulations, the results from which will contribute to the Intergovernmental Panel on Climate Change Fifth Assessment Report, present a first opportunity to combine these approaches. Ultimately, it is hoped that an improvement in the modeled atmospheric feedbacks will in turn lead to an improved model representation of ENSO and, hence, an increased ability to predict the evolution of this complex phenomenon over the coming decades.
We thank Sandrine Bony, Florent Brient, and Adam Scaife for useful discussions during the course of this work, as well as support from the European Union EUCLIPSE project (ENV/244067, FP7). JL acknowledges support by a CASE grant from the Met Office. We also acknowledge the modeling groups the Program for Climate Model Diagnosis and Intercomparison (PCMDI) and the WCRP’s Working Group on Coupled Modelling (WGCM) for their roles in making available the WCRP CMIP3 multimodel dataset. Support of this dataset is provided by the Office of Science, U.S. Department of Energy.