Intensity changes of Indian Ocean dipole mode in a carbon dioxide removal scenario

The Indian Ocean Dipole/Zonal mode (IOD) is an interannual phenomenon over the tropical Indian Ocean, causing a pronounced impact worldwide. Here, we investigate the mechanism of the change in IOD characteristics in a CO2 removal simulation for an earth system model (ESM). As the CO2 concentration increases, the intensity of IOD tends to increase, but at high CO2 concentrations, further increases decrease the IOD intensity. The minimum IOD amplitude was recorded during the early decrease in CO2. First, we developed a conceptual model for IOD that is composed of local air-sea coupled feedback, delayed ocean dynamics, El Niño impact, and noise forcing. Then, by adopting ESM results into this simple IOD model, we revealed that the local air–sea coupled feedback is a major factor for changing IOD amplitude, while El Niño does not exert a change in IOD amplitude. The local air–sea coupled feedback including thermocline feedback, wind-evaporation feedback, and Ekman feedback is strongly modified by the air–sea coupling strength during progression of a global warming. Consequently, under the higher CO2 concentrations, IOD amplitude is reduced due to the weakening of air-sea coupling over tropical Indian Ocean.


INTRODUCTION
The Indian Ocean Dipole/Zonal mode (IOD/IOZM) is a well-known phenomenon with an interannual variability over the Indian Ocean (IO) [1][2][3][4] , causing a pronounced impact over East Africa, western Indonesia, Australia, India, and other regions [5][6][7][8][9][10][11][12][13][14] . The positive IOD phase is characterized by cold sea surface temperature (SST) anomalies in the southeastern tropical IO and warm SST anomalies in the western tropical IO (Fig. 1a), and its negative phase is a mirror image. It tends to develop in boreal summer, peaks during fall, and decays rapidly in winter 1 , implying the seasonal-phase locking of IOD due to strong seasonal changes in the monsoon flow over the IO 4 . It is also known that El Niño-Southern Oscillation (ENSO; Fig. 1c), that is an interannual tropical Pacific SST anomaly changing phenomenon (Fig. 1c), can trigger IOD, which in turn influences the development and/or decay of ENSO 11,[15][16][17][18][19] .
IOD grows through wind-thermocline-SST (WTS) feedback 1,20 , wind-evaporation-SST (WES) feedback 21 , and cloud-radiation-SST feedback 22,23 , suggesting the importance of internal feedback processes. External forcing has also been proposed for IOD formation, such as ENSO 4,23,24 , the southern annular mode 25 , the Indo-Pacific warm pool 3,26 , and the Indonesian Throughflow 27 . Among them, the ENSO is likely to be the most influential. Most IOD events are accompanied by ENSO 11,[15][16][17] . That is, during a developing El Niño, a zonal shift in anomalous Walker cells due to central-to-eastern Pacific warming can trigger IOD events by modifying the surface winds over the tropical eastern IO 4,23 . IOD events are more predictable when they occur together with ENSO events because of the close connection between IOD and ENSO events than when they occur independently 28,29 . Therefore, the minimal but sufficient governing mechanisms of IOD include internal feedback processes, external driving factors (i.e., ENSO), and stochastic forcing such as weather noise 30 .
In recent decades, an increasing trend of positive IOD events has been observed [31][32][33] . Furthermore, a near-future IOD is projected to be stronger during boreal summer relative the current period 34 . A study 35 mentioned that the projected mean climate warming in boreal fall, displaying a positive IOD-like climate change, could lead to stronger easterly winds just south of the equator and a shoaling equatorial thermocline. Such a change in the mean climate may lead to stronger positive Bjerknes feedback that would lead to stronger IOD events 36 . However, the same study also indicated that the overall frequency of IOD events under a warmer Earth at the end of the 21st century would not change, while the difference in amplitude between positive and negative IOD events would decrease. These results suggest that the relationships between changes in the mean state over tropical IO associated with global warming and the projected changes in IOD have not yet been fully addressed. Therefore, we quantitatively diagnose the governing mechanisms of the changes in IOD characteristics in a changing CO 2 pathway (i.e., carbon dioxide removal scenario), and qualitatively investigate the possible role of changes in the mean state.

