Introduction

Predicting changes in the pattern of atmospheric circulation and associated atmospheric teleconnections over the tropics is critical for agriculture, water resources, and ecosystem management. The Walker circulation consists of the zonal atmospheric overturning cells over the tropics with surface westerlies (easterlies) over the equatorial Indian (Pacific) Ocean1. Variations in the structure of the Walker circulations are closely tied not only to precipitation and ocean temperature patterns2,3,4,5 but also to fishery stocks through variations in nutrient supply caused by the upwelling of cold water6. Recent studies of the Walker circulation over the Indo-Pacific region have conflictingly suggested either a weakening6,7,8,9,10,11,12,13,14,15,16 or strengthening1,3,4,17,18,19,20,21,22,23,24 pattern in response to global warming. Recent findings have also indicated changes in the Walker circulation using various paleo-proxy datasets and model simulations of the past climate25,26,27,28.

Model-simulated weakening of the Walker circulation in response to increasing greenhouse gas concentrations has been attributed to a differential increasing rate between precipitation (~3%/K) and atmospheric water vapor (~7%/K)9,10. However, the changes in the zonal overturning circulation in the tropics can also be explained by enhanced dry static stability over non-convective regions29, an influence of solar cycle maxima30, mean vertical advection of stratospheric change31, and trans-basin sea surface temperature (SST) differences12,21,32. In addition, several other factors also contribute to zonal circulation changes. For example, the meridional atmospheric circulation in the Pacific Ocean, which is controlled by the anomalously high-pressure centered at mid-latitude, acts to strengthen the ascending branch over the western Pacific Ocean and contributes to the strengthening of the Walker circulation20. However, no existing study addresses the contribution of meridional circulation to zonal circulation changes over the Indian Ocean (IO). Thus, it is essential to understand the underlying physical mechanism of how the local meridional circulation contributes to the zonal circulation changes over the IO in a warming world. Furthermore, it is important to find out if the SST warming pattern plays role in the weakening of zonal circulation or whether the weakening of zonal circulation causes the SST warming pattern.

The IO Walker circulation (IWC) is characterized by a rising motion over the eastern IO (EIO) and a sinking motion over the western IO (WIO)33. Presently, inter-annual variability of IWC is closely linked to IO Dipole (IOD) events6,33,34,35,36. A positive (negative) IOD, characterized by lower (higher) than normal SSTs in the tropical EIO and higher (lower) than normal SSTs in the WIO37,38, is often associated with a weakened33,39,40 (strengthened41) IWC, similar to how El-Niño (La-Niña) events weaken (strengthen) the zonal circulation in the Pacific Ocean42,43. Climate model simulations forced by a high greenhouse gas emission scenario (Representative Concentration Pathway 8.5 (RCP8.5)) indicate increased occurrences of extreme positive IOD (pIOD) events by almost a factor of three44,45, which is thought to imply a weakening of the IWC in a warmer world33,46.

However, the reliability of such a projection has been questioned due to common model biases in the IO47,48,49,50. In comparison to observation, climate models exhibit stronger easterly winds over the equatorial IO, a larger west-east SST gradient, and a shallower thermocline in the EIO47,48. Previously, an “emergent constraints” approach was used to investigate the relationship between present-day errors and future climate projections, with the conclusion that the model’s projection of pIOD-like warming pattern and the increasing frequency of extreme pIOD events are primarily artifacts of model errors51. A later study disputes the emergent constraints approach and demonstrates that other processes also contribute to the increased extreme pIOD events, which cannot be removed by the “emerging constraint” considering overestimating historical biases52.

The pIOD-like warming pattern has also been linked to a weakening of the Pacific Walker circulation53, but there is vigorous debate about whether the Pacific Walker Circulation projection is reliable. This suggests that the local component of the IO can play an important role in the projected IWC weakening in addition to the Pacific-induced changes. Therefore, we investigate the role of local component changes in the IWC weakening and the underlying physical mechanism. The present study is based on historical and high-emission coupled climate simulations from 24 models from the Coupled Model Intercomparison Project Phase 5 (CMIP5), 17 models from CMIP6, and a 40-member initial condition ensemble from the Community Earth System Model Large Ensemble (CESM-LE) climate simulations. We will treat the multimodel and initial condition ensembles somewhat differently in that it is only the latter for which the statistical significance of the forced response54 can be established, by evaluating statistics with respect to the ensemble-wise variability first55, because a multimodel ensemble does not represent an objective probability distribution56.

