Introduction

There is a debate about how the tropical atmospheric circulation (TAC) responds to external forcing, especially to rising CO2-concentrations that cause global warming. The two major atmospheric circulation systems over the tropical Pacific are the Walker Circulation, a zonal overturning circulation at the equator, and the Hadley Circulation (HC), a meridional overturning circulation in each hemisphere extending from the equator to the subtropics. It has been suggested that in a warming world, the overall exchange of mass between the boundary layer and the mid-troposphere must decrease1. Since much of this exchange occurs in tropical moist convection regions, there the convective mass flux must decrease, which would slow the TAC. The Pacific Walker Circulation (PWC) may have indeed slowed during the period 1861–19922. By performing experiments with an atmosphere model, it was found that SST warming patterns are the main cause of the weakened Walker Circulation over 1950–20093. Finally, the majority of current climate models predict that the tropical Pacific will exhibit an El Niño-like mean-state response to strongly elevated atmospheric CO2-concentrations, i.e. a reduced SST gradient across the equator, with more warming in the east than in the west, and slower PWC4,5,6,7.

Other observational analyses and climate model simulations, covering different time spans, provided evidence for a strengthening PWC. These cover the periods 1900–20078 and 1979-20089, the time since 195010 and the early 1990s11, the 20th century12, and the period 1970–199513. Further, the Held and Soden argument1 was disputed, noting that the connective mass flux argument doesn’t necessarily describe regional circulations12. In line with this, an unprecedented acceleration in the equatorial trade winds over the Pacific since the 1990s has been described14. For the period 1900–1991 an increase in the equatorial Pacific SST gradient was described15 that is consistent with the dynamical thermostat mechanism16, which would tend to strengthen the PWC. The dynamical thermostat involves vigorous upwelling in the eastern equatorial Pacific, which would counteract an externally forced SST warming in this region. The PWC response to external forcing eventually depends on the balance of the evaporative cooling17, cloud cover feedback18, and the ocean dynamical thermostat. The former two are more important over the warm pool in the west, while the latter is more important over the cold tongue in the east. Finally, an inter-basin thermostat mechanism has been suggested19. SST warming in the Indian Ocean drives stronger trade winds over the tropical Pacific that reduce the SST in the eastern equatorial Pacific.

Regarding the HC, a number of studies reported a strengthening or widening in the recent decades (e.g.,20,21,22,23,24). A poleward expansion of the Southern Hemisphere HC appears to be robust in centennial reanalyses25. An HC expansion is observed in state-of-the-art climate models, when the preindustrial CO2-concentration quadruples but with reduced amplitude relative to observations26. In agreement with the latest suite of climate models, it has been found that the Northern Hemisphere HC has considerably weakened over the recent decades27. By conducting simulations with an atmospheric general circulation model, it was shown that natural SST variability accounts for nearly all of the HC widening associated with recent SST evolution28. A possible connection of HC strength to the Atlantic Multidecadal Oscillation (AMO) has been suggested and that HC strengthening may thus have a contribution from internal variability29.

Here we investigate the trends in the SST and atmospheric circulation over the tropical Pacific during the last four decades by means of observations and reanalysis products and assess the performance of current climate models in simulating these trends. We find that the climate models fail to simulate major aspects of the observed trends in the SST and atmospheric circulation and argue that model bias leads to an erroneous response to external forcing.

Results

Sea surface temperature trends during 1981–2020

Tropical Pacific SST plays an important role in regional and also global climate variability, as, for example, during the so-called global warming hiatus around the turn of the century (e.g.,30,31,32). The SST over most ocean regions considerably warmed during 1981–2020, exemplifying the global-scale nature of the warming in an ensemble of six observational data sets (Fig. 1a; Supplementary Fig. 1). Nonetheless, a band with little warming is observed extending from western North America in a southwestward direction to the central equatorial Pacific. A similar band but with cooling in its center is observed south of the equator, whereas pronounced SST warming is observed over the western tropical Pacific. When subtracting the globally-averaged SST trend to derive relative SST trends the effects of global mean warming are reduced and regional differences emphasized (Fig. 1b), which also provides clues about the atmospheric circulation. The relative SST trends show a horseshoe-like warming pattern with warming in the interior North and South Pacific, and cooling in the tropical Pacific, reminiscent of the negative phase of the Pacific Decadal Oscillation (PDO;33,34) or related Interdecadal Pacific Oscillation (IPO;30,35,36).

