Reconciling opposite trends in the observed and simulated equatorial Pacific zonal sea surface temperature gradient

Abstract

The reasons for large discrepancies between observations and simulations, as well as for uncertainties in projections of the equatorial Pacific zonal sea surface temperature (SST) gradient, are controversial. We used CMIP6 models and large ensemble simulations to show that model bias and internal variabilities affected, i.e., strengthened, the SST gradient between 1981 and 2010. The underestimation of strengthened trends in the southeast trade wind belt, the insufficient cooling effect of eastern Pacific upwelling, and the excessive westward extension of the climatological cold tongue in models jointly caused a weaker SST gradient than the recent observations. The phase transformation of the Interdecadal Pacific Oscillation (IPO) could explain ~ 51% of the observed SST gradient strengthening. After adjusting the random IPO phase to the observed IPO change, the adjusted SST gradient trends were closer to observations. We further constrained the projection of SST gradient change by using climate models’ ability to reproduce the historical SST gradient intensification or the phase of the IPO. These models suggest a weakened SST gradient in the middle of the twenty-first century.

Introduction

The equatorial Pacific zonal sea surface temperature gradient (SST gradient) impacts the global mean temperature and is a pacemaker of global warming (Kosaka and Xie 2013). Observational records show a recent strengthening of the SST gradient, coupled with an intensification of the atmospheric Walker circulation and a strengthening of the Pacific trade winds (Cane et al. 2009; Heede and Fedorov 2023; McGregor et al. 2014), significantly influencing the El Niño-Southern Oscillation (Collins et al. 2010) and global ocean heat uptake (England et al. 2014; Kosaka and Xie 2013). However, climate models tend to simulate a weakening of the SST gradient during the twentieth century (Coats and Karnauskas 2017; Lee et al. 2022). It is exceedingly rare that the latest generation of models can produce the observations-based SST trends from 1958 to 2018 (Seager et al. 2022). This common failure may reduce the credibility of projected SST gradient changes.

The discrepancy between models and observations has been claimed to be related to model biases, such as the biases in the equatorial cold tongue (Li et al. 2016a; Seager et al. 2019; Zhou and Xie 2015), cloud-radiation feedback (Ying and Huang 2016), and effects of inter-basin warming contrast across the three tropical oceans. The negative local SST-cloud feedback and the surface net heat flux in the equatorial Pacific can jointly explain ~ 50% of the difference between observed and simulated SST gradient trends (Luo et al. 2017). It was found that SST gradient trends would be overestimated by ∼0.52 ℃ for every 1 °C global SST warming due to 13 common model biases (Tang et al. 2021).

In addition to the possible influence of model biases, coupled models' uncertainty of the simulated recent SST gradient trend can be attributed to internal variabilities, such as the Interdecadal Pacific Oscillation (IPO) (Coats and Karnauskas 2017; Olonscheck et al. 2020; Watanabe et al. 2020). The SST gradient trends in most models and observational datasets are insignificant relative to internal variability (Coats and Karnauskas 2017). In addition to internal variability, the enhanced tropical Indian Ocean warming (Luo et al. 2012) and the Atlantic Ocean warming (Kucharski et al. 2011; Li et al. 2016b; McGregor et al. 2014) in recent decades could strengthen the Pacific Walker Circulation and the SST gradient by inducing inter-basin interactions.

In future projections, the SST gradient will be weakened with a certain degree of uncertainty (Coats and Karnauskas 2017; Collins et al. 2010; IPCC 2021; Kim et al. 2014; Liu et al. 2017). Climate models project a strengthening of the SST gradient after removing the impacts of 13 common biases (Tang et al. 2021). The uncertainty of the projected SST gradient can be considerably reduced when the initial IPO phase (Bordbar et al. 2019) or the initial oceanic state of the Pacific Ocean (Watanabe et al. 2020) is known and well-represented in the model.