IOD and ENSO changes with changes in CO 2
To explore the changes in IOD as the CO 2 concentration changes, we computed the wavelet power spectrum of Indian Ocean dipole mode index (DMI 9 : difference between SST anomalies in IO west (50°E -70°E and 10°S-10°N) and IO east (90°E-110°E and 10°S-0°S)) from the 28-ensemble of CDR-reversibility experiments in CESM1.2 (see the "Methods" section), and then took their average. Therefore, IOD modulation by natural variabilities was likely removed, but the signal due to CO 2 forcing remained. As shown in Fig. 1b, during the rampup of the CO 2 concentration, the variability in IOD with around 4-year period increased with CO 2 , but as the CO 2 concentration increased further (more than double the initial concentration), the amplitude decreased. During the ramp-down of the CO 2 concentration, the amplitude of IOD recorded its minimum for the first 30 years throughout the entire integration period, and then started to slowly increase, but the overall amplitude of IOD was low. The frequency of the IOD became slightly faster during the ramp-down period. During the restoration period, the amplitude of IOD almost recovered, and its frequency was similar to that during the ramp-up period. During the restoring period, the amplitude of the IOD is very slowly recovered.
Similarly, we analyzed the changes in ENSO. During the rampup period, the amplitude of ENSO showed a weak increasing trend (Fig. 1d). Interestingly, the IOD amplitude decreased during the 100-140 years, but the ENSO amplitude increased. During the early ramp-down period, the amplitude of ENSO slightly decreased, and its period decreased. Later, the amplitude and frequency almost recovered. During the restoring period, the amplitude of ENSO was slightly lower than that during the rampup period. As seen in the wavelet power spectrums, the change in IOD characteristics is more distinct than that in ENSO. However, somewhat similar changing characteristics were observed between IOD and ENSO, yet a connection between them could not be justified from these results.

Reproduction of IOD using a simple IOD model and dynamics
To investigate the cause of the IOD changes, we used the following simple IOD model (see the "Methods" section).
where T is DMI; δ is a delayed time associated with ocean dynamic adjustment; E(t) is ENSO forcing (i.e., Nino3.4 index); and ξ t is a Gaussian white noise with a zero mean and unitary standard deviation. Physically, λ n , α n , β n , and σ n represent parameters for comprehensive local feedback, delayed feedback through the delayed oceanic wave response to the wind (see the "Methods" section), ENSO-driven impact, and amplitude of stochastic forcing, respectively.
The time series of the 30-year moving variance of the original and reproduced DMI indices for the PD experiment are presented in Fig. 2a. The DMI shows a multi-decadal amplitude modulation in both the original and reproduced DMI indices. The correlation between two time series is 0.66, which is statistically significant at the 99% confidence level. Therefore, this simple IOD model can be applied for IOD diagnosis. Note that some discrepancy between original DMI and reproduced DMI must be due to missing physics in IOD model. One of them may be nonlinear process, which is responsible for a skewed IOD.
The parameter λ n shows a dominant annual cycle (Fig. 2b), of which the maximum appears in July. As λ n is positive from July to September and then becomes negative in October, the peak IOD is expected to occur in September, as observed 1,2 . α n does not show any annual/semi-annual cycles and is negative throughout the year except in April, indicating that the main role of this term is the delayed negative feedback. The annual cycle of background conditions does not influence this delayed negative feedback (Fig. 2c). β n is also dominated by an annual cycle and also shows a weak semi-annual cycle (Fig. 2d). A distinct maximum is observed in May, while the values of β n in other calendar months are relatively low, indicating that the ENSO has a triggering effect on IOD rather than providing a driving force. σ n also shows a dominant annual cycle with the maximum in June, and it thus acts on the growth of IOD during the boreal summer, similar to λ n . The life cycle of IOD inferred from these parameters is likely as follows: ENSO and/or weather noise initiates IOD during the late spring; during the summer, IOD grows mainly by the local air-sea coupled feedback and stochastic forcing but is dampened by the delayed negative feedback, maturing during the early fall as λ n becomes negative, and owing to the weak influence of ENSO, IOD quickly dies out after the peak.

