The Subsurface and Surface Indian Ocean Dipoles and Their Association with ENSO in CMIP6 models

The leading Empirical Orthogonal Function (EOF) pattern of the subsurface temperature variability in the tropical Indian Ocean (TIO) closely resembles the Indian Ocean Dipole (IOD) pattern. This EOF pattern is represented by a quasi-biennial east-west dipole oscillation of the subsurface temperature anomalies, known as the subsurface Indian Ocean Dipole (Sub-IOD) (Rao et al., 2002; Rao and Behera, 2005; Sayantani and Gnanaseelan, 2015; Song et al., 2021). A positive Sub-IOD event is characterized by cool subsurface temperature anomalies off the coasts of Java and Sumatra and anomalous warming in the western TIO that peaks around November–January. The converse is true for the negative phase. The development of Sub-IOD events is related to significant zonal wind anomalies over the TIO, which act to gradually deepen or reduce the thermocline depth in the TIO by inducing zonally propagating equatorial and off-equatorial oceanic waves (Huang and Kinter, 2002; Rao and Behera, 2005; Yuan and Han, 2006; Yuan and Liu, 2009; Sun et al., 2010).

Most studies have shown that the Sub-IOD is closely associated with the surface IOD that usually occurs in boreal summer-to-autumn (Sayantani and Gnanaseelan, 2015; Kakatkar et al., 2020), and the Sub-IOD oscillation may have provided a baseline for the quasi-biennial variation of the surface IOD (Song et al., 2021). Although the Sub-IOD is well correlated with the surface IOD, the correlation between the Sub-IOD index and Niño-3.4 is quite weak (Shinoda et al., 2004; Song et al., 2021). Rao and Behera (2005) and Yu et al. (2005) indicated that the region north of 10°S in the TIO during a Sub-IOD event is mainly affected by ENSO. Sayantani and Gnanaseelan (2015) believe that the Sub-IOD mode is first triggered by the IOD and then governed by ENSO. Recent evidence has shown that during a Central Pacific (CP) ENSO event, characterized by maximum sea surface temperature anomalies (SSTA) over the central tropical Pacific (Ashok et al., 2009; Kao and Yu, 2009; Yeh et al., 2009; Xiang et al., 2013), a significant westwardly displaced SSTA center is present and peaks from winter (year +0) to spring (year +1). Such a event can promote the development of Sub-IOD events in (year +0) to spring (year +1), while the surface IOD (year +0) is not active (Song and Ren, 2022). Therefore, with the development of model simulations and the extension of reanalysis datasets, the mutual linkages among the Sub-IOD, IOD, and ENSO still require further investigation.

Numerous studies have evaluated the performance of coupled climate models from the third phase of the Coupled Model Inter-comparison Project (CMIP3) to the latest sixth phase (CMIP6) (Meehl et al., 2007; Eyring et al., 2016) in the simulations of ENSO and the surface IOD (Saji et al., 2006; Guilyardi et al., 2009; Yu and Kim, 2010; Kim and Yu, 2012; Weller and Cai, 2013; Bellenger et al., 2014; Planton et al., 2021). It has been revealed that most models suffer from common biases in their simulations of the intensity and diversity of the IOD and ENSO from CMIP3 to CMIP6. For example, the unrealistic thermocline tilt in models over a TIO accompanied by too strong a mean easterly wind and an overly strong east-to-west sea surface temperature (SST) gradient exerts a strong influence on the thermocline-SST feedback (Cai and Cowan, 2013; Wang et al., 2021). Therefore, the simulated positive IOD-like pattern is exaggerated in the TIO due to the overly strong Bjerknes feedback (Vecchi and Soden, 2007; Cai and Cowan, 2013; Li et al., 2015; Wang et al., 2021). The farther westward extent of the Pacific Ocean cold tongue in many CMIP3 and CMIP5 models can also affect the model's ability to reproduce ENSO (Li and Xie, 2014; Tian and Dong, 2020). However, few studies have evaluated the ability of the models to simulate a Sub-IOD closely linked to the IOD and ENSO. This study focuses primarily on the Sub-IOD and examines the reproducibility of CMIP6 models for the Sub-IOD, surface IOD, and their association with ENSO. The results can potentially lead to more reliable predictions and future projections of the regional and global climate.