There are several hypotheses of the dynamic mechanism for both negative and positive SST gradient changes. The strengthening of the SST gradient is supported by ocean thermostat mechanism (Cane et al. 2009; Clement et al. 1996), "iris effect" of cirrus contraction (Lee et al. 2022; Lindzen et al. 2001; Su et al. 2017), decreasing in ENSO skewness related to the increased upper ocean stability (Kohyama 2017), and the tropical ocean interaction (Kucharski et al. 2011; Luo et al. 2012; McGregor et al. 2014). The physical mechanisms for the weakening of the SST gradient include the energy balance theory of the water cycle (Held and Soden 2006; Vecchi and Soden 2007), the nonlinear correlation between evaporation and temperature (Knutson and Manabe 1995; Ying et al. 2016), and cloud radiation feedback mechanism (Meehl and Washington 1996; Ramanathan and Collins 1991; Song and Zhang 2014).

In summary, both model biases and the climate system's internal variability affect the SST gradient's simulation and projection. However, their respective contributions are still unclear. Therefore, based on the Coupled Model Intercomparison Project phase 6 (CMIP6) models and large ensemble simulations from the US CLIVAR Working Group on Large Ensembles, we investigated the possible impacts of model biases on the recent SST gradient trend simulation and underlying physical mechanisms. We also explored to what extent the observed SST gradient strengthening could be attributed to internal variability. Our results showed that both the model biases and phase changes of the IPO contributed to past and future SST gradient changes.

Data and methods

Observational data sets

SST observations were taken from the Extended Reconstructed SST v5 (ERSST) (Smith et al. 2008), Hadley Centre Sea Ice and SST (HadISST) (Rayner et al. 2003) v1.1 data set, and Centennial In Situ Observation Based Estimates SST v2 (COBESST) (Hirahara et al. 2014). We used the average of ERSST, HadISST, and COBESST SSTs as the optimal observational SST estimate (referred to as “SST observation”) to reduce uncertainties. Monthly gridded ocean temperature data were taken from COBE-SST2, having 24 vertical levels above 1500 m (Ishii and Kimoto 2009). SLP observations were taken from NCEP-NCAR Reanalysis 1 data set (Kalnay et al. 1996). We also used monthly wind data compiled with the ERA-Interim reanalysis data set (Berrisford et al. 2011).

CMIP6 models and large ensemble simulations

We used r1i1p1f1 outputs from 32 CMIP6 (Eyring et al. 2016) models (Additional file 1: Tables S1) driven by historical forcings (1850–2014) and extended by the shared socio-economic pathway SSP5-8.5 scenarios (2015–2100). To estimate the relative contributions of external forcing and internal variability to recent SST gradient trends, we used the outputs from 6 large ensembles with a total of 404 members (Additional file 1: Tables S2), including ACCESS-ESM1-5 with 40 members (Mackallah et al. 2022), CanESM5 with 50 members (Swart et al. 2019), EC-Earth3 with 55 members (Wyser et al. 2021), FGOALS-g3 with 110 members (Lin et al. 2022), MIROC6 with 50 members (Tatebe et al. 2019), MPI-ESM with 99 members (Maher et al. 2019).

The observational and historical simulations were selected from 1981 to 2010. The overall period (1981–2010) of the MPI-EMS large ensemble included historical runs for 1981–2005 and the RCP8.5 simulations for 2006–2010. All data were first compiled on a regular 1° × 1° grid through bilinear interpolation, and then the interdecadal signal was extracted using a 9-year running means.

Separate the internal variability and the external forcing

The same radiative forcing drives large ensemble simulations from a single model, the ensemble mean can be considered the response to the external forcing, and the ensemble spread of all simulations can be regarded as the effect of internal variability. Therefore, for the variable X of the member i in the large ensemble simulations, it can be expressed as:

$$X\left(i\right)=X_{forced}+X_{internal}(i){, i=1,2,\dots ,n-1, \, n,}$$

(1)

