CMIP6 models cannot capture long-term forced changes in the tropical Pacific sea surface temperature gradient

Domain: physics
Published: 2025-12-06
Editorial outcome: minor_revision

Abstract
The observed zonal sea surface temperature gradient in the tropical Pacific has strengthened over the last 150 years, but many CMIP6 models simulate a forced weakening of this gradient over the same period. This has spurred a multi-decade debate over whether models are correctly representing dynamics in the tropical Pacific and whether the observed strengthening is a forced response. We comprehensively assess all observed and modeled gradient trends of 20 years or longer from 1870-2024 across five observational datasets and 14 CMIP6 large ensembles and find that models are not able to match many long-term trends in the observed gradient, especially those that end more recently. Models that are able to match these trends do so through excessive internal variability that compensates for their gradient-weakening forced responses. We additionally find that trends in the observed gradient strengthen at an increasing rate with time, a forced response that is in contrast to the behavior of most models.

Introduction
In the equatorial Pacific Ocean, the sea surface temperature (SST) difference between the West Pacific Warm Pool (WPWP) and the cold tongue in the east plays a particularly outsized role in global climate.

The mean state of this gradient (ΔSST west-east) is established by the easterly trade winds that pile up warm water in the WPWP and deepen the thermocline in the west, while shoaling it in the east, where wind-driven upwelling of cold water causes the formation of the cold tongue. This mean state influences both regional and global climate. Atmospheric heating associated with deep convection over the WPWP contributes to the ascending branch of the Hadley circulation and drives the circumtropical Walker circulation, as well as Rossby waves that propagate poleward to help establish extratropical stationary waves that govern regional climates worldwide. The cold tongue in the east, as a consequence of upwelling, takes up more heat from the atmosphere than any other region on Earth1 and outgasses more CO2 than any other region2. The mean state is, however, not stable and varies on interannual to decadal timescales, with impacts on global and regional precipitation and temperature variability, air-sea CO2 fluxes3, and tropical cyclone activity in the Atlantic and Pacific4,5. It has also been argued that changes in the gradient affect climate sensitivity through low-level cloud feedbacks6, and the global warming rate via changes in ocean heat uptake7.

Observations show a statistically significant strengthening in the gradient since 18808–10. In contrast, many models participating in the Coupled Model Intercomparison Project (CMIP, most recently CMIP6) simulate a forced response that weakens the gradient. Arguments aimed at explaining why rising greenhouse gases (GHGs) should weaken the gradient rest on atmospheric and surface heat flux mechanisms. Knutson and Manabe11 and Betts12 argued that rising tropical tropospheric static stability, a consequence of rising moisture content and latent heat release in deep convection, will weaken overturning. Knutson and Manabe additionally argue that because of the nonlinearity of the Clausius–Clapeyron relation, regions with colder SST (the cold tongue) must warm more than regions with higher SST (the WPWP) to balance reduced surface longwave cooling with enhanced latent heat loss. Both effects would weaken the Walker circulation, further reducing the zonal SST gradient. Vecchi and Soden13, building on Held and Soden14, used climate models to show that, under increased radiative forcing, water vapor content in the lower troposphere increases faster with temperature than precipitation does, also requiring a reduction in Walker circulation strength that translates into a weakened equatorial Pacific zonal SST gradient. In contrast, arguments that support a strengthening gradient response to increasing GHGs9,15,16 focus on the Ocean Dynamical Thermostat mechanism17, which postulates that the shallow thermocline and upwelling in the eastern and central equatorial Pacific should restrict increases in SST in these regions, relative to the WPWP, leading to an increase in the gradient. In either case, the Bjerknes18,19 feedback works to amplify the tendency in the gradient response.

The difference between the observed gradient changes and forced response in models has spurred a debate over whether observed changes reflect internal variability, thus obscuring a forced response, or whether climate models are incorrectly representing forced gradient changes. Some studies have found that models can match observed trends, even when the forced response is a weakening gradient, via internal variability20,21. Other studies10,22,23 have argued that observed trends are at or beyond the limit of what climate models can generate, raising a question about whether the weakening forced response in models is due to model mean-state biases or misrepresentation of dynamics9,24–31. Nevertheless, previous studies assessing relationships between observed and modeled responses in the gradient have been performed over different time periods and include different models, and it is likely that their results are sensitive to these choices.

The apparently contradictory results of earlier studies that have compared observations and models raise the important questions of whether and, more importantly, why the choice of specific time periods and models may influence the ability of models to match observed trends. To comprehensively address these questions, we compare all observed trends of 20 years or longer in ΔSST west-east (1870–2024) to those from 14 CMIP6 large ensembles and consider the time-varying contributions of forced responses and internal variability to the ability of models to match observed trends in ΔSST west-east; collectively, our approach aligns with best practices recently articulated by Simpson et al.32. Because GHG-driven forcing has increased over the historical period33, we additionally analyze the time evolution of gradient trends in models and observations to identify whether their magnitudes are changing over time and are thus consistent with a response to intensifying GHG forcing.

Results
Model ability to capture trends is sensitive to trend choice, with models failing to capture recent trends

To characterize model performance over all time periods, we compare observed trends to the full 337-member grand ensemble comprising 14 CMIP6 large ensembles (Table 1, Fig. 1a) and each large ensemble individually (Fig. 1d). To rigorously address unavoidable observational uncertainties, we include five datasets in this analysis (HadISST, COBE, COBE2, Kaplan, ERSSTv5; Table 2), which collectively represent a variety of data sources, bias corrections, and interpolation techniques. The reliability of SST reconstructions in the early, sparsely sampled periods of these datasets has been validated using a number of techniques, including initializations of ENSO hindcasts34, reconstruction of verifiable SST fields from subsampled data35,36, and use of early reconstructed SST data as a boundary condition for atmospheric general circulation models, including in the 20th Century Reanalysis product37. An earlier study also demonstrates consistent positive trends for the period 1900–2013 in the tropical Pacific zonal SST gradient for five fully interpolated observational datasets from the same lineage as those included in this study, as well as all 100 members of the uninterpolated HadSST3 ensemble that samples a range of bias correction methods and their associated uncertainties22. Our analysis additionally quantifies observational uncertainty by representing each observed trend using a confidence interval derived from a probability distribution estimated from the five SST datasets (see “Methods” section). An observed trend is considered to have been matched by a given model ensemble if the 5–95% range in modeled trend values overlaps the observed trend probability distribution by >5% (see “Methods” section). This is a notably permissive test for models that is more objective and conservative than prior work that has addressed the question of whether observed trends in ΔSST west-east are encompassed by the model range of variability.

