Abstract

Studies of decadal-to-multidecadal ocean subsurface temperature variability are fundamental to improving the understanding of low-frequency climate signals. The present study uses the Simple Ocean Data Assimilation (SODA) version 2.2.4 product for the period 1950–2007 to identify decadal modes of variability that characterize the upper Indo-Pacific Ocean temperature structure (5–466-m depth). An empirical orthogonal function (EOF) analysis of the 10-yr low-pass filtered temperature field applied across four depths shows that the dominant mode is characterized by a long-term temperature trend, with warming at the surface and cooling at the thermocline depth connecting the tropical western Pacific with the southern Indian Ocean via the Indonesian Seas. EOF analysis of the detrended 10-yr filtered temperature data and correlation analyses of the EOF time series with established large-scale climate indices identified the interdecadal Pacific oscillation as EOF1, the North Pacific Gyre Oscillation as EOF2, and the decadal component of El Niño Modoki as EOF3 (respectively, modes 2, 3, and 4 of the nondetrended data). EOF2 identifies the Atlantic multidecadal oscillation when the analysis is applied to sea surface temperature anomalies only, suggesting that the surface is forced dominantly by fluxes associated with global-scale weather patterns, while the subsurface is dominantly forced by internal dynamics of the Pacific Ocean. This paper demonstrates that the decadal-to-interdecadal temperature variability in SODA has a pronounced vertical extension through the upper ocean. The upper thermocline accounts for most of the variance in the analysis. These results reinforce the importance of examining the subsurface ocean in climate dynamics studies that seek to understand the ocean’s role.

1. Introduction

Understanding the decadal-to-multidecadal variability of the upper-ocean thermal structure is essential to improve our understanding of the ocean’s role in coupled ocean–atmosphere climate dynamics, including climate variability modes and their mechanistic separation from long-term multidecadal trends. Importantly, better understanding of decadal climate modes of variability, for which the upper-ocean temperature (and density) structure plays a critical role in defining the time scale of the ocean and climate dynamics, is central to improved climate predictability. For example, it is known that the Pacific decadal oscillation (PDO) and the Atlantic multidecadal oscillation (AMO) influence climate by affecting the spatial and temporal variability in multidecadal drought frequency over the United States of America (McCabe et al. 2004), and the phase of the interdecadal Pacific oscillation (IPO)/PDO provides the background that underpins the potential interannual predictability of Australian rainfall, temperature, river flow, and wheat yield (Power et al. 1999; Arblaster et al. 2002; Power et al. 2006). Moreover, a more recent ocean climate feature named the North Pacific Gyre Oscillation (NPGO) explains the variations in salinity, nutrients, and surface chlorophyll-a (Chl-a) in the northeast Pacific Ocean remarkably well and is suggested to be a better indicator of upwelling strength, nutrient fluxes, and potential for ecosystem change in the northeast Pacific Ocean than the PDO (Di Lorenzo et al. 2008). The decadal climate prediction of those patterns is likely to be of great importance in the future (Hurrell et al. 2010).

Few studies have investigated the large-scale decadal temperature changes that characterize the upper Indo-Pacific Ocean, likely because of the sparsity of long-term observations of the vertical structure (Trenary and Han 2008). We know that the IPO (or the related PDO) dominates the low-frequency variability of sea surface temperature in the Pacific Ocean (Mantua et al. 1997; Power et al. 1999; Folland et al. 2002); however, its vertical expression is not well described yet, although various ocean modeling studies have attempted to define its subsurface signature and potential utility for ocean and climate predictability (McGregor et al. 2004, 2007, 2009; Power and Colman 2006; Holbrook et al. 2011).

The present study explores decadal temperature variability through the upper layers of the Indo-Pacific Ocean down to 466-m depth, its vertical expression in the upper ocean, and its relationship with known large-scale modes of climate variability, which are otherwise typically expressed as atmospheric or sea surface temperature patterns arising from EOF analyses or relatively simple climate indices. We use ocean-state estimates from version 2.2.4 of the Simple Ocean Data Assimilation (SODA) for the period 1950–2007 that have global spatial coverage and extend vertically to 5374-m depth. The representativeness of SODA for ocean decadal climate variability studies has been examined separately, where SODA was evaluated against independent measurements of sea level and other products (Vargas-Hernandez et al. 2014). The present study aims to describe and quantify decadal modes of Indo-Pacific Ocean temperature variability in the SODA reanalysis data characterized at specific levels within the mixed layer and thermocline (Carton and Giese 2008; SODA/TAMU Research Group 2008). By characterizing the decadal ocean subsurface temperature modes of variability in SODA, this informs our understanding of the expression of known climate modes in the upper ocean and provides a platform against which unconstrained ocean and/or coupled climate models might be evaluated.

2. Data and methods

This study uses version 2.2.4 of SODA, which applies an ocean general circulation model based on the Parallel Ocean Program (POP) physics version 2.0.1 numerics (Smith et al. 1992, 2010). SODA was generated for the purpose of using available ocean observations to improve the reanalysis of ocean variables with a simple assimilation algorithm (Carton et al. 2000; Carton and Giese 2008). It produces monthly averaged values across the globe on a uniform 0.5° × 0.5° × 40-level grid spanning from 1871 to 2008 at the date these data were downloaded. Here, we use only the period 1950–2007, as the year 2008 was incomplete when the data were processed, and, prior to 1950, the observations assimilated into SODA are extremely sparse. Many observations in SODA are from mechanical bathythermograph (MBT) and expendable bathythermograph (XBT) vertical profiles, with the MBTs only really increasing in density through the 1950s, while XBT profile collection only commenced in earnest in the 1970s. Prior to this, vertical profiles of temperature collected in the ocean are very limited.

To constrain the ocean model, SODA uses all available hydrographic temperature and salinity observations from the World Ocean Database 2009 (WOD09; Boyer et al. 2009) and sea surface temperature (SST) observations from the release 2.5 of the International Comprehensive Ocean–Atmosphere Dataset (ICOADS) (Woodruff et al. 2011). The POP model in SODA is driven by surface wind stress from the second version of the Twentieth Century Reanalysis (20Crv2) Project (Compo et al. 2008, 2011) carried out by the U.S. National Oceanic and Atmospheric Administration (NOAA) and the Cooperative Institute for Research in Environmental Science (CIRES) at the University of Colorado. The SODA data were downloaded from the SODA/TAMU website (http://soda.tamu.edu/). The 20Crv2 project uses reconstructed surface pressure observations and SST records for the period 1871–2008 to generate hemispheric weather patterns with the quality of a modern 3-day numerical forecast (Compo et al. 2011).

The representativeness of ocean-state estimates from SODA version 2.2.4 data for studying Indo-Pacific Ocean decadal temperature and sea level variability over the period 1950–2007 has been evaluated separately (Vargas-Hernandez et al. 2014). The SODA data were evaluated against independent sea level anomalies from long-record tide gauges at Midway Island and Fremantle, reconstructed sea surface height anomalies, and sea surface height anomalies from TOPEX/Poseidon satellite altimeter observations at the decadal time scale. From their analysis, Vargas-Hernandez et al. (2014) are able to show the utility of SODA as a potentially useful tool to examine ocean decadal climate variability across the Indo-Pacific Ocean.

In this study we applied empirical orthogonal function (EOF) analysis to identify patterns in the 1950–2007 SODA data. The EOF analysis was performed on the 10-yr low-pass Butterworth-filtered (LP10; Butterworth 1930) annual temperature anomalies simultaneously across four depth levels: 5, 97, 268, and 466 m. The four layers were chosen at the following: 5-m depth to explore decadal SST variability and its relationship with known decadal modes of climate variability typically expressed as atmospheric or sea surface temperature; and at depths of 97, 268, and 466 m to capture the vertical expression of decadal temperature variability below the mixed layer and within the thermocline. Those specific depths are related to vertical levels provided by the SODA product.

We have used the Butterworth method (Butterworth 1930) together with the function from MATLAB to filter the data. The filtfilt function performs zero-phase digital filtering by processing the input data, x, in both the forward and reverse directions. After filtering the data in the forward direction, filtfilt reverses the filtered sequence and runs it back through the filter. The filter has the following characteristics: 1) zero-phase distortion; 2) a filter transfer function, which equals the squared magnitude of the original filter transfer function; and 3) a filter order that is double the order of the filter specified by the numerator coefficients b and denominator coefficients a of the filter. The function filtfilt minimizes start-up and ending transients by matching initial conditions. Basically, this technique runs the filter twice—one time forward and the other backward—to rectify the edges of the time series. Therefore, there is no need to exclude data at the beginning or end of the time series.

