Introduction

Characterized by the Antarctic Circumpolar Current (ACC)—a dominant part of the global thermohaline circulation—the Southern Ocean (SO) plays a unique and key role in the global climate system1,2,3. It not only zonally connects three ocean basins—the Pacific, the Atlantic, and the Indian Ocean, but also vertically links the sea surface temperature (SST) with the ocean interior by bringing up deep water through the westerly wind-driven divergence. Normally, the deep water in the SO, consisting of the Northern Atlantic Deep Water (NADW) of high salinity and high oxygen concentration, upwells along the isopycnal1,4 to the sea surface5. One part of the surface flow travels northward to form the Subantarctic Mode Water and Antarctic Intermediate Water (AAIW), while the other moves southward and sinks to be part of the Antarctic Bottom Water6. The SO dynamics is critical in anthropogenic carbon and heat uptake and in redistributing global energy7,8,9. Therefore, understanding the basic state of the SO is among the most important tasks in climate change research. As sustained SO observations of are lacking10, numerical simulations based on coupled climate models are an indispensable and essential way to study the SO.

However, the SO has persistently suffered a systematic warm SST bias in several generations of climate models participating in the Coupled Model Inter-comparison Project (CMIP)11,12,13. This warm bias appears as one of the most prominent biases in global SST modeling12, precluding correct simulations of the current climate such as the location of the inter-tropical convergence zone14, the recent trend of the Southern Annual Mode (SAM)15, and the response of SO SST to SAM16. Moreover, it also lowers the reliability of model projections for future climate change17.

A few previous studies have proposed possible causes of the SO warm SST bias in CMIP models. For instance, ref. 12 found in a large set of CMIP5 models, that the models with a weaker Atlantic Meridional Overturning Circulation (AMOC) tend to have warmer SST in the SO, and thus the common SO warm SST bias can be attributed to an overly weak AMOC strength. In contrast, ref. 18 revealed an extremely weak inter-model correlation between the AMOC strength and the SO SST among CMIP5 models. Instead, they found the surface shortwave radiation is more responsible for the diversified simulations of SO SST and claimed that the warm SST bias in the SO is primarily caused by the overestimation of cloud-related shortwave radiation. These two contradictory explanations obscure a clear judgment of the true source of the common SO warm SST bias.

The warm SST bias anchored in the SO is still prevailing in the latest sixth phase of CMIP (CMIP6) models (Fig. 1a), and the magnitude of such a bias is comparable to those in the previous generations of CMIP models11,12, even though there are generally higher spatial resolutions and more reasonable physical parameterizations in CMIP6 models than in their predecessors19. Moreover, the bias of surface downward shortwave radiation has been revealed to be largely remediated20 and the AMOC has been shown to be reasonably simulated21,22 in CMIP6 models. All these raise a natural question: are the causes suggested to be responsible for the common SO warm SST bias among CMIP5 models still working in CMIP6 models? If not, what are the origins of this persistent common warm SST bias in CMIP6 models?

Fig. 1: Southern Ocean (SO) SST biases.
figure 1

a The multi-model mean (MMM) SST bias (units: °C) in the SO. Black lines are the observational northern (SAF, –0.035 dym) and southern (SACCF, 0.985 dym) boundaries of the ACC, which are defined by AVISO combined mean dynamic topography (MDT_CNES_CLS13) according to Kim and Orsi (2014). Stippling indicates areas where more than 80% of models agree on the sign of bias. The labels in (a) indicate prominent topographic features, including the Kerguelen Plateau (KP), the Macquarie Ridge (MR), the Pacific–Antarctic Ridge (PAR), the Drake Passage (DP), and the Mid-Atlantic Ridge (MAR). b The meridional averaged SST biases within SAF and SACCF for each model (gray lines) and CMIP6 MMM (red line). The colored error bars and diamond boxes denote the MMM SST bias and the standard deviations, from left to right, for the Indian Ocean sector (blue), the Pacific Ocean sector (green), and the Atlantic Ocean sector (purple).

Results

SST bias in the SO

