The primacy of dissolved organic matter to aquatic light variability

Absorption and scattering by optically active constituents (OACs) modify the sunlit aquatic light environment, facilitating the derivation of biogeochemical data products at scales spanning in situ to satellite observations. Excluding solar illumination, plus geometric and atmospheric effects, variability in an optical parameter arises from changing OAC concentrations, wherein observed patterns in the spectral evolution of data products are associated with the connectivity and spatiotemporal dynamics of OACs. In open-ocean waters far from terrestrial and riverine inputs, the content and mixture of OACs principally relates to the dynamics of phytoplankton and the microbial loop – a trophic pathway describing the cycling of microbial primary producers, remineralizers (e.g., bacteria and archaea), plus dissolved organic and inorganic materials (as applicable). Historical bio-optical models for the open ocean primarily invoke chlorophyll a concentration (C_a) – a commonly used proxy for phytoplankton biomass – as the ubiquitous independent variable governing optical data products such as the normalized water-leaving radiance, ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$. Formulation of ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ as a function of C_a invokes an idealized food chain, wherein phytoplankton are the dominant control of OACs, including the colored (or chromophoric, depending on the literature) portion of the dissolved organic matter (DOM) pool, hereafter CDOM. This prescription, in which C_a maximally explains oceanic light variability (hereafter primacy), is tested herein using eigenanalysis – e.g., an empirical orthogonal function analysis, principal component analysis, or other eigendecomposition depending on the literature. Analyses using three independent bio-optical datasets assess the shapes and associations of the principal and secondary eigenfunctions of aquatic ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ observations. The analyses reveal ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ variations to be more strongly associated with changes in CDOM rather than C_a – even for purely oceanic datasets – indicating that CDOM dynamics are more variable and exhibit greater independence from C_a than formerly ascribed. Blue and green band-ratio algorithms routinely used for remote sensing of C_a are found to be maximally sensitive to variability in CDOM rather than C_a based on validation tests of ocean chlorophyll (OC) algorithm performance (e.g., R² of 0.85 versus 0.78), plus partial correlation coefficients relating eigenfunction scalar amplitude functions to field or derived observations. Eigenanalyses applied to spectral subsets of the data indicate expansive spectral range observing improves the independence in retrieving CDOM absorption and C_a. The spectral subset comparisons indicate expanded spectral observations supported by recent domestic and international satellite missions constitute a new and unique opportunity to optically characterize surface ocean phytoplankton stocks without relying on explicit or implied empiricisms requiring CDOM and other OACs to covary with C_a. The shapes and associations of the eigenfunctions suggest a greater diversity of trophic pathways determine OAC dynamics – e.g., in addition to phytoplankton contributing CDOM via cellular lysis, excretion, and grazing – and are consistent with advancing knowledge of the microbial loop in the decades after bio-optical formulations based on C_a were proposed.

1 Introduction

More than seven decades ago, transitions in the visual perception of electromagnetic emission (hereafter “color”) from the deepest blue waters to the much shallower green coastal zone were attributed primarily to changes in the colored dissolved organic matter (CDOM) content (Jerlov, 1951; Yentsch, 1960; Kalle, 1966). Known absorption characteristics of photosynthetic pigments, principally absorption associated with the content of chlorophyll a, C_a, supported analogous dependencies of color on phytoplankton content, which motivated investigations into optical approaches for retrieving surface ocean phytoplankton abundance (Ryther and Yentsch, 1957; Tyler, 1961). Significant overlap in C_a and CDOM spectral absorption characteristics, combined with limited spectral resolution of early optical technologies, however, prevented the partition of C_a and CDOM signals (Yentsch, 1960). Contents of C_a and CDOM covary in oceanic environments away from the confounding effects of terrestrial and riverine inputs, and a generalized relationship between color and C_a content was demonstrated using airborne surveys (Clarke et al., 1970) – despite acknowledged confounding factors, including (but not exclusively) solar glint, sky state, and natural variability in the content of other optically active constituents (OACs) such as CDOM or suspended inorganic particles. Natural variability in OACs was temporarily managed by stipulating that OACs were adherent to stable empirical relationships in open-ocean environments, which were operationally defined as “case-1” waters (Morel and Prieur, 1977). The case-1 prescription circumvents uncertainty associated with natural variability in the relative content of different OACs so that the color of seawater could be related to individual biogeochemical parameters, namely C_a (Morel, 1980; Grew and Mayo, 1983; Gordon et al., 1983). The approach was consistent with contemporary knowledge of the grazing food chain – wherein phytoplankton growth and grazing regulate oceanic ecosystems (Kalle, 1966; Steele, 1974) – prior to growing recognition of the importance of microbial loop dynamics (Azam et al., 1983; Azam, 1998).

The launch of the Coastal Zone Color Scanner (CZCS) satellite mission in 1978 enabled remote assessment of electromagnetic flux emitted from aquatic surface waters in the visible (VIS) spectral domain – hereafter “ocean color” – at large spatial scales (Hovis, 1981). CZCS imagery corroborated linkages between ocean color and C_a, e.g., by demonstrating elevated oceanic C_a within eddy-like turbidity patterns (Gordon et al., 1980). Satellite remote sensing of oceanic apparent optical products (AOPs) relies principally on the normalized water-leaving radiance, ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$, or its conjugate product the remote-sensing reflectance, R_rs(λ), where λ denotes wavelength. Both were – and continue to be – derived for ocean color satellite observations by ascribing null oceanic flux at the longest signal-limited wavelength to facilitate the partitioning of oceanic and atmospheric observed signals (Gordon et al., 1983). As technology inexorably advanced, the wavelength associated with null flux was incrementally adjusted from the VIS (i.e., red) to the near-infrared (NIR), and then to the shortwave infrared domain (Gordon, 2021) – necessary, in pertinent part, because null flux approximations are inconsistent with the optical properties of the sunlit aquatic environment (Houskeeper and Hooker, 2023). Inversions relating ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ or R_rs(λ) waveband ratios to C_a – hereafter ocean chlorophyll (OC) algorithms – were defined for the CZCS mission (Clark, 1981) and subsequently extended to, for example, the Sea-viewing Wide Field-of-view Sensor (SeaWiFS), the MODerate resolution Imaging Spectroradiometer (MODIS), the Second-Generation Global Imager (SGLI), and other dedicated ocean color missions (Aiken et al., 1996; O'Reilly et al., 1998; Hirata et al., 2014; O'Reilly and Werdell, 2019; Isada et al., 2022).

As satellite missions progressed, pigment quantitation from water samples – requisite for developing and validating empirical OC algorithms – evolved from spectrophotometric and fluorometric techniques to high-performance liquid chromatography (HPLC) methods (Tyler, 1961; Smith et al., 1981; Hooker et al., 2000) and ultimately advanced to ultra HPLC (Suzuki et al., 2015). The trajectory in advancing laboratory methods allowed a concomitant improvement in OC algorithms by including chlorophyll allomers, epimers, and degradation products (Claustre et al., 2004), while lowering pigment uncertainties (Van Heukelem and Hooker, 2011; Hooker et al., 2012b). This era of development also included contemporaneous (i.e., occurring during the same time) advances in commercial-off-the-shelf (COTS) in situ optical instrumentation. The fixed-wavelength microradiometer (Morrow et al., 2010) reduced the size and weight of instruments while providing purely digital operations with improved dynamic range and accuracy. When coupled with a compact backplane equipped with hydrobaric buoyancy and small digital thrusters (Hooker et al., 2018a), oceanic, coastal, and inland (hereafter “global”) waters could be sampled with the same instrument suite. The new capabilities included autonomous at-sea operations (Hooker et al., 2018b) – plus onshore support measurements and semiautonomous airborne surveys (Guild et al., 2020) – to provide data products for the calibration and validation of the next generation of ocean color satellite missions (Hooker et al., 2018c).