Mechanism of IOD change with changes in CO 2
Here, to understand why IOD changes as CO 2 changes, each ensemble result from the CDR-reversibility experiment was adopted for the simple IOD model, and a total of 28 sets of parameters were computed for the simple IOD model. Using each parameter set, we conducted an ensemble simulation of the simple IOD model, and then, the ensemble mean was taken. All results are the ensemble mean from the simple IOD model simulation. As shown in Fig. 3a, the 30-year moving variance of DMI obtained from the simple IOD model follows that from CESM1.2, but it shows an overall underestimation bias. Before computing DMI from CESM1.2 simulation, trends in SST have been removed. Similar to the wavelet spectrum in Fig. 1b, an increasing tendency is presented in the IOD amplitude during the ramp-up, but it starts to decrease at around 90 model years. The minimum occurs during the early ramp-down period and then starts to recover. During the restoration period, the amplitude modulation of the DMI was quite weak. Overall, there are three interesting periods: the period with the largest amplitude during ramp-up (around 50 model years; LA), the period with the smallest amplitude during the early ramp-down (around 150 model years; SA), and the period with a recovered modest amplitude during the restoration period (around 350 model years; MA). As the IOD amplitude during the MA period is similar to that during the PD period, we mainly focus on the two distinct periods, LA and SA. Each period is defined as 30 years centered at 50 and 150 model years.
During the LA period, λ n is positive from May to September, indicating a much longer growing season than that for SA. α n was positive in April and June. β n is just slightly higher than that in SA except during August, indicating that the role of ENSO may not be important in influencing the IOD amplitude change. σ n is also slightly larger than that in SA except during February and March, but the difference in σ n between SA and LA is small compared to that in λ n . Therefore, the enhanced IOD amplitude during the LA period is likely due to an enhanced internal feedback (i.e., large λ n ).
During the SA period, the peak of λ n shifted to an earlier calendar month, and a positive λ n was observed for only two summer months (June and July). Thus, the shorter growing season and stronger damping result in a low amplitude for IOD, which are related to changes in the background state of the tropical IO (see next section). It has been mentioned that a new type of positive IOD is currently emerging, which develops and terminates earlier than a conventional IOD 35 . This type of IOD may be observed in our simulation, in the higher CO 2 concentration regime (i.e., SA period).
To further confirm the role of each parameter, we performed 1000-year integrations of the simple IOD model with the parameters for each LA and SA period, as exactly shown in Fig. 3b-e, which are referred to as the control simulations for LA and SA, respectively. As a sensitivity experiment, for example, λ n s for LA and SA were averaged, and then the original λ n in each control simulation was replaced by the averaged λ n . The same sensitivity experiment was performed for other parameters. Here, the ENSO forcing (that is, E(t) in Eq. (1)) was randomly selected from the 30-year time series of the monthly NINO3.4 index for each LA and SA period, but the seasonality of E(t) matched with that of the original Nino-3.4 timeseries. Table 1 shows the variance of DMI obtained from the control and sensitivity experiments. The smaller the difference in variance between the LA and SA sensitivity experiments, the stronger the impact of the corresponding parameter. Therefore, the most important factor for controlling the IOD amplitude is λ n , and the second most important factor is α n , while the influences of β n and σ n are negligible. These results further confirmed the importance of internal feedback for modifying the IOD amplitude in changing CO 2 forcing.