Fig. 1: Sea surface temperature (SST) trends 1981–2020 (°C per 40 yr).
figure 1

Left a, c, e: Full linear trends, right (b, d, f) relative linear trends (globally averaged trends removed). a, b Observations, (c, d) historical simulations extended with simulations employing the RCP8.5 scenario, (e, f) differences between (a, b) and (c, d). The observed trends are the averages over the trends calculated from 6 different SST datasets, which are shown in Supplementary Fig. 1.

The full (Fig. 1a) and relative SST trends (Fig. 1b) across the equatorial Pacific exhibit an enhanced SST gradient, similar to what is observed during short-term present-day La Niña events, which would tend to enhance the Walker Circulation. Both internal variability and external factors can influence the zonal SST gradient (e.g.,37,38). It has been argued that it is extremely unlikely that the strengthening SST gradient during 1958–2018 is consistent with modeled internal variability39. Eastern equatorial Pacific warming may have been delayed by both aerosols and dynamical thermostat response to CO2-forcing, driving a stronger SST gradient40. The SST cooling in the eastern equatorial Pacific during the last decades could also be a transient phenomenon. An El Niño-like response to CO2-forcing, with a reduced SST gradient, may prevail in the long run due to off-equatorial influences41,42. This would fit the idea of a two-timescale response to external forcing43.

We analyze the SST trends from a number of (extended) historical simulations (Methods, HIST hereafter) with state-of-the-art climate models from CMIP6 employing observed external forcing (Supplementary Table 1). There are major differences between the modeled SST trends (Fig. 1c, d) and the observations (Fig. 1a, b). For example, the trend pattern derived from the models lacks notable zonal gradients and spatial variations are mainly in the meridional direction. Further, the models simulate warming across the equatorial Pacific with a maximum in the east where cooling was observed. Moreover, the models lack the band of minimal warming (relative cooling) over the northern and the band with cooling in its center over the southeastern basin. The model-data differences calculated from the absolute (Fig. 1e) and relative SST trends (Fig. 1f) are similar, as expected, and resemble the negative phase of the PDO/IPO. As the multi-model mean calculated from HIST is an estimate of the externally forced SST signal, the model-data differences could suggest a major contribution of long-term internal variability to the tropical Pacific SST trends. On the other hand, the differences could also point to major deficiencies in the climate models. The present divergence between the historical simulations and the observed trends likely either reflects an error in the models’ forced response or an underestimate of the multidecadal internal variability by the models44.

Concurrent trends in atmospheric circulation

We investigate the trends in the PWC and Pacific HC (PHC) during 1981–2020 in ERA5 reanalysis45, a dataset combining vast amounts of historical observations into global estimates using advanced numerical modeling and data assimilation systems that serves here as validation data. The annual average zonal overturning streamfunction over the equatorial Pacific (Methods; Fig. 2a, contours), a measure of the mean-state PWC, is clockwise (positive streamfunction) with ascent over the western and subsidence over the eastern basin, implying easterly low-level winds across the basin. A PWC strengthening is observed during the last four decades (Fig. 2a, color) in all seasons (Supplementary Fig. 2), which is strongest in boreal winter (December-February, DJF) and fall (September-November, SON), the seasons when the equatorial SST gradient and large-scale equatorial atmosphere-ocean coupling are strong.

Fig. 2: Trends in Pacific Walker Circulation (PWC).
figure 2

ad Linear trend 1981-2020 in the PWC, as expressed by the zonal overturning streamfunction, in ERA5, AMIP (1979-2014), extended historical and 1%-CO2 runs. Stippling indicates in (a) a significant trend at the 90% confidence level, and in (bd) that at least two thirds of the models agree on the sign of the trend. eh Linear trend distribution of the PWC index (zonal overturning streamfunction averaged over the box 145°W – 180°, 700 hPa – 200 hPa) in CMIP6-AMIP, extended historical runs, 1%-CO2 runs, and in the preindustrial control runs (PICTL). The ensemble-mean trends are given by the light blue vertical lines and the trend from ERA5 by the yellow vertical lines.