where ${X}_{forced}$ is the ensemble mean of all simulations, representing the response to external forcing. ${X}_{internal}(i)$ is the difference between the original $X\left(i\right)$ and ${X}_{forced}$, which represents the residual part of $X\left(i\right)$ related to internal variability after subtracting the influence of external forcing.

Metrics

The equatorial Pacific zonal sea surface temperature gradient (SST gradient) is defined as the difference in SST anomalies between the equatorial eastern Pacific Ocean (5°S-5°N, 180°-80°W) and the equatorial western Pacific Ocean (5°S-5°N, 110°E -180°) (Watanabe et al. 2020). The mean seasonal cycle of SST gradient calculated for 1981–2010 is removed to calculate monthly SST gradient anomalies. We use the Tripole Index (TPI) for the Interdecadal Pacific Oscillation proposed by Henley et al. (2015). The TPI index is defined as the difference between the SST anomaly over the central tropical Pacific (10°S–10°N, 170°E–90°W) and the average of the SST anomalies over the Northwest (25°N–45°N, 140°E–145°W) and Southwest Pacific (50°S–15°S, 150°E–160°W) (Henley et al. 2015). The time series of the IPO index for each ensemble member is calculated by ${SST}_{internal}(i)$ (see Eq. (1)). The Atlantic Multi-decadal Oscillation (AMO) index is defined as the area-average SST anomalies over the North Atlantic Ocean (0°-65°N, 80°W-0°) after subtracting the global mean (80°S-80°N, 180°E-180°W) SST anomalies (Huang et al. 2020).

The contribution of IPO to SST gradient trend

On the decadal time scale, the contribution of IPO to the SST gradient for member i over the τ period (τ = 1981–2010) can be expressed as (Salzmann and Cherian 2015):

$$\frac{\partial {\Delta SST}_{IPO}\left(i\right)}{\partial t}={r}_{\Delta SST,IPO}\left(i\right)\cdot \frac{\partial IPO\left(i\right)}{\partial t}, i=\mathrm{1,2},\dots n-1, n,$$

(2)

$${r}_{\Delta SST,IPO}\left(i\right)=\frac{\partial \Delta SST\left(i\right)}{\partial IPO(i)},$$

(3)

where ${r}_{\Delta SST,IPO}\left(i\right)$ is the regression coefficient of the SST gradient index regressed onto the IPO index for member i over the τ period. $\frac{\partial IPO\left(i\right)}{\partial t}$ and $\frac{\partial {\Delta SST}_{IPO}\left(i\right)}{\partial t}$ are the IPO trend and the IPO-related SST gradient trend for member i over the τ period, respectively.

The IPO trends simulated by large ensembles are weaker than those in observations (Additional file 1: Fig. S1), linking to a weaker SST gradient trend. To quantitatively estimate the impact of the IPO phase transition on the SST gradient changes, the IPO phase transition for member i was adjusted to the observation. After the phase adjustment, each member can be regarded as affected by the same observational IPO phase transition. The adjusted SST gradient trend for the number i ($\frac{\partial {\Delta SST}_{adj}\left(i\right)}{\partial t}$) is the sum of the externally forced SST gradient trend ($\frac{\partial {\Delta SST}_{forced}}{\partial t}$) plus the internally adjusted SST gradient trend ($\frac{\partial {\Delta SST}_{internal\_adj}\left(i\right)}{\partial t}$)

$$\begin{aligned} & \frac{\partial {\Delta SST}_{adj}\left(i\right)}{\partial t}=\frac{\partial {\Delta SST}_{forced}}{\partial t}\\ &\quad+\frac{\partial {\Delta SST}_{internal\_adj}\left(i\right)}{\partial t}, i=\mathrm{1,2},\dots n-1, n.\end{aligned}$$

(4)

The internal adjusted SST gradient trend includes the internal component of the SST gradient trend ($\frac{\partial {\Delta SST}_{internal}\left(i\right)}{\partial t}$) and the adjustment term (${\alpha }_{internal}\left(i\right)$), which considers the observational IPO phase transition:

$$\begin{aligned} & \frac{\partial {\Delta SST}_{internal\_adj}\left(i\right)}{\partial t}=\frac{\partial {\Delta SST}_{internal}\left(i\right)}{\partial t}\\ &\quad+{\alpha }_{internal}\left(i\right), i=\mathrm{1,2},\dots n-1, n,\end{aligned}$$

(5)

$$\begin{aligned} & {\alpha }_{internal}\left(i\right)=-{r}_{\Delta SST,IPO}\left(i\right)\\ &\quad\times \left(\frac{\partial IPO\left(i\right)}{\partial t}-\frac{\partial {IPO}_{OBS}}{\partial t}\right), i=\mathrm{1,2},\dots n-1, n,\end{aligned}$$

(6)

where $\frac{\partial {IPO}_{OBS}}{\partial t}$ is the observed IPO trend for member i over the τ period. Combining Eq. (4), (5), and (6), the adjusted SST gradient trend is expressed as:

$$\begin{aligned} & \frac{\partial {\Delta SST}_{adj}\left(i\right)}{\partial t}=\frac{\partial {\Delta SST}_{forced}}{\partial t}\\ &\quad+\frac{\partial {\Delta SST}_{internal}\left(i\right)}{\partial t}-{r}_{\Delta SST,IPO}\left(i\right)\\ &\quad\times \left(\frac{\partial IPO\left(i\right)}{\partial t}-\frac{\partial {IPO}_{OBS}}{\partial t}\right), i=\mathrm{1,2},\dots n-1, n.\end{aligned}$$

(7)

The ensemble mean of $\frac{\partial {\Delta SST}_{adj}\left(i\right)}{\partial t}$ includes the effects of external forcing and observational IPO phase transition (Eq. 7). The relative contribution of the observed IPO phase transition is defined as the percentage of the IPO-related SST gradient trend to the observed SST gradient trend.

Results

SST gradient trends in observations and models

Observations consistently showed a strengthening of SST gradient trends (-1.59 ± 0.03 ℃/100 yr) between 1981 and 2010 (Fig. 1a). However, models commonly underestimated the observed SST gradient trends. The SST gradient trend in the CMIP6 ensemble mean was only − 0.02 ± 0.59 ℃/100 yr, and half of the CMIP6 models showed weakened SST gradient trends since 1981. We selected the six CMIP6 models with the largest SST gradient trend strengthening (S models) and weakening (W models). The mean SST gradient trends simulated by the S and W models were − 0.97 ± 0.19 and 0.79 ± 0.19 ℃/100 yr, respectively. The S models still underestimated the observed trend, even though CMIP6 models had a large inter-model spread. We also used six large ensembles with 404 members to check whether the internal variability influences SST gradient trends (Fig. 1b–h). The slight differences in the ensemble means of each large ensemble should imply that external forcing may not be the dominant driver of SST gradient trends in models, whereas the large spreads of ensemble simulations indicate that internal variabilities dominate the trends in models.

Fig. 1

Time series of the 9-year running mean of the equatorial Pacific zonal SST gradient (SST gradient, units: ℃) during 1981–2010 in CMIP6 (a), six large ensembles (b) composed by ACCESS-ESM1-5 (c), CanESM5 (d), EC-Earth3 (e), FGOALS-g3 (f), MIROC6 (g), and MPI-ESM (h). Black lines denote the COBESST, ERSST, and HadISST datasets. The equatorial Pacific zonal sea surface temperature gradient is defined as the difference in SST anomalies between the equatorial eastern Pacific Ocean (5°S-5°N, 180°-80°W) and the equatorial Western Pacific Ocean (5°S-5°N, 110°E -180°). Red lines represent models simulating strengthening trends consistent with the observations (S models), and blue lines represent models showing weakening trends opposite to the observations (W models). Grey shading represents the 5–95% range of internal variability. Numbers in parentheses indicate the sample size