The CMIP6 multi-model mean (MMM) SST bias in the SO displays a zonally nonuniform band-like warm pattern (Fig. 1a). The warm biases are more robust and consistent in the Indian and Atlantic Ocean sectors and weaker in the Pacific Ocean sector, which are similar to the zonal distributions of the common SO warm SST bias in the CMIP3 and CMIP5 models (see Fig. 5a in ref. 11 and Fig. 1a in ref. 12, respectively). Overlaying the observed circumpolar Subantarctic Front (SAF) and the southern Antarctic Circumpolar Current Front (SACCF) (see Methods), we find that the warm SST bias in the MMM, as well as in most individual models, is mainly confined within these two ACC fronts (black lines in Fig. 1a; Supplementary Fig. 1).

The zonally nonuniform SST bias pattern is to some extent associated with the complex topographies in the SO (Fig. 1a). For instance, extreme warm biases appear in the Kerguelen Plateau, the Macquarie Ridge, the Pacific–Antarctic Ridge, the Drake Passage, and the Mid-Atlantic Ridge. These biases coincide with the hot spots in the SO, including the steepening isopycnals, intensified eddy activities, enhanced cross-front eddy heat exchanges, and facilitated upwelling of deep water3,23,24. In addition, there are also some regional cold SST biases in the eastern boundary current areas, probably arising from model deficiencies in simulating these currents and associated local dynamics. The origins of these regional SST biases cannot be ascertained since variables essential to study the corresponding subgrid-scale processes are unavailable from climate models. Thus our focus here is to explore the origins of the common warm SST bias in the SO from the perspective of large-scale thermodynamics and dynamics.

By averaging the SST biases between the two observational ACC fronts (Fig. 1b), it again appears that the warm SST biases in the Atlantic and Indian sectors (reaching up to ~2 °C in MMM) are stronger than that in the Pacific sector (~1 °C in MMM). The inter-model standard deviation (STD) of SST bias in the Pacific sector is almost twice as much as those in the Indian and Atlantic sectors, indicating higher inter-model consistency of the warm SST bias in the latter two than in the former (Fig. 1b, colored error bars). The weak inter-model consistency in the Pacific sector is mainly caused by four models (EC-Earth3, EC-Earth3-Veg, MIRCO6, and MIRCO-ES2L) with extreme warm biases, which show larger warm biases in the whole SO than other models.

Effect of surface heat flux

Firstly we examine the biases of surface heat fluxes and their role in causing the common SO warm SST bias. Compared with the net surface heat flux (Qnet, positive downward) from ECMWF Reanalysis v5 (ERA5) dataset, the CMIP6 MMM Qnet displays an overall negative bias over the SO with high inter-model consistency (Fig. 2a, f), which generally opposes the positive MMM SST bias, with the spatial correlation between the two reaching up to −0.70 between 40°S and 60°S. This indicates that the Qnet bias does not contribute to the common warm SST bias in the SO but actually responds to it. For further verification, we choose six additional observational/reanalysis Qnet products to compute the MMM Qnet bias (Supplementary Fig. 2). The results show consistent negative Qnet biases, though the magnitude of the bias varies from one another due to observational uncertainties25.

Fig. 2: Surface heat flux biases.
figure 2

The left column is the MMM spatial patterns of (a) net heat flux (Qnet) bias, (b) shortwave radiation (SW) bias, (c) longwave radiation (LW) bias, (d) sensible heat flux (SH) bias, and (e) latent heat flux (LH) bias. The two black lines in (ae) are same as in Fig. 1. Stippling indicates areas where more than 80% of models agree on the sign of bias. The right column is the meridional-averaged (f) Qnet biases, (g) SW biases, (h) LW biases, (i) SH biases, (j) LH biases (units: W m−2) for individual models (gray lines) and MMM (bold red lines). The sign of surface heat flux terms is positive when going downward.