Results

IWC climatology and change in the historical period

To determine the climatological mean location of the ascending and descending branches of the IWC, we make use of the meridional average (10oN:10oS) of mid-tropospheric vertical velocity at 500 hPa (hereafter referred to as Omega500) in observation and climate model simulations during the period 1980–2015. We simply write ‘ensemble mean’ both with respect to members of a multimodel ensemble in CMIP5 and CMIP6, and an initial condition ensemble in CESM-LE.

Figure 1a shows the climatological mean location of the rising branch of the IWC over the EIO (94oE:104oE) and the sinking branch over the WIO (40oE:50oE). While the strength of the climatological descending branch over the WIO is consistent between observation and model simulations, climate models underestimate the strength of the ascending branch over the EIO. To determine the seasonal dependence of the inter-annual variability of the IWC during the period 1980–2015, we used the temporal standard deviation of the detrended seasonal mean IWC index (zonal Omega500 gradient), zonal SST gradient, zonal precipitation gradient, and surface zonal winds (SZW) over the central IO (CIO) (see “Methods” for more details). The observation and models in CMIP5 and CMIP6 show a pronounced seasonality with the peak occurring during the September–October–November (SON) season for the IWC index (Fig. 1b), zonal SST gradient (Fig. 1c), zonal precipitation gradient (Fig. 1d) and SZW over CIO (Fig. 1e). Thus, this study will focus on only the future change of the SON season. Note that, the seasonality is poorly represented in CESM-LE with the peak occurring during the December–January–February (DJF) season for the IWC index and SZW over CIO (Supplementary Fig. 1), which is caused in part by a bias in thermocline feedback57.

Fig. 1: IWC climatology based on the recent decades and its association with other variables.
figure 1

a Meridional average (10oN: 10oS) of mid-tropospheric vertical velocity (hereafter multiplied by the factor 100, Pa/s) at 500 hPa using ERA-5, CMIP5 (24 models), CMIP6 (17 models) and CESM-LE (40 members) simulations during the time period 1980–2015. The solid line represents the ensemble mean and the dashed line represents each individual model/member. Positive (negative) values represent the downward (upward) motion. The gray-shaded box shows the region of IWC. The vertical blue lines represent the IWC rising (94oE:104oE) and sinking (40oE:50oE) branches. be Temporal standard deviation of the seasonal mean b IWC index, c zonal SST gradient, d zonal precipitation gradient, and e CIO SZW in observation, CMIP5 and CMIP6 during the different seasons (DJF, MAM, JJA, SON) and annual mean for the time period 1980–2015, respectively. The stars (x) represent individual models/members and solid circles represent the ensemble mean temporal standard deviation (calculated upon detrending). fh Scatter plots represent the correlation of the observational IWC index anomalies with those of the f zonal SST gradient, g zonal precipitation gradient, and h CIO SZW. The years in green colors represent the extreme ±IOD year. The red dotted line represents the linear regression model. Correlation coefficients and corresponding P values are indicated.

In the observation and model simulations, the IWC index has a strong correlation with the zonal SST gradient, the zonal precipitation gradient, and SZW over CIO, suggesting a tight coupling between the oceanic and atmospheric anomalies12 (Fig. 1f–h and Supplementary Fig. 2). Regarding the ensemble mean trends over 1980–2015, model simulations exhibit a zonally oriented dipole pattern of Omega500 (more ascending in the west compared to the east), SST (more warming in the west compared to the east), and precipitation (more rainfall in the west compared to the east) overlaid by the anomalous surface easterlies over the equatorial IO (Supplementary Fig. 3a–c). A similar zonally oriented dipole pattern is also evident in the observations, indicating a weakening of the IWC over this time period (Supplementary Fig. 3d). However, our results or those from an earlier study58 cannot establish whether the dipole patterns of the forced response are related to the very similar dipole patterns of internal variability. Besides the spatial patterns, we also show in Supplementary Fig. 4 the changes in different IWC quantifiers are in loose association with Fig. 1f–h. In the case of the CESM-LE, only two out of the four IWC quantifiers indicate a statistically significant change.

Projected future weakening of the IWC