The SODA data were analyzed directly in the first instance (nondetrended EOF analysis) then detrended in the subsequent EOF analysis (detrended EOF analysis). Four separate EOF analyses were applied. These were undertaken on 1) nondetrended temperature anomalies across the four vertical levels to investigate the contribution of temperature trends connecting the upper Indo-Pacific Ocean; 2) detrended temperature anomalies across the four vertical levels in an attempt to connect the decadal and multidecadal modes of upper Indo-Pacific Ocean temperature variability through the water column and with known large-scale decadal climate signatures; 3) the nondetrended sea surface temperature (at 5-m depth; herein referred to as the SST); and 4) detrended SST. The last two analyses were done as a benchmark and validation of SODA surface temperatures against previous studies of SST and climate variability. We used the IPO index from the Hadley Centre Global Sea Ice and SST dataset (HadISST; Folland 2008) to compare with the results of the EOF analysis applied to SODA. The IPO, rather than PDO, time series was used for the analysis because the IPO is more representative of the entire Pacific Ocean (Power et al. 1999), while the PDO characterizes the North Pacific Ocean (inter-) decadal SST and climate variability (e.g., Mantua et al. 1997; Di Lorenzo et al. 2010).

In standard EOF analysis of the data covariance matrix, the units of the field (here, temperature) are formally carried by the principal component (PC) time series, while the EOF spatial patterns are dimensionless (Venegas 2001). Here, the PCs are also dimensionless, since the eigenvectors derived from the EOF analysis were scaled so that the norm of each is 1. Furthermore, the EOF analysis was applied to both the standardized and nonstandardized data. We found that both results were very similar. Nevertheless, the nonstandardized data provided clearer EOF patterns, and their PCs were better correlated with known climate indices. Hence, we only show the results for the nonstandardized data here.

To explore the relationship between the decadal temperature EOFs in SODA and large-scale climate signatures, correlation analyses were performed between the PC time series and established climate variability indices. A suite of climate indices was explored and utilized, including the following: IPO, El Niño Modoki (otherwise known as the central Pacific El Niño), NPGO, and AMO. A table is provided with all the correlations and significances across the various indices for the subsurface and surface analyses (Table 1).

Table 1.

Correlation coefficients and significance levels (in %) between the PC time series of the filtered and detrended temperature and climate variability indices. The results for the subsurface (PC1-D, PC2-D, and PC3-D) and surface (PC1-D5, PC2-D5, and PC3-D5) analyses are shown. Bold values indicate indices with the highest correlation coefficients.

Correlation coefficients and significance levels (in %) between the PC time series of the filtered and detrended temperature and climate variability indices. The results for the subsurface (PC1-D, PC2-D, and PC3-D) and surface (PC1-D5, PC2-D5, and PC3-D5) analyses are shown. Bold values indicate indices with the highest correlation coefficients.
Correlation coefficients and significance levels (in %) between the PC time series of the filtered and detrended temperature and climate variability indices. The results for the subsurface (PC1-D, PC2-D, and PC3-D) and surface (PC1-D5, PC2-D5, and PC3-D5) analyses are shown. Bold values indicate indices with the highest correlation coefficients.

El Niño Modoki is an important climate phenomenon of the Pacific Ocean, which has become more intense and frequent in the last decades (Lee and McPhaden 2010; Di Lorenzo et al. 2010). Studies of the El Niño Modoki have mostly focused on interannual time scales (Ashoket al. 2007; Ashok and Yamagata 2009; Li et al. 2010). Here, we explore the decadal component of El Niño Modoki using a time series of the improved El Niño Modoki index (IEMI) calculated from SODA according to Li et al.’s (2010) method for the period 1950–2007. The El Niño Modoki spatial pattern in SODA SST for the Indo-Pacific region (60°N–60°S, 20°E–60°W) was found by correlating the IEMI time series with the gridpoint SST time series across the entire domain. Both the Indo-Pacific El Niño Modoki SST pattern and the IEMI time series were compared with the results of the EOF analysis applied to SODA.

The NPGO, defined by Di Lorenzo et al. (2008) as the second PC of sea surface height anomalies (SSHA) in the North Pacific Ocean, is also an important decadal climate pattern of the Pacific Ocean, though it is still not very well understood. It seems to extend beyond the North Pacific to form part of a global-scale phenomenon (Di Lorenzo et al. 2008, 2010). In addition, decadal changes in the North Pacific Ocean are also influenced by the AMO through an atmospheric teleconnection (Zhang and Delworth 2007). Nevertheless, it is still not fully clear whether its influence also extends to the whole Indo-Pacific Ocean on vertical and horizontal scales. Here, we used the annual-averaged (based on calendar year) NPGO index generated by Di Lorenzo et al. (2008; 2015). We also used the unsmoothed annual-averaged (based on calendar year) AMO calculated from Kaplan’s SST dataset (NOAA 2015). As this paper concentrates on decadal time scales, we used the LP10 time series for all the climate indices correlated with the temperature PCs from the EOF analysis. Significances of the correlation coefficients and in all of the comparisons were assessed based on the effective number of degrees of freedom, taking account of serial correlation in the individual time series according to Davis (1976).

Power spectra of all the climate indices prior to LP10 filtering are shown in Fig. 1. The four unfiltered, but annually averaged, climate indices (IPO, NPGO, IEMI, and AMO) have significant low-frequency peak periodicities of >5 yr [values of log(f) < −0.6 in Fig. 1] when compared against the red noise spectrum. The peaks in IPO are strongest at interannual time scales, while the strongest peaks for the IEMI, NPGO, and AMO are at decadal-to-interdecadal time scales (Figs. 1b–d), at least for the annually averaged data. The IPO shows the maximum energy on interannual time scales likely associated with ENSO, since these two patterns have been shown to be linked (Di Lorenzo et al. 2010, 2013), although peaks on decadal-to-interdecadal time scales are also strong and significant (Fig. 1a).

Fig. 1.

Variance-preserving spectra of several large-scale climate indices: (a) IPO index, (b) IEMI, and (c),(d) NPGO and AMO indexes respectively. The spectra were calculated from the following record lengths: (a) 1901–2007 and (b)–(d) 1950–2007. Most of the low-frequency peaks are significant at the 95% level and above red noise. The periods with highest energy are indicated with numbers in years.

Fig. 1.

Variance-preserving spectra of several large-scale climate indices: (a) IPO index, (b) IEMI, and (c),(d) NPGO and AMO indexes respectively. The spectra were calculated from the following record lengths: (a) 1901–2007 and (b)–(d) 1950–2007. Most of the low-frequency peaks are significant at the 95% level and above red noise. The periods with highest energy are indicated with numbers in years.

3. Results

In the following subsections, we explain the main results for the subsurface EOF analysis performed on the LP10 annual temperature anomalies simultaneously across four depth levels: 5, 97, 268, and 466 m. These results include the spatial patterns of the EOFs, the PC time series compared with climate indices, and a coherence analysis performed for each mode and climate index.

a. Analysis 1: Nondetrended subsurface temperatures

Here, the EOF analysis was performed on the nondetrended subsurface temperatures. Only the first EOF is shown here because it captures the trend; the higher-order EOFs are represented equally well in the detrended time series.

EOF1-ND: Nondetrended subsurface temperature

The dominant EOF of the nondetrended temperature data (EOF1-ND) explains about 28% of the total upper-ocean temperature variance across the four vertical levels in the Indo-Pacific Ocean and is dominated by the long-term temperature trend (Fig. 2e). Broad-scale warming is identified at the 5-m level with cooling regions in the interior North Pacific subtropical gyre extending eastward from Japan and in the southern midlatitudes of the Indo-Pacific and Southern Ocean, defined within the domain of interest (Fig. 2a). At 97-m depth (Fig. 2b) there is a strong meridional gradient in the temperature anomaly pattern across the equator—from warming north of the equator to cooling south of the equator—which does not appear at the surface. At 268-m depth, we also observe a cooling pattern between the Pacific Ocean and the Indian Ocean, linked through the Indonesian Seas (Fig. 2c). In the Indian Ocean this cooling pattern is split into two paths: one across the interior Indian Ocean along the zonal band from ~10° to 19°S extending from western Australia to Madagascar and eastern Africa, and the other extending from the Indonesian Archipelago along the western Australian coastal waveguide and across southern Australia.

Fig. 2.

EOF1-ND (~28.6% of the variance of the LP10 filtered and nondetrended temperature anomalies) from analysis of the 1950–2007 SODA 2.2.4 data. The analysis is performed simultaneously across four vertical levels [(a) 5, (b) 97, (c) 268, and (d) 466 m] with the EOF1 spatial pattern loadings across the four levels shown in (a)–(d). (e) The corresponding PC1-ND. Units are nondimensional.

Fig. 2.

EOF1-ND (~28.6% of the variance of the LP10 filtered and nondetrended temperature anomalies) from analysis of the 1950–2007 SODA 2.2.4 data. The analysis is performed simultaneously across four vertical levels [(a) 5, (b) 97, (c) 268, and (d) 466 m] with the EOF1 spatial pattern loadings across the four levels shown in (a)–(d). (e) The corresponding PC1-ND. Units are nondimensional.

