Seasonal dynamics of major phytoplankton functional types in the coastal waters of the west coast of Canada derived from OLCI Sentinel 3A

1 Introduction

Phytoplankton represent only 0.2% of global primary producer biomass; however, they are responsible for more than 90% of the ocean’s primary productivity (Chassot et al., 2010), including half of the biosphere’s net primary productivity (Field et al., 1998). These ubiquitous single-celled autotrophs are the base of the marine food web (Falkowski, 2012), and are capable of fixing ~50-130 million tons of atmospheric carbon dioxide (CO₂) daily into the form of organic materials, playing a central role in the global biogeochemical cycles (Behrenfeld et al., 2006). Numerous in situ methods are available to quantify size fractions and community composition of phytoplankton in the ocean, yet these observations are time-consuming, labor-intensive, and sparse in space and time (IOCCG, 2014; Chase et al., 2017).

Space-based optical sensors measure water-leaving radiance in multiple wavebands, allowing the retrieval of the surface total chlorophyll-a (hereafter TChla), a proxy for phytoplankton biomass, with high spatial and temporal resolutions (McClain, 2009). Nevertheless, TChla alone is not sufficient to describe the ecological and biogeochemical processes in the ocean (Bracher et al., 2017). Extracting information about phytoplankton functional types (PFTs), and phytoplankton size classes (PSCs), is essential for understanding the linkage between phytoplankton composition and fisheries, and in determining the role of the ocean in a changing climate (e.g., Falkowski et al., 1998; Le Quere et al., 2005; IOCCG, 2014). Considering the highly diverse nature (e.g., taxonomical, morphological, and physiological) of marine phytoplankton, climatologists and marine biogeochemical modelers have defined PFTs according to their specific function or role in the biogeochemical cycles irrespective of their taxonomical affiliations (Le Quere et al., 2005; Nair et al., 2008). For example, diatoms are considered to be essential silicifiers, cyanobacteria as nitrogen fixers, Haptophytes as calcifiers, and Dinoflagellates/Haptophytes/Green Algae as dimethyl sulphide producers (DMSP) (IOCCG, 2014).

Given the importance of determining phytoplankton from space, several satellite-based algorithms have been developed to identify the significant PFTs (e.g., Alvain et al., 2005; Navarro et al., 2014; Correa-Ramirez et al., 2018), derive the chlorophyll-a concentration (hereafter Chla) of PFTs (e.g., Bracher et al., 2009; Sadeghi et al., 2012; Moisan et al., 2017), retrieve PSCs (e.g., Ciotti et al., 2002; Kostadinov et al., 2009; Mouw and Yoder, 2010; Brewin et al., 2010; Devred et al., 2011; Hirata et al., 2011; Brotas et al., 2013; Roy et al., 2013; Lamont et al., 2018; Moore and Brown, 2020; Liu et al., 2021), and retrieve accessory pigment concentrations (e.g., Pan et al., 2010; Bracher et al., 2015; Wang et al., 2018; El Hourany et al., 2019; Stock and Subramaniam, 2020; Kramer et al., 2022). These satellite-based algorithms explore specific signatures, such as anomalies in the reflectance spectra, inherent optical properties (IOPs), including TChla concentration (see IOCCG, 2014; Mouw et al., 2017; Bracher et al., 2017 for a detailed review). Of the many approaches developed, the Empirical Orthogonal Function (EOF) approach has been successfully used among these algorithms (e.g., Craig et al., 2012). The EOF approach is a statistical technique that reduces the high dimensionality of spectral data into dominant modes, explaining the variance of structures in the data. Specifically, the EOF approach allows the extraction of the spectral signature associated with specific PFTs, thus enabling the development of empirical models for various PFTs (e.g., Xi et al., 2020). This approach has been successful in retrieving the concentration of different accessory pigments from Medium Resolution Imaging Spectrometer (MERIS) and hyperspectral data (Bracher et al., 2015), Chla concentration of multiple PFTs using GlobColor merged and OLCI data (Xi et al., 2020), and cell abundance of Prochlorococcus, Synechococcus, and autotrophic picoeukaryotes based on multispectral and hyperspectral data (Lange et al., 2020). However, most of these studies were performed in open oceans or optically Case-1 type waters, with less application to coastal systems (e.g., Craig et al., 2012).

The coastal ocean of British Columbia, Canada, specifically the Strait of Georgia, is a spatially and temporally dynamic Case-2 water body (Phillips and Costa, 2017) crucial for Pacific salmon, herring, and hake fisheries (Beamish et al., 1994; Beamish et al., 2012; Thomson et al., 2012), which are of vital importance to the regional economy and First Nations (Nesbitt and Moore, 2016). In these waters, phytoplankton have been shown to be associated with the growth and survival of these species (Cushing, 1990) through their bottom-up control of key zooplankton species such as Neocalanus plumchrus (e.g., Sastri and Dower, 2009; Perry et al., 2021; Suchy et al., 2022). Hence, spatial-temporal satellite-derived PFTs will allow for a more coherent understanding of the function/role of multiple PFTs on the regional food web.

Here, we adapted an EOF-based algorithm developed by Xi et al. (2020) to the Strait of Georgia, British Columbia (BC), Canada, to retrieve TChla concentration and the Chla concentration of multiple PFTs. The model was developed for each PFT using the matchup between CHEMTAX-derived PFTs and EOF scores derived from the OLCI remote sensing reflectance (R_rs; sr^–1(λ)). The robustness of the fitted regression model was assessed statistically by the cross-validation procedure, and independent validation of the major derived PFTs was performed against an independent in situ CHEMTAX time-series. Further, the seasonal and spatial dynamics of major PFTs were determined by applying the regression model to daily OLCI imagery.

2 Materials and methods

2.1 Study area

The Strait of Georgia (SoG) is the largest semi-enclosed estuarine system in the Salish Sea, located on the west coast of Canada (Figure 1). This partially enclosed region has a total length of 222 km, a maximum width of 28 km, and a depth of about 350 m within its central region (Masson and Peña, 2009), which is generally divided into three marine areas: the southern, central, and northern Strait. The SoG is connected to the Pacific Ocean through the Johnstone Strait in the northern region and Juan de Fuca Strait in the southern boundary (Peña et al., 2016). An estuarine circulation typically drives the SoG regional dynamics, characterized by low saline surface waters from the adjacent Fraser River and highly saline nutrient-rich subsurface waters from the Pacific Ocean (Li et al., 2000; Masson and Cummins, 2004). The largest freshwater discharge to the SoG is from the Fraser River, which contributes 75% of the freshwater to the Strait, and is primarily responsible for the water column stratification and estuarine-like circulation in the central and southern Strait (Harrison et al., 1983; Johannessen et al., 2003; Pawlowicz et al., 2007). Incoming seawater dominantly enters the SoG through the southern Juan de Fuca Strait (Masson, 2002), bringing nutrient-rich water from the Pacific Ocean through the deep water inflow and is eventually upwelled into the surface waters by rigorous tidal mixing in the Haro Strait (Li et al., 2000). In the summer, the northern Strait becomes highly stratified, limiting surface nutrient renewal. This stratification is thought to be a function of high summer surface heating, low summer winds, and localized freshwater discharge, and can periodically be disrupted by strong wind events (Evans et al., 2019).

Figure 1 Map showing the study area. The red circles indicate the sampling locations in the central and northern Strait of Georgia (SoG). QU39 represents the Hakai time-series station. Bathymetry is shown in meters.

This region is one of the most productive estuaries in the Pacific North East with an average annual primary productivity of about 280 g C m^-2, which is the highest in the frontal regions and close to the Fraser River (Harrison et al., 1983; Sutton et al., 2013; Johannessen et al., 2021). Seasonally, SoG surface TChla greatly varies with winter concentrations <1 mg m^-3 and increases to > 15 mg m^-3 during spring (Harrison et al., 1983; Jackson et al., 2015; Peña et al., 2016; Suchy et al., 2019). As in many temperate coastal environments, spring blooms are dominated by diatoms with Skeletonema sp., Thalassiosira spp., and Chaetoceros spp. being the dominant bloom-forming species. By the end of spring, diatoms are usually replaced by high diversity communities constituted by flagellates such as cryptophytes, dinoflagellates, prasinophytes, and haptophytes (Harrison et al., 1983; Haigh and Taylor, 1991; Del Bel Belluz et al., 2021).

