Remote Sensing of Phytoplankton Size Class in Northwest Atlantic from 1998 to 2016: Bio-Optical Algorithms Comparison and Application

Abstract

Phytoplankton community structure and phytoplankton size class (PSC) are linked to ecological and biogeochemical changes in the oceanic environment. Many models developed to obtain the fraction of PSCs from satellite remote sensing have only been evaluated in open oceans, and very limited effort has been carried out to report on the performance of these PSC models in productive continental shelf waters. In this study, we evaluated the performance of nine PSC models in the coastal Northwest Atlantic (NWA) by comparison of in situ phytoplankton pigment measurements with coincidental satellite data from the Sea-Viewing Wide Field-of-View Sensor (SeaWiFS), Moderate-resolution Imaging Spectroradiometer (MODIS), and the Visible Infrared Imaging Radiometer Suite (VIIRS). Our results show that no PSC model retrieved all three phytoplankton size classes (pico-, nano-, and micro-phytoplankton) with reliable accuracy in the region of interest. In particular, these PSC models showed poor performance for retrieval of the picophytoplankton fraction of total phytoplankton in our study region, which could be related to the under-representation of pico-dominated samples in the productive waters of the NWA. For the accuracy of retrieved microphytoplankton and combined nano–pico phytoplankton fraction, the regional model developed by Devred et al. (2011) yielded the best result, followed by the model of Brewin et al. (2011). The model of Devred et al. (2011) was applied to satellite-derived chlorophyll-a concentration from the Ocean Color Climate Change Initiative (OC-CCI) archive in the NWA from 1998 to 2016. We report solely on the microphytoplankton biomass and fraction given the inverse relationship that exists with the nano–pico class. The multi-decadal trend along with the deseasonalized trend of microphytoplankton fraction was computed and analyzed for six biogeochemical provinces located in the NWA. Over the 19-year time series, there were significant, positive trends for four of the six provinces, with a slope of 0.36%·yr⁻¹ in the Northwest Continental Shelf (NWCS), 0.25%·yr⁻¹ in the Arctic Waters (ARCT), 0.12%·yr⁻¹ in the Slope Waters (SW) and 0.06%·yr⁻¹ in the Gulf Stream (GFST). Strong positive anomalies of microphytoplankton fraction were found in winter months in NWCS between 2009 and 2014, which could be associated with changes in environmental factors.

1. Introduction

Phytoplankton are the microscopic algae of size ranging from ~0.6 μm to individual and chain forming species greater than 200 μm [1]. Generally, this continuum of size is broken down into three size classes: pico- (0.2–2 μm), nano- (2–20 μm), and micro-phytoplankton (20–200 μm) [2]. Phytoplankton cell size has been considered a good indicator of functional role in many ecological and biogeochemical processes [2,3,4,5,6].

It is widely proposed that phytoplankton in different size classes occupy different physical and chemical niches based on their light-harvesting efficiency, nutrient-uptake ability, biogeochemical function, and mobility within the euphotic zone [1]. On a per volume or per biomass basis, smaller phytoplankton cells have higher metabolic rates [7], such that at constant intra-cellular pigment concentrations, the light-absorbing efficacy decreases with increasing phytoplankton cell size due to the package effect [8,9]. Conversely, large cells are associated with greater energy transfer efficiency to larger animals, and they are more likely to sink to deep waters where they are available to benthic species and contribute to the biological pump [1,10]. Therefore, the phytoplankton community size structure or phytoplankton size class (PSC) is linked to ecological and biogeochemical processes in the oceanic environment.

High-Performance Liquid Chromatography (HPLC) technique is currently the most extensively used non-microscopic method for characterizing in situ phytoplankton taxonomy because of its automatic operation, accuracy, and chemotaxonomy ability [11]. HPLC simultaneously provides a comprehensive description of the phytoplankton community composition over the entire cell size range and an accurate determination of the total Chlorophyll-a concentration (Chla) [12]. However, the HPLC-based method is not a direct measurement of phytoplankton size class, and it has several limitations, including few unambiguous marker pigments shared between several phytoplankton taxa [12] and existing assumptions about the ratios of Chla to diagnostic pigments [1,13,14]. Despite the limitations, we have used HPLC data in this study to derive the in situ PSC.

Information on PSC has been used in many large-scale marine biogeochemical models as a representation of ecological and biogeochemical functions [15,16,17,18]. In the past, validation and improvement of PSCs in biogeochemical models were based on comparisons with in situ measurements, but the elaboration of in situ databases to address spatial and temporal resolution requirements for models is time and labor consuming. In contrast, satellite-based remote sensing provides a synoptic view of optical properties in the upper ocean with high spatial and temporal resolution, making it an essential tool for the study of phytoplankton ecology. With advancements of optical remote sensing in recent decades, various PSC approaches have been developed, the most frequently used remote sensing PSC algorithms rely either on an abundance-based, spectral-based, or ecological-based approach [14,19,20,21]. Abundance-based approaches exploit the relationship between an indication of phytoplankton biomass (usually indicated by Chla or phytoplankton absorption coefficient) and cell size, following general ecological patterns wherein small phytoplankton size is generally associated with low phytoplankton biomass and vice versa [10]. The spectral-based approaches utilize optical properties of phytoplankton (or total particulate matter) that vary with phytoplankton size [1,9,22,23,24,25]. Based on the optical properties, the spectral-based approaches can be further partitioned into spectral absorption approaches and spectral backscattering approaches. For instance, spectral absorption approaches rely on the fact that small phytoplankton display higher specific absorption coefficients (i.e., absorption expressed in biomass unit) at blue wavelengths and steeper peaks when compared with larger phytoplankton [9,26,27]. These features can be exclusively tied to phytoplankton rather than all type of particles (phytoplankton, mineral particles, and detritus). Unlike the spectral absorption approach, the spectral backscattering approach provides size information on particulate assemblage rather than phytoplankton specifically, making them more suitable for application in the open ocean. Finally, ecological-based approaches retrieve the phytoplankton size using satellite-derived parameters of ecological indicators such as sea-surface temperature and mixed-layer depth [28,29]. Previous research has compared the performance of PSC algorithms at large oceanic scales; for example, Brewin, et al. [30] conducted an intercomparison between six satellite-derived PSC algorithms that were applied to Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) data and compared them with in situ observations collected in the global ocean. However, these in situ observations were taken mainly in the open oligotrophic ocean. Given that previous evaluations of PSC algorithms have focused on the open ocean [30,31,32], a comprehensive assessment of different PSC algorithms is needed for productive continental shelf waters, such as the Northwest Atlantic (NWA).

In this study, a total of 9 satellite-derived PSC algorithms were evaluated against in situ observations to assess their ability to detect dominant phytoplankton size classes (pico-, nano-, and micro-phytoplankton) from three sensors, namely the Sea-Viewing Wide Field-of-View Sensor (SeaWiFS), Moderate-resolution Imaging Spectroradiometer (MODIS), and Visible Infrared Imaging Radiometer Suite (VIIRS). The algorithms were selected based upon (a) their operability that has made them widely used by both experts and end-users; and (b) their different methodologies that provide information on the best approaches to retrieve PSC in the NWA. Based on an objective ranking, the best performing model was selected to study the variation of microphytoplankton between 1998 and 2016 using satellite data in six ecological provinces of the NWA as defined in Longhurst [33] and modified by Devred, et al. [34].

