1. Introduction

Changes in marine phytoplankton species compositions can be triggered by changes in the Earth’s climate [1–3]. Long-term changes in phytoplankton community composition contribute to changes in food webs and air–land–sea carbon cycles [4]. The northern Arabian Sea, for example, appears to be undergoing a shift in phytoplankton species composition from diatoms to Noctiluca miliaris during the annual winter (November to March) Northeast Monsoon (NEM) [5–7]. Blooms of N. miliaris are disrupting the traditional diatom-dominated food chain during the NEM and altering the magnitude of carbon export [8]. Ship- and aircraft-based sampling alone cannot produce sufficient data records to study such phytoplankton diversity shifts on large temporal and spatial scales. Following, the oceanographic community has invested in the development of methods to identify and discriminate between members of the phytoplankton community using remotely sensed data records [9]. Satellite ocean color instruments provide consistent and high-volume data records on scales that far exceed current ship and aircraft sampling strategies, with time-series of sufficient length to allow retrospective analysis of oceanographic trends. The imagery captured by the NASA Sea-Viewing Wide Field-of-View Sensor (SeaWiFS; 1997–2010) and the Moderate Resolution Imaging Spectroradiometer onboard Aqua (MODISA; 2002–present), for example, provide viable data records for observing almost two decades of changes in the biogeochemistry of both global and regional marine ecosystems [10–13]. The community considers ocean color satellite estimates of phytoplankton community composition to be sufficiently critical for advancing our understanding of biogeochemistry and carbon cycle science that all forthcoming ocean color satellite programs are now required to have the capability to discriminate between phytoplankton groups [14,15].

A variety of approaches exist for identifying phytoplankton communities from satellite ocean color data records. Briefly, satellite ocean color instruments measure the spectral radiance exiting the top of the atmosphere ($L_t(\lambda )$; $\mathrm{\mu W}{\mathrm{cm}}^{-2}{\mathrm{nm}}^{-1}{\mathrm{sr}}^{-1}$) at discrete visible and infrared wavelengths. Atmospheric correction algorithms are applied to $L_t(\lambda )$ to remove the contribution of the atmosphere from the total signal and produce estimates of remote-sensing reflectances $R_{\mathrm{rs}}(\lambda )$; ${\mathrm{sr}}^{-1}$), the light exiting the water mass was normalized to the hypothetical condition of an overhead Sun and no atmosphere [16,17]. Bio-optical algorithms are applied to the $R_{\mathrm{rs}}(\lambda )$ to produce estimates of additional marine geophysical properties, such as the near-surface concentration of the phytoplankton pigment chlorophyll-a ($C_{\varphi }$; $\mathrm{mg}{\mathrm{m}}^{-3}$) [18,19] and inherent optical properties (IOPs: the spectral absorption and scattering characteristics of ocean water and its dissolved and particulate constituents) [20,21]. Several methods for identifying phytoplankton communities do so using empirically derived thresholds on the magnitudes of the derived $C_{\varphi }$ or IOPs [22–27], whereas others interpret the retrieved radiometric spectral shapes [28–37]. In general, the former, abundance-based methods exploit observed relationships between the trophic status of the environment and the type of phytoplankton expected to be present, whereas the latter, spectral methods exploit differences in the optical signatures of specific size classes or functional groups to distinguish between phytoplankton types [9]. Corresponding studies related to changing ecosystems or the appearance of new species most often employ spectral methods, as abundance methods cannot capture the desired signals when the emerging environmental relationships do not appear in their in situ training data sets.

Ocean reflectance inversion models (ORMs) provide a common spectral method for inverting the “color” of the water observed by a satellite [e.g., $R_{\mathrm{rs}}(\lambda )$] into marine IOPs through a combination of empiricism and radiative transfer theory [20,21]. Generally speaking, ORMs attribute variations realized in the spectral shape of $R_{\mathrm{rs}}(\lambda )$ to varying marine populations of phytoplankton, nonalgal particles (NAP), and colored dissolved organic material (CDOM), all of which maintain unique optical signatures. They operate by assuming spectral shape functions of the constituent absorption and scattering components and retrieving the magnitudes of each constituent required to match the spectral distribution of $R_{\mathrm{rs}}(\lambda )$. Therefore, unlike abundance methods, ORMs can, in principle, discriminate between phytoplankton groups with common abundances, provided the groups present have contrasting optical signatures within the spectral bands detected. This ability to further deconvolve $a_{\varphi }(\lambda )$ into contributions by individual phytoplankton groups, however, remains inadequately demonstrated for a large range of environmental conditions [35,37]. When ORMs include only a single spectrum for a phytoplankton group, they remain confounded by natural variations in the spectral characteristics of that group due to growth stage, nutrient availability, and ambient light history. Furthermore, their solution can be statistically ambiguous [i.e., several combinations of IOPs can produce similar $R_{\mathrm{rs}}(\lambda )$] and dependent upon the wavelength suite (number and position) used in the inversion [38,39].

Here, we evaluate the ability of a common ORM to discriminate between two phytoplankton groups residing under varied biophysical conditions. Our goal was twofold: evaluate the ability of an ORM to quantify blooms with variable depths, thicknesses, and ages; and verify the spectral requirements for success using an ORM. For this case study, we selected a well-vetted ORM algorithmic form that we parameterized specifically to identify N. miliaris in a mixed phytoplankton community, using satellite ocean color data records collected in the Arabian Sea [40]. Briefly, we: (1) synthesized a series of vertical profiles of IOPs that represent a wide variety of bio-optical conditions for the northern Arabian Sea; (2) generated $R_{\mathrm{rs}}(\lambda )$ from these profiles using Hydrolight-Ecolight 5 (HE5) [41]; (3) applied the ORM to the synthesized $R_{\mathrm{rs}}(\lambda )$ to estimate the relative presence of diatoms and N. miliaris for each example; and (4) repeated the third step using $R_{\mathrm{rs}}(\lambda )$ with reduced spectral resolution. Comparing the estimates from the inversion model with those from the synthesized vertical profiles provided a mechanism for identifying those bio-optical conditions under which the ORM performed both well and poorly. While we focus on two specific phytoplankton populations, our ORM otherwise maintains the common form of many established approaches [20,21] and, thus, we expect our results to provide a case study that transfers well to other algorithms and phytoplankton populations.

2. Methods

A. Study Site and In Situ Data

Field experiments crafted to study the transition from N. miliaris to diatoms in the northern Arabian Sea were conducted during the 2011 NEM onboard the Indian Fishery Oceanographic Research Vessel Sagar Sampada, which traveled offshore following a northwest transect from Goa, India. Sampling dates ranged from March 7th–9th, 2011 and encompassed diatom-dominant, N. miliaris-dominant, and mixed population stations. We acquired measurements of $R_{\mathrm{rs}}(\lambda )$ and spectral absorption coefficients during this field campaign, using data collection and processing methods that appear in Roesler et al. [40] and Thibodeau et al. [7], respectively. With regard to the latter, we measured the spectrophotometric absorption of particles from discrete water samples using the quantitative filter technique, as modified by Mitchell [42] and Roesler [43]. We determined phytoplankton and non-algal absorption fractions by extraction [44] and identified and quantified phytoplankton species microscopically [5]. We assigned taxonomic dominance to each in situ sample using the microscopic species enumerations. We also collected discrete water samples via Niskin bottles for analysis of $C_{\varphi }$ and other pigments by the NASA high performance liquid chromatography facility at Goddard Space Flight Center. As in Thibodeau et al. [7], we used these measurements to derive average $C_{\varphi }$-specific absorption spectra for N. miliaris ($a_{\varphi N}^{*}(\lambda )$; ${\mathrm{m}}^2{\mathrm{mg}}^{-1}$) and diatoms ($a_{\varphi D}^{*}(\lambda )$; ${\mathrm{m}}^2{\mathrm{mg}}^{-1}$) (Fig. 1).

 

