Fortnightly Variability of Chl-a in the Indonesian Seas

. Twenty years of daily MODIS-Aqua ocean color observations (2002-2022) are used to identify periodic variability of near-surface chlorophyll (Chl-a ) in the Indonesian Seas. The frequency spectrum of Chl-a is dominated by the mean and low-frequency monsoonal variability; however, a prominent peak around the fortnightly tidal period, MS f , is present. Harmonic analysis is used to quantify and map the fortnightly Chl-a signal, which is discovered to be significant along the continental shelves of NW Australia, and at several sites associated with narrow passages between the Lesser Sunda Islands, within the 5 Sulu Archipelago, and at a few other sites in the Philippines Archipelago. Fortnightly variability at the shallow coastal sites is attributed to the spring-neap cycle of barotropic ocean currents, while we hypothesize that the variability in deeper water near the island passages is due to the modulation of vertical nutrient fluxes by baroclinic tidal mixing. The results document the significance of tidal mixing and highlight the heterogenous character of bio-physical processes within the Indonesian Seas.


10
The Indonesian Archipelago is a unique region of the world characterized by the warmest sea surface temperatures (SST) of the ocean, and it is a key region of the climate system (Nagai and Hibiya, 2015). The Indonesian Seas play a critical role in the climate system, with the Indonesian Throughflow (ITF) being the only tropical pathway connecting the global oceans (Sprintall et al., 2014). Pacific warm pool waters passing through the Indonesian Seas are cooled and freshened by strong air-sea fluxes and vertical mixing induced by internal tides to form a unique homohaline water mass (Indonesian Throughflow Water,or 15 ITW) that can be tracked across the Indian Ocean basin and beyond (Atmadipoera et al., 2022). The Indonesian Seas lie at the climatological center of the atmospheric deep convection associated with the ascending branch of the Walker Circulation (Clement et al., 2005). Regional SST variations cause changes in the surface winds that can shift the center of atmospheric deep convection, subsequently altering the precipitation and ocean circulation patterns within the entire Indo-Pacific region (Ropelewski and Halpert, 1987). The tidal vertical mixing induces an averaged cooling of 0.5 • C of the Indonesian warm SST by averaging Lomb-Scargle periograms based on Hamming-windowed data within 120-day segments. A narrow peak is evident at the MSf frequency (corresponding to 14.77 days, the M2-S2 sping-neap cycle), and broadband variance is present close to the first sub-and superharmonics of MSf (denoted Mm and 2MS,indicated). Harmonic analysis at the Mm and 2MS frequencies, and at nearby compound overtide frequencies, did not find significant amplitude. ticular matter, and estimating these quantities from a weigthed average of spectral data. The product provides daily estimates; although, the constituent measurements may occasionally fall outside of the calendar date of the product. For the purpose of 65 harmonic analysis, the measurements are treated as if they occur at 12:00 GMT of the product date. The spatial resolution of the gridded GSM data is approximately 5 km.
A map of the number of Chl-a observations per pixel illustrates the influence of clouds and other processes responsible for data gaps (Figure 1). The region just off the northwest coast of Australia contains the most data, the longest time series comprising about 10 years of observations; although, this is only about 1/2 the total possible if every daily data value were 70 populated. Individual pixel time series typically contain between 500 and 1000 data values.
The choice of frequencies used in the harmonic analysis was informed by previous studies (Susanto et al., 2006;Ray and Susanto, 2016;Xing et al., 2021) and trial-and-error. Figure 2 shows the power spectrum of Chl-a from a region off the northwest coast of Australia. The spectral estimate was computed by averaging the Lomb-Scargle periodogram over ovelapping blocks of 60 observations (120 day time series, typically). The spectrum shows a narrowband peak at MS f , and broadband 75 variability close to the sub-and super-harmonics of this frequency (labelled M m and 2SM, respectively), superimposed on a red broadband spectrum.
The choice of harmonic analysis' frequencies was also guided by experiments with analysis at the frequencies of S a (the annual period), S sa (the semi-annual period), M m , MS f , 2SM, and the annual modulates of MS f . Based on comparison of the harmonic constants with the estimated standard error, and bias errors discussed in the Appendix, only the signals at the S a 80 and MS f frequencies were of unambiguous significance. Thus, the plots below are restricted to the results for the S a and MS f variability. Experiments with harmonic analysis of log 10 (Chl-a) were conducted, but these did not reveal relevant differences in the goodness-of-fit or significance compared to Chl-a, so the latter is shown below.
The least-squares harmonic analysis of Chl-a was conducted using standard methods (Foreman et al., 2009). At each pixel a time series of Chl-a was assembled, y, and the design matrix, A, was constructed with columns consisting of a constant (mean) and the cosine and sine of the astronomical arguments for the annual (S a , Doodson number 0 565 555) and luni-solar fortnightly (MS f , Doodson number 0 735 555) cycles. The vector of in-phase and quadrature harmonic constants, x, was obtained as the least-squares solution of the linear system y = Ax.
Missing data in the ocean color observations are largely attributable to clouds. Because the clouds are associated with changes in weather related to the quasi-periodic monsoon cycle, there is concern that the Chl-a harmonic constants could be 90 biased by the correlations of data gaps with Chl-a variability. To examine this possibility, a simulation study was conducted using a gap-free time-series, namely, reanalysis winds from this region, as described in the Appendix. It was found that bias errors on the scale of 10 to 20% may be expected in the estimated mean and annual cycles of Chl-a over much of the domain.
In contrast, the S sa variability is strongly influenced by data gaps, and for this reason it is not shown. The data gaps can also cause significant error in MS f estimates due to crosstalk from the larger-amplitude low-frequencies. Therefore, MS f estimates 95 are not shown where they are smaller than 10% of the mean Chl-a. The reader is referred to the Appendix, below, for a detailed analysis.
In addition to the systematic error due to data gaps, it is necessary to quantify the random error related to measurement noise and non-tidal variability. Due to data gaps, the standard least-squares estimate computed from the spectrum of the residual is not applicable. Instead, the method of Matte et al. (2013) is used which is equivalent to the standard approach in the limit of 100 gap-free data. The method uses the unitary discrete Fourier transform, denoted with matrix F, to estimate a pseudo-spectrum of the residual. Given the residual, r = y − Ax, a pseudo-spectrum of r is defined as S(r) = Fr · (Fr) * , where · denotes the elementwise Shur product and super-script * denotes the complex conjugate. S(r) is referred to as a pseudo-spectrum because it is computed from unevenly-spaced data; however, it is equal to the periodogram in the case of an evenly-spaced (gap-free) time series. After smoothing S(r) to make the pseudo-spectrum continuous, the error variance of x is computed 105 with a Monte-Carlo estimate of the matrix, A −1 F T diag(S(r))FA −1 . This error estimate assumes that the covariance of the Fourier transform of the residual, Fr, is diagonal, i.e., that the noise in r is uncorrelated between pseudo-frequencies. This error estimate was compared to an approach mentioned in Ray and Susanto (2016), wherein the noise of the tidal harmonic constants was estimated by comparing it with the amplitude estimated at "fake" tidal frequencies, and found to agree well.
The fake frequencies used were S a plus 1 cycle-per-9-years (Doodson number 0 567 555) and MS f minus 3 cycles-per-year 110 (Doodson number 0 205 555). In the plots shown below, the amplitude and phase of the Chl-a harmonic constants are only plotted if the amplitude exceeds twice the estimated error level.  and non-significant amplitude compared to the noise estimate (light gray). Note the different colorscale used in each panel. Figure 3 shows the results of harmonic analysis for the mean Chl-a concentration and the S a and MS f signals. Comparing the mean and S a (panels a and b), it is evident that the seasonal variation of Chl-a is nearly 50% of the mean in many coastal areas, 115 specifically, in the northeast Bay of Bengal, in the northeast Andaman Sea, in the northeast Arafura Sea, and along the southern coast of Java. In contrast, the significant MS f variability of Chl-a is largely restricted to a few distinct sites and regions (panel d), including the northwest Australian shelf, the southern part of the Sulu Sea, between the Lesser Sunda Islands, at sites in the Molucca Seas, and in a few other coastal areas. Several of these regions will be enlarged, below.
Throughout the deep ocean, where its amplitude is small, S a is characterized by a smoothly-varying phase associated with  Figure 1). The largest amplitudes are found in the southern Sulu Sea (on the north side of the archipelago), in an area known for its large-amplitude tidal internal waves (Apel et al., 1985;Zhang et al., 2020;Zhao et al., 2021). The peak amplitude of MS f , 0.2mg/m 3 , is about 1/4 as large as the mean here, so it represents a significant modulation of the Chl-a concentration.