In observations, the positive SST trend in the western tropical Pacific was significantly greater than in the eastern tropical Pacific. This enhanced SST gradient is called a "La Niña-like" pattern (Fig. 2a). In contrast to observations, the CMIP6 ensemble mean showed a spatially uniform warming pattern (Fig. 2b), which disagreed with the observed warming pattern of the zonal SST in the tropical Pacific. S models tended to show a "La Niña-like" pattern, consistent with observations and opposite to the W models, in which the tropical Pacific zonal SST was an "El Niño-like" pattern (Fig. 2c, d). It is noted that even S models still miss the subtropical North Pacific cooling which provides for a wedge like pattern in the observations. The composite SST trend differences between S and W models were similar to the IPO pattern (Fig. 2e), suggesting that the IPO may be a prominent internal variability regulating SST gradient variabilities. A similar result was obtained from MPI-ESM (Additional file 1: Fig. S2e) and FGOALS-g3 (Additional file 1: Fig. S3e) supporting the premise that SST gradient change is modulated by internal variability.

Fig. 2

Spatial distributions of linear trends in annual mean sea surface temperatures (a–e, units: ℃/100 yr), sea level pressure (f–j, shading, units: Pa/100 yr), and wind at 850 hPa (f–j, vectors, units: m/(s·100 yr)) over the Indo-Pacific Ocean during 1981–2010 in observations (a, f), CMIP6 multi-model ensemble means (b, g), CMIP6 S model mean (c, h), CMIP6 W model mean (d, i), and S-W composite differences (e, j). Dotting denotes the trend is statistically significant at the 95% level

Possible causes of SST gradient trend biases

The strength of the SST gradient is closely related to the Pacific Walker circulation through the Bjerknes feedback (Bjerknes 1969). Therefore, we analyzed the sea level pressure (SLP) and low-level wind pattern from 1981 to 2010 to assess the dominant biases responsible for the SST gradient change. The "La Niña-like" SST pattern in observations was accompanied by an enhanced Walker circulation and easterly winds at 850 hPa near the equatorial Pacific (Fig. 2f). The enhanced warming over the Indian Ocean and the Western Pacific Ocean induced a Gill-type response (Gill 1980; Matsuno 1966), which drove the strengthening and westward shift of wind stress over the Pacific, indicating a stronger and westward shifted Walker circulation (Heede et al. 2021).

Compared with the spatial distribution of the observed SLP and low-level wind trends, the CMIP6 ensemble mean showed a weaker change in the Walker circulation and easterly winds off the equatorial Pacific (Fig. 2g). There was a strengthening trend in trade winds in the central equatorial and the southeast Pacific in the S models, with increased latent heat loss, coinciding with the cooling in the above region (Heede et al. 2021). The changes in the S models were more significant than those in the CMIP6 ensemble mean but still weaker than observed changes, indicating a bias in the simulated wind-evaporation-SST feedback (Fig. 2h and Additional file 1: Fig. S4b). The wind-evaporation-SST feedback bias was indistinctive in the W models, in which the strengthened trade wind trend was weaker and located around the dateline (Fig. 2i and Additional file 1: Fig. S4c), contributing positively to the El Niño-like SST warming (Luo et al. 2017). The difference between the S and W models may be partly due to internal variability, similar to the results of the MPI-ESM (Additional file 1: Fig. S2) and FGOALS-g3 large ensemble members (Additional file 1: Fig. S3).

