An Updated Estimate of the Indonesian Throughflow Geostrophic Transport: Interannual Variability and Salinity Effect

While the pivotal role of the Indonesian Throughflow (ITF) in the global ocean circulation and climate has been widely recognized, accurate estimation of its volume transport based on observational data remains challenging. This work provides an updated estimate of the monthly ITF geostrophic transport (ITFG) in the upper 700 m at the IX1 section between Indonesia and Australia. The ITFG estimate is better constrained by improved data correction and new salinity data products. The mean ITFG of 1993–2018 is 8.2 ± 0.2 Sv, in which contributions of the temperature and salinity components (ITFT and ITFS) are 5.5 ± 0.2 and 2.8 ± 0.1 Sv, respectively. The interannual variability (∼4.4 Sv in standard deviation) is dominated by ITFT, as the result of wind‐driven thermocline dynamics, and slightly attenuated in amplitude by ITFS. The strengthening trend of 1.33 Sv decade−1 during 1993–2018 primarily arises from ITFS and secondarily from ITFT.

• An updated estimate of the Indonesian Throughflow (ITF) is provided, with improvements in data correction, uncertainty estimation, and salinity effect • The ITF in the upper 700 m is of a mean geostrophic transport of 8.2 ± 0.2 Sv and interannual variability of ∼4.4 Sv • Salinity exerts a weak damping effect on ITF's interannual variability but fundamentally contributes to its strengthening trend since 1993 Supporting Information: Supporting Information may be found in the online version of this article.
sampling. These data capture the upper-layer geostrophic volume transport into the Indian Ocean, the leading component of ITF outflow transport (Potemra et al., 2002), with other components broadly including Ekman, barotropic, ageostrophic, and deep-layer transports. However, existing estimates of the ITF geostrophic transport (ITF G ) at the IX1 line ranged widely between 2 and 12 Sv (e.g., Liu et al., 2015;Meyers, 1996;Meyers et al., 1995b;Wainwright et al., 2008;Wijffels et al., 2008;Yuan et al., 2013). In addition to pronounced variability, uncertainties also arise from potential data errors and the approach. One major error source is the well-known XBT data bias identified soon after XBT manufacture began (Cheng et al., 2016).
Recently, a suite of correction schemes has been proposed to address this issue (Gouretski & Reseghetti, 2010;Levitus et al., 2009), with that of Cheng et al. (2014) the community-recommended method.
Another issue is the ocean salinity effect (Hu & Sprintall, 2017;Liu et al., 2005). Owing to the sparse historical sampling of ocean subsurface salinity, most previous estimates either adopted the climatological temperature-salinity relationship or ignored salinity effect in the geostrophic calculation (e.g., Liu et al., 2015;Meyers, 1996;Meyers et al., 1995a;Yuan et al., 2013). However, recent studies suggested that changes in ocean salinity have notable impacts on the upper-ocean circulation and sea level of ITF's exit region Llovel & Lee, 2015;Lu et al., 2022;Menezes et al., 2013). Hu and Sprintall (2016) suggested that the salinity component contributed ∼36% to the ITF interannual variability and dictated its strengthening since ∼2004. The recently available observational salinity data products Tian et al., 2022) may promote the understanding of the salinity effect.
This study updates the ITF G estimation at the IX1 line. By adopting the recently available data correction and quality control schemes and ocean salinity datasets, errors in the ITF G estimate are reduced. The uncertainty associated with data sampling is also provided through a Monte Carlo simulation approach, allowing the propagation of uncertainties from various data sources to the final ITF G estimate.

Data and Processing
Historical temperature (T) data near the IX1 line include profiles from XBTs, mechanical bathythermographs (MBT), conductivity-temperature-depths, bottles, moored buoys, gliders, and Argo floats. These data, in all 764,481 profiles with XBT being the dominant type, have undergone the XBT bias correction of Cheng et al. (2014), MBT correction by Gouretski and Cheng (2020), and quality control of Tan et al. (2023), yielding unified profiles with 1-m intervals between 0 and 100 m and 5-m intervals between 100 and 700 m. Further manual quality control was conducted to remove suspicious data (∼2.32% of the total). Then, T profiles are reorganized into grid boxes along the IX1 line (Figure 1a), following Liu et al. (2015). Our estimate of ITF G starts from 1984 when repeated XBT deploys induced a significant increase (Figure 1b). The monthly T climatology is computed by averaging all available data in each box bin, and then the T anomaly field relative to the climatology is obtained ( Figure S1 in Supporting Information S1). For each T anomaly member generated through the Monte Carlo simulation (Section 2.3), linear interpolation and 13-month low-pass filtering are applied to fill the data gap and reduce subseasonal fluctuations. Adding the filtered anomaly back to the climatology, the reconstructed T field ( Figure 1c) contains seasonal, interannual, and lower-frequency variations.
To account for the ocean salinity (S) effect, we used two versions of the Institute of Atmospheric Physics (IAP) monthly salinity gridded datasets Tian et al., 2022). One is the 25-member, 0.25° × 0.25° product of 1993-2018, and the other is the 26-member, 0.50° × 0.50° product of 1984-2017. These data are interpolated onto the same depth levels of the T profiles. In-situ salinity observations, sea surface temperature, surface winds, altimeter-based sea surface height, and a coarse-resolution gridded salinity product are merged to reconstruct the monthly S fields (Tian et al., 2022). The feed-forward neural network and the Monte Carlo dropout approach are used in member generation to ensure that the reconstructed members account for all major error sources including salinity instruments, representativeness, and methodology (Tian et al., 2022). These two S datasets are interpolated onto the IX1 line to match T (Figure 1d), forming two versions of the ITF G estimate: the 1993-2018 and 1984-2017 versions ( Figure S1 in Supporting Information S1). The former is of higher quality, owing to improved T data sampling (Figure 1b), and the following analysis is mainly based on this version.

