The ocean heat content (OHC) changes account for more than 90% of the changes in the globe’s energy budget (e.g., Trenberth and Fasullo 2010, 2012) and significantly modulate the pace at which the global mean surface temperature (GMST) is rising. The recent slowdown of GMST warming during the early twenty-first century, dubbed the “hiatus” (e.g., Easterling and Wehner 2009; Meehl et al. 2011; Trenberth and Fasullo 2013), has been attributed to the increased subsurface ocean heat uptake by the majority of the research community (e.g., Chen and Tung 2014; Roemmich et al. 2015; Xie 2016; Liu et al. 2016). Owing to its smaller area than the Pacific and the Atlantic, the importance of the Indian Ocean (IO) in the global ocean heat budget was not recognized until the recent hiatus period, during which the abrupt increase of upper-IO heat content served as a major sink of the excessive heat entering the Earth system (Cheng et al. 2015; Nieves et al. 2015; Lee et al. 2015).

In December 2015, the Paris Agreement of the United Nations Framework Convention on Climate Change (UNFCCC) agreed to pursue efforts to limit the GMST increase within 1.5 K above preindustrial levels (UNFCCC Conference of the Parties 2015). Scientific efforts in understanding OHC changes are helpful for reducing the uncertainties in policy making. In addition, the enhanced warming of the upper IO (Lee et al. 2015; Li et al. 2017) has brought many serious impacts to the regional climate and environment, such as weakening of the South Asian monsoon due to the decreasing land–sea thermal contrast (e.g., Bollasina et al. 2011; Roxy et al. 2015; Singh 2016), reduction in marine primary productivity (Roxy et al. 2016), frequent disastrous marine heatwave and coral bleaching events (Wernberg et al. 2013; Feng et al. 2015; Zinke et al. 2015), and increased stress of sea level rise in low-lying coastal areas and island nations (Han et al. 2010; Nicholls and Cazenave 2010). Consequently, quantifying and understanding the upper-IO heat content variations in the observed history are of scientific and practical significance.

Decadal (~10 yr)-to-multidecadal (>10 yr) variations of the upper IO are controlled by complex mechanisms involving local climate processes versus remote forcing from other ocean basins, dynamic versus thermodynamic processes, and natural climate variability versus anthropogenic effects (e.g., Allan et al. 1995; Schott et al. 2009; Han et al. 2014). Among the various factors, the interdecadal Pacific oscillation (IPO; Power et al. 1999) exerts a strong influence on the IO. Wind anomalies associated with the IPO drive oceanic waves that are able to transmit through the Indonesian channels (Clarke and Liu 1994; Meyers 1996), modulating the heat transport of the Indonesian Throughflow (ITF; Wijffels and Meyers 2004; Cai et al. 2005; Nidheesh et al. 2013; Sprintall et al. 2014) and causing upper-ocean variations in the south IO (Feng et al. 2004, 2013, 2015; Schwarzkopf and Böning 2011; Li et al. 2017; Ummenhofer et al. 2017). The negative IPO condition, with La Niña–like sea surface temperature (SST) anomalies and enhanced trade winds over the tropical Pacific, strengthens the ITF heat transport and results in rapid warming of the IO (England et al. 2014; Nieves et al. 2015; Lee et al. 2015), and meanwhile leads to decelerated GMST warming, such as for the cases of the “big hiatus” of the 1940–70s and the early twenty-first-century hiatus (Kosaka and Xie 2013, 2016; Meehl et al. 2014, 2016). The positive IPO condition, with frequently occurring strong El Niños, leads to cooling of the IO and accelerated GMST warming, such as for the period of the 1970s–90s.

Besides the oceanic connection, climate variability in other ocean basins can also affect the IO via the atmospheric bridge. El Niño–Southern Oscillation (ENSO) and IPO can induce large-scale anomalies of surface winds, heat fluxes, and rainfall over the IO (Klein et al. 1999; Alexander et al. 2002; Lee and McPhaden 2008). Meanwhile, climate modes of the Atlantic Ocean such as the Atlantic multidecadal oscillation (AMO) may also modulate the tropical IO climate via atmospheric teleconnection (Li et al. 2016). These external impacts are entangled with the IO local climate variability, with both causing variability in the upper IO (Lee 2004; Alory et al. 2007; Han et al. 2010; Trenary and Han 2013; Nieves et al. 2015). Regardless of different origins, in this study we do not distinguish the two forcing effects and collectively refer to them as “local atmospheric forcing effects within the IO.” Recent studies suggest that local atmospheric forcing effects, particularly those associated with natural internal climate variability, are responsible for most of the decadal variations of the upper IO (Han et al. 2006, 2010, 2017b; Trenary and Han 2013; Li and Han 2015). Stochastic atmospheric forcing, such as the tropical atmospheric intraseasonal oscillations (ISOs) and midlatitude storm tracks (Sura et al. 2006), can also cause decadal ocean variability. This effect, however, is shown to be important only in the subtropical southwest IO (Li and Han 2015). Anthropogenic forcing is another origin for the IO heat content variability, particularly for the multidecadal trend. The increased anthropogenic greenhouse gas concentration causes SST warming of the IO through surface longwave radiation (Du and Xie 2008; Dong et al. 2014; Zhang et al. 2018). Meanwhile, the increased anthropogenic aerosols over Asia in the past several decades and aerosols from volcanic eruptions (Gleckler et al. 2006; Smith et al. 2016) can cause SST cooling of the IO (Cai et al. 2006, 2007; Dong and Zhou 2014).