Mean state change and its impact on IOD
The previous subsection shows that comprehensive local feedback (i.e., λ n ) is the most influential factor in modifying the IOD intensity when the CO 2 concentration is changing. Here, we show the mean atmosphere and ocean states of the IO to investigate how the changes in the mean state modify IOD intensity. Here, we only conducted a qualitative approach because the overall quantitative analysis was conducted in the previous section, and focused on the analysis of the LA and SA periods because differences among PD, MA, and LA were not substantial. Figure 4a, b show the differences in the boreal summer mean SST and surface winds over the IO between LA (SA) and PD. The IO summer-mean SST during both LA and SA was warmer than that during PD. The basin-wide warming, removed in the SST pattern of Fig. 4a, b, during SA (3.65°C warmer than PD) was higher than that during LA (1.03°C warmer than PD). Furthermore, owing to the differential warming rate over the eastern and western IO, which is a common feature of the projected IO SST shown in most CGCMs 35,37 , the east-west contrast in the mean SST became stronger during the SA than that during the LA (Fig. 4a, b). This is presumably due to the ocean dynamic thermostat effect, which is stronger over the eastern IO. The eastern IO is the destination of the seasonal equatorial upwelling Kelvin waves and a shallow thermocline induced by strengthened equatorial easterlies 35 . Consequently, the easterly trade winds during SA were stronger than those during LA because of stronger global warming (Fig. 4a, b). Changes in trade winds modified the wind-evaporation-SST (WES), wind-upwelling-SST (WUS), and wind-thermocline-SST (WTS) feedbacks. Because WES feedback refers to a chain reaction referring to "wind-evaporative cooling-SST," its feedback strength is related to the extent to which the mean wind change modifies the anomalous latent heat flux associated with anomalous winds. A bulk formula for the anomalous latent heat flux can be represented as: where ρ a is the air density, L v is the latent heat of vaporization, C E is the turbulent exchange coefficient, V is the surface wind vector, u is the zonal component of the surface wind, and Δq is the air-sea humidity difference. The upper bar and prime indicate the mean and perturbation quantities, respectively. The approximation in the second row of Eq. (2) was applied because amplitude of zonal wind in an interannual timescale is much greater than that of meridional wind over the equatorial Indian Ocean 38 . For the boreal summer, we have an easterly trade wind over the equatorial IO such that u < 0 and usually u j j > u 0 j j. Therefore, the wind speed in the parentheses becomes −u′, where u′ can be linearly approximated to the east-west contrast of SST, that is, u′ = −μT. This is because the surface wind is proportional to T, i.e., DMI (east-west tropical SST gradient), and the timescale for the atmosphere to adjust to changes in SST is much faster than that for the DMI 39 . Finally, the anomalous latent heat flux becomes Q 0 LE ¼ ρ a L v C E μΔqT, indicating that the WES feedback in this case can be modified not by mean easterly trade winds over equatorial IO but by the mean differences in air-sea humidity (Δq). As surface air temperature and SST are closely linked, the change in Δq must be small.  As d½T dz is positive, WUS feedback is positive, and its intensity is proportional to the mean vertical ocean temperature gradient. The mean vertical gradient of the tropical IO between the surface and subsurface (Fig. 4d) for SA (~4.05°C/100 m) is smaller than that for LA (~4.26°C/100 m), indicating that the positive feedback by WUS for LA is greater than that for SA.
The WTS feedback (also known as thermocline feedback) is driven by anomalous vertical advection by mean upwelling such that Àw dT 0 dz % Àw ðT 0 s ÀT 0 sub Þ Δz , where T 0 s and T 0 sub indicate the ocean surface temperature and ocean subsurface temperature below the mixed layer, respectively, and Δz is the mixed layer depth. First, we neglect the mean upwelling difference between the eastern and western IO for the sake of simplicity, and thus the WTS feedback becomes % À½w ðT 0 s ÀT 0 sub Þ Δz . The DMI tendency by the first term is represented as For the second term, T 0 sub can be parameterized by a thermocline depth anomaly such that T 0 sub ¼ α h h 040,41 , where α h represents a conversion parameter from the thermocline depth anomaly to the subsurface temperature anomaly. The DMI tendency in the second term becomes αh With the aid of a dynamic balance between the east-west contrast of the ocean dynamic height and the equatorial zonal wind such that Δz T. Therefore, the WTS feedback becomes ðα h λμÀ1Þ Δz ½wT. As the thermocline feedback is positive, it must be α h λμ À 1 > 0. First, α h is inversely proportional to the mean thermocline depth 42 such that a deeper mean thermocline depth results in a smaller α h . As shown in Fig. 4d, the mean thermocline depth was shallower during SA than that during LA. Because of strong surface warming during SA, the mixed layer depth and thermocline depth were shoaling. The majority of CGCMs from the Coupled Model Intercomparison Project Phase 5 also agreed with the mean thermocline shoaling in their future projections 35 . Therefore, the  WTS associated with α h during SA was stronger than that during LA. Moreover, the stronger trade wind during SA induces a relatively stronger mean upwelling ½w, which also enhances the positive feedback by the WTS. Therefore, the WTS cannot explain the large (small) amplitude of the IOD during LA (SA).
Finally, we calculated the air-sea coupling factor, μ for the LA and SA periods using the regression coefficient between the DMI and surface zonal wind anomalies averaged over a domain of 5°S -5°N and 60-100°E. Here, μ from each ensemble for each period was computed, and 28 ensembles were then averaged. The ensemble-mean coupling factor for LA is 1.18 (ms −1°C−1 ) and that for SA is 0.83 (ms −1°C−1 ). It is unclear what drives the difference in air-sea coupling factor between LA and SA. One possible reason is related to the reduced gradient of tropospheric air temperature profile (Fig. 4c) such that ref. 35 argued that the anomalous atmospheric circulation response to the SST forcing under the global warming could be suppressed due to the relaxed stratification. Another possible reason is related to stronger near surface warming (Fig. 4c), which reduces the mean air density. When the air in planetary boundary layer is perturbed by SST change, the sensitivity of anomalous atmospheric circulation response to SSTA is reduced as the mean air density increases 43 . As seen above, all WES, WUS, and WTS feedbacks are directly related to the air-sea coupling factor; thus, the difference in the coupling factor between SA and LA must be the main cause of the difference in the DMI amplitude between the two periods as well as possibly whole period of a CO 2 removal simulation.