An ensemble of atmosphere models forced by the observed SSTs during 1979–2014 participating in the Atmosphere Model Intercomparison Project (Supplementary Table 1, AMIP hereafter) also simulate a PWC strengthening in the multi-model mean (Fig. 2b, color). A PWC index is defined (in ref. 46 Methods, box in Fig. 2a–d). The multi-model mean trend in the PWC index is in agreement with the trend from ERA5, demonstrating a strong influence of the SST on the PWC (Fig. 2e). In contrast, a PWC weakening is observed in HIST in the multi-model mean (Fig. 2c). Although there is a considerable spread among the models, the trend in the PWC index from ERA5 is well outside the trend distribution obtained from HIST (Fig. 2f), which is consistent with a previous study reporting that climate model simulations are inconsistent with the PWC strengthening observed during 1980–201247. A PWC weakening in the multi-model mean is also observed in an ensemble of global warming simulations with a number of CMIP6 models (Supplementary Table 1; Fig. 2d, g), in which the atmospheric CO2-concentration rises at a rate of 1% per year and quadruples at the end (1%-CO2 hereafter). This constitutes a radiative forcing that is much stronger than that in HIST. The probability that the trend in the ERA5 index is due to internal variability, as estimated from preindustrial control integrations of the CMIP6 models (Supplementary Table 1, PICTL hereafter), is very small (Fig. 2h). In summary, the CMIP6 models are virtually incapable of simulating a PWC trend which is consistent with that observed during the last four decades.

The mean-state PHC consists of two counter-rotating cells (Fig. 3a, contours), one centered north (positive streamfunction) and one south of the equator (negative streamfunction), with upward motion between them reflecting the Intertropical Convergence Zone (ITCZ). The PHC exhibited a strengthening of its southern cell, with the largest negative streamfunction trends in the lower equatorial atmosphere (Fig. 3a, color). The low-level negative streamfunction trends over the equatorial Pacific are most pronounced in June-August (JJA) and especially September-November (SON) (Supplementary Fig. 3c, d). Only a modest strengthening of the northern cell was observed.

Fig. 3: Trends in Pacific Hadley Circulation (PHC).
figure 3

ad Linear trend 1981–2020 in the PHC, as expressed by the meridional overturning streamfunction, in ERA5, AMIP (1979–2014), extended historical and 1%-CO2 runs. Stippling indicates in (a) a significant trend at the 90% confidence level and in (bd) that at least two third of the models agree on the sign of the trend. eh Linear trend distribution of the PHC index (meridional overturning streamfunction averaged over the box 5°S – 8°N, 925 hPa – 500 hPa) in CMIP6-AMIP, extended historical runs, 1%-CO2 runs, and in the preindustrial control runs (PICTL). The ensemble-mean trends are given by the light blue vertical lines and the trend from ERA5 by the yellow vertical lines.

In AMIP (Fig. 3b), HIST (Fig. 3c) and 1%-CO2 (Fig. 3d), the trend patterns (color) differ considerably from that in ERA5 (Fig. 3a) and also among the ensembles. The trend in the PHC is not captured in AMIP showing mostly opposite streamfunction trends (Fig. 3b) than ERA5 (Fig. 3a). This suggests that in contrast to the PWC, SST forcing does not account for the observed trends in the PHC. The trends in the PWC in HIST (Fig. 2c) and 1%-CO2 (Fig. 2d) are consistent with each other though inconsistent with ERA5, whereas the trends in the PHC differ markedly between HIST (Fig. 3c) and 1%-CO2 (Fig. 3d), exemplifying the complexity of the PHC variability. A PHC index is defined as the average of the streamfunction over the region 5°S – 8°N, 925 hPa – 500 hPa (box in Fig. 3a–d). The three corresponding trend distributions of the PHC index (Fig. 3e–g) have a mean close to zero, indicating an overall little sensitivity of the models’ PHC to SST or external forcing. Finally, the trend in the PHC index from ERA5 is outside the distribution from PICTL (Fig. 3h) and thus inconsistent with the models’ internal variability. In summary, the AMIP and CMIP6 models are virtually incapable of simulating the PHC trend observed during the last four decades.