Bias in the simulated Bjerknes feedback would also affect the SST gradient trend. Observations showed a subsurface cooling trend from the dateline to the eastern Pacific, indicating the upwelling of the thermocline to cool the SST in the central and east Pacific oceans (Fig. 3a). The enhanced equatorial trade wind reinforced the eastern Pacific Ocean SST cooling that drove the positive Bjerknes feedback. The tropical and subtropical cells could control the temperature of subsurface water upwelled in the eastern equatorial Pacific, which was strengthened and accompanied by accelerated equatorial undercurrents in reanalysis datasets such as the Ocean Reanalysis System 4 (Jayasankar et al. 2020). The subsurface cooling in the CMIP6 ensemble mean was weaker and located deeper (Fig. 3b). This discrepancy suggested that the upwelling did not effectively cool the eastern Pacific despite the rising of the thermocline. CMIP6 models underestimated the low-frequency variability of the subtropical cells' interior transport convergence and the subtropical wind stress (Graffino et al. 2021). The location of the subsurface cooling in the S models was consistent with observations, but the intensity was weaker (Fig. 3c). In the W models, the location of subsurface cooling was much deeper (Fig. 3d). The contrast in the location and intensity of subsurface cooling was evident in MPI-ESM (Fig. 3f–i) and FGOALS-g3 (Fig. 3j–m), which suggested that the difference between the S and W models was primarily due to internal variability.

Fig. 3

Linear trends of annual mean ocean potential temperature (units: ℃/100 yr) along the equator (5°S–5°N) during 1981–2010, in observations (a) and CMIP6 multi-model ensemble means (b), CMIP6 S model means (c), CMIP6 W model means (d), and S–W composite differences (e). Dots indicate the trend is statistically significant at the 95% level using Student's t-test; the climatological 20 isotherms are shown for three periods: 1981–1990 (green line),1991–2000 (yellow line), and 2001–2010 (purple lines). f–i, same as b–e, but for the MPI-ESM. j–m, same as b–e, but for the FGOALS-g3

Differences in the climatological cold tongue played a vital role in the SST gradient trend (Additional file 1: Fig. S5). The climatological cold tongue extended excessively westward, especially in the W models. This excessive cold tongue bias could induce a positive SST warming bias in the central Pacific by an overly negative shortwave–SST feedback (Ying et al. 2019). The excessive cold tongue bias suppressed deep atmospheric convection, increased the downward solar radiation flux, and reduced surface evaporation and ocean heat loss, which increased the downward net radiation flux and further led to a warm SST in the central Pacific. In addition, the excessive cold tongue bias might influence the SST gradient trend via other physical processes, such as an excessively weak net surface heat flux (Zheng et al. 2012), an excessively shallow thermocline depth (Li and Xie 2012), and an insufficient precipitation bias in the equatorial western Pacific (Du et al. 2015). Bjerknes' feedback could maintain or amplify the above processes (Li and Xie 2014).

Role of IPO in the SST gradient change

Internal variabilities, such as the IPO (Watanabe et al. 2020), modulate the recent strengthening of the SST gradient. We calculated the correlation between the IPO index and the SST gradient index to verify the role of the IPO. Significant linear correlations were identified in all 32 CMIP6 models and the 404 members of the six large ensembles. The mean correlation coefficients of the CMIP6 models and the six large ensembles were 0.39 (P < 0.01) and 0.4 (P < 0.01), respectively, indicating that SST gradient change was related to IPO phase transitions (Additional file 1: Fig. S6). The negative IPO phase was linked to a strengthened SST gradient and vice versa.

The AMO-related Atlantic SST anomalies and atmospheric teleconnection can also affect SST gradients on the multi-decadal time scale (Gan et al. 2023; Wang et al. 2013). However, the simulated equator Pacific zonal SST gradient trends were not significantly correlated with the AMO trends, indicating no robust, coherent relationships between the SST gradient trend and the AMO trend in a 30-year window (Additional file 1: Fig. S7). Hence, although the IPO and AMO could both affect SST gradient change, the IPO is likely the main factor regulating the recent 30-year period of SST gradient strengthening.

