Historical changes in wind-driven ocean circulation drive pattern of Pacific warming

The tropical Pacific warming pattern since the 1950s exhibits two warming centers in the western Pacific (WP) and eastern Pacific (EP), encompassing an equatorial central Pacific (CP) cooling and a hemispheric asymmetry in the subtropical EP. The underlying mechanisms of this warming pattern remain debated. Here, we conduct ocean heat decompositions of two coupled model large ensembles to unfold the role of wind-driven ocean circulation. When wind changes are suppressed, historical radiative forcing induces a subtropical northeastern Pacific warming, thus causing a hemispheric asymmetry that extends toward the tropical WP. The tropical EP warming is instead induced by the cross-equatorial winds associated with the hemispheric asymmetry, and its driving mechanism is southward warm Ekman advection due to the off-equatorial westerly wind anomalies around 5°N, not vertical thermocline adjustment. Climate models fail to capture the observed CP cooling, suggesting an urgent need to better simulate equatorial oceanic processes and thermal structures.

The tropical Pacific warming pattern since the 1950s exhibits two warming centers in the western Pacific (WP) and eastern Pacific (EP), encompassing an equatorial central Pacific (CP) cooling and a hemispheric asymmetry in the subtropical EP.The underlying mechanisms of this warming pattern remain debated.Here, we conduct ocean heat decompositions of two coupled model large ensembles to unfold the role of wind-driven ocean circulation.When wind changes are suppressed, historical radiative forcing induces a subtropical northeastern Pacific warming, thus causing a hemispheric asymmetry that extends toward the tropical WP.The tropical EP warming is instead induced by the cross-equatorial winds associated with the hemispheric asymmetry, and its driving mechanism is southward warm Ekman advection due to the offequatorial westerly wind anomalies around 5°N, not vertical thermocline adjustment.Climate models fail to capture the observed CP cooling, suggesting an urgent need to better simulate equatorial oceanic processes and thermal structures.
The tropical Pacific climate exhibits a prominent zonal asymmetry, consisting of a WP warm pool overlaid with vigorous convective rainfall and an EP cold tongue cooled by strong oceanic upwelling 1 .The pattern of tropical Pacific warming under climate change is crucial to climate sensitivity 2,3 , the global hydrological cycle [4][5][6] , and the El Niño-Southern Oscillation [7][8][9][10][11] .Nevertheless, how the tropical Pacific sea surface temperature (SST) has been changing in the past century and will change in the future has been under hot debate for decades (refer to the recent review by ref. 12).Climate model simulations under global warming often exhibit an equatorially enhanced warming in the eastern Pacific and therefore a weakened zonal SST gradient mimicking an El Niño-like state [13][14][15][16][17][18][19][20] .In contrast, observed SST trends in the 20 th century show a narrow band of equatorial CP cooling, sometimes referred to as a La Niña-like state [21][22][23] .Such model-observation discrepancies persisted through the beginning of the 21 st century 12,24,25 .The latest tropical Pacific SST trends since 1958 exhibit spatial structures that are more complex than a typical El Niño-or a La Niña-like pattern (Fig. 1a; also see in ref. 25).In particular, the tropical Pacific warming pattern consists of a warm-cold-warm tripolar structure along the equator and hemispheric asymmetry in the subtropical EP.
The formation mechanisms of the recent trend pattern remain elusive given the co-existence of various physical processes.From a thermodynamic perspective, in response to increased greenhouse gases, the WP warm pool will warm less than the cold EP due to its strong background evaporative cooling that is relatively more effective in offsetting the trapped longwave radiation 13,26 .This thermodynamic effect is expected to lead to an El Niño-like warming pattern as seen in slab ocean-atmosphere coupled models 19 .When the effect of ocean circulation is considered, the climatological equatorial upwelling in the CP-WP, where the thermocline is sufficiently shallow, can act to suppress the local surface warming, while this effect is rather weak in the WP 27 .Ocean circulation changes due to surface wind changes may further contribute to the warming pattern.For example, a strengthened Walker Circulation and an intensified equatorial easterly surface wind may enhance the zonal SST gradient, further amplified by the Bjerknes feedback 28 .Some theoretical studies suggest an enhanced Walker Circulation and therefore a strengthened zonal SST gradient due to the contraction of cloud area over the western Pacific 29,30 , while others imply the opposite due to energetic constraints 13,31 or low cloud-SST feedbacks 32 .
What role wind-driven ocean circulation changes play in shaping the historical SST trend pattern remains poorly understood as they are often hard to fully untangle from other mechanisms.To address this question, here we compare two Community Earth System Model version 2 (CESM2) large ensembles, first introduced in ref. 33 that allow us to isolate the impacts of externally forced changes in wind-driven ocean circulations on the tropical Pacific warming pattern over the historical period.As we will show later, this unique experimental set-up and our ocean heat decomposition analysis can lead us to identify a dynamic mechanism critical to the historical tropical Pacific warming pattern that has been overlooked in the existing literature.