Geostrophic Calculation
The dynamic height is calculated by integrating the specific volume anomaly upward from the 700 m reference level Yuan et al., 2013). Then, the cross-section geostrophic velocity of the ITF (V G ) is calculated through = 1 ( 2 − 1) , where L is the distance between neighboring box bins with the dynamic heights of D 1 and D 2 along the IX1 line, and f is the Coriolis parameter. The dynamic height is calculated by vertically integrated specific volume anomaly at pressure levels, possessing the unit of m 2 s −2 (Yang & He, 2014). Velocities toward the Indian Ocean (Indonesian Seas) are defined positive (negative). The calculated V G clearly shows multiple currents (Figure 1e; Text S1 in Supporting Information S1).
Integrating V G over 0-700 m of the IX1 section yields the total ITF G , which is further linearly decomposed into the temperature (ITF T ) and salinity components (ITF S ). Practically, in the calculation of ITF T , the dynamic height is computed using the average S profile of all along-section box bins, and similarly for the ITF S . In climatology, the temperature (V T ) and salinity components (V S ) are both essential in shaping the structure of V G ( Figure S3 in Supporting Information S1).

Generation of Ensemble Members and Uncertainty Estimation
The Monte Carlo approach is employed to quantify uncertainties, ensuring propagation of all data set errors into the final ITF G estimates, without making strong statistical assumptions on the error structure. The 100 realizations of ITF G are generated by randomly selecting one T anomaly profile for each box bin and selecting one member of the IAP S fields. Generating >100 realizations leads to similar results. Each member of the ITF G provides a consistent set of variability, with their ensemble mean representing the statistical average of the ITF G . The 90% confidence interval, given by the 5th and 95th percentiles, is the overall measure of uncertainties. This approach is also used for ITF T and ITF S . For ITF T , we randomly select T anomaly profiles for individual box bins to construct the T field, while for ITF S , we randomly select one member of the IAP S field (25 and 26 members for the 1993-2018 and 1984-2017 versions, respectively). Figure S1 in Supporting Information S1 provides the schematic view for the procedures described above.