Previous studies reporting on projected future changes of the Walker Circulation under global warming have primarily focused on the Pacific Ocean. In order to understand the zonal atmospheric circulation changes over the IO in a warmer climate, first, we examine the SST changes. It shows an increase of the zonal SST gradient (more warming over the WIO compared to EIO) in the future (lower panels in Fig. 2), which is consistent with the increase already seen during the historical period (Supplementary Fig. 4) and is also rather consistent in all the simulations (Supplementary Fig. 6b, f, j). Next, we examine the ensemble mean changes in vertical velocity (upper panel Fig. 2a–c; upper panel in Supplementary Fig. 5), and precipitation (middle panel Fig. 2a–c; middle panel in Supplementary Fig. 5), which show anomalous descending motion and drier conditions over the EIO rim countries, but ascending motion and wetter conditions over the eastern African countries33,34,35,44,45 (Supplementary Fig. 6c, g, k). This east–west contrast in the projected precipitation change is consistent with the “warmer gets wetter” and “colder gets drier” hypothesis2,10,59.

Fig. 2: Projected IWC weakening due to pIOD-like warming pattern and its implication.
figure 2

a The upper plot presents the RCP85 minus historical difference of the meridional average (10oN: 10oS) of the ensemble mean vertical velocity (multiplied by the factor 100, Pa/s) at different pressure levels (1000–100 hPa). The middle plot in (a) concern the meridionally averaged (10oN: 10oS) value of precipitation (mm/day, gray) and surface zonal wind (meridional average (5oN: 5oS), m/s, green). The shading marks the estimate of ±1 ensemble standard deviation which represents both internal and intermodel variability. The lower plot in (a) show the SST change (shading) (RCP85 − Historical period). The black dashes line from top to bottom represents the IO region in each plot (40E:110E and 30S:30N). b Same as in a but for CMIP6 (SSP585 minus Historical). c As in (a) but for CESM-LE (RCP85 minus Historical). Stippling in the upper and lower plot in (a, b) indicate the areas where at least 70% of CMIP5 and CMIP6 models analyzed in this study exhibit the same sign as the ensemble mean changes. Stippling in the upper and lower plot in (c) indicate the regions where the difference between future and historical changes is significant at 95% confidence level based on a two-sample t test on the time series of ensemble means.

The historical versus future values of the IWC index clearly indicate greenhouse gas-induced IWC weakening (Supplementary Fig. 6a, e, i) by ~26%, 32%, and 31% in CMIP5, CMIP6 and CESM-LE, respectively. The strengthening of anomalous easterlies over the CIO is a robust (CMIP5/6)/statistical significant (CESM-LE) signal in the simulations with greater magnitudes observed in the CESM-LE than in CMIP6 and CMIP5 (middle panels of Fig. 2; lower panel in Supplementary Fig. 5 and Supplementary Fig. 6d, h, l). The anomalous changes in the strength of ascending and descending branches of the IWC are found to be coherent with the SST changes. This shows that the slowdown of IWC is due to a pIOD-like warming pattern. The anomalous descending motion over the EIO coincides very well with the precipitation changes (upper and middle panels of Fig. 2). The precipitation changes show a greater magnitude in CMIP6 and CESM-LE compared to CMIP5, as the weakening of IWC (more strong anomalous descending branch over the EIO) is more pronounced in the former than in the latter (upper panels; Fig. 2a–c). Overall, the model-projected future IWC slowdown is robust/statistically significant.

We note, that Supplementary Fig. 6 reveals an emergent constraint60, namely, that despite model biases with respect to the IWC strength, evidenced by the directional stretching of the cloud of data points, they all tend to feature a comparable slowdown. Related emergent constraints have been shown also in previous studies35. Moreover, the pIOD-like warming pattern is expected to lead to a weakening of both ascending and descending branches of the IWC. However, the model-projected changes in the ascending branch over the WIO are distinctly weaker than those for the descending branch over the EIO (upper panel in Fig. 2), suggesting that other factors also play an important role in modulating IWC changes. This shows that the causality cannot be easily determined between the weakened IWC and the surface warming pattern12.

Concerning the potential impact of the slowdown of the IWC, we analyze the changes in surface diatoms chlorophyll (SDC) (a proxy for upwelling61,62) over the IO. Previous studies have discussed that the stronger easterly winds over the CIO lead to the shoaling of the thermocline over the EIO in the future climate33,34,35,63. The projected changes in SDC show positive values over the EIO (the region of strong anomalous descending motion) in all the simulations, which can be associated with the shallower thermocline and more upwelling (Supplementary Fig. 7a–c). Hence, the IWC weakening results in an increase of SDC in the EIO in all the simulations implying the effect of IWC on the marine productivity in the projected climate.