Complex tropical pattern warming pattern since 1958
Firstly, we present the tropical Pacific SST trends during 1958-2014, following ref. 25, averaged across five observational datasets (see Methods; Fig. 1a; Supplementary Fig. 1).The observed warming pattern is complex, characterized by two warming centers in the WP and EP, respectively, with a distinctive cooling in the CP extending towards the east along the narrow equatorial band.Off the equator, the WP exhibits enhanced warming in both hemispheres, whereas the EP exhibits a warm-north-cold-south hemispheric asymmetry.This warming pattern results from a combination of externally forced responses (i.e.greenhouse gases, aerosols, etc.) and internal variability 1,34,35 , making it hard to identify its formation mechanisms based on observations alone.
To better separate externally forced responses from internal variability, we analyze the CESM2 fully-coupled (FC), large ensemble historical simulations (see Methods).The simulations are integrated from 1850 to 2014, and the period from 1958-2014 is analyzed in our study to compare with observations.The ensemble-mean, FC CESM2simulated SST trends can reasonably capture some of the observed warming pattern with a spatial correlation of 0.4.In particular, the observed enhanced warming in the WP and EP along the equator can be reproduced by the ensemble-mean result, suggesting an external response, although the observed CP cooling is absent in the model simulations (Fig. 1c).For individual ensemble members, the spatial correlation of the tropical Pacific SST trend pattern with observations ranges from 0.44 to −0.60.This large inter-member spread suggests that internal variability can at least explain a fraction of the modelobservation discrepancies, among other factors like historical radiative forcing uncertainties and climate model biases 25 .A similar tripolar pattern is also evident in the multi-model mean historical SST trends from the Coupled Model Intercomparison Project Phase 6 (CMIP6) models, although the CP cooling is again absent (Fig. 1b).In the following analysis, we will focus particularly on the mechanisms of the enhanced WP and EP warming, and discuss the absence of CP cooling in climate models, both CESM2-LE and CMIP6, in the final discussion.

Mechanically decoupled simulations
To investigate the impact of wind-driven ocean circulation changes, we analyze another set of large ensemble members, mechanically decoupled (MD) historical simulations with CESM2, in which ocean surface wind stress is fixed to its seasonally varying pre-industrial climatological state 33 (see Methods) using the wind stress overriding technique [36][37][38][39][40] .Through this approach, the Bjerknes feedback is cut off, and therefore the effects of ENSO variability and Walker Circulation adjustment on SST trends are both eliminated 41 .We find that the MD ensemble-mean SST trends during 1958-2014 exhibit an enhanced warming spanning from the coast of Mexico to the tropical WP, resulting in a hemispheric asymmetry and an enhanced zonal SST gradient along the equator (Fig. 1d).In general, the MD ensemble mean trend agrees less well with observations than the FC ensemble mean trend with a spatial correlation of 0.25.The difference between the FC and MD represents the SST changes caused by externally-forced winddriven ocean circulation changes, and it is characterized by an equatorial EP warming resembling an El Niño-like pattern but confined to a much narrower band within 10°S-10°N (Fig. 1e).This EP warming, combined with the WP warming in the MD (Fig. 1d), results in the two warming centers in the FC (Fig. 1c) similar to observations (Fig. 1a).

