Increasing ENSO–rainfall variability due to changes in future tropical temperature–rainfall relationship

Intensification of El Niño-Southern Oscillation (ENSO)-rainfall variability in response to global warming is a robust feature across Coupled Model Intercomparison Project (CMIP) iterations, regardless of a lack of robust projected changes in ENSO-sea-surface temperature (SST) variability. Previous studies attributed this intensification to an increase in mean SST and moisture convergence over the central-to-eastern Pacific, without explicitly considering underlying nonlinear SST–rainfall relationship changes. Here, by analyzing changes of the tropical SST–rainfall relationship of CMIP6 models, we present a mechanism linking the mean SST rise to amplifying ENSO–rainfall variability. We show that the slope of the SST–rainfall function over Niño3 region becomes steeper in a warmer climate, ~42.1% increase in 2050–2099 relative to 1950–1999, due to the increase in Clausius–Clapeyron-driven moisture sensitivity, ~16.2%, and dynamic contributions, ~25.9%. A theoretical reconstruction of ENSO–rainfall variability further supports this mechanism. Our results imply ENSO’s hydrological impacts increase nonlinearly in response to global warming. Tropical rainfall variability associated with El Nino-Southern Oscillation increases nonlinearly with rising sea-surface temperatures in a warming climate, through increased moisture sensitivity as well as dynamic contributions, suggest analyses of the CMIP6 climate model simulations.

T he El Niño-Southern Oscillation (ENSO) is a key driver of global climate variability on interannual timescales 1 , influencing global monsoons 2,3 and extreme weather [4][5][6] . Thus, it is important to understand whether ENSO and/or its impacts will change in response to increasing concentrations of greenhouse gases (GHGs). However, projected changes in the amplitude of eastern equatorial Pacific sea-surface temperature (SST) variability exhibit little inter-model consensus over past Coupled Model Intercomparison Project (CMIP) iterations [7][8][9][10][11] . This large inter-model discrepancy is attributable to various causes, such as a low signal-to-noise ratio given large internal variability [12][13][14][15] as well as pronounced and divergent model biases in reproducing both the tropical mean-state and ENSO characteristics [16][17][18] .
Although the uncertainty in future ENSO SST amplitude changes is large, ENSO's rainfall changes are likely to intensify in a warmer world (see Supplementary Fig. 1 and Note 1), presumably due to the increase in mean-state atmospheric moisture 9,19,20 . Previous studies reported that changes in tropical rainfall are closely linked to the enhanced equatorial warming pattern 21 seen in most model projections: more warming in the eastern equatorial Pacific than the surrounding region enhances moisture convergence toward areas of greater surface ocean warming and thus rainfall increase over the eastern equatorial Pacific, according to the "warmer-gets-wetter" mechanism 21,22 . We emphasize here that uncertainties in projections of the SST mean-state warming pattern are also a source of inter-model spread in future ENSO projections 23,24 .
The relationship between tropical SST (T) and rainfall (P) is characterized by two important aspects 25 : (i) an SST threshold (~27.5°C in the present), which marks the occurrence of tropical deep convection, and (ii) a rainfall sensitivity to SST change. When considering only the effect of a convective threshold increase in response to global warming (due to increased tropical atmospheric stability) 26 , we expect a shift of the tropical rainfall probability density function (PDF) toward higher SST values. However, there is a growing evidence on the changing tropical T-P relationship in response to mean SST rise, such as nonlinear rainfall sensitivities to El Niño SST anomalies 27 , rainfall PDF shifts toward a more intense rainfall 28 , and rectified changes of mean rainfall by ENSO variability 29 . Although previous studies 9,[27][28][29] showed the contribution of the T-P relationship changes on the projected increase in ENSO-related tropical rainfall variability, a question still remains as to whether there are other fundamental changes in the underlying statistical relationship between tropical mean SST and rainfall. Furthermore, increases in mean SST and moisture that were previously suggested as important drivers 9,19,20 cannot fully explain how changes of the mean SST translate to an increase in tropical Pacific rainfall variability.
The T-P relationship changes in response to greenhouse warming are difficult to address due to limitations in separating contributions from different components (e.g., mean state, interannual variability, and nonlinear coupled feedbacks). To address these key questions, we apply a PDF transformation approach in which the rainfall PDF can be explained by a combination of the SST PDF distribution and the functional form of the T-P relationship (see Supplementary Note 2 for detail). We hypothesize that the intensification of ENSO-rainfall variability is primarily a consequence of the underlying nonlinear tropical T-P relationship changes, not resulting from the mean SST distribution changes. We re-examine the historical and projected changes of ENSOdriven tropical SST and rainfall variability, focusing in particular on the latest generation of CMIP6 models 30 (Supplementary  Table 1).

