How many independent quantities can be extracted from ocean color?

Products derived from remote sensing reflectances ( Rrsλ ), for example, chlorophyll, phytoplankton carbon, euphotic depth, or particle size, are widely used in oceanography. Problematically, Rrsλ may have fewer degrees of freedom (DoF) than measured wavebands or derived products. Here, we show that a global sea surface hyperspectral Rrsλ dataset has DoF = 4. MODIS‐like multispectral equivalent in situ data also have DoF = 4, while their SeaWiFS equivalent has DoF = 3. Both multispectral‐equivalent datasets predict individual hyperspectral wavelengths' Rrsλ within nominal uncertainties. Remotely sensed climatological multispectral Rrsλ have DoF = 2, as information is lost by atmospheric correction, shifting to larger spatiotemporal scales, and/or more open‐ocean measurements, but suites of Rrsλ ‐derived products have DoF = 1. These results suggest that remote sensing products based on existing satellites' Rrsλ are not independent and should not be treated as such, that existing Rrsλ measurements hold unutilized information, and that future multi‐ or especially hyper‐spectral algorithms must rigorously consider correlations between Rrsλ wavebands.

independent and should not be treated as such, that existing R rs λ ð Þ measurements hold unutilized information, and that future multi-or especially hyper-spectral algorithms must rigorously consider correlations between R rs λ ð Þ wavebands.
Ocean color satellites have revolutionized the study of ocean ecology and biogeochemistry in recent decades by providing a near-continuous global picture of surface ocean properties (Hovis et al. 1980;O'Reilly et al. 1998). Satellites measure the spectral radiance emanating from the ocean and atmosphere. Remote sensing reflectance (R rs λ ð Þ) is obtained following the removal of the contribution of atmospheric and surface effects and normalization to downwelling solar irradiance. Algorithms have been developed to estimate numerous biogeochemcally relevant surface variables from R rs λ ð Þ, such as chlorophyll concentration (Chl, [μg L À1 ]) (O'Reilly et al. 1998;Hu et al. 2012), the spectral slope of the particle size distribution (ξ) (Kostadinov et al. 2009), the concentrations of phytoplankton and particulate organic and inorganic carbon (C phyto , POC, and PIC, [μg L À1 ]) (Graff et al. 2015;Evers-King et al. 2017;Mitchell et al. 2017), euphotic layer depth (Z eu [m]) (Lee et al. 2007), and, using additional input variables, net primary production (NPP, [mg m À2 d À1 ]) (Behrenfeld and Falkowski 1997;Westberry et al. 2008;Silsbe et al. 2016). Such products are used in a wide variety of applications, such as validation of complex ocean ecosystem and biogeochemistry models (Dutkiewicz et al. 2020;Cael et al. 2021) or as inputs for simpler models that predict other variables such as vertical particulate organic carbon fluxes from ocean color (Siegel et al. 2014;Cael et al. 2017;DeVries and Weber 2017;Bisson et al. 2020;Nowicki et al. 2022).
Existing R rs λ ð Þ data are multispectral, measured at several wavebands. Derived products generally rely only on a subset of these wavebands and are commonly expressed as functions of band ratios between just two wavelengths (Hu et al. 2012). Some algorithms attempt to simultaneously estimate multiple products to match the full R rs λ ð Þ spectrum, for example, the generalized inherent optical properties approach (Werdell et al. 2013). However, the most widely used products, such as for Chl and POC, treat all outputs as independent quantities and are fully empirical.
Correlations between R rs λ ð Þ at different wavebands are strong (Huot and Antoine 2016) presenting multiple potential issues for both users and developers of derived products. If multiple products are used simultaneously and treated as independent when they are in fact not, this can lead to overconfidence in model skill or miscalculation of uncertainties. Adding different (yet correlated) satellite products to a model can result in model output redundancy (Bisson et al. 2020). These issues will be exacerbated by the hyperspectral resolution of the next generation of ocean color satellites, namely the Plankton, Clouds, Aerosols, and Ecosystems (PACE) satellite scheduled to launch January 2024 . In addition to the common suite of multispectral products, PACE also plans to enable characterizations of phytoplankton communities, for example, Chase et al. (2017), substantially increasing the number of products available from R rs λ ð Þ. The strong correlations among R rs λ ð Þ wavelengths can be framed in terms of the degrees of freedom (DoF) of R rs λ ð Þ measurements and suites of derived products. DoF represents the effective number of dimensions of a dataset after accounting for correlations and uncertainties between variables and is in essence the number of independent variables in that dataset. It has been shown that the DoF of globally distributed nearsurface measured hyperspectral absorption spectra is $ 5 (Cael et al. 2020). This could be considered a possible upper limit for the DoF of satellite-measured R rs λ ð Þ given higher uncertainties on satellite measurements-particularly associated with atmospheric correction (Cael et al. 2020;Bisson et al. 2021). The DoF of PACE's hyperspectral measurements might then be expected to be much lower than the number of wavelengths for which it will measure R rs λ ð Þ, which will appreciably affect how hyperspectral satellite R rs λ ð Þ products should be constructed. For both existing and future satellite R rs λ ð Þ, understanding the DoF of R rs λ ð Þ measurements and derived products is crucial for appropriate usage and optimal construction of such products.
Here, we investigate the DoF of R rs λ ð Þ. We show that a global sea surface hyperspectral R rs λ ð Þ database has four DoF. Coarsening hyperspectral R rs λ ð Þ to their moderate resolution imaging spectrometer (MODIS) equivalent retains four DoF, though the sea-viewing wide field of view sensor (SeaWiFS) equivalent only has three DoF. Both multispectral equivalents predict individual hyperspectral R rs λ ð Þ wavelengths within nominal uncertainties for satellite sensors. For climatological R rs λ ð Þ and derived products, both MODIS-Aqua and SeaWiFS R rs λ ð Þ have two DoF, suggesting R rs λ ð Þ complexity is lost either through atmospheric correction, relatively more inclusion of open-ocean data, or averaging over larger scales in space and time. Suites of derived products only retain one DoF. Therefore derived products should not be treated as independent by users. These findings have substantial implications for the construction and use of multispectrally and hyperspectrally derived ocean color products.

