Retrieving marine inherent optical properties from satellites using temperature and salinity-dependent backscattering by seawater.

Time-series of marine inherent optical properties (IOPs) from ocean color satellite instruments provide valuable data records for studying long-term time changes in ocean ecosystems. Semi-analytical algorithms (SAAs) provide a common method for estimating IOPs from radiometric measurements of the marine light field. Most SAAs assign constant spectral values for seawater absorption and backscattering, assume spectral shape functions of the remaining constituent absorption and scattering components (e.g., phytoplankton, non-algal particles, and colored dissolved organic matter), and retrieve the magnitudes of each remaining constituent required to match the spectral distribution of measured radiances. Here, we explore the use of temperature- and salinity-dependent values for seawater backscattering in lieu of the constant spectrum currently employed by most SAAs. Our results suggest that use of temperature- and salinity-dependent seawater spectra elevate the SAA-derived particle backscattering, reduce the non-algal particles plus colored dissolved organic matter absorption, and leave the derived absorption by phytoplankton unchanged.


Introduction
Satellite ocean color instruments provide consistent and high-density data records at sufficient temporal and spatial scales to allow retrospective analysis of long-term oceanographic trends. For example, the daily, synoptic images captured by the NASA Moderate Resolution Imaging Spectroradiometer onboard Aqua (MODISA) provide viable data records for observing decadal changes in biogeochemistry of both global and regional ecosystems [1]. Satellite ocean color instruments measure the spectral radiance emanating from the top of the atmosphere at discrete visible and infrared wavelengths. Atmospheric correction algorithms are applied to remove the contribution of the atmosphere from the total signal and produce estimates of remote sensing reflectances (R rs (Ȝ); sr −1 ), the light exiting the water normalized to a hypothetical condition of an overhead Sun and no atmosphere [2]. Bio-optical algorithms are applied to the R rs (Ȝ) to produce estimates of additional geophysical properties, such as spectral marine inherent optical properties (IOPs), namely the absorption and scattering properties of seawater and its particulate and dissolved constituents [3,4]. Time-series of these geophysical properties provide unparalleled resources for studying carbon stocks, phytoplankton population diversity and succession, and ecosystem responses to climatic disturbances on regional to global scales (e.g., [5][6][7][8]).
Semi-analytical algorithms (SAAs) provide one mechanism for estimating marine IOPs from R rs (Ȝ), through a combination of empiricism and radiative transfer theory. Many published SAAs attempt to simultaneously estimate the magnitudes of spectral backscattering by particles (b bp (Ȝ); m −1 ), absorption by phytoplankton (a ph (Ȝ); m −1 ), and the combined absorption by non-algal particles and colored dissolved organic material (a dg (Ȝ); m −1 ) [3,4 and references therein]. These SAAs generally assign constant spectral values for seawater absorption and backscattering (a w (Ȝ) and b bw (Ȝ), respectively; m −1 ) [9,10], assume spectral shape functions of the remaining constituent absorption and scattering components, and retrieve the magnitudes of each remaining constituent required to match the spectral distribution of R rs (Ȝ). As most SAAs differ only in the assumptions employed to define component spectral shapes and mathematical methods applied to calculate the magnitudes, the community typically focuses solely on refining these varied assumptions and methods [3].
While it is well known that a w (Ȝ) and b bw (Ȝ) in the visible spectrum vary with ocean temperature (T; °C) and salinity (S; g kg −1 ) [11][12][13][14][15][16], the impact of using constant values in lieu of temperature-and salinity-dependent values in an SAA has yet to be fully explored or quantified. Analyses conducted by Morel [11] revealed scattering increases of 30% in the spectral range 366-578 nm for S = 38.4 g kg −1 compared to those for pure water. This dependency of b bw (Ȝ) on S was most recently revisited by Zhang [16], whose theoretical model agrees to within 1% on average with Morel [12] at S = 38.4 g kg −1 . Despite this evolving understanding of how the spectral shape and magnitude of b bw (Ȝ) depend on T and S, the ocean color community continues to exclusively adopt the invariant Morel b bw (Ȝ) values [9,11,12] for use in SAAs and other bio-optical modeling activities [3,4], even though these values apply strictly to a singular temperature and salinity combination. The Zhang model suggests that b bw (Ȝ,T,S) can deviate from b bw (Ȝ,20,38.4) by up to 25% at 443 and 547 nm, respectively, for typical temperature and salinity ranges found in natural waters (−2 ≤ T ≤ 40°C and 0 ≤ S ≤ 40 g kg −1 ) (Fig. 1). Given that b bw (Ȝ) ≥ b bp (Ȝ) for some wavelengths in many offshore and open ocean waters [15], one might expect that including temperature and salinity-dependent b bw (Ȝ) in an SAA (in place of a constant spectrum) will spatially and temporally alter the magnitudes of its derived products.
We executed a series of analytical experiments to evaluate this expectation. In this paper, we present a preliminary quantification of the changes realized in SAA-derived products when using temperature-and salinity-dependent b bw (Ȝ) in lieu of a constant spectrum. We generated this quantification by applying the default SAA described in Werdell [3] (GIOP) to a wide dynamic range of in situ data and a global MODISA time-series using both b bw (Ȝ) from Morel [12] (which is the default for GIOP) and from Zhang [16]. For the latter, we calculated b bw (Ȝ,T,S) using the NOAA World Ocean Atlas sea surface salinity climatologies [17] and the NOAA optimally interpolated sea surface temperature climatologies [18] as inputs. While our use of climatologies may not have perfectly captured the potential differences realized using real-time data, they effectively enabled these proof-of-concept analyses and indicated that significant relative biases exist between b bp (Ȝ), a dg (Ȝ), and a ph (Ȝ) derived from the SAA runs with and without temperature-and salinity-dependent b bw (Ȝ). We conclude with a brief discussion regarding the use of temperature-and salinity-dependent a w (Ȝ) in an SAA, however, we do not quantify its impact.

