Decreased Indian Ocean Dipole variability under prolonged greenhouse warming

Domain: earth_science
Published: 2024-04-01
Editorial outcome: accept

Abstract
The Indian Ocean Dipole (IOD) is a major climate variability mode that substantially influences weather extremes and climate patterns worldwide. However, the response of IOD variability to anthropogenic global warming remains highly uncertain. The latest IPCC Sixth Assessment Report concluded that human influences on IOD variability are not robustly detected in observations and twenty-first century climate-model projections. Here, using millennial-length climate simulations, we disentangle forced response and internal variability in IOD change and show that greenhouse warming robustly suppresses IOD variability. On a century timescale, internal variability overwhelms the forced change in IOD, leading to a widespread response in IOD variability. This masking effect is mainly caused by a remote influence of the El Niño–Southern Oscillation. However, on a millennial timescale, nearly all climate models show a long-term weakening trend in IOD variability by greenhouse warming. Our results provide compelling evidence for a human influence on the IOD.

Introduction
The Indian Ocean Dipole (IOD) is a leading climate variability mode in the tropical Indian Ocean1–3. A positive IOD phase features anomalously warm sea surface temperature (SST) in the western Indian Ocean and cold SST in the east, and the zonal warm-cold pattern is reversed during its negative phase. IOD events alter ocean and atmospheric circulation patterns3–5 and the associated precipitation anomalies cause substantial impact on Indian Ocean-rim countries, such as droughts and floods in Australia6,7, Southeast Asia8–10, the Indian continent11,12, and East Africa13, wildfires in southeast Australia14 and Indonesia9, coral reef dieback in western Sumatra15, and malaria outbreaks in East Africa16.

Understanding the response of IOD to global warming has been one of the long-standing important problems in climate science4,5,17–20. The latest IPCC Sixth Assessment Report (AR6)21 has concluded that the human influence on IOD is not robustly detected during the observational period21,22 and the forced change in IOD by global warming in the future is uncertain due to a lack of robust evidence21,23. IOD variability, a fluctuation in climate variables induced by the IOD, is typically measured by the dipole mode index (DMI), which is defined as the difference between the western (50°E–70°E and 10°S–10°N) and eastern (90°E–110°E and 10°S–0°N) SST anomalies of the tropical Indian Ocean. This traditional IOD index serves as a fundamental physical metric for IOD dynamics and associated teleconnection impacts. So far, no significant trends in the DMI variability have been detected in either observations or future climate-model projections17,18, despite a strong warming signal being present in the tropical Indian Ocean18,24. Although some studies founded changes in the intensity and frequency of strong positive IOD events in the projected future climate using the principal component decomposition technique on SST5 and rainfall4,25, the response of DMI variability, which describes the warm-cold cycles of the tropical Indian Ocean SST anomalies, remains uncertain and no inter-model consensus has yet been found.

The lack of robust changes in DMI variability appears contradictory to our understanding of IOD dynamics. The mean state of the tropical Indian Ocean, which is known to physically shape the IOD feedback processes and variability, shows robust changes due to global warming: faster warming in the west than east (i.e., positive IOD-like warming), shoaling of the eastern thermocline, and easterly wind anomalies18. Such a clear warming pattern in the tropical Indian Ocean raises the question as to why the change in IOD amplitude is not as robust as the change in the background mean state.

One of the possible reasons for this problem is internal variability, a natural variation of IOD characteristics caused by intrinsic atmospheric and oceanic processes. Internal variability can mask forced changes in IOD variability due to global warming, resulting in a wide range of IOD amplitude projections and an obscured link with the changes in the mean background state of the tropical Indian Ocean. A large ensemble simulation using the Community Earth System Model shows that internal variability alone can generate widespread long-term trends in future IOD variability change; projected IOD amplitude in the twenty-first century either increases or decreases depending on the ensemble member with different internal variability, despite the same warming forcing26. Notably, a spectral analysis of the CMIP6 preindustrial control (piControl) experiment shows a surprisingly low-frequency peak in the inter-centennial IOD amplitude change, ranging from 1/400 year−1 to 1/50 year−1 (Supplementary Fig. 1). This implies that the internal variability has the potential to influence future CMIP6 IOD projections, which are typically assessed over an 86-year time frame (2015-2100). Taken together, this underscores the need for an unprecedented long length of warm climate simulations to reliably resolve forced response and internal variability.

In this Article, we utilize climate simulations from the Long Run Model Intercomparison Project27 (LongRunMIP), an archive of millennial-length scale fully coupled climate-model simulations, to clarify a forced response of IOD to greenhouse warming. The LongRunMIP provides global climate simulations subject to a wide range of atmospheric CO2 levels with a length of typically longer than 1000 years. The typical forcing scenario is an instantaneous doubling (abrupt2x), quadrupling (abrupt4x), and octupling (abrupt8x), and a gradual increase at a rate of 1% per year until the atmospheric CO2 level doubles (1pct2x) and quadruples (1pct4x). After the CO2 increase, the forcing is stabilized for typically longer than 1000 years. The LongRunMIP also provides a preindustrial control simulation with constant external forcings (control). Therefore, the LongRunMIP provides a unique opportunity to disentangle forced response and internal variability and to systematically study the response of IOD to different levels of greenhouse warming (e.g., ref. 28).

We evaluate the IOD simulation performance of all available models from the LongRunMIP archive, and select 9 models that can reasonably simulate the observed temporal and spatial characteristics of IOD (“Methods” section and Supplementary Discussion 1). Our analysis sample includes 9 control and 18 high-CO2 simulations (total 27) from 9 models (Supplementary Table 1). We use the DMI, a commonly used index for IOD, to describe IOD variability (“Methods” section). We apply two different measures for the strength of the IOD variability. The main metric is the IOD amplitude, defined as the standard deviation of the DMI. The other is IOD event intensity, defined as the peak DMI during the positive or negative IOD event (“Methods” section). Both metrics are highly correlated by definition (i.e., if the IOD amplitude is high, the IOD event intensity is very likely to be high, and vice versa), but we mainly use the IOD amplitude to examine the change in the IOD variability strength. The IOD amplitude is a fundamental metric that can show the overall response of the IOD variability to greenhouse gas warming. The simplicity of the IOD amplitude definition has advantages that facilitate the physical interpretation of the results.

Results
Equilibrium IOD response
We first analyze the equilibrium response of IOD variability to CO2 forcing. We use the last 500 years of variables where the global mean surface temperature (GMST) reaches nearly equilibrium (Supplementary Fig. 2) and define this time frame as the equilibrium period. The simulation array of different levels of equilibrium GMST allows us to examine the forced equilibrium response of IOD to greenhouse warming. We calculate the IOD amplitude for the equilibrium period of each control and high-CO2 simulation.

In almost all models, the IOD amplitude decreases monotonically with increasing CO2 (Fig. 1a). The exceptions are HadCM3L and MIROC3.2. HadCM3L simulates a decreased IOD amplitude in the abrupt4x, but an increased IOD amplitude in the abrupt2x and abrupt8x simulations. MIROC3.2 shows a slight increase in 1pct2x, but a large decrease in 1pct4x. Simulation-wise, the majority of the high-CO2 simulations (15 out of 18) show decreased IOD amplitude compared to their control simulations. The average intensity of positive and negative IOD events decreases in nearly all high-CO2 simulations (18 of 18 high-CO2 simulations for positive IOD and 15 of 18 for negative IOD), and the change is monotonic with CO2 in most models (Supplementary Fig. 3 and “Methods” section). Both the reduced SST variability in the eastern and western tropical Indian Ocean contribute almost equally to the decrease in IOD amplitude (Supplementary Fig. 4). The seasonal variability of IOD decreases in almost all seasons (Supplementary Fig. 5). These changes show that the weakening of IOD variability is robust under equilibrium greenhouse warming. The equilibrium sensitivity analysis (“Methods” section) suggests that the change in IOD amplitude per 1°C increase in GMST is estimated to range from –6.6% to 0.6%, and 8 out of 9 models show a negative sensitivity (Fig. 1b).

Notably, the CESM1.0.4 simulates a substantial weakening of IOD variability, that is a near-collapse response to CO2 forcing. The power spectra of the IOD variability in a broad frequency range (1/50 to 1/2 year−1) largely decreases with increasing CO2 levels (Supplementary Fig. 6). In particular, the spectral power of the IOD variability is almost indistinguishable from the climate noise level in the abrupt8x simulation. This suggests that the IOD mode is in a near-collapse state and that the tropical Indian Ocean SST variability is dominated by climate noise.

Such transition occurs not only in the amplitude but also in the seasonal phase locking and spatial pattern of the CESM1.0.4 simulations. The IOD commonly peaks in September in all CESM1.0.4 simulations, but a secondary peak appears in April with increasing CO2 levels (Supplementary Fig. 5), indicating an alteration of the typical seasonality of the IOD. This alteration in the seasonality also occurs in other models including HadCM3L and CCSM3II. HadCM3L simulations show a transition in seasonal peak from October (control and abrupt2x) to July (abrupt4x and abrupt8x). CCSM3II simulations show a gradual shift of seasonal peak from September (control) to August (abrupt700ppm) to June (abrupt1400ppm). Such a pronounced increase in IOD variability in boreal summer is consistent with the previous study19 shows an increased occurrence of early positive IOD events by global warming under the IPCC’s high-emissions scenario (although the warming scenario is different).

The spatial pattern of variability in CESM1.0.4 simulations also change from a dipole-like to a monopole-like pattern due to the weakening of the western Indian Ocean SST variability (Supplementary Fig. 7). The abrupt8x shows a clear monopole-like variability pattern, which is more appropriately referred to as Indian Ocean Monopole rather than an IOD. These results show that the main characteristics of the IOD in a significantly warm climate can substantially deviate from those of what are currently known about the IOD. This is somewhat analogous to the bifurcation or critical transition phenomenon29,30, which refer to qualitative changes in the characteristics of a system when it crosses the critical threshold.

Transient IOD response
Next, we analyze the transient response of the IOD amplitude to the CO2 forcing. We calculate the 100-year moving IOD amplitude for the entire simulation period. This window length corresponds to the typical time length for assessing future IOD change (e.g., comparing IOD variability for 1900–1999 and 2000–2099). The inter-simulation mean of the equilibrium GMST warming level of the quadrupling CO2 simulations (abrupt4x and 1pct4x) is 5.5°C (Supplementary Fig. 2). This is comparable to the projected warming level over 2081–2100 in the IPCC’s high-CO2 emissions scenario21, which ranges from 3.3°C to 5.7°C, although the detailed warming pathway is different.

The evolution of IOD amplitude against GMST level is shown in Fig. 2a. In nearly all high-CO2 simulations (17 out of 18), the IOD amplitude shows a long-term decreasing trend with increasing GMST. The exception is HadCM3L abrupt8x. This is consistent with the equilibrium response shown in Fig. 1, indicating that the long-term trend can be read as a forced response due to greenhouse warming towards the equilibrium state. However, within this long-term trend, the IOD amplitude largely fluctuates, not in line with the steadily increasing GMST. For example, the CESM1.0.4 abrupt4x simulation shows that the IOD amplitude increases in the GMST range between 18°C and 19°C, contrary to the long-term decreasing trend (Fig. 2a). This fluctuation is due to internal variability of the IOD amplitude, a deviation from the forced long-term trend. Throughout the paper, we refer to the anomalous short-term fluctuation in the variable that deviates from the long-term forced change as internal variability. This shows that the internal variability can mask the forced response to greenhouse warming by broadening the response range of the IOD amplitude.

To quantify the internal variability of IOD amplitude deviating from its forced long-term trend, we perform the transient sensitivity analysis (“Methods” section). Based on the results shown in Fig. 2a, we calculate the moving transient sensitivity of IOD amplitude against GMST and derive its distribution. This distribution represents an ensemble of all possible measurable transient sensitivities of IOD amplitude under a given warming forcing. If the change in IOD amplitude is driven solely by the forced response and not influenced by internal variability, the distribution would converge to the long-term trend sensitivity. In all high-CO2 simulations, the transient sensitivity shows widespread distribution that deviates from its long-term trend sensitivity (Fig. 2b and Supplementary Fig. 8). The distribution range of the transient sensitivity differs between the model and simulations. For example, the CCSM3II abrupt700ppm ranges from −92%/°C to 50%/°C and the IPSL-CM5A-LR abrupt4x ranges from −37%/°C to 13%/°C (mean minus/plus one standard deviation). However, in all high-CO2 simulations, the transient sensitivity ranges from negative to positive value, and its lower/upper bound is much larger than its long-term trend sensitivity.