The 97-m level (Fig. 2b) has the highest variance of temperature across the four depths. As expected, much of the variance explained appears in the depth range of the thermocline where gradients of temperature are largest. We found that the second and higher-order significant EOFs from the ND analysis can be directly relatable to known large-scale climate modes of variability. Here, modes 2–4 closely correspond to the first few EOFs of the detrended analysis and are instead discussed in the following subsection.

b. Analysis 2: Detrended subsurface temperatures

In this analysis, the linear trend was removed from the data before applying the EOF. The first three leading EOFs are shown.

1) EOF1-D: Detrended subsurface temperature

The PC1-D time series from the dominant EOF of the detrended temperature (EOF1-D; Fig. 3e) is significantly correlated with the IPO index (r = 0.8; significant at the >99.6% level). A spectral analysis shows significant coherence at the >95% level with little phase lag (from −10° to 0°) between PC1-D and the IPO for a broad range of periods longer than about 15 years (Fig. 4).

Fig. 3.

EOF1-D (~20.3% of the variance of the LP10 filtered and detrended temperature anomalies) from analysis of the 1950–2007 SODA 2.2.4 data. The analysis is performed simultaneously across four vertical levels [(a) 5, (b) 97, (c) 268, and (d) 466 m] with the EOF1 spatial pattern loadings across the four levels shown in (a)–(d). (e) The PC1-D (black line) and the IPO index with LP10 filtered (blue line with circles) and the annual IPO index (red line with squares). Units are nondimensional.

Fig. 3.

EOF1-D (~20.3% of the variance of the LP10 filtered and detrended temperature anomalies) from analysis of the 1950–2007 SODA 2.2.4 data. The analysis is performed simultaneously across four vertical levels [(a) 5, (b) 97, (c) 268, and (d) 466 m] with the EOF1 spatial pattern loadings across the four levels shown in (a)–(d). (e) The PC1-D (black line) and the IPO index with LP10 filtered (blue line with circles) and the annual IPO index (red line with squares). Units are nondimensional.

Fig. 4.