Optically, the SoG is a dynamic Case-2 water body characterized by high and spectrally dependent light attenuation with estuarine waters comprised of high loads of suspended inorganic sediment and dissolved organic matter from the Fraser River discharge (Komick et al., 2009; Loos and Costa, 2010; Phillips and Costa, 2017). For these highly attenuating waters, particulate scattering varies from approximately 0.2 to 16.0 m^-1 in the blue spectral region (Loos and Costa, 2010), and colored dissolved organic matter (CDOM) absorption varies from 0.007 to 3.07 m^-1 (Phillips and Costa, 2017). On the other hand, total attenuation in the northern SoG is predominantly controlled by absorption by CDOM and TChla (Loos and Costa, 2010; Loos et al., 2017; Phillips and Costa, 2017).

2.2 Data acquisition

2.2.1 In situ sata: water sampling and analysis

Surface water samples were collected at two regions of the SoG. For the first location, samples (N=108) were acquired along the BC ferry route aboard the Queen of Alberni (QA) crossing from Duke Point (Nanaimo) to Tsawassen (Vancouver; Stations 1-5, Figure 1) during the spring (March to May) and summer (June to September) of 2018 and 2019. Water samples were collected using a seawater intake pump installed on the ferry, representing a well-mixed layer of <3 m (Halverson and Pawlowicz, 2013; Wang et al., 2019; Travers-Smith et al., 2021). Triplicate samples from each station were collected for the different analyses (HPLC pigment concentration, CDOM absorption, and total suspended matter (TSM) concentration), kept in a dark environment, and immediately filtered onboard to minimize any potential compositional and concentration changes to phytoplankton pigments and dissolved organic matter (Mueller et al., 2003). For the second location, the northern SoG (Stations 6-7, Figure 1), surface water samples (~ 2 m depth) (N=19) were collected during summer (July) as part of an algal bloom monitoring component at aquaculture sites. Triplicate water samples were collected using a 5 L Niskin water sampler, transferred to a 10 L dark container, transported to the laboratory in the dark, and immediately filtered for HPLC pigments, CDOM absorption, and TSM. In addition, at the northern SoG location, single samples were collected for microscopic identification and enumeration of phytoplankton species (N = 8). Water for these samples was immediately transferred to 250 mL amber glass bottles and fixed using Lugol’s solution (1% concentration); samples were kept cool and dark until analysis.

Surface water samples for the TSM were filtered through precombusted and pre-weighed 47 mm Whatman GF/F 0.7 µm filter paper. Precombustion of the filter paper increases the efficiency by retaining more particulate matter as it shrinks the pore size to 0.3 μm during combustion (Nayar and Chou, 2003). Based on the color of the seawater, filtered volumes varied from 0.5 – 2 L. After the sample filtration, the filter paper was washed with 200 ml of distilled water to remove sea salt (Stavn et al., 2009), and immediately stored in dry ice, transported to the laboratory, and kept under -80°C until analysis. The filter paper was dried at 60°C for 24 hours in a clean oven and subsequently reweighed to 6 degrees of precision, at which point TSM was calculated (Röttgers et al., 2014). The final TSM concentration at each station was the average of the triplicates with a coefficient of variation (CV) < 25%.

Surface water samples to analyze a_CDOM (λ) were filtered through 47 mm Whatman 0.2 µm membrane filter paper. The filtrate was stored in acid-washed, precombusted (550°C for 1 hour) borosilicate glass bottles (Ferrari et al., 1996). Glass bottles were triple rinsed with sample water before the collection, and samples were kept in a cold environment, transferred to the laboratory and kept under 4°C until analysis. Prior to analysis, samples were allowed to reach room temperature to minimize temperature differences between the samples and the blank. The absorption of CDOM was measured using an Ocean Optics USB 4000 spectrophotometer (300-850 nm, 0.2 μm) against Milli-Q water as blank. The sample absorption was converted to CDOM absorption using Beer-Lambert’s Law equation (Ferrari et al., 1996).

A(λ) is the spectral absorption of the sample, and L is the cuvette path length in meters. The final a_CDOM (λ) for each station is the average of the triplicate with a CV < 25%.

Phytoplankton microscopy was performed by LCJL Marine Ecological Services using the Utermöhl method. Phytoplankton were identified and enumerated using 50 mL settling chambers and phase contrast microscopy with an inverted light microscope (Utermöhl, 1958). Count data were reported in cells L^-1. Counts from species within PFT groups used for analysis were summed to make comparisons with pigment-based estimates.