The remainder of this paper is organized as follows. Detailed descriptions of the models, reanalysis data, and analysis methods are given in section 2. Section 3 evaluates 31 coupled CMIP6 models regarding their ability to reproduce the basic properties of the Sub-IOD. In section 4, we focus on the association of the Sub-IOD and IOD with the canonical ENSO and the possible linkage of the Sub-IOD with CP-ENSO. A summary and discussion are given in section 5.

We analyze the historical simulations from 31 CGCMs participating in the CMIP6 (Table 1). Only one integration ensemble member is used for all models. Output of the monthly SST, wind, and ocean temperature from these models are remapped onto a 1° × 1° regular grid. The observational reanalysis of the monthly ocean temperature data is obtained from the Ocean Reanalysis System 4 (ORAS4) for 1958–2015 from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Balmaseda et al., 2013), which is also remapped onto a 1° × 1° regular grid. ORAS4 has been widely used to study variability in the equatorial Indian Ocean (Deepa et al., 2019; Kakatkar et al., 2020; Mohapatra et al., 2020) and was found to be one of the best reanalysis products over the Indian Ocean (Karmakar et al., 2018). The monthly wind from 1958 to 1978 and from 1979 to 2015 are obtained from the ECMWF 40 reanalysis (ERA-40) and ECMWF interim reanalysis (ERA-Interim) (Uppala et al., 2005; Dee et al., 2011), respectively. Oceanic heat content anomalies (HCA) are defined as the mean ocean temperature anomalies multiplied by the density (1025 kg m^–3) and the specific heat (3850 J kg^–1 °C^–1) of seawater for the upper 315 m of depth. The equatorial thermocline depth is defined by the depth of 20°C isotherm and the thermocline tilt is calculated as the difference in thermocline depth between the western (50°−70°E, 10°S−10°N) and eastern TIO (80°−95°E, 10°S−10°N).

An index representing the Sub-IOD (SDI) is defined as the normalized difference in the area-mean HCA between the western (50°–70°E, 5°S–5°N) and eastern (80°–100°E, 5°S–5°N) TIO (e.g., Rao et al. 2002 and their Fig. 3; Kakatkar et al. 2020 and their Fig.1; Song et al. 2021 and their Fig. 10). Following Saji et al. (1999), the IOD index is defined as the normalized SSTA difference between the western (50°–70°E, 10°S–10°N) and southeastern (90°–110°E, 10°S–0°) TIO. Sub-IOD (IOD) events are defined when the three-month running mean normalized SDI (DMI) exceeds its 0.9 (1.0) standard deviation. It should be pointed out that the IOD index defined from ORAS4 is highly consistent with that from the gridded SST products HadISST (correlation coefficient up to 0.81 and significant at the 99% confidence level). Therefore, we use the ORAS4 as the observational reanalysis dataset. As for ENSO, the Niño-3, Niño-3.4, and Niño-4 indices are defined as the area-mean SSTA respectively over the (5°N–5°S, 150°–90°W), the (5°N–5°S, 170°–120°W), and the (5°N–5°S, 160°E–150°W). In general, we use the Niño-3.4 index to describe the evolution and strength of ENSO and further define an index for CP ENSO (CPI) as:

where $ \alpha=0.4 $ when $\mathrm{N}\mathrm{i}\tilde{\mathrm{n}}\mathrm{o}4$×$\mathrm{N}\mathrm{i}\tilde{\mathrm{n}}\mathrm{o}3\gg 0$ and $ \alpha =0 $ when $\mathrm{N}\mathrm{i}\tilde{\mathrm{n}}\mathrm{o}4$×$\text{Ni}\tilde{\mathrm{n}}\mathrm{o}3\ll 0$ (Ren and Jin, 2011). Similarly, a positive (negative) three-month running mean, CPI normalized to be greater than 1.0 standard deviation indicates a CP El Niño (La Niña) event.

