Improvements and persistent biases in the southeast tropical Atlantic in CMIP models

State-of-the-art climate models simulate warmer than observed sea surface temperatures (SST) in eastern boundary upwelling systems (EBUS), generating biases with profound implications for the simulation of present-day climate and its future projections. Amongst all EBUS, the bias is largest in the southeastern tropical Atlantic (SETA). Here, we provide a comprehensive evaluation of the performance in the SETA of the Coupled Model Intercomparison Project phase 6 (CMIP6), including fine resolution (HighResMIP) and ocean-forced (OMIP) models. We show that biases in the SETA remain large in CMIP6 models but are reduced in HighResMIP, with OMIP models giving the best performance. The analysis suggests that, once local forcing errors have been reduced, the major source of the SETA biases lies in the equatorial Atlantic. This study shows that finer model resolution has helped reduce the local origin of the SETA SST bias but further developments of model physics schemes will be required to make progress.


INTRODUCTION
Eastern boundary upwelling systems (EBUS) are regions located on the eastern boundaries of subtropical gyres and characterized by coastal upwelling. There are 4 major EBUS: the California Current System (CCS), the Canary Current System (CaCS), the Humboldt Current system (HCS) and the Benguela Current System (BCS). These regions occupy merely 0.1% of the global ocean but with profound climatic 1 and ecological importance 2 . For example, EBUS contribute to the global carbon cycle and their role as either a source or sink of CO 2 is still uncertain 3 , they are the most biologically productive ecosystems in the ocean 4 , and are subject to basin-scale variability 5 and change 6,7 The large upwelling regions are dominated by two processes: offshore Ekman transport associated with equatorward alongshore wind-stress and Ekman pumping driven by positive (negative) wind stress curl in the Northern (Southern) Hemisphere. The two processes cause nutrient rich and cold subsurface water to rise to the surface influencing coastal sea surface temperature (SST). In coupled climate models, and in nature, SST is one of the most important ocean variables as it drives and controls many ocean-atmosphere feedbacks and modes of variability. However, most state-of-the-art climate models are characterized by large SST biases in EBUS regions, with their origin attributed to deficiencies in the atmospheric [8][9][10][11][12] , oceanic [13][14][15] , or both components of climate models [16][17][18][19] . SST biases, and in particular those in EBUS regions, are a long-lasting problem in coupled models 1,14 , and a major source of uncertainty in seasonal-to-decadal predictions and regional-to-global projections of climate change 20 . In fact, artificially correcting for those biases has shown to improve both the mean state and variability 1,21 .
Amongst all EBUS, the SST bias is warmest in the BCS, which extends from about 15 ∘ S (north of the Cunene upwelling cell) to 30 ∘ S (the Agulhas cell) along the coast of southern Africa 22 . Further, the BCS is the most peculiar EBUS as it is bounded by two opposing currents: the poleward Angola current and the equatorward Benguela current 23 . The former draws from the equatorial undercurrent and flows poleward. At about 16 ∘ S, the Angola current meets the cold northward Benguela current creating a sharp thermal front known as the Angola-Benguela Front (ABF) [24][25][26] (Fig. 1f). Such a dynamical feature is not observed in the southeastern Pacific, where downwind currents merge directly with equatorial waters, and thus the BCS is the most complex of the four boundary systems to model 18 . The SST bias in the southeastern tropical Atlantic (SETA; 10 ∘ S-30 ∘ S, 5 ∘ E-20 ∘ E), home to the BCS, has been mainly attributed to two causes: an incorrect representation of alongshore meridional wind stress 9,10,27 leading to weak coastal upwelling, and a poleward overshooting of the ABF together with a deep and diffusive thermocline along the eastern side of the Atlantic 14,15,21 creating a subsurface temperature bias that is carried south by the Angola current. Indeed, an improved representation of alongshore wind stress has been shown to play a crucial role in reducing the coastal bias, but with no significant change to the large scale bias in the SETA 9 . The SETA also suffers from a positive shortwave radiation bias in most climate models due to the lack of shallow stratocumulus clouds. Reducing this error has been shown to lead to a reduction in SST bias 8 but often with a minor impact 28 , supporting the hypothesis that ocean mechanisms are far more important 1,14 . Quantifying the role of atmospheric model biases and coupled feedbacks could be achieved through the analysis of net surface flux biases within the framework of equilibrium mixed layer heat budget theory 29 , however a large source of SETA SST bias has previously been identified within the equatorial Atlantic, suggesting coastal Kelvin waves 30 or horizontal advection 13,15,21 as possible mechanisms related to ocean dynamics.
Given the importance of the SETA warm SST bias in the accurate simulation of Atlantic modes of variability and global climate projections, it is of great importance to quantify these biases and sources of errors in state-of-the-art climate models. Phase 6 of the Coupled Model Intercomparison Project (CMIP6) has recently been released with different model inter-comparison projects (MIPs) 31 . It is therefore timely to evaluate the performance of CMIP6 climate models, assessing the role of increased resolution and coupling, with respect to their predecessors (CMIP5) and available observations in the SETA.

