1 Introduction

In spite of an overall performance improvement in simulating the world’s climate, coupled general circulation models (CGCMs) still display large biases in the tropical Atlantic (TA), both in the ocean and in the atmospheric circulation. Most CGCMs are unable to reproduce the observed cold tongue in the eastern part of the TA basin and there are biases in the slope of the equatorial thermocline (Richter and Xie 2008; Davey et al. 2002; Chang et al. 2007). Furthermore, there are biases in simulating tropical Atlantic variability (TAV).

A prerequsite of simulating TA climate correctly is that the mechanisms that govern TA climate and its variability are well represented in the models. While the large annual cycle causes most of the SST variablity in the TA, two other leading modes of interannual to decadal variability have been identified: the zonal and the meridional mode. The zonal mode with a period of two to four years is driven by the Bjerknes feedback (in the following abbreviated as BF) (Ruiz-Barradas et al. 2000; Keenlyside and Latif 2007; Janssen et al. 2008), whereas the meridional mode on decadal time scales (Zebiak 1993; Carton et al. 1996; Xie and Carton 2004, e. g.) arises through the wind evaporation sea surface temperature (WES) feedback (Chang et al. 1997). In this paper we will focus on the zonal mode.

Previous studies have identified the BF as driving mechanism of the zonal mode in reanalysis data. Keenlyside and Latif (2007) found the BF to be present in the TA. Janssen et al. (2008) found not only the BF, but also the recharge oscillator to be a dominant driver of TAV. Ding et al. (2010) confirmed the presence of the BF with principal oscillation pattern analysis, and investigated the seasonality of the upper equatorial Atlantic (Ding et al. 2009). A connection between the two modes of variability has been suggested by Foltz and McPhaden (2010), and a connection between the zonal mode and the Benguela Niño has been suggested by Lübbecke et al. (2010). Several studies have investigated the feedbacks in the TA with the help of general circulation models (GCMs) (DeWitt 2005; Mahajan et al. 2009; Muñoz et al. 2012; Ding et al. 2015).

Most state-of-the-art CGCMs simulate the mean state and variability in the TA poorly (Breugem et al. 2006). This has often been attributed to the westerly wind bias along the equator (Richter and Xie 2008). Hazeleger and Haarsma (2005) show that vertical mixing can play an important role and (Li and Xie 2012) propose that tropical-wide SST errors arise from incorrect simulation of cloud cover and thermocline depth. A recent study by Ding et al. (2015) shows that an improved mean state of the TA can lead to better representation of interannual variability. Our study aims to quantify the strength and phase of the air-sea interaction related to the BF in newest generation climate models, and possibly point to weaknesses in the model representation. We do so by considering different aspects of the feedback and perform statistical analysis on reanalysis and model data. Below, we shortly outline the BF.

In the warm phase of the zonal mode, positive sea surface temperature anomalies (SST’) in the eastern part of the basin cause the zonal winds \((\tau _u)\) in the western part of the TA to weaken. The associated weakening of the zonal pressure gradient allows the thermocline in the eastern basin to deepen and thereby increases the oceanic heat content (HC) in the region. This reduces entrainment of colder subsurface waters into the turbulently mixed surface layer. In addition, the weaker winds may lead to less atmospheric cooling by turbulent surface heat fluxes. The result is a strengthening of the positive SST’ in the region, closing the feedback loop. The opposite happens in the cold phase of the zonal mode (Zebiak 1993).

As a reference for the model results we use a set of reanalysis data which has, to our knowledge, not yet been used for the analysis of the zonal mode. We consider monthly as opposed to seasonal time series because of the pronounced annual cycle in the region. It has been suggested that the zonal mode can be understood as an interannual variation of the strength of the cold tongue (Burls et al. 2011, 2012). This leads to the hypothesis that if the annual cycle is simulated incorrectly, either in timing or strength, the representation and mechanism of the zonal mode will be flawed as well. Therefore, we compare the simulated annual cycles to annual cycles in the reanalysis before analysing the BF components. This allows us to a posteriori correct the model output for the errors in their annual cycles and to assess whether the mechanism driving the zonal mode is correctly represented in other months than in the reanalysis. Finally, we analyse the simulated subsurface ocean and its interaction with SST’ as compared to reanalysis.

This paper is structured as follows: in Sect. 2 we describe the data and methodology, in Sect. 3 we investigate the BF in reanalysis data, and in Sect. 4 we compare first the annual cycle (Sect. 4.1), and subsequently the three components of the BF (Sects. 4.2, 4.3) to reanalysis data. In Sect. 4.4 we examine subsurface temperature variance and its correlation to SST’ in the TA, followed by the conclusions of this study in Sect. 5.

2 Data and methodology

The reanalysis data used in this study is obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF). The variables \(\tau _u\) and SST are taken from the global atmospheric reanalysis ERA-Interim (Dee et al. 2011), and the potential ocean temperature from the Ocean Reanalysis System 4 (ORAS4) (Balmaseda et al. 2013). Note that the SST from ERA-Interim and ORAS4 do not differ significantly, since ERA-Interim was used to drive the ocean reanalysis for 20 years out of the 35 years we use in this study. We use time series from 1979 until 2013, due to availability (ERA-Interim) and trustworthiness of the data (ORAS4). For the first years of ORAS4, 1958 until 1979, before the satellite era, little data is available to constrain the reanalysis.

