The island rule with multiple islands and its application to the Indonesian Throughflow

The island rule theory in the case of complex geometry with multiple islands referring to the Indo-Pacific Maritime Continent is investigated on the basis of Godfrey’s island rule theory. The bottom friction and lateral friction of multiple channels are considered by employing the Munk and Stommel boundary layer models. Five idealized cases with various spatial distributions of islands are designed to examine the influence of shape and size of the islands. The analytical solutions of the streamfunctions of the through-flows among the islands are obtained and the volume transport through each channel is estimated. The Indonesian Throughflow (ITF) transport is then calculated using the analytical solutions with wind stress and compared with observations and previous theoretical results. The ITF transport from the multiple island rule is about 14.5 Sv during 2004–2006, which is close to the observed ITF transport (about 15.0 Sv) from the International Nusantara Stratification and Transport. We find that the multiple islands rule reproduces well the mean value and interannual variability of the observed ITF transport, and inclusion of wind stress in the North Pacific Ocean may improve the estimate of ITF transport. Sensitivity experiments indicate that frictional boundary layer thickness and channel size influence the estimated ITF transport under the multiple island framework. These results imply that the multiple island rule shows improvements in estimating the ITF transport relative to previous studies, and the multiple island rule can be used to produce long time series of ITF transport and might have implications for paleo-ITF study.


Introduction
The Indonesian Seas are the only channel that connects ocean basins in the tropics through the Indonesian Throughflow (ITF), which transfers tropical Pacific waters to the Indian Ocean (Gordon, 1986;Hu et al., 2015;Sprintall et al., 2019). The ITF acts to shape the distribution of Indo-Pacific ocean heat and fresh water, links to the global thermohaline circulation, and influences the Indo-Pacific and global climate system (Lee et al., 2002;Lee et al., 2015;Feng et al., 2018;Forget and Ferreira, 2019;Hu et al., 2019;Pujiana et al., 2019;Santoso et al., 2022).
The ITF and its variability are thought to be associated with large scale wind forcing that asserts on the Indo-Pacific Ocean with complex channels (e.g., Wyrtki, 1987;Wainwright et al., 2008). Godfrey (1989) (hereafter G1989) proposed the island rulethe theory derived from the Sverdrup model. This theory implies that, under idealized conditions, through the integration of the simplified equations of motion of the fluid, the flow of seawater around an island can be evaluated from wind stress (Godfrey, 1989). G1989 estimated the ITF transport using the island rule and concluded that the ITF transport was 16 ± 4 Sv. As the theory builds up the relationship between surface wind forcing, friction, geometry and ocean current, the theory is simple enough to be understood and one can understand how the ocean current response to wind forcing with this theory. Hence, the island rule theory has unique advantages comparing with other methodlogy like numerical modelling.
The Godfrey's island rule theory has been applied in many previous studies. For example, Firing et al. (1999) derived a timedependent island rule and investigated the North Hawaiian Ridge Current north of Oahu. Meng et al. (2004) used the island rule to examine the relationship between Pacific wind stress and ITF volume transport on an inter-decadal scale and found that the integral of zonal wind stress along the equator determined interdecadal change of the ITF. Cai (2006) utilized the Godfrey's island rule to investigate the linear trend of the Southern Ocean supergyre circulation. Liu et al. (2007) assessed the ITF transport and the Luzon Strait transport with the island rule and suggested that the westerly component of the equatorial Pacific winds accounted for the decrease of ITF transport after 1976. Joyce and Proshutinsky (2007) evaluated the streamfunction of flow around the Greenland by using the island rule. Feng et al. (2011) and Feng et al. (2018) suggested that the island rule theory is able to capture the decadal and multi-decadal variations of the ITF to some extent.
Although the Godfrey's island rule theory provides an invaluable approach for understanding observed features and changes in large scale ocean circulation like the ITF, further improvement of this theory is expected considering its extreme simplification and the nature of the complex geometry of the real oceans. Wajsowicz (1993) revised the Godfrey's island rule through including bottom topography and frictional effects along eastern boundaries, and for the first time extended the island rule to two islands. Pedlosky et al. (1997) re-derived Godfrey's island rule in a general form and discussed the role of dissipative boundary layers and inertial effects in estimating the net transport around an island. Pratt and Pedlosky (1998) took into account the dissipation on the northern, southern, and eastern boundaries of an island, and showed that lateral friction is a crucial factor in overestimating the ITF transport with Godfrey's island rule. Lian et al. (2017) proposed a parameterized scheme of friction-topography resistance on the basis of Wajsowicz (1993) and investigated the throughflow in the South China Sea. Wang et al. (2018) suggested that the island rule theory might be further improved through using an optimized path integral and friction parameterization, considering a more complex geometry and bathymetry, and/or adding the timedependent term.
Recently, Yang et al. (2019) and Yang et al. (2020) (hereafter Y2020) examined the streamfunction of each island under various geometries according to the different meridian lengths of the two islands, following Wajsowicz (1993) but considering both lateral and bottom friction. They found that the ITF transport, when considering both kinds of friction, was about 8.7 Sv, which was reduced by about 16.6% compared with that without considering friction . Although Y2020 considered the frictional effect, the ITF transport based on their new theory is significantly lower than the observed ITF transport (Sprintall et al., 2009).
Two key points can be drawn from previous research: 1) Importance of more realistic geometry in the island rule. The ITF consists of flows through several channels and/or straits, and considering multiple islands leads to including more channels and/or straits and hence may make the path of integration be more consistent with the real geometry. 2) Importance of including of wind forcing from the northern hemisphere in estimating the ITF with the island rule. Recent studies suggest that the north Pacific forcing is very important in determining the variability of ITF (e.g., Li et al., 2020), indicating that including the north Pacific wind forcing might be able to improve the estimate of ITF with island rule. Hence, in this study, we aim to further explore the island rule in the case of more complex geometry with multiple islands referring to the Indo-Pacific Maritime Continent on the basis of previous studies especially Y2020.
In the following, we will derive the analytical solutions of the streamfunctions of the through-flows among islands, and examine the ITF transport with the analytical solutions. The momentum equations and streamfunctions will be described in Section 2, followed by Section 3 presenting the estimated ITF transports and dynamic mechanism, and Section 4 examining the sensitivity of results to the frictional boundary layer thickness and channel size. We will discuss remaining issues and summarize the main results in Section 5.

