Decadal variability of the subsurface ocean heat content (OHC) in the Indian Ocean is investigated using a coupled climate model experiment, in which observed eastern tropical Pacific sea surface temperature (EPSST) anomalies are specified. This study intends to understand the contributions of external forcing relative to those of internal variability associated with EPSST, as well as the mechanisms by which the Pacific impacts Indian Ocean OHC. Internally generated variations associated with EPSST dominate decadal variations in the subsurface Indian Ocean. Consistent with ocean reanalyses, the coupled model reproduces a pronounced east–west dipole structure in the southern tropical Indian Ocean and discontinuities in westward-propagating signals in the central Indian Ocean around 100°E. This implies distinct mechanisms by which the Pacific impacts the eastern and western Indian Ocean on decadal time scales. Decadal variations of OHC in the eastern Indian Ocean are attributed to 1) western Pacific surface wind anomalies, which trigger oceanic Rossby waves propagating westward through the Indonesian Seas and influence Indonesian Throughflow transport, and 2) zonal wind anomalies over the central tropical Indian Ocean, which trigger eastward-propagating Kelvin waves. Decadal variations of OHC in the western Indian Ocean are linked to conditions in the Pacific via changes in the atmospheric Walker cell, which trigger anomalous wind stress curl and Ekman pumping in the central tropical Indian Ocean. Westward-propagating oceanic Rossby waves extend the influence of this anomalous Ekman pumping to the western Indian Ocean.
Global mean surface temperature (GMST) has increased unequivocally over the twentieth century, associated with ongoing increases in greenhouse gases emitted during anthropogenic activities (Meehl et al. 2005). Decadal periods of warming and cooling were superimposed on this long-term warming trend (England et al. 2014; Kosaka and Xie 2016) and the recent slowdown in GMST warming during the early 2000s has been primarily attributed to internal climate variability, with the interdecadal Pacific oscillation (IPO) making the largest contribution (Meehl et al. 2011; Huber and Knutti 2014; Dai et al. 2015; Trenberth 2015; Zhang 2016). A transition of the IPO toward its negative phase around 2000, characterized by cooling in the eastern tropical Pacific and strengthened Pacific trade winds, has been identified as a key factor in the recent hiatus event (Kosaka and Xie 2013; England et al. 2014; Watanabe et al. 2014). Furthermore, using a global coupled model with the eastern tropical Pacific sea surface temperature (EPSST) anomalies restored to observations, Kosaka and Xie (2016) highlighted the important role of EPSST variations associated with the IPO in modulating GMST. Variations in the Pacific also have been known to affect the Indian Ocean (Nidheesh et al. 2013; Dong et al. 2016; Krishnamurthy and Krishnamurthy 2016; Zhou et al. 2017). The focus of this paper is to examine how atmospheric circulation changes forced by EPSST anomalies influence subsurface ocean heat content (OHC) in the Indian Ocean on decadal time scales.
SST in the Pacific can influence the Indian Ocean through its role in modulating the atmospheric Walker circulation (e.g., Yu and Rienecker 1999; Alexander et al. 2002; Lee and McPhaden 2008; Zheng et al. 2010; Zheng et al. 2013; Yang et al. 2017). Driven by deep convection over the Indo-Pacific warm pool, the large-scale atmospheric circulations over the tropical Indian and the Pacific Oceans are dynamically linked through the equatorial zonal atmospheric overturning circulation. Prevailing easterly winds over the tropical Pacific form the low-level branch of the Pacific Walker cell, while annual-mean low-level westerly winds form that of the Indian Walker cell (Schott et al. 2009). The Walker cell response to variations in SST is associated with zonal shifts in the convective ascending branch (Philander 1985). Han et al. (2017) investigated the covariability of the Pacific and Indian Walker circulations on decadal time scales and found them to highly covary during winter associated with ENSO-like conditions. Dong et al. (2016) suggested that decadal variations in SST in the tropical Indian Ocean are strongly influenced by the IPO, which modulates conditions in the Indian Ocean through changes in surface heat fluxes and sea surface height (SSH) associated with changes in the zonal location and intensity of the Walker cell.
Winds over the Pacific can influence the Indian Ocean by altering the Indonesian Throughflow (ITF; Reason et al. 1996; Murtugudde et al. 1998; Sprintall and Révelard 2014). The ITF carries warm freshwater from the western Pacific into the Indian Ocean, playing a central role in the heat budget of the Indo-Pacific region (Hirst and Godfrey 1993; Godfrey 1996; Vranes et al. 2002). ITF transport tends to increase during La Niña events, responding to reinforced Pacific trade winds (Meyers 1996; England and Huang 2005). A pronounced increasing trend of approximately 1.0 Sv decade−1 (1 Sv ≡ 106 m3 s−1) in ITF transport over the past three decades has been linked to an intensification of the trade winds over the Pacific Ocean (Feng et al. 2011; Liu et al. 2015). Dong and McPhaden (2016) identified an interhemispheric asymmetry in SST trends during the recent global warming hiatus (2000–13), with accelerated warming in the southern Indian Ocean (south of 10°S) and little warming north of 10°S. Enhanced ITF transport induced by stronger Pacific trade winds is the main driver of this interhemispheric asymmetry in SST trends, accounting for most of the observed increase in southern Indian Ocean SST (Dong and McPhaden 2016, 2017b). Enhanced ITF transport also contributes to a pronounced warming in the Indian Ocean subsurface layer (100–300-m depth) during the early 2000s (Nieves et al. 2015; Liu et al. 2016), which in turn accounts for a large portion of the global ocean heat gain in the upper 700 m (Lee et al. 2015).
It is well known that remote signals from the western Pacific Ocean can be transmitted to the Indian Ocean via oceanic baroclinic waves (Clarke and Liu 1994; Chambers et al. 1999; Wijffels and Meyers 2004; Cai et al. 2008). Using observed sea level records and a simple model, Clarke and Liu (1994) suggested that variations in the northwestern Australia sea level and ITF transport are related with wind stress over the tropical Pacific via westward-propagating Rossby waves. Wijffels and Meyers (2004) provided further details on the wave pathways through the Indonesian region, where Rossby waves of Pacific origin become coastally trapped at the intersection of New Guinea and the equator before radiating into the southern Indian Ocean. Cai et al. (2005) identified a signal transmission from the subtropical North Pacific, a pathway that is associated with half of the thermocline variation off northwestern Australia. However, the transmission appears not to be stationary over time: few ENSO signals were transmitted into the Indian Ocean before the 1980s as the subtropical North Pacific pathway was not involved in the transmission, resulting from the weak and narrow meridional extent of ENSO events during that period (Shi et al. 2007).
Upper ocean thermal trends in the Indian Ocean have been discussed widely by previous studies (e.g., Han et al. 2006; Alory et al. 2007; Schwarzkopf and Böning 2011; Ummenhofer et al. 2017). In contrast to rapid surface warming, the subsurface Indian Ocean experienced a prominent cooling trend from the 1960s to the 1990s. This observed subsurface cooling was most pronounced within the tropical thermocline north of 20°S, at depths between 100 and 300 m. Han et al. (2006) attributed this cooling trend to changes in the local wind forcing that acted to shoal the thermocline. Remote contributions from slackened trade winds over the Pacific Ocean have also been identified as a driver for this subsurface cooling trend (Alory et al. 2007; Schwarzkopf and Böning 2011). In contrast, the southeastern Indian Ocean in the upper 400 m has warmed substantially during the early twenty-first century, which has been linked with enhanced ITF transport into the Indian Ocean, driven by strengthened Pacific trade winds associated with the negative phase of the IPO (Li et al. 2017). Despite the wealth of studies addressing cooling or warming trends in the Indian Ocean, it remains unclear how many of the OHC changes in the Indian Ocean arise from internal variability, particularly associated with the IPO, and how many are due to changes in external forcing, particularly that due to anthropogenic greenhouse gases. Previous studies (Dong and Zhou 2014; Han et al. 2014; Dong and McPhaden 2017a; Zhang et al. 2018) have suggested that external forcings are crucial for explaining observed variations of SST in the Indian Ocean. This study will assess the relative contributions of external forcing and internal variability in modulating decadal variations in subsurface OHC in the Indian Ocean based on coupled model experiments.
Using multiple ocean reanalyses, Jin et al. (2018) demonstrated that subsurface OHC decadal variability in the southern tropical Indian Ocean (10°–20°S) is characterized by an east–west dipole structure, with discontinuities in westward-propagating OHC anomalies often occurring around 90°–100°E. They suggested that the mechanisms by which the Pacific Ocean influences the southern tropical Indian Ocean on decadal time scales are distinct between the eastern and western parts of the Indian Ocean. In the eastern Indian Ocean, the IPO influence is largely communicated by the propagation of oceanic Rossby waves through the Indonesian Seas. Variations in the western Indian Ocean emerge from anomalous wind stress curl in the central Indian Ocean, forced remotely by conditions in the Pacific via the atmospheric bridge. These wind stress curl anomalies drive local downwelling or upwelling, which in turn excites westward-propagating Rossby waves that displace the thermocline and alter subsurface OHC in the western Indian Ocean.
Although Jin et al. (2018) suggested that the communication of decadal-scale anomalies in the eastern Indian Ocean with those in the western Pacific Ocean is mainly through Rossby waves, it remains unclear what triggers these oceanic Rossby waves in the western Pacific Ocean. In particular, applying the observed wind forcing to the Pacific Ocean alone in an ocean general circulation model, Schwarzkopf and Böning (2011) reproduced the subsurface cooling trends in the Indian Ocean, suggesting that winds in the Pacific sector play a key role in low-frequency subsurface variability in the Indian Ocean. However, the mechanism by which Pacific winds influence the eastern Indian Ocean was not discussed in their study. Moreover, while Jin et al. (2018) attributed the variations of subsurface OHC in the western Indian Ocean to anomalous wind stress curl over the central Indian Ocean associated with IPO, they did not explore the mechanistic connection between the IPO and this anomalous wind stress curl. Here we extend this previous work, highlighting the role of Pacific winds in the decadal variations in the eastern Indian Ocean and exploring the mechanisms responsible for decadal variations in wind stress curl anomalies in the central Indian Ocean. Our results are based mainly on analysis of the coupled model experiments presented by Kosaka and Xie (2016), which are used to investigate the above questions and to further verify the mechanism proposed by Jin et al. (2018) based on reanalysis products.
The rest of this paper is organized as follows. The coupled model, reanalysis products, and analysis method are described in section 2. The model’s ability to reproduce the timing and characteristics of decadal variations in subsurface OHC in the Indian Ocean is evaluated in section 3. Relative contributions of external forcing and internal variability to these decadal variations in the Indian Ocean are assessed in section 4. Section 5 presents the distinct mechanisms by which the Pacific Ocean influences the eastern Indian Ocean (section 5a) and western Indian Ocean (section 5b) on decadal time scales using the coupled model and then the same mechanisms are confirmed in the free-running coupled model simulations with historical radiative forcings (section 5c). A summary and discussion are provided in section 6.
2. Data and method
We analyze output from two sets of climate model experiments performed with the Geophysical Fluid Dynamics Laboratory Coupled Model version 2.1 (GFDL CM2.1; Delworth et al. 2006) in this study. The historical (HIST) experiment is based on the free-running coupled model subjected to historical radiative forcings prescribed by phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012) for 1861–2005, extended through 2014 using representative concentration pathway 4.5 (RCP4.5). The Pacific Ocean–Global Atmosphere experiment (POGA) is subjected to the same radiative forcing as HIST, but SST anomalies in the eastern tropical Pacific (15°S–15°N and from 180° to the eastern boundary, with a linearly tapering buffer zone of 5° in each direction) are nudged to observations (Kosaka and Xie 2016). Both HIST and POGA are performed with 10-member ensemble runs covering the period from 1861 to 2014, with slightly different initial conditions for each ensemble member. The atmospheric component, the GFDL Atmospheric Model version 2.1 (AM2.1), has 24 vertical levels with horizontal grid spacing of 2° in latitude and 2.5° in longitude. The ocean component, the Modular Ocean Model (MOM) version 4, has 50 vertical layers with a horizontal grid spacing of 1.0° refined to approximately ⅓° in the tropics. All simulations analyzed here are identical to those presented by Kosaka and Xie (2016), but our analyses focus on decadal variations in the Indian Ocean.
The POGA experiment reproduces global climate variability remarkably well and has been used for several purposes, including attributing the recent global warming hiatus (Kosaka and Xie 2013, 2016) and studying northwestern Pacific variability (Kosaka et al. 2013). Besides, Yang et al. (2015) used this POGA experiment to investigate the contribution of ENSO to interannual variability of Indian Ocean dipole events. Compared with observed data, it has been suggested that POGA can successfully capture much of the observed Indian Ocean climatology and interannual variability, including the leading modes of surface variability in the Indian Ocean (the Indian Ocean basin and dipole modes). In particular, the annual mean Indian Ocean SST in POGA is significantly correlated with that of the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) of r = 0.93 (Yang et al. 2015). These results suggest that tropical Pacific SST is the dominant driver of the internal climate variability in the Indo–Pacific sector. Since observed SST anomalies are prescribed only in the eastern tropical Pacific with the rest of the model’s coupled system free to evolve, we regard the POGA simulation as a suitable tool to evaluate the model’s response to observed tropical Pacific variability. In addition, POGA output has the advantage of allowing for direct comparison with observations, as the prescribed SST acts as a pacemaker. Therefore, we primarily present the results from POGA. To confirm the robustness of our results, we also make use of outputs from the HIST simulations subjected to the same external forcing but freely evolving in terms of the internal variability, which are briefly shown toward the end as ancillary examination.
To assess the representation of upper-ocean properties in POGA, the GECCO2 ocean reanalysis is used for comparison. The GECCO2 reanalysis (Köhl 2015) is produced by the German contingent of the Estimating the Circulation and Climate of the Ocean project (ECCO; www.ecco-group.org), based on the Massachusetts Institute of Technology ocean general circulation model (MITgcm). GECCO2 is provided on a 1° regular grid with 50 vertical levels for years between 1948 and 2014. The background atmospheric state is taken from the 6-hourly National Centers for Environmental Prediction–National Center for Atmospheric Research Reanalysis 1 (NCEP–NCAR R1), and data assimilation is conducted using an adjoint method.
To estimate changes due to external forcing, we calculate the time series of global-mean subsurface OHC from the HIST ensemble mean, assuming that internal variability has been averaged out. Then, the regional externally forced component is defined as the regional projection on this externally forced variation in global OHC, which is calculated as the linear regression of local OHC variations onto the time series of global-mean subsurface OHC based on the HIST. After removing the regional externally forced component from each POGA run at every grid point, the ensemble mean of the residual is assumed to represent internally generated variations associated with the tropical Pacific SST anomalies. Similar procedures have been widely used in climate detection and attribution studies based on output from coupled models (e.g., Santer et al. 2013; Dai et al. 2015; Dong and McPhaden 2017b).
Monthly anomalies of each variable are computed by subtracting their corresponding monthly climatology. A low-pass Butterworth filter (Butterworth 1930) has been used to isolate decadal variability by filtering out variability on periods shorter than 8 years. Considering the reduction in the effective degrees of freedom associated with the use of a low-pass filter, statistical significance is determined by using the nonparametric random phase method in the frequency domain to preserve serial correlations (Ebisuzaki 1997). The zonal atmospheric circulation is represented by the zonal streamfunction as defined by Yu et al. (2012):
where is the divergent component of the zonal wind, R is the radius of Earth, p is pressure, and g is the gravitational constant. This equation calculates the eastward mass flux between the top of the atmosphere and a given pressure level p as the vertical integral of the divergent component of zonal wind. The Walker cell is then defined as the zonal streamfunction along the equator averaged between 5°S and 5°N.
indicates the vertical integral from the nominal top of atmosphere (TOA; p = pt) to the surface (p = ps). Here, Qconv = LυP + S represents the column-integrated convective heating, and QR represents heating due to the absorption of longwave and shortwave radiation within the atmospheric column. The latter term is calculated based on net radiative fluxes into the atmosphere at the surface and at TOA. The product LυP represents the vertically integrated latent heating, with Lυ being the latent enthalpy of vaporization and P the surface precipitation, and S the (upward) sensible enthalpy flux at the surface.
3. Model evaluation
The output from the 10-member POGA coupled-model ensemble (in which tropical Pacific SST anomalies are nudged toward observations) is briefly assessed against the GECCO2 ocean reanalysis. To explore the role of tropical Pacific SST in modulating OHC in the Indian Ocean, Fig. 1 shows subsurface OHC anomalies regressed onto EPSST for GECCO2 and each of the 10 POGA ensemble members, along with the POGA ensemble mean. EPSST is defined as the area-mean SST anomalies over the eastern tropical Pacific (15°S–15°N, 180°–80°W), which is the same region as the SST anomalies restoring applied in POGA, and it has been 8-yr low-pass filtered to isolate the decadal variability. To avoid artifacts that arise from the low-pass filtering, the first and last four years are excluded in all our analyses. Subsurface OHC anomalies are calculated by integrating temperature anomalies over 50–300-m depth relative to the mean annual cycle during 1900–2014. In GECCO2, the spatial regression pattern in the Indian Ocean is asymmetric about the equator with the largest loading in the southern tropical Indian Ocean. A significant zonal dipole structure occurs in the southern Indian Ocean, especially in the latitude band of 10°–20°S (black box in Fig. 1), characterized by positive anomalies west of 100°E and negative anomalies to the east as noted by Jin et al. (2018). Connected by the Indonesian Archipelago, negative OHC anomalies in the western Pacific Ocean extend along the northwestern Australian coast and radiate into the eastern Indian Ocean in GECCO2. Compared with GECCO2, the POGA ensemble members capture the dipole structure in the southern Indian Ocean reasonably well, with a spatial pattern correlation of 0.69 (P < 0.05) between GECCO2 and the POGA ensemble mean. In particular, the largest loading is also located between 10° and 20°S in the POGA simulations. Moreover, an analogous connection between the western Pacific and the eastern Indian Ocean occurs via the Indonesian Seas among POGA members. Results based on the 10 HIST ensemble members are similar (not shown). As the maximum connection to the Pacific Ocean occurs in the latitude band of 10°–20°S in the Indian Ocean, we will focus on OHC variations over this band in the following analyses.
Figure 2 depicts longitude–time Hovmöller diagrams of subsurface OHC anomalies across the Indian Ocean from GECCO2, the individual POGA ensemble members, and the POGA ensemble mean. Subsurface OHC anomalies are averaged meridionally along the band of 10°–20°S where the largest variations occur. The Hovmöller plot based on GECCO2 suggests that decadal-scale OHC anomalies originating at the eastern boundary of the Indian Ocean propagate westward via oceanic Rossby waves to impact the southern Indian Ocean. However, the westward propagations were often blocked in the central Indian Ocean, especially during the 1980s. The proximate location of the discontinuity occurs near 100°E, with westward-propagating OHC anomalies in the western part of the Indian Ocean often appearing to originate from the central Indian Ocean rather than from the eastern boundary. The existence and possible origin of these discontinuities in ocean reanalyses were reported by Jin et al. (2018). Since POGA simulations allow for direct year-by-year comparisons against observational data due to the prescribed EPSST pacemaker, analogous Hovmöller diagrams based on POGA ensemble members are also shown in Fig. 2. The overall resemblance of OHC anomalies between POGA and GECCO2 implies that POGA reproduces decadal variations of subsurface OHC anomalies in the Indian Ocean, although the strong warming after 2005 in the western Indian Ocean is absent in all runs. The absence of this warming signal in POGA runs may arise from model bias. An apparent discontinuity also takes place at the same location (near 100°E) in each POGA run and even more prominent in their ensemble mean (Fig. 2l), consistent with that in GECCO2. Both these discontinuities and the dipole structure suggest that the mechanisms connecting decadal variability in the Pacific Ocean to Indian Ocean subsurface OHC are distinct between the eastern and western parts of the Indian Ocean. As such, we will separate the southern tropical Indian Ocean into eastern and western parts in the following analyses (with the dividing line at 100°E), and investigate how subsurface OHC in each part responds to decadal variations in the Pacific Ocean.
Time series of low-pass filtered subsurface OHC anomalies averaged over the eastern Indian Ocean (10°–20°S, 100°–120°E) based on the 10-member HIST and POGA ensembles are compared in Fig. 3. OHC anomalies in the HIST simulations show considerable diversity across ensemble members. By contrast, more consistent variations are exhibited among POGA members with specified EPSST, implying a critical role of the tropical Pacific in determining decadal variations in subsurface OHC in the eastern Indian Ocean. In addition, POGA ensembles do an excellent job in simulating the variation of OHC anomalies in the eastern Indian Ocean, as evidenced by their close correspondence to GECCO2 (solid black line). The correlation between the POGA ensemble mean and GECCO2 is 0.66 (P < 0.05), suggesting that around 43% of the decadal variance in the Indian Ocean based on GECCO2 is reproduced by EPSST and external forcing based on the POGA ensemble mean. In particular, the observed unprecedented warming in the eastern Indian Ocean after 2003 (Li et al. 2017) is also reproduced in the POGA members. Similarly, POGA also reproduces a realistic evolution of subsurface OHC in the western Indian Ocean (not shown). Given the close resemblance between POGA and GECCO2, POGA can be used to attribute and explain the mechanistic connections between the tropical Pacific Ocean and the Indian Ocean on decadal time scales.
4. Internal variability versus externally forced component
In this section, we examine the relative contributions of external forcing and internal variability to decadal variations in the Indian Ocean. The externally forced component and the internal variability are separated as explained in section 2. Note that the ensemble mean of the 10 POGA members can be interpreted as an estimate of the model’s forced response to external forcing plus observed variations in tropical Pacific SST. Assuming they are linearly additive, the internal variability associated with EPSST is obtained by subtracting the externally forced component from the POGA simulations. The competing effects of internal variability associated with EPSST and external forcing are then compared in Fig. 4 for the eastern and western parts of the southern tropical Indian Ocean.
In the eastern Indian Ocean, the overall trend of the externally forced variation is an upward curve, indicating that the general effect of external forcing is to increase OHC. The close resemblance of OHC anomalies between the total variation and the internally generated variation suggests that the internal variability associated with EPSST accounts for most of the decadal variability in the eastern Indian Ocean. To quantify the impact of EPSST on the decadal variability of subsurface OHC in the Indian Ocean, correlations between the OHC variation in each POGA member and the ensemble mean OHC in POGA (after removing the externally forced signal) are calculated. The averaged correlation across all POGA members is 0.76 ± 0.04 (range: 0.69–0.81) in the eastern Indian Ocean and 0.66 ± 0.06 (range: 0.58–0.78) in the western Indian Ocean, suggesting that around 58% (44%) of the decadal variation of OHC in the eastern (western) Indian Ocean is associated with fluctuations in tropical Pacific SST. The rapid warming in the eastern Indian Ocean since 2003 suggested by Li et al. (2017) was largely derived from internally generated variability, with external forcing making a secondary contribution to accelerate the warming trend (around 15%). Additionally, we note that the specification of EPSST in the POGA simulations results in much less spread, as shown in Fig. 3b. This reduction in the spread attests to the importance of EPSST-related internal variability in determining decadal variability in the Indian Ocean. The dominance of internal variability forced by EPSST is also exhibited in the western Indian Ocean (Fig. 4b). It is noteworthy that OHC anomalies in the western Indian Ocean have a variation that is out of phase with those in the eastern Indian Ocean with the simultaneous correlation being −0.80 (P < 0.05), consistent with the east–west dipole structure in the southern Indian Ocean.
As discussed in the introduction, many attempts have been made to explain OHC trends in the Indian Ocean; however, the extent to which OHC trends derive from external forcing as opposed to internal variability remains an open question. As the 30-yr cooling trend from the 1960s to 1990s has been discussed widely by many studies, we examine the relative importance of internal variability associated with EPSST and external forcing for generating 30-yr moving trends in subsurface OHC in the eastern and western parts of the Indian Ocean (Fig. 5). Our analysis indicates that external forcing cannot explain the observed cooling trend between the 1960s and 1990s. The contributions of external forcing acted to increase subsurface OHC in both the eastern and western Indian Ocean, with a stronger signal in the eastern part. By contrast, internally generated variations were responsible for the cooling trends in most of the POGA ensemble members over that period, with these trends appearing earlier in the western Indian Ocean. Indications of an internally generated cooling trend appeared first in the western Indian Ocean (1960–89; P < 0.1), before reverting to weak warming or neutral trends over 1965–94. Cooling trends appeared later in the eastern Indian Ocean, with weak cooling over 1965–94 becoming pronounced cooling over 1970–99, statistically significant at the 90% confidence level based on the ensemble spread. This cooling in the subsurface eastern Indian Ocean shifted back to warming for 1975–2004, with strong subsurface warming emerging during 1980–2009.
The 30-yr trends of subsurface OHC based on GECCO2 are shown by black dots in Fig. 5. Trends based on GECCO2 are within the range of internal variability-based trends from the POGA ensemble members, except for trends in the eastern Indian Ocean from the 1960s to 1990s. During this period, when trends in GECCO2 indicate cooling, trends based on the POGA ensemble mean are neutral, with some individual members even indicating warming. This difference may be due to the fact that a part of the internal variability is independent of the prescribed EPSST, but may also be due to model biases. It will therefore be necessary to systematically investigate the roles of internal variability and external forcing in the future using larger ensembles and different models to further resolve this issue.
5. Distinct mechanisms in the eastern and western Indian Ocean
The east–west dipole structure and frequent occurrence of discontinuities in westward-propagating signals in the central Indian Ocean suggest that the western and eastern parts of the Indian Ocean respond to the tropical Pacific in distinct ways on decadal time scales (Jin et al. 2018). Here, using POGA simulations, we diagnose the mechanisms by which the tropical Pacific Ocean impacts decadal variations in the eastern and western parts of the Indian Ocean, respectively. We then confirm that these responses arise from the same mechanisms in free-running HIST simulations. Since we focus on the relationship between the Indian Ocean and Pacific Ocean, the externally forced component of variations in OHC has been removed from each of the POGA and HIST runs before analysis.
a. Eastern Indian Ocean
As shown by the regression map of subsurface OHC anomalies onto EPSST (Fig. 1), negative OHC anomalies in the western Pacific are linked to negative OHC anomalies in the eastern Indian Ocean on decadal time scales via the Indonesian archipelago. Inspired by Schwarzkopf and Böning (2011), who suggested that Pacific wind is important to Indian Ocean decadal variation, we explore the role of Pacific wind by regressing the surface wind onto the time series of area-mean OHC in the eastern Indian Ocean using POGA simulations (Fig. 6). Positive OHC anomalies in the eastern Indian Ocean correspond to easterly wind anomalies over the equatorial western Pacific Ocean, indicating an intensification of the prevailing easterly trade winds. Enhanced easterly trade winds near the equator push more tropical water westward, piling up warm water in the western Pacific Ocean and deepening the thermocline along the western boundary. Note that easterly winds anomalies within the red box decrease northward, indicating the presence of strong negative wind stress curl anomalies over the western Pacific (0°–15°N, 125°–150°E) near the entrance to the Indonesian Sea. Zonal wind and wind stress curl are interdependent, with a simultaneous correlation of 0.57 (P < 0.05). Both factors contribute to a deepening of the thermocline. Given the highly coupled nature between the tropical wind and EPSST (simultaneous correlation is 0.58) in the Pacific Ocean through the Bjerknes feedback (Bjerknes 1969), we suggest that the tropical Pacific modulates the eastern Indian Ocean mainly through the western Pacific surface wind (zonal wind and wind stress curl) altering the depth of the western Pacific thermocline and thereby exciting oceanic Rossby waves that radiate into the southeastern Indian Ocean (Clarke and Liu 1994; Wijffels and Meyers 2004).
To further examine the association between the Pacific wind and eastern Indian Ocean OHC on decadal time scales, the 8-yr low-pass filtered time series of western Pacific zonal wind and wind stress curl averaged over the region of 0°–15°N, 125°–150°E (red box in Fig. 6) are shown in Fig. 7a, along with the 8-yr low-pass filtered time series of OHC anomalies in the eastern Indian Ocean (10°–20°S, 100°–120°E; black box in Fig. 6). All time series are based on the POGA ensemble mean and have been normalized to facilitate comparison. It is evident that temporal variation of OHC anomalies in the eastern Indian Ocean is anticorrelated with temporal variations of western Pacific zonal wind and wind stress curl. The lead–lag correlation between the OHC and zonal wind time series peaks at −0.75 (P < 0.05) when western Pacific wind leads OHC by approximately 5 months (Wijffels and Meyers 2004), implying the role of western Pacific surface wind in modulating decadal variations in subsurface OHC in the eastern Indian Ocean. Results are similar for lead–lag correlations against the western Pacific wind stress curl, which peak around −0.6 (P < 0.05) with 9-month leads (not shown).
Since the ITF delivers large amounts of heat into the Indian Ocean and has a large impact on the eastern Indian Ocean heat budget, the role of the ITF is examined. The ITF volume transport is integrated over the upper 700 m across the 113.5°E longitudinal transect (8°–25°S; see also England and Huang 2005; Dong and McPhaden 2016). It is estimated that the mean ITF transport in POGA is around 12 Sv with a standard deviation of 3.7 Sv, consistent with the estimate of 11.8 ± 2.9 Sv derived from other ocean reanalyses (Dong and McPhaden 2016), but slightly larger than the estimate of 8.9 ± 1.7 Sv based on expendable bathythermograph measurements (Wijffels et al. 2008; Liu et al. 2015). Westward transport by the ITF is defined as positive, so that positive ITF anomalies correspond to more warm water transported into the Indian Ocean. Variations in the ITF (yellow dashed line in Fig. 7a) are significantly correlated with OHC in the eastern Indian Ocean, peaking at 0.75 (P < 0.05) when ITF leads OHC by approximately 12 months. Responding to stronger trade winds, elevated sea level in the western Pacific Ocean tends to increase ITF transport, which in turn warms the eastern Indian Ocean several months later. The ITF transports are calculated at multiple zonal and meridional sections to confirm that this finding is not sensitive to the exact location of the section (not shown). Although western Pacific winds are important contributor to variations in ITF transport, other factors, such as the Indian Ocean dipole (Sprintall et al. 2009; Sprintall and Révelard 2014) and Asian monsoon (Gordon et al. 2003; Susanto et al. 2007) also exert large influences on the ITF transport.
We also note the existence of a significant correlation between zonal wind over the central tropical Indian Ocean (0°–5°S, 75°–95°E; blue box in Fig. 6) and OHC in the eastern Indian Ocean. Decadal variations of zonal wind in this region are highly correlated with decadal variations in OHC in the eastern Indian Ocean (Fig. 8a), suggesting that zonal wind over the central tropical Indian Ocean may also influence OHC in the eastern Indian Ocean via eastward-propagating equatorial Kelvin waves followed by southward-propagating coastal Kelvin waves. The lead–lag correlation supports this notion (Fig. 8b), with the maximum correlation occurring when zonal wind over the central tropical Indian Ocean leads OHC in the eastern Indian Ocean by 1 month. Furthermore, the zonal winds over the western tropical Pacific Ocean and the central tropical Indian Ocean are highly correlated through the modulation of the Walker cell as discussed in the next subsection.
b. Western Indian Ocean
Using reanalysis data, Jin et al. (2018) attributed variations of OHC in the western Indian Ocean to anomalous wind stress curl over the central Indian Ocean associated with the IPO. These wind stress curl anomalies induce anomalous Ekman pumping, which in turn generates westward-propagating Rossby waves that alter OHC in the western Indian Ocean. To test this mechanism in the coupled model experiments, the time series of area-mean wind stress curl over the central Indian Ocean (10°–20°S, 90°–100°E) and OHC in the western Indian Ocean (10°–20°S, 50°–90°E) from the POGA ensemble mean are shown in Fig. 9a. The wind stress curl is averaged over the region where the discontinuity emerges in westward-propagating subsurface OHC anomalies, as shown in Fig. 2. Obviously, OHC anomalies in the western Indian Ocean invariably follow wind stress curl anomalies in the central Indian Ocean, so we further examine the time lag relation between wind stress curl and OHC by using a lead–lag correlation analysis (Fig. 9b). The maximum correlation coefficient occurs when the wind stress curl leads the OHC by approximately 9 months. The results from the POGA coupled model experiment therefore confirm that wind stress curl anomalies over the central Indian Ocean associated with variations in EPSST drive decadal variations in subsurface OHC in the western Indian Ocean.
How could the tropical Pacific Ocean SST induce wind stress curl anomalies over the central Indian Ocean? As the tropical Pacific and Indian Oceans are coupled through the Walker cell, changes in this zonal atmospheric circulation would modulate the condition over the Indian Ocean. To investigate the response of the Walker cell to the tropical Pacific Ocean SST, Fig. 10 depicts the regression of zonal mass streamfunction anomalies in the atmosphere onto EPSST (shading), along with the climatological mean of the streamfunction (contours) based on the POGA ensemble mean. The Walker cell is defined as the zonal streamfunction along the equator, with positive values denoting clockwise circulation. In normal years, ascending motion prevails in the main convection region, including the western Pacific and eastern Indian Ocean. When positive SSTA occurs in the eastern tropical Pacific Ocean (corresponding to the positive phase of the IPO), the Walker cell shifts eastward (Deser and Wallace 1990). In response to the eastward shift in the ascending branch of the Walker cell, large-scale descending anomalies appear over the eastern Indian Ocean, suppressing the convection in the Maritime Continent (Lau and Nath 2003; Wang et al. 2003). The meridionally averaged (5°S–5°N) atmospheric diabatic heating anomalies regressed onto EPSST are also shown in Fig. 10. Corresponding to the suppression of convection, diabatic cooling appears over the Maritime Continent. These descending anomalies extend through most of the troposphere, with the most prominent anomalies located around 110°E. The surface circulation response to this deep layer of anomalous diabatic cooling manifests as the inverse of the classical Matsuno–Gill pattern (Matsuno 1966; Gill 1980), which includes an off-equatorial pair of anomalous surface anticyclones to the west of the anomalous descending branch centered around 100°E. The prominent anticyclonic anomaly over the southern Indian Ocean (and its weaker Northern Hemisphere counterpart; Fig. 11) represents the atmospheric Rossby wave response to anomalous diabatic cooling near 110°E, resulting from an eastward shift in the ascending branch of the Indo-Pacific Walker cell. This anomalous anticyclone gives rise to positive wind stress curl anomalies over the central Indian Ocean that trigger westward-propagating oceanic Rossby waves to impact the western Indian Ocean. The wind stress curl anomalies, and thus the connections between EPSST and OHC in the western Indian Ocean, arise primarily from EPSST-related changes in atmospheric heating associated with the Walker cell.
c. HIST simulations
The above analysis based on the POGA output identifies key oceanic and atmospheric mechanisms by which the influence of the tropical Pacific Ocean is transmitted to the eastern and western parts of the Indian Ocean on decadal time scales. These mechanisms should also operate in individual realizations of the HIST experiment, which allow for the free evolution of all modes of internal variability, including EPSST, although their time evolution is not synchronized as in the POGA pacemaker simulations. We therefore examine output from the HIST ensemble to further confirm the mechanisms proposed above and results are briefly shown here.
Analogous lead–lag correlations to Figs. 7b, 8b, and 9b based on HIST simulations are provided in Fig. 12. Recalling the large spread of subsurface OHC in the eastern Indian Ocean across HIST ensemble members (Fig. 3), a similar association between the western Pacific zonal wind and eastern Indian Ocean subsurface OHC highlights the role of western Pacific wind in modulating subsurface OHC in the eastern Indian Ocean on decadal time scales. The highest correlation with the eastern Indian Ocean OHC is around −0.7 when Pacific zonal winds lead by approximately 8 months (Fig. 12a). As in POGA, zonal wind over the central tropical Indian Ocean also exhibits significant correlation with OHC in the eastern Indian Ocean, with the coefficient peaking at 0.7 when zonal winds lead by approximately 3 months (Fig. 12b). As for the western Indian Ocean, similar results to POGA are obtained between wind stress curl over the central Indian Ocean and OHC in the western Indian Ocean (Fig. 12c), consistent with the notion that wind stress curl associated with the Pacific Ocean accounts for much of the decadal variability of OHC in the western Indian Ocean. The similar results from HIST confirm the robustness of findings based on POGA.
6. Summary and discussion
We have used output from coupled model experiments to examine the effects of internal variability associated with tropical Pacific SST relative to external forcing in generating decadal variations of ocean heat content (OHC) in the subsurface Indian Ocean, along with the mechanisms by which conditions in the Pacific Ocean influence the subsurface OHC in the southern tropical Indian Ocean on decadal time scales. This analysis is based on two sets of 10-member coupled model ensemble experiments performed with the GFDL CM2.1. In the historical (HIST) experiment, internal climate variability evolves freely, whereas in the Pacific Ocean–Global Atmosphere experiment (POGA) SST anomalies in the eastern Pacific Ocean are nudged to observations, leaving the rest of the climate system free to evolve.
We first assess the realism of POGA simulations by comparing the POGA output with ocean reanalysis data, namely GECCO2. The overall resemblance between POGA and GECCO2 confirms that the POGA ensemble members capture the historical evolution of conditions in the Indian Ocean on decadal time scales, and suggests that these simulations can be used to investigate decadal variations in subsurface OHC in the southern Indian Ocean. Comparing with HIST ensemble members, the decadal variations of subsurface OHC anomalies in the Indian Ocean are better reproduced after constraining the evolution of EPSST to closely match observed values in POGA, indicating that the tropical Pacific Ocean plays a key role in modulating conditions in the Indian Ocean on decadal time scales. GECCO2 and POGA both show the emergence of discontinuities in westward-propagating OHC anomalies in the central Indian Ocean around 100°E. Like GECCO2, prominent east–west dipole structures are also displayed in the regression maps of subsurface OHC anomalies onto EPSST among POGA ensemble members, consistent with the notion of distinct responses in the eastern and western parts of the Indian Ocean to the Pacific changes.
We quantified the relative roles of internal variability and external forcing to the variations of subsurface OHC in the tropical south Indian Ocean on the decadal and longer time scales based on coupled model simulations. Assuming that external forcing and internal variability are linearly superimposed, we separated these two components and found that the decadal variations of subsurface OHC in the Indian Ocean are dominated by internally generated variations associated with EPSST. In particular, the widely discussed subsurface cooling trend between the 1960s and 1990s is primarily a result of internally generated variability driven by the tropical Pacific changes, rather than external forcing. In terms of the rapid warming since 2003, internal variability is the main contributor to this warming with the externally forced component acting to augment the warming rate by around 15%.
Using the POGA simulations, we then explored how the Pacific Ocean impacts the eastern and western part of Indian Ocean on decadal time scales. In conjunction with the positive phase of the IPO (including the warmer eastern tropical Pacific SST), the surface winds change over both the western Pacific and central Indian Ocean owing to an eastward shift of the Walker cell. The weakened easterly winds in the western tropical Pacific drive negative sea level height anomalies and thermocline shoaling along the western boundary of the tropical Pacific, which in turn result in decreased OHC in the eastern Indian Ocean via oceanic Rossby wave propagation and reduced Indonesian Throughflow (ITF) transport. Wind changes in the central Indian Ocean play two roles. First, weakened westerly winds near the equator generate upwelling Kelvin waves that reduce subsurface OHC in the eastern Indian Ocean, amplifying the OHC anomalies that arise from the western Pacific wind anomalies. At the same time, the associated positive wind stress curl anomalies over the southern central tropical Indian Ocean generate downwelling Rossby waves, which propagate westward and increase OHC in the western Indian Ocean. The coupled model simulations therefore indicate that decadal changes in OHC in the eastern Indian Ocean are impacted by the tropical Pacific Ocean via oceanic processes through the Indonesian Seas, while opposite-signed OHC anomalies in the western Indian Ocean are linked to the tropical Pacific via atmospheric processes. In accordance with POGA, individual members of the HIST ensemble, the fully coupled free-running model with only external forcing prescribed, exhibit similar responses in the eastern and western parts of the Indian Ocean to decadal changes in tropical Pacific SST, supporting our findings based on POGA.
In this study, we have examined the relative contributions of external forcing and internal variability to decadal variations in subsurface OHC in the southern tropical Indian Ocean by using 10-member ensembles of coupled model simulations. Our findings may be impacted by the limited ensemble size. Indeed, we have examined the extent to which the HIST ensemble mean at each grid point can represent the externally forced component. However, when this method is used, the externally forced component does not exhibit an overall warming trend, which is inconsistent with our current knowledge. Thus, larger ensembles with the same model or other models will be needed to more robustly determine the relative roles of external forcing and internal variability in generating decadal variations in the Indian Ocean (e.g., Frankignoul et al. 2017).
In addition, while we have highlighted the role of conditions in the tropical Pacific Ocean in driving decadal variations in Indian Ocean subsurface OHC, the extent to which signals in the North Pacific and South Pacific influence decadal variations in the Indian Ocean remains to be clarified. Some studies (e.g., Krishnamurthy and Krishnamurthy 2016; Ummenhofer et al. 2017) have linked decadal variations in the Indian Ocean to the Pacific decadal oscillation (PDO), while others (e.g., Vargas-Hernandez et al. 2014) have suggested a potential role of the southern Pacific Ocean in modulating conditions in the Indian Ocean on decadal time scales. Decadal variability in the Pacific is characterized by a tripole pattern (Henley et al. 2015). The relative roles of these three poles of SST variability in driving decadal variability in the Indian Ocean remains an open question that needs to be quantified in future work.
This research was supported by the Independent Research and Development Program at WHOI to CCU, an NSF OCE PO grant (NSF OCE-1242989) to Young-Oh Kwon, NOAA CP CVP grants (NA15OAR4310176 and NA17OAR4310255) to Hyodae Seo, and a research grant from the Ministry of Science and Technology of the People’s Republic of China to Tsinghua University (2017YFA0603902).