ENSO feedbacks and their relationships with the mean state in a flux adjusted ensemble

The El Niño Southern Oscillation (ENSO) is governed by a combination of amplifying and damping ocean–atmosphere feedbacks in the equatorial Pacific. Here we quantify these feedbacks in a flux adjusted HadCM3 perturbed physics ensemble under present day conditions and a future emissions scenario using the Bjerknes Stability Index (BJ index). Relationships between feedbacks and both the present day biases and responses under climate change of the mean equatorial Pacific climate are investigated. Despite minimised mean sea surface temperature biases through flux adjustment, the important dominant ENSO feedbacks still show biases with respect to observed feedbacks and inter-ensemble diversity. The dominant positive thermocline and zonal advective feedbacks are found to be weaker in ensemble members with stronger mean zonal advection. This is due to a weaker sensitivity of the thermocline slope and zonal surface ocean currents in the east Pacific to surface wind stress anomalies. A drier west Pacific is also found to be linked to weakened shortwave and latent heat flux damping, suggesting a link between ENSO characteristics and the hydrological cycle. In contrast to previous studies using the BJ index that find positive relationships between the index and ENSO amplitude, here they are weakly or negatively correlated, both for present day conditions and for projected differences. This is caused by strong thermodynamic damping which dominates over positive feedbacks, which alone approximate ENSO amplitude well. While the BJ index proves useful for individual linear feedback analysis, we urge caution in using the total linear BJ index alone to assess the reasons for ENSO amplitude biases and its future change in models.