Momentum equations and streamfunctions
In order to investigate the island rule in the case of complex geometry, we consider five idealized cases with different spatial distributions of three islands, including one of which referring to the Indo-Pacific Maritime Continent. It is well known that there are many smaller islands within the Indonesian region and the ITF is composed of several throughflows, but for simplification, three islands [i.e., Kalimantan Island, Philippine Islands and Australia-Papua New Guinea (PNG)] are considered in this study. Previous studies indicate that the Makassar Strait Throughflow, the major component of ITF, has a transport occupies about 78% of the total ITF . So, the current across the straits between the three islands selected is the most component of ITF, and ignoring other small islands is expected to have neglecting contribution to the total ITF transport. Including more islands in this theory may further increase the accuracy, but also makes the analytical solutions be more complex.

Momentum equations
For simplification, we consider barotropic and steady flow with the time term and nonlinear term ignored, and the horizontal momentum equations is a simplified two-dimensional Navier-Stokes equations. In view of friction's importance, we consider the lateral and bottom frictions between islands and the friction is determined by the width of channel and the frictional boundary layer's thickness (see details in Section 2.3). Integrating the equations of motion from the bottom depth (i.e., z=-Z) to the sea surface (i.e., z=0) in the vertical direction with the Beta-plane approximation, the momentum equations are: udz represent the velocity component of the depth integral, f denotes the Coriolis parameter, r 0 is the fluid density, P is the pressure term of the vertical integration of the fluid, (F (x) ,F (y) ) represents the friction term of the vertical integration of the fluid, and (t (x) ,t (y) ) denotes the sea surface wind stress. Since the island rule calculates the volume transport in horizontal, the Ekman pumping caused by wind stress curl is not considered.

