<figure><img src="image_1.png" /></figure>

HAL
open science

# Fortnightly variability of Chl a in the Indonesian seas

Edward. D. D. Zaron, Tonia. A. A. Capuano, Ariane Koch-Larrouy

► To cite this version:

Edward. D. D. Zaron, Tonia. A. A. Capuano, Ariane Koch-Larrouy. Fortnightly variability of Chl a in the Indonesian seas. Ocean Science, 2023, 10.5194/os-19-43-2023. hal-04720390

HAL Id: hal-04720390

https://hal.science/hal-04720390v1

Submitted on 11 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Fortnightly variability of Chl *a* in the Indonesian seas

Edward D. Zaron<fnref n="1" />, Tonia A. Capuano<fnref n="2" />, and Ariane Koch-Larrouy<fnref n="2" />

**Correspondence:** Edward D. Zaron (edward.d.zaron@oregonstate.edu)

Received: 23 August 2022 – Discussion started: 31 August 2022

Revised: 17 December 2022 – Accepted: 20 December 2022 – Published: 17 January 2023

**Abstract.** 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, MSf, 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 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 heterogeneous character of biophysical processes within the Indonesian seas.

## 1 Introduction

The Indonesian Archipelago is a unique region of the world characterized by the warmest sea surface temperatures (SSTs) 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 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 Indonesian region (Ropelewski and Halpert, 1987). The tidal vertical mixing induces an averaged cooling of 0.5 °C of the Indonesian warm SST (Koch-Larrouy et al., 2008), which in turn reduces the precipitation by 20% (Koch-Larrouy et al., 2010) and regulates the amplitude and variability of El Niño–Southern Oscillation (ENSO), the Indian Ocean Dipole (IOD), and the Madden–Julian Oscillation (MJO) (Sprintall et al., 2019). Indeed, models that include tidal mixing have improved representation of ENSO, IOD, and MJO patterns (Susanto et al., 2005; Sprintall et al., 2019).

Previous studies have shown that SST patterns in the Indonesian seas are highly influenced by tidal mixing as well as monsoonal winds, surface heat fluxes, and the ITF (Koch-Larrouy et al., 2010; Kida and Wijffels, 2012; Sprintall et al., 2014, 2019). The spring-neap cycle of the tides has been measured in the SST using satellite data (Ray and Susanto, 2016) and simulated using models (Nugroho et al., 2018). In fact, only model simulations that include tides (Koch-Larrouy et al., 2007; Nugroho, 2017; Nagai et al., 2017, among others) are capable of reproducing the homohaline water (ITW).

However, the complexity of the atmospheric and ecosystem's response to tidal mixing remains poorly understood and is a topic of ongoing research. It is a challenge to disentangle the variability connected with long-term trends and monsoonal and seasonal processes, which are each heterogeneous.

<footnote n="1">¹College of Earth, Ocean and Atmospheric Science, Oregon State University, Corvallis, Oregon, USA</footnote>

<footnote n="2">²Laboratory for Studies in Geophysics and Spatial Oceanography (LEGOS), Université de Toulouse, Toulouse, France</footnote>
neous across the domain (Susanto and Ray, 2022). Recent joint analysis of chlorophyll *a* (Chl *a*), winds, ocean currents, and active tracers (Mandal et al., 2022); Chl *a*, winds, sea surface temperature, and sea level (Wirasatriya et al., 2021); Chl *a*, climate indices, precipitation, and river runoff (Siswanto et al., 2020); and sea level, winds, and stratification (Devlin et al., 2018) are all examples of multi-variable approaches to understanding the low-frequency dynamics of this complex region.

The goal of this study is to document the variability of upper-ocean Chl *a* concentration associated with the spring–neap tides throughout the Indonesian seas and nearby waters (from 23° S to 28° N and 87 to 148° E). Spring–neap variability arises when the vector sum of currents at the *M*₂ and *S*₂ tidal frequencies are rectified or otherwise interact to produce modulations of bottom stress or internal wave-driven mixing with a period of 14.77 d (denoted with the Darwin symbol, MSf), phase-locked with the astronomical forcing. At least three different mechanisms could explain the association of increased Chl *a* with spring–neap currents (Xing et al., 2021): a modulation of the tidally driven turbulent flux of nutrients into the euphotic zone, a modulation of suspended sediment concentration affecting the intensity of photosynthetically active radiation (which ought to decrease apparent Chl *a*), and a modulation of resuspended Chl *a* associated with the microphytobenthos. The present study uses daily averages of Chl *a* inferred from multi-spectral imagery acquired by the MODIS-Aqua mission from 2002 to 2022. It is motivated by several ongoing studies which seek to model the physical and biological dynamics of the ITF (Gutknecht et al., 2016; Nugroho et al., 2018). The results focus on observations of Chl *a* variability rather than its causes; however, satellite altimetry provides independent evidence for nonlinear interaction of *M*₂ and *S*₂ currents and suggests several sites where field investigations could identify the drivers and assess the dynamics of Chl *a* variability in detail.

