For the sensitivity analyses, we generated a synthesized dataset of coincident R_rs(λ), a_dg(λ), a_ph(λ), b_bp(λ), and C_a that closely mimics the IOCCG Ocean Color Algorithms Working Group synthetic dataset (IOCCG, 2006), but with expanded and finer spectral resolution. Appendix A describes its development (see the online article page to access the Appendix). This dataset maintains a spectral resolution of 1 nm intervals from 350 to 800 nm. We targeted R_rs(λ) and IOPs at the following center satellite wavelengths:

SeaWiFS: 412, 443, 490, 510, 555, 670

The latter three wavelength suites all represent possibilities for PACE, with the first of these three (15 wavelengths) indicating aggregated bands to be distributed by the mission (PACE Science Definition Team, 2018) and the second and third representing the full science resolution of the spectrometer (5 nm) spanning two ranges for inversion (410–600 nm vs. 410–670 nm). While it remains tempting to execute the inversion using the full spectrum of the synthesized dataset, we opted to limit our 5 nm examples to a SeaWiFS-like range to specifically compare differences in performance when using contiguous (so called hyperspectral) spectra vs. discrete heritage wavelengths. Furthermore, we feel the extension of SAA inversions into the ultraviolet remains sufficiently underdeveloped to warrant examination in separate works. We explored a more limited hyperspectral range—to only 600 nm in lieu of 670 nm—to evaluate the impact of very low red radiometric signals in the inversion. For comparative purposes with heritage instruments, we also estimated IOP retrievals for the two 5 nm cases (410–670 and 410–600 nm) at 443 nm using a cubic spline fit on the retrieved a_dg(λ), a_ph(λ), and b_bp(λ).