While the grand model ensemble captures many shorter observed trends and those trends that fall earlier in the observational interval (Fig. 1b,c), it largely fails to capture longer trends that end in recent decades (as recognized in a previous commentary6, although just for the period after 1950), with observed trend values increasingly approaching the limits of the model range for trends ending in more recent years (Fig. 1c). In particular, while the ensemble captures 84% of trends ending before 2000, it captures only 55% ending after 2000, and only 9% ending after 2020. 79% of observed trends shorter than 100 years fall within the range of trends across the model ensemble; only 50% of trends longer than 100 years are captured, and of those that end after 2000, only 29% are captured. 2% of 100-year or longer trends ending after 2020 are captured. These results remain largely unchanged when comparing observed trends to a bootstrapped grand ensemble that includes 10 randomly selected ensemble members per model so that each model is equally represented; they are also robust across all leave-one-out combinations of the five observational products (Supplementary Figs. 1, 2). While there is a broad spectrum of performance across the individual models – with models like MIROC1 capturing all but 21% of the trends and others, such as INM, failing to capture 64% of trends – the percentage of trends captured across all models decreases after 2000 (Fig. 1d).

Within the context of all calculated trends, our analyses clearly demonstrate that the conclusions from previous studies are sensitive to the chosen analysis period: studies concluding that models are unable to match observations focused on trend periods that end in more recent years or span longer intervals, while studies that have placed observed trends within the range of model variability focus on trends that either end less recently or are shorter (Fig. 1). Studies in the former group additionally include a greater number of models than those in the latter, indicating that the chosen model ensembles also influenced derived conclusions. Assessment over the largest possible ensembles are therefore necessary, given the wide range in model performance with regard to observed trends (Fig. 1d).

Models with oversized internal variability capture almost all trends despite their forced response

Observations and model ensembles include both internal variability and a forced response. Isolating these components is difficult in observations but straightforward in large ensembles. The forced response for each model is calculated as its ensemble mean; averaging across many simulations with common forcing removes the internal variability that is uncorrelated across ensemble members. Subtracting the forced response from all ensemble members yields an internal variability ensemble for each model in which trend ranges are centered about zero. The magnitude of this internal variability will influence the ability of each model to capture observed trends (see “Internal variability magnitude calculations” section), confirmed by the percentage of trends captured by each model (Fig. 2a), which is shown to positively and strongly relate to their estimated magnitude of internal variability (Fig. 2b).

The extent to which observed trends can be explained by internal variability has been characterized by the significance of observed trends relative to null hypotheses based on stationary statistics22 or relative to trends in unforced climate model control runs38. The climate model approach, however, relies on the assumption that models realistically reproduce observed internal variability; this assumption does not hold for all models at all trend lengths, with model internal variability lying above or below the observational range depending on model and trend length (Fig. 2b). It is possible that observed variability over the observational interval may represent an anomalously low period of variability. Nevertheless, proxy-based reconstructions provide insight on this point, and recent estimates of the gradient variance39,40 indicate that this variance (which, in models, positively relates to magnitude of internal variability as reflected in Fig. 2c) in the observational interval is actually high relative to the last 500–1000 years (Supplementary Fig. 3). Consequently, the 20th century appears to be a reasonable period to estimate internal variability in the observational gradient and suggests internal gradient variability in some models may be too large and thus allow more observed gradient trends to be captured for unrealistic reasons. Restricting the ensemble to only those models with internal variability similar to the observations yields an ensemble estimate that does not capture most long-term trend estimates ending after 2000 (Supplementary Fig. 4). This finding is consistent with other studies showing that recently-ending trends are not explained by internal variability and thus may be considered forced9,10,22,23.

The dominant role of ENSO in tropical Pacific SST variability, as well as models’ wide range in ENSO representation41, suggests modeled ENSO characteristics as a possible explanation for the large range of magnitudes in the modeled internal variability. These magnitudes, however, are not entirely explained by the modeled ENSO representation: ENSO variance explains just 45% of variance in gradient trends (r=0.67, p≪0.05; see Supplementary Fig. 5), implying influences from other factors.

For half of the models considered here in, the internal variability ensembles capture at least 5% more of the observed trends than their parent ensembles (compare inset percentages in Figs. 1d and 2a). We classify these models as Group A (Fig. 2 red triangles, Supplementary Table 1). Those models for which internal variability does not capture more trends than the parent ensembles are classified as Group B (Fig. 2, blue triangles, Supplementary Table 1). The implication is that the forced responses in the Group A models are toward a weakening in ΔSST west-east, which is confirmed across almost all timescales in Fig. 3a. While previous studies have highlighted this weakening gradient in models over specific periods10,20–23, the ubiquity of the weakening trends over almost all intervals in some models (e.g., MIROC1 and 2, EC-Earth3, MPI1, GISS) is strikingly in contrast to observations (Fig. 1a). On the other hand, Group B models demonstrate either muted strengthening or muted weakening trends throughout, which is qualitatively closer to the observed trends. Mean gradient trends for each trend length are universally and clearly strengthening in observations, modestly strengthening for almost all models in Group B, and universally and clearly weakening in Group A models (Fig. 3b).

Across this group of large ensembles, the mean forced response for each large ensemble is negatively correlated with the mean magnitude of internal variability (r=−0.83, p≪0.01, Supplementary Fig. 6), indicating that models in each group capture observed trends for different reasons (see “Methods” section). Group A models have oversized internal variability (Fig. 2b) that aids their ability to capture observed gradient trends, even as their forced responses weaken the gradient. In contrast, Group B models capture some observed trends because, despite undersized internal variability, their forced responses are neutral or weakly positive. While Group A large ensembles match, on average, a greater percentage of the observed trends than Group B, this is solely because oversized internal variability overwhelms a forced response that is opposite to the observed gradient strengthening. Both groups, it should be remembered, still struggle to capture the longest trends ending in the most recent years via internal variability alone, or combined with their forced response.