We first performed an EOF analysis on the HCA over the TIO for all 31 CMIP6 models to examine the ability of models to reproduce the leading east-west dipole mode of subsurface sea temperature variability in the observations. Figure 1 shows the leading EOF patterns of the TIO HCA in ORAS4 and the 31 models. All the CMIP6 models can simulate the Sub-IOD pattern over the TIO reasonably well except MCM-UA-1-0, which shows a basin-wide cooling of the HCA. For most models (18/31), the percentage of the variance explained by the leading EOF is much higher than that in ORAS4 (26.0%), implying an overestimated dominance of Sub-IOD variability in most models. For MCM-UA-1-0, the highest percentage of variability (53.7%) is attributed to its leading basin-wide mode. Among those models exhibiting lower variance loading than ORAS4, the HadGEM3-GC31-MM shows the lowest percentage of variability (14.7%) for its first Sub-IOD-like mode among the 31 models, suggesting the underestimated dominance of Sub-IOD variations.

By examining the intensity and seasonality of the Sub-IOD in CMIP6 models, it is found that the simulated Sub-IOD also tends to be phase-locked to certain seasons as in ORAS4, beginning during its development around August–October, peaking in NDJ (November–December), and decaying quickly thereafter, as manifested by the monthly standard deviation of SDI shown in Fig. 2a. Variability of the Sub-IOD in its peak season (NDJ) is overly strong in most (24/31) CMIP6 models, only 7/31 models simulate a weaker variability than that in ORAS4 (Fig. 2a).

Figure 1.  The first EOF modes of HCA over the TIO for the years 1958–2015 from ORAS4 and for 1850–2014 from the historical simulations of 31 CMIP6 models. The top right corner of each graph shows the percentage of variability explained by EOF1.

Figure 2.  (a) The standard deviation of SDI by month in 31 CMIP6 models sorted in ascending order from top to bottom; (b) Power spectrum of the SDI in 31 CMIP6 models. The multi-model means and that from ORAS4 are shown at the top for easy reference.

It is known that the Sub-IOD is an inherent east−west oscillation mode in the TIO, owning a quasi-periodicity of around 15−20 months (Song et al., 2021). We also perform a power spectrum analysis on the simulated SDI from CMIP6 models (Fig. 2b). It is seen that a larger part (21/31) of the CMIP6 models have a much longer dominant period of over three years, while a smaller part (10/31) including EC-Earth3-AerChem, EC-Earth3, E3SM-1-0, E3SM-1-1-ECA, ACCESS-CMS2, and BCC-ESM1 exhibit a period peak of 1.5–2-years (Fig. 2b). Aside from the overestimated dominant period of SDI, some models including FGOALS-f3-L, MRI-ESM2-0, and NESM3 also exhibit a much sharper spectral peak suggesting a highly oscillatory Sub-IOD in these models. The overly narrow spectral peak of the IOD in CMIP6 models was also found by McKenna et al. (2020).

As noted in the introduction, Sub-IOD events are usually triggered by internal dynamics in the TIO and are closely coupled with IOD events. By checking the lead-lag correlation between the Sub-IOD and IOD in individual models, it is found that most CMIP6 models can simulate the relationship between the Sub-IOD and IOD, showing a maximum correlation when the SDI lags DMI by ~1 month, except in the GISS-E2-1-H and GISS-E2-1-G, both of which show a much weaker correlation (Fig. 3a). However, the correlation in most (25/31) CMIP6 models is stronger than that in ORAS4 where the correlation coefficient is 0.48 (Fig. 3a). This is generally suggestive of a much stronger coupling between the Sub-IOD and IOD in CMIP6 models, or the processes that are responsible for their mutual interactions are simulated to a much stronger extent.