The model output used in this study is obtained from the Coupled Model Intercomparision Project Phase 5 (CMIP5) (Taylor et al. 2012). The CGCMs as well as the markers indicating them in the comparison plots, and the centers where the simulations were performed are shown in Table 1. We use the last 150 years of the pre-industrial control runs to investigate the natural variability in absence of anthropogenic forcing and with minimised model drift. It is noteworthy that the warm SST bias in those simulations is as large as it is. We are comparing the CGCM output to current reanalysis data which already includes the global warming in response to the heightened \(\hbox {CO}_2\) forcing. The CGCM output we use was produced with pre-industrial forcing and should therefore be colder than the reanalysis data.

All data are interpolated to a \(1^{o}\times 1^{o}\) grid and linearly detrended using monthly mean time series.

The upper ocean heat content is calculated grid-point-wise by integrating the potential ocean temperature upwards from the 293.15 K isotherm depth to the surface according to Eq. 1.

$$\begin{aligned} HC (lat, lon) = c_p \cdot \rho \int _{Z_{20}}^{0} T\left( z, lat, lon\right) \mathrm {d}z, \end{aligned}$$
(1)

where we assume a uniform density \(\rho\) of \(1024.75 \frac{kg}{m^3}\) and a specific heat capacity \(c_p\) of \(3993 \frac{J}{kg\cdot K}\) for sea water (Laboratory 2015).

Fig. 1
figure 1

Climatology obtained from ERAInterim reanalysis of SST [K] and \(\tau _u\,[\frac{N}{m^2}]\) in the TA throughout the year

Fig. 2
figure 2

The zonal mode obtained from ERA-Interim averaged over June, July and August. Left correlation between monthly SST in EA4 (right box), and SST elsewhere in the basin. The left box indicates the WA4 index used for other correlation analysis later on. Right linear regression analysis of monthly SST’\(_{EA4}\) on monthly SST’ elsewhere

We use the following index regions as indicators for the western TA (WA4, \(4^\circ \hbox {N}\)\(4^\circ \hbox {S}\); \(40^\circ \hbox {W}\)\(20^\circ \hbox {W}\)) and the eastern TA (EA4, \(4^\circ \hbox {N}\)\(4^\circ \hbox {S}\) and \(20^\circ \hbox {W}\)\(10^\circ \hbox {E}\)), cf. Fig. 2. The variables and correlations are spatially averaged over these index regions to obtain a time series or a mean value. Note that for the sample size used (length of the time series) and degrees of freedom of the data, the threshold value for a 95 % significant correlation is 0.35. The signal can be considered significant when the number of gridpoint correlations exceeding that threshold is larger than 5 %. Note that the variance that is explained by a correlation of 0.35 is still relatively small.

Table 1 Models used in this study, the markers denoting them in the graphics, and the centers where the simulations were run. For more information we refer to CMIP5 (2015)

3 Bjerknes feedback in reanalysis data

In this section we investigate characteristics of the zonal mode in the TA from reanalysis data. Because of the close connection of the zonal mode to the annual cycle mentioned above, we first examine the annual cycles of SST, \(\tau _u\), and HC, and subsequently the correlations between the variables.

Fig. 3
figure 3

Correlation pattern obtained from reanalysis (ERA-Interim) between SST anomaly in EA4 and wind stress anomalies for each month of the year, the first component of the BF \(\lambda _{SST \rightarrow \tau _u}\). Contours indicate significant correlation at the 90 % confidence level

Fig. 4
figure 4

Correlation pattern obtained from reanalysis (ERA-Interim) between \(\tau _u\) anomaly (ERA-Interim) in WA4 and upper ocean HC anomalies (derived from ORAS4) elswhere in the basin for each month of the year, the second component of the BF \(\lambda _{\tau _u \rightarrow HC}\). Contours indicate significant correlation at the 90 % confidence level

Fig. 5
figure 5

Correlation pattern obtained from reanalysis between upper ocean HC anomalies (derived from ORAS4) and SST anomalies (ERA-Interim), correlated pointwise, component three of the BF \(\lambda _{HC \rightarrow SST}\). Contours indicate significant correlation at the 90 % confidence level

The strong annual cycle of SST in the basin is evident from Fig. 1. During boreal winter direct solar radiation warms the ocean with high intensity, temperatures peak in boreal spring when they reach values larger than 300 K \((27~^\circ \hbox {C})\). From May onward the temperature distribution changes: the warm waters retreat north-westwards from the east and a cold tongue forms during boreal summer, spanning the region around the equator and covering the whole eastern to central part of the basin. This cooling is thought to be caused by upwelling of cold water in the eastern part of the basin, which is closely linked to the sudden onset of the West African Monsoon (Mitchell and Wallace 1992). With strengthening of the westward winds and cooling of the SST the BF becomes active (Burls et al. 2011). The cold tongue is most prominent in August. From September onward, upwelling is reduced through weakening of the southerly cross–equatorial winds, and the cold tongue retreats.

Fig. 6
figure 6

Seasonal cycles of the variables relevant to the BF: SST in EA4, \(\tau _u\) in WA4, and upper ocean HC in EA4. Reanalysis data is plotted in black, for reference of the individual models see Table 1

Fig. 7
figure 7

Seasonal cycles of variance of SST in EA4 \([\hbox {K}^2]\), \(\tau _u\) in WA4 \([\hbox {N}^2\hbox {m}]\), and upper ocean HC in EA4 \([\frac{kJ^2}{m}\cdot 10^9]\). Reanalysis values plotted in black