Results
Future changes in ENSO SST and rainfall variability. We first show the changes in ENSO amplitude (standard deviation σ T and σ P ; see "Methods") in terms of eastern equatorial Pacific T and P under the four tier-1 Shared Socioeconomic Pathway (SSP) scenarios using 30−34 CMIP6 models (depending on the scenario) and under the Representative Concentration Pathway (RCP) 8.5 scenario using 40 CMIP5 models (Fig. 1). The four SSPs are SSP1-2.6 (33) for sustainable development, SSP2-4.5 (34), SSP3-7.0 (30), and SSP5-8.5 (34 models) for fossil-fueled development. ENSO SST and rainfall indices are defined using monthly SST and rainfall, respectively, averaged over the Niño3 region [5°S-5°N, 150°W-90°W] where a strong T-P coupling occurs during extreme El Niño events in association with zonal and meridional shifts of the South Pacific Convergence Zone and Inter-Tropical Convergence Zone 9,31 . The fractional changes in ENSO SST and rainfall amplitude (Δσ T and Δσ P in % unit) between the period of 2050-2099 and the present-day period of 1950-1999 are correlated across CMIP6 SSPs, as well as in CMIP5 RCP 8.5 simulations with correlation coefficients attaining values of about 0.8. The strong relationship between Niño3 SST and rainfall changes is also found in time-averaged mean quantities ( Supplementary Fig. 2) due to ENSO-related air-sea-coupled processes 1,14 .
In spite of the strong relationship, increases in ENSO-rainfall variability are projected even with decreases in ENSO SST variability (~38% among all cases). The ENSO SST amplitude changes are indistinguishable between different SSP and RCP scenarios, except for the subset of models that show a large strengthening of ENSO SST amplitude in the future under highemission scenarios. We also see that the ENSO-rainfall amplitude changes are almost double the SST amplitude changes across CMIP models (ratio of Δσ P to Δσ T ), in agreement with the previous study 28 showing a minor impact of future ENSO properties (e.g., amplitude and frequency) changes on rainfall variability. This ratio Fig. 1 SST-based versus rainfall-based changes in ENSO amplitude due to increasing GHGs. Scatter plot of standard deviation change for Niño3 SST (Δσ T ) and Niño3 precipitation (Δσ P ) between 2050-2099 and 1950-1999, obtained from multi-model simulations of CMIP5 and CMIP6. Each filled symbol and ellipse indicate the multi-model mean and 5-95% range of the data. The ratio of Δσ P to Δσ T for each scenario is displayed. is 2.0 for CMIP6 SSP5-8.5 and 2.82 for CMIP5 RCP 8.5. Combining CMIP6 results with CMIP5 simulations it becomes more evident that the projected change in ENSO-rainfall amplitude cannot be fully explained by the ENSO SST amplitude change, in accordance with previous studies 9,19,20,28 . This result implies nonlinear atmospheric responses to mean SST warming.
Future change in the tropical mean SST-rainfall relationship. Our hypothesis is that the changes in tropical SST-rainfall relationship lead to amplification of ENSO-rainfall variability even in the presence of an ENSO SST variability decrease. To further elucidate the underlying statistical relations between SST and rainfall, we show a scatter plot for the CMIP6 multi-model mean of monthly SST and rainfall over the entire tropical ocean [20°S-20°N], time-averaged for the present-day 1950-1999 (blue dots) and the future 2050-2099 (SSP5-8.5; red dots) (Fig. 2a). This figure illustrates the non-monotonic T-P relationship as seen in observations 25 : rainfall increases sharply above a SST threshold (T c ,~26.1°C in the case of P > 2 mm day −1 ; see "Methods") and then decreases once SST exceeds~29.6°C in the present climate. The nonlinear T-P relationship has been suggested to result from the complex interactions between SST, convection, and the largescale circulation 32,33 and involve both thermodynamic and dynamic processes. There is a noticeable shift of the T-P relationship toward higher SST values from the present-day to the future (Fig. 2b). The mean SST and rainfall over the Niño3 region (star symbols, Fig. 2a) increase by~3.1°C (11.7%) and 0.8 mm day −1 (34.8%) in the future, associated with the enhanced equatorial warming pattern 9,20 . The frequency histograms of simulated SST aggregated over tropical ocean and Niño3 space (f(T), Fig. 2b; sample number is 2448 grid points × 12 months for tropics and 125 grid points × 12 months for Niño3) also indicate a GHG-forced robust shift in the mean SST.
The previous study 26 discussed this shift of the T-P relationship toward higher SST values, arguing that the rainfall PDF will not change in a warmer climate. However, the frequency histograms of rainfall (g(P), Fig. 2c) show a clear shift in the upper tail of the rainfall distribution. This disparity could be attributable to an overlooked aspect of the annual mean T-P relationship compared to a monthly mean resolved relationship (see also Supplementary Fig. 3 for the T-P seasonality). This might also be associated with different T-P mechanisms between annual mean and monthly mean-sampled data 34 . Moreover, we find an obvious change of the Niño3 rainfall distribution: an increase in the upper tail and a decrease in the lower tail, resulting in an increase of the mean rainfall. These PDF changes are also apparent in most of individual CMIP models ( Supplementary  Fig. 4). It was reported that the increase in mean rainfall is associated with the shift of mean SST distribution as a direct consequence to increasing GHGs 28 . Another explanation for the increase in mean rainfall is rectified changes in ENSO properties (e.g., amplitude and frequency) 29 . We also emphasize that the change in the rainfall PDF can be explained by a combination of two different factors: (i) shift of the SST distribution f(T) and (ii) change of the T-P relationship (i.e., estimated PDF ; see Supplementary Note 2). To delineate these nonlinear T-P changes, we next formulate rainfall as a function of SST (i.e., P T ), which can be expressed as: A Gaussian fit was applied for calculating P T over convective regions (T > T c ), whereas a third-order polynomial fit was used over non-convective regions (T < T c ) (see "Methods" section for details). In Eq. (1), the Gaussian fit parameters (P 0 , T 0 , σ) indicate the maximum precipitation, the SST at maximum precipitation, , and the pseudo future reconstruction of 2050-2099 based on T <20th> + ΔT and P <20th> (sky-blue; ΔT = T 0<21st> − T 0<20th> ). The colored shading indicates the slope of saturated vapor pressure change to SST change (dE s /dT). The star symbols show the mean position of Niño3 region. P T , rainfall as a function of SST, was calculated by a piecewise function (3rd-order polynomial fit, in T < T c ; Gaussian function fit in T > T c ; see "Methods" for details). The black, dark green, and light green lines show the fitted P T for present, pseudo future, and future climate, respectively. Simulated frequency histograms for b SST (f(T)) and c rainfall (g(P)) and d fitted rainfall histogram as a function of T (g(P T )) over the tropics (shading) and Niño3 region only [5°S-5°N, 150°W-90°W] (solid line) are shown as inserts. The fitted rainfall histogram of pseudo future reconstruction over the tropics (skyblue dashed line in d) is also displayed. and the SST-rainfall spread, respectively. There is a spread in the T-P scatter for the CMIP6 multi-model mean (Fig. 2a), partly due to inter-model diversity in the T-P peak (see Supplementary  Fig. 5). Despite this uncertainty in function fits (color lines in Fig. 2a), we can see that the fitted rainfall distributions over both tropics and Niño3 region (g(P T ), Fig. 2d) are in reasonable agreement with the simulated distributions (g(P), Fig. 2c).
The rainfall variability as a function of SST (i.e., P T ) is determined from the temperature T and the Gaussian fit parameters (i.e., T-P relationship). To examine the respective effect of changes in f(T) versus those in the T-P relationship ( dT dP ), we calculate the pseudo future rainfall distribution that would result from applying the present-day T-P relationship to the simulated future SST distribution (sky-blue dots, Fig. 2). The PDF for pseudo future P can be mathematically expressed as g E ðP <pseudo> Þ ¼ f T <20th> þ ΔT ð Þ j dT dP j <20th> (with ΔT = 2.5°C), whereas the PDFs for simulated present-day and future P are g E P <20th> ð Þ¼ f T <20th> ð Þj dT dP j <20th> and g E P <21st> ð Þ¼f T <21st> ð Þj dT dP j <21st> . By comparing the frequency distributions of pseudo future, present-day, and future P T (see sky-blue dashed line, blue and red shading in Fig. 2d), we find that the future change of the rainfall distribution over the entire tropics as well as over the Niño3 region can be largely explained by the shape change (i.e., dT dP ), rather than the SST PDF shift toward higher SST (i.e., f(T)).
The comparison of the future T-P relationships between the direct model simulation (Fig. 2a, red dots) and pseudo reconstruction (Fig. 2a, sky-blue dots) highlights an increase in the slope of the T-P curve in convective regions (i.e., T > T c ). The steeper slope of the T-P curve (i.e., smaller dT dP ) leads to the changes in mean rainfall distribution over the Niño3 region, characterized by an increase in convective rainfall (P > 2 mm day −1 ) and a decrease in non-convective rainfall (P < 2 mm day −1 ). Although the convective SST threshold over the entire tropics does not change much from present-day to future climate 26 , the convective region in the eastern Pacific expands in a warmer climate, particularly during the boreal winter ENSO mature season ( Supplementary Fig. 3). The increase in the extent of eastern Pacific convective rainfall is likely related to the eastward expansion of warm pool convection and an equatorward shift of the Inter-Tropical Convergence Zone, resulting from the enhanced equatorial warming pattern and a projected eastward shift of the Walker circulation 9,35 .
Linking the tropical SST-rainfall relationship to ENSO-rainfall variability. Considering the critical role of T-P slope on the rainfall PDF, the changes in mean-state T-P slope ( ∂P ∂T j T ; see "Methods") can reproduce the model-simulated changes in rainfall variability (σ P ). Based on the CMIP6 multi-model mean, the percentile changes in the T-P slope agree with overall changes in interannual rainfall variability (Fig. 3a). The T-P slope increases over the Niño3 region are statistically significant for the CMIP6 multi-model mean (p < 0.005) and across~75% of CMIP models (p < 0.025) (Supplementary Fig. 6). The change in mean-state T-P slope could be explained by both thermodynamic and dynamic processes 36,37 . The previous study 37 showed that the precipitation change can be reasonably decomposed into two dominant components: (i) the CC-related thermodynamic change and (ii) dynamic change associated with divergence feedback and spatial changes in circulation. In a similar way, we next decompose the T-P slope change into these two components. In a thermodynamic point of view, we emphasize an important increase in slope of saturated vapor pressure change to rising SST (i.e., larger dE s /dT; E s ¼ 6:1094e ½ 17:625T Tþ243:04 in the August-Roche-Magnus approximation; called the CC function), as shown in gradient shading in Fig. 2a. This implies a nonlinearity in future changes of the CC relationship. The CC acceleration to mean SST warming, which is related to steeper mean-state T-E s slope ( ∂E s ∂T j T ; see "Methods") rather than increasing mean E s , can intensify the mean-state T-P slope. We can see that the thermodynamic changes in T-P slope (Δ ∂P ∂T * T ; given as a function of ratio of 21st mean-state T-E s slope to 20th slope) are consistent with the enhanced equatorial warming pattern (see Fig. 3b and Supplementary Fig. 7).
Accompanied by the increasing thermodynamic T-P slope, the largest dynamic changes in T-P slope (as residual between total T-P slope changes and thermodynamic changes) also occur in the eastern tropical Pacific (Fig. 3c), indicating the moisture convergence in the eastern Pacific and overall divergence in the Indian and western Pacific. Based on 55 CMIP models (25 CMIP6 and 30 CMIP5 models) that have more reliable T-P fit estimation, the thermodynamic changes in T-P slope explaiñ 30% of inter-model spread in total T-P slope changes ( Supplementary Fig. 8). The thermodynamic changes in T-P slope increase, on average, by~16.2 ± 3.6% with a strong intermodel agreement, whereas the total T-P slope increases bỹ 42.1% ± 31.3% with a large inter-model spread. This implies that the dynamic T-P slope increases by~25.9% with a low confidence, which is likely to be dependent on the model physics and mean-state model bias. The stronger contribution of dynamic slope changes than thermodynamic ones is also consistent with the previous study demonstrating the importance of enhanced SST warming amplitude over the eastern equatorial Pacific on the nonlinear precipitation response 27 .
Previous studies documented the thermodynamic processes of increasing ENSO-rainfall variability, which was associated with an enhanced equatorial mean-state warming pattern and increase in tropical rainfall variability over the central-eastern Pacific 9,20,35,38 . We further elucidate the theoretical grounds linking the mean changes to increasing rainfall interannual variability. In light of this, we propose a simple rainfall variability equation as follows: where the overbar indicates the time average of SST over the present-day (20th; 1950-1999) and future climate (21st; 2050-2099) and the prime indicates the time fluctuation of SST anomalies (see also Supplementary Fig. 9 for illustration). This rainfall equation does not consider the rainfall variability due to atmospheric intrinsic noise which can explain a substantial portion of the total rainfall variability in some regions 39,40 .
Although the neglection of the atmospheric intrinsic variability can induce some discrepancies in σ P between simulation and reconstruction, we exclude the effect of atmospheric noise on σ P due to the small contribution of atmospheric intrinsic variability in the eastern Pacific 40 . Based on Eq.
(2), we re-calculate the rainfall standard deviations at each grid point and in individual 55 CMIP models, respectively, in the present-day and future climate under the highest emission scenario (SSP5-8.5 in CMIP6 and RCP 8.5 in CMIP5). The comparison between the reconstructed percentile changes in σ P and simulated changes reaffirms the fidelity in reconstructing the tropical rainfall variability (Fig. 4a; also see Supplementary Fig. 10 for individual models). We emphasize that the rainfall variability can be divided into two components: a contribution from the time-averaged T (PðTÞ with change in only space) and a contribution from the time-varying T ( ∂P ∂T j T T 0 with change in time and space). From this separation, the mean rainfall change due to the mean SST increase (PðTÞ) does not affect the rainfall standard deviation change. We next formulate the ratio of ENSO-rainfall variability change (i.e., σ P<21st> σ P<20th> ; R P ) as the combination of the ratio of mean-state T-P slope changes (S R ) and ratio of ENSO SST variability (R T ). These overall increases in T-P slope and ENSO-rainfall variability are coherent across the individual CMIP6 and even CMIP5 models (r~0.59 with a p value of 0.007; Fig. 4b). The reconstructed rainfall variability changes (R P ) as a function of only SST variability changes (R T ) cannot reproduce the simulated increase in R P (blue dots in Fig. 4c). We emphasize that the reconstructed rainfall variability changes taking into account both SST variability changes and T-P slope changes (green dots in Fig. 4c) can explain~58% of intermodel variance in the simulated R P and further the projected increases in ENSO-rainfall variability.