A general semi-analytical algorithm
Ocean color satellite instruments provide estimates of R rs (Ȝ), which can be converted to their subsurface values using the method presented in Lee [19]: Subsurface remote-sensing reflectances relate to marine IOPs via: where b b is the total backscattering coefficient (m −1 ), a is the total absorption coefficient (m −1 ), and g (sr −1 ) vary with illumination conditions, sea surface properties, and the shape of the marine volume scattering function. GIOP defaults to g 1 = 0.0949 and g 2 = 0.0794 sr −1 from Gordon [20]. The absorption coefficient can be expanded as the sum of all absorbing components. Further, each component can be expressed as the product of its concentrationspecific absorption spectrum (eigenvector; a*) and its magnitude (eigenvalue, M): GIOP expresses a* dg (Ȝ) as exp(-S dg Ȝ), where S dg describes a rate of exponential decay (nm −1 ), and assigns S dg = 0.018 nm −1 . The term M ph is equivalent to the concentration of the phytoplankton pigment chlorophyll-a (C a ; mg m −3 ) when a* ph (Ȝ) is expressed as a chlorophyll-a specific absorption spectrum (m 2 mg −1 ). GIOP uses the chlorophyll-specific a* ph (Ȝ) from Bricaud [21], seeded using an estimate of C a from O'Reilly [22] and normalized to a* ph (443) = 0.055 m 2 mg −1 . Similar to a(Ȝ), b b (Ȝ) can be expanded to: GIOP expresses b* bp (Ȝ) as Ȝ Sbp , where S bp defines the steepness of the power law (we acknowledge the validity of this power function remains debatable), and estimates S bp using the R rs (Ȝ) band-ratio approach described in Lee [19]. Werdell [3] describes in detail the rationale for the selection of the eigenvectors used in Eqs. (3) and (4).
Using R rs (Ȝ) and eigenvectors as input, the eigenvalues M can be estimated via linear or nonlinear least squares inversion of Eqs. (1) to (4). GIOP employs the nonlinear Levenberg-Marquardt method and uses all R rs (Ȝ) between 400 and 700 nm. We considered retrievals of b bp (Ȝ), a dg (Ȝ), and a ph (Ȝ) to be valid when they fell within predefined ranges (e.g., −0.05 a w (Ȝ) ≤ a dg (Ȝ) ≤ 5 m −1 ) and could reconstruct R rs (Ȝ) that differed from the input spectrum by less than 33% anywhere in the range 400-600 nm. Werdell [3] describes the additional algorithm metrics and quality control practices embedded within GIOP.