Local meridional circulation changes contribute to the IWC weakening

We next analyze the projected changes in the mean sea level pressure (MSLP), surface (1000 hPa), and upper-level (200 hPa) winds in CMIP5, CMIP6, and CESM-LE simulations. The ensemble mean changes in MSLP indicate that anomalous low pressure is centered over the Arabian Peninsula (AP) region, accompanied by anomalous high-pressure centered over the EIO (Fig. 3a–c). Consistent with the strong pressure gradient between the EIO and the AP, anomalous low-level convergence in the AP and divergence in the EIO regions are pronounced in all the simulations (Fig. 3a–c). Furthermore, the circulation changes show upper-level divergence in the AP and convergence in the CIO and the EIO (Fig. 3d–f). We also noticed a slight southward shift of upper-level divergent wind pattern over the AP. This shifting can be due to the southward movement of upper-level winds over the AP in the warming climate compared to the present-day climatology (Figure not shown). These changes in atmospheric circulation suggest an anomalous local Hadley circulation (LHC) between the EIO and AP.

Fig. 3: Contribution of anomalous local meridional circulation to the IWC weakening.
figure 3

Ensemble mean changes (RCP8.5 minus historical) in a MSLP (hPa, shaded) and divergent winds (m/s, vectors) at the surface level (1000 hPa), d Tsurf (oC, shaded), and divergent winds (m/s, vectors) at the upper level (200 hPa) for CMIP5 (24 models). b, c Same as in (a) but for CMIP6 (17 models) and CESM-LE (40 members), respectively. e, f as in (d), but for CMIP6 (17 models) and CESM-LE (40 members). Stippling in (a, b, d, e) indicates the areas where at least 70% of CMIP5 and CMIP6 models analyzed in this study exhibit the same sign as the ensemble mean changes. Stippling in panels c and f indicate the regions where the difference in MSLP between future and historical changes is significant at 95% confidence level based on a two-sample t test on the time series of ensemble means.

To determine the cause of anomalous local meridional circulation with a center of rising motion over the AP, we examine the surface temperature (Tsurf) changes in the projected future climate. The ensemble mean changes in Tsurf mirror the MSLP changes with more significant anomalous warming over the AP and anomalous lower warming over the whole EIO (Fig. 3d–f). The patterns of MSLP and Tsurf changes are consistent in all the simulations. More warming over the AP than EIO enhances the land-sea thermal contrast in the future climate in all the simulations. Larger warming over the land is attributed to a larger lapse rate over the land rather than a lower heat capacity of the land64. Because the troposphere over land is often drier than that over the ocean, the lapse rate over land is higher, and a higher lapse rate results in more surface warming, especially in dry areas. These characteristics explain the substantial warming, low-level convergence, and low-pressure region over the AP and the lesser warming, high-level convergence and high-pressure region over the whole EIO.

As an alternative way to identify the existence of the anomalous LHC, we examine the changes in the velocity potential and the associated divergent winds. It shows strong convergence and divergence at the upper and lower level over the EIO, respectively, and the associated divergent wind pattern is towards the AP (EIO) at the lower (upper) level (Supplementary Fig. 8). The projected changes in the meridional MSLP gradient (see “Methods” for more details) are generally similar across the seasons, with a slightly larger gradient in the SON season in CMIP5 and CMIP6 as compared to CESM-LE (Supplementary Fig. 9). Therefore, in a warming world, the strong pressure gradient due to enhanced land-sea thermal contrast between the EIO and the AP during the SON season results in the slanting LHC which feeds the anomalous descending branch of the IWC and contributes to the weakening of the IWC.

An important question arises then: does the anomalous LHC result from the pIOD-like warming pattern? To give an answer, we further analyze atmosphere-only model simulations (AMIP5 and AMIP6) forced by prescribed SST changes31,65,66 (see “Methods” for more details). We observe that the patterned SST warming experiment results in a zonal dipole pattern of an anomalous low (high) pressure and upper (lower) level divergence (convergence) of winds over the WIO (EIO) (Fig. 4e, f and Supplementary Fig. 10c, d). Also, the longitude-height cross-section of vertical velocity change shows a dipole pattern of anomalous ascending and descending motion over the WIO and EIO, respectively, (Fig. 4g, h and Supplementary Fig. 11c, d) with anomalous easterlies (westerlies) overlaid at the lower (upper) level over the IO.

