Upper-ocean stratification changes control ENSO amplitude shift under sustained global warming

Domain: chemistry
Published: 2026-02-24
Editorial outcome: major_revision

Abstract
As projected by climate models in the high emission scenarios, the El Niño–Southern Oscillation (ENSO) exhibits a non-monotonic amplitude shift. However, its key drivers remain poorly quantified. Here we introduce a framework using an intermediate coupled model (ICM) that coherently represents mean-state oceanic climatologies in the tropical Pacific, derived from eight selected climate models across three periods (1940–1990, 2040–2090 and 2240–2290). By applying vertical baroclinic mode decomposition to ocean density, we extract wind projection coefficients (p ; n is mode number) governing upper-ocean dynamical responses. The ICM with the explicitly prescribed climatological fields, including stratification and the thermocline structure, successfully reproduces the non-monotonic ENSO shifts, which is illustrated to be primarily driven by opposite changes in p1 and p2 post 2140. Sensitivity experiments further confirm stratification as the dominant modulator. This study establishes a coherent mechanistic framework for disentangling stratification impacts on ENSO in climate model projections under global warming.

Introduction
As the global climate continues to warm1–7, the El Niño–Southern Oscillation (ENSO)—the dominant mode of interannual variability—undergoes significant changes that have far-reaching impacts on global weather, ecosystems, and socio-economic systems8–14. Understanding how ENSO responds to sustained warming is a critical scientific challenge, especially given its influence on rainfall extremes, droughts, and marine productivity worldwide15–21. Recent advances in climate modeling have leveraged coupled general circulation models (CGCMs) to project ENSO behavior under future climate scenarios, including those represented in the Coupled Model Intercomparison Project Phase 6 (CMIP6)22–24.

Despite model-to-model differences, some of the CMIP6 models under the Shared Socioeconomic Pathway 585 (SSP585) scenario reveal a coherent yet non-monotonic response of ENSO amplitude to warming25–27. Specifically, ENSO variability is weak during 1940–1990, intensifies during 2040–2090, and then weakens again after 2100. However, the physical processes and mechanisms underlying these shifts remain poorly understood.

ENSO variability is tightly linked to the climatological mean state of the tropical Pacific Ocean28–30. Prior studies have pointed to key components of this mean state in the ocean—such as stratification, thermocline depth, surface-layer fields (currents and sea surface temperature (SST))—as important modulators of ENSO15,31–38. While CGCM simulations provide valuable datasets, the interactions among these components and their net effects on ENSO are complex and intertwined, making it difficult to rigorously isolate the individual contributions of each factor. As a result, many previous assessments have been largely qualitative and diagnostic in nature. Although often mentioned as a key influence in the literatures, stratification effects have rarely been represented and quantified in a consistent or physically explicit way26,27,39–43. In particular, the role played by upper-ocean stratification in modulating ENSO under global warming remains elusive.

This study addresses that gap by developing a modeling framework using an intermediate coupled model (ICM) in which three components of the mean ocean state—(1) upper-ocean stratification, (2) surface-layer fields (velocity and SST), and (3) thermocline structure—are explicitly prescribed based on eight CMIP6 model mean (8CMM) simulations across historical and future climate change scenarios during 1990–2300. This modeling framework enables us to systematically isolate and quantify individual and collective impacts of these mean oceanic components on ENSO responses under sustained global warming. In particular, we focus on the role of oceanic stratification, represented through vertical baroclinic mode decomposition, to explore how changes in ocean density structure alter the dynamical and thermodynamical response characteristics that govern ENSO variability. The mechanistic explanation resolves a key paradox in CMIP6 simulations: ENSO amplitude declines after 2100 despite continued increases in upper-ocean mean stratification.

Results
(1) Changes in ENSO amplitudes and the mean ocean state
Previous studies have comprehensively analyzed the ENSO responses and changes in the mean atmospheric and oceanic states using CMIP6 simulations under different climate change scenarios27,42. Here, we focus on CMIP6 historical and future simulations under SSP585 during 1990-2300 to illustrate ENSO responses (see the details in Method). To adequately represent ENSO variability, we use the standard deviation or magnitude of Niño3.4 or equatorial Pacific SST anomalies as proxies for ENSO amplitude, rather than the amplitudes of explicitly identified ENSO events. Also, we consider interannual variability of Niño3.4 or equatorial Pacific SST, which are used to reflect ENSO amplitude. So, SST anomalies used in the following analyses are based solely on the interannual component, which is obtained by subtracting the 51-year running mean climatology for each selected CMIP6 model.

Figure 1a shows the amplitude of Niño3.4 SST anomalies from eight individual CMIP6 models during 1900–2300 (Supplementary Table 1) and their ensemble mean; a 51-year running mean is performed to highlight a consistent change in ENSO amplitudes among the selected 8CMM simulations. The 8CMM in historical and future simulations reveal a non-monotonic evolution of ENSO amplitude: ENSO variability is weak during the historical simulation, increases during the future simulation up to around 2100, and is followed by a decline post-2100.

Specifically, SST fields from the 8CMM with a 51-year running mean show three distinct periods in terms of ENSO amplitudes: weak ENSO amplitude during 1940–1990 (denoted as P1), an enhancement during 2040–2090 (P2), and a decline during 2240–2290 (P3), respectively. The differences in ENSO amplitudes are quantified by the standard deviation of SST anomalies calculated separately from the 8CMM during these three different periods, which is shown in Fig. 2a.

These three periods also correspond to pronounced changes in the tropical Pacific mean ocean state, including upper-ocean stratification, surface wind and SST patterns, and thermocline structure; previous studies have analyzed these atmospheric and oceanic changes in great details25–27. Figure 1b, c illustrate examples for the pronounced changes in the stratification (as seen in upper-ocean density and the corresponding buoyancy frequency, N2) over the central equatorial Pacific during 1900-2300; these fields are obtained by performing a 51-year running mean based on monthly outputs. The stratification exhibits clear temporal changes: a stable stratification up to 2000, a gradual increase during 2000–2160, and strong stratification after 2160.

Also, we calculate the corresponding ocean climatology from the 8CMM simulations and highlight their differences for the three periods (Supplementary Figs. 1–3). There are apparent changes in SST and surface winds, characterized by reduced surface currents and zonal SST gradients. As a result, the equatorial upwelling is reduced in association with the changing mean surface velocity and SST fields. In addition, the thermocline shows pronounced changes in its depth, intensity and sharpness (Supplementary Fig. 4). For example, the mean thermocline continues to shoal in the equatorial Pacific, with its increased sharpness. In particular, the stratification exhibits non-uniform vertical changes in the upper ocean, with an increase in the upper 120 m, but a decrease below (Supplementary Fig. 5).

(2) ENSO amplitude changes in ICM simulations
Under sustained warming in the 8CMM, the mean ocean state in the tropical Pacific undergoes pronounced changes. From a process-based perspective, the mean ocean state modulating ENSO can be decomposed into three distinct components: stratification, surface fields (velocity and SST), and thermocline structure (Methods; Supplementary Table 2).

To quantify the influences of these mean ocean state on ENSO, we employ an intermediate coupled model (ICM)44,45 (see Method, Supplementary Text 1-3 and Supplementary Fig. 6). The ICM is an anomaly model featuring an empirically constructed atmospheric component to represent surface wind stress response to SST variability45, which is coupled to a simplified dynamic ocean model46,47. A key advantage of this framework is its capacity to directly produce interannual anomalies, while oceanic climatologies are allowed to be prescribed. Although the ICM has been widely applied in ENSO simulations and predictions44,45,48–52, the present version introduces several key advances (Supplementary Text 4). Notably, the thermocline component is now presented through an entrainment temperature (Te) parameterization, which is analytically formulated using a hyperbolic tangent (tanh) function53. This formulation explicitly links Te to mean thermocline properties (i.e., its depth, intensity, and sharpness), which collectively modulate the sensitivity of Te to thermocline depth anomalies. Furthermore, spatially-varying density and stratification fields are incorporated to derive corresponding modal parameters (e.g., wind projection coefficients, pn), enabling a more realistic representation of oceanic feedbacks. Critically, all three influencing components are coherently and explicitly integrated within the unified ICM framework. Then, this enhanced modeling framework is applied to the post-2100 period in CMIP6 simulations, where a counterintuitive decline in ENSO variability emerges despite continued warming, a phenomenon that existing theories struggle to explain. Through this approach, we systematically isolate the relative contributions of mean ocean state components to ENSO amplitude changes, offering a mechanistic interpretation of complex climate projections. Nevertheless, the current ICM framework and its modeling study for ENSO modulations in the global warming context also exhibits limitations in its representation of ENSO variability (Supplementary Text 4).

Considering the ENSO amplitude changes in the 8CMM simulations, three representative periods (1940–1990, 2040–2090, and 2240–2290) are separately considered to characterize changes in three components of the mean ocean state and their possible relationships with ENSO. We derive mean ocean fields, including stratification, surface velocity, and SST, and thermocline properties from the 51-year running mean of the 8CMM during the three periods. These related fields and parameters (including modal properties) are then explicitly incorporated into the ICM to assess ENSO responses under different mean ocean climatologies, which are fully three-dimensional and are all spatially varying (Method). Each experiment is freely run for a total of 10000 inter-realizations based on a Bootstrap method. The eight CMIP6 model mean (8CMM) upper-ocean (0–300 m) density in the central equatorial Pacific (CEP; 160°E–160°W, 2°S–2°N; contours with an interval of 1 kg m-3). The corresponding 8CMM upper-ocean (0-300 m) stratification (represented as buoyancy frequency, N2) in the CEP (contours with an interval of 5 × 10-4 s-2). The colored areas and vertical lines are added to highlight the three periods for 1940–1990 (P1), 2040-2090 (P2), and 2240–2290 (P3), respectively.

Figure 2c presents standard deviation of SST anomalies from the ICM simulations when the 8CMM ocean climatological fields estimated during the three periods are prescribed; these results can be compared with those directly from the 8CMM (Fig. 2a). Consistent with the 8CMM showing pronounced changes in ENSO amplitude during 1900–2300, ENSO amplitudes in the ICM simulations also exhibit apparent shifts (Fig. 2b): low amplitude during 1940-1990, enhanced one during 2040-2090, and low amplitude after 2100, respectively. Figure 3 illustrates examples for longitude-time sections of simulated SST anomalies over the equatorial Pacific from the ICM. Notably, weak ENSO variability during P1, enhanced ENSO variability during P2 and its later suppression during P3 are all simulated very well.

(3) Mechanism Analyses
The comparable SST amplitude changes and associated feedback structures between the 8CMM and the ICM-based experiments (Fig. 2) affirm the ICM’s capability to represent essential ENSO modulation-related processes under long-term global warming. Thus, the ICM experiments serve as a robust basis for mechanistic inquiry (Method; Supplementary Table 3). As indicated in Fig. 1a, variations in ENSO amplitude across the three distinct periods—P1, P2, and P3—coincide with notable changes in the mean ocean state (Fig. 1b, c). These ENSO modulations involve multiple oceanic parameters and fields, and three primary components are explicitly incorporated in the ICM framework. We begin with a focused analysis on the stratification component and its dynamical implications for ENSO variability.

(a) Baroclinic mode responses controlled by stratification and their relationship to ENSO variability. Whereas previous studies have offered qualitative assessments of stratification impacts, we here provide a quantitative treatment using vertical baroclinic mode decomposition32,37–39,46,47 (Methods). This approach is applied to upper-ocean density fields from the 8CMM to extract spatially varying vertical structure functions (ψn), wind projection coefficients (pn), and wave speeds (cn) for the first ten baroclinic modes (n = 1–10).

Figure 4 illustrates the spatial distributions of p1, p2, and p3 across the equatorial Pacific for each period. These parameters exhibit well-defined spatial structure across the equatorial Pacific. The amplitude of the first baroclinic mode (p1) decreases eastward along the equator, indicating a more surface-intensified response in the western-central equatorial Pacific. Conversely, p2 and p3 amplitudes increase toward the east, reflecting deeper thermocline sensitivities in the eastern equatorial Pacific. Using these modal parameters, reconstruction calculations for upper-ocean variability indicate that the first three baroclinic modes dominate the ocean’s dynamical responses to wind forcing.

As the climate becomes warm over the 1900–2300 period, substantial changes in the mean density field (Fig. 1b, c) yield corresponding variations in stratification and baroclinic mode parameters (Fig. 4 and Supplementary Figs. 4 and 5). One major finding from this study is the relatively changing amplitude of the first three vertical baroclinic modes across the three periods, as derived from the modal decomposition of stratification. While wave speeds (cn) remain relatively stable, the wind projection coefficients (pn) undergo pronounced temporal changes across P1, P2, and P3, which determine the vertical structure and spatial patterning of the ocean’s responses to wind forcing.

From P1 to P2, p1 increases notably in the western-central equatorial Pacific, signaling strengthened surface-layer responses. Simultaneously, p2 and p3 exhibit marked intensification in the eastern equatorial Pacific, indicating enhanced thermocline-layer effects. From P2 to P3, however, p1 decreases systematically, particularly over the western-central equatorial Pacific, whereas p2 exhibits a zonal dipole structure—increasing west of ~120°W and decreasing to the east. A similar pattern is observed in p3, with a weaker dipole centered near ~110°W.

These findings highlight a key transition in the wind projection coefficients that are determined by stratification: during P3, the amplitude of p1 decreases systematically, while p2 increases in the western-central equatorial Pacific. This shift reflects a redistribution of wind energy from surface layers to deeper thermocline layers, ultimately reducing the net upper-ocean response amplitudes. Consequently, given the same wind anomalies, the ocean’s dynamical response becomes attenuated in the vertical direction during P3, contributing to reduced ENSO amplitude.

The redistribution in the modal wind projections with opposite changes in p1 and p2 amplitudes from P2 to P3 corresponds to a broader regime shift in ENSO dynamics. During P2, increased stratification enhances the projection of wind forcing onto all three leading modes, reinforcing surface and thermocline anomalies and intensifying ENSO variability. During P3, in contrast, the diminished p1 amplitude suppresses surface-layer responses, while the enhancement of p2 in the western-central Pacific bolsters subsurface responses. These competing effects on upper-ocean responses during P3, with the decreased p1 amplitude but increased p2 amplitude, result in partial cancellation of their contributions to anomaly responses, causing reduced ENSO variability despite continued increases in stratification.

These are only qualitative arguments based on analyses from the 8CMM simulations. In addition, the stratification effects and changes are accompanied with other effects associated with changes in surface fields and the thermocline (Supplementary Figs. 1–5). So, there exist interactions among these influencing components, with their net effects on ENSO being complex and intertwined, making it difficult to rigorously isolate the individual contributions of each factor.