Heat budget decomposition
To determine how wind-driven ocean circulation changes induce enhanced warming in the equatorial EP, we conduct an ocean mixed layer heat budget analysis 26,42,43 (see Methods) and decompose the SST trend (T t s ) into 7 terms accounting for atmospheric adjustment and ocean dynamical processes (Fig. 2; Supplementary Fig. 2).The decomposition is summarized in Eq. ( 8) and the detailed methodology can be found in the Methods section.
First, the cumulative impact of these 7 factors (T t s ) can accurately reproduce the simulated SST trend in the FC-MD ensemble mean (Fig. 2b; c.f. Fig. 1e) with a spatial correlation of 0.94, thereby validating our methodology.Among the 7 terms, the dominant factor for the reduction in zonal SST gradient is the ocean heat transport change (T t Ocn ; see the quantitative decomposition in Fig. 2a), which is characterized by a strong warming confined to the eastern equatorial Pacific (Fig. 2c).It is partially offset by the shortwave radiative flux change (T t SW ), while other components related to atmospheric adjustment are generally weaker (Supplementary Fig. 2), including the changes in longwave radiative flux (T t LW ), sensible heat flux (T t SH ), and latent heat flux due to the adjustment of near-surface wind speed (T t LH,W ), near-surface relative humidity (T t LH,RH ) and air-sea temperature difference (T t LH,ΔT ).In summary, this decomposition confirms that dynamical oceanic processes dominate the wind-driven equatorial eastern Pacific warming.
Next, we further decompose the ocean mixed layer heat budget and compute 9 advective terms to elucidate specific ocean dynamical processes (see Eq. ( 16) in Methods).The cumulative effect of the 9 advective terms (2 linear terms and 1 nonlinear term in each direction; Supplementary Fig. 3) generally resembles the T t Ocn pattern (r = 0.80), indicating that the residual term (e.g.mixing, diffusion, etc.) is relatively small (Fig. 3a; Supplementary Fig. 3).Additionally, the three nonlinear terms are generally small.Among the 6 linear terms, -v t ∂ T ∂y , -u ∂T t ∂x , and -v ∂T t ∂y are the major contributing factors to the equatorial EP warming as are discussed further below (Fig. 3a).For -v t ∂ T ∂y , the offequatorial westerly wind anomalies in the EP around 5°N induce an anomalous southward Ekman flow, which transports the warm water beneath towards the equator, thereby warming the EP cold tongue (Fig. 3b, c).This meridional warm advection is particularly strong in this region due to the strong meridional SST gradient in the EP (Fig. 3c).The climatological divergence of surface ocean currents further spread the equatorial EP warm anomaly (i.e.-u ∂T t ∂x , and -v ∂T t ∂y ) both poleward and westward (Fig. 3d-g).
Surprisingly, vertical advection terms, and the thermocline adjustment term (-w ∂T t ∂z ), play only a minor role in the SST changes over the EP (Fig. 3a; Supplementary Fig. 3).The anomalous westerlies in the equatorial EP induce an anomalous downwelling that acts to warm the Fig. 2 | The decomposition of sea surface temperature (SST) trend differences between fully coupled ensemble mean (FC) and mechanically decoupled ensemble mean (MD) based on ocean mixed layer heat budget.a Bar chart showing the contribution of each component of the ocean mixed layer heat budget to FC-MD SST trends; see Methods for details.The region averaged terms consist of the western Pacific (5°S-5°N, 120°E-180°, blue box in Fig. 1a), the central Pacific (5°S-5°N, 180°−130°W, orange box in Fig. 1a), and the eastern Pacific (5°S-5°N, 130°W-80°W, red box in Fig. 1a).b Panel showing the SST spatial pattern of the sum of all components (T t s ) and spatial correlation coefficient (r) with the FC-MD SST trend (Fig. 1e).c Panel showing spatial pattern of the contribution of ocean dynamics processes to SST trends (T t Ocn ).Stippling indicates significance at the 95% confidence level.
local SST (-w t ∂ T ∂z ), but that effect is rather weak (Supplementary Fig. 3j), as the westerlies are found mostly off the equator (Fig. 1e).The climatological upwelling acts to damp the EP warming in FC-MD because of the enhanced vertical stratification (-w ∂T t ∂z ), but that effect is also small as it peaks around the thermocline depth (Supplementary Figs.3k and 4h).By applying a similar heat budget analysis to the MD simulations (Supplementary Fig. 5), we find that the dominant mechanism for the Fig. 3 | The decomposition of ocean dynamical terms differences between fully coupled ensemble mean (FC) and mechanically decoupled ensemble mean (MD) based on the oceanic mixed-layer heat budget analysis.a The diagnosed results of advective terms in the oceanic mixed-layer heat budget as described in Eqs. ( 16)- (25), with the same selected regions as in Fig. 2a; see details in Methods.enhanced SST gradient along the equator is the shortwave radiative flux changing (Fig. 4a; Fig. 4e).An increase (decrease) in low and middle clouds leads to a cooling (warming) in the eastern (western) equatorial Pacific (Supplementary Fig. 6).In contrast, longwave radiative flux changes (i.e.greenhouse effect) induce an El Niño-like warming due partly to the "evaporative cooling" mechanism (Fig. 4d) 13 , but this effect is relatively small compared to the shortwave effect.The hemispheric asymmetry of SST warming is dominated by the effect of near-surface wind speed changes (Fig. 4b, c).The cross-equatorial winds increase (decrease) the near-surface wind speed at the south (north) of the equator (Fig. 4c), thus leading to the northern hemisphere being relatively warmer.