Fig. 4: Circulation changes in the uniform and patterned SST warming atmospheric-only model experiment.
figure 4

Ensemble mean changes in MSLP (hPa, shaded) and the divergent wind at 200 hPa (m/s, vectors) in a AMIP5 uniform SST warming (AMIP4K minus AMIP) and b AMIP6 uniform SST warming (AMIPP4K minus AMIP) simulations forced by prescribed SST change. e, f Same as in (a, b) but for AMIP5 patterned SST warming (AMIPFuture minus AMIP4K) and AMIP6 patterned SST warming (AMIPFuture4K minus AMIPP4K). The yellow boxes in figure a i.e., (90°E–115°E, 10°S–10°N) and (30°E–55°E, 15°N–35°N) are used to calculate the meridional MSLP gradient which represents the manifestation of the LHC on the surface (see “Methods“ for more details). The longitude-height cross-section of ensemble mean changes in vertical velocity (Pa/s, shaded) and zonal wind-vertical velocity (vectors) from 1000–150 hPa in (c) AMIP5 and (d) AMIP6 uniform warming. g, h Same as in (c, d) but for AMIP5 and AMIP6 patterned warming experiment. Vertical velocity is multiplied by a factor 700. i Bar plot represents the ensemble mean changes in MSLP gradient (hPa) in CESM-LE (40 members, RCP8.5 minus historical), CMIP6 (17 models, SSP585 minus historical), and CMIP5 (24 models, RCP8.5 minus historical) and their comparison with uniform and pattern SST warming in AMIP5 and AMIP6. The error bars represent the ensemble standard deviation.

In contrast, the uniform SST warming experiment produces an anomalous low (high) pressure and upper (lower) level divergence (convergence) of winds over the AP (EIO) (Fig. 4a, b and Supplementary Fig. 10a, b), analogous to what is seen in the CMIP5, CMIP6, and CESM-LE simulations (Fig. 3a–c). Moreover, in the uniform SST warming experiments, the north-westerlies (upper LHC branch) are observed from the AP to the EIO (Supplementary Fig. 11a, b) and the longitude-height cross-section of vertical velocity change shows an anomalous sinking motion over the CIO and the EIO (Fig. 4c, d and Supplementary Fig. 11a, b). This implies that the anomalous LHC in the ocean-atmosphere coupled model simulations is unlikely to result from a pIOD-like warming pattern. Therefore, the anomalous LHC acts independently of the pIOD-like warming pattern, actually contributing to the weakening of the ascending branch of the IWC over the EIO.

Taken together, a comparison of the coupled model simulations (CMIP5, CMIP6, and CESM-LE) with the uniform and patterned SST warming atmospheric-only model experiment (Fig. 4i) implies that the changes seen in the coupled model simulations are unlikely to be driven solely by patterned SST warming, i.e., pIOD-like warming pattern.

Discussion

The pIOD-like warming pattern is expected to weaken both the ascending and descending branches of the IWC. However, model simulations indicate comparatively less pronounced weakening in the descending (or anomalous ascending) branch in the WIO. This suggests that the pIOD-like warming pattern is not the only process leading to the model-projected weakening of the IWC in a future warmer climate. In a warming climate, the ensemble mean changes in CMIP5, CMIP6 and CESM-LE display a strong pressure gradient due to enhanced land-sea thermal contrast between the AP and the EIO as well as convergence/divergence at the lower/upper-level, respectively, suggesting the existence of the anomalous LHC. The anomalous LHC acts to weaken the ascending branch of the IWC over the EIO and contributes to the IWC weakening. By analyzing the two AMIP experiments, we show that the anomalous LHC is closely linked to the enhanced land-sea thermal contrast between the AP and the EIO, rather than to the pIOD-like warming pattern. Thus, we conclude that although the pIOD-like warming pattern is responsible for some of the projected weakening of the IWC, the anomalous LHC—independent of the pIOD-like warming—also plays an important role in the IWC weakening. The spatially uniform SST forcing in an atmosphere-only simulation does produce an anomalous LHC and ascending motion over the AP, which further suggests that the changes seen in the coupled climate models are unlikely to be driven solely by pIOD-like warming pattern. A schematic diagram representing the studied mechanisms is shown in Supplementary Fig. 12. This projected slowdown of the IWC has far-reaching implications for climate and weather extremes over the nearby continents and marine productivity as discussed earlier. A considerable influence of the IWC changes on remote climate necessitates future studies.