Figure 3.  (a) Lead-lag correlations between the SDI and DMI in 31 CMIP6 models (sorted in ascending order from top to bottom) with the multi-model mean and that from ORAS4 shown at the top for easy reference. The abscissa is the lag month of SDI against DMI. (b) Model bias in the mean state of thermocline depth in the western (50°–70°E, 10°S–10°N; dashed black bar) and eastern TIO (80°–95°E, 10°S–10°N; solid black bar) and in thermocline tilt (red bar) as the difference in thermocline depth between the western and eastern TIO.

Since subsurface sea temperature variability can affect the SSTA over the TIO through Bjerknes feedback, which involves surface zonal winds, SST, and the state of the thermocline in the TIO (Bjerknes, 1969), the intensity of the thermocline-SST feedback should be closely related to the climatological mean ocean–atmosphere state, particularly the east-west tilt of the thermocline (Cai and Cowan, 2013; Wang et al., 2021). Figure 3b illustrates the model bias in thermocline depth in the western and eastern TIO, as well as the thermocline tilt represented by the depth difference between the western and eastern TIO. Though the thermocline depth in the eastern TIO in some CMIP6 models is deeper than that in ORAS4, most models generate an even deeper thermocline in the western TIO, which results in an even steeper thermocline tilt from the western to eastern TIO (Fig. 3b). The much steeper thermocline tilt favors an overly strong Bjerknes feedback over the TIO when anomalous easterlies appear over the TIO, which gives rise to an even stronger coupling between the surface and subsurface (Cai and Cowan, 2013; Prajeesh et al., 2022). The gradually reduced model bias in the thermocline tilt with the gradually reduced correlation between SDI and DMI can be clearly seen from the bottom to the top of Fig. 3 (Fig. 3a vs. 3b). The inter-model correlation between them is 0.35 (significant at 95% confidence level). Nevertheless, a few models (e.g., MCM-UA-1-0, GISS-E2-1-G) with a steeper thermocline tilt still show a relatively weaker correlation between DMI and SDI, which may be due to the overly deep thermocline in the eastern TIO. Aside from the tilt of the thermocline, the surface–subsurface coupling may be influenced by other factors, including the zonal wind, oceanic Rossby waves, and other oceanic processes (Chakravorty et al., 2014a, b).

Since IOD events that occur in autumn, as well as the accompanied Sub-IOD events, are usually influenced by the development of the canonical ENSO events that mature in winter, we first examine the models’ reproducibility of the seasonality of the canonical ENSO following McKenna et al. (2020). It is seen from Fig. 4 that most models can reproduce the winter-maturing ENSO. Among the 31 models, about 13 models are relatively poor in reproducing a proper ENSO seasonality. For example, a few models exhibit a dominant ENSO peak in other seasons (e.g., E3SM-1-1, BCC-ESM1, BCC-CSM2-MR, GISS-E2-1-G, MCM-UA-1-0, MPI-ESM1-2-HR, and ACCESS-ESM1-5), while a few other models exhibit double peaks of ENSO (GISS-E2-1-H, GISS-E2-2-H, CAMs-CSM1-0, CANESM5-CanOE, MPI-ESM1-2-LR, and CANESM5). In this regard, we exclude these 13 poor models when examining the relationship of the IOD and Sub-IOD with the canonical ENSO.

Figure 4.  Similar to Fig. 2a but for the standard deviation of Niño-3.4 by month.