Previous studies of Chl *a* variability employing satellite ocean color measurements have mentioned the challenges associated with data gaps due to cloudiness, sunglint, and water-type variability (Palacios, 2004; Susanto et al., 2006; Garcia and Garcia, 2008). A variety of approaches to overcoming these difficulties has been used to either fill data gaps (Yang et al., 2021) or phase-average the data according to days after full moon, for example, in the analysis of suspended sediment (Shi et al., 2011). Recently, another approach based on the wavelet analysis of 4.5 years of Himawari-8 imagery was used to study Chl *a*, which was capable of identifying seasonal and spring–neap variability (Xing et al., 2021). Given the long, 20-year, duration of the MODIS records considered here, we have found that least-squares harmonic analysis is adequate to quantify variability, even with the data gaps and noise, as detailed below. A similar approach highlighted spring–neap variability from 12 years of high-resolution SST data in the same region (Ray and Susanto, 2016).

<figure><img src="image_5.png" /><figcaption>**Figure 1.** The percentage of valid Chl *a* observations per pixel in the Level-3 GlobColour daily MODIS product. The rectangular boxes indicate regions enlarged in subsequent figures. The maximum possible observation count, if there were no data gaps, is 7264 for the time period of the observations (July 2002 to August 2022). The red marker at 21° S, 114° E indicates the location of the time series used for Fig. 2.</figcaption></figure>

## 2 Data and methods

