Any variable in our study can in general be denoted by

where ȷ = 1 , … , N is a serial number of the gridpoints at which the variable is defined, l = 1 , … , L is the time index, and m = 1 , … , M is the ensemble-member number. We sometimes use time t instead of index l and write, e.g., t = 1986 to denote the year 1986. The ensemble anomaly is naturally

We calculate the linear regression coefficient between the ensemble anomaly of the annual mean ITF transport, q ˜ ( l , m ) , and some variables, γ ˜ ( ȷ , l , m ) , to be named in Section 3. Details are found in Appendix A.1. Linear regression extracts two parts from γ ˜ as

at each ( l , m ) , where α ( ȷ ) is the regression coefficient at gridpoint ȷ and γ ˜ ′′ ( ȷ , l , m ) is the residual. The term α ( ȷ ) q ˜ ( l , m ) is that part of γ ˜ which is maximally correlated with q ˜ ( l , m ) and γ ˜ ′′ ( ȷ , l , m ) is the part which is uncorrelated. In Figures 7 , 8 below, we compare γ ˜ ( ȷ , l , m ) with α ( ȷ ) q ˜ ( l , m ) , not only in terms of the spatial pattern but also in terms of values. If the two fields are similar, that is an indication that the correlated part dominates at the particular ( l , m ) .

Note that this statistical calculation critically depends on the “degrees of freedom,” which are at most LM − L because there are L constraints that the ensemble average of the data is zero at each year (Walker, 1940). If each annual mean is statistically independent between years and between members, this is the actual degrees of freedom. To try to find other internal dependency than the L constraints, we have calculated the temporal autocorrelation of q ˜ ( l , m ) for all members combined ( Supplementary Figure S2 ; see also Section 3.1) and found that the correlation is small for lags of 1 year and larger. This suggests that ensemble anomalies may be independent year to year. This conclusion, however, is tentative, and we acknowledge that the degrees of freedom of LM − L can be an overestimation. [There is indeed some hint of spectral peaks (Section 3.1), which, if real, suggests an overestimation.] The statistical significance of the regression coefficients we show below is based on the assumption of degrees of freedom of LM − L, but the significance could be an overestimation as a result.

We also use empirical orthogonal functions (EOFs) to explore characteristic spatial patterns of ensemble anomalies. Details are found in Appendix A.2. EOFs are orthogonal to each other; in other words, the spatial correlation between any two EOFs vanishes. From this property, the original variable, γ ˜ ( ȷ , l , m ) , can be expanded in terms of EOFs as

where e n ( ȷ ) is the n th EOF, commonly termed “EOF n ”, and ϕ n ( l , m ) is the coefficient of expansion, often called the n th “principal component”.

From the orthogonality of EOFs, we can derive ( Appendix A.2 )

The EOFs are customarily numbered in the order of decreasing contribution to the overall variance:

by definition. On the other hand, the order of contribution, ϕ n 2 , is generally different at each ( l , m ) ; for example, see Figure 6 below.

In Figure 6 below, we compare γ ˜ ( ȷ , l , m ) with ϕ n ( l , m ) e n ( ȷ ) for a particular mode n or for the sum of two modes, not only in terms of the spatial pattern but also in terms of values. If the two fields are similar, that is an indication that the particular EOF mode(s) dominate(s) at the particular ( l , m ) .

Figure 1 shows the transects across which transports are calculated, and the definitions of the transects are listed in Table 1 .

These transects are designed to completely enclose the Indonesian Seas in the model, where some of the shallow, tiny gaps between some islands have already been closed. We mainly use the four transects (1 through 4) on the Indian Ocean side to calculate the total ITF transport. We have also calculated the total ITF transport on the Pacific side and found that, on annual average, the total transport differs by 0.02–0.05 Sv and sometimes up to 0.07 Sv (not shown), indicating that the barotropic flow has little divergence or convergence at this time scale, as expected. The small error must be because we ignore the variability of sea level in the calculation of transports: If there is a correlation u ′ η ′ ¯ between the horizontal and temporal variabilities of horizontal velocity u near the surface and sea level η , where the overline denotes the temporal average and average along the transects, this volume flux counts toward the error. For an order-of-magnitude estimate, let us suppose that η ∼ 0.5   m , u ∼ 0.1   m / s , and length of transect ∼ 1 , 000   km , which gives a transport of 0.05 Sv, not inconsistent with the volume imbalance. Since the model uses a Boussinesq approximation (Pacanowski and Griffies, 2000), volume is perfectly conserved (Section 2.1), and therefore the only other potential source of imbalance is the area-integrated sea-level change ( ∫ d x   d y   ∂ η / ∂ t ¯ ), which must be on the same order or less.