130
The region of maximum MS f amplitude is associated with a phase lag which varies from about 7 days in the east (x = 100 km) to 10 days (i.e, equivalent to -4 days on the colorscale) in the west (x = −100 km). The phase is measured relative to the tidal potential, but Ray and Susanto (2016) show that the barotropic semidiurnal tide at this location lags the potential by about 1.2 days. Thus, the fortnightly maximum in Chl-a occurs from about 6 to 9 days after the maximum spring-neap tides in this region.
135 Figure 5 shows MS f Chl-a in the area around the Lesser Sunda Islands, from Bali to Ombai Strait. The maximum in Chl-a occurs between the islands of Lombock and Sumbawa in the west, and within the northern part of the Ombai Strait in the east. The phase lag exhibits a relatively strong southward gradient, but close to the islands it is about 8 days lagged from the astronomical potential. The barotropic tide lags the potential by about 2 days here (Ray and Susanto, 2016), so the fortnightly Chl-a maximum again lags the spring tides by about 6 days. The phase propagates southward at a rate of roughly 5 days per 140 80 km, equivalent to 0.2m/s.

Discussion
The mean and annual variability of Chl-a shown above replicate the results of Susanto et al. (2006), providing an update with greater resolution based on the longer time series and different satellite missions now available. The use of harmonic analysis to explicitly separate the mean, annual, and fortnightly variability represents an alternative to the climatologies averaged into   calendar months or lunar phase used in previous works (Susanto et al., 2006;Shi et al., 2011). The harmonic analysis of Chl-a, together with its error estimate, appears to provide more reliable classification of significant and non-significant signals, for example, as compared with the wavelet-based approach of Xing et al. (2021), which indicated significant MS f variability in relatively large regions of the deep Indian and Pacific oceans.
The fortnightly maximum of Chl-a lags the spring tide by roughly 6 to 9 days for the examples cited previously in the 150 Southern Sulu Sea and along the Lesser Sunda Islands, and roughly the same phase lag occurs at other sites along the northwest Australian shelf and along the archipelagos in the Molluca Sea. The Chl-a phase lag thus exceeds the 1-to 3-day lag between the spring tide and the SST minimum estimated by Ray and Susanto (2016). The potential relationship between Chl-a and nutrient inputs due to tidal mixing is complex (Laws, 2013), but the timing of the fortnightly Chl-a maximum is logically consistent with the co-occurance of vertical nutrient and heat fluxes associated with tidal mixing, leading to a productivity 155 maximum and Chl-a peak a few days later.
In some respects the spatial distribution of MS f Chl-a variability coincides with the map of sub-surface tidal dissipation inferred in the numerical model of Nugroho et al. (2018) (their Figure 7), but it is notable that the MS f Chl-a signal is absent in many regions where the dissipation is large, e.g., in the Strait of Malacca, in the Makassar Strait, in the northern Arafura Sea, etc.