Besides atmospheric forcing, ocean internal instabilities, manifesting as mesoscale eddies and current meanders, can also cause significant decadal-to-multidecadal variability in some regions because of nonlinearity of the oceanic system. Modeling studies suggest that in eddy-rich areas of the IO—for example, the subtropical south IO, the Somali Current region, and the western Bay of Bengal—low-frequency oceanic variability induced by ocean internal instability is comparable in strength to that induced by atmospheric forcing (e.g., Wirth et al. 2002; Trenary and Han 2013; Li and Han 2015; Sérazin et al. 2016). Since ocean internal instabilities are generally unpredictable for ocean/climate models, they cause significant uncertainties for decadal ocean prediction.

Given the complexity of the problem, a thorough investigation is imperative for quantifying and understanding the upper-IO heat content variations in the observed history and thereby providing theoretical support for future prediction, policy making, and adaptation of human societies. This effort is pursued in the present study by analyzing historical ocean observational data and performing ocean general circulation model (OGCM) experiments that are capable of isolating effects from different processes. Our results show strikingly strong multidecadal changes of the upper-IO heat content. The decadal trends of the IO heat content in 0–400 m (OHC₄₀₀) reached as large as O(10²²) J decade⁻¹, on the same order of the global OHC changes (Chen and Tung 2014; Roemmich et al. 2015; Cheng et al. 2017) and are potentially influential for global climate change. The roles played by major processes, such as ITF heat transport due to Pacific wind forcing, local wind forcing, and radiation forcing within the IO, are examined.

The remainder of the paper is arranged as follows. In section 2, we introduce observational datasets, describe the OGCM utilized in this study, and design a hierarchy of OGCM experiments. In section 3, we describe the observed and simulated multidecadal OHC₄₀₀ variations and verify the performance of the OGCM. In section 4, we expound the mechanisms controlling the multidecadal OHC₄₀₀ variability by assessing the relative importance of different processes. In section 5, we provide a summary and discussion for our primary findings.

To estimate OHC₄₀₀ variability, four ocean subsurface temperature datasets are analyzed and compared: 1) the World Ocean Atlas (WOA; Levitus et al. 2009) during 1955–2016 provided by the National Oceanic and Atmospheric Administration (NOAA)’s National Centers for Environmental Information (NCEI) of the United States, 2) the Ishii data (Ishii et al. 2005) during 1955–2012 from the Frontier Research System for Global Change, the Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 3) the EN4.2.0 data (Good et al. 2013) during 1955–2016 provided by the Met Office, and 4) the Institute of Atmospheric Physics (IAP) ocean temperature analysis (Cheng et al. 2016, 2017) during 1955–2015 provided by the IAP, Chinese Academy of Sciences (CAS). All the four products provide 1° × 1°, standard-level, annual temperature estimates in the upper 700 m. Since large discrepancies exist among the four products in some regions (Wang et al. 2018), we consider the consensus of the four products as the “observational fact” and their spread as observational uncertainty. In addition, two versions of Simple Ocean Data Assimilation (SODA) products (Carton et al. 2005; Carton and Giese 2008) are also analyzed. SODA2.2.4 (SODA2) is forced with NOAA Twentieth Century Reanalysis, version 2 (20CR V2; Compo et al. 2011), and has an average resolution of 0.25° × 0.4° × 40 levels, covering 1871–2010, while SODA3.3.1 (SODA3) is forced by Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2; Rienecker et al. 2011), and has a resolution of 0.25° × 0.25° × 50 levels, covering the span of 1980–2015.

OHC₄₀₀ is used to quantify the upper-IO heat uptake, and for each dataset it is computed as

