Modulation of Bjerknes feedback on the decadal variations in ENSO predictability

Clear decadal variations exist in the predictability of the El Niño–Southern Oscillation (ENSO), with the most recent decade having the lowest ENSO predictability in the past six decades. The Bjerknes Feedback (BF) intensity, which dominates the development of ENSO, has been proposed to determine ENSO predictability. Here we demonstrate that decadal variations in BF intensity are largely a result of the sensitivity of the zonal winds to the zonal sea level pressure (SLP) gradient in the equatorial Pacific. Furthermore, the results show that during low‐ENSO predictability decades, zonal wind anomalies over the equatorial Pacific are more linked to SLP variations in the off‐equatorial Pacific, which can then transfer this information into surface temperature and precipitation fields through the BF, suggesting a weakening in the ocean‐atmosphere coupling in the tropical Pacific. This result indicates that more attention should be paid to off‐equatorial processes in the prediction of ENSO.


Introduction
Developing the ability to predict El Niño-Southern Oscillation (ENSO) events one or more seasons in advance is of great societal importance but remains a challenging task [Latif et al., 1998;Kirtman and Schopf, 1998;Chen et al., 2004;Chen and Cane, 2008;Jin et al., 2008;Barnston et al., 2012]. Extended-range forecasts of ENSO had met with encouraging successes during the 1990s [Latif et al., 1998;Kirtman and Schopf, 1998;Chen et al., 2004]. However, a number of analyses based on numerical model predictions indicate that there exist decadal variations in ENSO predictability [Chen et al., 2004;Chen and Cane, 2008;Jin et al., 2008;Barnston et al., 2012]. By investigating the ENSO prediction skill of 20 state-of-theart models, it was concluded that the forecasting reliability during the period 2002-2011 was relatively lower than that in the 1980s and 1990s; even the models used then were less advanced [Wang et al., 2010;Barnston et al., 2012;Xue et al., 2013]. This raises the question of why our ability to forecast ENSO events is declining even while our understanding of ENSO mechanisms as well as the ENSO models is significantly improving. This decline is deeply rooted in ENSO predictability itself and also possibly due to the imperfections of the forecast systems [Kirtman and Schopf, 1998;Chen et al., 2004;Chen and Cane, 2008;Barnston et al., 2012]. Based on signal-to-noise principles, it is generally believed that ENSO predictability is related to the magnitude of ENSO variability (i.e., the larger its variability, the less likely that the ENSO signal is disrupted by other physical processes, making it more easily simulated and forecast) [Kao and Yu, 2009;Wang et al., 2010;Barnston et al., 2012;Hu et al., 2013]. Therefore, an analysis of the processes that induce decadal variations in ENSO magnitude may lead to an understanding of the decadal variations in ENSO predictability.
It has been suggested that the amplitude, period, and sustainability of ENSO are determined by the oceanatmosphere coupling strength in the tropical Pacific [Capotondi et al., 2015;Yu and Kao, 2007]. Physically, this coupling strength can be largely interpreted as the strength of the so-called Bjerknes feedback (BF) [Bjerknes, 1969], which describes the sensitivity of the atmospheric response to ocean forcing and vice versa. In fact, the strength of the BF has been used as a criterion to measure the instability of the coupled ENSO mode in both observations as well as model simulations [Jin et al., 2006;Jin, 2011a, 2011b;Lübbecke and McPhaden, 2014]. In this study, a possible relationship between decadal variations in the BF intensity and ENSO predictability is examined by quantitatively characterizing the physical processes involved in the BF using observational and reanalysis data sets.

