Origin of Indian Ocean multidecadal climate variability: role of the North Atlantic Oscillation

The multidecadal variability of Indian Ocean sea surface temperature (IOSST) has an important impact on both the regional Indian Ocean climate and the global climate. Here, we explore multidecadal variability in the annual IOSST. Observational analysis shows that the annual IOSST multidecadal variability is not only related to the Pacific Decadal Oscillation (PDO), but also to the North Atlantic Oscillation (NAO). The NAO leads by 15–20 years the detrended annual IOSST in which the PDO signal of the same period has been removed. Further analysis reveals that the NAO leads the annual IOSST multidecadal variability through its leading effect on the Atlantic Multidecadal Oscillation (AMO). The AMO affects the vertical wind anomaly in the Indian Ocean region through the Atlantic–Indian Ocean multidecadal teleconnection (AIMT), which in turn affects the net longwave radiation in the Indian Ocean region, thus driving the annual IOSST multidecadal variability. A Hasselmann model based on NAO and PDO further verify the joint influence of the NAO and PDO on the multidecadal variability of the IOSST. A PDO-based linear model and a climate model that incorporates the NAO signal are also constructed for the annual IOSST. Results show that the climate model with the NAO signal can better simulate the annual IOSST. This again verifies that the NAO is part of the annual IOSST multidecadal variability source, indicating that the annual IOSST variability may be due to the combined influences of the NAO and PDO.