Distribution of islands and streamfunctions
The topography of Indo-Pacific Maritime Continent (from ETOPO1) and horizontal currents (vertically averaged over the upper 100 m layer) from the Simple Ocean Data Assimilation (SODA) are shown in Figure 1. In order to test the response of ITF to change of geometry (e.g., the latitudes of islands, spatial patterns, width and length of straits, etc) in Indonesian Seas under the framework of island rule, we simplify the real island geometries and design five idealized cases Case 1-Case 5. These cases are designed considering the nessary of univariate experiments with geometry changes and potential implications for paleo-ITF that was controlled by paleo-geometry. Case1 is directly simplified from the real geometry of the Indo-Pacific region, of which the western island represents the Kalimantan Island, the middle island represents the Philippine Islands, the eastern island represents the Australia-PNG, and the eastern boundary of the ocean is the western boundary of North and South America. Then, we design the other four cases Case 2-Case 5 by changing the spatial location or meridional extent of islands. From Case1 to Case5, only the latitude or size of one island is changed each time, which will cause the change of integral path. Below we present the streamfunctions in the five idealized cases.
Case 1: The geometry of Case 1 is shown in Figure 2. The western island partially overlaps with the middle and eastern islands in latitude, while the middle island and eastern island do not overlap with each other in latitude ( Figure 2). In the original Godfrey's island rule, the viscous effects of the interior ocean and the eastern boundary is ignored. But different from the single island, the frictional effects between the island channels for the multiple islands are considered in the present study. We integrate the wind stress and friction terms in Eq. (1), (2) along the closed curves C 0 , C 1 , and C 2 , and get: Topography of the Indo-Pacific Maritime Continent (shaded, unit in m) and mean horizontal currents averaged over the upper 100 m layer from SODA (black streamlines). Blue contour lines denote the 1000 m isobaths, and red contour lines denote the 5000 m isobaths. Then, Where Df i ,i=1,2,3,4 denotes the difference of Coriolis parameters between the north and south latitudes, and dl represents the differential vector along the direction of the integral path. And the streamfunctions of the island without considering the friction effect are: As shown in Figure 2, the integral path around the island is all along the western coastline of the island (similarly hereinafter). The path of the integral curve C 0 is denoted by FBGHF, the path of the integral curve C 1 is signified by IJCKI, the path of the integral curve C 2 is represented by KCDLMABFK, and the curves AB and CD stand for the right-hand side (RHS) boundaries of the island channel.
Case 2: In this case, latitudes of the western island partially overlaps with that of the middle and eastern islands, and latitudes of the middle island partially overlaps with that of the eastern island ( Figure 3). By integrating Eqs. (1) and (2) along the closed curves C 0 , C 1 , and C 2 , we get: where, and represent streamfunctions when omitting the friction term, the path of the integral curve C 0 is specified by FBGHF, the path of the integral curve C 1 is represented by IJDCBFI, the path of the integral curve C 2 is denoted by IJEKLABFI, and the curves AB, CB, and DE stand for the RHS boundaries of the island channel. Map of islands and path of integration in Case 1. Shi et al. 10.3389/fmars.2023.1145506 Frontiers in Marine Science frontiersin.org Case 3: The western island contains the middle island in latitude and partially overlaps with the eastern island, and the middle island and the eastern island moderately overlap in latitude ( Figure 4). Through integrating Eqs. (1) and (2) along the closed curves C 0 , C 1 , and C 2 , we get: In this formula, and denote streamfunctions when not considering friction. Further, the paths of the integral curve C 0 , C 1 , and C 2 are represented by FBGHF, IEDCBFI, and JKLABFJ, respectively, and the curves AB, CB, and DE denote the RHS boundary of the island channel.
Case 4: The western island partially overlaps with the middle island in meridional, the middle island partially overlaps with the eastern island in latitude, and the western island and eastern island do not overlap ( Figure 5). By integrating Eqs. (1) and (2) along the closed curves C 0 , C 1 , and C 2 , we have: Where, and are streamfunctions when omitting the friction term. Additionally, the paths of the integral curves C 0 , C 1 , and C 2 in order are specified by FBGHF, IDJABFI, and KLMCDIK, and the curves AB and CD represent the RHS boundaries of the island channel.
Case 5: The setup of this case is the same as Case 1 but the width of M-A is reduced by 100 km. As in Figure 2 but for Case 2. Shi et al. 10.3389/fmars.2023.1145506 Frontiers in Marine Science frontiersin.org