Data Sets, Reference Model, and Methodology
Several observational and reanalysis data sets are used in this study. The zonal wind at 850 hPa (hereinafter u_wind) and sea level pressure (SLP) data are from NCEP/NCAR (National Centers for Environmental Prediction/National Center For Atmospheric Research) Reanalysis 1 products [Kalnay et al., 1996], and the sea surface temperature (SST) data are from ERSST (NOAA Extended Reconstructed SST) v3b [Smith et al., 2008]. The analysis period for these data sets is from 1950 to 2014. A shorter period of 1979 to 2013 is used for the precipitation data from GPCP (Global Precipitation Climatology Project), v2.2 [Adler et al., 2003]. Anomalies are defined as the deviations from the mean seasonal cycle of each field. To confirm the major results of this study, other data sets were also used for calculation in the supporting information.
An ENSO ensemble prediction system (EPS, referred to as REF model in this article) is used to measure the practical ENSO predictability [Zheng et al., 2009]. The EPS used here consists of three main components: (1) an intermediate coupled model (ICM) [Zhang et al., 2005]; (2) an air-sea coupled data assimilation system that is used to minimize errors in both the initial atmospheric and oceanic conditions for the model through assimilation of available atmosphere and ocean observations simultaneously [Zheng and Zhu, 2010]; and (3) a linear first-order Markov stochastic model embedded within the ICM to simulate the time evolution of model uncertainties in the forecasted SST anomaly fields and prolong the period of ENSO prediction [Zheng et al., 2006[Zheng et al., , 2009. A 20 year retrospective forecast comparison demonstrates the skill of the EPS system in predictions with lead times of up to 1 year Zhu, 2010, 2016].
A persistence forecast is the simplest type of climate forecast, as it requires only knowledge of current conditions to forecast conditions in the coming months. In this study, we used the persistence forecasts to gauge the skill of the forecasts made by the REF model. We first documented the skill of the persistence forecast by applying a moving lagged autocorrelation analysis with the observed Niño3.4 SST index. A moving 120 month window is shifted each year from 1950 to 2014 to calculate the correlation coefficient. A detailed description can be found in Yu and Kao [2007]. The forecast skill of the REF model is measured in a similar way by calculating the 120 month running simultaneous correlation coefficients between the forecasted and observed Niño3.4 indices. As the ENSO persistence is used as an indication of the ENSO predictability, the correlation coefficients between model forecast skill and the persistence indicates that the model forecast skill is largely regulated by the ENSO predictability. In this work, we intended to measure the general ENSO forecast skill during all 120 months a decade, rather than focusing on the predictability during a particular season.
As data in geophysical time series are not always completely independent from each other, the persistence along with the finite size of the samples must be taken into account when computing statistics [Trenberth, 1984;Bretherton et al., 1999]. In this method, one can derive the degrees of freedom for significance tests using the correlation between two times series X i and Y i with different autocorrelation sequences ρ X τ and ρ Y τ as follows [Bretherton et al., 1999]: where T Ã XY is the number of degrees of freedom used in the significance calculation and T is the unadjusted number of degrees of freedom. Figure 1 displays the decadal variations in the prediction skill of persistence forecasts and the REF forecasts during the last six decades. The decadal variations in skill for these two forecast approaches are very similar, and the REF forecast skill has a high correlation with the persistence forecast skill at leads of 1-3 (0.86) and 4-6 months (0.82). This similarity suggests that persistence forecast skill can be used as a proxy for model forecast skill [McPhaden, 2003]. It is also in agreement with several previous studies [Kao and Yu, 2009;Wang et al., 2010;Barnston et al., 2012;McPhaden, 2012], which analyzed both physical and statistical models with vast differences in complexity, to suggest that an investigation of the decadal variations in the Geophysical Research Letters 10.1002/2016GL071636 persistence forecast skill can provide a reasonable understanding of the overall ENSO predictability. The striking decline of the persistence forecast skill in the recent decade is also clearly exhibited in Figure 1.

Consistency in the Decadal Variability of ENSO Predictability and BF Intensity
Recalling the simple signal-to-noise principle, a larger standard deviation (SD) in the Niño3.4 (5°S-5°N; 120°-170°W) SST index (i.e., signal) constitutes a larger signal-to-noise ratio, which implies that the ENSO anomalies cannot be easily disrupted by the noise and results in a more persistent signal [McPhaden, 2003;Yu and Kao, 2007]. To quantitatively verify this inference, the standardized decadal variations in the Niño3.4 SD and the ENSO persistence forecast skill are presented in Figure 2, which clearly shows that these forecasts have a high coherency across all six decades, with the correlation coefficient between the Niño3.4 SD and the averaged one to three leading months persistence forecast skill being as high as 0.84. With regard to the physical mechanism, the Niño3.4 SD relies heavily on the coupling strength of the tropical Pacific ocean-atmosphere system, which in turn depends primarily on the most important positive feedback in this region (i.e., the Bjerknes feedback).
The Bjerknes positive feedback describes how an initial positive SST anomaly in the equatorial eastern Pacific can reduce the east-west SST gradient and hence the strength of the Walker circulation to result in weaker trade winds along the tropical Pacific. These weaker trade winds can then further enhance the SST anomaly through a series of ocean dynamical processes [Bjerknes, 1969;Jin, 2011a, 2011b;Liu et al., 2014;Hu et al., 2015], including the zonal advective feedback (ZA), the local upwelling feedback (EK), and the thermocline feedback (TH). Here we choose to measure the intensity of the BF as the regression coefficient between the zonal gradient (difference between eastern and western Pacific regions) of SST anomaly (dSST/dx) and the zonal wind (u_wind) anomaly along the equatorial Pacific (i.e., dSST/dx → u_wind) [Zheng et al., 2014;Fang et al., 2015]. In the interpretation, the EP (eastern Pacific) region is defined over the Niño3.4 area, which can characterize the evolution of both types of El Niño and is used operationally to monitor the onset and define that amplitude of ENSO events by the climate  Figure S1 in the supporting information. As expected, the BF intensity correlates well with the Niño3.4 SD and persistence forecast skill during the entire period. Specifically, the correlation coefficient between the BF intensity and the persistence forecast skill at 1-3 months lead is as high as 0.63 (passing the 10% significance level with 8 degrees of freedom on the basis of a Student's t test), which supports the hypothesis that the decadal variations in ENSO predictability are strongly associated with the variations in the BF intensity in the tropical Pacific (Table S1 in the supporting information). To justify our method of quantifying the BF intensity (i.e., the overall wind response to the zonal SST gradient), we also explored quantifying the BF intensity using the ocean dynamical processes involved in the Bjerknes feedback (i.e., ZA, EK, and TH) (Text S2 in the supporting information). We found that these ocean dynamical processes show less synchronization with the ENSO predictability than the BF intensity defined in this study. The related comparisons and the detailed discussion of the methodology can be found in Tables S2 and S3 [Whitaker et al., 2004;Allan and Ansell, 2006;Compo et al., 2006;Carton and Giese, 2008;Jin, 2011a, 2011b;Tokinaga and Xie, 2011;Liu et al., 2014]. To further explore the atmospheric component of the BF processes, we choose to use wind variations at 850 hPa for the analysis, because it is known to be highly correlated with surface winds and at the same time reflect well the variation in the Walker circulation [e.g., Chung and Li, 2013;Wang et al., 2013;Xiang et al., 2013Xiang et al., , 2014. The BF intensity calculated using the 850 hPa wind is more closely related to the decadal variation of the ENSO predictability (i.e., the correlation coefficient is 0.63) than that calculated using surface wind stress or surface wind velocity (i.e., their respective correlation coefficients with the persistence forecast skill at leads 1-3 months are 0.47 and 0.49, and all pass the 10% significance level). The decrease in the correlation when adopting the surface wind/wind stress can be related to turbulent properties of the air-sea interface [e.g., Jenkins, 2007].
One way to identify the cause of the decadal variations in the ENSO predictability is to examine how the BF intensity changes from decade of high predictability to decade of low predictability. There are three major subprocesses that ultimately determine the BF intensity [Zheng et al., 2014;Fang et al., 2015]: (1) the response in total diabatic heating (Q) to ENSO SST anomalies (subprocess 1, dSST/dx → dQ/dx). The diabatic heating is composed of the latent heating from precipitation and the solar radiative heating related to clouds and water vapor. Very often the diabatic heating is represented only by the latent heating from precipitation (i.e., dQ/dx ≈ dPrec/dx), due to the lack of cloud data or the assumption that the contribution from the cloud radiative heating is small [Lin, 2007], (2) the sea level pressure (SLP) changes due to the heating (subprocess 2, dQ/dx → dSLP/dx), and (3) the zonal wind changes resulting from the SLP gradient changes (subprocess 3, dSLP/dx → u_wind). These subprocesses can be illustrated as follows [Bjerknes, 1969;Lin, 2007]:

Key Processes on Affecting the Decadal Variations in ENSO Predictability
In Figure 3, the monthly anomalies of each variable involved in the Bjerknes feedback are presented during the high-ENSO predictability decade (black circles) and the low-predictability decade (red asterisks) along with the corresponding linear fit lines, whose slopes represent the sensitivity between two specific variables (see also Figure S2 in the supporting information). As shown in Figure 3a, a significant asymmetry in the BF intensity appears between these two decades, with the coupling intensity being significantly larger in the high-ENSO predictability decade (a slope of 45.44 AE 3.19 (m/s)/(°C/deg) with 13 degrees of freedom) than that in the low predictability decade (a slope of 34.90 AE 3.35 (m/s)/(°C/deg) with 20 degrees of freedom). The slopes between the two decades are not significantly different for the first (dSST/dx → dPrec/dx; Figure 3b) and second (dPrec/dx → dSLP/dx; Figure 3c) subprocesses. However, the slopes are different for the third subprocess (dSLP/dx → u_wind, as shown in Figure 3d); that is, in response to a given SLP anomaly gradient, the basin-wide mean u_wind anomaly is stronger during the high-ENSO predictability decade (a slope of À45.47 AE 6.36 (m/s)/(hPa/deg) with 12 degrees of freedom) than during the low-predictability decade (a slope of À28.10 AE 6.41 (m/s)/(hPa/deg) with 27 degrees of freedom). This phenomenon is consistent with the intensity asymmetry in the BF. Figure 3 indicates that the different response of equatorial zonal wind to the SLP gradient (i.e., dSLP/dx → u_wind) is primarily responsible for the BF intensity variation between the decades of lowest and high ENSO predictability. It should be noted that there are also nonlinear processes in the BF processes, especially those involved in the SST-to-precipitation relationship for moderate and severe ENSO events. Most of the black dots in the top right part of Figure 3b (the major nonlinear component) came from the 1997-1998 ENSO event. This nonlinear component cannot be fully covered by the linear method and is not the focus of this work.
To further confirm the finding obtained from contrasting these two periods, Figure 4a shows the variations in the BF intensity and the wind-to-SLP-gradient sensitivity during the past six decades. It shows that these two quantities are indeed inversely correlated with each other throughout the analysis period, with a correlation coefficient as high as À0.89. To uncover the mechanism for the decadal variations in the wind-to-SLPgradient sensitivity, a comparison analysis is performed between high and low-ENSO predictability decades. In addition to the original pair of the decades, we added in the analysis another decade of low ENSO predictability (1959ENSO predictability ( -1968see Figure 1) and the neighboring decade (1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981) of high ENSO predictability (Figure 4 and Figure S3 in the supporting information). Figures 4b to 4e show the spatial distribution of the correlation coefficients between the zonal wind index and SLP anomalies for the four selected decades. For the positively correlated SLP anomalies, it can be seen that the maximum points are located in the equatorial western Pacific and with small differences among the four decades. However, for the negatively correlated SLP anomalies, there are significant differences between the high-and low-predictability decades. During both two decades of low predictability (Figures 4c and 4e), the maximum points of the negative correlations are separately located in the off-equatorial Pacific regions, implying that the zonal wind anomalies over the equatorial Pacific are more strongly related to the SLP variations off the equator. This is consistent with a closer relationship to the Hadley circulation than to the Walker circulation during these decades. This indicates that the weak ENSO variability during these two decades resulted from the interactions between the tropical Pacific Ocean and the Hadley circulation. However, during the two high-ENSO predictability decades (Figures 4b and 4d), the major negative correlation centers are more limited in the equatorial eastern Pacific regions, although there is also a minimum center over the south Pacific region in 1992-2001 period, indicating that the anomalous zonal wind over the equatorial Pacific is mostly regulated by the SLP gradient along the equatorial Pacific that is related to the Walker circulation, which suggests that the strong Geophysical Research Letters 10.1002/2016GL071636 ENSO variability in these two periods resulted from the strong coupling between the tropical Pacific Ocean and the Walker circulation. The comparison between the high and low ENSO predictability reveals that the ocean-atmosphere coupling within the equatorial Pacific is weaker in the low-ENSO predictability decades than in the high-predictability decades. This finding could explain why the ENSO variability becomes weaker during the periods of low ENSO predictability, and also suggests that more attention should be paid to processes outside the equator when attempting ENSO predictions during these periods.

Concluding Remarks
The relation between the extratropical variations and the equatorial physical fields can be potentially explained by three mechanisms. First, the mean SLP differences between high-and low-ENSO predictability decades ( Figure S4 in the supporting information) indicate that, during the low-ENSO predictability decades, there are two anomalous high-pressure centers in the extratropics of both the northern and southern Pacific. Note that there are differences among different data sets, which make the mean SLP differences we mentioned here not always obvious in every decade or every data set. For example, the anomalously high pressure center in the northeastern Pacific between the 1959-1968 (low) period and the 1972-1981 (high) period found in the HadSLP2 data set is not as strong as that found in the NCEP data set (Figures S4b and S4e). This difference may be caused by the different methods used in these two reanalysis products to reconstruct SLP information for early decades when SLP observations over oceans were limited. This possibility is supported by the fact that the mean SLP differences calculated from NCEP and HadSLP2 become similar for the recent decades (i.e., 1992-2001Figures 4a and 4d) when SLP observations become more available over oceans. The other possibility is that the mean SLP differences may be one of several possible mechanisms to cause the low-high ENSO predictability. The SLP mechanism may be particularly important for the change of ENSO predictability from 1992ENSO predictability from -2001ENSO predictability from to 2004ENSO predictability from -2013. Nevertheless, our results indicate that if there exist the anomalous high pressure centers in the extratropics, they can impact the zonal wind over the equatorial region by inducing local anticyclones and gradient forces between the equatorial and extratropical regions. Second, it is noted that the eastern Pacific (EP) ENSO variability is associated with a Southern Oscillation between eastern and western tropical Pacific [Bjerknes, 1969]. The SLP variations over the Maritime Continent are linked to those over the eastern equatorial Pacific through the Walker circulation. . BF and its subprocesses in the two representative decades. Scatter diagrams and the corresponding linear fitting lines (with the values of the slopes plus and minus the 95% confidence values from the Student's t test) between (a) wind index and zonal SST gradient (i.e., the BF intensity), (b) zonal precipitation gradient and zonal SST gradient (i.e., the first subprocess), (c) zonal SLP gradient and zonal precipitation gradient (i.e., the second subprocess), and (d) wind index and zonal SLP gradient (i.e., the third subprocess) during high (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001), black circles/lines) and low (2004-2013, red asterisks/ lines) ENSO predictability decades. Each dot in the figure represents the value for 1 month.

10.1002/2016GL071636
For the variability associated with the central Pacific (CP) ENSO [Kao and Yu, 2009], the SLP variations over the Maritime Continent are linked to those in the off-equatorial regions of both the northern and southern Pacific, suggesting that local Hadley circulation may be more important. This link associates ENSO predictability with the dominant ENSO type during a particular period [Yu et al., 2010] ( Figure S5 in the supporting information). Note that other mechanisms have also been proposed, but not invoked here, to explain the existence of the EP and CP types of ENSO, such as that considering these two ENSO types as manifestations of two unstable modes of the tropical Pacific Ocean-atmosphere [Bejarano and Jin, 2008;Yeh et al., 2014]. Third, because the Interdecadal Pacific Oscillation (IPO) is the most prominent decadal variation signal in the region, the correlation coefficients between the IPO index and SLP anomalies are also shown in Figure  S6 in the supporting information, where the IPO index was calculated based on the difference between the SST anomaly averaged over the central equatorial Pacific (170°E-90°W; 10°S-10°N) and the average of the SST anomalies in the northwest (140°E-145°W; 25°N-45°N) and southwest Pacific (150°E-160°W; 50°S-15°S) [Henley et al., 2015]. This figure also indicates two negative centers located in the extratropics of both northern and southern Pacific rather than in the equatorial eastern Pacific. Similar to Figures 4c and 4e, the coherency between the decadal variations in IPO and the persistence forecast skill indicates that variations in IPO also can influence the decadal variations in ENSO predictability to a certain extent. . Time series of BF intensity, its third subprocess, the SLP indices, and the wind index-related SLP distributions. (a) Decadal variations in intensities of BF (red curve) and its third subprocess (i.e., dSLP/dx → u_wind, black curve) and the corresponding error bars (shaded 95% confidence for a Student's t test). (b-e) Correlation coefficients (normalized by the maximum one in this region) between the wind index and SLP anomalies during two representative high (1992-2001 and 1972-1981) and low (2004-2013 and 1959-1968) ENSO predictability decades. The contour interval is 0.1.

10.1002/2016GL071636
It should be noted that our findings introduce an additional mechanism, rather than excluding other possible mechanisms, to explain the decadal ENSO variability. The possibility exists that the slow changes in ENSO activity in recent decades can also be caused by, for example, the positive feedback between ENSO and the ENSO-state-dependent atmospheric noise [Jin et al., 2007;Kug et al., 2008], the modulation effects from the North Atlantic Ocean [Ham et al., 2013;Kang et al., 2014;Yu et al., 2015], the Indian Ocean [Yu et al., 2002;Kug et al., 2006] via subtropical teleconnections or Walker circulation variations, the modulation effects from the Pacific mean state changes Chung and Li, 2013;England et al., 2014] and the weakening of Walker circulation [Bellomo and Clement, 2015], or the influence of westerly wind bursts [Fedorov et al., 2015;Chen et al., 2015].
In this article, we focused on exploring which atmospheric process plays the key role in inducing the decadal variation in the BF intensity and the ENSO predictability due to their high synchronization. This is a different way to investigate the ENSO predictability issues from that of Lübbecke and McPhaden [2014], in which they quantitatively compared the Bjerknes (BJ) index and its components between the two periods, i.e., 1980-1999 (P1) and 2000-2010 (P2), and concluded that the robust change in the total BJ index can be mainly attributed to the weakening of the thermocline feedback regarding both its contribution to the total BJ index and the robust change from P1 to P2. Besides, in this work, we used a simple linear method to measure the BF processes, because it is a common method that was widely used by many other studies. However, it should be noted that there are also nonlinear relationships in the BF processes [e.g., Hoerling et al., 1997;Choi et al., 2013;Levine and Jin, 2015;Takahashi and Dewitte, 2016] that should be studied in the future work.