Geostrophic Transport
The estimated ITF G shows prominent seasonal and interannual variabilities superimposed atop an increasing trend ( Figure 2a). The average ITF G of 1993-2018 is 8.2 ± 0.2 Sv (error representing the data uncertainty), close to Liu et al. (2015). This transport, even considering the ∼0.7 Sv Ekman transport (Figure 2d), is weaker than the mooring-based estimates of ∼14 Sv in the Makassar Strait (Gordon et al., 2010) and outlet straits (Sprintall et al., 2009). This discrepancy points to the limitation of our approach, given that deep-layer transport, ageostrophic flow, and the barotropic component are precluded. In addition, the data sampling at IX1 may lead to underestimation of the narrow currents along the Indonesian and Western Australian coasts.
The ITF T shows a mean value of 5.5 ± 0.2 Sv, closely resembling ITF G (Figure 2b). The ITF S shows a mean value of ∼2.8 ± 0.1 Sv (Figure 2c), explaining 34.2% ± 1.2% of the total ITF G . The uncertainties are typically one order of magnitude smaller than the transports. The ITF G shows an increasing trend of +1.33 Sv decade −1 during 1993-2018, significant at the 90% confidence level. This strengthening trend has been reported by previous studies based on IX1 data (Feng et al., 2018;Liu et al., 2015), proxies (Hu & Sprintall, 2017;Sprintall & Révelard, 2014), and model hindcasts (Feng et al., 2017;Lee et al., 2015;Li et al., 2017Li et al., , 2018 and attributed primarily to the enhanced Pacific Walker circulation (England et al., 2014). This trend occurs predominantly in ITF S (+0.81 Sv decade −1 ) and secondarily in ITF T (+0.48 Sv decade −1 ), broadly consistent with Hu and Sprintall (2017). The ITF S upward trend is probably induced by the strengthened precipitation around southeastern Indian Ocean and Indonesian Seas and ocean advection (Du et al., 2015;Guo et al., 2021;Hu & Sprintall, 2017). This trend is also examined in different density layers within 0-700 m ( Figure S4 in Supporting Information S1), and the significant strengthening of ITF S is uniformly seen in all layers ( Figure S4b and S4c in Supporting Information S1).
The seasonality of ITF G ( Figure S5a in Supporting Information S1) is characterized by a peak during the southeast monsoon (June-September; 12-14 Sv) and weakened transports during the transition seasons (April-May and October-November; 3-6 Sv), partly associated with the wind-driven SJC reversal (Meyers, 1996;Sprintall et al., 1999Sprintall et al., , 2010. This seasonality arises mainly from ITF T . The ITF S diminishes in January-February and generally fluctuates between 2 and 5 Sv in other months. The westward ITF E is enhanced to 2-3 Sv during the southeast monsoon, while an eastward ITF E of ∼1 Sv emerges in December-February (DJF) ( Figure S5b in Supporting Information S1). The ITF E from QuikSCAT is slightly stronger ( Figure S5b in Supporting Information S1). The analysis focuses on ITF G and its components.

The Temperature Component
The ITF T dominates the interannual variability of ITF G (Figure 3a). It shows a standard deviation of 4.7 Sv and a correlation of 0.97 with ITF G . Interannual variability in our estimation is weaker than Hu and Sprintall (2016) based on a proxy approach, likely arising from the differences in data source, methods, and time periods. Variations in ITF T are primarily the manifestation of the thermocline dynamics Meyers, 1996;Xie et al., 2002). In climatology, the along-section gradient of the thermocline depth (e.g., the depth of the 20°C isotherm; D20) is essential in shaping the overall westward V T (Figure 1c and Figure S3a in Supporting Information S1). We discern contrasting variability characteristics in D20 northern and southern portions of IX1 section, separated by ∼10 o S (Figure 3b). A depressed thermocline in the south enhances the thermocline tilt and leads to an enhanced westward ITF T , as in the La Niña episodes of 1999-2001 and 2010-2012, and vice versa for the El Niño episodes of 1993-1995, 1997-1998, 2003-2005 Figure S2 in Supporting Information S1 for Niño-3.4). The total ITF T shows a higher correlation with ENSO than with IOD ( Figure S6a in Supporting Information S1). The maximal correlation with the DJF Niño-3.4 is −0.74, occurring in the spring following the ENSO mature phase. This phase lag likely reflects the transmission of planetary and coastally trapped waves from the Pacific to Indian Ocean (Cai et al., 2005;Wijffels & Meyers, 2004).
D20 anomalies northern end of the IX1 are driven by IOD winds (Chen et al., 2016;Meyers, 1996;. The ITF T north of 10°S shows a significant correlation of −0.50 with the September-November (SON) DMI at a 2-month time lag, higher than the −0.39 correlation with Niño-3.4 (Figure 3c), while the ITF T south of 10°S is mainly dictated by ENSO (Figure 3d). These correlations likely suggest that the ITF T is generally enhanced during negative IOD and La Niña events and weakened during positive IOD and El Niño events. However, the negative correlation between IOD and ITF T arises partly from the relationship that the El Niño condition favors the development of positive IOD (e.g., Cai et al., 2009). Existing studies suggested that a positive IOD acts to enhance the ITF by driving upwelling coastal Kelvin waves near the Sumatra-Java coasts (Chen et al., 2016). However, this IOD effect is generally not discernible owing to the dominance of the ENSO effect.

The Salinity Component
Although ITF S displays sizable interannual variability with a standard deviation of ∼1.0 Sv, its phase is generally opposite to ITF G (Figure 4a). Overall, ocean salinity acts to slightly dampen ITF G anomalies, such as the 2000 maximum and the 1998 and 2013 minima. In the upper-layer (0-700 m average) salinity (Figure 4b), we again notice different variability characteristics north and south of ∼14.5°S. The northern sector exhibits weaker shortterm variations that persist typically for 1-2 years, while the southern sector shows stronger and lower-frequency variations. Salinities of both sectors can modulate the along-section dynamic height gradient and ITF S . The salinity contrast ∆S, as the south-minus-north salinity difference, is a useful proxy for understanding this issue, yielding a correlation of −0.73 with ITF S (Figure 4c). A salinity increase in the south or a freshening in the north causes a positive ∆S anomaly and an eastward ITF S anomaly by creating a northward dynamic height gradient along the section.
The total ITF S shows a significant correlation of +0.67 with ENSO and a low correlation of +0.21 with IOD ( Figure S6b in Supporting Information S1). These correlations are particularly tight for the salinity north of 14.5°S, which are +0.76 and +0.50 with Niño-3.4 and DMI, respectively (Figure 4d). Either an El Niño or a positive IOD acts to cause a salinity increase north of 14.5°S and enhance ITF S in ∼4-month time lag through both anomalous rainfall and wind-driven ocean dynamics (Cai et al., 2009;Chen et al., 2023;Masson et al., 2004;Zhang, Feng, et al., 2016;Zhang, Du, & Qu, 2016). For example, the freshening anomalies of 1999-2001 and 2008-2009 correspond to La Niña events, enhancing rainfall in the southeast tropical Indian Ocean (Zhang, Feng, et al., 2016); the positive IOD events in 1997 and 2006 caused local salinity increases via reduced rainfall, ocean advection, and coastal upwelling (Masson et al., 2004;Zhang, Du, & Qu, 2016). Note that the positive correlation with DMI partly arises from the covariance between ENSO and IOD. The positive correlation between ITF S and Niño-3.4 also explains the overall damping effect of salinity; an El Niño leads to an attenuated ITF T (Section 3.2) and an enhanced ITF S , with the former dominating the total ITF G change.
The Ningaloo Niño/Niña, the prominent surface temperature variability off the Australian coast, is seasonally phase-locked to DJF (Feng et al., 2013Kataoka et al., 2014). Some events develop without ENSO and IOD effects (Kataoka et al., 2018;Zhang et al., 2018) and represent an independent source for the ITF S variability.
The Ningaloo Niño/Niña shows a lower correlation with salinity north of 14.5°S (−0.41) than Niño-3.4 and DMI, implying limited influences (Figure 4d). The salinity south of 14.5°S shows a higher correlation (∼−0.5) with the NNI than Niño-3.4 and DMI (Figure 4e). The Ningaloo Niño (warm event) can induce large-scale freshening south of 10 o S primarily through enhanced local rainfall and secondarily through ocean current advection (Guo et al., 2021;Li et al., 2023).
The upward trend of ITF S of 0.81 Sv decade −1 indicates an overall increase of ∼72% of its climatological value (2.8 Sv) during 1993-2018. This notable trend is likely formed under multiple impacts. The interdecadal variations of the tropical Pacific climate (England et al., 2014) led to increases in the occurrence frequency and intensity of Ningaloo Niño Li et al., 2019;Tanuma & Tozuka, 2020) and a freshening trend near Australian coast during 1993-2015 (Du et al., 2015). Subsequently, the 2015-2016 super El Niño caused an abrupt salinification in the northern sector. These changes all contributed to the increasing trend of ITF S , particularly a strengthening after 2012. In comparison, a decreasing trend during 1999-2002 is responsible for a stronger freshening in the northern sector associated with La Niña and negative IOD events (Figure 4d). The signature of the 1997-1998 super El Niño was much weaker compared to that of 2015-2016, possibly owing to the insufficient data sampling prior to the Argo era. In addition to natural climate variability, externally forced change, such as amplification of the water cycle under greenhouse warming (Held & Soden, 2006), may also contribute to the ITF S trends (Hu & Sprintall, 2017). Moreover, we cannot ignore the caveat that the availability of Argo data after 2005 improves the sampling of subsurface salinity and better resolves the along-section slope of the halocline, favoring the ITF S trend since 1993. Further effort is required to disentangle the various impacts and better justify origin of the trend.