To evaluate the OGCM solution, we also used 0.25° × 0.25° monthly satellite sea level data based on the satellite altimeter sea surface height (SSH) product from the Archiving, Validation, and Interpretation of Satellite Oceanographic Data (AVISO) from January 1993 through August 2016 (Le Traon et al. 1998), and 1° × 1° monthly latent and sensible heat fluxes of the objectively analyzed air–sea fluxes (OAFlux) from 1958 through 2016 (Yu and Weller 2007). To quantify the IPO-relevant climate variability and examine its relationship with the Indian OHC₄₀₀ changes, the IPO index is computed with the monthly Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST) from 1871 through 2016 (Rayner et al. 2003), adopting the “tripole” index method of Henley et al. (2015), that is, the SST anomaly difference between the equatorial Pacific (10°S–10°N, 170°E–90°W) and the northwest plus southwest Pacific region (140°E–145°W, 25°–45°N plus 150°E–160°W, 50°–15°S). In addition to ocean datasets, five surface wind reanalysis products are analyzed and compared: 1) the 0.75° European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim) from 1979 to 2016 (Dee et al. 2011), 2) the 1° ECMWF Twentieth Century Reanalysis (ERA-20C) from 1955 to 2010 (Poli et al. 2016), 3) the 1.875° National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) from 1955 to 2016 (Kalnay et al. 1996), 4) the 0.5° × 0.625° MERRA-2 winds from 1980 to 2016, and 5) the 1° Climate Forecast System Reanalysis (CFSR) from 1979 to 2016 (Saha et al. 2010).

The OGCM used in this study is the Hybrid Coordinate Ocean Model (HYCOM), version 2.2.18 (Halliwell 2004), configured to the Indo-Pacific Ocean basin between 19°E–68°W and 55°S–50°N (see Fig. S1 in the online supplementary material). The zonal resolution is ⅓° between 70° and 170°E and gradually changes to 1° west of 40°E and east of 160°W, while the meridional resolution is ⅓° at low latitudes (25°S–25°N) and gradually changes to 1° at midlatitudes (poleward of ~44°). The model has 35 hybrid vertical layers, with a top layer thickness of 5 m. On the open-ocean boundaries of the model domain, 5° sponge layers are applied wherein model temperature and salinity are relaxed to World Ocean Atlas 2013 (WOA13) monthly climatology (Locarnini et al. 2013; Zweng et al. 2013). The surface atmospheric forcing fields include surface wind stress, 10-m wind speed, shortwave and longwave radiation (SWR and LWR), precipitation, and 2-m air temperature and humidity, which are taken from ERA-20C and ERA-Interim. Monthly climatologic river discharges of Dai et al. (2009) are used as lateral freshwater flux forcing. To suppress the drifting of model salinity, sea surface salinity is relaxed to WOA13 climatology with an e-folding time of ~135 days. Wind stress and wind speed are imposed separately on the model ocean and affect ocean temperature through different processes. Wind stress drives ocean dynamics (advection and upwelling) and mixing, while wind speed affects ocean temperature mainly through surface turbulent heat flux (latent plus sensible heat fluxes). In our model, turbulent heat fluxes and evaporation are calculated online with wind speed, air temperature, specific humidity, and modeled SST using the Coupled Ocean–Atmosphere Response Experiment, version 3.0 (COARE 3.0) algorithm (Fairall et al. 2003; Kara et al. 2005).

The model was spun up for 30 years, using monthly climatologic fields of ERA-20C as surface forcing. Subsequent to the spinup run, the model is integrated forward from 1940 to 2016 in three parallel experiments. The “main run” (MR) is forced by daily ERA-20C fields for 1940–2010 and daily ERA-Interim fields for 2011–16. ERA-20C and ERA-Interim are highly consistent with each other in surface winds during their overlapping period but show discrepancies in surface radiations (figures not shown). To reduce the effect of discontinuity in radiation forcing, LWR and SWR of ERA-Interim are tuned toward those of ERA-20C in annual mean climatology. MR contains the complete processes of the Indo-Pacific Ocean and is compared with observational data to evaluate the model performance.

The second experiment is the “Pacific run” (PAC), in which daily forcing is retained in the Pacific (the same as in MR), whereas monthly climatologic forcing (retaining only seasonal cycle) is imposed on the IO (Fig. S1b). There is a 5° transition zone where daily forcing outside the IO gradually changes to climatologic forcing. As such, the IO variability in the PAC is primarily induced by the Pacific Ocean variability through the ITF. The Southern Ocean may also affect the IO, but this effect is small in PAC because only a small part of the Antarctic Circumpolar Current is included in our model configuration. The difference between MR and PAC, that is, MR − PAC, can roughly represent the total effect of local atmospheric forcing within the IO (including both local climate variability and atmospheric teleconnections from other basins). It should be stated that the Pacific forcing effect and local forcing effect may not be clearly separable in such linear sense. As one of the reviewers pointed out, a more solid method is to perform an “Indian Ocean run” to assess the local forcing effect and analyze its results with those of MR and PAC to assess separately the two effects along their interaction. This has not been done because of the availability of computation resources.

The third experiment, the “wind stress run” (TAU), provides insights from another perspective. TAU is forced with daily wind stress, whereas other forcing fields, such as wind speed, radiation, and precipitation, are fixed to monthly climatology. As such, TAU roughly represents the effect of wind stress–driven ocean dynamics (hereafter “wind-forcing effect”), and the difference between MR and TAU, that is, MR − TAU, mainly assesses the effect of surface heat fluxes (SWR, LWR, and turbulent heat flux) on the IO heat content. Outputs of the three experiments during 1940–2016 were stored in monthly mean resolution and further averaged into annual mean fields for our analysis. Isopycnal-coordinate HYCOM data were interpolated to 10-m-depth intervals to facilitate the computation of OHC₄₀₀.

