Seasonal and interannual variations of MODIS Aqua chlorophyll-a (2003–2017) in the Upper Gulf of Thailand influenced by Asian monsoons

Seasonal and interannual variations of chlorophyll-a (chl-a) in the upper Gulf of Thailand (uGoT) were obtained using new regionally tuned algorithms applied to Moderate Resolution Imaging Spectroradiometer-Aqua. This long time-series (2003–2017) data were analyzed in the context of variations in environmental conditions associated with the Southeast Asian Monsoon. Chl-a distribution patterns were distinct for the non-monsoon (NOM), southwest-monsoon (SWM), and northeast-monsoon (NEM) seasons. During the SWM/NEM, high/low chl-a concentrations were associated with high/low precipitation and river discharge. During the NOM chl-a concentrations were generally low, because of low precipitation. In general, chl-a variability was tightly coupled to discharge from the Chao Phraya and Tha Chin rivers. Chl-a concentrations were generally higher in the north, but chl-a accumulation in the east/west of the uGoT could be linked to piling of freshwater to the east/west during the SWM/NEM caused by changes in wind direction and the reversal of currents. Interannual changes in chl-a were attributed to El Niño-Southern Oscillation (ENSO) rather than Indian Ocean Dipole (IOD) driven changes in precipitation, river discharge, and wind patterns. During the SWM, positive/negative chl-a anomalies coincided with high/low precipitation and river discharge during La Niña/El Niño. During the NEM, positive/negative chl-a anomaly coincided with high/low river discharge and strong/weak wind during La Niña/El Niño. Meanwhile, during NOM, positive chl-a anomaly could be attributed to anomalous high wind speed and precipitation during El Niño.