The near-surface (10 m) wind trends over 1981–2020 from ERA5 (Fig. 4) and normalized zonal and meridional components exhibit stronger equatorial easterlies (Supplementary Fig. 4a) and cross-equatorial flow (Supplementary Fig. 4b) that reflect the intensified PWC and low-level southern PHC cell (Figs. 2a3a, respectively). A slight strengthening of the northeasterly trade winds is observed in the east (Fig. 4), consistent with the minor strengthening of the PHC’s northern cell (Fig. 3a). The stronger equatorial easterlies and cross-equatorial winds are highly significant in the centers of action with respect to the internal variability derived from PICTL, exceeding several standard deviations of the models’ internal variability (Supplementary Fig. 4a, b).

Fig. 4: Trends in near-surface (10 m) winds.
figure 4

Shown are the linear trends 1981-2020 from ERA5. Arrows indicate wind direction and color wind speed-magnitude (ms-1). The box indicates the region over which the wind index is defined (9.5°S-5.5°S, 148°W-106°W).

The salient feature in Fig. 4 is the cross-equatorial wind trends over the central Pacific. A wind index is defined as the average of the meridional 10 m wind over the southeastern equatorial Pacific (box in Fig. 4 and Supplementary Fig. 5a). The 40-year trend 1981–2020 in the ERA5 index is outside the internal variability derived from PICT (Fig. 5a). The ERA5 trend is also outside the trend distribution derived from HIST (Supplementary Fig. 5b) and 1%-CO2 (Fig. 5b). Further, the distribution derived from 1%-CO2 is almost symmetric and does not differ much from that obtained from PICT. Thus, the trend in the ERA5-wind index is not captured by the CMIP6 models.

Fig. 5: Trend distribution of the wind index.
figure 5

Distribution of non-overlapping linear 40-year trends in the wind index calculated over the box shown in Fig. 4. a Distribution calculated from the preindustrial control integrations and (b) from the 1%-CO2 integrations with the CMIP6 models. The trend from ERA5 is given by the yellow vertical lines.

All considered AMIP models simulate positive trends in the wind index (Supplementary Fig. 5c), which, however, are considerably smaller than that from ERA5. Nevertheless, the AMIP simulations pose an advance upon the HIST simulations. As the CMIP6 models do not capture the strengthening of the cross-equatorial winds, irrespective of whether an external forcing is applied or not, while the AMIP simulations do capture some of the observed tendency, we argue that biases in the climate models distort their response to external forcing.

SST sensitivity to surface-wind trends

First, ensemble integrations are conducted with the Kiel Climate Model (KCM, Methods) forced by the pattern of global wind-stress trends 1981–2020 from ERA5. The ERA5 data is affected by ship measurement data and these suffer from a spurious trend48, but satellite data have been used during the period of interest to constrain the ERA5 winds. Therefore, the KCM’s response to the ERA5 winds should be interpreted with some caution. Local air-sea heat exchange is not changed by specifying the wind stresses, but coupled ocean-atmosphere interactions remain fully active. The KCM simulates a La Niña-like SST response, with enhanced equatorial easterlies over the central and western and stronger cross-equatorial winds over most of the equatorial Pacific (Fig. 6a) that reinforce the response to the specified wind stresses. The tropical Pacific mean state experiences a La Niña-like change when cross-equatorial wind strengthening is applied to a climate model49. In the KCM, it is the zonal wind-stress component accounting for most of the SST change (Supplementary Fig. 6a) as opposed to the meridional component (Supplementary Fig. 6b). The magnitude of the SST change in the center of the cold anomaly amounts to about 0.4 K (Fig. 6a). Wind stresses thus could have played an important role in counteracting externally forced SST warming over the equatorial Pacific, consistent with previous studies50,51.

Fig. 6: Influence of the wind trends on sea surface temperature (SST).
figure 6

a SST (K) and 10 m wind response (ms-1) simulated by the Kiel Climate Model (KCM) forced by the observed wind-stress trends 1981–2020. Shown is the ensemble-mean difference of the 40-year averages between the forced experiments and an unperturbed control integration over the same period. b SST anomalies (K) due to local wind-induced latent heat-flux change calculated with the GREB model forced by the ERA5 surface-wind trends 1981–2020. The global average was removed.