Discussion
Our study demonstrates that the complex warming structure in the tropical Pacific as observed since 1958 is a result of multiple physical processes.Historical radiative forcing can induce a subtropical northeastern Pacific warming extending towards the equatorial WP without invoking wind-driven ocean circulation changes.In contrast, the equatorial EP warming is driven by cross-equatorial wind anomalies in turn induced by the subtropical northeastern Pacific warming.
Interestingly, off-equatorial Ekman transport in the meridional direction turns out to be the dominant triggering mechanism for such warming.The warming effect of anomalous downwelling in the EP is relatively weak as the anomalous westerlies are located mostly off the equator, which emphasizes the recharge effect 44 due to the changes in wind stress curl instead of the Bjerknes feedback.This result highlights an overlooked role of wind-driven ocean circulation changes in the tropical Pacific warming pattern over the historical period.
Many previous studies have been devoted to understanding the cause of cross-equatorial winds in the EP, mostly related to the local or remote hemispheric asymmetry in warming as we see in our simulations 33 .The proposed mechanisms include asymmetric aerosol forcing 45,46 , tropical North Atlantic warming 47 , Atlantic meridional overturning circulation (AMOC) 48,49 , wind-evaporation-SST feedback (WES) 26,50,51 , and the Southern Ocean heat uptake [52][53][54] .A detailed attribution for the historical period and whether such a hemispheric asymmetry could continue into the future need to be addressed in follow-up studies.
Although the externally forced response in CESM2 can reasonably capture the equatorial WP and EP warming, it fails to reproduce the equatorial CP cooling seen in observations, like most CMIP6 models (Fig. 1a-c).Internal climate variability could potentially explain a part of the observation-model mismatch, but the mismatch may also originate from model biases in the tropical Pacific [55][56][57][58][59][60] .For example, our heat budget analysis reveals a CP cooling effect by the climatological upwelling in FC (supplementary Fig. 4b), consistent with the thermostat mechanism 27 , but this effect is rather weak in CESM2.A model bias in the magnitude of equatorial Pacific upwelling, in turn driven by surface winds, can potentially lead to a bias in the upwelling-induced CP cooling rate, and needs to be carefully assessed in future investigations.Other model biases, including the double ITCZ and the cold tongue biases, and their potential impacts on the tropical Pacific warming pattern also need to be addressed.

Observational data sets
We use five monthly observational SST datasets: (1) National Oceanic and Atmospheric Administration Extended Reconstructed Sea Surface Temperature Version 5 (ERSSTv5) with resolution 2°× 2°6  5) Kaplan Extended SST v2 with resolution 5°× 5°6 5 .In this study, the period from 1958 to 2014 has been chosen for calculating the SST trends.We chose 1958 as the starting year by following refs.1,25, which was determined according to the availability of ECMWF/ORAS4 ocean reanalysis data, and we chose 2014 as the ending year because our mechanically decoupled experiments followed the standard of the CMIP6 historical simulations that ended in the year 2014.

