This simple sounding question is quite complex to answer objectively because R _rs is a derived parameter and the uncertainties are quite challenging to quantify for each method under a wide array of environmental conditions. Spectra are not only influenced by the optical properties of a large diversity of dissolved and particulate components and the potential contribution from the benthos, but also by measurement artifacts related to the solar and viewing angles, sky and wave conditions, air−water interface, platform disturbance, spatial inhomogeneity in the water column, and distance from land masses (Voss et al., 2017; Bulgarelli and Zibordi, 2020; Shang et al., 2020). If we screen out field data to encompass only a set of predefined or “ideal” conditions with calm seas, clear skies, and constrained solar and viewing angles, for example, then the measurement uncertainty is decreased significantly. This can be justified for the generation of fiducial reference spectra for satellite valiation and calibration, such as data from the Aeronet-OC program (Zibordi et al., 2009). However, as noted by a recent intercomparison exercise, there can still be considerable uncertainty in R _rs even under ideal conditions and with well-calibrated sensors (Tilstone et al., 2020). Moreover, such ideal conditions are rarely met in the field and represent only a fraction of the diverse environmental conditions encountered in natural ecosystems across latitudes and seasons. Retrieval models applied to satellite imagery also contain a much broader set of conditions, in terms of water column, atmosphere and surface states, as well as observation and illumination geometries. As retrieval approaches are diversifying, the community is finding that even non-optimal spectral information can still be useful for algorithm development and parameter validation for a variety of applications. For example, algorithms with narrowband indices or relying on derivatives may be less influenced by surface reflected sunlight than full inversion-type algorithms (Dierssen et al., 2021).

Many standard options are available for processing radiometric data, including flags for low solar zenith angles and low light conditions and for identifying rain, outlined in processing software such as HyperInSpace (Aurin, 2022). Operationally, site-specific criteria and thresholds are often developed based on careful examination of the data from that region. For example, spurious ship-based radiometric data were identified using five different spectral shape metrics tuned for highly absorbing and weakly scattering conditions characteristic of the Baltic Sea (Qin et al., 2017). The thresholds in such metrics are set based on human evaluation and interpretation of the expected data for a region.

For many aquatic data, there are clear standards that can be used to assess the quality of data. Such metrics are still being defined in the aquatic optics community due to the large range in environmental conditions and the many corrections that must be considered for any measurement (e.g., most absorption measurements have an associated scattering correction and most scattering measurements have an associated absorption correction). Commonly, intercomparisons between instruments and methods are the primary means to assess the operational uncertainties (Tilstone et al., 2020) or through simulation and closure studies (Zaneveld, 1994; Tzortziou et al., 2006), although assessing the accuracy of a measurement still remains an outstanding problem (Ruddick et al., 2019; Zibordi et al., 2019). In this paper we do not assess the accuracy of R _rs, but rather develop a method to assess the quality of a water spectrum using data collected and filtered using community standard practices.

Experienced researchers tend to have a “gut” sense for what looks like a reasonable spectrum, but definitive approaches for assessing spectral shape are limited. Questionable spectral shapes are readily identifiable when reflectance values in the near-infrared wavelengths are much higher than anticipated for the type of water collected or when spectra exhibit an exponentially increasing tail from blue to ultraviolet, indicating that reflected diffuse skylight may not be removed fully from the spectrum (Mobley, 1999; Gould et al., 2001). A recent work by Wei et al. (2016) provides a means of scoring a spectrum compared to 23 different simulated optical water types at a maximum of nine wavebands (Wei et al., 2016). This spectral matching approach provides a score as to how close a spectrum matches one of the predefined water types. As noted below, this approach may be limited to standard simulated conditions and may not represent all water types encountered across the vast aquascape (Barnes et al., 2019).