Surface water samples to measure HPLC pigments were filtered through 25 mm Whatman GF/F 0.7 µm glass microfiber filters under a low vacuum (200mm Hg/vac). Based on the color of the seawater, filtered volumes varied from 0.5 – 2 L. Filters were immediately put in dry ice and kept under -80°C until laboratory analysis. TChla and phytoplankton accessory pigments were analyzed using HPLC analysis done at the University of South Carolina Baruch Institute (https://phytoninja.com/lab-protocols/), following the methods described in Pinckney (2010). HPLC samples were collected in duplicates, and the uncertainty in the HPLC pigment dataset was computed from the duplicate measurements for 20% of the samples, resulting in a CV < 20%.

The HPLC pigment dataset was used as the input in the CHEMTAX software (v1.95) to derive the Chla concentration of each PFT towards the TChla concentration of the sample (Mackey et al., 1996). The selection of PFT and input pigment ratio were largely based on those in Del Bel Belluz et al. (2021), which were defined via in-depth analysis of pigments and pigment ratios, microscopic observations, and literature review of phytoplankton species in the SoG.

Unlike Del Bel Belluz et al. (2021), a raphidophyte group was included in the analysis to account for the extensive 2018 Heterosigma akashiwo blooms observed via microscopy at comparable sampling locations and times as those investigated in this study (Esenkulova et al., 2021). Additionally, our 2018 observations showed high violaxanthin concentration events where relationships with TChla (generally driven by diatom blooms), fucoxanthin (diatoms), Chlb (green algae), and zeaxanthin (typically only present in low concentrations and related to green algae or cyanobacteria) diverged from those through the rest of the dataset (Supplementary Material 1). During these high concentration events, fucoxanthin, violaxanthin and zeaxanthin ratios to TChla were comparable to those found for raphidophytes in the literature (Supplementary Material 2, Lewitus et al., 2005; Higgins et al., 2011). It is unlikely that other violaxanthin containing species would bloom to such high violaxanthin and TChla concentrations in the SoG. For reference, through the four-year time-series presented in Del Bel Belluz et al. (2021), which spanned multiple high TChla bloom events, violaxanthin reached a maximum concentration of 0.41 mg m^-3, an order of magnitude lower than observed during this study. For the central SoG data, pigment ratios for the raphidophyte group were taken from Lewitus et al. (2005), derived from Heterosigma akashiwo. In the northern SoG, pigment and microscopic observations showed that raphidophytes were absent. As a result, this group was not included in the northern SoG analysis. With these considerations, the following groups were considered in CHEMTAX analysis: cyanobacteria-1 (Cyano), haptophytes (Hapto), green algae (GA), cryptophytes (Crypto), dinoflagellates (Dino), raphidophytes (Raphido, central SoG only), dictyochophytes (Dictyo), and diatoms (Diat). Input ratios for each taxonomic group are given in Supplementary Material 3 and Supplementary Material 4. While quantitive evaluation of CHEMTAX is often lacking, many studies have utilized similar groupings to describe phytoplankton community compositions in coastal waters (e.g., Lu et al., 2018; Ferreira et al., 2020). Here, it was decided to attempt to maximize the retrievable groups; however, it is important to note that following the methods of Catlett and Siegel (2018), clustering of the HPLC-derived pigment dataset suggests that only four PFT groups may be separatable and optically discernable (Supplementary Material 9). As a result, CHEMTAX model outputs are supported by additional information, including microscopy (northern SoG) and data from complimentary studies done at similar locations and times.

CHEMTAX analysis was run on the northern and central SoG data separately to account for regional differences in pigment ratios as a result of differing oceanographic conditions (i.e., stratification, nutrients, and light; Suchy et al., 2019). An initial seeding run was first performed on the full dataset from each region to optimize the literature-derived ratios to the SoG (Armbrecht et al., 2015; Del Bel Belluz et al., 2021). Then, the central SoG data were binned into different groups based on TSM concentrations to account for the high variability in light conditions along the ferry track and between surveys. The following TSM bins were used: Low (TSM< 5 mg/L), Medium (TSM 5 – 11 mg/L), and High (TSM > 11 mg/L), based on bio-optical classes derived by Phillips and Costa (2017) for the central SoG. CHEMTAX analysis on the binned data utilized the output ratios from the initial seeding run performed on the full central SoG dataset. For the northern SoG, TSM conditions were comparable across the surveys, and HPLC data were not binned for analysis, but a second run, using the initial seeding ratios from the full northern SoG dataset, was performed to ensure stabilized ratios and to reduce RMSE.

Beyond consideration of environmental conditions, we evaluated pigment colinearity, which can affect CHEMTAX results (Higgins et al., 2011; Catlett and Siegel, 2018). In this study, positive correlations did exist between pigments used for analysis (Supplementary Material 5). When standardized to TChla, many of these correlations were still strong (e.g., Chlb and prasinoxanthin, Supplementary Material 5); however, none were correlated with TChla. Careful consideration was taken in the CHEMTAX design to limit error; however, pigment colinearity can add bias to the CHEMTAX model output and complicate retrieval via satellite (Higgins et al., 2011; Catlett and Siegel, 2018). As independent validation of CHEMTAX model outputs was not performed here, it was difficult to quantify this error, and this is further discussed in section 4.1. Care should be taken to extend results until extensive validation has been performed.

Following the CHEMTAX analysis, we also performed a Shannon-Weaver Index (Lohrenz et al., 2003; Latasa et al., 2010) on the CHEMTAX model output to assess the diversity of phytoplankton community composition:

where, P_i is the proportion of each phytoplankton group’s Chla relative to the TChla concentration (%), and N is the number of groups.

2.2.2 Satellite data

High spatial (300 m) and temporal resolution (daily) OLCI data were used in this study as this sensor has improved signal-to-noise ratios and an off-nadir tilted view to minimize the effect of sun glint (Donlon et al., 2012). Therefore, OLCI Level-1 Full-Resolution data were downloaded from the Sentinel 3 Marine CODA (Copernicus Online Data Access) web service, and processed with POLYMER (POLYnomial-based algorithm applied to MERIS) version 4.10 (Steinmetz et al., 2011) to obtain R_rs. POLYMER is an atmospheric correction (AC) algorithm that separates the atmospheric and sunglint spectral reflectance from the water reflectance using a spectral matching technique (Steinmetz et al., 2011). To do so, it uses two models: a simple polynomial used to model the spectral reflectance from atmospheric and residual sunglint, and a model that delineates ocean water reflectance using TChla concentration and the particle backscattering coefficient (Park and Ruddick, 2005). The water reflectance model used in this algorithm is valid in coastal and open ocean waters, including the bidirectional reflectance, and it provides better spatial coverage retrievals because it accounts for sky glint and thin cloud conditions (Steinmetz et al., 2011). Quality flags such as “Cloud,” “Invalid,” “Negative Backscattering (BB),” “Out-of-bounds,” “Exception,” and “High Air Mass” (Steinmetz et al., 2016) were applied following Giannini et al. (2021) recommendations for this region. According to the technical document (Steinmetz et al., 2016), the “Inconsistency” flag should not invalidate the results; however, our matchup analysis showed artifacts with pixels associated with this flag. As a result, data with this flag were removed from analysis. Giannini et al. (2021) demonstrated that POLYMER OLCI-derived products such as TSM and TChla showed the best performance compared to the Case-2 Regional CoastColor (C2RCC) processor for the coastal waters of BC and Northeast Alaska.

2.3 Data analysis: PFT satellite-based method and validation

The Chla concentration of PFTs was derived from OLCI satellite imagery using the EOF-based algorithm by assessing the variance or spectral shape in the OLCI R_rs spectra (Bracher et al., 2015). This algorithm was recently tested globally by Xi et al. (2020), and we further adapted it considering the local matchups between OLCI R_rs and the CHEMTAX-derived PFTs Chla. A schematic workflow of the EOF-PFT algorithm is shown in Figure 2. The general steps of this approach are explained as follows:

Figure 2 Schematic flowchart illustrating the EOF-based algorithm for predicting multiple Phytoplankton Functional Types (PFTs). The left side demonstrates the model training part with matchup comprised of CHEMTAX model output and Ocean Land Color Instrument (OLCI) remote sensing reflectance (R_rs; sr^–1(λ)). The right side indicates the model application to OLCI imagery to derive the spatial distribution of PFTs.

Step 1: OLCI R_rs data at 10 bands (400, 412, 443, 490, 510, 560, 620, 665, 674, 681 nm) were standardized by subtracting the spectral average and normalized by the spectral standard deviation according to Xi et al. (2020). The standardized R_rs was then subjected to a singular value decomposition (SVD) to derive the EOF scores (U), singular values (Λ), and EOF loadings (V).

The dataset used in step 1 corresponds to matchup satellite data obtained from a 3×3 pixel window (900 × 900 m) centered around the in situ sampling region and meeting the following criteria: ± 3 hours time difference, valid pixels are ≥5/9, CV at 560 nm ≤15%, and the median of the 3×3-pixel window at each band was used to avoid outliers (Bailey and Werdell, 2006; Mograne et al., 2019). Furthermore, to minimize the spatial-temporal mismatch between the in situ water sampling and the satellite overpass in a dynamic waterbody such as the SoG, samples collected from the Fraser River plume were relocated according to the surface current direction derived from the Ocean Networks Canada CODAR system (Coastal Ocean Dynamic Application Radar) (Halverson and Pawlowicz, 2016), following Giannini et al. (2021).

Step 2: A generalized linear model (GLM) was built for each PFT, considering the log-transformed CHEMTAX-derived PFT Chla and the corresponding EOF scores subset. Following Xi et al. (2020), insignificant EOF modes were excluded, and a stepwise routine was then adapted to search for the most significant regression model through the minimization of the Akaike information criteria (AIC). AIC is an estimator of prediction error, which uses the number of scores to determine the best model among the total in a given dataset (Bracher et al., 2015). Through the AIC criteria, only significant EOF scores were used in the regression model. As a result, the regression model for the PFT Chla prediction was expressed as:

Where C_p is the Chla concentration of each PFT, a is the intercept, b_1,2….n are the regression coefficients, and u_1,2,…n are the n number of EOF scores.

Step 3: As part of a level-1 validation, a cross-validation was performed following Bracher et al. (2015) and Xi et al. (2020). The cross-validation technique was applied to assess the robustness of the fitted regression model by splitting the whole dataset into two subsets; 80% of the data were used for model training and 20% for model validation. This procedure was run for 500 permutations, and during each permutation, EOF predicted, observed PFT Chla, and the statistical metrics for each permutation were recorded for evaluation.

The regression model performance was assessed by calculating the coefficient of determination (R²), Mean-Absolute Error (MAE), Median Percentage Difference (MDPD), and Bias. R², slope (S), and intercept (a) were obtained using the log-scaled (in natural logarithm) predicted versus log-scaled (in natural logarithm) observed PFT Chla, while the MAE, MDPD, and Bias were based on the non-log transformed data. The equations of these statistical metrics were given in Xi et al. (2020), except MAE, expressed as follows:

where Chla_pi is the predicted (using the regression model) PFT Chla, Chla_oi is the CHEMTAX observed PFT Chla, and M is the number of observations. Moreover, the mean value of all the statistical metrics for the cross-validation was recorded (R²cv, MAEcv, MDPDcv) to evaluate the robustness of the fitted regression model. This step comprises the first level of evaluation of the EOF-based PFT retrievals.

Step 4: The EOF-based PFT algorithm was applied to the OLCI imagery to map the spatial distribution of PFTs. To derive PFT Chla from satellite imagery for which we do not have any CHEMTAX or HPLC pigment data, we projected the standardized OLCI R_rs data onto the EOF loading and derived a new EOF score set. This new set of EOF scores was used in Equation 4 to retrieve the spatial distribution of PFTs. The regression coefficient and intercepts were taken from Equation 4. Finally, spatial maps were produced for TChla and the following PFTs: diatoms, cryptophytes, green algae, and raphidophytes. PFT selection was based on the level-2 validation derived using a high-resolution independent CHEMTAX time-series from the Hakai Institute (called Hakai time-series hereafter) acquired between spring 2017 and summer 2020 (N=26) at station QU39, located in the northern part of the SoG (–125.0992⁰W,50.0307⁰N; Figure 1). The independent matchup extraction criteria performed in the level-2 validation were the same as in the model training. The CHEMTAX result used for this validation was analyzed the same way as the CHEMTAX used in the model development, following Del Bel Belluz et al. (2021). Beyond the level-1 and level-2 quantitative validation, we also conducted a qualitative evaluation of the seasonal PFT maps in comparison with the literature.

3 Results

3.1 Biogeochemical characterization

Surface TChla was highly variable, showing remarkable seasonal and spatial variation (Table 1, Figure 3A), ranging from 0.39 to 28.82 mg/m³ with an exceptionally high value of 120.11 mg/m³ recorded at Station 4 during the summer of 2018. Outside of this value, all of the highest TChla were observed in 2018 with Stn 4 (28.82 mg/m³) and Stn 2 (21.08 mg/m³) showing high magnitudes during the summer and Stn 1 (18.26 mg/m³), Stn 2 (21.56 mg/m³), and Stn 3 (12.80 mg/m³) during the spring. Similarly, high seasonal variability was noticeable for TSM (Table 1, Figure 3B), especially during spring 2018 with values ranging from 1.69 to 65.63 mg/L. Persistently high TSM was observed mainly at stations close to the Fraser River; for example, Stn 4 (45.90 mg/L), Stn 2 (17.01 mg/L), and Stn 5 (13.48 mg/L) showed high TSM during spring 2018. In agreement with the TSM, a_CDOM also demonstrated high seasonal variation (Table 1, Figure 3C) during spring 2018, ranging from 0.21 to 3.45 m^-1. Relatively higher (>1.14 m^-1) a_CDOM was noted at stations close to Fraser River, such as Stn 2, Stn 3, Stn 4, and Stn 5, during spring 2018 and 2019.

Table 1 Seasonal surface biogeochemical variables in the Strait of Georgia (SoG) during 2018 and 2019.

Figure 3 Seasonal variation of (A) total chlorophyll-a concentration (TChla, mg/m³), (B) total suspended matter concentration (TSM, mg/L), and (C) the absorption by colored dissolved organic matter [a_CDOM, 443 (m^-1] in the Strait of Georgia (SoG) during 2018 and 2019. The upper and lower edge of the box plot represents the first and third quartile, respectively. Outside of the boxes, the whiskers on both quartiles indicate standard error. Blue points represent data points, and points outside the third quartile indicate extreme values.

Surface phytoplankton community composition derived from CHEMTAX for the central SoG (Stn 1 to 5) during 2018 and 2019 is shown in Figure 4. Generally, during spring (March to May) 2018, diatoms dominated the community composition in almost all the stations (54 ± 32%; average ± SD). At this time of the year, the relative contribution of flagellates, including haptophytes, green algae, cryptophytes, dinoflagellates, and dictyochophytes, was low (<14%) and followed by cyanobacteria and raphidophytes (< 7%). A seasonal succession from diatoms to flagellates persisted across the sampling stations during the summer (June to August) of 2018 (Figure 4), when the phytoplankton community composition was highly dominated by raphidophytes (24 ± 23%), followed by green algae (23 ± 10%) and cryptophytes (21 ± 16%). In 2018, high TChla concentration events dominated by raphidophytes were observed towards the western portion of the transect (Station 2) in June and concentrated near the mouth of the Fraser River (Stations 4 and 5) in July. Outside of these events, raphidophyte contributions often persisted through low TChla conditions.

Figure 4 CHEMTAX-derived Phytoplankton Functional Types (PFTs) in the central (Station 1 to 5) Strait of Georgia (SoG) for 2018 and 2019. The primary y-axis shows the relative contribution (%), and the secondary y-axis (black line with yellow circles) indicates the total chlorophyll-a (TChla; mg/m³) concentration.

Unlike the spring of 2018, in 2019, we observed the predominance of cryptophytes (37 ± 13%), followed by diatoms (28 ± 17%). A very similar trend was evident during summer, when cryptophytes (23 ± 16%) and diatoms (22 ± 25%) dominated the community composition, followed by green algae (20 ± 11%). Nevertheless, the relative contribution of other flagellates, such as haptophytes, dinoflagellates, raphidophytes, and dictyochophytes, towards the total phytoplankton concentration was minimal (< 8%). Similar to 2018, raphidophyte contributions were small, but somewhat persistent through low TChla conditions.

In the northern SoG (Figure 5), diatoms were the predominant group at Stn 6 (66 ± 19%) and 7 (24 ± 20%). At station 6, the average contribution of all other groups was minimal (< 18%). In contrast, Stn 7 showed higher contributions from dictyochophytes (17 ± 14%), cryptophytes (17 ± 5%), haptophytes (17 ± 8%), and green algae (16 ± 2%).

Figure 5 CHEMTAX-derived Phytoplankton Functional Types (PFTs) in the northern (Station 6 to 7) Strait of Georgia (SoG) during 2019. The primary y-axis shows the relative contribution (%), and the secondary y-axis (black line with yellow circles) indicates the total chlorophyll-a (TChla; mg/m³) concentration.

After the spring bloom, a diverse phytoplankton community comprised of flagellates evolved and is summarized by the Shannon-Weaver Index (Supplementary Material 6). In general, the central SoG data showed increased diversity during the summer (1.5 ± 0.3) when diatom contributions were lowest (12 ± 21%). Similarily, in the northern SoG, diversity was higher at Stn 7 (1.6 ± 0.9), which had lower diatom contributions (24 ± 20%) (Supplementary Material 7).

The comparison of the CHEMTAX-derived PFTs with the phytoplankton abundance from microscopy showed significant relationships for diatoms (R² = 0.87, p< 0.001, n = 8) and dinoflagellates (R² = 0.52, p<0.05, n = 8). This is an important result that validates the CHEMTAX analysis; however, it is important to note that the sample size for these comparisons was small. Further, pico-sized species likely constituted many of the groups showing high contributions and these would have largely been missed by microscopy and many groups showed low variability in TChla contributions (small dynamic range).

3.2 OLCI derived PFTs

3.2.1 EOF-based analysis of the OLCI R_rs data

The spectral shapes and magnitude of OLCI R_rs and the corresponding standardized R_rs at 10 bands (400, 412, 443, 490, 510, 560, 620, 665, 674, 681 nm) are given in Figure 6. These OLCI R_rs matchup data were then decomposed using SVD to derive nine EOF modes (Table 2), in which the first four modes explained 99.48% of the variance. The spectral distribution of each EOF and the spectral features associated with several seawater constituents are shown in Figure 7. In general, EOF-1 explained 60.33% of the total variability, and exhibited a high reflectance peak at 560 nm and a 680 nm feature likely associated with phytoplankton and TChla fluorescence, respectively. EOF-2 explained 38.33% of the total variability and exhibited a sharp decrease from 500 to 400 nm, indicating the influence of high CDOM absorption in the blue region, and a broad peak from 560 to 700 nm, indicating high backscattering from inorganic suspended sediments. Moreover, EOF-3 explained 0.50% of the variability, showing a decrease in spectral structure from 400 to 450 nm and a spectral peak at 560 and 680 nm. EOF-4 explained only 0.32% of the total variance; however, a trough is evident from 443 to 510 nm likely associated with absorption by different accessory pigments, followed by several peaks at 412, 443, 560, and 681 nm. The spectral peaks observed at 560 and 681 nm could be associated with phytoplankton and TChla fluorescence, respectively. The higher EOFs, such as EOF-5, EOF6, EOF7, EOF8, and EOF9, only showed variance < 0.20%.

Figure 6 (A) Remote sensing reflectance (R_rs; sr^–1(λ)) and (B) corresponding standardized R_rs subjected to the Empirical Orthogonal Function (EOF) analysis. The grey lines show the median spectra from a 3×3 pixel window with the corresponding average spectra in red, and the standard deviation at each band is shown in black.

Table 2 Percentage of variance explained by Empirical Orthogonal Function (EOF) modes derived from the matchup between Ocean Land Color Instrument (OLCI) remote sensing reflectance [R_rs; sr^–1(λ)] and CHEMTAX-derived Phytoplankton Functional Types (PFTs).

Figure 7 Spectral loadings for the nine Empirical Orthogonal Function (EOF) modes derived from the Ocean Land Color Instrument (OLCI) remote sensing reflectance [R_rs; sr^–1(λ)]. This spectral loading was obtained from the matchup between OLCI R_rs and total chlorophyll-a (TChla; mg/m³) concentration.

The ∆AIC (Supplementary Material 8) showed that EOF-4 was the most significant for TChla, diatoms, and dictyochophytes, and higher EOFs were selected for most PFTs (Supplementary Material 8). At the same time, the second most important EOFs differed among PFTs (Supplementary Material 8). For cyanobacteria and dinoflagellates < 3 EOFs were selected to build the regression model (Supplementary Material 8).

3.2.2 Model evaluation

Figure 8 shows the regression between the CHEMTAX-derived PFTs versus TChla concentrations and the Chla concentrations of the eight PFTs retrieved from OLCI. Overall, the regression models showed moderate and highly significant retrieval (p<0.0001) for diatoms and raphidophytes (R² = 0.45, 0.37, respectively; Table 3). In addition, the regression models for dictyochophytes, haptophytes, green algae, and cryptophytes were also significant (p<0.001; Table 3); however, relatively lower R² (≤0.30) were observed. Poor performances were observed for dinoflagellates and cyanobacteria (R² = 0.07, 0.11, respectively; Table 3), with the best performance (R² = 0.51; p<0.0001) for TChla retrievals. The lowest MDPDs were observed for TChla, green algae, haptophytes, and cyanobacteria (34.26%, 41.18%, 44.18%, 44.29%, respectively), whereas the highest MDPDs were noted for raphidophytes, followed by dinoflagellates (66.78%, 60.27%, respectively; Table 3). MAE was relatively higher for diatoms (MAE = 0.23 mg/m³); however, it was <0.23 mg/m³ for all other PFTs.

Figure 8 Regression between CHEMTAX-derived (x-axis) versus OLCI-retrieved (y-axis) Chla concentration of (A) Cryptophytes, (B) Cyanobacteria, (C) Diatoms, (D) Dictyochophytes, (E) Dinoflagellates, (F) Haptophytes, (G) Green Algae, (H) Raphidophytes, and (I) total chlorophyll-a (TChla).

Table 3 Full-fit statistics and Cross-validation statistics derived from the EOF-based algorithm, developed using the matchup between Ocean Land Color Instrument (OLCI) remote sensing reflectance [R_rs; sr^–1(λ)] and CHEMTAX-derived Phytoplankton Functional Types (PFTs).

As part of the level-1 evaluation, the robustness of all PFT regression models was further assessed with the cross-validation approach (Table 3). The statistical parameters obtained from the cross-validation procedure were poorer than the entire dataset statistics. For example, a poor R²_CV was obtained for almost all PFTs, except TChla (0.43). In general, R²_CV ranged from 0 to 0.43 (Table 3). Whereas, MDPD_cv ranged from 39.12 to 74.17%, in which the highest MDPD_cv was obtained for raphidophytes (MDPD = 74.17%), followed by diatoms (65.77%). However, a low MDPD_cv (39.12%) was observed for TChla.

The level-2 evaluation (Figure 9) was conducted for OLCI-derived outputs against the independent Hakai time-series from station QU39 (N = 26, Figure 1). Of the considered groups, green algae, cryptophytes, and diatoms followed similar trends through the Hakai time-series, resulting in significant (p <0.05) retrieval (R² = 0.33, 0.23, 0.20, MDPD = 163, 46.86, 88.59%, respectively) (Figure 9; Table 4). The best retrieval was noted for TChla, which highly corresponded with the Hakai time-series (R² = 0.68, MDPD = 52.25%) (Figure 9; Table 4). Except for TChla, diatoms recorded the highest MAE (0.58 mg/m³; Table 4), and MAE for all other groups was < 0.58 mg/m³. The MDPD was the highest for green algae (163%) and raphidophytes (130%), whereas MDPD was relatively lower for cryptophytes (46.86%). Overall, OLCI-derived PFTs, including TChla, over-estimated the lower ranges and under-estimated higher ranges of in situ values. Raphidophytes showed poor correspondence (R² = 0.08, MDPD = 130%) with the Hakai time-series (Figure 9; Table 4); however, it is important to note that raphidophytes were not observed via microscopy on any of the considered dates, highlighting misclassifications by both CHEMTAX and OLCI-derived outputs under low TChla conditions.

Figure 9 Comparison of the performance of the OLCI-derived (A) Diatoms, (B) Cryptophytes, (C) Green Algae, (D) Raphidophytes, and (E) total chlorophyll-a (TChla) concentration obtained from the daily cloud-free images between spring 2017 to summer 2020 (blue diamonds) against the Hakai time-series data at station QU39 located in the northern Strait of Georgia (SoG) (red circles). The location of station QU39 is given in Figure 1. The vertical lines on each panel separate different years, and note that the matchup data points differ among Phytoplankton Functional Types (PFTs), mainly because all PFT may not be present at all stations. For example, Raphidophytes showed fewer points (19) out of the 26 matchup points, followed by Diatoms (22). Hence, the dates on the x-axis might not be similar for all PFTs.

Table 4 Validation statistics were obtained from the regression between OLCI-derived Phytoplankton Functional Types (PFTs), including total chlorophyll-a (TChla; mg/m³) versus Hakai time-series data.

3.3 PFT mapping in SoG from daily OLCI imagery

After the model evaluation, we applied the EOF-based regression algorithm to daily OLCI imagery acquired during cloud-free days between March and September 2018 to provide a qualitative evaluation based on the spatial distribution and seasonal dynamics of PFTs over the SoG. The spatial maps (Figure 10) were produced for diatoms, cryptophytes, green algae, and raphidophytes since they accounted for 74% of the total phytoplankton concentration. Overall, the PFT spatial maps showed a distinct spring and fall bloom in the SoG, which was mostly comparable with the Hakai time-series. During spring, diatoms showed the highest average concentration (avg 5.15 mg/m³) with high concentrations observed near the northern and central SoG, and the Juan de Fuca Strait (Figures 10A1-A2), similar to the Hakai time-series data (diatom Chla ranging from 0.23 to 6.00 mg/m³) in the northern SoG. In comparison, the spatial distribution of cryptophytes and raphidophytes for the spring followed a similar spatial distribution to diatoms, but with relatively lower concentrations (avg 0.70, 1.30 mg/m³, respectively; Figures 10B1-B2, Figures 10D1-D2). Furthermore, the spatial distribution of green algae for the spring showed a comparatively low average concentration (0.32 mg/m³) and tended to show inverse spatial trends to diatoms. Hakai’s data for this week of the year also showed comparable concentrations for cryptophytes, green algae, and raphidophytes in the northern SoG.

Figure 10 Spatial and seasonal distribution of Diatoms, Cryptophytes, Green Algae, Raphidophytes, and total chlorophyll-a (TChla), retrieved from Ocean Land Color Instrument (OLCI) over the surface waters of the Strait of Georgia (SoG). Only images with >80% coverage were selected for the final map. OLCI imagery was acquired on cloud-free days from spring to fall of 2018, and the specific days were given in the spatial map. White regions in the maps show the clouds.

Following the spring period, the summer average diatom concentration decreased (4.37 mg/m³), but strong diatom bloom events were still observed off of the southwestern coast of Vancouver Island, with patches of high contributions in the central and northern SoG near the beginning of summer (June 16^th – Figure10-A3). By July 25^th, diatom contributions had decreased considerably (Figure10-A4). Consistent with these observations, Hakai’s data for a similar period showed low diatom concentrations in the northern SoG, ranging from 0.12 to 0.14 mg/m³. Furthermore, comparable with the Hakai time-series, OLCI-derived green algae showed a relatively higher average concentration (0.37 mg/m³; Figure 10-C3) than spring with localized areas of high cryptophyte concentrations in the central SoG and Juan de Fuca Strait (Figure 10-B3). Interestingly, raphidophytes recorded the highest average concentration during the summer (2.24 mg/m³), and blooms were observed off the southwest coast of Vancouver Island and close to the Fraser River mouth (Figures 10-D3, 10-D4). Of note, on June 16^th, both diatoms and raphidophytes showed high concentrations off of the southwest coast of Vancouver Island; however, concentrations between the two groups was decoupled within the central and northern SoG, where diatom concentrations were high, and raphidophyte contributions were low. Finally, in fall, diatoms still constituted the highest average concentration (1.35 mg/m³), especially near the central and northern SoG and off the south of Vancouver Island (Figure 10-A6). During this season and, similar to the Hakai time-series, the spatial distribution of cryptophytes was similar to the distribution of diatoms but with a relatively low average concentration (0.49 mg/m³). In turn, green algae showed relatively high concentrations in the central and northern SoG and off of southern Vancouver Island (Figure 10-C6) resulting in the highest average seasonal concentration (0.41 mg/m³) for this PFT. Overall, TChla showed high concentrations throughout the season with high concentration regions localized in the central and northern SoG, near the mouth of the Fraser River, within the Juan de Fuca Strait, and off the south coast of Vancouver Island (Figures 10-E1, 10-E2, 10-E3, 10-E4, 10-E5, 10-E6).

4 Discussion

4.1 EOF-based PFT retrievals and source of uncertainties

EOF analysis has successfully been used to derive phytoplankton cell abundance and the concentration of accessory pigments and multiple PFTs from the global open ocean and regional waters (e.g., Craig et al., 2012; Bracher et al., 2015; Lange et al., 2020; Xi et al., 2020); however, few of these studies have addressed coastal Case-2 waters (e.g., Craig et al., 2012). In these waters, the relationship between EOFs to specific PFTs can be hindered by high concentrations of water constituents such as CDOM, detritus, and TSM, which can largely define the reflectance spectra and mask PFT signals (Mouw and Yoder, 2010; Craig et al., 2012; Mouw et al., 2017). In this study, high variability of these constituents were observed (Table 1), thus significantly contributing to the reflectance signal and EOF modes (Figure 7), as observed in Phillips and Costa (2017) and Craig et al. (2012). Furthermore, EOF modes from oceanic waters (Bracher et al., 2015; Xi et al., 2020) show high dissimilarity to those found in this study, likely due to the absence of high concentrations of non-phytoplankton constituents. As a result, the derivation of PFTs in this study was likely hindered by the high bio-optical complexity characteristic of the SoG (Loos and Costa, 2010; Loos et al., 2017; Phillips and Costa, 2017; Giannini et al., 2021). This difficulty would have been most pronounced for low biomass groups which showed weak signals.

Beyond the inherent role of the Case-2 waters’ different optical constituents on the EOF retrievals, the biomass contribution of a specific PFT to the total phytoplankton biomass is known to be relevant for accurate PFT retrievals from spectral data. Of the two levels of validation, the first validation showed good retrievals for diatoms and raphidophytes, and the second validation showed that green algae, cryptophytes, and diatoms followed the seasonal trends of their corresponding in situ values from a high temporal resolution Hakai time-series (Figure 9). This is likely associated with their larger contribution (62%) to the total concentration, as demonstrated in the CHEMTAX model output (Figure 4). These results are encouraging with regard to the remote sensing of PFTs over optically complex Case-2 waters, where previous studies have documented that PFT retrieval over Case-2 waters is challenging (Craig et al., 2012; IOCCG, 2014; Mouw et al., 2017). At the same time, the poorest retrievals obtained for haptophytes, dictyochophytes, cyanobacteria, and dinoflagellates was likely related to their low concentration (26%), resulting in low contributions to the R_rs signal. Similarly, Xi et al. (2020) also observed poor performance for low biomass prokaryotes (R²cv = 0.11; MDPDcv = 55.08%) and Prochlorococcus (R²cv = 0.18; MDPDcv = 42.68%), and Bracher et al. (2015) showed slightly poorer retrievals (R²cv = 0.28; MDPD = 31%) for the pigment zeaxanthin (found in low concentrations), a diagnostic pigment for cyanobacteria.

Other factors, such as the quantity and quality of in situ and satellite matchup data also impact the performance of satellite-derived PFTs. For the EOF-based PFT algorithm, the number of matchup data points and their representation of the range of phytoplankton groups in the region are crucial to enhance the accuracy of the retrievals (Bracher et al., 2015; Xi et al., 2020). Bracher et al. (2015) suggested that at least 50 valid matchup data points are necessary for statistically significant pigment estimation, specifically for open ocean waters. Xi et al. (2020) employed a larger number of matchup data points (52-394) for global open ocean waters, with an even higher number (483) in a more recent publication (Xi et al., 2021). Lange et al. (2020) employed 73-78 matchup points to develop an EOF-based algorithm to retrieve phytoplankton cell abundance from MODIS across the Atlantic Ocean. However, Craig et al. (2012) showed that a reliable EOF-hyperspectral-based model could be built with 15 valid data points, specifically for one season and at a local scale. In our study, the total number of matchup data points was between 61 and 102, which is within the range suggested by the above authors. However, given the optical complexity of our waters, a higher number of matchups would likely improve the PFT retrievals.

The spatial and temporal variability of the matchup samples, which has implications on data quality, is also an important factor in evaluating the EOF-based PFT retrievals, especially in Case-2 dynamic coastal waters such as the SoG. This is a common problem in validating satellite-derived products of coastal waters (e.g., Mahadevan and Campbell, 2002; Moses et al., 2016; Barnes et al., 2019; Tilstone et al., 2022). Specifically in the SoG, Nasiha et al. (2022) have shown large spatial variability of in situ above-water reflectance within the 300m spatial resolution of OLCI. Notably, most of this variability was observed at the frontal region of the Fraser River and within transitional oceanic waters and likely has considerable implications on the quality of in situ data for satellite validation purposes. Further, Giannini et al. (2021) have shown the importance of the temporal mismatch between in situ and OLCI data in the SoG due to the dynamic nature of this water system. Here, an attempt was made to reduce these spatial issues by rejecting matchups corresponding to transitional waters, and further, to minimize temporal mismatch by relocating the sample’s position using surface current maps generated using a CODAR system (Halverson and Pawlowicz, 2016). This method was successfully adopted by Giannini et al. (2021). Out of the 102 matchup samples utilized in this study, 20 were successfully relocated using the surface current map. Indeed, after the relocation, we observed improved performance in the R_rs matchup analysis (not shown). However, it is still challenging to detect and quantify any remaining source of uncertainties associated with this spatial and temporal water mass mismatch.

Uncertainties associated with CHEMTAX model outputs should also be considered in evaluating the EOF-based PFT retrievals. For example, CHEMTAX outputs are more prone to errors in environments with complex phytoplankton community composition, notably when local ratios are unknown, and many species have similar pigment composition (Lewitus et al., 2005); this resembles some of the conditions in the Salish Sea. For instance, the haptophyte, Phaeocystis pouchetti, has been shown to periodically lack 19′-hexanoyloxyfucoxanthin (19HF) and have relatively high 19′-butanoyloxyfucoxanthin (19BF) concentrations; in these conditions, the utilized CHEMTAX ratios may result in misclassification of the haptophyte as a dictyochophyte (Del Bel Belluz et al., 2021). In addition, raphidophytes, such as Heterosigma akashiwo, often bloom regionally and have similar pigment ratios to diatoms, with the separation of raphidophytes dependent on their high violaxanthin to TChla ratios. However, violaxanthin concentrations are variable and broadly present among groups (Lewitus et al., 2005; Laza-Martinez et al., 2007; Higgins et al., 2011). In this study, violaxanthin to TChla ratios were noticeably high during raphidophyte blooms, and CHEMTAX appeared to be able to distinguish these events. Yet, under lower TChla concentration conditions, raphidophyte contributions often persisted, likely representing misclassifications and underestimations of other groups. These misclassifications are evidenced by the small CHEMTAX and satellite-based raphidophyte contributions shown over the Hakai time-series comparison (Figure 9); no raphidophytes were observed via microscopy during much of this time-series (Del Bel Belluz et al., 2021). As such, errors from CHEMTAX may have been perpetuated to the satellite models.

To address this issue, novel “data driven” approaches have been developed to help account for PFTs present in low contributions and also, to deal with issues in CHEMTAX analysis such as misclassifications of groups with similar pigment profiles and pigment collinearity (Catlett and Siegel, 2018; Kramer and Siegel, 2019; Kramer et al., 2020). These approaches, which are not hindered by co-correlation between pigments, use both cluster and EOF analysis to group species with similar pigment ratios. Outputs from these methods typically have fewer PFT groups (< 5) that combine communities that covary in the environment. Using the same clustering methods, the pigment data in this study clustered into 4 pigment groups: diatoms, dinoflagellates, haptophytes/pelagophytes, and green algae/cryptophytes/cyanobacteria (Supplementary Material 9). The use of this more generalized approach in our study may have allowed for a more accurate retrieval of a broad “small flagellate” group during high diversity low TChla summer conditions, similar to those used effectively by regional particulate fatty acid studies (Costalago et al., 2020; McLaskey et al., 2022). However, the results of this clustering approach grouped violaxanthin with green algae pigments suggesting that this group drove its variability, which is generally the case under non-raphidophyte bloom conditions. As a result, this approach would not have allowed for the separation of raphidophytes, a toxic species with considerable regional impacts on fish (Esenkulova et al., 2021), but likely helps explain the misclassifications described above.

In order to fully assess the above error and decide on the best approaches, independent validation of in situ pigment-based PFT results are needed: this is a major shortcoming of many pigment-based studies (e.g., Chase et al., 2020; Kramer et al., 2020). In our study, regional CHEMTAX results corresponded well with expected trends in phytoplankton communities (discussed below), qualitative microscopy comparisons (Del Bel Belluz et al., 2021), and particulate matter fatty acids and size-fractionated Chla (SF-Chla) (Costalago et al., 2020; McLaskey et al., 2022); however, no direct cross comparisons have yet been performed. This lack of comparison is partially due to the difficulty in comparing pigment-based approaches with other measures of phytoplankton composition (Higgins et al., 2011). In particular, microscopy is skewed to larger species and cannot resolve pico-sized species known to often dominate summer conditions in the SoG (Del Bel Belluz et al., 2021). This trend is evidenced by the small and limited CHEMTAX validation dataset used here, which showed correspondence with only more easily preserved diatoms. Additional sources of validation, using flow cytometry, size-fractionated filtration (SFF), conventional imaging flow cytometry (Olson and Sosik, 2007), Imaging FlowCytobot, genomics, and next-generation sequencing (see IOCCG, 2014 and Lombard et al., 2019 for a detailed review) are needed to evaluate pigment-based community methods and satellite outputs.

Finally, the quality of atmospheric correction and its effect on the OLCI-measured signal is a consideration in evaluating the EOF-based PFT retrievals. For oligotrophic oceanic gyre waters, R_rs uncertainties in blue wavelengths are above 5% (Hu et al., 2013), and are even higher in optically complex waters (Moore et al., 2015). A recent study by Giannini et al. (2021) assessed the radiometric qualities of OLCI R_rs in northeast Pacific coastal waters, showing the highest errors associated with the blue bands (400, 412, 443 nm) when compared with autonomous above-water in situ reflectance measurements (Wang and Costa, 2022). A significant R_rs underestimate in the blue and red bands of OLCI was also reported in the European coastal and inland waters (Zibordi et al., 2018; Bi et al., 2018; Li et al., 2019; Alikas et al., 2020; Shen et al., 2020; Vanhellemont and Ruddick, 2021). Although POLYMER has been proven to outperform other methods to retrieve OLCI R_rs for the coastal waters of BC, Giannini et al. (2021) argue that uncertainties exist and are likely associated with the optimization of the model, which considers TChla and a coefficient (fb) that scales the particulate backscattering in the water column, but does not consider the contribution of CDOM (Steinmetz et al., 2016). As previously discussed, CDOM is a significant optical constituent affecting the light spectra in the SoG (Loos and Costa, 2010; Phillips and Costa, 2017, Table 1), with a_CDOM in our study ranging from 0.11 to 3.45 m^–1 (Table 1). Since the launch of the Sentinel 3A, the uncertainties found between the products derived from OLCI using different processing approaches, especially for coastal waters, have been under critical evaluation and are still an active area of research (e.g., Zibordi et al., 2018; Hieronymi, 2019; Mograne et al., 2019; Gossn et al., 2019; Giannini et al., 2021; Vanhellemont and Ruddick, 2021; Tilstone et al., 2021).

4.2 Dynamics of PFTs in the SoG

The OLCI-derived PFT maps showed strong spring diatom blooms with the highest concentrations observed in the central and northern SoG (Figure 10-A1). These findings are in-line with TChla concentrations and CHEMTAX, which showed that diatoms dominated the phytoplankton community in the spring of 2018 (Figure 4). For the same period, in situ observations (microscopy, CHEMTAX, particulate matter fatty acids, and SF-Chla) collected across the SoG showed a similar spring bloom timing with the community dominated by the centric diatoms Thalassiosira, Skeletonema, and Chaetoceros (Nemcek et al., 2018; Del Bel Belluz et al., 2021; Esenkulova et al., 2021; McLaskey et al., 2022). In the SoG, diatom-dominated spring blooms are well documented in the literature (e.g., Harrison et al., 1983; Allen and Wolfe, 2013; Suchy et al., 2019; Del Bel Belluz et al., 2021). As in many temperate systems, diatoms are the primary group that quickly respond to increased nutrient and light conditions resulting in rapid accumulation and bloom conditions (Allen and Wolfe, 2013). Both satellite-derived time-series data (Suchy et al., 2019) and coupled-biophysical model results (Collins et al., 2009) have shown that spring bloom timing and strength are related to early water column stratification due to positive anomalies in SST, Fraser River discharge, and reduced wind stress. These conditions promote the relatively high average surface TChla concentration in the SoG, which is one of the most productive regions on the west coast of North America (Jackson et al., 2015).

During summer, EOF-based PFT retrievals generally showed reduced TChla concentrations and diatom contributions, with exception to localized events off of the southwest coast of Vancouver Island and in the central and southern SoG. Furthermore, the northern SoG typically showed reduced TChla and diatom contributions compared to the central and southern regions (Figures 10-A4, 10-E4). This succession to lower concentration summer conditions is expected for the region. For instance, a monthly climatology of satellite-derived surface TChla over 13 years (2003 to 2016) showed relatively lower TChla during summer and reduced concentrations focused in the northern SoG (Suchy et al., 2019). Spatially, the northern SoG is more stable and stratified and prone to limiting nitrate concentrations, whereas the southern SoG is exposed to higher vertical mixing bringing nutrient-rich water from the subsurface to the sunlit layer (Peña et al., 2016). In contrast, higher phytoplankton concentrations found to the south of Vancouver Island have been linked to wind-induced seasonal upwelling, freshwater discharge from Fraser and Columbia Rivers, and the outflow of nutrient-rich low salinity water from the Juan de Fuca Strait (Masson and Cummins, 2000; Hickey and Banas, 2008; Thomson et al., 2014). Indeed, the brackish water outflow from the Juan de Fuca Strait is a significant source of nitrate to the continental shelf of Vancouver Island and northern California Current (Foreman et al., 2008; Hickey and Banas, 2008). These nutrient-rich waters supply adequate nitrate concentration to the mixed layer depth promoting increased phytoplankton biomass and primary productivity to the southwest of Vancouver Island during summer (MacFadyen and Hickey, 2010).

The summer succession to flagellate-dominated communities in the SoG was apparent in the in situ CHEMTAX model output and less so in the OLCI-derived spatial maps. In agreement with our observations, many authors using a variety of approaches have reported a similar seasonal succession from diatoms to small flagellates and noted that the northern SoG has a more diverse flagellate dominated community when compared to the central SoG (e.g., Harrison et al., 1983; Haigh and Taylor, 1991; Del Bel Belluz et al., 2021; McLaskey et al., 2022). The Shannon-Weaver Index calculated from our CHEMTAX model output supports the literature showing higher diversity in the summer when small flagellates were dominant (such as green algae and cryptophytes, followed by dinoflagellates, haptophytes, and dictyochophytes), and relatively higher diversity in the northern SoG than the central SoG (Supplementary Material 7).

The strong raphidophyte blooms observed in summer 2018 were an exception to the seasonal succession to a diverse range of small flagellates. These blooms were observed by multiple other studies in the region (Nemcek et al., 2018; Pawlowicz et al., 2020). In particular, Esenkulova et al. (2021) found June and July 2018 to be an exceptional year in terms of Heterosigma akashiwo, with abundance reaching 11,000 cells mL^-1, compared to data from 2015-2017 (150 cell mL^-1). Similar to the results of this study, these authors showed that in the SoG, Heterosigma akashiwo are typically observed from May to September, with the summer 2018 blooms having a distinct spatial distribution concentrated near the Fraser River mouth (Figure 10 D4). The location of these blooms near the Fraser River suggests that Heterosigma akashiwo is strongly associated with highly stratified low salinity waters with low nitrate and phosphate concentrations and high silicate concentrations (Esenkulova et al., 2021). In addition, and similar to the results of this study, Heterosigma akashiwo was observed off of the southwest coast of Vancouver Island (Figure 10-D3) on June 16^th, 2018 by Gower and King (2018) using the Maximum Chlorophyll Index (MCI) method on satellite imagery and these results were corroborated with microscopy. In the SoG, Heterosigma akashiwo are thought to have negative impacts on wild Pacific salmon populations and significantly impact the fin-fish and salmon aquaculture industries driving large fish kill events (Rensel et al., 2010; Haigh and Esenkulova, 2014). As a result, the ability to monitor these species using EOF-based PFT maps would provide a mechanism for better understanding and mitigating the negative effects of these harmful species.

Lastly, a fall bloom dominated by diatoms is evident in our spatial maps, which is typical for this region (Nemcek et al., 2018; Del Bel Belluz et al., 2021). This bloom was largely constituted by Pseudo-nitzschia and was observed spanning much of the northern SoG (Nemcek et al., 2018; Esenkulova et al., 2021; Del Bel Belluz et al., 2021). On a broader temporal scale, Jackson et al. (2015); Peña and Nemcek (2018), and Suchy et al. (2019) have reported fall blooms on the Vancouver Island coast. Generally, the primary driver for the 2018 fall bloom in the SoG was the reintroduction of nutrients into the sunlit layer following high summer stratification and grazing (Del Bel Belluz et al., 2021). In turn, autumn blooms to the south of Vancouver Island have been attributed to the seasonal wind-induced upwelling followed by freshwater discharge from the Columbia River or Fraser River and the nutrient-rich water outflow from the Juan de Fuca Strait (Jackson et al., 2015).

5 Conclusion and outlook

This is the first study to characterize the spatial-temporal dynamics of the major PFTs in the SoG, an optically complex Case-2 system on the west coast of Canada, with satellite data using an EOF-based algorithm combined with CHEMTAX model outputs. As a result of the complexity of these waters, many of the low concentration groups such as haptophytes, dictyochophytes, dinoflagellates, and cyanobacteria showed poor performance. However, the first level of validation demonstrated reliable performance for diatoms and raphidophytes, and the second level of validation with an independent Hakai time-series showed the best retrieval for TChla and green algae, cryptophytes, and diatoms showed similar seasonal trends as their corresponding in situ values. In this comparison, raphidophyte outputs by both CHEMTAX and satellite retrievals suggested misclassifications during the observed non-raphidophyte bloom conditions over the Hakai time-series when no raphidophytes were observed via microscopy.

Qualitative validations of the spatial-temporal 2018 EOF-based PFT generally agreed with expected patterns for this region with strong spring bloom highly dominated by diatoms, and primarily located in the central and northern SoG and Juan de Fuca Strait. During summer, diatoms showed high concentration off of the southwest coast of Vancouver Island and in patchy regions within the central SoG. As expected based on regional studies, summer TChla and diatom contributions decreased with the lowest contributions observed towards the northern SoG. In turn, patches of increased cryptophyte contributions were noticeable in the central SoG and Juan de Fuca Strait from spring through summer. Green algae, which generally shows its highest concentrations in the summer through autumn, showed low concentrations during spring throughout the central and northern SoG. Furthermore, similar to the CHEMTAX model output and literature, the EOF-based PFT maps revealed 2018 summer raphidophyte blooms (Heterosigma akashiwo) off of the southwest coast of Vancouver Island, within the Juan de Fuca Strait and close to the Fraser River mouth. Finally, a fall bloom was evident in our spatial maps in which diatoms displayed high concentration near the central and northern SoG and south of Vancouver Island. Overall, diatoms demonstrated high average concentrations throughout the study period, and relatively low average concentrations were observed from green algae, cryptophytes, and raphidophytes.

The results of this study suggest that retrieving PFTs in the SoG is highly challenging, likely as a result of the strong contributions of CDOM and particulate suspended matter to the satellite-measured reflectance signal. Further work is required to improve the algorithm and should rely on the following: better sensor vicarious calibration and atmospheric correction of the measured radiance signal; considerations of in situ data quality, including analytical errors; use of broader data-driven pigment methods and quantitative in situ validation utilizing a variety of methods; representation of the phytoplankton group’s dynamic range, and; improving the spatial-temporal quality of the matchups. Also, considering that the retrieval algorithm is an emerging technique, further research is necessary to investigate the over-estimation of the lower ranges and under-estimation of higher ranges of PFT Chla, and to promote the introduction of non-linear fitting in the regression model. Finally, another limitation of our EOF-based PFT analysis is the lack of estimation of per-pixel uncertainty, which to our knowledge has only been reported by Brewin et al. (2017) and Xi et al. (2021). In our study, we used OLCI R_rs as the spectral input data derived through the POLYMER algorithm, and the uncertainty for the POLYMER derived R_rs is currently under development (Steinmetz, personal communication). Hence, performing per-pixel uncertainty of our OLCI-derived PFTs is currently not feasible. Besides, the advance of POLYMER is still an ongoing field, so future studies that aim to elucidate PFTs using POLYMER R_rs should be able to justify the per-pixel uncertainty estimates of PFTs using an updated version of POLYMER, especially over coastal waters.

Despite the shortcomings and uncertainties of our satellite product, our study shows the potential of OLCI for deriving major PFTs over a large spatiotemporal scale, which are not achievable via ship-based measurements. With continued research, the methods deployed here will greatly help to further our understanding of phytoplankton dynamics in complex coastal regions such as SoG and aid in the management of critical species such as migrating Pacific Salmon.

Data availability statement

All analysis and model establishment were performed with R programming. The core functions used in the EOF model training part (Figure 2) in section 2.3 steps 1-2 are `SVD' (singular value decomposition) for EOF analysis, `GLM' (generalized linear models) for multilinear regression between EOF modes and PFTs, and `step' for the stepwise routine to remove insignificant predicting terms. The numeric matrices of V ∙ Λ-1 and the regression coefficients determined by the training part were used for the EOF model applications, and they are provided in Supplementary Materials 10 and 11. The R script on the model application procedure to estimate different PFT Chla is available at https://github.com/vishnu-star/Mapping-Phytoplankton-Functional-Types.git. The in situ and satellite data are presented in the article/Supplementary Materials, and further inquiries can be directed to the corresponding author.

Author contributions

PS: Conceptualization, data collection, analysis, methodology, figure preparation, and writing of the manuscript. HX: methodology, reviewing, and editing. JB: methodology, writing, reviewing, and editing. MH: methodology. AB: methodology, reviewing, and editing. MC: conceptualization, methodology, reviewing, and editing. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by NSERC NCE MEOPAR - Marine Environmental Observation, Prediction and Response Network; Canadian Space Agency (FAST 18FAVICB09). Additional support from CFI/BCKDF and NSERC-DG to Costa for instruments on the BC Ferry. Hongyan Xi’s and Astrid Bracher’s contributions were funded by the ESA S5P+Innovation Theme 7 Ocean Color (S5POC) project (No 4000127533/19/I-NS) and by the Copernicus Marine Service GLOPHYTS project (21036L05B-COP-INNO SCI-9000). Copernicus Marine Service is implemented by Mercator Ocean in the framework of a delegation agreement with the European Union.

Acknowledgments

The authors would like to thank the Spectral lab members, especially Dr. Giannini, who organized and participated in the fieldwork. We acknowledge the BC Ferries Corporation for their permission for their crew members’ technical support during the fieldwork. The authors also thank ONC (Ocean Network Canada) for their technical help with FOCOS (Ferry Ocean Color Observation Systems) and ferrybox data acquisition. We thank Hakai Institute for providing data for the validation of PFTs retrieved from OLCI imagery. We also acknowledge Cermaq, Canada, for providing the opportunity to conduct fieldwork in their aquaculture farms, logistical support, and help with fieldwork. Finally, the authors thank Keith Holmes for his help with graphics.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.1018510/full#supplementary-material

References