Spectral analysis for the low frequencies of EOF1-D (LP10 filtered and detrended temperatures across the four vertical levels) and the IPO index. (a) Power spectra of EOF1-D; (b) coherence between EOF1-D and the IPO index (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF1-D and the IPO index; and (d) power spectra of the IPO index.

Fig. 4.

Spectral analysis for the low frequencies of EOF1-D (LP10 filtered and detrended temperatures across the four vertical levels) and the IPO index. (a) Power spectra of EOF1-D; (b) coherence between EOF1-D and the IPO index (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF1-D and the IPO index; and (d) power spectra of the IPO index.

EOF1-D explained 20.3% of the total temperature variance across the four Indo-Pacific Ocean depth levels. Most of the variance explained is at the 97-m level (Fig. 3b), with the IPO signature being expressed strongly and a pronounced symmetric pattern around the equator in the Pacific Ocean and cooling patterns in the western North and southern Pacific Ocean. The Indonesian Seas and the tropical south Indian Ocean also showed a strong signal at 97-m depth that appears to be connected to the western tropical Pacific Ocean.

At the 268-m level (Fig. 3c), the IPO signature is reduced in amplitude, and the spatial pattern is very different from the corresponding pattern at the surface, with most of its variance in the western Pacific. The 466-m level (Fig. 3d) shows a very weak IPO signal.

2) EOF2-D: Detrended subsurface temperature

North of the equator, EOF2-D is similar in spatial and temporal structure to the NPGO (Fig. 5). The PC2-D time series correlates significantly (r = −0.57; significant at the >97% level) with the LP10 NPGO index. Spectral analysis shows that the correlation arises at the faster time scales, where PC2-D and NPGO are coherent at the >95% significance level and at frequencies between 0.08 and 0.1 cycles per year (i.e., periods of 12.5–10 yr; Fig. 6b) but not the lower frequencies, which explains the relatively low correlation. The coherent band includes the highest energy peak (around 0.08 cycles per year or ~12.8-yr period) of the NPGO spectral power for the period 1950–2007 (Fig. 1b). A relationship between the NPGO and temperature south of the equator has not, to our knowledge, been identified in previous studies.

Fig. 5.

EOF2-D (~17% of the variance of the LP10 filtered and detrended temperature anomalies) from analysis of the 1950–2007 SODA 2.2.4 data. The analysis is performed simultaneously across four vertical levels [(a) 5, (b) 97, (c) 268, and (d) 466 m] with the EOF2 spatial pattern loadings across the four levels shown in (a)–(d). (e) The PC2-D (black line) and the NPGO index with LP10 filtered (blue line with circles) and the annual NPGO index (red line with squares). Units are nondimensional.

Fig. 5.

EOF2-D (~17% of the variance of the LP10 filtered and detrended temperature anomalies) from analysis of the 1950–2007 SODA 2.2.4 data. The analysis is performed simultaneously across four vertical levels [(a) 5, (b) 97, (c) 268, and (d) 466 m] with the EOF2 spatial pattern loadings across the four levels shown in (a)–(d). (e) The PC2-D (black line) and the NPGO index with LP10 filtered (blue line with circles) and the annual NPGO index (red line with squares). Units are nondimensional.

Fig. 6.

Spectral analysis for the low frequencies of EOF2-D (LP10 filtered and detrended temperatures across the four vertical levels) and the NPGO index. (a) Power spectra of EOF2-D; (b) coherence between EOF2-D and the NPGO index (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF2-D and the NPGO index; and (d) power spectra of the NPGO index.

Fig. 6.

Spectral analysis for the low frequencies of EOF2-D (LP10 filtered and detrended temperatures across the four vertical levels) and the NPGO index. (a) Power spectra of EOF2-D; (b) coherence between EOF2-D and the NPGO index (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF2-D and the NPGO index; and (d) power spectra of the NPGO index.

EOF2-D explains 17% of the total decadal-to-multidecadal-scale temperature variability across the four vertical levels analyzed. The largest percentage of the variance explained by EOF2-D is at the 97-m (Fig. 5b) and 268-m levels (Fig. 5c). EOF2-D also shows a strong cooling at the 97-m level (Fig. 5b) connecting the western equatorial Pacific with the Indian Ocean. It is interesting to see that the SST response of the eastern Indian Ocean is very weak (as for EOF1-D).

3) EOF3-D: Detrended subsurface temperature

The EOF3-D surface spatial pattern (Fig. 7a) and PC3-D time series (Fig. 7e) clearly show the El Niño Modoki signal. The Modoki signature is apparent and strong at 97- and 268-m depth (Figs. 7b,c), but it does not appear to extend strongly to 466-m depth (Figs. 7d). We see that the highest energy spectral peak corresponds to a period of 12.8 yr (significant at the 95% level, Fig. 1c). This decadal-scale modulation of the El Niño Modoki signature appears to be captured well by EOF3-D (Fig. 7). Further, PC3-D is strongly and significantly correlated (r = 0.8, significant at the >99.9% level) with the IEMI (Fig. 7e). Moreover, spectral analysis indicates that the PC3-D is coherent with the IEMI for all frequencies <0.1 cycles per year (i.e., >10-yr period; significant at the >95% level; Fig. 8b).

Fig. 7.

EOF3-D (~11% of the variance of the LP10 filtered and detrended temperature anomalies) from analysis of the 1950–2007 SODA 2.2.4 data. The analysis is performed simultaneously across four vertical levels [(a) 5, (b) 97, (c) 268, and (d) 466 m] with the EOF2 spatial pattern loadings across the four levels shown in (a)–(d). (e) The PC3-D (black line) and the IEMI with LP10 filtered (blue line with circles) and the annual IEMI (red line with squares). Units are nondimensional.

Fig. 7.

EOF3-D (~11% of the variance of the LP10 filtered and detrended temperature anomalies) from analysis of the 1950–2007 SODA 2.2.4 data. The analysis is performed simultaneously across four vertical levels [(a) 5, (b) 97, (c) 268, and (d) 466 m] with the EOF2 spatial pattern loadings across the four levels shown in (a)–(d). (e) The PC3-D (black line) and the IEMI with LP10 filtered (blue line with circles) and the annual IEMI (red line with squares). Units are nondimensional.

Fig. 8.

Spectral analysis for the low frequencies of EOF3-D (LP10 filtered and detrended temperatures across the four vertical levels) and the IEMI. (a) Power spectra of EOF3-D; (b) coherence between EOF3-D and the IEMI (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF3-D and the IEMI; and (d) power spectra of the IEMI.

Fig. 8.

Spectral analysis for the low frequencies of EOF3-D (LP10 filtered and detrended temperatures across the four vertical levels) and the IEMI. (a) Power spectra of EOF3-D; (b) coherence between EOF3-D and the IEMI (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF3-D and the IEMI; and (d) power spectra of the IEMI.

EOF3-D explains 11% of the total decadal temperature variance of the Indo-Pacific Ocean across the four vertical levels. As before, the 97-m level corresponds to the highest variance explained across the four depth levels analyzed. The strength of the El Niño Modoki signature increased in the tropical Pacific over the last three decades or so (Ashok et al. 2007; Lee and McPhaden 2010). We note here the strongest correlation between PC3-D and the IEMI is evidenced from the 1970s onward (Fig. 7e). The strength of the similarity in spatial pattern of El Niño Modoki and EOF3-D is demonstrated in Fig. 9. A high percentage (84.9%) of the Indo-Pacific Ocean EOF3-D pattern loadings correspond in space (yellow) with the Modoki mode, and only a few regions (15.1%; red pattern) show discrepancies that are mostly localized in the west and south Indian Ocean, west Pacific Ocean, and southeast Pacific Ocean.

Fig. 9.

Spatial pattern of the decadal component of El Niño Modoki (a) generated by correlating the IEMI time series (Li et al. 2010) with the LP10 filtered and detrended annual sea surface temperature anomalies. (b) The spatial EOF3-D spatial pattern loadings, as in Fig. 7a (units are nondimensional). (c) The comparison between (a) and (b); the yellow color shows the regions where the patterns have the same sign (S; either positive or negative) and the red color shows where the patterns have opposite sign (D). Units are nondimensional.

Fig. 9.

Spatial pattern of the decadal component of El Niño Modoki (a) generated by correlating the IEMI time series (Li et al. 2010) with the LP10 filtered and detrended annual sea surface temperature anomalies. (b) The spatial EOF3-D spatial pattern loadings, as in Fig. 7a (units are nondimensional). (c) The comparison between (a) and (b); the yellow color shows the regions where the patterns have the same sign (S; either positive or negative) and the red color shows where the patterns have opposite sign (D). Units are nondimensional.

We have also correlated the NPGO index with PC3-D, given that Di Lorenzo et al. (2010) argue that El Niño Modoki (also known as central Pacific El Niño) forces an atmospheric response in the North Pacific, which in turn drives decadal variability of the NPGO. The correlation between the NPGO index and PC3-D at zero lag is r = −0.10 (insignificant), whereas the correlation between the IEMI and PC3-D at zero lag is r = 0.80 (significant at the >99.9% level, Table 1). We have also performed lagged correlation analysis between the NPGO index and PC3-D and found that the maximum correlation was r = 0.57 (significant at the >95%) at the 4-yr lag (PC3-D lagging the NPGO index). This suggests that PC3-D and the NPGO are related at some level, albeit that the strongest and most highly significant relationship is with the IEMI at zero lag.

c. Analysis 3: Nondetrended SSTA

We also performed an EOF analysis on the sea surface temperature anomalies (SSTA; Analyses 3 and 4), identified here as the 5-m level temperature anomalies in SODA, to explore whether there are any salient differences with respect to the surface pattern against the multidepth analysis. Large-scale climate indices were again used as per the previous section to identify relationships between the SSTA EOF modes (PC time series) and these proxies for known climate forcing.

EOF1-ND5: Nondetrended SSTA

EOF1-ND5 of the LP10 nondetrended SSTA (i.e., EOF1 of the 5-m depth data; Fig. 10) explained 32.7% of the total decadal-to-multidecadal SSTA variance and defines the long-term trend, as in the multidepth Analysis 1. However, the decadal modulation of the trend is different, in particular during the last part of the record from 2003 onward, where EOF1-ND (Fig. 2e) shows a minimum and then increasing tendency, whereas EOF1-ND5 simply continues to decrease over the same period. Analysis 1 fits the subsurface as well and explains subsurface variances, whereas Analysis 3 only considers SSTA variances explained by the EOF analysis. One explanation for the differences between the trends could be that the stored heat from subsurface levels has contributed to maintain the increasing tendency in subsurface temperature from 2003 onward. On the other hand, positive heat fluxes from the ocean to the atmosphere might be influencing the decrease in SSTA in those last years of the record.

Fig. 10.

EOF1-ND5 (~32.7% of the variance of the LP10 filtered and nondetrended SSTA) from analysis of the 1950–2007 SODA 2.2.4 data. The EOF is performed on the LP10 annual SSTA (at 5-m depth). (a) The spatial pattern loadings and (b) the PC1-ND5 are shown. Units are nondimensional.

Fig. 10.

EOF1-ND5 (~32.7% of the variance of the LP10 filtered and nondetrended SSTA) from analysis of the 1950–2007 SODA 2.2.4 data. The EOF is performed on the LP10 annual SSTA (at 5-m depth). (a) The spatial pattern loadings and (b) the PC1-ND5 are shown. Units are nondimensional.

d. Analysis 4: Detrended SSTA

1) EOF1-D5: Detrended SSTA

EOF1-D5 of the detrended SSTA (i.e., EOF1 of the 5-m depth data; Fig. 11a) explains 28% of the total decadal-to-multidecadal-scale variance of the detrended and LP10 SSTA in the Indo-Pacific Ocean. EOF1-D5 effectively defines the IPO SSTA mode signature for the Indo-Pacific Ocean domain (Fig. 12a). The correlation coefficient between PC1-D5 and the IPO index is large and highly significant (r = 0.98; significant at the >99.9% level), and the coherence between PC1-D5 and the IPO index is significant at the >95% level across the entire range of low frequencies (Fig. 13b).

Fig. 11.

Spatial patterns of (a) EOF1-D5, (b) EOF2-D5, and (c) EOF3-D5 of the LP10 filtered and detrended SSTA. Units are nondimensional.

Fig. 11.

Spatial patterns of (a) EOF1-D5, (b) EOF2-D5, and (c) EOF3-D5 of the LP10 filtered and detrended SSTA. Units are nondimensional.

Fig. 12.

Time series of (a) PC1-D5, (b) PC2-D5, and (c) PC3-D5 of the LP10 filtered and detrended SSTA. PC1-D5 is highly correlated with the IPO index, PC2-D5 with the AMO index, and PC3-D5 with the IEMI. Units are nondimensional.

Fig. 12.

Time series of (a) PC1-D5, (b) PC2-D5, and (c) PC3-D5 of the LP10 filtered and detrended SSTA. PC1-D5 is highly correlated with the IPO index, PC2-D5 with the AMO index, and PC3-D5 with the IEMI. Units are nondimensional.

Fig. 13.

Spectral analysis for the low frequencies of the EOF1-D5 (LP10 filtered and detrended SSTA) and the IPO index. (a) Power spectra of EOF1-D5; (b) coherence between EOF1-D5 and the IPO index (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF1-D5 and the IPO index; and (d) power spectra of the IPO index.

Fig. 13.

Spectral analysis for the low frequencies of the EOF1-D5 (LP10 filtered and detrended SSTA) and the IPO index. (a) Power spectra of EOF1-D5; (b) coherence between EOF1-D5 and the IPO index (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF1-D5 and the IPO index; and (d) power spectra of the IPO index.

2) EOF2-D5: Detrended SSTA

EOF2-D5 explains 19% of the total decadal-to-multidecadal variance of the detrended LP10 SSTA in the Indo-Pacific Ocean. We found that PC2-D5 (Figs. 11b and 12b) correlates more strongly with the AMO index (r = −0.87; significant at the >97.9% level) than with the NPGO index (r = 0.13; significant at the >40% level; see Table 1). Spectral analysis shows that PC2-D5 is coherent with the AMO index and significant at the >95% level for most of the frequencies <0.08 cycles per year (i.e., periods >12.5 yr), and both time series are around 170° out of phase (Fig. 14c).

Fig. 14.

Spectral analysis for the low frequencies of the EOF2-D5 (LP10 filtered and detrended SSTA) and the AMO index. (a) Power spectra of EOF2-D5; (b) coherence between EOF2-D5 and the AMO index (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF2-D5 and the AMO index; and (d) power spectra of the AMO index.

Fig. 14.

Spectral analysis for the low frequencies of the EOF2-D5 (LP10 filtered and detrended SSTA) and the AMO index. (a) Power spectra of EOF2-D5; (b) coherence between EOF2-D5 and the AMO index (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF2-D5 and the AMO index; and (d) power spectra of the AMO index.

EOF2-D5 has a meridional asymmetry in the Pacific Ocean. The North and central Pacific show negative anomalies in most of the locations (Fig. 11b), whereas the Southern Ocean shows positive anomalies, specifically south of Tasmania and New Zealand. EOF2-D5 has a weak signature in the Indian Ocean.

3) EOF3-D5: Detrended SSTA

EOF2-D5 explains 11% of the total decadal-to-multidecadal variance of the detrended LP10 SSTA in the Indo-Pacific Ocean. EOF3-D5 is very similar to EOF3-D (over four vertical levels) and is highly significantly correlated with El Niño Modoki (Figs. 11c and 12c). The correlation between the PC3-D5 and the IEMI time series is r = 0.75, significant at the >99.5% level (Fig. 11c). The time series are coherent at the 95% significance level across the decadal spectrum (Fig. 15b). Interestingly, we note that the phase relationship between the IEMI and mode 3 is stronger in the multilevel analysis case with PC3-D (r = 0.8; significant at the 99.9% level) than with PC3-D5 (r = 0.75; significant at the 99.9% level)—although their significance levels are comparable—even though the IEMI has been calculated from only SST. This demonstrates the strong connectivity of El Niño Modoki through the upper ocean.

Fig. 15.

Spectral analysis for the low frequencies of the EOF3-D5 (LP10 filtered and detrended SSTA) and the IEMI. (a) Power spectra of EOF3-D5; (b) coherence between EOF3-D5 and the IEMI (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF3-D5 and the IEMI; and (d) power spectra of the IEMI.

Fig. 15.

Spectral analysis for the low frequencies of the EOF3-D5 (LP10 filtered and detrended SSTA) and the IEMI. (a) Power spectra of EOF3-D5; (b) coherence between EOF3-D5 and the IEMI (dashed line is the 95% significance level corresponding to the red noise spectrum); (c) phase relationship between EOF3-D5 and the IEMI; and (d) power spectra of the IEMI.

4) Sampling errors in the EOF analysis

We have performed the North et al. (1982) test on the detrended decadal temperature anomalies and found that EOF1-D and EOF2-D are not separable (Fig. 16a). This is consistent with the similarities evident in the spatial patterns of EOF1-D and EOF2-D. It is notable, however, that Table 1 shows that the correlation of PC2-D with the NPGO index is higher and more significant than with the IPO. Conversely, PC1-D is not significantly correlated with the NPGO index but is instead highly correlated with the IPO index (Table 1). EOF3-D, associated here with the decadal component of El Niño Modoki, is clearly independent from EOF1-D and EOF2-D (Fig. 16a).

Fig. 16.

North et al. (1982) test for (a) the EOF analysis performed on the LP10 annual temperature anomalies simultaneously across four vertical levels (5, 97, 268, and 466 m); and (b) the EOF analysis performed on the LP10 annual SSTA (at 5-m depth).

Fig. 16.

North et al. (1982) test for (a) the EOF analysis performed on the LP10 annual temperature anomalies simultaneously across four vertical levels (5, 97, 268, and 466 m); and (b) the EOF analysis performed on the LP10 annual SSTA (at 5-m depth).

Interestingly, the North et al. (1982) test applied to the EOF analysis of the detrended decadally filtered sea surface temperature anomalies (5-m depth) shows that the first three EOFs are clearly independent. In this case, the corresponding PC time series for EOF1-D5, EOF2-D5, and EOF3-D5 characterize well indices of the IPO, AMO, and decadal component of El Niño Modoki, respectively, with Fig. 16b showing that those modes are clearly separated.

There appears to be a relationship between the IPO and NPGO, since PC1-D5 (first EOF of SSTA) is correlated with both the IPO and NPGO (Table 1), although PC1-D5 correlates better and more significantly with the IPO. Di Lorenzo et al. (2008, 2010) show that the IPO and the NPGO are separate EOFs for the North Pacific Ocean. Here, a much larger domain across the Indo-Pacific Ocean is included in the EOF analysis, which could generate more mixed results and explain the relationship between PC1-D5 (here associated with the IPO) and NPGO (Table 1). On the other hand, PC2-D (second EOF of multilevel analysis) is significantly correlated with the NPGO index but is also correlated with the IPO index, which is significant at the 93% level (Table 1).

One possible explanation for the apparent coupling between EOF1-D and EOF2-D in the multilevel EOF analysis applied here is that the North Pacific Gyre Oscillation has a strong expression through the vertical structure, and its spatial pattern at the ocean surface resembles the IPO when a larger domain across the Indo-Pacific Ocean is considered in the EOF analysis. It is important to keep in mind that the NPGO index used here was constructed using sea surface heights and therefore considers the upper-ocean variability of the North Pacific Ocean, whereas the IPO index is based on sea surface temperatures alone.

4. Discussion and conclusions

We found that the dominant mode of decadal-to-interdecadal temperature variability through the upper 450 m of the Indo-Pacific Ocean from 1950 to 2007 is identified as an increasing temperature trend. The spatial pattern of this temperature trend (Fig. 2a) at the sea surface is very similar to that shown in Ridgway (2007) for the period 1944–2005, using the Extended Reconstructed SST (ERSST) data of Smith and Reynolds (2003). Clearly, both the dominant SSTA (trend) mode using the SODA data and the trend pattern from the ERSST data show a warming trend at the sea surface across most of the Indo-Pacific Ocean and a strong cooling trend in the western North Pacific Ocean. This cooling pattern in the western North Pacific Ocean has a PDO-like shape. Meehl et al. (2009) argue that a pattern in SST similar to this one, based on results from a twentieth-century global climate model simulation, represents a superposition of the IPO and external forcing caused by increasing greenhouse gases.

We found pronounced differences in the patterns of temperature trends across different depths, in particular in the Indian Ocean. Our results show that the Indian Ocean has undergone a very homogeneous warming trend at the surface but has cooling patterns below the surface that are consistent with other reported observational studies by Alory et al. (2007), Trenary and Han (2008), and Alory and Meyers (2009)(Fig. 2). The trend mode from our analysis also shows a connection between the tropical west Pacific and the south Indian Ocean at 268-m depth, which splits into two paths that correspond to the confluence of the Indo-Pacific waveguides, as described by Wijffels and Meyers (2004).

One pathway extends from east to west across the Indian Ocean in the zonal band between ~10° and 19°S and then southward between mainland Africa and Madagascar, following a similar path to the Agulhas Current (Fig. 2c). This multidecadal-scale pathway identified here is consistent with the pathway of Trenary and Han (2008), who argue that remote forcing from the Pacific Ocean via the Indonesian Seas contributes to the cooling trend below 200-m depth in the south Indian Ocean. This also agrees with other similar findings reported by Wainwright et al. (2008), Timmermann et al. (2010), Schwarzkopf and Böning (2011), and Nidheesh et al. (2013), who suggest that there is indeed a dynamic connection in the trend between the western Pacific Ocean and Indian Ocean, based largely on ocean heat content and sea level data.

The other pathway in the trend mode is characterized by the cooling trend at 266-m depth from 1950 to 2007 along the coastal waveguide off Western Australia, which extends across southern Australia (Fig. 2c). This is in agreement with Feng et al. (2004), who show that multidecadal variations of coastal sea level at Fremantle in Western Australia are influenced by trade wind variations in the tropical Pacific Ocean via the equatorial and coastal waveguide. This pathway closely corresponds with that of the Leeuwin Current, where it has been previously shown that this warm surface ocean current flows strongly southward above 300-m depth and then eastward across southern Australia to as far as Tasmania (Feng et al. 2009). The connection between the western equatorial Pacific and the Leeuwin Current (Fig. 2c) is also in agreement with Feng et al. (2009, 2011), who suggest that the subsurface cooling transmitted from the tropical west Pacific to the southeast Indian Ocean also show a multidecadal weakening trend of the strength of Leeuwin Current between the mid-1970s and the mid-1990s. However, Feng et al. (2012) more recently show that the decreasing trend in the Leeuwin Current volume transport reversed after the mid-1990s so that the Leeuwin Current transport has been increasing overall ever since. This reversal in the volume transport of the Leeuwin Current is related to the change in phase of the IPO from positive to negative at the end of the 1990s and is captured in the spatial pattern and corresponding time series results from our EOF analysis of the detrended upper-ocean temperatures (see Fig. 3e). This finding highlights the need to estimate the trends in climate change from time series much longer than 50 years. It is important to recognize the difference between climate change and multidecadal variation, because we know that changes in the Leeuwin Current temperature can affect the recruitment of western rock lobster (Clarke and Li 2004), fish, and invertebrates, in particular in the Kalbarri region in Western Australia (Feng et al. 2012). Adaptation to climate change trends will be different than managing for interannual-to-decadal-scale variability.

EOF1-D (EOF1 of detrended and decadally filtered SODA upper-ocean temperatures) takes into account the temperature variability across four vertical levels simultaneously. We found that the sea surface spatial pattern of this multilevel EOF mode (Fig. 3a) is characteristic of the IPO, as compared with Fig. 2 of Power et al. (1999). Our analysis of the SODA data shows that the IPO signature extends into the subsurface layers of the upper ocean, and is indeed consistent with an ocean dynamic response of the upper ocean that may be considered independently from heat flux forcing across the sea surface (e.g., McGregor et al. 2004, 2007; Power and Colman 2006).The deepening of the thermocline associated with the surface signature of the IPO raises the important question: What role does the thermocline play in the surface layer heat budget?

McPhaden and Zhang (2002) provide an example of the important role of the pycnocline/thermocline on the surface heat budget. For instance, they suggest that an increase in sea surface temperature for the period 1970–90 in the equatorial Pacific cannot be explained by an increase in surface heating, because the heat fluxes into the ocean have decreased in the last 50 years in the equatorial Pacific. It is suggested that this might be associated with a decrease of the meridional overturning circulation between the tropics and subtropics, which was a manifestation of the phase shift associated with the Pacific decadal oscillation (PDO) in 1976/77 and possibly another one in the late 1990s in the North Pacific Ocean (McPhaden and Zhang 2002, 2004; Chavez et al. 2003; Peterson and Schwing 2003). The way that the overturning circulation works is that cold temperature anomalies of water masses created in the surface of the subtropics, by air–sea interactions, are subducted and transported equatorward through the upper pycnocline and are upwelled to the surface when they reach the Pacific equator. Further studies for the periods 1992–98 and 1998–2003 suggest that the cold pycnocline water within the meridional overturning circulation has increased in the latter period and thus has increased the upwelling to cool the surface in the eastern and central equatorial Pacific Ocean (McPhaden and Zhang 2004). Hence, the transport of the overturning circulation, through the pycnocline, oscillates on decadal time scales, producing sea surface decadal variability in the equatorial Pacific (McPhaden and Zhang 2002, 2004) and likely affecting the surface heat budget.

We found that the second detrended SSTA mode (EOF2-D5) of variability in the SODA data had a hemispherically asymmetric shape, and its PC2-D5 was significantly correlated with the Atlantic multidecadal oscillation (AMO) index. However, the second mode PC of the detrended decadal temperature anomalies across four depth levels (EOF2-D) was correlated with the North Pacific Gyre Oscillation (NPGO) index. PC2-D5 could be tracking the North Pacific aerosol loading as it increases rapidly after the 1960s; that is, it is possible that the SST mode is forced through an AMO teleconnection, similar to the “atmospheric bridge” discussed in previous studies (e.g., Alexander et al. 2002), while the corresponding subsurface mode is characterized by gyre-scale dynamics that appear to be dominated by the NPGO signature. There is consensus that anthropogenic and natural aerosols have likely played some role in forcing the North Atlantic SST, and it would be interesting to know whether AMO changes are anthropogenic aerosol or natural forced and whether this can affect the Pacific Ocean by atmospheric teleconnection (Chiang and Friedman 2012).

We have also performed an EOF analysis on the 10-yr low-pass filtered and detrended sea surface temperature anomalies across the global oceans (i.e., including the Atlantic sector), from the equator to 60° latitude (60°N–60°S) to determine whether the inclusion of the Atlantic Ocean in the EOF analysis might lead to an AMO-like pattern in the Atlantic Ocean. We found that the second EOF from this analysis (EOF2-D5-G) produces an AMO-like global spatial pattern (not shown) with the features in the Pacific Ocean being very similar to those found from the EOF analysis applied to the Indo-Pacific Ocean region (EOF2-D5; Fig. 11b). The pattern in the North Atlantic Ocean extends from 60°N to the equator and resembles the AMO pattern discussed in previous studies (e.g., Schlesinger and Ramankutty 1994; Enfield et al. 2001; Knight et al. 2005; Baines 2011). The PC time series from the EOF analysis correlates strongly with the AMO index [r = −0.58; significant at the >94% confidence level, according to Davis (1976)]. Interestingly, PC2-D5 from the EOF analysis applied to the Indo-Pacific Ocean is more strongly correlated with the AMO index than it is with PC2 from the global EOF analysis (Table 1). The AMO index is based on the Atlantic SSTA north of the equator (from 0° to 70°N) (Enfield et al. 2001). Both spatial patterns for the Pacific Ocean and North Atlantic Ocean during an AMO event show the same sign. However, this sign is opposite to that for the South Atlantic Ocean. Thus, including the entire Atlantic Ocean in the global EOF analysis appears to weaken the correlation between the PC time series and the AMO index.

Clearly, more research beyond the scope of this paper is needed to better understand the hemispheric, asymmetric shape of EOF2-D5 and the high correlation of PC2-D5 with the AMO index. Some questions that arise from these results are as follows: How much does the Atlantic Ocean influence the Pacific Ocean SST via atmospheric teleconnection? Is this influence anthropogenic aerosol forced or natural? A key part of the question here is whether the AMO is aerosol forced as well. The single forcing (historicalMisc AA) runs from phase 5 of the Coupled Model Intercomparison Project (CMIP5) may help to answer those questions. The historicalMisc runs include simulations forced by various combinations of anthropogenic and natural aerosols that are suitable for addressing the above questions in future work.

EOF3-D from the multilevel EOF analysis and EOF3-D5 from the SSTA-only analysis are both strongly correlated to the IEMI, suggesting that these decadal-scale upper-ocean temperature modes (here, EOF3) in the SODA data describe the low-frequency variability of El Niño Modoki through the upper-ocean temperature structure. Based on this analysis of SODA upper-ocean temperatures, we find that El Niño Modoki is largely confined to the upper 268 m, with weak signals just below the thermocline (Figs. 7a,b,c). We note that EOF3-D and EOF3-D5 imply that the decadally modulated El Niño Modoki signature has been more frequent and intense from the 1990s onward, appearing as two positive decadal modulations between 1990 and 2007 (Fig. 7e). Lee and McPhaden (2010) find that the warming trend in the Niño 4 (west-central Pacific) region on interannual time scales is mainly a result of the more frequent and intense central Pacific (or Modoki) El Niño events. The present results suggest that subsurface ocean helps to modulate El Niño Modoki on decadal time scales.

We have shown that the NPGO and the low-frequency component of El Niño Modoki have similar significant spectral peaks with maximum periods of 12.8 yr (Figs. 1b,c). Given that the domain of interest in this study is much broader than just the equatorial and North Pacific Ocean, this suggests that the decadal-scale connection reported between the NPGO and El Niño Modoki via the atmospheric bridge (Di Lorenzo et al. 2010; Chavez et al. 2011; Messié and Chavez 2011) is either expressed more broadly as an Indo-Pacific teleconnection or that there may be other oceanic interplays operating.

Our analysis has shown that the largest amplitude in decadal temperature variability in the SODA data is expressed near the base of the mixed layer/upper thermocline (at the 97-m depth level in SODA), with the largest amplitudes in the tropical Pacific Ocean and south Indian Ocean within the thermocline. The strong decadal variability across the tropical Pacific Ocean identified in the EOF modes of subsurface temperature is consistent with changes in the tropical Pacific winds, which also induce decadal-to-multidecadal variability of regional sea level, particularly in the western tropical Pacific, as discussed by Han et al. (2014b). Surface changes of dynamic variables in the south Indian Ocean are largely influenced by air–sea heat fluxes and decreased upwelling-related oceanic cooling (Alory and Meyers 2009); whereas the subsurface decadal variability is largely influenced by local wind stress (Nidheesh et al. 2013; Trenary and Han 2013; Han et al. 2014a), remote dynamics from the Pacific (Reason et al. 1996; Trenary and Han 2008; Schwarzkopf and Böning 2011; Trenary and Han 2013), and changes in salinity anomalies (Vargas-Hernández 2014).

In this paper, we have shown, based on reanalysis data from SODA 2.2.4 for the period 1950–2007, that the subsurface temperature structure is central to our understanding of the upper-ocean decadal climate variability in the Indo-Pacific Ocean region and is critical to the ocean dynamics and, hence, the ocean’s role in decadal climate variability. Particularly, we have shown that known large-scale modes of climate variability, including the IPO, NPGO, and El Niño Modoki, have a significant vertical extension through the upper ocean, including the thermocline, demonstrating that the vertical structure of the upper ocean is critically important to the decadal climate dynamics of the Indo-Pacific Ocean.

Acknowledgments

This work was conducted as part of the José Mauro Vargas Hernández Ph.D. in Quantitative Marine Sciences (QMS) Ph.D. Program, a joint initiative between the University of Tasmania and CSIRO. This research was partially supported by grants from the QMS Ph.D. Scholarship, the Postgraduate Studies Scholarship from Universidad Nacional of Costa Rica, and the CSIRO Wealth from Oceans Flagship Scholarship. The authors are thankful to Dr. Eric Oliver for advice on the spectral analysis used in this study. The SODA 2.2.4 data are provided by the SODA/TAMU research group, Texas A&M University (available at http://soda.tamu.edu/). The IPO index was obtained from the calculations of Chris Folland, Met Office Hadley Centre (available at http://www.iges.org/c20c/IPO_v2.doc). The NPGO index used here was generated by Di Lorenzo et al. (2008) (available at http://www.o3d.org/npgo/npgo.php). The AMO index used here is from Kaplan’s SST dataset (available from NOAA at http://www.esrl.noaa.gov/psd/data/timeseries/AMO/). [Kaplan SST V2 data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado (available at http://www.esrl.noaa.gov/psd/); see also Kaplan et al. (1998).] This paper makes a contribution to the goals of the Australian Research Council Centre of Excellence for Climate System Science and the Australian Climate Change Science Program.

REFERENCES

REFERENCES
Alexander
,
M. A.
,
I.
Bladé
,
M.
Newman
,
J. R.
Lanzante
,
N.-C.
Lau
, and
J. D.
Scott
,
2002
:
The atmospheric bridge: The influence of ENSO teleconnections on air–sea interaction over the global oceans
.
J. Climate
,
15
, 2205–2231, doi:.
Alory
,
G.
, and
G.
Meyers
,
2009
:
Warming of the upper equatorial Indian Ocean and changes in the heat budget (1960–99)
.
J. Climate
,
22
,
93
113
, doi:.
Alory
,
G.
,
S.
Wijffels
, and
G.
Meyers
,
2007
: Observed temperature trends in the Indian Ocean over 1960–1999 and associated mechanisms. Geophys. Res. Lett.,34, L02606, doi:.
Arblaster
,
J.
,
G.
Meehl
, and
A.
Moore
,
2002
:
Interdecadal modulation of Australian rainfall
.
Climate Dyn.
,
18
,
519
531
, doi:.
Ashok
,
K.
, and
T.
Yamagata
,
2009
:
Climate change: The El Niño with a difference
.
Nature
,
461
, 481–484, doi:.
Ashok
,
K.
,
S.
Behera
,
S.
Rao
,
H.
Weng
, and
T.
Yamagata
,
2007
: El Niño Modoki and its possible teleconnection. J. Geophys. Res.,112, C11007, doi:.
Baines
,
P. G.
,
2011
:
Patterns of decadal climate variability and their impact on global rainfall
.
Procedia Environ. Sci.
,
6
,
70
87
, doi:.
Boyer
,
T.
, and Coauthors
,
2009
: World Ocean Database 2009, S. Levitus, Ed., NOAA Atlas NESDIS 66, 216 pp. [Available online at ftp://ftp.nodc.noaa.gov/pub/WOD/DOC/wod09_intro.pdf.]
Butterworth
,
S.
,
1930
:
On the theory of filter amplifiers
.
Wireless Eng.
,
7
,
536
541
.
Carton
,
J.
, and
B.
Giese
,
2008
:
A reanalysis of ocean climate using Simple Ocean Data Assimilation (SODA)
.
Mon. Wea. Rev.
,
136
,
2999
3017
, doi:.
Carton
,
J.
,
G.
Chepurin
,
X.
Cao
, and
B.
Giese
,
2000
:
A Simple Ocean Data Assimilation analysis of the global upper ocean 1950–95. Part I: Methodology
.
J. Phys. Oceanogr.
,
30
,
294
309
, doi:.
Chavez
,
F. P.
,
J.
Ryan
,
S. E.
Lluch-Cota
, and
M.
Ñiquen
,
2003
:
From anchovies to sardines and back: Multidecadal change in the Pacific Ocean
.
Science
,
299
,
217
221
, doi:.
Chavez
,
F. P.
,
M.
Messié
, and
J. T.
Pennington
,
2011
:
Marine primary production in relation to climate variability and change
.
Annu. Rev. Mar. Sci.
,
3
,
227
260
, doi:.
Chiang
,
J. C.
, and
A. R.
Friedman
,
2012
:
Extratropical cooling, interhemispheric thermal gradients, and tropical climate change
.
Annu. Rev. Earth Planet. Sci.
,
40
,
383
412
, doi:.
Clarke
,
A.
, and
J.
Li
,
2004
:
El Niño/La Niña shelf edge flow and Australian western rock lobsters
.
Geophys. Res. Lett.
,
31
, L11301, doi:.
Compo
,
G.
, and Coauthors
,
2008
: The 20th Century Reanalysis Project. Extended Abstracts, Third WCRP Int. Conf. on Reanalysis, Tokyo, Japan, World Climate Research Programme, 5 pp. [Available online at http://wcrp.ipsl.jussieu.fr/Workshops/Reanalysis2008/Documents/V5-511_ea.pdf.]
Compo
,
G.
, and Coauthors
,
2011
:
The Twentieth Century Reanalysis Project
.
Quart. J. Roy. Meteor. Soc.
,
137
,
1
28
, doi:.
Davis
,
R.
,
1976
:
Predictability of sea surface temperature and sea level pressure anomalies over the North Pacific Ocean
.
J. Phys. Oceanogr.
,
6
,
249
266
, doi:.
Di Lorenzo
,
E.
, and Coauthors
,
2008
: North Pacific Gyre Oscillation links ocean climate and ecosystem change. Geophys. Res. Lett.,35, L08607, doi:.
Di Lorenzo
,
E.
,
K.
Cobb
,
J.
Furtado
,
N.
Schneider
,
B.
Anderson
,
A.
Bracco
,
M.
Alexander
, and
D.
Vimont
,
2010
:
Central Pacific El Niño and decadal climate change in the North Pacific Ocean
.
Nat. Geosci.
,
3
,
762
765
, doi:.
Di Lorenzo
,
E.
, and Coauthors
,
2013
: Synthesis of Pacific Ocean climate and ecosystem dynamics. Oceanography,26, 68–81, doi:.
Di Lorenzo
,
E.
, and Coauthors
,
2015
: North Pacific Gyre Oscillation (NPGO) index. Accessed 29 June 2015. [Available online at http://www.o3d.org/npgo/npgo.php.]
Enfield
,
D. B.
,
A. M.
Mestas-Nuñez
, and
P. J.
Trimble
,
2001
:
The Atlantic Multidecadal Oscillation and its relation to rainfall and river flows in the continental U.S
.
Geophys. Res. Lett.
,
28
,
2077
2080
, doi:.
Feng
,
M.
,
Y.
Li
, and
G.
Meyers
,
2004
: Multidecadal variations of Fremantle sea level: Footprint of climate variability in the tropical Pacific. Geophys. Res. Lett.,31, L16302, doi:.
Feng
,
M.
,
E.
Weller
, and
K.
Hill
,
2009
: The Leeuwin Current. A marine climate change impacts and adaptation report card for Australia 2009, E. S. Poloczanska, A. J. Hobday, and A. J. Richardson, Eds., NCCARF Tech. Rep., 11 pp. [Available online at http://www.oceanclimatechange.org.au/content/images/uploads/Leeuwin_Current.pdf.]
Feng
,
M.
,
C.
Böning
,
A.
Biastoch
,
E.
Behrens
,
E.
Weller
, and
Y.
Masumoto
,
2011
:
The reversal of the multi-decadal trends of the equatorial Pacific easterly winds, and the Indonesian Throughflow and Leeuwin Current transports
.
Geophys. Res. Lett.
,
38
, L11604, doi:.
Feng
,
M.
,
N.
Caputi
, and
A.
Pearce
,
2012
: The Leeuwin Current. A marine climate change impacts and adaptation report card for Australia 2012, E. S. Poloczanska, A. J. Hobday, and A. J. Richardson, Eds., NCCARF Tech. Rep., 61–80. [Available online at http://www.oceanclimatechange.org.au/content/images/uploads/2012_Leeuwin_current_marine_report_card_2.pdf.]
Folland
,
C.
,
2008
: Interdecadal Pacific Oscilation (IPO) index. Met Office Hadley Centre, accessed 3 October 2012. [Available online at http://www.iges.org/c20c/IPO_v2.doc.]
Folland
,
C.
,
J.
Renwick
,
M.
Salinger
, and
A.
Mullan
,
2002
:
Relative influences of the Interdecadal Pacific Oscillation and ENSO on the South Pacific Convergence Zone
.
Geophys. Res. Lett.
,
29
,
21-1
21-4
, doi:.
Han
,
W.
,
J.
Vialard
,
M. J.
McPhaden
,
T.
Lee
,
Y.
Masumoto
,
M.
Feng
, and
W. P.
de Ruijter
,
2014a
: Indian Ocean decadal variability: A review. Bull. Amer. Meteor. Soc.,95, 1679–1703, doi:.
Han
,
W.
, and Coauthors
,
2014b
:
Intensification of decadal and multi-decadal sea level variability in the western tropical Pacific during recent decades
.
Climate Dyn.
,
43
,
1357
1379
, doi:.
Holbrook
,
N. J.
,
I. D.
Goodwin
,
S.
McGregor
,
E.
Molina
, and
S. B.
Power
,
2011
:
ENSO to multi-decadal time scale changes in East Australian Current transports and Fort Denison sea level: Oceanic Rossby waves as the connecting mechanism
.
Deep-Sea Res. II
,
58
,
547
558
, doi:.
Hurrell
,
J.
, and Coauthors
,
2010
: Decadal climate prediction: Opportunities and challenges. Proc. OceanObs’09: Sustained Ocean Observations and Information for Society (Vol. 2), Venice, Italy, European Space Agency, WPP-306, doi:.
Kaplan
,
A.
,
M. A.
Cane
,
Y.
Kushnir
,
A. C.
Clement
,
M. B.
Blumenthal
, and
B.
Rajagopalan
,
1998
: Analyses of global sea surface temperature 1856–1991. J. Geophys. Res.,103, 18 567–18 589, doi:.
Knight
,
J. R.
,
R. J.
Allan
,
C. K.
Folland
,
M.
Vellinga
, and
M. E.
Mann
,
2005
:
A signature of persistent natural thermohaline circulation cycles in observed climate
.
Geophys. Res. Lett.
,
32
, L20708, doi:.
Lee
,
T.
, and
M.
McPhaden
,
2010
: Increasing intensity of El Niño in the central-equatorial Pacific. Geophys. Res. Lett.,37, L14603, doi:.
Li
,
G.
,
B.
Ren
,
C.
Yang
, and
J.
Zheng
,
2010
:
Indices of El Niño and El Niño Modoki: An improved El Niño Modoki index
.
Adv. Atmos. Sci.
,
27
,
1210
1220
, doi:.
Mantua
,
N.
,
S.
Hare
,
Y.
Zhang
,
J.
Wallace
, and
R.
Francis
,
1997
:
A Pacific interdecadal climate oscillation with impacts on salmon production
.
Bull. Amer. Meteor. Soc.
,
78
,
1069
1079
, doi:.
McCabe
,
G.
,
M.
Palecki
, and
J.
Betancourt
,
2004
:
Pacific and Atlantic Ocean influences on multidecadal drought frequency in the United States
.
Proc. Natl. Acad. Sci. USA
,
101
,
4136
4141
, doi:.
McGregor
,
S.
,
N. J.
Holbrook
, and
S. B.
Power
,
2004
: On the dynamics of interdecadal thermocline depth and sea surface temperature variability in the low to mid-latitude Pacific Ocean. Geophys. Res. Lett.,31, L24201, doi:.
McGregor
,
S.
,
N. J.
Holbrook
, and
S. B.
Power
,
2007
:
Interdecadal sea surface temperature variability in the equatorial Pacific Ocean. Part I: The role of off-equatorial wind stresses and oceanic Rossby waves
.
J. Climate
,
20
,
2643
2658
, doi:.
McGregor
,
S.
,
N. J.
Holbrook
, and
S. B.
Power
,
2009
:
The response of a stochastically forced ENSO model to observed off-equatorial wind stress forcing
.
J. Climate
,
22
,
2512
2525
, doi:.
McPhaden
,
M. J.
, and
D.
Zhang
,
2002
:
Slowdown of the meridional overturning circulation in the upper Pacific Ocean
.
Nature
,
415
,
603
608
, doi:.
McPhaden
,
M. J.
, and
D.
Zhang
,
2004
:
Pacific Ocean circulation rebounds
.
Geophys. Res. Lett.
,
31
, L18301, doi:.
Meehl
,
G. A.
,
A.
Hu
, and
B. D.
Santer
,
2009
:
The mid-1970s climate shift in the Pacific and the relative roles of forced versus inherent decadal variability
.
J. Climate
,
22
,
780
792
, doi:.
Messié
,
M.
, and
F.
Chavez
,
2011
:
Global modes of sea surface temperature variability in relation to regional climate indices
.
J. Climate
,
24
,
4314
4331
, doi:.
Nidheesh
,
A.
,
M.
Lengaigne
,
J.
Vialard
,
A.
Unnikrishnan
, and
H.
Dayan
,
2013
:
Decadal and long-term sea level variability in the tropical Indo-Pacific Ocean
.
Climate Dyn.
, 41,
381
402
, doi:.
NOAA
,
2015
: Atlantic Multidecadal Oscillation (AMO) index, Kaplan SST V2 data. NOAA/OAR/ESRL/PSD, Boulder, CO, accessed 29 June 2015. [Available online at http://www.esrl.noaa.gov/psd/data/timeseries/AMO/.]
North
,
G. R.
,
T. L.
Bell
,
R. F.
Cahalan
, and
F. J.
Moeng
,
1982
:
Sampling errors in the estimation of empirical orthogonal functions
.
Mon. Wea. Rev.
,
110
,
699
706
, doi:.
Peterson
,
W. T.
, and
F. B.
Schwing
,
2003
:
A new climate regime in northeast Pacific ecosystems
.
Geophys. Res. Lett.
,
30
, 1896, doi:.
Power
,
S.
, and
R.
Colman
,
2006
:
Multi-year predictability in a coupled general circulation model
.
Climate Dyn.
,
26
,
247
272
, doi:.
Power
,
S.
,
T.
Casey
,
C.
Folland
,
A.
Colman
, and
V.
Mehta
,
1999
:
Inter-decadal modulation of the impact of ENSO on Australia
.
Climate Dyn.
,
15
,
319
324
, doi:.
Power
,
S.
,
M.
Haylock
,
R.
Colman
, and
X.
Wang
,
2006
:
The predictability of interdecadal changes in ENSO activity and ENSO teleconnections
.
J. Climate
,
19
, 4755–4771, doi:.
Reason
,
C.
,
R.
Allan
, and
J.
Lindesay
,
1996
:
Evidence for the influence of remote forcing on interdecadal variability in the southern Indian Ocean
.
J. Geophys. Res.
,
101
,
11 867
11 882
, doi:.
Ridgway
,
K.
,
2007
: Long-term trend and decadal variability of the southward penetration of the East Australian Current. Geophys. Res. Lett.,34, L13613, doi:.
Schlesinger
,
M. E.
, and
N.
Ramankutty
,
1994
:
An oscillation in the global climate system of period 65–70 years
.
Nature
, 367,
723
726
, doi:.
Schwarzkopf
,
F. U.
, and
C. W.
Böning
,
2011
: Contribution of Pacific wind stress to multi-decadal variations in upper-ocean heat content and sea level in the tropical south Indian Ocean. Geophys. Res. Lett.,38, L12602, doi:.
Smith
,
R.
,
J.
Dukowicz
, and
R.
Malone
,
1992
:
Parallel ocean general circulation modeling
.
Physica D
,
60
,
38
61
, doi:.
Smith
,
R.
, and Coauthors
,
2010
: The Parallel Ocean Program (POP) reference manual: Ocean component of the Community Climate System Model (CCSM). Los Alamos National Laboratory Tech. Rep. LAUR-10-01853, 140 pp. [Available online at http://www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf.]
Smith
,
T.
, and
R.
Reynolds
,
2003
:
Extended reconstruction of global sea surface temperatures based on COADS data (1854–1997)
.
J. Climate
,
16
,
1495
1510
, doi:.
SODA/TAMU Research Group
,
2008
: Simple Ocean Data Assimilation (SODA), version 2.2.4. Texas A&M University, accessed 24 October 2010. [Available online at http://soda.tamu.edu/.]
Timmermann
,
A.
,
S.
McGregor
, and
F.-F.
Jin
,
2010
:
Wind effects on past and future regional sea level trends in the southern Indo-Pacific
.
J. Climate
,
23
,
4429
4437
, doi:.
Trenary
,
L. L.
, and
W.
Han
,
2008
:
Causes of decadal subsurface cooling in the tropical Indian Ocean during 1961–2000
.
Geophys. Res. Lett.
,
35
,
L17602
, doi:.
Trenary
,
L. L.
, and
W.
Han
,
2013
:
Local and remote forcing of decadal sea level and thermocline depth variability in the South Indian Ocean
.
J. Geophys. Res. Oceans
,
118
,
381
398
, doi:.
Vargas-Hernández
,
J. M.
,
2014
: Decadal climate variability in the Indo-Pacific upper ocean from state estimates. Ph.D. thesis, University of Tasmania, 198 pp.
Vargas-Hernández
,
J. M.
,
S.
Wijffels
,
G.
Meyers
, and
N. J.
Holbrook
,
2014
:
Evaluating SODA for Indo-Pacific Ocean decadal climate variability studies
.
J. Geophys. Res. Oceans
,
119
,
7854
7868
, doi:.
Venegas
,
S. A.
,
2001
:
Statistical methods for signal detection in climate
. Danish Center for Earth System Science Tech. Rep. 2, 96 pp.
Wainwright
,
L.
,
G.
Meyers
,
S.
Wijffels
, and
L.
Pigot
,
2008
: Change in the Indonesian Throughflow with the climatic shift of 1976/77. Geophys. Res. Lett.,35, L03604, doi:.
Wijffels
,
S.
, and
G.
Meyers
,
2004
:
An intersection of oceanic waveguides: Variability in the Indonesian throughflow region
.
J. Phys. Oceanogr.
,
34
,
1232
1253
, doi:.
Woodruff
,
S. D.
, and Coauthors
,
2011
:
ICOADS Release 2.5: Extensions and enhancements to the surface marine meteorological archive
.
Int. J. Climatol.
,
31
,
951
967
, doi:.
Zhang
,
R.
, and
T.
Delworth
,
2007
: Impact of the Atlantic Multidecadal Oscillation on North Pacific climate variability. Geophys. Res. Lett.,34, L23708, doi:.