Coastal phytoplankton blooms expand and intensify in the 21st century

Data sources

MODIS on the Aqua satellite provides a global coverage within 1–2 days. All images acquired by this satellite mission from January 2003 to December 2020 were used in our study to detect global coastal phytoplankton blooms, with a total of 0.76 million images. MODIS Level-1A images were downloaded from the Ocean Biology Distributed Active Archive Center (OB.DAAC) at NASA Goddard Space Flight Center (GSFC), and were subsequently processed with SeaDAS software (version 7.5) to obtain Rayleigh-corrected reflectance (R_rc (dimensionless), which was converted using the rhos (in sr⁻¹) product (rhos × π) from SeaDAS)⁴¹, remote sensing reflectance (R_rs (sr⁻¹)) and quality control flags (l2_flags). If a pixel was flagged by any of the following, it was then removed from phytoplankton bloom detection: straylight, cloud, land, high sunglint, high solar zenith angle and high sensor zenith angle (https://oceancolor.gsfc.nasa.gov/atbd/ocl2flags/). MODIS level-3 product for aerosol optical thicknesses (AOT) at 869 nm was also obtained from OB.DAAC NASA GSFC (version R2018.0), which was used to examine the impacts of aerosols on bloom trends.

We examined the algal blooms in the EEZs of 153 ocean-bordering countries (excluding the EEZs in the Caspian Sea or around the Antarctic), 126 of which were found with at least one bloom in the past two decades. The EEZ dataset is available at https://www.marineregions.org/download_file.php?name=World_EEZ_v11_20191118.zip. The EEZs are up to 200 nautical miles (or 370 km) away from coastlines, which include all continental shelf areas and offer the majority of marine resources available for human use. Regional statistics of algal blooms were also performed for LMEs. LMEs encompass global coastal oceans and outer edges of coastal currents areas, which are defined by various distinct features of the oceans, including hydrology, productivity, bathymetry and trophically dependent populations⁴². Of the 66 LMEs identified globally, we excluded the Arctic and Antarctic regions and examined 54 LMEs. The boundaries of LMEs were obtained from https://www.sciencebase.gov/catalog/item/55c77722e4b08400b1fd8244.

We used HAEDAT to validate our satellite-detected phytoplankton blooms in terms of presence or absence. The HAEDAT dataset (http://haedat.iode.org) is a collection of records of HAB events, maintained under the UNESCO Intergovernmental Oceanographic Commission and with data archives since 1985. For each HAB event, the HAEDAT records its bloom period (ranging from days to months) and geolocation. We merged duplicate entries when both the recorded locations and times of the HAEDAT events were very similar to one another, and a total number of 2,609 HAEDAT events were ultimately selected between 2003 and 2020.

We used the ¼° resolution National Oceanic and Atmospheric Administration Optimum Interpolated SST (v. 2.1) data to examine the potential simulating effects of warming on the global phytoplankton trends. We also estimated the SST gradients following the method of Martínez-Moreno³³. As detailed in ref. ³³, the SST gradient can be used as a proxy for the magnitude of oceanic mesoscale currents (EKE). We used the SST gradient to explore the effects of ocean circulation dynamics on algal blooms.

Fertilizer uses and aquaculture production for different countries was used to examine the potential effects of nutrient enrichment from humans on global phytoplankton bloom trends. Annual data between 2003 and 2019 on synthetic fertilizer use, including nitrogen and phosphorus, are available from https://ourworldindata.org/fertilizers. Annual aquaculture production includes cultivated fish and crustaceans in marine and inland waters, and sea tanks, and the data between 2003 and 2018 are available from https://ourworldindata.org/grapher/aquaculture-farmed-fish-production.

The MEI, which combines various oceanic and atmospheric variables³⁶, was used to examine the connections between El Niño–Southern Oscillation activities and marine phytoplankton blooms. The dataset is available from https://psl.noaa.gov/enso/mei/.

Development of an automated bloom detection method

A recent study by the UNESCO Intergovernmental Oceanographic Commission revealed that globally reported HAB events have increased⁶. However, such an overall increasing trend was found to be highly correlated with recently intensified sampling efforts⁶. Once this potential bias was accounted for by examining the ratio between HAB events to the number of samplings⁵, there was no significant global trend in HAB incidence, though there were increases in certain regions. With synoptic, frequent, and large-scale observations, satellite remote sensing has been extensively used to monitor algal blooms in oceanic environments^17,18,19. For example, chlorophyll a (Chla) concentrations, a proxy for phytoplankton biomass, has been provided as a standard product by NASA since the proof-of-concept Coastal Zone Color Scanner (1978–1986) era^43,44. The current default algorithm used to retrieve Chla products is based on the high absorption of Chla at the blue band^45,46, which often shows high accuracy in the clear open oceans but high uncertainties in coastal waters. This is because, in productive and dynamic coastal oceans, the absorption of Chla in the blue band can be obscured by the presence of suspended sediments and/or coloured dissolved organic matter (CDOM)⁴⁷. To address this problem, various regionalized Chla algorithms have been developed⁴⁸. Unfortunately, the concentrations of the water constituents (CDOM, sediment and Chla) can vary substantially across different coastal oceans. As a result, a universal Chla algorithm that can accurately estimate Chla concentrations in global coastal oceans is not currently available.

Alternatively, many spectral indices have been developed to identify phytoplankton blooms instead of quantifying their bloom biomass, including the normalized fluorescence line height²¹ (nFLH), red tide index⁴⁹ (RI), algal bloom index⁴⁷ (ABI), red–blue difference (RBD)⁵⁰, Karenia brevis bloom index⁵⁰ (KBBI) and red tide detection index⁵¹ (RDI). In practice, the most important task for these index-based algorithms is to determine their optimal thresholds for bloom classification. However, such optimal thresholds can be regional-or image-specific²⁰, due to the complexity of optical features in coastal waters and/or the contamination of unfavourable observational conditions (such as thick aerosols, thin clouds, and so on), making it difficult to apply spectral-index-based algorithms at a global scale.

To circumvent the difficulty in determining unified thresholds for various spectral indices across global coastal oceans, an approach from a recent study to classify algal blooms in freshwater lakes⁵² was adopted and modified here. In that study, the remotely sensed reflectance data in three visible bands (red, green and blue) were converted into two-dimensional colour space created by the Commission Internationale del’éclairage (CIE), in which the position on the CIE chromaticity diagram represented the colour perceived by human eyes (Extended Data Fig. 1a). As the algal blooms in freshwater lakes were manifested as greenish colours, the reflectance of bloom-containing pixels was expected to be distributed in the green gamut of the CIE chromaticity diagram; the stronger the bloom, the closer the distance to the upper border of the diagram (the greener the water).

Here, the colour of phytoplankton blooms in the coastal oceans can be greenish, yellowish, brownish, or even reddish⁵³, owing to the compositions of bloom species (diatoms or dinoflagellates) and the concentrations of different water constituents. Furthermore, the Chla concentrations of the coastal blooms are typically lower than those in inland waters, thus demanding more accurate classification algorithms. Thus, the algorithm proposed by Hou et al.⁵² was modified when using the CIE chromaticity space for bloom detection in marine environments. Specifically, we used the following coordinate conversion formulas to obtain the xy coordinate values in the CIE colour space:

where R, G and B represent the R_rc at 748 nm, 678 nm (fluorescence band) and 667 nm in the MODIS Aqua data, respectively. By contrast, the R, G and B channels used in Hou et al.⁵² were the red, green and blue bands. We used the fluorescence band for the G channel because, for a given region, the 678 nm signal increases monotonically with the Chla concentration for blooms of moderate intensity²¹, which is similar to the response of greenness to freshwater algal blooms. Thus, the converted y value in the CIE coordinate system represents the strength of the fluorescence. In practice, for pixels with phytoplankton blooms, the converted colours in the chromaticity diagram will be located within the green, yellow or orange–red gamut (see Extended Data Fig. 1a); the stronger the fluorescence signal is, the closer the distance to the upper border of the CIE diagram (larger y value). By contrast, for bloom-free pixels without a fluorescence signal, their converted xy coordinates will be located in the blue or purple gamut. Therefore, we can determine a lower boundary in the CIE two-dimensional coordinate system to separate bloom and non-bloom pixels, similar to the method proposed by Hou et al.⁵².

We selected 53,820 bloom-containing pixels from the MODIS R_rc data as training samples to determine the boundary of the CIE colour space. These sample points were selected from nearshore waters worldwide where frequent phytoplankton blooms have been reported (Extended Data Fig. 2); the algal species included various species of dinoflagellates and diatoms²⁰. A total of 80 images was used, which were acquired from different seasons and across various bloom magnitudes, to ensure that the samples used could almost exhaustively represent the different bloom conditions in the coastal oceans.

We combined the MODIS FLH_Rrc (fluorescence line height based on R_rc) and enhanced red–green–blue composite (ERGB) to delineate bloom pixels manually. The FLH_Rrc image was calculated as:

where R_rc667, R_rc678 and R_rc748 are the R_rc at 667, 678 and 748 nm, respectively, and F₆₆₇, F₆₇₈ and F₇₄₈ are the corresponding extraterrestrial solar irradiance. ERGB composite images were generated using R_rc of three bands at 555 (R), 488 (G) and 443 nm (B). Although phytoplankton-rich and sediment-rich waters have high FLH_Rrc values, they appear as darkish and bright features in the ERGB images (Extended Data Fig. 3), respectively²¹. In fact, visual examination with fluorescence signals and ERGB has been widely accepted as a practical way to delineate coastal algal blooms on a limited number of images^21,54,55. Note that the FLH_Rrc here was slightly different from the NASA standard nFLH product⁵⁶, as the latter is generated using R_rs (corrected for both Rayleigh and aerosol scattering) instead of R_rc (with residual effects of aerosols). However, when using the NASA standard algorithm to further perform aerosol scattering correction over R_rc, 20.7% of our selected bloom-containing pixels failed to obtain valid R_rs (without retrievals or flagged as low quality), especially for those with strong blooms (see examples in Extended Data Fig. 4). Likewise, we also found various nearshore regions with invalid R_rs retrievals. By contrast, R_rc had valid data for all selected samples and showed more coverage in nearshore coastal waters. The differences between R_rs and R_rc were because the assumptions for the standard atmospheric correction algorithm do not hold for bloom pixels or nearshore waters with complex optical properties⁵⁷. In fact, R_rc has been used as an alternative to R_rs in various applications in complex waters^58,59.

We converted the R_rc data of 53,820 selected sample pixels into the xy coordinates in the CIE colour space (Extended Data Fig. 1a). As expected, these samples of bloom-containing pixels were located in the upper half of the chromaticity diagram (the green, yellow and orange–red gamut) (Extended Data Fig. 1a). We determined the lower boundary of these sample points in the chromaticity diagram, which represents the lightest colour and thus the weakest phytoplankton blooms; any point that falls above this boundary represents stronger blooms. The method to determine the boundary was similar to Hou et al.⁵²: we first binned the sample points according to the x value in the chromaticity diagram and estimated the 1st percentile (Q₁%) of the corresponding Y for each bin; then, we fit the Q₁% using two-order polynomial regression. Sensitivity analysis with Q_0.3% (the three-sigma value) resulted in minor changes (<1%) in the resulting bloom areas for single images. Notably, sample points were rarely located near white points (x = 1/3 and y = 1/3, represent equal reflection from three RGB bands) in the diagram, and we used two polynomial regressions to determine the boundaries for x values greater and less than 1/3, which can be expressed as:

Based on the above, if a pixel’s xy coordinate (converted from R_rc spectrum) satisfies the conditions of (x < 1/3 AND y > y₁) or (x > 1/3 AND y > y₂), it is classified as a ‘bloom’ pixel.

Depending on the local region and application purpose, the meaning of ‘phytoplankton bloom’ may differ. Here, for a global application, the pixelwise bloom classification is based on the relationship (represented using the CIE colour space) between R_rc in the 667-, 678- and 754-nm bands derived from visual interpretation of the 80 pairs of FLH_Rrc and ERGB imagery. Instead of a simple threshold, we used a lower boundary of the sample points in the chromaticity diagram to define a bloom. In simple words, a pixel is classified as a bloom if its fluorescence signal is detectable (the associated xy coordinate in the CIE colour space located above the lower boundary). Histogram of the nFLH values from the 53,820 training pixels demonstrated the minimum value of ~0.02 mW cm⁻² μm⁻¹ (Extended Data Fig. 1a), which is in line with the lower-bound signal of K. brevis blooms on the West Florida shelf^21,47. Note that, such a minimum nFLH is determined from the global training pixels, and it does not necessarily represent a unified lower bound for phytoplankton blooms across the entire globe, especially considering that fluorescence efficiency may be a large variable across different regions. Different regions may have different lower bounds of nFLH to define a bloom, and such variability is represented by the predefined boundary in the CIE chromaticity diagram in our study. Correspondingly, although the accuracy of Chla retrievals may have large uncertainties in coastal waters, the histogram of the 53,820 training pixels shows a lower bound of ~1 mg m⁻³ (Extended Data Fig. 1a). Similarly to nFLH, such a lower bound may not be applicable to all coastal regions, as different regions may have different lower bounds of Chla for bloom definition.

Although the MODIS cloud (generated by SeaDAS with R_rc869 < 0.027) and associated straylight flags can be used to exclude most clouds, we found that residual errors from thin clouds or cloud shadows could affect the spectral shape and cause misclassification for bloom detections. Thus, we designed two spectral indices to remove such effects:

where Index1 and Index2 were used to remove cloud shadows and clouds, respectively. The nR_rc443, nR_rc488 and nR_rc555 in index1 are the normalized R_rc, obtained by normalizing R_rc488. Similar calculations were performed for index2. The purpose of normalizations is to eliminate the effect of the absolute magnitude of the reflectance, so that the thresholds of these two indices are influenced by only the relative magnitude (spectral shape). We determined thresholds for Index1 (>0.12) and Index2 (<0.012) through trial-and-error and ensured that the misclassifications caused by residual errors from clouds and cloud shadows could be effectively removed. After applying the cloud/cloud shadow and various other masks that are associated with l2_flags, we obtained an annual mean valid pixel observation (N_vobs) of ~2.0 × 10⁵ for global 1° × 1° grid cells, and the fluctuation patterns and trends of N_vobs, either annually or seasonally, are different from that of the global bloom frequency and affected areas (see Supplementary Fig. 1).

Assessments of the algorithm performance

In addition to phytoplankton blooms, macroalgal blooms (Sargassum and Ulva) frequently occur in many coastal oceans^60,61,62,63. To verify whether our CIE-fluorescence algorithm could eliminate such impacts, we compared the spectra between micro-and macroalgal blooms (see Extended Data Fig. 1b). We found that the spectral shapes of Sargassum and Ulva are substantially different from those of microalgae, particularly for the three bands used for CIE coordinate conversion. The converted xy coordinates for macroalgae were located in the purple–red gamut of the CIE diagram, which was far below the predefined boundary (Extended Data Fig. 1). Moreover, our algorithm is not affected by highly turbid waters for the following two reasons: first, extremely high turbidity tends to saturate the MODIS ocean bands⁶⁴, which can be easily excluded; second, without a fluorescence peak, the reflectance of unsaturated turbid waters, after conversion to CIE coordinates, will be located below the predefined boundary (see example in Extended Data Fig. 1b). We also confirmed that the spectral shapes of coccolithophore blooms are different from dinoflagellates and diatoms (see example in Extended Data Fig. 1b), and thus they are excluded from our algorithm.

Three different types of validation methods were adopted to demonstrate the reliability of the proposed CIE-fluorescence algorithm for phytoplankton bloom detection in global coastal oceans, including visual inspections of the RGB, ERGB and FLH_Rrc images, verifications using independent manually delineated algal blooms, and comparisons with the reported HAB events from the HAEDAT dataset.

First, we selected MODIS Aqua images from different locations where coastal phytoplankton blooms have been recorded in the published literature. We visually compared the RGB, ERGB, and FLH_Rrc images, and our algorithm detected bloom patterns (see examples in Extended Data Fig. 3). Comparisons from various images worldwide showed that our algorithm could successfully identify regions with high FLH_Rrc values and brownish-to-darkish features on the ERGB images, which can be considered phytoplankton blooms.

Second, we delineated additional 15,466 bloom-containing pixels from 35 images covering different coastal areas, using the same visual inspection and manual delineation method as for the training sample pixels. Moreover, we also selected 14,149 bloom-free pixels (bright or turquoise green colours on ERGB images or low FLH_Rrc values) within the same images as bloom-containing images. We applied our algorithm to all these pixels, and compared the algorithm-identified and manually delineated results. Our CIE-fluorescence algorithm showed high values in both producer and user accuracies (92.04% and 98.63%) (Supplementary Table 1), and appeared effective at identifying bloom pixels and excluding false negatives (blooms classified as non-blooms) and false positives (non-blooms classified as blooms).

Third, we validated the satellite-detected phytoplankton blooms using in situ reported HAB events from the HAEDAT dataset. For each HAB event in the HAEDAT dataset, we obtained all MODIS images over the reported bloom period (from days to months). Within each year, we estimated the ratio between the number of satellite images with ‘bloom detected’ against the number of valid images (see definition above) during the bloom periods across the entire globe (Supplementary Table 1). Moreover, we calculated the number of events with at least one successful satellite bloom detection (N_s), and then estimated the ratio between N_s and the total HAB events for each year. Results showed that substantial amounts (averaged at 51.2%) of satellite observations acquired during the HAB event periods were found with phytoplankton blooms by our algorithm. Overall, 79.3% of the global HAB events were successfully identified by satellite. The discrepancies between satellite and in situ observations could be explained by the following reasons: first, our study focused only on the phytoplankton blooms that are resolvable by satellite fluorescence signals; other types of HAB events in the HAEDAT dataset may not have been detectable by satellite observations, such as events with lower phytoplankton biomass but high toxicity, occurrences at the subsurface layers, or fluorescence signals overwhelmed by suspended sediments^65,66,67. Second, although the HAEDAT recorded HAB events could be sustained for long periods, high biomass of surface algae may not have occurred every day within this period due to the changes in stratification, precipitation, wind, vertical migration of cells, and many other factors⁶⁸. Third, the spatial scale of certain HAB events may have been too small to be identified using the 1-km resolution MODIS observations. Fourth, a reduced MODIS satellite observation frequency by the contaminations of clouds and land adjacency effects⁶⁹. Therefore, we believe the underestimations of satellite-detected blooms compared to the in situ reported HAB events were mainly due to inconsistencies between the two observations rather than uncertainties in our algorithm.

Because R_rc depends not only on water colour but also on aerosols (type and concentration) and solar and viewing geometry, a sensitivity analysis was used to determine whether such variables could impact bloom detection. Aerosol reflectance (ρ_a) with different AOTs at 869 nm was simulated using the NASA-recommended maritime aerosol model (r75f02, with a relative humidity of 75% and a fine-mode fraction of 2%). Then, ρ_a of each MODIS band was added to R_rc images, and the resulting bloom areas with and without added ρ_a were compared. Results showed that even with a change of 0.02 in AOT at 869 nm, the bloom areas showed minor changes (<2%) in the tested images; minor changes were also found when we used different aerosol models to conduct ρ_a simulations⁷⁰. Note that 0.02 represents the high end of the AOT intra-annual variability in coastal oceans (see Extended Data Fig. 5), and the associated interannual changes are much smaller. Thus, the use of R_rc instead of the fully atmospherically corrected reflectance R_rs could have limited impacts on our detected global bloom trend.

We also tried various index-based algorithms developed in previous studies. However, results showed that all these methods require image-specific thresholds to accurately determine algal bloom boundaries for different coastal regions (see Extended Data Fig. 6). By contrast, although our CIE-fluorescence algorithm may lead to different bloom thresholds for different regions, it can identify bloom pixels without adjusting the coefficients and, therefore, is more suitable for global-scale bloom assessment efforts.

We acknowledge that our satellite-detected algal blooms represent only high amounts of phytoplankton biomass on the ocean surfaces without distinguishing whether such blooms produce toxins or are harmful to marine environments. Furthermore, with only limited spectral information from MODIS, it is difficult to discriminate the phytoplankton species of algal blooms; such information could help to improve our understanding of the impacts of these phytoplankton blooms. However, we expect these challenges to be addressed soon with the scheduled launch of the Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission by NASA in 2024, where the hyperspectral measurements over a broad spectrum (350–885 nm) will make species-level classifications possible⁷¹.

Exploring the patterns and trends of global coastal phytoplankton blooms

We applied the CIE-fluorescence algorithm to all MODIS Aqua level-2 R_rc images, and a total number of 0.76 million images between 2003 and 2020 were processed. We mapped all detected blooms into 1-km daily scale level-3 composites. The number of bloom counts within a year for each location can be easily enumerated, and the long-term annual mean values were then estimated (Fig. 1a). We further calculated the total global bloom-affected area (the areas where algal blooms were detected at least once) for each year and examined their changes over time (Fig. 2b).

We defined bloom frequency (dimensionless) to represent the density of phytoplankton blooms for a year by integrating the bloom count and bloom-affected areas within 1°×1° grid cells within that year, which is expressed as:

where M_i is the enumerated bloom count for each 1-km resolution pixel in a year within one 1° × 1° grid cell, and n represents the associated number of bloom-affected pixels in the same cell (the number of pixels with M_i > 0), and N is the total number of 1-km MODIS pixels in this grid cell. We estimated the bloom frequency for each year between 2003 and 2020, and determined the long-term trend over global EEZs through a linear least-squares regression (see Fig. 2a).

Continental and country-level statistics were performed for bloom count, bloom-affected areas, and bloom frequency (Fig. 1b,c and Supplementary Table 2), using boundaries for the EEZs of different ocean-bordering countries (see above). Similar statistics were also conducted for 54 LMEs (Extended Data Fig. 7 and Supplementary Table 3).

Correlations with SST and SST gradient

To assess the impacts of climate change on long-term trends in coastal phytoplankton blooms, we correlated the annual mean bloom frequency and the associated SST and SST gradient in various coastal current systems for grid cells with significant changes in bloom frequency (Fig. 3c). The SST and SST gradient were averaged over the growth window within a year, assuming that the changes within the growth window, either in water temperatures or ocean circulations, play more important roles in the bloom trends compared to other seasons³².

We determined the growth window of phytoplankton blooms for each 1° × 1° grid cell (Extended Data Fig. 9a) using the following method: first, we estimated the proportion of cumulative bloom-affected pixels within the grid cells for a year. Second, a generalized additive model⁷² was used to determine the shape of the phenological curves (Extended Data Fig. 9b), where a log link function and a cubic cyclic regression spline smoother were applied^73,74. Third, the timing of maximum bloom-affected areas (TMBAA) was then determined by identifying the inflection point on the bloom growth curve (Extended Data Fig. 9c). To facilitate comparisons across Northern and Southern Hemispheres, the year in the Southern Hemisphere was shifted forward by 183 days (Extended Data Fig. 9c). We characterized the similarity of the bloom growth curve between different grid cells and grouped them into three distinct clusters using a fuzzy c-means cluster analysis method^75,76. We found uniform distributions of the clusters over large geographic areas. Cluster I is mainly distributed in mid-low latitudes (<45° N and <30° S), where the maximum bloom-affected areas were expected in the early period of the year. Cluster II was mostly found in higher latitudes, with bloom developments (quasi-) synchronized with increases in SST. Cluster III was detected along the coastlines, where the bloom-affected areas increase throughout the entire year. In practice, the growth window for clusters I and III was set as the entire year, and that for cluster II was set from day 150 to day 270 within the year. We further found that the TMBAA for cluster II showed small changes over the entire period (Extended Data Fig. 9d), indicating relatively stable phenological cycles for those phytoplankton blooms^32,77.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com