By decomposing the Qnet bias into biases of different flux terms, it becomes clear that the MMM shortwave radiation (SW) displays a totally negative bias over the SO (Fig. 2b), whereas the surface longwave radiation (LW) exhibits a positive bias that nearly counterbalances the bias of SW (Fig. 2c). In addition, both the MMM surface sensible heat flux (SH) and latent heat flux (LH) have negative biases between the two ACC fronts with similar spatial patterns but of different magnitudes (Fig. 2d, e). This indicates that the negative bias of Qnet is mainly due to turbulent heat flux biases, which release more heat from the ocean surface to balance the warm SST bias. Moreover, the negative bias of LH is stronger than that of SH, implying a dominant role of LH bias in responding to the warm SST bias. Most of the CMIP6 models yield the same results as the MMM (Fig. 2g–j), except for the two outliers with extreme warm SST bias in the SO, EC-Earth3-Veg and EC-Earth3, which show positive SW biases over the whole SO (Supplementary Fig. 3).

Effect of AMOC strength

Next, we turn to the effect of AMOC strength on the SO SST bias. Here and in the following analysis, the SO SST bias in each CMIP6 model is represented by the area-averaged SST bias between the observational SAF and SACCF, unless otherwise specified. As shown in Fig. 3a, there appears to be a significant negative inter-model relationship between the SO SST bias and the AMOC strength among 18 CMIP6 models, consistent with the result from CMIP5 multi-models12, with weaker AMOC strength tending to have stronger warm SST bias. However, more than 70% of the chosen models have positively biased AMOC strength, and the MMM bias of AMOC strength is also slightly positive (1.41 Sv, red diamond in Fig. 3a). This implies that the common SO warm SST bias cannot be caused by the common bias of AMOC strength.

Fig. 3: Inter-model relationships.
figure 3

a Inter-model relationship between AMOC strength biases (units: Sv) and SST biases in the SO. Here the AMOC strength is defined as the maximum value of the meridional overturning streamfunction in the Atlantic Ocean at 26°N. The red diamond box denotes the MMM AMOC strength bias, and the red error bar is one standard deviation of the 18 models. b Inter-model relationship between Northern Atlantic (30–60°N,70°W–20°E) deep ocean (1000–3000 m) temperature biases and SST biases in the SO for models with reasonable AMOC strength. c Inter-model relationship between AMOC strength biases and Northern Atlantic deep ocean warm biases. The solid black line denotes the linear regression, and the correlation coefficient and the p value determined by a two-side Student’s t test are also shown in each panel.

There are four models with their AMOC strength exceeding one inter-model STD. Two of them (FGOALS-g3 and GISS-E2-1-G) have their AMOC much stronger than the MMM while the other two (MIRCO-ES2L and MIRCO6) much weaker. It seems that only for these models with largely biased AMOC strength, there exists a strong negative correlation between the biases of AMOC strength and SO SST. If these four models are excluded from the model ensemble, the inter-model correlation coefficient between the bias of AMOC strength and that of SO SST reduces to as low as −0.04. Thus, the effect of AMOC strength on the SO SST is highly model-dependent, which may explain the conflicting arguments regarding the role of weakening AMOC strength in causing the SO warm SST bias12,18.

Effect of warm bias in the Northern Atlantic deep ocean

In the absence of a systematic AMOC strength bias that accounts for the common SO warm SST bias, together with a generally reasonable representation of surface heat flux and wind stress over the SO in CMIP6 models26 (Supplementary Fig. 4), we turn to remote factors and find that the common warm bias in the deep Northern Atlantic could be partly responsible for the common SO warm SST bias through AMOC transport. As shown in Fig. 4, the MMM zonally-averaged ocean temperatures in all three ocean sectors show overall warm biases throughout the water column with somewhat similar spatial patterns. Larger biases appear in the deep ocean and clime upward and southward to the upper ocean nearly along the isopycnals (Fig. 4), indicating a deep source of bias from the north. The warm biases in the Atlantic sector generally appear in the NADW, where salinity exceeds 34.7 psu (Supplementary Fig. 5), rather than in the AAMW/AAIW.

Fig. 4: Zonal mean temperature biases in three ocean sectors.
figure 4

Multi-model mean (MMM) zonally averaged ocean temperature biases (units: °C) in (a) Atlantic (70°W–20°E), (b) Indian (20°E–160°E), and (c) Pacific (160°E–70°W) ocean sectors. The black contours in each plot are the MMM potential density (units: kg m−3, with interval 0.2 kg m−3). Stippling indicates areas where more than 80% of models agree on the sign of bias.