Discussion
Due to the complex coupled feedbacks and nonlinear internal dynamics intrinsic to ENSO 1 , difficulties remain in constraining future changes in ENSO SST amplitude. There is also a large internal variability in ENSO SST amplitude change and divergent forced responses across models ( Supplementary Fig. 11). In contrast, ENSO-rainfall variability reveals a coherent forced response, showing a clear intensification in both CMIP5 and CMIP6 (Supplementary Fig. 1). Considering the disparity between SST and rainfall variability changes, mechanisms to account for the future changes in rainfall patterns, such as "warmer-gets-wetter" 21 and "wet-gets-wetter" 36 , cannot fully explain the pattern of changes in local tropical precipitation and its interannual variability. We propose the nonlinear thermodynamic rectification mechanism linking the mean SST warming to increasing rainfall variability, associated with the nonlinear change of CC relationship in response to global warming. The CC acceleration to mean SST warming (Δ ∂E s ∂T j T > 0), which is most pronounced in the region of the ) versus simulated rainfall variability changes (R P ; ratio of σ P<21st> to σ P<20th> ), c simulated rainfall variability changes versus reconstructed changes (R P ). Here, the R p is respectively reconstructed by different combinations of T-P slope changes (S R ) and SST variability changes (R T ; ratio of σ T<21st> to σ T<20th> ) based on Eq. (2). Green dots show the reconstructed rainfall variability changes by the combination of both S R and R T . Red and blue dots indicate the respective effects of S R and R T on the reconstructed rainfall variability changes. enhanced equatorial warming pattern, can intensify the tropical ENSO-rainfall sensitivity to mean SST increase (Δ ∂P ∂T j T > 0). ENSO-rainfall variability is coupled to the tropical T-P relationship. The models' fidelity in simulating the present-day T-P mean-state, particularly for tropical Pacific zonal T and P gradients ( Supplementary Fig. 12), might affect the projected future changes in ENSO property (e.g., amplitude and frequency) 7,17,24 . Although the projected changes in SST-based ENSO amplitude are unlikely related to the tropical T-P slope changes (see Supplementary  Fig. 8), the rectification of the mean T-P state by stronger ENSO SST variability 29 (r between S R and R T~0 .46 with p value of 0.01) may play a role in controlling the tropical rainfall variability. Understanding the model biases of the simulated present-day rainfall mean-state may be important in linking the mean-state T-P relationship to the projected ENSO-rainfall variability changes. We further highlight a novel aspect of the nonlinear changes of the CC relationship to rising mean SST and its thermodynamic and dynamic effects on the tropical T-P relationship. This newly proposed mechanism can provide additional insight into the fundamental properties of tropical T-P relationship in a changing climate, which is crucial for ENSO research and prediction. Our results further support the notion that ENSO-rainfall-driven teleconnections and extreme events, associated with more frequent occurrence of extreme El Niño events 9 , could further intensify in a warmer climate, irrespective of potential moderate changes in ENSO SST amplitude.