SST biases
We use monthly data from the Historical scenario and selected 15 models that participated in CMIP5, 17 models in CMIP6, 15 models in HighResMIP, 8 models in OMIP1 and 10 models in OMIP2 (Supplementary Tables 1-5). Only models providing all necessary variables for our analysis were considered in this study and we cover the period 1985-2004 in order to overlap with the optimally interpolated (OI) Reynolds SST observational data set (OISST; see "Methods"). Compared to CMIP5 models, the CMIP6 multi-model mean (MMM) is consistently colder in all EBUS (Fig. 1a, b). The maximum reduction in SST bias is observed in the HCS and in the BCS (see also Supplementary Fig. 1). In CMIP6, not all models agree on the representation of the HCS and the sign of the SST bias changes depending on the model. Out of 17 models, 9 have negative SST biases in the HCS, especially off the coast of southern Chile. This change in sign is responsible for a significant bias reduction in the southeast Pacific in CMIP6, and therefore explains the large improvement in SST bias in the region between CMIP6 and CMIP5 MMMs (Fig. 1a, b). On the contrary, all models have a positive SST bias in the SETA and, as in the case of CMIP5, the warmest SST bias in the CMIP6 MMM is still found in this region (Fig. 1b). With a maximum value of 6.45 ∘ C for CMIP6 and 7.6 ∘ C for CMIP5 MMMs, CMIP6 has only slightly improved the SETA SST bias, and the spatially averaged SST bias in the region is reduced by0 .45 ∘ C, with the largest signal located within 300 km from the coast as in CMIP5 1,14 . The CMIP5 and CMIP6 multi-model ensembles are made of different models (see Supplementary  Tables 1 and 2) and this might introduce a degree of uncertainty in our MMM comparison. However, if we select a subset of eight modeling centers that are present in both CMIP5 and CMIP6 ensembles, we find an even smaller improvement in the spatially averaged SST bias (~0.27 ∘ C).
Next, we consider the finer resolution models included in HighResMIP. As in previous MMMs, all HighResMIP models agree on a warm bias within the SETA (Fig. 1c), but with a much reduced maximum value of~4.9 ∘ C and with a spatial mean bias of~1.8 ∘ C (55% less than in CMIP5 MMM). Finally, the oceanforced OMIP1 and OMIP2 MMMs show a further reduction in SST bias, with maximum values of~3.9 ∘ C and~2.9 ∘ C, respectively (Fig. 1d, e). Consistent with previous studies, an improvement in atmospheric model resolution is responsible for reducing the SST coastal bias and the reduction is mainly achieved through an improvement of the atmospheric coastal jet representation in the SETA 9 . This result is confirmed by the ocean-only simulations as OMIP1, forced with the coarser CORE atmospheric state, is warmer in the BCS than OMIP2, forced with the finer JRA-55 atmospheric reanalysis (see Methods), especially at Cape Frio (17 ∘ S-18 ∘ S) and at L€ uderitz (26 ∘ S-27 ∘ S) ( Supplementary Fig. S1d). Also, the coastal bias produced by the OMIP MMMs is further reduced from HighResMIP, suggesting coupled feedbacks tend to exacerbate the bias 32 . We note that the SST bias within the equatorial Atlantic does not seem to improve significantly in any of the MMMs, particularly on the eastern side of the basin (see also Supplementary Fig. 3), suggesting that increased resolution in both the ocean and atmosphere is not sufficient in alleviating large scale biases both at the surface and in the interior, an aspect we will explore later.
The BCS upwelling is related to the strong southerly wind stress whose alongshore magnitude is highly variable. In austral summer, there is one maximum located at L€ uderitz and then a second maximum appears at Cape Frio in autumn and spring 33,34 . South of~30 ∘ S, upwelling is highly seasonal and less intense. It is therefore important to assess not only the time-mean SST bias but also its seasonality. Figure 2a shows the climatological annual cycle of the SETA SST bias, averaged over a 2 ∘ wide band along the coast, which can be directly related to coastal upwelling and alongshore wind stress in the MMMs. In both CMIP5 and CMIP6, the bias is maximum in austral autumn and winter, with some reduction in the median and in the maximum value in CMIP6 (the top extent of the whiskers). On the contrary, the maximum bias in HighResMIP MMM is observed in austral spring, the season with maximum bias improvement in CMIP6 and minimum bias in both CMIP5 and CMIP6 MMM, but with a significant reduction throughout the year. Both OMIP MMMs are colder than coupled MIPs and with a maximum bias during austral autumn and winter as in CMIP5 and CMIP6. The annual coastal SST bias (rightmost bars in Fig. 2a) is in line with the results shown in Fig. 1, with a gradual improvement from coarser to finer coupled models or atmospheric forcing. OMIP2 models are forced with the JRA55-do atmospheric dataset, which has double the temporal frequency and more than 3 times the spatial resolution of the CORE dataset used to force OMIP1 ("Methods").