Second, we use the globally resolved energy balance model (GREB, Methods) to study the influence of local wind-induced latent heat-flux changes on the SST. When applying the surface-wind trends 1981–2020 from ERA5 to GREB, two bands of SST cooling develop, one emanating from the northern and the other from the southern eastern subtropical Pacific, where the southern band is more pronounced. The two bands join in the western equatorial Pacific. In addition to wind stress-forced ocean dynamical processes acting mostly at the equator, enhanced wind-induced latent heat loss could have been another important factor opposing external SST warming especially off the equator.

Discussion

Large parts of the tropical Pacific SST warmed less than the global average SST or even cooled during the last four decades. Here we show that a strengthening atmospheric circulation over the tropical Pacific is a key to understand the SST trends. The associated changes in the trade winds counteracted externally forced SST warming by influencing ocean dynamical processes as well as air-sea heat exchange. Current climate models fail to simulate the observed changes in the atmospheric circulation, irrespective of whether an external forcing is applied or not, which is a major reason for the models’ inability to represent the observed SST trends.

Biases must be considered when assessing the climate models’ performance. There are two prominent tropical Pacific SST biases that influence the climate dynamics in that region. First, an equatorial cold bias and second, a warm bias in the southeastern tropical Pacific52. The major rainfall bands over the tropical Pacific are also biased in many climate models53, often exhibiting a double-ITCZ structure that can be traced back to the SST biases54. Major biases in tropical Pacific interannual variability have been described too55,56.

Mean-state biases may distort the climate models’ response to historical forcing57. By returning to the fundamental dynamics and thermodynamics, and avoiding sources of model bias, the erroneous warming in state-of-the-art climate models was shown to be a consequence of the cold bias of their equatorial cold tongues58. A sensitivity of the ocean thermostat to horizontal resolution in climate models was observed in CO2-response experiments42. A CO2-response sensitivity of the tropical Atlantic ocean-atmosphere system to the size of the SST bias over the southeastern basin has been observed in the KCM59. Analyses of seasonal forecast ensembles revealed that the SST trends in the central and eastern equatorial Pacific are more positive than those observed over the period 1982–2020, thereby reducing the cold bias60,61. An overestimated surface net heat flux and underestimated local SST-cloud negative feedback favor an El Niño-like warming bias in climate models62. Regarding the HC, the discrepancy between climate models and reanalyses may not stem from internal climate variability or biases in models but could also be related to artifacts in the representation of latent heating in the reanalyses63. Nevertheless, given the large mean-state biases and the flawed interannual variability in most climate models, it is reasonable to assume that their response to external forcing over the tropical Pacific is biased and erroneous.

Long-term internal variability, especially variability linked to the PDO, could be an alternative explanation for the discrepancies between the climate models and observations. The PDO is not a single phenomenon but is instead the result of a combination of different physical processes, tropical and extratropical64. Here we stick to the PDO as a simple measure of internal variability. The PDO index was in its negative phase for most of the time during the last four decades (Supplementary Fig. 7a), supporting SST cooling over large parts of the tropical Pacific (Supplementary Fig. 7b). The largest tropical SST changes linked to the PDO are typically centered at the equator whereas the largest SST cooling during the last four decades is observed south of the equator (Fig. 1a). Such discrepancies cast doubts about the PDO being the primary cause of the recent multidecadal SST trends. Finally, the observed trends in the PWC, PHC and trade winds are very unlikely due to internal variability as estimated by climate models, but the models may underestimate the variability65,66.

The failure of current climate models to represent some of the major tropical Pacific SST trends during the last four decades raises questions about the models’ ability to reliably project future climate over the tropical Pacific sector, which would also have global implications. The analyses presented in this study support the notion that model bias is the key factor hindering a realistic response to external forcing over the tropical Pacific. Enhancing confidence in 21st century tropical Pacific and global climate projections requires a better understanding of the atmospheric circulation over the tropical Pacific, its internal variability and response to external forcing, and interplay with SST.

Materials and Methods