This section explores the physical processes causing multidecadal OHC₄₀₀ variations. First, we examine the relative importance of remote forcing from the Pacific via the ITF versus local forcing over the IO (Figs. 4a,b). PAC (Figs. 4a,b, blue) explains a large portion of the decadal OHC₄₀₀ variations in MR, accounting for ~90% of the total OHC₄₀₀ increase of 1965–79, and showed even larger decreasing trend during 1980–96. However, during 2000–14, the Pacific-origin warming explained merely ~44% of the total, and it was slightly larger than that of 1965–79 (9.1 vs 8.1 × 10²¹ J decade⁻¹). On the other hand, the difference between MR and PAC (MR − PAC) quantifies local atmospheric forcing effect in the IO. This effect had little contribution to the first two periods but played the major role (~56%) in the post-2000 warming. This means that the mechanism is likely varying with time, with increasing importance of IO local processes in the twenty-first century.

View larger version (45K) Fig. 4. (a) Annual Indian OHC400 and (b) its linear trends for 1965–79, 1980–96, and 2000–14, derived from MR and PAC experiments of HYCOM and their difference MR − PAC. (c),(d) As in (a),(b), but derived from MR, TAU, and MR − TAU.

To assess the effects of wind versus heat flux forcing, we analyze the solutions of TAU and MR (Figs. 4c,d). The wind stress–driven ocean dynamics, as measured by TAU, play a fundamental role in modulating OHC₄₀₀ variability. Their contributions are large for all three periods. Particularly during 1965–79 and 2000–14, TAU explained more than half of the total OHC₄₀₀ increasing trend. MR − TAU roughly represents the heat flux forcing effect, which also made large contributions. The wind stress forcing effect exceeds that of heat flux during the two increasing periods, whereas heat flux forcing has a slightly larger contribution during the 1980–96 cooling period.

We further examine the contributions of different processes to the trend pattern (Fig. 5). During 1965–79, strong warming occurred in the southeast IO near the exit of the ITF, and weaker warming spread over the equatorial and north IO. Afterward, PAC exerted basinwide cooling and warming effects on the central and north IO during 1980–96 and 2000–14, respectively. All these trends were significant at 95% confidence level, whereas the subtropical south IO was exceptional. This region is largely affected by ocean internal variability and will be discussed later. Compared with the Pacific forcing, the IO local processes (MR − PAC) had a smaller contribution to the basin mean trend, partly because of the compensation between regional warming and cooling signals. However, local forcing induced strong warming over the western tropical IO and the southeast IO during 2000–14. Note that T₄₀₀ trends of MR − PAC during 2000–14 greatly resemble those of MR (Fig. 2f), confirming the dominant role of local processes during this period.

View larger version (137K) Fig. 5. T400 trends (K decade−1) during (left) 1965–79, (center) 1980–96, and (right) 2000–14 derived from PAC, MR − PAC, TAU, and MR − TAU. Stippling indicates exceeding 95% significance.

TAU includes both local wind-forcing effect within the IO and remote wind-forcing effect from the Pacific. Signatures of wind forcing include warming signals in the southwest tropical IO and southeast IO during 1965–79, cooling of the central-to-eastern IO during 1980–96, and enhanced warming of the western-to-central tropical IO after 2000. Heat flux forcing (MR − TAU) primarily affected the subtropical south IO during the first two periods, whereas after 2000, it induced enhanced warming throughout the IO Basin.

In almost all the HYCOM experiments, trends in the subtropical south IO tend to have large magnitudes but low statistical significance. These features conform to ocean internal variability. Existing research has pointed out the active mesoscale eddy variability in this region and its large impact on low-frequency ocean variations (e.g., Trenary and Han 2013; Li and Han 2015; Sérazin et al. 2016). This is proven by the results of the spinup run that used monthly climatologic forcing and therefore have no atmospheric forcing effect at interannual and lower frequencies (Fig. 6). The upper ocean has reached a quasi-steady state since the 10th year, and the standard deviation map of annual T₄₀₀ during the 11th–30th years (Fig. 6a) clearly suggests the ocean internal origin for the large interannual-to-decadal variability in the subtropical south IO. The decadal T₄₀₀ trend during the 21st–30th years shows strong signatures in the subtropical south IO (Fig. 6b), large enough to explain the trend signals in MR. These results raise caution for the interpretation of our HYCOM results, given that ocean internal variability is chaotic in nature and its phase cannot be reproduced by model simulation. In the subtropical south IO, MR showed strong warming trends during 1965–79 and 2000–14 and cooling trends in 1980–96, which were, however, not seen in PAC or TAU. It means that ocean internal variability happened to amplify the Indian OHC₄₀₀ trends of MR in all three periods, and this could be a major cause for the overestimated OHC₄₀₀ changes in MR relative to observation. It is, however, unknown how the ocean internal variability has affected the upper-IO variability in reality.