The cold tongue region displays high SST variability on both annual and interannual time scales. Correlating monthly SST time series of the chosen cold tongue index EA4 to the monthly time series of SST over the whole basin yields the pattern of the interannual mode of SST variability, the zonal mode. Performing regression analysis in the same manner delivers information about the strength of this mode of variability. Both are shown in Fig. 2.

The BF consists of three components: the influence of SST anomalies in the eastern part of the equatorial Atlantic basin on \(\tau _u\) (\(\lambda _{SST \rightarrow \tau _u}\)), the effect of \(\tau _u\) anomalies in the west on the east-equatorial HC (\(\lambda _{\tau _u\rightarrow HC}\)), and the local effect of the HC anomalies on overlying SSTs (\(\lambda _{HC \rightarrow SST}\)) in the cold tongue region. The connection between the two variables of the three components can be illustrated by correlating the anomaly time series to another according to Eq. 2, where \(\sigma _a\) and \(\sigma _b\) are the standard deviations of variable a and b, respectively, and primes denotes anomalies.

$$\begin{aligned} r(lat,lon)=\frac{\sum { a'(t)\cdot b'(t,lat,lon)}}{\sigma _a \cdot \sigma _b} \end{aligned}$$
(2)

Appropriate indices are chosen to represent the eastern and the western basin (see Fig. 2 and Sect. 2) for variable a.

Fig. 8
figure 8

Pattern correlation between models and reanalysis (MRA) plotted against the correlation value between the two variables of the respective component of the BF averaged over the area of interest. The red line denotes the reanalysis correlation value in the region, the pink line the multi model average. The three columns show the results for May, June and July individually (April on May, May on June, and June on July in the case of \(\lambda _{\tau _u \rightarrow HC}\))

Fig. 9
figure 9

Pattern correlation vs correlation strength plot, models lagging reanalysis by 1 month. The three colums show the results for \(\hbox {May}_{reanalysis}\) compared to \(\hbox {June}_{model}\), \(\hbox {June}_{reanalysis}\)\(\hbox {July}_{model}\), and \(\hbox {July}_{reanalysis}\)\(\hbox {August}_{model}\), respectively

Fig. 10
figure 10

Model reanalysis agreement plotted against the absolute error of the upper ocean HC in EA4 for May, June and July. This plot shows that the mean state of the HC does not give indication for the model performance regarding \(\lambda _{HC\rightarrow SST}\)

In order to detect the separate components of the BF, we adapt a method similar to Keenlyside and Latif (2007). For each component of the feedback we correlate the spatially averaged timeseries of anomalies in the region of interest of variable a point-wise to the time series of variable b. Each of the interactions is investigated seperately for every calendar month.

For the first component of the BF, \(\lambda _{SST\rightarrow \tau _u}\), we average SST’ over the box EA4 (see Fig. 2). The box was chosen to be slightly larger than reported in previous literature in order to avoid constraining the model output too much. In Fig. 3 the correlations between SST’ and \(\tau _u\)’ are shown for each month of the year. Because the response time of the atmosphere to SST’ is less than a month, we use a zero month lag for \(\lambda _{SST \rightarrow \tau _u}\).

Fig. 11
figure 11

Pattern correlation plotted against absolute error regression strength for May, June, and July. The bigger the error, the smaller is the model–reanalysis agreement. Hence, the stronger the influence of HC anomaly on SST anomalies, the closer the pattern becomes to a pure BF pattern

Fig. 12
figure 12

Left hand side Model ensemble mean response pattern of \(\lambda _{HC \rightarrow SST}\). Right hand side standard deviation of correlation values between the individual models for pattern \(\lambda _{HC \rightarrow SST}\)

Fig. 13
figure 13

Equatorial and \(30~^\circ \hbox {W}\) temperature variance cross sections in the month July of the series, upper 300 m of the TA, obtained from reanalysis. Black lines denote the \(20~^\circ \hbox {C}\) and \(22~^\circ \hbox {C}\) isotherms

A significant correlation associated with the first component of the BF appears in March, broadens in April, and is strongest in the western Atlantic along the equator in May. The equatorial response of \(\tau _u\) to SST’ is extends to a large area north and south of the equator. In June the area of highest correlation departs from the equator, indicating that another mechanism, possibly connected to the migrating ITCZ becomes more important. In July, the area with maximum correlations has moved northward together with the ITCZ.

The second component of the BF, \(\lambda _{\tau _u \rightarrow HC}\), describes the impact of wind stress anomalies in the western part of the basin on the eastern equatorial oceanic HC. \(\tau _u\)’ in the west generate Kelvin waves that propagate eastward along the equatorial waveguide and cause a deepening of the thermocline in the east, thereby increasing the upper ocean HC. Due to the time it takes a Kelvin wave to travel from east to west we introduce a lag of 1 month when correlating the two variables in order to exclude direct influences of \(\tau _u\) on HC via an atmospheric pathway. To identify \(\lambda _{\tau _u \rightarrow HC}\), we average \(\tau _u\) over the box WA4 (see Fig. 2) and correlate the anomalies to the HC’ in a point-wise manner. While testing the sensitivity of the results to the chosen lag, we observe that the highest correlation between \(\tau _u\)’ and HC’ with a BF like pattern occurs at zero month lag in the reanalysis data, although the results are very similar to a 1 month lag. We hypothesise that other mechanisms than the BF, such as Ekman pumping, can be responsible for this, and therefore retain the 1 month lag.

Fig. 14
figure 14

Equatorial temperature variance cross section in the upper 300 m of the TA for each of the models compared in this study, July. Black lines as in Fig. 13