The Chl *a* dataset used here consists of the Level-3 “Glob-Colour daily MODIS product”, version 2018.4, which is an estimate of the near-surface Chl *a* concentration in units of milligrams per cubic meter obtained from MODIS-Aqua data using a semi-analytical model (the Garver–Siegel–Maritorena (GSM) product; Maritorena and Siegel, 2005; Maritorena et al., 2010). The GSM product is based on representing the spectrum of outgoing radiation with contributions due to Chl *a*, dissolved organic matter, and suspended particular matter and estimating these quantities from a weighted 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 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 (Fig. 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 half the total possible if every daily data value were 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. Fig-
<figure><img src="image_6.png" /></figure>

**Figure 2.** Spectral estimate of Chl *a* concentration off the northwest coast of Australia (21° S, 114° E). The spectral estimate is computed by averaging Lomb–Scargle periodograms based on Hamming-windowed data within 120 d segments. A narrow peak is evident at the MS<sub>f</sub> frequency (corresponding to 14.77 d, the *M*<sub>2</sub>–*S*<sub>2</sub> spring–neap cycle), and broadband variance is present close to the first sub- and super-harmonics of MS<sub>f</sub> (denoted *M*<sub>m</sub> and 2SM). Harmonic analysis at the *M*<sub>m</sub> and 2SM frequencies, and at nearby compound overtide frequencies, did not find significant amplitude.

ure 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 overlapping blocks of 60 observations (120 d time series, typically). The spectrum shows a narrowband peak at MS<sub>f</sub>, and broadband variability close to the sub- and super-harmonics of this frequency (labeled *M*<sub>m</sub> 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), Ssa (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 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₁₀(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 lunisolar fortnightly (MS_f, Doodson number 0 735 555) cycles. The vector of in-phase and quadrature harmonic con-

stants, **x**, was obtained as the least-squares solution of the linear system **y** = **A****x**.

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 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 Ssa variability is strongly influenced by data gaps, and for this reason it is not shown. The data gaps can also cause significant error in MSf estimates due to crosstalk from the larger-amplitude low frequencies. Therefore, MSf estimates are not shown where they are smaller than 10% of the mean Chl *a*. The reader is referred to the Appendix 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 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** − **A****x**, a pseudo-spectrum of **r** is defined as S(**r**) = **F****r** · (**F****r**)*, where · denotes the element-wise Schur product and superscript * 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 with a Monte Carlo estimate of the matrix: A⁻¹**F**ᵀdiag(S(**r**))**F**A⁻¹. This error estimate assumes that the covariance of the Fourier transform of the residual, **F**r, 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 one cycle per 9 years (Doodson number 0 567 555) and MS_f minus three cycles per year (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.
## 3 Results

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, 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 Sea, and in a few other coastal areas. Two of these regions are enlarged in Figs. 4 and 5, below.

Throughout the deep ocean, where its amplitude is small, S_a is characterized by a smoothly varying phase associated with the winter season poleward of about 20° latitude in both hemispheres (Fig. 3c). Equatorward of 20°, the annual cycle peaks in June–September (Northern Hemisphere summer) throughout much of the Indonesian seas (Susanto et al., 2006); however, strong spatial gradients are associated with the northwest Australian shelf and other coastal areas. The seasonal cycle peaks earlier in the year (January–February) near the island of Borneo, in the Sulu Sea, and further eastward in the Pacific. The relationship between monsoonal and low-frequency variability of Chl a is described in detail in Mandal et al. (2021), Wirasatriya et al. (2021), and Mandal et al. (2022).

Enlarged maps of the MS_f amplitude and phase around the Sulu Archipelago are shown in Fig. 4 (region indicated in Fig. 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.2 mg m⁻³, is about one-quarter as large as the mean here, so it represents a significant modulation of the Chl a concentration. The region of maximum MS_f amplitude is associated with a phase lag which varies from about 7 d in the east (x = 100 km) to 10 d (i.e., equivalent to -4 d on the color scale) 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 d. Thus, the fortnightly maximum in Chl a occurs from about 6 to 9 d after the maximum spring-neap tides in this region.

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 Lombok 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 d lagged from the astronomical potential. The barotropic tide lags the potential by about 2 d here (Ray and Susanto, 2016), so the fortnightly Chl a maximum again lags the spring tides by about

6 d. The phase propagates southward at a rate of roughly 5 d per 80 km, equivalent to 0.2 m s⁻¹.

## 4 Discussion

The mean and annual variability of Chl a shown above replicates 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 d for the examples cited previously in the 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 Molucca Sea. The Chl a phase lag thus exceeds the 1 to 3 d 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-occurrence of vertical nutrient and heat fluxes associated with tidal mixing, leading to a productivity 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 Fig. 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, and in the northern Arafura Sea.

In order to understand the mechanisms of Chl a variability, Fig. 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 (Fig. 6a). Assuming that the energy source that drives the fortnightly mixing is proportional to the kinetic energy (KE) of the combined M_2 and S_2 tides,

$$KE = \frac{1}{2}(U_1 \cos(\omega_1 t - \phi_1) + U_2 \cos(\omega_2 t - \phi_2))^2, \quad (1)$$

where $U_i$ is the major axis of the tidal current, $\omega_i$ is the frequency, $\phi_i$ is the phase lag, and subscripts 1 and 2 correspond to the M_2 and S_2 tides, respectively, then the ampli-
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Results of harmonic analysis of Chl *a* concentration: (a) mean, (b) amplitude of annual cycle ($S_a$), (c) phase of $S_a$, measured relative to the vernal equinox (20 March; brown corresponds to a maximum in February–April, black to May–July, blue to August–October, and white to November–January), and (d) amplitude of the fortnightly cycle ($MS_f$). Two shades of gray are used to indicate land (dark gray) and non-significant amplitude compared to the noise estimate (light gray). Note the different color scale used in each panel.</figcaption></figure>

tude of the spring–neap cycle of KE is given by the prod-
uct ($\Delta U)^2 = U_1U_2$. In contrast, the baroclinic tide is gen-
erated by flow across topography, leading to localized verti-
cal displacements of isopycnals and baroclinic pressure gra-
dients. A metric of this quantity is $\Delta Z = w/\omega_1$, the verti-
cal displacement of a water parcel flowing along the bot-
tom (Fig. 6b), where $w = |w_1w_2|^{1/2}$ is given in terms of the
cross-isobath $M_2$ and $S_2$ currents, $w_i = -u_i \nabla H$, with $u_i$ the
major-axis current vector.

Figure 6a shows that ΔU is greatly enhanced in several
areas where the amplitude of Chl *a* MSf is large, espe-
cially along the northwest Australian shelf and in the Gulf
of Martaban at the north end of the Andaman Sea. In con-
trast, Fig. 6b shows that large values of ΔZ coincide with

large-amplitude MSf Chl *a* along the Lesser Sunda Islands and across sills adjacent to the Molucca and Sulu seas. Just as with the Nugroho et al. (2018) dissipation estimate, it is also notable that there are a number of sites where ΔZ is large, among the Andaman–Nicobar Islands and within Lushon Strait, where the MSf 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 MSf Chl *a* signal is not detected at these sites.

A more direct view of baroclinic tidal processes in the re-
gion is provided by satellite altimetry (Fig. 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
<figure><img src="image_8.png" /><figcaption>**Figure 4.** (a) Amplitude and (b) phase of the fortnightly MS<sub>f</sub> harmonic constants of Chl *a* along the Sulu Archipelago (indicated in Fig. 1).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** (a) Amplitude and (b) phase of the fortnightly MS<sub>f</sub> harmonic constants of Chl *a* along the Lesser Sunda Islands, from Bali to Ombai Strait (indicated in Fig. 1).</figcaption></figure>

the astronomical tide can be identified with sub-centimeter
precision (Zaron, 2019). Figure 7a and c illustrate that the
*M*₂ sea level anomaly in some places exceeds 5 cm (south of
Lombok in Fig. 7a and in the western Celebes Sea in Fig. 7c),
which is evidence of the strong and localized baroclinic tidal
generation in this region. Figure 7b and d show the presence of
small but clearly identifiable waves at the MS₄ overtide
frequency in the same locations just noted, as well as north
of Ombai Strait, in the western Sulu Sea (Fig. 7b), in the
Molucca Sea, and in the eastern Banda Sea (Fig. 7d). These
waves at the MS₄ frequency represent the high-frequency
counterpart of the nonlinearity which should also manifest
at the MSf frequency.

Although the altimeter ground tracks are spaced too far
apart to map the MS₄ 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 de-
cay of the waves southward from Lombok, northward from
Ombai, and southward from the Sulu Archipelago clearly im-
plicates narrow passages between the islands as their sites of
genesis.

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 co-
incide with regions where strong fortnightly modulation of
the cross-isobath flow is expected to occur (Fig. 6b). There is
evidence for nonlinear interactions at a subset of these sites,
from sparse altimetry, at the sum of the *M*₂ and *S*₂ frequen-
cies (Fig. 7b and d).
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Diagnostics of the barotropic tidal currents, based on the TPXO9-Atlas tide model (Egbert and Erofeeva, 2002). (a) The amplitude of the combined $M_2-S_2$ current at the MSf period. Solid contours show the coastline and 50, 200, and 2000 m isobaths. (b) The amplitude of the vertical particle displacement due to the combined $M_2-S_2$ bottom flow across isobaths.</figcaption></figure>

Fieldwork would be required to assess the validity of this hypothesis, since it presupposes the existence of a nutrient-limited 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₄ overtides have not been specifically studied, but the altimeter-derived estimates for this tide could be used to validate 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).

## 5 Conclusions

Harmonic analysis has been conducted of 20 years of daily MODIS-Aqua-derived Chl *a* measurements to identify regions of tidal influence on biological production in the Indonesian seas and the nearby Indian and Pacific oceans. The choice of analysis frequencies, $S_a$ and MSf, was dictated by the spectral resolution feasible with the time series, which suffer from extensive gaps due to clouds. The results agree with previous estimates of $S_a$ (Susanto et al., 2006) and reveal several sites and regions where the Chl *a* concentration is phase-locked with the spring-neap cycle, peaking around 6 to 9 d after the local maximum of tidal currents.

The locations of the largest MSf signals, off the north-west coast of Australia and in the vicinity of narrow channels between islands, generally correspond to regions where boundary-layer-driven and internal wave-driven mixing are expected to occur, based on comparisons with models. However, there are a number of regions where tidal dissipation is large or where large-amplitude internal waves are well-known but where a significant MSf Chl *a* amplitude is not found. This observation highlights the obvious point that tidal controls on Chl *a* occur only at sites where the environmental prerequisites for such control are suitable.

One limitation of the present estimates is the difficulty of identifying the small MSf Chl *a* signal from the gappy data. Experiments to improve the estimates were conducted using spatially coupled harmonic analysis (not shown); however, results 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 0 or 1 depending on whether the data value is valid (i.e., used in the harmonic analysis) or invalid (i.e., omitted as part of a data gap), respectively. Define the rectangular
<figure><img src="image_11.png" /><figcaption>Figure 7. Harmonic analysis of along-track altimeter data. Panels (a) and (c) show sea level anomaly associated with the M2 astronomical tide. Panels (b) and (d) show sea level anomaly associated with the MS4 nonlinear overtide. Ascending and descending orbit ground tracks are shown in the upper and lower panels, respectively. The background image is reproduced enlarged from Fig. 3d, with gray values plotted in white for clarity.</figcaption></figure>

diagonal matrix, $\hat{W}$, 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 $\hat{y} = \hat{W}y$, and the design matrix for the harmonic analysis of the gappy data shall be denoted $\hat{A} = \hat{W}\hat{A}$. For convenience, define the square diagonal weighting matrix, $\hat{W} = \hat{W}^T\hat{W}$, 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$, 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

$$ \hat{x} = (\hat{A}^T \hat{A})^{-1} \hat{A}^T \hat{y}, \qquad (A1) $$

$$ = (\hat{A}^T \hat{W} \hat{A})^{-1} \hat{A}^T \hat{W} y, \qquad (A2) $$
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**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 ($S_a$) and its harmonics (the semi-annual cycle, $S_{sa}$, is indicated). A peak at the fortnightly frequency ($MS_f$) is present, but it is almost 3 orders of magnitude smaller than the peak at the zero frequency (98 samples$^2$ cpd$^{-1}$ (cpd - cycles per day), marked with a dot). The spectrum is computed by averaging the Hamming-windowed periodograms computed from data indicator time series at each ocean pixel.</figcaption></figure>

$$= (A^T(I - W_c)A)^{-1}A^T(I - W_c)y. \qquad (A3)$$

A Neumann series can be used to approximate the inverse, above, in terms of powers of $W_c$:

$$(A^T(I - W_c)A)^{-1} = (A^T A - A^T W_c A)^{-1}, \qquad (A4)$$

$$= [(A^T A)(I - (A^T A)^{-1}A^T W_c A)]^{-1}, \qquad (A5)$$

$$= (I + (A^T A)^{-1}A^T W_c A + \dots)(A^T A)^{-1}, \qquad (A6)$$

where the dots denote quadratic and higher-order terms in $W_c$. Substituting this approximation in Eq. (A3), one obtains

$$\hat{x} = (I + (A^T A)^{-1}A^T W_c A + \dots)(A^T A)^{-1} \cdot A^T(I - W_c)y, \qquad (A7)$$

$$= (A^T A)^{-1}A^T(I - W_c)y + (A^T A)^{-1}A^T \cdot W_c A(A^T A)^{-1}A^T(I - W_c)y + \dots, \qquad (A8)$$

$$= x - (A^T A)^{-1}A^T W_c y + (A^T A)^{-1}A^T \cdot W_c A(A^T A)^{-1}A^T y + \dots, \qquad (A9)$$

$$= x - (A^T A)^{-1}A^T W_c (y - A(A^T A)^{-1}A^T y) + \dots, \qquad (A10)$$

$$= x - (A^T A)^{-1}A^T W_c (y - A(A^T A)^{-1}A^T y) + \dots, \qquad (A11)$$

$$= x - (A^T A)^{-1}A^T W_c (y - A(A^T A)^{-1}A^T y) + \dots, \qquad (A11)$$

$$\hat{x} - x = -(A^T A)^{-1}A^T W_c (y - A x) + \dots, \qquad (A11)$$

$$\hat{x} - x = -(A^T A)^{-1}A^T W_c (y - A x) + \dots, \qquad (A11)$$

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 insensitive 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 estimated (in which case $\hat{A}^T \hat{A}$ would not be invertible), the harmonic constants are uniquely determined, independent of the gaps.

<figure><img src="image_14.png" /><figcaption>**Figure A2.** The relative error in the estimated mean ($Z_0$) 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 850 mbar 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 $S_a$ 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.</figcaption></figure>

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 element-wise (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 a convolution in the frequency domain. The
<figure><img src="image_15.png" /><figcaption>**Figure A3.** The relative error due to data gaps in amplitude and phase for the harmonic constants at the S<sub>a</sub> (a, b) and S<sub>sa</sub> (c, d) frequencies, based on harmonic analysis of a representative time series of NCEP reanalysis meridional winds are shown, similar to Fig. A2. Assuming the low-frequency broadband variability of the wind is similar to Chl *a*, the S<sub>sa</sub> harmonic constants are judged as not reliable and are thus omitted from the presentation in the main text.</figcaption></figure>

effect of data gaps, then, can be regarded as redistributing
variance of the gap-free least-squares residual throughout the
frequency domain, according to the Fourier transform of the
data indicator function, *w*.

Figure A1a shows the count of valid data pixels versus the
day of year, averaged over the 2002–2022 time period. There
is a notable increase in the number of cloud-free pixels from
late March through August during the dry season in the re-
gion of the northern Australia and Indonesia monsoon (Mu-
rakami and Matsumoto, 1994), also called the dry southeast
monsoon (Aldrian and Susanto, 2003). Indeed, the power
spectrum of *w*<sub>c</sub> averaged over all the grid cells of the Chl *a*
dataset exhibits peaks at the S<sub>a</sub> and S<sub>sa</sub> periods, the larger

of which contains about 15% as much power as the zeroth
frequency (Fig. A1b). The power at the MSf frequency of the
spring–neap cycle is about 0.3% that at the zeroth frequency.
While the latter seems negligible, the expected crosstalk is
proportional to the square root of this quantity (about 0.05)
and the MSf amplitude is so much smaller than the mean
that this crosstalk could be significant. Thus, as mentioned in
the main text, MSf harmonic constants are considered 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<sub>f</sub> peak in the data indicator spectrum
suggests that the data flagging algorithm is influenced by
tides in this region. It is conceivable that fortnightly modu-
lations 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 as turbidity, might change the outlier flagging in the GSM algorithm.

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 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, 5.2, 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, 0.6, and 0.15 m s⁻¹, for the mean, $S_a$, and $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 in the mean wind due to the data gaps, $(\hat{x}_1/x_1 - 1) \times 100\%$, where $\hat{x}_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 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.

The crosstalk from the mean and $S_{sa}$ into the $S_a$ signal is illustrated in Fig. A3a and b; note the larger range of color scales compared to Fig. 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 negligible, 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, Fig. A3c and d indicate that the data gaps are likely to make the $S_{sa}$ variability of Chl *a* impossible to interpret. Indeed, an analysis of the formal (random) error in

the $S_{sa}$ harmonic constants suggests that the signal-to-noise ratio is below 1 in about half of the domain (not shown), and Fig. 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 main 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 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 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.

*Code and data availability.* Julia-language scripts for harmonic analysis and plotting the MODIS data are available from the lead author's software repository (https://ingria.ceoas.oregonstate.edu/fossil/CHL, Zaron, 2023).

*Author contributions.* EDZ drafted the paper and performed the data analysis. TAC provided the Chl *a* dataset and contributed to writing and editing the paper. AKL conceived of the objective of the research, guided the selection of significant geographic regions for study, and contributed to writing the paper.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

*Acknowledgements.* Support for Edward D. Zaron was provided by the NASA Ocean Surface Topography Science Team, award #80NSSC21K1189. The GlobColour daily MODIS product was downloaded from https://modis.gsfc.nasa.gov (last access: 11 January 2023) and was developed, validated, and distributed by ACRI-ST, France (http://globcolour.info, last access: 11 January 2023). The NCEP/DOE Reanalysis II data used in the Appendix were provided by the NOAA PSL, Boulder, Colorado, USA, from their website at https://psl.noaa.gov (last access: 11 January 2023). The TOPEX and Jason mission data used in Fig. 7 were obtained from the Radar Altimeter Database System (http://rads.tudelft.nl, last access)
cess: 11 January 2023). The comments of Raden Dwi Susanto and an anonymous reviewer, which helped to greatly improve the quality of this paper, are acknowledged and appreciated.

*Financial support.* This research has been supported by the National Aeronautics and Space Administration (grant no. 80NSSC21K1189).

*Review statement.* This paper was edited by Jochen Wollschlaeger and reviewed by Raden Dwi Susanto and one anonymous referee.

## References

Aldrian, E. and Susanto, R. D.: Identification of three dominant rainfall regions within Indonesia and their relationship to sea surface temperature, Int. J. Climatol., 23, 1435–1452, https://doi.org/10.1002/joc.950, 2003.

Apel, J. R., Holbrook, J. R., Liu, A. K., and Tsai, J. J.: The Sulu Sea Internal Soliton Experiment, J. Phys. Oceanogr., 15, 1625–1651, 1985.

Atmadipoera, A., Koch-Larrouy, A., Madec, G., Grelet, J., Baurend, F., Jaya, I., and Dadou, I.: Part I: hydrological properties within the eastern Indonesian Throughflow region during the INDOMIX experiment, Deep-Sea Res. Pt. I, 182, 103735, https://doi.org/10.1016/j.dsr.2022.103735, 2022.

Atmadipoera, A. S. and Suteja, Y.: Deep water masses exchange induced by internal tidal waves in Ombai Strait, IOP Conf. Ser.-Earth Environ. Sci., 176, 012017, https://doi.org/10.1088/1755-1315/176/1/012017, 2018.

Bouruet-Aubertot, P., Cuypers, Y., Ferron, B., Dause, D., Ménage, O., Atmadipoera, A., and Jaya, I.: Contrasted turbulence intensities in the Indonesian Throughflow: a challenge for parameterizing energy dissipation rate, Ocean Dynam., 68, 779–800, 2018.

Clement, A. C., Seager, R., and Murtugudde, R.: Why are there tropical warm pools?, J. Climate, 18, 5294–5311, 2005.

Devlin, A. T., Zaron, E. D., Jay, D. A., Talke, S. A., and Pan, J.: Seasonality of tides in Southeast Asian waters, J. Phys. Oceanogr., 48, 1169–1190, 2018.

Egbert, G. D. and Erofeeva, S. Y.: Efficient inverse modeling of barotropic ocean tides, J. Atmos. Ocean. Tech., 19, 183–204, 2002.

Egbert, G. D. and Erofeeva, S. Y.: An approach to empirical mapping of incoherent internal tides with altimetry data, Geophys. Res. Lett., 48, e2021GL095863, https://doi.org/10.1029/2021GL095863, 2021.

Foreman, M. G., Cherniawsky, J. Y., and Ballantyne, V. A.: Versatile Harmonic Tidal Analysis: Improvements and Applications, J. Atmos. Ocean. Tech., 26, 806–817, 2009.

Garcia, C. A. E. and Garcia, V. M. T.: Variability of chlorophyll-*a* from ocean color images in the La Plata continental shelf region, Cont. Shelf Res., 28, 1568–1578, https://doi.org/10.1016/j.csr.2007.08.010, 2008.

Gutknecht, E., Reffray, G., Gehlen, M., Triyulianti, I., Berlianty, D., and Gaspar, P.: Evaluation of an operational ocean model configuration at 1/12° spatial resolution for the Indonesian seas (NEMO2.3/INDO12) – Part 2: Biogeochemistry, Geosci. Model

Dev., 9, 1523–1543, https://doi.org/10.5194/gmd-9-1523-2016, 2016.

Hermansyah, H., Atmadipoera, A. S., Prartono, T., Jaya, I., and Syamsudin, F.: Energetics of internal tides over the Sangihe-Talaud ridge – Sulawesi sea, J. Phys.-Conf. Ser., 1341, 082001, https://doi.org/10.1088/1742-6596/1341/8/082001, 2019.

Kanamitsu, M., Ebisuzaki, W., Woollen, J., Yang, S.-K., Hnilo, J., Fiorino, M., and Potter, G. L.: NCEP-DOE AMIP-II Reanalysis (R-2), B. Am. Meteorol. Soc., 83, 1631–1643, 2002.

Kida, S. and Wijffels, S.: The impact of the Indonesian Throughflow and tidal mixing on the summertime sea surface temperature in the western Indonesian Seas, J. Geophys. Res., 117, C09007, https://doi.org/10.1029/2012JC008162, 2012.

Klymak, J. M., Pinkel, R., and Rainville, L.: Direct breaking of the internal tide near topography: Kaena Ridge, Hawaii, J. Phys. Oceanogr., 38, 380–399, 2008.

Koch-Larrouy, A., Madec, G., Bouruet-Aubertot, P., Gerkema, T., Bessieres, L., and Molcard, R.: On the transformation of Pacific Water into Indonesian Throughflow Water by internal tidal mixing, Geophys. Res. Lett., 34, L04604, https://doi.org/10.1029/2006GL028405, 2007.

Koch-Larrouy, A., Madec, G., Iudicone, D., Atmadipoera, A., and Molcard, R.: Physical processes contributing to the water mass transformation of the Indonesian Throughflow, Ocean. Dynam., 58, 275–288, https://doi.org/10.1007/s10236-008-0154-5, 2008.

Koch-Larrouy, A., Lengaigne, M., Terray, P., Madec, G., and Masson, S.: Tidal mixing in the Indonesian Seas and its effect on the tropical climate system, Clim. Dynam., 34, 891–904, https://doi.org/10.1007/s00382-009-0642-4, 2010.

Laws, E. A.: Evaluation of in situ phytoplankton growth rates: a synthesis of data from varied approaches, Annu. Rev. Mar. Sci., 5, 247–268, https://doi.org/10.1146/annurev-marine-121211-172258, 2013.

Levine, M. D. and Boyd, T. J.: Tidally-forced internal waves and overturns observed on a slope: Results from the HOME survey component, J. Phys. Oceanogr., 36, 1184–1201, 2006.

Mandal, S., Behera, N., Gangopadhyay, A., Susanto, R. D., and Pandey, P. C.: Evidence of a chlorophyll “tongue” in the Malacca Strait from satellite observations, J. Marine Syst., 223, 103610, https://doi.org/10.1016/j.jmarsys.2021.103610, 2021.

Mandal, S., Susanto, R. D., and Ramakrishnan, B.: Dynamical Factors Modulating Surface Chlorophyll-*a* Variability along South Java Coast, Remote Sens., 14, 1745, https://doi.org/10.3390/rs14071745, 2022.

Maritorena, S. and Siegel, D. A.: Consistent merging of satellite ocean color data sets using a bio-optical model, Remote Sens. Environ., 94, 429–440, 2005.

Maritorena, S., d'Andon, O. H. F., Mangin, A., and Siegel, D. A.: Merged satellite ocean color data products using a bio-optical model: Characteristics, benefits and issues, Remote Sens. Environ., 114, 1791–1804, https://doi.org/10.1016/j.rse.2010.04.002, 2010.

Matte, P., Jay, D. A., and Zaron, E. D.: Adaptation of Classical Tidal Harmonic Analysis to Non-Stationary Tides, with Application to River Tides, J. Atmos. Ocean. Tech., 30, 569–589, 2013.

Murakami, T. and Matsumoto, J.: Summer monsoon over the Asian Continent and western North Pacific, J. Meteorol. Soc. Jpn., 72, 719–745, 1994.
Nagai, T. and Hibiya, T.: Internal tides and associated vertical mixing in the Indonesian Archipelago, J. Geophys. Res., 120, 3373–3390, https://doi.org/10.1002/2014JC010592, 2015.

Nagai, T., Hibiya, T., and Bouruet-Aubertot, P.: Nonhydrostatic simulations of tide-induced mixing in the Halmahera Sea: A possible role in the transformation of the Indonesian Throughflow waters, J. Geophys. Res., 122, 8933–8943, https://doi.org/10.1002/2017JC013381, 2017.

Nugroho, D.: The tides in a general circulation model in the Indonesian Seas, PhD thesis, Paul Sabatier University, Toulouse, France, https://theses.hal.science/tel-01897523 (last access: 11 January 2023), 2017.

Nugroho, D., Koch-Larrouy, A., Gaspard, P., Lyard, F., RefrAY, G., and Tranchant, B.: Modelling explicit tides in the Indonesian seas: An important process for surface sea water properties, Mar. Pollut. Bull., 131, 7–18, https://doi.org/10.1016/j.marpolbul.2017.06.033, 2018.

Osborne, A. R. and Burch, T. L.: Internal solitons in the Andaman Sea, Science, 208, 451–460, https://doi.org/10.1126/science.208.4443.451, 1980.

Palacios, D. M.: Seasonal patterns of sea-surface temperature and ocean color around the Galapagos: regional and local influences, Deep-Sea Res., 51, 43–57, 2004.

Ray, R. D. and Susanto, R. D.: Tidal mixing signatures in the Indonesian Seas from high-resolution sea surface temperature data, Geophys. Res. Lett., 43, 8115–8123, https://doi.org/10.1002/2016GL069485, 2016.

Ray, R. D. and Susanto, R. D.: A fortnightly atmospheric “tide” at Bali caused by oceanic tidal mixing in Lombok Strait, Geosci. Lett., 6, 6, https://doi.org/10.1186/s40562-019-0135-1, 2019.

Robertson, R. and Ffield, A.: Baroclinic tides in the Indonesian Seas: Tidal fields and comparisons to observations, J. Geophys. Res., 113, C07031, https://doi.org/10.1029/2007JC004677, 2008.

Ropelewski, C. F. and Halpert, M. S.: Global and regional scale precipitation patterns associated with El Nino/Southern Oscillation, Mon. Weather Rev., 115, 1606–1626, 1987.

Shi, W., Wang, M., and Jiang, L.: Spring-neap tidal effects on satellite ocean color observations in the Bohai Sea, Yellow Sea, and East China Sea, J. Geophys. Res., 116, C12032, https://doi.org/10.1029/2011JC007234, 2011.

Siswanto, E., Horii, T., Iskandar, I., Gaol, J. L., Setiawan, R. Y., and Susanto, R. D.: Impacts of climate changes on the phytoplankton biomass of the Indonesian Maritime Continent, J. Marine Syst., 212, 103451, https://doi.org/10.1016/j.jmarsys.2020.103451, 2020.

Sprintall, J., Gordon, A. L., Koch-Larrouy, A., Lee, T., Potemra, J. T., Pujiana, K., and Wijffels, S.: The Indonesian Seas and their impact on the Coupled Ocean Climate System, Nat. Geosci., 7, 487–492, https://doi.org/10.1038/NGEO2188, 2014.

Sprintall, J., Gordon, A. L., Wijffels, S. E., Feng, M., Hu, S., Koch-Larrouy, A., Phillips, H., Nugroho, D., Napitu, A., Pujiana, K., Susanto, R. D., Sloyan, B., Pena-Molino, B., Yuan, D., Riama, N. F., Siswanto, S., Kuswardani, A., Arifin, Z., Wahyudi, A. J., Zhou, H., Nagai, T., Ansong, J. K., Bourdalle-Badié, R., Chanut, J., Lyard, F., Arbic, B. K., Ramdhani, A., and Setiawan, A.: Detecting change in the Indonesian Seas, Front. Mar. Sci., 6, 257, https://doi.org/10.3389/fmars.2019.00257, 2019.

Susanto, R. D. and Ray, R. D.: Seasonal and interannual variability of tidal mixing signatures in Indonesian seas from high-resolution sea surface temperature, Remote Sens., 14, 1934, https://doi.org/10.3390/rs14081934, 2022.

Susanto, R. D., Mitnik, L., and Zheng, Q.: Ocean internal waves observed in the Lombok strait, Oceanography, 18, 80–87, https://doi.org/10.5670/oceanog.2005.08, 2005.

Susanto, R. D., Moore II, T. S., and Marra, J.: Ocean color variability in the Indonesian Seas during the SeaWiFS era, Geochem. Geophy. Geosy., 7, Q05021, https://doi.org/10.1029/2005GC001009, 2006.

Wirasatriya, A., Susanto, R. D., Kunarso, K., Jalil, A. R., Ramdani, F., and Puryajati, A. D.: Northwest monsoon upwelling within the Indonesian Seas, Int. J. Remote. Sens., 42, 5437–5458, https://doi.org/10.1080/01431161.2021.1918790, 2021.

Xing, Q., Yu, H., Yu, H., Wang, H., Ito, S., and Yuan, C.: Evaluating the spring-neap tidal effects on chlorophyll-*a* variations based on the geostationary satellite, Front. Mar. Sci., 8, 758538, https://doi.org/10.3389/fmars.2021.758538, 2021.

Yang, M., Khan, F. A., Tian, H., and Liu, Q.: Analysis of the monthly and spring-neap tidal variability of satellite chlorophyll-*a* and total suspended matter in a turbid coastal ocean using the DINEOF method, Remote Sensing, 13, 632, https://doi.org/10.3390/rs13040632, 2021.

Zaron, E. D.: Baroclinic tidal sea level from exact-repeat mission altimetry, J. Phys. Oceanogr., 49, 193–210, 2019.

Zaron, E. D.: Scripts related to the paper, Zaron Oregon State University [code] https://ingria.ceoas.oregonstate.edu/fossil/CHL, last access: 11 January 2023.

Zaron, E. D., Musgrave, R., and Egbert, G. D.: Baroclinic tidal energetics inferred from satellite altimetry, J. Phys. Oceanogr., 52, 1015–1032, 2022.

Zhang, X., Li, X., and Zhang, T.: Characteristics and generations of internal wave in the Sulu Sea inferred from optical satellite images, J. Ocean. Limnol., 38, 1435–1444, https://doi.org/10.1007/s00343-020-0046-1, 2020.

Zhao, X., Xu, Z., Feng, M., Li, Q., Zhang, P., You, J., Gao, S., and Yin, B.: Satellite investigation of semidiurnal internal tides in the Sulu-Sulawesi Seas, Remote Sensing, 3, 2530, https://doi.org/10.3390/rs13132530, 2021.

Zhao, Z. and Alford, M. H.: Source and propagation of internal solitary waves in the northeastern South China Sea, J. Geophys. Res., 111, C11012, https://doi.org/10.1029/2006JC003644, 2006.