View larger version (48K) Fig. 6. (a) Standard deviation of annual T400 from the 11th to the 30th year of the HYCOM spinup run. (b) Linear trend of annual T400 from the 21st to the 30th year of HYCOM spinup run.

The above analysis broadly assesses the relative importance of different processes, and hereafter we further examine how each of them affects the upper IO, and their relationship with well-known climate phenomena such as the IPO and anthropogenic warming. First, we look at the Pacific forcing through the ocean channel (the ITF). Figure 7a compares ITF heat transport and temporal tendency of OHC₄₀₀ (∂OHC₄₀₀/∂t) in PAC. The two variables show good agreement, with a linear correlation of r = 0.71, confirming the dominance of the ITF in the Indian OHC₄₀₀ variability in PAC. Strong and weak ITF years (Fig. 7, red and blue dots) were identified with a 0.5 standard deviation criterion. There were more strong ITF years in the two warming periods and more weak ITF years in the 1980–96 cooling period. To better understand the ITF variability, Fig. 7b further compares the ITF heat transport estimates based on MR, PAC, and TAU. It is found that the three are highly consistent, showing merely detailed differences. A similar situation is seen in volume transport and mean temperature of the ITF (Fig. S3). These results indicate that the ITF variability is primarily dictated by the Pacific wind forcing. The mean temperature of the ITF water also exhibited pronounced decadal variations. The ITF water was evidently warmer during the two warming periods and thus contributed to the heat transport changes of the ITF, as suggested by Lee et al. (2015).

View larger version (73K) Fig. 7. (a) Normalized Indian OHC400 tendency (∂OHC400/∂t) and the ITF heat transport in PAC. ITF heat transport values higher (lower) than the +0.5 (−0.5) standard deviation value are plotted as red (blue) dots, representing the strong (weak) ITF years. (b) Heat transport (1 PW = 1015 W) of the ITF derived from MR, PAC, and TAU. (c) Normalized tropical Pacific τx (150°E–120°W, 10°S–10°N), IPO index, and ITF heat transport. The mean values of τx (IPO index) for the three periods are plotted as dashed black (solid blue/red) straight lines. Here the ITF transport is computed at 5°S by integrating the meridional current over 106°–136°E and 0–400 m, and southward transport (toward the IO) is defined positive. Each variable is shown as the anomaly relative to the 1980–2010 mean value.

The results of Fig. 7 point to the need to understand the Pacific wind variations. The mean zonal wind stress τ^x anomalies of the three periods are displayed in Fig. S4. The tropical Pacific was dominated by anomalous easterly winds (negative τ^x anomalies) in the two warming periods, indicative of a strengthening of the tropical trade winds, whereas during 1980–96 anomalous westerly winds suggested a weakening of the trade winds. These changes in Pacific trade winds were closely associated with the IPO (Fig. 7c). The correlation between the IPO index and the tropical Pacific–averaged τ^x is 0.86. The two warming periods were generally under negative IPO conditions with more La Niña events and enhanced trade winds, while the 1980–96 period was characterized by a positive IPO condition with frequently occurring El Niños and weakened trade winds. The enhanced (weakened) Pacific trade winds during the negative (positive) IPO condition lead to strengthened (attenuated) ITF heat transport toward the IO during 1965–79 and 2000–14 (1980–96). The correlation between tropical Pacific τ^x and ITF heat transport is r = −0.63, also exceeding 95% significance.

Existing studies have demonstrated that the ITF variability mainly affects the south IO and has much weaker signatures in the north (Feng et al. 2004; Zhuang et al. 2013; Trenary and Han 2013; Dong and McPhaden 2016; Li et al. 2017). This is confirmed by composite maps of ∂T₄₀₀/∂t for the strong and weak ITF years (Figs. 8a,b). A stronger (weaker) ITF directly induces warming (cooling) of the south IO, with hotspots near the exit of the ITF and in the subtropical southeast IO. Moreover, annual ∂OHC₄₀₀/∂t of the south IO (south of 5°S) shows consistent variability with the ITF heat transport in PAC (Fig. 9a), showing a correlation of r = 0.73. However, Fig. 5 shows that in PAC, the ITF is able to induce basinwide decadal warming/cooling over the north IO, although the trend signals are weaker than those in the south IO. The multidecadal warming/cooling and the corresponding sea level rise/fall in the north IO have been observed and reported by several recent studies but primarily attributed to local wind forcing, such as equatorial and monsoonal winds (Thompson et al. 2016; Srinivasu et al. 2017; Swapna et al. 2017). Therefore, a significant impact of the ITF over the north IO is quite new for our existing context of knowledge and requires careful examination.

View larger version (107K) Fig. 8. (a),(b) Composite ∂T400/∂t maps for the (a) weak and (b) strong ITF years derived from PAC. Weak and strong ITF years are shown in Fig. 7a. (c),(d) Composite ∂T400/∂t for two years after the weak and strong years, also derived from PAC. Stippling indicates exceeding 95% significance.