Integral of lateral and bottom frictions
As shown in Section 2.2, the streamfunctions can be expressed as the sum of wind stress term and friction term. The wind stress term can be calculated from observed data of sea surface wind, but the friction force terms are not readily available, and the vorticity equations need to be solved firstly to establish the relationship between transport and friction term. Wajsowicz (1993) considered the situation of lateral and bottom frictions separately and modified the island rule. Yang et al. (2020) suggested that considering both the bottom and lateral friction may be of much importance. Hence in this study, we adopt the Munk-Stommel model (e.g., Stommel, 1948;Munk, 1950) following Yang et al. (2020), and the bottom and lateral frictions are substituted into Eqs. (1) and (2) at the same time: Then the vorticity relation is given by: where A S = D f Z and A H represent the coefficients of the bottom friction and the lateral friction respectively, and b = df dy denotes the gradient of Coriolis parameter. We take into account a combination As in Figure 2 but for Case 3. As in Figure 2 but for Case 4. Shi et al. 10.3389/fmars.2023.1145506 of friction coefficients and geometric/integral paths, although the bottom and lateral frictions used are simplified and idealized frictions and may have potential influences. This simplification may affect the exact value of ITF shipments, but does not change the main result in terms of ITF variability. Eq. (32) is a fourth-order partial differential equation, and it is difficult to obtain its analytical solution. Hence, the island channel is assumed to be a rectangular channel with width and length denoted by W and L ( Figure 6). For the case of W L ≪ 1, we can rationally state that ∂ ∂ x ≫ ∂ ∂ y , and thereby the vorticity equation is simplified to: We first nondimensionalize Eq. (33), which is particularly important for the subsequent discussion of the contribution of bottom friction and lateral friction to transport and helps to diagnose the differences between the model and previous ones. We assume that there is a Sverdrup balance in the channel, then the last term on the left and the first term on the right of Eq. (33) are balanced. If U is the characteristic scale of horizontal velocity, W is the characteristic scale of horizontal length, and t 0 is the characteristic scale of t , then there is U=t 0 r 0 −1 b −1 W −1 . If y is normalized by UW, we get the following formula: where, d M is the Munk boundary layer thickness and d s corresponds to the Stommel boundary layer thickness. The lateral friction coefficient is set as A H =10 4 m 2 s −1 (Wajsowicz, 1996), the bottom friction coefficient is set as A S = 2 3 Â 10 −5 s -1 (The water depth is taken as 300m) following Yang et al. (2020), and the gradient of Coriolis parameter is set as b=1.62×10 −11 m −1 s −1 . Therefore, the Munk boundary layer is d M =85 km, the Stommel boundary layer is d s = 412 km, and d M is less than d s . For simplification, friction is considered only when the strait width is less than the boundary layer thickness, and hence we choose different models according to the width of channel.
1) When 0<W≤d M <d s , both the bottom friction and lateral friction are considered and the model is similar to the Munk-Stommel model. Eq. (33) is then expressed as follows: with the following boundary conditions: By introducing boundary conditions, we get the general solution of Eq. (35): y (x, y) = l 0 + l 1 e r 1 x + l 2 e r 2 x + l 3 e r 3 x + d b x, where, r 1 , r 2 and r 3 represent the roots of the characteristic equation, l 0 , l 1 , l 2 and l 3 are constants. Please refer to the appendix for detailed derivation process.
Integrating the friction expression along the length and we have: ðhere; k = ( − 1) j i + jÞ Schematic map of the passage between the two nearby islands. Shi et al. 10.3389/fmars.2023.1145506 Frontiers in Marine Science frontiersin.org 2) When d M <W≤d s , the bottom friction is considered and the lateral friction coefficient is zero, and in this circumstances, the model is similar to the Stommel model. Eq. (33) is expressed as follows: with the following boundary conditions: y(0,y)=y 1 , y (W, y) = y 0 By introducing boundary conditions, we get the general solution to Eq. (39): Where, r 1 represents the root of the characteristic equation, l 0 and l 1 are constants. Please refer to the appendix for detailed derivation process.
Integrate the friction expression along the length; therefore, we have: In the formula, m = A S Lb a ( Wr 1 1−e −r 1 W − 1), n = − A S Lr 1 1−e −r 1 W . 3) When d M <d s <W, the bottom friction and lateral friction are not considered, so the integral of friction is zero.