By checking the lead-lag correlation between the Sub-IOD and canonical ENSO in the remaining 18 models, it is seen that almost all of them can simulate the lead-lag relationship between the Sub-IOD and ENSO, generally showing a maximum positive correlation when the Niño-3.4 lags the SDI by 0~1 month (Fig. 5a). However, the correlations in most of the models are significantly stronger than that in the ORAS4 (Fig. 5a), which is consistent with the overly regular Sub-IOD events in models [our Fig. 2 and also Fig. 3b in McKenna et al. (2020)]. Comparing Fig. 3a and Fig. 5a, it is evident that the stronger (weaker) coupling between the IOD and Sub-IOD largely corresponds to a stronger (weaker) correlation between ENSO and the Sub-IOD across models. This may be because the stronger (weaker) surface-subsurface coupling yields more (less) Sub-IOD events in autumn-to-winter, which thus favors a stronger (weaker) linkage with the canonical ENSO. Of course, a stronger (weaker) linkage between the Sub-IOD and canonical ENSO can, in turn, produce feedback to support a stronger (weaker) surface-subsurface coupling. In addition, the overestimated correlation between the Sub-IOD and canonical ENSO may also explain the simulated overall longer timescale of the Sub-IOD indicated in Fig. 2b.

Figure 5.  Similar to Fig. 3a, but for the lead-lag correlations (a) between the SDI and Niño-3.4, and (b) between the DMI and Niño-3.4 in ORAS4 and the 18 CMIP6 models. The abscissa is the lag month of Niño-3.4 against the SDI (DMI). The sequence of the ordinates is consistent with that in Fig. 3, except that the 13 models which poorly reproduced the canonical ENSO are excluded.

Accordingly, the lead-lag correlation between the IOD and canonical ENSO is displayed in Fig. 5b. Though the lag time varies across models, most of them generally show the correlation peak when the canonical ENSO lags the DMI by about 1−3 months, which is similar to that in ORAS4. Particularly, the lagged correlation also gets stronger (weaker) when the coupling between the IOD and Sub-IOD becomes stronger (weaker) across models (Fig. 5b), again manifesting as the linkage between the strength of surface-subsurface coupling and the associations of the Sub-IOD and IOD with the canonical ENSO.

To further demonstrate the linkage between the strength of surface-subsurface coupling and the associations of the Sub-IOD and IOD with the canonical ENSO, Fig. 6 displays the scatter diagrams of the Sub-IOD vs. IOD, Sub-IOD vs. ENSO, and IOD vs. ENSO events in ORAS4 and the three selected models that exhibit a strong (FGOALS-f3-L), a medium (EC-Earth3) and a weak (NESM3) surface-subsurface coupling, respectively. For the ORAS4 from 1958–2015, the IOD and the accompanied Sub-IOD were significantly correlated at 0.52, and correlations between the Sub-IOD and ENSO and between the IOD and ENSO are also significant at 0.39 and 0.32, respectively (Figs. 6a–c). This corresponds to a situation where some Sub-IOD and IOD events have occurred under opposing ENSO conditions. For FGOALS-f3-L, with a much stronger surface-subsurface coupling, the correlation coefficients are 0.78, 0.62, and 0.51, respectively, which are much larger than that in ORAS4 (Figs. 6d–f). Consequently, there are fewer numbers of Sub-IOD and IOD events (1 or 2) that occurred under opposing ENSO conditions, indicating a highly consistent relationship with ENSO. EC-Earth3 shows slightly lower correlation coefficients (0.68, 0.51, 0.31) and an increased number of Sub-IOD and IOD events (more than 5) occurred under opposing ENSO conditions (Figs. 6g–i). As expected, associated with the weaker surface-subsurface coupling (0.46) than ORAS4, NESM3 also shows a smaller correlation (0.39, 0.27) with ENSO and a much more scattered distribution of the events in the bottom panel (Figs. 6j–l). In general, the surface-subsurface coupling is stronger in most models, and the IOD and Sub-IOD are more closely associated with the canonical ENSO.

Figure 6.  Scatter diagrams of the Sub-IOD vs. IOD (left, with Sub-IOD lagging IOD by one month), Sub-IOD vs. ENSO (middle, with ENSO lagging Sub-IOD by one month), and IOD vs. ENSO (right, with ENSO lagging IOD by two months) and their linear fitting line for (a–c) ORAS4, for (d–f) FGOALS-f3-L, for (g–i) EC-Earth3 and for (j–l) NESM3. Correlation coefficients are shown in each panel, and the superscript “*” indicates that the correlation coefficient is statistically significant at or above the 95% confidence level. All indices are normalized.