The native synthesized dataset defines its IOPs using four orders of magnitude of C_a (20 values within the range 0.03 to 30 mg m⁻³, N = 500) (IOCCG, 2006). In most cases (specifically noted), we limited our analyses to 0.03 ≤ C_a ≤ 5 mg m⁻³, as the SAA became sufficiently unstable from instrument-to-instrument above 5 mg m⁻³ such that low retrieval sample sizes confounded meaningful reporting of results (N ≅ 350). We discuss the strengths and weaknesses of these synthesized data for our analyses in section Discussion. For the satellite comparisons, we acquired the mission-mean (1997–2010) SeaWiFS Level-3 R_rs(λ) bin file at 9 km spatial resolution from the NASA Ocean Biology Processing Group (OBPG; https://oceancolor.gsfc.nasa.gov; reprocessing version 2018.0). For several analyses, we stratified the results into three C_a ranges—C_a ≤ 0.1 mg m⁻³ (intended to represent oligotrophic water), 0.1 < C_a ≤ 1 mg m⁻³ (mesotrophic), and C_a > 1 mg m⁻³ (eutrophic)—to broadly evaluate performance differences as a function of water type. We note that shape and magnitude of synthesized R_rs(λ) spectra shift from being blue-dominated to red-dominated as the water type transitions from oligotrophic to eutrophic as shown in Figure 1.

Regarding our use of single center wavelengths for this exercise, we acknowledge that different satellite instruments maintain inherently different radiometric bandpasses. Standard NASA processing applies an out-of-band correction using instrument-specific spectral response functions that adjusts the retrieved R_rs(λ) to the nominal band centers after atmospheric correction (Wang et al., 2001; Patt et al., 2003). Given that this is done prior to execution of an SAA, our use of center wavelengths sufficiently mimics the satellite input R_rs(λ)—acknowledging, of course, that the out-of-band adjustment carries uncertainties that inevitably impact the final IOP retrievals.

Our default SAA adopts the default generalized IOP (GIOP) form described in Werdell et al. (2013a). Briefly, R_rs(λ) are first converted to their subsurface values following (Lee et al., 2002):

Subsurface remote-sensing reflectances relate to marine IOPs following (Gordon et al., 1988):

where b_b(λ) is the total backscattering coefficient (m⁻¹) and a(λ) is the total absorption coefficient (m⁻¹). Total backscattering can be expanded as the sum of all backscattering components. Further, each component can be expressed as the product of its mass-specific backscattering spectrum (spectral shape; b^*) and its magnitude or concentration (M):

Similarly, total absorption can be expanded to:

where the subscript w indicates contributions by water. Both b_bw(λ) and a_w(λ) are known (Pope and Fry, 1997; Zhang et al., 2009). Werdell et al. (2013a) describes in the detail the rationale for the default choices of bbp*(λ), adg*(λ), and aph*(λ). By default, we expressed bbp*(λ) as λ^Sbp, where S_bp defines the steepness of the power law and typically varies between −2 and 0 in natural waters (Stramski et al., 2004). We dynamically assigned S_bp to each station/pixel using the R_rs(λ) band-ratio method described in Lee et al. (2002). We expressed adg*(λ) as exp(-S_dg λ), where S_dg describes the rate of exponential decay and typically varies between 0.01 and 0.02 nm⁻¹ (Roesler et al., 1989), and adopted a fixed S_dg of 0.018 nm⁻¹. Finally, we expressed aph*(λ) as the C_a-specific absorption spectrum of Bricaud et al. (1998) seeded using each station/pixel-specific C_a and normalized to 0.055 m² mg⁻¹ at 443 nm. Broadly speaking, many alternatives for such parameterizations exist (as explored in Werdell et al., 2013a), but our choices in parameterization remain of secondary importance for this paper given that we apply them consistently across analyses.

Three unknowns remain in Equations (1–4) after assignment of the eigenvectors, namely M_bp, M_dg, and M_ph. Using R_rs(λ) as input, the SAA estimates these eigenvalues via non-linear least squares (Levenberg-Marquardt) inversion of Equation (2) with minimization of a χ² cost function [see Equation (11) of Werdell et al., 2013a]. We retained only those solutions with viable estimates of M_bp, M_dg, and M_ph [for example, −0.005 ≤ a_dg(λ) ≤ 5 m⁻¹] that resulted in reconstructed R_rs(λ) that differed from the input R_rs(λ) by < 33% in the range 400 to 600 nm. We reconstructed R_rs(λ) using the retrieved eigenvalues as input into Equations (1–4) and defined failure as non-convergence in the inversion. Werdell et al. (2013a) outlines exclusion criteria used to determine viable estimates, as well as the similarities in accuracy shared by this SAA configuration and other common approaches.

We ran the inversion model using the synthesized R_rs(λ) for every satellite input wavelength suite listed in section Data Acquisition. We then statistically compared the model-derived IOPs (hereafter referred to by satellite to identify the input wavelength suite) with the source synthesized IOPs. Our reported statistics include the sample size (N), the Type I linear regression slope, the median absolute error (MAE), and the bias. We calculated MAE as:

and bias as:

Note that observations are log-transformed such that MAE and bias present multiplicative metrics (that is, a MAE of 1.5 indicates that the satellite observations are 1.5 × (50% greater) on average than the synthesized observations and a bias less than unity indicates a negative bias). Note also that Table 1 reports bias and MAE as calculated in Equations (5) and (6), but we discuss these performance metrics using their more intuitive percentage counterparts throughout the Results and Discussion sections.

Next, we compared derived IOPs from each satellite wavelength suite with those from every other satellite wavelength suite. This comparison included determination of unbiased percent differences (UPD). We calculated UPD for two satellites (satellite₁ and satellite₂) as:

and report the median UPD for the complete population of match-ups. A negative UPD indicates that the IOP from satellite₂ is biased low relative to satellite₁.

To evaluate if the results varied with the configuration of the inversion model (that is, if they varied with assignment of spectral shapes), we ran the inversion model for each satellite input R_rs(λ) wavelength suite an additional seven times, sequentially re-assigning one variable per run as: (1) default S_bp from Lee et al. (2002) + 33%; (2) default S_bp from Lee et al. (2002)−33%; (3) default S_dg + 33% (= 0.024 nm⁻¹); (4) default S_dg−33% (= 0.012 nm⁻¹); (5) seed C_a to calculate aph*(λ) using Bricaud et al. (1998) + 50%; (6) seed C_a to calculate aph*(λ) using Bricaud et al. (1998)−50%; and, (7) replacing 0.0949 and 0.0794 in Equation (2a) with 0.0895 and 0.1247 as in Lee et al. (2002).

Finally, we twice ran the default inversion model (as defined in section Inversion Modeling) on the mission-long mean Level-3 SeaWiFS R_rs(λ) to determine if results from our controlled modeling environment were generally representative of real satellite data. To do so, we used OBPG-supported satellite data processing software (l3 gen) to run the inversion model via the GIOP software framework, which allows assignment of wavelengths to be used in the inversion, on the Level-3 SeaWIFS R_rs(λ). The NASA SeaWiFS Data Analysis Software (SeaDAS; https://seadas.gsfc.nasa.gov) package includes this software. The first run used all six SeaWiFS wavelengths. The second run served to mimic VIIRS and, as such, we only assigned the five common SeaWiFS-VIIRS wavelengths. Given the known differences between SeaWiFS and VIIRS R_rs(λ) (Eplee et al., 2015; Barnes and Hu, 2016), using a narrower SeaWiFS band set as a surrogate for VIIRS provided a more suitable controlled environment to explore differences in band selection.

Broadly speaking, the satellite-to-synthesized IOP validation analyses demonstrated equivalent performance across all wavelength suites (Figure 2, Table 1), with performance metrics that mirror those from previous studies that used equivalent simulated datasets (e.g., IOCCG, 2006; Lee et al., 2010; Werdell et al., 2013a). All runs reported very similar MAE, ranging from 5 to 10% for a(443) and 10 to 16% for b_bp(443) and from 5 to 8% for a(555) and 8 to 14% for b_bp(555). In general, we achieved comparable results across blue-green wavelengths for all analyses and, thus, our tables and figures focus only on 443-nm for clarity in presentation of results. MAEs for the absorption subcomponents ranged from 25 to 27% for a_dg(443) and 47 to 54% for a_ph(443). Every wavelength suite carried its lowest MAEs for a(443) and the highest MAE for a_ph(443). The regression slopes behaved similarly. All runs reported very similar values for a(λ), a_dg(λ), and b_bp(λ), ranging from 1.00 to 1.07 for a(443), 1.03 to 1.06 for a_dg(443) and 0.94 to 1.01 for b_bp(443). Across IOPs for each wavelength suite, a_ph(443) deviated most substantially from unity. It also exhibited the broadest range across wavelength suites, ranging from 0.72 to 1.05. Interestingly, the runs with the most wavelengths (PACE, 410–670 nm, and 410–600 nm) reported the a_ph(443) regression slopes that deviated most from unity. The satellite-to-synthesized biases diverged more across runs than either MAE or regression slope. While a(443) biases spanned across unity, all a_dg(443) biases fell below unity (recall that a bias less than unity indicates a negative bias), ranging from 7% for SeaWiFS to 15% for the 410–670 nm run, and all a_ph(443) biases fell either very close to unity or were positive, ranging from exactly unity for the 410–670 nm run to 24% for SeaWiFS and VIIRS. All b_bp(443) biases exceeded unity (were positive), ranging from 1% for the two 5 nm runs to 10% for SeaWiFS.

These results predominantly encompass simulations with C_a ≤ 5 mg m⁻³, as 85% of SAA failures occurred for those simulations with higher C_a values (reducing the viable N to ≅350 as previously mentioned). Few failures resulted from non-convergence in the inversion. Rather, most occurred because of quality assurance metrics applied to the IOP retrievals, including requiring them to fall within predefined acceptable ranges (for example, −0.005 ≤ a_dg(λ) ≤ 5 m⁻¹) and requiring that reconstructed R_rs(λ) differ from input R_rs(λ) by < 33% at every wavelength in the range 400 to 600 nm. While relaxing or modifying these quality assurance metrics would enable extending our analyses into more complex water, we opted to not do so to: (a) ensure that our results encompassed only the highest quality SAA retrievals; (b) limit the introduction of results biases that emerge from global SAA parameterizations that do not properly represent complex waters; and, (c) acknowledge that the simulated dataset itself may be the least representative of natural conditions in this C_a range. Moreover, as simulations with C_a > 5 mg m⁻³ realized so few results (low sample sizes), we opted to exclude them from subsequent analyses and data product intercomparisons. Follow-on analyses will benefit from additional attention to SAAs and datasets that better represent optically complex water.

IOP retrievals for the synthetic dataset varied across wavelength suites and as a function of water type as defined by C_a (Figures 3–6, Tables 2–4). UPDs in a(443), a_dg(443), a_ph(443), and b_bp(443) across all wavelength suites ranged from −3.83 to 2.91%, −3.14 to 5.27%, −11.57 to 6.33%, and −6.09 to 4.74%, respectively, for the lowest C_a waters (Table 2). For comparison, UPDs in a(555), a_dg(555), a_ph(555), and b_bp(555) for SeaWiFS, PACE, and the two 5 nm wavelength suites (the only four suites with 555 nm) ranged from −0.51 to 0.35%, −0.74 to 3.13%, −6.65 to 1.47%, and −2.78 to 1.52%, respectively. Generally speaking, the differences between any two wavelength suites remained constant for all C_a ≤ 0.1 mg m⁻³ (Figures 3–6). Retrievals for the wavelength suites we consider most similar, SeaWiFS, MODIS, and MERIS, generally differed the least—by ≤ 0.92, 0.40, 1.11, and 1.25% for a(443), a_dg(443), a_ph(443), and b_bp(443), respectively, in absolute terms. OLCI retrieved low a_dg(443) and high a_ph(443) relative to all other wavelength suites, perhaps due to the inclusion of a 400 nm band that altered the separation of the two absorption subcomponents. The two 5 nm wavelength suites retrieved high a_dg(443), low a_ph(443), and low b_bp(443) relative to all other wavelength suites. With only several exceptions, they also realized the highest magnitude UPDs.

For the mesotrophic C_a range, UPDs in a(443), a_dg(443), a_ph(443), and b_bp(443) across all wavelength suites ranged from −1.66 to 1.80%, −3.18 to 4.79%, −10.46 to 4.50%, and −2.27 to 2.60%, respectively (Table 3). These ranges for a_dg(443) and a_ph(443) generally match those for the oligotrophic C_a range, whereas these ranges for a(443) and b_bp(443) compress closer to unity. For comparison, UPDs in a(555), a_dg(555), a_ph(555), and b_bp(555) for SeaWiFS, PACE, and the two 5 nm wavelength suites ranged from −1.50 to 0.21%, −1.74 to 3.99%, −7.47 to 1.34%, and −1.32 to 0.24%, respectively. Generally speaking, the differences between any two wavelength suites remained invariant for 0.1 ≤ C_a ≤ 0.5 mg m⁻³, above which they started to diverge (Figures 3–6). SeaWiFS, MODIS, and MERIS differed by ≤ 1.50, 2.32, 0.75, and 2.59% for a(443), a_dg(443), a_ph(443), and b_bp(443), respectively, in absolute terms. Like the oligotrophic subset, OLCI retrieved low a_dg(443) relative to all other wavelength suites. Unlike that subset, however, its a_ph(443) retrievals exceeded all other wavelength suites. As for the oligotrophic subset, the two 5 nm wavelength suites retrieved high a_dg(443), low a_ph(443), and low b_bp(443) relative to all other wavelength suites. Their UPD magnitudes, however, fell more in line with the others.

The eutrophic C_a range (capped at 5 mg m⁻³) realized the most diversity and highest magnitudes in UPDs across all wavelength suites [noting that: (a) SAA retrievals often realize instabilities in such ranges (Werdell et al., 2009, 2018); and (b) the synthesized dataset represents the widest diversity in water in this range (IOCCG, 2006)]. UPDs in a(443), a_dg(443), a_ph(443), and b_bp(443) across all wavelength suites ranged from −16.07 to 12.70%, −19.69 to 17.67%, −53.92 to 12.01%, and −19.50 to 12.99%, respectively (Table 4). For comparison, UPDs in a(555), a_dg(555), a_ph(555), and b_bp(555) for SeaWiFS, PACE, and the two 5 nm wavelength suites ranged from −17.00 to 5.33%, −0.24 to 11.23%, −42.26 to 9.31%, and −10.60 to 3.54%, respectively. The differences between any two wavelength suites varied considerably for C_a > 1 mg m⁻³ (Figures 3–6). Even SeaWiFS, MODIS, and MERIS differed more substantially than in other C_a ranges—by ≤ 15.74, 15.03, 18.64, and 17.12% for a(443), a_dg(443), a_ph(443), and b_bp(443), respectively, in absolute terms. Like the oligotrophic subset, OLCI retrieved low a_dg(443) and high a_ph(443) relative to all other wavelength suites, with the exception of MERIS. As for the other two subsets, the two 5 nm wavelength suites retrieved high a_dg(443), low a_ph(443), and low b_bp(443) relative to all other wavelength suites. Their UPDs for a_ph(443) far exceeded those for the other suites.

Differences in IOPs from SeaWiFS and SeaWiFS-modified-to-look-like-VIIRS largely demonstrated similar patterns. For context, Figure 7 shows SeaWiFS-derived C_a (Hu et al., 2012) and the spatial distribution of olig-, meso-, and eutrophic waters, and Figure 8 shows the SeaWiFS-derived a_dg(443), a_ph(443), and b_bp(443). Naturally, the cumulative distribution of observations falling into each of the trophic levels differs widely in comparison to the synthesized data set, which might be considered in the interpretation of these results. The median SeaWiFS-derived C_a for its olig-, meso-, and eutrophic regions were 0.066, 0.19, and 1.6 mg m⁻³, respectively. The corresponding median UPDs for a_dg(443) for these regions were 1.7, 1.8, and 4.2%, respectively, which approximately track the magnitude of the UPDs at these C_a values in Figure 4 (upper left panel, green line). Note, however, that VIIRS a_dg(443) falls below that for SeaWiFS by ~2% on average across the synthesized mesotrophic subset. Note also that the mesotrophic subset dominates globally and represents a great diversity in water type (Figure 7). The median UPDs for a_ph(443) for the three trophic regions were 9.8, 7.7, and −3.4%, respectively. These also roughly track the UPDs at the corresponding C_a shown in Figure 5, with the exception of the eutrophic subset, which also has a sign difference. Finally, the median UPDs for b_bp(443) for these regions were 8.9, 8.0, and 2.3%, respectively. These UPDs match those shown in Figure 6 for the oligo- and eutrophic subsets. Geographically, a_dg(443) differ the least between SeaWiFS and surrogate VIIRS, with UPDs that hover around unity (Figure 8B). Both a_ph(443) and b_bp(443) show some latitudinal dependence in their UPDs and surrogate VIIRS-derived values almost unequivocally exceed those from SeaWiFS (Figures 8D,F).

As an additional thought experiment to loosely begin to tie our results to marine particle size distributions (PSDs) and, thus, phytoplankton functional types (PFTs), we applied the approach of Hirata et al. (2011) to locate pixels in this SeaWiFS imagery dominated by pico-, nano-, and microplankton (with dominance defined as >50%). On average, SeaWiFS and surrogate VIIRS-derived a_dg(443) differed by 1.7, 1.7, and 2.1% for these size classes, whereas a_ph(443) and b_bp(443) differed by 8.4, 7.6, and 7.8% and 8.0, 7.4, and 8.0%, respectively. These values do not differ substantially from those of the previous exercise as expected, given that Hirata et al. (2011) seeds with C_a, but they offer a modest confirmation of the spatial patterns realized in our analyses. While evaluating the impact of our findings on PSD and PFT retrievals extends beyond the scope of this paper, we feel it remains an obvious follow-on to this work. Many options for estimation of PSDs and PFTs exist and follow-on activities could: (a) build upon our brief example using model-derived distributions of phytoplankton (e.g., Dutkiewicz et al., 2009; Rousseaux and Gregg, 2015); (b) quantify the impacts of varied SAA retrievals within PSD and PFT algorithms that use IOPs as input (see, e.g., IOCCG, 2014; Bracher et al., 2017; Mouw et al., 2017); or; (c) explore in more detail how differences in SAA retrievals for multiple sensors compare to spectral changes in marine IOPs that occur across natural variations in phytoplankton community composition (Alvain et al., 2012).

Our choice in SAA parameterization made little difference in results achieved from the synthesized dataset. Figures 9–12 follow Figures 3–6 with the exception that the reported UPDs represent the average for eight SAA parameterizations. At C_a ≤ 0.5 mg m⁻³, the patterns in satellite-to-satellite UPDs for the ensemble of runs almost exactly mimicked those for the default SAA runs. This C_a range represents roughly 86% of the retrieved ocean in Figure 7. The most obvious difference in UPD patterns emerge at C_a > 2 mg m⁻³, which represents roughly 1.5% of the retrieved ocean. Interestingly, when C_a exceeded 2 mg m⁻³, UPD values for b_bp(443) remained reasonably steady (Figure 12) for wavelength suites that included one or more orange/red bands placed between 550 and 670 nm (e.g., 620 nm on MERIS, OLCI, PACE, and the two 5 nm suites). This may be explained, in part, as due to eutrophic R_rs(λ) spectral features being more distinct in the red and less so in the blue. While differences in SAA retrievals almost certainly increase when directly comparing two parameterizations (vs. comparing one SAA to an ensemble of SAAs), this analysis serves to confirm that our default parameterization sufficiently represents a global average to reinforce that, in general, our results apply universally and do not depend on our choice in SAA parameterization.

LM was employed by company Go2Q Pty Ltd. The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.