(b) Sensitivity experiments and component-specific contributions. Given the interdependence of mean stratification, thermocline structure, and surface fields, it is non-trivial to isolate their individual contributions of influencing components to ENSO amplitude changes using CMIP6 outputs alone. Therefore, a suite of controlled sensitivity experiments is conducted using the ICM, in which each component of the mean ocean state (stratification, surface velocity/SST, and thermocline structure) is independently varied across the three periods. Each experiment involves a 50-year free run using the ICM, during which a single mean state component from one period is substituted into the baseline conditions of another period. This allows us to isolate and quantify the individual and combined effects of each component on ENSO amplitude.

As shown in Fig. 5a for considering the two periods (P1 and P2), replacing the P1 stratification with that from P2 (denoted as ICM-P1 P2-N2) leads to a significant enhancement of ENSO amplitude relative to the P1 control run. In contrast, substituting the surface velocity/SST (ICM-P1 P2-U) and thermocline components (ICM-P1 P2-Te) respectively from P2 into P1 both act to dampen ENSO amplitude. Among these, the stratification effect is the most influential, confirming its primary amplifying role in driving the P1-to-P2 increase in ENSO variability. Here, the enhanced stratification increases the wind projection onto the first three baroclinic modes (Fig. 4), thereby amplifying both surface-layer and thermocline effects through stronger surface-layer and thermocline anomalies.

The P2-to-P3 transition presents a more complex and puzzling picture (Fig. 5b). Notably, when the P3 stratification (all pn fields) is substituted into the P2 baseline (ICM-P2 P3-N2), ENSO amplitude does not decrease as expected. To resolve this apparent paradox, we then perform additional mode-specific experiments: replacing only p1 or p2 from P3 into the P2 stratification configuration. The results (Fig. 5c) reveal that the decline in p1 during P3 leads to reduced ENSO amplitude, while the rise in p2 leads to an increase. Thus, the effects due to the decreased p1 but increased p2 are opposing in terms of their contributions to ENSO amplitude; their nearly canceling effects on upper-ocean responses explain why ENSO amplitude remains largely unchanged between the ICM-P2 baseline and ICM-P2 P3-N2 as shown in Fig. 5b.

Thus, the observed ENSO amplitude reduction from P2 to P3 is attributable to the spatial redistribution of baroclinic mode contributions: the decreased p1 weakens surface responses primarily in the western-central equatorial Pacific, while the increased p2 enhances subsurface effects. Meanwhile, in the eastern equatorial Pacific, both p2 and p3 decrease during P3, acting to weaken the subsurface effects. Then, despite continued stratification intensification during P3, the wind-induced upper-ocean responses become increasingly inefficient, leading to a weakened ENSO amplitude. Notably, the intensification of stratification during P3 does not translate into stronger ENSO variability. Instead, wind projections onto the surface-layer (p1) and eastern equatorial Pacific thermocline layers (p2) are diminished, while higher-mode contributions in the western-central equatorial Pacific increase. These opposing effects result in weakened surface-layer currents and thermocline fluctuations, particularly in the eastern equatorial Pacific, thereby suppressing ENSO amplitude.

In addition to stratification, changes in the other two components—the mean surface velocity–SST and the thermocline—also play important roles in damping ENSO amplitude (Fig. 5a, b). For example, replacing the P1 surface velocity–SST component with that from P2 (denoted as ICM-P1 P2-U) leads to a reduction in ENSO amplitude relative to the P1 control run (Fig. 5a). Similarly, substituting the P1 thermocline component with that from P2 (ICM-P1 P2-Te) leads to a reduction in ENSO amplitude (Fig. 5a). In a parallel manner, replacing respectively the surface velocity-SST component and thermocline component from P2 with those from P3 (ICM-P2 P3-U and ICM-P2 P3-Te) leads to a reduction in ENSO amplitude relative to the P2 control run (Fig. 5b).

It is worth noting that the influence of the mean thermocline changes on Te is not apparent by directly examining the spatial distributions of the thermocline-related mean fields alone (Supplementary Figs. 2–4). However, its effects can be illustrated and quantified from ICM-based simulations, as clearly demonstrated in Fig. 5a, b. In this study, the tanh-based parameterization scheme for Te, coherently incorporated into the ICM, allow to estimate upper-ocean temperature anomalies (the equation (3.4) in Supplementary Text 3) and changes in the relationship between the mean thermocline properties and Te. The ICM-based sensitivity simulations reveal period-dependent changes in ENSO amplitude, resulting from coherent interplays among upper-ocean stratification, mean thermocline properties and surface fields, which is consistent with the changes of ENSO amplitude seen in the 8CMM. These results provide another evidence for using the ICM as an appropriate modeling tool to study the physical processes that drive future changes in ENSO amplitude.

(c) Further mechanistic insights into post-2100 ENSO weakening.
The proposed mechanism invokes changes in the wind projection coefficients (pn) derived from stratification and their amplifying effects on ENSO, which are accompanied by damping effects from the thermocline and surface field components. The increases in the amplitudes of p1 and p2 offer a plausible explanation for the amplification of ENSO until the late-21st century (during P2). However, its ability to account for the subsequent weakening beyond 2100 (e.g., Fig. 5) appears less compelling. A particular challenge lies in explaining why continued stratification strengthening from period P2 to P3 differentially projects onto ocean vertical baroclinic modes (primarily p1 and p2). This sub-section aims to clarify this issue through additional analyses and discussions, thereby strengthening the physical interpretation for ENSO changes under prolonged warming.

As shown in Fig. 1c, the stratification in the upper ocean undergoes marked changes with progressive warming. Closer examination reveals that these changes are non-uniform in both vertical and horizontal directions. For instance, the vertical structure of stratification (Supplementary Fig. 5) indicates intensification in the upper 120 m, but a reduction at subsurface depths between 120–150 m. This vertical reconstruction of the stratification alters the relative contributions of wind projections to the baroclinic mode amplitudes: the amplitude of p1 decreases while that of p2 increases during P3 (Supplementary Fig. 4), reflecting an amplifying of surface-layer stratification but vertically a gradual weakening of subsurface stratification. From P2 to P3, this shift in the vertical distribution of stratification modifies the respective influences of p1 and p2 on upper-ocean responses. During P3, the reduced amplifying effects on ocean anomalies from the decreased p1 amplitude are insufficient to offset the increased amplifying effects from the increased p2 amplitude, resulting in reduced stratification-related net amplifying effects.

It should be noted that other processes associated with mean ocean fields also contribute to ENSO modulation in the coupled system. As indicated in Fig. 5, the Te modulation due to changes in thermocline depth, strength, and sharpness exerts a damping influence on ENSO variability. Similarly, changes in the mean surface ocean fields also play a suppressive role (Fig. 5). Thus, in the ICM framework, three key components coexist that represent their competing effects on ENSO: the thermocline and surface field components both dampen ENSO, whereas the stratification generally amplifies it, which is determined by the respective amplitudes of p1 and p2 (i.e., when both p1 and p2 increase, ENSO intensifies). Thus, the overall ENSO amplitude is governed by the combined and synergistic effects of these influencing components, and the resultant ENSO amplitude reflects a balance between amplifying and damping effects in the coupled system. In particular, during P3, p1 decreases while p2 increases, resulting in a partial cancellation of their amplifying contribution to ENSO amplitudes. Concurrently, damping mechanisms remain active. When the thermocline- and surface fields-related damping effects dominate, the amplification due to changes in stratification is insufficient to offset the damping, and the net effect leads to a reduction in ENSO amplitude—as is the case during P3.

(4) Budget analyses for ENSO feedback processes
The underlying mechanisms are further analyzed using the BJ-index and its contributing terms54 (Method; Supplementary Table 4; Supplementary Text 6), which are calculated from the 8CMM and ICM simulations (Fig. 2b, d), respectively. Various feedback processes responsible for ENSO modulations show coherent relationships with SST variabilities, including the thermocline feedback (TH), dynamical damping (DD), and thermodynamical damping (TD), respectively. These budget analyses of the mixed-layer SST further corroborate the changes in the three key influencing components of the mean ocean state and dominant feedbacks. In general, when the first baroclinic mode (p1) dominates in the upper-ocean responses and its amplitude is large, surface-layer responses are intensified, leading to enhanced zonal advection feedback. In contrast, when higher-order baroclinic modes (p2 and p3) dominate and their amplitudes are large, thermocline anomalies strengthen, corresponding to enhanced thermocline effect.

These arguments indicate how changes in mean stratification influence the strength of major ENSO feedbacks in the tropical Pacific. During P2, for example, the increased stratification strengthens both zonal advection and thermocline feedbacks, leading to increased ENSO amplitude. The stratification favors pronounced thermocline effects (reflected in increased p2 and p3), allowing thermocline variability to exert a strong influence on SST through the thermocline feedback. Simultaneously, surface-layer responses (indicated by p1) are also amplified, enhancing zonal advection feedback through surface currents. The combined effect leads to strong ENSO variability during this period.

During P3, structural and amplitude shifts occur in the mean stratification and associated baroclinic modal properties. The reduced amplitude of p1 over the western-central equatorial Pacific reflects weakened surface-layer responses, resulting in diminished zonal advection feedback. Meanwhile, thermocline feedback weakens, particularly in the eastern equatorial Pacific, where the amplitudes of p2 and p3 decrease. These simultaneous reductions in thermocline and surface zonal advection feedbacks collectively dampen SST variability. Furthermore, weakened thermocline variability and associated vertical advection in the eastern equatorial Pacific under continued warming also contribute to the overall decline in thermocline feedback strength. Consequently, ENSO amplitude is reduced during P3.

The changes in ENSO amplitude across the three periods align with findings from the baroclinic mode analyses. For instance, the enhanced ENSO amplitude during P2 is associated with stronger mean stratification and corresponding increases in both zonal advection and thermocline feedbacks. Conversely, the weakened ENSO amplitude during P3 results from structural shifts in mean stratification that act to reduce the efficacy of these feedbacks, especially in key regions of the eastern equatorial Pacific.

This mean stratification-induced shifts in baroclinic mode properties—both in structure and amplitude—drive changes in the dominant ENSO feedback mechanisms. This mechanistic explanation resolves a key paradox in CMIP6 simulations: ENSO amplitude declines after 2100 despite continued increases in upper-ocean mean stratification. Our findings underscore the importance of not only the magnitude of stratification, but also its vertical structure in modulating ENSO variability. The ICM-based simulations (Fig. 2b, d) provides another evidence for using the ICM as an appropriate modeling tool for studying the physical feedbacks that drive future changes in ENSO amplitude.

Methods
CMIP6 simulations and analysis periods
To investigate ENSO responses to global warming, we analyze CMIP6 historical and future simulations under the high-emission SSP585 scenario, following methodologies established in previous studies27. These simulations are further used to derive climatological mean oceanic fields, which are subsequently prescribed in an ICM to quantitatively assess their influences on ENSO amplitude.

We select eight CMIP6 models (Supplementary Table 1) based on their consistent projection of a distinct non-monotonic ENSO amplitude changes under sustained warming, a feature that is not widely present across the full CMIP6 ensemble. The eight CMIP6 model mean (8CMM) is characterized by three phases in terms of ENSO amplitudes: a weak period (1940–1990, P1), an enhanced period (2040–2090, P2), and a later declining period (2240–2290, P3), respectively. It should be noted that the mechanisms identified using this subset of models require validation across a broader ensemble, given the considerable intermodel spread in ENSO projections. A detailed discussion of sampling considerations is further provided in Supplementary Text 5.

The selection of these three periods (1940-1990 (P1), 2040-2090 (P2), and 2240-2290 (P3)), while consistent with prior work27, reflects key differentially evolving behaviors in both ENSO amplitude and oceanic mean state. Notably, the changes in mean stratification associated modal parameters (pn) and ENSO amplitude are not synchronous—each peaks and declines with distinct timing (Fig. 1 and Supplementary data Fig. 5). The P3 period (2240-2290), for example, represents a later phase during which ENSO amplitude has declined considerably (Fig. 1). Should an early period (say 2100–2150) be considered, the ENSO amplitude is just seen to begin to decline under intensified stable stratification, with the first baroclinic mode (p1) amplitude beginning to decrease. In this case, the quantitative relationship is just seen to begin to decline under intensified stable stratification with the increased p1 amplitude. The results suggest that the physical processes identified are robust based on the earlier selected periods (1940–2150), but many changes during the stratification intensity and ENSO amplitude are asynchronous under global warming. Further analyses and modeling efforts are clearly needed to address these issues as the mean ocean state continues to evolve.