Over the time span of advancing satellite, laboratory, and field capabilities, the architecture of OC algorithms remained largely unchanged, i.e., OC algorithms continued to rely on VIS band ratios (O'Reilly and Werdell, 2019), although improvements were obtained, e.g., in oligotrophic waters (Hu et al., 2012) and by simultaneous (VIS) retrieval of multiple OACs (e.g., Garver and Siegel, 1997). Continuity of satellite OC algorithm retrievals supports global monitoring of oceanic C_a needed to assess planetary change (McClain et al., 2006), but correspondence between the OC algorithm data product and oceanic C_a content varies significantly in space and time (Dierssen, 2010; Sauer et al., 2012). For example, natural variability in OAC relationships (Dierssen and Smith, 2000; Siegel et al., 2005) injects regional biases into OC algorithmic products, even for large oligotrophic (purportedly case-1) water bodies, such as the Mediterranean Sea (Claustre et al., 2002; Morel and Gentili, 2009). Evolving climate and biological conditions further alter OAC relationships and derivation of optical data products, thereby modifying OC algorithm performance (Sauer et al., 2012). Variability in OAC relationships confounds detection of long-term C_a trends (Dierssen, 2010), although trends in the fundamental remote sensing optical data products, e.g., ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ or R_rs(λ), are detectable (Dutkiewicz et al., 2019). Relating changes in oceanic optical properties to specific OACs – e.g., robust estimation of C_a – would help to assess changes in oceanic ecosystem functionality and provide context for cause-and-effect relationships in aquatic environments coinciding with a changing global climate.

Interpreting changes in oceanic AOPs remains challenging – as noted in earlier historical studies, e.g., Yentsch (1960) – due to confounding similarities in the spectral absorption of OACs within the VIS domain relative to the spectral range, resolution, and signal-to-noise ratio (SNR) of legacy ocean color datasets combined with non-adherence to an absolute radiometric scale (Houskeeper and Hooker, 2023). An eigenanalysis for modeled aquatic spectra demonstrated optical data products require the shortest available wavelength (modeled to 400 nm at the time), i.e., the spectral end member, to support discrimination of independently varying CDOM and C_a (Sathyendranath et al., 1987). As field technologies advanced, spectrally expansive above- and in-water observations were shown to provide optimal estimation of CDOM across global conservative waters with maximal independence from C_a (Hooker et al., 2020; Houskeeper et al., 2020) – with spectral expansivity referring to the spectral range of compliant data products, with the most spectrally expansive observations characterizing the greatest spectral range. Conservative waters – wherein the inflow and outflow of properties constrain the range in the gradient of a constituent (Hooker et al., 2020) – include all case-1 or otherwise oligotrophic waters, while also including significant representation from optically complex, hereafter case-2 (e.g., including coastal and inland), waters. The inclusion of case-2 waters means that regional (or global average) empirical OAC relationships are not invoked, and global algorithms are, thus, more capable of partitioning OACs and maintaining consistent performance under changing environmental conditions (Houskeeper et al., 2020). New algorithmic approaches included using ratios of the most spectrally separated optical data products (Hooker et al., 2013), especially those from the UV and NIR spectral domains, an approach hereafter termed end-member analysis (EMA) following Hooker et al. (2020). EMA analyses were subsequently extended to single-channel UV algorithms (Hooker et al., 2021b), indicating the importance of accurately deriving UV data products for robust CDOM estimation. VIS algorithms, conversely, were shown to exhibit decreased robustness to variability in OAC relationships relative to non-visible – or invisible (INV) – algorithms (Hooker et al., 2020, 2021b; Houskeeper et al., 2020, 2021, 2022).

In February 2024, the launch of the Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) satellite mission provided state-of-the-art hyperspectral and spectrally extended oceanographic observing using the Ocean Color Instrument (OCI). Hyperspectral observations obtained using OCI are anticipated to reveal new opportunities for satellite observing of the microbial community (Cetinić et al., 2024), but the fidelity of inversions relating hyperspectral features to biological parameters is mostly unknown. The challenges stem, in part, from the use of legacy in situ instrumentation, e.g., wherein COTS spectrometers lacking paired microradiometers – i.e., hybrid-spectral configuration – to mitigate low sampling rate and SNR (Hooker et al., 2022; Houskeeper et al., 2024; Kudela et al., 2024) preclude most INV data products and do not support adherence to an absolute radiometric scale. For example, non-hybrid-spectral data products frequently rely on a null long-wavelength bias. Reporting of observations at resolutions finer than the nominal resolution of the spectrograph – often in addition to spectral smoothing – further challenges assessments of spectral content and independence of adjacent waveband observations.

Nonetheless, assessments of optical information content of the sunlit aquatic environment are ongoing and indicate adequate derivation of independent OACs is plausible – but not guaranteed. Hooker et al. (2021a) assessed spectral variability of global conservative plus non-conservative water bodies, with the latter modified by internal processes rather than linear mixing of parent or source waters, and demonstrated that a finite number of spectral modes adequately describe the majority of spectral variability across globally representative (and including case-2) waters. Cael et al. (2023) applied an eigenanalysis to legacy smoothed and interpolated hyperspectral observations and demonstrated two components capture over 80 % of the variance – with a similar information constraint also previously shown for hyperspectral observations of particulate absorption (Cael et al., 2020). Prochaska and Frouin (2025) applied a Bayesian framework testing whether hyperspectral observations could overcome an acknowledged degeneracy in radiative transfer (wherein multiple combinations of spectral absorption and backscattering are possible and reasonable) and reasoned that hyperspectral ocean color inversions are ill-posed.

Quantifying information content using eigenanalyses can be challenging: dimensionality bias can produce an incomplete perspective of information content by decreasing the variance captured by spectrally broadband versus spectrally narrow features. The latter – when not spectrally diluted – often provide informative and exploitable information (Houskeeper et al., 2020). Furthermore, despite known challenges in finite information content and spectral degeneracy, new information sources have steadily been discovered in aquatic optics throughout the field's long history. Legacy satellite observations of phytoplankton fluorescence added new biological observables complementary to (but distinct from) OC observations, while also minimizing vulnerabilities to atmospheric correction (Letelier and Abbott, 1996). Hyperspectral airborne investigations leveraged spatially cohesive optical signatures associated with harmful algal bloom events (e.g., Lee and Carder, 2004; Kudela et al., 2015) and retrieved macroalgal physiological condition (Bell and Siegel, 2022). Mixture density networks (MDNs) and inverse models that incorporate prior knowledge (e.g., leveraging a Bayesian framework) improve the management of degeneracy in radiative transfer and are forthcoming for PACE science objectives (O'Shea et al., 2021; Erickson et al., 2020, 2023).

A fundamental and timely question remains how to best leverage PACE observations within the trajectory of ocean color satellite missions to improve robust quantification of phytoplankton biomass and support characterization of biogeochemical changes in the sunlit ocean. Advancing this trajectory – wherein confounding effects of CDOM have been reported for at least over seven decades (e.g., Jerlov, 1951; Kalle, 1966; Yentsch, 1960; Sathyendranath et al., 1987) and recent studies suggest high spectral autocorrelation and degeneracy (Hooker et al., 2021a; Cael et al., 2023; Prochaska and Frouin, 2025) – requires improvements in understanding the primary sources of marine spectral variability. Aquatic spectral variability is assessed herein for natural aquatic environments using three independent bio-optical datasets, and analyses are replicated using the datasets' oligotrophic subsets, as applicable. Associations relating spectral variability to contemporaneous OAC field observations indicate the relative importance of individual constituents in modifying AOP spectral shapes. The key findings document, as follows: continuing challenges in retrieving C_a using observations constrained to the VIS domain; opportunities to advance independent OAC retrieval using spectrally expansive observations; greater variability in OAC relationships for oceanic environments than formerly ascribed; and primacy of CDOM absorption in spectrally modifying AOPs. The outcomes reflect a greater complexity of OAC dynamics consistent with advancing knowledge of the microbial loop, wherein new trophic pathways have been discovered contemporaneously with the aforementioned advances in optical oceanography (Azam et al., 1983; Azam, 1998; Fenchel, 2008).

2 Methods

Aquatic observations presented herein correspond to published datasets of above- and in-water observations of the emergent aquatic light field. The datasets rely on instrumentation, including legacy and state-of-the-art (SOTA) radiometers and spectrometers, corresponding to a technological trajectory spanning decades of improvements applicable across the hardware, software, and processing domains. For example, legacy oceanic radiometric observations were obtained using in-water optical instrument suites consisting of primarily analog radiometers mounted on rapidly descending, so-called rocket-shaped, profilers (Morel et al., 2007). Rocket-shaped profilers confer coarse vertical sampling resolution (VSR), deep extrapolation intervals, and a limited spectral range of observations. These deficiencies prevent observations in shallow, highly turbid, or high-flow-rate waters, meanwhile introducing uncertainties including, but not limited to, those arising from depth aliasing, expanded extrapolation distances, and low signal (Hooker et al., 2002). The Compact-Optical Profiling System (C-OPS) introduced a novel kite-shaped backplane (Morrow et al., 2010) to slow the rate of descent of the profiler while ensuring planar stabilities, along with corresponding improvements in microradiometer technology to advance sampling rates from usually less than 1 Hz for legacy spectrometers to routinely 15 Hz, with improved SNR and dynamic range. The addition of digital thrusters, i.e., the Compact-Propulsion Option for Profiling Systems (C-PrOPS), ensures planar stability at the initiation of a cast and mitigates adjacency effects by supporting navigation away from a ship or dock (Hooker et al., 2018a). The C-OPS with C-PrOPS technologies greatly expand the spectral range of in-water observations and reduce uncertainties in optical data products relative to legacy observations, in part, by supporting improved VSR (as fine as 0.9 mm) to enable accurate derivation of the central tendency of the wave-focusing field (e.g., Zaneveld et al., 2001) and reducing the initial depths of upwelling radiance observations, L_u(z,λ), to as shallow as 0.3 m, i.e., approximately the length of the downward-pointing radiance radiometer (Hooker et al., 2020). Improvements in sampling rates support management of high-frequency, non-Gaussian variability in flux observed by an above- or in-water instrument (due to glint and wave focusing, respectively), thereby expanding the spectral range of optical data products to preserve INV information and supporting necessary corrections to retain information associated with spectral signal amplitudes or brightness (Houskeeper and Hooker, 2023; Houskeeper et al., 2024).