160
In order to understand the mechanisms of Chl-a variability, Figure 6 illustrates two diagnostics computed from the barotropic currents predicted by the TPXO9 tidal model (Egbert and Erofeeva, 2002). The first diagnostic is denoted ∆U , and it is the barotropic tidal current speed associated with the spring-neap cycle (Figure 6a). Assuming that the energy source that drives  the fortnightly mixing is proportional to the kinetic energy of the combined M 2 and S 2 tides, 165 where U i is the major axis of the tidal current, ω i is the frequency, ϕ i is the phase lag, and subscripts 1 and 2 correspond to the M 2 and S 2 tides, respectively, then the amplitude of the spring-neap cycle of KE is given by the product, (∆U ) 2 = U 1 U 2 . In contrast, the baroclinic tide is generated by flow across topography, leading to localized vertical displacements of isopycnals and baroclinic pressure gradients. A metric of this quantity is ∆Z = w/ω 1 , the vertical displacement of a water parcel flowing along the bottom (Figure 6b), where w = |w 1 w 2 | 1/2 is given in terms of the cross-isobath M 2 and S 2 currents, w i = −u i ∇H, 170 with u i the major-axis current vector. where the MS f Chl-a amplitude is not significant. Indeed many of these same regions are well-known for the presence of nonlinear internal waves (Osborne and Burch, 1980;Zhao and Alford, 2006), and it is puzzling that a robust MS f Chl-a signal is not detected at these sites.
A more direct view of baroclinic tidal processes in the region is provided by satellite altimetry (Figure 7). The nearly 30- year record of data along the TOPEX/Jason ground tracks is long enough that the sea level anomaly phase-locked with the 180 astronomical tide can be identified with sub-centimeter precision (Zaron, 2019). The left panels illustrate that the M 2 sea level anomaly in some places exceeds 5 cm (south of Lombock in Figure  frequency represent the high-frequency counterpart of the nonlinearity which should also manifest at the MS f frequency. Although the altimeter ground tracks are spaced too far apart to map the MS 4 waves in two dimensions, their along-track wavelength is consistent with the interpretation that they are mode-1 internal waves (59 to 66 km, based on the World Ocean Atlas stratification; Zaron et al. 2022). The decay of the waves southward from Lombock, northward from Ombai, and southward from the Sulu Archipelago clearly implicates narrow passages between the islands as their sites of genesis. 190 We hypothesize that the fortnightly modulation of surface Chl-a identified from ocean color is the result of tidally-driven mixing, which supports a fortnightly modulation of the vertical nutrient transport. The locations with significant fortnightly Chl-a variability, aside from the northwest shelf of Australian and the northern Andaman Sea, nearly all coincide with regions where strong fortnightly modulation of the cross-isobath flow is expected to occur (Figure 6b). There is evidence for nonlinear interactions at a subset of these sites, from sparse altimetry, at the sum of the M 2 and S 2 frequencies (Figure 7b and d).

195
Field work would be required to assess the validity of this hypothesis, since it presupposes the existence of a nutrientlimited standing stock of phytoplankton capable of assimilating the periodic influx of nutrients, a vertical nutrient gradient, and a marginally-stable water column susceptible to destabilization by tidal currents. Among the high-resolution tide-resolving numerical models of the Indonesian Seas (e.g., Robertson and Ffield 2008;Nagai and Hibiya 2015;Hermansyah et al. 2019), the MS 4 overtides have not been specifically studied, but the altimeter-derived estimates for this tide could be used to validate 200 the nonlinear dynamics resolved in these models. One plausible mechanism for tidally-induced mixing involves near-bottom overturning initiated where cross-isobath particle displacement is comparable to the vertical motion of the linear internal tide (Levine and Boyd, 2006;Klymak et al., 2008). Another mechanism, which would lead to mixing higher in the water column, is the convective instability of the massive internal waves which have been observed in locations of strong tidal currents in this region (Atmadipoera and Suteja, 2018;Bouruet-Aubertot et al., 2018). One limitation of the present estimates is the difficulty of identifying the small MS f Chl-a signal from the gappy data.
Experiments to improve the estimates were conducted using spatially-coupled harmonic analysis (not shown); however, results 220 were found to be overly sensitive to details of the spatial basis functions used. Ocean models with the capability of simulating the tidal-biological interactions are now being implemented (e.g., Gutknecht et al. 2016), and it may be fruitful to use the surface Chl-a fields from these models as spatial basis functions within a spatially-coupled harmonic analysis, as was recently done for sea surface-height data (Egbert and Erofeeva, 2021). While imperfect, such an approach may provide the capability to utilize the gappy data more effectively to map periodic components of the Chl-a fields.
Appendix A: The impact of data gaps on harmonic analysis Building on the notation of the main text, consider a data indicator vector, w, the same size as the data vector, y, with entries zero or one depending on whether the data value is valid (i.e., used in the harmonic analysis) or invalid (i.e., omitted as part 230 of a data gap), respectively. Define the rectangular diagonal matrix,Ŵ, which maps from the original data vector, y, to the gappy data vector by picking out those entries which correspond to nonzero entries in the w vector. Given the data vector, y, the corresponding gappy data vector shall be denotedŷ =Ŵy, and the design matrix for the harmonic analysis of the gappy data shall be denotedÂ =ŴA. For convenience, define the square diagonal weighting matrix, W =Ŵ TŴ , which contains the elements of w along its diagonal. Finally, let W c = I − W denote the complement of the weight matrix with diagonal w c ,

235
where w c = 1 − w is the complement of the data indicator vector.
With the above extensions of the original notation, it is possible to approximately compute the sensitivity of the harmonic constants to the data gaps. The harmonic constants based on the gappy data are given by, A Neumann series can be used to approximate the inverse, above, in terms of powers of W c :

245
where the dots denote quadratic and higher-order terms in W c . Substituting this approximation in equation (A3), one obtains, Note that y − Ax is the residual, r, data minus predicted signal, obtained from the harmonic analysis of the gap-free data.
Equation (A11) shows that if the gap-free time series is described exactly by the harmonic components (i.e., r = 0), then the harmonic constants are insensitve to data gaps. The amplitude and phase of each harmonic component is described by exactly two parameters, therefore, as long as the gaps do not reduce the quantity of data below the number of parameters to be 255 estimated (in which caseÂ TÂ would not be invertable) the harmonic constants are uniquely determined, independent of the gaps. (b) Power-spectral density of data indicator function Figure A1. (a) The count of good data, i.e., cloud-free pixels, exhibits a notable annual cycle associated with monsoon variability. Data counts are about 30% greater during the dry season than during the monsoon season (roughly September through March). (b) The power spectral density of the data indicator function, w, contains significant peaks at the frequency of the annual cycle (Sa) and its harmonics (the semi-annual cycle, Ssa, is indicated). A peak at the fortnightly frequency (MSf) is present, but it is almost 3 orders of magnitude smaller than the peak at the zero frequency (98 samples 2 /cpd, marked with a dot). The spectrum is computed by averaging the Hamming-windowed periodograms computed from data-indicator time series at each ocean pixel.
Of course, with real data, the residual will be nonzero whether or not the time series contains gaps. In this case it is useful to note that W c r is equal to the elementwise (Schur) product of w c and r, a time-domain product of the complement of the data-gap indicator and the residual time series. Using the convolution theorem, this time-domain product can be interpreted as   (Murakami and Matsumoto, 1994), also called the dry southeast monsoon (Aldrian 265 and Susanto, 2003). Indeed, the power spectrum of w c averaged over all the grid cells of the Chl-a dataset exhibits peaks at the S a and S sa periods, the larger of which contains about 15% as much power as the zeroth frequency ( Figure A1b). The power at the MS f frequency of the spring-neap cycle is about 0.3% that at the zeroth frequency. While the latter seems negligable, the expected crosstalk is proportional to the square root of this quantity, about 0.05, and the MS f amplitude is so much smaller than the mean that this crosstalk could be sigificant. Thus, as mentioned in the text, MS f harmonic constants are considered 270 non-significant and not plotted if they are smaller than 0.1 of the mean Chl-a value in a given pixel.
As an aside, the MS f peak in the data indicator spectrum suggests that the data flagging algorithm is influenced by tides in this region. It is conceivable that fortnightly modulations in SST could affect atmospheric stability and cloudiness (Ray and Susanto, 2019); however, it is also possible that some other process affecting ocean color, such turbidity, might change the outlier flagging in the GSM algorithm.

275
Error of mean, Z0 Figure A2. The relative error in the estimated mean (Z0) of the meridional wind shows the effect of Chl-a data-gaps at each pixel. For this example, the time series of the NCEP-reanalysis 850mb meridional wind at (7.5 • S,110 • E) is taken as the reference gap-free "truth", and the mean is compared with the mean estimated from the gappy time series at each oceanic pixel, according to the Chl-a data indicator. The error is a measure of the crosstalk of non-sinusoidal (broadband) variability near the Sa frequency which is aliased back to the mean by the gappy sampling. This map shows where data gaps caused by cloudiness are likely to bias the mean estimated for Chl-a.
The concern that data gaps could lead to errors in the S a and S sa harmonic constants is also warranted, but it is difficult to estimate the magnitude of this effect for Chl-a since the gap-free residual is unknown. In order to assess the potential crosstalk between these frequencies, we have made harmonic analyses of a known signal, in this case, the vector winds at a location in the Java Sea, where the broadband (non-sinusoidal) signals associated with the monsoon cycle are expected to be reasonably similar to the broadband monsoon signals in Chl-a, assuming the latter are largely the result of wind-driven variability. The 280 wind time series is taken from the NCEP reanalysis product (Kanamitsu et al., 2002), so the gap-free time series and its harmonic analysis are available for comparison.
The (gap-free) time series of zonal wind has values of -1.5 m/s, 5.2 m/s, and 2.2 m/s for its mean and amplitudes of the S a and S sa harmonic constants, respectively. It thus differs from the Chl-a time series in which the mean is about twice as large as the S a amplitude. In contrast, the (gap-free) meridional wind has values of 0.8 m/s, 0.6 m/s, 0.15 m/s, for the mean, S a , and 285 S sa amplitudes, respectively, which are in the same order, by size, as the Chl-a harmonic constants. Results from the analysis of the meridional wind component are shown below, since these should be more representative of the errors due to data gaps in Chl-a. Figure A2 shows the relative error of the mean wind due to the data gaps, (x 1 /x 1 − 1) × 100%, wherex 1 and x 1 are means estimated from the gappy and gap-free time series, respectively. The typical magnitude of the error, about 15%, is what could 290 be expected from the crosstalk of the S a and S sa signals with the mean. As a consequence of the spatial correlations of the data gaps, the errors are also spatially correlated. The largest errors occur in the northern Gulf of Thailand, in the northern Andaman Sea, and near the coastline in the Bay of Bengal. The relative sizes of the S a and S sa amplitudes compared to the mean are smaller for Chl-a than for the wind, so the bias errors in the estimated Chl-a mean field are likely to be smaller than shown, but it is nonetheless clear that the errors due to data gaps could be a significant fraction of the mean.

295
The crosstalk from the mean and S sa into the S a signal is illustrated in Figures A3a and A3b; note the larger range of colorscales compared to Figure A2. These panels show, again, a geographic pattern with maximal errors in many of the same regions as above. While errors in the phase are not negligable, they are generally small enough and spatially noisy enough, that one would not likely misrepresent the seasonal phasing of the signal, assuming the Chl-a errors are comparable.
In contrast, Figures A3c and A3d indicate that the data gaps are likely to make the S sa variability of Chl-a impossible to 300 interpret. Indeed, an analysis of the formal (random error) in the S sa harmonic constants suggests that the signal-to-noise is below 1 in about half of the domain (not shown), and Figure A3c indicates that the systematic bias due to data gaps is likely to be severe. For this reason, the S sa analysis is not shown in the text; although, the harmonic constants were computed and are available from the authors.
To conclude, this appendix provides an analysis of the effect of data gaps on the harmonic analysis of Chl-a. Because the 305 cloudiness is correlated with the seasonal and semi-annual cycles, it has the potential to create systematic errors in the seasonal variability estimated via harmonic analysis. Because the magnitude of the effect depends on the residual, i.e., on the component of the signal which is not phase-locked with the components included in the harmonic analysis, it is difficult to analyze without a comparable gap-free dataset. The approach taken here has involved the harmonic analysis of a gap-free dataset consisting of reanalysis winds, which should have low-frequency broadband properties similar to that of Chl-a. Based on the present 310 analysis, the mean and annual cycles of Chl-a are presented in the main text, for comparison with previous works, but the S sa variability is not shown since it is unlikely to be reliably determined.
Author contributions. Zaron drafted the manuscript and performed the data analysis. Capuano provided the Chl-a dataset and contributed to revising the manuscript. Koch-Larrouy conceived of the objective of the research, guided the selection of significant geographic regions for study, and contributed to revising the manuscript.  Figure A3. The relative error due to data gaps in amplitude and phase for the harmonic constants at the Sa (panels a and b) and Ssa (panels c and d) frequencies, based on harmonic analysis of a representative time series of NCEP reanalysis meridional winds are shown, similar to Figure A2. Assuming the low frequency broadband variability of the wind is similar to Chl-a, the Ssa harmonic constants are judged as not reliable, and are thus omitted from the presentation in the main text.