We first integrate along the transect the velocity component normal to it

depending on the orientation of the transect, where x t or y t is its longitude or latitude. For transects 54 and 55, each of which consists of a meridional and a zonal segment, V is the sum of the zonal and meridional transports.

This integral gives the volume transport per unit depth across the transect as a function of depth. Upon vertical integration, it gives the total volume of transport across the transect. The sum of the transports for transects 1–4 is the total ITF transport on the Indian side and the sum over transects 51–56 is the total ITF transport on the Pacific side.

In addition to the transports across the transects, we calculate a barotropic streamfunction, ψ , defined as

where U and V are the horizontal velocity components integrated from the bottom to the sea surface, and the subscripts indicate partial derivatives. Here, we use the Cartesian coordinates for simplicity, but we use the spherical coordinates for the actual calculation. This definition implies that when ψ x > 0 and ψ y = 0 , for example, the barotropic current is northward: V > 0 and U = 0 . When ψ x = 0 and ψ y > 0 , the current is westward. In general, the barotropic circulation is clockwise around a local maximum of the streamfunction.

This definition implies that U x + V y = ( − ψ y ) x + ( ψ x ) y = 0 . Therefore, strictly speaking, the barotropic streamfunction does not exist because of small divergence or convergence. As described above, however, divergence and convergence within the Indonesian Seas are negligible on an annual average. It should also be small everywhere because barotropic adjustment should be fast.

For convenience, then, we use a “pseudo-streamfunction”: we first plug in zeros to the velocity variable at land grid points and then integrate V in the zonal direction starting from a point in Australia, which determines the values of the streamfunction along this latitude circle. From each gridpoint on the circle, U is integrated meridionally, determining the values everywhere on the sphere. Since we give zero velocity to land, streamfunction values are formally defined on land. As will be seen later, this arrangement is convenient to visualize transport between two landmasses. In particular, a western boundary current leads to a rapid zonal change in ψ , resulting in a maximum or minimum value right at the coast. Since the integration of ψ is continued inland, the landmass acquires the same value as along its coast. If the landmass were masked out, the maximum or minimum value at the coast would be hard to see on the map of ψ .

Figure 4 shows the vertical profiles of horizontally integrated, cross-sectional anomalous velocity, V ˜ ( z ; t , m ) , for transects 4 and 54 (Section 2.2.4); the transport values quoted earlier are the vertical integrations of this quantity.

In both transects, the profiles are very noisy above 200 and 300 m on the Indian and Pacific sides. The high-wavenumber structure is particularly prominent in transect 54. We have examined several instances of the meridional–depth section of annual-mean zonal velocity along this transect (not shown) and found that the noisy feature above 300 m is located near the northern end of the transect, that is, near the coast of Mindanao ( Figure 1 ). It is interesting that such a high-wavenumber feature remains on an annual average.

In contrast, no particularly noisy features are visible in the meridional–depth plots we examined (not shown) along transect 4; the wiggliness in the upper ocean is due to the horizontal integration of features with various vertical scales and structures. This difference could be because transect 4 is quite far from the actual narrow channels (Ombai and Timor).

Below this surface layer, the vertical structure appears more systematic, especially in transect 4. We here calculate EOFs (Section 2.2.3) of V ˜ ( z ; t , m ) .

Figure 5 shows the gravest five EOF modes for transects 4 and 54. The noisy near-surface layer is excluded from the calculation. For transect 54, the top 320 m is excluded. When the threshold was 300 m (not shown), one of the gravest modes peaked at 300 m, likely because it still caught the variability of the near-surface layer. The lower bound is set at 2,200 m, below which the EOF modes tend to vanish for both transects.