DISCUSSION
In this study, we analyzed the change in the IOD characteristics in a CO 2 removal simulation by using CESM1.2. Since our analysis is based on one model result, a generalization has to be done with a caution because climate models showed the amount of diversity in current and future IOD simulations 10,[35][36][37] . Nevertheless, the enhance IOD activity that was observed in the late 20th century 32 and projected in near future 34 and the suppressed moderated IOD activity projected for 21st century 44 are somewhat consistent with the IOD feature during LA and SA periods in our simulation, respectively. Therefore, the model result in this study is not so different from other models. Another caution to be considered is that the simulated IOD by CESM1.2 tends to be overestimated, which is very common in other climate models 36 . Nevertheless, since our main purpose is to compare IOD characteristics under two different climate states forced by ramp-up and ramp-down CO 2 forcing, a bias effect may be small via focusing on their relative changes.
During the ramp-up period, the intensity of IOD tended to increase in early period, like the observed in recent decades [31][32][33] and a near future projection of IOD 34 , but as the CO 2 concentration increasement passed a certain point, the IOD intensity decreased. The minimum IOD amplitude is recorded during the early ramp-down period and the amplitude recovered afterward. A simple model for IOD was developed to investigate the process leading to changes in the intensity of IOD during changes in CO 2 concentrations. This simple model includes comprehensive local air-sea coupled feedback, delayed negative feedback, ENSO forcing, and stochastic forcing. Among them, we found that the change in the local air-sea coupled feedback due to CO 2 forcing mainly changes the IOD amplitude. For the large amplitude period, the growing season of IOD controlled by the local air-sea coupled feedback was longer, while it was shorter for the small amplitude period, and the peak appeared during earlier calendar months. ENSO characteristics also showed some changes during changes in CO 2 concentrations, but the impact of ENSO on IOD during both small and large amplitude periods was similar.
Finally, we showed the mean state changes in the tropical Indian Ocean owing to the changes in local feedback. The air-sea coupling strength and mean vertical ocean temperature gradient during LA were larger than those during SA. The former enhances WES, WUS, and WTS feedback, and the latter enhances WUS during LA. The deeper mean thermocline depth during LA compared to that during SA reduces the WTS feedback. Therefore, the air-sea coupling strength and mean vertical ocean temperature gradient are responsible for the amplitude change of the IOD during changes in CO 2 concentrations. Although this approach to seeking specific feedbacks is qualitative, the simple IOD model quantitatively revealed the IOD changing mechanism. Experimental design of CESM1.2