One of the most important modes of natural climate variability is the El Niño Southern Oscillation (ENSO). Many climate models can now simulate an ENSO cycle with a reasonable level of fidelity, but errors persist in both the mean climate of the tropical Pacific and in the structure, frequency and magnitude of ENSO (Bellenger et al. 2013). The CMIP3 models (Coupled Model Intercomparison Project, version 3), and more recently the CMIP5 models, have been used in many studies to assess ENSO processes and to investigate the impact anthropogenic climate change may have on ENSO, often with inconclusive results . Studies that consider alternative definitions of ENSO events (and extreme ENSO events) in terms of variations in precipitation, seemingly show more robust changes in the future, but may also be influenced by mean biases in models 2015). Errors and biases in models can influence the characteristics of the modelled ENSO and its projected response under climate change. Developing a greater understanding of ENSO feedback processes and their relation to model biases and the uncertainty of long-term ENSO projections is an ongoing priority area of research (e.g. the CLIVAR ENSO in a Changing Climate Research Focus).
ENSO is governed by a suite of air-sea feedback processes in the equatorial Pacific. As described by Bjerknes (1969), the development of El Niño involves weakening of the easterly trade Winds as sea surface temperature (SST) warms in the eastern equatorial Pacific, rendering a reduced zonal SST gradient. The westerly wind anomalies drive weaker surface currents and upwelling, causing the thermocline to tilt eastward. The deepened thermocline in the east leads to increased SSTs which in turn reinforce the westerly wind anomaly, and so forth. These reinforcing processes are termed the Bjerknes feedback. Conversely, negative SSTAs associated with La Niña, the cool ENSO phase, cause an opposite ocean response resulting in further SST cooling. The development of the Bjerknes Stability Index (BJ Index, Jin et al. 2006), allows for a relatively simple method of quantifying these equatorial Pacific feedback strengths. The BJ index includes the damping feedback of the mean ocean currents and damping via reduced atmosphere-ocean heat flux, termed the thermodynamic damping. The positive feedbacks to a warm SST anomaly resulting from anomalous eastward advection in the ocean, weakening east Pacific upwelling and a flattening thermocline, referred to as the zonal advective, Ekman and thermocline feedbacks respectively, are also included in the BJ index. The BJ Index is derived from a linear stability analysis of a simplified model of the upper ocean and atmosphere in the Pacific (see Section 3 for details).
Coupled Atmosphere-Ocean models still struggle to accurately represent ENSO and the feedbacks that control ENSO characteristics. ENSO amplitude and period are often diverse in multi-model ensembles. Here the focus is on the dominant positive and negative ENSO feedbacks, which are found in a number of studies assessing ENSO feedbacks in models to be areas for improvement (Guilyardi 2006;Lin 2007;Lloyd et al. 2009;Kim and Jin 2010a;Bellenger et al. 2013;Kim et al. 2013).
Thermodynamic damping, the strength of the response of surface heat flux to sea surface temperature anomalies in the east equatorial Pacific, is a major source of ENSO diversity. This feedback is usually too weak in models and shows large variations in both CMIP3 and CMIP5 (Lloyd et al. 2009;Kim and Jin 2010a;Kim et al. 2013). This is often primarily related to an underestimated shortwave damping response caused by weak atmospheric ascent in response to east Pacific sea surface temperature anomalies during El Niños as well as a weak response of clouds (Lloyd et al. 2012). Guilyardi et al. (2009) and Bony and Dufresne (2005) also find that the distribution of atmospheric convection is important for an accurate thermodynamic damping.
Another important component of thermodynamic damping is the cooling of SSTs via latent heat flux, which is also found to be weak in many CGCMs and can be linked to biases in wind speed or near surface humidity sensitivity in the equatorial Pacific (Lin 2007;Lloyd et al. 2010).
wind stress to sea surface temperature anomalies, is important. This coupling has been found to be generally weak in models and varies in strength between them (Guilyardi 2006;Lloyd et al. 2009; Kim and Jin 2010a;Kim et al. 2013). The strength of this coupling is slightly improved in atmosphere-only models and CMIP5 relative to CMIP3 (Lloyd et al. 2010;Bellenger et al. 2013;Kim et al. 2013) though the reasons behind this improvement are unclear. Weak positive ENSO feedbacks in models are also related to other ocean-atmosphere feedback loops such as the response of ocean currents to wind stress anomalies, important to the zonal advective feedback, and the response of the thermocline slope to wind stress anomalies which is a key component of the thermocline feedback (Lübbecke and McPhaden 2013;Graham et al. 2014). The biases in these feedbacks can often be linked to errors in the mean equatorial Pacific climate. In particular, the cold tongue bias and weak upper ocean stratification common in models can coincide with weak air-sea couplings (Kim et al. 2013).
Multi-model studies into the response of ENSO to climate change generally do not indicate robust responses, particularly when investigating ENSO amplitude and period Guilyardi 2006;Merryfield 2006;Yeh et al. 2006;Vecchi and Wittenberg 2010;DiNezio et al. 2012). However, recent studies using alternative metrics (e. g. precipitation as an indicator of extreme El Niño events) find a projected increase in strong events 2015). Santoso et al. (2013) also find a projected increase in extreme events when considering zonal SST propagation as an indicator of strong El Niños. These studies highlight the need to also include process-based metrics in ENSO analysis. It has been suggested that initial mean state or the relative dominance of different feedbacks may have an impact on the ENSO amplitude response to climate change (Zelle and Oldenborgh 2005;Merryfield 2006;Yeh and Kirtman 2007).
A stability analysis of CMIP3 models by Kim and Jin (2010)  Similar results have also been found in analyses by DiNezio et al. (2012) and Philip and van Oldenborgh (2006). Despite this, Kim and Jin (2010) find that ENSO stability is a good predictor of ENSO amplitude in CMIP3. In more recent models, ENSO stability climate change responses may be less diverse; Kim et al. (2014) find that a subset of the CMIP5 ensemble, which they judge to most accurately simulate the observed relative importance of feedback terms, do agree on an ENSO amplitude and stability response which rises then falls in later decades. However, the inter-ensemble link between stability and amplitude is less clear in CMIP5 when compared to CMIP3 (Kim et al. 2013).
While the BJ index has proved to be a powerful tool in the assessment of ENSO feedback strength, an ongoing discussion involves the reliability and accuracy of the BJ index as a quantitative predictor of ENSO amplitude. Graham et al. (2014) compare the BJ index with a calculation of the full non-linear mixed-layer heat budget and find that positive feedbacks, associated with ocean processes, tend to be misrepresented by the BJ index. The linear calculation of atmospheric feedbacks is also known to result in underestimated feedback strength in the case of the shortwave damping and the Bjerknes feedback (Jin et al. 2003;Timmermann et al. 2003;An and Jin 2004;Lloyd et al. 2012;Bellenger et al. 2013).
One potential issue is in the area-averaging used in the BJ Index calculation. This involves the choice of various boxes, over which different variables are averaged in order to estimate ENSO feedbacks. CMIP models still feature systematic errors and biases in e.g. the distribution of SSTs, with the cold tongue being too cold and extending too far into the west. Biases are different in different models and, ideally, one should adjust the locations and sizes of the different box-averages in order to accurately capture the different feedback processes (see Kim and Jim 2010a;b). One way of artificially minimising biases in models is to use flux adjustments. While not a way of formally correcting models, the approach can be useful in understanding the role of model biases on e.g modelled variability or the climate change response. Here we exploit a 'perturbed physics' ensemble of models in which flux adjustments are employed to prevent model drift due to radiation imbalances caused by the parameter perturbations (Collins et al. 2011). This has the side effect that the SST distribution in the different model versions is relatively close to that observed, thus alleviating some of the issues of model biases. However, flux adjustment does not solve all problems and there are still some biases in models and some differences in mean state. These are exploited to relate errors in feedbacks to errors in that mean state. The use of a perturbed physics approach results in an ensemble that covers a range of ENSO variability allowing for inter-ensemble relationships to suggest possible causes of feedback biases. While it does not represent the diversity of different mean-state errors, we consider it a stepping-stone to the understanding of the multi-model ensemble. If we can relate feedbacks to mean state errors in an ensemble in which model mean states are close to each other, but in which ENSO characteristics are quite different, this should aid us in the much harder multi-model problem.
Section 2 outlines the perturbed physics ensemble (PPE) used here as well as the reanalysis data used for comparison. The BJ index and its calculation is introduced in section 3. In section 4 the stability analysis results for the HadCM3 PPE during 1985-2015 are given in comparison to the reanalysis data. ENSO feedbacks are also related to the equatorial mean state of the PPE. The ENSO climate change response of the PPE is shown in section 5 and possible reasons for the responses found are suggested. Section 6 briefly assesses the relationship in the PPE between ENSO amplitude and atmospheric noise. Results are summarised and discussed in section 7 and suggestions for future work are given. The analysis carried out here uses a HadCM3 perturbed physics ensemble. Perturbed physics ensembles are generated from a single model and are a relatively easy way of producing a large ensemble, in comparison to a multi-model ensemble, which uses a number of different models.