Influence of geometry on channel transport
To examine the influence of geometry on the channel transport, we then estimate the transport y 0 −y 2 , i.e., the through flow in the channel between the western island and the eastern island, in the five cases. The monthly average wind speed data at 10 meters above the sea surface based on the CCMP (Cross-Calibrated Multi-Platform) are employed. The CCMP wind speed data spans from 2004 to 2006 with a spatial resolution of 0.25°×0.25°, and the sea surface wind stress is obtained by using the empirical formula. It should be noted that, previous studies have suggested that wind stress calculated by actual wind and relative wind might have an impact on ocean current simulations (e.g., Wu et al., 2012;Sun et al., 2021), but in order to contrast with the previous island rule, only actual wind is considered in this paper. The drag coefficient of wind stress has been investigated in numerous previous studies (e.g., Garratt, 1977;Powell et al., 2003). Here, we employ the empirical formula associated with the wind stress and drag coefficient following Yelland and Taylor (1996): where t x and t y represent the zonal and meridional wind stresses, r a =1.29 kg m −3 is the air density, C D denotes the drag coefficient, jṼj is the modulus of the wind vector, and u and v are the components of the wind speed in the east-west and north-south directions, respectively. We assumed that the width of the channel is constant in each case. The width (W 2 ) of the channel between the western island and the middle island is set as 100km, the width (W 1 ) of the channel between the western island and the eastern island is set as 200km (this assumption is valid only when there is a channel between the two islands), and the climatological wind stress curl between channels is adopted. Table 1 presents the latitude ranges of the islands and corresponding channel transports in the five cases. Note that the total transport in Case 1 equals to the channel transport through the section M-A, in Case 2 equals to the channel transport through the section L-A, in Case 3 equals to the channel transport through the section L-A, and in Case 4 equals to the channel transport through the section M-C plus J-A. As shown in Table 1, the channel transport decreases from 14.47 Sv in Case 1 to 8.70 Sv in Case 2. The difference of islands' setup between Case 1 and Case 2 is the latitude of the middle island, which leads to changes of channel friction and integral paths. The main changes of geometry in Case 3 from Case 2 are a reduction in length and southward shift of the middle island (Table 1 and Figure 4). As a result, the western island streamfunction in Case 3 would be no longer affected by the friction of the western boundary of the middle island, and the transport rises to 10.12 Sv. In Case 4, the western island is shifted to the north by 6°relative to Case 3 and does not coincide with the eastern island in latitude. Although the geometry changes in case 4 lead to an increase of friction, the streamfunction of the western island in Case 4 is affected by the middle island compared to Case 3, in which the transport is further increased to 12.75 Sv. The Case 5 is similar to Case 1. Compared with Case 1, the latitude range of islands in Case 5 remains unchanged, but the width of the channel between the middle and eastern islands is reduced by 100 km. As a result, the transport is decreased from 14.47 Sv in Case 1 to 9.02 Sv in Case 5, due to the narrowing down of the channel and a larger contribution of bottom friction. To summarize, the above comparison indicates that the channel's transport is very sensitive to geometry through the effect of friction and integral paths.

