Decadal variabilities in Indian Ocean subsurface ocean heat content (OHC; 50–300 m) since the 1950s are examined using ocean reanalyses. This study elaborates on how Pacific variability modulates the Indian Ocean on decadal time scales through both oceanic and atmospheric pathways. High correlations between OHC and thermocline depth variations across the entire Indian Ocean Basin suggest that OHC variability is primarily driven by thermocline fluctuations. The spatial pattern of the leading mode of decadal Indian Ocean OHC variability closely matches the regression pattern of OHC on the interdecadal Pacific oscillation (IPO), emphasizing the role of the Pacific Ocean in determining Indian Ocean OHC decadal variability. Further analyses identify different mechanisms by which the Pacific influences the eastern and western Indian Ocean. IPO-related anomalies from the Pacific propagate mainly through oceanic pathways in the Maritime Continent to impact the eastern Indian Ocean. By contrast, in the western Indian Ocean, the IPO induces wind-driven Ekman pumping in the central Indian Ocean via the atmospheric bridge, which in turn modifies conditions in the southwestern Indian Ocean via westward-propagating Rossby waves. To confirm this, a linear Rossby wave model is forced with wind stresses and eastern boundary conditions based on reanalyses. This linear model skillfully reproduces observed sea surface height anomalies and highlights both the oceanic connection in the eastern Indian Ocean and the role of wind-driven Ekman pumping in the west. These findings are also reproduced by OGCM hindcast experiments forced by interannual atmospheric boundary conditions applied only over the Pacific and Indian Oceans, respectively.
Recent work has demonstrated the important role of the Indian Ocean in modulating global climate variability (SanchezGomez et al. 2008; Schott et al. 2009; Luo et al. 2012) and regional rainfall (Ashok et al. 2001; Ummenhofer et al. 2009). In particular, the role of upper-ocean heat content (OHC) in the Indian Ocean has been highlighted in recent discussions of the so-called global warming hiatus (Lee et al. 2015; Nieves et al. 2015). Several studies have linked the hiatus of the global-mean surface warming during the early 2000s to the heat redistribution in the upper 700 m of the Pacific and Indian Oceans (Liu et al. 2016). Although heat uptake in the Pacific Ocean increased during the hiatus, the OHC in the Pacific upper layer has decreased (Meehl et al. 2011; Balmaseda et al. 2013a). Strengthened easterly trade winds associated with a series of long-lasting La Niña events (England et al. 2014) resulted in an anomalously strong Pacific–Indian Ocean pressure gradient, contributing to an increase in Indonesian Throughflow (ITF) heat transport (Lee et al. 2015). Through this increased ITF heat transport, cooling in the top 100 m of the Pacific Ocean was mainly compensated by warming in the subsurface Indian Ocean (100–300 m) post 2000 (Nieves et al. 2015). Observations suggest that the rapid increase in Indian Ocean OHC accounted for more than 70% of the global ocean heat gain in the upper 700 m during 2003–12 (Lee et al. 2015).
In addition to OHC redistribution in the Indo-Pacific region, subsurface OHC is important to climate variability, both as a background state for the sea surface temperature (SST) and as an important factor in the development of Indian Ocean dipole events (Shinoda et al. 2004; Annamalai et al. 2005). Noting a high correlation between thermocline depth and SST in the southwestern tropical Indian Ocean, Xie et al. (2002) proposed that subsurface ocean dynamics can impact the SST variability in the southwestern tropical Indian Ocean and alter the state of the overlying atmosphere. Therefore, subsurface variability in the Indian Ocean is increasingly being recognized as an important factor in the climate system.
While the tropical Indian Ocean SST has exhibited faster warming since the 1950s than the tropical Atlantic or Pacific (Han et al. 2014), the subsurface tropical Indian Ocean has displayed a prominent cooling trend from the 1960s through the 1990s (Han et al. 2006; Alory et al. 2007). Surface warming has been primarily trapped above the 20°C isotherm in the tropical Indian Ocean from 15°S to 5°N, whereas a strong subsurface cooling trend has been observed within the tropical thermocline between 100- and 300-m depth north of 20°S. Previous studies have suggested various candidates as the drivers for this observed multidecadal cooling trend in the subsurface Indian Ocean. Han et al. (2006) proposed that the upper-thermocline cooling trend was primarily caused by local wind forcing that resulted in an enhancement of upward Ekman pumping. The increase in upward Ekman pumping velocities shoals the thermocline, inducing the cooling trend in the subsurface. By contrast, others suggested that the cooling in the southern Indian Ocean resulted from changes in the strength of the Pacific trade winds (Alory et al. 2007; Schwarzkopf and Böning 2011). Thermocline depth anomalies in the equatorial Pacific induced by changes in the trade winds are transmitted to the southern Indian Ocean via the Indonesian region in the form of baroclinic waves or ITF heat transport. Schwarzkopf and Böning (2011) reproduced the subsurface cooling trends in the Indian Ocean by applying observed wind forcing to the Pacific Ocean alone in an ocean general circulation model (OGCM), suggesting that remote forcing in the Pacific is an important contributor to low-frequency variations in the subsurface Indian Ocean. Zhou et al. (2017) further investigated contributions of local and remote forcing from the Pacific to the Indian Ocean by opening and closing the Indonesian passages in an OGCM. They suggested that the Pacific exerts large influence on subsurface variations in the Indian Ocean via oceanic baroclinic Rossby waves.
It is well known that Pacific Ocean variability modulates Indian Ocean conditions on interannual time scales (Du et al. 2009; Xie et al. 2009). El Niño–Southern Oscillation (ENSO), the most prominent interannual mode of climate variability, influences the Indian Ocean via both the Walker circulation (Klein et al. 1999; Alexander et al. 2002; Rao et al. 2002) and westward-propagating oceanic Rossby waves (Cai et al. 2005; Shi et al. 2007; Cai et al. 2008). The ascending branch of the Walker circulation shifts eastward during El Niño events, resulting in anomalous easterlies over the eastern Indian Ocean. These easterly anomalies influence latent heat fluxes and thermocline depths in the eastern Indian Ocean. In the paradigm of the oceanic connection, westward-propagating Rossby waves generated by zonal wind anomalies over the Pacific Ocean become coastally trapped waves at the intersection of New Guinea and the equator. These waves then propagate poleward along the western coast of Australia and radiate into the southern Indian Ocean (Li and Clarke 2004; Wijffels and Meyers 2004). The ITF, which connects the Pacific Ocean to the Indian Ocean, delivers large amounts of heat from the western Pacific to the eastern Indian Ocean. Variations in the ITF can thus change the thermal properties of the upper layer of the Indian Ocean (Godfrey 1996; Zhou et al. 2008). However, the relative influence of atmospheric and oceanic pathways from the Pacific Ocean to the Indian Ocean remains uncertain, particularly on decadal time scales.
Furthermore, the eastern and western Indian Ocean are known to be forced by different mechanisms on interannual time scales (Klein et al. 1999; Murtugudde and Busalacchi 1999). The net heat flux, which responds primarily to changes in cloud cover and surface latent heat flux, is the primary factor in SST variations in the eastern Indian Ocean during ENSO events. However, correlations between heat flux anomalies and SST tendencies are weak in the western Indian Ocean, suggesting that oceanic mechanisms are the primary contributors to SST variability in that part of the basin (Klein et al. 1999). Xie et al. (2002) identified a key ocean dynamic process in the western Indian Ocean where wind stress curl associated with anomalous easterlies induces downwelling Rossby waves on interannual time scales. However, it is still unclear whether the drivers of the variability on decadal time scales are distinct between the eastern and western parts of the tropical Indian Ocean. Previous studies of decadal variability have typically treated the tropical Indian Ocean as a whole (Han et al. 2006; Trenary and Han 2008; Dong et al. 2016) and focus on the zonal-mean trend in the Indian Ocean, despite anomalies of opposite sign in the eastern and western parts being reported in some studies (Vargas-Hernandez et al. 2015). Indeed, performing a set of sensitivity experiments in a regional ocean model, Trenary and Han (2013) suggested that thermocline variations are largely determined by wind stress curl imposed on the western Indian Ocean, while the ITF dominates variability in the eastern Indian Ocean after the1990s. Therefore, unlike previous studies that primarily focused on the zonal-mean trends, here we suggest the existence of distinct mechanisms for OHC variations in the eastern and western Indian Ocean on decadal time scales. We synthesize observations, theoretical considerations, and output from a high-resolution OGCM to elaborate on how the Pacific Ocean impacts OHC in the southern tropical Indian Ocean through oceanic and atmospheric pathways.
The interdecadal Pacific oscillation (IPO), defined as the leading empirical orthogonal function (EOF) mode of low-pass-filtered SST anomalies over the Pacific Ocean, is the dominant mode of variability in the Pacific on decadal time scales (Meehl and Arblaster 2002). The IPO is highly correlated with the Pacific decadal oscillation (PDO; Mantua et al. 1997), with some studies suggesting that the PDO should be regarded as an expression of the IPO in the North Pacific (Folland et al. 2002). However, the IPO and PDO are not identical, as there are notable differences in their spatial distributions. The northern and southern centers of activity have roughly comparable amplitudes in the IPO pattern, while the northern signal is preeminent in the PDO pattern (Newman et al. 2016). Lead–lag correlations of sea surface height anomalies (SSHA) in the Indo–Pacific Ocean region from Simple Ocean Data Assimilation products against the IPO index suggest a meridionally asymmetric dynamical connection between the western Pacific Ocean and the southern tropical Indian Ocean, whereby signals from the Pacific propagate into the Indian Ocean through the Indonesian Seas (Vargas-Hernandez et al. 2014).
Noting that most previous studies focus on long-term trends in the Indian Ocean from the early 1960s to late 1990s with the available dataset ending in the early 2000s (Han et al. 2006; Alory et al. 2007; Han et al. 2014), it is therefore timely to investigate the decadal variability based on more up-to-date datasets. We will focus particularly on assessing the processes through which the Pacific influences the Indian Ocean on decadal time scales. In this study, we suggest that decadal variations in the eastern Indian Ocean are mainly caused by oceanic signals entering from the Pacific in the form of baroclinic Rossby waves, while those in the western Indian Ocean are dominated by the local wind stress curl response to conditions in the Pacific via the atmospheric bridge.
The remainder of this paper is organized as follows. The data and model experiments used in this study are described in section 2. Decadal variations in the Indian Ocean subsurface OHC and their connection with IPO variability are examined using reanalysis products in section 3. A linear baroclinic Rossby wave model is then used to illustrate the relative roles of local wind stress curl and remote Pacific Ocean forcing to the Indian Ocean (section 3c). OGCM sensitivity experiments provide additional support for the hypotheses derived from the ocean reanalyses and linear baroclinic Rossby wave model (section 3d). Summary and discussion are provided in section 4.
2. Data and method
a. Observational data and reanalysis products
Several observational and reanalysis products are used in this study. The GECCO2 reanalysis (Köhl 2015), produced by the German contingent of the Estimating the Circulation and Climate of the Ocean project (ECCO; www.ecco-group.org), is used as the primary dataset to investigate the decadal variability across the Indian and Pacific Oceans during the period of 1948–2014. GECCO2 has been produced using the Massachusetts Institute of Technology General Circulation Model (MITgcm) with 50 vertical levels and a horizontal grid spacing of 1°. An adjoint method has been used to adjust the model outputs for consistency with available observational data. The background atmospheric state is taken from the 6-hourly National Centers for Environmental Prediction (NCEP)–National Center for Atmospheric Research (NCAR) reanalysis (R-1).
In addition, the Hadley Centre EN4.0.2 dataset, the Simple Ocean Data Assimilation version 2.2.4 reanalysis (SODA2.2.4), and the European Centre for Medium-Range Weather Forecasts (ECMWF) Ocean Reanalysis System, version 4 (ORAS4) are also used as ancillary datasets to assess the robustness of findings based on GECCO2. EN4, an objective gridded hydrographic dataset produced by the Met Office Hadley Centre (Good et al. 2013), covers the period from 1900 to the present with a 1° horizontal grid spacing, and its main data source is the 2009 version of the World Ocean Database (WOD09). Extensive quality control procedures are used to generate monthly objective analyses of ocean temperature and salinity, along with uncertainty estimates. The gridded data are relaxed to the 1971–2000 climatology when observational data are unavailable. SODA2.2.4 is provided on a 0.5° regular grid for the period of 1871–2008 (Carton et al. 2000; Carton and Giese 2008). This reanalysis was constructed using version 2.0.1 of the Parallel Ocean Program model (POP2.0.1) forced by surface wind stresses from the NOAA–Cooperative Institute for Research in Environmental Sciences (CIRES) Twentieth Century Reanalysis (20CRv2). Surface heat fluxes are calculated using bulk formulas, and data assimilation in SODA is conducted using the sequential estimation method. ORAS4 (Balmaseda et al. 2013b) uses version 3.0 of the Nucleus for European Modelling of the Ocean (NEMO3) ocean model with a standard horizontal grid spacing of 1°. Wind stresses and other atmospheric forcings are taken from the 40-yr ECMWF Re-Analysis (ERA-40) from September 1957 to December 1989, the ECMWF interim reanalysis (ERA-Interim) from 1989 to 2009, and the ECMWF operational analyses from 2010 onward. The NEMO variational data assimilation (NEMOVAR), which uses a three-dimensional variational data assimilation system (3DVAR) in a first guess at appropriate time configuration, has been introduced in ORAS4. ORAS4 uses a model bias-correction scheme, which is a latitudinal-dependent method. Key details of the ocean reanalyses used in this study are listed in Table 1. Despite discrepancies among these different datasets, the results obtained in this analysis are robust. We therefore only present those results based on GECCO2 in the following section.
Monthly multisatellite merged SSHA for 1993–2014 on a 0.25° horizontal grid are obtained from the Collecte Localisation Satellites (CLS) Archiving, Validation, and Interpretation of Satellite Oceanographic Data (AVISO). Monthly SST data are taken from the NOAA Extended Reconstructed SST (ERSST) version 3 dataset, which is available starting from 1854 (Smith et al. 2008). Surface wind stress are from ERA-40 for 1958–2001 (Uppala et al. 2005). Twentieth Century Reanalysis (20CR) wind stress data for 1871–2010 (Compo et al. 2011) and R-1 data from 1948 to the present (Kalnay et al. 1996) are also used to assess the robustness of the results. Discrepancies among reanalysis datasets may arise from differences in physical parameterizations, model boundary conditions, bias-correction procedures, assimilated observational data, and other factors (Fujiwara et al. 2017). Comparisons among these datasets help to confirm the robustness of our results regardless of different processing techniques used in each dataset, and sensitivities to the choice of dataset are discussed in section 4.
Monthly anomalies of each variable are computed by subtracting the corresponding monthly climatology. Bilinear interpolation is used to map all variables onto a common 1° × 1° regular latitude–longitude grid. To isolate variability at decadal and longer time scales, the long-term linear trend is first removed, and then an 8-yr low-pass Butterworth filter (Butterworth 1930) is applied to all variables. The first and last four years of data are excluded in our analysis to avoid artifacts that arise from the low-pass filtering. The results are similar when 8–30-yr bandpass filters are applied, indicating that decadal variability dominates variations at time scales longer than 8 yr in these datasets. EOF is calculated based on the covariance matrix of the corresponding variable with the weight proportional to the square root of the area. Principal component (PC) time series are normalized, so that the spatial EOF patterns show typical magnitudes associated with one standard deviation of the corresponding PCs.
The IPO is defined as the leading EOF of low-pass-filtered monthly SST anomalies (based on ERSST) over the Pacific Ocean (60°S–60°N) and its corresponding first PC (PC1) time series is herein referred to as the IPO index. The positive phase of the IPO is characterized by anomalously warm SST in the eastern tropical Pacific Ocean. We use the IPO index rather than the PDO index for this analysis because the IPO better represents variability in the entire Pacific Ocean (Power et al. 1999; Meehl and Arblaster 2002; Meehl and Hu 2006), whereas the PDO primarily characterizes SST variability in the North Pacific Ocean (Alexander 2010). Because of the high degree of serial correlation in the low-pass-filtered time series, statistical significance is determined by the nonparametric random phase method (Ebisuzaki 1997). Unless specified otherwise, all correlation coefficients are statistically significant at the 95% confidence level. The 20°C isothermal depth (D20) is used as a proxy for thermocline depth (Wyrtki and Kendall 1967; Xie et al. 2002; Annamalai et al. 2003).
b. Ocean model simulation
The OGCM used in this study is the numerical ocean–sea ice model NEMO, version 2.3, (Madec et al. 2015) in its global tripolar configuration at 0.5° horizontal resolution (ORCA05 configuration) with 46 vertical levels ranging in thickness from 6 m near the surface to 250 m in the deepest layers, while bottom topography is represented by partial steps (Barnier et al. 2006). The experiments are driven by prescribed atmospheric boundary conditions as given by the Co-ordinated Ocean–Ice Reference Experiments (CORE; Griffies et al. 2009, based on Large and Yeager 2004). The reference simulation (REF) is a hindcast forced by interannually varying fields for the period 1958–2004, preceded by a 20-yr-long spinup. Additionally, two sensitivity experiments have been conducted. In the PAC experiment, the interannual wind stress and heat fluxes are applied only over the Pacific Ocean, while the rest of the ocean is forced with repeated climatological fields; and in the IND experiment, the interannual forcing is only applied in the Indian Ocean. To remove spurious model drift, the trend in an experiment driven by climatological forcings (CLIM) was subtracted from the interannually forced cases prior to further analysis.
The experiments and their performances have been described and evaluated in previous studies (Schwarzkopf and Böning 2011; Ummenhofer et al. 2013), showing that the Indo-Pacific upper-ocean properties are simulated remarkably well, including the main features of the observed mean seasonal cycle and its associated variance. Compared with the ORAS4 reanalysis, REF also reproduces the linear trend during the 1960s–1990s, with subsurface cooling at 50–300-m depth in the 10°–20°S band (Ummenhofer et al. 2017). Good agreement between observed and modeled SST anomalies in the eastern Indian Ocean indicates the model’s credibility. However, the amplitudes of the Indian Ocean dipole (IOD) are overestimated. which can be attributed to the shallower thermocline off the coast of Sumatra–Java with respect to observations (Ummenhofer et al. 2013).
a. Decadal variability of the subsurface Indian Ocean
To investigate decadal variations in the subsurface Indian Ocean, EOF analysis is performed on detrended 8-yr low-pass-filtered subsurface (50–300 m) OHC. The first two EOF modes and their corresponding PC time series based on the GECCO2 reanalysis are shown in Fig. 1. The first mode, accounting for around 30% of the variance, is characterized by a meridional dipole between the northern and southern Indian Ocean and manifests on the multidecadal time scale (Fig. 1c). This multidecadal seesaw of OHC between the northern and southern parts of the tropical Indian Ocean has previously been linked with variations in the meridional overturning circulation of the upper Indian Ocean (Lee 2004) and ITF heat transport (Dong and McPhaden 2016). The second EOF mode (Fig. 1b), which accounts for around 20% of the OHC total variance, is expressed primarily on decadal time scales (Fig. 1d). The maximum anomalies are found between 10° and 20°S in the southern Indian Ocean, with zonal dipole structure characterized by positive anomalies west of 100°E and negative anomalies to the east. The first two EOF patterns are robust and also exhibited in other datasets, although the order of the two leading modes is sensitive to the choice of dataset or time period (figure not shown). As our intention is to address decadal variability in this study, we therefore focus on the latitude band with the maximum signal in EOF2 (10°–20°S) in the following analysis.
The climatological thermocline depth (here represented by D20) in the Indian Ocean is typically between 100 and 200 m below the surface (Fig. 2a), with a thermocline dome in the band of 5°–10°S in the western Indian Ocean. Temporal variations in thermocline depth are closely associated with those in subsurface OHC as illustrated by the strong correlation between D20 and subsurface OHC (Fig. 2b), where the latter is integrated between 50- and 300-m depth. The correlation between D20 and subsurface OHC exceeds 0.9 for most of the Indian Ocean. Moreover, the correlation map between D20 and PC2 of OHC anomalies (Fig. 2c) displays a similar east–west dipole structure as the EOF2 (Fig. 1b), indicating that the decadal variation of OHC is dominated by that in D20 where vertical gradients of temperature are largest and hence, the variation of OHC reflects the dynamical fluctuation in the thermocline. In our subsequent analysis, we will therefore focus especially on decadal variations in subsurface OHC associated with fluctuations in the thermocline.
b. Connections between the Indian and Pacific Oceans
In this section, we investigate how the Pacific Ocean influences the Indian Ocean on decadal time scales. The IPO index is used to represent the Pacific Ocean decadal variability, and the regression map of the Indian Ocean subsurface OHC anomalies onto the IPO is shown in Fig. 3. A zonal dipole structure is evident in the band of 10°–20°S in the southern Indian Ocean (black box in Fig. 3). The regression map bears a remarkable resemblance to the OHC EOF2 spatial pattern (Fig. 1b). The spatial correlation between Fig. 1b and Fig. 3 is 0.85 (P < 0.05), suggesting that Indian Ocean subsurface decadal variability is associated with the state of the Pacific Ocean.
Figure 4 shows the longitude–time Hovmöller diagram of 8-yr low-pass-filtered anomalies in subsurface OHC based on the GECCO2 reanalysis for 1948–2014. OHC anomalies are averaged along 10°–20°S where the largest variance occurs in the southern Indian Ocean (Fig. 4a), and along 5°–10°N in the Pacific Ocean (Fig. 4b). The 8-yr low-pass-filtered IPO index based on ERSST is also shown for reference (Fig. 4c). Previous studies (Wijffels and Meyers 2004; Cai et al. 2005) have suggested that signals in the tropical North Pacific Ocean (here represent by 5°–10°N; red box in Fig. 3) could propagate into the southern Indian Ocean through the ITF region. Although the South Pacific Ocean also contributes to Indian Ocean variations on decadal time scales (Vargas-Hernandez et al. 2014), IPO-related anomalies in the Pacific are almost symmetric. We therefore use the signal in the North Pacific to represent variations in both the North and South Pacific. As shown in the regression map (Fig. 3), OHC anomalies in the western Pacific link with the eastern Indian Ocean through the Indonesian Seas, suggesting an oceanic connection between the western Pacific and eastern Indian Ocean on decadal time scales. Furthermore, anomalies in the eastern Indian Ocean vary consistently with those in the western Pacific Ocean at a lag of several months (Fig. 4), implying the oceanic pathway for the Pacific Ocean to influence the eastern Indian Ocean decadal variability. The OHC anomalies entering at the eastern boundary of the southern Indian Ocean originate from the western Pacific Ocean and then propagate westward via oceanic Rossby waves or advection to impact the western Indian Ocean on decadal time scales (Li and Clarke 2004; Wijffels and Meyers 2004; Zhou et al. 2017). However, discontinuities often appear in these westward-propagating OHC anomalies around 90°–100°E (Fig. 4a). Indeed, many of the westward-propagating OHC anomalies in the western part of the Indian Ocean appear to originate from the central Indian Ocean rather than from the eastern boundary, particularly during the 1970s and 1980s. In addition, the OHC anomalies during that period were organized into dipole structures, with opposing anomalies in the eastern and western parts of the south Indian Ocean. A similar discontinuity emerged near the same location after 2005. The frequent discontinuities in signal propagation in the central Indian Ocean and the prominent east–west dipole structure in OHC anomalies lead us to hypothesize that the mechanisms by which the Pacific Ocean influences the southern Indian Ocean on decadal time scales are distinct between the eastern and western parts of the Indian Ocean.
To further explore the influence of Pacific variability on the Indian Ocean, we calculate lead–lag correlations between subsurface OHC anomalies and the IPO index (Fig. 5). Correlations are predominantly negative in the Indian Ocean east of 100°E, with the maximum amplitude of the correlation coefficient lagging the IPO index by about one year. This negative correlation propagates westward from the west coast of Australia, reaching 100°E in the central Indian Ocean after approximately 1.5 years. The relatively slow propagation of this signal suggests that the influence of the Pacific Ocean on the eastern Indian Ocean operates primarily through oceanic pathways. The estimated propagation speed of these decadal-scale anomalies in the eastern Indian Ocean is approximately 5.5 cm s−1, consistent with previous estimates of Rossby wave propagation in this latitude band (White 2000).
By contrast, significant positive correlations west of 100°E peak at zero lag. The nearly instantaneous response of the central–western Indian Ocean to IPO variability suggests that the Pacific Ocean influences conditions in the western Indian Ocean mainly via the atmospheric bridge. Regressing Ekman pumping onto the simultaneous IPO index reveals a striking negative anomaly in Ekman pumping in the central Indian Ocean (Fig. 6). Positive wind stress curl over the central Indian Ocean during the positive phase of the IPO induces Ekman downwelling, which in turn triggers westward-propagating downwelling Rossby waves. The westward propagation of these downwelling Rossby waves deepens the thermocline and increases OHC in the western Indian Ocean. Note that the maximum OHC anomalies (around 80°E) are located slightly west of the region with the largest regression of Ekman pumping onto IPO (90°–100°E), where the discontinuities appear for those westward-propagating OHC anomalies from the eastern boundary origin. This mismatch can be reconciled by the fact that Ekman pumping impacts the OHC not only locally but also remotely via the cumulative signals of westward-propagating Rossby waves. To further examine the mechanisms described above, low-pass-filtered time series of wind stress curl in the central Indian Ocean and subsurface OHC in the western Indian Ocean are shown in Fig. 7a. Variations in OHC averaged over the western Indian Ocean follow variations in wind stress curl averaged over the central Indian Ocean (Fig. 7a). The lead–lag correlation between these two time series (Fig. 7b) peaks at a value of 0.9, with central Indian Ocean wind stress curl leading western Indian Ocean OHC by about 12 months. This result suggests that decadal-scale variations in subsurface OHC in the western Indian Ocean primarily emerge as a response to westward-propagating Rossby waves induced by wind-driven Ekman pumping in the central Indian Ocean. Variations in the latter are linked to the IPO via the atmospheric bridge. This relationship contrasts with that in the eastern Indian Ocean, where the IPO influence is largely communicated by the propagation of Rossby waves carrying oceanic signals from the Pacific westward through the eastern boundary.
c. Linear Rossby wave model
The statistical analysis described above suggests that the Pacific Ocean affects decadal variations in the eastern and western parts of the Indian Ocean through different mechanisms. The primary role in the western Indian Ocean is played by the surface wind forcing in the central Indian Ocean via the atmospheric bridge, whereas the primary role in the eastern Indian Ocean is played by ocean transport across the eastern boundary. In this section, this interpretation is further evaluated using a simple reduced gravity model (Woodberry et al. 1989; Périgaud and Delecluse 1992; Masumoto and Meyers 1998; Birol and Morrow 2001; Potemra 2001). Comparison with observations indicates that this simple linear model captures variations in SSHA well (Périgaud and Delecluse 1992). As SSHA are tightly correlated with vertical displacements in the thermocline, variations in SSHA reflect corresponding variations in D20 and subsurface OHC. The advent of satellite retrievals of SSHA has enabled global coverage and dramatically expanded data availability, especially in the Indian Ocean where observational records are sparse. Hindcast SSHA derived from the simple reduced gravity model are compared with both reanalysis (GECCO2) and satellite (AVISO) estimates.
Here, the simple linear 1.5-layer reduced gravity model is employed to separate the relative role of wind stress curl and the eastern boundary condition. Under the long-wave approximation, the large-scale linear vorticity equation of baroclinic ocean response to wind stress curl can be written as
where is the SSHA, cR is the observed speed of the first baroclinic-mode Rossby wave, g′ is the reduced gravity, ρ0 is the reference seawater density, f is the Coriolis parameter, k is the unit vector in the vertical direction, and τ is the wind stress. The solution can be obtained by integrating Eq. (1) from the eastern boundary (x = xe) along the baroclinic Rossby wave characteristic:
To get the hindcast SSHA field in the southern Indian Ocean, monthly wind stress from ERA-40 (1958–2001) is used and the Rossby wave speed at each latitude is estimated from the observed SSHA. The eastern boundary is set at 120.5°E, and h(xe, y, t) is obtained from the GECCO2 monthly SSHA at 120.5°E.
Hovmöller diagrams in Fig. 8 show part of the simulated SSHA averaged over the zonal band of 10°–20°S during 1993–2001 to visually validate against the observed satellite SSH. For comparison, SSHA from GECCO2 (Fig. 8d) and AVISO (Fig. 8e) during the same period are also shown. The simulated total SSHA bears an evident resemblance to the observed SSHA. The negative SSHA, for example, originated from the eastern boundary in 1997 (Figs. 8d,e) and propagated across the entire basin before reaching the western coast of the Indian Ocean in late 1999. Both reanalysis and satellite estimates likewise show positive anomalies emerging in the eastern Indian Ocean as the negative anomaly propagated westward across the basin. The linear model (Fig. 8c) captures both features well. Relative to the eastern Indian Ocean, wind-driven Ekman pumping causes the larger magnitude of SSHA variation in the western Indian Ocean, especially west of 100°E (Fig. 8a).
To evaluate the hindcast skill of the linear Rossby wave model, the longitudinal distribution of correlations between SSHA from the Rossby wave model and GECCO2 is calculated during the period of 1960–2001 (Fig. 9). All estimates have been low-pass filtered to remove interannual variability. Correlations between the SSHA from the different components of the Rossby wave model, that is, those generated by wind-driven Ekman pumping (Fig. 9, blue line), eastern boundary condition (Fig. 9, green line), and total SSHA (Fig. 9, red line; the sum of SSHA caused by Ekman pumping and eastern boundary condition) against GECCO2 are shown. The component of the model-simulated SSHA generated by wind-driven Ekman pumping (Fig. 9, blue line) is significantly correlated with GECCO2 estimates in the western Indian Ocean but not in the eastern Indian Ocean. Conversely, the correlation between the component of model-simulated SSHA due to the eastern boundary condition and GECCO2 (Fig. 9, green line) is strong and significant in the eastern Indian Ocean but gradually diminishes as distance to the eastern boundary increases. Zonal variations in the correlation between the simulated total SSHA and GECCO2 (Fig. 9, red line) are likewise dominated by the eastern boundary condition east of 100°E, but by the wind-driven Ekman pumping west of 90°E. These results support our contention that oceanic signals from the western Pacific Ocean dominate decadal-scale thermocline and OHC variations in the eastern Indian Ocean, while wind-driven Ekman pumping associated with the IPO (Fig. 6) dominates these variations in the western Indian Ocean.
d. OGCM experiments
A REF simulation and two OGCM sensitivity experiments (PAC and IND) are used to evaluate the role of the Pacific Ocean in Indian Ocean decadal variability (see section 2b). Figure 10 shows the zonal propagation of 8-yr low-pass-filtered subsurface OHC averaged over the 10°–20°S band in the southern Indian Ocean for 1960–2004. The output from the REF simulation (Fig. 10b) broadly resembles the GECCO2 reanalysis (Fig. 10a). During the 1960s, anomalies that entered the southern Indian Ocean at the eastern boundary then propagated across the Indian Ocean to the western Indian Ocean. Both the reanalysis and the model output indicate that signal propagation from the eastern boundary was often blocked at other times, with discontinuities consistently emerging near 100°E. During these periods, OHC anomalies in the southern Indian Ocean were organized into a zonal dipole with opposite signs in the eastern and western parts of the basin. Both the dipole structure and the discontinuity near 100°E were more prominent after the mid-1970s.
Two parallel sensitivity experiments were conducted by imposing interannually varying surface heat fluxes and wind stresses only in the Pacific Ocean (PAC experiment) or Indian Ocean (IND experiment), respectively, as described in section 2. In the PAC experiment (Fig. 10d), anomalies from the western Pacific Ocean entered the eastern Indian Ocean through the Indonesian region and then propagated westward across the basin. This experiment reproduces much of the decadal-scale OHC variations in the eastern Indian Ocean simulated by REF, but captures neither the magnitude nor the temporal evolution of OHC variations in the western Indian Ocean. Relative to the PAC experiment, the IND experiment reproduces much more of the OHC variance in the Indian Ocean (Fig. 10e). It is noteworthy that the IND experiment mainly produces OHC anomalies west of 100°E in the Indian Ocean, while the magnitude of OHC anomalies east of 100°E is typically smaller than those produced by the PAC experiment. This result suggests that variations in the eastern Indian Ocean are weakened substantially when variations in the Pacific Ocean are suppressed, and supports our contention that the eastern Indian Ocean is linked to interdecadal variability in the Pacific Ocean primarily via oceanic pathways. The sum of OHC anomalies produced by the PAC and IND experiments (Fig. 10c) closely matches the output from the REF simulation, indicating that the total OHC anomalies in the Indian Ocean (REF) can be regarded as a linear superposition of a locally driven signal (IND) and a remote signal that propagates though oceanic pathways from the Pacific (PAC).
Regression of Ekman pumping onto the IPO index using outputs from the REF simulation (Fig. 11) shows prominent negative correlations over the central south Indian Ocean, with a spatial distribution quite similar to that based on GECCO2 (Fig. 6). This spatial distribution further supports the idea that conditions in the Pacific Ocean can induce Ekman pumping in the central Indian Ocean through the atmospheric bridge. When the IPO is in its positive phase, positive wind stress curl over the central Indian Ocean induces downward Ekman pumping that deepens the thermocline and excites downwelling signals that propagate westward as Rossby waves. The influence of these downwelling Rossby waves is to warm the subsurface layer of the western Indian Ocean. The results of these OCGM sensitivity experiments are consistent with the results based on GECCO2 and the linear Rossby wave model: the eastern boundary condition dominates decadal-scale OHC variability in the eastern Indian Ocean, while the wind-driven Ekman pumping occurring in the central Indian Ocean impacts the western part by triggering westward Rossby waves.
4. Summary and discussion
Decadal variations in the subsurface (50–300 m) OHC in the southern Indian Ocean and their links to the IPO are examined based on reanalysis products and OGCM hindcast simulations. Changes in OHC are closely linked to changes in thermocline depth (D20) on decadal time scales, as shown by correlations exceeding 0.9 over most of the Indian Ocean. Decadal variations in Indian Ocean subsurface OHC can thus be mainly attributed to dynamical fluctuations in the thermocline. The spatial distribution of OHC regressed onto the IPO index closely matches the second EOF mode of OHC, which in turn dominates variations in the Indian Ocean OHC on decadal time scales. This close correspondence suggests that decadal signals in the subsurface OHC of the Indian Ocean are associated with the IPO. As described above, previous studies have already suggested that distinct mechanisms are responsible for the interannual variability of the western and eastern Indian Ocean. Here we add to this body of work by detailing the existence of these distinct mechanisms but on the decadal time scale. In contrast to previous studies of decadal variability in the Indian Ocean that have treated the Indian Ocean as a whole basin, results shown in this study highlight the different mechanisms by which the Pacific Ocean influences the eastern and western parts of the Indian Ocean on decadal time scales.
Examination of reanalysis data indicates that the Pacific influence on decadal variability in the eastern Indian Ocean operates mainly via oceanic pathways. This conclusion is supported by results based on a linear Rossby wave model and OGCM sensitivity experiments, which indicate that decadal variations in SSHA and OHC in the south Indian Ocean east of 100°E mainly enter via the eastern boundary. In the OGCM sensitivity experiment of PAC, which imposes observed surface forcing over the Pacific only, westward-propagating OHC anomalies entering at the eastern boundary of the Indian Ocean dominate OHC variability, especially in the region east of 100°E. These results implicate the Pacific as the source of these anomalies and confirm that anomalies enter the eastern Indian Ocean primarily through oceanic pathways rather than atmospheric forcing. By contrast, the Pacific influence on the western Indian Ocean operates through the atmospheric bridge. Remote variations in the Pacific Ocean induce anomalous wind-driven Ekman pumping in the central Indian Ocean. Rossby waves triggered by this anomalous Ekman pumping transmit these signals westward, displacing the thermocline and altering subsurface OHC in the western Indian Ocean. The linear Rossby wave model indicates that variations in SSHA west of 100°E arise mainly from Ekman pumping, and much of the OHC variability in the western Indian Ocean can be reproduced in an OGCM by specifying the atmospheric forcing over the Indian Ocean alone. In summary, the Pacific Ocean drives decadal variations in the eastern Indian Ocean subsurface OHC via oceanic pathways, while wind-driven Ekman pumping, modulated via the atmospheric bridge from the Pacific, is the key driver of decadal variations in subsurface OHC in the western Indian Ocean.
Although we have emphasized the oceanic connection between the western Pacific and eastern Indian Ocean, local winds also contribute to the decadal variations of OHC in the eastern Indian Ocean. As shown in the IND simulation (Fig. 10e) in which atmospheric forcing is only applied in the Indian Ocean, OHC anomalies in the eastern Indian Ocean are affected by local winds. Comparison of OHC anomalies from the PAC and IND experiments reveals that these wind-driven anomalies have smaller amplitudes than those transmitted from the Pacific via oceanic pathways in the PAC simulations. This result further bolsters the hypothesis that the oceanic pathway plays the dominant role in determining subsurface OHC variation in the eastern Indian Ocean on decadal time scales. It is worth noting that the Indian Ocean local winds also contain the signal from the Pacific, as the tropical Pacific and Indian Ocean are highly coupled through the atmospheric Walker circulation.
The results shown in this paper are based primarily on the GECCO2 reanalysis. Other oceanic gridded observation and reanalyses have been used to assess the robustness of our findings, including the EN4 analysis and the ORAS4 and SODA2.2.4 reanalyses. Figure 12 shows Hovmöller diagrams of 8-yr low-pass-filtered Indian Ocean subsurface OHC anomalies averaged over 10°–20°S during 1960–2004. Five different estimates are shown, including the OGCM REF simulation. Notable discrepancies are evident among these datasets, particularly during the early part of the record when observational data in the Indian Ocean was especially sparse (Nidheesh et al. 2017). Noting the limited number of observed temperature profiles between 100 and 2000 m in the World Ocean Atlas 2001, Harrison and Carson (2007) suggested that sparse data coverage in the Indian Ocean (and especially in the southern Indian Ocean) would pose difficult challenges for decadal climate research in this region. Reliable, long-term oceanic observations are crucial for advancing our understanding of decadal variability in the climate system. Significant strides have been made to improve the availability of hydrographic measurements in the Indian Ocean, including the deployment of the Indian Ocean Observing System (IndOOS), the Argo network, and satellite observations of surface winds and SSH. Despite the discrepancies shown in Fig. 12, many of the key features are evident among all five datasets. These robust features include the discontinuities in subsurface OHC anomalies at 100°E, which we attribute to wind stress curl anomalies in the central Indian Ocean.
To further evaluate the robustness of our results regardless of the choice of dataset, we have forced the linear Rossby wave model using eastern boundary conditions derived from the other two ocean reanalyses (ORAS4 and SODA2.2.4; Figs. 13a,b). The atmospheric forcing for all three simulations is based on ERA-40. Correlation patterns produced by using ORAS4 and SODA are qualitatively consistent with that produced by using GECCO2, with the dominance of the eastern boundary condition in the eastern Indian Ocean and wind stress forcing in the western Indian Ocean. Furthermore, we test the sensitivity of our results to the atmospheric forcing in a similar way, by combining wind stress forcing based on 20CR and R-1, respectively, with eastern boundary conditions derived from GECCO2. Both results show similar correlation patterns: the contribution of wind stress curl dominates variations in SSHA west of 90°E but has little impact in the eastern Indian Ocean. This sensitivity analysis indicates that the conclusions drawn from the linear Rossby wave model results are robust to the choice of the dataset used to force the model. Despite the similarity in correlation results based on different wind data, we note the presence of substantial discrepancies between the three wind stress curl fields. Figure 14 shows meridionally averaged climatological wind stress curl from the three atmospheric reanalyses along with their standard deviations. Relative to R-1 and 20CR, ERA-40 wind stress curl has both larger mean values and larger standard deviations across the basin. Similar discrepancies exist in the zonal wind stress (not shown). This discrepancy in the amplitude of wind stress curl has been previously noted by McGregor et al. (2012) and reasons that account for this discrepancy should be investigated in future studies.
We have primarily used the IPO index to represent the Pacific decadal variability in our analysis. The IPO is characterized by a tripole pattern of SSTA, with three large centers of action in the Pacific Ocean on decadal time scales, namely, the North Pacific, the eastern tropical Pacific, and the South Pacific (Henley et al. 2015). Ummenhofer et al. (2017) suggested that multidecadal variations in the Indian Ocean are associated with the PDO, and Krishnamurthy and Krishnamurthy (2016) highlighted the link of decadal variations of Indian Ocean dipole with the PDO. In addition, the connection between the western South Pacific Ocean associated with IPO and the southeastern Indian Ocean along the western coast of Australia was identified by Vargas-Hernandez et al. (2014). It, therefore, remains unclear which part of the Pacific plays the predominant role in driving Indian Ocean variability. We will further examine the relative influence from these three key regions of the Pacific Ocean on decadal variability in the Indian Ocean in future work.
Use of the following datasets is gratefully acknowledged: GECCO2 from ECCO, ERA-40 and ORAS4 from ECMWF, EN4 from the Met Office, SODA from the SODA/TAMU Research Group, and 20CR and NCEP–NCAR reanalysis from NOAA. OGCM simulations were performed at the Computing Centre in Kiel University. We are grateful for the discussions with Allan Clarke and Susan Wijffels. Comments from three anonymous reviewers and the editor Rong Zhang are gratefully acknowledged. This research was supported by a scholarship from the China Scholarship Council (CSC) to X. J., a research fellowship by the Alexander von Humboldt Foundation to C. C. U., an NSF OCE PO Grant (OCE-1242989) to Y.-O. K., the ONR Young Investigator Award (N00014-15-1-2588) to H. S., and a research grant from the Ministry of Science and Technology of the People’s Republic of China to Tsinghua University (2017YFA0603902).