The EOF modes tend to have smooth vertical profiles. EOF1 does not change signs and has a vertical transport, whereas the other modes are baroclinic, reminiscent of the vertical dynamical modes associated with baroclinic waves. For transect 4, EOF1 is clearly dominant, whereas for transect 54, power goes down more slowly with increasing mode number.

Potemra et al. (2003) obtained EOF modes for the vertical profiles of the volume transports across Ombai Strait and Timor Passage in a 20-year integration of an OGCM. One of their gravest two modes (their Figure 6) is surface-intensified and nearly zero below 400 m. The other mode peaks at 400 m and vanishes below 800 m in Ombai or peaks at 200 m and vanishes below 600 m in Timor, with a weaker reverse flow in the top 100 m. The time series (principal components) of both modes include large interannual variability (their Figures 6e, f). These modes stay essentially the same after the removal of annual and semiannual harmonics. Our EOFs for transect 4 ( Figure 5 , left panel) miss these modes, all of which have multiple peaks below 400 m, either because we omit the top 200 m or because our variability is of different nature. It would not be surprising if the latter is the case because our data is a deviation from the ensemble average, not a temporal variability originating directly from external forcing.

The upper panels of Figure 6 compare V ˜ ( z , t ; m ) with the vertical profiles, [ ϕ n ( t , m )   e n ( z ) ] of EOF modes, for t = 1986 , for m = 1 and m = 3 , and for transect 4. In both members, EOF1 dominates while the second contribution, which is below 10%, comes from different modes. For all the top 5 years with large transport differences (1986, 1992, 1993, 1996, and 2013; see Figure 2B ), EOF1 is the dominant mode with a 70% contribution or larger except for t = 1993 and m = 3 , where EOF1 has only 42% and EOF1 + EOF2 has 76% ( Supplementary Figure S4 ). That is, EOF1 has a large amplitude and explains a large part of the vertical profile when the ITF ensemble anomaly is large, at least for the 5 years we have examined.

On the Pacific side, however, EOFs 1, 3, 4, and 9 have the largest contributions (see Supplementary Figure S4 for an example where EOF4 is the largest) at one of the top 5 years for the maximum and minimum members, and even the sum of the top 2 modes sometimes explains only approximately 35% of the variance (not shown). EOF1 is baroclinic, without much horizontal transport. There are large intrinsic variabilities at the entrance to the Indonesian Seas that are not directly related to large ITF transport anomalies. The lack of particularly dominant modes on the Pacific side might be because of the direct impacts of the western boundary current or other coastal currents. It is interesting that even an annual mean and meridionally integrated velocity (blue curves in Figure 6 ) includes extremely high-wavenumber features.

While EOFs describe intrinsic variability itself, they do not necessarily describe the intrinsic variability associated with the intrinsic variability of the ITF transport as in the case of transect 54. For transect 4, EOF1 does seem to be correlated with the ITF transport anomaly, but we looked at only the top 5 years of large ensemble spread. Here, we calculate the linear regression coefficient which extracts the component of the intrinsic variability which is linearly correlated with the total transport of the ITF, q ˜ ( t , m ) .

Figure 7 shows the part of V ˜ ( z ; t , m ) that is linearly correlated with q ˜ ( t , m ) for t = 1986 and m = 1 and m = 3 . See Section 2.2.2 for the method. The linear model approximately matches the actual profiles, particularly well below 600 m depth on the Indian side, whereas the actual flow is much more complicated on the Pacific side than the linear model, again suggesting that variability not related to the ITF variability contributes much more at the entrance into the Indonesian Seas than after the exit. The regression coefficient, α ( z ) , has a similar vertical structure to EOF 1 ( z ) on the Indian side with a double peak at 300 and 700 m. This result strongly suggests that the dominant variability ( EOF 1 ( z ) ) is the one associated with the variability of the total ITF transport. It is interesting to note that on the Pacific side, the component of V ˜ ( z ; t , m ) that is correlated with the total ITF transport has a smooth vertical profile below the near-surface noisy layer.