SOTA advances in hyperspectral instrumentation include hybrid-spectral sensing configurations (Hooker et al., 2018c, 2022), wherein microradiometer and spectrograph observations are obtained in concert, with the multispectral radiometer providing necessary quality control for the hyperspectral spectrograph. The quality control is desirable, in part, because COTS spectrographs suffer from slower integration times, narrower dynamic range, and a degraded SNR relative to COTS microradiometers (Houskeeper et al., 2024; Kudela et al., 2019, 2024). In addition, a radiance control arm positions the L_u(z,λ) aperture near the water surface and approximately aligned with the upward-pointing irradiance radiometer (Hooker et al., 2018a). Improvements in the number of spectrograph pixels to as high as 2048, denoted λ₂₀₄₈, with 3 nm resolution is in keeping with PACE mission requirements for vicarious calibration exercises (Hooker et al., 2012a). For comparison, spectral resolution of legacy spectrometers include λ₂₅₆ spectrographs sampling the visible domain as often as every 3 nm, although bandwidths can be on the order of 10 nm resolution (Sea-Bird Scientific, 2024), i.e., comparable to COTS microradiometers but conferring lower SNR, dynamic range, and sampling rates. Legacy spectrometer observations are often interpolated to synthetic 1 nm intervals – which eases archiving requirements for dissimilar spectral configurations but over-reports the spectral resolution – and smoothed, e.g., using a 5 nm moving mean, which removes spectral information. Spectral smoothing and interpolation creates challenges for comparisons of spectral resolution between instruments.

2.1 Datasets

Three AOP datasets are considered herein, with biogeochemical parameters included when available. No dataset includes SOTA hybrid-spectral observations, because hybrid-spectral observations are presently too sparse in terms of the number of observations. Naming conventions for the datasets correspond to the applicable published literature, and their descriptions plus quality control are briefly summarized, as follows:

• RSE2007. Oceanic ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ observations obtained using primarily analog free-falling instrument suites during open-ocean field campaigns, including a minority of stations in the coastal zone, are described by Morel et al. (2007). Three stations (0.5 % of the dataset) corresponding to surface mucilage and coccolithophore blooms were removed, and observations were thinned so that only one spectrum was retained per station, although at least three spectra were generally obtained per station. The dataset is defined as RSE2007 and consists of contemporaneous observations of C_a and ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ spanning the VIS domain.

• RSE2021. Global ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ observations obtained using the digital C-OPS with C-PrOPS instrument suite deployed at oceanic (31.2 %), coastal (36.6 %), and inland (32.2 %) water sites are presented by Houskeeper et al. (2021). Optical observations featured average VSR of 6.0 mm (maximum VSR of 0.9 mm) and initial L_u(z,λ) observations routinely obtained at 0.3 m, i.e., the depth corresponding to the length of the L_u(z,λ) radiance radiometer. Contemporaneous observations of C_a and the absorption coefficient of CDOM at 440 nm, a_CDOM(440), supported global algorithm development (Hooker et al., 2020, 2021b; Houskeeper et al., 2021, 2022) and confirmed expansive ranges for OACs. For example, observed a_CDOM(440) spanned 0.001–2.305 m⁻¹, and the spectral slope of CDOM spanned 0.0095–0.0410 nm⁻¹. Conservative versus non-conservative designations were objectively determined prior to sampling (Hooker et al., 2020), and only the conservative observations are used herein. The dataset is defined as RSE2021 and consists of contemporaneous observations of C_a, a_CDOM(440), and ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ spanning the INV plus VIS domains.

• RSE2022. Oceanic hyperspectral observations with concurrent samples of phytoplankton pigments are described in Kramer et al. (2022) and accessed via Kramer et al. (2021). Observations span the open and coastal oceans, with C_a ranging from 0.019–4.150 mg m³. Optical observations were obtained using rocket-shaped profilers and buoys equipped with legacy COTS spectrometers, including systems with 10 nm bandwidths (Sea-Bird Scientific, 2024). Vertical profiles of measurements obtained deep in the water column were averaged to 2 m or coarser depth bins and extrapolated to the surface (Taylor et al., 2011). The compiled dataset overcomes differences in spectral configuration between various COTS systems by spectral interpolation to 1 nm resolution and mitigates optical data product artifacts associated with legacy profiling systems by subsetting spectra to VIS wavelengths, smoothing using a 5 nm moving mean, and additionally removing 19 % of the observations by visual inspection (Kramer et al., 2022). The dataset is defined as RSE2022 and consists of contemporaneous observations of C_a and ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ spanning the VIS domain.

The field campaign supporting the RSE2021 dataset (Hooker et al., 2020; Houskeeper et al., 2021) was designed to capture maximal range in global OAC variability, including the clearest waters to the most sedimented, eutrophic, or humic-rich waters. The clearest waters were observed, e.g., near Kawaihae, HI, and in Crater Lake, OR, both of which corresponded to an a_CDOM(440) value of 0.001 m⁻¹, approaching so-called pure seawater. The turbidity limit is inherently undefined, e.g., there is no known natural maximum to turbidity, but the RSE2021 dataset greatly exceeds the turbidity range of the RSE2007 and RSE2022 datasets because it contains significant representation from inland waters. The absolute radiometric data products from RSE2021, therefore, constitute a reasonably expansive range in optical properties such that open and coastal ocean optical data products (collected with legacy systems and with degraded SNR and VSR relative to the RSE2021 dataset) that lie substantially outside the range of global conservative data products indicate likely non-physical or biased observations. Quality control based on the threshold removal of non-physical radiometric brightness was thus applied, albeit only using VIS data products (given that RSE2007 and RSE2022 are restricted to the VIS domain), following sensitivity testing. Briefly, more stringent quality control metrics (e.g., requiring observations to comply with the brightness ranges of the globally representative RSE2021 dataset) produced similar eigenfunction shapes but reduced the size of RSE2022 by approximately 50 %. Alternately, relaxed quality control metrics permitting radiometric values up to 50 times darker than the global dataset did not alter the primary findings of RSE2022, and no threshold relaxation affected the primary findings of RSE2007. Considering acknowledged uncertainties relating to hardware, processing, data acquisition, and water mass differences, the quality control thresholds were performed within this insensitivity range by filtering: RSE2022 observations radiometrically darker than the darkest RSE2021 observations (by at least approximately 50 %); and RSE2007 observations radiometrically darker or brighter than RSE2021 observations (approximately) at corresponding green (555 nm) or red (683 nm) wavebands, respectively. The brightness filtering resulted in the removal of 18 of 590 observations (3 %) from RSE2007 and 22 of 144 observations (15 %) from RSE2022. Because both RSE2007 and RSE2022 include observations obtained using rocket-shaped profilers, the difference in compliance is likely due to differences in the optical sensing technologies – i.e., COTS spectrometers (RSE2022) are more degraded by slow sampling rates, reduced dynamic range, and low SNR (Houskeeper et al., 2024; Kudela et al., 2024). The final number of compliant observations for the RSE2007, RSE2021, and RSE2022 datasets are 572, 612, and 122, respectively. The quartile summary of the biochemical quantities corresponding to each dataset is shown in the Supplement.

2.2 Biogeochemical algorithms