ITF transport estimate
The Indonesian Seas where the ITF passes by have a very complex geometry and bathymetry. Here we optimize the integral path for simplification. Figure 7 demonstrates the climatological sea surface winds in the Pacific Ocean and the integral path we used. Arrows represent the climatological winds from the CCMP product. The distribution of the three islands and the streamfunctions are represented by the idealized Case 1 discussed in Section 2.2. The Makassar Strait is the main inflow channel for the ITF (Wajsowicz, 1996;Susanto and Gordon, 2005;Vranes and Gordon, 2005), accounting for 75%-80% of the ITF water transport (Gordon et al., 2008;Li et al., 2018;Gordon et al., 2019). Case 1 is an idealized geometry, and it needs to be adjusted according to the actual geometry in the estimation of ITF. The channel's width W 1 shown in Figure 7 includes the width of the Makassar Strait and the Lifamatola Strait. Although the Lifamatola Strait is also an entrance to the ITF, its contribution to the ITF is small compared to the overall transport and its width is negligible compared to the Makassar Strait. Therefore, we only take into account the width of the Makassar Strait. Suppose that the width and length of the Makassar Strait are represented by W 1 and L 1 , and the width and length of the Mindoro Strait are denoted by W 2 and L 2 , respectively. Then the transports through the Mindoro Strait and Makassar Strait specified by T 1 and T 2 are: T 2 = y 0 − y 2 = y 0,wind − y 2,wind + R 3 m 1 + R 4 m 2 + R 3 n 1 (y 0 − y 2 ) + R 4 n 2 (y 1 − y 2 ): (44) Df 1 Df 2 : m 1 and n 1 are formally the same as the expression m and n in Section 2.3, except that L 1 and W 1 replace L and W in the formula. Same thing with m 2 and n 2 .
By simultaneous solving the relations in Eq. (43) and Eq. (44), the transport of the Makassar Strait is obtained as: T 2 = R 4 n 2 (y 1,wind − y 2,wind + R 1 m 1 + R 2 m 2 ) + (1 − R 2 n 2 ) (y 0,wind − y 2,wind + R 3 m 1 + R 4 m 2 ) 1 − R 2 n 2 − R 3 n 1 + R 2 R 3 n 1 n 2 − R 1 R 4 n 1 n 2 : (45) To calculate the ITF transports over a long period of time, we extend the time span of the monthly sea surface wind product of CCMP to 1988-2016. The geometry data of the Makassar Strait and Mindoro Strait are given by: W 1 =200km, L 1 =1200km, W 2 =100km, L 2 =700km, and the average wind stress curl between channels is set as rot( t * ) = −5:42 Â 10 −8 N m −3 . Because d M <W 1 ≤d s , the width .70 Sv (Section: L-A) Case 5 As Case 1 but width of M-A is reduced by 100 km. 9.02 Sv (Section: M-A) of the channel is much larger than the Munk boundary layer thickness, and the lateral friction is not included in the calculation, so the Stommel model is adopted. The time series of ITF from 1988 to 2016 can be obtained by inserting parameters and wind stress data into the Eq (45). Figure 8A illustrates the yearly mean transport of ITF from the INSTANT observations and that calculated with island rule during 2004-2006. The INSTANT program deployed 11 moorings in the major inflow and outflow passages of the ITF during 2004-2006, and the total transport through three passages (Lombok Strait, Ombai Strait, and Timor Passage) during the INSTANT period was 15.0 Sv (Sprintall et al., 2009). As shown in Figure 8A, the ITF transport predicted by G1989 and Y2020 are smaller than the observed value, and the ITF transport derived from the multiple islands rule of this study has the smallest difference from the observations. Figure 8B Figure 8B, the ITF transport from the triple island rule theory has a significant correlation coefficient of about 0.9 with the INSTANT ITF transport, while the correlation coefficient between original island rule transport and the INSTANT ITF transport is only 0.25, implying that the multiple islands rule of this study is able to infer a more reasonable time series and shows significant improvement in estimating the ITF transport comparing with the original island rule theory.
To further verify the multiple islands rule in a longer time period, we plot in Figure 9 the yearly transports during 1988-2016 from G1989, Y2020 and the present study, with a 6-month running mean of the data in advance. The annual and monthly average ITF transport by INSTANT during 2004INSTANT during -2006   Integral paths (red, blue and black lines) and climatological sea surface winds (vectors). Color shading shows the topography (unit in m). Shi et al. 10.3389/fmars.2023.1145506 Frontiers in Marine Science frontiersin.org ITF transport calculated using the MITF velocity are also added to the picture. Note that the monthly average has been subjected to a 13-month running mean to remove seasonal signals. The MITF program deployed two moorings in the Makassar Strait, on which two acoustic Doppler current profilers (ADCPs) measured the Makassar throughflow from surface to the 680 m sill depth (Gordon et al., 2003;Gordon et al., 2019). The time span of MITF velocity dataset is from November 2006 to December 2016 (Only eight months of data in 2017 and hence are not used). By splicing it with the INSTANT observation dataset, the time span is expanded to January 2004 -December 2016, and it was preprocessed into a time interval of one day. Since MITF program only conducted mooring observation in the Makassar Strait and not for the entire ITF, we use the INSTANT time series to correct the MITF transport through a formula to get a MITF-based transport of the total ITF. First of all, we averaged the single point velocity measured by ADCP in the vertical direction to obtain the average velocity of each day, and then averaged to get the average speed of the whole time period. The ratio of the two was multiplied by the average flow of INSTANT to obtain the MITF transport, and the corresponding expression is MITF ¼V vertical =V time Â 15, V vertical represents the vertical average of the horizonal velocity, the number 15 represents the average ITF transport during the INSTANT period, and V time represents the average over the time dimension of V vertical .

and the annual mean
As shown in Figure 9, the average transports from G1989, Y2020 and the present study are 13.03 ± 3.30 Sv, 9.94 ± 2.52 Sv, and13.00 ± 3.94 Sv during 1988-2016, respectively. Observations show that the ITF transport during 2004-2015 was about 15.64 Sv (Li et al., 2018), and for comparison, the average transports during 2004-2015 from G1989, Y2020 and the present study are 13.60 ± 0.87 Sv, 10.38 ± 0.51 Sv and 14.10 ± 1.02 Sv, respectively. Hence, in general, the multiple islands rule of this study predicts a reasonable mean transport that is consistent with observations. Both the G1989 and the multiple islands rule predict a reasonable mean transport of about 13.03 and 13.00 Sv, but the G1989 transport show significant differences in interannual peaks and valleys comparing with observations. For example, geostrophic calculations indicate that the ITF in about 1988-1989 was smaller than that of the 1990s and 2000s , but in contrast, the transport from G1989 and Y2020's island rule is greater in 1988 than that in 1990s and 2000s (Figure 9).
We then calculated the average errors of the three (G1989, Y2020, and this study) time series and their correlation coefficients with the MITF time series. Here, the average error is defined as the average of the difference between a time series and the MITF time series. The average errors of the former two (G1989 and Y2020) are 1.32 Sv and 4.48 Sv, respectively, which are much greater than that of the present study that is only 0.87 Sv. The correlation coefficient between the G1989 (Y2020) time series and the MITF is both 0.37 (0.37), while the time series from the multiple island rule and the MITF has a correlation coefficient of 0.81. Therefore, the comparison implies that the multiple islands rule might be improved from previous studies of island rule theory. It should be noted that, since the time series from the island rule shown in Figure 9 were calculated using the same wind stress data set, the significant differences between them should be attributed to from the optimizing of the fraction and path of integration related to the geometry.
As we mentioned, Pratt and Pedlosky (1998) examined the role of lateral friction in the island rule theory, and pointed out that lateral friction is an important factor in the overestimation of ITF in the original island rule. But it seems that the major straits (e.g., Makassar Strait) of the ITF passing by have a width greater than the thickness of Munk boundary layer except a few narrow passages, and the lateral friction from these narrow passages may have limited influence on the overall transport. As an alternative, the interchannel bottom friction, which can significantly reduce the estimate of ITF, may play an important role. As we tested, the ITF transport can reach 27.6 Sv, which is obviously unreasonable, if the inter- Yearly ITF transport during 1988-2016 derived from G1989 (red), Y2020 (blue) and the solution of this study (black). Yearly and monthly INSTANT observations during 2004-2006 (cyan and magenta respectively) and yearly MITF transport during 2004-2016 (gray) are also shown for comparison. channel bottom friction dissipation in the Indonesian sea is not taken into account.