Data.
Perturbed physics allows for the experiment to be controlled while exploring uncertainties in processes or feedbacks. A single model is used while values of uncertain parameters are perturbed.
Different physical schemes can also be switched in and out as well as perturbing parameters. This approach enables sources of uncertainty in projections due to uncertain parameters (rather than model structure) to be found. The PPE features a fully coupled version of HadCM3 (Collins et al. 2011) with an interactive sulphur cycle. The ensemble consists of a total of thirty-three members, which can be split into two sixteen-member sub-ensembles with one ensemble member having standard HadCM3 parameters. The first sub-ensemble features multiple perturbations to only parameters in the atmosphere component simultaneously (as opposed to single parameter perturbations). The second has perturbations to only schemes and parameters in the ocean component but features standard atmosphere settings. The PPE also uses flux adjustment in order to avoid model drift caused by top of atmosphere imbalances due to perturbations and to improve regional climate change and feedback simulation. The flux adjustment varies spatially and seasonally and is computed during a 'spin-up' phase and fixed (with a seasonal cycle) throughout the integration of all the experiments used here. Thus flux adjustments do not damp SST anomalies, as they do not depend on the value of the SST anomaly.
We analyse output from experiments spanning two hundred years beginning in 1896 with historical forcings followed by the SRES A1B emissions scenario. This scenario describes a future of rapid economic growth and a peak in global population mid-century. Eighteen overlapping de-trended thirty year time periods with start years staggered at ten year intervals (1896-1925, 1906-1835, …, 2066-2095) are used when responses under climate change are analysed. Fields used are sea temperatures to depth of 300m, surface zonal wind stress, surface zonal ocean current, surface meridional ocean current, upwelling ocean current, precipitation, surface sensible heat flux, surface latent heat flux, shortwave radiation and longwave radiation.