OC algorithms are derived for each dataset using the OC4 formulation (O'Reilly and Werdell, 2019), a fourth-order polynomial with coefficients configured to match the spectral configurations of the SeaWiFS mission and others, as follows:

$$\begin{array}{}\text{(1)}& C_{\mathrm{a}}=a_{\mathrm{0}}+\sum _{n=\mathrm{1}}^{\mathrm{4}}{a}_n{\left({\text{log}}_{\mathrm{10}}{\displaystyle \frac{R_{\mathrm{rs}}\left({\mathit{\lambda }}_{\mathrm{b}}\right)}{R_{\mathrm{rs}}\left(\mathrm{555}\right)}}\right)}^n\left[\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}\right],\end{array}$$

where the a terms are fitting constants, n is the order index, and λ_b indicates the wavelength corresponding to the radiometrically brightest R_rs(λ) observation at applicable blue SeaWiFS wavelengths (i.e., 443, 489, and 510 nm). The OC4 formulation was selected to ensure maximal applicability to RSE2007, which supported the fewest spectral channels, while maintaining a consistent formulation across datasets. Nonetheless, other OC formulations are anticipated to produce similar results.

The RSE2021 dataset contains contemporaneous observations of a_CDOM(440), but contemporaneous a_CDOM(440) observations were not contained in the other datasets. Consequently, an EMA algorithm derived using both case-1 and case-2 waters to minimize biases in algorithm fitting corresponding to regional OAC covariances (Houskeeper et al., 2021) is applied to the RSE2007 and RSE2022 datasets to estimate a_CDOM(440). The approach leverages ratios of ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$, hereafter ${\mathrm{\Gamma }}_{{\mathit{\lambda }}_{\mathrm{2}}}^{{\mathit{\lambda }}_{\mathrm{1}}}$, using the most spectrally separated wavelengths, which mitigates confounding signals from internal VIS wavelengths to improve the independence of a_CDOM(440) retrieval from C_a variability (Houskeeper, 2020; Houskeeper et al., 2021). The EMA formulation for above-water, spectrally expansive observations is modeled as a power-law relationship, as follows:

$$\begin{array}{}\text{(2)}& a_{\mathrm{CDOM}}\left(\mathrm{440}\right)=A{\left[{\mathrm{\Gamma }}_{{\mathit{\lambda }}_{\mathrm{2}}}^{{\mathit{\lambda }}_{\mathrm{1}}}\right]}^B\left[{\mathrm{m}}^{-\mathrm{1}}\right].\end{array}$$

The λ₁ and λ₂ values correspond to the shortest and longest wavelengths, respectively, supported by the given dataset. For the RSE2007 and RSE2022 datasets, wherein INV observations were not supported, the applicable wavebands are 412 and 670 nm, and the A and B terms are 0.242 and −0.961, respectively. EMA algorithm performance improves with expanding spectral range (Hooker et al., 2020; Houskeeper et al., 2021), but robust performance of VIS algorithms has been documented for ocean color observations (Houskeeper et al., 2021; Hooker et al., 2021b), and simulated geostationary weather satellite observations (Houskeeper et al., 2022). The uncertainty of EMA, however, is anticipated to increase when applied to optical data products obtained using spectrometers and rocket-shaped profilers (e.g., RSE2022) due to increasing extrapolation depths combined with strong attenuation of signal, particularly within the end-member spectral domains (Kudela et al., 2019, 2024; Hooker et al., 2020; Houskeeper et al., 2021, 2024).

2.3 Spectral comparisons

The dimensionality for each of the RSE2007, RSE2021, and RSE2022 datasets is reduced using an eigenanalysis of normalized observations of the square root of ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$, hereafter ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\text{N}}^{{\mathrm{0.5}}^{\prime }}$, because the square-root transformation was found to improve normality in the blue spectral domain based on the Shapiro–Wilk test statistic, W. Nonetheless, sensitivity testing was performed by reproducing the eigenanalyses using natural-log transformation, as well as no transformation. The eigenanalysis captures linear patterns of variability in the datasets and is performed as follows:

$$\begin{array}{}\text{(3)}& \mathbf{C}\mathit{\psi }=\mathbf{\Lambda }\mathit{\psi },\end{array}$$

in which C is the covariance matrix of the normalized ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\text{N}}^{{\mathrm{0.5}}^{\prime }}$ values for each dataset, with square dimensions corresponding to the number of wavelengths in each dataset. The ψ value denotes the eigenfunction matrix (matrices are denoted by bold symbols). The diagonal matrix Λ contains scalar information corresponding to the variance captured by each eigenfunction. The mean spectra from each dataset, plus reconstructed spectra indicating the effects of the leading and secondary eigenfunctions, are shown in the Supplement.

Scalar amplitude functions quantifying the stretching and compressing necessary to represent the ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\text{N}}^{{\mathrm{0.5}}^{\prime }}$ values in the coordinates defined by ψ are derived by projecting the normalized optical observations onto ψ, as follows:

$$\begin{array}{}\text{(4)}& \mathbf{S}={\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\text{N}}^{{\mathrm{0.5}}^{\prime }}\mathit{\psi },\end{array}$$

in which S is the scalar amplitude function matrix with columns, S_i, corresponding to scalar amplitude functions for individual eigenfunctions, ψ_i (where i is the eigenfunction index). The length of each S_i column is the number of observations in the underlying dataset.

Comparisons between available biogeochemical parameters – e.g., field plus algorithm values of C_a and a_CDOM(440) – and the S matrices are derived using Pearson's correlation coefficient, ρ_i,x, where i is the eigenfunction index and x represents the biogeochemical quantity under comparison with the corresponding column i of the S matrix. Because biogeochemical parameters covary, e.g., C_a and a_CDOM(440) are strongly correlated (Morel and Prieur, 1977), partial correlation coefficients (Fisher, 1924) are derived using the notation ${\mathit{\rho }}_{i,x|y}$ to indicate correlation of S_i with biogeochemical quantity x after adjusting for covariance with the biogeochemical quantity y. The subscript representations used herein are g (for the legacy linkage to gelbstoff absorption), C_a, EMA, and OC4, which correspond to field observations of a_CDOM(440), field observations of C_a, EMA algorithm observations of a_CDOM(440), and OC algorithm observations of C_a, respectively. The R² value is also utilized when quantifying variance in biogeochemical field observations captured by a remote sensing algorithm. Pearson's correlation coefficients comparing C_a and a_CDOM(440) – measured or derived algorithmically – for the RSE2021, RSE2022, and RSE2007 datasets are 0.812, 0.912, and 0.909, respectively.

3 Results

Eigenanalyses indicate that the first, second, and third eigenfunctions capture 60 %, 32 %, and 5 % (± 1 %), respectively, of the variance within each dataset. The decomposition of the majority of the variance in the aquatic spectral dataset into three functions is consistent with the findings of Hooker et al. (2021a) and Cael et al. (2023), as follows: Hooker et al. (2021a) demonstrated that spectral changes corresponding to increasing optical complexity were represented using only five spectral modifications – two of which are similar but in opposing directions and would therefore correspond to a single eigenfunction, plus one which is characterized by the emergence of multiple peaks and corresponds to severely complex waters; and Cael et al. (2023) showed the aquatic optical spectra could be compressed into three (or four) principal components with minimal information loss. The Cael et al. (2023) dataset includes significant representation by RSE2022, and similar results obtained herein using RSE2022 provide support for the normalizations and methodologies being comparable.

Figure 1 Spectral representations of the eigenfunction matrices corresponding to the RSE2021, RSE2022, and RSE2007 datasets are represented as blue, orange, and green lines, respectively. The first three eigenfunctions ψ₁–ψ₃ are organized vertically based on descending variance explained, Λ. The y axis label superscripts indicate the phase, i.e., either positive (+) or negative (−), in which the eigenvectors are presented, and y axis label subscripts indicate the ordering of the eigenvector. Correlation coefficients relating S for each eigenfunction to recorded environmental (g and C_a) or algorithm (EMA and OC4) values are indicated using ρ notation, wherein the left subscript indicates the parameter under comparison with S, given the parameter indicated by the right subscript, when applicable. Eigenfunctions for oligotrophic subsets (C_a less than 0.5 mg m⁻³) of the datasets are plotted on the same axes using gray lines. For clarity, all statistics presented correspond to the complete datasets (colored lines), and the sign for ρ corresponds to the phase in which ψ is plotted. The nominal locations of the blue and green wavebands used by the OC4 algorithm are indicated by blue and green vertical dashed lines, respectively.

3.1 Eigenfunction characteristics

The three eigenfunctions capturing maximal variance for the RSE2021, RSE2022, and RSE2007 datasets are presented as a function of λ in Fig. 1. Similarities in the eigenfunction spectral shapes are apparent across datasets for each of ψ₁, ψ₂, and ψ₃, and similar shapes are likewise produced by the eigenanalyses for the oligotrophic (C_a less than 0.5 mg m⁻³) datasets (gray). Although the tertiary eigenfunction explains only approximately 5 % of the variance for each dataset, the spectral shapes represented by ψ₃ are consistent across all datasets. The similarities in ψ are summarized as follows:

• ψ₁: reversed sign for the anomaly between shorter (UV/blue) and longer (red/NIR) wavelengths;

• ψ₂: maximum amplitude for the anomaly in the domain (i.e., a peak), plus the emergence of spectral features in the spectral vicinity of the C_a fluorescence peak (except for RSE2007, which lacks the requisite spectral resolution in the applicable wavelength domain); and

• ψ₃: reversed sign for the anomaly between interior VIS and exterior VIS plus INV wavelengths.

Individual eigenfunctions are not anticipated to be singularly associated with variability in a specific OAC (e.g., because OACs covary), but the spectral shapes of the eigenfunctions may nonetheless be assessed qualitatively in light of known spectral properties of individual OACs to identify similarities to spectral properties of OACs. The spectral shapes of the eigenfunctions are primarily considered as a function of absorption processes, which generally – notwithstanding absorption band effects (Zaneveld and Kitchen, 1995) – confer stronger spectral dependencies than scattering processes (Kirk, 2011). Differences in brightness or peak height, however, could produce reversals in the VIS and INV anomalies, possibly relevant to ψ₃, but are not assessed herein.

Spectral features associated with C_a absorption – e.g., a relative maximum or minimum in the blue and red domains (Bridigare et al., 1990), or a fluorescence peak (Letelier and Abbott, 1996) – are most clearly visible in ψ₂. For example, although both ψ₁ and ψ₂ indicate spectral dependencies in the vicinity of 670 nm for RSE2022, the amplitude of the anomaly (as a percentage of the range expressed in ψ) is on the order of 10 % versus 50 % for ψ₁ and ψ₂, respectively. Similarly, all datasets indicate a blue- or green-domain peak (or trough) in ψ₂ of greater magnitude than in ψ₁. Conversely, spectral darkening of shorter wavelengths with minimal spectral features is most consistent with CDOM absorption (Kalle, 1966; Jerlov, 1968; Siegel et al., 2013), which approximately describes ψ₁. The similarity of spectral shapes for ψ between the complete and oligotrophic (gray) datasets – as well as the differences in the range and mean properties of the water bodies sampled across the three datasets – indicates that spectral variability is consistent across datasets and not primarily a function of differences in sampling within each dataset.

The spectral shape comparisons, while qualitative, can be verified using partial correlation coefficients comparing the S amplitude functions and corresponding values of C_a and a_CDOM(440), while controlling for a_CDOM(440) and C_a, respectively. RSE2021 includes field observations of a_CDOM(440) and C_a, whereas RSE2022 and RSE2007 include field observations of C_a but not a_CDOM(440). Consequently, the RSE2021 results are considered first. Briefly, partial correlation coefficients for S₁ indicate maximal correlation to a_CDOM(440) values (${\mathit{\rho }}_{\mathrm{1},g|C_{\mathrm{a}}}=\mathrm{0.80}$) compared to C_a values (${\mathit{\rho }}_{\mathrm{1},C_{\mathrm{a}}|g}=\mathrm{0.11}$), whereas partial correlation coefficients for S₂ indicate maximal correlation to C_a values (${\mathit{\rho }}_{\mathrm{2},C_{\mathrm{a}}|g}=-\mathrm{0.22}$) compared to a_CDOM(440) values (${\mathit{\rho }}_{\mathrm{2},g|C_{\mathrm{a}}}=\mathrm{0.02}$). Both of the maximal correlation results are highly significant (P<0.01), whereas the non-maximal comparisons are not (P>0.18). These results are in keeping with the qualitative assessments of the spectral shapes of ψ, in which ψ₁ was shown to produce spectral dependencies more consistent with that of CDOM as documented in the literature, whereas ψ₂ was shown to produce spectral dependencies with features most closely matching the spectral domains of C_a absorption or fluorescence properties. The sign of the significant partial correlation coefficients is also consistent with the ψ phase and the spectral shape of absorption for a_CDOM(440) and C_a. For example, negative anomalies in the blue and UV indicated by S₁ produce a positive correlation with a_CDOM(440) values, whereas positive (peak-like) anomalies in blue wavelengths indicated by S₂ produce a negative correlation with C_a values.

Estimates of a_CDOM(440) are also obtained for RSE2022 and RSE2007 using EMA. The partial correlation values derived using the EMA product for the RSE2022 dataset are very similar to those derived using the RSE2021 dataset. Briefly, partial correlation coefficients for S₁ indicate maximal correlation to a_CDOM(440) values, whereas partial correlation coefficients for S₂ indicate maximal correlation to C_a values. Both of the maximal correlation results are highly significant (P≤0.01), whereas the other comparisons are not (P>0.38). The results from the RSE2007 dataset are different from those of RSE2022 and RSE20221 as follows: partial correlation coefficients for S₁ still indicate maximal correlation to a_CDOM(440) values, but correlation to C_a is also significant, albeit lower; and partial correlation coefficients for S₂ indicate slightly higher correlation to a_CDOM(440) values compared to C_a values, and both are highly significant (P<0.01). The differences are consistent with decreasing independence of EMA to C_a variability with decreasing spectral range (Houskeeper et al., 2021), plus maximal representation of open-ocean waters in RSE2007.

Correlation tests comparing the S functions with OC algorithm values find the strongest relationships for all datasets correspond to S₁, ranging from 0.89 to 0.96. Less strong, albeit significant, relationships are found for S₂, ranging from −0.38 to 0.16. The significant results for comparing OC algorithm values with S₂ are likely due to the test not accounting for covariance in C_a and a_CDOM(440). For example, the correlation coefficients using RSE2021 – the highest-quality dataset and the only dataset to contain contemporaneous C_a and a_CDOM(440) observations – still indicate strong significant relationships when comparing OC algorithm values and S₁ while controlling for C_a, but indicate weak and insignificant relationships when comparing OC algorithm values and S₂ while controlling for a_CDOM(440).

In summary, analyzing the shapes of the eigenfunctions and comparing correlation matrices using the eigenfunction amplitude functions and contemporaneous OAC field observations yields five observations, as follows: (a) spectral shapes of the eigenfunctions for an individual dataset are distinct from one another and characterized by maxima and minima consistent with OAC absorption properties described in the applicable literature; (b) the transitions in spectral shapes for the first, second, and third eigenfunctions are similar across all datasets and the oligotrophic subsets of those same datasets; (c) the interpretation of the eigenfunction shapes based on the spectral positioning of observed minima and maxima is consistent with the results of the partial correlation tests (e.g., ψ₁ confers spectral similarities with CDOM absorption and the corresponding amplitude function produces highly significant correlation to CDOM absorption); (d) the OC algorithm produces stronger correlation to S₁ than S₂; and (e) S₁ produces stronger correlation to field observations of a_CDOM(440) than C_a values. Although RSE2007 indicates relatively more OAC covariance than RSE2021 and RSE2022, as anticipated, the key findings are nonetheless consistent.

Figure 2 Algorithmic relationships for the RSE2021, RSE2022, and RSE2007 datasets are represented as blue, orange, and green open circles, respectively. Panels (a)–(c) indicate relationships between OC algorithms and field C_a, and a 1:1 correspondence line is shown in black. Panels (d)–(f) indicate relationships between OC algorithms and field (d) or algorithmic (e, f) observations of a_CDOM(440), with R² values shown. Algorithmic products (using an EMA relationship following Houskeeper et al., 2021) were required for the RSE2022 and RSE2007 datasets because only RSE2021 included routine field sampling of a_CDOM(440).

3.2 OC algorithm retrievals

Section 3.1 documents that the primary mode of variability in ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ for all datasets is maximally correlated with a_CDOM(440) and OC algorithm retrievals (Fig. 1). The observed linkage between OC algorithm sensitivity and a_CDOM(440) variability is further investigated by comparing OC algorithm retrievals with field observations of C_a, plus field or algorithmic a_CDOM(440) observations (Fig. 2), as applicable. The coherence of a relationship between the OC algorithm and C_a observations is assessed using R², and differences across datasets in the coherences observed are consistent with differences in the oceanic contributions of each dataset, as follows: (a) the least coherence is associated with RSE2021, which contains similar representation of oceanic, coastal, and inland waters, and therefore consists of a high percentage of optically complex water types; (b) intermediate coherence is associated with RSE2022, which contains coastal and inland water observations and also consists of data products obtained with coarse VSR, low SNR, and deep sampling intervals far from the surface; and (c) the most coherence is associated with RSE2007, which corresponds to primarily open-ocean waters and was obtained with technology most similar to that used in the earliest OC algorithm derivations.