Methods

Coupled model simulations

This study used coupled model simulation outputs from the CMIP5 (24 models) archive for the historical and RCP8.5 scenario67. The CMIP5 climate models analyzed in this study are ACCESS1–0, ACCESS1–3, BCC-CSM1–1, CSIRO-Mk3–6–0, CanESM2, CMCC-CESM, GFDL-CM3, GFDL-ESM2M, GFDL-ESM2G, GISS-E2-H, GISS-E2-R, HadGEM2-ES, HadGEM2-AO, INMCM4, IPSL-CM5A-LR, IPSL-CM5A-MR, MIROC-ESM, MIROC-ESM-CHEM, MIROC5, MPI-ESM-LR, MPI-ESM-MR, MRI-CGCM3, MRI-ESM1, and NorESM1-M. In addition to CMIP5 simulations, coupled model simulations from the CMIP6 archive are used to further examine changes in the IWC for the historical and SSP585 scenario68. The CMIP6 climate model simulations used in this study are ACCESS-CM2, ACCESS-ESM1–5, BCC-CSM2-MR, CAMS-CSM1–0, CNRM-CM6–1, CNRM-ESM2–1, CanESM5, EC-Earth3, EC-Earth3-Veg, GFDL-ESM4, GISS-E2–1-G, HadGEM3-GC31-LL, IPSL-CM6A-LR, MIROC6, MRI-ESM2–0, NESM3, and UKESM1–0-LL. For each model in CMIP5 and CMIP6, only the first realization (i.e., r1i1p1 for CMIP5 and r1i1p1f1 for CMIP6) is used in this study. Supplementary Fig. 13 shows that the models are not sensitive to the number of realizations selected. Also, the statistically significant difference between the initial realization and the ensemble of realizations is negligible, in terms of the long-term mean climate. Furthermore, the number of available ensemble members is different among climate models. This means that if all available ensemble members are used to compute the ensemble mean change, more weighting is given to models that have a larger number of ensemble members. Because of this, for all models, the first ensemble was chosen.

Furthermore, we also analyzed the coupled model simulation output from the CESM-LE project, which consists of 40 members with slightly different initial conditions in the atmospheric state69. All the members of the CESM-LE have the same emission scenario as CMIP567, in which historical forcing is applied from 1920 to 2005 and RCP8.5 forcing from 2006 to 2100. The periods analyzed span 55 years each, 1950–2005 and 2035–2090 for the historical and RCP8.5 (SSP585 for CMIP6) simulations, respectively. We calculate the changes in the projected climate by taking the difference between the future (2035–2090) and the historical (1950–2005) periods in all the simulations. All CMIP5, CMIP6, and CEMS-LE coupled model simulations are re-grided to the same grid of 2.5 × 2.5 using bilinear interpolation.

Atmosphere-only model simulations

We analyzed atmosphere-only simulations from the CMIP567 and CMIP668 archives and name them AMIP5 and AMIP6, respectively. The atmosphere-only simulations in AMIP5 used in the study are (1) AMIP simulations which run from 1979 to 2008 forced with observed monthly mean SST and sea ice concentration; (2) the AMIP plus patterned anomaly simulations (AMIPFuture) which are the same as AMIP except for adding SST anomalies taken from the coupled model CMIP3 experiment; (3) the uniform SST increase simulations (AMIP4K), which are the same as AMIP except for adding a uniform +4 K SST “anomaly”. The effect of the pattern of SST change (pattern warming) is estimated by subtracting changes in AMIPFuture from AMIP4K, furthermore, the influence of uniform SST change (uniform warming) is calculated by subtracting AMIP4K from AMIP. The soundness of this methodology hinges on the fact that in a coupled model the SST change acts approximately as a forcing to the atmosphere. It is verified by Fig. 366, showing that a realized SST change by the coupled model taken as a forcing for the atmosphere-only simulation does result approximately in the atmospheric response of the coupled model. The AMIP5 models considered in this study are IPSL-CM5A-LR, MPI-ESM-LR, MPI-ESM-M, MRI-CGCM3, CanAM4, HadGEM2-A, and IPSL-CM5B-LR.