Also, it is interesting that some of the peaks in α ( z ) in the near-surface noisy layer lie above the significance threshold (lower panels of Figure 7 ). The anomalous velocity itself (purple curves) includes very high-wavenumber features down to ~1,000 m. The depths of some of the peaks in the purple curve agree with those of the peaks in the regression (green curve) in the upper 300 m. It would be surprising if this noisy feature, also in Figure 4 , is really correlated with the ITF transport. If it is real, it might be due to the western boundary current or other coastal currents acting on the coast of Mindanao. This issue may be an interesting subject of future study.

Figure 9 shows the regression coefficient between q ˜ and ψ ˜ this time using ψ ˜ mapped onto a 0.5° × 0.5° grid to save computational efforts. Otherwise, the coefficient is the same as that plotted in Figure 8C . Again, we multiply the coefficient α ( x , y ) by the value of the ITF transport at t = 1986 and m = 1 to make the physical interpretation of the coefficient more convenient because with this multiplication, we can describe how the 0.5 Sv increase in the ITF transport above the ensemble average flows around the global ocean as below. If the ITF transport is 0.5 Sv below the average as in m = 3 of t = 1986 , the signs should simply be flipped. If the amplitude of the increase or decrease is smaller, the streamfunction values below should simply be reduced by the same factor.

The ITF transport anomaly has a robust correlation with streamfunction values on the continents, which is not surprising at all as streamfunction values on landmasses represent transports between the landmasses (Section 2.2.5). For this reason, estimates of statistical significance are not shown in this plot. Correlation is perfect on the continents and neighboring, relatively quiescent regions and is weak in other parts of the ocean only because of background noise. Moreover, the regression coefficient is always zero on Australia by definition and is very small in the region surrounding Australia, either (1) because the barotropic velocity is weak (and therefore the streamfunction values do not change much from the zero value of Australia) or (2) because of the background noise. Since small regression values are formally classified as nonsignificant, it is difficult to distinguish these two cases. For this reason, the statistical significance of simultaneous correlation makes sense only for values on landmasses.

The Eurasian continent has a value of about 0.5 Sv, which is equal to the ITF transport anomaly for t = 1986 and m = 1 ( Figure 2B ), which manifests itself as a westward transport between Australia and Eurasia. This transport, moreover, takes the form of a narrow flow along 10° S; that is, it takes the form of an increase in the South Equatorial Current in the Indian Ocean. This narrow flow is represented by the rapid northward increase in streamfunction values (equation 6), which agrees with the linear barotropic response (McCreary et al., 2007).

It is very interesting that the streamfunction value is negative (about 0.4 Sv) on Antarctica, which means that the transport anomaly is westward between Australia and Antarctica. Even though the correlation of q ˜ with the Antarctic value of ψ ˜ is well above the significance threshold (not shown), no clear pathway is visible in the Southern Ocean as the streamfunction values vary wildly at small spatial scales and the regression coefficient is under the 99% significance level. To obtain a clear mean pathway, we would need many more samples (ensemble members) to cancel out the noise.

In any case, combined with the 0.5 Sv increase in ITF transport, there is a westward transport anomaly of 0.9 Sv in the Indian sector, which flows westward south of Africa into the Atlantic sector of the Southern Ocean and then through the Drake Passage into the Pacific.

In the Pacific, the flow pattern is not clear because of the noisy streamfunction field, but the following description should be correct to satisfy the mass conservation. Most of this anomalous transport appears to bend northwestward to the west of South America and then flows back southwestward, circling the purple patch at ∼100° W, 45° S. It then continues westward again. Upon hitting New Zealand, 0.5 Sv bends northward along the Tonga–Kermadec arc and then bends northwestward to join the ITF; the rest flows south and joins the westward transport south of Australia.