Consistent with the IOD amplitude change, the intensities of both positive and negative IOD events also show long-term decreasing trends with increasing GMST (15 of 18 high-CO2 simulations for positive IOD and 15 of 18 for negative IOD), but show large fluctuations and widespread transient sensitivity to GMST (Supplementary Figs. 9 and 10).

This suggests that internal variability significantly influences the inter-centennial IOD variability change to GMST increase and masks the long-term response. The masking effect of internal variability becomes stronger with increasing strength and decreasing frequency of internal variability (Supplementary Discussion 2).

Physical sources of the IOD change
To investigate the physical sources of the forced change and internal variability of IOD amplitude, we utilize a simple model for the IOD20,31,32 (“Methods” section). The simple IOD model includes three key parameters/terms that characterize IOD variability, namely, local feedback processes (λ), external El Niño–Southern Oscillation (ENSO) forcing (β×ENSO), and weather noise (σ). The simple IOD model, which reasonably reproduces the IOD variability (blue lines in Fig. 3f, g), disentangles the key physical sources of IOD as each parameter (“Methods” section and Supplementary Figs. 11–14) and allows us to quantify the contribution of each physical source to the change in IOD amplitude. All three parameters/terms, λ, β×ENSO, and σ are robustly decreased in the equilibrium period of high-CO2 simulations (Fig. 3a and Supplementary Fig. 11) and exhibit considerable internal variability (Fig. 3e and Supplementary Figs. 12–14) (“Methods” section).

We perform two different sets of experiments with the simple IOD model. First, we perform the parameter perturbation experiment (“Methods” section) to quantify how much each forced change in λ, β×ENSO, and σ contributes to the forced change in the IOD amplitude. Second, we perform the low-pass parameter experiment (“Methods” section) to quantify how much each of the internal variability of λ, β×ENSO, and σ contributes to the internal variability of the IOD amplitude.

The results of the parameter perturbation experiment (Fig. 3f) show that the decrease in λ is the most dominant factor for the forced IOD amplitude change, and the decrease in σ is a secondary contributor. The decrease in β×ENSO has a minimal contribution. This shows that greenhouse warming suppresses the local feedback processes in the tropical Indian Ocean, leading to the decrease in IOD amplitude. However, the weakened ENSO forcing makes virtually no contribution to the IOD change, even though ENSO is known to be a major forcing of IOD31,33. The change in local feedback processes is associated with the changes in mean climatological state of the tropical Indian Ocean17,18,20 (detailed discussion can be found in Supplementary Discussion 3).

In contrast, the results of the low-pass experiments (Fig. 3g) show that the internal variability in β×ENSO is the most dominant factor for the internal variability of the IOD amplitude, and λ and σ play secondary roles. This indicates that the inter-centennial change in ENSO forcing is a primary source of the internal variability of IOD amplitude, while the local feedback processes have a lesser influence. Thus, the masking effect of internal variability, which leads to the widespread response of IOD amplitude, is mainly attributed to the remote influence of ENSO.

The remote ENSO forcing is a combined effect of ENSO variability (ENSO) and sensitivity of IOD to ENSO variability (β). We perform additional experiments to separately quantify their contribution to the forced change and the internal variability of the IOD amplitude (“Methods” section). The changes in both ENSO and β make virtually no contribution to the forced decrease in IOD amplitude (Supplementary Fig. 15). The decrease in ENSO variability makes a larger contribution than β, but the contribution is very small compared to λ and σ. The internal variability of ENSO alone makes a considerable contribution to the internal variability of the IOD amplitude, similar in magnitude to that of λ and σ, and β makes a slightly smaller contribution (Supplementary Fig. 16). This shows that the large internal variability of β×ENSO and its influence on IOD is attributed to combined internal variability of ENSO and β.

In summary, these results from the simple IOD model experiments show that the internal variability of the IOD amplitude is primarily caused by the remote influence of ENSO, while the forced weakening of the IOD amplitude is mainly induced by the suppression of local feedback processes in the tropical Indian Ocean.

Remote ENSO influence on the IOD change
The quantification by the simple IOD model experiments shows that the remote ENSO forcing is the dominant source of the internal variability of the IOD amplitude, while it has a very small contribution to the forced IOD amplitude change. This contrasting role of the ENSO forcing can be understood in physical terms.

The remote ENSO forcing on the IOD occurs by altering the atmospheric circulation in the tropical Indian Ocean basin17,31,34,35. The example ENSO regression pattern for wind stress and SST anomaly in the tropical Indian Ocean is shown in Supplementary Fig. 17. The warm (cold) phase of ENSO induces an anomalous easterly (westerly) wind stress in the eastern and central tropical Indian Ocean. This anomalous easterly (westerly) wind stress strengthens the wind-thermocline-SST feedback which acts as positive feedback, and consequently facilitates the development of the positive (negative) phase of the IOD. Thus, the ENSO-induced IOD variability is linked to the wind stress variability in the eastern and central tropical Indian Ocean.

In the high-CO2 simulations, the eastern and central wind stress variability in the tropical Indian Ocean is robustly weakened (Supplementary Fig. 18). This forced weakening of the wind stress variability is correlated with the decrease in the remote ENSO forcing (Supplementary Fig. 18) and also with the ENSO-induced decrease in IOD amplitude (Supplementary Fig. 19). This suggests that the forced weakening of the wind stress variability and has a subsequent effect on the decrease in IOD amplitude. However, the reduced ENSO forcing actually has a very small contribution to the decrease in IOD amplitude as shown in the simple IOD model experiments. This is because the IOD feedback processes are related not only to the wind stress variability, but also to multiple mean states of the tropical Indian Ocean, including thermocline depth, zonal SST gradient, and surface-subsurface coupling strength (e.g., sensitivity of SST to thermocline depth perturbation) (Supplementary Discussion 3). The forced change in these multiple factors would also contribute to the change in the feedback process and the associated IOD amplitude. Thus, the ENSO-induced weakening of the wind stress variability and associated wind-thermocline-SST feedback alone has a limited impact on the net feedback process, and is therefore unlikely to induce a significant change in IOD amplitude.

In contrast, the eastern and central wind stress variability can considerably contribute to the internal variability of IOD amplitude. The eastern and central wind stress variability shows a very high level of internal variability (Supplementary Fig. 20). The internal variability of eastern and central wind stress is associated with the internal variability of ENSO and its remote forcing on IOD (Supplementary Fig. 21), and is closely linked with the internal variability of IOD amplitude (Supplementary Figs. 22 and 23). This suggests that the internal variability of ENSO influences the eastern and central wind stress variability in the tropical Indian Ocean, and contributes to the internal variability of IOD amplitude. The strong internal variability of ENSO and the associated remote forcing by the anomalous wind stress variability make the ENSO forcing the largest contributor to the internal variability of IOD amplitude.

Methods
LongRunMIP data
We use simulation output from the LongRunMIP archive. The LongRunMIP archive includes a total of 11 models that provide monthly surface temperature (or monthly surface air temperature) variables longer than 500 years. We evaluate the IOD simulation performance of the 11 models, and select 9 models that can reasonably simulate the observed spatial and temporal characteristics of IOD. The model evaluation process and selection criteria are explained in Supplementary Discussion 1. The selected models are listed Supplementary Table 1.

The selected analysis sample includes a total of 27 simulations from 9 models. The simulation sample consists of 9 control and 18 high-CO2 simulations (each model has a single number of control experiments). The control simulation is preindustrial run with a constant 1850 forcing. The high-CO2 simulation includes the following forcing scenarios: instantaneous doubling (abrupt2x), quadrupling (abrupt4x), and octupling (abrupt8x); and gradual increase of 1% per year until doubling (1pct2x) and quadrupling (1pct4x) of atmospheric CO2 levels relative to the control simulation. As an exception, CCSM3II includes the abrupt700ppm and abrupt1400ppm experiments. These experiments abruptly increased CO2 levels to 700 ppm and 1400 ppm from the control level, respectively. After the increase in the CO2 level, the CO2 level remains unchanged until the end of the simulation. The length of the simulation differs by simulation (Supplementary Table 1), but all simulations are longer than 500 years. The analysis sample of the high-CO2 simulations does not include their control simulation period. Details of the LongRunMIP protocol can be found in ref. 27.

Monthly surface temperature (variable name: ts) is used as a proxy for the SST since monthly SST is not consistently available (only a few models offer monthly SST data). Ref. 28 also applied a similar method to study ENSO amplitude changes using the LongRunMIP archive. For GISS-E2-R, monthly surface air temperature (variable name: tas) is utilized instead of monthly surface temperature, which is unavailable.

We use annual ocean potential temperature (variable name: thetao) and zonal wind stress (variable name: tauuo) variables. The monthly ocean potential temperature and zonal wind stress are not consistently available (only a few models offer monthly data). The ocean potential temperature is used for the calculation of the thermocline depth. The ocean potential temperature is unavailable for the CCSM3II model, and the zonal wind stress is unavailable for the CCSM3II and GISS-E2-R models. Therefore, these models are excluded from the results related to the thermocline depth and zonal wind stress.

CMIP6 data
We use the piControl experiment from the CMIP6 archive36. A total of 36 models are used. The model name, ensemble, and simulation length are listed in Supplementary Table 2. We use monthly SST (variable name: tos) for the DMI calculation.

Calculation of the DMI and Niño 3.4 index
The DMI is calculated for the LongRunMIP and CMIP6 simulation output. For the LongRunMIP, we applied quadratic detrending to the monthly SST to remove the forced long-term change (i.e., secular mode). The detrending is performed for the entire period of the simulation (see Supplementary Table 1 for the length of each simulation). The anomaly is defined as a deviation from this long-term trend. The area-averaged monthly SST anomaly is calculated for the western (50°E–70°E and 10°S–10°N) and eastern (90°E–110°E and 10°S–0°N) tropical Indian Ocean. The DMI is calculated as the west minus east monthly SST anomaly. The Niño 3.4 index is calculated as the area-averaged monthly SST anomaly for the Niño 3.4 region (170°W–120°W and 5°S–5°N). As a result, for each simulation, we obtain the DMI and Niño 3.4 index time series, the length of which is equal to the length of the entire simulation. Using these DMI and Niño 3.4 index time series, we perform analyses including the calculation of IOD amplitude and event intensity.

IOD and ENSO amplitude
The amplitude of IOD and ENSO is defined as the standard deviation of the DMI and Niño 3.4 index, respectively.

Detection of IOD events
The positive (negative) IOD event is defined as a three-month running mean of DMI exceeding (below) the one standard deviation of the time series for at least three consecutive months. The intensity of a positive (negative) IOD event is the highest (lowest) DMI level during the event.

Definition of the equilibrium and transient response
The equilibrium IOD response is defined as the response of IOD variability to an increased equilibrium GMST level. The equilibrium response reflects the forced equilibrium response due to greenhouse gas warming. The transient response is defined as the evolution of IOD variability in response to increasing GMST with time. The transient response reflects both the forced response and internal variability. Specifically, we interpret the long-term trend in IOD amplitude (or IOD event intensity) as the forced response, and the fluctuation in IOD amplitude (or IOD event intensity) that deviates from the long-term forced trend as the internal variability.

Spectrum analysis
The spectrum analysis is performed using the Multi-Taper Method (MTM), a widely used technique for spectrum estimation. The parameters of the MTM analysis are chosen by following the standard guidelines of ref. 37,38. We set the time-frequency bandwidth parameter (NW) to 2, and the number of Slepian Tapers (K) to 3. The significance of the power spectrum is estimated based on the red noise, which is a typical null hypothesis for climate variability. The 99%, 95%, and 50% significance levels are calculated using the Monte Carlo method with 100 ensembles.

Equilibrium sensitivity
We calculate the equilibrium sensitivity of the IOD amplitude to GMST. We use the IOD amplitude and GMST in the equilibrium period (i.e., the last 500 years of the simulation). For each model, we calculate the linear regression coefficients of {DMI}eq = aeq {GMST}eq + b where {DMI}eq and {GMST}eq are an array of IOD amplitude and GMST for the control and high-CO2 simulations, respectively. The aeq is the equilibrium sensitivity. For {DMI}eq, we use percentage change from the control simulation, not the absolute value. For example, in the equilibrium period, the CCSM3 model exhibits IOD amplitudes of 0.26°C, 0.17°C, 0.14°C, and 0.13°C for the control, abrupt2x, abrupt4x, and abrupt8x simulations, respectively. The corresponding percentage changes in the IOD amplitude from the control simulations are 0%, −34%, −47%, and −50%, respectively. The corresponding GMSTs are 14.3°C, 16.6°C, 18.9°C, and 21.5°C, respectively. We perform linear regression on this array of IOD amplitude ({DMI}eq = {0%, −34%, −47%, −50%}) and GMST ({GMST}eq = {14.3°C, 16.6°C, 18.9°C, 21.5°C}) and estimate the aeq. The unit of aeq is %/°C. This equilibrium sensitivity quantifies the percentage change in IOD amplitude per 1°C increase in GMST under an equilibrium climate condition.