2. Materials and Methods

2.1. In Situ Measurements

Water samples used in this study were collected in the Northwest Atlantic region (39°N~62.5°N and −71°W~−42°W) between April and October from 1998 to 2016 as part of the Atlantic Zone Monitoring Program led by the Department of Fisheries and Oceans Canada [35]. The distribution of the sampling sites is shown in Figure 1. All the samples used in this study were collected in the surface layer of the ocean with a depth of 10 m or less. The total number of samples is listed in

Table 1. List of datasets used in this study.

Dataset	Sensor/Source	Quantity	Usage
Coincident Rrs and in situ pigments a (Figure 1a)	SeaWiFS	268	Validation of the satellite algorithms against in situ data Evaluation of the OC4 Chla products
	MODIS	644
	VIIRS	319
Coincident Rrs and in situ phytoplankton absorption (Figure 1b)	SeaWiFS	138	Evaluation of the IOP products
	MODIS	318
	VIIRS	130
Daily satellite image	OC-CCI v3.1	1998 to 2016	Generate the PSC in NWA

^a In this study, pigment includes: chlorophyll-a, fucoxanthin, peridinin, 19′-hexanoyloxyfucoxanthin, 19′-butanoyloxyfucoxanthin, alloxanthin, chlorophyll-b, and zeaxanthin.

. Note that from an initial database of 3117 HPLC samples, we only retained those that were coincidental (i.e., collected the same day) with satellite imagery from SeaWiFS, MODIS, and VIIRS.

Following collection at sea, water samples were filtered through GF/F glass fiber filters, which were then frozen in liquid nitrogen (−80 °C) for later analysis in the laboratory. Total phytoplankton absorption coefficients between 400 and 750 nm were measured using a Shimadzu UV-2101 spectrophotometer following the method of Yentsch [36] as modified by Mitchell and Kiefer [37], Hoepffner and Sathyendranath [38] and Kyewalyanga, et al. [39]. Absorption by detrital materials was measured after extraction of the pigments using hot methanol. Phytoplankton absorption was then obtained by subtracting the detrital absorption from the total phytoplankton absorption. Concentrations of phytoplankton pigments were measured using HPLC following Stuart and Head [40]. Pigments were extracted from the frozen filters in 90% acetone and were identified according to their retention times and by comparison with those pigment’s standards. The measured pigments include: chlorophyll-a, fucoxanthin, peridinin, 19′-hexanoyloxyfucoxanthin, 19′-butanoyloxyfucoxanthin, alloxanthin, chlorophyll-b, and zeaxanthin.

2.2. Satellite Data

2.2.1. Satellite Dataset

●  Satellite Level 2 ocean-color data

L2 daily products from SeaWiFS (1999 to 2010), MODIS (2002 to 2016) and VIIRS (2012 to 2016) sensors were downloaded from the National Aeronautics and Space Administration (NASA) ocean-color website (https://oceancolor.gsfc.nasa.gov/) and matched to in situ HPLC data using a window of 3 × 3 pixels centered on the location of the sample at the nominal spatial resolution of each sensor (1.1, 1.0, and 0.75 km for SeaWiFS, MODIS, and VIIRS respectively). The associated standard deviation of remote sensing reflectance (R_rs) of each matchup was calculated across the 9 pixels and used to account for uncertainties in the PSC-models (see Section 2.3.2). Chla, which was calculated following the standard OC_x band ratio algorithm (https://oceancolor.gsfc.nasa.gov/atbd/chlor_a/#sec_2), were also downloaded from the NASA website. Here, the x subscript denotes the sensor specific algorithm to derive Chla, which accounts for waveband shifts between sensors (see https://oceancolor.gsfc.nasa.gov/atbd/chlor_a/ for details). In addition, the MODIS SST L3 products at 9 km resolution were obtained from the NASA website (https://oceancolor.gsfc.nasa.gov/cgi/l3) as input of the ecological-based model I (

Table 2. Description of models used in this study.

Model	Reference	Type	Size Classes	Satellite Input Variables
				Chla	aph(λ) 1	aph(λ) 2	SST
A	Hirata et al. (2008)	Abundance-based	3		√
B	Hirata et al. (2008)	Abundance-based	3			√
C	Hirata et al. (2008)	Abundance-based	3	√
D	Devred et al. (2006)	Abundance-based	2	√
E	Devred et al. (2011)	Abundance-based	3	√
F	Devred et al. (2011)	Spectral-based	3		√
G	Devred et al. (2011)	Spectral-based	3			√
H	Brewin et al. (2010)	Abundance-based	3	√
I	Barnes et al. (2011)	Ecological-based	3	√			√

¹ Refers to a_ph(λ) calculated by reparameterization in NWA following the procedure described in Craig et al. (2012) [41]; this IOP model is termed as EOF_NWA thereafter. ² refers to a_ph(λ) calculated according to Lee et al. (2002) [42]; QAA v5 was selected in order to be consistent with Brewin et al. (2011) [30].

).

●  Matchups datasets

Matchups are defined as in situ samples that were collected on the same day of the satellite pass. We did not follow the stringent criteria used by NASA for satellite validation (i.e., collection of in situ data within 3 h of satellite pass in open waters) given that the current study is not a validation exercise of satellite Chla but rather an evaluation of the performance of PSC algorithms for operational applications. Information on the remote sensing retrieval of Chla, however, will provide a useful insight in the understanding of PSC-model performance. The total number of matchups between SeaWiFS, MODIS, VIIRS R_rs spectra (or Chla) and in situ HPLC are 268, 644, and 319, respectively (Figure 1a and Table 1). Although the samples were collected between the months of April and October, the samples distribution is biased toward the months of April, September, and October, corresponding to the Atlantic Zone Monitoring Program cruises (Figure 2). Chla concentration ranged from 0.09 to 13.21 mg·m⁻³ for SeaWIFS, and 0.07 to 13.94 mg·m⁻³ for both MODIS and VIIRS (Figure 3). The frequency distribution of Chla for the global SeaBASS archive (NOMAD, 0.02 to 70.21 mg·m⁻³) is also shown in Figure 3 for reference.

In addition to R_rs and Chla concentration, two inherent optical properties (IOP) models were used to derive the phytoplankton absorption coefficient (a_ph(λ)) from spectral R_rs. These two IOP models rely on either advanced statistical methods (principal component analysis, PCA, as in Craig, et al. [41]) or on analytical approximation of the radiative transfer equation (quasi-analytical algorithm, QAA, [42]). The coefficients of the advanced statistical method were recomputed based on an empirical orthogonal function (EOF) analysis using the MATLAB function princomp following the procedure described in Craig, et al. [41] for the three sets of matchups (i.e., SeaWiFS, MODIS, and VIIRS). This IOP model is then termed as EOF_NWA thereafter. QAA version 5 was employed to remain consistent with Brewin, et al. [30], so this IOP model is termed as QAA_v5 thereafter. The rationale for using two different approaches, a global and a regional, was to assess the impact of the absorption algorithms on the performance of PSC methods.