Here, we present a simple quantitative index approach to conduct quality control of an R _rs spectrum and thereby determine whether it appears similar to other ocean color spectra from a wide variety of blue, green and brown water types. The Quality Water Index Polynomial (QWIP) is founded on the Apparent Visible Wavelength (AVW) a one-dimensional geophysical metric of “color,” calculated as the weighted harmonic mean of the reflectance spectrum across a range of wavelengths (Vandermeulen et al., 2020). This metric reduces the hyperspectral or multispectral (after applying sensor-specific correction factors) data to a continuous variable representing the mean “color” in wavelength (expressed in nm). As noted in Figure 10B of Vandermeulen et al. (2020), a good correlation is found between AVW and chlorophyll a concentrations using a global dataset. Building upon this, we developed a relationship between AVW and standard multi-channel waveband indices to identify spectra that fall outside the general trends observed in aquatic optics for optically deep waters. The approach was developed with a large global dataset representing blue, green, and brown waters and was further tested extensively with field and satellite datasets. This simple approach can be rapidly implemented in ocean color processing chains to provide a level of uncertainty about the spectrum and flag questionable or unusual spectra for further analysis.

A simple decision tree was implemented to determine whether the spectral shape of R _rs would be classified as blue-green, green, or brown in color following from the approach presented in (Balasubramanian et al., 2020). This simplified classification (Figure 2) was used to indicate where different types of waters fall within the QWIP schema.

Quality of all spectra was assessed according to the Wei et al. (2016) approach, using the Matlab-based code provided therein. Specifically, a multispectral subset was extracted for each spectrum (at wavelengths (nm): 412, 443, 488, 510, 531, 547, 555, 667, and 678), which was then normalized and compared to each of 23 reference spectra representing disparate optical water types. A spectral angle mapper was used to identify the most spectrally similar reference spectrum. Normalized reflectance data at each wavelength of the subset were then compared to the corresponding reference spectrum and associated boundaries. The final QA score was determined as the number wavelengths for which the reflectance datum fit within the reference boundaries, divided by the total number of wavelengths assessed (nine for this work). Thus QA scores ranged from 0 (at no wavelength does the target spectrum fit within boundaries of the identified water type) to 1 (target spectrum is fully within identified water type boundaries). Wei et al. (2016) qualitatively discusses “very high scores (>0.9)” and “very low QA scores (<0.5),” and reports that 90% of the evaluation spectra had “high QA scores” of >0.8. Based on this, a Wei score of >0.5 was considered as “Passing” and ≤0.5 was considered “Failing.”

The equations used in the QWIP calculation include the Apparent Visible Wavelength (AVW) calculated from 400 to 700 nm at 1 nm intervals and the Normalized Difference Index (NDI) at two wavelengths as formulated below. Any negative values of R _rs are included in the calculations. We note that the QWIP and NDI acronyms are used in the equations below for clarity, such that: A V W = ( ∑ i = 1 n R r s ( λ i ) ∑ i = 1 n R r s ( λ i ) λ i ) (2) N D I =   ( R r s ( λ 2 ) − R r s ( λ 1 ) ) ( R r s ( λ 2 ) + R r s ( λ 1 ) ) (3)

Statistical tests were conducted in Matlab (The Mathworks, Inc.). Model II regression analysis was used presuming similar uncertainty magnitudes in the involved variables.

The algorithm is first developed using a global dataset covering a diverse set of (mostly) optically deep water types. A training dataset was compiled (n = 1,629) that included the CASCK data (see Supplemental Material) and 300 random points selected from each of the PANTHYR datasets, dubbed here the “CASCK-P” dataset (see Figure 1). We evaluated the relationships between AVW (the mean “color” metric, Eq. 2) and other indices to find trends across the wide range of spectral shapes in the dataset. The objective was to find a simple metric using common spectral bands found in a variety of multi-spectral datasets (i.e., wavebands from the VIIRS sensor wavelengths were selected here) where the central tendency of data formed a well-described continuum across the wide range of AVW values and the amount of deviations from the central tendency could be easily scored. To better visualize the data trends, the water color was further differentiated into the categories of blue-green (blue dots), green (green dots), and brown (red dots) following from the decision tree (see Figure 2). Blue-green waters have AVW values ranging from 400 to 510 nm and green waters have AVW from 510 to 590 nm. The AVW for brown waters overlaps with the upper end of the green waters, ranging from 555 to 575 nm.