McCreary et al. (2007) used an OGCM in a global domain without the Arctic Ocean and obtained a standard climatological solution. In a sensitivity solution, the Indonesian passages are closed at 8.5° S. The difference in the barotropic streamfunction between the two solutions shows a pure ITF. It flows westward along 10° S, and at Africa, it bends southward along the coast of Africa. Most of this flow bends eastward at the southern tip of Africa, forming a narrow eastward flow, and then bends northward at the southern tip of Tasmania, flowing along the east coast of Australia. Some of the transport, however, flows westward from the southern tip of Africa and eventually crosses the Drake Passage into the Pacific and eventually joins the northward flow along the east coast of Australia. This flow pattern is broadly similar to the ensemble anomaly of the streamfunction described above, except that the narrow eastward return flow that reaches Tasmania does not exist in our ensemble anomaly. The reason for the westward flow in the Southern Ocean is not clear, either in our result or in the result of McCreary et al. (2007). This potential connection of the net Southern-Ocean transport all around Antarctica with the ITF transport might be through some changes in overturning circulation similar to the change explored by Sen Gupta et al. (2016, their Supplementary Figure S10).

Lead-lag regression coefficients (Section 2.2.2) are also calculated ( Figure 10 ). The figure shows only the North Pacific since there is no systematic and significant signal in other regions ( Supplementary Figure S5 ) perhaps except for the potentially weaker stripes in the mid-latitude southern hemisphere. At lag 5 ( Δ t = − 5 , ψ ˜ leading q ˜ ), there is a negative band extending southeastward from about 170° E, 30° N, and other positive and negative stripes south and southwest of it. There is another weaker negative band extending eastward from the south of Taiwan. As time goes on, this negative anomaly moves, together with the positive anomaly on its southern flank, westward and equatorward. Other lesser stripes appear to move similarly, but their motions are not as clear. At lag 1 ( Δ t = − 1 ), the strongest negative band reaches Okinawa and the associated positive band reaches Taiwan. At the same time, there is a negative band extending eastward from the southern tip of the Philippines. At lag zero ( Figure 9 ), the major positive and negative bands are still visible, and the other lesser stripes are buried in the large-scale signal described above. At lag −1 ( Δ t = + 1 , ψ ˜ lagging q ˜ , Supplementary Figure S6 ), the negative band is still visible but much weaker extending from Taiwan, and there are no significant and systematic signals elsewhere, either.

The same analysis is carried out for sea-surface height (SSH). The simultaneous linear regression (not shown) is noisier than the corresponding linear regression on ψ ˜ and provides similar information: SSH tends to be lower north of 10° S in the Indian Ocean, consistent with the westward barotropic transport there; SSH is higher near Antarctica, consistent with the westward transport in the Antarctic Circumpolar Current; and west of South America, there is a region of higher SSH coincident with the lower ψ ˜ region there, consistent with the counterclockwise barotropic circulation there. On the other hand, SSH does not show the northward barotropic flow along the Tonga–Kermadec ridge ( Figure 9 ). This difference may be an indication that the barotropic flow there is deeper, with less SSH signal, or that SSH “noise” is larger there.

Figure 11 shows a set of lead-lag regressions of the SSH anomaly on q ˜ corresponding to that of ψ ˜ . The spatial pattern of the anomaly is very similar to that of ψ ˜ . The strong negative and positive bands of SSH that correspond to those of ψ ˜ are particularly clear, and they migrate westward and equatorward. Other, lesser, stripes equatorward of this pair are also sometimes visible. The potential relation of these stripes in streamfunction and SSH to the ITF transport will be discussed below.

The stripes are indeed the most dominant feature in the anomalous barotropic streamfunction itself (not shown) and in its EOFs (not shown either), except for the obvious strong variability in the western boundary currents like the Kuroshio and the Gulf Stream and in the Antarctic Circumpolar Current.

These stripes are a manifestation of “zonal jets” (see the review in the Introduction of Furue et al., 2021 and references cited there). The ocean basins are populated with quasi-barotropic alternating zonal flows. The jets are tilted in the northwest–southeast direction and migrate equatorward and westward ( Figures 10 , 11 ). The equatorial migration has also been found in the previous version of OFES (Richards et al., 2006) and in the same OFES2 dataset (Furue et al., 2021) as we use here in the present study. Both the tilt and equatorial migration are consistent with the general westward deepening of the sea floor in the Pacific or the westward thickening of the main pycnocline due to the Subtropical Gyre according to the idealized theoretical studies by Boland et al. (2012) and Khatri and Berloff (2018, 2019).