Wind stress biases
Wind stress affects coastal SST through two processes: offshore Ekman transport and Ekman pumping 35 . The two processes are associated with the rising of cold subsurface waters to compensate for the divergent transport. Due to the coarse grid of most climate ocean models, offshore transport is underestimated due to unresolved eddy transport exacerbating coastal SST biases 18 . The alongshore meridional wind stress (taken at the grid cells closest to the coast) remains weaker than observations in the CMIP6 MMM (Fig. 2b). However, the magnitude of the alongshore wind stress HighResMIP is much stronger than in both CMIP5 and CMIP6 MMMs at both L€ uderitz (~26 ∘ S) and Cape Frio (~17 ∘ S) and approaching observational values. The strengthening of alongshore wind stress in HighResMIP can be attributed to the improved representation of orography situated along the Namibian coast due to the increased atmospheric model resolution 9,11 . A larger meridional wind stress and associated offshore transport triggers stronger upwelling 10,36 and consequently reduces the SST bias, as observed in Fig. 2.
Other than intensifying the alongshore wind stress, increased atmospheric resolution also brings the core of the coastal jet, the maximum wind stress, closer to the coast. Figure 3a-d shows the  location of the coastal jet core in the different MMMs and reanalysis, as well as the difference in meridional wind stress among MMMs. Away from L€ uderitz and Cape Frio, the meridional wind stress close to the coast is more intense in CMIP5 than in CMIP6 (Fig. 3a), however, just offshore, the CMIP6 MMM coastal jet is more intense, especially in the southern Benguela. In this region, the upwelling is highly seasonal and is separated from the northern Benguela, which exhibits stronger upwelling all year round, by the prominent L€ uderitz upwelling cell acting as a physical boundary 22,37 . Compared to both HighResMIP and JRA-55, CMIP6 MMM has a large negative bias in meridional wind stress along the coast and its core is located too offshore in the main upwelling cells of L€ uderitz and Cape Frio (Fig. 3b, c). Relative to CMIP5, the core remains nearly unchanged except for the minor variations in the two cells mentioned above. Because of its finer resolution, HighResMIP locks the maximum wind stress closer to the coast, almost matching observational estimates (yellow and black lines in Fig. 3d), and its magnitude also presents weaker biases in the SETA, as suggested by Fig. 2b.
If the core of the atmospheric coastal jet is biased offshore, it creates a band of cyclonic (i.e., negative) wind stress curl close to the coast. The wind stress curl is predominantly dependent on the distance between the coast and the core of the meridional wind stress and on the zonal gradient of the meridional wind stress 10 , the former being too large in both CMIP5 and CMIP6 (yellow and black lines in Fig. 3a). The MMMs wind stress curl is shown in Fig. 3e-g along with the one from JRA-55 (Fig. 3h). The wind stress curl is too broad and with an alongshore intensification south of 20 ∘ S in both CMIP5 and CMIP6. HighResMIP produces a more narrow, coastal-intensified, cyclonic wind stress curl, but the two positive wind stress curl cells located offshore L€ uderitz and Cape Frio are poorly reproduced in HighResMIP. A direct consequence of biases in wind stress curl is a misrepresentation of the induced Ekman pumping. CMIP5 and CMIP6 exhibits a too-wide and intense positive (upward) Ekman pumping, particularly in the southern Benguela (Fig. 3i, l), whereas HighResMIP seems to reproduce the main upwelling cells with more fidelity (Fig. 3m, n). A larger than observed simulated Ekman pumping would seem at odds with a large SST bias, given that the larger Ekman-induced vertical flux in CMIP5 and CMIP6 would be expected to bring subsurface colder waters and cool the surface. However, upwelled waters in both MMMs also have large temperature biases, partially contributing to the resulting SST biases in those models.