The coherence of the relationship between the OC algorithm and a_CDOM(440) observations is also assessed using R², and improvements in R² are observed for the a_CDOM(440) comparison relative to the C_a comparison. The magnitude of the improvements are likewise consistent with the differences in the oceanic contributions of each dataset, as follows: (a) maximal improvement in coherence is indicated by RSE2021; (b) intermediate improvement in coherence is indicated by RSE2022; and (c) minimal improvement in coherence is indicated by RSE2007. Although the differences observed are slight for RSE2007 (indicating the dataset adheres well to case-1 waters), all datasets nonetheless produce stronger R² values when comparing the OC algorithm values to a_CDOM(440) versus C_a. The findings are consistent with those of Sect. 3.1, in which the OC algorithm products produced maximal correlation to ψ₁, which in turn was maximally associated with a_CDOM(440) – and not with C_a.

Figure 3 Variance explained by the first six eigenfunctions (ψ₁–ψ₆) for eigenanalyses using spectral subsets of the RSE2021 dataset. The full spectral range (320–875 nm) is shown in blue, matching the presentations of Fig. 1 and 2, and denoted UVN21; the visible (412–683 nm) spectral subset results are shown in yellow and denoted VIS21; and the invisible (313–395 and 710–875 nm) spectral subset results are shown in red and denoted INV21. Pearson's correlation coefficients for a_CDOM(440) and C_a (ρ_1,g and ${\mathit{\rho }}_{\mathrm{1},C_{\mathrm{a}}}$, respectively) are denoted using arrows for the VIS21 and INV21 analyses, with bold font indicating significance (P<0.05). The inlay panel shows that RSE2007, RSE2021, and RSE2022 datasets have similar maximal variance explained in ψ₁ and ψ₂.

3.3 Information content of invisible and visible data products

Information content as a function of spectral range is further investigated using the RSE2021 dataset, which confers the most expansive spectral range (313–865 nm). Scree plots indicating variance captured for eigenfunctions ψ₁–ψ₆ are shown for three spectral partitions of the RSE2021 dataset in Fig. 3, as follows: UVN21 corresponds to the full spectral range of RSE2021 data products (i.e., UVN21 is identical to RSE2021); INV21 corresponds to the invisible data products of RSE2021 (i.e., wavelengths shorter than 400 nm or longer than 700 nm); and VIS21 corresponds to the visible data products of RSE2021 (i.e., wavelengths within 400–700 nm). Pearson's correlation coefficients in Fig. 3 compare the S₁ and S₂ amplitude functions with C_a and a_CDOM(440) values. The comparisons explore the results of Fig. 1, wherein ψ₁ and ψ₂ are found to confer information likely most related to a_CDOM(440) and C_a, respectively.

Variances captured by the UVN21, INV21, and VIS21 datasets decrease rapidly with increasing eigenfunction number. Considering the UVN21 scree plot as a baseline, ψ₁ captures similar variance and ψ₂ captures more variance for the VIS21 dataset, whereas ψ₁ captures more variance and ψ₂ captures less variance for the INV21 dataset. Information content of INV21 is thus compressed more fully into a singular eigenfunction, which captures nearly 78 % of the variance in the dataset. The INV21 leading eigenfunction is very strongly associated with variability in a_CDOM(440), and the correlation coefficient comparing S₁ and a_CDOM(440) for the INV21 dataset, −0.96, approaches that of a perfect inverse correlation. The INV21 secondary eigenfunction indicates minimal adherence to either a_CDOM(440) or C_a, with correlation values of −0.07 and −0.09, respectively. Very strong correlation of ψ₁ with a_CDOM(440), no correlation of ψ₂ with C_a, and very low explanatory power of the subsequent eigenfunctions for INV21 indicate photosynthetic pigmentation confers minimal spectral modifications in the INV domain besides those covarying with a_CDOM(440), consistent with previous EMA assessments (Hooker et al., 2020; Houskeeper et al., 2021). Inorganic particles, not measured in any of the datasets assessed herein, modify INV data products (e.g., Doron et al., 2011) and may contribute to the variance captured in INV21 ψ₂. Inorganic particles are not assessed herein due to the sparsity of independent particle field observations contemporaneous and coregistered with the AOP observations.

Information content of VIS21 is less compressed into a singular eigenfunction (although the leading eigenfunction still dominates), and the correlation coefficients comparing ψ₁ and ψ₂ to a_CDOM(440) and C_a are more similar than for the INV21 dataset. Similar correlation coefficient values indicate spectral modifications associated with a_CDOM(440) and C_a are less readily partitioned using only VIS data products, despite the high SNR of the RSE2021 dataset. The result suggests independent derivation of C_a is more challenging using VIS observations in the absence of UV and NIR information. Spectrally expansive (UV, VIS, and NIR) observations thus support more accurate and independent retrievals of a_CDOM(440) compared to VIS observations wherein a_CDOM(440) and C_a signals are maximally confounded.

4 Discussion

4.1 Historical formulations invoke primacy of chlorophyll a