Does the lagged correlation ( Figures 10 , 11 ) mean that these zonal jets force the ensemble anomaly of ITF transport? If so, this variability of ITF is random (that is, not directly determined by external forcing) because the phases of the zonal jets are random (except near the equator) despite the regularity of their propagation (Furue et al., 2021). The stripes in the ensemble anomaly of barotropic streamfunction appear blocked by prominent sea-floor ridges such as Tonga–Kermadec arc and Ninety East Ridge [not shown, but the blockage is visible in the linear regression ( Figure 9 ), clearly for Tonga–Kermadec and barely for Ninety East Ridge], which suggests that the zonal jets are affected by bottom topography. There could then be systematic JEBAR stress that could alter the ITF transport similarly to the way that wind stress in the Pacific determines the overall ITF transport (Godfrey, 1989). According to (Godfrey’s, 1989) island rule, wind stress along latitudes just north of Halmahera and along latitudes just south of the southern tip of Australia (if we ignore the blockage by New Zealand) is key. When particularly large jets arrive north of Halmahera, the ITF transport would increase or decrease, and this anomaly would run across the three ocean basins in the southern hemisphere at a speed of barotropic adjustment ( Figure 9 ). The weakness of the lagged correlation with the southern hemisphere field (Section 3.5.2; Supplementary Figure S5 ) might be explained by the blockage of the stripes by the Tonga–Kermadec ridge (anonymous reviewer, private communication).

This hypothesis could also explain the apparent weak ∼4-year peak in the autocorrelation ( Supplementary Figure S2 ) through the potential periodicity in the migration of zonal jets. For example, a weak yellow (positive) anomaly is located in southern Japan, and another is located in Taiwan at lag 5 ( Figure 10 ). At lag 4, stronger positive anomalies extend eastward or southeastward from southern Japan and from Taiwan. At lag 1, similar positive anomalies are located at similar latitudes. This suggests a potential periodicity of 3–4 years.

The strong negative simultaneous correlation extending eastward or southeastward from Taiwan in Figure 9 is puzzling. One hypothesis to explain this is that the entire set of positive and negative stripes between Taiwan and Halmahera, which systematically migrates toward the equator, is actually correlated with the ITF transport anomaly. Another hypothesis is that this negative anomaly in Taiwan is related to the magnitude of the transport through the South China Sea to the Indian Ocean. This is an interesting subject for future investigation.

Calil (2023) showed that zonal jets at depth are much more robust in the 1/30°-resolution version of their Atlantic model than in their 1/10°-resolution version. Furue et al. (2021) also found that OFES2’s zonal jets have roughly half the amplitude of the observation by Cravatte et al. (2017) in the central-to-eastern tropical Pacific. If these findings apply to the western tropical and subtropical Pacific and if the zonal jets really have something to do with the intrinsic variability of the ITF transport, the real amplitude of the intrinsic variability might be larger than found in the present study.

In any case, this is speculation, and the evidence we have provided is hardly sufficient. We just propose this hypothesis for future studies to test.

Sasaki et al. (2018) found an increase in the ITF transport when a tidal-mixing parameterization is introduced to OFES2. Considering the size of the ensemble spread found in the present study, one might wonder whether the increase they found may have been just an ensemble spread triggered by the introduction of tidal mixing. The increase, however, is a long-term trend from 1960 to 2014 (their Figure 1) and therefore is not likely a stochastic (random) variability, as the latter is not persistent in time ( Figure 2B ). The interannual variability they found, on the other hand, should include the variability we have found in the present study.

Intrinsic variability should depend on resolution, model parameters, forcing, and boundary conditions. Since chaoticness depends on nonlinearity, coarser models are expected to be generally less chaotic. Damping of variability would naturally weaken intrinsic variability, and therefore surface heat flux would dampen it and mixing would, too. The so-called "mixed boundary conditions" are, however, known to produce chaotic oscillation in coarse-resolution models with large vertical diffusion (Huang and Chou, 1994).

Finally, the present study has explored mainly barotropic circulation. The vertical structure, however, has a number of curious features, as shown in Sections 3.3.2 and 3.3.3. Since the speed of horizontal propagation of baroclinic waves strongly depends on their vertical structure, much more sophisticated analyses would be necessary than those employed in the present study. The vertical structure of the ITF’s intrinsic variability will be an interesting subject of future studies.