The Chl-a data used in this study were obtained from OC-CCI database, one of fourteen CCI projects by the European Space Agency (https://www.oceancolour.org). This dataset provides consistent, stable, long-term error-characterized and merged time series of global ocean color data products from multi-sensor data archives by integrating remote sensing reflectance (Rrs) values using selected algorithms (Brewin et al., 2015; Sathyendranath et al., 2019). The dataset features a fine spatial resolution of 0.04°×0.04°, enabling precise global mapping of Chl-a concentration. For this study, monthly Chl-a composites (Version 6.0, Sathyendranath et al., 2023) specific to the Gulf were extracted, utilizing the extensive temporal range of the dataset, which spans from January 2003 to December 2023. This extended coverage supports the analysis of both historical and recent trends with better validation, offering valuable insights into the long-term dynamics and recent changes in Chl-a concentrations in the region (Jackson et al., 2022). Although the OC-CCI dataset spans from September 1997 onward, we restricted our analysis to the post-2003 period due to known quality issues in the earlier segment of the dataset. During the initial phase (up to April 2002), SeaWiFS was the single contributing sensor, resulting in a relatively sparse sampling frequency and increased noise in the data, including isolated anomalies and spikes (Sathyendranath et al., 2019). From mid-2002 onward, the inclusion of additional sensors such as MODIS-Aqua and MERIS led to a substantial improvement in temporal coverage and data quality, with a manifold increase in the number of observations per day. This multi-sensor merging enhances the robustness of the product by minimizing gaps and random noise while maintaining a stable long-term signal that is essential for climate related studies. Moreover, previous studies on both the Global Ocean and the Gulf have verified that OC-CCI products are more consistent than other available ocean color satellite products with proper validation of in-situ measurements (Belo Couto et al., 2016; Al-Naimi et al., 2017; Moradi, 2021). Moreover, version 6.0 of the OC-CCI uses blended chlorophyll algorithm, which improves performance in Case-2 waters compared to earlier versions that were mostly focused in open-ocean waters (Jackson et al., 2022).

The SST data of Operational SST and Ice Analysis (OSTIA) with 0.05°×0.05° spatial resolution were obtained from the Marine Data Store (MDS) of Copernicus Marine Environment Monitoring Service (CMEMS). OSTIA data provide daily gap-free maps of foundation SST and ice concentration using in-situ and satellite data from both infrared and microwave radiometers (Good et al., 2020). The OSTIA system is run by the Met Office of the United Kingdom and delivered by the French Institute for Ocean Science (IFREMER). It uses satellite data provided by the Group for High Resolution (GHRSST) project together with in-situ observations to determine the SST (https://data.marine.copernicus.eu/product/SST_GLO_SST_L4_REP_OBSERVATIONS_010_011).

Global analyzed sea surface salinity data were accessed from CMEMS’ Multi Observation Global Ocean Sea Surface Salinity and Sea Surface Density (https://data.marine.copernicus.eu/product/MULTIOBS_GLO_PHY_S_SURFACE_MYNRT_015_013). This product consists of daily global gap-free Level-4 (L4) analyses of the SAL and sea surface density at 0.125°×0.125° of resolution, obtained through a multivariate optimal interpolation algorithm that combines sea surface salinity images from multiple satellite sources as NASA’s Soil Moisture Active Passive (SMAP) and ESA’s Soil Moisture Ocean Salinity (SMOS) satellites (Nardelli et al., 2016; Sammartino et al., 2022).

SLA, surface geostrophic northward sea water velocity and surface geostrophic eastward sea water velocity data were obtained from MDS-CMEMS (https://doi.org/10.48670/moi-00148). The data are generated by the processing system combining data from multi-satellite altimetry missions that include all altimeter Copernicus missions (Sentinel-6A, Sentinel-3A/B) and other collaborative or opportunity missions (e.g. Jason-3, Saral/AltiKa, Cryosat-2, OSTM/Jason-2, Jason-1, Topex/Poseidon, Envisat, GFO, ERS-1/2, Haiyang-2A/B). All normal geophysical and environmental corrections have been applied to the data, including sensor errors and tidal impacts, as well as dynamic atmospheric correction. With a horizontal resolution of 0.125°×0.125°, the dataset supports detailed spatial analyses of global sea level changes. Additionally, it offers temporal resolutions ranging from daily to monthly, providing continuous monitoring of sea level anomalies from 1993 to the present. Total CSP was calculated using northward and eastward velocity components.

The wind data used to estimate wind speed and direction were sourced from the ERA5 reanalysis dataset, the fifth-generation global climate and weather reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). ERA5 wind components at a height of 10m above the surface of the Earth with spatial resolution 0.25°×0.25° were obtained from the Climate Data Store of Copernicus Climate Change Service (Hersbach et al., 2023; https://cds.climate.copernicus.eu). ERA5 provides extensive temporal coverage from 1940 to the present, making it an invaluable resource for climate and atmospheric research. As the successor to the ERA-Interim reanalysis, ERA5 offers enhanced accuracy and higher resolution of atmospheric variables. By blending model outputs with a wide range of observational data, ERA5 delivers a globally complete and consistent dataset. ERA5 also supports multiple temporal resolutions, from hourly to monthly, encompassing a broad spectrum of atmospheric, ocean-wave, and land-surface variables, making it a comprehensive dataset for both historical and real-time analysis.

PAR is defined as the quantum energy flux from the Sun from 0.4 µm to 0.7 µm that is absorbed, transferred and stored within ecosystems. PAR data from SeaWiFS and Aqua-MODIS with spatial resolution of 0.04°×0.04° were achieved from NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group (http://oceancolor.gsfc.nasa.gov). For ocean color applications, PAR is a common input used in modeling marine primary productivity. Net primary productivity (NPP) is a measure of the rate at which phytoplankton convert carbon dioxide into organic matter (Westberry et al., 2023). NPP data with spatial resolution 0.083°×0.083° used in this study was computed from satellite-derived Chl-a, SST and PAR data mostly from SeaWiFS, MODIS and VIIRS satellites using the Vertically Generalized Production Model (VGPM, Behrenfeld and Falkowski, 1997), available at ocean productivity website (https://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php).

The standard deviation (SD) of Chl-a in the Gulf provides critical insights into the temporal and spatial fluctuations of phytoplankton biomass. The results were presented across three different temporal scales: annual mean, summer and winter. This SD analysis aided in identifying regions of high Chl-a variability, which could then be used to define study areas (0.5^o×0.5^o boxes) for further investigation to understand the relationship between the distribution of Chl-a and environmental factors. This analysis aligns with previous studies on Chl-a variability in semi-enclosed seas, where nutrient dynamics, physical forcing, and climatic variability strongly regulate phytoplankton distributions (Gregg and Casey, 2004; Behrenfeld et al., 2006). A high Chl-a SD value implies regions where phytoplankton blooms are strongly influenced by any of the seasonal shifts, such as winds, stratification, and nutrient availability mostly from river discharge. In contrast, regions with consistently low SD indicate lower or stable Chl-a throughout the year, often found in deeper or more oligotrophic waters. We have also performed spatial correlation between Chl-a and major environmental variables to confirm the selection of specific regions and thereby to investigate the relationship between Chl-a and other environmental variables. To calculate the spatial correlation, we used data from 2003–2023 and regridded all datasets to a common spatial resolution of 0.04°×0.04°.

Time-series analysis of Chl-a was carried out for the selected boxes on a monthly, seasonal and annual timescales along with a trendline. We applied the modified Mann-Kendall trend test (Hamed and Rao, 1998) to evaluate long-term trends in Chl-a concentration. This method adjusts the variance of the Mann-Kendall statistic to account for autocorrelation in the time series, providing a more reliable assessment of trend significance. The trend calculated for NPP data were normalized to the range of –1 to +1 to ensure comparability with other variables, as this scaling transforms the data to a dimensionless interval where –1 represents the minimum observed value and +1 represents the maximum, preventing its larger magnitude from biasing the analysis.

Empirical Orthogonal Function (EOF) analysis were performed spatially and temporally in the Gulf domain (Longitude: 48–56^oE; Latitude: 24–32^oN) for the period 2003 to 2023. EOF is a robust statistical method used to identify patterns and extract significant features from complex datasets. The dominant spatial patterns of variabilities are presented as different modes of EOF and the temporal variability has shown as time series of principal components. Here, the eigenvalue λk corresponding to k^th EOF gives a measure of the explained variance by ak, where k = 1,2,3,…., p. The following Equations 1-3 are used to calculate the major terms in this analysis. The complete theory of EOF analysis is available at Hannachi (2004).

The explained variance is often presented in percentage as:

The k^th principal component is the projection of the anomaly field X′ on to the k^th EOF ak, i.e.:

where t, s represent the temporal and spatial dimensions. The leading EOFs were analyzed to reveal significant spatial patterns of Chl-a distribution, while the corresponding principal components provided insights into temporal changes. We removed the season and trend from the data for this analysis.

Mutivariate analyses such as correlogram and Principal Component Analysis (PCA) were also done to quantify the relationship between Chl-a and other environmental variables by identifying the key Chl-a patterns of spatial variability (Jackson, 1991; Jolliffe, 2002). These analyses were also aimed to reveal the dominant factor for Chl-a variability at selected locations. Generally, these are useful for reducing the dimensionality of large number of data into a smaller number of orthogonal factors called Principal Components (PCs). PCs represent the major part of the overall data variance along with the correlation amongst the variables without losing much information (Vega et al., 1998; Helena et al., 2000; Lins et al., 2018). In simple terms, these methods break the complex data into a few main patterns, showing where the highest changes occur in the Gulf and how they evolve over months and years. A two-tailed student’s t-test was applied to assess whether the correlation coefficients (r) between Chl-a and other environmental parameters were statistically significant (Kreyszig, 1970). Interannual variability of Chl-a and other associated variables was also analyzed by removing the seasonal cycle from each variable.

The monthly climatology spatial maps of Chl-a concentration provide an overview of the spatiotemporal distribution of phytoplankton biomass in the Gulf in different months (Supplementary Figure S1). Strong seasonal and spatial discrepancy in the distribution of Chl-a was clearly observed in the climatology map, with clear distinctions between coastal and offshore waters. Elevated Chl-a values (>2 mg/m³) were predominantly observed in coastal regions, particularly near the north and northwestern coasts throughout the year. In contrast, the eastern coasts and offshore regions of the Gulf exhibit lower Chl-a. The map also revealed marked seasonal variability in Chl-a concentration all over the Gulf. During winter (December–March), higher Chl-a exhibited (>2 mg/m³) along the central (Iran coast) and southeastern coast adjacent to the Strait of Hormuz and moderate Chl-a values along north and western coasts. During summer (June–September), nutrient availability in the surface layer is limited due to the high SST (>35°C) and strong stratification (Polikarpov et al., 2016). As a result, Chl-a concentrations were generally lower (<2 mg/m³) during this period. In contrast with the winter, higher Chl-a (~3 mg/m³) was exhibited along the central (Saudi and Qatar coast) and southwestern coasts (U.A.E.) of the Gulf during summer (Supplementary Figure S1). The transitional periods of Spring (April–May) and Autumn (October–November) months exhibit moderate to lower Chl-a levels (~1 mg/m³), reflecting the gradual shifts in environmental conditions. The observed patterns of monthly climatology provide insights into the modulation of Chl-a in the Gulf and its response to seasonal changes.

The monthly climatology maps of SST overlaid with winds (Supplementary Figure S2) and SLA overlaid with surface currents (Supplementary Figure S3) provide valuable insights into major physical processes driving Chl-a variability in the Gulf. The SST climatology maps reveal pronounced seasonal variability, with temperature ranging from approximately 17°C in winter to 35°C in summer. The superimposed wind vectors highlighted the role of atmospheric forcing in modulating SST and, consequently, Chl-a concentration. During winter, cooler SST (17–25°C) dominated the Gulf waters, with the lowest SST observed in the northwestern region maybe due to the enhanced vertical mixing. In contrast, summer was characterized by extreme SST (29–35°C), with the highest values observed in the southern and central Gulf (Supplementary Figure S2). During this period, strong thermal stratification limits vertical mixing, reducing the upward transport of nutrients to the surface layer. The transitional periods (spring and autumn) exhibited moderate SST (21–29°C), reflecting the gradual warming in spring and cooling in autumn. Wind patterns during these periods are less intense and more variable and hence influence the mixing and nutrient dynamics. These transitional seasons were marked by moderate Chl-a levels, since these seasons are in between the winter bloom and summer stratification phases.

The SLA climatology maps, overlaid with surface currents, provided further insights into the circulation patterns and their influence on nutrient transport and Chl-a distribution in the Gulf (Supplementary Figure S3). SLA reflects changes in sea surface height due to variations in SST, salinity, and wind forcing, while surface currents play a critical role in redistributing nutrients and phytoplankton blooms. During winter, positive SLA (up to 10 cm) were observed in the southern Gulf, indicating higher sea levels due to cooler, denser water and wind-driven convergence (Al-Subhi and Abdulla, 2021). In summer, SLA varies from 3–8 cm with the presence of strong cyclonic eddies along central and southeastern coasts. During the transitional periods (spring and autumn), SLA ranges varied from –1 to 10 cm, reflecting the gradual transition between winter and summer conditions. Surface currents during these periods were weaker, with moderate flows regulating nutrient distribution. The interaction between SLA and circulation patterns during these seasons creates a dynamic environment for phytoplankton, with moderate nutrient availability and transitional Chl-a levels. The overall counterclockwise current pattern plays a critical role in explaining the observed distribution and transport of the algal blooms in the Gulf. The current enters the Gulf through the Strait of Hormuz, flows northwestward along the eastern Gulf, and then loops back southward along the western Gulf before exiting through the Strait (Mussa et al., 2024). The climatology analyses (Supplementary Figures S1–S3) highlight the overall distribution of Chl-a with unique oceanographic and environmental characteristics of the Gulf, which strongly influence Chl-a dynamics. The observed patterns in SST, winds, SLA, and surface currents provide insights into the complex relationship between physical processes and phytoplankton dynamics. In order to select specific locations for further analyze the spatiotemporal variability of Chl-a and its major environmental drivers, we performed SD of Chl-a and spatial correlation between Chl-a and major environmental variables.

We analyzed the SD of monthly Chl-a concentrations to identify regions with high seasonal or interannual fluctuations. The annual mean SD Chl-a map represents the mean spatial variability, providing the total variability of annual Chl-a concentration during the study period (Figure 1b). Higher SD was exhibited along the coastal areas, particularly in the northern Gulf near the Shatt al-Arab is mainly attributed to the nutrient-rich river discharge (Nezlin et al., 2010; Al-Naimi et al., 2017; Al-Thani et al., 2023). On the other hand, offshore regions experience relatively lower SD, reflecting a more stable environment with less pronounced productivity and seasonal changes. The summer SD highlights spatial variations in Chl-a during warm months when the Gulf experiences high SST and intensified stratification (Figure 1c). Phytoplankton productivity during this period is often limited by nutrient availability, as thermal stratification restricts the upward transport of nutrients from deeper waters (Doney, 2006; Al-Thani et al., 2023). However, localized high SD values in certain coastal and shallow regions, particularly along northern Gulf indicate wind-induced mixing, sediment resuspension, or nutrient input from riverine sources. Mesoscale features such as eddies and upwelling events may also contribute to localized Chl-a variability during this period (Thoppil and Hogan, 2010a, Thoppil and Hogan, 2010b). In contrast, the winter SD reveals Chl-a variability during the cooler months (Figure 1d), a period typically associated with enhanced vertical mixing and nutrient replenishment in the Gulf (Al-Thani et al., 2023). The presence of high SD values along the southeastern coast towards the Strait of Hormuz in winter suggested that phytoplankton biomass is highly responsive to seasonal environmental changes, particularly nutrient availability driven by convective mixing and wind-driven turbulence. Northern coastal regions exhibited greater variability due to fluctuations in salinity and nutrient influx influenced by freshwater discharge. These high-variability zones provided insights for selecting specific regions to further explore the relationships between Chl-a and environmental drivers in the Gulf.

The spatial correlation between Chl-a and SST revealed a clear regional gradient. There were moderate negative correlations along the entire eastern and northern coasts, as well as offshore regions (Figure 2). In contrast, moderate to strong positive correlations were exhibited along western and southern coastal areas. This spatial variation highlights how SST control on Chl-a is regionally dependent within the Gulf. SAL demonstrated moderate to strong positive correlations throughout the basin, particularly in the northern half of the Gulf. However, some weak negative correlations were also found near the southern and western coasts. WSP displayed positive correlations in the eastern and central Gulf, especially along the Iranian coasts. Conversely, weak to moderate negative correlations were evident along the western and southern coast. The correlation with SLA was generally weak or no correlation across the Gulf, though some negative patches appear in the eastern and northern coastal regions; whereas positive patches along the western and southern coasts. Similar to SLA, CSP also showed weak scattered positive and negative correlations. Notably, the correlation between Chl-a and PAR was strongly negative throughout most of the Gulf, while narrow patches of moderate negative correlation was found along the westren and southern coasts. Based on the SD map and spatial correlation map, four representative zones (0.5^o×0.5^o boxes), with consistently high SD and correlation values were selected for further analyses.

The detailed analysis of long-term Chl-a over the four selected regions in the Gulf provides more insights into spatial and temporal variability, seasonal dynamics, and long-term trends. Even though all the boxes displayed decreasing trend, each region exhibited distinct characteristics both in spatial and temporal scales (Figure 3). Box 1 (central-eastern Gulf, along Iran coast) exhibited significant declining trend, with a total decrease of –0.31 mg/m³ and values ranging from 3.3 to 0.7 mg/m³ throughout the study period. Seasonal patterns were evident, with winter peaks of 2.3 mg/m³ and summer lows below <1 mg/m³. Interestingly two peaks during 2005–2006 and 2008–2009 were identified at this region with the values of Chl-a >3 mg/m³. Box 2 (northwestern Gulf, Kuwait coast) was also experienced a significant decline with a total reduction of –0.66 mg/m³. However, highest Chl-a variability (from 8 to 1.5 mg/m³) was displayed at this region during the study period. Winter peaks of 4.3 mg/m³ reflect strong productivity, while summer lows below 1.8 mg/m³ pointed to seasonal nutrient limitations.

Interestingly, unlike from boxes 1 and 2, the declining trend in both Chl-a and seasonal peaks were weak for the remaining boxes (3 and 4); in which box 4 showed non-significant correlation. In contrast, winter dominance of higher Chl-a observed along the eastern and northern Gulf (boxes 1and 2), seasonal shift in Chl-a maximum from winter to summer were observed for boxes 3 and 4 (western and southern Gulf). Box 3 (central western coast of the Gulf, along the Saudi Arabia) illustrated a decreasing trend, with Chl-a values vary from 3.8 mg/m³ to 1.3 mg/m³ and net variability of –0.18 mg/m³. Summer peaks reach up to 3 mg/m³, while winter lows remain above 1.6 mg/m³. Finally, Box 4 (Southern Gulf, UAE coast) showed a decrease in Chl-a concentration, ranging from 4.5 to 1.4 mg/m³ and a total decrease of –0.1 mg/m³. Summer peaks were reached up to 4 mg/m³, while average winter Chl-a showed <2 mg/m³. In summary, time series analysis reveals an overall decreasing trend with clear regional contrasts in both seasonal and long-term Chl-a. The eastern and northern coasts (boxes 1 and 2) show severe decline and winter dominants, while western and southern coasts (boxes 3 and 4) exhibited comparatively moderate/ weak waning in Chl-a with summer peaks. The overall decreasing Chl-a trend in the time series analysis is also consistent with recent studies (Moradi, 2020; Yang et al., 2024). The seasonal dominance patterns along the eastern and western regions of the Gulf are shaped by the interaction of physical and regional factors. Winter dominance is primarily driven by nutrient mixing and favorable temperature, while summer dominance is linked to sustained nutrient inputs and higher light availability in coastal regions (Moradi and Kabiri, 2015). These findings underscore the importance of localized environmental controls, which are further examined in the following sections.

The time series analysis of Chl-a reveals discrete peaks and declines across various boxes in the Gulf, corresponding to different periods. Previous studies confirmed the indirect influence of large-scale climatic oscillations such as El Niño-Southern Oscillation (ENSO), Indian Ocean Dipole (IOD), Atlantic Multi-decadal Oscillation (AMO) and North Atlantic Oscillation (NAO) on phytoplankton dynamics in the Gulf through changes in SST, winds, mixing, nutrient supply, and precipitation patterns (Nezlin et al., 2007; Nezlin et al., 2010, Al Senafi, 2022; Hamdeno et al., 2022). Notably, box 1 exhibits significant peaks during 2008 to 2009, highlight the influence of episodic algal blooms or extreme weather events, which suggest further investigation to understand their ecological impacts (Richlen et al., 2010). Nezlin et al. (2010) observed that the bloom occurrence in 2008 was associated with variability in NAO phases. Al Senafi (2022) showed that ENSO, IOD and AMO significantly influence the SST and circulation, thereby affecting marine productivity in the Gulf. Moreover, recent work has also highlighted the impact of marine heatwaves on Chl-a, and ENSO/ IOD variability in the region (Hamdeno et al., 2022; Shaltout and Eladawy, 2024). Long-term warming trends further amplify these dynamics, contributing to stronger stratification and reduced nutrient replenishment (Roxy et al., 2016; Al Senafi, 2022). Overall, while this study does not directly focus on extreme events, past literatures indicate that climate oscillations can modulate Chl-a variability in the Gulf through indirect pathways, emphasizing the importance of future targeted analyses of ecosystem responses to these events.

We applied Correlogram, PCA and EOF analyses to explore how Chl-a varies over space and time and to identify dominant patterns deriving this variability in the Gulf. Correlogram and PCA helps to quantify the relationship between Chl-a and environmental factors and identify the most influencing environmental drivers. EOF, in turn, decomposes the Chl-a dataset into spatial modes that capture the most significant variance, while the corresponding time series (principal components) extracts temporal evolution associated with these modes. Together, these approaches allow us to distinguish large-scale, recurring patterns such as seasonal cycles, localized or interannual variations, making the complex dynamics of Chl-a easier to interpret. The insights from these analyses are detailed in the following two subsections.

Correlogram and PCA were performed by including PAR as another forcing factor along with the variables used in the previous analysis as causative factors for Chl-a modulation. The correlation analysis of Chl-a with SLA, SST, SAL, WSP, CSP, and PAR (Figure 4) across four spatial boxes (Figure 1) in the Gulf revealed complex, region-specific relationships, highlighting strong spatial differences in the environmental drivers of phytoplankton dynamics. All variables showed significant relationship with Chl-a except CSP in box 1. Here, a negative correlation between SST and Chl-a was evident, with a strong correlation (r= –0.50), which may suppress phytoplankton growth due to high SST. A notable positive correlation between Chl-a and WSP (0.46) emphasized the role of wind-driven mixing in enhancing nutrient availability, which supports phytoplankton growth. A moderate positive correlation with SAL (0.42) suggested that higher salinity waters in this box may coincide with nutrient-rich conditions, potentially influenced by regional circulation patterns or saline water influxes. Further, the strong negative correlation between Chl-a and PAR (–0.60) underscored the influence of photoinhibition in this region, where excessive light exposure inhibits phytoplankton productivity (Behrenfeld et al., 2001).

In box 2, Chl-a exhibited a very weak non-significant negative and positive correlation with SLA (–0.13), WSP (–0.03), and CSP (0.02) respectively (Figure 4). So, these variables have a limited role in nutrient dynamics and phytoplankton growth in this region. A weak negative correlation with SST (–0.26) highlighted the minimal impact of thermal stratification in reducing phytoplankton productivity. Notably, the negative correlation with PAR (–0.41) indicated severe photoinhibition similar to previous location (box 1), where high light levels hinder photosynthesis, particularly in nutrient-limited conditions (Behrenfeld et al., 2001; Falkowski and Raven, 2013). In box 3, all variables had a significant correlation with Chl-a. SLA, SST and PAR exhibited weak to strong positive correlation with Chl-a (r = 0.13, 0.79 and 0.60) respectively, which indicated their roles in promoting phytoplankton growth. However, the strong negative correlation with SAL (–0.50) suggested that lower salinity conditions, possibly caused by freshwater inputs, may limit phytoplankton productivity by increasing stratification. A strong negative correlation with WSP (–0.58) may lead to the disruption of phytoplankton stability in the water column due to excessive wind-driven mixing, thence diluting localized Chl-a concentrations. Whereas, a strong positive correlation with PAR supports photosynthesis in this region with high light availability. Finally, box 4 revealed a very strong positive correlation between Chl-a and SST (0.85), indicating that seasonal warming supports phytoplankton growth in this area, likely by enhancing nutrient cycling and biological processes. A moderate negative correlation with WSP (–0.52) implied that wind-induced turbulence may dilute phytoplankton concentrations. Interestingly, PAR showed strongest positive correlation (0.72), suggesting that in this region, light availability supports photosynthesis without significant photoinhibition.

While correlation analysis identifies pairwise relationships, PCA provides a broader multivariate perspective by revealing dominant environmental regimes that drive Chl-a variability. We considered the first two dominant principal components (PC1 and PC2), which explain most of the variance (>60%) in Chl-a and highlight the strongest environmental controls in the Gulf. The key insights obtained from the PCA (Figure 5) are as follows. The analysis showed that high SST and PAR were associated with opposite phases and lower Chl-a for boxes 1 and 2. This indicated that enhanced surface warming and photoinhibition reduce phytoplankton growth in these regions. For instance, the dominant drivers in box 1 include a negative relationship of Chl-a with SST and PAR and a positive influence of wind (WSP). Here, enhanced Chl-a was promoted by strong wind-induced upwelling with relatively cooler surface conditions (Al-Yamani et al., 2020; Ismail and Al Shehhi, 2022). Similarly, in box 2, the negative correlations with SST and PAR, combined with the positive contributions from salinity (SAL) and WSP, underscore a regime where upwelling and nutrient enrichment through river discharge play critical roles in driving Chl-a variability (Ershadifar et al., 2020; Ismail and Al Shehhi, 2022). Similar to the correlogram analysis, the PCA also indicated a shift in environmental drivers in boxes 3 and 4 (Figure 5). In these regions, SST, SAL and WSP emerge as major drivers, pointing to a stratification‐dominated regime where thermal conditions and elevated water levels enhance phytoplankton productivity. Here, PAR also exhibits a positive influence, implying that ample light in conjunction with warmer, more stratified waters supports augmented phytoplankton biomass. Furthermore, the negative associations with SAL and WSP in these boxes indicate that lower salinity and reduced wind mixing favor phytoplankton bloom formation, likely by minimizing dilution and promoting water column stability.

The correlogram and PCA analyses further highlights the complex and regionally variable relationship between Chl-a and environmental parameters in the Gulf. The results further confirm the crucial interplay of physical and biological factors governing phytoplankton productivity in this ecologically and economically vital region. It underlines the importance of analyzing the spatial heterogeneity of major environmental drivers when studying marine ecosystems. The distinct correlation patterns across the four selected boxes indicated that local environmental conditions such as hydrodynamics, light and nutrient availability, or temperature gradients significantly influence the relationships between Chl-a and the selected variables. SST, SAL, WSP and PAR emerge as key factors influencing nutrient dynamics and phytoplankton productivity, with varying impacts across different regions. SST often showed a significant correlation with Chl-a, though the direction (positive or negative) varies across boxes, reflecting potential regional influences of temperature on phytoplankton dynamics. While SST is a strong positive driver for the regions along the southern and western Gulf (boxes 3 and 4), it negatively impacts Chl-a in the eastern and northern regions (boxes 1 and 2). The weak or moderate positive correlation of SLA and SAL with Chl-a, indicating that changes in sea level anomaly and salinity slightly influence nutrient availability or vertical mixing. PAR generally exhibits moderate or strong positive (negative) correlations with Chl-a in western and southern (eastern and northern) boxes, suggesting light availability may (not) be the primary limiting factor in these regions. The consistent negative correlations with PAR, particularly in boxes 1 and 2, highlight the persistent role of photoinhibition in limiting phytoplankton growth under high light conditions, a well-documented phenomenon common in shallow or clear waters with high solar insolation (Behrenfeld et al., 2001). Additionally, correlation between current speed and Chl-a were weaker in most boxes, indicating less direct influence. Beyond the pairwise relationship from correlogram analysis, PCA further confirmed and highlights how multiple factors interact to form distinct environmental regimes and complements the findings by decomposing the variance in Chl-a into key multivariate components. It also demonstrates that while individual parameters have distinct pairwise correlations with Chl-a, their combined effects vary regionally. These findings are essential for refining ecological models and developing effective management strategies in the Gulf, particularly in the context of ongoing environmental changes (Al-Yamani et al., 2020). The overall results of the correlogram and PCA are summarized in Table 1.

EOF decomposition was applied to Chl-a in the Gulf to identify the primary modes of variability and their potential driving mechanisms. The results exhibited four dominant modes: each captures distinct patterns of variability and contributes to understand the Chl-a dynamics in the region (Figure 6). The first mode of EOF accounted for 47% of the total variance, making it the dominant pattern of Chl-a variability. The spatial plot indicated a homogeneous distribution across the Gulf, with stronger variability observed along the northern and eastern coastal regions, particularly near the strait of Hormuz and along the Iran coast (Figure 6a). The corresponding temporal series showed sharp peaks during 2008–2009 (Figure 6b), indicating the influence of extreme climate anomalies apart from episodic algal bloom events (Chinta et al., 2024). These climatic events impact the Indian Ocean monsoon system and oceanic productivity in adjacent regions, including the Gulf (Al-Yamani et al., 2004; Shafeeque et al., 2021a; Seelanki et al., 2022). Other contributing processes may include seasonal upwelling, mesoscale eddies and anthropogenic nutrient inputs, which enhance nutrient availability and drive Chl-a fluctuations (Kämpf and Sadrinasab, 2006; Thoppil and Hogan, 2010a; Ismail and Al Shehhi, 2022). The second EOF mode explained 35% of the variance and is characterized by a localized pattern, with strong positive loadings near the northern Gulf, particularly along the Kuwait and Iraq coasts (Figure 6a). Temporal evolution suggested a distinct seasonal cycle, indicating that this mode is strongly influenced by regional hydrodynamic processes. A key driver of this pattern is likely the river discharges, especially from the Shatt al-Arab river, which injects freshwater and nutrients into the Gulf, stimulating phytoplankton blooms (Abaychi et al., 1988; Sheppard et al., 2010). Additionally, wind-driven circulation and coastal currents play a role in shaping this variability by modulating nutrient transport and stratification dynamics (Mussa et al., 2024).

The third EOF mode captured 10% of the variance, which highlights Chl-a variability mainly along the central and southeastern Gulf. The time series indicates interannual variability, possibly linked to fluctuations in wind-driven mixing and atmospheric dust deposition (Nezlin et al., 2010; Asgari and Soleimany, 2023). The interplay between wind mixing and nutrient availability appears to be a key factor in driving this mode. The fourth EOF mode accounts for 8% of the variance and exhibits a more localized pattern along the northwest (Kuwait coast) and near to the Strait of Hormuz (north U.A.E. coast) of the Gulf. The associated time series suggests short-term fluctuations (Figure 6b), indicating that this mode is influenced by regional oceanographic processes rather than long-term climatic patterns. Additionally, extreme weather events, such as dust storms and cyclones, may contribute to the observed variability by altering light penetration and nutrient fluxes (Nezlin et al., 2010; Jish Prakash et al., 2015; Madhusoodhanan et al., 2024). Overall, the EOF results reveal that Chl-a dynamics in the Gulf are not uniformly driven by seasonal processes, but controlled by a combination of large-scale climatic influences, regional hydrodynamic processes, and localized oceanographic phenomena. These modes offer valuable baselines for long-term monitoring and ecosystem forecasting. The dominance of the first two Modes (82%) highlights the predominant role of climatic and seasonal forcing, whereas the other modes capture more localized and transient processes. This analysis strengthens the importance of understanding multiple scales of variability to improve predictions of phytoplankton dynamics and ecosystem responses in this ecologically and economically significant region.

The interannual variability of Chl-a in the Gulf was briefly revealed in the time series (section 3.3) and EOF analyses (section 3.5), which indicated that Chl-a concentrations have undergone significant fluctuations during the study period. Based on the correlogram and PCA results presented in Table 1, we selected two regions (boxes 1 and 3) representing the eastern and western coasts of the Gulf to further explore interannual changes in Chl-a along with dominant environmental factors (SST, SAL, SLA, WSP and PAR). The time series of these variables are plotted for the period 2003 to 2023 (Figure 7). Moreover, a lag correlation analysis for these boxes was also done and presented as supplementary figures (Supplementary Figures S4, S5). Prolonged high values of Chl-a concentration were recorded in box 1 (along the Iran Coast) from 2006 to 2009 (Figure 7a). The years 2008 and 2009 showed the highest Chl-a anomaly (1.5 mg/m³), indicating an intense phytoplankton bloom; conversely, the lowest Chl-a occurred in 2010, 2013 and 2020, indicating a period of markedly reduced biological productivity. The central-eastern coast of the Gulf, particularly along the Iran coast is characterized by frequent upwelling events that significantly influence Chl-a variability (Ismail and Al Shehhi, 2022). Upwelling brings nutrient-rich deep waters to the surface, supporting phytoplankton blooms and leading to higher Chl-a concentration (Moradi and Kabiri, 2015; Shafeeque et al., 2021b). The period of lower Chl-a could be attributed to weaker winds and higher SST, limiting vertical mixing and reducing nutrient availability. Recent studies have highlighted the increasing frequency of marine heatwaves over this region, which can disrupt upwelling and reduce phytoplankton productivity (Chaidez et al., 2017; Hamdeno et al., 2022; Kashkooli et al., 2022; Shaltout and Eladawy, 2024). Salinity in this region is influenced by the influx of nutrient-rich waters from upwelling and the exchange with the Indian Ocean through the Strait of Hormuz (Lorenz et al., 2020; Ghaemi et al., 2021). Earlier studies have emphasized the role of salinity gradients in shaping phytoplankton communities in upwelling zones (Wabnitz et al., 2018; Ben-Hasan and Christensen, 2019; Alnafissa et al., 2021). The lag correlation analysis for box 1 highlights the temporal shifts in the relationship between Chl-a and environmental parameters, reflecting the influence of seasonal progression, mixing processes, and hydrodynamic conditions. SST exhibited strong negative correlation (–0.62) at lag +1 month, consistent with nutrient depletion during periods of sustained warming (Supplementary Figure S4). SAL exhibited positive correlations (+0.45) peaking at lag +1 month. Wind speed showed maximum correlation (+0.46) at a lag 0, highlighting the direct impact of wind-driven nutrient entrainment on Chl-a production. CSP had moderate positive correlations at long negative lags (+0.30 at –6 months). PAR displayed strong negative correlations near lag 0 (–0.60), likely reflecting that the immediate impact of irradiance on Chl-a. Overall, these results underscore that Chl-a variability in box 1 is shaped by cumulative warming and light availability months in advance, coupled with immediate influences of wind-driven mixing and dynamic changes in salinity and current regimes.

Unlike box 1, peak Chl-a values in box 3 (Saudi Coast) mostly occurred during summer months, with highest Chl-a (3.9 mg/m³) was recorded in August 2017 (Figure 7b). The seasonal peak during summer is associated with intensified thermal stratification, coupled with wind-driven mixing that redistributes nutrients in the upper water column (Reynolds, 1993; Doney, 2006). The year 2017 marked the highest Chl-a concentration, potentially influenced by anomalous weather conditions, including stronger than usual Shamal winds and regional upwelling events (Al Senafi and Anis, 2020). The dominance of SST as the key driver was evident, as this region follows a summer Chl-a peak rather than a winter peak (Figure 7b). Salinity in this region is normally high due to limited freshwater input and high evaporation rates. However, localized reductions in salinity can enhance nutrient availability and support phytoplankton growth (Lattemann and Höpner, 2008). Recent research has highlighted the role of salinity gradients in shaping phytoplankton communities (Al-Thani et al., 2023). In box 3, SST showed its highest correlation (+0.79) with Chl-a without any lag (Supplementary Figure S5). However, at a lag of –6 months, the same correlation also existed with opposite direction (–0.79). Salinity peaked at a lag of –1 to –2 months, exhibiting a strong negative correlation of –0.72. Wind speed and CSP presented their maximum correlations (–0.58 and –0.19, respectively) with 0 lags. PAR displayed a high positive correlation (+0.76) at a lag of –2 months. Despite a lag or lead of a few months, box 3 also revealed more or less similar correlation pattern to the results of correlogram analysis.

The interannual patterns across the two selected regions highlight the disparity in the role of environmental drivers in Chl-a variability. A key distinction between the regions was the seasonality of peak Chl-a concentration. Box 1 exhibited winter-dominated Chl-a peaks, reflecting SST-driven mixing and enhanced nutrient supply during colder months. Box 3 followed a summer peak, suggesting a different nutrient cycling mechanism, possibly linked to temperature-induced stratification dynamics. These findings align with prior sections and studies also indicated that winter blooms in the Gulf are primarily regulated by nutrient enrichment via convective mixing and riverine input, while summer productivity in the western Gulf is influenced by local wind-driven upwelling and water column stability (Wiggert et al., 2005; Thoppil and Hogan, 2010a, 2010b). While interannual variability is controlled by oceanographic and atmospheric forcing, potential impacts of climate change and human activities cannot be ignored. Warming trends in the Gulf intensify stratification, reduce nutrient replenishment, and shifting phytoplankton community structures (Al-Yamani et al., 2004; Wabnitz et al., 2018). Major climatic events have influenced productivity, reinforcing the role of global climate variability in regional ocean dynamics (Doney, 2006; Park et al., 2018; Shafeeque et al., 2021b; Hamdeno et al., 2022; Kashani et al., 2025). The lag correlation analysis of Chl-a with environmental variables further revealed complex distinct temporal patterns shaped by seasonal dynamics, hydrodynamic conditions, and mixing processes. Overall, these patterns demonstrate that Chl-a dynamics in the Gulf are governed by a combination of cumulative seasonal warming and light availability, modulated by lagged effects of mixing, circulation, and salinity changes which need to be further explored in the future work. The varying timing and magnitude of these relationships across boxes further confirmed the spatial heterogeneity of phytoplankton responses to environmental forcing in this highly dynamic marine ecosystem. In addition to the seasonal and interannual variability of Chl-a with contributions from environmental factors, we analyzed the spatial trends of NPP along with Chl-a and SST to enhance our findings in the following section.

A two-decade (2003– 2023) cumulative spatial trend for Chl-a, NPP and SST in the Gulf was calculated and presented in Figure 8. The result showed an overall decrease in Chl-a and NPP trend (Figure 8a) in the Gulf, with small patches of increasing trend along the coasts of western and southern Gulf. Whereas SST showed a warming trend in the northern and central Gulf. The declining trend in Chl-a and NPP is due to the high SST and stratification, which limit nutrient availability in the surface layer (Polikarpov et al., 2016; Al-Shehhi et al., 2021). The oligotrophic nature of these regions makes them particularly sensitive to changes in environmental conditions, such as reduced vertical mixing and nutrient upwelling (Reynolds, 1993; Thoppil and Hogan, 2010a). The declining Chl-a trend was also aligned with the observed decreasing trend in NPP (Figure 8b). This decline is likely driven by the combined effects of warming SST, reduced nutrient availability, and increased stratification (Reynolds, 1993; Behrenfeld et al., 2006; Polikarpov et al., 2016). Recent studies have also linked these trends to the intensification of climate change impacts, such as more frequent marine heatwaves and reduced oxygen levels (Chaidez et al., 2017; Wabnitz et al., 2018; Shaltout and Eladawy, 2024). Continued warming and increasing frequency of marine heatwaves may further weaken offshore productivity, with potential implications for carbon cycling and ecosystem resilience (Wabnitz et al., 2018; Kashkooli et al., 2022).

The observed long-term trends in Chl-a, NPP and SST have significant implications for the marine ecosystem of the Gulf. Coastal regions with decreasing Chl-a and NPP trend could negatively impact pelagic species and commercial fisheries that rely on open-water resources (Mieszkowska et al., 2009; Gamito et al., 2015; Wabnitz et al., 2018). High SST and stratification reduce vertical mixing and nutrient supply to the surface layer (Chaidez et al., 2017). The Gulf experiences a more pronounced SST shift due to its shallow depths and limited water circulation, which result in faster and more intense heating. The SST trend across the Gulf exhibited significant escalation over the past two decades which leads to extensive implications. The escalating SST levels, ranging from 0.5°C in the Strait of Hormuz to 1.5°C in the Northern Gulf, underscore the severity of this thermal shift (Figure 8c). The central Gulf also experienced a ~0.5^oC SST increase during the past two decades. These changes are pronounced in the Gulf as observed in the previous sections, where dominant environmental factors vary regionally. Human activities, such as coastal development, desalination and oil exploration also impact the marine ecosystem. Desalination is the primary source of freshwater in the region, discharges highly saline and chemically altered brine back into the coastal waters. This can alter the salinity and nutrient dynamics and eventually enhance Chl-a and NPP in the coastal areas (Peterson et al., 2003; Purnama et al., 2003; Lattemann and Höpner, 2008). In parallel, large-scale coastal development, land reclamation, industrial expansion, dredging activities, overfishing, and other population-driven stressors have reshaped nearshore habitats of the Gulf and modified circulation patterns, further impacting water quality and primary productivity (Burt, 2014; Naser, 2014; Al-Thani et al., 2023). Further, these activities can also lead to eutrophication and harmful algal blooms, which may disrupt ecosystem balance (Glibert et al., 2014; Madhusoodhanan et al., 2025). Understanding such dynamics is crucial for sustainable marine management, particularly in light of rapid coastal development and climate-driven changes. The ecological impacts of rising SST in the Gulf are manifold. Elevated SST can disrupt marine ecosystems, leading to coral bleaching, changes in species distribution, and alterations in food webs (Riegl et al., 2011; Cavalcante et al., 2016; Harrison et al., 2017). This, in turn, can have cascading effects on biodiversity and fisheries productivity, threatening the livelihoods of communities reliant on marine resources. Moreover, the socioeconomic consequences of escalating SST in the Gulf are significant. Fisheries, tourism, and other industries dependent on the marine environment may suffer due to declining fish stocks, damaged coral reefs, and reduced attractiveness of coastal destinations. The cumulative impacts of the anthropogenic pressures add to the natural variability driven by climate and oceanographic factors, making it difficult to separate human-induced signals from natural ones. Nonetheless, integrating such drivers into the interpretation of Chl-a trends is essential for a more holistic understanding of ecosystem changes in this semi-enclosed and heavily developed marine basin.

Global satellite-derived surface Chl-a data have become an essential tool for studying ecosystem and biogeochemical processes. However, the potential offered by these data sets can be fulfilled only if their uncertainties are quantified, and indeed, this assessment is considered an integral part of any ocean color mission (McClain et al., 2004; Mélin, 2009). Even though the ocean color datasets have been extensively validated, certain limitations remain, particularly in optically complex Case-2 waters such as the Gulf waters. Shallow depth, high turbidity, and elevated concentration of optically active substances can reduce algorithm accuracy and lead to greater uncertainties in nearshore Chl-a retrievals. The lack of long-term in-situ measurements further limits validation at finer spatial and temporal scales. The global distribution of random uncertainty in satellite-derived Chl-a varies significantly by region, with higher values typically observed in coastal and turbid waters due to challenges in atmospheric correction and the separating water-leaving radiance from atmospheric signals. While the general global patterns of Chl-a are relatively consistent across different satellite products, regional uncertainty differences remain substantial, highlighting the need for region-specific approaches to uncertainty assessment.

Nevertheless, global ocean color products, including the OC-CCI dataset used in this study, have been undergoing continuous improvements in merging techniques, atmospheric correction, and algorithm refinements (Sathyendranath et al., 2019). The OC-CCI dataset, in particular, has been shown in both global and Gulf-specific studies to provide consistent and reliable estimates, outperforming earlier satellite products in terms of stability and cross-sensor consistency (Belo Couto et al., 2016; Al-Naimi et al., 2017; Moradi, 2021). Uncertainty in OC-CCI Chl-a data arises from factors such as sensor noise, optical classification limits, differences in atmospheric correction and algorithm performance across various water types and missions, and spatiotemporal scale mismatches (Mélin et al., 2016; Brewin et al., 2017; IOCCG, 2019). To address this, OC-CCI provides per-pixel uncertainty estimates based on an improved optical classification scheme and blending algorithms. The development of a new, improved optical classification scheme based on a larger dataset of satellite-derived spectra, and the use of class membership to assess the performance of multiple algorithms within each water class and to blend algorithms have further enhanced the assignment of uncertainties, and improved the quality of the final product, especially for Chl-a (Jackson et al., 2022). These developments make OC-CCI a robust resource for long-term Chl-a studies, while the acknowledged limitations highlight the need for complementary in-situ observations and tailored regional algorithms. Moreover, the use of latest version (Version 6.0) in this study also ensures improvements to the algorithm blending and error correction to provide more reliable results.