Improving the Remote Sensing Retrieval of Phytoplankton Functional Types (PFT) Using Empirical Orthogonal Functions: A Case Study in a Coastal Upwelling Region

Abstract

An approach that improves the spectral-based PHYSAT method for identifying phytoplankton functional types (PFT) in satellite ocean-color imagery is developed and applied to one study case. This new approach, called PHYSTWO, relies on the assumption that the dominant effect of chlorophyll-a (Chl-a) in the normalized water-leaving radiance (nLw) spectrum can be effectively isolated from the signal of accessory pigment biomarkers of different PFT by using Empirical Orthogonal Function (EOF) decomposition. PHYSTWO operates in the dimensionless plane composed by the first two EOF modes generated through the decomposition of a space–nLw matrix at seven wavelengths (412, 443, 469, 488, 531, 547, and 555 nm). PFT determination is performed using orthogonal models derived from the acceptable ranges of anomalies proposed by PHYSAT but adjusted with the available regional and global data. In applying PHYSTWO to study phytoplankton community structures in the coastal upwelling system off central Chile, we find that this method increases the accuracy of PFT identification, extends the application of this tool to waters with high Chl-a concentration, and significantly decreases (~60%) the undetermined retrievals when compared with PHYSAT. The improved accuracy of PHYSTWO and its applicability for the identification of new PFT are discussed.

1. Introduction

Phytoplankton Functional Types (PFTs) represent an operational division that combines highly diverse taxa (species level) into groups which share traits (morphological, physiological, behavioral, and/or life-history) associated with ecological and/or biogeochemical functions (resource acquisition, predator avoidance, and metabolite production) [1]. The purpose of this division is to simplify community analyses and aid in the model building associated with climate change impacts on the biogeochemical and ecological components of oceans [2,3,4]. PFTs are usually defined using a combination of high phylogenetic grouping and functionality (for example, silicifying diatoms, mixotrophic dinoflagellates, nitrogen fixing and non-nitrogen fixing cyanobacteria, and calcifying coccolithophorids) or cell size [1,5,6,7]. Although size classification is not based on functional criteria, there is an alignment of the functional roles of phytoplankton with the size categories and the environmental niches (biogeochemical provinces) they occupy [1,8,9]. Altogether, the use of the PFT concept varies according to the scientific objectives of a given study and the observational capabilities available or required to tackle them [1].

Over the past two decades, the development of remote sensing techniques to identify PFTs has increased due to the wider spatial and temporal scales over which they can be used compared to in situ measurements. The most common satellite PFT approaches can be classified into abundance (size classes) and spectral (bio-optical) techniques [1,10,11]. Abundance-based approaches use satellite chlorophyll-a (Chl-a) as an input and, by using that, they exploit the largest signal in water-leaving radiance to extract the variability due to PFTs out of Chl-a. Spectral approaches are based on optical properties such as radiance, absorption, and/or backscattering spectra of phytoplankton, the variation of which is linked to community structure and pigment composition [10]. Different algorithms for the spectral-based discrimination of PTFs include the decomposition of absorption spectra, derivative analysis, inversion modeling, or Empirical Orthogonal Function (EOF) decomposition [11,12,13,14,15,16]. The validation of all these algorithms is a critical task that depends heavily on in situ measurements with sufficient spatial and temporal coverage [17]. Overall, the development and validation of PFT methods rely on HPLC (High-Performance Liquid Chromatography) pigment-based proxies of taxonomic composition or size structure and there is a clear need to complement these validations with supplemental datasets, including flow-cytometry, microscopy, and size-fractionated estimates [10,18,19].

Most algorithms perform similarly well over large gradients of co-varying bio-optical properties and reproduce expected trends in the global distribution of PFTs but not at smaller scales of variability [8]. Since different PTF algorithms use distinct approaches, datasets, and validation metrics, the evaluation of their performances requires comprehensive inter-comparisons using the same validation data and needs to consider errors associated with each one of them [10,20,21]. Thus, an appropriate choice of a PFT algorithm will depend on the scientific objectives and the observational capabilities available to validate its performance. PHYSAT is a spectral-based method that relies on the principle that spectral changes in anomalies of normalized water-leaving radiance (nLw) are due to changes in the dominant PFT [22]. This method was first developed for Case 1 waters using an empirical comparison of global SeaWiFS imagery with in situ data collected by the GEP&CO cruises (Geochemistry, Phytoplankton and Color of the Ocean [23]) along the North Atlantic, the Caribbean Sea, and the South Pacific. In its latest version, PHYSAT uses MODIS-Aqua data and discriminates between six phytoplankton groups (Nanoeukaryotes, Prochlorococcus, Synechococcus, Diatoms, Phaeocystis-like, and coccolithophorids) [24]. The current calibration of the PHYSAT algorithm has been updated using additional in situ measurements [25].

PHYSAT removes the dominant effect of Chl-a in the nLw spectrum to highlight the reflectance spectral anomaly generated by accessory pigments (that is, it exploits second-order anomalies of reflectance spectra) [22]. This proceeding is performed using an empirical reference model of nLw (nLw^ref) containing the average nLw at different Chl-a concentrations, based on a matchup between satellite records and in situ measurements. Since most of GEP&CO data come from open ocean sites with low Chl-a concentrations (<3 mg·m⁻³), the application of nLw^ref may not be appropriate in high Chl-a environments such as coastal and upwelling systems. In addition, coastal environments usually show a higher content of suspended inorganic matter and colored dissolved organic matter (CDOM), which could bias the PFT discrimination [10,25]. In the case of mid-latitude upwelling systems, most of the observed CDOM is transported from the subsurface by the same vertical Ekman fluxes that supply nutrients and increase biological production in the surface layer; therefore, CDOM could increase and covary on similar spatial and temporal scales as phytoplankton biomass [26]. To minimize the effect of CDOM, PHYSAT can be adapted by the development of regional nLw^ref models. Recently, a regionalized version of PHYSAT was developed for the Mediterranean Sea (PHYSAT-Med) based on the optical retrievals of the phytoplankton assemblages and their succession patterns in this basin [27]. Thus, it is expected that the development of new regional bio-optical models will improve the capabilities of PHYSAT.