Biases in coastal currents and the Angola-Benguela Front
The wind stress curl field located north of Cape Frio is associated with the upwind Angola current 38,39 . In observations, the southward Angola current meets the northward Benguela current at about 16 ∘ S creating a sharp thermal front, the Angola-Benguela Front (ABF) [24][25][26] . The ABF latitudinal position is said to be controlled by the two opposing currents 38 : the upwind Angola current fueled by the wind stress curl and the Benguela current that is formed within a Rossby radius to the coast and driven by the pressure gradient generated by the wind-driven coastal divergence 14,15 . Hence, the wind stress, and associated wind stress curl, not only affects coastal upwelling but also controls the dynamics of the coastal currents. Supplementary Fig. 2 shows the meridional coastal currents along the African coast (from the equator to 30 ∘ S) averaged within 2 ∘ off the coast in all MMMs. The thick line denotes zero velocity, indicating the latitude where the two meridional currents meet. The latitude where this line intersects the surface is a measure of the position of the ABF. The intensity of the upwind Angola current core is similar in both CMIP5 and CMIP6 and located too far south (~17-18 ∘ S) because the excessively large wind stress curl causes the Angola current to overshoot. On the other side of the front, the downwind Benguela current remains fairly weak in CMIP6, although deeper than in CMIP5 due to an improvement in simulated alongshore wind stress (Fig. 3a). The subsurface velocity structure improves considerably in HighResMIP with an equatorward shift of the Angola current core, which shows signatures of a subsurface intrusion below a much deeper and stronger Benguela current ( Supplementary Fig. 2c), as the southward Angola current is forced to become deeper when colliding with the Benguela current in order to conserve potential vorticity 14 . Reduction in overshooting of the Angola current can be attributed to improvements in the wind stress curl (Fig. 3g, h) and the intensified Benguela current to a stronger and coastal-intensified meridional wind stress (Fig. 3c, d). A similar picture is found in OMIP MMMs, with a progressively reduction in Angola current overshooting and stronger and deeper simulated Benguela current. Both an intense Benguela current and a more equatorward Angola current core should then lead to a less spurious southward displacement of the ABF, as the ABF position is tied to the intensity of the two opposing currents 38 .
To quantify the ABF latitudinal position in models and be able to compare with observations, we use the definition of maximum SST gradient within the SETA 11,19 . The observed SST (OISST; see Methods) shows little seasonal variations in ABF, with an equatorward displacement in austral winter coincident with a stronger intrusion from the downwind Benguela current 40 (Fig. 4a). The coupled MIPs all show a strong seasonal dependence of the ABF position and the SST gradient does not have a coherent structure throughout the year. The equatorward propagation of the Benguela current, highlighted by the 20 ∘ C isotherm, is also limited compared to observations, although the HighResMIP MMM is closer in reproducing the observed extent of the seasonal variations (Fig. 4d). As inferred from the analysis of the Angola and Benguela current dynamics, uncoupled OMIP MMMs both show an equatorward shifted maximum SST gradient with a very weak seasonal migration and large Benguela current meridional intrusion during austral winter (Fig. 4e, f). CMIP5 and CMIP6 place the annual-mean ABF at~21 ∘ S (Fig. 4g), which is far too south compared to 15-16 ∘ S from OISST. OMIP models have a mean ABF at~17.5 ∘ S and HighResMIP, although with an incorrect representation of the ABF in austral summer, locks the annual-mean ABF at the same latitude. The misrepresentation of the ABF position in coupled MIPs is well correlated with the coastal SST bias seasonality shown in Fig. 2a. In austral autumn, when the ABF is furthest south in CMIP5 and CMIP6, it attains its northmost position in HighResMIP. The situation is reversed in austral spring, when the ABF in CMIP5 and CMIP6 is equatorward of HighResMIP. This explains why the bias in HighResMIP is maximum in austral spring and winter, and minimum in austral autumn, as opposed to CMIP6 MMM that shows a maximum bias in austral autumn.