To examine whether the linkage between the Sub-IOD and westward-centered CP-ENSO events occurring in winter–spring (Song and Ren, 2022) can be reproduced by CMIP6 models, we first display the seasonality of the Sub-IOD (Fig. 7a) and IOD (Fig. 7b) events in ORAS4 and 31 models. In ORAS4, most Sub-IOD events occur from autumn into winter (November–February) (Fig. 7a), while most IOD events are observed strictly during the autumn months (September–November) (Fig. 7b). This is also largely true in most models; however, this phase-locking feature is much more prominent in models than that in ORAS4. This is consistent with the more regular IOD vs. Sub-IOD events due to stronger surface-subsurface coupling and a tighter association of IOD vs. Sub-IOD events with ENSO, as indicated in the above section.

Figure 7.  Similar to Fig. 3, but for the number per century of (a) Sub-IOD and (b) IOD events peaking in each calendar month for ORAS4 and the CMIP6 models.

Figure 8 further shows the percentages of Sub-IOD events which co-occur with the IOD in autumn–winter as well as those which co-occur with CP-ENSO in winter–spring. It can be seen that the percentage of Sub-IOD events co-occurring with IOD events (or the regular events) generally decreases along with the weakening of the surface-subsurface coupling moving from left to right in Fig. 8. Meanwhile, the percentage of independent Sub-IOD events that co-occur with neither the IOD nor CP-ENSO significantly increases along with the weakening of surface-subsurface coupling, which is understandable because more and more Sub-IOD events are losing their connection with the surface due to the weakened surface-subsurface coupling. The percentage of Sub-IOD events co-occurring with CP-ENSO is steady around 10% in the multi-model mean. Next, we choose three models that exhibit the highest (18.2%, FGOALS-f3-L), intermediate (10.4%, GFDL-ESM4), and lowest percentages (2.3%, GISS-E2-1-H) of these winter–spring Sub-IOD events to examine the existence of their linkage with CP-ENSO.

Figure 8.  Percentage of the Sub-IOD events (orange bars) co-occurring with the IOD in autumn (year +0) to the total number of Sub-IOD events, and that of the Sub-IOD events (blue bars) co-occurring with CP-ENSO during winter (year +0) into spring (year +1) and also independent of the IOD in autumn (year +0) in ORAS4 and 31 CMIP6 models. Ratios of the second (blue bars) to the first type (orange bars) of Sub-IOD events are marked above the corresponding bars. The sequence of the 31 models along the abscissa is consistent with that in the ordinate of Fig. 3.

Figure 9 displays the composite subsurface temperature patterns for these winter–spring Sub-IOD events in ORAS4 and the three selected models. An east-west, seesaw-like SST dipole in ORAS4 is well reproduced by all three models. Associated with the stronger, medium, and weaker surface-subsurface coupling, the thermocline tilt in the TIO is stronger in FGOALS-f3-L (Fig. 9b) and GFDL-ESM4 (Fig. 9c) and weaker in GISS-E2-1-H (Fig. 9d) than that in ORAS4 (Fig. 9a). Associated with this, a CP-ENSO SSTA pattern exists over the central tropical Pacific in the three models, similar to ORAS4 (Fig. 10a vs. Figs. 10b-d). However, the intensity of the main SSTA center is slightly stronger, especially in FGOALS-f3-L, and weaker in GISS-E2-1-H. The SSTA center is also obviously westward centered, as in ORAS4, except for that in FGOALS-f3-L, where the SSTAs are much stronger over the central and eastern Pacific (Fig. 10b). In general, whether there are more or less winter–spring Sub-IOD events in models depends upon the westward-centered CP-ENSO over the equatorial Pacific, or the linkage between the Sub-IOD and CP-ENSO events identified in Song and Ren (2022). This teleconnection can be reproduced by most CMIP6 models.