Fig. 15
figure 15

Temperature variance cross section along \(30~^\circ \hbox {W}\) in the upper 300 m of the TA obtained from model output, July. Black lines as in Fig. 13

Fig. 16
figure 16

Annual cycle of equatorial \(\hbox {SST}_{EA4}\)\(\varTheta _o\) correlation, reanalysis data (ORAS4). Contours indicate statistically significant correlation at the 90 % confidence level

The monthly results are shown in Fig. 4. We notice a build up of high correlation centered on the equator in March–April and a further increase in April–May. The correlation weakens in May–June, and fades into a pattern that is reminiscent of a discharge process in June–July. The timing of the correlation and its maximum strength in May coincides with that of \(\lambda _{SST\rightarrow \tau _u}\), suggesting that the two processes are closely connected. In July, the pattern clearly deviates from the expected \(\lambda _{\tau _u \rightarrow HC}\) pattern associated with the BF. This is also the case for \(\lambda _{SST\rightarrow \tau _u}\), and indicates the end of the season when the BF is active.

The third component of the BF describes the influence of oceanic HC on overlying SST’, \(\lambda _{HC \rightarrow SST}\); mainly to be seen in the cold tongue region where the thermocline is shallow and allows cold subsurface waters to mix with the warm surface waters (Fig. 5). This is the least well-defined of the three components, displaying correlations between subsurface anomalies and SST’ in several regions of the basin. Focussion on the equator, the correlation increases as early as April, albeit very weakly at this stage. The signal peaks in May, and persisits throughout June. In the latter month, the signal is still located almost symmetrically along the equator, which is hardly the case for \(\lambda _{SST\rightarrow \tau _u}\). The signal quickly vanishes afterward, consistent with the timing of the other two components which have also begun to cease by that time. There is also a significant correlation signal in November up to January, indicative of the second peak of the zonal mode described by Okumura and Xie (2006). Signs of this can also be found in the other two components of the BF loop shown above.

Fig. 17
figure 17

Model equatorial \(\hbox {SST}_{EA4}\)\(\varTheta _o\) correlation in June. Contours as Fig. 16. The observerd subsurface–surface coupling is simulated by only a few models

Fig. 18
figure 18

As Fig. 17, but for August. The CMIP5 models display an interaction between subsurface- and surface temperature later in the year

From the monthly stratified correlation analysis of reanalysis data we conclude that the BF mechanism indeed exists in the TA, confirming earlier studies. The activity of the three components coincides in May and June, when distinct correlation patterns are visible. Zonal wind stress anomalies originate from SST’ early in the year (March, April). In these months, \(\tau _u\)’ induce oceanic perturbations in the western part of the TA that travel to the east of the basin. There, the ocean stratification reaches a state in May in which the \(\lambda _{HC \rightarrow SST}\) coupling is possible. The resulting SST’ again affect the atmospheric circulation along the equator, impacting \(\tau _u\)’. During these months the feedback loop is active. The equatorial atmosperic response decreases in June, when the atmosphere is influenced by other factors, most importantly the migration of the ITCZ. Remnants of the feedback loop can be observed in July, although other mechasnisms seem to more dominant at this point.

In contrast to other studies about the BF, we show a monthly stratified picture of the BF. It becomes clear that the signal is not always centered on the equator throughout boreal summer. Other processes seem to mix with the BF leading to correlation patterns that deviate from the one found by Keenlyside and Latif (2007). Their study also analyses the seasonality of the BF. Our results are in line with the authors’ finding that \(\lambda _{SST \rightarrow \tau _u}\) is mostly active in boreal spring and early boreal summer, and that \(\lambda _{HC \rightarrow SST}\) is mostly active in late boreal spring and early summer. The regression analysis in Keenlyside and Latif (2007) indicates a coupling until as late as August, which we cannot confirm, as the signal is no longer centered on the equator in our analysis. While the tropical Atlantic atmosphere is still sensitive to SST’ in the east, the BF loop is not completed by an equatorial Kelvin wave. In the reanalysis data set analysed here the oceanic response of HC to \(\tau _u\)’ is also subject to seasonal modulation. The signal on the equator is most prominent in boreal spring until early boreal summer.

We have demonstrated that it is important to look at the monthly stratified picture of the tropical Atlantic when investigating the BF, as the correlation patterns of the three components vary strongly between the months. We continue to examine the CGCMs performance in the months May, June, and July. During May and June the correlation patterns are most distinct, especially for \(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{HC\rightarrow SST}\). In July the correlation strength weakens and the pattern changes, indicating the end of the BF.

4 Comparison between reanalysis and model output

4.1 Annual cycle

In the former section we have established that a BF exists in the equatorial Atlantic ocean, confirming other studies (Keenlyside and Latif 2007; Okumura and Xie 2006). It is strongest in the early boreal summer months May and June and fades in July. In the following we will compare the model output to reanalysis data in order to investigate whether the models capture the BF. Because of the strong seasonality and the phaselocking to the annual cycle we will first compare the annual cycle of the reanalysis to those of the models, as shown in Fig. 6. The three variables that are related to the BF, SST, \(\tau _u\), and HC, are examined in the regions of interest for the BF, i. e., EA4 for SST and HC, and WA4 for \(\tau _u\).

