Projected future changes in the contribution of Indo-Pacific sea surface height variability to the Indonesian throughflow

The Indonesian throughflow (ITF) transports a significant amount of warm freshwater from the Pacific to the Indian Ocean, making it critical to the global climate system. This study examines decadal ITF variations using ocean reanalysis data as well as climate model simulations from the Coupled Model Inter-comparison Project Phase 5 (CMIP5). While the observed annual cycle of ITF transport is known to be correlated with the annual cycle of sea surface height (SSH) difference between the Pacific and Indian Oceans, ocean reanalysis data (1959–2015) show that the Pacific Ocean SSH variability controls more than 85% of ITF variation on decadal timescales. In contrast, the Indian Ocean SSH variability contributes less than 15%. While those observed contributions are mostly reproduced in the CMIP5 historical simulations, an analysis of future climate projections shows a 25–30% increase in the Indian Ocean SSH variability to decadal ITF variations and a corresponding decrease in the Pacific contribution. These projected changes in the Indian Ocean SSH variability are associated with a 23% increase in the amplitudes of negative zonal wind stress anomalies over the equatorial Indian Ocean, along with a 12º eastward shift in the center of action in these anomalies. This combined effect of the increased amplitude and eastward shift in the zonal wind stress increases the SSHA variance over the Indian Ocean, increasing its contribution to the ITF variation. The decadal ITF changes discussed in this study will be crucial in understanding the future global climate variability, strongly coupled to Indo-Pacific interactions.