The correlation coefficients between the SST gradient and the IPO index varied among the six large ensembles, ranging from − 0.1 to 0.64. Since the simulated IPO phase transitions in the large ensemble members were random, ranging from − 3.7 to 2.8 ℃/100 yr, we adjusted the IPO trends in each member according to the observed IPO trend during 1981–2010 to quantitatively estimate the contribution of the IPO phase evolutions to the recent SST gradient strengthening. We first eliminated the random IPO-related SST gradient trends through linear regression and then superimposed the observed IPO phase transition. After adjustment, the ensemble mean included the response to both the observed IPO phase transition and external forcing (Wu et al. 2021). The SST gradient trends attributed to external forcing ranged from -0.31 to 0.67 ℃/100 yr in the six large ensembles, with an average of 0.26 ℃/100 yr. After adjustment, the SST gradient showed a slightly increasing trend of − 0.56 ℃/100 yr (− 1.04–− 0.02 ℃/100 yr), closer to the observed mean of − 1.60 ℃/100 yr. The IPO-related SST gradient trend was − 0.81 ℃/100 yr (− 1.0–− 0.67 ℃/100 yr), accounting for 51% (range: 42–62%) of the observed trend (Fig. 4). The IPO phase transition was essential to reproduce the observed SST gradient trend successfully.

Fig. 4

SST gradient trends (units: ℃/100 yr) from 1981 to 2010 before and after the observational IPO phase transition adjustments. White, brick red, purple, and blue bars represent observed SST gradient trends, externally forced SST gradient trends, IPO-related SST gradient trends, and total SST gradient trends after adjustments obtained from ACCESS-ESM1-5 (a), CanESM5 (b), EC-Earth3 (c), FGOALS-g3 (d), MIROC6 (e), and MPI-ESM (f) separately. Error bars represent one standard deviation of model intervals. Numbers in parentheses indicate the percentage of the IPO-related SST gradient trend relative to the observed SST gradient trend

Future SST gradient projection

The SST gradient trends showed significant uncertainties among the CMIP6 models and the six large ensembles under the RCP8.5 scenario. Model biases and the IPO phase transition could partly explain these uncertainties. We used two criteria to select the optimal models for future projection constraints. The first was to choose the S and W models according to the largest magnitude of SST gradient strengthening and weakening trends, respectively. The second was to select the S and the W models based on the highest and lowest correlations between the member-simulated and the observed IPO time series, respectively. We chose six members from CMIP6 and twenty members (the top 5%) from the six large ensembles. The S models showed a positive-to-negative phase transition of IPO during 1981–2010, and the W models a negative-to-positive transition. To assess the reliability of this method, we select S models from 1956 to 1985 to explore whether these models could capture the observed SST gradient change during 1981–2010. The S models (Additional file 1: Fig. S8b, c) reproduced the IPO phase shift from 1951 to 2010 and performed better than the ensemble mean (Additional file 1: Fig. S8a). Furthermore, the SST gradient strengthening was well reproduced in the S models, which agreed with the SST gradient observed in 1956–1985 (Additional file 1: Fig. S8b).

In the S models, the SST gradient trends calculated based on the 30-year sliding window changed their signs from negative to positive by the mid-twenty-first century. In contrast, the sign change was opposite in the W models in the same period (Fig. 5a). A similar result was obtained from the large ensembles, indicating that the near-term SST gradient trend was significantly modulated by internal variability (Fig. 5b). A shift toward a positive phase of the IPO in the S models influenced the near-term SST gradient trend (Additional file 1: Fig. S9b, c). The negative-to-positive phase shift of the IPO, superimposed on the externally forced 'El Niño-like' pattern response to projected CO₂ emissions, likely weakened the SST gradient in the mid-twenty-first century in the S models but strengthened it in the W models (Watanabe et al. 2020).

Fig. 5

30-year sliding linear trends of SST gradients (units: ℃/100 yr) simulated by CMIP6 (a, c) and six large ensembles (b, d) from 1981 to 2100. The dots and dashed curves indicate SST gradient trends obtained from historical and SSP5-8.5 simulations by all models (black), S models (red), and W models (blue). Vertical error bars indicate the 66% error range, and the horizontal axis indicates the end year for the 30-year segment. The criterion of S/W models in the first column (a, b) is based on the strengthening/weakening SST gradient trends. The criterion of S/W models in the second column (c, d) is based on the highest/lowest correlation between the model IPO trend and the observed IPO trend