Large ensemble historical simulations
In this study, we analyze two large ensemble coupled climate model simulations using the Community Earth System Model version 2 (CESM2) 66,67 .Two CESM2 ensembles both have nominally 1°horizontal resolution and are forced with realistic, time-varying 1850-2014 external forcing, including greenhouse gasses, natural and anthropogenic aerosol emissions, and solar radiation.Smoothed biomass burning forcing was used for these simulations 68 .
The first CESM2 ensemble, referred to as FC for "Fully Coupled Model", is created by branching 50 ensemble members from a preindustrial model run for over 1000 years by NCAR.Ten of the ensemble members are branched from different pre-industrial control simulation dates (macro ensemble).The other 40 ensemble members are created by branching from four pre-industrial control simulation starting dates and adding a random perturbation to the atmospheric potential temperature field to create 10 micro ensemble members per branch date 68 .These runs are publicly available as part of the CESM2 Large ensemble.In the FC ensemble, the ocean and atmosphere exchange buoyancy fluxes and wind stress.
Wind stress overriding experiments have been demonstrated to be helpful in isolating the influence of wind-driven ocean circulation on the ocean-atmosphere coupled system [36][37][38][39][40] .The second CESM2 ensemble, MD for "Mechanically Decoupled Model", is created by branching 20 members from a pre-industrial MD model run 33 .Each ensemble member is branched from a different pre-industrial control date, making it a macro ensemble.The 6-hourly climatological wind stress forcing in the MD ensemble is calculated from 50 years of hourly wind stress output from an FC pre-industrial run.The ocean and atmosphere exchange time-varying buoyancy fluxes, but the ocean is forced by a seasonally-varying fixed wind stress climatology, calculated from pre-industrial conditions.More details can be found in ref. 33.

Ocean mixed layer heat budget analysis
To further understand the mechanisms of wind-driven ocean circulation changes that induce an eastern equatorial Pacific warming, we analyzed the ocean mixed layer heat budget following previous studies 42,43 .The ocean surface mixed-layer budget equation is: Here, the left-hand side represents the mixed-layer heat storage term, where ρ s is the ocean density, c p is the specific ocean heat, H is the time-varying mixed layer depth from the model, and SST is the mixedlayer temperature.The right-hand side consists of net surface shortwave (F SW ), longwave (F LW ) fluxes, sensible heat fluxes (SH) and latent heat fluxes (LH), and heat flux due to ocean dynamics (Ocn).
Here we define downward as positive.
When we take the linear trend of Eq. ( 1), the left-hand side which represents the trend of heat storage is negligible 26,69 , which is written as follows: the superscript t denotes linear trends from 1958 to 2014.The latent heat flux term is directly related to SST via saturation vapor pressure: where q a is the specific humidity of the air above the sea surface, L v is the latent heat of vaporization, ρ air is the density of the air and W is near-surface wind speed.The q a term can be expressed as: where RH 0 is the relative humidity at the sea surface, ΔT is the temperature gradient near the sea surface, which is defined as (SST -Ta).
Using the Clausius-Clapeyron equation, Eq. ( 4) can be expressed as: where α = L v R v T 2 ≈ 0:06 K À1 .We can substitute Eq. ( 5) into Eq.( 3) and obtain a new expression for the latent heat flux: Then, we use linear regression to get the linear trend of latent heat flux from Eq. ( 6): We represent SST t as T t s in the equation, the latent heat term can be further broken down into contributions from trends in near-surface wind speed (T t LH,W ), near-surface relative humidity (T t LH,RH ), and air-sea temperature gradient (T t LH,ΔT ), where the right-hand side terms can be calculated as follows Eqs. ( 9)- (15).
Here, the cumulative of these 7 factors (T t s ) may not exactly match with the model-simulated SST due to the assumptions made during the derivation process.