The global deep ocean temperature bias averaged from 1000 m to 3000 m reveals that the largest warm bias is located in the Atlantic Ocean, and the warm biases in the SO gradually decrease from offshore of Argentina in the Atlantic Ocean sector to the east (Fig. 5). Such a spatial distribution of the global deep ocean warm bias resembles that of the dissolved oxygen concentration (black contours in Fig. 5), which is often used as an indicator for the transport of NADW starting from the deep convective region in the Northern Atlantic Ocean27. This implies that the warm bias in the North Atlantic deep ocean can be advected almost adiabatically to the surface of the SO, contributing to the warm SST bias there. Indeed, as shown in Fig. 3b, all of the models have a warm bias in the Northern Atlantic deep ocean (30–60°N, 70°W–20°E, 1000–3000 m), and the models with larger warm bias in the Northern Atlantic deep ocean tend to have larger SO warm SST bias (Fig. 3b), with a significant correlation of 0.54 (two-sided Student’s t test, p < 0.05). The warm bias in the Northern Atlantic deep ocean is not linked to the bias of AMOC strength, since there is no significant correlation between the two (Fig. 3c).

Fig. 5: Deep ocean temperature biases.
figure 5

Multi-model mean (MMM) ocean temperature biases (units: °C, shading) averaged between the depths of 1000 m and 3000 m. The black contours are the MMM climatological dissolved oxygen concentration (µmol kg−1) at the same depth. Stippling indicates areas where more than 80% of models agree on the same sign of bias.

To further demonstrate the processes through which the warm bias in the Northern Atlantic deep ocean impacts the SO warm SST bias, we perform a set of idealized model experiments with passive-tracer (see Methods). After tracers are released through the entire water column in the northern Atlantic Ocean (north of 60°N), they travel with the deep flow and arrive at the SO in about 200 years (Supplementary Fig. 6), and finally upwell to the surface (Fig. 6). The experimental results show that high concentrations at the surface are primarily confined to the SO (Fig. 6a–c), whereas higher concentrations in the deep ocean are in the Atlantic Ocean (Fig. 6d–f), which is very similar to the distribution of warm bias in the deep ocean shown in Fig. 5. Additionally, the zonal mean vertical structures of the tracer concentration also bear resemblance to that of the warm bias shown in Fig. 4, with the largest concentration located at depths about 1000–3000 m declining from the north to south (Fig. 7). These idealized passive-tracer experiments confirm that the common Northern Atlantic deep ocean warm bias can significantly contribute to the common warm SST bias in the SO.

Fig. 6: Tracer concentration in idealized MITgcm experiments.
figure 6

Tracer concentration at the surface (ac) and in the deep ocean (averaged between 1000–3000 m) (df). The concentration of each tracer is relaxed to one (top), two (middle), and three (bottom).

Fig. 7: Zonal-mean tracer concentration in idealized MITgcm experiment.
figure 7

Zonally averaged tracer concentration in the Atlantic Ocean sector for three MITgcm tracer experiments, with tracer concentration relaxed to (a) one, (b) two, (c) three.

Discussion

Our analyses reveal that the biases of the net surface heat flux and the AMOC strength, which have been suggested to be the primary causes for the warm SST bias in the SO for CMIP5 models, are not the reasons for the same persistent warm bias in CMIP6 models. Instead, the remote warm bias in the Northern Atlantic deep ocean is found to be conducive for the common warm SST bias in the SO through adiabatic transport by AMOC, based on a large ensemble of climate model outputs, various observation/reanalysis datasets, as well as a series of idealized model experiments with passive-tracers. This remote transport of warm bias to the SO is similar, to some extent, to the southward propagation of excessive anthropogenic warming from the deep layer of Nordic Sea to the Southern Atlantic28. The cause of the NADW warm bias is still unclear, but might be tied to the bias in ocean dynamics, such as unrealistic simulation of deep ocean convection, which will be explored in the future.