Subsurface biases
Next, we look at the subsurface temperature structure in the SETA. We have assessed and quantified the extent of the SST bias in the SETA and how it is related to the incorrect atmospheric forcing, in particular the wind stress curl in the northern part of the SETA and the alongshore meridional winds responsible for the overshooting of the Angola current and a weak Benguela current, respectively. What remains to be understood is the origin of the remaining and significant SST bias in the ocean-forced OMIP and HighResMIP MMMs. If the SST bias cannot be attributed mainly to a deficient atmospheric forcing within the SETA, then it will have to reside in ocean dynamics or be of remote origin, or both.
In the observed temperature vertical structure (WOA18; see "Methods"), the thermocline north of the ABF is sharp and shallow, with the 20 ∘ C isotherm above 50 m (Fig. 5f). Southward of the 20 ∘ C outcropping latitude, observations are characterized by vigorous upwelling portrayed by the doming of isotherms towards the surface. This is mostly around 25 ∘ S where water parcels with temperatures as low as 16 ∘ C reach the surface. Both CMIP5 and CMIP6 present a deep and diffusive thermocline, with the 20 ∘ C isotherm outcropping far too south, related to the location of the ABF, and a similar vertical structure develops in the OMIP MMMs (Fig. 5a, b, d, e). Due to the finer oceanic resolution and improved wind stress, HighResMIP simulates a shallower and sharper thermocline and extensive upwelling in the southern Benguela, as in observations. However, the surface water entering the SETA is as warm as in CMIP5 and CMIP6, suggesting horizontal advection of warm remote waters (Fig. 5c). The existence of a subsurface bias in temperature implies that, even with an intensified surface wind stress, the SST bias in the SETA would still persist due to its remote origin 15,30 . As shown by previous targeted sensitivity model simulations 14,15,21,27 , the temperature bias, located within the top 50 m and above the 20 ∘ C isotherm, is advected southward by the upwind Angola current which submerges under the Benguela current south of the ABF. Within the upwelling region, this subsurface bias is raised to the surface and then spreads offshore, resulting in one of the main contributors to the SST bias south of the ABF 14,15 . This mechanism, shown here in its final stage by MMM time-mean temperature biases, seems at play in all CMIP coupled and uncoupled models (Fig. 6).
Our analysis supports the hypothesis that the failure of state-ofthe-art coupled models to simulate a shallow and sharp thermocline results in a subsurface bias whose origin can be traced in the eastern and central equatorial Atlantic 15 , as coupled models tend to fail to simulate the cold tongue there 8,41,42 . This deficiency is attributed mostly to a wind stress bias in austral autumn that creates a deep (shallow) thermocline in the eastern (western) equatorial Atlantic 27 , preventing the upwelling of cold water and creating an SST bias that is maximum in austral winter. These equatorial Atlantic biases have been suggested to be related to the SETA SST bias through coastal kelvin waves or lateral advection 13,15,43 . Supplementary Fig. 3 shows the SST zonal gradient at the equator, averaged between 2 ∘ N and 2 ∘ S, and we find that the zonal gradient west of 10 ∘ W is reversed in all CMIPs. In the western tropical Atlantic the simulated SST is much colder than observed, and only OMIP MMMs simulate warmer than observed SSTs. The mean zonal gradient in HighResMIP is flatter than in CMIP5 and CMIP6 MMMs. This is because, although the western equatorial Atlantic is warmer in HighResMIP, the magnitude of the cold-tongue remains similar to CMIP5 and CMIP6. A plausible reason is the lack of significant improvement in zonal wind stress biases within the eastern equatorial Atlantic sector ( Supplementary Fig. 4). In the uncoupled ocean simulations, both MMMs are able to simulate the observed zonal gradient and a realistic cold-tongue, but the gradient remains weak and biased warm. Given that the two OMIPs are forced with atmospheric datasets at both coarse and fine resolutions, and that participating ocean models have relatively coarse resolutions (see Supplementary Table 3), the more realistic equatorial Atlantic cold tongue can be attributed to the absence of coupled feedbacks exacerbating biases and the warm bias to deficiencies in vertical mixing and eddy flux parameterizations 44 .
In fact, although wind stress and its curl have largely eliminated their biases in the HighResMIP MMM over the SETA, large biases persist in the equatorial sector, resulting in limited changes in SST there ( Supplementary Fig. 1). The persistent zonal wind stress bias is also responsible for a large negative subsurface temperature bias in the western tropical Atlantic in all coupled models (Fig.  6a-c), which is absent in both OMIP MMMs. However, HighResMIP, OMIP1 and OMIP2 all suffer from an eastern equatorial subsurface temperature bias of similar intensity to the subsurface bias found along the African coast (Fig. 6c-e and h-l). Indeed, full-field initialized seasonal hindcast simulations have shown the horizontal advection of warm anomalies propagating along the equatorial thermocline within the equatorial undercurrent and later along the African coast, eventually upwelling in the southern Benguela region 21,27 . Our MMM analysis seems to support this Fig. 6 Equatorial and SETA subsurface temperature biases. a-e Time-mean subsurface equatorial temperature biases (zonally averaged between 2 ∘ S and 2 ∘ N) in CMIP5, CMIP6, HighResMIP, OMIP1 and OMIP2 MMMs relative to WOA18. f-l Time-mean subsurface temperature biases along the African coast, zonally averaged over a 2 ∘ wide band along the coast, for all MMMs relative to WOA18. Contour interval is 1 ∘ C. The solid and dashed yellow contour indicate the depth of the 20 ∘ C isotherm in the MMMs and WOA18, respectively. mechanism, resembling the final equilibrated state reached after only a few months and ascribed to an equatorial westerly wind bias. If a large fraction of the SETA SST bias in all MMMs is originated remotely then a clear relationship between temperature biases at the equator and SETA should hold. We choose to correlate SST biases in the ATL3 region (5 ∘ S-5 ∘ S, 10 ∘ W-0 ∘ E) 21 and SETA for each model considered in this study (Fig. 7a). CMIP5 shows a weak correlation of 0.14, and this is explained by the large local contribution from wind stress curl and alongshore wind stress in generating the SETA bias. CMIP6 has a 0.52 correlation (significant at the 95% level) between ATL3 and SETA SST biases, due to minor but significant improvements in surface forcing over the SETA region and particularly the southern Benguela (Fig. 3). HighResMIP, OMIP1 and OMIP2 show high correlations of 0.72, 0.75 and 0.81 (all significant at the 95% level). These correlations imply that a more realistic wind stress forcing (HighResMIP) or forcing the ocean with the observed atmospheric state, largely eliminates the local source of SETA SST bias and the remote equatorial origin remains the main contribution. As in the case of SST ATL3 biases, there is a weak correlation between SETA SST biases and equatorial subsurface temperature biases in CMIP5 (0.1) and CMIP6 (0.36) (Fig. 7b). When local sources of SETA SST biases are greatly reduced, however, the correlation amounts to 0.59, 0.73 and 0.78 in HighResMIP, OMIP1 and OMIP2, respectively, suggesting the origin of SST biases in the SETA resides in the equatorial thermocline.