Methods
CMIP5 and CMIP6 models. This study analyzes multi-model simulations from CMIP5 41 Table 1). In CMIP6, a new range of SSP scenarios representing different future socioeconomic developments was developed 42,43 . The SSPs represent a range of challenges for mitigating and adapting to climate change. For example, the SSP5-8.5 scenario has an identical radiative forcing level to RCP 8.5 (i.e., 8.5 W m −2 at 2100) but the scenario assumes accelerated globalization and rapid economic and social development of developing countries coupled with the exploitation of abundant fossil fuel resources. The SSP narratives lead to higher levels of warming from SSP1-2.6, SSP2-4.5, SSP3-7.0, to SSP5-8.5. We primarily compared the second half of twentieth century (1950-1999 of historical run) with the second half of twenty-first century (2050-2099 of SSP5-8.5 run for CMIP6 and RCP 8.5 for CMIP5).
ENSO amplitude and its change due to GHG forcing. ENSO amplitude was measured as the standard deviation of Niño3 [averaged over 5°S-5°N, 150°W-90°W] SST and rainfall anomalies 1,9,17 . Monthly anomalies were calculated by removing the long-term mean seasonal cycle and the linear trend in each period. The linear trend was estimated by an ordinary least-squares fit of monthly SST and rainfall anomalies. Detrending isolates the interannual variability by removing the externally forced low-frequency variability. To calculate the change in ENSO amplitude, we calculated the fractional change between twenty-first century (2050-2099) Niño3 SST/rainfall standard deviation and twentieth century  standard deviation in each model.
Tropical SST-rainfall relationship. Following the method in ref. 26 , we determined the SST threshold (T c ) by finding the minimum SST for which the mean rainfall exceeds 2 mm day −1 , based on the entire tropical Ocean grids [20°S-20°N]. We attained similar results with different rainfall criteria (e.g., lower 25% of P), consistent with the previous study 26 . Here, the rainfall rate was calculated using the mean rainfall rate corresponding to SST binned at intervals of 0.1°C. For the rainfall rate curves in Fig. 2, we similarly calculated the mean rainfall rate for SST bins.
While rainfall increases almost linearly below T c , rainfall and convection above T c increases sharply 25 . Based on this process, P T was formulated by a piecewise function (Eq. (1)): The Gaussian fit was applied for calculating P T , if T is greater than the SST threshold (T c ; convective regions); a third-order polynomial fit was used for T < T c (non-convective regions) (see also Supplementary Fig. 13a). Although we quantify the rainfall variability as a function of SST, the rainfall variability itself can also affect the SST change. It was reported that the effect of precipitation on SST is particularly pronounced in warm pool regions 39,44 . Because of the two-way T-P interaction, the use of P T can underestimate the rainfall variability over the warm pool region. Moreover, in vicinity of high SSTs (i.e., warm pool region), a Gaussian fit function has a limitation owing to the non-monotonous increase of T-P relationship. This calls for a caution interpreting the warm pool T-P relationship, which will be carefully examined in future work.
The frequency histogram for simulated P (i.e., g(P)) and histogram for P T calculated by Eq. (1) (i.e., g(P T )) are compared in Supplementary Fig. 13b. To assess the fidelity of the Gaussian fit to the data, the chi square statistic (χ 2 ) was calculated (χ 2 ¼ P n P i ÀP Ti ð Þ 2 i =σ 2 i , where σ 2 is variance of data with n number). Although the future T-P scatter shows larger spread than the present-day scatter (χ 2 <20th> = 196.55 and χ 2 <21st> = 276.56), the overall results in both present and future climate are statistically significant with a p value < 0.05.
Mean-state T-P slope and T-E s slope. The mean-state T-P slope ( ∂P ∂T T ) over entire tropical ocean grids was calculated using P T formula: e.g., ∂P ∂T j T ¼ T 0 ÀT ð Þ σ 2 P T T > T c À Á . The mean-state T-E s slope was also calculated in the same way: ∂E s ∂T j T ¼ ð17:625 243:04ÞE s ðTþ243:04Þ 2 . Thermodynamic changes in T-P slope (Δ ∂P ∂T j * T ) between the present-day and future climate were defined as Δ ∂P ∂T j * T ¼ ∂P ∂T j T<20th> S T À ∂P ∂T j T<20th> . Where S T is the ratio of 21st mean-state T-E s slope to 20th slope. The dynamic changes in T-P slope were simply measured by residual anomalies between total T-P slope change and thermodynamic change (i.e., Δ ∂P ∂T T À Δ ∂P ∂T * T ).

Data availability
The CMIP data can be obtained from https://esgf-node.llnl.gov.