Sea-surface temperature observations

We use an ensemble of six SST data sets, namely COBE67, ERSST68, GISS69, HadISST70, KAPLAN71, and OISST72 to calculate the observed trends for 1981–2020. Each data set was first bi-linearly remapped to a common 1° × 1° horizontal grid. The monthly anomalies were calculated and averaged to provide annual anomalies. Then, the linear trends were calculated using the least-squares method and tested for statistical significance using the Student’s t-test after accounting for autocorrelation in the time series73. Finally, the trends and significance were averaged over all six SST datasets. A global SST average was calculated from yearly anomalies by weighting with the cosine of the latitude. The global average was subtracted at each grid point before calculating linear trends and statistical significance to obtain the relative trends.

CMIP6 control simulations

We analyze preindustrial control integrations from the Coupled Model Intercomparison Project phase 6 (CMIP674,). The models used in this study are listed in Supplementary Table 1.

CMIP6 historical simulations

We analyze historical simulations with CMIP6 models employing observed external forcing 1850–2014 and extended to 2020 with data from simulations employing the IPCC-SSP585 scenario, roughly representing the observed radiative forcing from 2015 onward. The models are listed in Supplementary Table 1.

CMIP6 1%-CO2 simulations

We analyze CMIP6 models which were forced by a CO2-concentration that rises at a rate of 1% per year and quadruples at the end. The models used in this study are listed in Supplementary Table 1.

AMIP simulations

We analyze atmosphere models from the Atmosphere Model Intercomparison Project (AMIP) as part of CMIP6 that were forced by observed SSTs and sea-ice concentrations for 1979–2014. The models used in this study are listed in Supplementary Table 1.

Influence of surface-wind trends on SST

The surface winds influence the SST in two ways, through the wind stress and the air-sea heat exchange. The wind stress trends can be directly specified to a climate model. However, the trends in the air-sea heat exchange cannot be easily specified to a climate model due to drift problems. We therefore proceed in two steps, as described below.

Kiel Climate Model

First, we use the Kiel Climate Model (KCM75,) to investigate the influence of the wind-stress trends on the SST. The KCM version used here employs ECHAM5 as an atmospheric component, with a horizontal resolution of T42 (2.8° × 2.8°) and 31 vertical levels. The ocean-sea ice component is NEMO on a 2° Mercator mesh (ORCA2) horizontally, with increased meridional resolution of 0.5° near the equator and 31 vertical levels. We use for reference a multi-millennial control simulation with preindustrial CO2-concentration amounting to 286.2 ppm. An ensemble of 20 wind stress-forced simulations is conducted. Initial conditions for the wind-forced simulations are taken every 50 years from the preindustrial control integration (after allowing for an about 1000 years long spin-up phase). The wind stress-forced simulations are performed by adding the wind-stress trends for 1981–2020 from ERA5 to the model and they have a length of 40 years. Ensemble-mean differences between the 40-year averages of the forced integrations and of the control run over the same period are shown in Fig. 6a.

GREB model

Second, we use the Globally Resolved Energy Balance (GREB) model to investigate the influence of the surface-wind trends on the SST. GREB allows to isolate the SST change solely due to local wind-induced latent heat-flux change. GREB is a flux-corrected energy balance climate model that realistically simulates the mean climate state and exhibits a global warming SST response consistent with current climate models76. GREB consists of prognostic equations for the temperature in 3 layers (mixed layer ocean, surface, lower atmosphere) and atmospheric humidity. The wind climatology is prescribed, as are cloud cover and ocean mixed layer depth. The surface-wind climatology was changed by adding the trends for 1981–2020 from ERA5.

Zonal and meridional streamfunction

The zonal and meridional streamfunction are used as measures for the Pacific Walker and Hadley cells, respectively:

$$\Psi =2\pi a\, \mathop {\int} \nolimits_{0}^{p} {u}_{D}\frac{{dp}}{g}$$

Here uD is the divergent component of the zonal or meridional wind, a the Earth’ radius, p the pressure and g the constant gravity. The zonal streamfunction is averaged over the latitude band 5°N-5°S6, and the meridional stream function over the Pacific is averaged over the longitude band 120°W-70°E.