In all three panels of Fig. 6 there are obvious deviations of the model output from reanalysis, and there is a large spread between the models. For SST (left panel of Fig. 6) the models simulate the annual cycle reasonably well during boreal winter and early spring. The model ensemble mean almost coincides with reanalysis from November to April. This implies an overall warm bias in the CGCMs, as pre-industrial control simulations are used. From May onward the reanalysis SST cools strongly, generating a slope that none of the models reproduces. In the following months the models’ SST lag behind reanalysis SST, reaching their minimum SST one to 2 months later than in the reanalysis data. Another feature in the CGCMs’ SST annual cycle is their relatively weak amplitude, especially during boreal summer and fall. The weak cooling of eastern equatorial SST is associated with the models’ inability to reproduce the boreal summer cold tongue.

The annual cycle of \(\tau _u\) in the western part of the basin displays similar shortcomings (center panel of Fig. 6). Overall, \(\tau _u\) is markedly underestimated by almost all CGCMs under investigation. None of the examined CMIP5 models simulate the weakening of the winds in September and the following strengthening in November associated with the second zonal mode (second Niño, Okumura and Xie 2006). In the reanalysis \(\tau _u\) reaches its minimum in April, while several models simulate \(\tau _u\) at its weakest in May, some even in June. \(\tau _u\) in the reanalysis reaches its first local maximum in August, while most of the models simulate strengthening until September. This is in agreement with the lag we also find for the minimum of \(\tau _u\), but is also connected to the missing weakening of \(\tau _u\) in September as described above.

The oceanic HC displays even larger differences in seasonality with respect to the reanalysis than \(\tau _u\) and SST. In addition to the large spread between the models, the absolute value of the HC is up to twice as large as that of the ocean reanalysis and there is a distinct lag in annual cycle of the CGCMs with respect to the reanalysis. The seasonal minimum of the HC is clearly misplaced by up to 3 months and the models fail to capture the details of the biannual characteristic. Note that the similarities in the annual cycle of HC shown here and the annual cycle of \(\hbox {Z}_{20}\) as shown by Richter [Fig. 6 in Richter and Xie (2008)] demonstrate that HC is a valid measure for the oceanic subsurface condition as indicated by \(\hbox {Z}_{20}\).

In this section, large errors in the annual cycle for each of the three variables that are involved in the BF have been shown, confirming earlier studies, e. g., Richter and Xie (2008), Richter et al. (2014). We continue to show the annual cycle of variances for each of the three variables.

In Fig. 7 shortcomings of the models are distinctly visible. For example, for the variance of \(\tau _u\) (center panel) we find that the reanalysis data displays a plateau of high variance from April to May, when the winds strengthen. This is also the time when the atmosphere at the equator reacts most strongly to SST’ in the east. Most of the models display a single peak which occurs about one to 2 months later in the year than observed. This lagged relationship also appears in the annual cycles of variance for SST and HC. Apart from the lag, there are also large errors in the amplitude of variance in CGCMs compared to the variance in reanalysis. These errors are both positive and negative.

For the reanalysis data, the time frame when the BF is active coincides with the maximum variance of the relevant variables, especially in case of \(\tau _u\). SST variance steeply increases in May and peaks in June, when the BF loop is closed. We have found biases in the CGCMs’ annual cycles, but does this also imply errors in the simulated BF? This could be expected fom the result that the BF is only active in a particular season in the reanalysis: it needs a specific background state to operate. We investigate this question by correlating the three variable pairs obtained from model output in the same manner as done for the reanalysis, and subsequently compare the spatial response patterns to the pattern obtained from reanalysis by pattern correlation analysis.

4.2 Pattern correlation analysis

To investigate whether a model reproduces the BF we correlate the time series of variable a and b with each other to obtain their coupling pattern as shown for the reanalysis (see Eq. 2). Subsequently, we perform pattern correlation between the model and the reanalysis correlation field, according to Eq. 3. Here, \(r_{ra}\) and \(r_{m}\) are the correlation coefficient fields from correlating the variables a and b from the reanalysis and model output, respectively. \(\sigma _{r_{ra}} \cdot \sigma _{r_{m}}\) is the product of the standard deviations of the reanalysis and model correlation fields. Prime denotes anomalies as in Eq. 2.

$$\begin{aligned} MRA = \frac{\sum {r_{ra}'(lat,lon) \cdot r_{m}'(lat,lon)}}{\sigma _{r_{ra}} \cdot \sigma _{r_{m}}} \end{aligned}$$
(3)

A model that reproduces the coupling pattern yields a high pattern correlation value. We call this value “model-reanalysis agreement” (MRA), it is plotted on the y-axis of the multi-model comparison plots (Figs. 8, 9). Note that the MRA depends only on the spatial pattern of the response, and not on its correlation strength. Hence much weaker ab correlation values in the model can still score highly on the MRA axis. To take the correlation between the variable pair into account, we plot the correlation value spatially averaged over the region of interest (i. e., WA4 for \(\lambda _{SST\rightarrow \tau _u}\), and EA4 for \(\lambda _{\tau _u\rightarrow HC}\) and \(\lambda _{HC\rightarrow SST}\)) on the x-axis. The original correlation value obtained from the reanalysis is marked as a red line. The ensemble mean correlation strength between the two variables under investigation is marked by a pink line.