The QWIP is formally a mathematical model relating an optical index to the AVW. Different examples of the relationships between AVW and other indices are provided in Figure 3 for illustrative purposes. Blue-to-green band ratios, like those used in chlorophyll algorithms and Vandermeulen et al. (2020), showed a reasonable overall trend, but resulted in the values being spread across a log scale with clustering of the green and red spectra (e.g., Figure 3A). Band difference algorithms provided lower predictive power, especially for the brown waters which had much higher differences than the blue and green waters (Figure 3B). Relationships between AVW and maximum wavelength in visible wavelengths (data not shown) were effective for blue and green waters, but proved to have low predictive power for intense green algal blooms with high red edge values and turbid brown waters. Brown waters can have variable peak wavelengths from green to red to near infrared wavelengths (see Figure 2).

The NDI (Eq. 3) provided a means to highlight the variability of logarithmically distributed data on a linear scale such that the distance either above or below the central tendency was scored with a positive or negative value. Different combinations of NDI were systematically evaluated using wavebands found on historic ocean color sensors. For example, the wavebands used in the standard chlorophyll a concentration algorithm (492 and 560 nm) (O’Reilly et al., 1998) did not have high rank correlation with the AVW for green and brown waters (Figure 3C). The best relationship was found using the NDI calculated with blue/green (λ ₁ = 492 nm) and red (λ ₂ = 665 nm) bands (Figure 4A). A 4th degree polynomial fit between NDI (492,665) and AVW described the variability across the AVW range (R ² = 0.974). This QWIP relationship followed the overall objectives in that the central tendency was clearly outlined across a wide range of data and distance was easily scored. To implement, the user would calculate their AVW and NDI (492,665) for a spectrum following from Eqs. 2, 3 and then calculate the QWIP Score as the difference between the measured and the predicted NDI based on the QWIP. Note that the acronyms NDI, AVW, and QWIP are used in the equations below for clarity and the five coefficients in Eq. 4 correspond to five variables p provided below the equation: QWIP = p 1 AVW 4 + p 2 AVW 3 + p 3 AVW 2 + p 4 AVW + p 5 p = ( − 8.399885 x 10 − 9 ,   1.715532 x 10 − 5 , − 1.301670 x 10 − 2 , 4.357838 x 10 0 , − 5.449532 x 10 2 ) (4) QWIP score = NDI ( 492,665 ) − QWIP (5)

As shown in Figure 4A, 97.5% of the CASCK-P data used in the calibration of the QWIP fell within ± 0.1 of the QWIP and 99.4% within ± 0.2 (Figure 4B). Outlier data with QWIP scores greater than ± 0.2 were subject to additional screening to determine any evident spectral anomalies. The points with more negative QWIP scores (Figure 4C) appeared to be related to some types of optically shallow water where the red wavelengths were much lower than anticipated for highly green peaked waters due to water absorption between the scattering seafloor and the sea surface (e.g., Dierssen et al., 2003) suggesting the QWIP score may prove a useful metric for identifying certain types of optically shallow water. Data with high positive QWIP scores often had rising tails in the blue end of the spectrum consistent with data containing residual surface reflected skylight (Gould et al., 2001) (Figure 4C). In such cases, the AVW would be weighted incorrectly towards the blue end of the spectrum resulting in a positive QWIP score.

QWIP scores were calculated for two different above water field datasets for evaluation across a variety of water types with different instruments.

The QWIP procedure was applied to a hyperspectral HICO image (H2012236112610), covering the complex waters of the Nile Delta, and extending north into the clear waters of the Mediterranean Sea (Figure 9A). Pseudo-true color composites of the HICO imagery can be found in Supplementary Figure S1 (see Supplementary_README.pdf). For this demonstration, a binary exclusion flag was set to distinguish spectra that pass/fail a nominal threshold value, as defined below. Using rigid exclusion criteria, 99.55% of the spectra from the atmospherically-corrected HICO image passed quality control with QWIP score < ± 0.2), while 0.29% fell above the QWIP score threshold, and 0.16% fell below the QWIP score threshold (Figure 9B). The average gradient of integral-normalized R _rs(λ) spectra (as a function of AVW) are plotted for these three scenarios, showcasing instances in which data were above/failed (Figure 9C), between/passed (Figure 9D), and below/failed (Figure 9E) the threshold QWIP score. Most flagged spectra were characterized by continuous negative reflectance either below 500 nm, or above 600 nm, and sometimes both. Note, that anomalous negative data should have been screened out in the SeaDAS processing prior to this analysis.