Introduction
The Indonesian throughflow (ITF) transports a huge volume of seawater from the warmer and less saline western tropical Pacific to the colder and more saline southern tropical Indian Ocean via the straits of the Indonesian Archipelago. The observed net mass transport through the straits of the Indonesian Archipelago is approximately 15 Sverdrup (Sv, 1 Sv = 10 6 m 3 s −1 ), with seasonal variations ranging from 0.5 to 2.5 Sv (Gordon et al. 2008;Sprintall et al. 2009). Sizable interannual variability in the total ITF, with a typical amplitude of 3-3.7 Sv, has also been reported based on hydrographic observations and numerical modeling studies (Susanto et al. 2012;Gordon et al. 2008;Meyers 1996;Sprintall and Revelard 2014). These ITF fluctuations increase basin-scale temperature and salinity variations in the ocean surface and subsurface layers in the Indian Ocean. They modulate climate variability over the countries along the Indian Ocean coast via the ocean-atmosphere interaction (Lee and McPhaden 2008;Murtugudde et al. 1998;Potemra and Schneider 2007;Schneider 1998;Song et al. 2007;Tokinaga et al. 2012). Lee et al. (2015) suggested that decadal ITF variability played an essential role in redistributing the excess heat from the tropical Pacific into the Indian Ocean during the early 2000s. This may be linked to the significant ocean heat uptake over the tropical Pacific during the same period which slowed global-mean surface temperature warming (Kosaka et al. 2013;Watanabe et al. 2013).
Wyrtiki (1987) proposed that the sea surface height (SSH) difference between the northwestern tropical Pacific (NWP) and the southeastern Indian Ocean (SEI) is an excellent measure to represent the ITF transport. Large-scale surface winds mostly maintain this SSH difference over the tropical Pacific and Indian Oceans. Climatological surface easterlies over the tropical Pacific associated with the Hadley and Walker circulations act to pile up surface waters, increasing SSH in the NWP region than the relatively lower SSH in the SEI region. Since the SSH in the SEI region is controlled by the climatological Indian Monsoon (Han et al. 2014), seasonal reversal of surface winds can contribute to seasonal ITF variations. Besides the climatological viewpoint, modulating these surface winds on interannual and decadal time scales can induce ITF fluctuations by changing the SSH difference between the NWP and SEI regions.
While basin-scale surface wind anomalies with the El Niño-Southern Oscillation (ENSO) and the Indian Ocean Dipole (IOD) events can induce SSH anomalies (SSHAs) on either side of the Indonesian Archipelago, respectively (Clark and Liu 1994;Gordon et al. 1999;Meyers 1996;Sprintall et al. 2009), frequent co-occurrences of the El Niño-positive IOD and La Niña-negative IOD events (Roxy et al. 2011) scarcely influence the SSH difference between the NWP and SEI regions (Murtugudde et al. 1998;Potemra and Schneider 2007;Sprintall and Revelard 2014). Along with this subtle effect of the interannual wind variability on the SSHAs, Lee and McPhaden (2008) mentioned the negative correlation of the surface winds over the two basins also occurred in the recent interdecadal change between 1993-2000 and 2000-2006, indicating little effect on the SSH difference across the Indonesian Archipelago straits during that interdecadal change. In contrast to the argument that the interannual and decadal time scales are similar, Nidheesh et al. (2013) demonstrated that decadal surface wind variations are relatively independent over the two basins for a longer period .
The inconsistent results among previous studies may have arisen from the following reasons, such as the varying relationship between the Pacific and Indian Oceans in different decades (Han et al. 2014), limited duration of available satellite observations (Han et al. 2014), and/or different representations of long-term variability over the Indian Ocean in atmospheric reanalysis datasets (Nidheesh et al. 2013).
They also demonstrated that dominant spatial patterns of the simulated SSHA in an ocean general circulation model are different between interannual and decadal time scales (see Fig. 6 in Nidheesh et al. 2013). The leading SSHA mode on the interannual time scale shows comparable loadings with the same polarity across the straits of the Indonesian Archipelago, resulting in an unchanged SSH difference between the NWP and SEI regions. In contrast, that on the decadal time scale depicts significant loadings in the NWP region, with subtle loadings in the SEI region. Hence, the SSH difference across the Indonesian Archipelago straits associated with this decadal leading mode can be the major contributor to decadal ITF variations, implying that the NWP and SEI regions contribute differently to decadal ITF variations.
It is well known that the Pacific Decadal Oscillation (PDO, Mantua et al. 1997) and/or Interdecadal Pacific Oscillation (IPO, Power et al. 1999 are the most dominant decadal variability all over the North and South Pacific via the Bjerknes feedback (Bjerknes 1969) over the tropical Pacific. The PDO/IPO defined in the SST field accompanies the decadal modulation of the easterly trade winds. Abe et al. (2016) showed that the surface wind variability over the tropical Pacific during 1993-2014 induced the decadal SSHA variability in the off-equatorial Pacific. A recent study by Li et al. (2018), based on the observed and reanalysis of surface winds only over the North Pacific, also showed the substantial contribution of anomalies in surface winds over the tropical Pacific to decadal ITF variations. In addition, the PDO/IPO can involve a trans-basin atmospheric varialibity via the modulation of the Walker Circulation (Dong et al. 2016). These decadal variations of the surface winds are associated with the possible changes in the SSHA in the Pacific and Indian Ocean (Deepa et al. 2019).
Due to global warming, the decadal ITF variation, and relative contributions of decadal SSH variations over the NWP and SEI regions to the ITF, will change, because a change in the fundamental basement in the climate system, including the surface wind field, affects the behaviors of perturbations around the mean state. Future climate projections conducted in Coupled Model Intercomparison Project Phase 5 (CMIP5) consistently show increased sea surface temperature (SST) variances of the SEI region due to the shallowing of the thermocline in the eastern equatorial Indian Ocean (Cai et al. 2013). Hence, this study aimed to understand the relative contribution of decadal SSH variability over the Pacific and Indian Oceans to decadal ITF variations and assess the dynamics behind a change in the contribution in a future climate by analyzing the SSH and surface wind fields. Samanta et al. (2021) investigated the vetical structure of the volume transport in the Indonesian stratis using high resolution regional ocean models and global climate models, and then pointed out the diversity in these vertical structures in the Indonesian straits, which is indicative of the role of widely varying bathymetry. Indeed, CMIP5 models cannot resolve the Indonesian Archipelago limiting the direct calculation of ITF transport based on the simulated velocity in the ocean. To overcome this inability in the velocity field, we employed a new method to estimate the anomalous SSH difference using only observed and reanalysis atmospheric wind forcings of the Pacific and Indian Oceans to obtain proxy ITF variations. Apart from overcoming the inability of CMIP5 models to simulate ITF, this is new method of estimating ITF variability which is independent of complex ocean processes and yet proved to be very effective in capturing decadal ITF variability.
Section 2 describes the data and methodology used in this study. Section 3 presents estimation techniques used to develop proxy ITF time-series and the relative contribution from the Pacific and Indian Oceans. Section 4 explains the challenges in using CMIP5 data and the diagnosis we performed for selecting CMIP5 models for analysis. Section 5 discusses future changes in the contribution. A discussion in Sect. 6 demonstrates how these changes are linked to largescale wind forcings. We conclude the study in Sect. 7.

Data and methodology
Observational and numerical studies on the ITF show that Makassar and Lifamatola are the inflow straits of the ITF from the NWP region, whereas Lombok, Ombai, and Timor are the outflow straits of the ITF into the SEI region (Gordon et al. 2008;Sprintall et al. 2004Sprintall et al. , 2009. While they showed that the volume transport in the Makassar strait accounts for 80% of the total volume transport of the ITF, observations conducted in the Makassar strait are not long enough to examine decadal ITF variations. Hence, monthly mean velocity is integrated from the surface to the depth above the bottom on a line from Australia to Indonesia between 9º S and 20º S along 115º E, based on the Ocean Reanalysis System 4 (ORAS4) dataset (Balmaseda et al. 2013), which depicts the long-term variability in ITF volume transport from 1959 to 2015 (Fig. 1). While the standard deviation of total ITF variation is 1.70 Sv, the one on the decadal variation is 0.50 Sv, which is around 29% of the total ITF variation. This means that the decadal ITF variation is significant enough to obtain measurable anomalies from the vertically integrating velocity of the ORAS4 reanalysis data.
The monthly mean SSH of the same ORAS4 dataset is used to obtain the SSH difference between the NWP and SEI region for 1959-2015. We also use monthly mean SSH observed by satellite altimeters provided by the Archiving, Validation, and Interpretation of Satellite Oceanographic data (AVISO) from 1992 to 2015 to validate the monthly SSHA in the ORAS4 reanalysis dataset. While the SSH in the ORAS4 and AVISO are not independent, because the observed data including satellite measurements are assimilated in the reanalysis, we made an inter-comparison to validate usefulness of the monthly SSHA in the ORAS4 reanalysis.
We will attempt to estimate departures from the mean SSH difference (i.e., SSHA difference) between the NWP and SEI regions by applying only surface wind data to a 1.5layer reduced gravity model which has successfully reproduced Rossby wave features in many previous studies (Abe et al. 2014;Meyers et al. 1979;Kessler 1990;Holbrook and Bindoff 1999). Details in the model settings in the present study are identical as described in Abe et al. (2014). This model experiment was conducted over the domains of the North Pacific (125ºE-90ºW, 6-16ºN) and the Indian Ocean (85-115º E, 6-16º S), respectively, with 1º × 1º resolution. As coastal boundary conditions, the SSH of the ORAS4 along the eastern border of each domain is used. As the forcing product, we used the monthly zonal and meridional surface wind stress in the previous version of the product, ORA-S3 (Balmaseda et al. 2008(Balmaseda et al. ), for 1959(Balmaseda et al. -2005. While the forcing product of the ORA-S3 was ERA40 (1959ERA40 ( -2002 and ECMWF operation data (2002)(2003)(2004)(2005), the one of the ORAS4 was from three different products of ERA40 , ERA-interim (1989-2009), and ECMWF operational data (2010-2015. We prioritize the longer Fig. 1 Time-series of the ITF volume transport (1 Sv = 10 m 6 s −1 ) based on the ORAS4. The transport is calculated by vertically integrating the monthly mean velocity from the surface to the depth above the bottom on a line from Australia to Indonesia between 9º S and 20º S along 115º E. A 13-month (85-month) running mean is applied in the thin curve (thick curve) 1 3 uniform quality of wind stress data over the 4 decades to avoid artificial errors in the estimated SSHA. We note that our inter-comparison of the wind stress decadal variability over the North Pacific are mostly identical among the ORA-S3 and the ORAS4.
We used the CMIP5 dataset to analyze future changes in the behavior of decadal ITF variability. Among 23 models archived in the CMIP5 datasets, listed in Table 1, we eventually chose only 13 models with bold names. As we will describe in Sect. 4, the simulated SSH difference between the NWP and SEI regions in these selected models can represent a certain amount of variance in decadal ITF variability in historical simulations of CMIP5 models. We also used the Representative Concentration Pathways 8.5 (RCP8.5) scenario simulations (Moss et al. 2010). Future changes are evaluated as differences between 2049-2095 and 1959-2005. The monthly mean SSH and wind stress products in all employed datasets are re-gridded onto a 1° longitude × 1° latitude grid. Monthly anomalies as departures from the respective climatologies were smoothed by applying an 85-month running mean to highlight decadal variations.

ITF volume transport anomalies calculated from direct and estimated SSHA differences
We will describe procedures to estimate ITF volume transport anomalies in the reanalysis datasets and coupled models. Although these assimilated and simulated datasets contain velocity profiles in ocean surface layers as output variables, the procedures in this study only need a largescale SSHA distribution rather than anomalous velocity profiles in individual Indonesian straits.
Following previous studies, we employed the same technique of estimating the ITF variation using an anomalous SSH difference between the Pacific and Indian Oceans on a decadal timescale. Based on the assimilated and simulated datasets, we can obtain two different sorts of SSHA fields, such as "the direct SSHA" and "the estimated SSHA" in the following ways. As a direct SSHA field, we used the SSH field provided as the output variables in each dataset. While the direct SSHAs in the coupled models are the result of two-way communications between the ocean and atmosphere with considering not only wind forcing but also heat exchanges at the surface as well as vertical thermal structures in the ocean, the estimated SSHA is the result of one-way communication where the ocean is forced only by wind stress with a certain parameter setting of the ocean. By comparing the direct SSHAs in the coupled model and estimated SSH in the reduced gravity model, we can evaluate the contribution of wind forcing to the SSH variability in the NWP and SEI regions. Sections 3.1 and 3.2 explain the ITF estimation methods from the direct and estimated SSHAs, and in Sect. 4, we select a set of CMIP5 models after validating them by similar reproducibility in simulating the direct and estimated SSHAs in the NWP and SEI regions.

Decadal ITF variations by the direct SSHA difference
As discussed in the introduction, the SSHA difference across the Indonesian Archipelago is a good indicator of ITF variations on the annual and interannual time scales. The regression coefficients between the vertically integrated meridional velocity and SSHA variations on a decadal timescale were calculated using the ORAS4 dataset to determine whether it is still a good indicator in representing the decadal ITF variations (Fig. 2). A pair of significant positive regressed anomalies in the NWP region and negative ones north of 20º S in the Indian Ocean are associated with the intense ITF transport, vice versa, showing that Wyrtki's theory is applicable to the decadal time scale.
In the entrance and exit regions of the ITF, the SSHA is significantly correlated with ITF variations. Hence, the NWP (6-16º N, 125-155º E) and SEI (6-16º S, 85-115º E) regions are selected to calculate the SSHA difference as a proxy index ITF variability. In addition to the SEI region, the EEIO region (5º S-5º N, 65-95º E) is also selected for calculating SSHA differences. This additional selection is because the SSHA in the EEIO region is highly correlated with the decadal variability of the ITF transport, although the regressed anomalies in the EEIO region are not as much as in the SEI region. The selected regions are represented as the inset black rectangles in Fig. 2. SSHA variations from 1993 to 2015 based on ORAS4 and AVISO are mostly identical in these three regions, respectively (Fig. 3), although there is a slight systematic bias in mean values. This confirms that the reanalysis data can reproduce SSHA variations observed by satellites, apparently owing to the successful as simulation of observed data in the ORAS4. Figure 4 shows that the SSHA difference between the NWP and SEI regions correlated well with decadal ITF volume transport (correlation coefficient, R = 0.95) compared to the SSHA difference between the NWP and EEIO regions (R = 0.93). Although no significant difference between the two Rs, in subsequent sections, we selected the SEI region for the exit region to calculate the SSHA difference. Following Wyrtki's (1987) idea, the SSHA difference was linearly regressed onto the ITF variation. The resulting linear regression relation is expressed as follows:

Decadal ITF variations by the estimated SSHA difference
Besides using the direct SSHA as an output variable in the reanalysis datasets, the estimated SSHA fields are calculated by a 1.5-layer reduced gravity model forced by wind stress with the boundary condition along the coast to the east. This method for estimating the SSHA difference using only atmospheric data has a substantial advantage in overcoming the challenges of lower spatial resolution in the ocean component of CMIP5 models. Since previous studies (Abe et al. 2014;Vivier et al. 1999;Watanabe et al. 2016) confirm the reduced gravity model can reproduce observed SSHA, including the feature of the westward propagating baroclinic Rossby waves indicated by the satellite measurements over the tropical Pacific, we used the same model as used in Abe et al. (2014). A baroclinic Rossby wave in this model propagates across the ocean basin within 1-2 years. Hence, the temporal lags at the western boundary after the formation by wind Fig. 3 a Time-series of the decadal SSHA (cm) over the NWP region based on the ORAS4 (black thick curve) and AVISO (blue thin curve). b As in a, but for the SEI region (solid curves) and the EEIO region (dashed curves) (a) (b) Fig. 4 Time-series of the decadal ITF transport anomalies (Sv, black line) and the SSHA difference (cm) across the NWP-SEI regions (blue solid curve), and across the NWP-EEIO regions (blue dashed curve) on a basis of the ORAS4 reanalysis. The estimated SSHA dif-ference across the NWP-SEI regions is also superimposed (light blue solid curve). Note that the axis for the ITF transport is at the left of the panel, but the axis for the SSHA difference is at the right forcing can be negligible, because we will examine decadal time scales.
In the present study, we use monthly wind stress anomalies in the reanalysis as well as in the CMIP5 models to force the 1.5-layer reduced gravity model all over the Pacific basin with the boundary condition along the eastern coast in the 6-16º N latitude band to reproduce the SSHA in the NWP region (6-16º N, 125-155º E). Our inter-comparison of the estimated SSHAs by the reduced gravity model (blue curve of Fig. 5a) with the SSHAs using satellite observations (red curve of Fig. 5a) and direct output of the ORAS4 (black curve in Fig. 5a) shows that the decadal SSH variations of the NWP region are successfully reproduced in the reduced gravity model. Due to the long distance of wave propagation toward the NWP region, the boundary condition along the eastern coast has no substantial effect on the estimated SSHA in the NWP region (Fig. 5b). This means that the SSHA in the NWP region can be estimated without the coastal boundary condition as the ocean component output.
By contrast, the SEI region is close to the eastern coastal boundary. Therefore, the reduced gravity model will require the wind stress curl and the boundary condition for a better estimation of the SSHA. As expected, both the wind forcing and the boundary condition contribute to reproducing the observed SSHA in the SEI region (Fig. 6a). In contrast, either alone is unable to do so (Fig. 6b).
However, our intention for the estimated SSHA is to seek the appropriate procedure without using the SSHA information for the calculations of the reduced gravity model. Schouten et al. (2002) show that an eastward propagating Kelvin wave caused by equatorial zonal wind stress over the Indian Ocean controls SSH variations along the eastern coast in the SEI region. Indeed, the SSH variations of the SEI region are significantly correlated with zonal wind stress anomalies over the central portion of the equatorial Indian Ocean (5º N-10º S, 65-105º E) (Fig. 7). Hence, rather than using the boundary condition directly from the satellite observation and/or reanalysis dataset, we calculated a linear regression of the equatorial zonal wind stress onto the contribution of the boundary condition to the SSHA in the SEI region using decadal anomalies of those variables during a period of 1959-2005, as shown below where h b (black curve of Fig. 8) is a contribution of the coastal boundary condition onto the estimated SSHA in the SEI region, and x is 85-month running mean anomalies of zonal wind stress over the central portion (5º N-10º S, 65-105º E) of the Indian Ocean.
Consequently, there is a two-step procedure to obtain the estimated SSHA in the SEI region. First, the reduced gravity model calculates the wind-forced component. Then, as the second step, the contribution of the coastal boundary condition from the zonal wind stress anomalies using the regression analysis is added to the wind-forced component. Our inter-comparison shows that the estimated SSHA by this two-step procedure (blue curve in Fig. 6a) is mostly identical to the direct SSHA (black curve in Fig. 6a).
The differences in the estimated SSHA between the NWP and SEI regions are used to represent the ITF variations, as discussed in Sect. 3.1. Figure 4 compares the total ITF transport of the ORAS4 reanalysis data (see also Fig. 1) with direct and estimated SSHA differences between the NWP and SEI regions. The significant correlation coefficients between those time-series show that the estimated SSHA proposed in this subsection is useful to represent the decadal variations in the ITF volume transport based solely on atmospheric variables. We will apply this understanding and procedure to the CMIP5 dataset, as documented in later sections.

CMIP5 model selection
One of this study's objectives is to understand future changes in the relative contribution from the Pacific and Indian Oceans to decadal ITF variation. For this purpose, we used historical and RCP8.5 scenario simulations as the current and future climate realizations, respectively. However, the climate model simulations for this objective are challenging, because model reproducibility of the ocean components in the marginal seas, including the SSH variations around the Indonesian Archipelago, is not always guaranteed, primarily due to insufficient horizontal resolution and associated unrealistic bottom topography. Hence, we should select appropriate models for analyzing decadal ITF variations.
Our metric for model selection is based on the mean SSH difference between the NWP and SEI regions. As listed in  Fig. 7 A map of correlation coefficients on a decadal timescale (contoured for every 0.2) in zonal wind stress anomalies in ORA-S3 onto the effect of boundary condition SSHAs averaged over the SEI region on a basis of the ORAS4. The coloring convention is shown on the right of the panel Table 1, there is a large inter-model spread in reproducing the mean SSH difference between the 23 CMIP models. While 21 models represented a downward SSH gradient from the Pacific to Indian oceans, the other two showed an upward gradient. A multi-model ensemble mean using 21 models with a realistic downward gradient is 16.7 cm, which is close to the observed mean (18.4 cm). An inter-model standard deviation among 21 models is 6.0 cm. We eventually retained thirteen (13) of the CMIP5 models screened by the mean SSH difference in a range of one standard deviation to the multi-model ensemble mean (10.7-22.6 cm).
These selected models are then evaluated for their ability to reproduce comparable amounts of variance in both direct and estimated SSHA differences. This direct SSHA in the simulations is based on the model outputs. The estimated SSHA in the simulation is calculated with the reduced gravity model forced by the simulated surface winds, as conducted with the observed datasets in Sect. 3. Figure 9 shows the decadal time-series of the SSHA differences between the NWP and SEI regions using both direct and estimated SSHAs from the selected thirteen (13) CMIP5 models. While the decadal variability of the individual time-series has amplitude modulations on the interannual time scale, displaying the significant decadal covariation of direct and estimated SSHA differences with small root-mean-square errors in the individual selected models indicates that the simulated surface winds mostly determine the simulated SSHA difference across the Indonesian Islands over the Indo-Pacific regions.

Future changes in the relative contributions from the Pacific and Indian Ocean SSHA variability onto the decadal ITF variations
To examine relative contributions of decadal SSH variability in the NWP and SEI regions to their SSHA difference, respectively, as the indicator of the ITF variation, we defined the positive ( Figure 10a shows a difference composite map of SSHAs between the positive and negative phases for the ORAS4 dataset. As expected from the definition, we found a coherent pattern with higher SSHA in the western North Pacific than in the Indian Ocean. While this large SSHA difference across the Indonesian Islands is mostly attributed to the SSHA over the NWP region, the relative contribution of the SEI region is approximately one-sixth of the contribution from the NWP region. Thus, the NWP region contributes to pushing water from the Pacific into the Indian Ocean. Meanwhile, the smaller contribution of the SEI region indicates that drawing the water from the Pacific is weak. Composite SSHAs over 6.0 cm in the NWP region (Fig. 10a) are about 40-50% of the regressed SSHA in the same region (Fig. 2), which contributes to 0.4-0.5 Sv of the decadal ITF variability, although their statistical significance is not so high owing to the small degree of freedom. We note that the positive SSHA in the NWP region accompanies with zonally elongated significant positive SSHA over the zonal band at around 5-10ºN across the North Pacific basin.
A similar SSHA pattern is found in the difference composite map based on the historical simulations of the selected CMIP5 models (Fig. 10b), even though positive anomalies over the NWP and negative ones over the SEI regions are weaker than those of the ORAS4 dataset (Fig. 10a). While the spatial correlation coefficient between Fig. 10a and b is 0.55, the correlation becomes significant (0.62) when the domain is limited in 20ºS-30ºN (upper row in Table 2). This significant spatial correlation indicates that the simulated SSHA pattern Fig. 8 Time-series of the boundary condition component (blue dashed curve) in the estimated SSHA by the reduced gravity model (same as in Fig. 6b) and the calculated boundary condition contribution to the SSHA in the SEI region using a regression analysis with the equatorial (5º N-10º S, 65-105º E) zonal wind stress of the ORA-S3 (solid black curve) 1 3 associated with decadal ITF variations based on the selected CMIP5 model is adequately reproduced. As in the ORAS4, the large SSHA difference across the Indonesian Islands is mostly attributable to the SSHA over the NWP region. The SEI region's relative contribution is approximately one-sixth of the contribution from the NWP region, which fairly agrees with the result of the ORAS4. As shown in Fig. 10a and b, the direct and estimated SSHAs over the NWP region in the CMIP5 historical simulation are underestimated than in the ORAS4 (first and second rows of Table 3). As indicated in Table 1, the climatological mean SSH difference is smaller in ten models among the selected 13 CMIP5 models than in the ORAS4. This implies that the contribution from the NWP tend to be underestimated in most of the selected models. This underestimation in the historical simulations is mostly attributed to the weaker atmospheric wind forcing in the simulated mean states (Lyu et al. 2016). The spatial correlation in the difference composite maps of the SSHA is still high (second row in Table 2) between the RCP8.5 (Fig. 10c) of the selected CMIP5 models and the ORAS4 as well as between the historical and RCP8.5 simulations (bottom row in Table 2). A closer look at the SSHA magnitudes, specifically over the NWP and SEI regions, shows decreased changes in the relative contributions of the direct SSHA over the NWP region from 3.34 cm in the historical simulations to 2.50 cm in the RCP8.5 simulation (a 25% decrease) as well as the increased changes over the SEI region from − 0.51 to − 0.75 cm (a 47% increase). Comparable values of the changes are also found in the estimated SSHA from 2.26 to 1.67 cm in the NWP region (a 26% decrease) and from − 0.71 to − 1.00 cm in the SEI region (a 41% increase). Table 3 summarizes these changes in the relative contributions in the multi-model ensemble mean and individually selected models. We should note that these comparable results in the direct and estimated SSHAs in the individual models and the model ensemble mean strongly indicate that future changes in the relative contributions of the SSHAs across the Indonesian islands are attributable to changes in the atmospheric forcing.

Future changes in atmospheric forcing
As displayed in Fig. 9, identical features in the decadal variability of the direct and estimated SSHAs indicate that the decadal surface wind variability is a major driver to induce the direct SSHAs in Fig. 10. In addition, the comparable values of the direct and estimated SSHA in Table 3 additionally prove this notion. Hence, to discuss the causes of future changes in the relative contributions of the SSHAs on either side of the Indonesian Islands, we show the difference composite maps of the multi-model ensemble mean wind stress curl anomalies over the North Pacific (Fig. 11) and zonal wind stress anomalies over the Indian Ocean (Fig. 12). Negative wind stress curl anomalies are significant over the zonal band in 6-16º N over the North Pacific in the historical simulations (Fig. 11a), which is associated with anticyclonic surface wind anomalies over the subtropical North Pacific (not shown here). Indeed, the weak anticyclonic surface winds over the subtropical Pacific were illustrated during the positive phase of the PDO (Mantua et al. 1997). These negative wind stress curl anomalies become much weaker in the RCP8.5 simulation than in the historical climate simulation (Fig. 11). This reduced magnitude of wind stress curl anomalies over the NWP region will reduce the magnitudes of wind-forced Rossby waves, which mostly explain the reduced relative contribution of direct and estimated SSHAs over the NWP region during the future climate.
As displayed in Fig. 7, negative zonal wind stress anomalies over the central equatorial Indian Ocean significantly contribute to the negative SSHA over the SEI region through the wind-forced equatorial Kelvin wave. The difference composite map of the zonal wind stress anomalies over the Indian Ocean shows enhanced easterly wind stress anomalies over the central equatorial Indian Ocean in the RCP8.5 simulation than in the historical simulation by + 23% over the SEI region, which results in stronger negative SSHA. Compared to historical simulations, the analysis also shows the eastward shift of zonal wind stress anomalies (Fig. 12) in the future climate scenario. While the zonal wind stress anomalies center is around 70º E longitude in the historical climate simulation, it is shifted eastward around 82º E in the RCP8.5 simulation.
The combined effect of increased amplitude of zonal wind stress anomalies and eastward shift in the anomalies can explain the increased SSHA over the SEI region by + 47% in the direct SSHA and + 41% in the estimated SSHA (Table 3) and, eventually, increased contribution from the SEI region to the decadal variation of ITF transport.

Conclusion
CMIP5 models project consistent changes in the climate system due to global warming, which may affect the decadal perturbations from the mean state. Hence, we examined the decadal ITF variations and relative contribution of decadal SSHA variations over the NWP and SEI regions to the ITF in the current and future climate conditions. Fig. 9 Time-series of the decadal SSHA difference (cm) across the NWP-SEI regions, based on the simulated direct SSHA in the selected CMIP5 models (black curve), and the estimated SSHA (red curve) using the simulated wind stress of those models. Simultaneous correlation coefficients between those direct and estimated SSHAs are indicated on the upper right corner in each of the panels ◂ Decadal timescale analysis of ITF requires long-term datasets, which are unavailable. Hence, a new method is required for estimating the long-term variability of ITF. Here, the ITF volume transport was calculated based on monthly mean velocity integrated from the surface to the depth above the bottom on a line from Australia to Indonesia between 9º S and 20º S along 115º E, based on the ORAS4 reanalysis data.
The decadal variability of ITF is captured by calculating 85-month running mean anomalies that show significantly high correlation (r > 0.9) with the SSHA difference between the NWP and SEI regions. This confirms that the SSHA difference is capable of reproducing the ITF variation on decadal timescales. Further, we reconstruct ITF's decadal variability based solely on wind stress data over the Pacific and Indian Oceans, avoiding the dependence on a Fig. 10 a Difference composite maps of the decadal SSHA (contoured for every 1.0 cm) based on the ORAS4 dataset. b As in a, but for the multi-model ensemble mean decadal SSHA in the historical simulations of the selected CMIP5 models. c As in a, but for the decadal SSHA in the RCP8.5 simulations of those models. The coloring convention is shown on the right of the panels, respectively (a) (b) (c) complex climate model. The SSHAs were simulated using a simple 1.5 reduced gravity model and 46 years of monthly wind stress data from the Pacific and Indian Oceans. The connection with the ITF volume transport in the Indonesia straits was built using statistical analysis. The correlation coefficient between reanalysis and estimated ITF based on wind stress in this period is greater than 0.8. However, the SSHAs' relative contribution to the ITF variability on either side of the Indonesia strait is not the same due to the large difference in the SSHA variance. Hence, the contribution of the decadal variability of the Indo-Pacific SSHA on ITF was investigated using ORAS4 reanalysis data from 1959 to 2005.
The contributions from the NWP and SEI regions were estimated based on the difference composite anomalies of SSHAs over the NWP and SEI regions between positive and negative phases of the ITF. Our analysis of the historical period  concludes that the SSHA over the NWP region drives ITF variability on a decadal timescale. The SEI region plays a negligible role. Furthermore, the changes in the contribution from the NWP and SEI regions were studied using thirteen (13) CMIP5 models in historical simulations  and RCP8.5 future scenarios . While the analysis of SSHA suggests that anomalies over the SEI region should increase significantly by + 47% (direct SSHA) and + 41% (estimated SSHA) in the future climate, the anomalies over the NWP region show a significant decrease by − 25% (direct SSHA) and − 26% (estimated SSHA). These changes in the anomalies are equivalent to the changes in the contribution from the NWP and SEI regions.
These changes in the contribution from direct/estimated SSHAs over the NWP and SEI regions during future climate scenarios  are driven by changes in the wind stress curl over the Pacific and equatorial zonal wind stress over the Indian Ocean. A decreased amplitude of wind stress curl anomalies over the off-equatorial Pacific decreases the SSHA and, therefore, a contribution from the NWP region. While the increased anomalies of zonal wind stress over the equatorial Indian Ocean contribute to the increased contribution from the SEI region, the eastward shift of the anomalies over the equatorial region also contributes to the increased contribution changes from the SEI region to decadal ITF transport.
These changes in the contribution from the NWP and SEI regions to the decadal variation of ITF have significant importance. Previous studies and our analysis have also shown that in the current/historical climate, ITF can be estimated only with the SSHA over the NWP region without considering the variation in the SEI region. However, with global warming and climate change, the projected contribution changes from the NWP and SEI regions to ITF variability indicate that estimating ITF is difficult without considering future climate variation in the Indian Ocean. While this implication based on the analysis only for the 13 selected models, there is a substantial diversity of the simulated SSHA difference between the NWP and SEI regions among all of the CMIP5 models (Table 1). A rationale on this diversity in the simulations including a reserved SSHA gradient should be investigated to further improve the performance in representing the future climatological mean and decadal variability in the climate models. Table 2 Spatial correlation coefficients of difference composite SSHA maps among the three datasets such as the ORAS4 reanalysis dataset (Fig. 10a), the multi-model ensemble mean of the historical simulations (Fig. 10b), and one of the RCP8.5 simulations (Fig. 10c) from the selected CMIP5 models While values in the middle column indicate the spatial coefficients over the whole domain of the figures (30º S-30º N), the ones in the right column are over a domain only in 20º S-30º N  . 11 a Difference composite maps of the decadal wind stress curl anomalies (contoured for every 0.2 N m −3 ) based on the multi-model ensemble mean of selected CMIP5 historical simulations. The coloring convention is shown on the right of the panels. b As in a, but for the simulations under the RCP8.5 scenario. c is the difference between b and a (a)

(b) (c)
program for analysis and graphics in this manuscript. We also would like to thank Vineet Singh and Sophia Jacob from Indian Institute of Tropical Meteorology for helping download data whenever needed.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.