For the first component of the feedback, \(\lambda _{SST\rightarrow \tau _u}\) (first row), we see an improvement in model performance of \(\lambda _{SST \rightarrow \tau _u}\) from May to June (Fig. 8). In May (first row, left panel) the ensemble mean does not simulate the BF, only a few models are located in the upper right of the plot, which means that those CGCMs show some agreement with the reanalysis in both pattern and \(\lambda _{SST \rightarrow \tau _u}\) correlation strength. However, none of the models reaches the inter-variable correlation strength obtained from reanalysis. In June (first row, center panel) the overall pattern correlation value increases compared to May, along with an increase in the correlation strength of \(\lambda _{SST\rightarrow \tau _u}\). The ensemble mean correlation strength moves closer to the observed one and more models display large pattern agreement with the reanalysis. In July (upper right) the picture is improved again, now most of the models under investigation agree with the reanalysis. Note that in this month the reanalysis pattern of \(\lambda _{SST \rightarrow \tau _u}\) does not consist of the pure BF pattern anymore, because the atmosperic response to SST’ has moved towards the north.

The second component of the feedback, \(\lambda _{\tau _u \rightarrow HC}\), is simulated better by the CGCMs. Throughout the 3 months in boreal summer we investigate, high model–reanalysis agreement (MRA) is obtained. The correlation between \(\tau _u\) and HC is only slightly underestimated and especially in June nearly all models simulate this component of the feedback loop correctly. In July we observe higher \(\tau _u \rightarrow HC\) correlations in the EA4 region than in the reanalysis.

The influence of HC’ on overlying SST’ at the equator (\(\lambda _{HC\rightarrow SST}\)) is simulated less accurately. The spatial pattern of this BF component is poorly represented throughout the whole season we investigate. Similarly to \(\lambda _{SST \rightarrow \tau _u}\), the MRA scores lowest in May and improves slightly with increasing \(\lambda _{HC \rightarrow SST}\) correlation strength in June. Only very few models reach high MRA, even though the correlation strength between HC’ and SST’ almost reach reanalysis values. As for the other components, the multi-model mean variable correlation strength exceeds the reanalysis value in July, but in this case MRA has not increased significantly. The low MRA shows that, while the CGCMs display an interaction almost as strong as in the reanalysis in the eastern part of the basin, the response pattern is different from the one in the reanalysis. Strong correlations between HC’ and SST’ are simulated in locations that are absent in the reanalysis and vice versa. This could be due to the fact that we compare basinwide patterns. However, when we restrict the area we use for the pattern correlation to the eastern basin the model–reanalysis agreement is only marginally increased. This indicates that the low MRA also arises from the discrepancy in the spatial representation of \(\lambda _{HC \rightarrow SST}\) in the cold tongue region. The same holds for the \(HC \rightarrow SST\) correlation strength on the x-axis if we choose a smaller index region (\(2^\circ \hbox {N}\)\(2^\circ \hbox {S}\), \(20^\circ \hbox {W}\)\(10^\circ \hbox {E}\), not shown). Our conclusion is therefore that the third component of the BF is the one represented worst with respect to the spatial structure out of the three components of the BF.

In the following we investigate whether the incorrectly simulated annual cycle is responsible for errors of the simulated BF. Most models display lags of variance by 1–2 months, but the variance does not uniformly lag for all variables. We compare the model correlation patterns to the ones from reanalysis, 1 and 2 months earlier each, that is \(\hbox {May}_{Reanalysis}\) to \(\hbox {June}_{CGCM}\), and so on. Regarding the size of the basin and the time scale on which the BF is active, a longer lag time between reanalysis and CGCM would not be physical. This is determined by the atmospheric and upper ocean adjustment time scales.