View larger version (75K) Fig. 9. (a) Annual ∂OHC400/∂t of the south IO (south of 5°S) and the ITF heat transport in PAC. (b) Annual ∂OHC400/∂t of the south IO and the north IO (north of 5°S). (c) Annual ∂OHC400/∂t of the north IO and meridional heat transport at 5°S (northward transport defined positive) in the IO for 0–400 m. Each time series is normalized by its standard deviation.

For the first step, we answer the question as to how the ITF affects the north IO in PAC. Figure 9b compares the ∂OHC₄₀₀/∂t time series of the north and south IO in PAC, divided by the 5°S latitude. A lagged relationship between them is clearly discerned. The two show low concurrent correlation (r = −0.11; insignificant at 95% confidence level) but high correlation when the south IO leads by two years (r = 0.86). This means that the north IO shows a lagged response to the ITF variability. The ITF is dynamically linked to the ocean circulation of the upper IO, particularly the shallow overturning circulation (Schott et al. 2002). Liu et al. (2016) showed that the IPO can modulate the cross-equatorial overturning cell (their Fig. 6b), which may cause anomalous upper-ocean heat exchanges between the south and north IO. To test this hypothesis, we compute the heat transport at 5°S over the IO Basin (Fig. 9c), which shows rather consistent variations with the north IO ∂OHC₄₀₀/∂t (r = 0.56; significant at 95% confidence level). This effect led to strong north IO warmings in 1976, 2001/02, and 2008–12, and coolings in 1966–68, 1984–93, and 2004/05. A further check over the 5°S section reveals that the heat transport mainly occurs in the western boundary region rather than in the interior (Fig. S5). The western boundary current at 5°S, namely, the East African Coastal Current (Schott et al. 2009), is not significantly enhanced by the ITF. Instead, the water transported by the western boundary current is warmer, particularly in the thermocline (Fig. S5c). This effect leads to a basinwide cooling and warming signatures north of 5°S, as shown by the composite ∂T₄₀₀/∂t for two years after the weak and strong ITF years, respectively (Figs. 8c,d). In the south IO, the tendencies that occurred two years later can partly offset those that occurred in strong and weak ITF years when computing the decadal trend. The 2-yr phase lag is short in comparison with the three periods considered here (>14 years), and therefore basinwide north IO warming and cooling are seen in trend maps of PAC (Fig. 5). This also explains why the trend pattern of 2000–14 shows fundamental differences from the 2003–12 one revealed by Li et al. (2017) in which the southeast IO has much stronger warming than other areas. The 2000–14 period includes the lagged effects of the 1998–2000 and 2010/11 La Niña conditions that induce north IO warming and southeast IO cooling as in Fig. 8d.

A further question naturally arises as to whether the processes elucidated above are realistic or artificial. Here we are able to gain some insights by looking at MR results and observational data. Figure S6 demonstrates that in MR, despite local atmospheric forcing effects, the lagged correlation between the south and north IO still holds true (r = 0.67), and the heat transport at 5°S strongly affects the north IO (r = 0.41; significant at 95% confidence level). In the composite maps of MR, the overall warming/cooling of the north IO are also seen two years after the ITF extremes (Figs. S7a,b), although the spatial patterns are not quasi-uniform because of the modifications of local forcing. These results indicate that the ITF’s effects on the north IO can be tested with observations. Figures S7c and S7d show the observed composite ∂T₄₀₀/∂t maps two years after the weak and strong ITF years, respectively, which are basically consistent with those of MR, albeit with prevailingly weaker magnitudes. Furthermore, the lagged correlation between the south and north IO also holds in observation (Fig. S8a), showing a correlation of r = 0.32 when the south IO leads by two years (exceeding 95% significance, assuming each year has independent degree of freedom). The observed north IO OHC₄₀₀ agrees quite well with that in MR (Fig. S8b; r = 0.68, significant at 99%). Recalling Fig. 3d suggests that the HYCOM-simulated SSH in the Bay of Bengal is highly consistent with altimeter observation. If the ITF’s large effect is an artifact of the model, these favorable comparisons with observation cannot be achieved. With these regards, the notable effect of the ITF on the north IO is likely realistic.

Given that the ITF heat transport variability is mainly controlled by the Pacific wind forcing, the difference between TAU and PAC (TAU − PAC) can approximately represent the effect of local wind stress forcing on the upper-IO heat content. Note that here the “local wind stress forcing” contains the atmospheric bridge effect from the Pacific. A weak cooling over the equatorial and north IO and a warming of the southwest tropical IO are seen during 1965–79 (Fig. 10a). The trend patterns were quite similar in the subsequent two periods (Figs. 10b,c), both showing widespread warming over the western tropical IO and weak cooling in the eastern equatorial IO and the eastern Bay of Bengal. During 2000–14, enhanced and significant warming trends were seen in the Arabian Sea and the southwest tropical IO. To better understand these trend structures, we also computed SSH trends (Figs. 10d–f) that show similar distributions to T₄₀₀ trends over the tropical IO in all three periods. The pattern correlations between T₄₀₀ and SSH trends are 0.69, 0.64, and 0.67 for the three periods. The consistency between SSH and T₄₀₀ indicates that the upper-ocean warming/cooling trends are a quasi-adiabatic ocean response to surface wind stress forcing, in the form of upper-ocean convergence/divergence.