Notably, a series of quality flags already exist within the SeaDAS framework to identify pixels that fail various levels of quality control (https://oceancolor.gsfc.nasa.gov/atbd/ocl2flags/). Pixels identified as land (LAND flag), radiance saturation (HILT flag), and clouds or ice (CLDICE flag) are masked during processing (blue in Figure 10), thus no R _rs data were returned for these pixels. The only other standard flags that were indicated over water for this scene were the ATMWARN and MAXAERITER, which indicate a warning in the atmospheric correction procedure, and that the aerosol model reached the maximum amount of iterations, respectively. The ATMWARN flag included all instances of MAXAERITER, and this flag comprises <1% of the ocean pixels (Figure 10) in the HICO scene. The QWIP identified questionable pixels from regions along the edge of the scan line, inland, and even a few offshore patches that were not flagged by ATMWARN (brown regions). In addition, some of the pixels flagged with ATMWARN had spectral shapes with passing QWIP scores (yellow regions).

A series of additional HICO scenes representing a diverse range of optical water types (H2012237230813, Columbia River outflow, OR, USA; H2014191103614, Danube River Outflow, Romania; H2009344060219, QLD, Australia; see Supplementary Figures S1B–D) were sub-sampled by a range of incremental QWIP score threshold values, to illustrate the connection between QWIP score and spectral shape. Figures 11A−F identify more satellite-derived spectra that failed the nominal quality control criteria, as color coded by discretized QWIP scores. The spectra are often characterized by sharp increases or decreases in the blue range of R _rs(λ) and/or contain negative values. In most cases, as spectral data increasingly deviate from the polynomial relationship between AVW and NDI (492,665), the anomalous spectral features become more prominent.

While the QWIP was developed for hyperspectral measurements, the approach can be applied to multi-spectral data using sensor-specific coefficients to derive the hyperspectral-equivalent AVW as per Vandermeulen et al. (2020). Here, we test the algorithm efficacy using multi-spectral validation measurements (in this case, for SNPP-VIIRS) retrieved from NASA’s SeaBASS. In order to use the QWIP as defined in this manuscript, the AVW derived from multi-spectral measurements must first be translated to a hyperspectral-equivalent AVW through the use of sensor-specific polynomial offsets (Vandermeulen et al., 2020). The most recent updates to the coefficients for all satellite sensors processed by OBPG have been developed and published (Vandermeulen, 2022) (also see Matlab scripts in Supplemental Material).

Applying the same procedure as the HICO analysis, 5.2% and 12.1% of the VIIRS satellite data were flagged as falling above/below a QWIP threshold of 0.2, respectively, and 0.2% and 4.2% of in situ data were flagged as falling above/below the QWIP threshold of 0.2. Note, given that SNPP-VIIRS has the fewest number of spectral channels relative to the other validation data streams, it exhibits a higher uncertainty in the calibrated AVW values relative to many other heritage sensors, with a mean absolute error (MAE) of 1.21 nm, and bias of −0.15 nm. Such uncertainty is expected when converting a product derived from five bands into an approximation of its hyperspectral equivalent. To account for this additional uncertainty, we chose here to relax the nominal threshold value to 0.3 for subsequent reporting, as the 0.2 threshold appeared too stringent for practical application.

Using a nominal threshold of 0.3, we found that 2.1% and 8.2% of the VIIRS satellite data were flagged as falling above/below the QWIP threshold, respectively, and 0.0% and 1.0% of in situ data were flagged as falling above/below the QWIP threshold. The flagged spectra were mostly characterized by elevated blue reflectance when the QWIP scores are above the threshold (Figure 12A), and negative/depressed blue reflectance when the QWIP scores are below the threshold (Figure 12B). A scatter plot comparison between in situ (410 nm) and VIIRS (410 nm) R _rs matchups (Figure 12C) shows that the flagged outliers tend to form tightly grouped clusters with independent trends that deviate significantly from the overall linear regression fit. A modest reduction can be seen in the MAE, and bias between the in situ v. satellite matchups for R _rs (412) (Figure 12C), while differences at all other wavelengths were negligible. Notably, the removal of flagged pixels improves the matching of data frequency distributions at R _rs (412) (Figure 12E) relative to the native dataset (Figure 12D).