The Northern Bering and Chukchi seas, located in the Pacific Arctic, are productive ocean areas supporting important commercial and subsistence fisheries and serving as a crucial habitat for sea-icedependent species such as polar bears and walruses (National Research Council, 1996). This region is experiencing rapid environmental shifts due to climate change, affecting both physical and biological components of the ecosystem, with warming temperatures leading to reductions in sea-ice extent and alterations in the marine food web (Huntington et al., 2020;Grebmeier et al., 2015), with significant implications for marine and human lives (Stabeno et al., 2019).
Ocean warming, and associated ecosystem changes, are leading to an increase in frequency and spatial extent of harmful algal bloom (HAB) events, particularly of the dinoflagellate Alexandrium catenella (Anderson et al., 2021(Anderson et al., , 2022;;Natsuike et al., 2017). Note that this species has also been previously identified as Gonyaulax catenella, G. tamarensis, A. tamarense, A. tamarense Group I, and A. fundyense (Lilly et al., 2007;Brosnahan et al., 2020;John et al., 2014;Litaker et al., 2018). Here, we follow the taxonomic revision proposed by John et al. (2014) that discriminates this species based on modern understanding its molecular phylogeny, toxicity, and mating barriers rather than the flawed morphological criteria that were applied historically. As A. catenella blooms terminate, cysts are formed that settle to bottom sediments. These cysts then initiate new blooms through germination, which happens earlier and progresses more rapidly at warmer bottom temperatures (Fischer and Brosnahan, 2022;Anderson et al., 2021). Moreover, the warming of surface waters drives temperature-driven water column stratification, and also favours earlier ice melt, further enhancing salinity-driven stratification. This stratification facilitates the development, maintenance and expansion of dinoflagellate blooms (Margalef, 1978(Margalef, , 1997)).
A. catenella produces paralytic shellfish toxins (including saxitoxin and analogues) which bind to voltage-gated sodium channels of neurons and muscle cells, blocking synaptic transmission and causing paralysis (Gad, 2014). These toxins accumulate in shellfish, zooplankton and small fishes, and have been found widely in marine birds, invertebrates, and mammals throughout the Pacific Arctic region, posing a threat to marine ecosystems and human health (Lefebvre et al., 2016).
In the Bering and Chukchi seas, blooms of Alexandrium catenella (A. catenella) exhibit complex dynamics influenced by a combination of environmental conditions such as temperature and salinity (which affect stratification and light availability to phytoplankton cells), and nutrient availability and stoichiometry (Anderson et al., 2021;Tobin et al., 2019). In other regions, blooms of Alexandrium species commonly follow or co-occur with diatom blooms, benefiting from post-bloom nitrogen depletion (Collos et al., 2014;Love et al., 2005), are associated with water column stratification driven by the warming of surface waters and reduced salinity from ice melt at higher latitudes (Klemm et al., 2022;Devred et al., 2018), or by terrestrial freshwater input and coastal upwelling at lower latitudes (Xie et al., 2022;Barría et al., 2022;Crawford et al., 2021).
As a consequence of the present trend in ocean warming, Alexandrium spp. blooms are expanding in the Arctic (Fachon et al., 2025;Anderson et al., 2021;Boivin-Rioux et al., 2021) and in other ecosystems across the globe (Sacilotto Detoni et al., 2023;Barría et al., 2022;Klemm et al., 2022;Xie et al., 2022;Crawford et al., 2021;Brosnahan et al., 2020;Devred et al., 2018;Collos et al., 2014;Gettings et al., 2014;Townsend et al., 2005). This reveals a need for improved monitoring tools and continued research to better understand and mitigate the impacts of Alexandrium blooms (Fachon et al., 2025;Anderson et al., 2022).
At present, HAB detection in the Northern Bering and Chukchi seas is largely based on in situ sampling and analysis of phytoplankton community composition and abundance in research surveys, with some sporadic collection of marine invertebrates for toxin testing. Some regional coastal communities implement additional nearshore sampling for early warning on a primarily volunteer basis. Data are also obtained from ships of opportunity transiting the region for other research or commercial purposes. Such direct measurements are labor intensive and intermittent in time and space, making it difficult to detect and track the rapid development and movement of HABs over larger spatial areas. These challenges are further exacerbated by the unique geography of the region, where large stretches of sparsely inhabited coastline are nearly impossible to access for sampling with in situ methods. In that sense, remote sensing tools such as satellite imagery can be an extremely helpful complement, providing observations with higher frequency, synoptic assessments of near-surface phytoplankton communities, and insights on their transport over vast geographic areas.
Remote sensing models estimate the biomass of phytoplankton taxonomic groups from apparent optical properties such as remotesensing reflectance (R rs (λ)), which is measured using satellite radiometric sensors. These models leverage the unique optical properties (absorption and backscattering) of phytoplankton by aligning features in the R rs (λ) spectrum with the specific spectral signature of electromagnetic energy absorption by different phytoplankton pigments (a ph (λ)), associated to their particulate backscatter properties (b bp (λ)). Yet, despite the fact that Alexandrium and similar peridinin-rich dinoflagellates exhibit specific spectral features due to their high pigment content, and low b bp (λ) (Lomas et al., 2024), early efforts to detect and track their blooms with satellites did not leverage ocean color data in its entirety. Instead, satellite-derived chlorophyll concentration estimates were associated with ship-based cell abundance measurements to track blooms and/or predict their development from satellite sea surface temperature measurements (Li et al., 2020;Condie et al., 2019;Guallar et al., 2017;Luerssen et al., 2005). Recently, Sacilotto Detoni et al. ( 2023) developed a band-ratio model to use ocean color data to detect Alexandrium blooms in the Galician Rías. Although other band-ratio models such as the Normalized difference chlorophyll index -NDCI (Mishra and Mishra, 2012) and fluorescence line height -FLH (Behrenfeld et al., 2008) are often used to detect phytoplankton blooms in general, specific satellite algorithms to enhance Alexandrium bloom monitoring in the Bering and Chukchi Seas have never been parameterized or evaluated.
It is also the case that other peridinin dinoflagellates (Scrippsiella, Pentapharsodinium, and others) can be abundant across this region as evidenced by the presence their cysts in underlying sediments (Fachon and Anderson, 2021;Gu et al., 2013;Orlova and Morozova, 2013). Because these cells present similar optical properties to that of Alexandrium (Zapata et al., 2012), an assessment of bloom taxonomic composition is essential for interpreting outputs of satellite algorithms. Moreover, satellites observe the ocean color signal generated from a depth of only a few metres below the ocean surface, and are unable to detect subsurface phytoplankton blooms in most cases (IOCCG et al., 2021) except when there is a clear relationship between subsurface and surface phytoplankton communities (Brewin et al., 2022;Lange et al., 2020). Therefore, the combination of in situ and satellite observations of phytoplankton blooms dramatically enhances the spatial and temporal coverage of HAB monitoring, providing insights on bloom initiation, transport, and fate (e.g., Carvalho et al., 2010;Fick et al., 2024).
This study leverages the optical properties of Alexandrium spp. to develop novel combinations of ocean color algorithms, regionally-tuned to the Bering and Chukchi seas, for detecting and quantifying Alexandrium catenella blooms using R rs (λ) data from Sentinel-3's Ocean and Land Colour Instrument (OLCI). Specifically, our objectives are: 1) to develop and validate specific color index models for detecting dinoflagellate blooms and differentiating these cells from other phytoplankton groups that form blooms in the region; 2) to define the regional parameterization of the color index models and a multivariate spectral analysis model (vector decomposition) to quantify A. catenella in dinoflagellate blooms of the Bering and Chukchi seas; 3) to compare the performance of novel color indexes and vector decomposition algorithms and legacy algorithms for detecting and quantifying A. catenella blooms in the Bering/Chukchi region; 4) to implement the novel and legacy dinoflagellate algorithms to generate maps of A. catenella cell abundance using Sentinel-3 remote sensing reflectance data of several phytoplankton blooms in the Bering and Chukchi seas. Finally, using shipboard measurements of A. catenella abundances collected during summers 2023 and 2024, we demonstrate that the developed dinoflagellate algorithms can provide near-real time detection of dinoflagellate blooms in the Northern Bering Sea, Chukchi Sea, and Gulf of Alaska regions.
The Northern Bering Sea and Chukchi Sea (Fig. 1a) are located between the US State of Alaska and Eastern Russia, and are characterized by broad shallow continental shelves influenced by the dynamic interplay of various currents. The water flow on the Northern Bering and Chukchi shelves is generally northward. There are three well-defined currents on the eastern Bering Sea shelf: the Alaskan Coastal Current (ACC); the Bering Shelf Current (BSC); and the Anadyr Current (AC). The ACC is a buoyancy driven coastal flow that generally follows the 50-m isobath, although it shoals north of 60 • N. It originates in shallow eastern Aleutian Passes, in particular Unimak Pass. The BSC is a northward flow that follows the 100-m isobath between the outer and middle shelf domains and transports Bering shelf water northward. The source of this water is on shelf flow of slope water primarily but not exclusively through canyons (Stabeno and Van Meurs, 1999;Stabeno et al., 2016). The AC originates in the north, and is the strongest of the currents. These three currents combine to form Bering Strait throughflow (Woodgate and Aagaard, 2005), with the AC providing saltiest, nutrient-rich waters and the ACC the freshest. This system transports significant heat and nutrients northward into the Chukchi Sea. Within the Chukchi Sea this inflow intersects the Siberian Coastal Current (SCC) which transports cold, low salinity water eastward. The SCC is also known as the East Siberian Current. Although the SCC enters the Bering Strait along Chukotka, it typically turns northward and joins the flow through Bering Strait (Weingartner et al., 1999;Stabeno et al., 2018).
Alexandrium catenella abundance data used for initial algorithm development were collected by a shipboard Imaging FlowCytobot (IFCB, McLane Laboratories) deployed aboard the R/V Norseman II during a two-part HAB-focused research cruise that took place during July-September 2022 (NRS2022_01S, NRS2022_02S; Fachon et al., 2025). This instrument is an automated imaging flow cytometer that continuously captures submicron resolution images of individual phytoplankton and other particles ranging between 5 and 150 μm in length (Olson and Sosik, 2007), as shown in Fig. S2 of the Supplementary Material. For these cruises, the IFCB was configured to collect samples from an underway seawater line (intake depth ~3 m), with image capture triggered by chlorophyll fluorescence, processing 5 mL every 20-30 min throughout the two-month sampling effort. During this time, a bloom of A. catenella developed in the study region, ultimately extending more than 600 km between the northern Bering Sea and Chukchi Shelf. This event was observed multiple times as the vessel travelled back and forth through the region, providing the opportunity to track bloom progression between late July and early September (Fachon et al., 2025).
Images were classified using a convolutional neural network (CNN) classifier, which identified 94 taxonomic groups, mostly to genus level, including Alexandrium (https://github.com/WHOIGit/ifcb_classifier). Following this initial step, Alexandrium classifications were reviewed and corrected manually to produce final measures of cell concentration at each sampling point (n = 2680). The resulting dataset is available at https://doi.org/10.18739/A2804XM7S. Additional seawater samples collected during the cruise were used for microscopy analyses that confirmed the species ID as A. catenella.
Additional shipboard observations were conducted in the summers of 2023 and 2024 and used to test and refine model outputs. For this, we look at underway IFCB imagery collected during cruises of opportunity aboard the R/V Sikuliaq (SKQ2023_10S, June 16 -July 30, 2023; SKQ2024_14S, August 6-27 2024). Once again, the IFCB was configured to sample from the underway seawater line (intake depth ~6m). Images were again classified through initial processing with a CNN classifier model, then manually reviewed to produce final A. catenella cell concentration measurements. In late July 2023 during SKQ2023_10S, an A. catenella bloom was detected in the northern Bering Sea extending into the Bering Strait. In 2024, no Alexandrium blooms were detected by shipboard observing efforts in the northern Bering Sea, Bering Strait, or Chukchi Shelf, however a dense diatom bloom was recorded in the northern Bering Sea in early August.
The Sentinel-3's Ocean and Land Colour Instrument (OLCI) has a spatial resolution (pixel size) of 300 m and provides global coverage on a daily or bi-daily basis depending on location. At high latitudes during summer, up to 8 images can be provided for the same location in one day. For all sampling dates, radiometric data (radiance at the top of atmosphere (TOA) captured by Sentinel-3's OLCI were downloaded as Level-1 Earth Full Resolution (OL_1_EFR) from the former Copernicus Open Access Hub (now Copernicus Data Space Ecosystem: https://data space.copernicus.eu), hosted by the European Space Agency. Level-1 TOA radiance images were then converted into sea-level remote sensing reflectance (R rs (λ)) images using the ACOLITE processor through a Python routine. ACOLITE is freely available software for atmospheric correction and processing of satellite imagery, with a set of radiative transfer calculations developed for marine, coastal and inland waters by the Royal Belgian Institute of Natural Sciences (RBINS) (Vanhellemont and Ruddick, 2021), where settings can be adjusted for the specific area of interest.
Once calculated for the whole data file, R rs (λ) at 12 wavelength bands (400, 412, 443, 490, 510, 560, 620, 665, 674, 682, 709 and 754 nm) were extracted from images acquired at the same dates as the IFCB sampling dates. All Sentinel-3 pixels within a maximum distance of 1 km from the exact sampling location (latitude and longitude) were extracted. Within this distance radius, there were 3401 satellite match-ups at 218 in situ sampling sites (more than one satellite image matches one in situ sampling site). When more than one image or more than one pixel matched with the sampling location and date, the mean R rs (λ) value for each wavelength (mean spectrum) was calculated, generating one match up for each of the 218 sampling locations. A. catenella was present in samples at 108 of the 218 match-up sites, with cell abundances higher than 1000 cells/L at 64 locations and a maximum match-up cell abundance of 94,946 cells/L. The number of observations used for model parameterization was reduced to 45 or less after image masks using bloom criteria (FLH > 0.0001, see section 2.4.4 and Table 2 for details) were applied.
Remote sensing reflectance (R rs (λ)) varies as a function of the absorption and backscatter of optically-active water constituents, including phytoplankton. Phytoplankton influences the particulate absorption (a p (λ)) and particulate backscatter (b bp (λ)) of radiation, in addition to producing colored dissolved organic matter (CDOM) which also absorbs light (a g (λ)). a p (λ) is a representation of the absorption of all particulate materials, including algal (a ph (λ)) and non-algal particles (a NAP (λ)). In a simplified expression of the radiative transfer model, R rs (λ) can be represented in equation ( 1),
where a(λ) is the sum of a ph (λ), a NAP (λ), a g (λ) and a w (λ) (absorption by water molecules), and b b (λ) is the sum of b bp (λ) and b bw (λ) (backscatter by water). Therefore, R rs (λ) varies spectrally as a function of a ph (λ), a g (λ), a NAP (λ) and b bp (λ). The terms f and Q vary as a function of wavelength, solar zenith angle, aerosol optical thickness and surface roughness (Tzortziou et al., 2007). The ratio f/Q varies between 0.08 and 0.15 sr -1 (Morel and Mueller, 2002) and is assumed to be independent of wavelength and solar zenith angle in some algorithms (Tzortziou et al., 2007). Alexandrium spp. and similar large, peridinin-rich dinoflagellates have distinct R rs (λ) spectral signatures (Zapata et al., 2012). Their high intracellular pigment content leads to high phytoplankton absorption (a ph (λ)) throughout the visible spectrum, especially at blue and green wavelengths (Fig. S1 of the Supplementary Material). As a consequence, in Alexandrium blooms, R rs (λ) is low at visible bands (400-700 nm), except for the regions that comprise the red edge (670-700 nm) and near-infrared (NIR: 700-5000 nm). In the NIR region, R rs (λ) is enhanced by the cells' high chlorophyll fluorescence signal in combination with high b bp (λ), which contrasts their very low b bp (λ) at visible wavelengths (see Fig. 9.1 of Bernard et al., 2014).
Dinoflagellates, including A. catenella, exhibit high fluorescence compared to other phytoplankton groups, generating a marked red emission signal in the R rs (λ) spectrum when combined to their low b bp (λ) (Soto et al., 2016;Amin et al., 2009). Chlorophyll fluorescence can be estimated from R rs (λ) using the fluorescence line height (FLH, Behrenfeld et al., 2008), which is often calculated following equation ( 2):
] .
Eq. 2
Different phytoplankton groups such as cyanobacteria (Mishra et al., 2019), flagellates (Lavigne et al., 2022), diatoms, and dinoflagellates (Amin et al., 2009) also exhibit high FLH signals. For example, bright green blooms with high fluorescence signals and blooms of diatoms are often observed in the Bering/Chukchi region. Therefore, some color indexes can be used in combination in order to differentiate between A. catenella and other phytoplankton blooms.
Considering that immense diatom blooms are recurrent in the Bering and Chukchi seas, and that A. catenella blooms often succeed diatom blooms in other areas of the globe (Devred et al., 2018), we utilized the FLH peak difference as an important ocean color filter to differentiate between diatom and A. catenella blooms.
Marked differences in the spectral shape and magnitude of a ph (λ) are observed between dinoflagellates and diatoms (Fig. S1 of the Supplementary Material), with dinoflagellates absorbing a significant amount of radiation at the ultraviolet bands (300-400 nm) (Kahru et al., 2021;Orlova et al., 2022). However, this spectral range is not covered by Sentinel-3 and other legacy satellites. There are also differences in the a ph (λ) at green bands, between 500 and 550 nm (Fig. S1 of the Supplementary Material), but the effect of this a ph (λ) difference in the R rs (λ) spectra could not be observed using Sentinel-3's OLCI when comparing the 2022 Alexandrium catenella bloom and the 2024 dark diatom bloom (Fig. 2).
An interesting feature in the chlorophyll fluorescence emission signal in the R rs (λ) spectra, is the high gaussian R rs (λ) signal at the red edge to near infrared (NIR) region of the spectrum. When considering the OLCI
Factors of equations to calculate the abundance of A. catenella in cells/L from R rs (λ). Models were parameterized with complete datasets from 2022, considering sameday Sentinel-3 match-ups (mean R rs (λ) of all pixels within a 1 km radius from the sampling location).
Model Equation Abundance threshold a b c d e M1 Eq. 10 >10,000 cells/L 4.133 0.017 M2 Eq. 12 >3000 cells/L 4.009 5350 2.27E+06 PCR1 Eq. 15 >1000 cells/L 3.998 -0.304 2.712 -2.718 -4.016 PCR3 Eq. 16 >3000 cells/L -0.767 -1.574 -1.573 3.027 PCR10 Eq. 17 >10,000 cells/L 4.549 -1.055 -1.329 -3.777 8.770 Table 2 Uncertainty of the models to calculate the abundance of A. catenella (cells/L) from R rs (λ). Models were parameterized and validated using the complete dataset (selfvalidation), considering same-day Sentinel-3 match-ups, considering the mean R rs (λ) of all pixels within a 1 km radius from the sampling location.
Model Equation Abundance threshold R 2 RMSE MAE bias n M1 (DINI) Eq. 10 >10,000 cells/L 0.53 0.19 0.15 3.7 × 10 -6 23 M2 (EBI) Eq. 12 >3000 cells/L 0.67 0.25 0.20 8.1 × 10 -5 31 PCR1 Eq. 15 >1000 cells/L 0.68 0.37 0.32 1.1 × 10 -6 45 PCR3 Eq. 16 >3000 cells/L 0.67 0.25 0.20 7.8 × 10 -4 31 PCR10 Eq. 17 >10,000 cells/L 0.71 0.15 0.12 -4.4 × 10 -5 23 KBBI Eq. 6 >3000 cells/L 0.18 4.36 4.34 -4.3 31 KBBI Eq. 6 >10,000 cells/L 0.24 4.55 4.54 -4.5 23 NDNI Eq. 13 >3000 cells/L 0.19 5.00 4.99 -5.0 31 NDNI Eq. 13 >10,000 cells/L 0.12 5.18 5.18 -5.2 23 ( Tao et al., 2013).
Here, we used this spectral difference in R rs (λ) to differentiate A. catenella and diatom blooms. We calculated the fluorescence line height at 682 and 709 nm using Equation (3),
where λ is the wavelength where the peak of the fluorescence line height lies, which can be attributed to 682 or 709 nm considering the Sentinel-3 bands. Then, blooms were flagged as A. catenella if FLH 709 > FLH 682 . Finally, we used a quantitative model (PCR, NDNI, EBI, etc.) to calculate the abundance of A. catenella cells in the valid pixels. This masking method effectively differentiated blooms of A. catenella from diatom blooms, as shown in Fig. 2c-h. To differentiate A. catenella from algae that form green blooms (e.g., green algae, cyanobacteria), we used two alternative color indexes: the green line height (GLH) and the blue line height (BLH). First, GLH was used to differentiate between A. catenella blooms and bright green blooms. The GLH is expressed in equation ( 4 ] .
Eq. 4
Based on these color indexes, green blooms can be distinguished from A. catenella blooms using the green bloom index (GBI), expressed in equation ( 5 GBI = FLH*GLH*10 6 , Eq. 5
where green blooms are considered dominant where GBI >1, as shown in Fig. 3.
2.4.2.1. Color index models. Because of their high chlorophyll fluorescence emission compared to other phytoplankton groups (Fig. 2a and b), satellite detection of dinoflagellates can be performed using the Karenia brevis bloom index (KBBI, calculated as the normalized index between the bands 682 and 665 nm, equation ( 6) following Amin et al., 2009), where
Eq. 6
However, high backscatter levels from suspended sediments and non-algal particles, sun glint and backscatter of atmospheric particles can significantly increase the R rs (λ) at red and near-infrared wavebands, increasing the KBBI signal. To avoid this contamination of the R rs (λ) signal, we opted to use the fluorescence line height (FLH), which is minimally influenced by the high absorption by water (thus low R rs (λ)) in the NIR.
As the concentration of dinoflagellate cells increases during a bloom event, the R rs (λ) peak at the red edge of the visible spectrum (generated by chlorophyll fluorescence emission) is enhanced and moves towards the NIR portion of the spectrum due to the increase in backscatter by phytoplankton cells (IOCCG, 2021), especially when these cells are concentrated in the near surface. As a consequence, there is a spectral shift in the peak of the FLH (682 nm) towards the NIR (710 nm). To minimize this spectral shift effect, we replaced FLH with FLH max , calculated in equation ( 7),
] ,
Eq. 7
where FLH max is the fluorescence line height at the wavelength of maximum emission peak between R rs (682) and R rs (709), R rs (λ max ) is the highest value between R rs (682) and R rs (709), and λ max is wavelength where R rs (λ) is the highest, considering 682 or 709 nm. FLH max expresses an emitted spectral signal that often travels through thin clouds, improving bloom diagnosis when compared to other band ratio models. However, it can be significantly enhanced in certain cloud conditions (shadow, thin clouds) (Fig. S3 in Supplementary Material).
The BLH (adapted from Tomlinson et al., 2009), is a proxy for chlorophyll-a absorption and less subject to changes in CDOM and NAP when compared to the classic blue/green ratio (R rs (490)/R rs ( 560)). The BLH is calculated using equation ( 8),
] .
Eq. 8
Here, for line height reference, we replaced R rs (560) (R rs (555) for MODIS as used in Tomlinson et al., 2009) with R rs (665) because the latter is less subject to sun glint and atmospheric and ocean backscatter effects that enhance gain.
The color indexes FLH max , GLH and BLH extracted from Sentinel-3 OLCI R rs (λ) observations showed significant correlation (p > 0.05) with A. catenella cell abundances, despite the relatively low uncertainty metrics, especially for FLH max (Fig. S6 and Table S1 in Supplementary Material). These color indexes were combined in two models. The first model (M1) is based on the dinoflagellate bloom index (DINI) which assumes that A. catenella blooms present high FLH max and low GLH (compared to other blooms), and that R rs ( 490) is reduced due to the radiation absorption by chlorophyll. The DINI is calculated using equation ( 9),
FLH max GLH* R rs (490) , Eq. 9
and model 1 (M1) calculates the abundance of A. catenella using equation (10),
Eq. 10
The second model (M2) is based on the enhanced bloom index (EBI) which detects dinoflagellate blooms in developing to senescent state, where cells are expected to be under physiological stress, with high chlorophyll-specific fluorescence (Tomlinson et al., 2009;Alexanin et al., 2023). The EBI is calculated using equation ( 11),
Eq. 11
In M2, A. catenella cell abundances are calculated following equation ( 12):
Eq. 12
In equations ( 10) and ( 12), x is the abundance of A. catenella in cells/L, and parameters (a, b and c) of both equations are listed in Table 1.
In addition to the above-mentioned color indexes (FLH max , GLH, BLH, and KBBI), we also evaluated empirical relationships between A. catenella cell abundances and the normalized difference Noctiluca model (NDNI), a method that has been used to estimate cell abundances of Alexandrium tamarense in the Galician Rías (Spain) by Sacilotto Detoni et al. ( 2023). The NDNI is calculated using equation ( 13),
Eq. 13
Next, we used a Principal Component Analysis (PCA) to decompose the R rs (λ) spectra into dominant modes (loadings) with specific spectral shapes. These spectral modes show specific characteristics that represent the contribution of optical properties of distinct water constituents (absorption by different pigments and CDOM as well as particulate backscatter). The sum of the scores (magnitudes) of the spectral modes result in the observed reflectance spectra (Bracher et al., 2015;Lange et al., 2020).
Considering that specific phytoplankton taxa contain unique combinations of pigments and backscatter properties, the scores of some specific spectral modes can be used as predictors of the abundance or biomass of target phytoplankton species. Following equation ( 14), the cell concentration of a target phytoplankton group (x, cells/L) is a function of the scores of a R rs (λ) spectrum (u {1, 2, …, p} ) in the pth PC from the PCA, the regression coefficients (b {1, 2, …, p} ) and the intercept a:
Eq. 14
The principal component regression models were parameterized using subdatasets which were extracted from the complete dataset based on three thresholds for A. catenella cell abundances: >1,000, >3000 and > 10,000 cells/L. This was done to allow the comparison of PCR models with the color index models, which were also parameterized using these subdatasets. For the color index models, the thresholds were established based on a sensitivity analysis where, for each model (including the color index models), the dataset was subsampled using various thresholds (ranging from 10 3 to 10 4 cells/L) and the subdataset with the best model fitting (highest R 2 , at p < 0.05) was selected for each model.
For each of these subdatasets, different principal components (PCs) were selected as explanatory variables using a stepwise selection based on the Akaike Information Criteria (stepAIC function in R). At cell abundances above 1000 cells/L, A. catenella is estimated using equation (15) (PCR1),
Eq. 15
where PC2, PC4, PC5 and PC6 are the scores of that R rs (λ) spectrum at principal components 2, 4, 5 and 6 of a principal component analysis performed with the reference R rs (λ) spectra library. Similarly, principal component regression models parameterized with the dataset restricted to cell abundances above 3000 and 10,000 cells/L are described in Equation ( 16) (PCR3) and 17 (PCR10), Eq.17
The parameters for equations ( 15)-( 17) are listed in Table 1, whereas the spectral distribution of the loadings (PCs) used in equations ( 15)-( 17) are shown in Fig. S4 of the Supplementary Material.
The spectral distribution of PC loadings resemble spectral features of optical properties of the optically-active constituents in the ocean (Craig et al., 2012;Bracher et al., 2015;Lange et al., 2020). The spectral variability of the PC1 is driven mainly by the particulate backscattering of the in-water constituents and the absorption of water molecules, and explained more than 96 % of the data covariance R rs (λ) dataset with A. catenella cell abundances >3000 cells/L (PCR3). The second PC (PC2), used in models with cell abundances >1000 and > 10,000 cells/L (PCR1 and PCR10) highlights spectral features related to the absorption by chlorophyll at the ocean surface, explaining ~6.9 and 2.3 % of the covariance of the PCR1 and PCR10 R rs (λ) datasets. The third and fourth PCs (PC3 and PC4) are driven by the spectral variation of R rs (λ) due to the absorption of accessory pigments and explained less than 0.3 % of the dataset covariance.
All model performances were evaluated by comparing in situ (observed) and modeled (predicted) A. catenella cell abundances. From this comparison, performance metrics such as the coefficient of determination (R 2 ), root mean square error (RMSE), mean absolute error (MAE), and bias were calculated following equations described in Lange et al. (2020). These metrics were calculated in log 10 space, thus MAE and RMSE values can be interpreted as relative or percentage uncertainties (Seegers et al., 2018). Moreover, the comparisons between different models were based on the Akaike Information Criterion (AIC), where the model with the best performance presents the lowest AIC value.
Considering the reduced number of observations (Table 2) for model parameterization and validation, model uncertainties were calculated in two dataset arrangements: 1) Full-fit model parametrization and selfvalidation, where models were trained and validated using the complete dataset, and; 2) Partial model parametrization with crossvalidation, where models were trained with randomly sub-sampled training datasets (80 % of the complete dataset), then used to compute cell abundances using the remaining 20 % of the dataset, and finally these predictions were compared with observations from this 20 % sub-dataset. This process was repeated (2000 Monte-Carlo permutations) and the average performance metrics were computed. These arrangements (1 and 2) were used to evaluate model performance and robustness, respectively. Similar results of the full-fit and the crossvalidation uncertainties indicate that the model is not overtrained to this regional dataset. Finally, the uncertainty metrics of the different models were compared in order to choose the best model to estimate cell abundances of A. catenella in the Bering and Chukchi seas.
The statistical analyses, model parameterizations, graphic figures and maps of model implementations to satellite data were computed in R (R Core Team, 2017) using the packages stats (Venables and Ripley, 2002), scales, and ncdf4 (Pierce, 2015). Figures of model output and true color maps were generated using the software SeaDAS (NASA).
In the color index models, high model outputs are generated when FLH max is high and GLH is low (Equation ( 7)). However, high outputs are also generated when FLH max is low and GLH is even lower (which is the case of blue waters with low phytoplankton concentrations). To avoid clear water pixels (where a phytoplankton bloom event is not present) and cloud cover, image masks were applied prior to the implementation of all A. catenella models. In addition, bloom identification masks were used to differentiate dinoflagellates from diatoms or green algae. In summary, the following criteria were used to discard (mask) pixels with bad R rs (λ) spectra: 1) FLH < 0.0001 (not a bloom); 2) GLH < 0.0005 (insufficient phytoplankton concentration); 3) R rs (665) < R rs (754) (high NIR, not water); 4) R rs (865) > 0.005 (cloud mask); 5) R rs (865) < 0 (mask for negative R rs (λ) spectra due to cloud shade or inefficient atmospheric correction); 6) GBI > 1 (green bloom); 7) FLH 709 -FLH 682 < 0 (diatom bloom).
In summary, to implement the algorithms (color index models or PCR models) to quantify A. catenella cells in surface waters of the Bering and Chukchi seas, one must do the following: A) download the Sentinel-3 image of interest; B) process the image for atmospheric correction using ACOLITE; C) Use the image masks to discard pixels with clouds and other image artifacts; D) Use image masks ( 6) and ( 7) to differentiate between A. catenella and other phytoplankton blooms; E) Run the A. catenella cell quantification algorithm of choice using R rs (λ) data as input, and finally; F) project the estimated A. catenella cell abundances in maps. Using the image masks suggested above is key to avoid wrong taxonomic identification of the bloom. The workflow of actions for algorithm parametrization, validation and implementation is shown in Fig. 4.
Despite the reduced number of match-ups between samples and satellite data (45 match-ups from 2680 shipboard measurements), model self-correlations (Table 2) demonstrated similar uncertainties to model cross-correlations (Table 3). Regarding the color index models, M1 showed its best performance when A. catenella cells were present at concentrations above 10,000 cells/L, whereas M2 was efficient for the detection of A. catenella with cell concentrations above 3000 cells/L (Table 2).
The PCR model showed good performance for predicting A. catenella cell abundances from PC scores (Table 2, R 2 > 0.5, Fig. 5). In models parameterized with cell abundances higher than 10,000 cells/L, the performance of the PCR (PCR10) model (self-correlation R 2 = 0.71, cross-correlation R 2 = 0.63) was superior to the color index models (R 2 = 0.53 and 0.67 for M1 and M2, respectively) (Fig. 5, Table 2, Table 3). However, for cell abundances higher than 3000 cells/L, the PCR (PCR3) and M2 models demonstrated similar performances (R 2 = 0.67), Fig. 5, Table 2, Table 3). For abundances below this threshold (and above 1000 cells/L), only the PCR (PCR1) model generated satisfactory outputs (R 2 = 0.68, Fig. 5c-Table 2). For reference, the comparison between measured and modeled A. catenella cell abundances using model PCR1, considering the full dataset which includes samples with cell abundances lower than 1000 cells/L, is shown in Fig. S3d of the Supplementary Material.
The principal component regression (PCR) model to estimate A. catenella abundance and color index models to detect Alexandrium blooms were implemented in Sentinel-3 images of blooms with similar spectral signatures, but which occurred in early Spring, when A. catenella blooms are not currently known to occur in this region (high likelihood of false positives of A. catenella bloom detection).
Model projections using Sentinel-3 images of the Bering and Chukchi seas on three bloom dates (Fig. 6) show that, among the color index models, M1 is able to distinguish between A. catenella blooms and green blooms with high FLH max (Figs. 6e and 3), despite being efficient only for cell concentrations above 10,000 cells/L. Conversely, M2 detects A. catenella and other blooms (green bloom -Fig. 6f, and possible diatom blooms, Fig. 6i). As a final step, the PCR1 model was implemented in Sentinel-3 data shown in Fig. 6d and g to test its ability in distinguishing between A. catenella blooms and other bloom events.
The color index algorithms presented here effectively capture the distinct spectral characteristics of Alexandrium-like dinoflagellates, allowing for rapid identification, quantification and mapping of potential HAB events. Taxonomic confirmation that blooms are in fact A. catenella and not other peridinin-rich dinoflagellates can be made more efficient by using derived maps to direct ship-based sampling and other in situ bloom measurements in a regionally integrated observing system to protect human health from paralytic shellfish toxin poisonings. Once confirmed, A. catenella blooms can be more confidently mapped and cell abundances quantified over vast temporal and spatial windows using the satellite algorithms presented in this manuscript. After color masks are implemented to differentiate between dinoflagellates and other phytoplankton groups, model M1 (based on index DINI (FLH and GLH)) effectively quantified A. catenella cells in the dinoflagellate bloom when cell concentrations exceed 10,000 cells/L, while model M2 (based on EBI (FLH and BLH)) was effective for concentrations higher than 3000 cells/L. In turn, the principal component regression (PCR) model, which utilizes all available bands provided by Sentinel-3, quantifies A. catenella cells at abundances as low as 1000 cells/L.
A key to the success of these algorithms is the use of ocean color filters for prior differentiation between bloom populations dominated by dinoflagellates and those dominated by other phytoplankton group that in some cases have similar optical characteristics. Examples of other potentially confused groups include ice-algae and benthic/tychoplanktonic blooms which occur near the ice edge and in shallow turbulent waters. Because they are adapted to low light levels, they thus emit high chlorophyll fluorescence when exposed to high light conditions at the sea surface, increasing their FLH max . In addition, benthic (or tychoplanktonic) diatoms often form agglomerates, increasing their apparent particle size to a range where geometric scattering prevails, reducing b bp (λ) (low particulate backscatterb bp (λ) hence low remote sensing reflectance -R rs at 560 nm (R rs ( 560)). This combination of high FLH max and low b bp (λ) is similar to the optical characteristics of a dinoflagellate bloom observed by multispectral sensors such as Sentinel-3's OLCI. Another example is that some blooms of flagellates and cyanobacteria, which are bright green due to high backscatter, present similarly high FLH values and bloom indexes as A. catenella blooms (see Fig. 6d-f), in which case other color indexes (e.g. GBI) must be used for differentiation. Most importantly, the difference between FLH709 and FLH682, which we use as a filter, is essential for differentiating between diatoms and dinoflagellates (Tao et al., 2013). Additionally, bloom seasonality and timing are important for interpreting maps of bloom distribution generated by these models, since various phytoplankton groups have different optimal growing seasons.
Hereafter we examine the model shortcomings, suggested Fig. 4. Workflow of actions to parameterize and implement models to estimate cell abundances of Alexandrium catenella from R rs (λ) provided by Sentinel-3 OLCI: IFCB -Imaging FlowCytobot; ACOLITE is a method for atmospheric correction (Vanhellemont and Ruddick, 2021).
Uncertainty metrics of the models to calculate the abundance of A. catenella (cells/L) from R rs (λ), with models parameterized and validated using cross-correlation (80 % for parameterization, 20 % for validation, 2000 bootstrap) considering same-day Sentinel-3 match-ups within a 1 km radius from the sampling location.
Model Equation Abundance threshold R 2 RMSE MAE bias M1 (DINI) Eq. 10 >10,000 cells/L 0.54 ± 0.33 0.20 ± 0.05 0.163 ± 0.05 0.3 ± 10.3 × 10 -2 M2 (EBI)
Eq. 12 >3000 cells/L 0.37 ± 0.24 0.76 ± 0.16 0.652 ± 0.15 6.5 ± 1.5 × 10 -2 PCR1 Eq. 15 >1000 cells/L 0.48 ± 0.25 0.53 ± 0.11 0.440 ± 0.10 0.8 ± 2.0 × 10 -2 PCR3 Eq. 16 >3000 cells/L 0.53 ± 0.25 0.29 ± 0.06 0.248 ± 0.06 0.5 ± 1.2 × 10 -2 PCR10 Eq. 17 >10,000 cells/L 0.63 ± 0.28 0.24 ± 0.11 0.204 ± 0.09 3.3 ± 1.2 × 10 -2 applications, and ecological context for HAB event monitoring. For bloom monitoring purposes, the fact that Alexandrium bloom events can still present considerable toxin risk at cell abundances lower than 1000 cells/L, which is lower than those detected by our satellite models (Vandersea et al., 2018;Lefebvre et al., 2016), must be taken into consideration. Moreover, based on mapped model outputs (see Fig. show results for models parameterized with the dataset where cell abundances are higher than 3000 cells/L, and; in panels g,h,i) models were parameterized for cell abundances higher than 10,000 cells/L. Validation and cross-validation uncertainty metrics are listed in Tables 2 and 3, respectively. In panels a-c, note that the Alexandrium patch was detected east of the Diomede Islands when sampled in situ (panel a). By comparison, the satellite imagery predicts higher bloom densities in the west, outside the ship's survey area. The period of the survey (July 22 to 31) preceded the satellite image acquisition (August 2, panel b). Given the difference in observation times, this satellite image was not used in the match-up dataset, however, it is shown here for illustration purposes since this was the image where the bloom could best be observed without cloud cover obstruction.
observed in other open ocean areas (McGillicuddy et al., 2014) but not in the North Bering or Chukchi seas. Thus, at extreme cell abundances, bloom presence can be robustly detected, but there is considerable uncertainty in their quantification of absolute cell densities.
The PCR model provides more accurate quantitative A. catenella abundances at lower cell concentrations, which is valuable for understanding bloom evolution and dynamics. Previous PCR models demonstrated good performance for retrieving group-specific phytoplankton biomass of cyanobacteria (Craig et al., 2012), cell abundances of picophytoplankton cyanobacteria and eukaryotes (Lange et al., 2020) and other phytoplankton groups (Bracher et al., 2015). PCR models are expected to improve significantly when parameterized and implemented using hyperspectral R rs (λ) imagery (Lange et al., 2020) because unique optical pigment absorption signals occur at lower spectral resolution than that provided by legacy multispectral satellites (IOCCG et al., 2021). In the future, hyperspectral models can be implemented using data from new satellites such as the Plankton, Aerosol, Clouds, ocean Ecosystem (PACE) (Cetinić et al., 2024). These models are expected to improve the detection of dinoflagellate blooms, increasing taxonomic discrimination, but will still benefit from application of the spectral masks presented in this manuscript as a way streamline their computation.
Color index models similar to those presented here (DINI and EBI) have been shown to detect blooms of other dinoflagellates such as Karenia spp., which present some similar optical spectral features as A. catenella such as the high chlorophyll fluorescence emission at the red edge and low particulate backscatter (Amin et al., 2009;Soto et al., 2016). The DINI (M1) and EBI (M2) indexes, as well as the PCR model, were more efficient for detecting and quantifying A. catenella blooms than the KBBI (Amin et al., 2009) and NDNI (Soto et al., 2016) (Figs. S6 and S7), which were developed to detect Alexandrium blooms in other areas of the globe. Despite the lower accuracy (Table 2, Fig. 6 and A. catenella (Figs. S6 and S7 in Supplementary Material). However, KBBI, NDNI and EBI (M2) also identify blooms of other highly fluorescent phytoplankton such as the green bloom (possibly Phaeocystis), that co-occurred with an Alexandrium-like bloom in the Bering Strait on July 6th, 2023 . In this bloom event, high FLH is observed in both blooms (Fig. 3d), and the GLH (or GBI) was required to differentiate between the possible A. catenella signal and the green bloom (Fig. 3b and c).
Most color indexes (Table 4) are influenced by atmospheric components such as aerosols, clouds, fog, cloud shadows and sun glint (see Fig. S3 in the Supplementary Material). For example, KBBI (normalized difference between 682 and 665 nm) is high in spectra where there is an exponential increase in R rs (λ) from the visible towards the near-infrared (NIR), which is often observed in fog. However, the same high KBBI value may be observed in a R rs (λ) spectrum with low NIR and a peak at 682 nm. To address this issue, we applied image masks prior to model implementation. These masks use portions of the R rs (λ) spectrum (e.g. NIR) where the water-leaving R rs (λ) signal is significantly lower than that from cloud tops, land, and floating material (Huber et al., 2023). Masks also consider spectral features that are expected in phytoplankton blooms such as some level of blue R rs (λ) decrease due to the electromagnetic energy absorption by chlorophyll (expressed in the GLH index), and a minimum chlorophyll fluorescence emission (FLH) by the bloom (IOCCG, 2018). Therefore, if masks are not implemented, false bloom positives may be detected when clouds and fog are present.
The models presented here are parameterized and validated with satellite R rs (λ) data matched to in situ biological measurements, serving as proof-of-concept since regional inversion models must be parameterized and validated using in situ radiometric (R rs (λ)) and biological (phytoplankton biomass) measurements. The use of these data avoids the influence of the uncertainty embedded in satellite R rs (λ) retrievals and ensures similar spatial scales and minimal temporal mismatch between radiometric and biological measurements (IOCCG, 2020). Therefore, further reparameterization with in situ data is a priority for continued efforts. However, these models are still useful for fast and remote detection of HABs, especially in regions such as the Bering and Chukchi Seas where combined in situ phytoplankton and radiometric measurements are extremely scarce.
The limited number of match-ups between satellite and ship-based observations for model parameterization and validation is another constraint of the models presented in this study: only 218 (or 45 after masks are applied) of the 2680 in-situ measurements (8 %) had Sentinel-3 match-ups. This constraint is expected in the Pacific Arctic, where cloud and ice cover are prevalent most days of the year and sampling is restricted primarily to the summer ice free months. The reliance on satellite R rs (λ) data for model parameterizations significantly restricted the number of observations considered in these models, highlighting the importance of collecting in situ radiometric measurements concurrent to the phytoplankton samples.
Satellite observations of phytoplankton bloom events are a powerful monitoring tool but are not sufficient by themselves for effective A. catenella monitoring. Their limitations include the fact that these products poorly resolve inshore areas that may harbor localized A. catenella blooms (Anderson et al., 1994;Brosnahan et al., 2020). Satellite observations are also restricted to only the first few metres from the surface of the water column in the open water areas for which they are most adept (IOCCG, 2021). Given that A. catenella blooms originate from cyst beds located on the ocean bottom, they may cover vast subsurface areas before reaching the optically-observable ocean surface, especially in turbid waters such as the Alaskan Coastal Current (see Fig. 1). Therefore, to be effective, A. catenella bloom monitoring must rely on a mixture of approaches and practices that detect A. catenella cells and toxins with sufficient sensitivity to protect human health. Satellite observations augment these systems by alerting communities to the presence of large offshore blooms that pose some of the greatest paralytic shellfish toxin poisoning risks. If blooms are far from shore, satellite detection can also provide substantial advance warning, making it possible for communities to enhance and/or better focus their nearshore in situ observing efforts.
Overall, our analysis suggest that the developed models are a valuable tool for 1) locating dinoflagellate blooms and targeting in situ sampling by research vessels or coastal communities in the vicinity of the bloom, which can provide more information about bloom taxonomic composition and toxicity, and; 2) quantifying the concentration of A. catenella cells and diagnose the extent and transport of blooms that are initially detected (and confirmed) in situ during monitoring or scientific research efforts. As an example, during Summer 2023, fauna intoxication events in the Gulf of Alaska were reported by the Knik Tribe (Martínezcuello, 2023), and another A. catenella bloom event was detected during a cruise of opportunity passing through the Bering Strait region. In response to these findings, dinoflagellate distribution maps Table 4 Advantages and limitations of the predictive dinoflagellate models presented in this manuscript: Principal Component Regression (PCR) and the color indexes dinoflagellate index -DINI (which uses the fluorescence line height -FLH and the green line height -GLH and is the basis of model M1), enhanced bloom index -EBI (combination of FLH with blue line height -BLH, which is the basis of M2), Karenia brevis bloom index -KBBI from Amin et al. (2009) and Normalized difference Noctiluca index (NDNI) from Sacilotto Detoni et al. (2023).
Advantages Limitations
• Utilizes all available R rs bands, optimizing phytoplankton taxon-specific quantification using ocean color information; • Differentiates bloom signal from clouds, fog and sun glint.
• Possibly overestimates cell abundances higher than 200,000 cells/L a ; • Computationally demanding for generating maps.
• The reduced number of wavebands (6) facilitates model implementation using satellite imagery. • Color signal is visible through fog, thin clouds and mild sun glint; • Distinguishes dinoflagellate blooms from blooms of green algae and cyanobacteria.
• Possibly overestimates cell abundances higher than 200,000 cells/L a ;
• High uncertainty in the quantification of cells when abundances are lower than 10,000 cells/L (low index signal is similar to that of atmospheric artifacts and suspended sediments in the ocean).
• The reduced number of wavebands (6) facilitates model implementation using satellite imagery; • Color signal is visible through fog, thin clouds and mild sun glint.
• Possibly overestimates cell abundances higher than 200,000 cells/L a ;
• Unable to distinguish between dinoflagellates and other bloom types such as green algae and cyanobacteria.
• The reduced number of wavebands (2) facilitates model implementation using satellite imagery.
• Overestimates cell abundances;
• Defines clouds and cloud shadows as phytoplankton;
• Highly affected by the presence of thin clouds, fog and sun glint. NDNI
• The reduced number of wavebands (2) facilitates model implementation using satellite imagery.
• Overestimates cell abundances;
• Defines clouds and cloud shadows as phytoplankton;
• Highly affected by the presence of thin clouds, fog and sun glint.
a Model uncertainty cannot be quantified for cell abundances above 200,000 cells/L due to in-situ (ground-truth) data limitation.
P.K. Lange et al. Journal of Environmental Management 380 (2025) 125042
(Fig. S10 of the Supplementary Material) were produced using the index DINI as an indicator of A. catenella. These maps provided information about potential bloom extent, location and transport, and were included in reports transmitted to the scientific community and Alaskan stakeholders.
This study demonstrates the efficacy of satellite remote sensing as a valuable tool for detection and monitoring of HABs in the Pacific Arctic region. The regionally-parameterized HAB models are able to rapidly detect A. catenella blooms in the Bering and Chukchi seas. The developed HAB algorithms can be a helpful tool aiding northern Bering and Chukchi Sea coastal communities that rely on the comprehensive use of the marine ecosystem as food, with data for early warning systems. Given the limited coverage of sampling across this dynamic and transboundary region, satellite bloom detection offers a promising method for augmenting in-situ observations or directing in-situ efforts. Moreover, remote observations conducted by satellites are capable of detecting blooms in regions where it is not possible to collect in situ data (ex. Russian Federation). These observations are essential to identify regions where blooms frequently originate (which can be associated with the occurrence of cyst beds), aiding our understanding about bloom origin, transport, fate, and phenology.
The combination of a FLH diff maks with the PCR model provided efficient quantification of A. catenella cells and discrimination from blooms of diatoms and green blooms. We advise the combined use of masking layers, using GLH (present work) and FLH diff to differentiate A. catenella blooms from other recurrent blooms in the Bering/Chukchi region: GLH separates A. catenella from blooms of green algae, whereas FLH diff segregates A. catenella from diatom blooms.
The model frameworks presented here can also be parameterized to detect dinoflagellate blooms, especially those dominated by species of Alexandrium and potentially Karenia, in temperate, subtropical, and tropical regions. Moreover, the A. catenella algorithms can be adapted to other satellite sensors such as MODIS and VIIRS to assess longer timeseries of past images, contributing significantly to understanding, monitoring, and managing the increasing threat of HABs in the Pacific Arctic. By providing timely and accurate information on bloom distribution, these techniques can aid in mitigating the risks associated with HABs in the Pacific Arctic. These satellite-derived products can also guide targeted testing on the ground, improving the effectiveness and efficiency of the larger HAB monitoring and surveillance system that is needed to protect all peoples in northern Bering Sea and Chukchi Sea coastal communities, to whom HAB events are a novel and emerging threat to food security, food safety, human health, as well as a potential conservation concern for diverse marine taxa.
Journalof Environmental Management 380 (2025) 125042