Remarks
This study provides an updated estimate of ITF G at the IX1 section by adopting advanced data bias correction, more realistically considering the ocean salinity effect, and uncertainty quantification. The mean ITF G is estimated to be 8.2 ± 0.2 Sv with an evident seasonality ranging from 3 to 14 Sv. The interannual variability of ITF G shows a standard deviation of 4.4 Sv and is dominated by ITF T , with the thermocline dynamics driven by ENSO and IOD winds being the central mechanism. The interannual variability of ITF S , caused mutually by ENSO, IOD, and Ningaloo Niño/Niña, is considerable in amplitude (∼1.0 Sv) and overall acts to slightly dampen ITF G variability. The 1.33 Sv decade −1 increasing trend of ITF G during 1993-2018 arises primarily from ITF S , dictated by salinity trends in the southeast Indian Ocean.
The 1984-2017 version ( Figure S7 in Supporting Information S1) agrees with 1993-2018 version in their overlapping period but shows larger uncertainty in ITF S prior to 1993. Our future effort is to refine the quality of the 1984-2017 version in analysis method to provide a more reliable long-term trend. Other components of the ITF transport, particularly the deeper-layer flow, should also be estimated. There is also a need to examine the heat and freshwater transports of the ITF and evaluate how these transports are simulated by numerical models.