Introduction
The influence of the seasonally reversing Asian monsoons on the distribution of phytoplankton in the upper euphotic column has been reported for several regions within southeast Asia, including the South China Sea (e.g., Shen et al. 2018;Tang et al. 2003;Yuan-Jian et al. 2012), the eastern and western coasts of Sabah (e.g., Abdul-Hadi et al. 2013), the Strait of Malacca (e.g., Siswanto and Tanaka 2014), and the south-eastern Arabian Sea (e.g., Shafeeque et al. 2019). Several of these studies have shown that the strength of monsoonal winds can enhance phytoplankton biomass by stimulating oceanic physical processes (e.g., coastal upwelling and vertical mixing), that facilitate the transport of nutrient-rich waters from the deep and enhance primary productivity within the entire euphotic layer (Pinet 2014). The seasonal cycle of the monsoons also impacts rainfall patterns in southeast Asia, and consequently phytoplankton production via the influx of riverine nutrients (Simpson and Sharples 2012).
Often, the annual cycle of the Asian monsoons can be interrupted by global scale climate phenomena, such as El Niño-Southern Oscillation (ENSO) and the Indian Ocean Dipole (IOD). These climatic events contribute to variabilities in monsoonal wind strength and precipitation patterns that can have cascading impacts on regional phytoplankton growth and bloom formation in the region. Within the South China Sea and the region south of Malacca Strait for instance, reductions in phytoplankton biomass during an El Niño event, were attributed to weaker than normal winds and significant reductions in wind-induced upwelling (Zhao and Tang 2007;Yuan-Jian et al. 2012;Siswanto and Tanaka 2014;Zhao et al. 2018). Unfortunately, there is no comparable information for the Gulf of Thailand which also comes under the influence of both the Asian monsoons and global scale climate phenomenon.
Our present study is focused on the upper Gulf of Thailand (hereafter uGoT, Fig. 1), the largest (approximately 100 × 100 km 2 ) shallow (average depth = 15 m), semienclosed basin of Thailand. Eutrophication is often reported in the coastal area due to high nutrient inputs from the major rivers namely Phetchaburi, Mae Klong, Tha Chin, Chao Phraya, and Bang Pakong Rivers, located from the west to east along the upper coast, that flow into the uGoT. These rivers, especially the Chao Phraya River, carry with them large volumes of freshwater along with land-based contaminants and nutrients from agriculture, aquaculture, and human sewage into the uGoT (Wattayakorn 2006). This freshwater discharge along with its constituents have a large influence on the seawater properties of uGoT, such as salinity, nutrients, turbidity, etc. (Buranapratheprat et al. 2002;Kobayashi et al. 2010a, b). Furthermore, phytoplankton blooms, which generally follow peaks of freshwater discharge, have often been reported several times a year at different locations, some of which cause mass fish mortality. Although it is assumed that the bloom is related to anthropogenic impacts, little is known about the evolution of these blooms or their relationship with precipitation-driven changes in freshwater discharge or to seasonal to interannual Asian Monsoon variabilities.
The climate of Thailand is dominated by the southwest (SWM) and northeast (NEM) monsoons, that are responsible for fluctuations in precipitation (Thai Meteorological Fig. 1 Map of the uGoT and sampling stations. Circles, stars, pentagons, triangles, and squares on the bathymetric map of the uGoT represent data collected during four different surveys, KB, KU, CU, and AK listed in Table 1. Hollow squares of W, CE, C, CW, E, and O represent the average areas (15 × 15 pixels) for analysis in this study. The region north of 12° 35.4′ (dotted line) is referred to as the whole area used for obtaining spatial averages for this study 1 3 Department 2014). The SWM (mid-May to mid-October) brings a stream of warm and moist air from the Indian Ocean, that drives the rainy season. On the other hand, the NEM (mid-October to mid-February) brings cold and dry air from the China mainland to Thailand and causes the weather to be cold and dry, referred to as the winter of Thailand. The period in between (mid-February to mid-May) that is not influenced by the monsoon is considered as the summer season. The monsoonal winds, combined with the unique topography of the uGoT, generate two specific hydrographic circulation patterns along the northern coast; i.e., eastward during the SWM and westward transport during the NEM (Buranapratheprat et al. 2002(Buranapratheprat et al. , 2006. It is believed that the variations in monsoonal rainfall, river runoff and monsooninduced changes in circulation, can therefore, potentially influence seasonal phytoplankton variability in the uGoT. Using Level-2 (L2) Medium Resolution Imaging Spectrometer Instrument (MERIS) data with a local in-water chl-a algorithm of Matsumura et al. (2006) and  revealed seasonal distribution of surface phytoplankton in the uGoT with six daily data during [2003][2004][2005]. The high chlorophyll-a (chl-a) in the northeastern/western coastal areas during the SWM/NEM season corresponded with the monsoon winds. The result was consistent with modeling results of chl-a distribution, simulated from monthly local wind, water movement, and river discharge data (Buranapratheprat et al. 2008b). High phytoplankton biomass was also observed in relation to nutrient (dissolved inorganic nitrogen and phosphorus) inputs into the uGoT during the rainy season in estuaries near the northwest coast (Thongdonphum et al. 2011(Thongdonphum et al. , 2014. In most cases, high phytoplankton biomass was caused by high cell densities of green Noctiluca scintillans, dinoflagellate Ceratium furca, and diatoms, frequently in the northeastern coastal area during the SWM Sriwoon et al. 2008). Although variations in phytoplankton biomass seem to be related to seasonal monsoonal winds and to precipitation and circulation patterns, these relationships can best be described as cursory, lacking details of mechanisms that contribute to this seasonality.
Impacts of ENSO and IOD on Asian monsoons is well known that affect the monsoonal rainfall variability over the mainland of Thailand (Chaowiwat and Koontanakulvong 2016). Low/high rainfall is during El Niño/La Niña (Singhrattna et al. 2005;Kirtphaiboon et al. 2014) and strong positive/negative IOD (Chansaengkrachang et al. 2011). Recently, Kim et al. (2020) revealed that the rainfall variability in South China and the Indochina Peninsula including Thailand, is amplified by combined and/or independent effects of the ENSO and IOD. The ENSO-induced variability in precipitation, resulted in the great flood of Thailand in 2011 due to high precipitation during strong La Niña (Gale and Saunders 2013). In contrast, during El Niño, the positive phase of the IOD (pIOD), and the ENSO neutral condition, dramatic declines in precipitation and droughts have been reported (Kim et al. 2020;Power et al. 2020).
Although the impacts of ENSO and IOD on local weather and rainfall patterns are well known, the downstream impacts of precipitation on the uGoT ecosystem are not known. One of the only reports by Intacharoen et al. (2018) available describes a reduction of high chl-a area in the uGoT by the standard MODerate-resolution Imaging Spectroradiometer-Aqua (MODIS) product during SWM of the 2011 when Thailand experienced significant flooding. Although global climate variability (e.g., ENSO and IOD) is expected to affect plankton dynamics, no studies have reported long-term interannual variability and the climate variability effects on phytoplankton in Thai Waters. One of major impediments has been the paucity of field observations over the space and time scales required to understand the response of the uGoT lower trophic ecosystem to seasonal and interannual changes in monsoonal cycles.
To fill this knowledge gap in our understanding of the seasonal to inter-annual variability in chl-a distribution as well as the formation of blooms in the uGoT in response to Asian monsoon-driven changes in circulation and precipitation, we relied on a time-series of MODIS ocean color data sets between 2003 and 2017. Prior to doing so, we improved the accuracy of the standard National Aeronautics and Space Administration (NASA) chl-a product of MODIS for the uGoT using a new empirical algorithm for the uGoT that was developed using in situ optical and chl-a datasets. After this step, we composited the improved MODIS chl-a into a coherent time-series, to investigate seasonal and interannual chl-a variations. Finally, we analyzed these improved chl-a satellite time-series dataset in the context of monsoon-driven variability in sea surface winds, precipitation and river discharge, as well as with the variations in the periods of ENSO and IOD.

In situ optical and chl-a data
In situ chl-a and remote sensing reflectance ( R rs (λ)) datasets were collected at several locations in the uGoT and used in this study (Table 1, Fig. 1). In situ chl-a data measured by the fluorescence method (Strickland and Parsons 1972;Suzuki and Ishimaru 1990) varied between 0.04 and 84.5 mg m −3 . The R rs (λ) were obtained by using a Profiling Reflectance Radiometer model 600 (PRR-600) for 2003 and 2004 measuring downward irradiance, upward radiance, and sky irradiance just below the sea surface (Matsumura et al. 2006). Starting in 2009, water-leaving radiances and downwelling irradiance at the sea surface (Kobayashi et al. 2011) and just above the sea surface (Kobayashi et al. 2010a, b) were undertaken using a hyperspectral radiometer RAMSES (TriOS, Germany). In situ R rs and chl-a were used to develop a new uGoT specific 'switching' algorithm for chl-a, for use with R rs , derived from MODIS (see in Supplementary Appendix 1). The global chl-a product for MODIS is based on the OCx chl-a, which is a fourth-order polynomial function of R rs blue-to-green (443-555 nm) (O'Reilly et al. 1998). A subset of the in situ chl-a listed in Table 1 data was used for validating the chl-a.
MODIS data were matched up to in situ data using WIM Automation Module (WAM; http:// www. wimso ft. com/) which uses the nearest location and time of the sampling point for validating satellite products against in situ data. Satellite match-ups were based on averages 3 × 3 pixels from the center of sampling point and data obtained within ± 24 h time difference from the in situ sampling. MODIS pixels (1) Chl a = 10 (0.2424−2.7423R+1.8017R 2 +0.0015R 3 −1.2280R 4 ) , R = log 10 max R rs 443 R rs 547 , R rs 488 R rs 547 .
selected for the match-ups were quality controlled by the masks of LAND, HIGLINT, HILT, CLDICE, HISOLZEN, LOWLW, CHLFAIL, and NAVFAIL (https:// www. ocean color. gsfc. nasa. gov/ atbd/ ocl2fl ags/). We used the matched-up data to verify the accuracy of the MODIS R rs values by comparing with the in situ R rs data before applying chl-a algorithm to MODIS data to examine the performance of the MODIS atmospheric correction algorithm. Also, the matched-up R rs and chl-a data were used validate the performance of our switching chl-a algorithm (Eqs. 2, 3) to evaluate the reliability of the new chl-a algorithm and compared with the standard MODIS chl-a algorithm (Eq. 1).
Statistical analyses (Campbell and O'Reilly 2006) of the percent root-mean-square error (RMSE, %) and the percent bias (Bias, %) were used to evaluate the performance of R rs data from satellite and in situ sensors. The statistical equations are as the following: (2) Chl a = 10 (−3.5076R−0.1132) for R rs 547 =< 0.005sr −1 (3) Chl a = 10 (−4.0317R+0.1734) for R rs 547 > 0.005sr −1 where X sat and X obs are satellite R rs and in situ R rs , respectively, and N is the number of data.

Satellite data processing
The new chl-a algorithms (Eqs. 2, 3) were used to generate daily L2 MODIS data for the period from 2003 to 2017 after geometric correction. The masks mentioned in 2.2 were utilized to control the pixel quality. Because of the limitation of the dataset used for algorithm development, the maximum chl-a was set at 100 mg m −3 in each pixel. The daily MODIS chl-a after log-transformation, were binned to generate monthly chl-a composites, for the data from 2003 to 2017. The number of the daily images was shown in Supplementary Appendix 2. Monthly climatological chl-a and seasonal percentage chl-a anomalies were calculated using these regional tuned products. The seasonal percentage chl-a anomaly was difference between the seasonal composite of each year and seasonal climatological composite normalized by the seasonal climatology (e.g., Kahru et al. 2010). The seasonal composites were calculated for three seasons; (1) non-monsoon (NOM; February-April), (2) southwest-monsoon (SWM; May-September), and (3) northeast-monsoon (NEM; October-January of the following year).
Finally, chl-a averages of the uGoT were calculated after log-transformation and then reversed to chl-a values for the whole uGoT (12°35.40′ N-13°32.59′ N and 99°57.26′ E-100°58.27′ E). They were then spatially averaged (15 × 15 pixels) at five locations near the northern coast and for one offshore location to the south (Fig. 1).

Other related environmental variables
River discharge data (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) for the principal rivers along the coast were provided by Royal Irrigation Department of Thailand. Because there are no hydrographic stations at the river mouths, the river discharge data of the central and east coasts was retrieved from upstream stations, C.2 (Ban Phai Lom; Nakorn Sawan, representing Chao Phraya and Tha Chin rivers) and Kgt.3 (Ban Kabin Buri; Prachin Buri, representing Bang Pakong River), respectively. The total of B.10 (Ban Tha Yang, Phetchaburi, representing Phetchaburi River) and K.35A (Ban Nong Bua; Kanchanaburi, representing Mae Klong River) stations were used as representative river discharge from the west coast. Total river discharge was calculated by the sum of monthly river discharge from all rivers.
The bi-monthly NCEP-NCAR Multivariate ENSO Index (MEI) used in our study comes from the Webberweather website (https:// www. webbe rweat her. com/ multi varia teenso-index. html). It is calculated based on the six observed atmospheric and oceanic variables over the Equatorial Pacific Ocean (i.e., sea surface temperature, surface air temperature, sea level pressure, zonal winds, surface outgoing longwave radiation, and precipitable water). The MEI was used to represent the phase of ENSO. The three-month running mean Dipole Mode Index (DMI, monthly) was obtained from the Japanese Agency for Marine-Earth Science and Technology website (http:// www. jamst ec. go. jp/ virtu alear th/ gener al/ en/ index. html), calculated as the differences in the sea surface temperatures between the west and southeast of Equatorial Indian Ocean, was used to present the phase of IOD. We further averaged the MEI and DMI values in each season for investigating their role in interannual chl-a variation.
Multivariate correlation analysis was performed to investigate relationships among the monthly climatological and seasonal percentage anomaly data of chl-a and other variables (i.e., river discharge, winds, and precipitation) as well as the seasonal ENSO and IOD indexes. Pearson's correlation (r) with p value < 0.05 was used to indicate the significance of the correlation.

Improvement of satellite MODIS chl-a
Comparisons of in situ and MODIS R rs 488 and R rs 547 estimates used for the OCx (Eq. 1) and in our switching algorithms (Eqs. 2, 3) are shown in Fig. 2a, b. Both satellite R rs 547 and R rs 488 were underestimated, and this implied the problem of atmospheric corrections. However, the errors were reduced when used as the ratio of R rs 488∕R rs 547 (R 2 = 0.93, RMS = 16%, bias = − 3%) (Fig. 2c). Overestimation of the standard satellite chl-a algorithm (Eq. 1), which has been reported in and Intacharoen et al. (2018), due to underestimates in R rs 547 and R rs 488 were also observed in our study (Fig. 3a). With the local log-log R rs 488∕R rs 547 algorithm, and switching at R rs 547 = 0.005 (Eqs. 2, 3), the accuracy of MODIS derived chl-a showed a marked improvement (Fig. 3b) over those estimated by the standard OC3M algorithm (R 2 from 0.77 to 0.83, bias from 0.174 to 0.058). These results indicate that the local log-log R rs 488∕R rs 547 algorithm with switching at R rs 547 = 0.005 is more reliable than the standard MODIS algorithm for estimating chl-a in the uGoT.

uGoT wide chl-a variation
The 15-year monthly climatological data of MODIS chla revealed three surface chl-a distribution patterns corresponding with surface wind direction (Fig. 4). During NOM (February-April), chl-a was uniformly distributed along the entire coast, north of 13°12′ N, and in general higher than 10 mg m −3 . Most of this region is shallow and under the influence of riverine discharge and southerly wind (4.1-5.1 m s −1 ). During SWM (May to September), a period of persistent southwesterly winds (4.8-6.1 m s −1 ), chl-a was high in the northeastern area and characterized by a broader offshore spread than during NOM. In October, the shift in the high chl-a from the east to west coincided with the change of the southwesterly wind to the southward. The high chl-a area, however, tended to reduce with the dominance of northeasterly winds (2.4-5.7 m s −1 ) during the rest of NEM (November-January). Based on these chl-a distribution patterns, we therefore divided the seasonality of chl-a into three seasons (i.e., NOM, SWM, and NEM) as mentioned in Sect. 2.3.
Area averaged chl-a for the whole uGoT increased from January to September i.e., through NOM and SWM and decreased during NEM (Fig. 5b). Chl-a was high during SWM (June to September, 3.2 mg m −3 ) while the minimum was in January (1.7 mg m −3 ).
The total river discharge ( Fig. 5b) was low and nearly constant (524-608 m 3 s −1 ) during NOM, and then increased and decreased during SWM and NEM, respectively, with the peak in October (2427 m 3 s −1 ). Meanwhile, precipitation was low and gradually increase during NOM, and steadily high during SWM. It was also high in September and October (> 230 mm) and decreased during NEM with the lowest in December (11 mm) (Fig. 5a).
In general, the area average chl-a for uGoT increased during SWM corresponding with the increase in total river discharge as well as high precipitation and u-wind speed, while their variations reversed during NEM (Fig. 5a, b). During NOM, the increase in chl-a coincided with the increase in Comparison of in situ and estimated chl-a by the standard OC3M in-water chl-a algorithm (a) and the local switching chl-a algorithm (b) using satellite R rs . The dashed line is a one-to-one line. Solid lines indicate the ratio of estimated chl-a to in situ chl-a is 2 and 0.5. Cross and triangle symbols indicate R rs 547 = < 0.005 and R rs 547 > 0.005, respectively precipitation during low river discharge. Correlation analysis (Table 2) suggested that the area averaged chl-a was significantly correlated with both total river discharge (r = 0.52, p < 0.0001) and precipitation (r = 0.58, p < 0.0001), as well as u-wind speed (r = 0.40, p < 0.0001).

Regional variation
In the west (Fig. 5c), chl-a was nearly constant (8.4-11.7 mg m −3 ) throughout the year. Chl-a in this area was the highest of the three coastal regions during NEM (November-January) when chl-a of other areas were low, while it was the lower than the other areas during NOM and SWM. In the central area (Fig. 5d), the chl-a was higher than 10 mg m −3 from NOM to early NEM in October, and highest in the three coastal regions during the NOM. In the east (Fig. 5e), chl-a was the highest (18 mg m −3 ) in the three coastal regions during SWM, and rapidly decreased during NEM to the minimum in January (5.1 mg m −3 ). In the offshore area (Fig. 5d), chl-a was lower (< 3 mg m −3 ) than coastal area and revealed a similar seasonal variability over the whole area.
River discharge from the west via Mae Klong and Phetchaburi rivers was low with little variation (118-204 m 3 s −1 ) throughout the year (Fig. 5c). River discharge via Tha Chin and Chao Phraya rivers in the central area was the highest in three regions (337-1882 m 3 s −1 ) and corresponded 74% of the total river discharge (Fig. 5d). Variability of the total river discharge was largely contributed by variations in discharge from these two rivers in the central area. Variation of the river discharge via Bang Pakong River in the east was similar to the variations of rivers in the central area, but the amount was lower even with the peak (352 m 3 s −1 ) in September (Fig. 5e). Precipitation showed significant correlations with the river discharges except in the west, and it was more correlated to the regional chl-a than the river discharge ( Fig. 5; Table 2).
Correlation between the regional chl-a and both regional and total river discharge was significant, except the chl-a in the central, and similar to the whole area in the offshore area (Table 2). In the west (Fig. 5c), both chl-a and river discharge were lowest in the coastal areas, and seasonal variations were also small. In the central area (Fig. 5d), the river discharge was highest, and chl-a was also high. The chl-a was constantly high during NOM and SWM, while peak of the river discharge was in late SWM to early NEM. In the east (Fig. 5e), chl-a was high during SWM, although the river discharge was low and peaked in late SWM to early NEM. This discrepancy may be indicative of the importance of wind in as a driver of chl-a in the coastal area.
Comparison between chl-a along the coast and wind showed the monthly peak of chl-a corresponding with the  (Figs. 5a, 6). Chl-a in the central area between the mouths of the Tha Chin and Chao Phraya rivers was persistently higher especially when the southerly wind steadily blew during NOM. The chl-a minimum, which was in March in the west, occurred when the u-wind velocity was nearly zero. The chl-a peak moved to the east of center during SWM, associated with the strong southwesterly wind. The high chl-a from June to September at the east of the Chao Phraya River coincided with the wind speed maximum. In October, the chl-a peak returned to the center, while the wind speed dropped and changed direction from northeast to south. From November, with the onset of NEM, the high chl-a waters moved to the west simultaneously with the reduction of chl-a in the east, especially from the mouth of the Chao Phraya River. Correlation analysis also showed significant correlation of the chl-a in the central and east area with u and v component of the wind (Table 2).

Interannual variations
Seasonal percentage anomalies of the MODIS chl-a, total river discharge, precipitation, and wind speed were  Fig. 1. Dashed lines and error bars in b-e are a reference of high chl-a of 10 mg m −3 and temporal mean ± standard deviation of chl-a in each month, respectively examined in relation to ENSO and IOD events. Over the entire period of our study, bi-monthly MEI (− 2.3 to 2.8) and DMI (− 1.4 to 1.4) (Fig. 7a) indicated five EN and six LN years of ENSO, while positive and negative IOD often coincided with the ENSO events. Variations in seasonal anomalies of chl-a (− 30 to 41%), total river discharge (− 66 to 115%), precipitation (− 64 to 107%), and wind speed (− 31 to 30%) were with standard deviations (SD) of 17%, 42%, 37%, and 10%, respectively (Fig. 7b-d).
Chl-a anomaly showed no correlation with DMI; whereas it showed significant and non-significant negative correlation with MEI during NEM (r = − 0.54, p < 0.05) and SWM (r = − 0.41, p > 0.05), respectively (Fig. 8a-f). The chl-a anomaly was mostly positive (negative) in LN (EN) during NEM and SWM (Fig. 8d, e). The positive chl-a anomaly exceeded 30% of the annual variation in SWM of 2007 (34%) and 2010 (41%) with LN, and the negative chl-a anomaly below -30% occurred in SWM of 2012. The chla anomaly was abnormally positive in SWM and NEM of 2009 with EN. On the other hand, chl-a anomaly showed no correlation with MEI during NOM, and positive chl-a anomaly occurred only in 2005, 2010, and 2015 with EN (Fig. 8f).
Anomaly of river discharge in all seasons were negatively correlated (r = − 0.49 to − 0.60) with MEI, although it was only significant at p < 0.05 during NOM (Fig. 8g-i). The correlation was significant during SWM (r = − 0.58, p < 0.05) and NEM (r = − 0.67, p < 0.01) when the flood of 2006 under EN was excluded. The positive river discharge anomaly exceeded 45% from SWM and NEM of 2011 to NOM of 2012 under LN (89-115%) just after the precipitation anomaly (107%) in NOM of 2011 (the late super LN2010-11) (Fig. 7). The negative river discharge anomaly below − 45% from NEM of 2014 to NOM of 2016 under EN with the minimum in NEM of 2015 (− 66%) (the middle super EN2015-16). Negative anomaly occurred during the early strong LN in SWM of 2010 just after EN.
Anomaly of precipitation was negatively correlated with MEI in SWM (r = − 0.55, p < 0.05) and NOM (r = − 0.81, p < 0.0005), and not correlated in NEM (Fig. 8j-l). The positive (negative) anomaly peaked in NOM of 2011 (2016) with super LN (EN). The precipitation anomaly was abnormally positive during SWM of EN in 2009 just after LN, which coincided with the oddly positive chl-a anomaly during EN. Anomaly of wind speed and MEI was correlated positively in SWM (r = 0.47) and NOM (r = 0.35), and negatively in NEM (r = − 0.43), but these relationships were not statistically significant at the 0.05 level (Fig. 8m-o). The wind speed anomaly was up to 30% during NOM of 2005 with EN.
The anomalies of chl-a and total river discharge were positively correlated over the entire area, particularly in the west and east of Chao Phraya River toward the offshore during Table 2 Correlation analysis results of monthly data of chl-a, river discharge, precipitation, and wind speeds *Correlation is significant at the 0.05 level (2-tailed)  Fig. 1, regional river discharges (bars), and wind speeds and directions over the uGoT (arrows). Dash lines are chl-a of 10 mg m −3 the SWM and NEM, respectively, and only along the central and west coast during the NOM (Fig. 9a-c). Using the area average of chl-a for the whole Gulf ( Fig. 9d-f), the correlation was distinctly significant during NEM. The correlation was also significant for the SWM and NOM, but only if the years impacted by ENSO events during SWM (i.e., 2010 and 2011 with LN, and 2012 just after LN) and NOM (i.e., 2005, 2010, and 2015 with EN and 2012 with LN) were excluded.
Besides river discharge, chl-a anomaly was also significantly correlated to anomalies of precipitation during the SWM and wind speed during the NEM (Table 3). The positive correlation of chl-a with precipitation anomaly during the SWM and with wind speed anomaly during the NEM was more significantly or equally well correlated as the relationship of chl-a with river discharge anomalies. The low correlations with river discharge could be explained by the variation of chla anomaly in the specific years (Fig. 10). During the SWM, the chl-a anomaly was positive almost the entire area during 2010 and in the northern half of the uGoT during 2011 with LN coincided with positive precipitation during low and high river discharge, respectively. It was negative in almost the whole area during 2012 just after LN when both river discharge and precipitation were low. During the NEM in 2011 with LN, chl-a anomaly was positively correlated with the strong wind and high river discharge for the entire study area, except in the west. During the NOM in 2005 and 2010 with EN, the positive chl-a anomaly was high in the east toward offshore area with positive anomaly of wind speed during low river discharge and precipitation. On the other hand, the positive chl-a showed only nearshore when the wind was weak during the positive anomalies of river discharge in 2012 with LN and precipitation in 2015 with EN.

Improvement of satellite MODIS chl-a
The local switching R rs 488∕R rs 547 algorithms developed by us helped improve the accuracy of chl-a estimation over that possible by the standard MODIS chl-a algorithm. Our switching algorithm especially helped in reducing the overestimates generated by the standard algorithm when chl-a values were lower than 3 mg m −3 (Fig. 3a, b). One of the major reasons why the standard NASA chl-a algorithm does not perform well in coastal waters is because most of in situ data (i.e., the NASA bio-Optical Marine Algorithm Dataset) that has been used in developing it, come from oceanic waters, where the influence of CDOM and suspended sediments on the optical properties of seawater are minimal (Werdell and Bailey 2005). Similar to the studies of Matsumura et al. (2006), our R rs band ratio algorithms developed for the MODIS clearly reduced the overestimation of chl-a for the low chl-a water, which in reality is probably the result of high inorganic suspended solids and CDOM (Kobayashi et al. 2011). Although MODIS R rs 488 and R rs 547 Fig. 9 Spatial correlation of seasonal anomaly of chl-a (a-c) and scatter plots of the chl-a averaged in the whole area (d-f) with total river discharge during SWM (a, d), NEM (b, e), and NOM (c, f). Colors, blackline, and grey strip are the data years, linear fit, and 95% prediction interval of the significant relationship. Hollow dots are years causing the insignificant of the correlation. Stars indicate the significance of the correlation (p < 0.05) (color figure online) were underestimated because of errors in the atmospheric correction, comparison of in situ and MODIS showed no applicable difference for R rs 488∕R rs 547 (Fig. 2c). Because our new local switching algorithm is developed using a database comprised of several previously collected bio-optical datasets, it is more valid spatially covering both clear to turbid waters.

Seasonal chl-a variation
The seasonal chl-a variation that we report here is based on the monthly climatological data which were generated using improved MODIS chl-a data. The spatial and temporal continuity of the climatological data (Fig. 4) improved our understanding of the seasonal cycle of chl-a variation and the controlling mechanisms. Our results demonstrated that the monsoons crucially drive seasonal chl-a variation in the uGoT. First, we described the monthly climatology images of MODIS showing high chl-a near the northern coast, and the westward and eastward shift of the high chl-a during SWM and NEM, respectively. Similar seasonal variation of the distribution of high chl-a water and correspondence to the monsoon was also reported in the previous studies by ship survey (Matsumura et al. 2006) and satellite observations Intacharoen et al. 2018). Time series of the whole area average of monthly climatological chl-a, showed gradual increase during the NOM and SWM, and rapid decrease during NEM (Fig. 5b). This chl-a variation was totally consistent with chl-a observed in the uGoT during October 2003-May 2004 by Matsumura et al. (2006). Meanwhile, the chl-a increase during SWM, particularly in the east (Fig. 5e), was consistent with the increase in cell abundance of green Noctiluca as reported by Sriwoon et al. (2008). Furthermore, we demonstrated that this seasonal variation of chl-a generally corresponded to changes in precipitation and river discharge patterns during the respective seasons. The correspondence was evident from the significant correlation between whole area chl-a, river discharge and precipitation as well as the u-wind (Table 2) indicating monsoonal control of precipitation and 1 3 river discharge. Those time series and the correlations are probably indicating that the story of seasonal increase of precipitation makes increase of nutrient load from river discharge and high average chl-a during NOM to SWM, and decrease of precipitation makes decrease of nutrient load from river discharge and low average chl-a during NEM in the uGoT. Likewise, the chl-a increase owing to increased river discharge induced by high precipitation during SWM is also reported in other areas such as the Bay of Bengal (Amol et al. 2020).
The time series and correlations in whole area average chl-a were similar with the local chl-a average in offshore area (Fig. 5b, d, Table 2). The offshore area of the study was the edge of the high chl-a water and indicating the expansion of the high chl-a water to offshore. The variation of the chla in this offshore area probably corresponds to the offshore extent of the nutrient-rich discharge water (cf. Amol et al. 2020;Sun 2017;Yamaguchi et al. 2012).
On the other hand, regional chl-a along the northern coast showed little correlation with river discharge (Figs. 5, 6, Table 2). Chl-a variation of the three areas (and five areas of Fig. 6) near the coast was strongly influenced by wind as expected from the satellite images (cf. Fig. 4). The high chla during SWM/NEM in the east/west, when river discharge was low, indicated that phytoplankton growth in those areas was supported by the transport of river discharge from the central coast (i.e., Chao Phraya and Tha Chin River). These results are consistent with the results of wind-driven transport of discharged water along the northern coast obtained using a hydro-dynamic model (Buranapratheprat et al. 2002(Buranapratheprat et al. , 2006. Using a coupled hydrodynamic ecosystem with nutrient inputs from the major rivers and monsoon winds, Buranapratheprat et al. (2008b) were able to clearly simulate high chl-a near the central river mouths and the variation of east/west coast during SWM/NEM. Furthermore, we clarified for the first time that the chl-a distribution patterns during the NOM were highest between Tha Chin and Chao Phraya river mouths (Fig. 6), where small canals and coastal aquaculture farms are located [Fisheries Development Policy and Strategy Division (Thailand) 2019]. The results suggested that during the NOM of low river discharge and moderately high precipitation, the major source of nutrients for the large blooms observed in the central region of the uGoT be from small canals and aquaculture farms. It is also conceivable that the strong southerly wind could have pushed the seawater to the north and limited the distribution of nutrients only near the coast.
Based on our studies and the previous studies on uGoT hydrodynamics, seasonal chl-a variation in the uGoT can be summarized as follows. During the SWM, high precipitation enhanced river discharge, mostly in the central uGoT via the Tha Chin and Chao Phraya rivers. Chl-a concentration increased with the increase in river discharge. However, the southwest monsoon wind caused the piling of discharged water towards east resulting in higher chl-a concentrations in the east. During the NEM, the decrease in precipitation and the reverse of wind caused the reduction of chl-a in the east, and transport of the discharged water from central area to the west, resulting in higher chl-a in the western area. During the NOM, increased precipitation and strong southerly wind led to an increase in chl-a in the coastal as well as the offshore areas.
Our analysis showed a higher correlation of precipitation, than river discharge with chl-a in both whole area and most of the coastal areas ( Table 2). The correlation was caused by the increase of chl-a during NOM under increase of precipitation, but a constant river discharge (Fig. 5a, b). Although precipitation has been suggested as a source of nutrients to enhance phytoplankton production in oligotrophic areas (Paerl 1985;Zou et al. 2000;Kim et al. 2014;Sedwick et al. 2018), its influence in a high nutrient coastal area like the uGoT may be minimum. Those higher correlations between chl-a in the coastal area of the uGoT and precipitation can be attributed to drainage of freshwater with nutrients from the surrounding land following precipitation and not due to river discharge from major rivers. This lack of a robust correlation between river discharge and chl-a might the result of a lag time between river discharge and precipitation caused by the long distance between the river-discharge measuring station, and the river mouths, and also the regulated release of water through dams all along the river (Kotsuki and Tanaka 2013).
Our analysis did not show any significant correlation between coastal chl-a and wind speed (Table 2). On the other hand, chl-a showed positive correlation with wind speed in offshore area during strong wind in EN events (Fig. 10). This implied that wind-induced vertical mixing in the water column might also be a factor to enhance chl-a during NOM. These results are in contrast with those obtained from modeling studies which showed chl-a enhancement during weak winds in October and during low periods of weak vertical diffusivity in the water column (Buranapratheprat et al. 2008b) as was also shown in March 2009 using in situ data (Buranapratheprat et al. 2010). More detailed in situ studies will be required to better understand the role of winds in regulating chl-a distribution in the shorter time scale events than seasonal scale in the uGoT.

Interannual chl-a variation
The influences of ENSO and IOD on the interannual chla variation between 2003 and 2017 were studied because of the known impacts of these events on monsoon variability and precipitation. Our results indicate that ENSO events have a far greater influence on uGoT chl-a than IOD (Fig. 8a-f), confirmed from low correlation between IOD and precipitation ( Fig. A3g-i, Supplementary Appendix 3). Kim et al. (2020) illustrated the impact of IOD on rainfall variability over the uGoT is not as intense as in Myanmar (the west of Thailand) and the southern part of Thailand, except when it happens with ENSO events. We expect that the lower IOD impact in the uGoT might be caused by the mountain ranges (200-1000 m) between Thailand and Myanmar separating the Indian and Pacific oceans. Thus, IOD might be modulated by ENSO resulting in the low correlation between IOD and the variabilities of precipitation and chl-a.
Chl-a concentrations increased (decreased) during LN (EN), similar to those reported in the southern region of the strait of Malacca (Siswanto and Tanaka 2014) and the northern, middle, and western South China Sea (Tang et al. 2011;Yuan-Jian et al. 2012;Zhao et al. 2018;Zhao and Tang 2007). Our results; however, differ with those from the Indonesian Seas, such as the Bali Strait, Halmahera Sea, and Gulf of Tomini (Sambah et al. 2021;Sari et al. 2018;Setiawan et al. 2020;Wijaya et al. 2020) when are increase/ decrease in chl-a was observed during EN/LN events. In contrast to our observations in the uGoT, strong influence of IOD on chl-a variation has been reported in the eastern and southeastern tropical Indian Ocean (Iskandar et al. 2009(Iskandar et al. , 2017Mashita and Lumban-Gaol 2019;Sari et al. 2020). It is worthy to note that in most of the other regions, chl-a variation during ENSO or IOD events is associated with strong upwelling and upward Ekman pumping, which is in contrast to the uGoT, where we demonstrated that variability in chl-a is the result of ENSO driven variations in river discharge, precipitation, and wind speed (Fig. 8d-o).
Precipitation, induced by moisture-laden southwest winds during the SWM and by the south wind during NOM, anomalously increased/decreased in association with weak/strong winds during LN/EN (Fig. 8j-o). The anomalously high/ low precipitation patterns resulted in anomalously high/ low river discharge in all seasons during LN/EN events ( Fig. 8g-i). The onset of an EN, induced low river discharge, which at times continued into the following season well past the event, such as SWM of 2010 under LN just after EN (Fig. 7b, c). This underlined the time lag between the two as discussed in Sect. 4.2, and partly might be related to water management to meet agricultural and industrial needs during the EN dry condition (Yamsiri 2014).
The LN-induced high river discharge during 2011-2012 is referred to as the historical great flood of Thailand (Gale and Saunders 2013;Promchote et al. 2016). The flood occurred over a wide area, but intensively in Chao Phraya River basin causing considerable damage and negative impacts on human life and the economy of Thailand (Aon Benfield 2012). Flooding also occurred during EN (e.g., 2006), but it was caused by the crossing of low-pressure troughs and tropical cyclones [National Hydroinformatics and Climate (THAIWATER) 2006]. The intensity of precipitation and the flooding was not as large as the LNinduced flood of 2011-2012.
In general, chl-a variation in the uGoT was significantly correlated to variations in river discharge (Fig. 9), emphasizing the importance of the interannual variation of river discharge on phytoplankton variability in the uGoT. However, this strong correlation was observable when river discharge data for years associated with LN/EN during SWM/NOM and the 2011-12 flood were excluded. These observations suggest that the chl-a variation might also be affected by other factors besides river discharge, such as the other significant environmental factors (i.e., precipitation and wind speed, Table 3). Given that the focus of this study was on the monsoon-influenced environmental factors, hereafter, the impact of ENSO events on the variability of chl-a in the uGoT is discussed seasonally.

Southwest monsoon season (SWM)
As discussed in Sect. 4.2, SWM is the season of enhanced precipitation and river discharge that typically leads to an increase of chl-a. Chl-a anomaly, river discharge and precipitation were negatively correlated with MEI (i.e., increase/ decrease with LN/EN) and positively correlated with wind although only the correlation with only precipitation was significant (Fig. 8d, g, j, m, Table 3). Meanwhile, while the correlation between chl-a anomaly and precipitation was significant, the correlation between chl-a anomaly and river discharge was not significant unless the data for 2010-2012, and for the super and moderate LN were excluded (Fig. 9d).
In 2010, the chl-a anomaly was higher, despite river discharge being lower than expected during a LN event (Fig. 7). Precipitation and wind speed were high and low, respectively, as expected as super LN. The higher chl-a anomaly was clearly caused by a uGoT wide increase in chl-a (Fig. 10). We believe that this increase could be the result of land-based nutrient-rich freshwater drainage into the uGoT, combined with low winds that could have enhanced water column stability resulting in the great chl-a accumulation near the sea surface. The chl-a enhancement during the low wind condition is consistent with simulated high surface chl-a when weakening the winds during May 2004 and July 2005 in Buranapratheprat et al. (2008b).
In 2011 (the flood year), chl-a anomaly was closer to the average anomaly for the SWM and not very high as LN, although the river discharge was extremely high (Fig. 7), as reported in Intacharoen et al. (2018). High river discharge was the result of heavy precipitation, but the wind speed was low. The distribution of chl-a anomaly indicated that chl-a was high in most of the northern area and low in the southern area (Fig. 10). The southwesterly wind, which controls the surface uGoT circulation in this season (Buranapratheprat et al. 2002(Buranapratheprat et al. , 2006, could have reduced the spread of freshwater to the south leading to the accumulation of chl-a in the north. In 2012, chl-a anomaly was very low, and river discharge, precipitation and wind speed were fair as a normal ENSO year (Fig. 7). The distribution of chl-a anomaly indicated that the chl-a was low in the whole uGoT (Fig. 10). Even though the chl-a reduction can be the combined effect of low river discharge and precipitation, it also seems to relate to after the LN. Since limitation of the data, it was not possible to clarify in this study.

Northeast monsoon season (NEM)
As discussed in Sect. 4.2, NEM is the season of decrease of chl-a after the maximum of precipitation and river discharge. Anomalies of chl-a, river discharge, and wind speed were negatively correlated with MEI (i.e., increase/decrease with LN/EN), although only chl-a was significant (Fig. 8e, h, k, n). The chl-a anomaly was significantly correlated with river discharge and wind speed (Fig. 9e, Table 3).
The distribution of chl-a anomaly during 2011 flood with the strong wind as LN was intensively positive in the eastern and offshore area, while negative in the west (Fig. 10). The positive chl-a anomaly, which was more extended to the south than during the SWM, might be due to the northeasterly wind induced southwestward flow (Buranapratheprat et al. 2002) pushing the high freshwater discharge into the south. The large expanse of chl-a was consistent with the spread of freshwater in the southern part of the Gulf due to the stimulated high river discharge of October 2014 (Yu et al. 2018). However, negative chl-a anomaly found in the west and south of the Chao Phraya River is possibly due to short flushing time resulted from strong southward transport during NEM. This was induced by strong northeasterly wind and estuarine circulation during this intense LN year.

Nonmonsoon season (NOM)
As discussed in Sect. 4.2, NOM is the season of low and slight increase of chl-a after the minimum with slight increase of precipitation and constant river discharge. Anomalies of chl-a and wind speed were positively correlated (i.e., decreased/increased with LN/EN), and river discharge and precipitation were negatively correlated with MEI (i.e., increase/decrease with LN/EN), although only river discharge and precipitation were significant (Fig. 8f, i, l, o, Table 3). The correlation between chl-a anomaly and river discharge was only significant when the data of 2005, 2010, 2012, and 2015 was excluded (Fig. 9f).
In 2005 and 2010, chl-a anomaly was high, although river discharge and precipitation were low as EN year. Wind speed of those years was relatively higher as the EN year (Fig. 7). The distribution of chl-a anomaly indicated that the chl-a was high in the eastern and southern areas which corresponded to lower chl-a areas (Figs. 4,10). Low river discharge during intense surface heating in NOM (Yanagi et al. 2001;Buranapratheprat et al. 2008a) of EN can lead to thermal stratification. The strengthening of wind; thus, potentially generated stronger vertical mixing to break down the stratifications and facilitate chl-a enhancement at the sea surface as shown in the chl-a enhancement in the offshore area during the events. Because this area was expected to be not influenced by river discharge, we suspect vertical mixing transported nutrients or high chl-a, which was accumulated in subsurface or near-bottom during the high river discharge seasons, although we cannot find any nutrient or subsurface chl-a data to support this hypothesis at this moment.
In 2015, chl-a anomaly was also high, and river discharge and wind speed were low as super EN (Fig. 7). However, the precipitation was slightly high. Distribution of the high chl-a anomaly indicated that the enhanced chl-a was mainly nearshore area (Fig. 10). This may indicate that the major sources of nutrients for the increase with high precipitation were from the coastal area not through the river discharge stations as discussed in Sect. 4.2.
In 2012, chl-a was low, although the river discharge was very high as the extension of the flood of 2011 (Fig. 7). The precipitation and wind speed were fairly low as LN. The distribution of chl-a anomaly indicated that the low chl-a area was in the middle of uGoT, and chl-a was higher in the shallow coastal area (Fig. 10). Although the percentage anomaly of river discharge was high, the absolute value was not as high as during the SWM and NEM. Since no dominant circulation along the northern coast in this season generated by the southerly wind (Buranapratheprat et al. 2002), we suspect that the northerly wind could have limited the flowing out of freshwater from the major rivers to the uGoT, resulting in the low chl-a.

Conclusions
We improved the accuracy of the MODIS chl-a data for the uGoT using new local empirical chl-a algorithms capable of switching between regions of low and high R rs 547. By applying these algorithms, we investigated the seasonal variations of surface chl-a for the period between 2003 and 2017. We also investigated the interannual variations of chl-a with respect to the Asian monsoon variations due to ENSO-and IOD-induced changes in precipitation, river discharge, and sea surface winds. We demonstrated that chl-a concentrations slightly increased with an increase in precipitation during NOM, increased with increases of precipitation and river discharge during SWM, and decreased with decrease of precipitation and river discharge during NEM. At the same time, monsoonal winds were responsible for transport of surface chl-a and river discharge from the central coast to either east (SWM) or west (NEM), causing a unique pattern of surface chl-a distribution in each season. ENSO rather than IOD had a more significant impact on interannual chl-a variation, by inducing anomalies in river discharge, precipitation, and surface winds. The interannual chl-a variation showed strong correspondence with variations in river discharge and precipitation each season. There are some years with weaker correlation with river discharge and precipitation. In cases like these, it seems that the intensity of winds had a greater role. In near future, it is necessary to study the dynamics of red-tide species in the uGoT with high-resolution satellite data, such as the Second-generation Global Imager (SGLI).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.