Methods.
The BJ index is derived from the linear equation for sea temperature anomalies averaged over the mixed layer (e.g. An et al. 1999), which are then area-averaged over the central and east equatorial Pacific, the areas which are most relevant to ENSO-related variability. Here the mixed layer is taken to have a fixed depth of 50m. This is based on a PPE mean mixed layer depth of 50.5m +/-4.9 m, calculated as the depth of SST-0.5°C (e.g. Philip and van Oldenborgh 2006). A number of linear approximations can then be used to obtain equations which form a simple linear coupled recharge oscillator system with growth rate I BJ such that R represents the collective strength of various ENSO feedbacks. As the BJ index ( I BJ ) represents the growth rate of the system, this means that for I BJ >0 the leading mode of the system is linearly unstable, and when I BJ <0 , the leading ENSO mode is damped.
The first two feedbacks represent the damping feedbacks of mean current damping (CD) and where 〈 〉 E represents area-averaging over the east equatorial Pacific ( The Bjerknes feedback, μ a , is the linear regression coefficient of equatorial Pacific wind stress anomalies against east Pacific SSTAs. β u is found by regressing east Pacific zonal surface ocean current anomalies against equatorial Pacific surface wind stress anomalies, β w is the regression coefficient of east Pacific upwelling ocean current against equatorial Pacific surface wind stress anomalies and β h is obtained by the regression of thermocline slope (represented as the difference between east and west Pacific sea temperatures averaged from the surface to a depth of 300m) anomalies against equatorial Pacific surface wind stress anomalies. Once calculated, the feedbacks are summed together to obtain the total BJ index (example shown in figure 2).
Clearly, the use of area-averaging in the calculation of the BJ index is important and while a method for choosing the longitudinal extents of the area-averaging boxes exists (given in Kim and Jin 2010b) this is not used here as the use of flux adjustment in the ensemble minimises spatial differences between ensemble members. Instead, a fixed division between the east and west Pacific at longitude 180° is chosen for all time periods and ensemble members with an east boundary at 80 °W for east Pacific area-averaging. Averaging over the full basin (e.g. for wind stress in calculation of μ a ) has the same east boundary but a west boundary at 120 °E. Sensitivity tests to the definitions of the dividing boundary of the two regions show that the main conclusions of the paper are not affected by these choices. The role of atmospheric noise is also estimated in section 6 as the residual of wind stress as a function of equatorial Pacific SST (Philip and van Oldenborgh 2009; given as: Where methods of these two studies that result in slight differences in the feedback strengths found. This study mainly follows the original BJ index formula given in Jin et al. (2006), though we use the calculation of μ a (the sensitivity of zonal wind stress to SSTA) given in Kim et al (2013) as we feel this to be more consistent with recent studies (e.g. Graham et al. 2014;Lubbecke and McPhaden 2013). Our analysis is consistent with Graham et al. (2014), however Kim et al. (2013) follows the formula in Kim and Jin (2010a;2010b), which may give a different result. The main difference in the formula used here and the one in Kim et al. (2013) is in the calculation of the mean current damping which includes regression of the SSTAs at the longitudinal and latitudinal boundaries of the averaging box against the full box averaged SSTA as well as consideration of the boundary ocean current, rather than the full box averaged ocean current. The influence of mean upwelling is also not included in the alternative mean current damping calculation.
A large discrepancy between the analysis presented here and that of Kim et al. (2013) lies in the thermodynamic damping, found by Kim et al. (2013) to be approximately 0.15-0.2 y r −1 for the CMIP5 version of HadCM3. Thermodynamic damping for the standard PPE member in this study is significantly larger in strength. Aside from the formula differences outlined above, one of the main differences between these two studies is that the model experiment design differs between them. Kim et al. (2013) use the CMIP5 historical non-flux-adjusted HadCM3 run. Here we use a perturbed physics ensemble of HadCM3 with flux adjustment. Repeating the methods used here on the CMIP5 historical HadCM3 run gives a thermodynamic damping of -0.28 y r −1 ± 0.06 y r −1 , much weaker than that found for the standard parameter HadCM3 PPE member (table 1) in this study.
Other differences include the calculation of anomalies; Kim et al. (2013) uses a 7 year smoothing as opposed to seasonal anomalies used here. Finally Kim et al. (2013) find the separating boundary between east and west Pacific averaging areas by using a method based on SST EOFs, whereas here the boundary here is taken to be fixed at 180°. Therefore the feedbacks will be Positive feedback discrepancies are smaller. Kim et al. (2013)     Similarly, the thermocline feedback can also be linked to the mean equatorial Pacific climate. Figure 3c shows the positive relationship between the thermocline feedback and mean Niño 4' zonal surface ocean current. By decomposing the thermocline feedback (equation (8)) into the sensitivity of the equatorial Pacific thermocline slope to surface wind stress anomalies ( β h ), the sensitivity of equatorial Pacific wind stress to east Pacific SSTAs ( μ a ) and the mean upwelling current (w), the cause of the weak feedback can be found in the sensitivity of the thermocline slope to surface wind stress anomalies, β h , which is significantly weaker in the PPE than the reanalysis. Inter-ensemble variations in the remaining positive feedback, the Ekman feedback (equation (7)) is similarly found to be linked to the zonal surface ocean current (correlation of 0.61, figure 3e).
This can be found to be caused by the dominating component, the upwelling ocean current response to zonal wind stress anomalies ( β w ). PPE members with weaker zonal ocean currents precipitation than the reanalysis in Niño 4 (though the difference is only significant for eleven PPE members; ten of which have atmosphere perturbations), which may be linked to the weak bias in thermodynamic damping as suggested by the correlation of 0.72 between the two. By decomposing thermodynamic damping into individual heat flux damping rates (equation (5)), it is found that the bias is largely caused by weak shortwave (-α SW ) and latent heat flux (-α LH ) damping terms, the dominant components of thermodynamic damping.
Shortwave damping represents the damping due to decreased incoming shortwave radiation in response to warming SSTs during El Niño. This feedback has two regimes in the Tropical Pacific depending on the large-scale circulation. In regions of subsidence, a warm SST anomaly acts to reduce static stability, which breaks up marine stratiform clouds in the area, leading to an increase of solar radiation reaching the surface. This is a positive feedback that typically occurs at higher latitudes (Klein and Hartmann 1993). However, in warmer tropical areas of ascent, such as over the warm pool in the west equatorial Pacific, a positive SST anomaly increases convective cloud cover causing less shortwave radiation to reach the surface and creating a damping feedback (Ramanathan and Collins 1991).
The shortwave radiation damping in the ensemble, which shows the largest bias, is found to be linked to mean SST and precipitation ( the decreasing wind stress sensitivity, μ a , resulting in the reducing thermocline feedback. This is in contrast to a studies by Kim et al. (2015) and Borlace et al. (2013)  and span less of the equatorial Pacific, resulting in a basin-wide decrease in wind stress sensitivity.
The response of this feedback is found to vary in MMEs (Kim and Jin 2011). Philip and van Oldenborgh (2006) suggest that while warmer mean SST can increase wind sensitivity to SSTAs, the more stable atmosphere can counteract this which may explain the variation in its projected response.
Similarly perturbation PPE members and also shows larger variation when atmosphere parameters are perturbed, particularly for shortwave damping. For the positive zonal advective and thermocline feedbacks, the main causes of ensemble bias are weak responses of ocean currents and the thermocline slope to wind stress anomalies in keeping with previous studies (Kim et al. 2013).
Positive feedbacks also tend to be weaker for atmosphere perturbations. This is caused by weaker mean state contributions to the feedbacks (reduced zonal temperature gradient and reduced east Pacific upwelling) and reduced ocean sensitivity. Many of these feedback biases can be related to biases in the mean state equatorial Pacific which then impacts variability. Most importantly, the feedbacks, an addition to the mean climate links to ENSO feedbacks suggested by Kim et al. (2013). PPE members with strong mean zonal advection have an ocean that is less sensitive to El Niño induced wind stress anomalies than in observations, resulting in reduced zonal advective and thermocline feedbacks. Weak shortwave damping and latent heat flux damping, also found by Lloyd et al. (2010) to be the biggest cause of thermodynamic damping biases in AMIP models, are found to be linked to reduced precipitation. It is possible this may signify convection biases causing less cloud cover formation during El Niño (reducing shortwave damping) and may alter the wind speed response (reducing latent heat flux damping).
Understanding the impact of mean state biases on dominant feedbacks also helps in controlled by the thermocline feedback (Kim et al. 2015;Borlace et al. 2013), which in turn is governed by the thermocline response to equatorial zonal surface wind stress ( β h ). For the thermocline feedback here the initial mean state (as governed by perturbation type) has an impact on the climate change response such that ensemble members featuring atmosphere perturbation (initially the ensemble members with least ocean sensitivity) show significant increases in the sensitivity of the thermocline to surface wind stress anomalies under climate change but ensemble members with ocean perturbations tend not to. This results in less weakening of the thermocline feedback for the atmosphere perturbation ensemble members as the increase in thermocline response to zonal surface wind stress ( β h ) tends to balance out the decreasing wind stress response to SSTA ( μ a ) for these ensemble members. This demonstrates an importance of the mean state of the climate to ENSO feedback projection.
In contrast to previous studies which find a positive relationship between ENSO amplitude and the BJ index (Kim and Jin 2010a;2010b;Borlace et al. 2013), a key finding here is that both inter-ensemble differences and climate change response show a weakly or negatively correlated BJ index and ENSO amplitude -both for present-day conditions and for future changes. (Note that for CMIP5 models this relationship also breaks down as outlier models weaken the relationship between BJ index and ENSO amplitude (Kim et al. 2013)). There are a number of reasons why this may be the case.
The first of these relates to the assumption that the amplitude of components of the BJ index can be directly related and therefore can be added together to provide a measure of ENSO stability.  Secondly, the derivation of the BJ Index relies heavily on assumptions of linearity and neglects processes that lie outside of the framework from which the BJ index is derived. For example, we find here significant relationships between ENSO amplitude and atmospheric noise which will not be accounted for in the BJ index calculation and should be considered in future work. It is known that a number of these linearity assumptions are not complete, such as the linearity of the wind stress response (Kang and Kug 2002;Philip and van Oldenborgh 2009;Choi et al. 2013) or the approximation of thermodynamic damping (Lloyd et al. 2010;Bellenger et al. 2013) which may cause inaccuracies in the BJ index (Graham et al. 2014). It is also important to note that at the height of an El Niño event when SSTAs peak is the point at which non-linearities are more likely to occur and the assumptions of linearity break down. The strength of this peak is essentially what is measured by ENSO amplitude, whereas the BJ index aims to quantify the growth rate of the El Niño event in the lead up to this peak, when the linearity assumptions are more feasible, meaning that the two measures may not be as strongly linked as suggested in previous studies.
Bearing these caveats in mind and considering the results found here, the use of the BJ index as a measure of ENSO amplitude, either on an inter-ensemble basis or as a tool for assessing projected ENSO response, remains questionable. Despite this, the BJ index still allows for assessment of linear dynamics of ENSO and proves useful when attempting to relate feedback bias and response to the background climate. Future work should focus on improvement of the understanding of the feedbacks that the BJ index attempts to approximate and ways in which these Results here suggest that even when SST biases are improved (here by the use of flux adjustment) and analysis is confined to a single model framework (by use of a PPE) common CGCM biases (e.g. strong equatorial Pacific zonal ocean currents) have an impact on ENSO feedbacks causing persistent feedback biases, such as weak thermocline and zonal advective feedbacks and thermodynamic damping, and are a cause of uncertainty in projections. The representation of mean zonal advection and precipitation are found to be particularly important.
Links between heat flux feedbacks and precipitation are suggestive of links to the hydrological