DISCUSSION
In this study, we use models from phase 6 of the Coupled Model Intercomparison Project (CMIP6) to assess their performance in simulating the southeastern tropical Atlantic (SETA), a region that experiences the largest SST bias in climate models. We focus on the ocean structure and its biases, since errors in shortwave radiation have been shown not to play a leading role in setting the warm SST bias 1,8,14 . We show that the CMIP6 ensemble, compared to CMIP5, shows limited improvements in representing the surface and subsurface dynamics in the Benguela Current System (BCS), and the overall bias along the SETA region remains large across all seasons. The SST bias in the BCS and the large SETA region is attributed to two main factors: an incorrect location of the Angola-Benguela Front (ABF) 14,15 that sits too far south and a deep and diffusive thermocline simulated by climate models along the eastern side of the Atlantic 21,27 . Both problems persist in CMIP6. Our analysis of the finer resolution models included in HighResMIP shows a substantial reduction in coastal SST bias and a shallower thermocline, mostly because of an improvement in the simulation of the coastal jet leading to a strengthening in coastal upwelling and also due to a weaker and narrower negative wind stress curl preventing the ABF to overshoot. Further, a similar picture arises when considering uncoupled ocean models (OMIP1 and OMIP2), presenting a further reduction in SST bias within the SETA and a better representation of both the Angola and Benguela currents and associated ABF, particularly for those models forced with the finer atmospheric datasets (OMIP2). However, similar surface and subsurface biases persist in both HighResMIP and OMIP simulations, suggesting a common source independent of refined model resolution.
We find that equatorial temperature biases, which are mainly subsurface rather than SST biases, are not improved in CMIP6 and HighResMIP, and they are only partially alleviated within forced ocean models, suggesting a role for coupled feedbacks in the misrepresentation of the equatorial cold tongue. In fact, refining the resolution in both oceanic and atmospheric components does not help in cooling the eastern equatorial Atlantic, and both biases in SST and zonal wind stress remain large in all multi-model means. The equatorial subsurface temperature bias is shown to be advected along the African coast generating a warm bias at the latitude of the Angola-Benguela front where it is upwelled to the surface. In the case of HighResMIP and OMIP models, where local sources of biases are largely removed, most of the SETA SST bias can be traced to an equatorial origin. We argue that finer resolution in both oceanic and atmospheric components of climate models can reduce the local source of the SETA SST bias but that large biases in the equatorial Atlantic sector still persist and they constitute the main reason for the persevering warm bias. Our multi-model intercomparison results suggest that only through addressing equatorial Atlantic biases and physical processes, with efforts that go beyond increasing model resolution 17,21,45 , could result in eliminating SETA biases and yielding a more reliable representation of the coupled climate system, with implications for modes of variability and future projections in the Atlantic sector.