Observations and models exhibit contrasting changes in trends consistent with opposing forced responses

We next examine whether trends in ΔSST west-east are themselves showing systematic changes over time. A positive second-order trend (see “Methods” section and dashed lines in Fig. 4b) over the observational interval indicates that ΔSST west-east is increasing at a growing rate or decreasing at a slowing rate; a negative value indicates that ΔSST west-east is decreasing at a growing rate, or increasing at a slowing rate (Supplementary Fig. 7). Radiative forcing has been increasing at a growing rate over this full analysis interval42. If zonal gradient trends are sensitive to this forcing, forced responses should grow larger over time and thus induce non-zero second-order trends in the gradient.

For both models and observations, significant positive or negative second-order trends in ΔSST west-east, either increasing or decreasing at a growing rate, are therefore interpreted as indicative of a forced response. On the other hand, gradient changes attributable solely to internal variability should yield no long-term significant changes in ΔSST west-east trends. For observations, the second-order trend is positive for almost all trend lengths (with a small range of marginally insignificant trends with lengths between 74 and 78 years, all with 0.05<p-value<0.08), reflecting that ΔSST west-east is increasing at a growing rate (Fig. 4b,c). This finding is robust even if observed trends ending after 1997 are excluded to eliminate the influence of the recent decades of cool tropical Pacific SSTs following the 1997/1998 El Niño43. In contrast, five of the seven Group A models show consistent negative second-order trends across almost all trend lengths, indicating that ΔSST west-east is decreasing at a growing rate. The remaining two models (MPI(LR) and MPI(HR)) exhibit insignificant second-order trends for many trend lengths. Group B models show less consistent behavior, with two models (CanESM5 and ACCESS) exhibiting uniformly negative second-order trends, mirroring the evolution of the forced response observed in most Group A models, while the remaining models yield combinations of positive and negative trends.

This collection of results indicates that observations display increasingly positive trends as radiative forcing has intensified at a growing rate, while most Group A models and two Group B models show increasingly negative trends. Not only does this support the conclusion that many models exhibit forced ΔSST west-east changes that are opposite to those in observations, but it also indicates that ΔSST west-east trends in these models are diverging from observed trends with time. This divergence relates directly to changes in the mean state ΔSST west-east, rather than changes in ENSO. Observations demonstrate positive trends in both skewness (a measure of ENSO nonlinearity in which positive values indicate stronger or more frequent El Niño events relative to La Niña events) and variance of Niño3 SST anomaly (SSTA; see Supplementary Fig. 8) over this period and, hence, observed El Niños are getting stronger and/or more frequent, consistent with Cai et al.44. This alone would act to weaken ΔSST west-east, implying that the observed strengthening must be related to other processes and is more likely a change in the mean state. Group A models, on the other hand, demonstrate increasing ENSO variance and either a reduction or very little change in skewness. This implies an increase in amplitude of both El Niño and La Niña events, which would not be reflected in consistent and growing gradient-weakening trends. Therefore, these forced trends in the models also likely represent a change in the mean state rather than a change in ENSO, though one that is in the opposite direction to observations.

Methods
CMIP6 large ensembles

To compare observed gradient trends to those simulated by CMIP6 large ensembles (also referred to hereafter as ‘models’), this study uses the surface temperature (SST) output from 14 large ensembles (337 total ensemble members). Ensembles were selected on the basis of having ten or more historical runs and at least one SSP245 scenario realization, all with uniform initial conditions, physics, and forcing schemes. The selected models are: CanESM5, MIROC6, MIROC-ES2L, MPI-ESM1-2-LR, MPI-ESM1-2-HR, GISS-E2-1-G, IPSL-CM6A-LR, CNRM-CM6-1, CNRM-ESM2-1, CESM2, INM-CM5-0, UKESM1-0-LL, ACCESS-ESM1-5, and EC-Earth3. To encompass our entire period of interest, we concatenate historical runs with SSP245 scenario runs. For historical runs for which scenario runs with identical variant labels are available (i.e., runs that are continuous from historical to future scenarios), runs are concatenated as such; in cases for which matching pairs cannot be identified, historical runs are randomly paired with one of the scenario runs available from the same model.

Observed SST datasets