Data acquisition
We acquired coincident observations of in situ R rs (Ȝ), b b (Ȝ), a dg (Ȝ), and a ph (Ȝ) from the NASA bio-Optical Marine Algorithm Data set (NOMAD) [23]. Sample sizes for b b (Ȝ) (N = 217) differ from a dg (Ȝ) and a ph (Ȝ) (N = 609) in NOMAD because of instrumental variability in field campaigns included in its compilation. We obtained MODISA Level-3 R rs (Ȝ) monthly composites for 2010 from the NASA Ocean Biology Processing Group (OBPG; http://oceancolor.gsfc.nasa.gov). We also acquired NOAA World Ocean Atlas 2009 sea surface salinity monthly climatologies (WOA-SSS) [17] and NOAA optimally interpolated sea surface temperature weekly climatologies (OISST) [18] from the OBPG, as they are available through the NASA SeaWiFS Data Analysis System (SeaDAS; http://seadas.gsfc.nasa.gov). The NOAA National Oceanographic Data Center and NOAA Earth System Research Laboratory provide and maintain the official source WOA-SSS and OISST data records, respectively.

Backscattering of seawater
As suggested above, most SAAs employ the b bw (Ȝ) provided by Morel [12], as tabulated via a best-fit expression by Smith and Baker [9]. This fit is commonly expressed as a power-law, namely b bw (Ȝ) = 0.0038 (400 / Ȝ) -4.32 , where b bw (400) = 0.0038 m −1 (e.g., [24]). This expression provides the constant b bw (Ȝ) used in this activity, which we hereafter refer to as b bw c (Ȝ). Using the WOA-SSS and OISST data records as input into Zhang [16] (via MathWorks MATLAB software kindly provided by X. Zhang), we calculated temperature-and salinity-dependent b bw (Ȝ) for each station in NOMAD and pixel in the MODISA monthly composites. We hereafter refer to temperature-and salinity-dependent b bw (Ȝ) as b bw TS (Ȝ).

Data processing and analysis
We applied GIOP to each station in NOMAD and each pixel in the MODISA monthly composites twice -first using b bw c (Ȝ) and second using b bw TS (Ȝ), with all other parameters held constant. We hereafter refer to these two instances of GIOP as GIOP-C and GIOP-TS, respectively. We focused on two biogeochemically distinct regions in the MODISA monthly composites -the North Atlantic (40 to 50°N, −50 to −20°W) and South Pacific Gyre (−20 to −30°S, 100 to 130°E). We used OBPG-supported satellite data processing software (l3gen; distributed with SeaDAS) to run GIOP on the Level-3 MODISA R rs (Ȝ). Currently, l3gen accepts WOA-SSS and OISST data records as standard ancillary input. For NOMAD, we compared modeled and in situ b bp (Ȝ), a dg (Ȝ), and a ph (Ȝ) for the two GIOP runs using ordinary least squares regression. Our statistics of interest included the coefficient of determination (r 2 ), the regression slope, the bias, and the absolute median percent difference (MPD; see Table 1 for its calculation). We executed these comparisons using

Results
Direct comparisons of modeled and in situ b bp (Ȝ), a dg (Ȝ), and a ph (Ȝ) provided estimates of the accuracy of GIOP-C and GIOP-TS. Visually, the model-versus-in situ regression analyses for the two runs showed only subtle differences (Fig. 2). For MODISA wavelengths, GIOP-TS realized r 2 that were unchanged for b bp (Ȝ), slightly degraded for a dg (Ȝ) (~-3%), and modestly improved for a ph (Ȝ) (~ + 10%) relative to GIOP-C (Table 1). GIOP-C showed spectral dependence in its regression slopes for b bp (Ȝ) (rising from 0.88 to 1.00 for 412 to 667 nm), but GIOP-TS did not, hovering around unity (0.97-0.98) for all wavelengths. GIOP-C and GIOP-TS yielded similar spectral dependence in the regression slopes for a dg (Ȝ) (falling from 1.11 to 0.90 and 1.13 to 0.92 for 412 to 667 nm, respectively). But, the regression slopes for a ph (Ȝ) fell closer to unity for GIOP-TS (average of 1.07) than for GIOP-C (average of 1.11). Finally, GIOP-TS realized improved (reduced) biases and MPDs for b bp (Ȝ) relative to GIOP-C, with the exception of b bp (667). GIOP-TS yielded slightly degraded biases and MPDs for a dg (Ȝ), but comparable statistics for a ph (Ȝ), relative to GIOP-C. On average, estimates of b bp (Ȝ) from GIOP-TS rose 2.7% compared to those from GIOP-C for the full population of data in NOMAD (Fig. 3). In contrast, estimates of a dg (Ȝ) and a ph (Ȝ) fell 6.1 and 1.0%, respectively. These departures showed dependence on the trophic state of the station (e.g., low C a , oligotrophic water versus high C a , eutrophic water). In particular, GIOP-TS yielded much higher b bp (Ȝ) than GIOP-C (3-10%) for the clearest NOMAD stations (C a < 1 mg m −3 ), whereas a dg (Ȝ) from GIOP-TS fell by 6% in this range (Fig. 4). Both GIOP-TS and GIOP-C produced comparable a ph (Ȝ) over the full dynamic range of NOMAD C a , in particular when C a > 0.2 mg m −3 . GIOP-TS yielded slightly lower a ph (443) than GIOP-C in the clearest waters, for example, where C a ~0.1 mg m −3 .
MODISA global imagery from May 2010 showed geographic features that mimicked this "trophic-level" behavior (Fig. 5). In the oligotrophic sub-tropical gyres, for example, b bp (443) from GIOP-TS exceeded that from GIOP-C by 10-20%, whereas a dg (443) and a ph (443) fell 5-10% and 1-5%, respectively. By far, b bp (443) realized the largest spatial differences in magnitude, as most open ocean waters showed significantly elevated GIOP-TS to GIOP-C ratios. The Southern Ocean (circumpolar Antarctica) and near-shore waters presented the only exceptions, where b bp (443) was largely equivalent for both GIOP runs. The modest differences in a ph (443) followed a similar and equally broad spatial pattern to b bp (443), differing by a few percent in the open ocean and approaching unity in the Southern Ocean and along the coasts. In contrast, a dg (443) from GIOP-TS and GIOP-C only truly differed in the centers of the sub-tropical gyres (and occasionally along the coast). While the Southern Ocean and coasts often realize the largest departures in T and S from the 20°C and 38.4 g kg −1 inherent to b bw c (Ȝ) (Fig. 6), the differences between b bw TS (Ȝ) and b bw c (Ȝ) remain minimal for the T and S combinations of these regions (Fig. 1).
Time-series of MODISA monthly averages in the North Atlantic showed seasonal patterns in the ratios of GIOP-TS to GIOP-C b bp (443) (Fig. 7). GIOP-TS produced b bp (443) that exceeded that from GIOP-C by 3-5% in April-July and by 7-12% in January-March and August-December. Ratios of GIOP-TS to GIOP-C a ph (443) followed a subtle, inverse seasonal pattern, although the two retrievals never departed by more than 3% (in boreal winter). The b bp (443) and a ph (443) time-series tracked each other through a natural, seasonal progression of phytoplankton biomass [25], with higher b bp (443) generally occurring in the presence of higher a ph (443) (e.g., during and following the spring bloom). The highest GIOP-TS to GIOP-C ratios for b bp (443) appeared with the lowest b bp (443), while the highest ratios for a ph (443) appeared with the highest a ph (443) in the late boreal spring. Interestingly, despite a spike in magnitude in boreal spring and a dip in boreal winter, the ratios of a dg (443) showed no temporal dependence, but rather remained fixed with GIOP-C exceeding GIOP-TS by 6%.
MODISA monthly averages in the South Pacific Gyre showed little seasonal pattern in the ratios of GIOP-TS to GIOP-C (Fig. 7). This region encompasses the clearest known ocean water [15] and the annual patterns in IOP magnitudes followed expectations, with minor variability in b bp (443), a peak in a ph (443) in June and July (austral winter), and a subsequent rise in a dg (443) [26,27]. Overall, GIOP-TS produced b bp (443) that exceeded GIOP-C by 11%, a dg (443) that fell below by 6%, and a ph (443) that fell below by 4% consistently over the full annual cycle. In general, the average magnitudes of these departures follow those observed for the clearest stations in NOMAD (e.g., C a < 0.1 mg m −3 ; Fig. 4).

Discussion
We demonstrated that using b bw TS (Ȝ) (e.g., [16]) in an SAA in lieu of the community standard b bw c (Ȝ) [12] impacts the derived products in a sufficiently significant way as to merit further consideration. Theoretically and pragmatically, our results suggest that an SAA configuration such as GIOP-TS provides a step towards more accurate estimation of marine IOPs. Quantitatively, we expect our results could be refined with more contemporaneous values of T and S -for example, using coincident sea surface temperature data records from MODISA and complimentary sea surface salinity records from the Aquarius instrument onboard the Argentinean SAC-D spacecraft. Our case study cumulatively suggests that using b bw TS (Ȝ) results in elevated b bp (Ȝ) (3-10%), reduced a dg (Ȝ) (1-6%), and negligibly changed a ph (Ȝ).
Given that b bw TS (Ȝ) typically falls below b bw c (Ȝ) for most S<38.4 g kg −1 (Fig. 1), we expected GIOP-TS estimates of b bp (Ȝ) to exceed those from GIOP-C (Figs. 3-6). Also as expected, the largest departures in b bp (Ȝ) occurred in the open ocean where b bw (Ȝ) ≥ b bp (Ȝ) [15] (in contrast to highly turbid or productive environments where b bw (Ȝ) << b bp (Ȝ)). Of interest, however, which we elaborate upon later, is that b b (Ȝ) from GIOP-TS showed a consistently low bias relative to that from GIOP-C (Fig. 8). The NOMAD regression results indicate that b bp (Ȝ) from GIOP-TS realized improved agreement with in situ values (Fig. 2, Table 1). For the ocean color remote sensing paradigm, this implies that SAAs currently underestimate b bp (Ȝ). Following, the switch to b bp TS (Ȝ) could refine our future scientific understanding of marine particle dynamics and primary productivity (e.g., [7]). More compellingly, the MODISA time-series indicate that using b bw TS (Ȝ) in place of b bw c (Ȝ) imparts an amplified seasonal signal in b bp (Ȝ), which could further alter our understanding of and expectations for these phenomena. Finally, the spectral slope for b bw TS (Ȝ) also varies with salinity (Zhang [16] reported a range of −4.286 to −4.306 from 0 to 40 g kg −1 ), whereas the slope of b bw c (Ȝ) does not, which will impact SAAs such as Loisel [5] that produce estimates of the spectral slope of b bp (Ȝ), S bp , to infer marine particle sizes. Using b bw TS (Ȝ) yielded reduced estimates of both a dg (Ȝ) and a ph (Ȝ), although departures in the latter were largely negligible (Figs. [3][4][5][6]. The NOMAD regression results indicate that a ph (Ȝ) from GIOP-TS improves upon that from GIOP-C, but that a dg (Ȝ) degrades slightly (Fig.  2, Table 1). Werdell [3] demonstrated that SAA retrievals are highly sensitive to the parameterization of the spectral slope of a dg (Ȝ), S dg , and that changes in a dg (Ȝ) often accompany changes in b bp (Ȝ) (a classic case of transfer of variance, as a dg (Ȝ) and b bp (Ȝ) have similar spectral shapes). However, repeating our experiment with S dg = 0.015 and = 0.02 nm −1 , and with S dg dynamically estimated using the R rs (Ȝ) band-ratio approach described in Lee [19], produced equivalent results (not shown). Globally, the most significant departures occurred in subtropical gyres (Fig. 5). While this maintains some significance for studies relating to carbon budgets (particularly, the optically-relevant dissolved component, e.g., [6,8]), the bias between a dg (Ȝ) from GIOP-TS and GIOP-C remained fairly constant over most trophic levels (Fig. 4) and seasons (Fig. 7). Following, the switch to b bw TS (Ȝ) may produce a dg (Ȝ) with reduced magnitudes, but not with significantly altered spatial distributions. We recommend that future studies continue to focus on methods to dynamically assign S dg to each station in lieu of using spatially and temporally fixed values, with the purpose of reducing differences in SAA-derived and in situ a dg (Ȝ).
Intuitively, the elevated estimates of b bp (Ȝ) realized with b bw TS (Ȝ) appear in contradiction with the reduced estimates of a dg (Ȝ) and a ph (Ȝ). Given fixed a w (Ȝ) [10], these reduced estimates yielded a(Ȝ) from GIOP-TS that fell below those from GIOP-C by 3.2% on average (Fig. 8). The contradiction lies within Eq. (2) -that is, we expect common u(Ȝ) from both GIOP runs, assuming equivalent least-squares fitting performances, which suggests that b b (Ȝ) from GIOP-TS should be reduced in conjunction with a(Ȝ). With b bw TS (Ȝ) typically less than b bw c (Ȝ), and b bp (Ȝ) from GIOP-TS holistically greater than that from GIOP-C, we initially anticipated b b (Ȝ) from GIOP-TS to be roughly equivalent to that from GIOP-C. In other words, we expected a corresponding rise in magnitude in derived b bp (Ȝ) that accounted for the decline in magnitude in b bw TS (Ȝ). In practice, however, estimates of b b (Ȝ) from GIOP-TS fell 3.5% below those from GIOP-C (Fig. 8). This decline in b b (Ȝ) ultimately yielded equivalent u(Ȝ) for the two runs (0.15% different on average). In the end, using b bw TS (Ȝ) in lieu of b bw c (Ȝ) yielded reduced b b (Ȝ) and a(Ȝ), despite producing elevated b bp (Ȝ).
Changes in the parameterization of an SAA (e.g., the component spectral shapes) change the derived products to varying degrees [3]. Our results do not imply that GIOP-TS provides the optimal SAA configuration for all water masses at all times [28]. As discussed in Werdell [3], significant work remains to refine SAAs for both global and regional applications. Westberry [29], for example, recently demonstrated that accounting for Raman effects can additionally improve SAA-estimates of b bp (Ȝ). Despite the work that remains to be accomplished, however, we felt it prudent to repeat several of our experiments with an alternative form of SAA, namely the Quasi-Analytical Algorithm (QAA [19];). A description of the differences in algorithm form relative to GIOP exceed the scope of this paper, but can be found in Werdell [3]. For the NOMAD R rs (Ȝ), using b bw TS (Ȝ) in lieu of b bw c (Ȝ) as input into QAA yielded elevated b bp (Ȝ) (3.2%), reduced a dg (Ȝ) (5.5%), equivalent a ph (Ȝ) (0.5%), reduced b b (Ȝ) (2.8%), and reduced a(Ȝ) (2.9%) on average. These departures effectively match those from GIOP in both magnitude and direction. Cumulatively, with the assumption that global climatologies of T and S provided sufficiently reliable values for this case study, our results imply that using temperature and salinity-dependent b bw (Ȝ) [16] in place of the (most commonly adopted) constant spectrum [12] in an SAA results in elevated b bp (Ȝ) and reduced a dg (Ȝ). Spatially and temporally, a ph (Ȝ) from GIOP-TS and GIOP-C showed little variability, with the exception of the highest latitudes. We acknowledge that the departures in b bp (Ȝ) and a dg (Ȝ) may be sufficiently small to also be considered insignificant, however, it appears that using b bw TS (Ȝ) produces improved comparisons with field data, at least for b bp (Ȝ) and a ph (Ȝ). Using b bw TS (Ȝ) in an SAA also appears to change the spatial and temporal distributions of b bp (Ȝ), which should be sufficiently compelling itself to stimulate subsequent investigations, despite the degraded comparisons between modeled and in situ a dg (Ȝ). We limited our activity to the use of b bw TS (Ȝ) in an SAA. In practice, b bw (Ȝ) appears elsewhere within the satellite ocean color paradigm, such as within the atmospheric correction process [2]. For example, the standard NASA algorithm applied to account for non-negligible near-infrared (NIR) radiances uses b bw (Ȝ) and a w (Ȝ) to model R rs (NIR) [30]. Although this specific algorithm employs an iterative process, using b bw TS (Ȝ) in lieu of b bw c (Ȝ) has potential to alter R rs (Nǿȇ) and, thus, the aerosol models selected with the atmospheric correction process. Likewise, a w (Ȝ) appears throughout the ocean color paradigm. Unlike b bw (Ȝ), the temperature-and salinity-dependence of a w (Ȝ) most strongly effects the NIR portion of the spectrum [13,14]. We, therefore, do not anticipate that using a w TS (Ȝ) in place of a w c (Ȝ) will significantly alter SAA-derived IOPs. With regards to the atmospheric correction component presented above, however, this replacement will almost certainly alter the modeled R rs (Nǿȇ) derived following the method Bailey [30]. Accounting for non-negligible R rs (Nǿȇ) remains a critical step within the atmospheric correction process in turbid waters ( [31] and references therein). We recommend that similar, subsequent studies emerge to evaluate the use of temperature-and salinity-dependent a w (Ȝ) in standard satellite ocean color data processing.