The generalized transition from bluer oceanic to greener coastal waters is routinely mapped with OC algorithms to increases in the near-surface abundance of oceanic phytoplankton (McClain et al., 2006). The theoretical basis for OC algorithms is the absorption of blue light by biogenic compounds (Gordon et al., 1980; Gordon and Morel, 1983), including C_a but also other covarying OACs, e.g., a_CDOM(440), various secondary photopigments, and detritus. OC algorithms and the ocean color activities supporting their development and validation historically pose C_a as the predominant determinant of oceanic optical properties, demonstrated in pertinent part, as follows: (a) OC algorithms regress ratios of oceanic reflectances to C_a only and do not include secondary photopigments or dissolved compounds (O'Reilly et al., 1998); (b) C_a field observations are more ubiquitous than a_CDOM(440) field observations (Werdell and Bailey, 2002); (c) spatiotemporal differences in OC algorithm performance are routinely mitigated using regional tunings, in which OC algorithms are refitted for a subset of global waters using regional field observations of R_rs and C_a – but not other parameters such as a_CDOM(440); and (d) primacy of C_a is consistent with the published literature regarding the modeling of optical properties of case-1 water masses. For example, Morel (1988) leverages C_a for quantifying oceanic living and detrital algal material, as well as the depth of the euphotic layer.

Bio-optical formulations to derive optical properties as a function of C_a do not explicitly require primacy of C_a. Rather, algorithms parameterize the evolution in optical properties – notwithstanding observational challenges – based on empirically observed relationships (Morel and Prieur, 1977; Morel, 2009). Nonetheless, variability in OAC relationships exists for in situ water bodies (Hansell and Orellana, 2021). Empirical approaches such as the OC family of algorithms must mitigate observational artifacts (Uitz et al., 2006; Kudela et al., 2019, 2024), plus regional variability and non-linearity in OAC relationships (Morel, 2009). Regional variability in OAC relationships manifests with the need to tune VIS algorithms for specific water bodies. For example, tunings for arctic waters account for higher a_CDOM(440) relative to C_a content (Matsuoka et al., 2014; Lewis and Arrigo, 2020), and, conversely, tunings for Antarctic waters account for lower a_CDOM(440) relative to C_a (Dierssen and Smith, 2000). Non-linear OAC relationships are in keeping with highly non-linear formulations of OC algorithms, e.g., the number of power terms used in the polynomial model to fit C_a to ratios of R_rs is routinely four (Morel, 2009; O'Reilly and Werdell, 2019). Linearity and log-linearity of CDOM algorithms were shown previously to be a function of the spectral separation of algorithm end members (Hooker et al., 2020, 2021a; Houskeeper et al., 2021). The results herein suggest that differences in linearity correspond, in part, to differences in the separability of signals associated with the spectral range.

4.2 Observations indicate primacy of CDOM over C_a

The eigenanalyses herein assess the primacy of C_a using three independent datasets, plus the oligotrophic subsets of each dataset. Phytoplankton absorption was not tested due to the limitations of the datasets. Although phytoplankton absorption is generally considered to correspond to C_a, processes such as pigment packaging can add variability relevant to the correlations presented herein (Bricaud et al., 2004). Nonetheless, despite significant differences in the spectral range and uncertainties of the datasets due to differences in hardware, software, and processing implemented, as well as the water bodies assessed plus variability associated with pigment packaging, the eigenanalyses show qualitative and quantitative similarities for all datasets and oligotrophic subsets (Fig. 1). The similarities favor primacy of a_CDOM(440) rather than C_a as the predominant driver of optical variability for waters, including those of the open ocean, summarized briefly, as follows: (a) the leading eigenvectors from each dataset capture approximately 60 % of the variance and indicate opposing anomalies for longer versus shorter wavelengths with minimal amplitude in the blue–green transition domain. The secondary eigenvectors from each dataset capture approximately 32 % of the variance and indicate internal VIS spectral dependencies characteristic of C_a absorption, including a maximum near the blue–green transition; (b) the S₁ term representing the stretching and compressing necessary to best represent the data using the leading eigenfunction is always more strongly associated with variability in a_CDOM(440) than C_a values. The S₂ term more strongly associates with field observations of C_a than a_CDOM(440) for two of the three datasets, likely due to reduced applicability of the EMA algorithm for the less spectrally expansive dataset; and (c) the signs of the partial correlation coefficients derived comparing the leading and secondary eigenfunctions to C_a and a_CDOM(440) variability are in all cases consistent considering the short wavelength anomaly in ψ₁ and the blue–green anomaly in ψ₂ compared to the spectral absorption profiles of CDOM and C_a. Observations of detrital absorption were not available, and detritus and CDOM share some general similarities in the spectral shape of absorption, i.e., absorption increases for both at shorter wavelengths. Some key factors relevant to the focus herein on primacy of a_CDOM(440) rather than detritus include the leading eigenfunction capturing approximately (within 1 % of) 60 % of the variance for each dataset, whereas each dataset corresponds to dissimilar contributions from inland, coastal, and open-ocean waters. Absorption by CDOM and detritus also exhibit variable correlation patterns (Babin et al., 2003; Twardowski et al., 2004), which would degrade the correlation between the leading eigenfunction and a_CDOM(440) if the leading eigenfunctions were most strongly associated with detritus. The comparisons between S₁ and a_CDOM(440) indicate log-linear relationships, which are shown – along with Spearman's rank correlation coefficient values – in the Supplement.

Although primacy of a_CDOM(440) is conceptually distinct from formulations in the published literature posing C_a as an independent variable primarily determining optical properties (e.g., Gordon, 1988; Morel, 1988), the results do not necessarily conflict with published observations in the bio-optical literature. For example, the inadequacy of multispectral VIS observations to partition C_a and a_CDOM(440) signals was noted at least six decades ago (Yentsch, 1960), and the case-1 prescription does not explicitly require primacy of C_a (Morel and Prieur, 1977). Rather, primacy of C_a is invoked by models representing the oceanic light field and its OACs as a function of C_a (Morel, 1988), plus the implementation of OC algorithms wherein regional tunings refit R_rs ratios and C_a to overcome regional differences in relationships between OACs. The principal finding wherein CDOM – rather than C_a – variability is most directly associated with optical data products emerging from the aquatic surface, e.g., ${\left[L_{\mathrm{W}}\left(\mathit{\lambda }\right)\right]}_{\mathrm{N}}$ or R_rs, may reflect differences in temporal and vertical dynamics of particulate versus dissolved constituents, including via the complex pathways dissolved materials cycle in the microbial loop (Azam et al., 1983) that would contribute to decoupling CDOM dynamics from C_a.

For example, the connectivity, directionality, and rates of fluxes connecting phytoplankton, dissolved organics, nutrients, and other parameters determine the degree to which one parameter could feasibly predict another. The case-1 model, wherein C_a predominantly modifies the aquatic light environment and predicts CDOM and other OACs, is perhaps most logically consistent with historical understandings of the so-called grazing food chain, a simplified oceanic trophic cascade wherein bacteria were considered negligible in the cycling of organic materials (Steele, 1974) and phytoplankton degradation produces CDOM (Kalle, 1966). Increasingly complex microbial loop dynamics – wherein phytoplankton contribute dissolved organic matter (DOM) directly through cellular lysis, excretion, and grazing, and wherein bacteria regulate the remineralization of DOM that fuels phytoplankton growth – support an increasing diversity of pathways representing the flow of energy and materials through marine systems (Azam et al., 1983; Azam, 1998). Additional discoveries, including new functional groups (e.g., viruses) and pathways (e.g., direct phytoplankton uptake of DOM), further complicate oceanic food webs (Granéli et al., 1999; Fenchel, 2008) relative to those of the grazing food chain.

Stratification and mixing modify vertical time scales of sinking and have been related to global-scale patterns in surface concentrations of DOM, e.g., wherein stratified mid-latitude waters are characterized by higher concentrations of DOM and well-mixed waters of higher latitudes – with the exception of the arctic – contain relatively lower concentrations of DOM (Hansell and Orellana, 2021). The optically active DOM subset, i.e., CDOM, is further influenced by the kinematic processes governing photobleaching and degradation rates, affected by mixed layer depth, seasonality, latitude, sea surface temperature, and sky condition (Gonsior et al., 2014; Del Vecchio and Blough, 2002). Changing ocean temperatures and acidity also relate to temporal dynamics in DOM via the restructuring of microbial food webs through altering stratification (and therefore nutrient availability), microbial metabolic rates plus nutrient demands, remineralization of DOM, and perhaps, species distributions (Sala et al., 2016). Horizontal spatial comparisons are not assessed, but the results are, perhaps, in keeping with recent findings that patchiness in oceanic CDOM is not well explained by patchiness in C_a and exhibits a distribution more similar to sea surface temperature (Gray et al., 2025). Vertical differences may also be relevant, although neither the depth dimension nor vertical flux are investigated here. Negative buoyancy of particles – e.g., sinking velocity increases proportionally with density differences and the squared particle radius (Bach et al., 2012) – or nutricline structure may drive particulate organic distributions away from the surface, i.e., away from where remote sensing signals are weighted (e.g., Morel and Berthon, 1989). DOM – including a refractory pool that also persists with lower sensitivity to the dynamics of microbial food webs (Jiao et al., 2010) – is transported primarily through advection with water masses.

The eigenanalyses presented herein additionally suggest case-1 conditions are not representative of global water bodies nor selective oceanic subsets, because the partial correlation coefficients indicate the eigenfunctions associate (to a greater extent than anticipated with individual OACs, and not with bulk combinations. For example, considering an idealized case-1 dataset wherein all OACs are accurately derived from C_a, the leading eigenfunction would be anticipated to produce similar rank correlation coefficients for each OAC. In contrast to expectations that water mass OACs should be modeled from C_a, the eigenfunction analyses indicate that the majority of the variance (approximately 60 %) in oceanic optical data products most strongly corresponds to a_CDOM(440), and an orthogonal mode capturing less variance (approximately 32 %) may correspond to residual variability in C_a. Finally, although the results suggest a case-1 model with a_CDOM(440) as the fundamental variable rather than C_a might better represent the three datasets and their oligotrophic subsets, the underlying assumption that a single OAC captures oceanic light variability would limit the model to representing only 60 % of the variance observed. The result suggests oceanic trophic cascades are characterized by greater complexity than can be represented by a univariate model.

Overall, the findings demonstrate high variance captured by the leading eigenfunction and the strong association of the leading eigenfunction with a_CDOM(440) variability, but the findings would be strengthened by adding field sampling for additional biochemical or bio-optical parameters. For example, RSE2021 demonstrates the importance of adding observations of a_CDOM(440) plus spectrally expansive radiometry to observations of C_a. Also adding observations of non-algal particles (NAPs) and detrital absorption would improve clarity on the eigenfunction interpretations and provide another option for investigating subsequent modes. The findings presented herein correspond to datasets spanning globally representative or oceanic waters and do not rank the importance of constituents in any specific region. For example, regional scenarios wherein high variability of NAPs is associated with riverine flux or resuspension would not be anticipated to reflect the findings herein supporting a_CDOM(440) primacy.

4.3 OC algorithm retrievals

High sensitivity of OC algorithms to regional differences in bio-optics is likewise consistent with the findings herein of S₁ producing maximal correlation to a_CDOM(440) and the OC algorithm products. For example, regional tunings of OC algorithms have proven necessary even at spatial scales as large as currents (Kahru and Mitchell, 1999), seas (Morel and Gentili, 2009), and oceans (Johnson et al., 2013), revealing an ever diminishing fraction of oceanic surface environments wherein OC algorithms are robust to natural variability in the mixtures and dynamics of OACs.

Sensitivity of the OC algorithm to a_CDOM(440) is supported by the R² metrics, which indicate the OC algorithm captures greater variance in a_CDOM(440) compared to C_a, consistent with the results of the eigenanalyses. The assessment warrants caution, as follows: a_CDOM(440) is less commonly obtained in bio-optical field sampling than C_a, which is why two of the datasets herein do not include contemporaneous a_CDOM(440) observations; and the RSE2021 dataset, which does include contemporaneous observations of both a_CDOM(440) and C_a, also includes representation from coastal and inland waters. Regarding the former, the high incidence of optical measurements with contemporaneous field observations of C_a but not a_CDOM(440) poses a persistent difficulty in aquatic optics, and the results herein indicate the importance of routine a_CDOM(440) field observations in bio-optical field activities. Regarding the latter, the water mass representation of RSE2021 is addressed to the extent practicable by considering the scatter of the data points in the separate quadrants of the Fig. 2 validation comparisons. For example, a hypothesis wherein the OC algorithm would produce maximal association to C_a in open-ocean waters would correspond to reduced scatter of the data points in the lower left quadrant in the C_a comparison (panel a) relative to the a_CDOM(440) comparison (panel d). Conversely, the Fig. 2 results show tighter association of the OC algorithm to a_CDOM(440) relative to C_a in the lower left quadrant.

Comparing between datasets, the RSE2007 observations produce the most similar R² values when comparing OC algorithm products to C_a and a_CDOM(440), consistent with the RSE2007 dataset being maximally represented by open-ocean waters far from coastal inputs and RSE2007 conferring the least expansive spectral range, which challenges the applicability of the EMA algorithm for retrieving a_CDOM(440) independently of C_a. Regarding the former, nearly identical agreement in R² comparing the OC algorithm products to C_a and a_CDOM(440) suggests RSE2007 waters are most consistent with the case-1 scenario wherein OACs covary. Nonetheless, a_CDOM(440) always produces higher R² values than C_a, with differences being insignificant for RSE2007 but increasing for RSE2022 and RSE2021.

The results of the algorithmic comparisons (Fig. 2) are consistent with those of the eigenanalyses (Fig. 1) in that both suggest that the OC algorithm more strongly associates with variability in a_CDOM(440) rather than C_a. The OC algorithm's applicability, therefore, depends strongly upon the relationships between OACs. The finding regarding the applicability of the OC algorithm is not the primary focus herein, and is consistent with previous studies demonstrating vulnerabilities in band-ratio retrievals of C_a to changing bio-optical relationships (Dierssen, 2010). The comparisons further suggest estimation of a_CDOM(440) is more readily obtained than estimation of C_a using remote sensing. For example, the OC algorithm produces a_CDOM(440) estimation with an R² of 0.85 for the RSE2021 dataset wherein optically complex waters are included, although applying the OC algorithm for a_CDOM(440) estimation is not recommended. Improved R² values greater than 0.90 have been shown for remote sensing of a_CDOM(440) using end-member analysis with VIS or INV wavebands (Houskeeper et al., 2021).

4.4 Spectral information content of aquatic spectra

Maximal association of the leading eigenfunction for each dataset with a_CDOM(440) rather than C_a indicates that accurate retrieval of a_CDOM(440) is fundamental to the objectives of the ocean color community to quantify oceanic productivity. Information extraction for subsequent eigenfunctions – including the second eigenfunction observed herein for most datasets to correlate most strongly with C_a – demands increasing radiometric capabilities (e.g., higher accuracy, more expansive spectral range, and perhaps finer spectral resolution), but also accurate retrieval of the information associated with the preceding eigenfunctions. Regarding the latter, small uncertainties or residual errors in retrievals of information corresponding to a leading eigenfunction produce so-called noise amplification in retrieval of information corresponding to subsequent eigenfunctions, just as the mean spread of S_i functions decreases with increasing i. Accurate retrieval of a_CDOM(440) is, therefore, a prerequisite to – not separate from – accurate retrieval of phytoplankton properties, in keeping with the trajectory of research into simultaneous inherent optical property inversions (e.g., Maritorena et al., 2002; Lee et al., 2002; Werdell et al., 2013), machine learning approaches (O'Shea et al., 2021), and degeneracy (Prochaska and Frouin, 2025).

Fine-scale spectral dependencies have previously been shown to confer useful and exploitable information – e.g., fluorescence (Letelier and Abbott, 1996), absorption band effects (Houskeeper et al., 2020), and macroalgal physiological state (Bell and Siegel, 2022) – which in some cases may not emerge with strong explanatory power in an eigenanalysis using the full spectral range, because Λ is a metric of variance captured across the full spectral range of the analysis. Eigenfunction analyses such as those performed here – as well as in Sathyendranath et al. (1987), Cael et al. (2023), and others – are useful tools for exploring the major patterns of a dataset. However, caution is also warranted because information only impacting a subset of the wavebands in an eigenanalysis may average out or suffer spectral dilution (Houskeeper et al., 2020).

The eigenanalyses recreated for the RSE2021 dataset but using VIS or INV spectral subsets (denoted VIS21 and INV21, respectively) also support previous findings demonstrating improvements in a_CDOM(440) algorithmic relationships with increasing spectral separation of the algorithm wavebands from the VIS domain (Hooker et al., 2020, 2021b; Houskeeper et al., 2020, 2021). For example, the results indicate INV information is more readily represented by a single mode of variability than VIS information, as follows: (a) the leading eigenfunction of VIS21 captures similar variance to that of the spectrally expansive UVN21; (b) the leading eigenfunction of INV21 captures approximately 17 % more variance than that of either UVN21 or VIS21; and (c) the second eigenfunction of VIS21 captures approximately 15 % more variance than that of INV21. Pearson's correlation coefficients relating S₁ to a_CDOM(440) increase strongly for the INV21 subset, approaching that of a perfectly coherent relationship (−0.96). The results indicate INV21 captures information more associated with a_CDOM(440) and less associated with C_a. Conversely, the leading eigenfunction of the VIS21 subset captures information more similarly associated with both a_CDOM(440) and C_a. The findings confirm previous literature based on theoretical formulations and observations, wherein spectrally expansive observations mitigate the confounding effects of spectral overlap in C_a and CDOM absorption (Houskeeper et al., 2021).

The recent launch of the PACE satellite mission supports novel spectrally expansive and hyperspectral observations of surface ocean environs. The spectral subset eigenanalyses indicate opportunities to leverage spectrally expansive PACE observations to advance OAC retrievals. For example, accurate retrieval of INV reflectances support retrieval of a_CDOM(440) with improved independence from C_a, and knowledge of a_CDOM(440) could then inform C_a retrievals using VIS information. Alternately, the results indicate that simultaneous retrieval of OACs, e.g., using inversions (Garver and Siegel, 1997) or MDNs (O'Shea et al., 2021), may improve if (quality assured) input reflectances are spectrally expansive and include INV domains. However, spectrally expansive in situ datasets are presently sparse and corresponding improvements in atmospheric correction are likewise required to support spectrally expansive satellite ocean observing.

5 Conclusions

Spectral AOP data products of aquatic environments produce leading eigenfunctions more strongly correlated to variability in CDOM absorption than C_a, consistent with advancing knowledge of microbial loop dynamics and an increasing diversity of trophic pathways represented therein (Azam, 1998). The similarity of results using three independent datasets strengthens support for CDOM primacy and suggests that accurate and independent estimation of CDOM is requisite to the retrieval of C_a and other phytoplankton parameters. Acknowledged limitations include the eigenanalyses capturing linear – but not non-linear – patterns of variability, although the relationships between the leading eigenfunction and log-transformed a_CDOM(440) indicated linearity. The eigenanalyses results are insensitive to moderate changes in the quality control filtering, but not to complete removal of quality control filtering for the legacy data, which contain significant proportions of likely non-physical data products. Confounding signals identified for CDOM and C_a are consistent with early investigations into the drivers of color variability (Yentsch, 1960), as well as subsequent work assessing vulnerabilities in band-ratio algorithms for characterizing C_a (Dierssen, 2010; Sauer et al., 2012). Spectrally expansive data products have been shown to improve retrieval of a_CDOM(440) independent of C_a (Sathyendranath et al., 1987; Hooker et al., 2020; Houskeeper et al., 2021), attributed, in part, to improved robustness to variability in the relationships between OACs. The comparisons herein of eigenanalyses using the INV21 and VIS21 spectral subsets further support the hypothesis that an expansive spectral range for data products improves the separability of signals. Spectrally expansive data products have been demonstrated in situ spanning the UV to shortwave infrared wavelength domain, with the latter formerly ascribed as null (Houskeeper and Hooker, 2023). The findings herein support opportunities associated with leveraging the spectrally expansive capabilities of PACE OCI, although further advances in image processing, atmospheric correction, algorithm design, and in situ datasets would be necessary.