The influence of including North Pacific winds
The island rule reflects the relationship between the ITF transport and the Pacific wind field. Previous studies suggest that the decadal variability of ITF is related to the North Pacific trade winds and the ITF is primarily drawn from North Pacific thermocline waters (e.g., Li et al., 2018;Li et al., 2020). The North Pacific trade winds has significant impact on the Northern Equatorial Current (NEC) bifurcation latitude and Mindanao Current (MC) and hence the ITF (e.g., Hu et al., 2020;Hu et al., 2021). One of changes in the present multiple islands rule from previous theories is the inclusion of wind stress of northern Pacific Ocean in the path of integration, and this change might play an important role in producing a better estimate of ITF transport.
The original island rule suggests that the major forcing for the ITF transport (T=y 0, wind ) is the wind field in the eastern ocean of Australia-PNG. In other words, the magnitude and variability of the ITF would depend only on the zonal wind stress across the Pacific at the northernmost and southernmost latitudes of Australia-PNG and the alongshore wind stress along the western coasts of Australia and South America (Wajsowicz, 1993). The wind in the southern hemisphere definitely plays an important role in determining the ITF transport. However, the upstream of ITF include not only the southern Pacific ocean circulation, but also the major currents in the norther Pacific Ocean which are mainly controlled by wind forcing from the norther Pacific Ocean. For example, the MC, which comes from the wind-driven current the NEC (McCreary and Lu, 1994), is an important source of the ITF but controlled by the northern Pacific winds (e.g., Hu et al., 2021). Godfrey et al. (1993) proposed that the ITF is supplied by the MC, and Gordon (1995) suggested that the ITF mainly originates from the North Pacific Ocean, and the seawater from the MC is mainly transported southward through Makassar Strait.
In addition, the northern Pacific Ocean also feeds the ITF through the South China Sea (e.g., Gordon et al., 2019). The Kuroshio flows into the South China Sea through the Luzon Strait in winter under the influence of the Northeast monsoon (Wyrtki, 1961). Lebedev and Yaremchuk (2000) pointed out that the inflow of the Luzon Strait contributes significantly to the ITF, the seawater from the North Pacific flows into the South China Sea through the Luzon Strait, and part of it flows south through the Mindoro Strait, and then farther on until it joins a powerful current, which can be seen as an extension of the MC. So, changes of the North Pacific wind field cause the response of MC and Kuroshio and further affect the ITF.
As can be seen from Eq. (45), considering that the ITF volume transport of multiple islands is determined by the wind stress field in the South Pacific (0°-45°S) and the North Pacific (0°-20°N), the transport value includes the flow through Mindoro Strait (i.e., T=y 1, wind −y 2, wind , with a source of the South China Sea) and Makassar Strait (i.e., T=y 0, wind −y 2, wind , with a source of the MC). So, the optimizing of the path of integration related to the geometry (i.e., including the northern Paicific Ocean) expectedly acts to improve the ITF estimates.