Fig. 1. Absorption spectra for N. miliaris and diatoms, normalized to 440 nm [7]. Thin blue lines show spectra from in situ stations dominated by diatoms, with the thick blue line indicating the mean spectrum. Thin orange lines show spectra from in situ stations dominated by N. miliaris, with the thick red line indicating the mean spectrum. The average $a_{\varphi D}^{*}(440)$ and $a_{\varphi N}^{*}(440)$ were 0.052 ($\pm 0.022$) and 0.066 ($\pm 0.015$) ${\mathrm{m}}^2{\mathrm{mg}}^{-1}$, respectively. Vertical solid lines indicate center MODISA wavelengths. Vertical dotted lines show the additional wavelengths considered in this study.

Download Full Size | PDF

The species N. miliaris appearing in the Arabian Sea is a mixotrophic dinoflagellate with green prasinoxanthin-containing symbionts (Pedinomonas noctilucae), which is physiologically different than the pink, bioluminescent variety more commonly found in temperate coastal waters [45] and optically different than most other resident diatom species, including mixtures of Rhizosolenia spp., Thallassiothrix spp., Skeletonema spp., and Chaetoceros spp. We acknowledge that N. miliaris and diatoms alone cannot universally represent the phytoplankton community of the northern Arabian Sea at all times and that single $a_{\varphi }^{*}(\lambda )$ cannot universally represent a phytoplankton group at all times. However, the contrasting optical signatures of N. miliaris and diatoms and their known colocation under many circumstances provided an ideal “real world” scenario for our sensitivity analyses.

B. Ocean Reflectance Inversion Model

Our ORM adopts the generalized IOP (GIOP) form described in Werdell et al. [21]. Briefly, ocean color satellite instruments provide estimates of $R_{\mathrm{rs}}(\lambda )$, which we convert to their subsurface values using the method presented in Lee et al. [46]:

Subsurface remote-sensing reflectances relate to marine IOPs, following Gordon et al. [47]:

Similar to absorption, total backscattering can be expanded to:

Four unknowns remain in Eqs. (1)–(5) after we assign the eigenvectors: $M_{dg}$, $M_{bp}$, $M_{\varphi D}$ (hereafter $C_{\varphi D}$), and $M_{\varphi N}$ (hereafter $C_{\varphi N}$). Using $R_{\mathrm{rs}}(\lambda )$ as input, we estimate these eigenvalues via nonlinear least squares (Levenberg–Marquardt) inversion of Eq. (2) [21]. We retained only those solutions with viable estimates of $M_{dg}$, $M_{bp}$, $C_{\varphi D}$, and $C_{\varphi N}$ (e.g., $-0.05a_w(\lambda )\le a_{dg}(\lambda )\le 5{\mathrm{m}}^{-1}$) that resulted in reconstructed $R_{\mathrm{rs}}(\lambda )$, which differed from the input $R_{\mathrm{rs}}(\lambda )$ by less than 33% for all wavelengths between 400 and 600 nm. We reconstructed $R_{\mathrm{rs}}(\lambda )$ using the retrieved eigenvalues as input into Eqs. (1)–(5) and defined failure as nonconvergence in the inversion. Werdell et al. [21] outlined the similarities in form and accuracy shared by this ORM configuration and other common approaches (e.g., [19,20,31,37,46]).

C. Synthesis of IOP Profiles

We synthesized a series of vertical profiles of spectral IOPs, representing a wide variety of bio-optical conditions for the northern Arabian Sea, for use as input into the “IOP Data” model of HE5. Specifically, we constructed vertical profiles of absorption and attenuation by particles plus CDOM ($a_{pg}(\lambda )$ and $c_{pg}(\lambda )$, respectively; ${\mathrm{m}}^{-1}$):

Unlike our ORM, we deconvolved $a_{dg}$($\lambda ,z$) into its two components:

Construction of $a_{pg}(\lambda ,z)$ and $c_{pg}(\lambda ,z)$ using Eqs. (6)–(10) required defining 14 eigenvectors and eigenvalues: $C_{\varphi D}(z)$, $C_{\varphi N}(z)$, $a_{\varphi D}^{*}(\lambda )$, $a_{\varphi N}^{*}(\lambda )$, $b_d(555,z)$, $b_{\varphi D}^{*}(555)$, $b_{\varphi N}^{*}(555)$, ${\eta }_d$, ${\eta }_{\varphi D}$, ${\eta }_{\varphi N}$, $a_d(443,z)$, $a_g(443,z)$, $S_g$, and $S_d$ (Table 1). We only assigned depth dependence to a single term $C_{\varphi N}(z)$, which we constructed as a Gaussian expression defined by a null background signal, a full width at half-maximum ($N^W$), and a maximum value ($N^{\max }$) at an assigned depth ($N^Z$). $C_{\varphi N}(z)$ alone controlled the vertical structure of our simulated profiles (Fig. 2).

 

Fig. 2. Example simulated profiles for $C_{\varphi }$, absorption coefficients at 443 nm, and scattering coefficients at 443 nm. The simulation parameters for these profiles are: $C_{\varphi D}=1\mathrm{mg}{\mathrm{m}}^{-3}$, $N^Z=5\mathrm{m}$, $N^W=5\mathrm{m}$, $N^{\max }=6\mathrm{mg}{\mathrm{m}}^{-3}$, $a_d(443)=0.02{\mathrm{m}}^{-1}$, $a_g(443)=0.05{\mathrm{m}}^{-1}$, and $b_d(555)=0.2{\mathrm{m}}^{-1}$ (Table 1). Panel (A) shows $C_{\varphi }$ and $C_{\varphi D}$ and indicates the three parameters that describe the Gaussian-shaped $C_{\varphi N}$ ($N^Z$, $N^W$, and $N^{\max }$). Panel (B) shows the corresponding $a_{\varphi N}(443)$, $a_{\varphi D}(443)$, $a_d(443)$, and $a_g(443)$. Panel (C) shows the corresponding $b_{\varphi N}(443)$, $b_{\varphi D}(443)$, and $b_d(443)$.

Download Full Size | PDF

Table 1. Values Used to Simulate Vertical IOP Profiles for Input in HE5

View Table | View all tables in this article

We carefully selected appropriate ranges of simulation parameters, based on extensive literature reviews and local knowledge from field sampling (Table 1). In addition to $a_{pg}(\lambda ,z)$ and $c_{pg}(\lambda ,z)$, running the HE5 “IOP Data” model also required generation of vertical profiles of $a_g(\lambda ,z)$ [see Eq. (10)] and $C_{\varphi }(z)$ $(=C_{\varphi N}(z)+C_{\varphi D}(z))$, plus spectra of $a_{\varphi }^{*}(\lambda )$. For the latter, we generated $a_{\varphi }^{*}(\lambda )$ for each simulation following:

Using all combinations of the values presented in Table 1 resulted in 1,920 simulated stations. We generated vertical profiles and spectra at 5 nm intervals from 400 to 750 nm to loosely mimic the output of a WET Labs, Inc. AC-S (a hyperspectral absorption and attenuation meter), which HE5 readily accepts as input into its “IOP Data” model. Figure 2 provides an example suite of profiles for one simulated station.

D. Modeling and Data Analysis

Using our simulated data as input into HE5, we generated simulated remote-sensing reflectances (Fig. 3) and depth profiles of spectral upwelling radiance and downwelling irradiance. We configured HE5 as follows: spectral output from 400 to 700 nm at 5 nm intervals; vertical output from 0 to 30 m at 3 m intervals; an infinitely deep bottom; cloud cover of 25%; wind speed of $5\mathrm{m}{\mathrm{s}}^{-1}$; solar geometry specific to year day 45, latitude 20°N, and longitude 60°E at local noon; the RADTRAN sky model; all inelastic scattering options enabled; and a spectrally constant and depth-independent $b_b/b$ of 0.01 [53]. Although some uncertainty accompanied our assignment of the latter, we repeated the analysis with other values and found our overall results to be unchanged.

 