Transient sensitivity
We calculate the transient sensitivity of IOD amplitude to GMST. We use times series of 100-year moving IOD amplitude and GMST. For each high-CO2 simulation, we perform moving window linear regression for {DMI}ts = ats {GMST}ts + b where {DMI}ts and {GMST}ts are time series of 100-year moving IOD amplitude and GMST. The ats is the transient sensitivity. For {DMI}ts, we use percentage change from the initial value, not the absolute value. The unit of ats is %/°C. The size of moving window is 1°C. The statistics of the ats distribution (mean and standard deviation) is presented in Fig. 2b. The full probability distribution of the ats is presented in Supplementary Fig. 8. To derive the probability distribution of the ats, we perform kernel density estimation on the ats. The bandwidth is set as the optimal value for normal densities.

Calculation of thermocline depth
The thermocline depth is calculated as the maximum gradient depth of the ocean potential temperature profile. For each horizontal grid point, we perform a cubic interpolation on the vertical ocean potential temperature profile. The maximum temperature gradient depth is determined in the interpolated ocean potential temperature profile.

Simple IOD model
The simple IOD model20,31 is written as

dT/dt = λT + βENSO(t) + σξ

where T is the DMI and ENSO is the Niño 3.4 index. The simple IOD model includes three terms representing local feedback processes (λT), external ENSO forcing (βENSO(t)), and stochastic forcing (σξ). λ characterizes local feedback processes associated with the IOD. β characterizes the sensitivity of the DMI to external ENSO forcing. σ is the noise amplitude and ξ is white noise with zero mean and unit standard deviation. We consider the seasonality of λ, β, and σ. Thus, these three parameters vary seasonally. The model can consider the lag between ENSO and T because the model equation relates the ENSO with the tendency of T, not with the T itself.

We note that there are several versions of the simple IOD model in the current literature depending on the types of feedback processes included. The simple IOD model used in this study is the basic type model introduced in ref. 31, which considers the essential physical process of IOD in the simplest way. Hence, the model does not explicitly account for the time-delayed Rossby wave process as done in ref. 20, and implicitly considers it by λT, which comprehensively represents local feedback processes in the tropical Indian Ocean. The simplicity of the model used in this study has advantages that facilitate the physical understanding of the model results, especially in determining whether the IOD amplitude change is caused by a local process or a remote influence.

Fitting the parameters of the simple IOD model
We perform parameter fitting of the simple IOD model for the given monthly DMI and Niño 3.4 index time series. We estimate the values of three parameters (λ, β, σ) for each month of the year from January to December. We discretize the simple IOD model as follows:

ΔTi/Δt = λTi + β(ENSO)i

where i denote time step (i = 1, 2…N) and ΔTi = Ti+1 − Ti. Δt is 1 month. We perform the linear regression and estimate λ and β using yearly time series of ΔTi, Ti, and (ENSO)i for each month of the year. The σ is estimated as the standard deviation of the residual of the linear regression.

Solving the simple IOD model
The simple IOD model is numerically solved using the Euler–Maruyama method with a time step of 0.1 months. The Euler–Huen method is a numerical scheme for solving the Stochastic Differential Equation. The input data is seasonally varying parameter and the Niño 3.4 index. We performed the linear interpolation for the parameters and the Niño 3.4 index, which have a time resolution of 1.0 month, and downscaled the time resolution to 0.1 month to match the time resolution with the numerical time step.

Parameter perturbation experiment
We perform the parameter perturbation experiment using the simple IOD model. The parameters of the simple IOD model are fitted using the DMI and Niño 3.4 index time series in the equilibrium period. The fitted parameters are shown in Supplementary Fig. 11 and Fig. 3a. The statistical relationships between the fitted parameters and the GMST are shown in Supplementary Fig. 24. We reproduce the DMI time series by running the simple IOD model using the fitted parameters (EXP-all-high-CO2) for each LongRunMIP simulation. We then run the simple model with the control simulation parameters but with the replaced λ (EXP-λ-high-CO2), β (EXP-β-high-CO2), ENSO (EXP-ENSO-high-CO2), β and ENSO (EXP-β*ENSO-high-CO2), and σ (EXP-σ-high-CO2) values with the ones from the high-CO2 simulations. For example, the EXP-λ-high-CO2 uses λ from the high-CO2 simulation and the β, ENSO, and σ from its control simulation. The parameter perturbation experiment is performed for each high-CO2 simulation. Additionally, for the baseline experiment, we run the simple model using the control simulation parameters (EXP-all-control).