View larger version (118K) Fig. 10. T400 trends during (a) 1965–79, (b) 1980–96, and (c) 2000–14, derived from TAU − PAC. (d)–(f) SSH trends derived from TAU −PAC. (g)–(i) Surface wind stress trends (black vectors; 10−4 N m−2 decade−1) and Ekman pumping velocity wE trends (color shading; cm day−1 decade−1) derived from ECMWF winds. Stippling in (a)–(g) indicates exceeding 95% significance.

Equatorial zonal wind stress and off-equatorial Ekman pumping are both important in driving oceanic adjustments, and their trends are shown in Figs. 10g–i. Here, Ekman pumping velocity w_E = curl(τ/f)ρ₀⁻¹, where f is the Coriolis parameter and ρ₀ = 1025 kg m⁻³ is the mean seawater density of the Ekman layer. During 1965–79, westerly and easterly wind trends were observed in the western and eastern equatorial basins, respectively. Together with the westerly wind trends in the south IO, these surface wind changes induced strong upwelling (positive w_E trend) in the off-equatorial region in the central-to-eastern IO, which was the primary cause for the tropical IO cooling. In the subsequent two periods, easterly wind trends persisted near the equator, causing the eastern IO cooling through driving a shoaling thermocline (upwelling), and causing the western IO warming by deepening the thermocline, and therefore setting up the east–west temperature and SSH gradients. On decadal time scale, the east–west SSH gradient is in quasi equilibrium with the zonal wind forcing (Han et al. 2017a), given that the decadal period is much longer than the equatorial adjustment time over the Indian Ocean, which is the time it takes for the equatorial Kelvin wave to cross the basin and the first meridional mode Rossby wave to return. It is approximately six months for the second baroclinic mode and even shorter for the first baroclinic mode (e.g., Han et al. 1999, 2011). The decadal period is also much longer than the adjustment time for the north IO, which is about 17 months (Shankar et al. 2010). This quasi-equilibrium balance between equatorial winds and zonal SSH gradient is set up by radiating transient equatorial Kelvin waves eastward [and subsequently poleward via coastal Kelvin waves to the Bay of Bengal and Sumatra–Java coasts, with part of their energy dissipated by mixing (e.g., Clarke and Liu 1994; Han and Webster 2002)] and westward-propagating Rossby waves. The easterly winds were also partly responsible for the western IO warming, through exciting downwelling Rossby waves. Propagation of equatorial waves and their close relationship with equatorial winds are shown in Fig. 11a.

View larger version (73K) Fig. 11. (a) Time–longitude plot of 8-yr low-pass-filtered SSH anomaly (cm) derived from TAU − PAC in the equatorial IO (3°S–3°N). Black (green) contours denote the positive (negative) 8-yr low-passed τx anomaly (10−3 N m−2) derived from ECMWF winds. (b) As in (a), but in the tropical south IO (15°–5°S), and black (green) contours denote positive (negative) 8-yr low-passed wE (cm day−1).

The equatorial zonal winds were meridionally bounded by off-equatorial Ekman pumping. Positive (negative) w_E trends drove westward-propagating upwelling (downwelling) Rossby waves in the tropical south IO and led to persistent SSH fall (rise) and upper-ocean cooling (warming) in the southwest tropical IO (Fig. 11b). Coastal winds were important in forming some regional trend structures through coastal downwelling/upwelling. During 2000–14, the northeasterly winds along the west Asia–East Africa coast were favorable for the generation of local warming maximum in the western Arabian Sea, while in the southeast IO, the northerly winds induced strong warming along the West Australian coast (Li et al. 2017).

It is necessary to examine the relationship between the IO surface winds and the IPO. The composite for negative IPO years shows significant westerly anomalies in the equatorial region (Fig. S9a), as in the cases of 1965–79 and 2000–14 (Fig. S4). The IPO also induces westerly winds in the south IO, which were seen in 1965–79 (Fig. S4a) but not in 2000–14 (Fig. S4c). These zonal wind anomalies also induce strong w_E anomalies in the off-equatorial tropics (Fig. S9b) that are critical for the redistribution of OHC₄₀₀ changes.