Fig. 3. $R_{\mathrm{rs}}(\lambda )$ collected in situ (thick black lines) and synthesized using HE5 (colored thin lines). Panel (A) shows all synthesized and in situ $R_{\mathrm{rs}}(\lambda )$. Panel (B) shows only synthesized $R_{\mathrm{rs}}(\lambda )$ with optically weighted $C_{\varphi }$ ranging from 0.95 to $1.05\mathrm{mg}{\mathrm{m}}^{-3}$. Colors are only used to visually distinguish between spectra.

Download Full Size | PDF

We applied the ORM to each simulated $R_{\mathrm{rs}}(\lambda )$ using two different suites of input wavelengths with relevance to current and planned satellite ocean color science. First, we used only six visible MODISA wavelengths to represent the present state of satellite ocean color (412, 443, 488, 531, 547, and 667 nm). Second, we considered a suite of 16 visible wavelengths in consideration for future ocean color satellite instruments (the six MODISA wavelengths plus 400, 425, 460, 475, 510, 583, 617, 640, 655, and 665 nm) [15,54]. Hereafter, we refer to this expanded suite of wavelengths as “all wavelengths” in figures and tables. While we acknowledge the value in evaluating alternate combinations of wavelengths to study phytoplankton species composition [55], our analysis of these two suites also provides practical results for existing and forthcoming satellite instruments. These two $R_{\mathrm{rs}}(\lambda )$ suites serve well to explore what can be accomplished today (e.g., using MODISA) and what we might be able to accomplish in the near future {e.g., via the upcoming NASA Pre-Aerosols, Clouds, and Ocean Ecosystems (PACE) mission [15]}.

At this stage, our synthesized $C_{\varphi }(z)$ and $\mathrm{IOP}(z)$ profiles provide the “ground truth” for comparison with the modeled values from the ORM. Approximately 90% of what a satellite ocean color “sees” includes weighted contributions from all water column constituents shallower than the first optical depth ($z_{90}$; ${\mathrm{m}}^{-1}$); that is, the e-folding depth for diffuse attenuation coefficients for downwelling irradiance ($K_d(\lambda )$; ${\mathrm{m}}^{-1}$) [56]. To properly account for the subsurface $C_{\varphi N}(z)$ maxima, we optically weighted the $C_{\varphi }(z)$ and $\mathrm{IOP}(z)$s using the HE5-generated radiance and irradiance profiles, following the method of Zaneveld et al. [57], which produced the corresponding pseudo-depth-integrated values for these variables detected by the satellite instrument. We hereafter imply that values are optically weighted, unless depth dependence is specifically indicated [e.g., by “($z$)”]. As will be elaborated upon later, this weighting process highlights an inherent ambiguity within the ocean color paradigm; that is, identical $C_{\varphi }$ can be produced by larger, deeper blooms and smaller, shallower blooms, simply because of their relative vertical location in the water column.

3. Results

A. Synthesized ${\mathsf{R}}_{\mathsf{rs}}(\mathsf{\lambda })$

The synthesized $R_{\mathrm{rs}}(\lambda )$ reasonably reproduced the dynamic range of in situ values acquired in the Arabian Sea in March, 2011 (Fig. 3). Unusual $R_{\mathrm{rs}}(\lambda )$ can exist in any suite of synthesized values when input parameters combine in ways that do not occur naturally. We identified unnatural simulations for our study area as those with single scattering albedos at 443 nm less than 0.7 [occurring for stations with combinations of the lowest $C_{\varphi }$ and highest $a_{dg}(\lambda )$] or $R_{\mathrm{rs}}$ (490) greater than $0.008{\mathrm{sr}}^{-1}$ [occurring for stations with combinations of the lowest $C_{\varphi }$ and lowest $a_{dg}(\lambda )$]. We considered simulations with low single scattering albedos to be overly atypical, as previous radiative transfer and in situ studies indicated it exceeds 0.8 in most natural waters [58,59]. Similarly, we considered simulations with $R_{\mathrm{rs}}(490)>0.008{\mathrm{sr}}^{-1}$ to be overly unusual, as this threshold exceeds our in situ measurements by $>40\%$ and an extreme value for open ocean models by $>10\%$ [60]. Elimination of these unnatural simulations left 1437 stations for our analyses. As these remaining synthesized values fall within the envelope of spectra measured in the field, we propose that they adequately represent a natural range of environmental conditions for the northern Arabian Sea, at least for the purposes of this theoretical analysis.

B. Retrieval of Total IOPs

Retrievals of $b_{bp}(\lambda )$, $a_{pg}(\lambda )$, and $a_{\varphi }(\lambda )$ compared favorably with the ground-truth-synthesized values (Fig. 4, Table 2). It is clear by difference [$a_{pg}(443)-a_{\varphi }(443)$] that retrievals of $a_{dg}(443)$ also compared favorably. Using all available wavelengths yielded better results than only considering MODISA wavelengths. When using all wavelengths in the inversion, the comparisons maintained coefficients of determination ($r^2$), least squares regression slopes, and median ORM-to-synthesized ratios that ranged from 0.77 to 0.97, 0.65 to 0.95, and 0.85 to 1.21 respectively, indicating good performance over the full dynamic ranges of retrievals. The mean $r^2$, regression slopes, and ratios were 0.85, 0.83, and 1.07, respectively. The results changed somewhat when using only MODISA wavelengths, with $r^2$, regression slopes, and median ratios ranging from 0.67 to 0.94, 0.74 to 1.02, and 0.88 to 1.22, respectively. The corresponding mean $r^2$, regression slopes, and ratios degraded slightly to 0.82, 0.84, and 1.09. The most notable statistical difference between the two wavelength suites was an increase in root mean square error (RMSE) of 69%, 45%, and 35% for $b_{bp}(\lambda )$, $a_{pg}(\lambda )$, and $a_{\varphi }(\lambda )$, respectively, when using only MODISA wavelengths (Table 2).

 

Fig. 4. Comparison of ground-truth (synthesized) and ORM-derived $b_{bp}(443)$ (A) and (B), $a_{pg}(443)$ (C) and (D), and $a_{\varphi }(443)$ (E) and (F) using all available wavelengths in the inversion (left column) and only MODISA visible wavelengths in the inversion (right column). Colors show the numerical density of the retrievals, with red to purple indicating high to low volumes of sample sizes. We present comparative statistics in Table 2.

Download Full Size | PDF

Table 2. Ordinary Least Squares Regression Statistics for Ground-Truth-Synthesized Versus ORM-Derived IOPs^a

View Table | View all tables in this article

Our results suggest that ORM performance increases with finer $R_{\mathrm{rs}}(\lambda )$ spectral resolution, most notably with regard to retrieval variability (via RMSE). However, the remaining statistics ($r^2$, regression slopes, and ratios) suggest that data records from multispectral satellite instruments, such as MODISA, provide reliable ORM-derived IOPs on average for a significant dynamic range of water types. In general, these validation statistics fall well within the range of those presented in previous ORM studies [20,21]. To verify this, we re-ran our analyses using all available wavelengths as input into the quasi-analytical algorithm (QAA) of Lee [46]. The $r^2$, regression slopes, ratios, and RMSEs for QAA were: 0.93, 1.22, 1.18, and 0.0051 for $b_{bp}(443)$; 0.95, 1.23, 1.09, and 0.01617 for $a_{pg}(443)$; and, 0.93, 0.83, 0.97, and 0.01220 for $a_{\varphi }(443)$. More often than not, the performance of our ORM matched or exceeded that of QAA, which confirms to a first-order that our ORM yields results comparable with those of alternate configurations [21] (Table 2).