Low-pass parameter experiment
We perform the low-pass parameter experiment using the simple IOD model. The model parameters are fitted for a moving 100-year window using the DMI and Niño 3.4 index time series. The annual mean value of the fitted parameters is shown in Supplementary Figs. 12–14 (note that the parameter fitting is performed for each calendar month as shown in Supplementary Fig. 11, but we only present their annual mean value for simplicity). For each segment of a window, we reproduced the DMI time series using the fitted parameters (EXP-all-varying). To remove the fluctuations of the parameters, we low-pass filter λ, βENSO, and σ and only leave their quadratic trend. We run the EXP-all-varying but with the replaced original parameter from the low-pass filtered parameters for λ (EXP-λ-varying), β (EXP-β-varying), ENSO (EXP-ENSO-varying), β and ENSO (EXP-β*ENSO-varying), and σ (EXP-σ-varying). For example, the EXP-λ-varying uses the original λ and low-pass filtered β, ENSO, and σ. The strength of the internal variability of IOD amplitude is quantified as the standard deviation of the IOD amplitude. The low-pass parameter experiment is performed for each high-CO2 simulation. Additionally, for the baseline experiment, we run the simple model solely using only the low-pass parameters (EXP-all-fixed).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors invesfigated the variability of IOD amplitude under prolonged greenhouse warming based\non model simulafions, and found that IOD variability weakens under prolonged greenhouse warming due\nto the suppressed local feedback processes in tropical Indian Ocean. On centennial fime scale, the\nweakening trend is disturbed by the internal variability associated with ENSO. Although these results are\ninteresfing and show the possible response of IOD to global warming on millennial-length scale, some\nmajor conclusions are not convincing and further evidences on physical processes are needed.\nneeds to be carefully verified. However, the models are not validated and selected before they were\nused for analysis in the manuscript. This may affect the credibility of the results. The authors should\nselect the models based on the fidelity of model results on the temporal and spafial characterisfics of\nIOD. Objecfive criteria are needed for model selecfion.\nsimulafions (Supplementary Fig. 7) seems to be different from that in observafions. In observafions, IOD\naccounts for the second EOF mode rather than the first EOF mode. In addifion, the SST variafion is\nstrongest in the off-Sumatra region in observafions, while in simulafions there are two regions with\nstrong SST variafions. This may deserve further discussion.\nmechanisms. The authors concluded that (L221) internal variability of the IOD amplitude is primarily\ncaused by the remote influence of ENSO, while the forced weakening of the IOD amplitude is mainly\ninduced by the suppression of local feedback processes in tropical Indian Ocean. Why is ENSO important\non centennial fime scale but not important on millennial fime scale? What are the underlying\nmechanisms?\n(L216) The internal variability in β×ENSO is the most dominant factor for the internal variability of IOD\namplitude. How does ENSO influence the IOD amplitude? Is it due to variafions in ENSO intensity or the\nvariafions in the connecfion between ENSO and IOD?\nlike change which weakens the local feedback processes. (Supplementary Discussion 2). But previous\nstudies show that such IOD-like paftern could enhance the thermocline–SST feedback. What is the cause\nfor such difference? Bjerknes feedback is the major feedback process in IOD development. Are the\nvariafions of the three components of Bjerknes feedback consistent with the conclusion?\ne.g. L252 Furthermore, our study suggests the possibility of dramafic changes in the characterisfics of\nIOD in a very warm climate environment, such as the occurrence of double peaks in the boreal spring\nand fall, and a zonally uniform oscillafion paftern.\nAre these results significant?\nresponse, need to be further clarified, make them easier for readers to understand.", "reviewer_2": "Review for submission “Decreased Indian Ocean Dipole variability under prolonged\ngreenhouse warming”\nThis study ufilized LongrunMIP to address the forced long-term changes in the Indian Ocean Dipole (IOD)\namplitude arising from natural variability. The manuscript is well-wriften, and the arguments are\nsupported. Nevertheless, there are some concerns that the authors should address in the revision.\n1. About the natural variability. Line 80: \"IOD amplitude in the twenty-first century either increases or\ndecreases depending on the pathways of internal variability, despite the same warming forcing.\" The\nterm 'pathways' may be more closely associated with forcing rather than natural variability. Please\nconsider revising for clarity.\"\n2. About the simple IOD model. In paper 20, four terms are discussed: local feedback, delayed feedback,\nENSO impact, and stochasfic forcing. In this study, the authors only consider three terms, excluding the\ndelayed feedback. Please provide more discussions or explanafions on this considerafion. Line 381: “We\nconsider the seasonality of λ, β, and σ. Thus, these three parameters vary seasonally. The typical value of\nthe annual mean λ, esfimated from observafions and climate models in previous studies, is negafive.' It\nis assumed here that the model is fifted for each month with ENSO and DMI index in phase (lag-0).\nWhether you consider a lead-lag relafionship? Also, for the annual mean model, does the equafion sfill\nhold?\"\n3. About the seasonal changes. Paper 19 addressed that the early IOD is more frequently observed\nunder warming condifions. In this study, Line 124 menfions, 'The seasonal variability of IOD decreases in\nalmost all seasons (Supplementary Fig. 5).' Moreover, Figure S5 illustrates that HadCM3L and a few other\nmodels peak early in JJA. It would be beneficial to delve into a more detailed discussion of this change.\nLines 135-146: You might consider separafing the discussions on seasonal changes and spafial pafterns\ninto two paragraphs for clarity.\n4. About the relafionship between changes in IOD and ENSO. Paper 28 addressed long-term ENSO\nchanges based on LongrunMIP. It would be valuable to explore whether models agree on changes in\nboth IOD and ENSO and if the dominafing mechanisms are similar. Authors should engage in discussions\non this aspect.\n5. About the definifion of IOD amplitude. Line 102: 'We measure the IOD amplitude as the standard\ndeviafion of the DMI.' Typically, the community uses the peak season (SON) DMI to gauge IOD changes.\nHere, monthly DMI is employed; the authors may wish to explain this choice.\n6. About transient response. Line 95: For the 1pct2x and 1pct4x runs, 'After the CO2 increase, the forcing\nis stabilized for typically longer than 1,000 years.' While the first 200 years may be considered the\ntransient period, the runs predominantly reflect the equilibrated response after the forcing stabilizes.\nLine 148: 'We analyze the transient response of the IOD amplitude to the CO2 forcing. We calculate the\n100-year moving IOD amplitude for the enfire simulafion period.' It is noted that the enfire period (1000\nyears or more) covers both the transient and equilibrated periods. The authors may need to consider\ntaking the first 100-200 years as the transient period when the forcing is changing. After the forcing\nstabilizes, it may not be considered part of the transient period.\n7. About the period of analysis. Line 107: 'last 500 years.' Are the enfire 500 years used for detrending\nthrough quadrafic fifting? Line 312: 'For the LongRunMIP, we applied quadrafic detrending to the\nmonthly SST to remove the forced long-term change (i.e., secular mode).' Could the authors clarify\nwhether they used a 100-year or 500-year period for this detrending process?\n8. Figure S13. The upper panel x-y axis should be swapped to be consistent with the lower panel."}, "author_response": "Author’s Reply to the Reviewer #1\nReviewer #1 (Remarks to the Author):\nThe authors investigated the variability of IOD amplitude under prolonged\ngreenhouse warming based on model simulations, and found that IOD variability\nweakens under prolonged greenhouse warming due to the suppressed local\nfeedback processes in tropical Indian Ocean. On centennial time scale, the\nweakening trend is disturbed by the internal variability associated with ENSO.\nAlthough these results are interesting and show the possible response of IOD to\nglobal warming on millennial-length scale, some major conclusions are not\nconvincing and further evidences on physical processes are needed.\nReply: We appreciate your recognition and valuable comments on our paper.\nThanks to your constructive comments, we have made a substantial improvement\nto the manuscript. The revised manuscript contains some major updates as follows:\n1. We evaluated the performance of the LongRunMIP model and selected models\nthat can reasonably simulate the IOD. We updated the analyses in the paper with\nthese selected models (9 out of 11 models). For the model selection process and\ncriteria, please see the newly added section Supplementary Discussion 1.\n2. We performed the physical mechanism analysis on how ENSO influences IOD.\nPlease see the newly added section Remote ENSO influence on the IOD change.\nThe reply to your comment is attached below. Please note that we will refer to the\ntext line and figure/table number based on the revised version of the manuscript.\nmodels in reproducing IOD needs to be carefully verified. However, the models are\nnot validated and selected before they were used for analysis in the manuscript.\nThis may affect the credibility of the results. The authors should select the models\nbased on the fidelity of model results on the temporal and spatial characteristics\nof IOD. Objective criteria are needed for model selection.\nReply: In the original manuscript, we used all 11 available LongRunMIP models for\nthe IOD analysis. Following your comment, we performed a model performance\nevaluation and identified the model that can reasonably simulate IOD. We selected\n9 out of 11 models that can reasonably simulate the observed characteristics of\nIOD and used them for the analysis.\nThe detailed model evaluation and selection process is fully documented in\nSupplementary Discussion 1. We evaluate model performance based on the DMI,\nthe main analysis variable of this study. We examine three basic characteristics of\nthe IOD: (1) the peak season of the DMI standard deviation, (2) the DMI skewness\nsign, and (3) the spatial regression pattern of the DMI. The peak DMI season is\nrelated to the seasonal phase locking of the IOD, which is associated with various\nphysical processes including the local feedback process and ENSO. Therefore, it is\na good test metric for whether the model correctly simulates the IOD feedback\nprocess and associated variability. The DMI skewness characterizes the nonlinearity\nof the IOD, and is related to the response of the IOD to greenhouse gas warming.\nThe DMI spatial regression pattern shows the spatial pattern of the IOD. We\nevaluate the ability of the model to simulate these three metrics. Based on the\nevaluation, we score the model performance level and classify it into four classes:\nhigh, medium, low, and very low. The summary of the model evaluation is shown\nin Table R1.\nWe select the high and medium performance models for the IOD analysis, and\nexclude the low performance models. Finally, out of 11 models, the following 9\nmodels are selected: CCSM3, CCSM3II, CESM1.0.4, GFDL-CM3, GFDL-ESM2M, GISS-\nE2-R, HadCM3L, IPSL-CM5A and MIROC3.2. The following 2 models are excluded\nfrom the analysis: CNRMCM61 and MPIESM12.\nWe accordingly updated all related material in the manuscript. First, we removed\nthe results of CNRMC61 and MPIESM12, and updated all related figures, tables,\nand metrics. Second, we added an explanation of the model evaluation and\nselection process. We newly added Supplementary Discussion 1 with\nSupplementary Figs. 25 and 26 and Supplementary Table 3. We also accordingly\nadded explanations in L100 and L349.\nOur model selection process does not affect the main results of the paper, but\nrather strengthens them. Table R2 shows the inter-model (or inter-simulation)\nconsensus ratio of the IOD response for three different model groups, all models\n(11 models; model samples in the original manuscript), selected high and medium\nperformance models (9 models; model samples in the revised manuscript), and\nselected high performance models (6 models). As you can see, the selected high\nand medium performance model group shows a higher inter-model consensus on\nthe IOD response than the all models group. In particular, for the selected high\nperformance model group, the inter-model consensus ratio reaches almost 1. This\nincreasing inter-model consensus ratio adds confidence to the results of the paper.\nThe main results of the parameter perturbation experiment (Fig. 3f) and the low-\npass parameter experiment (Fig. 3g) which show the physical source of the IOD\namplitude change, also remain almost unchanged. Therefore, the main results and\nthe associated conclusion of the paper remain the same after applying the model\nselection.\nTable R1. Model evaluation summary. Same as Supplementary Table 3.\nModel\nDMI\nPeak DMI Standard DMI Regression Pattern Performance\nType Product Skewness\nDeviation Season (Peak Positions) Level\nSign\n(Score)\nDipole\nERSSTv5 October Positive (Sumatra-Java coast / -\nEast African coast)\nDipole\nObservation\nCOBE October Positive (Sumatra-Java coast / -\n(1900-2022)\nEast African coast)\nDipole\nHadISST September Positive (Sumatra-Java coast / -\nEast African coast)\nDipole\nCCSM3 September Positive (Sumatra-Java coast / High (3)\nEast African coast)\nDipole\nCCSM3II September Positive (Sumatra-Java coast / High (3)\nEast African coast)\nDipole\nCESM104 October Positive (Sumatra-Java coast / High (3)\nEast African coast)\nDipole\nCNRMCM61 December Negative (Sumatra-Java coast / Low (1)\nEast African coast)\nDipole\nGFDLCM3 September Positive (Sumatra-Java coast / High (3)\nEast African coast)\nDipole\nLongRunMIP\nGFDLESM2M September Negative (Sumatra-Java coast / Medium (2)\n(control)\nEast African coast)\nDipole\nGISSE2R October Positive (Sumatra-Java coast / Medium (2)\nArabian Sea)\nDipole\nHadCM3L October Positive (Sumatra-Java coast / Medium (2)\nArabian Sea)\nDipole\nIPSLCM5A October Positive (Sumatra-Java coast / High (3)\nEast African coast)\nDipole\nMIROC32 September Positive (Sumatra-Java coast / High (3)\nEast African coast)\nDipole\nMPIESM12 August Negative (Sumatra-Java coast / Low (1)\nEast African coast)\nTable R2. The inter-model (or inter-simulation) consensus ratio of the IOD response\nfor different model groups. The consensus ratio for all models (11 models; model\nsamples in the original manuscript), selected high and medium performance\nmodels (9 models; model samples in the revised manuscript), and selected high\nperformance model (6 models).\nSelected Model\nResponse Response in Selected Model\nVariable All Model (High and Medium\nType the Majority (High Performance)\nPerformance)\n17 out of 23 15 out of 18 11 out of 12\nIOD Amplitude Decrease\n(74%) (83%) (92%)\nIOD Event\n20 out of 23 18 out of 18 12 out of 12\nIntensity Decrease\n(87%) (100%) (100%)\n(Positive)\nEquilibrium\nResponse\nIOD Event\n17 out of 23 15 out of 18 11 out of 12\nIntensity Decrease\n(74%) (83%) (92%)\n(Negative)\nIOD Amplitude 9 out of 11 8 out of 9 6 out of 6\nNegative\nSensitivity (82%) (89%) (100%)\nIOD Amplitude\n21 out of 23 17 out of 18 12 out of 12\nSensitivity Negative\n(91%) (94%) (100%)\n(Mean)\nTransient\nResponse\nIOD Amplitude\n20 out of 23 17 out of 18 12 out of 12\nSensitivity Negative\n(87%) (94%) (100%)\n(Long-term Trend)\nOcean in CESM1.0.4 simulations (Supplementary Fig. 7) seems to be different from\nthat in observations. In observations, IOD accounts for the second EOF mode rather\nthan the first EOF mode. In addition, the SST variation is strongest in the off-\nSumatra region in observations, while in simulations there are two regions with\nstrong SST variations. This may deserve further discussion.\nReply: Following your comment, we carefully reexamined the EOF analysis results\nof the CESM1.0.4 simulations.\n(1) The first and second EOF mode\nAs we explain in the Methods section, we use surface temperature (variable name:\nts) as a proxy for the SST (variable name: tos) because monthly SST is not\nconsistently available (only a few models provide monthly SST data). It has come\nto our attention that we did not apply the land masking to the surface temperature\nvariable when performing the EOF analysis on it. Thus, the original EOF analysis\nresults included both land and sea surface temperature anomalies (we applied the\nland masking for all other calculations, but the CESM1.0.4 EOF analysis was the\nonly exception).\nTo correct this, we applied the land masking to the surface temperature field and\nre-performed the EOF analysis. Please see Fig. R1.\nFig. R1. EOF analysis of the CESM1.0.4 surface temperature anomaly. The first\nand second EOF modes of the monthly surface temperature anomaly over the\ntropical Indian Ocean (40° E-120° E and 20° S-20° N). The results for the control\n(a), abrupt2x (b), abrupt4x (c), and abrupt8x (d). The upper panel is the first mode\nand the lower panel is the second mode. Variance explained by the mode is marked\nin the title of each panel.\nAs you can see, the IOD now accounts for the second EOF mode, consistent with\nthe observation.\nWe accordingly revised the figure and its caption with the land mask applied result\nshown in Fig. R1. Please see the updated Supplementary Fig. 7. We slightly\nextended the eastern edge of the analysis domain from 110° E to 120° E to more\nclearly show and discuss the secondary strong SST location, that you pointed out\nin your comment. We also removed the EOF analysis result for the SON surface\ntemperature anomaly for a fair comparison between the simulations, as the\nseasonal phase locking pattern is significantly different in the abrupt8x (secondary\nseasonal peak emerges in March-April-May).\n(2) Difference in the spatial pattern\nAs you pointed out, there is a difference in the EOF pattern between the CESM and\nthe observation. We added a discussion of this difference in the caption of\nSupplementary Fig. 7 as (we also considered adding it in the main text, but given\nthe flow of the paper, we decided to add it in the figure caption):\n“We note that the spatial pattern of the second EOF mode of the CESM1.0.4 shows\nsome difference from the typical observation. Both the CESM1.0.4 and the historical\nobservational records show a strong SST EOF pattern in the off-Sumatra region,\nbut the CESM1.0.4 also shows a strong SST in the Sumatra coastal region. The\ndifference can be partly attributed to the bias of the CESM1.0.4 in simulating the\nmean climatological state of the tropical Indian Ocean, including the wind stress\nfield, the SST pattern, and the ocean potential temperature profile. A detailed\ndiscussion of this subject is beyond the scope of this study which focuses primarily\non the response of the tropical Indian Ocean to greenhouse gas warming. Despite\nsuch bias in the mean state, our conclusion in this paper is well supported by the\nstrong inter-model consensus in the IOD response.”\nThank you.\nsupported by physical mechanisms. The authors concluded that (L221) internal\nvariability of the IOD amplitude is primarily caused by the remote influence of ENSO,\nwhile the forced weakening of the IOD amplitude is mainly induced by the\nsuppression of local feedback processes in tropical Indian Ocean. Why is ENSO\nimportant on centennial time scale but not important on millennial time scale?\nWhat are the underlying mechanisms?\n(L216) The internal variability in β×ENSO is the most dominant factor for the\ninternal variability of IOD amplitude. How does ENSO influence the IOD amplitude?\nIs it due to variations in ENSO intensity or the variations in the connection between\nENSO and IOD?\nReply: Following your comment, we performed a series of additional experiments\nand analyses.\n(1) Disentangling the role of ENSO intensity and IOD sensitivity to ENSO\nIn the original version of the paper, we quantified the contribution of remote ENSO\nforcing (β×ENSO) to the forced change and internal variability of IOD amplitude.\nTo separate the contribution of ENSO variability (ENSO) and the IOD sensitivity to\nENSO variability (β), we performed additional experiments. We repeated the\nprevious parameter perturbation and low-pass parameter experiment, but changed\nENSO and β separately.\nPlease see the newly added experiment result in L245 and associated\nSupplementary Figs. 15 and 16. The added explanation of the result is read as\nfollows:\n“The remote ENSO forcing is a combined effect of ENSO variability (ENSO) and\nsensitivity of IOD to ENSO variability (β). We perform additional experiments to\nseparately quantify their contribution to the forced change and the internal\nvariability of the IOD amplitude (Methods). The changes in both ENSO and β make\nvirtually no contribution to the forced decrease in IOD amplitude (Supplementary\nFig. 15). The decrease in ENSO variability makes a larger contribution than β, but\nthe contribution is very small compared to λ and σ. The internal variability of ENSO\nalone makes a considerable contribution to the internal variability of the IOD\namplitude, similar in magnitude to that of λ and σ, and β makes a slightly smaller\ncontribution (Supplementary Fig. 16). This shows that the large internal variability\nof β×ENSO and its influence on IOD is attributed to combined internal variability\nof ENSO and β.”\nIn answer to your question, the additional experiment shows that the internal\nvariability of β×ENSO is caused by contribution from both internal variability of\nENSO amplitude (i.e., the standard deviation of ENSO) and internal variability of the\nsensitivity of IOD to ENSO variability (β). The internal variability of ENSO amplitude\nis the primary contributor, and β is the secondary effect, however, they have a\nsimilar magnitude of contribution.\n(2) Physical mechanism of the remote ENSO forcing\nWe performed the analysis on the physical mechanism of the remote ENSO forcing\non the IOD amplitude change. Please see the newly added section, “Remote ENSO\ninfluence on the IOD change” in L260 and associated Supplementary Figs. 17-\n23. The added section is read as follows.\n“Remote ENSO influence on the IOD change\nThe quantification by the simple IOD model experiments shows that the remote\nENSO forcing is the dominant source of the internal variability of the IOD amplitude,\nwhile it has a very small contribution to the forced IOD amplitude change.\nThis contrasting role of the ENSO forcing can be understood in physical terms.\nThe remote ENSO forcing on the IOD occurs by altering the atmospheric circulation\nin the tropical Indian Ocean basin. The example ENSO regression pattern for wind\nstress and SST anomaly in the tropical Indian Ocean is shown in Supplementary\nFig. 17. The warm (cold) phase of ENSO induces an anomalous easterly (westerly)\nwind stress in the eastern and central tropical Indian Ocean. This anomalous\neasterly (westerly) wind stress strengthens the wind-thermocline-SST feedback\nwhich acts as positive feedback, and consequently facilitates the development of\nthe positive (negative) phase of the IOD. Thus, the ENSO-induced IOD variability is\nlinked to the wind stress variability in the eastern and central tropical Indian Ocean.\nIn the high-CO2 simulations, the eastern and central wind stress variability in the\ntropical Indian Ocean is robustly weakened (Supplementary Fig. 18). This forced\nweakening of the wind stress variability is correlated with the decrease in the\nremote ENSO forcing (Supplementary Fig. 18) and also with the ENSO-induced\ndecrease in IOD amplitude (Supplementary Fig. 19). This suggests that the forced\ndecrease in ENSO forcing weakens the wind stress variability in the tropical Indian\nOcean and has a subsequent effect on the decrease in IOD amplitude. However,\nthe reduced ENSO forcing actually has a very small contribution to the decrease in\nIOD amplitude as shown in the simple IOD model experiments. This is because the\nIOD feedback processes are related not only to the wind stress variability, but also\nto multiple mean states of the tropical Indian Ocean, including thermocline depth,\nzonal SST gradient, and surface-subsurface coupling strength (e.g., sensitivity of\nSST to thermocline depth perturbation) (Supplementary Discussion 3). The forced\nchange in these multiple factors would also contribute to the change in the\nfeedback process and the associated IOD amplitude. Thus, the ENSO-induced\nweakening of the wind stress variability and associated wind-thermocline-SST\nfeedback alone has a limited impact on the net feedback process, and is therefore\nunlikely to induce a significant change in IOD amplitude.\nIn contrast, the eastern and central wind stress variability can considerably\ncontribute to the internal variability of IOD amplitude. The eastern and central wind\nstress variability shows a very high level of internal variability (Supplementary Fig.\n20). The internal variability of eastern and central wind stress is associated with the\ninternal variability of ENSO and its remote forcing on IOD (Supplementary Fig. 21),\nand is closely linked with the internal variability of IOD amplitude (Supplementary\nFigs. 22 and 23). This suggests that the internal variability of ENSO influences the\neastern and central wind stress variability in the tropical Indian Ocean, and\ncontributes to the internal variability of IOD amplitude. The strong internal\nvariability of ENSO and the associated remote forcing by the anomalous wind stress\nvariability make the ENSO forcing the largest contributor to the internal variability\nof IOD amplitude.”\nIn addition to the newly added section, we would like to specifically answer your\nquestions.\n- How does ENSO influence the IOD amplitude?\nThe ENSO forcing influences the eastern and central wind stress in the tropical\nIndian Ocean by altering the Walker circulation pattern. This anomalous wind in the\ntropical Indian Ocean facilitates the development of IOD. Specifically, the warm\n(cold) phase of ENSO induces an anomalous easterly (westerly) wind stress in the\neastern and central tropical Indian Ocean. This anomalous easterly (westerly) wind\nstress strengthens the wind-thermocline-SST feedback which acts as positive\nfeedback, and facilitates the development of the positive (negative) phase of the\nIOD. Consequently, the ENSO forcing increases the IOD amplitude.\n- Why is ENSO important on centennial time scale but not important on millennial\ntime scale? What are the underlying mechanisms?\nOn the centennial time scale, internal variability dominates the change in IOD\namplitude as shown in the “Transient IOD response” section of the paper. The\n100-year moving IOD amplitude fluctuates strongly with time under greenhouse\ngas warming. In contrast, on the millennial time scale, the forced change in IOD\namplitude is detectable and shapes the IOD amplitude change. The forced change\nin IOD amplitude is consistently detectable with the long-term trend in the 100-\nyear moving IOD amplitude (as shown in the “Transient IOD response” section)\nand the change in IOD amplitude under a different equilibrium warming level (as\nshown in the “Equilibrium IOD response” section). Therefore, the change in IOD\namplitude over the centennial and millennial timescales can be translated as\ninternal variability and forced change of IOD amplitude, respectively.\nENSO-induced IOD variability is linked to wind stress variability in the eastern and\ncentral tropical Indian Ocean. As we explained in the newly added section, the\nforced change in wind stress variability alone has a limited role in altering the net\nIOD feedback process, which involves with multiple physical factors. The decrease\nin wind stress variability possibly contributes to the decrease in IOD amplitude, but\nthe series of analyses and experiments (explained in “Physical sources of the IOD\nchange” and “Remote ENSO influence on the IOD change”) show that the\nmagnitude is actually very small. Thus, ENSO does not have a significant influence\non the forced IOD change that occurs on a millennial time scale.\nIn contrast, the wind stress variability shows a very strong degree of internal\nvariability compared to other variables, which means that its centennial-scale\ntemporal fluctuation is very large (we define internal variability as anomalous short-\nterm fluctuation in the variable that deviates from the long-term forced change).\nThe strength of the temporal fluctuation of the IOD amplitude is determined by\nthe strength of the temporal fluctuation of the IOD feedback processes. Thus, the\nwind stress variability, which exhibits a considerable degree of fluctuation, can\nconsiderably contribute to the fluctuation in the IOD feedback processes and the\nIOD amplitude. Consequently, internal variability of ENSO forcing via the wind stress\nin the tropical Indian Ocean, contributes significantly to the internal variability of\nIOD amplitude, that occurs on a centennial timescale. On this internal variability of\nENSO forcing, the internal variability of ENSO amplitude is the primary contributor,\nand the internal variability of IOD sensitivity to ENSO is the secondary effect.\nshows a positive IOD-like change which weakens the local feedback processes.\n(Supplementary Discussion 2). But previous studies show that such IOD-like pattern\ncould enhance the thermocline–SST feedback. What is the cause for such difference?\nBjerknes feedback is the major feedback process in IOD development. Are the\nvariations of the three components of Bjerknes feedback consistent with the\nconclusion?\nReply: We assume that you are referring to the line “Conversely, the reduced east-\nto-west SST gradient (i.e., positive IOD-like change) (Fig. 3c) would weaken the\nsurface zonal advective feedback” in Supplementary Discussion 3 (previously\nSupplementary Discussion 2, but renumbered). Our intention was to discuss the\nindividual effect of a reduced mean zonal SST gradient and associated zonal\nadvective feedback on the net local feedback process. The mean east-to-west SST\ngradient amplifies the development of IOD through the zonal advection process.\nAnomalous easterly (westerly) wind during the positive (negative) IOD development\ncarries anomalous warm water from east to west (west to east), and amplifies the\npositive (negative) IOD development. Thus, the decrease in the mean zonal SST\ngradient weakens this zonal advective feedback, and contributes to the weakening\nof the Bjerkness feedback processes.\nHowever, as you pointed out, the term – ‘positive IOD-like change’ includes not\nonly the change in the zonal SST gradient, but also the change in the thermocline\ndepth and the wind pattern. Therefore, the ‘positive IOD-like change’ is related to\nvarious feedback processes, not only with the zonal advective process, but also with\nsuch as the thermocline-SST feedback. Such comprehensive terminology would\ncause confusion to the readers, especially when discussing the related feedback\nprocesses. To clarify the meaning, we removed the term ‘positive IOD-like change’\nin Supplementary Discussion 3.\nWe carefully reexamined Supplementary Discussion 3 and found that our\ndescription of the Bjerkness feedback process coincides with the previous studies.\n• In response to your concern, the shoaling thermocline, which is part of the\n‘positive IOD-like change’, enhances the thermocline-SST feedback, and contributes\nto the strengthening Bjerkness feedback. This role of thermocline shoaling is\ndescribed in Supplementary Discussion 3 as “Shoaling of the thermocline depth\nincreases the sensitivity of the thermocline depth to the zonal wind stress anomaly\n(Supplementary Fig. 29a). This change would intensify the thermocline-wind\nfeedback”. We revised the term ‘thermocline-wind feedback’ to ‘wind-thermocline-\nSST feedback’ to clarify the feedback loop chain in its term.\nThere is another description of the thermocline-SST feedback in Supplementary\nDiscussion 3. The increased easterly wind, which is part of the ‘positive IOD-like\nchange’ and associated thermocline depth gradient increase strengthens the\nthermocline-SST feedback and contributes to the increase in the Bjerkness feedback.\nIt is described as “The increased climatological west-to-east thermocline depth\ngradient (Fig. 3b) due to the strengthened easterly wind would also enhance the\nthermocline-SST feedback as well.”\n• The ‘positive IOD-like change’ involves changes in multiple mean states including\nthe zonal SST gradient, thermocline depth, and wind stress anomaly. These changes\nin mean states are related to the various parts of Bjerkness feedback processes.\nSome of them suppress the Bjerkness feedback, and some of them strengthen the\nBjerkness feedback. For example, the decreased zonal SST gradient reduces the\nzonal advective feedback and contributes to the weakening of the Bjerkness\nfeedback. In contrast, the decreased thermocline depth intensifies thermocline-\nwind-SST and thermocline-SST feedback processes and contributes to the\nstrengthening of the Bjerkness feedback. The net changes in Bjerkness feedback\nare determined by the sum of these positive/negative effects. Therefore, the change\nin each Bjerkness feedback does not necessarily follow the total net change in the\nIOD feedback process. Such offsetting effect was discussed in Supplementary\nDiscussion 3 and the previous studies.\nneeds to be ascertained.\ne.g. L252 Furthermore, our study suggests the possibility of dramatic changes in\nthe characteristics of IOD in a very warm climate environment, such as the\noccurrence of double peaks in the boreal spring and fall, and a zonally uniform\noscillation pattern.\nAre these results significant?\nReply: First, the L252 refers specifically to the CESM1.0.4 results, which are\nexplained in the “Equilibrium IOD response” section of the paper. The aim of this\nline, which is placed in the discussion section, is to address the implication of this\nstudy rather than to explain the detailed result. In the high-CO simulations of the\nspatial pattern shifts from the dipole- to the monopole-like pattern. Although such\nspecific transitions occur in only one model, it raises the possibility of a dramatic\ntransition of IOD (i.e., critical transition). We believe that the dramatic transition in\nCESM1.0.4 is worth mentioning in the discussion section where an implication of\nthis study is addressed. Given its motivation and purpose, we stated in the L252\n“…suggest the possibility…”.\nTo clarify that the results come specifically from the single model, we revised the\nL252 as “Furthermore, our study suggests the possibility of dramatic changes in\nthe characteristics of IOD in a very warm climate environment, such as the\noccurrence of double peaks in the boreal spring and fall, and a zonally uniform\noscillation pattern, as seen in the CESM1.0.4 simulations.”\nWe assume that your concerns about the L252 are partly based on a\nmisunderstanding, as if it represents common results of the LongRunMIP models.\nHowever, as we explained, the L252 is only related to the single model (which is\nexplained in the “Equilibrium IOD response” section). It is therefore not subject to\nthe inter-model consensus test.\nFollowing your comment, we performed a series of additional statistical significance\ntests on the LongRunMIP model results. This additional statistical significance test\nalso covers the CESM1.0.4 results (i.e., the results in the L252) which you raised\nconcerns about. We also performed additional analysis on the inter-model\nconsensus ratio. The revised parts are as follows:\n(1) Statistical significance\nTo measure the statistical significance of the equilibrium IOD amplitude change, we\nperformed a two-sample F-test of equal variance. For each high-CO simulation\n• IOD amplitude: In all high-CO simulations, the equilibrium IOD amplitude change\n2\nis statistically significant (P<0.05). We accordingly added the result in the caption\nof Fig. 1 as “In all high CO2 simulations, the IOD amplitude change compared to\nthe control simulation is statistically significant (P<0.05). The statistical significance\nis tested with the F-test for a pair of DMI time series of the control and high CO2\nsimulations.”\n• IOD seasonal amplitude: We marked a seasonally significant point (P<0.05) in\nSupplementary Fig. 5. The majority of the high-CO simulation shows the\n2\nsignificant seasonal IOD amplitude change in the peak season. We accordingly\nadded the result in the caption of Supplementary Fig. 5 as “The statistical\nsignificance of the change in the seasonal DMI standard deviation in the high-CO2\nsimulation compared to the control simulation is tested with the F-test. The F-test\nis performed for a pair of DMI time series of the control and high-CO2 simulations\nfor each calendar month. The statistically significant (P<0.05) point is marked as a\nblack circle.”\nTo measure the statistical significance of the CESM1.0.4 EOF pattern changes in the\nhigh-CO simulations, we performed a two-sample F-test. The F-test is performed\n2\non a pair of the reconstructed second mode surface temperature anomalies of the\ncontrol and high-CO2 simulations for each grid point.\n• EOF patterns of CESM1.0.4 simulations: We marked a statistically significant point\n(P<0.05) in Supplementary Fig. 7. In abrupt2x, abrupt4, and abrupt8x simulations,\nthe change is significant in the vast majority of the grid points. We accordingly\nadded the result in the caption of Supplementary Fig. 7 as “The statistical\nsignificance of the change in the second EOF mode surface temperature anomaly\nin the high-CO2 simulation compared to the control simulation is tested with the\nF-test. The F-test is performed for a pair of the reconstructed second mode surface\ntemperature anomalies of the control and high-CO2 simulations for each grid point.\nThe statistically significant (P<0.05) point is marked as a black circle.”\nThe significance test on the seasonal IOD amplitude and CESM1.0.4 EOF patterns\ntogether resolve the concern you raised on the L252.\nAdditionally, to test the statistical significance of the long-term trend in the\ntransient IOD amplitude, we calculated the p-value for the linear trend. We\ncalculated the p-value for the transient long-term sensitivity of IOD amplitude,\npositive IOD event intensity, and negative IOD event intensity. The results are as\nfollows:\n• IOD amplitude: In all high-CO simulations, the linear trends are statistically\n2\nsignificant (P<0.05). We accordingly added the result in the caption of Fig. 2 as “All\nlinear trends are statistically significant (P<0.05).”\n• Positive IOD event intensity: In all high-CO simulations, the linear trends are\n2\nsignificant (P<0.05) except for GISS-E2-R abrupt4x. GISS-E2-R abrupt4x shows a\nnegative linear trend, but is not statistically significant (P=0.06). When counting\nthe number of high-CO simulations with a decreasing trend, GISS-E2-R\n2\nabrupt4x is excluded. We accordingly added the result in the caption of\nSupplementary Fig. 9 as “All linear trends are statistically significant (P<0.05)\nexcept for GISS-E2-R abrupt4x. GISS-E2-R abrupt4x shows a negative linear trend,\nbut is not statistically significant (P=0.06). When counting the number of high-CO\n2\nsimulations with a decreasing trend, GISS-E2-R abrupt4x is excluded.”\n• Negative IOD event intensity: In all high-CO simulations, the linear trends are\n2\nsignificant (P<0.05) except for CCSM3 abrupt8x. CCSM3 abrupt8x shows a positive\nlinear trend, but is not statistically significant (P=0.06). We accordingly added the\nresult in the caption of Supplementary Fig. 10 as “All linear trends are statistically\nsignificant (P<0.05) except for CCSM3 abrupt8x. CCSM3 abrupt8x shows a positive\nlinear trend, but is not statistically significant (P=0.06).”\n(2) Inter-model consensus\nIn the main text of the original paper, we explicitly provided the inter-model\nconsensus ratio for the equilibrium response for IOD amplitude (L127), positive\nIOD event intensity (L129), negative IOD event intensity (L130), IOD amplitude\nsensitivity (L137). We found that inter-model consensus ratios for the transient\nresponse are not explicitly provided in the original paper (although they can be\nfound from the corresponding figure). We accordingly added them.\n• We added the inter-model consensus ratio for the long-term sensitivity of IOD\namplitude in L175. The added result is read as “In nearly all high-CO2 simulations\n(17 out of 18), the IOD amplitude shows a long-term decreasing trend with\nincreasing GMST. The exception is HadCM3L abrupt8x.”\n• We additionally performed analysis on the transient IOD event intensity changes\nfollowed by Reviewer #2’s comment (Supplementary Figs. 9 and 10). We\naccordingly added the inter-model consensus ratio for the long-term sensitivity of\npositive/negative IOD event intensity in L202. The added result is read as:\n“Consistent with the IOD amplitude change, the intensity of positive and negative\nIOD events also shows a long-term decreasing trend with increasing GMST (15 of\n18 high-CO2 simulations for positive IOD and 15 of 18 for negative IOD), but shows\nlarge fluctuations and widespread transient sensitivity to GMST (Supplementary\nFigs. 9 and 10).”\nvariability, transient IOD response, need to be further clarified, make them easier\nfor readers to understand.\nReply: Following your suggestion, we clarified the following terms:\n(1) IOD variability\nWe added a definition of IOD variability, where the term ‘IOD variability’ first\nappears in the manuscript in L54. The full sentence with the added definition reads\nas follows:\n“IOD variability, a fluctuation in climate variables induced by the IOD, is typically\nmeasured by the dipole mode index (DMI), which is defined as the difference\nbetween the western (50° E-70° E and 10° S-10° N) and eastern (90° E-110° E and\n10° S-0° N) SST anomalies of the tropical Indian Ocean.”\n(2) IOD amplitude and IOD event intensity\nWe added a clear explanation of ‘IOD amplitude’ and ‘IOD event intensity’ in L105.\nThe added explanation reads as follows:\n“We apply two different measures for the strength of the IOD variability. The main\nmetric is the IOD amplitude, defined as the standard deviation of the DMI. The\nother is IOD event intensity, defined as the peak DMI during the positive or negative\nIOD event (Methods). Both metrics are highly correlated by definition (i.e., if the\nIOD amplitude is high, the IOD event intensity is very likely to be high, and vice\nversa), but we mainly use the IOD amplitude to examine the change in the IOD\nvariability strength. The IOD amplitude is a fundamental metric that can show the\noverall response of the IOD variability to greenhouse gas warming. The simplicity\nof the IOD amplitude definition has advantages that facilitate the physical\ninterpretation of the results.”\n(3) Equilibrium and Transient IOD response\nWe added a new section in Methods, which explains the exact definition of\nequilibrium and transient IOD response. As we briefly explain the equilibrium and\ntransient response in the main text, we decided to place their exact definitions in\nMethods, considering the flow of the main text. Please see a newly added section\nentitled, ‘Definition of the equilibrium and transient response’ in L409. The\nsection is read as follows:\n“The equilibrium IOD response is defined as the response of IOD variability to an\nincreased equilibrium GMST level. The equilibrium response reflects the forced\nequilibrium response due to greenhouse gas warming. The transient response is\ndefined as the evolution of IOD variability in response to increasing GMST with\ntime. The transient response reflects both the forced response and internal\nvariability. Specifically, we interpret the long-term trend in IOD amplitude (or IOD\nevent intensity) as the forced response, and the fluctuation in IOD amplitude (or\nIOD event intensity) that deviates from the long-term forced trend as the internal\nvariability.”\n(4) Internal variability\nWe added the exact meaning of ‘internal variability’ in the newly added section,\nentitled ‘Definition of the equilibrium and transient response’ in L409. The\nexplanation in this part read as follows:\n“Specifically, we interpret the long-term trend in IOD amplitude (or IOD event\nintensity) as the forced response, and the fluctuation in IOD amplitude (or IOD\nevent intensity) that deviates from the long-term forced trend as the internal\nvariability.”\nWe also added the exact definition of the internal variability in the main text in\nL184. The added explanation read as follows:\n“Throughout the paper, we refer to the anomalous short-term fluctuation in the\nvariable that deviates from the long-term forced change as internal variability.”\nIn addition to this newly added exact definition of IOD internal variability, we refer\nto the general meaning of the internal variability in the main text, such as in L73\n(general meaning of internal variability). We believe that this newly added\nexplanation can make it easier to understand.\nAuthor’s Reply to the Reviewer #2\nReviewer #2 (Remarks to the Author):\nReview for submission “Decreased Indian Ocean Dipole variability under prolonged\ngreenhouse warming”\nThis study utilized LongrunMIP to address the forced long-term changes in the\nIndian Ocean Dipole (IOD) amplitude arising from natural variability. The manuscript\nis well-written, and the arguments are supported. Nevertheless, there are some\nconcerns that the authors should address in the revision.\nReply: We appreciate your recognition and valuable comments on our paper.\nThanks to your constructive comments, we have made a substantial improvement\nto the manuscript. The revised manuscript contains some major updates as follows:\n1. We evaluated the performance of the LongRunMIP model and selected models\nthat can reasonably simulate the IOD. We updated the analyses in the paper with\nthese selected models (9 out of 11 models). For the model selection process and\ncriteria, please see the newly added section Supplementary Discussion 1.\n2. We performed the physical mechanism analysis on how ENSO influences IOD.\nPlease see the newly added section Remote ENSO influence on the IOD change.\nThe reply to your comment is attached below. Please note that we will refer to the\ntext line and figure/table number based on the revised version of the manuscript.\n1. About the natural variability. Line 80: \"IOD amplitude in the twenty-first century\neither increases or decreases depending on the pathways of internal variability,\ndespite the same warming forcing.\" The term 'pathways' may be more closely\nassociated with forcing rather than natural variability. Please consider revising for\nclarity.\"\nReply: We agree that the term ‘pathways’ is not clear for this sentence. Following\nyour suggestion, we remove the term ‘pathways’ and clarified the sentence as\nfollows:\n“A large ensemble simulation using the Community Earth System Model shows that\ninternal variability alone can generate widespread long-term trends in future IOD\nvariability change; projected IOD amplitude in the twenty-first century either\nincreases or decreases depending on the ensemble member with different internal\nvariability, despite the same warming forcing.”\nThank you.\n2. About the simple IOD model. In paper 20, four terms are discussed: local feedback,\ndelayed feedback, ENSO impact, and stochastic forcing. In this study, the\nauthors only consider three terms, excluding the delayed feedback. Please provide\nmore discussions or explanations on this consideration. Line 381: “We consider the\nseasonality of λ, β, and σ. Thus, these three parameters vary seasonally. The typical\nvalue of the annual mean λ, which was estimated from observations and climate models in\nprevious studies, is negative.' It is assumed here that the model is fitted for each\nmonth with ENSO and DMI index in phase (lag-0). Whether you consider a lead-\nlag relationship? Also, for the annual mean model, does the equation still hold?\"\nReply: We would like to explain about the simple IOD model as follows:\n(1) Delayed feedback term\nAs you pointed out, we did not explicitly consider the delayed feedback term. There\nare several versions of the simple IOD model. For example, Malte et al. (2017) and\nour study consider three terms, local feedback (λT), remote ENSO impact (β×ENSO),\nand stochastic forcing (σξ). An et al. (2022) (which is written by our research group)\nadditionally considers the delayed feedback term (αT(t-τ)) as you mentioned. Both\nversions of the model are known to reproduce the DMI time series very well, and\nthe difference comes from the consideration of the delayed feedback process.\nTo explain the difference between the two versions of the model, we would like to\nexplain the physical meaning of αT(t-τ). As the wind forcing is exerted in the\nIndian Ocean, the response of the thermocline to the wind forcing is achieved by the\nocean waves. There are two waves involved in this process, the Kelvin wave and the\nRossby wave. The Kelvin wave is fast; therefore, the Kelvin wave response of the\nthermocline and the associated SST are almost simultaneous, and this process can\nbe included in λT. In contrast, the Rossby wave is three times slower than the Kelvin\nwave (time delay is typically 2 months), thus it is a delayed process. To explicitly\naccount for this delayed Rossby wave response, An et al. (2022) introduced a new\nterm, αT(t-τ).\nAs we explained, Malte et al. (2017) and our study did not include this term in the\nsimple IOD model, and only used λT to represent the feedback process occurring\nin the tropical Indian Ocean. The reason for this is that λT can implicitly account\nfor the delayed Rossby wave response.\nIf we fit the model parameter for the simple IOD model with αT(t-τ), the effect of\nthe delayed Rossby wave response would be explicitly represented by αT(t-τ). In\ncontrast, if we fit the model parameter for the simple IOD model without αT(t-τ),\nthe effect of delayed Rossby wave response cannot be explicitly represented by\nαT(t-τ). Instead, it would be implicitly included in the fitted seasonally varying λ. As\nan idealized example, we assume that the Rossby wave initiation activity peaks in\nOctober and the lag time is 2 months. For the model with αT(t-τ), this process\nwould be explicitly represented as αT(t-2) with α peaking in October. In contrast,\nfor the model without αT(t-τ), this process would be included in λ. Specifically, it\nwould increase the value of λ in December. Therefore, λT can implicitly account for\nthe Rossby wave response effect (at least partially).\nThe very high simulation skill of the simple IOD model demonstrated in our study\nand in Malte et al. (2017), shows that the exclusion of αT(t-τ) is a reasonable\nsimplification that doesn’t hinder the actual model performance. Moreover, the\nsimplicity of the model without αT(t-τ) also brings advantages that facilitate the\nphysical understanding of the model results, especially in determining whether the\nIOD amplitude change is caused by a local process or a remote process.\nTo clarify this, we added the explanation in the main text. Please see L474. The\nadded explanation is read as follows:\n“We note that there are several versions of the simple IOD model in the current\nliterature depending on the types of feedback processes included. The simple IOD\nmodel used in this study is the basic type model introduced in ref. 31, which\nconsiders the essential physical process of IOD in the simplest way. Hence, the\nmodel does not explicitly account for the time-delayed Rossby wave process as\ndone in ref. 20, and implicitly considers it by λT, which comprehensively represents\nlocal feedback processes in the tropical Indian Ocean. The simplicity of the model\nused in this study has advantages that facilitate the physical understanding of the\nmodel results, especially in determining whether the IOD amplitude change is\ncaused by a local process or a remote influence.”\n(2) Lag between ENSO and DMI\nThe simple IOD model can consider the lag between the Niño 3.4 index and DMI.\nAs you can read from the equation, the model relates the Niño 3.4 index (ENSO)\nwith the tendency of DMI (dT/dt), not with the DMI itself (T). Therefore, Niño 3.4\nindex and DMI are not in-phase and would show an out-of-phase relationship.\nWe mathematically demonstrate that the model can consider the lag between Niño\n3.4 and DMI. For simplicity, we consider only the ENSO forcing in the simple IOD\nmodel, dT/dt=β×ENSO(t). We set the Niño 3.4 index to vary sinusoidally with\namplitude A and frequency w, ENSO(t)=Acos(wt). The analytical solution of T is\nT(t)=(βA/w)sin(wt). As you can see, ENSO and T show a phase lag of π/2. This is a\nhighly idealized example showing that the basic form of the model (which relates\nthe Niño 3.4 index with the tendency of the DMI) can consider a lag between the\nNiño 3.4 index and the DMI. For the full model with a real Niño 3.4 index, the lag\nwould depend on the λ, β, and characteristics of the Niño 3.4 index time series.\nThe model parameter fitting is the process of finding the optimal set of parameters\nthat can best explain the given time series under the constraints of the model\nequation. Since the model itself can naturally consider the lag relationship, it is not\nnecessary to modify the equation and the associated parameter fitting process to\nconsider the lag. It is okay to just fit the model parameters to the original equation\nwithout modification.\nTo clarify that the model can consider the lag, we added an explanation in L470.\nThe added explanation is read as follows:\n“The model can consider the lag between ENSO and T because the model equation\nrelates the ENSO with the tendency of T, not with the T itself”.\n(3) Annual mean model\nFirst of all, we would like to clarify the meaning of the ‘annual mean’ in this sentence\n“The typical value of the annual mean λ…”. The ‘annual mean’ here, means the\nannual mean value of the seasonally fitted parameter, not the annually fitted\nparameter. We would like to note that the two are different. The former is obtained\nby (1) fitting the parameter for each season (12 calendar month value for each\nparameter), and (2) taking the annual mean value of the parameter. The latter is\nobtained by fitting the parameters for the entire period without considering\nseasonality (single value for each parameter).\nIn this paper, we only used and discussed the seasonally-varying model. We also\nonly cite the reference that uses the seasonally-varying model, including the\nsentence you mentioned (as far as we know, all literature using the simple IOD\nmodel uses the seasonally-varying model to consider the seasonality of the physical\nprocess). The annual model (i.e., model with seasonally-fixed parameters) is neither\nused nor mentioned. Thus, the discussion of the annual mean model is out of the\nsubject of this study, because we neither use it, nor explain about it, nor cite the\nrelated material.\nWe assume that your concern stems from the misunderstanding that we mentioned\nthe annual model in this sentence, whereas our intention was to explain the annual\nmean value of the seasonally-fitted parameter. We found the sentence has a dual\nmeaning and can cause confusion to readers. We would like to apologize for the\nunclarity and any confusion you may have experienced.\nTo explain the exact meaning of this sentence, background knowledge of the non-\nautonomous linear ordinary differential equation needs to be provided. λ\ncharacterizes the linear stability of T. For the case where λ is constant (i.e.,\nseasonally-fixed annual mean model, which was not used in this study), if λ is\nnegative (positive), T decays (grows). For the case where λ is seasonally-varying (i.e.,\nthe seasonal model used in this study), the annual mean value of the seasonally\nvarying parameter, characterizes the stability of T. If the annual mean value of\nseasonally-varying λ is negative (positive), T decays (grows). The first case is\nstraightforward, can be understood by a simple eigenvalue analysis, and does not\nrequire additional explanation. However, the second case (the model used and\nexplained in this study) is not straightforward and requires heavy mathematical\ncalculations to prove. The mathematical proof of this argument can be found in\nKim and An (2021) (Journal of Climate) (https://doi.org/10.1175/JCLI-D-20-0495.1)\n(although it is for ENSO and the model is two-dimensional, it presents an essentially\nmathematically identical system).\nAfter careful consideration, we decided to remove the sentence, “The typical value\nof the annual mean λ, which was estimated from the observations and climate\nmodels in the previous studies, is negative20,31,32. Hence, without external ENSO and\nstochastic forcing, the simple IOD model decays because of a negative growth rate.”.\nBecause, it requires a heavy and complex explanation of the mathematical\nbackground to properly deliver the information to readers, and it is not such an\nimportant discussion that affects understanding results in our paper.\nThank you.\n3. About the seasonal changes. Paper 19 addressed that the early IOD is more\nfrequently observed under warming conditions. In this study, Line 124 mentions,\n'The seasonal variability of IOD decreases in almost all seasons (Supplementary Fig.\n5).' Moreover, Figure S5 illustrates that HadCM3L and a few other models peak\nearly in JJA. It would be beneficial to delve into a more detailed discussion of this\nchange. Lines 135-146: You might consider separating the discussions on seasonal\nchanges and spatial patterns into two paragraphs for clarity.\nReply: We appreciate your insightful suggestion. As you noticed, some of the\nmodels including HadCM3L and CCSM3II, also show a shift of the seasonal peak\nto an earlier season, which is consistent with the findings from Sun et al. (2022)\n(although the warming scenario is different). In the main text, we added a discussion\nof this seasonal shift in L148. The added discussion read as follows:\n“This alteration in the seasonality also occurs in other models including HadCM3L\nand CCSM3II. HadCM3L simulations show a transition in seasonal peak from\nOctober (control and abrupt2x) to July (abrupt4x and abrupt8x). CCSM3II\nsimulations show a gradual shift of seasonal peak from September (control) to\nAugust (abrupt700ppm) to June (abrupt1400ppm). Such a pronounced increase in\nIOD variability in boreal summer is consistent with the previous study19 shows an\nincreased occurrence of early positive IOD events by global warming under the\nIPCC’s high-emissions scenario (although the warming scenario is different).”\nWe also separated the paragraph. Thank you.\n4. About the relationship between changes in IOD and ENSO. Paper 28 addressed\nlong-term ENSO changes based on LongrunMIP. It would be valuable to explore\nwhether models agree on changes in both IOD and ENSO and if the dominating\nmechanisms are similar. Authors should engage in discussions on this aspect.\nReply: We performed a series of additional experiments and analyses on the role\nof ENSO on the IOD. First, we performed the additional simple IOD model\nexperiment to separate the contribution of the ENSO amplitude change to the IOD\namplitude change (previously we quantified the contribution of the combined effect\nof ENSO variability and sensitivity of IOD to ENSO). Please see the newly added\nexperiment result in L245. Second, we performed an analysis of the physical\nmechanism of the remote ENSO forcing on the IOD amplitude change. Please see\nthe newly added section, “Remote ENSO influence on the IOD change” in L260.\nThese two analyses deepen the understanding of the relationship between ENSO\nand IOD.\nThe ENSO amplitude decreases in all of the high-CO simulations, consistent with\n2\nthe paper 28 you mention (Supplementary Fig. 15a and Fig. R1). The additional\nexperiment shows that the forced decrease in ENSO amplitude has a very small\ncontribution to the forced decrease in IOD amplitude. Therefore, the forced changes\nin IOD and ENSO amplitude are independently arising from greenhouse gas\nwarming. The forced change in IOD amplitude is attributed to the change in the\ntropical Indian Ocean, and the forced change in ENSO amplitude is attributed to\nthe change in the tropical Pacific Ocean (paper 28). Both ocean basins show similar\nfeatures, positive IOD-like change and El Niño-like change, respectively. The inter-\nsimulation correlation between forced ENSO and IOD amplitude change is 0.56 (Fig.\nR1). The correlation shows that there is some model agreement on the changes,\nbut does not imply causality between them.\nFig. R1. Relationship between change in ENSO amplitude and IOD amplitude in\nthe equilibrium period. The change is a percentage change from the control\nsimulation value (colored dots). Each colored dot represents the results for each\nhigh-CO2 simulation. The linear regression line is shown as a black line. The Pearson\ncorrelation coefficient is 0.56.\nWe have carefully considered including this part in the discussion section, in\naddition to the newly added analyses (L245 and “Remote ENSO influence on the\nIOD change”), but we concluded that it is somewhat out of the scope of the paper\nwhich focuses on the IOD. We feel that the analysis of the ENSO itself and the\ntropical Pacific Ocean should be done consistently with the IOD in order to\nconfidently argue it, although the paper 28 provides a very relevant reference. We\nbelieve that such a detailed analysis of ENSO and the tropical Indian Ocean is\nbeyond the scope of this paper. Furthermore, the relationship between ENSO and\nIOD is already elucidated in the newly added analyses. Therefore, we decided not\nto include this part in this paper. Thank you for your suggestion.\n5. About the definition of IOD amplitude. Line 102: 'We measure the IOD amplitude\nas the standard deviation of the DMI.' Typically, the community uses the peak\nseason (SON) DMI to gauge IOD changes. Here, monthly DMI is employed; the\nauthors may wish to explain this choice.\nReply: As you pointed out, the IOD change can be examined with the peak season\n(SON) DMI. In this study, we use two different metrics to examine IOD change. One\nis “IOD amplitude” which is defined as the standard deviation of DMI. The other\none is “IOD event intensity” which is defined as the highest (or lowest) DMI level\nduring the IOD event. Here, the positive (negative) IOD event is defined as a three-\nmonth running mean of DMI exceeding (below) the one standard deviation of the\ntime series for at least three consecutive months. Therefore, the “IOD event intensity”\nis close to a typical metric for the IOD change (a peak season DMI).\nThese two metrics are highly correlated by their definition (i.e., if the IOD event\nintensity is high, the IOD amplitude is likely to be high, and vice versa). In this study,\nwe analyze both metrics, but mainly focus on the “IOD amplitude”.\nThere are several reasons for this. Most of all, the IOD amplitude is a more\nfundamental metric than the IOD event intensity. The IOD amplitude measures the\noverall variability of a warm-cold cycle of DMI, considering the entire SST cycle of\nthe tropical Indian Ocean. Meanwhile, the IOD event intensity measures peak DMI,\nemphasizing the maximum SST variability. Therefore, the IOD amplitude is a more\nfundamental metric than the IOD event intensity, which can show more overall\nresponse of tropical Indian Ocean SST anomaly.\nThis facilitates the physical interpretation of the IOD change. We utilize the simple\nIOD model to disentangle the physical source of the IOD change. Here, the main\nvariable of this model is DMI (T) (including this model, theory and model of IOD is\ntypically based on DMI). Therefore, IOD amplitude, which is a standard deviation of\nDMI, can be directly linked with the simple IOD model, given by ((<T2>-<T>2)1/2)\nwhere <> denotes ensemble mean (even an analytical solution of DMI standard\ndeviation can be directly retrieved under very idealized parameter settings).\nMeanwhile, the IOD event intensity cannot be directly expressed with T as done for\nthe IOD amplitude due to the complexity of its definition. Also, using IOD amplitude\nis easier to directly attribute which region (eastern or western tropical Indian Ocean)\ncontributed to the IOD change by simply measuring the standard deviation of SST\nof these two regions.\nIn summary, these two metrics are highly correlated and both show the response\nof IOD to greenhouse gas warming. We analyzed both metrics, but more focus on\nthe “IOD amplitude” throughout the paper by its usefulness for the physical\ninterpretation.\nTo clarify this point we made two changes in the manuscript:\n(1) We added the additional explanation of the IOD amplitude definition. Please\nsee L105. The added explanation is read as follows:\n“We apply two different measures for the strength of the IOD variability. The main\nmetric is the IOD amplitude, defined as the standard deviation of the DMI. The\nother is IOD event intensity, defined as the peak DMI during the positive or negative\nIOD event (Methods). Both metrics are highly correlated by definition (i.e., if the\nIOD amplitude is high, the IOD event intensity is very likely to be high, and vice\nversa), but we mainly use the IOD amplitude to examine the change in the IOD\nvariability strength. The IOD amplitude is a fundamental metric that can show the\noverall response of the IOD variability to greenhouse gas warming. The simplicity\nof the IOD amplitude definition has advantages that facilitate the physical\ninterpretation of the results.”\n(2) We additionally performed analysis on the transient IOD event intensity changes\n(previously there was only analysis for the equilibrium IOD event intensity change\nwhich was shown in Supplementary Fig. 3). Please see newly added\nSupplementary Figs. 9 and 10. As you can see, the change in positive and negative\nIOD intensity is consistent with the change in IOD amplitude. The majority of high-\nCO simulation shows a decreasing trend in IOD event intensity against GMST (16\n2\nout of 18 for the positive IOD event and 15 out of 18 for the negative IOD event)\nand shows large fluctuation. We also accordingly added an explanation in the main\ntext. Please see L202. The added explanation is read as follows:\n“Consistent with the IOD amplitude change, the intensities of both positive and\nnegative IOD events also show long-term decreasing trends with increasing GMST\n(15 of 18 high-CO2 simulations for positive IOD and 15 of 18 for negative IOD),\nbut show large fluctuations and widespread transient sensitivity to GMST\n(Supplementary Figs. 9 and 10).”\n6. About transient response. Line 95: For the 1pct2x and 1pct4x runs, 'After the\nCO2 increase, the forcing is stabilized for typically longer than 1,000 years.' While\nthe first 200 years may be considered the transient period, the runs predominantly\nreflect the equilibrated response after the forcing stabilizes. Line 148: 'We analyze\nthe transient response of the IOD amplitude to the CO2 forcing. We calculate the\n100-year moving IOD amplitude for the entire simulation period.' It is noted that\nthe entire period (1000 years or more) covers both the transient and equilibrated\nperiods. The authors may need to consider taking the first 100-200 years as the\ntransient period when the forcing is changing. After the forcing stabilizes, it may\nnot be considered part of the transient period.\nReply: We carefully considered your suggestion and would like to explain in great\ndetail why we use the entire simulation period to examine the transient response.\nAs you can see in Supplementary Fig. 2, the GMST, which is the representative\nindicator of planetary warming level, increases over time due to the CO forcing.\nthere is a common GMST evolution pattern in the high-CO simulations. It increases\n2\nrapidly in the early period, and then gradually adjusts to the equilibrium level. In\nthe final period of the simulation, GMST increases only slightly and reaches the\n(quasi-) equilibrium level.\nFor the equilibrium response analysis, we defined the last 500 years of the\nsimulation as the “equilibrium period’. We examined the equilibrium response of\nIOD by performing an inter-experiment comparison for the IOD amplitude (and\nevent intensity) in the equilibrium period.\nIn contrast, for the transient response analysis, we used the entire simulation period\nand did not specifically define the ‘transient period’ (we do not use the term\n‘transient period’ throughout the paper). Here, transient response refers to the\ntemporal evolution of IOD variability (the strength of which is quantified as IOD\namplitude and IOD event intensity) in response to greenhouse gas warming.\nTherefore, by definition, the entire simulation period shows the transient response\nof IOD to greenhouse gas warming. The continuous evolution of IOD variability\nover the whole simulation period shows the transient response of IOD during the\nperiod where both large (e.g., 0-200 years) and small (e.g., 2000-3000 years)\nincreases in GMST occur. We would like to note that ‘the time evolution of the IOD\nvariability during the equilibrium period’ and ‘the response of IOD variability to\nincreasing equilibrium GMST’ are different. The former is the transient response of\nIOD by definition, and the latter is the equilibrium response of IOD that we analyzed\nin the first part of the paper.\nAs you suggested, there is an alternative approach to understanding the transient\nresponse for LongRunMIP simulations. We can define the transient period as the\nfirst 100-200 years (or the first 0-200 years) as you suggested. This type of approach\nconsiders the IOD amplitude change when the GMST (or forcing) change is large.\nThe difference between these two approaches lies fundamentally in how we define\n‘transient’ and related terms. The approach of your suggestion implicitly defines a\n‘transient period’ as the period during which the significant change in GMST (or\nCO forcing) occurs. Therefore, it focuses on the change in IOD variability during\n2\nthe period when a large increase in GMST occurs. In contrast, our approach\nconsiders the ‘transient response’ as the temporal evolution itself. Therefore, our\napproach covers the change in IOD variability for the entire period, including the\nperiod when GMST increases small.\nAs you can see in Fig. 2, we analyzed the transient IOD response on the basis of\nGMST level, not time. This GMST-based analysis can naturally consider the\ndifference in GMST level depending on the time period (the concern you raised).\nSince the changes in IOD-related variables are analyzed on the basis of GMST\nchange, the transient IOD response during the period of large GMST increase\nperiod (e.g., 0-200 years), will make a large contribution to the result, while the\ntransient IOD response during the period of small GMST increase period (e.g., 2000-\n3000 years), will make a relatively small contribution. Therefore, we can expect that,\nthere would be no significant difference in the main results even if we perform the\nanalysis for the selected transient period of large GMST increase.\nTo examine this, we calculate the IOD amplitude for the first 0-500 years for all\nsimulations. The period 0-500 years is roughly a median value of the period of\nlarge GMST increase (e.g., GFDLCM3 1pct2x shows large GMST increase roughly up\nto year 1000 and MIROC 3.2 1pct2x shows up to year 200) (please note that there\nis a large difference in the GMST increase rate between simulations). Please see Fig.\nR2. As you can see, there is practically no difference in the main result when we\nchoose the first 0-500 years for the transient response. We reach the same\nconclusion that the decreasing forced long-term trend in IOD amplitude occurs in\nthe vast majority of the simulations.\nFig. R2. Same as Fig. 2a in the main paper, but for the transient period (years 0-\n500).\nAs we demonstrated, both approaches to ‘transient response’ analysis are\nreasonable and bring no practical difference to the main conclusion of the study.\nHowever, our approach has several advantages. First of all, we can avoid the\nproblem of period selection dependence. If we specifically choose the transient\nperiod, we have to specifically set the period criteria (e.g., whether we choose the\n50-150 years or the 10-200 years for the transient period). The result would depend\non how we choose the transient period (although the difference may be small).\nAlso, the transient period approach can potentially hinder the fair comparison due\nto the large inter-model spread in the rate of GMST increase. For example, MIROC32\n1pct2x shows fast adjustment to equilibrium GMST level (proper transient period\nwould be 0-200 years), while the GFDLCM3 1pct2x shows very slow adjustment to\nequilibrium level (proper transient period would be 0-1000 years). Therefore, it is\nnot fair to compare them in the same line by applying uniform transient period\ncriteria. We also considered setting the ‘transient period’ criteria based on the rate\nof GMST increase, but this would be more complex and would hinder the\ninterpretation of the results (somehow, we have to set specific criteria for the\nwarming rate as well). In addition, our approach provides a complete description\nof the IOD amplitude evolution spanning the entire simulation period (the evolution\nduring the transient period is part of our result using the entire period).\nTherefore, our approach, which considers the entire simulation period, is the most\nnatural approach to analyzing the time evolution of IOD in response to greenhouse\ngas warming. It allows a fair comparison between simulations with different rates\nof GMST increase and facilitates the understanding of the results in as consistent a\nframework as possible.\nIn addition, we would like to note that it would make more sense to use the GMST\nas the basis for analysis rather than the atmospheric CO concentration level.\n2\nplanetary warming and alteration of the hydrological cycle. Therefore, this warming\nand associated changes are key physical effects of the CO forcing. As there is a\n2\nlarge time lag between CO and GMST due to the inertia effect, it would be more\n2\nphysically reasonable to use GMST, which brings practical impact to the climate\nsystem, as the basis of the analysis.\nFor these reasons, after serious and careful consideration, we decided to retain the\noriginal analysis using the entire simulation period. To clarify our approach, we\nadded a new section in Methods, which explains the exact definition of transient\nIOD response. Please see a newly added section entitled, ‘Definition of the\nequilibrium and transient response’ in L409. The section is read as follows:\n“The equilibrium IOD response is defined as the response of IOD variability to an\nincreased equilibrium GMST level. The equilibrium response reflects the forced\nequilibrium response due to greenhouse gas warming. The transient response is\ndefined as the evolution of IOD variability in response to increasing GMST with\ntime. The transient response reflects both the forced response and internal\nvariability. Specifically, we interpret the long-term trend in IOD amplitude (or IOD\nevent intensity) as the forced response, and the fluctuation in IOD amplitude (or\nIOD event intensity) that deviates from the long-term forced trend as the internal\nvariability.”\nThank you for the suggestion.\n7. About the period of analysis. Line 107: 'last 500 years.' Are the entire 500 years\nused for detrending through quadratic fitting? Line 312: 'For the LongRunMIP, we\napplied quadratic detrending to the monthly SST to remove the forced long-term\nchange (i.e., secular mode).' Could the authors clarify whether they used a 100-year\nor 500-year period for this detrending process?\nReply: We apologize for the inconvenience you might have experienced due to the\nunclarity of the description. The detrending is performed for the entire period of\nthe simulation (see Supplementary Table 1 for the length of each simulation).\nTherefore, we did not use either the 100-year or 500-year period for the detrending\nas you questioned. This detrending removes the long-term warming trend induced\nby CO forcing (please see Supplementary Fig. 2 as an example of the long-term\n2\nwarming trend).\nAs a result, for each simulation, we obtain the eastern/western SST anomaly and\nDMI time series, the length of which is equal to the length of the entire simulation.\nFor the equilibrium response analysis, we select the last 500 years of the\nDMI time series and calculate its standard deviation (results shown in Fig. 1a). For the\ntransient response analysis, we calculate the 100-year moving DMI standard\ndeviation for the entire period of DMI time series (results shown in Fig. 2a).\nTo clarify this process, we added explanations in L388 and L394. The added\nexplanations are read as follows:\n“The detrending is performed for the entire period of the simulation (see\nSupplementary Table 1 for the length of each simulation).”\n“As a result, for each simulation, we obtain the DMI and Niño 3.4 index time series,\nthe length of which is equal to the length of the entire simulation. Using these DMI\nand Niño 3.4 index time series, we perform analyses including the calculation of\nIOD amplitude and event intensity.”\nThank you.\n8. Figure S13. The upper panel x-y axis should be swapped to be consistent with\nthe lower panel.\nReply: We appreciate your suggestion and carefully considered the swapping the\nx-y axis in the upper panel.\nUsually for a scatter plot, if the variable ‘A’ potentially causes variable ‘B’, ‘A’ is\nplotted on the x-axis and ‘B’ is plotted on the y-axis. In our case, the upper panel\nshows the relationship between change in GMST and change in the model\nparameter (λ, β, ENSO amplitude, σ). If there is no potential causal relationship\nbetween these two variables, it would be reasonable to swap the x-y axis and make\nit consistent with the lower panel, following your suggestion. However, there is a\nvery strong potential causal relationship between them. Physically, the change in\nmodel parameters would be caused by the change in GMST (i.e., forced warming\nby the CO increase). Therefore, considering the strong causal relationship between\n2\nthem, it is more reasonable to plot the GMST change on the x-axis and the model\nparameter change on the y-axis.\nThe same plot style is applied to the lower panel. The lower panel shows the\nrelationship between the change in the model parameter (λ, β, ENSO amplitude, σ)\nand the change in the DMI standard deviation. As we can read from the simple\nIOD model equation, the change in the model parameter directly causes changes\nin the DMI and its standard deviation. It is therefore reasonable to plot the GMST\nchange on the x-axis and the model parameter change on the y-axis.\nIn conclusion, the current upper and lower panel shows a causal relationship\nbetween two variables in a consistent way (cause variable on the x-axis and effect\nvariable on the y-axis). We believe maintaining the plotting style can deliver the\nresults of this plot to readers emphasizing the physical relationship between them.\nTo clarify the causal relationship between the two variables, we added the following\nexplanations in the figure legend:\n“The GMST change by CO forcing would cause the change in the model\n2\nparameters and ENSO, contributing to the IOD amplitude change.”\nThank you.\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have responded to all my concerns and revised their manuscript accordingly. I have no\nfurther concerns or comments.\nReviewer #2 (Remarks to the Author):\nThe authors have effecfively addressed my concerns by providing addifional evidence and enhancing the\nclarity of the study. At this point, I have no further comments to add."}]

Decision letter