Similar to HPLC-derived Chla, satellite data coincident with in situ a_ph(λ) were extracted from daily images to assess the accuracy of the a_ph(λ) products (λ = 443, 490, and 555 nm). The number of matchups for SeaWiFS, MODIS, VIIRS were 138, 318, and 130, respectively (Figure 1b and Table 1).

2.2.2. Biogeochemical Provinces and Climate Change Compatible Time Series of Ocean Color in the NWA

The best performing algorithm was applied to the Ocean Color Climate Change Initiative (OC-CCI) v3.1 daily Chla product, which spans from 1998 to 2016 and used merged data from the Medium Resolution Imaging Spectrometer (MERIS), Aqua-MODIS, SeaWiFS, and VIIRS sensors. Daily R_rs and Chla concentration were downloaded from the OC-CCI website (https://www.oceancolour.org) and were extracted for the region of interest (see Section 2.1). This dataset has the advantage of being quality controlled and inter-sensor calibrated, which ensure climate-compatible and continuous ocean-color observations at the global scale [43]. Information on Chla was converted to percentage of microphytoplankton using Devred, et al. [25].

In order to account for different ecological regimes and detect decadal trends of PSC in the NWA, the region was divided into six ecological provinces based on physical and biological characteristics [33] (Figure 1b), following the method of Devred, et al. [34]. The six ecological provinces are the Gulf Stream (GFST), Slope Waters (SW), North Atlantic Drift (NADR), Northwest Continental Shelf (NWCS), Arctic Waters (ARCT), and Polar Boreal Current (BPLR) (Figure 1b). In the present study, we used the climatology of the provinces boundaries that was generated using the time series of MODIS Chla and SST data from 2003 to 2015.

The GFST is an oligotrophic region with high salinity and high temperature relative to the other five regions; it is subject to strong mesoscale variability. Two to ten rings a year, that break off the main stream and move towards the slope water may be observed from satellite products of thermal, surface elevation or Chla concentration [33]. The SW is a province that was added to the original scheme of Longhurst [33] by Devred, et al. [34], this province represents a surface water mass located on the continental slope between the high-primary-production waters (on the continental shelf) and the poor-primary-production waters (oligotrophic water of the Gulf Stream and subtropical gyre). This is an intermediate water mass with mean temperature close to the GFST province but mean Chla concentration close to the NWCS province (Devred, et al. [34]). The NADR corresponds to a mixed of two contrasted water masses, the cold nutrient-rich waters from the north (BPLR and ARCT) and warmer oligotorophic waters from the south (GFST) in that respect, this province shows the highest seasonal variability in terms of temperature and Chla concentration. Note that our definition of the Northwest Atlantic encompasses only a small fraction of the NADR, which extends across most of the north Atlantic Basin. The NWCS spans from the Gulf of Maine (South) to the northern tip of Newfoundland. It is a productive water mass with intense mixing given the shallow waters (< 200m) and the strong atmospheric forcing. It is influenced by the BPLR and SW waters. The NADR comprises two cyclonic Subpolar Gyres of the North Atlantic Ocean, within which the water is relatively cold and saline [44]. The BPLR on the western side of the region is a cold arctic water mass that is subject to seasonal ice cover [33]. On the eastern side, it is influenced by cold Atlantic water that flow northwards along the Greenland shelf. Both BPLR and ARCT cold provinces (<10 °C) that sustain a phytoplankton bloom that occurs late in the year compared to the southern provinces.

Due to the extensive and frequent coverage of sea-ice and cloud in winter, low numbers of valid pixels are obtained in these months in the high-latitude ecological provinces (e.g., NADR, ARCT, and BPLR). To avoid large biases due to the lack of data in a given month in a given province, monthly-composites of PSC with more than 40% missing pixels in that province were excluded during the time-series analysis. Thus, for NADR, ARCT, and BPLR, only the months of March to October, March to October, and May to October were included in the time-series and anomaly analysis, respectively.

For each province, monthly time series of PSCs were calculated by averaging all valid pixels within each month for each individual year. The time-series of monthly PSC in each ecological province were deseasonalized by applying “unstable” seasonal filters (i.e., using a multiplicative decomposition, provided by MATLAB, https://www.mathworks.com/help/econ/seasonal-adjustment-using-snxd7m-seasonal-filters.html). After removing the seasonal component from the time-series, a ‘LOESS’ filter (Local regression using weighted linear least squares and a second-degree polynomial model, with 25% span) was applied to the deseasonalized time-series trend signal to characterize the ‘long-term trend’.

2.3. PSC Algorithms and Ranking Method

2.3.1. PSC Algorithms

Four abundance-based (model A, B, C, D, and E), two spectral-based (model F and G) and one ecological-based PSC algorithm (model H) were selected to carry out a performance assessment (Table 2). Given that two different IOP models were used to derive phytoplankton absorption (i.e., EOF_NWA and QAA_v5), we tested a total of 9 models. A brief description of each model is given below.

●  Model A B C: Hirata et al. (2008) [45] Abundance-based method

Hirata, et al. [45] is a threshold-based method that uses a single variable, a_ph(443) or Chla, to quantify the dominant phytoplankton size classes. Models A and B use a_ph(443) as input, and the three populations of phytoplankton are classified as follows: pico, a_ph(443) < 0.023 m⁻¹; nano, 0.023 ≤ a_ph(443) < 0.069 m⁻¹; micro, a_ph(443) ≥ 0.069 m⁻¹. Model A uses satellite-derived a_ph(443) according to Craig, et al. [41] (EOF_NWA) and model B uses satellite-derived a_ph(443) according to Lee, et al. [42] (QAA_v5). Model C uses satellite-derived OC_x Chla product as an input. The thresholds of Chla used in model C are: pico, Chla < 0.25 mg·m⁻³; nano, 0.25 ≤ Chla < 1.8 mg·m⁻³; micro, Chla ≥ 1.8 mg·m⁻³.

●  Model D: Devred et al. (2006) [26] Abundance-based method

Devred, et al. [26] developed a two-population model based on Sathyendranath, et al. [19]. This model relates the absorption coefficient of phytoplankton to Chla concentration assuming the presence of two populations of phytoplankton in a seawater sample with different absorption properties related to the mean size of each population. The two populations of phytoplankton classified in this study are assumed to be (a) microphytoplankton, and (b) mixed nano–pico phytoplankton. The core formulations of this model are:

$$\mathrm{Chl}{a}_{nano-pico}=\mathrm{Chl}{a}_{nano-pico}^{max}·\left(1-e^{-S·\mathrm{Chl}a}\right)$$

(1)

and

$$\mathrm{Chl}{a}_{micro}=\mathrm{Chl}a-\mathrm{Chl}{a}_{nano-pico}$$

(2)

where Chla is the total chlorophyll-a concentration, $\mathrm{Chl}{a}_{micro}$ is the Chla of microphytoplankton, $\mathrm{Chl}{a}_{nano-pico}$ is the Chla of mixed nano–pico phytoplankton, $\mathrm{Chl}{a}_{nano-pico}^{max}$ is the maximum concentration of mixed nano–pico phytoplankton, S is the increasing slope of $\mathrm{Chl}{a}_{nano-pico}$ with Chla.

This model can be applied in two ways: (a) Either by representing the satellite-derived a_ph as a linear combination of the specific absorption of the two populations; (b) or by using it as an abundance-based approach by using Equations (1) and (2). In this study, we used it as an abundance-based approach (e.g., Equations (1) and (2)) by applying the parameters derived from the total dataset provided in Devred, et al. [26].

Since model D is a two-population model, it was excluded in the comparison of nanophytoplankton and picophytoplankton models to in situ measurements.

●  Model E: Devred et al. (2011) [25] Abundance-based method

In this model, the three populations of phytoplankton are obtained by successive application of the two-population model D (i.e., first classify the total phytoplankton as micro- and mixed nano–pico phytoplankton, and then further classify the mixed nano–pico phytoplankton as nano- and pico- phytoplankton). As model D, this model can be applied either using a spectral- or abundance-based approach. Here, Model E was used as an abundance-based method, the parameters simulated from the NWA and NOMAD (as shown in [25]) were adopted in this comparison.

●  Model F and G: Devred et al. (2011) [25] Spectral-based method

Model F and G used Devred, et al. [25] as a spectral-based method, with satellite-derived a_ph as inputs. The a_ph used in model F and G were calculated from R_rs using a method based on PCA analysis from Craig, et al. [41] (EOF_NWA) and following the scheme of [42] (QAA_v5), respectively.

●  Model H: Brewin et al. (2010) [46] Abundance-based method

This method estimates the Chla of the pico-, nano-, and micro-phytoplankton as a function of total Chla. The ecological basis of this method is similar to that of Devred, et al. [26] and Devred, et al. [25], based on observations that small-sized phytoplankton tends to dominate in oligotrophic waters and large-sized phytoplankton tends to dominate in waters where Chla concentrations are high. The main difference between the Devred, et al. [25] and Brewin, et al. [46] approaches is that the coefficients used in the Brewin, et al. [46] formulation were derived using HPLC data, whereas the parameters in Devred, et al. [25] were derived by fitting a model on absorption data against Chla. The coefficients of the parameterization in Brewin, et al. [46] computed for the Atlantic Ocean were used in the present study.

●  Model I: Barnes et al. (2011) [28] Ecological-based method

In this model, the size composition of the phytoplankton communities was predicted using two remotely-sensed variables: Chla and sea-surface temperature (SST). This model described the relationship between the environment and the size composition of phytoplankton communities, using a collation of empirical measurements of size composition from sites that include polar, tropical and upwelling environments. Empirical relationships were established between phytoplankton biomass (calculated as the product of numerical abundance and cell mass) and satellite-derived Chla and SST. Using these relationships, the biomass in each group (i.e., pico-, nano-, and micro-phytoplankton) can be determined by integrating the biomass within the corresponding cell mass ranges.

2.3.2. Algorithm Assessment Method

Size class fraction of all HPLC samples were derived using diagnostic pigment analysis (DPA) according to the method proposed by Vidussi, et al. [13], extended by Uitz, et al. [14] and further adjusted by Brewin, et al. [46] and Devred, et al. [25]. The ranges of Chla for micro-, nano- and pico-phytoplankton derived by this method are 0~13.53, 0~2.74 and 0~1.05 mg·m⁻³, respectively.

We also tested the HPLC-based classification method of Hirata, et al. [45], which gave similar results to the Uitz, et al. [14] method (result not shown). The dominant size class was established based on whether a size class (pico-, nano- or micro-phytoplankton) had a diagnostic pigment to total Chla ratio of greater than 0.45.

In order to make the nine models comparable, a scoring method similar to Hirata, et al. [45] and Brewin, et al. [30] was used to evaluate the performance of the satellite-derived PSCs when compared with in situ matchups (Figure 4). A “correct”, “near-correct”, and “incorrect” classification result was scored as 2, 1, and 0 points, respectively.

First, the standard deviation (SD) of Chla, a_ph, and R_rs for the 9 pixels (3 × 3 window) of the satellite input was calculated for each matchup. The 9 PSC models were then run on the mean, mean plus SD and mean minus SD of the satellite-derived property of interest (i.e., either Chla, a_ph or R_rs) giving three, or two in the case of Model D, dominant classes for each model. If any of the three satellite-derived dominant classes matched the dominant size class calculated from the in situ pigments, a “correct” classification of 2 points was assigned. There were rare situations where two dominant size classes were identified in the in situ dataset, e.g., both nano- and micro-phytoplankton had a diagnostic pigment to total Chla ratio of over 0.45. In this case, if one of the three satellite-derived dominant size classes matched the dominant in situ size class, a correct classification of two points was assigned.

If the three satellite-derived size classes mismatched with the HPLC size-class assignment, the results were further assessed to test if a second dominant size class in the in situ sample had a DPA ratio equal to or greater than 0.4 (based on an uncertainty estimate of 9.3% for the DPA ratio; [47]), these cases were labeled as a co-dominance situation. If any of the three results matched the second dominant size class, a “near-correct” classification of 1 point was assigned. For matchups where the three satellite-derived size classes mismatched with the in situ classification, no points were assigned to the model. For a given model, the percentage accuracy (%) was calculated by dividing the number of points accumulated using this model by the maximum possible number of points and multiplying by 100.

Note that in some cases, the diagnostic pigment to total Chla ratio did not return a dominant size class (e.g., pico 33.3%, nano 33.3%, micro 33.4%) or returned two dominant size classes in one sample (e.g., pico 46%, nano 48%, micro 6%). This explains why the sum of micro-, nano-, and pico-phytoplankton dominated samples was not necessarily equal to the total number of samples.

In this study, models A, B, and C are designed to obtain the dominant phytoplankton size class, while models D, E, F, G, H, and I derive fractional contributions. This evaluation method was developed to assess the dominance, so those approaches that derive fractional contributions (models D, E, F, G, H, and I) have been adapted in order to make the models comparable. Therefore, these methodological differences between models should be accounted for when discussing the performance of the models.

2.4. Accuracy Assessment

The accuracy of satellite derived Chla and a_ph were assessed by calculating the percentage error (PE), mean absolute percentage deviation (MAPD) and mean percentage deviation (MPD) between measured and derived values, using the following equations:

$$\mathrm{PE}=100\%·\frac{x_d-x_m}{x_m}$$

(3)

$$\mathrm{MAPD}=100\%·\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\left|\frac{x_d-x_m}{x_m}\right|$$

(4)

$$\mathrm{MPD}=100\%·\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\frac{x_d-x_m}{x_m}$$

(5)

where x_m is the measured value, x_d is the derived value and n is the total number of samples.

3. Results

3.1. PSC Algorithm Comparison Results

3.1.1. SeaWiFS

The nine models were applied to SeaWIFS, MODIS, and VIIRS datasets, and a comparison of the percentage accuracy was carried out (

Table 3. Comparison results showing the percentage accuracy of different models with SeaWiFS, MODIS, and VIIRS products as model inputs. The comparison results of models B, C, and D in the Open Ocean conducted by Brewin, et al. [30] are also presented.

	SeaWiFS					MODIS				VIIRS
Dataset	Method	Micro %	Nano + Pico %	Nano %	Pico %	Micro %	Nano + Pico %	Nano %	Pico %	Micro %	Nano + Pico %	Nano %	Pico %
Northwest Atlantic	A	24.2 ± 6.5	75.3 ± 7.4	94.8 ± 5.5	13.0 ± 12.9	21.2 ± 4.5	56.2 ± 4.8	75.1± 5.4	52.1 ± 8.2	40.8 ± 8.8	48 ± 6.2	69.3 ± 8.5	50 ± 8.6
	B	17.5 ± 5.8	24.7 ± 7.3	46.0 ± 9.3	69.6 ± 19.4	1.0 ± 1.1	0.6 ± 0.7	0	90.7 ± 5.5	1 ± 1.8	0.42 ± 0.8	0	97.3 ± 4.1
	C	47.5 ± 7.7	69.3 ± 7.9	95.4 ± 5.3	52.2 ± 18.4	40.2 ± 5.5	71.4 ± 4.4	94.8 ± 3.5	24.6 ± 6.8	19.4 ± 7.0	76.8 ± 5.3	95.0 ± 5.0	30.4 ± 7.7
	D	80.6 ± 6.2	88.0 ± 5.7	\	\	74.0 ± 5.0	83.1 ±3.7	\	\	60.0 ± 8.8	89.4 ± 3.9	\	\
	E	83.1 ± 5.9	81.7 ± 6.7	0	52.2 ± 18.4	79.5 ± 4.6	78.0 ± 4.1	0	28.0 ± 7.0	70.1 ± 8.4	86.0 ± 4.4	0	30.0 ± 7.7
	F	86.6 ± 5.3	16.7 ± 6.4	21.3 ± 7.7	0	87.9 ± 3.8	40.0 ± 4.7	43.6 ± 6.1	0	89.1 ± 6.3	64.3 ± 6.0	69.8 ± 8.5	0
	G	55.7 ± 7.6	91.2 ± 5.0	3.4 ± 2.9	78.3 ± 18.3	52.8 ± 5.6	85.4 ± 3.5	0	57.2 ± 7.8	59.2 ± 9.0	96.2 ± 2.5	0	65.6 ± 8.3
	H	63.7 ± 7.5	94.4 ± 4.2	92.0 ± 6.2	30.4 ± 17.3	53.0 ± 5.6	94.4 ± 2.4	89.3 ± 4.2	11.4 ± 4.9	33.3 ± 8.2	99.2 ± 1.4	93.0 ± 5.5	13.0 ± 5.6
	I	0	100 ± 1.7	100 ± 3.9	47.8 ± 18.8	0	100.0 ± 0.9	100 ± 2.5	2.1 ± 1.9	0	100 ± 0.9	96.0 ± 4.7	21.4 ± 6.9
	No. of samples	165	121	81	14	305	403	187	133	110	240	87	122
Open Ocean (Brewin et al. 2011)	B	36.5 ± 7.6	96 ± 1.8	22.1 ±5.4	93.1 ±3.2	\	\	\	\	\	\	\	\
	C	90.1 ± 4.7	95.1 ± 1.9	37.1 ± 6.4	87.4 ±4.1	\	\	\	\	\	\	\	\
	D	91 ± 4.5	94 ± 2.1	\	\	\	\	\	\	\	\	\	\
	No. of samples	92	285	112	173	\	\	\	\	\	\	\	\

and Figure 5). Among the models applied to SeaWIFS data, D, E, and F showed relatively high accuracy (>80%) for the microphytoplankton size class, followed by G and H (55.7~63.7%). The performance of models A, B, C, and I was found to be relatively low (0~47.5%).

None of the models performed well for both nanophytoplankton and picophytoplankton (details below), and therefore a comparison of model performance for mixed nano–pico phytoplankton was made by treating the above models as two-population models (microphytoplankton and mixed nano–pico phytoplankton). Models D, E, G, H, and I were found to perform with higher accuracy (81.7~100%) than models A, B, C, and F (16.7~75.3%) although models A and C were not greatly different from model E.

For the nanophytoplankton size class, the performing accuracy of models A, C, H, and I (92.0~100%) were much higher than those of models B, E, F, and G (0~46.0%). For the picophytoplankton size class, the performing accuracy for all of the models was low, among which models B, C, E, G, and I performed with relatively higher accuracy (52.2~78.3%) than the rest of the models (13.0~47.8%).

The performance of models B, C, and D was evaluated using a global dataset in Brewin, et al. [30] (Table 3). The performance of the models was higher with the global dataset than with the dataset used in the present study, except for identification and estimation of the nanophytoplankton class. This could be explained by the fact that our in situ dataset is biased towards high Chla concentration (mesotrophic to eutrophic waters, Figure 3).

3.1.2. MODIS

For microphytoplankton-dominated populations, the accuracy of the nine algorithms applied to MODIS data was generally the same as found when using SeaWiFS: Models D, E, and F exhibited higher accuracy (74.0~87.9%) than models G and H (52.8~53.0%), followed by models A, B, C, and I (0~40.2%). Regarding the mixed nano–pico phytoplankton size classes, models C, D, E, G, H, and I performed with higher accuracy (71.4~100%) than models A, B, and F (0.6~56.2%).

For nanophytoplankton dominated populations, the accuracy of models A, C, H, and I were relatively high (75.1~100%) among the nine models, followed by model F (43.6%). The accuracy of models B, E, and G at retrieving dominance of nanophytoplankton was null. Regarding picophytoplankton dominated populations, model B exhibited the highest accuracy among all of the models (90.7%), followed by models G (57.2%) and A (52.1%), models C, E, H, and I exhibited relatively low accuracies (0~28%).

3.1.3. VIIRS

Among the models applied to VIIRS data in microphytoplankton dominated samples, models D, E, F, and G (59.2~89.1%) exhibited higher accuracy than models A, B, C, H, and I (0~40.8%). For mixed nano–pico phytoplankton populations, models C, D, E, G, H, and I (76.8~100%) performed with higher accuracy than models A, B, and F (0.42~64.3%). For nanophytoplankton dominated populations, the accuracy of models C, H, and I were highest (93~96%) among the nine models, followed by model A and F (69.3% and 69.8%, respectively). The accuracy of models B, E, and G were all found to be null. For picophytoplankton dominated samples, model B performed with the highest accuracy (97.3%), followed by models A and F (50.0~65.6%), and then models C, E, G, H, and I (0~30.4%).

Taking all of the satellite sensors and three phytoplankton size classes into consideration, none of the models performed with reasonable accuracies (e.g., >60%) at retrieving all three of the micro- nano-, and pico- phytoplankton size classes together. In other words, retrieving one size class accurately is achieved at the expenses of accurately retrieving another class(es). If we consider SeaWiFS and two populations (i.e., microphytoplankton and mixed nano–pico phytoplankton), models D, E, and H achieve greater than 60% accuracy for both size classes in the NWA. When the three sensors and two-population models were accounted for, model E was the only one that performed with a greater than 60% accuracy for both size classes.

4. Discussion

4.1. Uncertainties Associated with Satellite-Derived Inputs

4.1.1. OC_x Chla Products

Models C, D, E, H, and I are all based on Chla, and in the case of model I, also SST, and therefore the accuracy of satellite Chla products is critical for phytoplankton size classification. The dominant PSC of Model C is based on thresholds of Chla. For models D, E, and H, the variation of each phytoplankton size class was defined as a non-linear function of Chla. However, for a given model and set of parameters, dominance of one group over others is always set by a given Chla (Figure 6), and thus dominance of a given phytoplankton class in these models is assessed similarly to a threshold algorithm such as model C.

For model D, samples with Chla less than 1 mg·m⁻³ were classified as mixed nano–pico phytoplankton dominated; whereas those with Chla greater than 1 mg·m⁻³ were classified as dominated by microphytoplankton (Figure 6a). For model E with the parameters used here, the maximum fraction of the nanophytoplankton size class that can be achieved is less than 45% (~40%, Figure 6b). This means that the criterion of dominant size class defined as greater than 45% is never realized and therefore model E was unable to retrieve nanophytoplankton as a dominant size class. In model E, samples with Chla values less than 0.3 mg·mg⁻³ were classified as picophytoplankton dominated, while those with Chla values greater than 1 mg·m⁻³ were classified as microphytoplankton dominated. Model H (Figure 6 c) classified the dominant PSC was classified as picophytoplankton when Chla was lower 0.15 mg·m⁻³, nanophytoplankton when Chla ranged between 0.15 and 1.5 mg·m⁻³, and microphytoplnakton when Chla > 1.5 mg·m⁻³.

In order to evaluate the impact of Chla on the estimation of PSC, a sensitivity analysis was also conducted (Figure 6). Here, we assumed that the Chla varies within a range of ±30%, which is the targeted accuracy for satellite-derived Chla in open oceans. When Chla had a positive bias, the microphytoplankton was overestimated in all three models (models D, E, and H), but the mixed nano–pico phytoplankton (model D) and picophytoplankton (models E and H) were underestimated. It is worth pointing out that when Chla had a positive bias, the nanophytoplankton was overestimated for the low Chla (<~0.5 mg·m⁻³) but underestimated for high Chla (~0.5 mg·m⁻³) (Figure 6b,c). The opposite was true when the Chla had a negative bias.

The error of satellite-derived Chla products is biased such that Chla is mainly overestimated when in situ Chla is low and underestimated when it is high (Figure 7). As expected, similar patterns of bias were observed for all three sensors. SeaWiFS OC_x Chla was generally overestimated when in situ Chla concentration was less than 0.4 mg·m⁻³, while SeaWiFS OC_x Chla was underestimated when in situ Chla was over 1 mg·m⁻³. For MODIS OC_x Chla, the PE was generally positive when in situ Chla was less than 0.3 mg·m⁻³, and a systematic underestimation occurred when in situ Chla was greater than 2 mg·m⁻³. For the VIIRS OC_x product, a systematic overestimation occurred when in situ Chla was less than 0.3 mg·m⁻³, and a systematic underestimation was observed when in situ Chla was greater than 1 mg·m⁻³.

Both the relationships between PSC and Chla and the uncertainties associated with satellite-derived Chla explain the poor performance of the abundance-based models to infer nano- and pico-phytoplankton on the NWA continental shelf (Figure 6 and Figure 7). On one hand, given the systematic overestimation of Chla by the OC_x algorithm for Chla below 0.3 mg·m⁻³, in situ samples identified as picophytoplankton dominant can be misclassified as nanophytoplankton or even microphytoplankton dominated (Figure 6 and Figure 7). On the other hand, despite the underestimation of satellite-derived Chla using the OC_x algorithm for in situ Chla above 2 mg·m⁻³, the outcome of the abundance-based models remained accurate for the microphytoplankton dominated size class, even considering the vast range of Chla (>1 mg·m⁻³, Figure 6) for the microphytoplankton dominated region. The in situ samples identified as microphytoplankton dominant but vulnerable to be misclassificated as nano- or pico- phytoplankton dominated (in situ Chla >2 mg·m⁻³ and satellite Chla <1 mg·m⁻³) were identified using blue bars in Figure 7.

4.1.2. Phytoplankton Absorption Coefficient Products

Models A, B, F, and G were based on absorption coefficient products. Two IOP models were selected to calculate a_ph from R_rs, EOF_NWA and QAA_v5. The PE of satellite-derived a_ph at 443, 490 and 555 nm compared to in situ data was further evaluated to help understand the impact of absorption retrieval on the performance of the algorithms. For all three satellite sensors (i.e., SeaWiFS, MODIS, and MERIS), the EOF-NWA model provided better statistical metrics for a_ph retrieval than the QAA model (Figure 8), specifically, a narrower distribution of PE and lower values of MAPD ranging from ~35% to ~47%. The a_ph derived from QAA_v5 shows a relatively wide distribution of PE, with most of MAPD higher than 76%. In particular, the EOF_NWA model showed positive MPD for all of three bands of the three satellite sensors (9.4~25.15%), demonstrating a slight overestimation with this model (Figure 8). For the QAA_v5 model, systematic negative MPD values were found for all three wavebands for all three satellite sensors (−29.6~−78.7%), with the exception of 490 and 555 nm wavebands of SeaWiFS. This demonstrated an overall underestimation on a_ph with the QAA_v5 model (Figure 8).

The original version of QAA was used in our analyses and no regional parameterization was performed, which may explain the relatively poor performance of this method when tested on a regional dataset. Similar to our findings, Shanmugam, et al. [48] showed that phytoplankton absorption coefficients retrieved by QAA were generally underestimated for the selected wavebands (412, 443, 490, 510, and 555 nm) in coastal regions. The underestimation revealed by Shanmugam, et al. [48] (as well as in our study) was attributed to the spectral slope of colored dissolved organic matter absorption (a_CDOM) and detritus (a_d), which resulted in an overestimation of a_CDOM and a_d, and corresponding underestimation of a_ph at these wavelengths. In agreement with previous works [25,26], our results suggest that a regional parameterization of IOP models by a fully suited dataset is required in our region of study.

4.2. Model Selection for Monitoring PSC in the NWA

Our results showed that dominance of microphytoplankton and mixed nano–pico phytoplankton populations (i.e., two-population model) in the NWA were more accurately detected than three-population classifications (Figure 5). This might be due to (i) the abundance-based models’ reliance on low Chla values to detect nano- and picophytoplankton and the OC_x algorithms’ systematic overestimation of Chla at low concentration (i.e., Chla < ~0.3/0.4 mg·m⁻³) (Figures 7), and (ii) productive waters in the study areas, in which oligotrophic conditions are under-represented in the in situ dataset (Figure 3).

Among the three different approaches, abundance-based and spectral-based models performed better than the ecological-based model (model I), which tended to over classify pixels in the nanophytoplankton size class in the region of interest. The ecological-based model, based on phytoplankton biomass, was developed to retrieve PSCs over large spatial scales, limiting the application of these large-scale empirical ecological relationships, which were rather scattered with low determination coefficients, as is shown in Barnes, et al. [28], to the relatively small shelf region.

The performance of the absorption-based models (e.g., A, B, F, and G) depends on the IOP model selected. Models A and F, which used the EOF_NWA approach to derive phytoplankton absorption, performed better at detecting microphytoplankton and nanophytoplankton dominated pixels than the picophytoplankton dominated pixels. The derived phytoplankton absorption by EOF_NWA is, on average, overestimated, which results in an underestimation of picophytoplankton. Models B and G that had inputs calculated according to QAA_v5 [42], worked better at detecting picophytoplankton dominated pixels than those dominated by the other two size classes. This is a result of the systematic underestimation of a_ph by QAA_v5, compared to estimates by EOF_NWA (Figure 8).

Compared to open ocean, the 3-class PSC models showed relatively poor performance in coastal waters (Scotian Shelf, Northwest Atlantic), in particular for the picophytoplankton class ([30] and Table 3 in this study). This phenomenon is mainly a result of: (a) The larger uncertainty of the Chla product in productive waters compared to the open ocean, and the possible contamination of yellow substance absorption to the total water leaving signal, known to artificially increase Chla in band-ratio algorithms; and (b) under-representation of samples dominated by picophytoplankton in coastal regions (Figure 3) since they generally dominate in oligotrophic waters. Therefore, the application of the models to productive waters must be handled cautiously and would benefit from regional adjustment.

Regarding the use of HPLC-based PSC as a reference, several limitations exist that have to be accounted for when interpreting the results. Some marker pigments are shared among several phytoplankton taxa, making it difficult to identify different groups unambiguously. For example, Chlb is the marker pigment for green flagellates and prochlorophytes [49,50,51,52]. However, it is also present in Euglenophyta, which is classified as nano- or micro-phytoplankton [12]. In addition, the HPLC-based phytoplankton size class determination is sensitive to assumptions about the Chla to diagnostic pigment ratios [1,13,14]. Therefore, uncertainties of the in situ data are not negligible. None of the current in situ methods to identify phytoplankton taxa or functions (such as microscopy, flow cytometry, and filtration etc.) are considered ideal for detection and quantification of PSC [1]. A compromise has to be achieved between the accuracy at retrieving information of phytoplankton against the range covered by the method used. For instance, while microscopic counts and flow cytometry have high accuracy of measurements, these methods are constrained by the phytoplankton cell size they can resolve (~5 to 300 $\mathsf{\mu }$m and 0.2 to ~20 $\mathsf{\mu }$m for light microscopy and flow cytometry respectively). Conversely, HPLC-derived pigment composition and concentration resolve the entire size range of the bulk phytoplankton population but have less capacity to identifying phytoplankton groups accurately. Here we have relied on HPLC pigment information, one of the most widely used measurements in the field of marine optics/biogeochemistry. It is also noteworthy that PSC models that were developed using HPLC-pigment database are the ones that provide the best results, which advocate for an independent measurement to truly assess the performance of these algorithms. In addition, it is expected that models developed using regional datasets will perform better than those developed using global datasets. However, some global models, such as Brewin et al. (2010) showed that they could still be applied on the NWA continental shelf with meaningful results.

4.3. PSC in Northwest Atlantic

Our results revealed that none of the PSC models tested were sufficiently accurate in retrieving the three PSCs simultaneously (micro-, nano-, and pico- phytoplankton). It appears that for any given algorithm, successfully retrieving information on a particular size class is achieved at the expense of other classes, meaning that an algorithm that performs well to retrieving the fraction of microphytoplankton will do so at the expense of retrieving nano- or pico-phytoplankton fractions accurately. When we applied the stringent criteria of at least 60% accuracy in retrieving two size class fractions (microphytoplankton and mixed nano–pico phytoplankton) for all three sensors, only model E performed successfully. The above average performance of model E is not a surprise since: (a) simpler models tends to be more robust; and (b) it was developed using a dataset from this NWA region (spanning 1996–2003). Given the inverse relationship between the microphytoplankton and mixed nano–pico phytoplankton classes (see Equation (2)), only the time series analysis of CCI-derived microphytoplankton fraction percentage is reported here.

4.3.1. Temporal Variation of Microphytoplankton in NWA

Time-series of microphytoplankton fractions in the six ecological provinces were decomposed into three components (seasonal cycle, trend, and residual) to detect possible long-term changes in phytoplankton community structure (Figure 9, Figure 10 and Figure 11). For the three southern provinces, the highest average microphytoplankton fraction from 1998 to 2016 was showed in most productive region NWCS (40.49%), then followed by SW (21.73%) and GFST (18.83%) (Figure 9a,b,d). The seasonal cycles in microphytoplankton fraction in NWCS, SW, and GFST had pronounced peaks in April (NWCS, GFST) or April–May (SW), corresponding to the spring bloom timing. This is followed by an abrupt decrease in the fraction of microphytoplankton, which last the entire summer and early fall, and a less intense secondary bloom occurring around November–December (late fall, early winter) (Figure 9 and Figure 10). In the GFST and SW, the microphytoplankton cycle was less pronounced than in the NWCS, with the microphytoplankton fraction reaching at most 30% in the spring compared to a value of 60% in the NWCS (Figure 9 and Figure 10). The seasonal cycles in NADR, ARCT, and BPLR were truncated due to the presence of ice cover in winter months that prevent recording of satellite ocean-colour measurements. The peak of the microphytoplankton fraction occurred between May and June in the three provinces. In the NADR and ARCT provinces, the mean values between March and October from 1998 to 2016 were 35.51% and 37.59% respectively (Figure 9c,e). In BPLR, the mean value between June and October in the 19 years was 41.99% (Figure 9f).

For four of the six provinces (GFST, SW, NWCS, and BPLR) the trend in microplankton fraction showed an increase from the late 1990s to approximately 2010 followed by a sharp decline (Figure 11a,b,d,f). The two remaining provinces (ARCT and NADR) exhibited a less pronounced decrease after 2010, which lasted 2 to 4 years, before oscillating positively again (Figure 11c,e). Therefore, in addition to the temporal analysis over the entire period of observation (19 years, 1998–2016), we studied the trends of the time series for each province for two distinct periods corresponding to 1998–2010 and 2011–2016 (Figure 11). In general, the trend for the microphytoplankton fraction in GFST, SW, and NWCS were similar, with a significant increase until 2010s followed by a decrease, although there was an overall increasing trend over the full 19 year time-series (Figure 11a,b,d). Over the 19 years of observation, the slope of the microphytoplankton fraction versus time (in year) was positive and significant (p < 0.001) for four provinces, including NWCS (0.36%·yr⁻¹), ARCT (0.25%·yr⁻¹), SW (0.12%·yr⁻¹), and GFST (0.06%·yr⁻¹), and positive but not significant for BPLR (0.06%·yr⁻¹) and NADR (0.11%·yr⁻¹). From 1998 to 2010 the trends in all six provinces showed a positive slope with higher values found for the coastal provinces: NWCS (0.68%·yr⁻¹), BPLR (0.43%·yr⁻¹), and ARCT (0.41%·yr⁻¹) and smaller positive slopes for the pelagic provinces in descending order from the least to the most oligotrophic. After 2010, the relative abundance trends of microphytoplankton were negative in NWCS (−1.33%·yr⁻¹), BPLR (−0.84%·yr⁻¹), GFST (−0.64%·yr⁻¹), and SW (−0.53%·yr⁻¹). However, a significant positive slope was found in the NADR (0.87%·yr⁻¹) and a positive slope in ARCT (0.23%·yr⁻¹, p > 0.001) for the second truncation of the time-series.

4.3.2. Anomaly Analysis of Microphytoplankton in NWA

Overall annual trends of microphytoplankton fraction in the six provinces conceal large variation in the seasonal signal. Percentage anomalies of microphytoplankton in each ecological province are shown in Figure 12. The anomaly was calculated by subtracting the overall mean of each month (1998 to 2016) from each cell. In GFST, the anomalies from March to May showed the largest variations, positive values were found from March to May between 2010 and 2014 (Figure 12a). Generally, the anomaly in SW was negative before 2003, predominantly positive anomalies were found in 2010 and 2011(Figure 12b). In NADR, the anomalies showed a vast variation (Figure 12c). The NADR province represents a boundary area influenced by BPLR, NWCS, and SW waters [34], resulting in a highly dynamic climatological character. In NWCS, the anomalies were negative before 2006, while between 2008 and 2013 positive anomalies were observed in winter months (Figure 12d). It is noteworthy that the anomalies were all positive in 2011, which is consistent with the pronounced peak in 2011 for the time-series analysis results in NWCS (Figure 11d and Figure 12d). In ARCT, the anomalies showed positive values in 2010 and 2016, while negative anomalies were generally found before 2003, which is consistent with the trend shown in Figure 11e. In BPLR, the anomaly of microphytoplankton in September and October was positive from 2009 to 2011 (Figure 12f), which corresponds to the peak of the time-series trend within these years (Figure 11f).

4.3.3. Discussion on the Microphytoplankton Variation in NWA

Model E, was selected following the performance assessment study of nine PSC models. However, this study focused only on the ability of the models to derive the dominant PSCs, and not the accuracy of the retrieval of the PSCs. Therefore, we further compared the in situ and CCI-derived microphytoplankton (both concentration and fraction), to test the performance of model E (Figure 13). The in situ and CCI-derived microphytoplankton are consistent with each other, with a highly significant linear relationship for both concentration and fraction (r² = 0.73 and 0.72, respectively, p-value < 0.01). As observed with Chla concentration derived from ocean colour satellite, a systematic overestimation for microphytoplankton concentration and fraction was observed when in situ Chla was lower than ~0.5 mg·m⁻³, while CCI-derived values were generally underestimated when in situ Chla was higher than ~ 0.5 mg·m⁻³. Given that model E is an abundance-based model, the overestimation in the lower-end and the underestimation in the higher-end of the Chla concentration spectrum result from uncertainties in the OC_x Chla algorithm (Figure 6, Figure 7 and Figure 13). Indeed, this uncertainty on Chla directly propagate to microphytoplankton estimates. However, the effect of this systematic bias will not change the sign of long-term trends in each ecological province, but perhaps decrease its magnitude.

Our findings on microphytoplankton variation in NWA are consistent with the results of a remote-sensing-based study of phytoplankton phenology from 1998 to 2012 that identified a shift toward earlier spring bloom timing in the NWCS in recent years as well as increasing magnitude of fall/winter blooms in the entire NWA [53]. A comprehensive analysis of factors influencing phytoplankton phenology in the North Atlantic (including SST, photosynthetic active radiation, wind stress, and eddy kinetic) revealed that SST was the main driver of the deviance in inter-annual winter bloom magnitude in the NWA [53]. Between 2012 and 2014, SST anomalies (based on a 1985–2010 climatology) from July to November remained positive for the majority of the NWA region [54,55,56]. In contrast, in 2015 and 2016, SST anomalies from July to November in the ARCT and the Northwest of NWCS were negative [57,58]. This is consistent with the pronounced “fall-winter peak” in the fraction of microphytoplankton recorded in the early 2010s for the NWCS and ARCT provinces, and the relatively low microphytoplankton anomalies recorded in the winter months of 2015 and 2016 there (Figure 12c–e).

A previous study based on Continuous Plankton Recorder (CPR) data collected from 1958–2006 characterized variability in surface-water diatom and dinoflagellate abundance in an area corresponding to the NWCS [59]. Diatom abundance was higher in the 1990s and early 2000s than in the 1960s and 1970s in regions corresponding to the NWCS, due mainly to differences in abundance in the winter months [59]. The abundance of dinoflagellates was also significantly higher in the 1990s and early 2000s than in previous decades, but the differences were driven by shifts in abundance mainly in the first half of the year on the Scotian Shelf and throughout the year on the Newfoundland Shelf [59]. In our study, an increase of microphytoplankton fraction in the NWCS was observed in winter months in 2010s relative to the 2000s (Figure 12). In this context, although the methodologies and spatial coverage employed by these two studies were different, a general conclusion can be drawn that the fraction of microphytoplankton had been increasing in the last three or four decades, with a reversal in the trend around 2010. The increase in CPR-based diatom abundance in the 1990s and 2000s was attributed to higher stratification in the winter, which would enhance phytoplankton production by stabilizing the water column at a time when strong mixing would normally inhibit production [59]. In our study, persistent higher than average fall–winter microphytoplankton fractions were observed from 2009 to 2014, even as the trend in microphytoplankton fraction became negative due to below-average fractions observed in the summer, when stronger stratification would reduce the supply of nutrients to the euphotic zone and inhibit production of large phytoplankton. A decline in deep-water silicate concentrations that started in 2013 may have further inhibited production of diatoms [35].

5. Conclusions

Nine Phytoplankton Size Class (PSC) models were evaluated against in situ observations in the Northwest Atlantic (NWA). The PSC models’ performance was lower in productive waters (Scotian Shelf, Northwest Atlantic; within 4 degrees of the coastline) compared to the results obtained from the global open ocean. Model performance was especially poor for picophytoplankton-dominated communities (Table 3). This could be related to large errors in satellite-derived model inputs (e.g., Chla concentration and absorption coefficient of phytoplankton), and a general under-representation of pico-dominated samples in productive waters. Considering the accuracy of retrieving micro- and mixed nano–pico phytoplankton, Devred et al. 2011 [25] (model E) showed the best performance among the nine models for the three satellite products (i.e., SeaWiFS, MODIS, and VIIRS), which is not surprising given that this model was developed including a regional dataset.

The best performance model was then applied to the NWA daily OC-CCI images from 1998 to 2016, to demonstrate the spatial and temporal variations of microphytoplankton fraction. For the three southern ecological provinces, there was a significant increase in the microphytoplankton fraction until 2010s followed by a decrease, although there was an overall increasing trend over the full 19 year time-series. Strong anomalies in winter months between 2008 and 2013 were found in Northwest Continental Shelf, which could be driven by the changes in environmental factors (e.g., SST, stratification, and nutrient supply).

In a changing climate, long-term dynamics on the PSC variations would aid in the monitoring and investigating of water system changes. In addition, in our studied area, the selected model can be also applied at a smaller scale in representative areas and compared with changes in environmental drivers at more relevant scales, to further reveal the influence of environmental factors on the changes of PSC. Besides, our study indicates that the application of PSC models to productive waters has its limitations. As an active area of research, new ocean-color-based PSC models are under development. Ongoing efforts strive to include more data, more diverse conditions, and to improve model formulation and parameterization. Improved (both methodologically and quantitatively) in situ observations and satellite products in productive waters together will help to advance PSC research in the future.