Sea surface R rs : Hyperspectral vs. multispectral
We first analyze a global sea surface hyperspectral R rs λ ð Þ dataset to determine its DoF and how the DoF depends on spectral resolution (Chase et al. 2017;Kramer et al. 2022). The dataset includes R rs λ ð Þ data at 191 locations at an effective 3.35 nm resolution (Chase et al. 2017) from 400 to 800 nm, linearly interpolated to 1 nm (Fig. 1). We trimmed spectra to Cael et al. Quantites extracted from ocean color 700 nm due to the large fraction of missing values and zeros > 700 nm; our conclusions are not affected by using a lower or higher maximum wavelength. The dataset includes measurements taken from 2004 to 2018 evenly distributed across months of the year, and from all major ocean basins ranging in latitude from 41 S to 74 N. We also compare these data to their MODIS-Aqua and SeaWiFS multispectral equivalents by convoluting the hyperspectral R rs λ ð Þ with the MODIS-Aqua and SeaWiFS spectral response functions (available at https:// oceancolor.gsfc.nasa.gov/docs/rsr/HMODISA_RSRs.txt and https://oceancolor.gsfc.nasa.gov/docs/rsr/SeaWiFS_RSRs.txt) to generate 10-waveband and 6-waveband datasets which correspond to what each instrument would have measured from the same optical input that the radiometer received when generating the hyperspectral R rs λ ð Þ data. We then apply principal component analysis (PCA) (Wold et al. 1987) to these 301-, 10-, and 6-dimensional R rs λ ð Þ datasets. PCA is a widely used method to reduce the dimensionality of datasets by identifying orthogonal vectors that explain the most variance in the data. PCA is linear in nature, which may result in an overestimation of effective dimensions by poorly approximating nonlinear relationships between variables (e.g., a PCA performed on the pair (x, y) where y = x 2 will yield two DoF). Nonlinear generalizations do exist (Weinberger et al. 2004;Scholz et al. 2008), though these are less widely applied due to their additional complexity and computational requirements that make interpretation challenging. One may therefore consider the DoF we report to be upper bounds. We perform a PCA on each R rs λ ð Þ dataset, after standardizing each waveband. We use the broken-stick rule to choose the DoF, which states that the DoF is equal to the  number of components that explain more variance than would be expected by randomly distributed data; this method was shown to be more consistent than a suite of others in a comparison (Jackson 1993; note that this threshold b i ð Þ ¼ 1 n P k¼i…n 1=k is a function of the dimensionality of the dataset, so the significance cutoffs are different in Figs. 2, 4). These results can be shown visually as a "scree" plot, which plots the percentage of variance explained by each component and for randomly distributed data; the DoF is the number of components with a higher percentage of variance explained than would be expected for randomly distributed data. Our figures also visibly demonstrate that one would get the same results from using the scree plot rule, which states that the DoF is equal to the number of components not sitting on the straight line made by the higher-order components, and was found to consistently capture the correct DoF plus one when the first point on this straight line was included (Jackson 1993).
PCA reveals the hyperspectral R rs λ ð Þ dataset has four DoF (Fig. 2); the first four components explain 54%, 33%, 8%, and 2%, totalling 97%, of the variance. The first four MODIS-Aqua equivalent R rs λ ð Þ principal components have very similar percentages of variance explained: 49%, 37%, 10%, and 2%, totalling 99% of the total variance. In contrast, the first three SeaWiFS equivalent R rs λ ð Þ principal components explain 63%, 28%, and 8%, totalling 99%, of the variance. This suggests that the hyperspectral R rs λ ð Þ have four DoF, or four independent variables within the data, and that these four variables are effectively captured when reducing spectral resolution to the 10 MODIS-Aqua wavebands, but not to the 6 SeaWiFS wavebands.
Note that this difference of 1 DoF between the SeaWiFSand MODIS-equivalent data is because MODIS includes a band centered at 645 nm which captures spectral variation in the range of $ 610-650 nm, a spectral region which is not covered by any SeaWiFS waveband. When this 645 nmcentered waveband is excluded from the MODIS analysis, we find three DoF for the remaining nine MODIS wavebands (this is not true of other wavebands, e.g., the fluorescence waveband at 673-683 nm). Irrespective of the exact number of DoF, Fig. 2 demonstrates that the bulk of the variance in these data, whether hyperspectral or multispectral, collapses along a few modes of variation, with the first two modes containing almost all ($ 90%) of the explanatory power. This suggests that only two independent quantities can be estimated accurately with these data. Even if one or two additional quantities can be estimated independently (i.e., DoF = 3 or 4), regardless of the exact number of quantities, these will necessarily be estimated with low signal-to-noise ratios.
The ability of coarsened, MODIS-equivalent data to obtain the same number of DoF as the hyperspectral dataset is further supported by predictions of hyperspectral R rs λ ð Þ from multispectral equivalents. To illustrate this, for each hyperspectral wavelength we perform a multivariate linear regression of R rs λ ð Þ at that wavelength regressed against R rs λ ð Þ at each waveband of both the MODIS-Aqua and SeaWiFS equivalent R rs λ ð Þ. For all wavelengths < 578nm in the SeaWiFS case and 582 nm in the MODIS-Aqua case, the root-meansquare-error (RMSE) is smaller than 5% of the mean R rs λ ð Þ at that wavelength, where 5% is a nominal relative uncertainty for satellite R rs λ ð Þ (Fig. 3). Even for wavelengths greater than this, the RMSE is still very small in absolute terms, < 0.00007sr À1 , far smaller than the nominal 0.0003 sr À1 absolute error for 1 km-by-1 km pixels for PACE (Gordon and Wang 1994). This underscores the extent to which different wavelengths' R rs λ ð Þ are correlated and explains the ability of MODIS-Aqua equivalent multispectral R rs λ ð Þ to preserve the dimensionality of hyperspectral R rs λ ð Þ. The fact that SeaWIFSlike R rs λ ð Þ can accurately predict hyperspectral R rs λ ð Þ to within PACE uncertainties but has fewer DoF than the in situ hyperspectral dataset is a reflection of the lower uncertainty on the in situ dataset than the expected PACE R rs λ ð Þ, and suggests that PACE R rs λ ð Þ may have fewer DoF than the in situ hyperspectral dataset.

Climatologies: R rs vs. products
The analysis above is based on instantaneous, sea surface R rs λ ð Þ values. The power of satellite R rs λ ð Þ and derived products, however, lies in their near-continuous global spatial coverage, and many users are primarily interested in climatological data, that is, the the coarsest spatial and temporal scales. In this section we therefore analyze climatological R rs λ ð Þ and derived products, again via PCA to determine DoF.
We generated a 1 Â 1 monthly climatology for SeaWiFS R rs λ ð Þ, excluding 2009-2010 due to known instrument issues (Siegel et al. 2014), using data downloaded from https:// oceancolor.gsfc.nasa.gov/. Note that as we are using climatological data, this is a global scale analysis. We did the same for MODIS-Aqua, spanning July 2002-June 2022. We generated analogous climatologies, from each satellite over the same periods and spatiotemporal resolutions, for Chl, C phyto , POC, PIC, Z eu , ξ, the fraction of biovolume in the microplankton size class f micro calculated from ξ as described in (Kostadinov et al. 2009), the particulate backscatter to chlorophyll ratio b bp : Chl, and NPP as estimated by the CAFE (Silsbe et al. 2016) and CbPMv2 (Westberry et al. 2008) models. Chl, POC, and PIC were downloaded from https://oceancolor.gsfc. nasa.gov/, as was b bp to calculate C phyto according to Graff et al. (2015) and b bp : Chl and the diffuse attenuation coefficent at 490 nm to calculate Z eu according to Lee et al. (2007); SeaWiFS ξ and f micro were derived as in Kostadinov et al. (2009); and NPP products were downloaded from http:// sites.science.oregonstate.edu/ocean.productivity/index.php. In total we then have climatologies for MODIS-Aqua, SeaWiFS R rs λ ð Þ, and 10 derived products. We consider the six products Chl, C phyto , POC, PIC, ξ, and Z eu , to be core products and f micro , b bp : Chl, CAFE NPP, and CbPMv2 NPP to be ancillary products as these are either derived from the core products or rely on ancillary data other than R rs λ ð Þ.
We note that a PCA on the MODIS-Aqua climatologies of R rs λ ð Þ and products other than ξ and f micro yields the same results as those for SeaWiFS below, so we focus here only on the SeaWiFS climatologies because ξ and f micro are not readily available for MODIS-Aqua. We find two DoF for climatological SeaWiFS R rs λ ð Þ, but only one for the products (Fig. 4). This result is not sensitive to which combination of products is used; for example, including all the ancillary products as well still results in one DoF for the products. This result is also not sensitive to log-transformations of the variables that are log-normally (e.g., Chl, POC, PIC, C phyto ; Campbell 1995) or log-skew-normally (e.g., NPP; Cael 2021; Cael et al. 2018) distributed, or removal of outliers, zeros, or negative values.
That R rs λ ð Þ have more DoF for the data in the previous section than for satellite-derived climatologies suggests that some reduction of complexity of the data occurs via some combination of increased sensor noise relative to ship-based data, atmospheric correction, or averaging over large space and time scales (Scott and Werdell 2019). Two DoF remain in satellite climatological R rs λ ð Þ for both SeaWiFS and MODIS-Aqua, indicating the possibility of generating two independent products from these data. The suite of products tested above, however, has one DoF. This is likely due to derived products' appreciable uncertainties and/or strong correlations with chlorophyll. POC, ξ, and Z eu , for instance, have Spearman rank correlations (across all months and 1 grid cells) of > 0.9 with Chl. C phyto 's rank correlation with Chl is still fairly high, at 0.61, and is lower largely due to small fluctuations when both are small; a simple spline fit of log(C phyto ) against log(Chl) yields an r 2 of 0.7.
The exception is PIC, which has a rank correlation with Chl of 0.11. For typical R rs λ ð Þ values, however, PIC is highly uncertain-that is, PIC estimates are very sensitive to small variations in R rs λ ð Þ-as substantiated by the following analysis. We performed a simple sensitivity analysis with the standard two-band PIC algorithm used by NASA for all but the most optically bright waters (see https://oceancolor.gsfc.nasa. gov/atbd/pic/). We calculated PIC for the climatological median R rs λ ð Þ at 443 and 555 nm and for 5% variations, converting to normalized water-leaving radiance by multiplying by the global mean extraterrestrial solar irradiance. We then perturbed these R rs λ ð Þ values with Gaussian noise at the 5% level, corresponding to the nominal uncertainty in R rs λ ð Þ. This noise at 443 nm results in 68% noise in PIC. By contrast, POC only varies 5% with these 5% variations in R rs λ ð Þ at either wavelength. This indicates that in the bulk of cases, satellitederived PIC is highly uncertain, on the order of 70% (and note the PIC uncertainty will be magnified more when considering documented uncertainties for R rs λ ð Þ of 15-40% in some regions; Bisson et al. 2021). In contrast, for relatively bright waters, the same exercise resulted in PIC variations of < 10%, indicating that this algorithm performs well in instances when PIC values are high. Nonetheless, the high sensitivity to typical uncertainty in R rs λ ð Þ for median waters explains why we find one DoF for the products even though PIC and Chl are not strongly correlated: derived PIC is noisy most of the time.
These results have two key implications. One is that there is additional information in climatological R rs λ ð Þ that is not included in current derived products, because climatological R rs λ ð Þ has more DoF than a suite of climatological products. The other implication is that these products are not at all independent, because a suite of them only has one DoF. For instance, a numerical ecosystem model that reproduces the satellite-derived climatology of chlorophyll and of the particle size distribution's spectral slope should not be considered to be capturing two independent properties of the Earth system. When using satellite products as inputs to other models, these products and their propagated uncertainties must be treated simultaneously rather than independently.
The results presented here are appropriate for global and hence primarily open ocean analyses, composed primarily of Case 1 waters where optical variability is dominated by chlorophyll (Morel and Prieur 1977). It is therefore arguably unsurprising that the suite of R rs λ ð Þ-derived products produced one DoF. Coastal and inland waters' optical variability is influenced by other constituents, such as colored dissolved organic material (CDOM), inorganic particles, and other pigments (Brown et al. 2008;Nelson and Siegel 2013). Analyses focused on these waters is likely to reveal a higher number of DoF from both R rs λ ð Þ and derived products. However, we note that the in situ dataset used here ( Fig. 1) represents waters with R rs λ ð Þ variability similar to that of the ocean as a whole, which can be seen by comparing the variation in R rs λ ð Þ at each MODIS-Aqua wavelength from global satellite data with the same satellite data subsampled to the locations with in situ measurements (or the closest non-cloudy location). Subsampled satellite measurements have similar, and slightly lower, R rs λ ð Þ in bluer wavelengths, indicating that the in situ dataset is oriented more toward optically complex coastal waters with substantial CDOM. This suggests that part of the explanation for the drop in DoF in satellite-derived climatologies comes from the fact that the in situ dataset sampled, as a whole, more optically complex waters. For future work, it would be valuable to perform a similar analysis to what we have done here for coastal waters, regional or shelf seas, or other more optically complex environments.
We find that both R rs λ ð Þ and variables derived from R rs λ ð Þ are highly inter-correlated, reducing the number of DoF associated with each, with a greater reduction in DoF in the derived products. This becomes a problem when products are derived using empirical relationships with R rs λ ð Þ, and especially when the same wavelengths are used for the products that are assumed to be independent of each other; for example, over much of the ocean PIC, POC, and chlorophyll all are functions only of R rs λ ð Þ at two wavelengths, at (or near, depending on the sensor) 443 and 555 nm. Certain combinations of PIC, POC, and chlorophyll, which may occur in the surface ocean, are therefore impossible to find using these algorithms. This is distinct from algorithms, typically called "quasi-analytical" or "semi-empirical", that use known or assumed spectral shapes for absorption and scattering properties of optical constituents that can be related to the same derived products, such as PIC, POC, and chlorophyll (Werdell et al. 2013). These approaches may result in similar correlations and DoF between derived products, but do not inherently have the same problems as empirical approaches. We note that PACE will have, in addition to hyperspectral visible bands, UV bands from 350 nm as well as spectral polarized bands. These measurements are expected to both improve the atmospheric correction (hence reduce the R rs λ ð Þ uncertainties) as well as provide their own ocean signals, both of which may increase the DoF compared to those found here. UV data in particular is potentially rich with information about phytoplankton physiology and community structure, and independent of variability at other wavelengths, though there will be challenges associated with simultaneously using it for atmospheric correction and extracting biogeochemical information, and not enough data exist at present for reliable statistical analysis. In addition, it has been shown that adding other environmental variables such as SST can add useful information to inversions of phytoplantkon groups, for example, Chase et al. (2022) and thus another approach to increase DoF for inversions by adding relevant and independent information (e.g., mixed-layer depth and nutrients from BGC-Argo assimilating models). It may also be fruitful to include spatiotemporal information to specify, for example, where blooms of a particular plankton type are expected.

Conclusion
The results presented here highlight the high degree of codependence between remote sensing reflectances at different wavelengths and of the products derived from these reflectances. For users of products based on existing reflectances, this primarily means factoring in the relationships between products when using more than one simultaneously. For the algorithms that generate these products from existing reflectances, these results indicate a potential to improve the suite of available products to be more accurate and precise, and to account for the relationships between products and R rs λ ð Þ wavebands. One way to do this, consistent with the findings above, would be to derive a single product such as chlorophyll as a function of all reflectance wavebands, derive an anomaly from chlorophyll-based expectations of a secondary product, then specify all other products explicitly as a function of these two, along the lines of Alvain et al. (2005).
These findings are most relevant for algorithms that will generate products from hyperspectral reflectances in the future. The small number of DoF in hyperspectral reflectances indicates that only a few quantities can be estimated independently, and that different wavelengths' reflectances as measured from space will be strongly correlated. Complex algorithms that utilize the full spectrum of reflectance will need to factor in these correlations in order to generate reliable products. Crucially, if more than a few products are generated from hyperspectral reflectances, as is likely the case, such algorithms will also need to output the covariance information encoding the uncertainty in each product and the relationships between them. The fact that hyperspectral reflectances can be predicted within nominal uncertainties by their multispectral equivalents suggests that hyperspectral resolution can play a role in improving ocean color products, but that it will be challenging to provide a substantially finer-grained picture of surface ocean ecosystems and biogeochemical cycles. In particular, these results present a fundamental challenge to (or at least ceiling on the ecological resolution of ) algorithms that attempt to extract the abundance of different phytoplankton functional types from remote sensing reflectance. Here by relying on PCA we have focused on broad, first-order variations, but where such resolution may be most useful and generate novel insights is in investigating outliers and rare events, such as blooms or binning data over coherent features like eddies, where e.g., monospecific signatures may be resolved with spectral precision.