METHODS Earth System Model configurations
We used two ideal CO 2 scenarios. One was a constant CO 2 concentration level in the atmosphere, and the other had a varying CO 2 concentration 49 . These scenarios are almost the same as the protocol of the carbon dioxide removal model intercomparison project (CDR-reversibility) 50 except with a different initial CO 2 level. In this study, the present-day CO 2 level was used instead of the pre-industrial level. First, the constant CO 2 scenario was conducted for 900 years, with 367 ppm as the present-day climate simulation (PD period). Then, under the varying CO 2 scenario with 28 ensemble members, the CO 2 concentration increased by 1% per year for 140 years until it quadrupled (i.e., 1468 ppm) (ramp-up period), followed by a decrease in a symmetric pathway (about 1% per year) for 140 years until the CO 2 concentration reached the initial state (367 ppm) (ramp-down period). Subsequently, a constant CO 2 (367 ppm) simulation was conducted for 220 years (restoration period). The 28 ensemble members were identical except for the initial oceanic conditions.

A simple model for IOD and dynamics
A simple model for IOD was adopted from the previous study 30 and modified by including the delayed feedback process. In addition to a delayed feedback term, λ n and β n are not fitted to sinusoid function. The ref. 30 demonstrated that the original IOD model well reproduces the observed IOD. In Eq. (1), parameters λ n , α n , and β n are computed using the CESM1.2 produced DMI and NINO3.4 indices by applying a least-square method at each calendar month, and σ n is computed as the residual using the computed λ n , α n , and β n . The subscript 'n' indicates the calendar month (n = 1, …, 12), and the parameters for each calendar month were separately calculated. Using the computed parameters, the DMI was reproduced by integrating the simple IOD model.
Comprehensive local feedback includes the wind-thermocline-SST (WTS) feedback, wind-evaporation-SST (WES) feedback, cloud-radiation-SST (CRS) feedback 21,23 , and wind-upwelling-SST (WUS) feedback. The WES, CRS, and WUS feedbacks operate simultaneously during an IOD event 21,22 ; thus, these can be parameterized as simply λ n T. However, the WTS feedback involves dynamic ocean processes. For the ENSO case, the thermocline response to the ENSO-induced surface wind stress anomalies is not simultaneous but lags because of oceanic wave propagation and subsurface adjustment. On the one hand, a much smaller basin scale of IO can possibly make oceanic wave response, especially that of equatorial Kelvin waves, to be almost simultaneous (time for half-basin crossing by an equatorial Kelvin wave is~0.5 months). Thus, the WTS feedback associated with a forced equatorial Kelvin wave can be approximated as λ n T. On the other hand, the Rossby wave response cannot be simultaneous because its speed is three times slower than that of a Kelvin wave. Thus, WTS feedback in the western IO associated with Rossby waves and WTS feedback in the eastern IO associated with reflected Kelvin waves originating from Rossby waves involve a delay. The delayed time estimated by wave speeds and ocean basin size is about 2 months 42 . Therefore, δ in Eq. (1) is two months. These wave dynamics, which can be analogous to the delayed oscillator theory of ENSO in the Pacific Ocean 42 , also work in the IO. It has been also argued that changes in heat content along the equator occurred before the onset of IOD events, but its role in generating IOD variability is relatively weaker than that in the Pacific Ocean 51 . Interestingly, ref. 52 proposed that the equatorial Indian Ocean subsurface can be pre-conditioned such that the water from the sub-surface of Southern Ocean, which starts about 8-10 years prior to an IOD, reaches the equatorial Indian ocean 2-3 years prior to the event, and awaits a trigger. Note that one important missing process must be nonlinear process including nonlinear advection, which influences extreme IOD event 52 .

CODE AVAILABILITY
Any codes used in the manuscript are available upon request from H.-J. Park, hjpark1021@yonsei.ac.kr.