We quantified the contribution of IPO-related internal variability to the SST gradient projection by adjusting the random IPO trends in all 404 members (Additional file 1: Fig. S9a) based on the averaged IPO change projected by the S models (see Eq. (7)). The ensemble mean showed a change from a strengthened to weakened an SST gradient after using the IPO constraint method. The original and the adjusted SST gradient trends based on the phase adjustment method were still different in the late twenty-first century, mainly due to various simulations of the IPO. Differences between the S and W models in the late twenty-first century were also significant. The impact of the IPO on the SST gradient could last until the end of the twenty-first century on a 30-year timescale.

Discussion

Observations indicate a strengthening of the tropical Pacific SST gradient since 1981. However, most state-of-the-art coupled models cannot reproduce this recent strengthening trend. Understanding this discrepancy between models and observations is vital to understanding the present and projected climate in the tropical Pacific and the globe through teleconnections. We explored several model biases in the tropics that might influence the simulation of SST gradient change. We quantified the contributions of internal variability to the observed SST gradient strengthening to understand the causes of past changes and future projections using 32 CMIP6 models and six large ensemble simulations.

The externally forced SST gradient change, represented by the ensemble mean of a single large ensemble, showed a slight change from 1981 to 2010. In contrast, the uncertainty induced by internal variability was significant. For example, under the same external forcing, the ensemble mean of the MPI-ESM showed an average increasing trend of − 0.31 ℃/100 yr, ranging from − 2.00 to 1.61 ℃/100 yr. On the other hand, the ensemble mean of the FGOALS-g3 showed a decreasing trend of 0.38 ℃/100 yr, ranging from − 0.50 to 1.17 ℃/100 yr. The inter-model spread in the MPI-ESM was about twice as large as in the FGOALS-g3. In addition, the magnitude of the ensemble mean trend varied among the six large ensembles, ranging from − 0.31 to 0.67 ℃/100 yr, indicating model uncertainty in response to external forcings, such as sulfate aerosols (Takahashi and Watanabe 2016) or Atlantic warming (McGregor et al. 2014).

Models less biased in easterly trends in the central-eastern Pacific, the subsurface cooling trends along the equator, and the excessive cold tongue exhibit more realistic changes in SST gradients (Additional file 1: Fig. S10). Underestimation of the strengthened trade wind trend in the central equatorial and southeast Pacific, insufficient subsurface cooling in the eastern Pacific, and the excessive westward extension of the climatological cold tongue might have favored a weaker SST gradient trend in the models compared to the observations. Large ensembles could identify the difference in trade wind trends and subsurface cooling effects between S and W models. However, they could not distinguish between S and W models' extreme cold tongue biases. The difference between the S and W models might have been partly due to internal variability, but systematic errors could not be excluded.

With the positive-to-negative phase shift of the IPO, the SST gradient strengthened during 1981–2010, implying that the IPO might be the main factor for the observed SST gradient changes. IPO accounted for ~ 51% (42–62%) of the observed SST gradient strengthening. Based on the S models, our results showed a weakened SST gradient by the mid-twenty-first century. The magnitude of the SST gradient trend varied among the six large ensembles, but they all showed a weakened SST gradient in the coming decades (Additional file 1: Fig. S11). The weakened SST gradient in the projection indicated the slowdown of the Walker circulation, which was consistent with previous studies, in which the Walker circulation was projected to be weakened after adjustments (Huang and Ying 2015; Wu et al. 2021). Our findings suggest that internal variability can influence the SST gradient on decadal timescales. However, the relative importance of internal variability will decline as atmospheric greenhouse gas concentrations continue to increase into the future (Bordbar et al. 2019).