C. Optically Weighted ${\mathsf{C}}_{\varphi }$

The ORM demonstrated dependence on the depth of the subsurface maxima in its ability to estimate absolute magnitudes of N. miliaris. Before proceeding with our interpretation of this dependency, we feel it worthwhile to remind the reader how vertical structure in a water column constituent shapes its optically weighted value and the corresponding $R_{\mathrm{rs}}(\lambda )$ [38,39]. Multiple variants of $C_{\varphi N}(z)$ yielded common $C_{\varphi }$, the optically integrated values measured by an ocean color satellite instrument. Figure 5 presents $C_{\varphi }$ as a function of bloom depth ($N^Z$) and magnitude ($N^{\max }$) for a subset of our synthesized stations with constant bloom thickness ($N^W$), $C_{\varphi D}$, $a_d(443)$, $a_g(443)$, and $b_d(555)$. The deepest variants of $N_{\max }$ (3 and 6 mg ${\mathrm{m}}^{-3}$), for example, produced $C_{\varphi }$ that matched the shallower variants of $N_{\max }$ (0.5 and $1\mathrm{mg}{\mathrm{m}}^{-3}$) (Table 3). Generally speaking, $C_{\varphi }$ for N. miliaris blooms of all magnitudes converged to the background diatom signal as $N^Z$ increased. The weighted $C_{\varphi }$ increasingly underestimated $N^{\max }$ as the subsurface maximum increased in magnitude and deepened. For the smallest $N^{\max }$, $C_{\varphi }$ depended most significantly on the background $C_{\varphi D}$ and showed very little dependence on the depth and thickness of the subsurface N. miliaris maxima (Table 3). Ultimately, the magnitude of the peak subsurface population exceeded its remote-sensing value (e.g., $N^{\max }>C_{\varphi }$) for most examples.

 

Fig. 5. Optically weighted $C_{\varphi }$ and corresponding $R_{\mathrm{rs}}(\lambda )$ as a function of N. miliaris bloom depth ($N^Z$) and magnitude ($N^{\max }$). This subset of simulations has constant $C_{\varphi D}$ ($=0.5\mathrm{mg}{\mathrm{m}}^{-3}$), $N^W$ $(=5\mathrm{m})$, $a_g(443)$ $(=0.005{\mathrm{m}}^{-1})$, $a_d(443)$ $(=0.002{\mathrm{m}}^{-1})$, and $b_p(555)$ $(=0.1{\mathrm{m}}^{-1})$. Panel (A) shows optically weighted $C_{\varphi }$ versus bloom depth, with colors representing different bloom magnitudes. Panel (B) shows $R_{\mathrm{rs}}(\lambda )$ for a bloom magnitude of $0.5\mathrm{mg}{\mathrm{m}}^{-3}$, with different line styles representing different bloom depths. Panels (C)–(E) follow Panel (B), but for bloom magnitudes of 1, 3, and $6\mathrm{mg}{\mathrm{m}}^{-3}$. Colors in Panels (B)–(E) follow Panel (A) for clarity.

Download Full Size | PDF

Table 3. Optically weighted ${\mathsf{C}}_{\varphi }$, ${\mathsf{K}}_{\mathsf{d}}(\mathsf{490})$, and ${\mathsf{z}}_{\mathsf{90}}$ for the Example Simulations Presented in Fig. 5

View Table | View all tables in this article

Changing the vertical position of a subsurface bloom of fixed size and thickness yielded different $C_{\varphi }$. Consider, for example, the profiles with $N_{\max }=6\mathrm{mg}{\mathrm{m}}^{-3}$ in Figure 5(A) (red circles). The corresponding optically weighted $C_{\varphi }$ dropped from 5.58 to $0.75\mathrm{mg}{\mathrm{m}}^{-3}$ as this synthesized bloom progressed from the surface to a depth of 15 m (Table 3). Given that standard ocean color algorithms are tuned to optically weighted constituent stocks [18,61], it stands to reason that satellite retrievals of these stocks do not independently provide information about the vertical distribution of the constituents or the thicknesses of their layers. To circumvent this, previous studies used external hydrographic and environmental information to make assumptions about vertical distributions [27,62]. Uitz et al. [27], for example, suggested that the vertical structure of $C_{\varphi }$$(z)$ could be inferred on average in the open ocean from its absolute magnitude and some knowledge of the local average mixed layer depth. Not surprisingly, the rate of change of $C_{\varphi }$ with bloom depth ($N^Z$) varied strongly with the magnitude of the subsurface maximum ($N^{\max }$).