We directly compare surface wind anomalies with the IPO index in three key regions (Fig. 12). Equatorial zonal winds and tropical south IO w_E from ECMWF data show moderate correlations with IPO (−0.40 and −0.25). Most of IPO maxima (minima), indicating El Niño (La Niña) events, can cause equatorial easterly (westerly) winds and Ekman downwelling (upwelling) in the tropical south IO. These results are in line with existing understanding that ENSO and IPO have large atmospheric teleconnection imprints on the tropical IO (e.g., Klein et al. 1999; Alexander et al. 2002; Lee and McPhaden 2008; Han et al. 2017a). There are also variations not explained by IPO, probably arising from local variability in the IO such as the Indian Ocean dipole events independent of ENSO (Saji et al. 1999). In comparison, the alongshore winds in the southeast IO (along the West Australian coast) show larger uncertainties in the relationship with IPO (r = 0.08). Recent studies have pointed out a strong control of ENSO/IPO in this region, with El Niño/positive IPO (La Niña/negative IPO) driving anticyclonic (cyclonic) anomalous winds over the southeast IO, which involve southerly (northerly) winds along the West Australian coast (Tozuka et al. 2014; Feng et al. 2015; Li et al. 2017). This relationship is, in fact, approximately valid in some periods, particularly after 1990. The low correlation might be attributed to the time-varying relationship between ENSO/IPO and the IO atmospheric circulation, or else to the large impact of atmospheric internal variability. NCEP–NCAR winds are highly consistent with ECMWF winds in the equatorial IO (Fig. 12a) but show large discrepancies in the south IO (Figs. 12b,c), particularly for the 1960s–70s when satellite data were unavailable for constraining the assimilation systems. Different wind products (ERA-Interim, ERA-20C, MERRA-2, NCEP–NCAR, and CFSR) show broadly consistent variations after 1979 (Fig. S10), confirming the relationships revealed above.

View larger version (73K) Fig. 12. (a) Zonal wind stress τx in the equatorial IO (60°–95°E, 5°S–5°N), (b) Ekman pumping velocity wE in the tropical south IO (50°–90°E, 15°–8°S), and (c) meridional wind stress τy along the West Australian coast (110°–120°E, 30°S–15°N), derived from ECMWF and NCEP–NCAR wind data. The IPO index is also plotted as gray bars.

In addition to wind forcing, surface heat flux is another important process causing OHC₄₀₀ variability. In HYCOM simulations, surface SWR and LWR are taken from ECMWF data, while surface turbulent heat flux Q_T is computed internally in the model with bulk formula. As referenced to OAFlux data, Q_T in MR shows realistic spatial distribution and temporal variability over the Indian Ocean (Fig. S11). We first examine the heat flux forcing effect on the basin-integrated OHC₄₀₀. The net surface heat flux Q_net can well explain the Indian ∂OHC₄₀₀/∂t in MR − TAU, showing a correlation of r = 0.71 (Fig. 13a). The three components of Q_net, that is, SWR, LWR, and Q_T, played distinct roles in the three periods (Fig. 13b). The high Q_net during 1965–75 was mainly induced by large SWR toward the ocean. In the subsequent two periods, SWR was greatly reduced, probably owing to the increasing anthropogenic aerosols over Asia (e.g., Cai et al. 2006, 2007); eruptions of El Chichón in 1982 and Mount Pinatubo in 1991 also contributed to the low SWR during 1980–96. In sharp contrast, because of the rise in greenhouse gas concentration, surface LWR exhibited a monotonic increase over the three periods and acted as the major driver of the enhanced warming during 2000–14. This is also consistent with the larger IO local forcing effect in this period. Although the effect of Q_T was overall much weaker than SWR and LWR, it contributed to the cooling of 1980–96.

View larger version (76K) Fig. 13. (a) Normalized ∂OHC400/∂t and surface net heat flux Qnet integrated over the IO, derived from MR − TAU; dashed lines denote their mean values for the three periods. (b) Anomalies of the IO-integrated Qnet, SWR, LWR, and turbulent heat flux QT (1 TW = 1012 W) averaged over the three periods. SWR and LWR are based on merged ERA-20C and ERA-Interim data that are used to force HYCOM, while Qnet and QT are derived from MR − TAU. (c)–(e) Distributions of Qnet anomaly for 1965–79, 1980–96, and 2000–14 derived from MR − TAU. Stippling indicates exceeding 95% significance.

Surface heat fluxes are also important in determining the warming and cooling patterns. For example, the obvious and significant warming in the south IO during 1965–79 in observation and MR (Fig. 2) is seen in the MR − TAU solution (Fig. 5). Figure 13a confirms increased Q_net in the south IO during this period. Checking the distributions of the three components suggests that this increased Q_net arose primarily from SWR and secondarily from Q_T (Fig. S12). The Q_net had basinwide cooling and warming effects on the tropical IO for 1980–96 and 2000–14, respectively, and LWR was the leading contributor; Q_T also contributed to the enhanced warming of the tropical IO during 2000–14. A comparison with IPO composite maps (Fig. S13) indicates that Q_T anomalies in the three periods contain evident signatures of the IPO, whereas SWR and LWR anomalies are not likely associated with the IPO.