It should be noted, however, that the remote common warm bias in the North Atlantic deep ocean can only explain about 30% of the total inter-model variance of the SO SST bias, implying that there must be other causes of the SO warm SST bias remaining to be revealed. For example, the incorrect model representation of ocean eddy fluxes due to coarse spatial resolution could be important for the warm SST bias in the eddy-rich SO. Moreover, regional factors such as complicated local topography that strongly modulate the spatial pattern of the SO warm SST bias, are not considered within the current large-scale framework. Nevertheless, the warm bias in the Northern Atlantic deep ocean, which has been persistently present in almost all CMIP5 and CMIP6 models22,29, as well as in climate models with higher spatial resolution22, have to be taken into account to reduce the SO bias, which is needed not only for simulating the current climate state, but also for projecting the future climate change.

Methods

CMIP6 models

Monthly outputs from 18 CMIP6 models in the historical run are used in the present study. The monthly CMIP6 outputs are SST, net shortwave radiation (SW), net longwave radiation (LW), surface latent heat flux (LH), sensible heat flux (SH), sea water potential temperature, salinity, and overturning mass streamfunction. The net surface heat radiation (Qnet) is the sum of SW, LW, LH, and SH. For each CMIP6 model, we used only one ensemble member. We prioritized the “r1i1p1f1” ensemble member and selected an alternative member when the preferred member was unavailable. The model names, ensemble members, and institutions are listed in Table 1. The current climatology in CMIP6 models is defined as the 33-year long-term mean during the 1982–2014 period.

Table 1 List of models and associated information.

Observational datasets

The model outputs are compared against several observational/reanalysis datasets (both called observations for simplicity in this study). The SST observation is obtained from the Optimum Interpolation SST (OISSTv2)30. The dynamic sea surface height is extracted from MDT_CNES_CLS1331 and the surface heat fluxes (including the LH, SH, LW, and SW) are obtained from seven datasets: the ECMWF Reanalysis v5 (ERA5)32, ECMWF Reanalysis − 40 years (ERA40)33, NCAR/NCEP Reanalysis (NCEP‐NCAR)34, NCEP-DOE AMIP-II reanalysis (NCEP-DOE)35, NCEP Climate Forecast System Reanalysis (NCEP-CFSR)36, the Japanese 55-year Reanalysis (JRA‐55)37, and 20th Century Reanalysis (20Cv2c)38. The three-dimensional ocean temperature is extracted from the Simple Ocean Data Assimilation (SODA3.4.2)39, and the AMOC index was obtained from RAPID40. The observational datasets are detailed in Table 2. Note that the time periods of some observational datasets (such as ERA40, JRA55, and NCEP-CFSR) are slightly different from each other (see Table 2). The main conclusions are not sensitive to the observational datasets being used. In this study, the ACC fronts are defined based on the observed dynamic sea-surface height according to Kim and Orsi41, with the circumpolar SAF associated with 0.985 dym isoline and the SACCF with −0.035 dym isoline, respectively. The sign of surface heat flux terms is positive when going downward. All model outputs and observational datasets are interpolated into a 2.5° × 2.5° grid before analyses.

Table 2 List of observational/reanalysis datasets.

Idealized experiments

To verify whether the SO warm SST bias can be attributed to the warm bias in the Northern Atlantic, we employ an ocean model with idealized-geometry based on the Massachusetts Institute of Technology general circulation model (MITgcm)42. The simulation domain, with 1° resolution, spans from 71.5°S to 71.5°N in latitude and 240° in longitude. Two continents divide the domain into the Atlantic, Indian, and Pacific Oceans. The SO is represented by a reentrant channel at the southern edge of the domain. In the vertical direction, the ocean is 4000 m deep everywhere except in the Drake Passage (the submarine sill in the SO), which is 2500 m deep. Detailed information on this model setting refers to ref. 43. In this study, three passive-tracer experiments were conducted. In the northern Atlantic Ocean (>60°N), the tracer concentrations are relaxed to one, two, and three through the entire water column. In each case, the simulation runs for 1000 years, sufficient to reach equilibrium. The model outputs are averaged over the last 30 years of simulation for analysis.