An intermediate coupled model (ICM)
We adopt an ICM to examine ENSO responses to changes in mean ocean state (Supplementary Text 1–3; Supplementary Fig. 6). The ICM is an anomaly model of intermediate complexity, which directly takes interannual anomaly fields as predictive variables from its governing equation, while mean ocean climatology is prescribed from observations or model simulations47,48. Its atmospheric component is represented by surface wind stress anomalies that are formulated by singular value decomposition (SVD) method, with anomaly coupling to the ocean. Symbolically, this coupling can be expressed as τ=ατ × f(SST), where ατ is a scalar parameter representing coupling intensity, and f denotes the spatial structure derived from SVD modes. This formulation allows ατ to modulate the strength of ocean–atmosphere coupling, offering a convenient means to tune model sensitivity. In all ICM experiments, the same value is prescribed for ατ (ατ=1.0). The oceanic component was developed by Keenlyside and Kleeman47 (2000); its dynamical ocean model is solved by the vertical mode decomposition method (Supplementary Text 1). An SST anomaly model is embedded into this dynamical ocean model to depict SST variability (Supplementary Text 2). In order to close the equation of SST anomalies, the temperature of subsurface waters entrained into the mixed layer (Te) need to be determined; previously, an SVD-based statistical model was developed to parameterize interannual variability of Te in terms of sea level anomalies for SST A simulations48,49. In this work, an analytical parameterization of Te has been developed to link Te to thermocline depth (or sea level) anomalies, as well as mean thermocline properties53 (Yuan et al. 2020; Supplementary Text 3). The ICM has been widely used for ENSO modeling studies48–52, including routine applications for real-time predictions since 2003 (see a summary of the model ENSO forecasts at the International Research Institute for Climate and Society (IRI) website: http://iri.columbia.edu/climate/ENSO/currentinfo/update.html; Supplementary Fig. 7).

To investigate the role of the mean ocean state in modulating ENSO under global warming, we develop a novel modeling framework based on the ICM (Supplementary Text 4), in which three mean ocean components in the tropical Pacific are explicitly incorporated into the ICM (Supplementary Text 2). This approach enables us to explicitly prescribe and systematically evaluate the effects of three mean components of the upper ocean state in the tropical Pacific. In this study, the ICM is freely run for 50 years to represent ENSO responses to different global warming conditions during the three periods. ENSO simulations from targeted ICM experiments can be compared with those from the 8CMM and can be further analyzed to quantify the effects of various components of the mean ocean state.

Vertical baroclinic mode decomposition performed using upper-ocean density
Given a density field, we perform vertical decomposition into baroclinic modes46,47. The relevant variables of the model in physical space (Equation 1.1 in Supplementary Text 1) can be transformed into modal space (Equation 1.3 in Supplementary Text 1) as a function of vertical structure; its n-th baroclinic mode is written to be ψn, which, given the background ocean density field (N2), is determined by

(N2 ψn)z = - c-2n ψn; the parameter cn denotes the characteristic speed of Kelvin waves associated with different baroclinic modes. This equation, together with boundary conditions, can be solved as a value problem to produce ψn; its integral in the vertical direction is given by

pn = H-1 ∫0H ψn dz, which govern the coupling strength between baroclinic modes and wind forcing; note that the lowercase letter is used for pn, whereas a capital letter is used in Supplementary Text 1 (that is, Pn in Equation 1.6). These coefficients determine how wind energy is distributed among different baroclinic modes, thereby shaping the characteristics of upper-ocean responses and influencing the sensitivity of ocean responses to atmospheric wind forcing.

In this study, we compute the corresponding wind projection coefficients (pn) and phase speeds (cn) for the first ten modes (n is the mode number). Different from previous similar studies32,33,35–39, these modal parameters in our modeling study are calculated in a spatially varying manner so that spatial structure of climatological mean density can be adequately taken into account. Computational analyses reveal pronounced spatial changes in pn across the equatorial Pacific, with p1 having a large amplitude in the western-central equatorial Pacific; in contrast, the second and third modes (p2 and p3) tend to have large amplitudes in the eastern equatorial Pacific. As these modal parameters determine the characteristics of upper ocean responses to wind, their spatial structures control how surface wind forcing projects onto the ocean’s vertical modes, thereby shaping the amplitudes and characteristics of upper-ocean responses, with oceanic anomalies being dominantly seen at surface layers or in the thermocline. Then, changes in these modal parameters act to modulate characteristics and vertical structure of the upper-ocean dynamical responses to wind forcing.

Representation of mean ocean state in the ICM
It is well recognized that the mean ocean state in the tropical Pacific is a main factor affecting ENSO characteristics. The ICM is formulated in ways that climatological mean ocean fields can be explicitly represented (Supplementary Fig. 6), including stratification, surface velocity and SST, and thermocline properties. Correspondingly, mean ocean fields are explicitly considered to comprise the mean stratification component, the velocity_SST component, and the thermocline component, respectively.

(1) Stratification component. Upper-ocean stratification, represented by the climatological mean buoyancy frequency (N2), is used to perform vertical mode decomposition analyses. It is quantified through vertical structure ψn, projection coefficients (pn), and wave speeds (cn). This component governs the vertical structure of the ocean response to wind forcing, with wind projections onto baroclinic modes (pn) being particularly relevant. The relationship among these parameters can be denoted as N2 ↔ f(ψn, pn, cn). Note that the effects of this stratification component on dynamical ocean responses are not apparent when directly examining the distributions of the stratification-related parameters and fields (ψn and pn), but can be revealed through ICM-based simulations.

(2) Surface velocity-SST component. In the SST anomaly equation (Supplementary Text 2), mean ocean fields of surface velocity (U; its components taken as u, v, w) and SST directly affect the surface-layer heat budget via zonal advection, vertical advection, Ekman transport, and other processes. This mean surface field component is symbolically denoted as U ↔ f(u, v, w, T).

(3) Thermocline component associated with Te. The thermocline component is represented through Te, which is analytically expressed in a hyperbolic tangent (tanh) function (Yuan et al., 2020; Supplementary Text 3) as follows.

Te = (Tr - Ttc) × tanh(Dtc + Dtc' - Hm) / h* - tanh(Dtc - Hm) / h*

in which variables with bar indicate mean fields and with primes are for interannual anomalies; Tr represents the radiative-convective equilibrium temperature (i.e., the upper bound of the warm pool SST in °C); Hm is the mixed layer depth; Dtc is the thermocline depth that is defined by the depth of the maximum temperature vertical gradient; Ttc represents the temperature at the climatological thermocline depth; h* indicates the thermocline sharpness53 (Yuan et al., 2020). This component represents the modulating effects of mean thermocline properties on Te and further on SST. Thermocline properties include its depth, intensity, and sharpness, which collectively modulate the sensitivity of Te to thermocline depth anomalies. So, Te depends on not only thermocline depth anomalies, but also the mean thermocline properties (Dtc, Ttc, and h*). For convenience, the thermocline component associated with mean thermocline effects on Te is symbolically denoted as Te ↔ f(Dtc, Ttc, h*). Examples for the distributions and changes associated with this component are shown in Supplementary Figs. 3 and 4. While the direct effects of this component on subsurface thermal anomalies are not obvious when directly examining the mean fields, they can be clearly illustrated in ICM-based simulations.

Experimental design using the ICM
Mean ocean fields and parameters are estimated from eight selected CMIP6 models spanning 1900–2300 under the SSP585 scenario. Monthly-mean fields from the 8CMM during 1900–2300 are subjected to a 51-year running mean26,27. Mean ocean states (12-month seasonally varying climatological fields and annual mean) are then separately calculated during three representative periods for 1940–1990 (historical), 2040–2090 (mid-to-late 21st century), and 2240–2290 (mid-to-late 23rd century), which are used in the ICM. For the stratification and thermocline components, we use annual-mean fields to calculate the corresponding parameters for each period in ICM simulations, including pn, Dtc, Ttc, and h*. For the surface velocity–SST component, we use seasonally varying climatological wind stress fields from the 8CMM to force the ICM’s dynamical ocean model and derive climatological mixed-layer currents (u, v, w); climatological SST fields (T) are directly taken from the 8CMM.

These mean annual and climatological fields are incorporated into the ICM in a modular manner, enabling targeted sensitivity experiments performed to isolate the effects of each component (Supplementary Fig. 6 and Table 3). The ICM has been freely run for 50 years in each case to capture ENSO responses to different global warming conditions during the three periods. ENSO simulations from targeted ICM experiments can be compared with those from the 8CMM and can be further analyzed to quantify the effects of various components of the mean ocean state. Importantly, these experiments reproduce key features of ENSO amplitude changes seen in CMIP6 simulations, namely weak amplitude in the mid-20th century, enhanced amplitude during the mid-to-late 21st century, and a later decline in the mid-to-late 23rd century, respectively.

Analyses of the Bjerknes stability index (BJ-index)
The Bjerknes stability index (BJ-index) provides a concise framework for assessing ENSO-related feedback processes controlling ENSO amplitude increase or decrease. Developed by Jin et al. (2006) within the recharge–discharge oscillator paradigm and the Zebiak–Cane (ZC) model57, the BJ-index decomposes ENSO growth into distinct positive and negative feedback processes (Supplementary Text 6).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this study, the authors use an intermediate climate model (ICM) of the coupled atmosphere–ocean system to address a paradox in future ENSO projections. Specifically, they show that continuous upper-ocean stratification in the equatorial Pacific tends to increase ENSO amplitude in the early stages of a warming world but suppresses it later on. Their model provides a plausible explanation for this paradox by emphasizing how ongoing stratification affects different vertical modes of ocean response. This manuscript offers valuable insights into future ENSO projections and the key oceanic processes that shape them under changing mean states. From this perspective, the work could be considered for publication in Nature Communications after major revisions. Two areas require substantial improvement: (1) the presentation needs to be clearer and more accessible so readers can easily follow and appreciate the findings, and (2) the suitability of the ICM for this study must be more thoroughly justified, with further clarification of some physical interpretations.\n\nThe presentation of this manuscript is overly complicated and difficult to follow. Much of the critical information is buried in the supplementary files, without being summarized clearly and concisely in the main text. As a result, I found myself having to study supplementary figures, tables, and descriptions just to extract the key points. In many places, the authors vaguely state that “something changes across the three periods” but do not describe what those changes are, instead pointing only to supplementary figure numbers. Essential information should be brought into the main text so that readers can grasp the findings without having to consult the supplementary files in detail. I provide several examples of such omissions below. In addition, figure captions in the main text are poorly connected to the figures themselves, which makes the discussion even harder to follow. Overall, the clarity and presentation of the manuscript need major improvement.\n\nAnother major concern is the lack of detailed discussion about whether the ICM is truly an appropriate tool for studying the physical processes that drive future changes in ENSO amplitude. The atmospheric component of the model is empirically constructed, which raises questions. If the atmospheric model was trained only on present-day climate data, then it may not be valid for future projections. However, if the atmospheric model was separately trained using empirical SVD modes derived from each of the three periods (P1, P2, and P3) based on the mean of the eight selected CMIP6 models, then this concern would be less serious. Unfortunately, the manuscript does not clarify which approach was taken. In addition, the manuscript does not adequately explain why the continued strengthening of stratification from the P2 to P3 periods projects differently onto the vertical ocean modes.\n\nThe specific comments that follow highlight these presentation issues in more detail and also raise questions related to the physical interpretation of the results.\n(1) This study relies heavily on the ICM for its findings and conclusions, which makes the assumptions behind the model particularly important. A key concern lies in the atmospheric component of the ICM, which is based on empirical relationships between surface wind anomalies and SST anomalies derived through SVD analysis. In its current form, the air–sea coupling is fixed using modes constructed from present-day climate and then applied directly to future climate states. This raises questions about whether such couplings are valid under altered mean climate conditions. I suggest that the authors provide a clear physical justification for why they believe these couplings remain appropriate for future projections.\n(2) Related to the previous point, lines 132–137 of the manuscript seem to suggest that the authors may have re-trained or re-constructed all parameters separately for the three periods (P1, P2, and P3). If the atmospheric SVD models were indeed recalculated for each period, meaning that there are three versions of the anomalous atmospheric model in the ICM, then my earlier concern would not apply. If this is the case, I recommend that the authors make this point clear in the manuscript, so that readers understand the basis of the model setup.\n(3) The phrase “CMIP6 MME” in the manuscript appears to refer specifically to the eight CMIP6 models used in this study. However, this wording may mislead readers into thinking it represents the mean of all, or at least most, CMIP6 models. To avoid confusion, I suggest that the authors use a clearer term to indicate that this refers to the mean of the eight selected models rather than the full CMIP6 ensemble.\n(4) It would be helpful if the authors could provide more information on how the eight models used in this study were selected. In the Methods section, the manuscript only states that “eight models that exhibit consistent ENSO responses” were chosen (Supplementary Table 1). However, there is no detailed explanation of what is meant by “consistent ENSO responses” or what specific criteria were applied. I recommend that the authors elaborate on these selection criteria, at least in the Supplementary Information, so readers can better understand the basis for the model choice.\n(5) L96-99: It is not clear why the authors chose the 1940–1990 period to represent a “weak ENSO amplitude” phase. Figure 1a appears to show that ENSO amplitude in the 1900–1940 period is actually weaker than in the 1940–1990 period across all eight models and their mean.\n(6) The term “ENSO amplitude” is not precisely defined in the manuscript. Throughout the text, it is often used interchangeably with the Niño3.4 index or the magnitude of SST anomalies, which are not exactly equivalent to ENSO amplitude. While it is reasonable to use large Niño3.4 values or SST anomalies as proxies for ENSO events, it is important to be explicit about this assumption. For instance, several figure captions begin with “ENSO amplitude,” when in fact they are showing the standard deviation or magnitude of Niño3.4 or equatorial Pacific SST anomalies. To improve accuracy and clarity, I suggest that the authors clearly state near the beginning of the manuscript that they use Niño3.4 and equatorial SST anomalies as proxies for ENSO amplitude, rather than amplitudes of explicitly identified ENSO events.\n(7) Related to the previous comment, it is really the interannual variability of equatorial Pacific SST anomalies or Niño3.4 that should be used to reflect ENSO amplitude. However, the manuscript and Supplementary Information do not mention whether the SST anomalies used in the analysis (e.g., in Figure 2) are based solely on the interannual component, nor do they describe the procedures used to isolate this component. While I think that filtering may not substantially change the overall results, it is scientifically important to state whether such filtering was performed and, if so, how it was done.\n(8) L104-107: This statement is somewhat unclear. I assume it is meant to describe the results shown in Figures 2a and 2b, but it is difficult to understand what it is conveying about the relative contributions of the BJ index terms to the different amplitudes simulated in the P1, P2, and P3 periods. I recommend that the authors make this description more explicit by stating which physical terms of the BJ index primarily determine the ENSO amplitude in each of the three periods. While a detailed discussion of the changing importance of these terms is provided later in subsection (3) of the “Mechanism Analyses,” the description here should at least summarize the key features in Figure 2b and include a note that more details are discussed later in the manuscript.\n(9) L117-127: Please provide brief but specific information on the mean surface velocity, SST, upwelling, and thermocline changes shown in Figure S5 across the three periods. In several places, the manuscript simply refers readers to the supplementary figures without summarizing the key findings they contain. For example, in lines 111–113, the authors note that they use a 51-year running mean from 1900–2300 based on the ensemble mean of CMIP6 simulations and highlight the differences across the three periods (Supplementary Figs. 1–4). However, the text does not explicitly describe what those key differences are. Relying on readers to interpret supplementary figures on their own does not effectively convey the main scientific findings. I suggest that the authors briefly summarize the essential results from the supplementary figures directly in the main texts.\n(10) Figure 3: The ICM appears to produce an unusually regular ENSO cycle. Could the authors explain why this occurs? Is a similarly regular cycle also present in the eight selected CMIP6 models? If not, does this overly regular behavior in the ICM raise concerns that its ENSO dynamics may differ from those in the CMIP6 models? It would strengthen the manuscript if the authors briefly addressed this potential model performance issue.\n(11) L177-181: It is not clear which figures are being referenced to support the findings described in this part of the manuscript. I was unable to locate any figure numbers cited in the text.\n(12) Figure 4 and the related discussion in L188-212: The authors use this figure to suggest that ENSO dynamics shifted across the P1, P2, and P3 periods, and that this shift helps explain the differences in ENSO amplitude among the three periods. The changing amplitudes of the first three vertical modes shown in Figure 4 are particularly interesting and represent an important finding of the study. Could the authors also identify and name the specific ENSO regimes associated with each of these three periods?\n(13) As mentioned earlier, the manuscript does not provide a clear explanation of why the continued strengthening of stratification in the P3 period changes the relative importance of the first three vertical ocean modes. Clarifying this point would make the physical interpretation of the results much stronger.", "reviewer_2": "This study investigates the non-monotonic change in ENSO amplitude under the SSP5-8.5 scenario using a novel intermediate coupled model (ICM) framework that explicitly incorporates mean-state changes derived from CMIP6 simulations. By decomposing upper-ocean density into vertical baroclinic modes and quantifying their wind projection coefficients, the authors show that changes in ocean stratification—particularly in the first two modes—play a dominant role in modulating ENSO variability. The study provides an elegant mechanistic interpretation of how enhanced stratification can lead to the mid-century strengthening and late-century weakening of ENSO amplitude seen in CMIP6 models.\n\nOverall, this is a valuable and timely contribution. The approach of combining CMIP6-based mean-state information with an ICM and baroclinic mode analysis is innovative and effective in elucidating the physical drivers of future ENSO changes. The manuscript is generally well written, and the results are clearly presented. However, there are several issues that should be addressed before the paper can be considered for publication.\n\nMajor Comments\n1: Limited CMIP6 sampling:\nThe analysis appears to be based on only seven CMIP6 models, which raises concerns about representativeness. Given the large inter-model spread in projected ENSO responses, it would be important to verify whether the identified mechanisms hold across a broader model ensemble. At minimum, the authors should justify their model selection and discuss potential biases arising from this limitation.\n2: Mechanism of ENSO weakening after 2100:\nThe proposed mechanism involving the offsetting contributions of the first and second baroclinic modes provides a reasonable explanation for the mid-21st-century ENSO amplification. However, the subsequent weakening after 2100 (e.g., Fig. 5) seems less convincingly captured by this interpretation. Additional discussion or sensitivity tests focusing on other relevant processes—such as the entrainment temperature (Te) effect and changes in mixed-layer depth—would help clarify this aspect.\n\nMinor Comments\n1: Please clarify whether the prescribed mean states are zonally symmetric or fully three-dimensional.\n2: The authors may wish to comment on whether similar stratification-driven mode changes are found in the Indian or Atlantic basins, as a possible avenue for future work.\n\nRecommendation:\nMajor revisions are required to strengthen the robustness and completeness of the results. Addressing the above points will considerably enhance the paper’s impact and clarity.", "reviewer_3": "The present paper analyzes ENSO mean-state interactions using a refined ICM (with baroclinic modes). It aims in particular at explaining the ENSO amplitude increase then reversal in CMIP6 SSP5-8.5 scenario from the strong ocean stratification increase. To assess this, sensitivity experiments are carried out where the ICM mean state is prescribed from a CMIP6 MME on distinct periods. The authors also make use of the BJ-index to quantify ENSO feedback changes on each period. Changes in wind stress projection on baroclinic modes (as driven by a strongly increasing ocean stratification) are found to be the dominant explainer for the ENSO amplitude increase then reversal.\n\nUnfortunately, I recommend the paper to be rejected. I do like this paper and I think the ocean stratification approach is well grounded and meaningful for understanding ENSO in a climate change context. But the approach is not really novel and has already been applied in similar context, the only novelty being, maybe, that it is now applied post 2100 to a decreased period of ENSO variability. This just my understanding (and I could be wrong) that Nature Communications emphasizes more groundbreaking or timely advances, albeit I believe the paper would fit well in other journals. In addition to this, there are some technical issues in the paper (e.g. wrong figures, lack of details) but these seem that they could be addressed.\nDespite my negative recommendation I want to acknowledge to the authors some solid work and an interesting read.\n\nMajor Comments:\n1) The paper needs to state more clearly what is novel in their approach:\n- The stated novelty seems to be a new modeling framework for mean-state ENSO interactions (l24, 42, 128, 323, etc). However the present framework is already well established. For instance the ICM with its singular baroclinic mode approach is from Zhang et al. (2003, 2020) and other ICMs of this kind also exist as cited. The modified Te parameterization may also be slightly new but seems non-essential to the results. Prescribing the ICM mean state from CGCMs is also pretty common.\n- Applying such a framework in a climate change context is also well established. As one notable example, Cai et al. 2018 analyzed in CMIP5 RCP8.5 scenario the interplay between increased stratification, baroclinic wind stress projection (their Fig. 3) and increased ENSO variability for 2000-2099. What is different from their approach here? Is it that ENSO variability is now found to decrease post 2100?\n2) Some Figures are just wrong or uncommented, this needs polish:\n- If I compare Fig. 4 and SI Fig. 9, it doesnt match: p1 decreases over P3 in Fig. 4 but increases in SI Fig. 9.\n- If I compare Fig. 1a and Fig. 2a, it doesnt match. Nino 3.4 SST s.d. is about the same between P1 and P3 in Fig. 1a, yet the SSTA s.d decreases strongly on P3 in Fig 2a (from which I would expect Nino 3.4 SST s.d. to decrease too).\n- Between SI Fig. 2 and 4 the thermocline depths lines are repeated but they are not the same.\n- SI Fig. 5. shows thermocline sharpness, but that sharpness is not defined anywhere. Also how does this relate the mean thermocline depth properties (mean depth Dtc, temperature at mean depth Ttc, and steepness h*, l463)?\n- SI Fig. 8 and 9 are not referenced anywhere in the manuscript or SI. SI. Fig 9 should also show p3.\n- SI Fig 10 I do not understand how this is obtained or and what is the point. This shows (I guess) T(x,z) regressed on central Pacific sea level in CMIP6 vs the ICM. First how do you obtain T(x,z,t) in the ICM? The caption also says T(x,z,t) is obtained using the Te parameterization scheme (SI eq. 3.5) but that is impossible as the scheme is vertically averaged. Second what is point of this diagnosis and the main takeaway?\n\nMinor Comments:\n1) l84 and l361: You select only 8 models out of CMIP6 SSP5-8.5 for the MME. I recommended discussing why you use these models in particular and why not more out of the 40+ available? Is the stratification increase and the ENSO amplitude increase and reversal also found in the other models?\n2) The choice of periods P1 (1940-1990),P2 (2040-2090), P3 (2100-2149) puzzles me. Period P3 is transient with stratification not fully increased and ENSO amplitude just starting to decay (Fig 1). In fact during P3 the ENSO amplitude is still above or similar to the ENSO amplitude during P1 (around 1 s.d. in Fig. 1a), not lower. So why choose P3? Why not consider a period (say P4) centred somewhere on 2200-2300 with more pronounced and representative changes?\n3) The current ICM framework has some limitations as regards its representation of ENSO variabilty, which is fine but should at least be acknowledged somewhere. Looking at Fig. 3 the ICM produces way to regular ENSO oscillations compared to nature, e.g. there is no ENSO asymmetry, irregularity, diversity (e.g. central Pacific CP El Nino) nor extreme events. These features are important in a climate change context too (e.g. Cai et al. 2014 or 2018, Yeh et al. 2009 and other cited references).\n4) The details of stratification changes in SSP5-8.5 during P1, P2, P3 could be more thoroughly discussed paragraph l108, as it motivates the present study. First, you may simply state why stratification increases under climate change. Second, briefly summarize the main take aways from previous studies (beyond just l110: \"previous studies have analyzed these...changes in great details\"). Third, briefly give the main takeaways from SI Figures 1-5 (these figures are not commented anywhere).\n5) The formulation of the ICM lacks some important details:\n- In the Supplementary Information you omitted to relate the SST budget terms (u,v,w in eq 2.1, Dtc in SI eq. 3.4) to the baroclinic modes (un,vn,wn,pn, rho_n in eq. SI 1.6). Note SI eq 2.1 is vertically averaged (within the SST mixed layer) so the u,v,w here are not the same as in SI eq 1.3 (and it is a bit misleading in terms of notations).\n- The thermocline Te component (l451) depends on thermocline mean depth, intensity and sharpness parameters. These mean state parameters are varied between periods P1,P2,P3 (Fig 5.), but it is not stated how they are determined on each period.\n6) Stratification changes modify the projection of wind energy to the (baroclinic) ocean, as considered here, but beyond this it also influences the ocean feedbacks on SST in the mixed-layer budget. For example the amplitude of thermocline displacements for a given energy is modulated by stratification (see e.g. the scln coefficients in Dewitte 2000), and so is the thermocline feebdack (e.g. thermocline mean depth, intensity and sharpness in the present model's Te parameterization l451). Why is this other sensitivity to stratification omitted here? It has been shown to be important including in climate change context (e.g. Carreric et al. 2020).\n7) It is unclear how the ICM mean state parameters (pn, V-SST or Te, SI Table 2) intervene in the BJ index. When interpreting Fig 2d(section l276), only pn contributions are invoked but not V-SST and Te (despite that they do play a role according to Fig. 5). I could suggest (but this is up to you) to provide the BJ-index deduced from each sensitivity experiments (SI Table 3 or Fig. 5), using the same modular approach as Fig. 5.\n\nText:\n- l53: Maybe briefly motivate why you choose the SSP5-8.5 scenario. It is noteworthy the highest emission scenario within CMIP6. Also is the non-monotonic increase (i.e. increase then decrease) in ENSO amplitude unique to this scenario?\n- l120: in its depth and intensity -> in its depth, intensity and sharpness. Because Fig 5 shows sharpness.\n- l188: Figure 4 illustrates... This sentence should be moved to l177 or 178, no?\n- l316 \"This mechanistic explanation resolves a key paradox...\". This is a very nice two sentences that summarizes some important goal of the paper. You chould put these sentences in the introduction too.\n- l468: Extended Data -> Supplementary Information\n- Fig. 1: add vertical lines or colored areas to highlight periods P1, P2, P3\n- SI l110: Just to be clear, are you using exactly Zhang et al. 2020 ICM or some differences? The text suggests that the Te scheme is modified (Supplementary text 3 l333, e.g. the Te-SVD scheme has been replaced by Te-tanh scheme). Maybe add a sentence to summarize what is modified or not. Same for l24 and l128 of main text (the model is described as novel, but what is novel?).\n- SI l157. Equation 5 -> 1.5\n- SI l278: 'h* is related to the vertical gradient of temperature field\"\". Please specify. Steepness is the amplitude of the dT/dz (or N^2) profile at the thermocline depth?\n- SI l314: maximum temperature gradient -> maximum temperature vertical gradient\n- SI l 349: The notation ... (doesnt appear, I guess it is the brackets)\n- SI l366: BJ coefficients in SI Table 4 not 3.\n- SI Fig. 2: Specify the meridional average used to derive this figure. Is it 5N-5S? In panel c it would be more consistent to show thermocline depth P2 (red) and P3 (blue).\n- SI Fig. 2 is the same information at Fig. 2 (although showing totals instead of differences). These figures could be merged (or at least shown one after the other). The thermocline depths should match but they dont between SI Fig. 2 and 4\n- SI Fig. 9, panel a: N->N^2\n- SI Table 2: Some missing parenthesis in the row \"Denotations\".\n- SI Table 4: You could add a column that indicates to which feedback (TH, EK, ZA...) each parameter is associated.", "reviewer_author_response": "Responses to reviewer comments\nReviewer #1 (Remarks to the Author):\nThank you very much for your comments, which are all valuable and very helpful\nfor revising and improving our paper. We have studied the comments very carefully\nand have revised the paper accordingly. The major changes are in red font in the\nrevised manuscript. The main revisions and the responses to specific suggestions are\nas follows.\nIn this study, the authors use an intermediate climate model (ICM) of the coupled\natmosphere–ocean system to address a paradox in future ENSO projections.\nSpecifically, they show that continuous upper-ocean stratification in the equatorial\nPacific tends to increase ENSO amplitude in the early stages of a warming world but\nsuppresses it later on. Their model provides a plausible explanation for this paradox by\nemphasizing how ongoing stratification affects different vertical modes of ocean\nresponse. This manuscript offers valuable insights into future ENSO projections and\nthe key oceanic processes that shape them under changing mean states. From this\nperspective, the work could be considered for publication in Nature Communications\nafter major revisions. Two areas require substantial improvement: (1) the presentation\nneeds to be clearer and more accessible so readers can easily follow and appreciate the\nfindings, and (2) the suitability of the ICM for this study must be more thoroughly\njustified, with further clarification of some physical interpretations.\nThe presentation of this manuscript is overly complicated and difficult to follow.\nMuch of the critical information is buried in the supplementary files, without being\nsummarized clearly and concisely in the main text. As a result, I found myself having\nto study supplementary figures, tables, and descriptions just to extract the key points.\nIn many places, the authors vaguely state that “something changes across the three\nperiods” but do not describe what those changes are, instead pointing only to\nsupplementary figure numbers. Essential information should be brought into the main\ntext so that readers can grasp the findings without having to consult the supplementary\n1\nfiles in detail. I provide several examples of such omissions below. In addition, figure\ncaptions in the main text are poorly connected to the figures themselves, which makes\nthe discussion even harder to follow. Overall, the clarity and presentation of the\nmanuscript need major improvement.\nAnother major concern is the lack of detailed discussion about whether the ICM\nis truly an appropriate tool for studying the physical processes that drive future changes\nin ENSO amplitude. The atmospheric component of the model is empirically\nconstructed, which raises questions. If the atmospheric model was trained only on\npresent-day climate data, then it may not be valid for future projections. However, if\nthe atmospheric model was separately trained using empirical SVD modes derived from\neach of the three periods (P1, P2, and P3) based on the mean of the eight selected\nCMIP6 models, then this concern would be less serious. Unfortunately, the manuscript\ndoes not clarify which approach was taken. In addition, the manuscript does not\nadequately explain why the continued strengthening of stratification from the P2 to P3\nperiods projects differently onto the vertical ocean modes.\nAccording to the two concerns, we have revised the ms very carefully. In particular,\nwe add two new supplementary text 4 & 5 to provide detailed discussion about whether\nthe ICM is truly an appropriate tool for studying the physical processes that drive future\nchanges in ENSO amplitude.\nThe specific comments that follow highlight these presentation issues in more detail\nand also raise questions related to the physical interpretation of the results.\n(1) This study relies heavily on the ICM for its findings and conclusions, which makes\nthe assumptions behind the model particularly important. A key concern lies in the\natmospheric component of the ICM, which is based on empirical relationships between\nsurface wind anomalies and SST anomalies derived through SVD analysis. In its current\nform, the airsea coupling is fixed using modes constructed from present-day climate\nand then applied directly to future climate states. This raises questions about whether\n2\nsuch couplings are valid under altered mean climate conditions. I suggest that the\nauthors provide a clear physical justification for why they believe these couplings\nremain appropriate for future projections.\nAccordingly, we add a new subsection “5.3 Considerations for the suitability\nof the ICM for this study” within the supplementary text 5 (Considerations for\nlimitations in this study) to provide detailed discussion about whether the ICM is truly\nan appropriate tool for studying the physical processes that drive future changes in\nENSO amplitude.\n(2) Related to the previous point, lines 132–137 of the manuscript seem to suggest\nthat the authors may have re-trained or re-constructed all parameters separately for the\nthree periods (P1, P2, and P3). If the atmospheric SVD models were indeed recalculated\nfor each period, meaning that there are three versions of the anomalous atmospheric\nmodel in the ICM, then my earlier concern would not apply. If this is the case, I\nrecommend that the authors make this point clear in the manuscript, so that readers\nunderstand the basis of the model setup.\nAccordingly, we add some texts in the revised version as follows:\nThe ICM is an anomaly model featuring an empirically constructed atmospheric\ncomponent to represent surface wind stress response to SST variability45, which is\ncoupled to a simplified dynamic model 46-47. A key advantage of this framework is its\ncapacity to directly produce interannual anomalies, while oceanic climatologies are\nallowed to be prescribed. Although the ICM has been widely applied in ENSO\nsimulations and predictions44-45,48-49, the present version introduces several key\nadvances (Supplementary Text 4). Notably, the thermocline component is now\npresented through an entrainment temperature (T ) parameterization, analytically\ne\nformulated using a hyperbolic tangent (tanh) function. This formulation explicitly links\nT to mean thermocline properties (i. e., its depth, intensity and sharpness),\ne\nwhich\ncollectively modulate the sensitivity of T to thermocline depth anomalies. Furthermore,\ne\nspatially-varying density and stratification fields are incorporated to derive\ncorresponding modal parameters, enabling a more realistic representation of oceanic\nfeedbacks. Critically, all three influencing components are coherently and explicitly\nintegrated within the unified ICM framework. Then, this enhanced modeling\nframework is innovatively applied to the post-2100 period in CMIP6 simulations,\nwhere a counterintuitive decline in ENSO variability emerges despite continued\nwarming, a phenomenon that existing theories struggle to explain. Through this\napproach, we systematically isolate the relative contributions of mean ocean state\ncomponents to ENSO amplitude changes, offering a mechanistic interpretation of\ncomplex climate projections. Nevertheless, the current ICM framework and its\nmodeling study for ENSO modulation in the global warming context also exhibits\nlimitations in its representation of ENSO variability (Supplementary Text 4).\n(3) The phrase “CMIP6 MME” in the manuscript appears to refer specifically to the\neight CMIP6 models used in this study. However, this wording may mislead readers\ninto thinking it represents the mean of all, or at least most, CMIP6 models. To avoid\nconfusion, I suggest that the authors use a clearer term to indicate that this refers to the\nmean of the eight selected models rather than the full CMIP6 ensemble.\nAccordingly, we have revised the ms carefully, including text revision in the\nMethod section. Furthermore, we add a new subsection “5.1 Issues regarding limited\nsampling of CMIP6 models ” within the supplementary text 5 to provide detailed\ndiscussion.\n(4) It would be helpful if the authors could provide more information on how the\neight models used in this study were selected. In the Methods section, the manuscript\nonly states that “eight models that exhibit consistent ENSO responses” were chosen\n(Supplementary Table 1). However, there is no detailed explanation of what is meant\n4\nby “consistent ENSO responses” or what specific criteria were applied. I recommend\nthat the authors elaborate on these selection criteria, at least in the Supplementary\nInformation, so readers can better understand the basis for the model choice.\nAccordingly, we have revised the ms carefully, including text revision in the\nMethod section. Furthermore, we add a new subsection “5.1 Issues regarding limited\nsampling of CMIP6 models ” within the supplementary text 5 to provide detailed\ndiscussion.\n(5) L96-99: It is not clear why the authors chose the 1940–1990 period to represent a\n“weak ENSO amplitude” phase. Figure 1a appears to show that ENSO amplitude in the\n1900–1940 period is actually weaker than in the 1940–1990 period across all eight\nmodels and their mean.\nAccordingly, we have revised the ms carefully, including text revision in the\nMethod section. Furthermore, we add a new subsection “5.1 Issues regarding limited\nsampling of CMIP6 models ” within the supplementary text 5 to provide detailed\ndiscussion.\n(6) The term “ENSO amplitude” is not precisely defined in the manuscript. Throughout\nthe text, it is often used interchangeably with the Niño3.4 index or the magnitude of\nSST anomalies, which are not exactly equivalent to ENSO amplitude. While it is\nreasonable to use large Niño3.4 values or SST anomalies as proxies for ENSO events,\nit is important to be explicit about this assumption. For instance, several figure captions\nbegin with “ENSO amplitude,” when in fact they are showing the standard deviation or\nmagnitude of Niño3.4 or equatorial Pacific SST anomalies. To improve accuracy and\nclarity, I suggest that the authors clearly state near the beginning of the manuscript that\n5\nthey use Niño3.4 and equatorial SST anomalies as proxies for ENSO amplitude, rather\nthan amplitudes of explicitly identified ENSO events.\nThank you very much for the points. Accordingly, we have revised the ms carefully\nand add the following text:\nTo adequately represent ENSO variability, we use the standard deviation or\nmagnitude of Niño3.4 or equatorial Pacific SST anomalies as proxies for ENSO\namplitude, rather than the amplitudes of explicitly identified ENSO events. Also, we\nconsider interannual variability of Niño3.4 or equatorial Pacific SST that are used to\nreflect ENSO amplitude. So, SST anomalies used in the following analyses are based\nsolely on the interannual component, which is obtained by subtracting the 51-year\nrunning mean climatology for each selected CMIP6 model.\n(7) Related to the previous comment, it is really the interannual variability of equatorial\nPacific SST anomalies or Niño3.4 that should be used to reflect ENSO amplitude.\nHowever, the manuscript and Supplementary Information do not mention whether the\nSST anomalies used in the analysis (e.g., in Figure 2) are based solely on the interannual\ncomponent, nor do they describe the procedures used to isolate this component. While\nI think that filtering may not substantially change the overall results, it is scientifically\nimportant to state whether such filtering was performed and, if so, how it was done.\nAccordingly, we have revised the ms carefully and add the following text:\nTo adequately represent ENSO variability, we use the standard deviation or\nmagnitude of Niño3.4 or equatorial Pacific SST anomalies as proxies for ENSO\namplitude, rather than the amplitudes of explicitly identified ENSO events. Also, we\nconsider interannual variability of Niño3.4 or equatorial Pacific SST that are used to\nreflect ENSO amplitude. So, SST anomalies used in the following analyses are based\n6\nsolely on the interannual component, which is obtained by subtracting the 51-year\nrunning mean climatology for each selected CMIP6 model.\n(8) L104-107: This statement is somewhat unclear. I assume it is meant to describe the\nresults shown in Figures 2a and 2b, but it is difficult to understand what it is conveying\nabout the relative contributions of the BJ index terms to the different amplitudes\nsimulated in the P1, P2, and P3 periods. I recommend that the authors make this\ndescription more explicit by stating which physical terms of the BJ index primarily\ndetermine the ENSO amplitude in each of the three periods. While a detailed discussion\nof the changing importance of these terms is provided later in subsection (3) of the\n“Mechanism Analyses,” the description here should at least summarize the key features\nin Figure 2b and include a note that more details are discussed later in the manuscript.\nThank you very much for your suggestions. Accordingly, we have revised the ms\ncarefully and place the related texts to the “Mechanism Analyses” section, where\nbudget analyses for CMIP6 and ICM simulations are analyzed.\n(9) L117-127: Please provide brief but specific information on the mean surface\nvelocity, SST, upwelling, and thermocline changes shown in Figure S5 across the three\nperiods. In several places, the manuscript simply refers readers to the supplementary\nfigures without summarizing the key findings they contain. For example, in lines 111–\n113, the authors note that they use a 51-year running mean from 1900–2300 based on\nthe ensemble mean of CMIP6 simulations and highlight the differences across the three\nperiods (Supplementary Figs. 1–4). However, the text does not explicitly describe what\nthose key differences are. Relying on readers to interpret supplementary figures on their\nown does not effectively convey the main scientific findings. I suggest that the authors\nbriefly summarize the essential results from the supplementary figures directly in the\nmain texts.\nThank you very much for your suggestions. Accordingly, we have revised the ms\ncarefully and add more texts as follows.\nsimulations and highlight their differences for the three periods (Supplementary Figs.\n1-3). There are apparent changes in SST and surface winds, characterized by reduced\nsurface currents and zonal SST gradients. As a result, the equatorial upwelling is\nreduced in association with the changing mean surface velocity and SST fields. In\naddition, the thermocline shows pronounced changes in its depth, intensity and\nsharpness (Supplementary Fig. 4). For example, the thermocline continues to shoal in\nthe equatorial Pacific, with its increased sharpness. In particular, the stratification\nexhibits non-uniform vertical changes in the upper ocean, with an increase in the upper\n120 m, but a decrease below (Supplementary Fig. 5).\n(10) Figure 3: The ICM appears to produce an unusually regular ENSO cycle. Could\nthe authors explain why this occurs? Is a similarly regular cycle also present in the eight\nselected CMIP6 models? If not, does this overly regular behavior in the ICM raise\nconcerns that its ENSO dynamics may differ from those in the CMIP6 models? It would\nstrengthen the manuscript if the authors briefly addressed this potential model\nperformance issue.\nThank you very much for your comments. Accordingly, we add a new subsection\n“5.2 Considerations for model performance” within the supplementary text 5\n(Considerations for limitations in this study) to provide detailed discussion about\nwhether the ICM is truly an appropriate tool for studying the physical processes that\ndrive future changes in ENSO amplitude.\n(11) L177-181: It is not clear which figures are being referenced to support the findings\ndescribed in this part of the manuscript. I was unable to locate any figure numbers cited\nin the text.\nThank you very much for pointing this. Accordingly, we have revised the ms\ncarefully and move some texts around. Now, the texts read as\nFigure 4 illustrates the spatial distributions of p₁, p₂, and p₃ across the equatorial\nPacific for each period. These parameters exhibit well-defined spatial structure across\n8\nthe equatorial Pacific. The amplitude of the first baroclinic mode (p₁) decreases\neastward along the equator, indicating a more surface-intensified response in the\nwestern-central equatorial Pacific. Conversely, p₂ and p₃ amplitudes increase toward the\neast, reflecting deeper thermocline sensitivities in the eastern equatorial Pacific. Using\nthese modal parameters, reconstruction calculations for upper-ocean variability indicate\nthat the first three baroclinic modes dominate the ocean’s dynamical responses to wind\nforcing.\n(12) Figure 4 and the related discussion in L188-212: The authors use this figure to\nsuggest that ENSO dynamics shifted across the P1, P2, and P3 periods, and that this\nshift helps explain the differences in ENSO amplitude among the three periods. The\nchanging amplitudes of the first three vertical modes shown in Figure 4 are particularly\ninteresting and represent an important finding of the study. Could the authors also\nidentify and name the specific ENSO regimes associated with each of these three\nperiods?\nThank you very much for pointing this. Accordingly, we have revised the ms\ncarefully, by adding a new subsection “(3) Further mechanistic insights into post-\n5), with corresponding texts. Some additional texts including the following:\nOne major finding from this study is the relatively changing amplitude of the first\nthree vertical baroclinic modes across the three periods, as derived from modal\ndecomposition of stratification. While wave speeds (cₙ) remain relatively stable, the\nwind projection coefficients (pₙ) undergo pronounced temporal changes across P1, P2,\nand P3, which determine the vertical structure and spatial patterning of the ocean’s\nresponses to wind forcing.\n(13) As mentioned earlier, the manuscript does not provide a clear explanation of why\nthe continued strengthening of stratification in the P3 period changes the relative\nimportance of the first three vertical ocean modes. Clarifying this point would make the\n9\nphysical interpretation of the results much stronger.\nThank you very much for pointing this. Accordingly, we have revised the ms\ncarefully, by adding a new subsection “(3) Further mechanistic insights into post-\n5), with corresponding texts. Some additional texts including the following:\n(3) Further mechanistic insights into post-2100 ENSO weakening\nThe proposed mechanism invokes changes in the wind projection coefficients\nderived from stratification and their amplifying effects on ENSO, which are\naccompanied with damping effects from the thermocline and surface field components.\nThe increases in the amplitudes of p and p offer a plausible explanation for the\n1 2\namplification of ENSO until the mid-21st century (during period P2). However, its\nability to account for the subsequent weakening beyond 2100 (e.g., Fig. 5) appears less\ncompelling. A particular challenge lies in explaining why continued stratification\nstrengthening from period P2 to P3 differentially projects onto ocean vertical baroclinic\nmodes (primarily p₁ and p₂). This subsection aims to clarify this issue through additional\nanalyses and discussions, thereby strengthening the physical interpretation of ENSO\nchanges under prolonged warming.\nAs shown in Fig. 1c, oceanic stratification in the upper ocean undergoes marked\nchanges with progressive warming. Closer examination reveals that these changes are\nnon-uniform in both vertical and horizontal directions. For instance, the vertical\nstructure of stratification (Supplementary Fig. 5) indicates intensification in the upper\n1\n10\n120 m, but a reduction at subsurface depths between 120–150 m. This vertical\nreconstruction of the stratification alters the relative contributions of wind projections\nto the baroclinic mode amplitudes: the amplitude of p₁ decreases while that of p₂\nincreases during P3 (Supplementary Fig. 4), reflecting an amplifying of surface-layer\nstratification but vertically a gradual weakening of subsurface stratification. From P2\nto P3, this shift in the vertical distribution of stratification modifies the respective\ninfluences of p₁ and p₂ on upper-ocean responses. During P3, the reduced amplifying\neffects on ocean anomalies due to the decreasced p amplitude are not sufficiently\n1\nstrong to offset the increased amplifying effects due to increased p amplitude, making\n2\nthe stratification-related net amplifying effects reduced.\nIt should be noted that other processes associated with mean ocean fields also\ncontribute to ENSO modulation in the coupled system. As indicated in Fig. 5, the T\ne\nmodulation due to changes in thermocline depth, strength, and sharpness exerts a\ndamping influence on ENSO variability. Similarly, changes in the mean surface ocean\nfields also play a suppressive role (Fig. 5). Thus, in the ICM framework, three key\ncomponents coexist that represent their competing effects on ENSO: the thermocline\nand surface field components both dampen ENSO, whereas the stratification generally\namplifies it, which is determined by the respective amplitudes of p₁ and p₂. (that is,\nwhen both p₁ and p₂ increase, ENSO intensifies). Thus, the overall ENSO amplitude is\ngoverned by the combined and synergistic effects of these influencing components, and\nthe resultant ENSO amplitude reflects a balance between amplifying and damping\neffects in the coupled system. In particular, during P3, p₁ decreases while p₂ increases,\nresulting in a partial cancellation of their amplifying contributions. Concurrently,\ndamping mechanisms remain active. When the thermocline- and surface fields-related\ndamping effects dominate, the amplification due to changes in stratification is\ninsufficient to offset the damping, the net effect leads to a reduction in ENSO\namplitude—as is the case during P3.\nThank you very much for your comments, which are all valuable and very helpful\nfor revising and improving the paper. We have studied the comments very carefully\nand have revised our paper accordingly. The major changes are in red font in the\nrevised manuscript. The main revisions and the responses to the suggestions are as\nfollows.\nThis study investigates the non-monotonic change in ENSO amplitude under the\nSSP5-8.5 scenario using a novel intermediate coupled model (ICM) framework that\nexplicitly incorporates mean-state changes derived from CMIP6 simulations. By\ndecomposing upper-ocean density into vertical baroclinic modes and quantifying their\nwind projection coefficients, the authors show that changes in ocean stratification—\nparticularly in the first two modes—play a dominant role in modulating ENSO\nvariability. The study provides an elegant mechanistic interpretation of how enhanced\nstratification can lead to the mid-century strengthening and late-century weakening of\nENSO amplitude seen in CMIP6 models.\nOverall, this is a valuable and timely contribution. The approach of combining\nCMIP6-based mean-state information with an ICM and baroclinic mode analysis is\ninnovative and effective in elucidating the physical drivers of future ENSO changes.\nThe manuscript is generally well written, and the results are clearly presented. However,\nthere are several issues that should be addressed before the paper can be considered for\npublication.\nMajor Comments\n1: Limited CMIP6 sampling:\nThe analysis appears to be based on only seven CMIP6 models, which raises concerns\nabout representativeness. Given the large inter-model spread in projected ENSO\nresponses, it would be important to verify whether the identified mechanisms hold\nacross a broader model ensemble. At minimum, the authors should justify their model\nselection and discuss potential biases arising from this limitation.\nThank you very much for pointing out this. Accordingly, we have revised the ms\n12\ncarefully, including text revision in the Method section. Furthermore, we add a new\nsubsection “5.1 Issues regarding limited sampling of CMIP6 models ” within the\nsupplementary text 5 to provide detailed discussion. Some texts are included as follows.\nTo investigate ENSO responses to global warming, we analyze CMIP6 historical\nand future simulations under the high-emission SSP5-8.5 scenario, following\nmethodologies established in previous studies27. These simulations are further used to\nderive climatological mean oceanic fields, which are subsequently prescribed in an\nICM to quantitatively assess their influence on ENSO amplitude.\nWe select eight CMIP6 models (Supplementary Table 1) based on their consistent\nprojection of a distinct non-monotonic ENSO amplitude changes under sustained\nwarming, a feature not widely present across the full CMIP6 ensemble. The eight\nCMIP6 model mean (8CMM) is characterized by three phases in terms of ENSO\namplitudes: a weak period (1940–1990, P1), an enhanced period (2040–2090, P2), and\na subsequent declining period (2100–2149, P3), respectively. It should be noted that\nthe mechanisms identified using this subset of models require validation across a\nbroader ensemble, given the considerable intermodel spread in ENSO projections. A\ndetailed discussion of sampling considerations is further provided in Supplementary\nText 5.\nThe selection of these periods, while consistent with prior work27, reflects key\ntransitional behaviors in both ENSO amplitude and oceanic mean state. Notably, the\nchanges in mean stratification associated modal parameters (pₙ) and ENSO amplitude\nare not synchronous—each peaks and declines with distinct timing (Fig. 1 and\nSupplementary data Fig. 5). The P3 period, in particular, represents a transitional phase\nduring which ENSO amplitude begins to decline around 2100 (Fig. 1), while\nstratification continues to intensify and the first baroclinic mode (p₁) peaks near 2150\nbefore decreasing. Should a later period (say 2200–2250) be considered, ENSO\namplitude further declines under intensified stable stratification and reduced p₁, yet the\nqualitative relationship among these variables remains consistent. These results suggest\n13\nthat the physical processes identified based on the currently selected periods are robust\nacross post-2100 climates, even as the mean state continues to evolve.\n2: Mechanism of ENSO weakening after 2100:\nThe proposed mechanism involving the offsetting contributions of the first and second\nbaroclinic modes provides a reasonable explanation for the mid-21st-century ENSO\namplification. However, the subsequent weakening after 2100 (e.g., Fig. 5) seems less\nconvincingly captured by this interpretation. Additional discussion or sensitivity tests\nfocusing on other relevant processes—such as the entrainment temperature (Te) effect\nand changes in mixed-layer depth—would help clarify this aspect.\nThank you very much for pointing this. Accordingly, we have revised the ms\ncarefully, by adding a new subsection “(3) Further mechanistic insights into post-\n5), with corresponding texts. Some additional texts including the following:\n(3) Further mechanistic insights into post-2100 ENSO weakening\nThe proposed mechanism invokes changes in the wind projection coefficients\nderived from stratification and their amplifying effects on ENSO, which are\naccompanied with damping effects from the thermocline and surface field components.\nThe increases in the amplitudes of p and p offer a plausible explanation for the\n1 2\namplification of ENSO until the mid-21st century (during period P2). However, its\nability to account for the subsequent weakening beyond 2100 (e.g., Fig. 5) appears less\ncompelling. A particular challenge lies in explaining why continued stratification\nstrengthening from period P2 to P3 differentially projects onto ocean vertical baroclinic\nmodes (primarily p₁ and p₂). This subsection aims to clarify this issue through additional\nanalyses and discussions, thereby strengthening the physical interpretation of ENSO\nchanges under prolonged warming.\nAs shown in Fig. 1c, oceanic stratification in the upper ocean undergoes marked\nchanges with progressive warming. Closer examination reveals that these changes are\nnon-uniform in both vertical and horizontal directions. For instance, the vertical\nstructure of stratification (Supplementary Fig. 5) indicates intensification in the upper\n14\n120 m, but a reduction at subsurface depths between 120–150 m. This vertical\nreconstruction of the stratification alters the relative contributions of wind projections\nto the baroclinic mode amplitudes: the amplitude of p₁ decreases while that of p₂\nincreases during P3 (Supplementary Fig. 4), reflecting an amplifying of surface-layer\nstratification but vertically a gradual weakening of subsurface stratification. From P2\nto P3, this shift in the vertical distribution of stratification modifies the respective\ninfluences of p₁ and p₂ on upper-ocean responses. During P3, the reduced amplifying\neffects on ocean anomalies due to the decreasced p amplitude are not sufficiently\n1\nstrong to offset the increased amplifying effects due to increased p amplitude, making\n2\nthe stratification-related net amplifying effects reduced.\nIt should be noted that other processes associated with mean ocean fields also\ncontribute to ENSO modulation in the coupled system. As indicated in Fig. 5, the T\ne\nmodulation due to changes in thermocline depth, strength, and sharpness exerts a\ndamping influence on ENSO variability. Similarly, changes in the mean surface ocean\nfields also play a suppressive role (Fig. 5). Thus, in the ICM framework, three key\ncomponents coexist that represent their competing effects on ENSO: the thermocline\nand surface field components both dampen ENSO, whereas the stratification generally\namplifies it, which is determined by the respective amplitudes of p₁ and p₂. (that is,\nwhen both p₁ and p₂ increase, ENSO intensifies). Thus, the overall ENSO amplitude is\ngoverned by the combined and synergistic effects of these influencing components, and\nthe resultant ENSO amplitude reflects a balance between amplifying and damping\neffects in the coupled system. In particular, during P3, p₁ decreases while p₂ increases,\nresulting in a partial cancellation of their amplifying contributions. Concurrently,\ndamping mechanisms remain active. When the thermocline- and surface fields-related\ndamping effects dominate, the amplification due to changes in stratification is\ninsufficient to offset the damping, the net effect leads to a reduction in ENSO\namplitude—as is the case during P3.\n1: Please clarify whether the prescribed mean states are zonally symmetric or fully\nthree-dimensional.\nRevised as suggested, indicating that the prescribed mean states are fully three-\ndimensional.\n2: The authors may wish to comment on whether similar stratification-driven mode\nchanges are found in the Indian or Atlantic basins, as a possible avenue for future work.\nRevised as suggested. Some related texts are included as follows.\nLooking forward, this mechanistic framework opens several promising research\navenues. Future work could incorporate transient background climate states and\nstochastic atmospheric component to assess robustness across more realistic climate\nregimes. Moreover, while this study focuses exclusively on ENSO modulations in the\ntropical Pacific, examining analogous processes in the Indian and Atlantic Oceans may\noffer valuable insights for understanding climate variability in these basins,\nrepresenting a promising research direction that could advance understanding of global\nclimate variability patterns and potentially guiding the development of next-generation\nregional climate projections.\nRecommendation:\nMajor revisions are required to strengthen the robustness and completeness of the\nresults. Addressing the above points will considerably enhance the paper’s impact and\nclarity.\nThe present paper analyzes ENSO mean-state interactions using a refined ICM (with\nbaroclinic modes). It aims in particular at explaining the ENSO amplitude increase then\nreversal in CMIP6 SSP5-8.5 scenario from the strong ocean stratification increase. To\nassess this, sensitivity experiments are carried out where the ICM mean state is\nprescribed from a CMIP6 MME on distinct periods. The authors also make use of the\nBJ-index to quantify ENSO feedback changes on each period. Changes in wind stress\nprojection on baroclinic modes (as driven by a strongly increasing ocean stratification)\nare found to be the dominant explainer for the ENSO amplitude increase then reversal.\nUnfortunately, I recommend the paper to be rejected. I do like this paper and I think\nthe ocean stratification approach is well grounded and meaningful for understanding\nENSO in a climate change context. But the approach is not really novel and has already\nbeen applied in similar context, the only novelty being, maybe, that it is now applied\npost 2100 to a decreased period of ENSO variability. This just my understanding (and\nI could be wrong) that Nature Communications emphasizes more groundbreaking or\ntimely advances, albeit I believe the paper would fit well in other journals. In addition\nto this, there are some technical issues in the paper (e.g. wrong figures, lack of details)\nbut these seem that they could be addressed. Despite my negative recommendation I\nwant to acknowledge to the authors some solid work and an interesting read.\nThank you very much for your comments, which are all valuable and very helpful\nfor revising and improving the paper. We have studied the comments very carefully\nand have revised our paper accordingly. In particular, we add a new subsection\n“Supplementary Text 4: Methodological and scientific innovations in this study ”\nto provide detailed discussions in terms of technical and scientific advances that we\nbelieve they are new. The major changes are in red font in the revised manuscript.\nThe main revisions and the responses to the suggestions are as follows.\nMajor Comments:\n1) The paper needs to state more clearly what is novel in their approach:\n17\n- The stated novelty seems to be a new modeling framework for mean-state ENSO\ninteractions (l24, 42, 128, 323, etc). However the present framework is already well\nestablished. For instance the ICM with its singular baroclinic mode approach is from\nZhang et al. (2003, 2020) and other ICMs of this kind also exist as cited. The modified\nTe parameterization may also be slightly new but seems non-essential to the results.\nPrescribing the ICM mean state from CGCMs is also pretty common.\n- Applying such a framework in a climate change context is also well established. As\none notable example, Cai et al. 2018 analyzed in CMIP5 RCP8.5 scenario the\ninterplay between increased stratification, baroclinic wind stress projection (their Fig.\n3) and increased ENSO variability for 2000-2099. What is different from their\napproach here? Is it that ENSO variability is now found to decrease post 2100?\nThank you very much for your comments. Accordingly, we add a new subsection\n“Supplementary Text 4: Methodological and scientific innovations in this study ”\nto provide detailed discussions in terms of technical and scientific advances. Also,\nsome texts in the introductions are revised to highlight the motivation of this study.\nSome examples include the following:\nAs a result, many previous assessments have been largely qualitative and\ndiagnostic in nature. Although often mentioned as a key influence in the literatures,\nstratification effects have rarely been represented and quantified in a consistent or\nphysically explicit way26, 27, 39-42. In particular, the role played by upper-ocean\nstratification in modulating ENSO under global warming remains elusive.\n“Supplementary Text 4: Methodological and scientific innovations in this study ” is\nattached at the end of this response for your convenience.\n2) Some Figures are just wrong or uncommented, this needs polish:\n- If I compare Fig. 4 and SI Fig. 9, it doesnt match: p1 decreases over P3 in Fig. 4 but\nincreases in SI Fig. 9.\n- If I compare Fig. 1a and Fig. 2a, it doesnt match. Nino 3.4 SST s.d. is about the same\nbetween P1 and P3 in Fig. 1a, yet the SSTA s.d decreases strongly on P3 in Fig 2a (from\n18\nwhich I would expect Nino 3.4 SST s.d. to decrease too).\n- Between SI Fig. 2 and 4 the thermocline depths lines are repeated but they are not the\nsame.\n- SI Fig. 5. shows thermocline sharpness, but that sharpness is not defined anywhere.\nAlso how does this relate the mean thermocline depth properties (mean depth Dtc,\ntemperature at mean depth Ttc, and steepness h*, l463)?\n- SI Fig. 8 and 9 are not referenced anywhere in the manuscript or SI. SI. Fig 9 should\nalso show p3.\n- SI Fig 10 I do not understand how this is obtained or and what is the point. This shows\n(I guess) T(x,z) regressed on central Pacific sea level in CMIP6 vs the ICM. First how\ndo you obtain T(x,z,t) in the ICM? The caption also says T(x,z,t) is obtained using the\nTe parameterization scheme (SI eq. 3.5) but that is impossible as the scheme is vertically\naveraged. Second what is point of this diagnosis and the main takeaway?\nThank you very much for pointing these. Accordingly, we have carefully check these\nand revised the ms. The main reasons for the figure mismatch are due to differences in\nthe averaged regions: the Nino3.4 region is defined as (170°W - 120°W and 5°S - 5°N),\nwhereas some fields along the equator are defined as be averaged in 2°S - 2°N and the\ncentral equatorial Pacific (CEP) is referred to be averaged in the region (160°E - 160°W,\n2° S - 2° N). We have re-organized the related figures (including the deletion of the\noriginal Figs. 8-10 and replacements of something new) and add some texts to explain\nhow the mean thermocline depth properties (mean depth Dtc, temperature at mean\ndepth Ttc, and steepness h*) are defined and calculated. For example, some texts are\nincluded as follows.\nIn the SST equation, subsurface entrainment temperature into the ML (T ) affects\ne\nSST through mean vertical upwelling and mixing, which is an essential component of\nSST variability. Previously, a statistical T model is developed based on empirical\ne\nrelationships between anomalies of T and thermocline depth using SVD analysis32,\ne\nwhich has been widely used for ENSO modulation studies48,49,57. However, the\nmodulating effects of the mean thermocline parameters (its mean depth and sharpness)\non T have not been explicitly represented in the SVD-based T parameterizations. It is\ne e\n19\ndesirable to have explicit relationships of T with the thermocline depth anomaly, as\ne\nwell as its mean parameters (including depth and sharpness). In this work, , following\nthe work of Yuan et al. (2020), we take T as an analytical function linking it with\ne\nthermocline depth anomalies and mean thermocline properties in a hyperbolic tangent\n(tanh) form.\nIn our ICM simulations, mean fields for 𝐷‾ , ℎ̅̅̅∗ and 𝑇‾ are prescribed from\nt𝑐 𝑡𝑐\nthe 8CMM data; ocean anomaly fields are produced from the dynamical model using\nthe governing equation during time integration, including sea level anomalies (h´) and\ndynamical pressure anomalies. To relate h´to 𝐷 ， for uses in Equation (3.5) their\nt𝑐\nempirical relationship is derived from historical ocean reanalysis data using a regression\n，\nmethod; the derived 𝐷 fields are used in Equation (3.3) to derive T e . During the three\nt𝑐\n，\nperiods, the regression-based relationship between h´ and 𝐷 are kept the same in\nt𝑐\nICM experiments. In addition, 𝑇‾ and 𝐷‾ are estimated from the 8CMM, and ℎ̅̅̅∗\nt𝑐 𝑡𝑐\nis analytically calculated using Equation (3.2) from the 8CMM data; these mean fields\nact to modulate the sensitivity of subsurface temperature to thermocline depth\nvariations in the ICM, accompanied with stratification effects represented.\nMinor Comments:\n1) l84 and l361: You select only 8 models out of CMIP6 SSP5-8.5 for the MME. I\nrecommended discussing why you use these models in particular and why not more out\nof the 40+ available? Is the stratification increase and the ENSO amplitude increase and\nreversal also found in the other models?\nThank you very much for pointing this. Accordingly, we have revised the Method\nsection. Furthermore, we add a new subsection “5.1 Issues regarding limited\nsampling of CMIP6 models ” within the supplementary text 5 to provide detailed\ndiscussion. Some texts are revised as follows.\nWe select eight CMIP6 models (Supplementary Table 1) based on their consistent\nprojection of a distinct non-monotonic ENSO amplitude changes under sustained\n20\nwarming, a feature not widely present across the full CMIP6 ensemble. The eight\nCMIP6 model mean (8CMM) is characterized by three phases in terms of ENSO\namplitudes: a weak period (1940–1990, P1), an enhanced period (2040–2090, P2), and\na subsequent declining period (2100–2149, P3), respectively. It should be noted that\nthe mechanisms identified using this subset of models require validation across a\nbroader ensemble, given the considerable intermodel spread in ENSO projections. A\ndetailed discussion of sampling considerations is further provided in Supplementary\nText 5.\n2) The choice of periods P1 (1940-1990),P2 (2040-2090), P3 (2100-2149) puzzles me.\nPeriod P3 is transient with stratification not fully increased and ENSO amplitude just\nstarting to decay (Fig 1). In fact during P3 the ENSO amplitude is still above or similar\nto the ENSO amplitude during P1 (around 1 s.d. in Fig. 1a), not lower. So why choose\nP3? Why not consider a period (say P4) centred somewhere on 2200-2300 with more\npronounced and representative changes?\nThank you very much for pointing this. Accordingly, we have revised the Method\nsection. Furthermore, we add a new subsection “5.1 Issues regarding limited\nsampling of CMIP6 models ” within the supplementary text 5 to provide detailed\ndiscussion. Some texts are revised as follows.\nThe selection of these periods, while consistent with prior work27, reflects key\ntransitional behaviors in both ENSO amplitude and oceanic mean state. Notably, the\nchanges in mean stratification associated modal parameters (pₙ) and ENSO amplitude\nare not synchronous—each peaks and declines with distinct timing (Fig. 1 and\nSupplementary data Fig. 5). The P3 period, in particular, represents a transitional phase\nduring which ENSO amplitude begins to decline around 2100 (Fig. 1), while\nstratification continues to intensify and the first baroclinic mode (p₁) peaks near 2150\nbefore decreasing. Should a later period (say 2200–2250) be considered, ENSO\namplitude further declines under intensified stable stratification and reduced p₁, yet the\nqualitative relationship among these variables remains consistent. These results suggest\n21\nthat the physical processes identified based on the currently selected periods are robust\nacross post-2100 climates, even as the mean state continues to evolve.\n3) The current ICM framework has some limitations as regards its representation of\nENSO variabilty, which is fine but should at least be acknowledged somewhere.\nLooking at Fig. 3 the ICM produces way to regular ENSO oscillations compared to\nnature, e.g. there is no ENSO asymmetry, irregularity, diversity (e.g. central Pacific CP\nEl Nino) nor extreme events. These features are important in a climate change context\ntoo (e.g. Cai et al. 2014 or 2018, Yeh et al. 2009 and other cited references).\nThank you for your comments. Accordingly, we have revised the ms very carefully.\nIn particular, we add a new supplementary text 5 “Considerations for limitations in this\nstudy”, with the subsection named as “5.2 Considerations for model performance ”.\n4) The details of stratification changes in SSP5-8.5 during P1, P2, P3 could be more\nthoroughly discussed paragraph l108, as it motivates the present study. First, you may\nsimply state why stratification increases under climate change. Second, briefly\nsummarize the main take aways from previous studies (beyond just l110: \"previous\nstudies have analyzed these...changes in great details\"). Third, briefly give the main\ntakeaways from SI Figures 1-5 (these figures are not commented anywhere).\nThank you for your comments. Accordingly, we have revised the ms very carefully.\nIn particular, we add a new subsection “(3) Further mechanistic insights into post-\n5), with corresponding texts. Some additional texts including the following:\nOne major finding from this study is the relatively changing amplitude of the first\nthree vertical baroclinic modes across the three periods, as derived from modal\ndecomposition of stratification. While wave speeds (cₙ) remain relatively stable, the\nwind projection coefficients (pₙ) undergo pronounced temporal changes across P1, P2,\nand P3, which determine the vertical structure and spatial patterning of the ocean’s\nresponses to wind forcing.\n22\n5) The formulation of the ICM lacks some important details:\n- In the Supplementary Information you omitted to relate the SST budget terms (u,v,w\nin eq 2.1, Dtc in SI eq. 3.4) to the baroclinic modes (un,vn,wn,pn, rho_n in eq. SI 1.6).\nNote SI eq 2.1 is vertically averaged (within the SST mixed layer) so the u,v,w here are\nnot the same as in SI eq 1.3 (and it is a bit misleading in terms of notations).\n- The thermocline Te component (l451) depends on thermocline mean depth, intensity\nand sharpness parameters. These mean state parameters are varied between periods\nP1,P2,P3 (Fig 5.), but it is not stated how they are determined on each period.\nThank you for your comments. Accordingly, we have revised the ms very carefully.\nIn particular, we add some new texts as follows:\nFrom equations (1.3) and (1.6), current anomalies (u´, v´, w´) are calculated, which\nare then averaged within the mixed layer for uses in calculating the SST anomalies\nusing Equation (2.1). Note that u´, v´, w´ in Equation (2.1) are indicated for vertically\naveraged fields within the mixed layer and thus are not the same as in Equation (1.3).\nThe same notations adopted in Equations (1.3) and (2.1) are a bit misleading but simply\nare for conciseness and convenience.\nMean ocean fields and parameters are estimated from eight selected CMIP6\nmodels spanning 1900–2300 under the SSP5-8.5 scenario. Monthly-mean fields from\nthe 8CMM during 1900–2300 are subjected to a 51-year running mean26,27. Mean ocean\nstates (12-month seasonally varying climatological fields and annual mean) are then\nseparately calculated during three representative periods for 1940–1990 (historical),\n2040–2090 (mid-to-late 21st century), and 2100–2149 (early 22nd century), which are\nused in the ICM. For the stratification and thermocline components, we use annual-\nmean fields to calculate the corresponding parameters for each period in ICM\nsimulations, including p , D̅ , T̅ , and h̅*. For the surface velocity–SST component, we\nn tc tc\nuse seasonally varying climatological wind stress fields from the 8CMM to force the\nICM’s dynamical ocean model and derive climatological mixed-layer currents (u̅, v̅, w̅ );\nclimatological SST fields (T̅ ) are directly taken from the 8CMM.\n23\n6) Stratification changes modify the projection of wind energy to the (baroclinic)\nocean, as considered here, but beyond this it also influences the ocean feedbacks on SST in the\nmixed-layer budget. For example the amplitude of thermocline displacements for a\ngiven energy is modulated by stratification (see e.g. the scln coefficients in Dewitte\n2000), and so is the thermocline feebdack (e.g. thermocline mean depth, intensity and\nsharpness in the present model's Te parameterization l451). Why is this other sensitivity\nto stratification omitted here? It has been shown to be important including in climate\nchange context (e.g. Carreric et al. 2020).\nThank you for your comments. Accordingly, we have revised the ms very carefully.\nIn particular, we add some new texts as follows:\nIn addition, the stratification effects and changes are accompanied with other effects\nassociated with changes in surface fields and the thermocline (Supplementary Figs. 1-\n5). So, there exist interactions among these components, with their net effects on ENSO\nbeing complex and intertwined, making it difficult to rigorously isolate the individual\ncontributions of each factor.\n7) It is unclear how the ICM mean state parameters (pn, V-SST or Te, SI Table 2)\nintervene in the BJ index. When interpreting Fig 2d(section l276), only pn contributions\nare invoked but not V-SST and Te (despite that they do play a role according to Fig. 5).\nI could suggest (but this is up to you) to provide the BJ-index deduced from each\nsensitivity experiments (SI Table 3 or Fig. 5), using the same modular approach as Fig.\n5.\nThank you for your comments. Accordingly, we have revised the ms very carefully.\nIn particular, we add some new texts as follows:\nMean ocean fields and parameters are estimated from eight selected CMIP6\nmodels spanning 1900–2300 under the SSP5-8.5 scenario. Monthly-mean fields from\nthe 8CMM during 1900–2300 are subjected to a 51-year running mean26,27. Mean ocean\nstates (12-month seasonally varying climatological fields and annual mean) are then\n24\nseparately calculated during three representative periods for 1940–1990 (historical),\n2040–2090 (mid-to-late 21st century), and 2100–2149 (early 22nd century), which are\nused in the ICM. For the stratification and thermocline components, we use annual-\nmean fields to calculate the corresponding parameters for each period in ICM\nsimulations, including p , D̅ , T̅ , and h̅*. For the surface velocity–SST component, we\nn tc tc\nuse seasonally varying climatological wind stress fields from the 8CMM to force the\nICM’s dynamical ocean model and derive climatological mixed-layer currents (u̅, v̅, w̅ );\nclimatological SST fields (T̅ ) are directly taken from the 8CMM.\nText:\n- l53: Maybe briefly motivate why you choose the SSP5-8.5 scenario. It is noteworthy\nthe highest emission scenario within CMIP6. Also is the non-monotonic increase (i.e.\nincrease then decrease) in ENSO amplitude unique to this scenario?\nRevised as suggested as indicated in the revised method section.\n- l120: in its depth and intensity -> in its depth, intensity and sharpness. Because Fig 5\nshows sharpness.\n- l188: Figure 4 illustrates... This sentence should be moved to l177 or 178, no?\n- l316 \"This mechanistic explanation resolves a key paradox...\". This is a very nice two\nsentences that summarizes some important goal of the paper. You chould put these\nsentences in the introduction too.\n- l468: Extended Data -> Supplementary Information\n- Fig. 1: add vertical lines or colored areas to highlight periods P1, P2, P3\nAll are revised as suggested.\n- SI l110: Just to be clear, are you using exactly Zhang et al. 2020 ICM or some\ndifferences? The text suggests that the Te scheme is modified (Supplementary text 3\n25\nl333, e.g. the Te-SVD scheme has been replaced by Te-tanh scheme). Maybe add a\nsentence to summarize what is modified or not. Same for l24 and l128 of main text (the\nmodel is described as novel, but what is novel?).\nRevised as suggested. In particular, a paragraph is added in the revision as follows.\nIn the SST equation, subsurface entrainment temperature into the ML (T ) affects\ne\nSST through mean vertical upwelling and mixing, which is an essential component of\nSST variability. Previously, a statistical T model is developed based on empirical\ne\nrelationships between anomalies of T and thermocline depth using SVD analysis32,\ne\nwhich has been widely used for ENSO modulation studies48,49,57. However, the\nmodulating effects of the mean thermocline parameters (its mean depth and sharpness)\non T have not been explicitly represented in the SVD-based T parameterizations. It is\ne e\ndesirable to have explicit relationships of T with the thermocline depth anomaly, as\ne\nwell as its mean parameters (including depth and sharpness). In this work, , following\nthe work of Yuan et al. (2020), we take T as an analytical function linking it with\ne\nthermocline depth anomalies and mean thermocline properties in a hyperbolic tangent\n(tanh) form.\n- SI l157. Equation 5 -> 1.5\nRevised as suggested.\n- SI l278: 'h* is related to the vertical gradient of temperature field\"\". Please specify.\nSteepness is the amplitude of the dT/dz (or N^2) profile at the thermocline depth?\nRevised as suggested, making it more clearly in the Supplementary Text 3:\nAn analytical parameterization of T as follows:\ne\nobservations or ocean reanalysis data as follows.\nT 20C\nh* (D z)/arctanh( sub )\nr\n(3.2)\nBy inputting monthly temperature fields and 20°C isotherm depth data into Equation 3.2, the\nvalue of ℎ∗ can be derived. Given that interannual variability of ℎ∗ is relatively small, ℎ∗ can be\napproximated by its climatological mean.\n- SI l314: maximum temperature gradient -> maximum temperature vertical gradient\n- SI l 349: The notation ... (doesnt appear, I guess it is the brackets)\n- SI l366: BJ coefficients in SI Table 4 not 3.\n- SI Fig. 2: Specify the meridional average used to derive this figure. Is it 5N-5S?\nIn\npanel c it would be more consistent to show thermocline depth P2 (red) and P3 (blue).\n- SI Fig. 2 is the same information at Fig. 2 (although showing totals instead of\ndifferences). These figures could be merged (or at least shown one after the other). The\nthermocline depths should match but they dont between SI Fig. 2 and 4\n- SI Fig. 9, panel a: N->N^2\n- SI Table 2: Some missing parenthesis in the row \"Denotations\".\n- SI Table 4: You could add a column that indicates to which feedback (TH, EK, ZA...)\neach parameter is associated.\nAll revised as suggested and thank you for your details points!\nSupplementary Text 4: Methodological and scientific innovations in\nthis study\nThis study introduces a novel modeling framework that integrates three key\nphysical components into an ICM, enabling a systematic identification of the dominant\nfactors governing ENSO amplitude changes under sustained global warming. By\nisolating and quantifying the relative contributions of these components, we provide\nnew mechanistic insights into their combined effects on ENSO—particularly the\namplifying role of stratification and the damping roles of thermocline and surface field\nprocesses, respectively. To substantiate the advances made in this work, we outline the\ntechnical and scientific innovations below.\n4.1 Technical advances\n(a) Spatially explicit representation of stratification effects\nObservations reveal significant spatial and vertical heterogeneity in density and\nstratification across the tropical Pacific. While stratification is widely recognized as a\nkey modulator of ENSO under climate change, most prior studies have treated it\nqualitatively or diagnostically. In contrast, our approach incorporates spatially varying\ndensity fields to derive corresponding modal parameters (e.g., p₁ and p₂) for uses in the\nICM. This allows us to realistically capture how stratification changes under\nwarming—both horizontally and vertically—affect ENSO variability. For instance, the\namplitude of p₁ dominates over the western-central equatorial Pacific, influencing\nsurface responses, whereas p₂ is more prominent in the eastern equatorial Pacific,\naffecting subsurface processes. These spatially resolved modal structure and changes\nare critical for accurately depicting ENSO amplitude changes.\n(b) Coherent representation of the thermocline component\ne\nformulated using a hyperbolic tangent function that depends on mean thermocline depth,\nintensity, and sharpness50. This analytical representation ensures physical consistency\nwith the stratification component in the same ICM framework and allows T to\ne\ndynamically respond to changes in the mean ocean state.\n(c) Integrated framework with three interactive components\nAlthough stratification has been previously incorporated in certain modeling\ncontexts, this study is the first to coherently unify three key components—stratification,\nthermocline structure, and surface fields—within a single ICM framework. This\nintegration enables a systematic separation and evaluation of their individual and\nsynergistic effects on ENSO, through a set of carefully designed sensitivity experiments.\n(d) Application to a puzzling climatic projection period\nWe apply this framework to the post-2100 period, during which ENSO amplitude\ndeclines despite continued intensification of stratification, a paradox that current\nunderstanding fails to explain. Our model reconciles this counterintuitive behavior by\nrevealing how shifts in baroclinic mode structure alter the balance between amplifying\nand damping processes associated with the three influencing components of the mean\nocean state.\n4.2 Scientific insights\nUnlike previous largely qualitative or diagnostic assessments, our study quantifies\nthe effects of stratification and other components in a physically consistent and explicit\nmanner. The three mean-state components—upper-ocean stratification, surface fields,\nand thermocline structure—are derived from a 51-year running mean of eight CMIP6\nmodels under the SSP5-8.5 scenario (1990–2300). This setup supports several key\nfindings:\n(a) Non-monotonic ENSO responses linked to baroclinic mode shifts\nfuture climate projections, tracing it to changes in the amplitudes of the first two\nbaroclinic modes. These modal shifts alter the dominant feedback mechanisms\ngoverning ENSO growth.\n(b) A fruitful resolution of the post-2100 ENSO to stratification paradox\nWe provide a mechanistic resolution to the paradox of weakening ENSO amid\nstrengthening stratification. The answer lies in the vertical redistribution of\nstratification: surface-intensified warming weakens p₁ but enhances p₂, reducing the net\namplifying effect. When combined with persistent damping from thermocline and\nsurface field components, the result is a decline in ENSO amplitude.\n(c) Synergistic view of competing processes\nENSO amplitude emerges from a delicate balance between competing processes.\nThe damping effects from the thermocline and surface field components counteract the\namplify effects from stratification. During P3, the opposing changes in p₁ and p₂ further\ndiminish the net amplifying effect, leading to an overall damping of ENSO.\n4.3 Research positioning of this study relative to previous work\nWhile the use of ICMs and baroclinic mode theory is established (e.g., Zhang et\nal., 2003, 2020), our integration of a new T parameterization with stratification and\ne\nsurface components in a unified framework is also novel. This allows us to explicitly\nisolate and quantify competing effects, an apparent advance over previous studies such\nas Cai et al. (2018), which only identified a correlation between stratification and ENSO\namplitude, but could not explain its breakdown post-2100.\nIn this study, by coherently prescribing CMIP6-derived mean ocean states and\nsystematically conducting targeted sensitivity experiments, we not only replicate\ncomplex ENSO responses but also quantify the roles of individual components. This\n30\napproach offers a coherent and synergistic perspective on how mean ocean state\nchanges under prolonged warming shape the future of ENSO.\n31"}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have properly addressed my comments, and the revised manuscript is acceptable for publication.", "reviewer_2": "The authors have addressed all of my concerns and questions in a comprehensive and reasonable manner. The inclusion of the new subsection (5.1) and the supplementary discussion regarding the CMIP6 model sampling provides a clear and honest justification for their methodology. Furthermore, the additional mechanistic insights into the post-2100 ENSO weakening—specifically the analysis of the first two baroclinic modes (p_1 and p_2) and the competing effects of stratification and damping—significantly strengthen the physical interpretation of the results.The revised manuscript is much more robust, and the clarifications regarding the three-dimensional mean states and future research directions are well-taken. I am satisfied with the revisions and recommend the paper for publication in its current form.", "reviewer_3": "Second review:\nThe authors have addressed most of my previous comments and the revised version looks nice. My previous main comment was on some lack of novelty but given the editorial decision (which is fine by me) I wont comment further on that. There are two inconsistencies from my previous review, repeated below, that have not been fully addressed: I would strongly recommended checking them as they impact results that are central to the manuscript. Apart from this the manuscript is well revised and suitable for publication, with no further comments from me.\nComments:\n1) I have the same comment than in previous review: according to some figures, p1 continously increases from P1 to P3 which goes against the main point of the paper that p1 decreases over P3 leading to decreased ENSO amplitude. For reference your time periods are P1 (1940-1990), P2 (2040-2090) and P3 (2100-2149) (l107 or Fig 1a). Your baroclinic wind projection coefficients are p1, p2, p3 that depend on longitude x and time t.\n- According to Fig 4, p1(x,t) is p1(P3)<p1(P1)<p1(P2): p1 is strongest during P2 and weakest during P3. But according to SI Fig 4a, p1 in central Pacific is p1(P1)<p1(P2)<p1(P3): it continously increases from P1 to P3 then starts decreasing AFTER P3 i.e. after 2050. Thus the figures still mismatch. I also note both figures use same 2N-2S averaging area according to text (thus meridional averaging is not a reason for mismatch).\n- Could you also provide SI Fig 9 from previous version showing p1(x,t), p2(x,t) and also p3(x,t) to resolve the above? It is up to you but I would suggest keeping this figure in the revised manuscript or SI, it is crucial to your ICM model with zonally varying stratification.\n2) In my previous review, I had noted another inconsistency on SST standard deviation (std). In Fig1a, SST std in Nino 3.4 region shows std(P1)<std(P3)<std(P2). In Fig 2a, SST std along the equator shows std(P3)<std(P2)<std(P1), with std(P3) being dramatically lower. As you noted in your response, Fig1a uses 5N-5S average while Fig2a uses 2N-2S, which may explain the differences.\n- This however still raises the question of why the 2N-2S std decreases so dramatically (by a factor 2) while the 5N-5S std increases slightly comparing P1 to P3 (or P2 to P3). In my opinion there is still something important to check here, as 2N-2S vs 5N-5S SST averages usually dont show such drastic differences.\n- Could you briefly check the reason for this discrepancy? For example you may check how the meridional profile evolves by plotting SST std(y,t) as a function of latitude (range 10N-10S) and time (range 1900-2300), when averaged 170W - 120W on x as in Fig 1a. Maybe the std profile narrows or shifts off equator explaining the discrepancy. I could also comment that baroclinic modes have different Rossby radius (higher modes between more confined to the equator), thus modified baroclinic contributions could slightly alter the SST profile, but frankly I don't know if that is a sufficient explanation for such drastic 2N-2S vs 5N-5S discrepancies.\nText:\n- SI Fig 8-10: specify in captions that the winds are regressed on Nino 3.4 SST (as stated SI l640)"}, "author_response": "Thank you very much again for your comments, which are all valuable and extremely\nhelpful for revising and improving our paper. We have studied the comments very\ncarefully and have revised the paper accordingly. The major changes are in red font\nin the revised manuscript. The main revisions and the responses to specific\nsuggestions are as follows.\nSecond review:\nThe authors have addressed most of my previous comments and the revised version\nlooks nice. My previous main comment was on some lack of novelty but given the\neditorial decision (which is fine by me) I wont comment further on that. There are two\ninconsistencies from my previous review, repeated below, that have not been fully\naddressed: I would strongly recommended checking them as they impact results that\nare central to the manuscript. Apart from this, the manuscript is well revised and\nsuitable for publication, with no further comments from me.\nComments:\n1) I have the same comment than in previous review: according to some figures, p1\ncontinously increases from P1 to P3 which goes against the main point of the paper that\np1 decreases over P3 leading to decreased ENSO amplitude. For reference your time\nperiods are P1 (1940-1990), P2 (2040-2090) and P3 (2100-2149) (l107 or Fig 1a). Your\nbaroclinic wind projection coefficients are p1, p2, p3 that depend on longitude x and\ntime t.\n- According to Fig 4, p1(x,t) is p1(P3)<p1(P1)<p1(P2): p1 is strongest during P2 and\nweakest during P3. But according to SI Fig 4a, p1 in central Pacific is\np1(P1)<p1(P2)<p1(P3): it continously increases from P1 to P3 then starts decreasing\nAFTER P3 i.e. after 2050. Thus the figures still mismatch. I also note both figures use\nsame 2N-2S averaging area according to text (thus meridional averaging is not a reason\nfor mismatch).\n- Could you also provide SI Fig 9 from previous version showing p1(x,t), p2(x,t) and\nalso p3(x,t) to resolve the above? It is up to you but I would suggest keeping this figure\nin the revised manuscript or SI, it is crucial to your ICM model with zonally varying\nstratification.\nResponse: According to your comments, we have carefully check all the texts and\nexperimental designs, and we found a typo mistake that is made in the text writing: the\nthird period (P3) defined in the paper should be 2240-2290, but is accidently written as\n2100-2149. We apologize for this mistake that has caused the misleading and\ninconvenience.\nAs a result, all these issues would disappear; the related texts and figures have been\nupdated accordingly. For example, Fig. 1 looks like this (the third period is for 2240-\n2290).\nFig. 4 and Supplementary Figure 4 are all consistent, as indicated below.\nFigure 4 Zonal structure of wind projection coefficients (p_1,p_2, p_3) over the\nequatorial Pacific (averaged between 2°S and 2°N). a-c, The first, second and third\nbaroclinic modes are calculated separately from the 8CMM during three periods of\n1940-1990 (denoted as P1; red line), 2040-2090 (P2; green line) and 2240-2290 (P3;\nblack line), representing their historical and future distributions under the SSP585\nscenario.\nSupplementary Figure 4. Time evolutions during 1900-2300 for the 8CMM-derived\nocean fields averaged over the central equatorial Pacific (CEP; 160°E - 160°W, 2°S -\n2°N). (a) Stratification-related wind projection coefficients for the first three baroclinic\nmodes; (b) Thermocline and pycnocline depth (unit: m); (c) Thermocline sharpness at\ndepths of 60 m and 80 m (unit: m).\n2) In my previous review, I had noted another inconsistency on SST standard\ndeviation (std). In Fig1a, SST std in Nino 3.4 region shows std(P1)<std(P3)<std(P2).\nIn Fig 2a, SST std along the equator shows std(P3)<std(P2)<std(P1), with std(P3)\nbeing dramatically lower. As you noted in your response, Fig1a uses 5N-5S average\nwhile Fig2a uses 2N-2S, which may explain the differences.\n- This however still raises the question of why the 2N-2S std decreases so dramatically\n(by a factor 2) while the 5N-5S std increases slightly comparing P1 to P3 (or P2 to P3).\nIn my opinion there is still something important to check here, as 2N-2S vs 5N-5S SST\naverages usually dont show such drastic differences.\n- Could you briefly check the reason for this discrepancy? For example you may check\nhow the meridional profile evolves by plotting SST std(y,t) as a function of latitude\n(range 10N-10S) and time (range 1900-2300), when averaged 170W - 120W on x as in\nFig 1a. Maybe the std profile narrows or shifts off equator explaining the discrepancy.\nI could also comment that baroclinic modes have different Rossby radius (higher modes\nbetween more confined to the equator), thus modified baroclinic contributions could\nslightly alter the SST profile, but frankly I don't know if that is a sufficient explanation\nfor such drastic 2N-2S vs 5N-5S discrepancies.\nResponse: According to your comments, we have carefully check all the texts and\nexperimental designs, and we found a typo mistake that is made in the text writing: the\nthird period (P3) defined in the paper should be 2240-2290, but is accidently written as\n2100-2149. We apologize for this mistake that has caused the misleading and\ninconvenience.\nAs a result, all these issues would disappear; the related texts and figures have been\nupdated accordingly. Additionally, all these comments are still very helpful for our\nfuture analyses and modeling efforts.\nText: - SI Fig 8-10: specify in captions that the winds are regressed on Nino 3.4 SST\n(as stated SI l640)\nResponse: Accordingly, we have revised it. Thank you very much!"}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have addressed my last comments on figure mismatches, it looks good. I have no further comments and recommend for publication."}, "author_response": ""}]

Decision letter