By introducing a lag of 1 month, the MRA is drastically improved for \(\lambda _{SST \rightarrow \tau _u}\) in all months (cf. Fig. 9). The majority of the models is now situated in the upper right corner of the plot revealing that the pattern response of \({\tau _u}'\) to SST’ is similar to the reanalysis response pattern if this lag is taken into account. In August, the CGCMs produce a pattern that is close to the reanalysis pattern in July. This indicates that the atmosperic response move from pure BF to another process, such as following the migration of the ITCZ, as has been suggested earlier.

For \(\lambda _{\tau _u \rightarrow HC}\) similar behaviour is observed, even though the agreement between model output and reanalysis is reasonable without a lag as well. This component of the BF is similar to the one in reanalysis in almost all of the models.

The simulation of the third component of the BF, \(\lambda _{HC\rightarrow SST}\), is also improved by introducing a lag of 1 month between reanalysis and CGCMs, but here the model output does not reach similar agreement with the reanalysis as for \(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\). For all three components of the BF a lag of 2 months yields similar results as a lag of 1 month (not shown).

From the results of this section we conclude that \(\lambda _{SST\rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\) of the BF are reasonably well represented by the models, with \(\lambda _{\tau _u \rightarrow HC}\) simulated better than \(\lambda _{SST \rightarrow \tau _u}\). The pattern of the correlation between HC’ and overlying SST’, \(\lambda _{HC \rightarrow SST}\), is not represented as well. One hypothesis is that the systematic error in the mean state of the thermocline Richter et al. (2014), and hence the HC, is responsible. If this is the case we expect to see a connection between the mean state error of HC in EA4 and the agreement between models and reanalysis \(\lambda _{HC \rightarrow SST}\) in such a way that a model which displays a reduced mean state error would score a higher MRA value. However, Fig. 10 reveals that no significant linear connection can be found. The misrepresentation of \(\lambda _{HC \rightarrow SST}\) can, hence, not simply be explained by the erroneous mean state (of HC) in the TA basin. Neither is there a connection between the mean state of SST in EA4 and the model performance with respect to simulating \(\lambda _{HC \rightarrow SST}\) (not shown). Note that the missing connection might indicate that the CGCMs’ mean state is too far removed from the reanalysis to obtain a linear relation.

4.3 Regression values

Linear regression analysis (not shown) reveals that for \(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\), which are reasonably well simulated in the models, there is no consistent connection between regression strength and MRA. The regression strength is under- and overestimated by the individual models, but without significant connection to the spatial pattern obtained from the correlation analysis. This is consistent with the results of the pattern correlation analysis. However, for \(\lambda _{HC \rightarrow SST}\) we notice a relationship between MRA and regression strength. This relationship becomes more obvious when the MRA value is plotted against the absolute error in regression values, as done in Fig. 11. In general, the models underestimate the regression strength of \(\lambda _{HC \rightarrow SST}\), especially in May and June. Most models with large absolute errors in regression strength strongly underestimate the actual influence of HC’ on SST’. We deduce that a model that shows a larger influence of HC’ on overlying SST’ is also likely to have a response pattern more similar to the actual response as obtained from reanalysis. In other words, the larger the influence of HC’ on SST’, the closer the \(\lambda _{HC \rightarrow SST}\) pattern is to reanalysis, i. e., the better the BF is simulated.

4.4 Subsurface structure

As we have demonstrated above, the influence of HC’ on overlying SST’ for the individual models cannot readily be linked to a typical BF pattern in CMIP5 models. In reanalysis, on the other hand, there is a distinct pattern on the equator in May and June. When averaging the response pattern of all models we obtain a pattern (Fig. 12) that is similar to the one obtained from reanalysis (cf. Fig. 5). The correlation strength of \(\lambda _{HC \rightarrow SST}\) is, however, much lower than in reanalysis, due to different correlation patterns in the different CGCMs. Note that the standard deviation is large all over the basin, which indicates that individual models differ greatly in their representation of \(\lambda _{HC \rightarrow SST}\).

We investigate the origin of the inaccurate spatial representation of \(\lambda _{HC \rightarrow SST}\) by analysing the subsurface temperatures of the equatorial Atlantic. Two temperature variance cross sections (along the equator and \(30^\circ \hbox {W}\)) are shown for the reanalysis and for each of the models under investigation (Figs. 1314, and 15). In regions where the temperature variance extends to the surface, coupling between subsurface (and hence HC) anomalies and SST’ is present. In the reanalysis it occurs predominantly at the equator around \(10~^\circ \hbox {W}\) in the EA4 region (Fig. 13). This is vastly misrepresented by the CGCMs. The model output shows large deviations from the reanalysis; coupling from subsurface to surface takes place at different locations, if at all, throughout the months when the BF is active (only July shown). This is also illustrated by the longitude cross section in Fig. 15. Here, off-equatorial temperature anomalies couple with the surface much more than observed in reanalysis.

In the Pacific, Zelle et al. (2004) have identified two mechanisms through which SST’ are influenced by HC’, the upwelling and the wind coupling pathway. In the former, the subsurface anomaly propagates vertically through mixing and upwelling. In the latter, it propagates along the equatorial thermocline, causes an SST anomaly on the eastern side of the basin (through the upwelling pathway), which causes a \(\tau _u\) anomaly over the basin. This \(\tau _u\)’ influences the SST’ overlying the original subsurface anomaly. The wind coupling pathway is considered to be active in the western part of the Pacific, where the thermocline reaches depths of more than 120 m. Where this pathway is active, the thermocline–SST’ correlation peaks at lags of 4–8 months. Considering the size of the Atlantic ocean compared to the Pacific, the lag would be considerably shorter if this pathway was active. In the reanalysis the \(\lambda _{HC \rightarrow SST}\) correlation peaks at zero month lag, and the thermocline in the eastern part of the Atlantic is very shallow (ca. 50 m). This indicates that the direct pathway is dominant in the eastern part of the TA, where it forms the third component of the BF, \(\lambda _{HC\rightarrow SST}\). Some models display growing correlation between HC’ and SST’ with increasing lag in the western part of the basin, suggesting the presence of the wind coupling pathway (not shown). However, this feature is absent in the reanalysis.

We can obtain more insight into the oceanic origin of the SST’ by correlating SST’s in the cold tongue region with subsurface temperature anomalies \(\varTheta _o\)’. In the reanalysis a clear picture of the BF and its seasonality is found (Fig. 16). In February and March, SST’s in the cold tongue are mostly correlated with SST’ in the well mixed surface layer stretched along the equator. In April, when both \(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\) become active, correlation with the underlying ocean temperature \(\varTheta _o\) becomes significant. The coupling strengthens and reaches greater depths in May and June, and retreats in July. In August it has almost disappeared. The correlation between \(\varTheta _o\)’ and SST’ is significant only in the months when the BF is active. The second Niño in boral winter can also be seen. It is noteworthy though, that here the correlation between \(\hbox {SST}_{EA4}\)’ and \(\varTheta _o\)’ is more confined to the coast, while in boreal summer it is present in the whole eastern TA. We repeated the previous analysis for each of the individual CGCMs for June and August. In June (Fig. 17), very few CGCMs display a \(\hbox {SST}_{EA4}\)’ – \(\varTheta _o\)’ correlation pattern similar to the one in reanalysis. The coupling between subsurface and surface anomalies does not take place adequately. In August (Fig. 18), most models simulate the coupling, although some models display spurious subsurface patterns that are not connected with the surface (e. g. panels five, twelve, and fourteen). These might be associated with meridional current patterns in the CGCMs that are not present in the reanalysis. Whether this coupling is associated with the BF is questionable. In August the BF loop is no longer closed. The \(\hbox {SST}_{EA4}\)’ is not centered on the equator anymore, suggesting that other mechasnisms, such as the WES feedback, dominate. Note that we compared the correlation patterns of \(\lambda _{SST \rightarrow \tau _u}\) with a lag of 1 month. August in the CGCMs is similar to July in the reanalysis, when the pattern has departed from the equator and the BF coupling is weakening. Some models have already established an \(\hbox {SST}_{EA4}\)’ – \(\varTheta _o\)’ correlation in July. When that is the case, the models most likely simulate the BF completely. The lack of high MRA values can then be explained by (misrepresented) off-equatorial variability.

5 Conclusion

In this study we investigate the BF in the tropical Atlantic ocean from both reanalysis data and CMIP5 model output. The reanalysis data clearly shows the presence of a zonal mode which is driven by the BF, confirming previous studies. The monthly stratified analysis shows that the BF is not equally active throughout the year, it is strongest in May and June and begins to fade again in July.

\(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\) are reasonably well simulated by the models based on correlation strength and spatial pattern. While \(\tau _u\)’ in the western part of the basin reacts to SST’\(_{EA4}\) within a month, we introduce a lag between \(\tau _u\)’ and HC’ in order to exclude direct influences of \(\tau _u\)’ on SST’ in the eastern part of the basin. By evaluating the eastern equatorial HC’ response to \(\tau _u\)\(_{WA4}\) with a lag of 1 month we allow for a Kelvin wave to propagate from west to east.

The agreement of \(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\) between reanalysis and model output improves further when introducing a lag of 1 month between them, correcting for the lag of the annual cycle (of variance) in the CGCMs. We conclude that the physical mechanisms behind \(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\) are simulated well by the models, but that the lag in the annual cycle introduces a lag in the timing of the BF.

\(\lambda _{HC \rightarrow SST}\), on the other hand, is generally not well represented by the models. The patterns of the responses of CGCM output and reanalysis agree to less than 25 % in almost all cases. We note that there is a significant correlation between HC’ and SST’ in the eastern part of the basin, but its location and pattern are significantly different from the one detected in reanalysis. No significant linear correlation between mean state error of the HC in the eastern part of the basin and the model–reanalysis agreement could be found, so the reason for the error in the third component of the loop could not conclusively be linked to the mean state bias of SST and HC. This does not exclude the possibility that individual models perform better when their mean state is corrected. In fact, this was recently shown to be the case for the Kiel Climate Model, which does not belong to the CGCM ensemble investigated here (Ding et al. 2015). We speculate that differences between the mean state of the CGCMs and the reanalysis mean state may be too large to display a linear relationship between the mean state error and the performance with respect to the BF.

Examining the subsurface temperature variance provides insight into the origin of the largely varying \(\lambda _{HC \rightarrow SST}\) response fields. This component of the BF describes an interaction between subsurface and surface anomalies. In the equatorial cross section of temperature variance in the reanalysis there is a clear connection between subsurface and surface anomalies. Here, mixing takes place and SST is influenced by subsurface temperatures. None of the models is able to reproduce this structure. Hence, while there is some interaction indicated by the \(\lambda _{HC \rightarrow SST}\) correlation, the location of this interaction is different from reanalysis. This is also demonstrated by longitude–depth cross sections that show \(SST_{EA4}\)’–\(\varTheta _o\)’ correlations. The seasonality displayed by reanalysis is modelled by hardly any of the CGCMs. When the correlation decreases in reanalysis, it increases in the models and even persisits throughout the year (not shown).

The timing between \(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\) and the subsurface–surface coupling is misplaced in such a way that the full BF loop is most likely not closed for most of the models. \(\lambda _{HC \rightarrow SST}\) in the CGCMs becomes active too late in the year. The center of atmospheric sensitivity to \(\hbox {SST}_{EA4}\)’ has moved north of the equator, indicating that other processes play a more important role than the BF along the equator. The lack of seasonality displayed in \(SST_{EA4}\)’–\(\varTheta _o\)’ in the CGCMs hints on other processes inducing the \(SST_{EA4}\)’–\(\varTheta _o\)’ coupling, rather than being initiated by the BF.

Based on the results of this study, we conclude that the erroneous vertical oceanic stratification influences the interannual variability in the region by hindering \(\lambda _{HC \rightarrow SST}\) from being represented correctly.

Even though the models display an SST annual cycle much weaker than in reanalysis, \(\lambda _{SST \rightarrow \tau _u}\) and \(\lambda _{\tau _u \rightarrow HC}\) are reasonably well simulated, albeit with an error in timing because of the lags in seasonality. This disproves our hypothesis that the functioning of the BF depends on the correctly simulated annual cycle. The mean state of the TA is misrepresented in the CGCMs, but the influence of perturbations of variable a (SST, \(\tau _u\)) on variable b (\(\tau _u\), HC) is still present. This is in line with earlier findings that even though the mean state of the TA is heavily biased, interannual variability can still be reasonably represented Richter et al. (2014). However, biases in the subsurface temperature variance affect the spatial pattern of \(\lambda _{HC\rightarrow SST}\). \(\lambda _{HC\rightarrow SST}\) explains a large part of SST variability in the region, and exactly this part of the feedback is not simulated well.

The errors in the simulated subsurface thermal structure can be due to both ocean mixing and/or arise from shortcomings of the atmospheric forcing which sets the structure of the ventilated thermocline. This should be the focus of further investigation and a starting point for work on model improvement in order to improve the HC–SST relationship in the eastern tropical Atlantic.