Sensitivity to boundary layer thickness and channel size
In the case of multiple islands, several factors, including bottom and lateral friction, channel width and length, are able to influence the estimated throughflow. The frictional boundary layer thickness is a monotonically increasing function of the frictional coefficient. To test the sensitivity of ITF transport to frictional coefficient and channel size, we examined the response of ITF transport from the multiple island rule to these factors. For simplification, we assume that the frictional coefficient can be represented by the boundary layer thickness in exploring the relationship between frictional coefficients and ITF transport. Figure 10A shows the relationship between ITF transport and Stommel boundary layer thickness when the channel's width is 200 km. It can be seen from the graph that the transport decreases with the increase of d S , but the decreasing rate gradually reduces. When the channel's width is set to 50 km, Munk-type friction is also taken into account, and the flow decreased to varying degrees with the increase of d S and d M ( Figure 10B). The contribution of the two kinds of friction is different. When the passage is relatively narrow, the Munk type friction dominates and the blockage of the current is stronger. It is interesting that, with the increase of Munk and Stommel boundary layer thickness, the lateral friction's blocking effect on the ocean current is enhanced, while the bottom friction's blocking effect on the ocean current is weakened, which is different from the conclusion when the passage is wider.
We plot the derived ITF transport as a function of the channel's width and length in the bottom panels of Figure 10. When d M <W 1 ≤d s , there is only bottom friction, ITF transport increases with the increase of channel's width, and the growth rate gradually increases ( Figure 10C). When 0<W 1 ≤d M <d s , both bottom friction and lateral friction play a role and the ITF transport still increases with the increase of channel's width, but the rate of increase is accelerated ( Figure 10D).
The reason why the rate of change is inconsistent is that when the channel is relatively narrow (less than d M ), the dominant lateral friction nonlinearity decreases rapidly with the increase of the channel's width, and when the critical width (d M ) is reached, the lateral friction decreases to zero. When the width continues to increase, the current is completely controlled by the bottom friction, and the bottom friction obstructs the current more and more until the current rate changes stabilize. In contrast, the ITF transport is generally decreasing with an increasing of the channel's length, which is not affected by the thickness of the boundary layer, and the rate of decrease gradually reduces ( Figures 10C, D).

Discussion and conclusions
In this study, the island rule theory in the case of complex geometry with multiple islands referring to the Indo-Pacific Maritime Continent is investigated based on Godfrey's island rule theory. To this end, the horizontal momentum equations were vertically integrated, the streamfunctions of islands under five cases with different geometry were calculated, and the transport through the channels was predicted. By comparing the transport in the five cases, we find that the geometry has a substantial influence on the transport of throughflow, due to changes in bottom and lateral friction and integral path associated with the geometry. The transport of the throughflow is negatively correlated with the frictional boundary layer thickness and channel length, but positively correlated with channel width, and the flow varies at different rates in different cases.
The ITF transport is estimated using this multiple islands rule with yearly wind stress data. During 2004During -2006, the derived ITF transport is about 14.5 Sv, which is close to the observed ITF transport from INSTANT in the same duration (about 15.0 Sv). We compared the yearly ITF transport during 1988-2016 with that using G1989 and Y2020, and find that the multiple islands rule in this study is able to infer an improved estimate of ITF transport time series in terms of mean value and interannual variability. The ITF transport derived from the multiple islands rule is 13.00 ± 3.94 Sv, which is very close to direct observations. In addition, the island rule theory is applicable for low-frequency variability on which time scale the ITF adjustment to basin scale wind forcing is basically finished, so the theory is also applicable for variability of lower frequency like decadal variability.
Through a series of comparison between the multiple island rule, previous studies and observations, and sensitivity experiments, we suggest that the optimizing of the fraction and path of integration related to the geometry makes the multiple island rule be able to produce a time series of ITF transport that is in well agreement with observations. The comparison indicates that it is possible to produce a long time series of ITF transport from wind products using the multiple island rule, which is of much importance in monitoring the ITF and understanding the ocean circulation and climate system. It is worth mentioning that this study may have a potential value to the estimation of paleo-ITF transport.
Previous studies suggest that changes, including closures, opening and widening of the Indonesian seaway have happened in history (e.g., Kuhnt et al., 2013). The spatial features of the Indonesian straits might be similar with the paleo-geometry, hence the triple island rule theory is expected to be a good tool to estimate the paleo-ITF transport. However, there is still room for further improvement in the multiple islands rule theory. For instance, in solving the vorticity equation, we assume that the width of the channel is much smaller than the length of the channel for simplification, and the problem has been not solved in the current investigation. In addition, time dependence may be important but is not explored here. In the near future, the theory might be further improved by adding the baroclinic term, time term, and other aspects.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

Author contributions
SH designed the study. WS and SH conducted the analysis and wrote the initial version of the paper. All authors contributed to the article and approved the submitted version.