METHODS
The CMIP models 17 historical simulations from the CMIP6 ensemble 31 Fig. 7 Relation between SETA and equatorial SST biases. a SETA vs ATL3 SST biases in CMIP models. The ATL3 region is defined between 5 ∘ S-5 ∘ S and 10 ∘ W-0 ∘ E. b Same as in (a) but for SETA SST and the equatorial subsurface temperature biases (averaged between 2 ∘ S-2 ∘ N, 10 ∘ W-10 ∘ E, and 50-150 m).  Table 5). HighResMIP is a newly endorsed MIP within the CMIP6 framework and is meant to act as a basis to tackle the effect of horizontal resolution on model biases. It includes models with resolution of at least 50 km in the atmosphere and 0.25 ∘ in the ocean 46 . OMIP has two versions: OMIP1 forced with the CORE.v2 atmospheric state (112 × 94 grid cells) 47 and OMIP2 forced with the JRA-55 reanalysis with a finer spatial resolution (640 × 320 grid cells) 48 . Coupled models with coarse ocean horizontal resolution (greater than 1.5 ∘ ) have not been considered in this study because they are known to underestimate upwelling dynamics. The analysis is based on historical simulations. Historical simulations use prescribed CO 2 atmospheric concentrations to simulate the recent past climate conditions, from the mid-19th century to present day, and hence can be directly compared with available observations. Only the first ensemble member (r1i1p1 & r1i1p1f1 for CMIP5 and CMIP6 respectively) is considered. All models are remapped onto a regular 1 ∘ × 1 ∘ grid before the analysis and the comparison is for the time period spanning the years 1985 to 2004. In the subsurface, variables are interpolated to common depth levels (6, 15, 25, 35, 45, 55,

The observational data sets
For surface temperature, we use the daily optimally interpolated (OI) Reynolds SST (OISST) data set 49 . This data set spans from 1982 to 2019 and is used as the observed SST against which models are evaluated. OISST is based on the second version of the National Oceanic and Atmospheric Administration daily Optimum Interpolation SST with a horizontal resolution of 0.25 ∘50 . For the subsurface, monthly temperatures from the World Ocean Atlas 2018 (WOA18) 51 , with a time span from 1985 to 2004, are used. The data is on a regular 1 ∘ horizontal resolution with 57 vertical levels from the surface to 1500 m. OMIP1 is forced with the CORE dataset, which consists of surface atmospheric states based on the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) atmospheric reanalysis 52 , covering the period from 1948 to 2009. Instead, OMIP2 uses the JRA55-do surface-atmospheric dataset, which is based on the Japanese 55-year atmospheric reanalysis (JRA-55) 53 , covering the period from 1958 to 2018. Relative to CORE, the JRA55-do forcing has an increased temporal frequency (from 6 h to 3 h) and refined horizontal resolution (from 1.875 ∘ to 0.5625 ∘ ). Moreover, JRA55-do has the advantage of being based on a single reanalysis product and hence producing fields that are more selfconsistent 54 . We use JRA-55 when validating wind stress from model simulations. The choice is dictated mainly for consistency with the OMIP2 MMM but also because JRA-55 wind stress and QuickSCAT observations have been found to be very similar 55 , especially in the BCS region. JRA-55 is therefore used as the wind stress reference since the QuickSCAT temporal resolution does not cover the period of study (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004).

DATA AVAILABILITY
The model output from CMIP5 and CMIP6 simulations used in this study are distributed through the Earth System Grid Federation (ESGF). Model output is freely accessible through ESGF data portals after registration. OISST data sets can be downloaded from https://www.ncei.noaa.gov/products/optimum-interpolation-sst, WOA18 data sets from https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/, and the JRA-55 atmospheric reanalysis data sets from https://jra.kishou.go.jp/JRA-55/ index_en.html. Post-processed data used in this study can be obtained from the corresponding author upon request.