We use five observational SST analysis products: the Hadley Center Global Sea Ice and Sea Surface Temperature dataset (HadISSTv1.155, available at https://www.metoffice.gov.uk/hadobs/hadisst/), the National Oceanic and Atmospheric Administration Extended Reconstructed SST V5 (ERSSTv556, available at https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html), COBE Sea Surface Temperature (COBE57, available at https://psl.noaa.gov/data/gridded/data.cobe.html), COBE-SST2 and Sea Ice (COBE235, available at https://psl.noaa.gov/data/gridded/data.cobe2.html), Kaplan Extended SST V3 (Kaplan36, available at https://iridl.ldeo.columbia.edu/SOURCES/.KAPLAN/.EXTENDED/.v3/.ssta/). HadISST incorporates both in situ and satellite data into analysis; COBE2 and ERSSTv5 rely primarily on in situ observations, supplementing with satellite data to reconstruct SST variability in data-sparse regions; COBE and Kaplan rely solely on in situ measurements. All data products use Optimal Interpolation to fill coverage gaps; Kaplan additionally applies Optimal Smoothing. Because each dataset includes errors and uncertainties, the inclusion of all five provides a sample of the observational uncertainty and, therefore, a fairer test of the ability of models to capture observed trends.

SST zonal gradient trends

We define the Pacific zonal gradient as the difference in area-weighted SST anomaly (SSTA) between a previously defined western tropical Pacific box (140°–170°E, 3°S–3°N) and eastern tropical Pacific box (170°–90°W, 3°S–3°N)10. Positive values indicate a relatively warmer western Pacific and a relatively cooler eastern Pacific, while negative values indicate a reduced SSTA difference between the two regions.

We calculate SSTA by removing the mean monthly SST from each grid cell over the entire period defined for this study (1870–2024). This was performed for each observational dataset and each ensemble member of the 14 large ensembles.

On both observational and model gradient timeseries, we perform a linear least-squares regression using monthly data for all 9180 trends longer than 20 years between years 1870 and 2024. Trends were calculated in sliding windows, incremented annually according to the calendar year. Slopes are computed in units of K per decade.

Calculating the combined observed trend

Combined error-weighted average trend values are calculated for each start-end date combination across all five observational datasets as follows58:

x wav = (P 5 w_i x_i) / (P 5 w_i) (1)

where w_i = 1 / σ_i^2 (2)

where x_wav is the weighted average trend, x_i is the trend value from each observational dataset, and σ_i is the standard error associated with the linear trend. Uncertainties in the weighted average trend (‘combined uncertainty’) were calculated as

σ_wav = 1 / sqrt(P 5 w_i) (3)

The standard error is a representation of the uncertainty of the linear least-squares regression fit arising from variance in the underlying data. In this case, this variance includes both the variability of the physical system as well as potential data errors and uncertainties. For each of these datasets, the ratio between the standard error in the slope and the gradient variance is very similar, indicating that the proportion of the standard slope error attributable to physical variability versus data errors is likely comparable across all datasets. Furthermore, the standard errors for each trend across all datasets are similar, justifying the choice of the slope standard error as a metric of observational uncertainty.

Identification of trends captured by models

To identify the observed trends that are ‘captured’ by models, the following test is applied: if a trend value within the 5–95% range of the model realizations has a greater than 5% likelihood of falling within the range of observational trend uncertainty, that trend is considered captured by the model. The uncertainty range of the combined observed trend is defined by a probability density function with a mean of the weighted average observed trend and a standard deviation of the combined uncertainty. Traditional parametric ‘difference of means’ tests are not applicable in this context because assessment of each trend requires the comparison of a parametric (observational) and a non-parametric (model) distribution. As such, by treating the model distribution, with its discrete points, as a known range and the observations as a probability distribution, this test enables the calculation of the likelihood that these two distributions overlap significantly (i.e., p-value > 5%).

Isolation of model forced response and internal variability

For each model, the forced response in the gradient is calculated as the mean gradient value across all ensemble members because averaging across many simulations with common forcing removes the internal variability that is uncorrelated across ensemble members. Trends in this mean gradient are then calculated as described in the “SST zonal gradient trends” section and are considered to be the forced trends in the gradient. To calculate the internal variability of each model, the forced response in the gradient is removed from each ensemble member’s gradient timeseries. Trends in this gradient are then calculated for each ensemble member to create an internal variability ensemble for each model.

Internal variability magnitude calculations

Calculating the true magnitude of internal variability in observations is complicated by the fact that only one realization of the real world is available. As shown in one of the results of this study, no long-term changes in trends appear to be consistent with a forced response. As such, by removing the long-term linear trend in observed trend values for each trend length, we aim to remove the influence of external forcing and, to the degree possible, isolate the internal variability in the observational gradient. We then estimate the magnitude of this internal variability by calculating the standard deviation of the observational internal variability trends for each trend length. While it is possible to isolate the internal variability in models by removing the forced response, in order to directly compare a measure of internal variability in models to that in observations, we instead apply the same metric used to calculate the observed internal variability to each ensemble member in every model. For each trend length, we linearly detrend each ensemble member’s trend values and then calculate the combined detrended standard deviation across the full ensemble for each model. The use of intermember standard deviation to quantify internal variability magnitude has also been used in other studies (e.g., Olonscheck et al.20; Wang et al.59). The additional step of detrending at each trend length allows for a comparison to be made between modeled and observed internal variability. Note that this method for estimating the magnitude of internal variability is specific to Fig. 2b.

Correlation between forced response and internal variability

The Pearson correlation coefficient between models’ forced response and internal variability is calculated based on the mean forced response and mean magnitude of internal variability for each ensemble. Mean forced response is calculated as the average of all 9180 trends, while mean internal variability magnitude is calculated as the average internal variability magnitude (as described in the “Internal variability magnitude calculations” section) across all trend lengths.

Second-order linear trends

To compare long-term changes in trends between models and observations, we calculate second-order trends. For observations, all trends are calculated for a given trend length using a sliding analysis window; the weighted least-squares regression of these trends is then calculated, with weights of 1/σ_wav^2. Only significant trends are included (p-value < 0.05). For forced responses in the model ensembles, a least-squares regression is applied to all trends for each trend length for each model, with only significant values included.

ENSO analysis

Niño3 SSTA is calculated as the area-weighted mean SSTA in the Niño3 region (150°W–90°W, 5°S–5°N). Skewness and standard deviation are calculated on the detrended annual mean model and observational data. Annual mean values are calculated with the year being defined from May to April to better represent the ENSO cycle. Model SSTA is detrended by subtracting the ensemble mean SSTA, while observations are detrended by subtracting the 60-year low-pass filtered SSTA. For models, ensemble mean values are shown, where the mean across the ensemble is calculated after skewness and standard deviation have been calculated for each ensemble member.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This new study is a nice extension of the senior authors' prior work on the same topic that the observed tropical Pacific SST trend pattern is outside or at the margins of the internal variability of the climate models. In this updated work, the authors show that this real world - vs - model simulation gap is accelerating. They further show that the models that are relatively successful at simulating the SST gradient do so for a wrong reason. These are important results which deserve to be published.\nThe manuscript is well written in that the ideas are clearly presented. The analyses are excellent. I do appreciate seeing the 'location' of the prior work by various authors in Fig. 1. This figure is particularly informative.", "reviewer_2": "This study compares tropical Pacific zonal SST gradient trends between CMIP6 historical simulations and observations since 1870. Observations consistently show long-term strengthening, whereas models tend to fail to reproduce this, with discrepancies especially large after 2000. The authors classify models into two groups: Group A matches observed trends by overestimating internal variability but fails to capture the forced response, while Group B underestimates internal variability but better represents the forced tendency. They also show that models with a weakening gradient tend to have larger internal variability. Finally, second-order linear trends analysis reveals that the observed strengthening of SST gradient accelerates over time, consistent with increasing radiative forcing but not reproduced by models, suggesting a major model–observation discrepancy in the long-term forced response. The first part of the results extends the trend analysis of Rugenstein et al. (2023), which went back to 1950, further back to 1870, and the findings are consistent with theirs. The inverse correlation they showed between the signal of zonal gradient strengthening in each model and the magnitude of internal variability is intriguing, and I would like to see a more detailed physical explanation of this relationship. Furthermore, the final section of the results, which shows that the observed trend in the zonal gradient has been strengthening at an accelerating rate over time, presents a forced response that contrasts with the behavior of most models. This suggests that the discrepancy between models and observations will persist or even worsen in the future. This finding is particularly interesting because it contradicts the optimistic view that mechanisms such as the ocean thermostat—which strengthen the SST gradient—are fast responses to forcing and that, eventually, both models and observations would converge toward an El Niño–like SST change. I want authors to address several comments before publication.\nMajor Comments:\n1. It has been pointed out that tropical SST observations before 1950 are sparse and thus uncertain. For example, Chung et al. (2025) noted that instrumental observations prior to 1950 have large uncertainties, and therefore used proxies in a paleoclimate context to show that, over the very long term, the SST gradient decreases in response to rising CO₂. How would you respond to the concern that SST data prior to 1950 may not be suitable for analyzing tropical SST gradients?\n2. L166-176: Can you explain why models with larger internal variability tend to exhibit a weakening SST gradient? Rugenstein et al. (2023) demonstrated a correlation between internal variability of the SST gradient and ECS. They hypothesized that it possibly due to an incorrect balance of cloud feedbacks driven by changing inversion strength versus surface warming. Can you provide a physical explanation for the correlation you found?\n3. L166-176: Does the anti-correlation between the strength of internal variability and the forced response (r=-0.87) hold across all trend lengths, or is it more robust for long-term trends? Placing a scatterplot in the Supplementary that shows this anticorrelation, along with an explanation of how the correlation coefficient is derived, would make it clearer what analysis was performed.\nSpecific comments:\n1. A 100-year trend line should be added in Fig. 1b.\n2. How is the standard deviation of observations in Fig. 2b calculated?\n3. The method for calculating the second-order linear trends is unclear. Is my understanding correct that you fixed the trend length, slid the analysis window, and then calculated the trend of those trends?\n4. L389–405: You state that internal variability in the observations is extracted using detrending, and in that explanation you also mention applying the same operation to the models. Please specify in the same paragraph for which analyses or figures (only Fig. 2c?) this operation was applied to the models. Otherwise it becomes confusing with the explanation that model internal variability is defined as the deviation from the ensemble mean.\n5. L95: “with observed trend values increasingly approaching the limits of the model range as time progresses” — I don’t understand this. What does “as time progresses” mean?\n6. In the caption of Fig. S4, you should clarify how you selected the models whose internal variability is consistent with the observations (does this correspond to the gray shading in Fig. 2c?).", "reviewer_3": "It is now well established that zonal SST gradients in the Equatorial Pacific have strengthened in recent decades. But this trend is not captured by most climate models. This discrepancy between the observed and modeled response to global warming of zonal SST trends in the Equatorial Pacific, a region known to remotely affect various aspects of Earth’s climate, has been a key issue in recent years, since it potentially points to a critical deficiency of climate models. Here the authors compare trends in the zonal temperature gradient in the Equatorial Pacific between ensembles of 14 CMIP6 models and five observational datasets. They find diverging responses in models and observations (i.e., a growing discrepancy) and that larger than observed variability in said trends in CMIP6 models, can mask their overall inability to capture the observed forced response. The paper is very well written, timely, and the results are well supported. I only have one significant comment:\n- As mentioned in the discussion, there is some consistency between Groups A and B and between the way CMIP6 models capture ENSO variability. The changes in the zonal SST trends can result from changes in the distribution of ENSO flavors, for example, toward more prevalent, longer, or stronger Niña episodes (and potentially similarly, for the negative phase of the PDO mode). This can affect several aspects of the results. First, Niño events tend to be stronger than Niña events, which suggests that the trend variance about the mean zonal SST gradients is not normally distributed. Second, are the models diverging from observations in the manner in which they represent the distribution of ENSO flavors or in the mean zonal SST gradient? Third, are the variations in zonal SST gradient trend variability (Figure 2) linked to the way models represent ENSO (e.g., the zonal SST gradient is likely strongly correlated with ENSO indices; are the statistics of the zonal SST gradient index similar to those of ENSO related indices?). These seem to be an important distinctions, that can be addressed within the scope of the present analysis.\nI would therefore recommend accepting the work, after addressing the above comment and some additional minor comments listed below.\nComments by line number\n145-147 (I hope I’m not misinterpreting Supplementary Figure 4.) Would this statement be true if the heat maps are constructed such that the end year is fixed at 2024 and the begin year is walked back? In this case, would long term trends that begin earlier than ~1900 not be captured? And further, if there is such a difference, can we rule out the possibility that part of the discrepancy is due to the larger uncertainty in earlier years (i.e., observed uncertainty likely reduces with more modern data)?\n229-231 Note however, that Adam et al. (2023) find that mean-state evaporation biases alone cannot account for the biased equatorial warming trend.\nFigure 3a Colors are too pale"}, "author_response": "Response to Reviewers Document for Nature Communications Submission NCOMMS-25-41165 by Hannah Byrne, Richard Seager and Jason E. Smerdon\nWe thank the Reviewers for their insightful and constructive feedback.\nThis file provides a complete documentation of the changes that have been made in response to each individual Reviewer comment, with formatting indicating the following: Reviewer comments are shown in blue; responses are in plain text, with quotations from the revised manuscript in italics; line numbers in our responses refer to locations in the revised manuscript. New figures relevant to responses are included at the end of this document, as well as in the revised Supplemental Information.\nReviewer #1 Evaluations:\nThis new study is a nice extension of the senior authors' prior work on the same topic that the observed tropical Pacific SST trend pattern is outside or at the margins of the internal variability of the climate models. In this updated work, the authors show that this real world - vs - model simulation gap is accelerating. They further show that the models that are relatively successful at simulating the SST gradient do so for a wrong reason. These are important results which deserve to be published.\nThe manuscript is well written in that the ideas are clearly presented. The analyses are excellent. I do appreciate seeing the 'location' of the prior work by various authors in Fig. 1. This figure is particularly informative.\nWe thank the Reviewer for their endorsement of this work and are pleased that they view the results as valuable and informative.\nReviewer #2 Evaluations:\n[...Reviewer 2's comments repeated...]\nWe thank the Reviewer for this positive assessment.\nMajor comments\n1. It has been pointed out that tropical SST observations before 1950 are sparse and thus uncertain. For example, Chung et al. (2025) noted that instrumental observations prior to 1950 have large uncertainties, and therefore used proxies in a paleoclimate context to show that, over the very long term, the SST gradient decreases in response to rising CO2. How would you respond to the concern that SST data prior to 1950 may not be suitable for analyzing tropical SST gradients?\nWe agree that observational uncertainties are an important consideration, which is the rationale for both the inclusion of five observational datasets in our study, as well as the treatment of observed trends as a probability distribution rather than a singular value. We additionally demonstrate that our conclusions are insensitive to the exclusion of any observed dataset (see Fig. S2). Despite this multifaceted approach to uncertainty assessments, it was only briefly mentioned in the body of the original manuscript. We therefore better highlight this rationale and methodology in the revised manuscript.\nAll of the employed SST reconstructions have been demonstrated to be reliable for use in earlier periods through a variety of studies, which we now discuss and include in the revised manuscript. In one previous study (Chen et al., 2004), three of the datasets (Kaplan, HadISST and ERSST) were evaluated by conducting ENSO hindcasts from 1856 onwards (with a 6-month lead time) using only reconstructed SST data for initialization of the Zebiak-Cane (ZC) geophysical model. In the forecast system, observed prior SSTs are used to force the atmosphere component of the ZC model and the modeled winds are used to drive the ocean model up to the forecast start time. Remarkably, the SST-initialized forecasts were able to hindcast, at a 6-month lead, almost all of the El Niño-events in the observed record from 1856 on. Because the distribution of ocean heat content is key to ENSO prediction, the SST data must be sufficiently accurate to lead to modeled winds and upper ocean heat content that enable predictability of subsequent ENSO evolution. This is an extremely demanding test of the realism of the SST data. While it is true the test does not necessarily assess the ability of the SST products to capture the zonal SST gradient, it does provide considerable confidence that the SST products have a high degree of realism over their entire period of record.\nThe extensive validation carried out for each of the individual SST datasets should also not be dismissed. The Kaplan dataset, one of three used in the Chen et al. (2004) study, underwent rigorous validation as described in Kaplan et al. (1998). In one test, observational data from a more recent period of higher data availability were subsampled to match the sparseness of data available in earlier years (e.g., 1868, 1877, 1918) and then used to reconstruct the complete SST field during recent decades. The striking consistency between these reconstructions and the actual SST data demonstrates the effectiveness of their reconstruction techniques even in periods of data sparsity. The COBE2 dataset underwent a similar validation procedure, through which the authors conclude that Niño 3 SST data are verifiably reliable from the 1880s onwards, while noting that the data remain very similar to HadISST even before this period (Hirahara et al., 2014). HadISST itself was independently validated through its use as the boundary condition for the production of the 20th Century Reanalysis product (Compo et al., 2011), which has successfully reproduced key atmospheric variables from 1871 onwards. The reliability of HadISST, verified through both the Chen et al. (2004) study and its use in 20th Century Reanalysis, lends credibility to COBE, which reports on the similarity between its reconstruction and HadISST as a means of validation (Ishii et al., 2005). In another study, Coats and Karnauskas (2017) find positive zonal SST gradient trends from 1900 to 2013 in five fully interpolated observational datasets, all of which are within the same lineages as the datasets that we use (ERSSTv2, ERSSTv3b, ERSSTv4, Kaplan SSTv2 and COBE2). They additionally find positive trends for this period in all 100 ensemble members of the uninterpolated HadSST3 ensemble (Kennedy et al., 2011a, 2011b), which samples a range of bias correction methods and associated uncertainties, demonstrating that the long-term strengthening of the gradient is unlikely to be an artefact of bias correction and interpolation choices. We also note that nothing in our trend analysis (or that of others) suggests unreliable data before the mid 20th century: the long-term gradient strengthening in the more data-rich post 1950s period is consistent with that found using earlier start periods.\nThese arguments have been added to the revised manuscript as follows:\nLine 90: To rigorously address unavoidable observational uncertainties, we include five datasets in this analysis (HadISST, COBE, COBE2, Kaplan, ERSSTv5), which collectively represent a variety of data sources, bias corrections, and interpolation techniques. The reliability of SST reconstructions in the early, sparsely sampled periods of these datasets has been validated using a number of techniques including initializations of ENSO hindcasts (Chen et al., 2004), reconstruction of verifiable SST fields from subsampled data (Hirahara et al., 2014; Kaplan et al., 1998) and use of early reconstructed SST data as a boundary condition for atmospheric general circulation models, including in the 20th Century Reanalysis product (Compo et al., 2011). An earlier study also demonstrates consistent positive trends for the period 1900-2013 in the tropical Pacific zonal SST gradient for five fully interpolated observational datasets from the same lineage as those included in this study, as well as all 100 members of the uninterpolated HadSST3 ensemble that samples a range of bias correction methods and their associated uncertainties (Coats & Karnauskas, 2017). Our analysis additionally quantifies observational uncertainty by representing each observed trend using a confidence interval derived from a probability distribution estimated from the five SST datasets (see Methods). An observed trend is considered to have been matched by a given model ensemble if the 5-95% range in modeled trend values overlaps the observed trend probability distribution by >5% (see Methods). This is a notably permissive test for models that is more objective and conservative than prior work that has addressed the question of whether observed trends in ΔSSTwest-east are encompassed by the model range of variability.\nBottom line: The inclusion of this new discussion on validation of SST reconstructions over early periods and the improved articulation of our uncertainty assessment across the five SST datasets now serves to better inform the reader about the data uncertainties and instill confidence that they are robust, as assessed in our study and in extensive prior work.\n2. L166-176: Can you explain why models with larger internal variability tend to exhibit a weakening SST gradient? Rugenstein et al. (2023) demonstrated a correlation between internal variability of the SST gradient and ECS. They hypothesized that it possibly due to an incorrect balance of cloud feedbacks driven by changing inversion strength versus surface warming. Can you provide a physical explanation for the correlation you found?\nWhile a full model-by-model physical diagnosis of this relationship is outside of the scope of our study, we have conducted additional analyses into ENSO characteristics in both models and observations to interrogate this relationship. These analyses paint a more complete picture of the potential nonlinear ENSO-mean state rectification explanation that we briefly covered in the Discussion section of the original manuscript (Hayashi et al., 2020; Kohyama et al., 2017). The results of this analysis have added valuable support for our initial arguments. We address the subject in the revised manuscript as follows:\nLine 294: It is possible that this nonlinear rectification of the mean state partially explains the significant anti-correlation between forced response trends and the magnitude of internal variability in models. While this effect likely plays no role in Group B models, for which Niño 3 SSTA skewness and thus ENSO nonlinearity is underestimated (Supplementary Fig. 8), for Group A models, high and increasing ENSO variance and more realistic ENSO nonlinearity are consistent with, respectively, a higher magnitude of internal variability and a ΔSSTwest-east weakening forced response through nonlinear rectification onto the mean state. In agreement with this hypothesis, we note that this anti-correlation is not significant for either individual model group (Supplementary Fig. 6), but is larger for Group A than Group B, implying a more direct connection between the magnitude of internal variability and the sign of forced response for these models. Given that the correlation is strongest and significant only for the full model group, it is likely that this nonlinear rectification only partially explains the relationship between forced responses and the magnitudes of internal variability. Another potential contributing effect is the role that simulated cloud feedbacks play in SST variability (and change) within models (Rugenstein et al., 2023). We additionally acknowledge the lack of consensus over whether changes in ENSO variance rectify the mean state, or vice versa (Fedorov & Philander, 2000; Jin et al., 2003; Sun & Zhang, 2006).\n3. L166-176: Does the anti-correlation between the strength of internal variability and the forced response (r=-0.87) hold across all trend lengths, or is it more robust for long-term trends? Placing a scatterplot in the Supplementary that shows this anticorrelation, along with an explanation of how the correlation coefficient is derived, would make it clearer what analysis was performed.\nThis anti-correlation does hold across all trend lengths, as shown in the new Fig. S6. While there was an explanation for how this coefficient is derived in the Methods section, it was not previously referenced in the body of the manuscript; Line 195 of the revised manuscript now points the reader to the relevant section of the Methods.\nSpecific comments\n1. A 100-year trend line should be added in Fig. 1b.\nDone.\n2. How is the standard deviation of observations in Fig. 2b calculated?\nFor Fig. 2b, the standard deviation of the observations is calculated as described in Methods (under “Internal variability magnitude calculations”). To ensure this is clear for the reader, we now include a direct reference to the Methods section in Line 146 of the revised manuscript.\n3. The method for calculating the second-order linear trends is unclear. Is my understanding correct that you fixed the trend length, slid the analysis window, and then calculated the trend of those trends?\nThis is the correct interpretation of how we calculate the second-order linear trend. For clarity, we have updated the “Second order linear trends” section under Methods as follows:\nLine 468: For observations, all trends are calculated for a given trend length using a sliding analysis window; the weighted least-squares regression of these trends is then calculated, with weights of 1 / σ2wav.\n4. L389–405: You state that internal variability in the observations is extracted using detrending, and in that explanation you also mention applying the same operation to the models. Please specify in the same paragraph for which analyses or figures (only Fig. 2c?) this operation was applied to the models. Otherwise it becomes confusing with the explanation that model internal variability is defined as the deviation from the ensemble mean.\nWe have amended this section for clarity as follows:\nLine 456: Note that this method for estimating the magnitude of internal variability is specific to Figure 2b.\n5. L95: “with observed trend values increasingly approaching the limits of the model range as time progresses” — I don’t understand this. What does “as time progresses” mean?\nWe meant that models increasingly struggle to match observed trend values for trends that end in more recent periods. We have clarified the revised manuscript as follows:\nLine 113: …with observed trend values increasingly approaching the limits of the model range for trends ending in more recent years.\n6. In the caption of Fig. S4, you should clarify how you selected the models whose internal variability is consistent with the observations (does this correspond to the gray shading in Fig. 2c?).\nThis does not relate to the gray shading in Fig. 2c. We have clarified this in the Supplemental Information as follows:\nFig. S4: Hatching indicates observed trends captured by the two models (EC-Earth3 and MPI-ESM1-2-HR) with internal variability that is most similar to observations, thus removing the effects of models with unrealistically high internal variability. This selection was performed by ranking all models and observations according to magnitude of internal variability for each trend length. Average rank across all trend lengths was calculated and the two models selected are those immediately above and below the rank of the observations.\nReviewer #3 Evaluations:\n[...Reviewer 3's comments repeated...]\nWe thank the Reviewer for these positive remarks.\nMajor comments\n- As mentioned in the discussion, there is some consistency between Groups A and B and between the way CMIP6 models capture ENSO variability. The changes in the zonal SST trends can result from changes in the distribution of ENSO flavors, for example, toward more prevalent, longer, or stronger Niña episodes (and potentially similarly, for the negative phase of the PDO mode). This can affect several aspects of the results.\n1) First, Niño events tend to be stronger than Niña events, which suggests that the trend variance about the mean zonal SST gradients is not normally distributed.\nWe have extended our analysis to draw out some of the distinctions between changes in ENSO and the mean state by considering how ENSO flavors have changed over the analysis period. As the Reviewer suggests, the nonlinearity in ENSO amplitude with larger El Niño events and smaller La Niña events would create a positive skewness in the distribution of SSTA in the Eastern Equatorial Pacific (EEP). Increases in this skewness would be reflected in negative skewness of ΔSSTwest-east trends (because relatively larger SSTAs in the EEP are associated with a weakened ΔSSTwest-east). This suggests that if models demonstrated higher increases in ENSO skewness than observations, the behavior might explain excessive model ΔSSTwest-east weakening relative to observations. In the real world, skewness has been increasing, which we would expect to be reducing the gradient through nonlinear rectification; however, the observed gradient has been strengthening, reflecting more of a direct mean state change. All but two models, on the other hand, underestimate Niño 3 SSTA skewness relative to observations as noted in previous studies (Hayashi et al., 2020; Kohyama et al., 2017), and for the model group as a whole there is no clear increase in skewness. Hence, model skewness and changes in skewness over time cannot explain the gradient-weakening trends in models or, therefore, resolve the observation-model discrepancy.\nWe summarize this analysis in a new figure Fig. S8, which will replace the similar but limited original Fig. S6.\n2) Second, are the models diverging from observations in the manner in which they represent the distribution of ENSO flavors or in the mean zonal SST gradient?\nWe find that observations and models do show different changes in ENSO characteristics over this period, but those changes do not appear to explain the divergence in ΔSSTwest-east trends that we highlight in the paper and rather exemplify another area in which models and observations differ. We have amended the manuscript as follows:\nLine 236: This divergence relates directly to changes in the mean state ΔSSTwest-east, rather than changes in ENSO. Observations demonstrate positive trends in both skewness (a measure of ENSO nonlinearity in which positive values indicate stronger or more frequent El Niño events relative to La Niña events) and variance of Niño 3 SST anomaly (SSTA; see Supplementary Fig. 8) over this period and, hence, observed El Niños are getting stronger and/or more frequent, consistent with Cai et al., (2021). This alone would act to weaken ΔSSTwest-east, implying that the observed strengthening must be related to other processes and is more likely a change in the mean state. Group A models, on the other hand, demonstrate increasing ENSO variance and either a reduction, or very little change in skewness. This implies an increase in amplitude of both El Niño and La Niña events, which would not be reflected in consistent and growing gradient-weakening trends. Therefore, these forced trends in the models also likely represent a change in the mean state rather than a change in ENSO, though one that is in the opposite direction to observations.\n3) Third, are the variations in zonal SST gradient trend variability (Figure 2) linked to the way models represent ENSO (e.g., the zonal SST gradient is likely strongly correlated with ENSO indices; are the statistics of the zonal SST gradient index similar to those of ENSO related indices?). These seem to be an important distinctions, that can be addressed within the scope of the present analysis.\nWhile models with the highest ENSO variance also have the highest magnitude of internal variability, the Pearson correlation between model gradient trend variance (Fig. 2b) and Niño 3 SSTA variance (Fig. S8b) is 0.67 (Fig. S5), implying that less than half of the magnitude of internal variability of ΔSSTwest-east trends is explained by ENSO magnitude. We have amended the manuscript as follows:\nLine 170: The dominant role of ENSO in tropical Pacific SST variability, as well as models’ wide range in ENSO representation (Bayr et al., 2019), suggests modeled characteristics as a possible explanation for the large range of magnitudes in the modeled internal variability. These magnitudes, however, are not entirely explained by the modeled ENSO representation: ENSO variance explains just 45% of variance in gradient trends (r=0.67, p<<0.05; see Supplementary Fig. 5), implying influences from other factors.\nFig S5: This new figure is now included in the supplemental material, which shows standard deviation of trends in ΔSSTwest-east versus standard deviation of Niño 3 SSTA for all models.\nSpecific comments\nThe revised manuscript now reads:\nLine 28: Atmospheric heating associated with deep convection over the WPWP contributes to the ascending branch of the Hadley circulation and drives the circumtropical Walker circulation, as well as Rossby waves that propagate poleward to help establish extratropical stationary waves that govern regional climates worldwide.\n145-147 (I hope I’m not misinterpreting Supplementary Figure 4.) Would this statement be true if the heat maps are constructed such that the end year is fixed at 2024 and the begin year is walked back? In this case, would long term trends that begin earlier than ~1900 not be captured? And further, if there is such a difference, can we rule out the possibility that part of the discrepancy is due to the larger uncertainty in earlier years (i.e., observed uncertainty likely reduces with more modern data)?\nWe meant that sufficiently long trends (e.g., longer than 100 years) with end dates after 2000 are not captured. To clarify this point, we have added a 100-year trend reference line to Fig. S4 and have amended the previous statement as follows:\nLine 164: Restricting the ensemble to only those models with ΔSSTwest-east internal variability similar to the observations yields an ensemble estimate that still does not capture most of the long-term trend estimates ending after 2000 (Supplementary Fig. 4).\nReviewer #2 also raised a related question about observational uncertainty. Please consult our response to Reviewer #2’s Major comment 1 where we explain our approach to incorporating this uncertainty and our confidence in the early data.\n229-231 Note however, that Adam et al. (2023) find that mean-state evaporation biases alone cannot account for the biased equatorial warming trend.\nThis line has been updated as follows:\nLine 269: For example, Seager et al. (2019) attribute model biases to excessively high relative humidity and wind speeds that are too low over the eastern Pacific cold tongue, a feature related to the double ITCZ bias in models, although further studies have shown that thermodynamic effects alone cannot account for this warming trend (Adam et al., 2023).\nFigure 3a Colors are too pale\nThe figure has been revised.\nAdditional comments from Authors:\nFor the sake of consistency with other observational datasets, we have removed the climatology from ERSSTv5 (which was accidentally retained in the analysis in the initial submission). The effects of this are negligible and do not affect any of our conclusions. We additionally made a modest update to the method for calculating the correlation between trends in the gradient forced response and internal variability magnitude in order to weight each forced response trend equally, which changed the r value from -0.82 to -0.83. This too does not affect any conclusions of the study."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have adequately addressed my comments. I recommend publication in its current form.", "reviewer_3": "The authors have thoroughly addressed my comments. I have no further comments on the revised manuscript, and I recommend the paper for publication."}, "author_response": ""}]

Decision letter