Figure 9.  Longitude–depth cross-sections of the positive-minus-negative composite of upper-ocean temperature anomalies (averaged over 5°S–5°N; units: °C) in the TIO during JFM for the Sub-IOD events co-occurring with CP ENSO in (a) ORAS4, (b) FGOALS-f3-L, (c) GFDL-ESM4, and (d) GISS-E2-1-H. Dotted areas indicate composites that are significant at or above the 80% confidence level based on a Student’s two-tailed t-test.

Figure 10.  As in Fig. 9 but for the composite SSTAs (shading; units: °C) and 10-m wind anomalies (vectors; units: m s⁻¹) over the tropical Indo-Pacific region. Dotted areas and vectors indicate that the composite anomalies are statistically significant at or above the 80% confidence level based on a two-tailed Student’s t-test.

Sub-IOD, or the east–west dipole oscillation of the subsurface temperature in the TIO, is the first leading mode of the interannual variability of ocean temperature. Many previous studies have confirmed the presence of this mode in different datasets (e.g., Kakatkar et al., 2020). Here we use 31 CMIP6 historical simulations to examine the fidelity of the models in representing the Sub-IOD characteristics and their association with IOD and ENSO. It is found that most CMIP6 models can reproduce the Sub-IOD mode and the main features of its spatial pattern, occurrence seasonality, and its association with the IOD and canonical ENSO. However, the intensity and dominant period of the Sub-IOD are overestimated in most CMIP6 models. The overly exaggerated lead-lag correlation between the Sub-IOD and IOD events indicates a bias towards a stronger coupling between the surface and subsurface, which is associated with the bias of the much steeper west-to-east thermocline tilt in the TIO in most CMIP6 models.

Although the simulated lead-lag relationship between the Sub-IOD (IOD) and canonical ENSO is poorly represented in some models, the general diversity of the relationship between the Sub-IOD (IOD) and canonical ENSO can also be reproduced by a larger part of the 31 CMIP6 models, as long as the model can properly simulate the seasonality of ENSO. Most models generally simulate a tighter association with the canonical ENSO; thus, a more regular pattern of Sub-IOD and IOD events occurs in autumn compared to ORAS4.

A stronger (weaker) relationship among these indices is associated with a stronger (weaker) surface-subsurface coupling. The stronger (weaker) surface-subsurface coupling also yields a decrease (increase) of the independent Sub-IOD events that co-occur with neither the IOD nor CP-ENSO. Nevertheless, the dynamical linkage between Sub-IOD and CP-ENSO is captured by CMIP6 models, but the number of such events and the intensity of the SSTA anomalies varies among models.

In general, most CMIP6 models simulate a much stronger surface-subsurface coupling in the TIO, thought to be related to a much steeper east-west thermocline tilt in the TIO. This yields a much more regular overall pattern of Sub-IOD and IOD events in autumn associated with ENSO. Therefore, a stronger (weaker) surface-subsurface coupling also corresponds to a tighter (looser) association of Sub-IOD and IOD events with the canonical ENSO. Though the Sub-IOD events which co-occur with CP-ENSO are relatively few among the 31 models, the dynamical linkage between them can be well reproduced. The concrete evidence provided by this study can support future model improvements related to ENSO and IOD reproducibility and predictability.

Acknowledgements. This work was supported by the National Key R&D Program of China (Grant No. 2019YFA0606701) and the Guangdong Major Project of Basic and Applied Basic Research (Grant No. 2020B0301030004). We acknowledge the World Climate Research Programme’s Working Group on Coupled Modeling, responsible for CMIP, and thank the climate modeling groups for producing and making the model output available. We also acknowledge the European Centre for Medium-Range Weather Forecasts (ECMWF) for providing the Ocean Reanalysis System 4 (ORAS4), ERA-40, and ERA-Interim Reanalysis datasets.