The atmosphere-only simulations in AMIP6 are (1) AMIP run with historical observed SST and external forcing from 1979 to 2014; (2) AMIPFuture4K simulations are the same as the AMIP experiment, but SSTs are subject to a composite of SST warming pattern derived from coupled models, scaled to an ice-free ocean mean of 4K. (3) AMIPP4K simulations involving AMIP experiments where SSTs are subject to uniform warming of 4K. We subtract AMIPP4K from AMIPFuture4K to get the response of IWC to the patterned SST change (pattern warming) and AMIPP4K from AMIP to get the response of IWC to the uniform SST change (uniform warming) in AMIP6 simulations. The AMIP6 models considered in this study are BCC-CSM2-MR, CanESM5, IPSL-CM6A-LR, MRI-ESM2–0, GFDL-CM4, and HadGEM3-GC31-LL.

In both AMIP5 and AMIP6 uniform SST warming experiments, the SST uniformly increases over the ocean, it does not affect the mean west-east gradient in the base climate (Supplementary Fig. 14). We choose only those models in the AMIP5 and AMIP6 that have runs for all the variables used in this study. Only the multimodel ensemble mean variations for the AMIP5 and AMIP6 simulations are shown. All AMIP5 and AMIP6 simulations are re-grided to the same grid of 2.5 × 2.5 using bilinear interpolation.

Model simulations for SDC

To determine the implication of the IWC on the ocean ecosystem, we used the SDC (a proxy for ocean upwelling region/shoaling thermocline zone) simulations from seven CMIP5 models67 (CMCC-CESM, CNRM-CM5, GISS-E2-H-CC, HadGEM2-CC, HadGEM2-ES, IPSL-CM5A-LR, and IPSL-CM5A-MR), four CMIP6 models68 (CESM2-WACCM, CESM2, CanESM5-CanOE, and GFDL-ESM4) and ten CESM-LE members69 (002, 009–017 members). The periods span 55 years each, 1950–2005 and 2035–2090 for historical and RCP8.5 (SSP585 for CMIP6) simulations, respectively. Only a limited number of models in CMIP5 and CMIP6 provide information related to surface diatoms chlorophyll. The difference between the two time periods is used to study the associated impact of the IWC on SDC in the projected climate.

Reanalysis and observational datasets

To define the criteria for the IWC rising and sinking branches, we used mid-tropospheric vertical velocity from the reanalysis dataset (ERA-5)70 from 1980 to 2015. SST observations are acquired from HadISST71. Precipitation observations are acquired from GPCP, version 2.372, which is a combined precipitation dataset (in situ observations over land merged with satellite-based precipitation) during the same period.

Estimating the IWC index and other variables

The meridional average (10oS–10oN) of the omega500 over the IO is used to define the IWC rising branch over the EIO (10°S–10°N, 94°E–104°E) and the IWC sinking branch over the WIO (10°S–10°N, 40°E–50°E). Positive (negative) values of omega500 represent the sinking (rising) motion. The IWC index is defined as the difference of area- average Omega500 between the EIO and the WIO, i.e., IWC index = (EIOomega500) – (WIOomega500). The SZW averaged over the central equatorial IO (5°S–5°N, 60°E–90°E) is used to determine the strength of the surface branch of IWC. The Zonal SST gradient is calculated as the difference of area-averaged SST anomaly in the WIO and the EIO (Zonal SST gradient = WIOSST – EIOSST)37. Also, the zonal precipitation gradient is calculated as the difference of area average precipitation between the EIO and the WIO, i.e., zonal precipitation gradient = (EIOprecip) – (WIOprecip). Furthermore, the MSLP/Tsurf gradient is calculated firstly as the difference of area average MSLP/Tsurf between the EIO (90°E–115°E, 10°S–10°N) and AP (30°E–55°E, 15°N–35°N), i.e., MSLP/Tsurf gradient = (EIOmslp) – (APmslp). As the gradient is not entirely meridional, we modified the index by removing the zonal component such as meridional MSLP/Tsurf gradient = MSLP/Tsurf gradient – (10S:25N area average of MSLP/Tsurf over EIO – 10S:25N area average of MSLP/Tsurf over AP). For all of the analyses in the manuscript, we used the LHC as the meridional MSLP gradient. The most prominent changes in the IWC occur during the IOD years, which peaks during the SON season (Fig. 1b–e). Therefore, all the calculations in this study are performed over the austral spring (SON) season.