Vertically displacing the subsurface N. miliaris maxima resulted in different $R_{\mathrm{rs}}(\lambda )$ [Figs. 5(B)–5(E)]. The depths of subsurface features and the bulk attenuating (absorbing and scattering) properties of the water mass modulate the average contributions of each constituent to $R_{\mathrm{rs}}(\lambda )$ [56,57]. Consider again the $N^{\max }=6\mathrm{mg}{\mathrm{m}}^{-3}$ in Fig. 5 (red circles). The color of the water appeared green under a near-surface bloom and progressively bluer as this bloom sank deeper into the water column [Fig. 5(E)]. All $N^{\max }$ demonstrated this, but the effect lessened as it decreased. For the case study presented in Fig. 5, $R_{\mathrm{rs}}(\lambda )$ appeared similar for the multiple scenarios with $C_{\varphi }$ between 0.6 and $0.75\mathrm{mg}{\mathrm{m}}^{-3}$ (compare, e.g., the dashed lines in Figs. 5(B) and 5(C) with the dashed–dotted lines in Figs. 5(D) and 5(E). Within the full population of simulations, however, $R_{\mathrm{rs}}(\lambda )$ differed for many instances of common $C_{\varphi }$. Fig. 3(B) shows $R_{\mathrm{rs}}(\lambda )$ for simulations with $C_{\varphi }$ near $\sim 1\mathrm{mg}{\mathrm{m}}^{-3}$, many of which vary significantly in the blue–green part of the spectrum. This variability in $R_{\mathrm{rs}}(\lambda )$ for common $C_{\varphi }$ contributes to the overall variability seen in algorithm development data sets used to develop common bio-optical models (see, e.g., the vertical scatter of the radiometric data for each given optically weighted $C_{\varphi }$ presented in Figs. 3 and 6 of O’Reilly [18] and Fig. 5(A) of Werdell and Bailey [61]).

D. Retrieval of Multiple Phytoplankton Species

The ORM demonstrated decent skill in identifying N. miliaris and diatoms, but with sufficient variability to indicate that the inversion cannot provide absolute magnitudes for both species at all times. All ORMs most effectively retrieve total absorption (e.g., [20,21]). Deconstructing total absorption into its $\mathrm{CDOM}+\mathrm{NAP}$ and phytoplankton components remains more uncertain (Table 2) and, not surprisingly, decomposing the total phytoplankton signature into multiple components carries additional uncertainties [23,24,29]. In general, the ORM tended to underestimate N. miliaris and to overestimate diatoms. When using all wavelengths in the inversion, the ORM-to-synthesized bias, RMSE, and regression slope changed from $0.093\mathrm{mg}{\mathrm{m}}^{-3}$, $0.261\mathrm{mg}{\mathrm{m}}^{-3}$, and 0.82 for total $C_{\varphi }$ to $0.097\mathrm{mg}{\mathrm{m}}^{-3}$, $0.215\mathrm{mg}{\mathrm{m}}^{-3}$, and 0.93 for $C_{\varphi D}$ and $-0.005\mathrm{mg}{\mathrm{m}}^{-3}$, $0.310\mathrm{mg}{\mathrm{m}}^{-3}$, and 0.64 for $C_{\varphi N}$ (Fig. 6). We calculated bias as the average difference between modeled and measured values. The differences amplified when using only MODISA wavelengths, with biases, RMSEs, and slopes changing from $0.094\mathrm{mg}{\mathrm{m}}^{-3}$, $0.313\mathrm{mg}{\mathrm{m}}^{-3}$, and 0.82 for $C_{\varphi }$ to $0.113\mathrm{mg}{\mathrm{m}}^{-3}$, $0.402\mathrm{mg}{\mathrm{m}}^{-3}$, and 1.02 for $C_{\varphi D}$ and $-0.019\mathrm{mg}{\mathrm{m}}^{-3}$, $0.345\mathrm{mg}{\mathrm{m}}^{-3}$, and 0.55 for $C_{\varphi N}$. The ORM exclusively underestimated N. miliaris at concentrations greater than $2\mathrm{mg}{\mathrm{m}}^{-3}$.

 

Fig. 6. Comparison of ground-truth (synthesized) and ORM-derived $C_{\varphi }$ (A) and (B), $C_{\varphi D}$ (C) and (D), and $C_{\varphi N}$ (E) and (F) using all available wavelengths in the inversion (left column) and only MODISA visible wavelengths in the inversion (right column). Colors as in Fig. 4.

Download Full Size | PDF

The ORM frequently identified monospecific phytoplankton populations, where only N. miliaris or only diatoms were present. For the former, the ORM routinely detected N. miliaris, but often also returned false-positives for diatoms [Figs. 7(A) and 7(B)]. In general, the ORM underestimated $C_{\varphi N}$ (most significantly when the N. miliaris bloom resided at the surface), which resulted in an overestimation of $C_{\varphi D}$. The lower cluster of $C_{\varphi N}\ge 2\mathrm{mg}{\mathrm{m}}^{-3}$ in Figs. 7(A) and 7(B) corresponds to the $C_{\varphi D}\ge 0.5\mathrm{mg}{\mathrm{m}}^{-3}$ in these panels. Reducing the spectral resolution of $R_{\mathrm{rs}}(\lambda )$ to MODISA wavelengths amplified the underestimation of $C_{\varphi N}$ for concentrations greater than $4\mathrm{mg}{\mathrm{m}}^{-3}$ and increased the magnitudes of false positives for diatoms. When only diatoms were present, however, the ORM never significantly detected N. miliaris (always $<0.2\mathrm{mg}{\mathrm{m}}^{-3}$) and successfully estimated $C_{\varphi D}$ over its full dynamic range [Figs. 7(C) and 7(D)]. Using only MODISA wavelengths did not significantly degrade these results. In general, the ORM did not identify N. miliaris when it was absent and often identified N. miliaris when it was present. Although the quantification of $C_{\varphi N}$ and $C_{\varphi D}$ may be imperfect, this suggests that a positive ORM-derived $C_{\varphi N}$ reliably indicates the presence of N. miliaris somewhere in the upper water column.

 

Fig. 7. Comparison of ground-truth (synthesized) and ORM-derived $C_{\varphi D}$ and $C_{\varphi N}$ for a mono-species subset of simulations using all available wavelengths in the inversion (left column) and only MODISA visible wavelengths in the inversion (right column). We considered only synthesized $C_{\varphi D}=0.02$ and $C_{\varphi N}\ge 0.5\mathrm{mg}{\mathrm{m}}^{-3}$ in panels (A) and (B). We considered only synthesized $C_{\varphi D}\ge 0.1$ and $C_{\varphi N}=0\mathrm{mg}{\mathrm{m}}^{-3}$ in panels (C) and (D). Crosses and filled circles show $C_{\varphi N}$ and $C_{\varphi D}$, respectively.

Download Full Size | PDF

Equally mixed populations of N. miliaris and diatoms produced $R_{\mathrm{rs}}(\lambda )$ that occasionally challenged the ORM. Consider the following scenarios presented in Fig. 8: monospecific $C_{\varphi D}=0.5\mathrm{mg}{\mathrm{m}}^{-3}$, monospecific $C_{\varphi N}=0.5\mathrm{mg}{\mathrm{m}}^{-3}$, monospecific $C_{\varphi D}=1\mathrm{mg}{\mathrm{m}}^{-3}$, monospecific $C_{\varphi N}=1\mathrm{mg}{\mathrm{m}}^{-3}$, and mixed $C_{\varphi D}=C_{\varphi N}=0.5\mathrm{mg}{\mathrm{m}}^{-3}$ (total $C_{\varphi }=1\mathrm{mg}{\mathrm{m}}^{-3}$). The four monospecific populations produced sufficiently unique $R_{\mathrm{rs}}(\lambda )$ such that the ORM effectively estimated $C_{\varphi }$, $C_{\varphi N}$, and $C_{\varphi D}$ without false positives for the absent species. In the mixed population, where N. miliaris and diatoms both contributed equally to $C_{\varphi }$, the ORM achieved accurate $C_{\varphi }$, but less effectively separated the two phytoplankton populations. The ORM attributed much of the $C_{\varphi }$ to $C_{\varphi D}$ ($0.75\mathrm{mg}{\mathrm{m}}^{-3}$), which is consistent with our previous general observation that the ORM tends to slightly over- and underestimate diatoms and N. miliaris, respectively. Equal contributions by N. miliaris and diatoms did not produce $R_{\mathrm{rs}}(\lambda )$ that fell equally between the two monospecific spectra, suggesting to a first-order that the contributions of the two species to $R_{\mathrm{rs}}(\lambda )$ were not strictly additive. The case studies with equivalent $C_{\varphi }$ presented in Figs. 5(B) and 5(C) show a clean progression in $R_{\mathrm{rs}}(\lambda )$ as a bloom deepened in the water column and the color of the water transitioned from green to blue. In contrast, the spectra in Fig. 8(B) do not demonstrate a clean transition as species mix.

 

Fig. 8. Five example simulated $R_{\mathrm{rs}}(\lambda )$ and their corresponding ORM-derived $C_{\varphi }$. Each simulation was created using $a_d(443)=0.002$, $a_g(443)=0.005$, and $b_d(555)=0.1{\mathrm{m}}^{-1}$. Panel (A) presents spectra for $C_{\varphi D}\sim 0.5\mathrm{mg}{\mathrm{m}}^{-3}$ without N. miliaris (black) and $C_{\varphi N}\sim 0.5\mathrm{mg}{\mathrm{m}}^{-3}$ without diatoms (orange). Panel (B) presents spectra for $C_{\varphi D}\sim 1\mathrm{mg}{\mathrm{m}}^{-3}$ without N. miliaris (blue), $C_{\varphi N}\sim 1\mathrm{mg}{\mathrm{m}}^{-3}$ without diatoms (green), and a mixed population with $C_{\varphi D}$ and $C_{\varphi N}\sim 0.5\mathrm{mg}{\mathrm{m}}^{-3}$ each (red). Panel (C) presents ground-truth (synthesized) versus ORM comparisons, with colors referring to the corresponding $R_{\mathrm{rs}}(\lambda )$, and circles, crosses, and asterisks indicating $C_{\varphi }$, $C_{\varphi D}$, and $C_{\varphi N}$, respectively.

Download Full Size | PDF

The accuracy of the ORM showed dependency on the depth of the N. miliaris bloom when considering all available wavelengths and only MODISA wavelengths. On average, this dependency did not vary with the total magnitude of $C_{\varphi }$ (Table 4). For the full population of synthesized data, the ORM consistently under- and over-estimated $C_{\varphi N}$ and $C_{\varphi D}$, respectively, when $N^Z=1\mathrm{m}$ (Fig. 9). This under- and overestimation of $C_{\varphi N}$ and $C_{\varphi D}$ generally reversed when $N^Z$ deepened. In all cases, the differences in $C_{\varphi N}$ exceeded those for $C_{\varphi D}$. With one exception (the $C_{\varphi }>1\mathrm{mg}{\mathrm{m}}^{-3}$ subset at 15 m), the largest positive differences in $C_{\varphi N}$ appeared at 5 m and decreased as the subsurface maxima deepened. Overall, the ORM overestimated N. miliaris by 23.8% to 77.0%, on average, for $N^Z\ge 5\mathrm{m}$. Absolute biases in $C_{\varphi N}$ decreased with increasing bloom depth, which is not surprising given that the (optically weighted) total $C_{\varphi }$ also decreased with increasing bloom depth (see, e.g., Fig. 5). In contrast to N. miliaris, the largest differences in $C_{\varphi D}$ occurred near the surface. At depth, ORM retrievals of $C_{\varphi D}$ differed from ground truth by $<8\%$, with only one exception (the $C_{\varphi }<1\mathrm{mg}{\mathrm{m}}^{-3}$ subset at 15 m), and maintained null absolute biases ($<0.07\mathrm{mg}{\mathrm{m}}^{-3}$).

Table 4. Median Relative Percent Differences (MPD) and Absolute Biases between the Ground-Truth-Synthesized and ORM-Derived ${\mathsf{C}}_{\varphi \mathsf{N}}$ and ${\mathsf{C}}_{\varphi \mathsf{D}}$^a

View Table | View all tables in this article

4. Discussion

A. Interpreting the ORM Retrievals

Ocean color satellite instruments provide rich data streams for studying phytoplankton dynamics. The community considers these data streams to be sufficiently critical to require that upcoming ocean color satellite programs include capabilities to discriminate between phytoplankton groups in support of advanced carbon studies [14,15]. Our study focused on a class of algorithm known as ORMs (also known as semi-analytic algorithms). Few studies have demonstrated the ability of an ORM to successfully deconvolve $a_{\varphi }(\lambda )$ into contributions by individual phytoplankton groups over a large range of environmental conditions [37]. In practice, ORMs are currently applied to satellite measurements of $R_{\mathrm{rs}}(\lambda )$ without exact a priori knowledge of phytoplankton physiological state or environmental conditions for a given pixel. They operate assuming that their adopted eigenvectors adequately represent local conditions at the moment the measurement is made. We constructed a controlled environment where the eigenvectors used to parameterize the ORM were the same as those used to build the input $R_{\mathrm{rs}}(\lambda )$. Doing so generated scenarios for which we had exact (simulated) knowledge of phytoplankton physiological state and environmental conditions. Comparisons of ORM retrievals and simulation inputs within this ideal environment should reveal weaknesses in ORM operation and interpretation that merit subsequent review. To ensure broad applicability of our results, we conducted our analyses using a very common ORM algorithmic form (built using the GIOP framework [21] and analogous in form to that of, e.g., Maritorena et al. [19], Roesler et al. [35], and Westberry et al. [37]).

Our results reaffirm previous conclusions that ORMs effectively separate the total signal of the collective phytoplankton community from other water column constituents, such as CDOM and NAP [20,21] (Fig. 4, Table 2). Our results suggest, however, that the absolute accuracy of ORM retrievals degrades when further deconstructing the derived $C_{\varphi }$ and IOPs into additional subcomponents, even under perfectly controlled conditions (Fig. 6). Despite this, we view our results as positive. Although the absolute magnitudes of $C_{\varphi N}$ and $C_{\varphi D}$ maintained biases, the ORM successfully detected whether or not N. miliaris appeared in the simulated water column for both the limited (MODISA) and expanded (PACE-like) wavelength suites (Figs. 6–8). This quantitatively calls for caution when interpreting the absolute magnitudes of the retrievals, but qualitatively suggests that the ORM provides a robust mechanism for identifying the presence or absence of N. miliaris.

Let us briefly explore this suggestion further. A positive retrieval of $C_{\varphi N}$ reliably indicated the presence of N. miliaris at some near surface ($\le 15\mathrm{m}$) position in our simulations. The ORM effectively identified diatoms and often identified N. miliaris in mono- and mixed-populations, but often with imperfect absolute magnitudes (e.g., $C_{\varphi N}$ showed negative biases when the bloom maxima resided near the surface and positive biases when it deepened). Despite this, N. miliaris was always present when the ORM retrieved positive $C_{\varphi N}$. The ORM did fail to identify N. miliaris on occasion and produced corresponding false positives for diatoms, but it never substantially identified N. miliaris for simulations with only diatoms. In other words, when N. miliaris was present, the ORM often (but not always) estimated $C_{\varphi N}>0\mathrm{mg}{\mathrm{m}}^{-3}$, but when N. miliaris was not present, the ORM never estimated $C_{\varphi N}\gg 0\mathrm{mg}{\mathrm{m}}^{-3}$. At a minimum then, this suggests that the ORM can be used to indicate whether or not N. miliaris is present in a satellite pixel. Indeed, we explored the qualitative use of satellite-derived $C_{\varphi N}$ in a companion study via application of the ORM to MODISA to simply identify and catalog pixels with substantially positive $C_{\varphi N}$. To a first-order, applying this (binary absence or presence) identification scheme to the full MODISA data record permitted us to study the spatial and temporal distribution of N. miliaris appearances over a decade in the northern Arabian Sea [63].

B. Phytoplankton Absorption Eigenvectors

Both spectral resolutions yielded comparable regression statistics for the retrieved $b_{bp}(\lambda )$, $a_{pg}(\lambda )$, and $a_{\varphi }(\lambda )$, as evidenced by $r^2$, regression slopes, and median ORM-to-synthesized ratios (Fig. 4, Table 2). Using all available wavelengths enabled somewhat improved separation of $C_{\varphi D}$ and $C_{\varphi N}$ (Fig. 6). The variability in all IOP and $C_{\varphi }$ retrievals (not just $C_{\varphi D}$ and $C_{\varphi N}$), however, increased significantly when using $R_{\mathrm{rs}}(\lambda )$ with reduced spectral resolution (see, e.g., the RMSE in Table 2). The slopes from 400–440 nm and from 440–510 nm differ substantially for $a_{\varphi D}^{*}(\lambda )$ and $a_{\varphi N}^{*}(\lambda )$ (Fig. 1). The MODISA wavelength suite maintains less spectral information in these spectral ranges (and, thus, fewer degrees of freedom in the inversion) with which to discriminate between N. miliaris and diatoms. In other words, the optical signatures of $a_{\varphi D}^{*}(\lambda )$ and $a_{\varphi N}^{*}(\lambda )$ contrast less when viewed at only MODISA wavelengths versus at PACE-like wavelengths. This has implications for planning new satellite missions with mandates to deliver data records with improved capabilities for discriminating between phytoplankton species [14,15]. Certainly, the SeaWiFS and MODISA data records (and others) stimulated innovation and progress toward the remote identification of phytoplankton groups [22–31,34,36,37,40]. Given our results, however, we expect that new satellite instruments with finer spectral resolution will enable more quantitative resolution of multiple species with fewer uncertainties [32,35] (Figs. 6 and 7).

Single $a_{\varphi }^{*}(\lambda )$ cannot perfectly represent any phytoplankton species under all conditions of growth, nutrient availability, and ambient light [29,64]. This presents problems for ocean color remote sensing, as the phytoplankton population under observation and its physiological state cannot be known a priori. Previous studies adopted varied methods of expressing phytoplankton absorption: (1) single $a_{\varphi }^{*}(\lambda )$ derived from global assemblies of in situ measurements [19,46]; (2) multiple $a_{\varphi }^{*}(\lambda )$ to simultaneously represent several phytoplankton size classes [23,24,29,35]; (3) dynamic assignment of $a_{\varphi }^{*}(\lambda )$ based on estimates of $C_{\varphi }$ [21]; and (4) iterations on ranges of $a_{\varphi }^{*}(\lambda )$ [65,66]. Our ORM generically falls into method (2), but might realize future benefits from being refined to the forms of (3) or (4). Additional in situ measurements from various stages of diatom and N. miliaris blooms during the NEM could be used to call the ORM in an iterative fashion using ranges of observed eigenvectors. More sophisticated parameterizations of $a_{\varphi }(\lambda )$ might also be developed to replace the linear expression we adopted $(=C_{\varphi }{a}_{\varphi }^{*}(\lambda ))$ [64,67]. Nevertheless, these approaches all require some a priori assumption of phytoplankton spectral shape(s). The need to make such an assumption argues in favor of retrospective analysis when interpreting satellite time-series of $C_{\varphi }$ and $a_{\varphi }(\lambda )$ to verify that local hydrodynamics, environmental conditions, and biology support the a priori selection of eigenvectors.

With respect to the mechanics of our analyses, our choice in $a_{\varphi }^{*}(\lambda )$ remains irrelevant. However, they merit additional discussion in the context of interpretation of satellite data records. We generated our simulated IOP profiles and configured the ORM using consistent $a_{\varphi }^{*}(\lambda )$. Despite this, the ORM retrievals did not always exactly reproduce the inputs $C_{\varphi N}$ and $C_{\varphi D}$. We chose diatoms and N. miliaris because they best represent significantly different end-members in the Arabian Sea, and because the ability to detect N. miliaris using long-term satellite data records presents an opportunity to quantify its emerging presence over the past decade. We know, however, that additional phytoplankton species coexist and flourish in the Arabian Sea, such as dinoflagellates and Trichodesmium spp. Unfortunately, the ORM became unstable when we added a third phytoplankton component (raising the unknowns to five), particularly when limiting the runs to MODISA wavelengths. In the two-component system, however, certain combinations of our $a_{\varphi N}^{*}(\lambda )$ and $a_{\varphi D}^{*}(\lambda )$ yielded spectra that resembled that of dinoflagellates (see Thibodeau et al. [7]). Following, in scenarios where $a_{\varphi }^{*}(\lambda )$ within an ORM do not represent the population under observation, the ORM defaults to combinations of its native phytoplankton components and misidentifies the species present. This supports our recommendation of caution when interpreting the absolute magnitudes of multiple phytoplankton eigenvalues. It also argues that, for targeted application of ORMs, a configuration for alternate phytoplankton species would produce more reliable and appropriate IOPs outside of the NEM, when we would not expect an abundance of N. miliaris.

Our choice in $a_{\varphi D}^{*}(\lambda )$ and $a_{\varphi N}^{*}(\lambda )$ contributed in part to the vertical biases seen in $C_{\varphi N}$. Concentrations of absorbing constituents control the shape of $R_{\mathrm{rs}}(\lambda )$, whereas the magnitude of scattering relative to absorption modulates the brightness. Near the surface, pure seawater absorption dominates the total absorption signal at wavelengths 600 nm [51]. As CDOM and NAP contribute very little to total absorption in this spectral range [48], phytoplankton begin to shape $R_{\mathrm{rs}}(\lambda )$ in the absence of these other water column constituents as blooms near the surface. Blooms of sufficient magnitude also enhance backscattering and can notably fluoresce [68], which elevates red $R_{\mathrm{rs}}(\lambda )$, as demonstrated in Figs. 5(D) and 5(E). In these figures, the ratio of $R_{\mathrm{rs}}(443)$-to-$R_{\mathrm{rs}}(670)$ decreases toward unity as the simulated bloom maxima shallows. The $a_{\varphi D}^{*}(\lambda )$ and $a_{\varphi N}^{*}(\lambda )$ we adopted also maintain substantially different blue-to-red ratios, with diatoms more closely approaching a ratio of one (Fig. 1). We expect that these near-unity blue-to-red relationships favored the retrieval of $C_{\varphi D}$ when $N^Z=1\mathrm{m}$. This favoritism switched to $C_{\varphi N}$ as the bloom maxima deepened and the red $R_{\mathrm{rs}}(\lambda )$ became water-dominated instead of phytoplankton-dominated.

Several of our processing choices complicate this interpretation of vertical biases. First, as clearly demonstrated near 683 nm in Figs. 5(D) and 5(E), we enabled the fluorescence capabilities of HE5. A discussion of the fluorescence model incorporated into HE5 exceeds the scope of this manuscript. In contrast to HE5, however, our ORM (as well as most others) does not account for phytoplankton fluorescence [20,21]. To our knowledge, most (quasi-single scattering) relationships between apparent and IOPs, e.g., the coefficients in Eq. (2), or simplify contributions by fluorescence to $R_{\mathrm{rs}}(\lambda )$. Second, the Levenberg–Marquardt inversion scheme employs a chi-squared cost function that minimizes relative differences between the fit and measured $R_{\mathrm{rs}}(\lambda )$ with equal consideration of all wavelengths [21]. Given the weakness of the red $R_{\mathrm{rs}}(\lambda )$ signal in most open ocean waters, small absolute differences between fit and measured spectra amplify into substantial relative differences, which can overweight this region of the spectrum. Previous studies partially addressed this by constraining their inversion to 412–555 nm [19,65]. Doing so using our simulations increased the variability of our bulk-derived IOPs. For example, the RMSEs for $b_{bp}(443)$ and $a_{pg}(443)$ rose from 0.00031 to 0.00056 and 0.0125 to 0.0167, respectively. We recommend that subsequent studies explore revised cost functions that consider both absolute and relative goodness-of-fits [15] and/or alternate spectral weighting schemes (see e.g., the discussion of this in Werdell et al. [21]).

 

Fig. 9. Comparisons of ground-truth (synthesized) and ORM-derived $C_{\varphi N}$ (left column) and $C_{\varphi D}$ (center column) using all available wavelengths in the inversion. Results are stratified by depth of the N. miliaris subsurface maxima, with $N^Z=1$, 5, 10, and 15 m shown from top to bottom. The right column presents frequency distributions of relative percent differences, calculated as 100%* (ORM/synthesized-1). Solid and dotted lines indicate differences for $C_{\varphi N}$ and $C_{\varphi D}$, respectively.

Download Full Size | PDF

C. Future Directions

We conclude with a discussion of several general paths forward: (1) refining the ORM parameterization and evaluating alternate configurations (e.g., following Werdell et al. [21]); (2) using other ORM products, namely $b_{bp}(\lambda )$, to help detect N. miliaris; and (3) re-exploring the role of vertical structure in interpretations of ocean color data records. With regard to the ORM itself, we reaffirmed that a GIOP-like form can produce accurate estimates of bulk IOPs [$b_{bp}(\lambda )$, $a_{dg}(\lambda )$, and $a_{\varphi }(\lambda )$] as has been repeatedly demonstrated (e.g., [19–21,35,37]). The estimates showed biases, whereas our ORM effectively discriminated between $C_{\varphi }$ from diatoms and N. miliaris. With regard to improving the latter, reconsidering our treatment of phytoplankton absorption remains an obvious next step as discussed in Section 4.B. Other ORM parameterizations could also be refined, several of which we briefly explored. Using linear matrix inversion [65,66] in lieu of the nonlinear Levenberg–Marquardt inversion scheme yielded almost identical results, but with fewer overall viable retrievals. Alternate forms of $b_{bp}^{*}(\lambda )$ and $a_{dg}^{*}(\lambda )$ could also be adopted [31,46]. That said, ORM-retrieved $C_{\varphi N}$ demonstrated no significant dependence on the magnitudes of $a_g(\lambda )$, $a_d(\lambda )$, or $b_d(\lambda )$ (Fig. 10). Higher $a_g(443)$ and $a_d(443)$ corresponded to slightly elevated ORM-derived $C_{\varphi N}$, but these differences fell well below variances seen elsewhere in our analyses, such as by $N^Z$. We expect to pursue incorporating pixel classification [69] and ensemble (iterative or bootstrapping) [65,66] schemes into future studies, both of which provide vehicles for considering ranges of eigenvectors within the ORM.

 

Fig. 10. Comparisons of ground-truth (synthesized) and ORM-derived $C_{\varphi N}$ using all available wavelengths in the inversion. Results are stratified by the varied $a_g(443)$ (left), $a_d(443)$ (center), and $b_d(555)$ (right) used to generate the synthesized IOP profiles and $R_{\mathrm{rs}}(\lambda )$.

Download Full Size | PDF

Our analyses attempted to discriminate between diatoms and N. miliaris using phytoplankton absorption coefficients and $C_{\varphi }$. However, differences in diatom and N. miliaris backscattering properties could also be exploited to identify these species. Changes in the $b_{bp}(\lambda )$-to-$C_{\varphi }$ ratio over time, for example, would indicate changes in the optics (e.g., average particle size) or biology (e.g., phytoplankton or nonpigmented grazers). Figure 11 presents the relationship between $b_{bp}(443)$ and $C_{\varphi }$ as a function of the abundance of N. miliaris for both the simulated data set (panel A) and the ORM (panel B). The linear pattern of $b_{bp}(443)$ versus $C_{\varphi }$ results from our parameterization of the simulated profiles (Table 1). Two features stand out. First, the simulated and ORM relationships exhibit consistent patterns in $C_{\varphi N}$, dynamic ranges, and slopes. The simulated relationships represent a forward solution for the exact radiative transfer equations (RTE) while the ORM relationships represent an inversion solution for a simplified (single scattering) approximation of the RTE. The similarities in the patterns suggest that the ORM provides an adequate approximation of the RTE, with the smearing of points [relative to the clean slopes in Fig. 11(A)] resulting from the approximations made. Second, $b_{bp}(443)$ changes with $C_{\varphi N}/C_{\varphi }$; although the direction of the change appears to reverse several times as $C_{\varphi }$ rises from 0 to $6\mathrm{mg}{\mathrm{m}}^{-3}$. In practice, such a pattern could provide a useful indicator of N. miliaris if robust patterns in $b_{bp}(\lambda )$ versus $C_{\varphi N}/C_{\varphi }$ emerge. However, as no consensus covariance emerged from our synthesized data set, similar analyses of in situ data (as they become available) are warranted.

 

Fig. 11. $C_{\varphi }$ versus $b_{bp}(443)$. Panel (A) shows $C_{\varphi }$ from our simulated profiles versus backscattering from HE5. Panel (B) shows $C_{\varphi }$ from the ORM versus backscattering from the ORM. Colors indicate the magnitude of $C_{\varphi N}$ in mg ${\mathrm{m}}^{-3}$.

Download Full Size | PDF

Synthesized data cannot represent all conditions at all times. We could endlessly refine the parameterizations used to generate the simulated IOP profiles (Table 1); however, the values we selected correspond closely to field measurements and produce complementary $R_{\mathrm{rs}}(\lambda )$ to those measured in situ (Fig. 2) [7,40]. Regarding the role of scattering, we assumed that the N. miliaris and diatom species present during the NEM are large relative to visible wavelengths (thus, our choice of $\eta =0$) [70,71]. Optically, N. miliaris remains poorly studied relative to other phytoplankton [49]; however, preliminary studies suggest inefficient light scattering, consistent with its cell-specific properties [40] (thus, our choice of $b^{*}(555)=0.15{\mathrm{m}}^2{\mathrm{mg}}^{-1}$). Although our constant $b_b(\lambda )/b(\lambda )$ of 0.01 remains a vetted global average, we also acknowledge that this ratio ranges from 0.005 to 0.02 over a range of oceanographic conditions [53,72]. Currently, we have no knowledge of the scattering phase function of N. miliaris, but expect it might differ from that embedded in the Monte Carlo simulations used to generate the apparent-to-IOP relationship expressed in Eq. (2) [47]. An optical closure study following Tzortziou et al. [73] would be prudent as subsequent work to better constrain $b_b(\lambda )/b(\lambda )$ (and, thus, the scattering phase function) of species present during the NEM.

Finally, it remains well-known within the oceanographic community that ocean color estimates of water column properties represent optically weighted contributions that do not independently provide information on vertical structure [27,56,62]. While we did not originally intend to revisit this, our simulations provided an opportunity to pose a reminder and offer two related suggestions with potential for advancing our collective ability to remotely sense phytoplankton community composition. First, if the community wants true estimates of vertical structure from space, alternate technologies need to be developed [74] and/or methods that make use of complementary environmental information need to be nurtured [27,62]. The case studies presented in Fig. 5 demonstrate that optically weighted values cannot independently describe the vertical position, subsurface magnitude and thickness, or age of a bloom. Furthermore, our results provided little evidence that the ORM better detects N. miliaris as its bloom expands in abundance or thickness. While the known biases with $N^Z$ were quantifiable in this controlled study, the vertical position of N. miliaris detected from space will not be known a priori and cannot currently be independently inferred in all water masses at all times using passive ocean color radiometry [27,62]. In practice, this simply means that absolute abundances for multiple ORM-derive phytoplankton communities can only be compared with the understanding that they represent pseudo-depth-integrated values.

Second, reconstructing the optically weighted $C_{\varphi }$ seen by a satellite instrument requires collecting biogeochemical and radiometric data with sufficient vertical resolution to capture heterogeneity over the layer defined by $z_{90}$ [57]. The vertical composition and structure of the water column shapes what the satellite measures and, therefore, drives how in situ measurements should be collected and prepared for use in both ORM development (parameterization) and ocean color satellite data product validation. For a heterogeneous water column with deep $C_{\varphi }$ maxima appearing within the first e-folding depth, near surface measurements alone cannot be accurately compared with satellite-derived data products [61]. In practice, this argues for specialized in situ data collection and processing to support ocean color satellite validation activities.

5. Conclusions

We pursued this work for two primary reasons: (1) to contribute to the growing foundation upon which advanced methods for remote phytoplankton detection are being built and (2) to begin developing a method to identify N. miliaris in the northern Arabian Sea. We explored the capabilities of an ORM to discriminate between two distinct phytoplankton communities under a diverse array of biophysical conditions. As in previous studies, the ORM successfully separated the total phytoplankton signal from those of other water column constituents. The ORM effectively separated the individual contributions of N. miliaris and diatoms; however, the absolute estimates showed biases even under our perfectly controlled conditions. The vertical structure of our synthesized blooms highly influenced these biases, analysis of which reaffirmed that ORM retrievals alone cannot provide information on the vertical structure of a phytoplankton bloom. In the end, our results quantitatively call for caution when interpreting the absolute magnitudes of the retrievals, but qualitatively suggest that the ORM provides a robust mechanism for identifying the presence or absence of N. miliaris. Not surprisingly, incorporating $R_{\mathrm{rs}}(\lambda )$ at additional wavelengths improved the quality of the ORM retrievals, underscoring the benefit of additional spectral information on forthcoming satellite instruments with mandates to produce phytoplankton community data products. In the short-term, we propose that our ORM (and other ORMs of this form) can adequately support studies that require only identifications of phytoplankton groups. In the long-term, we anticipate the quantitative skills of ORMs to improve as additional technologies, relevant in situ data, and environmental data get routinely included in modeling and data analysis activities.