Oceanic mixed-layer heat budget Equations
To understand the roles of ocean dynamical processes in causing the eastern equatorial Pacific warming, we further decompose the influence of ocean heat transport changes as follows: where the overbar denotes the climatological mean state, the prime denotes anomaly, and the superscript 't' denotes long-term trend.Each zonal, meridional, and vertical ocean 3-dimensional advection term consists of three components, two linear terms, and one nonlinear term.T t advection is the sum of the terms on the right-hand side of Eq. ( 16), excluding the residual terms (R), which means that the residual terms are not as important for the heat budget at the mixed layer depth.We use T t advection for comparison with T t Ocn .T is the mixedlayer average temperature, T b is the bottom of mixed-layer temperature, u, v, and w represent 3D mixed-layer current velocity, w b is the bottom of mixed-layer temperature and current velocity, and H denotes the climatological seasonally varying mixed-layer depth in the equatorial Pacific, which is located at a 0.5 °C difference from the surface SST.R is the residual term.The term ρ s c p H αLH varies in space much less than those advective terms.The climatological annual cycle is calculated based on the period from 1958-2014.The anomaly field is calculated by subtracting the monthly mean field from its climatological annual cycle at each period respectively.
For clarity, the terms on the right-hand side of Eq. ( 16) are abbreviated in the figures as follows, Trends and significance All the linear trends are calculated by applying a linear least-squares regression model to the spatially integrated time series.To determine the regions and time periods where trends are significantly different, we consider the spread of the ensemble members as a normal distribution.To test if the distributions are significantly different, we calculate the Z statistic as follows Eq. ( 26) and use 95% significance (Z ≥ 1.96) as a threshold, where: X is the ensemble mean from each experiment.with shading showing the 95% confidence interval.σ is the standard deviation of each ensemble member divided by the square root of the number of ensemble members.cobe2.html;Kaplan Extended SST v2 at https://psl.noaa.gov/data/gridded/data.kaplan_sst.html;CMIP6 data at https://esgf-node.llnl.gov/search/cmip6.CESM2 Fully-Coupled (FC) Large ensemble members at https://www.cesm.ucar.edu/projects/community-projects/LENS2/.Select CESM2 Mechanically Decoupled (MD) pre-industrial and used Large ensemble data is found at https://www.earthsystemgrid.org/dataset/ucar.cgd.cesm2.mdpc.htmland https:// zenodo.org/records/10484207.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.

ρ
s c p H αLH ∂ðT 0 Þ t ∂t is the mixed layer temperature tendency, here abbreviated as is ∂T t ∂t .b The spatial pattern of -v t ∂ T ∂y term (unit: °C/decade).c Decomposing -v t ∂ T ∂y into the trend of meridional advective currents term, denoted as v t (colors; cm/s/decade) and the sea surface temperature (SST) of the climatology term, represented as T (black line; unit: °C).d The spatial pattern of -v ∂T t ∂y term (unit: °C/decade).e Decomposition of -v ∂T t ∂y into the meridional advective currents of the climatology term, denoted as v (colors; unit: cm/s) and the trend of SST term, represented as T t (the black solid line represents the positive SST trend, the black dashed line represents the negative SST trend; unit: °C/decade).f The spatial pattern of -u ∂T t ∂y term (unit: °C/decade).g Decomposition of -u ∂T t ∂y into the zonal advective currents of the climatology term, denoted as u (colors; unit: cm/s) and the trend of SST term, represented as T t (the black solid line represents the positive SST trend, the black dashed line represents the negative SST trend; unit: C/decade).Stippling indicates significance at the 95% confidence level.

Fig. 4 |
Fig. 4 | Mechanically decoupled ensemble mean (MD) sea surface temperature (SST) trend decomposition based on the ocean mixed layer heat budget analysis.a Same bar chart as Fig. 2a, but for MD.b Panel showing the SST spatial pattern of the sum of all components (T t s ) and spatial correlation coefficient (r) with the MD SST trend (Fig. 1d).c Panel showing the spatial pattern of the contributions from trends in near-surface wind speed (T t LH,W ) (colors; unit: °C/decade), with surface wind trend overlaid (arrows; unit: m/s/decade).d, e Panels showing the spatial pattern of the contribution of longwave radiative flux (T t LW ) and surface shortwave radiative flux (T t SW ) to SST trends (colors; unit:°C/decade), respectively.Stippling indicates significance at the 95% confidence level.