Introduction
The climate system in the Indian Ocean is complex and can play an important role in both regional and global climate variability (Ashok et al. 2001;Annamalai and Murtugudde 2004;Du and Xie 2008;Schott et al. 2009;Reverdin et al. 2010;Dong and McPhaden 2018;Zhang et al. 2019). Meteorological and climatic disasters often occur in the Indian Ocean region and surrounding areas, e.g., the record-breaking precipitation occurred in East Africa during October and November 1997, killing thousands and leaving hundreds of thousands of people displaced, and at the same time, Indonesia suffered from severe drought with uncontrolled fires on several of its islands (Saji et al. 1999;Webster et al. 1999;Annamalai and Murtugudde 2004;Schott et al. 2009). In March 2019, the Indian Ocean tropical cyclone Idai swept through East Africa, causing one of the worst weather related disasters in the Southern Hemisphere (Wikipedia 2019). The Indian Ocean sea surface temperature (IOSST) also has an impact on the El Niño-Southern Oscillation (ENSO; Yu et al. 2002;Wu and Kirtman 2007;Izumo et al. 2010;Cai et al. 2019), the East Asian climate (Guan and Yamagata 2003;Annamalai and Sperber 2005;Yang et al. 2007), the North Atlantic climate (Bader and Latif 2003), and etc.
The IOSST has multiscale characteristics ranging from seasonal through interannual and interdecadal to multidecadal or even longer timescale (Huang and Shukla 2007;Lee and McPhaden 2008;Schott et al. 2009;Trenary and Han 2013;Han et al. 2014;Dong et al. 2016;Zhang et al. 2018;Cai et al. 2019). Due to its strong seasonal variability, previous studies have focused on the intra-seasonal to interannual variability of the IOSST, but less on its multidecadal variability (Allan et al. 1995;Dong et al. 2016). Significant multidecadal variabilities in annual IOSST have been recorded since the twentieth century ( Fig. 2; Allan et al. 1995;Dong et al. 2016), and have been determined from the seventeenth century onward based on coral isotope records (Damassa et al. 2006). Multidecadal temperature variabilities are closely related to human life (Yoshino 1978;Li et al. 2013;Robeson et al. 2014;Keenlyside et al. 2015;Stolpe et al. 2017) and are also a focus of current international climate research (Taylor et al. 2012;Stocker et al. 2014;Stolpe et al. 2017). Thus, recognizing the source of Indian Ocean multidecadal variability is very important.
Some studies have found that ocean-atmosphere coupling in the Indian Ocean (Lee 2004;Yamagata et al. 2004;Schott et al. 2009;Rintoul 2013;Nagura and McPhaden 2018), ENSO (Baquero-Bernal et al. 2002;Gualdi et al. 2003;Lau and Nath 2004;Yu and Lau 2005;Luo et al. 2010), and Indonesian throughflow from the Pacific (Godfrey and Weaver 1991;Hughes et al. 1992;Qu et al. 1994;Godfrey 1996;Potemra and Schneider 2005;Halkides and Lee 2009;Feng et al. 2018) can affect the IOSST. However, few studies have examined their relationship on multidecadal time scales. Krishnamurthy and Krishnamurthy (2016) found that the Pacific Decadal Oscillation (PDO) has an effect on the interdecadal characteristics of Indian Ocean Dipole. Dong et al. (2016) found that the Pacific SST multidecadal variability has a great driving effect on the IOSST based on observation and simulation experiments. However, the relationship between the annual IOSST and PDO may be not stable. As shown in Fig. 1, sliding correlations on a 25-year moving window between them was significantly positive before 1990, but sharply weakened and turned to negative after 1990, implying their relationship has broken down in recent decades, and this is also consistent with previous studies (Dong and McPhaden 2018;Zhang et al. 2018). Therefore, this result suggests that there may be factors other than PDO that can affect the IOSST. Some studies pointed out that the change in the relationship between the IOSST and PDO is affected by external forcing effect as anthropogenic greenhouse gas forcing and volcanic eruptions (Dong and McPhaden 2018;Zhang et al. 2018). The influence of internal variability on the change in the relationship between the IOSST and PDO is still unclear. Sun et al. (2019) found that the SST multidecadal variability in the Arabian Sea is affected by the Atlantic Multidecadal Oscillation (AMO) and the Western Pacific SST multidecadal variability. This study finds that the AMO can also affect the IOSST multidecadal variability, not just the Arabian Sea   . The solid green line is the 11-year running mean of annual IOSST anomalies. b Power spectrum of annual IOSST anomalies for the period 1900-2017. The blue and red dashed lines show the 95% confidence level and the reference red noise spectrum, respectively region. Moreover, the North Atlantic Oscillation (NAO), as an early signal of the AMO (Li et al. 2013(Li et al. , 2018Sun et al. 2015;Wang et al. 2017;Xie et al. 2019), may have an important role in the source of the IOSST multidecadal variability, besides the PDO.
In this paper, we offer preliminary evidence of a lagged correlation between the annual IOSST anomalies and the NAO 16 years earlier. The mechanism by which the NAO influences the annual IOSST multidecadal variability is also investigated. The arrangement of the text is as follows. Section 2 describes the data and methodology usesd in this study. In Sect. 3, we investigate multidecadal variability in the annual IOSST and their relationship to the NAO. In Sect. 4, we investigate the mechanism by which the NAO influences multidecadal variability in the annual IOSST. In Sect. 5, we build a multidecadal climate model of the annual IOSST through combining the signals of both the NAO and PDO. A summary and discussion are provided in Sect. 6.

Data
For sea level pressure (SLP) data we used the National Center for Atmospheric Research (NCAR) SLP observational dataset (Trenberth and Paolino 1980) (Kaplan et al. 1998) data from the National Oceanic and Atmospheric Administration (NOAA), at a horizontal resolution of 1° × 1°, 2° × 2° and 5° × 5°, respectively. As the relationship between the NAO and annual IOSST is not sensitive to the selection of dataset for SST in this study, we only show results obtained using HadISST data. We also used Twentieth Century Reanalysis Version 2 (20CRv2) data (Compo et al. 2011) from NOAA for heat flux, radiation and wind data, at a resolution of 2° × 2°. For all of these datasets, the period 1900-2017 are analyzed except the 20CRv2 data, which only covers the period 1900-2012.

Statistical methods and indices
In this paper, the NAO Index (NAOI) is calculated as the difference between the normalized zonal mean SLP at middle (35° N) and high latitudes (65° N) of the North Atlantic (80° W-30° E) following Li and Wang (2003). The AMO index is calculated following Enfield et al. (2001) as the area-weighted mean SST anomalies of the North Atlantic (0°-60° N, 7.5°-75° W). Annual IOSST is defined as the area-weighted mean SST of the Indian Ocean region (20° S-10° N, 60°-100° E). The PDO index (PDOI) is defined as the first leading principal component of North Pacific (north of 20° N) monthly sea surface temperature variability (Mantua et al. 1997). In exploring the multidecadal teleconnection between the Atlantic and Indian Ocean, the tropical Atlantic-Indian Ocean multidecadal teleconnection (AIMT) index is defined as follows: The analysis of this study focuses on multidecadal scale, so an 11-year running mean is applied to low-pass filtering. For determining the significance levels of correlation and regression, a two-tailed t-test is used. To consider the influence of low-pass filtering on degrees of freedom, the number of effective degrees of freedom (N eff ) is thus calculated as following approximation (Pyper and Peterman 1998;Li et al. 2013;Xie et al. 2019): where N is the sample size and XX (j) and YY (j) are the autocorrelations of two sampled time series X and Y at time lag j , respectively.
In research, the Hasselmann climate model was established to further illustrate the multidecadal change mechanism of the annual IOSST. For the integral solution of linear ordinary differential equations in the Hasselmann climate model, the classical fourth-order Runge-Kutta formula is used (Alexander 1977). The fourth-order Runge-Kutta calculation formula is: where y i+1 is the approximation value, h is integration stepsize and (x i , y i ) are initial values, k 1 , k 2 , k 3 and k 4 are the slopes at points (

Multidecadal variability in annual IOSST and relationship to NAO
Since the beginning of the twentieth century, the annual IOSST has shown a trend of warming with fluctuations and multidecadal variability, similar to the change in global temperature (Fig. 2a). The annual IOSST showed a trend of fluctuating warming between 1905 and 1940, reaching a maximum around 1940, and then gradual cooling to a minimum around 1950. From 1950 to around 2010, annual IOSST showed a trend of continuous warming with two turning points in the warming trend: around 1976 and 1987. The trend of warming from 1976 to 1987 was larger than that between 1950 and 1976. Since around 1987, the warming trend has slightly slowed down. This shows that although human activities in this period have contributed a lot to the IOSST (Sabine et al. 2004;Cowan et al. 2013;Roxy et al. 2015;Srivastava et al. 2016), there are still internal variability effects. After 2012, the 11-year running mean annual IOSST shows a downward trend. In the continuous power spectrum analysis as shown in Fig. 2b, we can find that the annual IOSST has a significant 70-year cycle. In addition, there are peaks over interannual scales at ~ 3 to ~ 4 years. Although the continuous power spectrum analysis only covers one and a half cycles, that the annual IOSST has a multidecadal cycle of 50-70 years has been found in several studies based on observed (Allan et al. 1995;Dong et al. 2016) and reconstructed (Damassa et al. 2006) data. Both the NAOI and annual IOSST anomalies show significant multidecadal variability. As shown in Fig. 3, the NAOI has a period of 50-70 years, consistent with some previous studies (e.g., Schlesinger 1994;Xie et al. 2019). The NAO leads the detrended annual IOSST in phase, but the relationship is not obvious (Fig. 4a). Studies have found that the IOSST can be driven by the PDO of the same period (Fig. 4c, d;Mantua and Hare 2002;Krishnan and Sugi 2003;Krishnamurthy and Krishnamurthy 2016). The phase relationship between the NAOI and annual IOSST with the  Lead-lag correlation between NAOI, PDOI and IOSST. a Lead-lag correlation between annual mean NAOI and detrended annual IOSST anomalies . The red (blue) line is for the annual mean (11-year running mean) time series. Negative (positive) lags indicate that the NAOI leads (lags) detrended annual IOSST, and the red (blue) dashed lines denote the 98% confidence levels for unsmoothed (smoothed) time series using the effective number of degrees of freedom. b As in (a), but for IOSST anomalies after removing the PDO signal. c As in (a), but for PDOI and IOSST. d As in (c), but for IOSST anomalies after removing the NAO signal PDO signal removed (IOSST_rm_PDO) is significant, in where the partial regression method is used to remove the PDO-related signals in the IOSST. This is also seen in the lead-lag correlation between the annual mean NAOI and the detrended annual IOSST_rm_PDO (Fig. 4b). The NAOI leads IOSST by 15-20 years on the multidecadal scale and the correlation coefficient reaches its maximum of 0.76 (significant at the 98% confidence level) when the NAOI leads the annual IOSST by 16 years (Fig. 4b). These results illustrate that the annual mean NAOI is ahead of the annual IOSST, and may be part of the annual IOSST multidecadal variability source.
After removing the PDO signal in the 11-year running mean NAOI with a lead time of 16 years has a significant correlation with the spatial field of the IOSST on the decadal scale. The correlation between the 11-year running mean NAOI and IOSST is generally significant in the Indian Ocean region, and the minimum value appears in the eastern equatorial region of the Indian Ocean region (Fig. 5a). After removing the PDO signal in the IOSST, the correlation coefficients increase and the significance is also substantially enhanced (Fig. 5b). The lead correlation coefficients exceed 0.8 in many areas and pass the 98% significance level. At the same, the 11-year running mean NAOI with a lead time of 16 years correlates well with the Arabian sea region SST multidecadal variability. This is consistent with the study by Sun et al. (2019). As similar results were obtained through multiple datasets for SST in these analyses, we only show results obtained using the HadISST data. This further illustrates that the NAO can be part of the annual IOSST multidecadal variability source.

Mechanism by which the NAO influences multidecadal variability in annual IOSST
It remains unclear how the NAO affects the annual IOSST with a phase lag of 16 years. According to NAO's delayed oscillator model, the NAO signal can be converted to the AMO (Sun et al. 2015). We propose a hypothesis that the NAO affects the annual IOSST via the AMO.

Relationship between the NAO and AMO
There have been many studies of the relationship between the NAO and AMO. The NAO can influence the North Atlantic SST (Hu and Huang 2006;Li et al. 2013;Sun et al. 2015;Stolpe et al. 2018). The NAO can drive deep ocean convection by affecting anomalous turbulent heat fluxes, and then drive the Atlantic meridional overturning circulation (AMOC). The AMOC then gives rise to SST anomalies as AMO-type (Sun et al. 2015;Stolpe et al. 2018). Some studies have found that the NAO leads the AMO by 15-20 years (Li et al. 2013;Xie et al. 2019). The correlation between the annual NAO and the annual SST variability of the North Atlantic is significant when the NAO leads by 16 years (Li et al. 2013). Wang et al. (2017) and Xie et al. (2019) further confirmed the relationship between the NAO and AMO based on CMIP5 historical simulation data.

Relationship between the AMO and IOSST and its interdecadal changes
The relationship between the PDO, AMO and annual IOSST shows interdecadal changes (Fig. 1). As shown in Fig. 6, the PDO and IOSST are relatively consistent in the period of 1900-1980. After 1980, the trends of the PDO and annual IOSST changed significantly. The PDO show a downward trend during the period of 1982-2008, while the annual IOSST had an upward trend during this period. The PDO showed a significant upward trend during the period of 2008-2017, while the annual IOSST has a downward trend during this period. In these periods when the relationship between the PDO and annual IOSST out of phase, the AMO and annual IOSST have good consistency, and the contributions of the AMO and PDO to annual IOSST are complementary. In other words, the IOSST shows three quasi-steady states in the periods of 1905-1940, 1950-1975 and 1985-2015, and two quick transitions in the periods of 1940-1950 and 1975-1985. The first one triggered by PDO alone but the second one by both. The first one is triggered by the PDO, while the second one is triggered by both the PDO and AMO, with the AMO having a greater impact. In the sliding correlation between the PDO, AMO, and IOSST ( Fig. 1), this phenomenon is more clear. The PDO and annual IOSST are significantly correlated during the period 1920-1990. The AMO and annual IOSST are significantly correlated during the periods 1915-1933, 1940-1954and 1980-2012. During 1934-1939and 1955-1975, the correlations between the AMO and annual IOSST are not significant, but the correlations between the PDO and annual IOSST are significant. During 1990-2012, the correlations between the PDO and annual IOSST are not significant, but the correlations between the AMO and annual IOSST are significant. These results indicate that rather than just the PDO driving the annual IOSST variability, the AMO can also have a driving effect on the annual IOSST variability, which means that the AMO and PDO may jointly drive the annual IOSST variability.

Mechanism by which the AMO affects IOSST multidecadal variability
In order to explore the physical processes that drive the multidecadal variability in IOSST, the surface heat budget over the Indian Ocean is analyzed by regressing net surface heat flux and radiation flux components onto the decadal IOSST and AMO index. As seen from Fig. 7, the warming in IOSST is driven by positive net longwave radiation into the ocean, whereas anomalous sensible heat flux, latent heat flux and net shortwave radiation contribute less to the IOSST. As seen from Fig. 8, the net longwave radiation in the Indian Ocean is related to the AMO, whereas the AMO is not strongly correlated with the sensible heat flux, the latent heat flux nor the net shortwave radiation. These indicate that the AMO may drive the IOSST multidecadal variability via net longwave radiation in the Indian Ocean. The correlation in some regions is not very strong, which is consistent with that the relationship between the AMO and IOSST has interdecadal changes. This result can also be obtained in time series of sensible heat flux, latent heat flux, net shortwave radiation and net longwave radiation in the Indian Ocean region. As shown in Fig. 9, the changes in the annual IOSST are in good agreement with the net longwave radiation in the Indian Ocean. For example, the volatility increased before 1940, the rapid decline between 1940 and 1951, the volatility increased between 1952 and 1962, the volatility decreased between 1963 and 1975, and the consistency of the two increased after 1975. For other heat fluxes, the consistency of the sensible heat flux in Indian Ocean region with the annual IOSST is good before 1960, but the sensible heat flux in Indian Ocean region is on the upward trend from 1961 to 1975, which is inconsistent with the downward trend of the annual IOSST during this period, and since 1987, the sensible heat flux is on a downward trend, while the annual IOSST is a fluctuating warming trend. The latent heat flux in Indian Ocean region was in good agreement with the annual IOSST before 1987, but after 1987, the latent heat flux showed a downward trend, which was inconsistent with the fluctuation of the annual IOSST during this period. This also corresponds to the results of the Fig. 7. For net shortwave radiation in Indian Ocean region, its volatility decreased from 1910 to 1940, which is inconsistent with the warming trend of the annual IOSST during this period, and from 1952 to 1968, the trend of the net shortwave radiation in Indian Ocean region first decreased and then increased, which was inconsistent with the trend of the annual IOSST warming first and then cooling during this period. These correspond to the previous results of heat fluxes regressing to the annual IOSST, indicating that the annual IOSST is mainly affected by net long wave radiation in Indian Ocean region. For net longwave radiation in Indian Ocean region and AMO index, the two are generally in good agreement, with a consistent decline before 1910, a consistent increase between 1910 and 1940, and the consistent trend of first volatility decreased after 1961 and then continued to rise. These also correspond to the relationship between the annual IOSST and the AMO index, which has interdecadal changes. This further shows that the annual IOSST is affected by net longwave radiation in Indian Ocean region, and the net longwave radiation in the Indian Ocean region may be affected by the AMO.
To further explore the relationship between the AMO and annual IOSST, we regress the surface heat flux and radiation flux onto the AMO for the period 1955-1975 when the annual IOSST is not closely related to the AMO (Fig. 10), and the period 1980-2012 when annual IOSST is significantly correlated with the AMO (Fig. 11). The net longwave radiation in the Indian Ocean region is strongly correlated with AMO in the period 1980-2012. The correlation for 1955-1975 is significantly different from that in 1980-2012. These results are consistent with interdecadal changes in the relationship between the AMO and annual IOSST. When the AMO has strong relationship with the net longwave radiation in the Indian Ocean region, the AMO is also strongly correlated with the annual IOSST. When the AMO has weak relationship with the net longwave radiation in the Indian Ocean region, the correlation between the AMO and the annul IOSST is weak. These further confirm that the AMO can drive the IOSST multidecadal variability via net longwave radiation in the Indian Ocean.
To investigate how the AMO affects net longwave radiation in the Indian Ocean region, we explore the differences heat flux is defined to be positive downward. The dotted area denotes regressions significant at the 90% confidence level. b-d As in (a) but for latent heat flux, shortwave radiation and longwave radiation, respectively between 1980-2012 and 1955-1975 in meridional wind, zonal wind and vertical wind anomaly averaged between 20° S and 10° N (Fig. 12). During the period when the annual IOSST and the AMO have a strong relationship, the vertical wind in the Atlantic region is abnormally enhanced, as is the vertical wind in the Indian Ocean region. At the same time there are two nearly parallel wave trains from the Atlantic Ocean to the Indian Ocean. Further analysis of this wave train shows that it has significant vertical structure and is defined as the tropical Atlantic-Indian Ocean multidecadal teleconnection (AIMT; Fig. 13). From the Atlantic Ocean to the Indian Ocean, it presents three positive and two negative structures, and the positive centers of the AIMT are distributed over the Atlantic Ocean and the Indian Ocean. Previous studies have found that the warm AMO phase can cause the tropical Atlantic anomalous vertical winds in West Africa, forming anomalous upper-level divergence and exciting Rossby wave (Watanabe 2004;Hodson et al. 2010;Alexander et al. 2014;Li et al. 2016). Then, enhance the vertical wind anomaly in the Indian Ocean region through the AIMT, which is consistent with the wind field over the Indian Ocean region during the warm AMO phase (Fig. 12) and the previous research that the wave source excited by the warm AMO phase can affect the wind field in the downstream Indian Ocean region . The vertical wind anomaly can increase lower atmospheric humidity and clouds. As shown in Fig. 14, in the high correlation period between the annual IOSST and AMO, the AMO can significantly enhance the relative humidity, precipitable water content, cloud water content and low cloud cover anomalies in the Indian Ocean. In the low correlation period between the annual IOSST and AMO, the correlations between the AMO and the relative humidity, precise water content, cloud water content and low cloud cover in the Indian Ocean are significantly different from those in high correlation period between the annual IOSST and AMO. In conclusion, these results suggest that the AMO may enhance the vertical wind anomaly in the Indian Ocean region through the AIMT, affecting the net longwave radiation by enhancing lower atmospheric lower atmospheric humidity and clouds in the Indian Ocean region, and thus affecting the IOSST.
The SLP field over the entire period 1900-2017 in the Atlantic and Indian Ocean regions is significantly negatively correlated with the AMO (Fig. 15a). This is consistent with the AMO corresponding to low pressure over the Atlantic Ocean region. At the same time, the corresponding low pressure in the Indian Ocean region indicates the presence of updrafts in the Indian Ocean region, which is consistent with results of previous studies of vertical wind anomalies and net longwave radiation. Compared with the period 1955-1975, the AMO positive position corresponds to the North Atlantic low pressure in the period 1980-2012, which is also consistent with the low pressure anomaly in the Indian Ocean region (Fig. 15b). These results further validate the impact of the AMO on the annual IOSST multidecadal variability.
Based on these findings, we propose that the AMO first enhances vertical wind anomalies in the Indian Ocean region through the AIMT, thereby increasing the net longwave radiation by enhancing lower atmospheric lower atmospheric humidity and clouds in the Indian Ocean region, which in Fig. 10 Regression patterns of surface sensible heat flux, latent heat flux, shortwave radiation and longwave radiation onto AMO. As in Fig. 7, but for the period 1955-1975 Fig. 11 Regression patterns of surface sensible heat flux, latent heat flux, shortwave radiation and longwave radiation onto AMO. As in Fig. 7, but for the period 1980-2012  1980-2012and 1955-1975. Differences between 1980-2012and 1955-1975 in the annual mean climatologies of the zonal-vertical wind circulation (vectors, m s −1 ) and meridional wind (shading, m s −1 ) averaged over 20° S-10° N turn drives the IOSST. When these relationships are strong, the AMO contributes more to the IOSST. When these relationships are weak, the AMO's contribution to the IOSST is also reduced.

Multidecadal climate model of annual IOSST
To further illustrate the multidecadal change mechanism of the annual IOSST, analogous to the approach proposed by Hasselmann (1976), we model the annual IOSST response to NAO and PDO variability (T NAO+PDO ) following the Hasselmann climate model: where T NAO+PDO denotes the simulated response of the annual IOSST to annual NAO and PDO variability, the estimated effective heat capacity C of climate system is between 13 and 35 W year m −2 K −1 (Knutti et al. 2008), and β is the linear damping coefficient, in the range from 0.3 to 1.2 K (W m −2 ) −1 (Foster et al. 2008). Based on the previous studies, the effective heat capacity actually represents the thermal inertia of the Earth's climate system, which is broadly defined as a long-scale ocean process (Foster et al. 2008;Knutti et al. 2008). Here, C is set to 25 W year m -2 K -1 . Based on the difference between the NAO and PDO for the IOSST forcing, β 1 and β 2 are set to 0.7 K (W m −2 ) −1 and 0.3 K (W m −2 ) −1 , respectively. In fact, as long as the values falls within a physically reasonable range, the results are not sensitive to the selection of β or C. The parameter α 1 and α 2 is determined empirically so that the root-mean-squared error between T NAO+PDO and the NAO and PDO-related IOSST is minimized. In this study, α 1 and α 2 are 0.6 and 0.9, respectively. This also corresponds to the correlation coefficients between the NAO, PDO and IOSST. The simulated response of the annual IOSST to the unfiltered original year-to-year NAO and PDO variability is shown in Fig. 16a. T NAO+PDO has similar phase and amplitude to the smoothed IOSST (Fig. 2) and the NAO and PDO-related IOSST which is obtained by the linear regression of the annual mean IOSST onto the NAO and PDO index (Fig. 16a). The correlation coefficient between T NAO+PDO and the NAO and PDO-related IOSST is 0.77 and pass the 98% significance level (t-test), providing further evidence that the multidecadal variability of the IOSST has received the joint effect of the NAO and PDO.
To further demonstrate the importance of the NAO for IOSST multidecadal variability, we construct two linear models, one relating IOSST to the PDO and the other to both the PDO and NAO. As shown in Fig. 17a, a PDObased linear climate model for the annual IOSST is established as follows: where t is time in years and the coefficients a 1 (0.18), b 1 (0.006), and c 1 (− 11.94) are determined empirically by linear regression based on the observed data, so that the root-mean-squared error of the model results and the observations is minimized.
(5) IÕSST(t) = a 1 PDO(t) + b 1 t + c 1 , Fig. 13 Meridional wind regressed onto the AIMT index. Regression patterns of annual meridional wind (shading, m s −1 ) at 20° S-10° N onto the annual AAIT index at decadal timescales . Dots denotes the regressions significant at the 95% confidence level, using the effective number of degrees of freedom This PDO-based decadal annual IOSST climate model reflects some of the characteristics of the IOSST. The correlation between the observations and the model's IOSST is 0.94, and the root-mean-squared error is 0.077. However, it is obvious that the model is better for the IOSST in the early stage (before 1950), but the simulation of the subsequent IOSST is not very good (after 1950). This is consistent with previous results indicating that the relationship between the PDO and IOSST is significant in early stages, but is not significant in later stages (Fig. 1). In addition, this can also be seen from Fig. 6 that from 1980 to 2010, the PDOI shows a downward trend, while the annual IOSST shows a warming trend. Therefore, an annual IOSST climate model that also incorporates the 16-year leading NAO signal is established as follows (Fig. 17b): where the coefficients a 2 (0.15), b 2 (0.22), c 2 (0.007) and d 2 (− 14.31) are determined using the same method as the previous Fig. 17a. Moreover, combined with the amplitudes of the NAOI and PDOI, it can be seen from the coefficients a 2 (0.15) and b 2 (0.22) that the PDO and NAO have equivalent contributions to the annual IOSST in this annual IOSST climate model. As shown in Fig. 17b, after adding the NAO signal, the correlation coefficient between the model and the observed IOSST increases to 0.98, and the root-meansquared error decreases to 0.041. Furthermore, the warming trend of the annual IOSST between 1985 and 2010 can be well simulated. The model (Formula 6) is able to simulate the annual mean IOSST well over the entire period.
In order to analyze the influence of the linear trend on the results of the model, the models are constructed without trends, respectively. The formulas are as follows: where the IÕSST* denotes the model IOSST constructed after detrended and the coefficients a 11   as the previous Fig. 17a. As shown in Fig. 17c and d, the model (Eq. 8) with the NAO signal added is better than that relies solely on the PDO. The coefficient a 11 (0.18) has little change compared with a 1 , and b 11 (0.000026) is very small, which indicates that although there is a linear trend term in the model (Eq. 7), it can be almost ignored. Also in Fig. 17d, the coefficients a 22 (0.15) and b 22 (0.22) did not change compared with a 2 and b 2 , and the linear term c 22 (0.001) was very small. Thus, the empirical models (Eqs. 6, 8) that includes the NAO signal is better than the empirical models (Eqs. 5, 7) that relies solely on the PDO, respectively. This is consistent with previous results that the NAO may be part of the annual IOSST multidecadal variability source. This also indicates that the annual IOSST multidecadal variability is co-driven by the PDO and NAO.

Summary and discussion
We study here the changing characteristics of the annual IOSST and find that IOSST undergoes significant multidecadal variability. Based on previous research of the relationship between the annual IOSST and PDO, the PDO signal from the same period as the annual IOSST was removed. A clear lead-lag relationship became evident between the NAO and the annual IOSST with the PDO signal removed. The NAO leads the annual IOSST, with the PDO signal removed, by 16 years, as shown in a lead-lag correlation diagram (Fig. 4). This suggests that the NAO may be part of the annual IOSST multidecadal variability source, and results for the spatial fields confirm this hypothesis. After removing the PDO signal, the correlation between the 16-year leading NAO and the IOSST multidecadal variability strengthens.
Based on the NAO's delayed oscillator model (Sun et al. 2015), and the relationship between the NAO and AMO, we propose that the AMO acts as a "bridge" by which the NAO affects the annual IOSST. The relationship between the AMO and IOSST is studied. The relationship between the PDO, AMO and annual IOSST varies interdecadally. When the PDO is not significantly related to the annual IOSST, the relationship between the AMO and the annual IOSST is significant. When the relationship between the AMO and the annual IOSST is not significant, the PDO is significantly related to the annual IOSST. These results indicate that it may be more than the PDO alone driving the annual IOSST variability, but the AMO also has a driving effect on the annual IOSST, which means that the AMO and PDO jointly drive the IOSST.
The mechanism by which the AMO drives the annual IOSST multidecadal variability is also explored. We find that net longwave radiation in the Indian Ocean region can warm IOSST. The net longwave radiation in the Indian Ocean region is related to the AMO, and the relationship Fig. 14 Regression patterns of relative humidity, precipitable water content, cloud water content and low cloud cover onto AMO. a, c, e and g are regression patterns of relative humidity (%), precipitable water content (kg m −2 ), cloud water content (kg m −2 ) and low cloud cover (%) onto the normalized AMO index over decadal timescales for the low correlation period  between the IOSST and AMO based on 20CRv2 and HadISST dataset, respectively. The dotted area denotes regressions significant at the 90% confidence level. b, d, f and h As in (a), (c), (e) and (g), but for the high correlation period  between the IOSST and AMO, respectively ◂ between the Indian Ocean region net longwave radiation and the AMO has interdecadal changes. The differences in meridional wind, zonal wind and vertical wind anomaly over a 20° S-10° N meridional mean between the periods 1980-2012 and 1955-1975 and the vertical structure of the AIMT show that the AMO can enhance the vertical wind anomaly in the Indian Ocean region through the AIMT, as in the SLP field. This affects the net longwave radiation by enhancing lower atmospheric lower atmospheric humidity and clouds in the Indian Ocean region, and thus the IOSST.
A Hasselmann model based on NAO and PDO was constructed to further verify the joint influence of the NAO and PDO on the multidecadal variability of IOSST. An empirical model models for the annual-mean IOSST were also established, one based on the PDO and another incorporating both the PDO and the leading 16-year NAO signal. The empirical model that includes the NAO signal performs better in predicting the IOSST than does the empirical model that relies solely on the PDO. This further verifies that the NAO is part of the annual IOSST multidecadal variability source, and is also consistent with previous results that indicate that the annual IOSST multidecadal variability is co-driven by the PDO and NAO (or AMO).
Note that this work only considers the NAO as part of the annual IOSST multidecadal variability source and the mechanism by which the NAO influences the annual IOSST multidecadal variability. Sun et al. (2017) found that the AMO can drive the multidecadal sea temperature variability in the western tropical Pacific. Therefore, the Pacific Ocean's contribution to the annual IOSST multidecadal variability may also come from the AMO. In other words, the direct effect of the AMO on the annual IOSST multidecadal variability and the effect of Pacific Ocean SST multidecadal variability on the annual IOSST multidecadal variability may all originate from the NAO. This warrants further investigation in future work.

Fig. 15
Regression patterns of SLP onto AMO and the differences between two periods. a Regression patterns of SLP (hPa) onto the normalized annual AMO index over decadal timescales based on 20CRv2 and HadISST dataset . The dotted area denotes regressions significant at the 90% confidence level. b As in (a), but for the regression pattern differences between 1980-2012 and 1955-1975  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.