In the present study, we propose and test an alternative method to improve the identification of PFTs done using the PHYSAT approach. In theory, this new approach does not require the construction of regional reference models because it allows for an adequate PFT identification in both low and high Chl-a environments. The improved method, called PHYSTWO (PHYSAT based on second-order spectral modes), relies on the assumption that the dominant effect of Chl-a in the nLw spectrum can be effectively isolated from the signal of accessory pigments—biomarkers of different PFTs—by using EOF decomposition. To test and validate this approach, the PHYSTWO performance is compared against that of PHYSAT based on a regional nLw model for upwelling conditions using global datasets, as well as in-situ observations, from the coastal upwelling system (coastal and coastal transition zones) off central Chile.

2. Data and Methodology

2.1. Satellite Data and In Situ Measurements

Daily level-2 products of the MODIS-Aqua mission with a spatial resolution of 1 km were obtained from the OceanColor Web (http://oceancolor.gsfc.nasa.gov) for the region off central Chile (72–76°W and 35–38°S). The MODIS-Aqua products used were the Chl-a concentration estimated by the standard algorithm OC3, the aerosol optical thickness (AOT) at 869 nm, and the remote sensing reflectances (Rrs, in sr⁻¹) at 412, 443, 469, 488, 531, 547, and 555 nm. The Rrs were converted to normalized water-leaving radiances (nLw), adding the respective MODIS-Aqua spectrally convolved values of nominal solar irradiance (Fo, in W·m⁻²·μm⁻¹) at each wavelength:

$$nLw=Rrs(\lambda )\ast Fo(\lambda )$$

(1)

The accuracy of PFT retrievals from both methods (PHYSAT and PHYSTWO) was assessed by comparing these with in situ data collected from different sources. Some of these data were obtained from in vivo spectral fluorescence profiles using a submersible spectro-fluorometer—the FluoroProbe (bbe Moldaenke GmbH, Kiel, Germany). This profiling instrument measures fluorescence emission at ~680 nm in response to excitation by light emitting diodes (LEDs) centered at ~370, 470, 525, 570, 590, and 610 nm. Its software provides estimates of CDOM or “yellow substances” and of the Chl-a concentration associated with the four PFTs derived from their fluorescence excitation spectra: green algae (Chlorophyta and Euglenophyta), brown algae (Bacillariophyta, Chrysophyta, and Dinophyta), blue algae (Cyanophyta), and red algae (Cryptophyta) [28]. In our study, the instrument did not register Cyanophyta because of a sensitivity problem (too high of a threshold level for its detection). Eighteen FluoroProbe profiles were obtained off Concepción (~36.5°S) as part of a cruise performed on an upwelling frontal area (PHYTOFRONT cruise) during 4–6 February 2014 (Figure 1). Six additional FluoroProbe profiles came from monthly samplings at a shelf time-series station (St. 18; 73°W and 36.3°S) during 2013. The CDOM data from the FluoroProbe profiles were used to assess its vertical distribution in the water column and its potential interference with the PFT retrievals.

Since the PFTs obtained with the FluoroProbe are not the same as those defined by PHYSAT, we used additional data collected during the PHYTOFRONT cruise to develop equivalences. These data were derived from the analyses of discrete water samples for PFT identification and enumeration using epifluorescence and inverted microscopy, as well as size-fractionated Chl-a and estimations of carbon biomass [29]. The FluoroProbe data were considered to represent two PHYSAT groups: the Nanoeukaryotes (including mostly flagellated green and red algae in pico-, nano-, and micro-size fractions) and the Diatoms (brown, non-flagellated algae, including nano- and micro-size fractions). PHYSAT originally defined Nanoeukaryotes as a group of small, flagellated algae containing characteristic carotenoids [24]. In the coastal region off Concepción, photosynthetic flagellates (including Chlorophyta and Cryptophyta) have been recognized as the largest contributors to the biomass of small phytoeukaryotes [29,30,31,32]. Hence, we conclude that the FluoroProbe Chl-a contributed by Chlorophyta and Cryptophyta (GA+C) represents the PHYSAT Nanoeukaryotes in this study. In addition, the FluoroProbe and PHYSAT Diatoms (Bacillariophyta) data can be considered as equivalent since this group is usually a dominant component of the microphytoplankton in the coastal zone during the spring–summer period [33]. FluoroProbe data from the upper 20 m layer for these two PHYSAT PFT (Nanoeukaryotes and Diatoms) were averaged and expressed as a percentage of the total Chl-a to compare their distribution with that derived from satellite observations (Figure 2).

For the remaining four PHYSAT PFTs, a direct comparison with in-situ data was not possible. During the PHYTOFRONT sampling, Synechococcus had low abundance/biomass and a lack of Prochlorococcus cells was presented [29]. Furthermore, no data on Coccolithophorids and Phaeocystis were obtained in this cruise. However, satellite estimates of Phaeocystis were compared using the available databases (BODC: British Oceanographic Data Centre; OBIS: Ocean Biogeographic Information System; OCB DMO: Ocean Carbon and Biogeochemistry Coordination and Data Management Office; Pangea: Data Publisher for Earth and Environmental Science; WOD09: World Ocean Database 2009; and US JGOFS: US Joint Global Ocean Flux Study) compiled by the MARine Ecosystem DATa (MAREDAT) initiative [34]. The dataset used here (141 samples) was restricted to samples collected during cloud-free days following the launch of the MODIS-Aqua mission (July 2002). Most of these correspond to coastal time-series stations. The Phaeocystis biomass data were transformed to Chl-a concentration values assuming a C:Chl-a ratio equal to 60:1 [35].

2.2. The PHYSAT Regional Model for the Coastal Upwelling System Off of Central Chile

To compare the PHYSTWO and PHYSAT retrievals, an adapted version of PHYSAT was developed using a new regional reference model (nLw^upw) for the coastal upwelling region off central Chile. This model was used instead of the standard nLw^ref to remove the Chl-a spectral contribution from the nLw retrieved by the satellite since Chl-a has the dominant (first order) contribution of the retrieved signal from color sensors in open oceans [36,37]. The regional nLw^upw model is based on the mean of the daily MODIS-A nLw during the three months of the main upwelling activity (January–March 2014), using 1490 Chl-a concentrations in the range between 0.1 and 15 mg·m⁻³ (Figure 3), following the approach described in Reference [27].

In the PHYSAT approach, the spectral removal of the Chl-a signal is achieved by dividing the observed nLw by the nLw^ref (or nLw^upw) to obtain the radiance anomalies (Ra) at each wavelength:

$$Ra(\lambda )=nLw(\lambda )/nLw{(\lambda )}^{upw}.$$

(2)

Alvain [22] showed that every PFT group can be associated with a specific Ra produced by the signal of their non-chlorophyll pigments. The obtained Ra are then compared with the acceptable ranges of Ra (Figure 4) for each PFT [24] to assign a dominant PFT at each pixel. Figure 5 displays the PFT fields obtained using nLw^ref and nLw^upw. The differences between the estimates of both models, as well as their accuracy, are discussed in Section 3.1.

2.3. Singular Value Decomposition (SVD) in Spectral Orthogonal Modes: The PHYSTWO Approach

The PHYSTWO method assumes that the Chl-a spectral signal, beyond being the dominant signal in the nLw retrieved by color satellites, is also independent (or not-correlated) from the signal of the other non-chlorophyll pigments (such as the ones used by PHYSAT to assign PFTs) or even from the signal of other substances such as CDOM in the surface ocean. From this perspective, it is expected that the Chl-a spectral signal can be separated from that of other pigments in the form of a dominant variability mode in an EOF analysis. An EOF analysis is a decomposition of a signal or dataset in terms of orthogonal modes which represent different (and, in theory, uncorrelated) fractions of data variability [38]. Usually, the first mode accounts for the largest fraction of variability. In nLw data from Case 1 (offshore) waters, this variability is expected to be associated to the Chl-a signal [36,37]. The remaining modes account for lower and uncorrelated fractions of variability and, therefore, they could contain the spectral signal of non-chlorophyll pigments that can be used to identify PFT assemblages.

PHYSTWO performs the EOF decomposition in a two-dimensional space–nLw matrix N (Equation (A2)). This matrix is a normalized version of the R matrix (Equation (A1)) which is built by reordering the satellite nLw fields, as is detailed in Appendix A. The N matrix is decomposed into three matrices using an SVD approach:

$$N=U\ast S\ast V^T.$$

(3)

The new matrix U has dimensions P × 7 (P rows × 7 columns) and contains the spatial information of the N matrix. S is a diagonal matrix of dimensions 7 × 7 and V is a square matrix of dimensions 7 × 7 containing the spectral information. As in a typical SVD, the singular values on the diagonal S matrix are proportional to the explained variance contained in each mode m:

$$\%VarEx{p}_m=\frac{S_m^2}{{\displaystyle \sum _{i=1}^m{S}_i^2}}\ast 100.$$

(4)

The SVD operates in a similar way on the R matrices constructed with the nLw data from different dates or geographic regions. In general, the first two resulting modes always account for ~99% of the explained variance (m1 ~ 60% and m2 ~ 39%) and the remaining 1% is scattered into the other five variability modes. The spatial pattern of the first two modes (U₁ and U₂) is recovered by rearranging the dimensionless values of the first two columns of the U matrix in an original longitude–latitude matrix (Figure 6). A detailed inspection shows a close relationship between U₁ and the Chl-a field (Figure 1), whereas U₂ shows a good spatial agreement with the PFT retrievals obtained by PHYSAT (Figure 5). Therefore, U₂ may contain the spectral signal of non-chlorophyll pigments useful for PFT identification (see Section 3.2 and Section 3.3 ). However, U₂ by itself does not provide information on the dominant PFT. Hence, it is not possible to establish reference values for each PFT (such as the PHYSAT acceptable ranges in Figure 4) in the U₂ field since dimensionless values obtained in the U, S, and V matrices can change depending on the numerical properties of the R matrix, the architecture of the computer hardware, and the numerical precision of the operating systems and programs used. For these reasons, SVD does not have a unique numerical solution. Any reference values must be incorporated before the SVD together with the values of the matrix R so that they undergo the same process of normalization and SVD that the R values have to ensure that they have comparable magnitudes in the resulting U-fields. The reference values must be organized within a synthetic matrix that can be attached to the R matrix without affecting the SVD.

2.4. The Synthetic Matrix of Typical nLw Values for PFT Categories and PFT Estimation from Orthogonal Models

Obtaining a reference nLw matrix useful for PFT identification is not an easy task since it requires observational information about abundance, biomass, and the spectral signal generated by different concentrations of several phytoplankton groups in waters with different optical properties. This type of information is generally scarce because it requires intense field measurements and laboratory analyses [34]. A first version of a synthetic matrix Rs, containing reference values for PFT, can be obtained from the PHYSAT acceptable ranges of radiance anomalies (Figure 4) [22,24,27] following the procedure detailed in Appendix B. This Rs matrix contains the expected radiances of each dominant PFT under different Chl-a concentrations. Therefore, it is a concentration–nLw matrix whose values have the same units (W·m⁻²·μm⁻¹·sr⁻¹) of the R matrix built with satellite nLw data. Although R is a space–nLw matrix and does not have the same size of Rs, both matrices have the same number of columns as they share the same wavelength frequencies. Therefore, it is mathematically possible to concatenate the Rs matrix after the R matrix to obtain a new combined matrix RRs:

$$RRs=\left({}_{R_S}^R\right).$$

(5)

The concatenation is a common variation of the EOF method used in the analysis of the co-variability (or joint variability) of two or more fields at a time [38]. Here, the concatenation was performed as a way to provide reference values before the SVD. As was done for the R matrix, the RRs matrix is also standardized by removing the mean and standard deviation (Equation (A2)) to obtain the NNs matrix. The SVD of the standardized matrix is achieved through:

$$NNs=U\ast S\ast V^T.$$

(6)

This decomposition results in a U matrix, containing in its last 1800 rows, the dimensionless reference values of PFT at different Chl-a concentrations. As for the decomposition of the N matrix, the first two modes of the NNs decomposition accounts for most of the explained variance (M1 = 58.7% and M2 = 39.8%).

In the dimensionless plane formed by the first two spatial modes U₁ and U₂, we can locate the orthogonal values resulting from the matrix R decomposition and the typical values for the PFT derived from the decomposition of the synthetic matrix Rs (Figure 7). The typical PFT values appear distributed from the lower right to the upper left corner in Figure 7. In this same direction, the Chl-a concentration increases and therefore, the typical values at low (high) concentrations are located in the lower-right (upper-left) sector. It is noticeable that the typical values of all PFTs are distributed almost parallel to each other, without any overlap between them. This is an inherited characteristic of the acceptable ranges of the Ra matrix which helps to avoid a double assignment of PFTs.

The lines formed by typical values in the U₁–U₂ plane can be considered as PFT orthomodels. The assignment of the PFTs corresponding to each point in the R matrix is easily achieved using the Euclidean distance to the orthomodels values. The PFTs assigned to any point will be that of the orthomodel with the closest typical value in the U₁–U₂ plane, as shown in Figure 8a. By rearranging the assignments in a longitude–latitude space, the first and unadjusted PHYSTWO estimation field is obtained (Figure 8b). A schematic diagram summarizing the PHYSTWO method is shown in Figure 9.

The unadjusted PFT estimation of PHYSTWO does not fully agree with the in-situ observations made during the PHYTOFRONT cruise (see Section 3.4). Some erroneous assignments are due to the fact that current orthomodels have no values in the region of high Chl-a concentration in the U₁–U₂ plane (since the typical values in the Rs matrix come from waters with relatively low Chl-a concentrations <3 mg·m⁻³; GEP&CO cruises). Additionally, some PFT reference values are out of the range of their known environmental preferences of background Chl-a concentration. It is possible to improve the PHYSTWO estimations by directly fitting the orthomodels in the U₁–U₂ plane by using additional information (field measurements or available databases) as described in the next subsection.

2.5. The Adjustment of PFT Orthomodels with In Situ Data

2.5.1. Diatoms and Nanoeukaryotes

Figure 10a shows the U₁–U₂ plane obtained from the SVD of an Rs matrix together with an R matrix composed by the matchups with the satellite nLw values corresponding to the place and date of the samples during the PHYTOFRONT cruise (16 stations) and St. 18 (four stations). Each matchup was obtained as the mean nLw value in a 3 × 3 pixels box (~9 km²) around the latitude and longitude of any station for the same date. The size of the circles represents the relative abundance of Diatoms and Nanoeukaryotes based on FluoroProbe profiles. In this plane, Diatoms are most abundant at stations located in the upper left corner, associated with higher Chl-a concentrations (3.5–10.5 mg·m⁻³). These stations are located beyond the Diatom orthomodel that considers Chl-a concentrations <3 mg·m⁻³. To fit the diatom orthomodel, an extension towards the region of higher Chl-a concentrations is required. This is achieved by incorporating new typical values following the path determined by Diatom proportion values of up to 95% in the interpolated relative abundance field, as represented by the orange dots in Figure 10b.

Figure 10 also indicates that the Nanoeukaryotes were relatively more important in stations located towards the right central and upper sectors in the U₁–U₂ plane. However, the upper end of the Nanoeukaryotes orthomodel extends to the region where our observations suggest a greater dominance of Diatoms. Therefore, the fit for Nanoeukaryotes was done by excluding the top points of the orthomodel where the Nanoeukaryotes proportion was <45% (white dots, Figure 10c). Unfortunately, the PHYTOFRONT data were not sufficient to modify the orthomodels in regions of moderate to low Chl-a concentrations. These adjustments should be addressed in future studies covering information from the mesotrophic and oligotrophic regimes. However, considering that most of GeP&CO observations come from open ocean regions, we expect that the current values of the orthomodels in this range of Chl-a concentrations do not require significant changes.

2.5.2. Phaeocystis

The worldwide distribution of Phaeocystis is poorly understood but the locations that it is abundant in have been reported in the Southern Ocean southward of 70°S and in the North Atlantic northward of 50°N [34]. In the coastal upwelling system off Peru, the presence of Phaeocystis globosa has been reported during iron-fertilization experiments [39]. However, further in situ data on Phaeocystis in the coastal upwelling area off Chile are not available [40]. Despite this lack of information, retrievals from PHYSAT and PHYSTWO (Figure 5 and Figure 9) suggest that this PFT may be present in our study region.

The SVD of the R matrix corresponding to the MAREDAT stations shows a core of high Phaeocystis abundance that coincides with the middle region (moderate Chl-a concentrations) of the Phaeocystis orthomodel. Since PHYSAT lacks direct evidence to constrain the acceptable ranges for Phaeocystis, these estimates were initially categorized as Phaeocystis-like [24]. The spatial coincidence between the high concentrations of Phaeocystis observed from MAREDAT data and the orthomodel suggests that the PHYSAT method made an adequate assignment of this PFT, at least in regions with intermediate background concentrations of Chl-a. Under this premise, the Phaeocystis orthomodel was fitted to just exclude values corresponding to pixels with high and low Chl-a concentrations in the area below 30% of the Phaeocystis abundance (dark red dots, Figure 11).

2.5.3. Prochlorococcus and Synechococcus

Photosynthetic cyanobacteria have a wide distribution and are dominant in the ocean regions with very low Chl-a concentrations where low nutrient availability does not favor the dominance of other large-size PFTs. For this reason, the orthomodel extension towards the upper left corner of high background Chl-a concentration is not consistent with what we know about the ecological niches of the dominant cyanobacteria genera, Prochlorococcus and Synechococcus [41]. Using molecular techniques, Bibby et al. [42] have shown an alternation in the dominance of Prochlorococcus and Synechococcus in open waters related to the background concentration of Chl-a and nutrients. Prochlorococcus only appears dominant at Chl-a concentrations <0.35 mg·m⁻³, while Synechococcus dominates at values >0.35 and <0.5 mg·m⁻³. However, in regions with iron limitation, Synechococcus can extend its dominance to the 0.26 and 0.5 mg·m⁻³ range thanks to the greater efficiency of its photosynthetic system under stressed conditions compared to that of Prochlorococcus. Based on these results, the Prochlorococcus orthomodel was restricted to the typical values associated with Chl-a concentrations between 0.1 and 0.3 mg·m⁻³ and the Synechococcus orthomodel to values between 0.26 and 0.5 mg·m⁻³ (Figure 12).

2.6. PFT Estimation with Adjusted Orthomodels

A new matrix of typical values (Rt) can be obtained with a reconstruction of the adjusted orthomodels. The fitted values in U₁ and U₂ were included in a new matrix Ut_1,2 and multiplied by their corresponding values in the S and V matrices to obtain the new synthetic Rt adjusted matrix:

$$Rt=U{t}_{1,2}\ast S_{1,2}\ast V_{1,2}^T$$

(7)

After the mean and standard deviation (which were previously removed before the SVD) are added back, Rt becomes a matrix of nLw radiances expected for each PFT at different background Chl-a concentrations. Hereafter, this new matrix can be used to concatenate with other R matrices and, in this way, will allow for an adjusted PFT estimation. In our study case, PFT estimations using Rt (Figure 13) show a large agreement with the in situ measurements rather than the PHYSAT estimations (see Section 3.4). This new PHYSTWO method and the Rt matrix of typical adjusted values were implemented as a Matlab function and are freely available at https://www.dropbox.com/sh/1j7sn3uo7fa60r5/AABwub8DEg1YOqsICY_PAnkEa?dl=0.

3. Results

3.1. Adaptation of PHYSAT to Upwelling Conditions

To obtain more accurate PFT estimates that could be used to compare and measure the performance of PHYSTWO, the performance of PHYSAT was improved by the development of a regional model nLw^upw for the coastal upwelling region off central Chile (35–38°S and 72-76°W). The new nLw^upw greatly expands the original PHYSAT Chl-a upper limit from 3 to 15 mg·m⁻³ and it allows PFT estimations at higher Chl-a concentrations than it is possible using nLw^ref. This is consistent with observations in the region where Chl-a concentrations >10 mg·m⁻³ are common in coastal waters during the upwelling season [33]. The resulting nLw^upw values were lower (about half) than the standard nLw^ref and displayed a dome pattern due to a central maximum at 469 nm (Figure 3b). This pattern is not unexpected since the shape and amplitude of the spectra are not unique to the global ocean but are highly dependent on the bio-optical environment of each region [25]. Furthermore, it is not possible to compare directly both reference models at a wavelength of 469 nm since this wavelength does not form part of the nLw^ref model originally proposed by Alvain et al. [22].

The distributions of the dominant PFTs based on the standard PHYSAT reference model (nLw^ref) and the regional model for the coastal upwelling and transition conditions (nLw^upw) during the dates of the PHYTOFRONT cruise are presented in Figure 5. The nLw^ref model (Figure 5a) produces PFT retrievals for only 26.9% of the valid pixels (pixels with valid nLw data and free of cloud interference), while over 70% of pixels remain without an assigned PFT, especially in the coastal zone where Chl-a exceeds 3 mg·m⁻³ (

Table 1. The percentage of phytoplankton functional types (PFT) retrievals obtained from the PHYSAT and PHYSTWO methods off central Chile for 4–6 February 2014.

Method	Retrievals	% Ret. 1	NEU (%)	PRO (%)	SLC (%)	DIA (%)	PHA (%)	COB (%)
PHYSAT nLwref	18,944	26.9	16.1	3.1	3.3	3.0	1.4	0.0
PHYSAT nLwupw	27,874	39.5	11.9	10.1	7.1	10.3	0.2	0.0
PHYSTWO unadj.	70,009	99.3	18.6	35.4	22.6	10.1	11.5	0.0
PHYSTWO adj.	70,009	99.3	37.7	8.9	8.6	38.4	5.7	0.0

¹ The percentage of retrievals in relation to the total number of pixels with valid nLw data and where satellite Chl-a estimation was possible.

). The retrievals in the coastal upwelling zone suggest a dominance of Nanoeukaryotes. Diatoms, on the other hand, dominated a high Chl-a filament and an anticyclonic subsurface mesoscale eddy, structures which were observed in the transition from the coastal to the oceanic zone. Phaeocystis was observed along the filament’s central axis. In the rest of the coastal and offshore areas, Synechococcus and Prochlorococcus were the dominant groups.

The PFT distribution using the nLw^upw model (Figure 5b) is similar to the one derived from nLw^ref. However, it provides a slight improvement regarding PFT retrievals (39.5%) in high Chl-a areas. One of the differences between the models was that the nLw^upw model showed a higher spatial coverage of Diatoms (10.3%) than the nLw^ref model (3.0%). Both PFT models (nLw^ref and nLw^upw) yielded a greater dominance of Nanoeukaryotes in the coastal zone, an expected result during the relaxation of the upwelling conditions [30]. During the PHYTOFRONT cruise, a relaxation of upwelling was detected, and this might have led to a slightly higher dominance of the nanoplankton Chl-a fraction compared with the microplankton in the upper layer [29]. The coastal retrievals of pico-prokaryotic phytoplankton by PHYSAT are in clear disagreement with the environmental preference of Prochlorococcus for Chl-a < 0.35 mg·m⁻³ [42]. Besides, during the PHYTOFRONT cruise, low abundances of Synechococcus and no cells of Prochlorococcus were detected by flow-cytometry [29]. These results suggest that PHYSAT based on either the standard or the regional reference models may not be fully appropriate for the estimation of PFT distributions in waters with high Chl-a concentrations as those sampled during the PHYTOFRONT cruise.

3.2. First Spectral Mode and Satellite Chl-a

The PHYSTWO approximation relies on the assumption that Chl-a is the dominant and independent signal in the nLw retrieved by ocean color satellites and, therefore, it is susceptible to be separated as the first variability mode through an SVD. The SVD of the space–nLw matrix composed by MODIS-A nLw data off central Chile, on the dates of the PHYTOFRONT cruise (4–6 February 2014), resulted in seven variability modes. The first two modes together accounted for most of the explained variance (>99%). The spatial pattern of the first mode (U₁, VarExp = 79.9%, Figure 6a) displays a close similarity with that of the Chl-a field (Figure 1), with higher U₁ values in the coastal zone and in the filament north of the anticyclonic eddy. Indeed, Chl-a displays a high correlation with U₁ (r = 0.71) and does not show a significant correlation with U₂ (r = 0.05; Figure 7a,b). This high correlation suggests that the first mode contained the major part of the Chl-a spectral signal and, therefore, the back reconstruction of the first-order matrix M₁ must contain nLw values associated to chlorophyll reflectance. This reconstruction is easily achieved using Equation (8), where U₁, V_1, and S₁ are the first column of these matrices:

$$M_1=U_1\ast S_1\ast V_1^T.$$

(8)

Since the M₁ matrix contains the Chl-a signal of the phytoplankton community, this could be used to obtain a first-order reference model nLw^M1 (Figure 14c) following the same proceeding as that used to obtain the regional upwelling model nLw^upw in Section 2.2. When these models are compared, it is possible to observe that nLw^M1 has a close similarity to nLw^upw, in both shape and magnitude, showing the same spectral threshold at the 469 frequency, which determines a high correlation between both models (r = 0.91; Figure 14d). This spectral similarity corroborates that the M₁ mode—where the nLw^M1 comes from—contains the spectral signal of Chl-a.

The small differences between nLw^M1 and nLw^upw may be linked to other substances with a strong effect on the spectral nLw radiances such as suspended sediments or CDOM. The FluoroProbe profiles show an increase of CDOM with depth (Figure 2), leading to a low but significant inverse correlation with the vertical distribution of Chl-a (

Table 2. The correlation coefficients of the relationship between Colored Dissolved Organic Matter (CDOM), measured from FluoroProbe profiles (CDOM-FP, 0–60 m depth) and integrated in the upper 20 m layer (CDOM-20m), and the in-situ measured FluoroProbe Chl-a (Chl-FP), the satellite MODIS-A chlorophyll-a (Chl-sat), and the first two spatial modes obtained from the singular value decomposition (SVD) of a R matrix composed by MODIS-A nLw data in the region off central Chile (35–38°S, 72–76°W).

	Chl-FP	Chl-sat	U1	U2
CDOM-FP	−0.29 S	-	-	-
CDOM-20m	-	−0.14 NS	−0.11 NS	0.22 NS

^S The correlation is significant at the 0.05 level; ^NS Not significant.

). This pattern has been previously reported as a typical condition in upwelling systems [26]. However, the low CDOM concentrations in the upper 20 m do not have significant relationships either with the MODIS-A Chl-a or with the first two spatial modes. This suggests that the CDOM signal in the study region may represent a negligible contamination on the MODIS-A retrievals and in the subsequent PFT estimation based on the SVD of the R matrix, as performed by PHYSTWO.

3.3. Use of the Second Spectral Mode for Analyzing PFT Spatial Succession

The variability associated with the second spatial mode (U₂, VarExp = 39.1%) is, by definition, orthogonal and independent from U₁ variability. Since U₁ contains the variability associated with Chl-a (see the previous section), we expect U₂ to contain the signal associated with non-chlorophyll accessory pigments that can be used as indicators of the PFT. Evidence of this can be found in the spatial pattern of U₂ (Figure 6b), which closely resembles the PFT field derived from PHYSAT (Figure 5). Comparing these two fields, we observe that high U₂ values in the filament zone correspond to regions dominated by Diatoms, while low U₂ values fall in regions dominated by Nanoeukaryotes and intermediate values by smaller phytoplankton. This spatial coherence is expected and relies on the fact that PHYSAT is also designed to highlight the effect of secondary pigments that are overshadowed by the Chl-a signal [22]. Unfortunately, there is no way to quantify the degree of agreement between U₂ and PFTs because U₂ is a continuous field (that is, not categorical like the PFT estimations of PHYSAT) and because our results show that PHYSAT estimations of PFTs that could be used in this comparison may not be appropriated for the region. Indirect evidence of this agreement is obtained from the fact that an accurate PFT identification based in the U₂ field is possible, as performed by PHYSTWO.

3.4. PFT Accuracy for the PHYSAT and PHYSTWO Methods

PFT retrievals obtained by the PHYSAT (using nLw^ref and nLw^uwp reference models) and PHYSTWO (using unadjusted and adjusted orthomodels) were compared with the relative abundance of Nanoeukaryotes and Diatoms observed along the PHYTOFRONT transects off central Chile (

Table 3. The comparison between PFT retrievals with in situ and published data. Ret = the total number of retrievals on PHYTOFRONT transects. % Agr. = The percentage of retrievals on PHYTOFRONT transects that are in agreement with the in situ dominance of Diatoms and Nanoeukaryotes measured in the FluoroProbe profiles. % PAgr. = the percentage of retrievals that are in agreement with the areas of likely dominance by Prochlorococcus (where Chl-a concentration is <0.35 mg·m⁻³) and Synechococcus (where Chl-a is between 0.26 and 0.5 mg·m⁻³).

Method	NEU		DIA		PRO		SLC
	Ret.(1)	% Agr.	Ret.	% Agr	Ret.	%PAgr.	Ret.	%PAgr.
PHYSAT nLwref	64	6.1	0	0.0	2169	10.3	2305	9.7
PHYSAT nLwuwp	56	0.0	9	1.9	7111	13.8	5028	8.0
PHYSTWO unadj.	73	0.0	93	19.7	24,975	55.7	15,955	11.9
PHYSTWO adj.	174	81.6	327	68.9	6289	29.9	6051	22.1

). The Synechococcus and Prochlorococcus groups were not dominant or present in these transects but they are expected to be present offshore, out of the direct influence of upwelling waters and where low nutrient conditions coincide with their known environmental preferences. Therefore, PFT retrievals for these groups were compared with the potential dominance they may have in areas defined by their ecological preferences in terms of Chl-a concentration (as a proxy of nutrient concentration) as reported by Bibby et al. [42]. For the remaining groups (Phaeocystis and Coccolithophorids), we were not able to perform any validation since neither in situ nor published data were available for the study region. However, considering that retrievals for Phaeocystis are obtained by both PHYSAT and PHYSTWO estimations, besides the fact that previous studies have reported its presence in upwelling system, it is expected that this group could be present in the region.

Using the nLw^ref reference model, PHYSAT does not produce any Diatom retrievals for the PHYTOFRONT transects, although this group was the dominant component in 90.6% of the FluoroProbe profiles (Figure 2). There were 64 Nanoeukaryote retrievals for the transects but only 6.1% of these were observed as dominant in the FluoroProbe profiles. When the regional model nLw^upw is used, PHYSAT slightly increases the Diatoms retrievals and the percentage of agreement with FluoroProbe profiles to 1.9% but it does not produce any coincidences for Nanoeukaryote retrievals. The two PHYSAT models produced an important number of retrievals for Prochlorococcus and Synechococcus. However, only a low percentage (13.8% and 8.0%, respectively) of such retrievals coincided with their environmental preference, suggesting a misassignment for these groups.

PHYSTWO estimations produce retrievals over 99.3% of the area, which is >60% greater than that achieved by PHYSAT (Table 1). This large number of PFT retrievals is due to the SVD performed by PHYSTWO, which does not produce undetermined values as is the case with PHYSAT as an effect of removing the Chl-a through the division of the satellite nLw observations by the reference models (Equation (2)). However, PHYSTWO retrievals with unadjusted orthomodels are similar to PHYSAT estimations and, therefore, they still do not agree with the in-situ observations. The percentage of agreement in the case of Diatoms slightly increases with regard to PHYSAT but remains low (~19.7%), whereas there is no agreement with the retrievals for Nanoeukaryotes (Table 3). Only the Prochlorococcus agreement increased over 50% but, as observed in Figure 8, there were also an important number of retrievals close the coast in waters with high Chl-a concentration, where this group is unlikely to be found.

When adjusted orthomodels are used, PHYSTWO retrievals of Diatoms and Nanoeukaryotes contribute to the increase in the percentage of agreement to 81.6% and 68.9%, respectively. In the case of Synechococcus, the agreements remain similar, while those of Prochlorococcus decrease, although they are still higher than those reached by PHYSAT. Thus, PHYSTWO produces a PFT distribution closer to that expected from the in situ measurements during the PHYTOFRONT cruise and from the known environmental preferences of the PFTs (Figure 13). This type of PFT estimation displays a dominance of Diatoms followed by Nanoeukaryotes in coastal areas. The extension of coastal waters through the filament is also dominated by Diatoms with small patches of Phaeocystis. Finally, the most offshore areas are dominated by Synechococcus and Prochlorococcus, including the intrusions of oceanic waters in the southern section of the eddy, an effect of its circulation pattern [29].

4. Discussion

Since the first satellite-color missions, satellite data in the visible spectrum have been used to provide valuable information about the optical properties and the main compounds found in the surface ocean. Chl-a concentration as a proxy of phytoplankton biomass has been by far the most used product from ocean-color satellites and, as the first biological variable measured from space, it triggered an important breakthrough in our knowledge and perception of the ocean as the largest ecosystem on Earth. An evolving interest in retrieving information on other properties, including the composition of phytoplanktonic communities, has emerged in recent years [1]. This interest relates to the need of characterizing the spatial and temporal variability of primary producers, which in turn influence biogeochemical cycles, ecosystem structures, trophic transfers, and the climate.

Satellite Chl-a measurements are currently straightforward and obtained through empirical algorithms [43,44] since this is the main photosynthetic pigment in the upper ocean [36]. The Chl-a signal in the nLw radiances retrieved by satellite dominates most of the ocean landscape and it overshadows the signals of other phytoplanktonic (accessory) pigments, which are useful in PTF identification. The PHYSAT approach removes the effect of Chl-a by normalizing nLw data using a reference Chl-a spectrum between 412 and 555 nm as a way to highlight the signal of the accessory pigments [22]. However, the normalization process makes the information on the environmental Chl-a concentration unavailable for further PFT identification. Chl-a concentration is a very important piece of information on the environmental regime of the phytoplanktonic community under study since, in general, the dominance of a given PFT is associated to particular nutrients and Chl-a concentrations (that is, diatoms are likely the dominant PFT in high nutrient and Chl-a conditions while cyanobacteria are expected to dominate under low nutrient and Chl-a concentrations) [45,46].

In the present study, we have applied PHYSAT to an area ranging from nutrient-rich nearshore waters to offshore nutrient-poor waters in a temperate coastal upwelling region. We have shown that the main errors in the PHYSAT estimates are caused by PFT assignments to environments where these PFTs are not dominant or less likely to live. We have also provided evidence supporting the idea that the signal of the accessory pigments is spectrally uncorrelated with that of Chl-a so that both signals (Chl-a and non-Chl-a) conform independent modes of spectral variation in the nLw data. Based on this, the proposed PHYSTWO approximation isolates the Chl-a signal into a dominant (first) spectral orthogonal mode (U₁) and together with a second mode (U₂), containing the signal of the accessory pigments, performs a PFT estimation in the U₁–U₂ plane. Thus, the main advantage of PHYSTWO over PHYSAT is that it does not require the removal of the Chl-a spectral effect from the nLw data. This allows PFT estimations over a wide scope of environmental conditions with different Chl-a concentration. In theory, the U₁–U₂ plane contains all possible spectral combinations of the seven nLw wavelengths, but it is also a bi-dimensional arrangement for all possible combinations of the spectral signal coming from different phytoplankton communities under the diverse ecological regiments that characterize the oceans. Thus, PHYSTWO appears to be more suitable than PHYSAT to perform PFT estimations in complex regions like coastal upwelling regions where the optical properties change quickly from high to low Chl-a concentrations in a coastal to offshore gradient.

The presence of other optically significant agents in surface waters, such as detritus and colored dissolved organic matter (CDOM) may bias the remote detection of phytoplankton groups [25]. In the case of the upwelling system off central Chile, our observations indicated an increase of CDOM with depth and had a low correlation with the satellite surface Chl-a. Despite this, the low CDOM levels in surface waters may account for the spectral maximum at 469 nm observed in the regional reference model nLw^upw, since the absorption by CDOM is large in the blue part of the spectrum. However, the PFT estimation using this regional model produces a similar estimation to that using the standard nLw reference model, which suggests that the CDOM represents only a low interference for PFT estimation in this region. It is likely, however, that this is not the case for all coastal waters, especially for those influenced by river discharges, where higher CDOM concentrations are expected and their optical characteristics could be more similar to the case-2 waters. Unfortunately, the few CDOM observations included in our study do not allow an appropriated analysis of the CDOM interference but we expect that further investigations will allow the inclusion of CDOM as an additional PHYSTWO orthomodel and, with it, improve the PFT estimations.

The PFT estimation through PHYSTWO relies on the use of orthomodels in the dimensionless plane U₁–U₂. This plane is a continuous field which provides information on the spatial changes associated with the transition between different spectral regions. As such, this approximation could be used to identify the relative abundance over a greater number of phytoplankton assemblages compared to the categorical dominance obtained by PHYSAT. Hence, PHYSTWO can potentially provide information on the spatial succession or co-dominance in the phytoplankton communities. However, the orthomodels used in this study to perform the PFT estimation are in effect an adjusted version of the PHYSAT acceptable ranges [22,24,27], so they cannot identify the co-dominance regions in their current state. Future PHYSTWO versions could be evolved towards the estimation of mixed PFT dominance in different assemblages when detailed information on in situ phytoplankton community structure, or on accessory pigments, becomes available to build new orthomodels.

The adjustments and validation of the PHYSTWO orthomodels in a given region rely on the availability of in-situ information on the phytoplankton community structure, which in general, is scarce for the oceans as a whole. In the present study, the orthomodel adjustments were done employing different sources of data. The Phaeocystis orthomodel was adjusted using the global databases collected in the last 13 years; those of Cyanobacteria orthomodels were adjusted using genomic studies performed at the basin scale; and in the case of Diatoms and Nanoeukaryotes, the adjustments were done using in vivo fluorescence obtained from a spectro-fluorometer (FluoroProbe profiler; bbe Moldaenke GmbH, Kiel, Germany ). Altogether, the improvement of PHYSTWO over PHYSAT in PFT identification in the region of study and using local, as well as global databases, suggest that it has the potential for further improvements with larger datasets and that it can be applied to other regions of the oceans, even in its present format. On the one hand, our results indicated that spectro-fluorometers provide suitable information for this purpose because it reports the relative abundance of the phytoplankton groups in terms of their chlorophyll content. This is a desirable characteristic that allows a direct estimation of the degree of local spectral dominance without the requirement of bio-volume transformations from biomass information. Previous studies on the comparison of different techniques for assessing phytoplankton community structures have shown that FluoroProbe is potentially capable of determining its general characteristics when compared to HPLC analyses [47,48,49]. However, its application improves greatly when species-specific calibrations under different environmental conditions are considered [49]. In the case of the FluoroProbe data used in this study, a comparison of the Chl-a estimates with those obtained from fluorometric analyses revealed very similar patterns of distribution, though the magnitudes were higher in the first case. Additionally, the surface satellite and in situ Chl-a data presented a similar distribution [29]. Differences in the magnitude of the estimates do not affect PFT assignation since orthomodels represent the relative abundances of the different groups with respect to the total.

The PHYSTWO approach incorporates a synthetic matrix, initially based on the PHYSAT acceptable ranges, which after it undergoes an SVD, produces orthomodels with reference values that are useful to perform PFT identification. These orthomodels operate in the U₁–U₂ plane and they can be easily adjusted in this plane using additional new in-situ observations and published data. Because the U₁–U₂ plane has a finite space, the changes in the scope of any orthomodel (that is, with the addition or removal of orthomodel points) affect the scope of the neighbor orthomodels, even though the latter is not modified. After adjustment, we have shown that the orthomodels are suitable to produce accurate PFT estimations in environments with high Chl-a concentrations. However, these adjustments do not constrain PFT identification to this kind of environment. The orthomodels were modified mostly in their top-left ends to ensure an adequate PFT identification in high Chl-a waters (>3 mg·m⁻³), leaving their bottom-right values unaltered, where the reference values for clear waters (Chl-a concentrations between 0.1 and 3 mg·m⁻³) are located. Thus, PFT identification by PHYSTWO and PHYSAT remained similar in environments with low Chl-a concentrations. Altogether, the adjustments performed to orthomodels make PHYSTWO suitable to operate over a wider range of environments than PHYSAT and it does not require the use of regional reference models to acknowledge the particular optical differences between ocean regions (Figure 13c). However, further adjustments to PHYSTWO orthomodels beyond those performed in the present study may be required to increase the improvement in PFT estimation. Potential users of PHYSTWO need to determine if the error inherent in this technique is acceptable for the type of research being conducted. Matlab codes for PHYSTWO and PHYSAT are freely-available at https://www.dropbox.com/sh/1j7sn3uo7fa60r5/AABwub8DEg1YOqsICY_PankEa?dl=0.

5. Conclusions

Knowledge on the phytoplankton community structure in the oceans and on its spatial and temporal scales of variability directly contribute to more accurate assessments of the state of pelagic ecosystems, the extent of carbon cycling, and the impacts of climate change. PHYSAT is a pioneering method in the identification of PFTs using remote sensing retrievals of radiance reflectance spectra. PHYSAT uses reference data from open ocean environments and it fails to produce appropriate estimates in waters that are optically more complex, such as those in highly productive coastal waters. An adaption of the PHYSAT method to more optically complex waters through the construction of regional reference models does not substantially improve PFT estimation in the case of coastal upwelling regions, as assessed during the present study. The proposed PHYSTWO method represents an alternative approach applicable to systems with a larger range of Chl-a concentration and without requiring the construction of additional regional reference models. Since the signal of Chl-a and non-chlorophyll pigments conform to different modes of the spectral variability retrieved by satellites, PHYSTWO separates these signal by employing an orthogonal decomposition of space–nLw matrices constructed with satellite nLw information. With the use of adjusted orthomodels in the dimensionless plane formed by the first two modes. PHYSTWO produces PFT retrievals with a better agreement regarding in-situ observations than those obtained by PHYSAT in the coastal upwelling region under study. Furthermore, PHYSTWO reduces the number of indeterminate assignments generated by PHYSAT by more than 60%. PHYSTWO orthomodels can be easily adjusted to improve PFT determination and can be complemented by the incorporation of new PTFs beyond the six categories defined by PHYSAT if additional in situ data are employed. These aspects will require further assessments of the PHYSTWO method in different regions and environmental conditions.