Determining the Primary Sources of Uncertainty in Retrieval of Marine Remote Sensing Reflectance From Satellite Ocean Color Sensors

Uncertainties in the retrieval of the remote sensing reflectance, R_rs, from Ocean Color (OC) satellite sensors have a strong impact on the performance of algorithms for the estimation of chlorophyll-a, mineral concentrations, and inherent optical properties (IOPs). The uncertainties are highest in the blue bands. The total radiance measured at the top of the atmosphere captures the instantaneous state of the atmosphere-ocean system: the in-water conditions, sky and Sun glint reflected from the wind-roughened ocean surface, as well as light scattered from molecules and aerosols in the atmosphere. Each of these components has associated uncertainties, and when combined with the additional uncertainties from the instrument noise and the atmospheric correction process, they contribute to the total uncertainty budget for the retrieved R_rs. We analyzed the contribution of each component uncertainties to the total R_rs uncertainties in SNPP-VIIRS level 2 products, taking advantage of the spectral differences between the components. We examined multiple scenes in the open ocean and coastal waters at spatial resolutions ranging from 2250 to 5250 m by comparing the retrieved R_rs to in situ measurements made at several AERONET-OC sites and at the MOBY site. It was shown that uncertainties associated with the molecular (Rayleigh) scattering play the most significant role, while the contributions of other components are usually smaller. Uncertainties in Rayleigh scattering are primarily attributed to the variability of Rayleigh optical thickness (ROT) with a standard deviation of approximately 1.5% of ROT, which can largely explain the frequency of negative R_rs retrievals as observed using the current standard atmospheric correction process employed by NASA. Variability of the sky light reflected from the ocean surface in some conditions also contributed to uncertainties in the blue; water variability proportional to R_rs had a very pronounced peak in the green at coastal sites.

Introduction

Ocean Color (OC) is indicative of ocean health and biochemistry, and for that reason is listed as an essential climate variable (ECV) (IOCCG, 2008). The color of a water body is determined by scattering and absorption of pure water and its natural constituents, such as phytoplankton, non-algal particles, and colored dissolved organic matter (CDOM) (Mobley, 1994). Some phytoplankton species form harmful algal blooms (HABs), which can be toxic and affect human and marine life, and more generally the health of the ecosystem, fishing industry, and recreation activities (IOCCG, 2021). Typically, less than 10 percent of the top of the atmosphere (TOA) radiance is due to the water signal at sea level (Gordon and Morel, 1983), with the remainder originating from scattering processes in the atmosphere and reflections of the Sun and sky on the wave-roughened water surface. It is paramount to accurately estimate radiances at the surface level from the ones at the TOA, as uncertainties propagate into the retrieval of water parameters, characteristics of in-water particulates, concentrations of chlorophyll-a, and detection of algal blooms (IOCCG, 2010; IOCCG, 2019).

Atmospheric correction uncertainties stem at least partially from the estimation of aerosol models, and air–water interface effects due to sky and Sun light reflections at the wind-roughened air–water interface (Gordon and Wang, 1992; Gordon and Wang, 1994; Frouin et al., 1996; Wang and Bailey, 2001; Ahmad et al., 2010; Frouin et al., 2019). Atmospheric correction uncertainties often have a stronger impact on retrievals in coastal waters with low water reflectance values in the blue bands (Carrizo et al., 2019; Groetsch et al., 2020; Wei et al., 2020). Currently, high uncertainties in the blue reflectance observations are widely acknowledged (IOCCG, 2019; Wei et al., 2020; Li et al., 2019; Herrera-Estrella et al., 2020). However, little is known about the specific dependencies of uncertainties concerning their spectral and scaling behavior in various water areas, and their dependence on meteorological conditions.

The spectra of aerosol radiances in coastal areas are significantly affected by the presence of absorbing aerosols (Gordon et al., 1997; Ransibrahmanakul and Stumpf, 2006; Shi and Wang, 2007). Atmospheric correction processing schemes for the current satellite sensors do not account for this effect due to a lack of information about aerosol parameters, and it is assumed that this leads to negative values of R_rs in blue bands (Frouin et al., 2019). This makes the estimation of chlorophyll concentration and water parameters inaccurate in such waters. Several partial solutions to this problem are found in the existing literature and include more complex processing in the atmospheric correction (Gordon et al., 1997; Oo et al., 2008), removal of the uncertainty as a power law-like “artifact” with exponent -6 (Ransibrahmanakul and Stumpf, 2006), neural network approaches (Fan et al., 2021), utilization of atmospheric correction algorithms based on the fitting of Rayleigh spectra (Steinmetz et al., 2011; Zhang et al., 2019), and simply avoiding blue bands in algorithms for the retrieval of water parameters (El-Habashi et al., 2019; Gilerson et al., 2021). The upcoming NASA PACE mission (Werdell et al., 2019) will have the hyperspectral Ocean Color Instrument (OCI) and two polarimeters on board, which are expected to provide broader information on aerosol parameters. With these, atmospheric correction processing is expected to be significantly improved. However, for the current sensors, it is important to have a more accurate understanding of the impact of absorbing aerosols on R_rs retrievals.

There are differences in the terminology regarding “uncertainties.” In the Guide to Uncertainty in Measurement (Sayer et al., 2020), uncertainty is defined as an expression of the dispersion of the measurand (in our case, R_rs retrieved from satellite observations), and it is often represented as one standard deviation around the retrieved value. In our case, the standard deviation would represent R_rs uncertainties due to spatial variability (Herrera-Estrella et al., 2021), where the difference between the mean R_rs value and “true” R_rs value is up to 4–5 times greater than one standard deviation. In this work, uncertainties of R_rs will be described as the root mean square difference (RMSD) between the mean R_rs value determined from satellite data and in situ R_rs, which is considered as a “true” value, similar to the approach in IOCCG (2019).

The estimation of the uncertainties can be carried out by the comparison of the parameters determined from the satellite imagery with the “true” values. These comparisons can be made in clear waters, where all the water parameters can be connected to the concentration of chlorophyll-a, [Chl] (Hu et al., 2013). Another approach is to compare data from satellite sensors with field measurements from offshore platforms, autonomous systems such as the AERONET-OC network (Zibordi et al., 2009, 21), buoys like Marine Optical BuoY (MOBY) (Clark et al., 1997), and from ships (Moore et al., 2015). Specifically, the (Moore et al., 2015) uncertainties in R_rs were estimated for seven optical water types (OWT) using data from the SeaBASS optical database (Werdell et al., 2003) for the MOBY site, the BOUSSOLE mooring site in the Mediterranean Sea (Antoine et al., 2008), and AERONET-OC sites within 5 × 5 pixel boxes, plus/minus 3 h differences. It was found that R_rs uncertainties are generally the highest in the blue part of the spectrum in both clear and coastal waters.

In yet another approach, using Monte Carlo (MC) simulations for Sea-viewing Wide Field-of-view Sensor (SeaWiFS) observations (Franz et al., 2016), the retrieval process for R_rs was repeated 1000 times, and uncertainties in R_rs were then estimated as the “standard deviation of the 1000 perturbed R_rs retrievals in each band.” This derived uncertainty was interpreted “as the precision of the R_rs retrieval due to instrument noise.” It was about 4 times smaller than the observed R_rs uncertainties based on in situ validation (Moore et al., 2015).

Furthermore, in the study by Herrera-Estrella et al. (2021), a model was developed to evaluate the spectral composition of R_rs uncertainties, which was applied to characterize uncertainties due to the R_rs spatial distribution in images from the Visible Infrared Imaging Radiometer Suite (VIIRS) sensor on the SNPP platform and the Landsat-8 Operational Land Imager (OLI) at different spatial resolutions. Most of these uncertainties were attributed to the surface effects and water variability conditions.

In this work, a similar model is applied to estimate the spectral components of the uncertainties in R_rs retrieval by comparing SNPP VIIRS satellite data and in situ data from the MOBY site and eight AERONET-OC stations in US and European waters.

In Theoretical Considerations in the Estimation of Uncertainties, the model (Herrera-Estrella et al., 2021) is replicated with several modifications, in VIIRS Satellite and AERONET-OC, data satellite and in situ data are described, and the results are presented in Results. Discussion and conclusions are provided in Discussion and Conclusion.

Theoretical Considerations in the Estimation of Uncertainties

Main Relationships

The main radiometric quantity in the processing of satellite data is the remote sensing reflectance, R_rs, which is defined as the ratio of the water-leaving radiance to the downwelling irradiance at the sea surface, $R_{rs}\left(\lambda \right)=L_w\left(\lambda \right)/E_d\left(\lambda \right)$, where $L_w\left(\lambda \right)$ is the water-leaving radiance, $E_d\left(\lambda \right)$ is the downwelling irradiance, and λ is the wavelength. At the top of the atmosphere (TOA), the total radiance, $L_t^0\left(\lambda \right)$, can be represented as (Gordon and Wang, 1994; Mobley, 2022)

${\mathit{L}}_{\mathit{t}}^0\left(\lambda \right)={\mathit{L}}_{\mathit{r}}^0\left(\lambda \right)+{\mathit{L}}_{\mathit{a}}^0\left(\lambda \right)+{\mathit{L}}_{\mathit{g}}^0\left(\lambda \right)+{\mathit{t}}^0\left(\lambda \right){\mathit{L}}_{\mathit{w}}^0\left(\lambda \right),(\mathrm{1a})$

where $L_r^0\left(\lambda \right)$ is the total Rayleigh radiance at the TOA, which includes Rayleigh scattering and surface effects, $L_a^0\left(\lambda \right)$ is the total aerosol radiance, $L_g^0\left(\lambda \right)$ is the direct sun glint radiance from the water surface at TOA, $L_w^0\left(\lambda \right)$ is the water-leaving radiance just above the surface, and $t^0\left(\lambda \right)$ is the diffuse transmittance of light from the water surface to the TOA in the viewing direction. Superscript “0” denotes the actual parameters in the water and atmosphere. $L_t\left(\lambda \right)$, measured at the satellite sensor at 3 × 3 or more pixels, has uncertainties due to all of these components and to vicarious calibration and sensor noise. In the process of retrieval of the water-leaving radiance $L_r\left(\lambda \right)$, $L_a\left(\lambda \right)$, $L_g\left(\lambda \right)$ , and $t\left(\lambda \right)$ are modeled, and radiances are subtracted from the measured radiance, $L_t\left(\lambda \right)$, which introduces another set of uncertainties between actual and modeled radiances and transmittance coefficients:

${\mathit{L}}_{\mathit{w}}\left(\lambda \right)=({\mathit{L}}_{\mathit{t}}\left(\lambda \right)-{\mathit{L}}_{\mathit{r}}\left(\lambda \right)-{\mathit{L}}_{\mathit{a}}\left(\lambda \right)-{\mathit{L}}_{\mathit{g}}\left(\lambda \right))/\mathit{t}\left(\lambda \right),(\mathrm{1b})$

where radiances $L_r\left(\lambda \right)$, $L_a\left(\lambda \right)$, $L_g\left(\lambda \right)$, and $L_w\left(\lambda \right)$ are modeled radiances.

In addition, $L_r^0\left(\lambda \right)$ and $L_r\left(\lambda \right)$ can be divided into the radiance from the Rayleigh scattering in the atmosphere and reflectance from the ocean surface:

${\mathit{L}}_{\mathit{r}}^0\left(\lambda \right)={\mathit{L}}_{\mathit{R}}^0\left(\lambda \right)+\mathit{t}\left(\lambda \right){\mathit{L}}_{\mathit{surf}}^0\left(\lambda \right)\text{,}(\mathrm{2a})$

$\mathit{ }{\mathit{L}}_{\mathit{r}}\left(\lambda \right)={\mathit{L}}_{\mathit{R}}\left(\lambda \right)+\mathit{t}\left(\lambda \right){\mathit{L}}_{\mathit{surf}}\left(\lambda \right)\text{,}(\mathrm{2b})$

where $L_{surf}\left(\lambda \right)=L_{sky}\left(\lambda \right)\ast \rho$, $L_{sky}\left(\lambda \right)$ is the sky radiance, and ρ is the reflectance coefficient from the water surface; this is similar for $L_{surf}^0\left(\lambda \right)$. In the satellite atmospheric correction procedure, averaged surface effects are considered in the vector radiative transfer (VRT) equations which are based on Cox–Munk distributions (Gordon and Wang, 1992; Cox and Munk, 1954; Hu and Carder, 2002) together with the Rayleigh scattering, and thus, $L_r\left(\lambda \right)$ in Eq. 1 is not separated into its components. In this work, one of the goals is to estimate separately uncertainties from the Rayleigh scattering and surface effects, and that is why both components are considered separately, and a very small term of Rayleigh–surface interactions is not considered. In the atmospheric correction procedure, surface effects are also included in a similar manner in the modeling of aerosol $L_a\left(\lambda \right)$ radiance (Gordon and Wang, 1992). Here, such effects are not considered separately as well. Regarding the VRT estimation of surface effects, each satellite image captures a specific snapshot of the ocean, where the actual spatial average of the light field reflected from the wave facets may not exactly match the average predicted by the VRT model. The actual signal may have its own features due to the instantaneous water and atmospheric conditions, spatial scales in the area, or due to simplifying assumptions made within the VRT model, as the Cox–Munk model is not necessarily valid for waters in coastal areas.

Then, from Eqs 1, 2 for $L_w\left(\lambda \right)$ assuming $t\left(\lambda \right)=t^0\left(\lambda \right)$,

${\mathit{L}}_{\mathit{w}}=({\mathit{L}}_{\mathit{t}}-{\mathit{L}}_{\mathit{t}}^0+{\mathit{L}}_{\mathit{R}}^0-{\mathit{L}}_{\mathit{R}}+{\mathit{L}}_{\mathit{a}}^0-{\mathit{L}}_{\mathit{a}}+{\mathit{L}}_{\mathit{g}}^0-{\mathit{L}}_{\mathit{g}})/\mathit{t}+{\mathit{L}}_{\mathit{surf}}^0-{\mathit{L}}_{\mathit{surf}}+{\mathit{L}}_{\mathit{w}}^0.(3)$

Uncertainties from all the components included in Eqs 1–3 in the recording of the signal and in the retrieval process need to be taken into account. Normalizing by the downwelling irradiance, $E_d\left(\lambda \right)$, the uncertainty in remote sensing reflectance σ in sr⁻¹ can be determined from

${\mathit{\sigma }}^2=({\mathit{\sigma }}_{\mathit{t}}^2+{\mathit{\sigma }}_{\mathit{R}}^2+{\mathit{\sigma }}_{\mathit{a}}^2+{\mathit{\sigma }}_{\mathit{g}}^2)/\mathit{ }{\mathit{t}}^2+{\mathit{\sigma }}_{\mathit{surf}}^2+{\mathit{\sigma }}_{\mathit{water}}^2+{\mathit{\sigma }}_{\mathit{noise}}^2.(4)$

Variances for the quantities at TOA ${\sigma }_t^2$, ${\sigma }_R^2$, ${\sigma }_a^2$, and ${\sigma }_g^2$ are divided by $t^2$ in accordance with Eqs 1, 2; ${\sigma }_{noise}^2$ includes $1/t^2$ in its definition (Qi et al., 2017) and characterizes the impact of sensor noise, which affects $L_t\left(\lambda \right)$ and through it R_rs(λ). In this work, the transmittance coefficient spectrum is considered not to be dependent on the aerosol optical thickness ${\tau }_a\left(\lambda \right)$, as this dependence is usually small (Wang, 1999). A possible small variability of the Rayleigh optical thickness values, ${\tau }_R\left(\lambda \right)$, which is discussed as a result of this work, is also not considered in the transmittance coefficient. Following Eq. 3, ${\sigma }_R^2$, ${\sigma }_a^2$, ${\sigma }_g^2$, and ${\sigma }_{surf}^2$ contain uncertainties due to both natural variability inside the set of pixels and uncertainties due to inaccuracies of modeling, while ${\sigma }_t^2$ is at least partially due to the vicarious calibration. It can also include other systematic errors due to detectors, polarization effects, and stray light, but these errors are not included in the model.

It was shown in the study by Herrera-Estrella et al. (2021) that estimated ${\sigma }_{noise}$ for the VIIRS sensor (Qi, et al., 2017; Xiong, et al., 2020) is significantly smaller than the total uncertainties $\sigma \left(\lambda \right)$ (Moore et al., 2015). With the representation of data averaged over 3 × 3, 5 × 5, and 7 × 7 pixels considered in this article, noise contribution is further reduced in a way that is inversely proportional to the square root of the number of pixels, in these cases, by 3, 5, and 7 times, respectively, and was therefore not considered further.

All other standard deviation components in Eq. 4, except ${\sigma }_t$, as a first approximation, were considered proportional to the corresponding mean values of the normalized radiances with k as proportionality coefficients:

${\mathit{\sigma }}^2=({\mathit{\sigma }}_{\mathit{vc}}^2+{({\mathit{k}}_{\mathit{R}}{\mathit{L}}_{\mathit{R}}/{\mathit{E}}_{\mathit{d}})}^2+{({\mathit{k}}_{\mathit{a}}{\mathit{L}}_{\mathit{a}}/{\mathit{E}}_{\mathit{d}})}^2+{({\mathit{k}}_{\mathit{g}}{\mathit{L}}_{\mathit{g}}/{\mathit{E}}_{\mathit{d}})}^2)/{\mathit{t}}^2+{\left({\mathit{k}}_{\mathit{S}}\mathit{S\ast \rho }\right)}^2+{\left({\mathit{k}}_{\mathit{Rrs}}\mathit{Rrs}\right)}^2,(5)$

where ${\mathit{\sigma }}_{\mathit{vc}}\left(\mathit{\lambda }\right)={\mathit{\sigma }}_{\mathit{gains}}\left(\mathit{\lambda }\right){\mathit{L}}_{\mathit{t}}\left(\mathit{\lambda }\right)/{\mathit{E}}_{\mathit{d}}\left(\mathit{\lambda }\right)$, ${\mathit{\sigma }}_{\mathit{gains}}\left(\mathit{\lambda }\right)$ is the standard deviation of gains (unitless) for the VIIRS from NASA processing.

${\mathit{L}}_{\mathit{R}}\left(\mathit{\lambda }\right)={\mathit{F}}_0\left(\mathit{\lambda }\right){\mathit{\tau }}_{\mathit{R}}\left(\mathit{\lambda }\right)\mathit{\ast }0.75\mathit{\ast }\frac{1+{\cos }^2\mathit{\Theta }}{4\mathit{\pi }\cos \mathit{\theta }},(\mathrm{6a})$

${\mathit{L}}_{\mathit{a}}\left(\mathit{\lambda }\right)=\frac{{\mathit{\omega }}_0\left(\mathit{\lambda }\right){\mathit{F}}_0\left(\mathit{\lambda }\right){\mathit{\tau }}_{\mathit{a}}\left(\mathit{\lambda }\right){\mathit{P}}_{\mathit{a}}}{4\mathit{\pi }\cos \mathit{\theta }},(\mathrm{6b})$

${\mathit{L}}_{\mathit{g}}\left(\mathit{\lambda }\right)={\mathit{F}}_0\left(\mathit{\lambda }\right){\mathit{T}}_0\left(\mathit{\lambda }\right)\mathit{T}\left(\mathit{\lambda }\right)\mathit{\ast }\mathbf{0.005},(\mathrm{6c})$

${\mathit{E}}_{\mathit{d}}\left(\mathit{\lambda }\right)={\mathit{F}}_0\left(\mathit{\lambda }\right){\mathit{t}}_0\left(\mathit{\lambda }\right)\cos {\mathit{\theta }}_0,(\mathrm{6d})$

${\mathit{L}}_{\mathit{surf}}\left(\mathit{\lambda }\right)=\mathit{\rho }{\mathit{L}}_{\mathit{sky}}\left(\mathit{\lambda }\right)\text{; }\mathit{S}\left(\mathit{\lambda }\right)=\frac{{\mathit{L}}_{\mathit{sky}}\left(\mathit{\lambda }\right)}{{\mathit{E}}_{\mathit{d}}\left(\mathit{\lambda }\right)}.(6e)$

In Eq. 6a, $F_0\left(\lambda \right)$ is the extraterrestrial irradiance, $\theta$ is the sensor zenith angle, and $\Theta$ is the scattering angle, the angle between the solar and viewing directions; in Eq. 6b, ${\omega }_0\left(\lambda \right)$ is the single scattering albedo and $P_a$ is the scattering function for aerosols; in Eq. 6c, $T_0\left(\lambda \right)$ and $T\left(\lambda \right)$ are the direct transmittance coefficients for TOA to surface and surface to TOA, respectively, and 0.005 is the threshold for glint detection $L_{GN}$, in sr⁻¹ (Wang and Bailey, 2001); in Eq. 6d, ${\theta }_0$ and $t_0\left(\lambda \right)$ are the Sun zenith angle and the corresponding diffuse transmittance, respectively; and in Eq. 6e, a representative normalized sky reflectance, $S=L_{sky}/E_d$, was simulated by the VRT RayXP code (Tynes et al., 2001) for the Sun zenith angle ${\theta }_0=42°$ at a viewing zenith angle of 40°, with ${\tau }_a\left(443\right)$ and Angstrom coefficient average values for each specific area based on the numbers from the satellite processing given below in Table 2. While the sky radiance, $L_{sky}$, differs from scene to scene, the normalized sky radiance, $L_{sky}/E_d$, has much less variability; reflectance coefficient was considered as ρ = 0.025, which is the typical reflectance coefficient at 40° viewing angle. Water component ${\sigma }_{water}$ was expressed directly proportional to the remote sensing reflectance, R_rs, recalling its definition as $R_{rs}\left(\lambda \right)=L_w/E_d$.

The uncertainty due to the aerosol component was estimated in two ways, both based on aerosol radiance. The first way is analogous to that of the other radiance components, using normalized radiance with a proportionality coefficient ($k_a{L}_a/E_d$) and $L_a\left(\lambda \right)$ determined from Eq. 6b with ${\omega }_0=1$ for all wavelengths, with aerosol optical thickness values and the Angstrom coefficient derived from the satellite imagery and with the phase function (PF) assumed equal to 0.3 corresponding to the scattering angle around 120°. The second technique used the AERONET data, with the differences between satellite and in situ radiances calculated as differences between the aerosol radiances for the VIIRS from NASA SeaDAS processing software package and radiances calculated from Eq. 6b, where ${\omega }_0$ spectra, aerosol optical thicknesses, Angstrom coefficient, and phase function (PF) values were all derived from the AERONET data. The representative term $\Delta L_a\left(\lambda \right)$ was determined as

$\mathrm{\Delta }{\mathit{L}}_{\mathit{a}}\left(\mathit{\lambda }\right)=\sqrt{\frac{{{\displaystyle \sum }}_{\mathit{i}=1}^{\mathit{N}}{\left({\mathit{L}}_{{\mathit{a}}_{\mathit{SeaDAS}}}^{\mathit{i}}-{\mathit{L}}_{{\mathit{a}}_{\mathit{model}}}^{\mathit{i}}\right)}^2}{\mathit{N}}},(7)$

with ${\sigma }_a=(k_a\Delta L_a/E_d)$ and N being the number of available measurements. The Rayleigh–aerosol interactions were not considered in the model since VRT simulations showed their potential contributions to be very small, 2-5% of the total uncertainties due to aerosol radiances.

Eq. 5 includes VIIRS vicarious calibration uncertainties, where sensor gains are determined by the comparison of the water-leaving radiance from VIIRS after atmospheric correction with in situ measurements at the MOBY site (Franz et al., 2007). ${\sigma }_{gains}\left(\lambda \right)$ is included as a constant spectrum with values ${\sigma }_{gains}\left(\lambda \right)=\left[0.01220,0.01125,0.01157,0.00904,0.00580\right]$ for the corresponding wavelengths 410, 443, 489, 551, and 671 nm, respectively. Data processing was carried out first without ${\sigma }_{vc}\left(\lambda \right)$ in Eq. 5 to understand the contributions of other components to the total $\sigma \left(\lambda \right)$ and then with ${\sigma }_{vc}\left(\lambda \right)$ included.

For each available matchup between the satellite and AERONET-OC measurements, all radiance spectra in Eqs 6a–d were calculated, then spectra were averaged over the total number of available measurements. Mean spectra were then used in the fitting procedure based on Eq. 5, together with $L_{surf}\left(\lambda \right)$ calculated based on the spectrum of the sky radiance, S(λ), for each station. There are 5 bands on the VIIRS sensor, and there are 5 unknown k coefficients in Eq. 5, which are determined from the fitting procedure, and which define the contribution of each term to the total uncertainty.

It should be emphasized that the spectral shapes of the main components in Eqs 6a–e are of primary interest: changes in values, which were assumed constant in the model, do not affect these shapes nor the contribution of the corresponding uncertainties from these components to the total, ${\sigma }^2\left(\lambda \right)$; only the values of the coefficients k will be affected.

$\sigma \left(\lambda \right)$ was calculated from the comparison of the VIIRS data with the corresponding AERONET-OC station data in all the bands and in four wind speed W brackets (W < 3 m/s, 3 < W < 5 m/s, W>5 m/s, and all wind speeds together) as

$\mathit{RMSD}=\sqrt{\frac{{{\displaystyle \sum }}_{\mathit{i}=1}^{\mathit{N}}{\left({\mathit{R}}_{\mathit{r}{\mathit{s}}_{\mathit{sat}}}^{\mathit{i}}-{\mathit{R}}_{\mathit{r}{\mathit{s}}_{\mathit{in}-\mathit{situ}}}^{\mathit{i}}\right)}^2}{\mathit{N}}},(8)$

with $\sigma \left(\lambda \right)=RMSD$.

Biases were also calculated as

$\mathit{bias}=\frac{{{\displaystyle \sum }}_{\mathit{i}=1}^{\mathit{N}}\left({\mathit{R}}_{\mathit{r}{\mathit{s}}_{\mathit{sat}}}^{\mathit{i}}-{\mathit{R}}_{\mathit{r}{\mathit{s}}_{\mathit{in}-\mathit{situ}}}^{\mathit{i}}\right)}{\mathit{N}}.(9)$

Following Eqs 5, 6, coefficients of variation (CVs) can be further determined by normalizing the $\sigma \left(\lambda \right)$ components by $R_{rs}^2\left(\lambda \right)$

$\mathit{C}{\mathit{V}}^2=(\mathit{C}{\mathit{V}}_{\mathit{R}}^2+\mathit{C}{\mathit{V}}_{\mathit{a}}^2+\mathit{C}{\mathit{V}}_{\mathit{g}}^2)/{\mathit{t}}^2+\mathit{C}{\mathit{V}}_{\mathit{surf}}^2+{\mathit{k}}_{\mathit{Rrs}}^2,(10)$

where the total CV on the left side represents $\sigma \left(\lambda \right)/R_{rs}\left(\lambda \right)$.

Optimization Procedure

With all $\sigma \left(\lambda \right)$ components calculated for a scene, a non-linear least-squares fit optimization was carried out in the MATLAB using the default trust-region-reflective algorithm (Coleman and Li, 1994; Coleman and Li, 1996) to determine the respective values of the k coefficients for MOBY and each AERONET-OC site based on the spectra of the $\sigma \left(\lambda \right)$ components and their individual contribution to the total observed R_rs variance, ${\sigma }^2\left(\lambda \right)$, as described in Eq. 5. For each set of the components (defined by a specific site, pixel averaging resolution, and wind speed bracket, for a total of 120 sets), a corresponding set of k coefficients was determined based on VIIRS spectra from 5 bands, that is, 410, 443, 489, 551, and 671 nm. To do so, the components in each set as well as the total observed uncertainty were all normalized to 1 with a simple division by their maximum value across the available spectrum. The normalized spectra were then fed to the optimization algorithm, which was run 100 times for each set to test the robustness of the solution against the possible local minima. For every run, a new randomized set of initial conditions was employed, and the normalized k coefficients were then let to vary between an upper bound of 1 and a lower bound of 0.02, to encourage the avoidance of non-physical solutions. The robustness of each solution was evaluated by the ratio of the standard deviation to the mean of the normalized k coefficients across all 100 outputs, with the mean being used as the final value of the coefficients. For all 5 k coefficients within each set, and across all 120 sets in the dataset, the median of the ratio was found to be of the order of 10^–6 or smaller, and its maximum to be of the order of 10^–3 or smaller, saved for one single set out of 120, where the ratio corresponding to the Rayleigh k coefficient was of the order of 10^–1. Overall, these statistics indicate strong convergence in the optimization results. No further constraints were applied to the possible solutions, and in particular, no expectations in terms of typical relative magnitudes of the components were used to direct the algorithm, in accordance with the exploratory nature of our study. Once a solution was reached, all normalized k coefficients were then scaled back to their true scale values, using the formula

${\mathit{k}}_{\mathit{i}=}\frac{{\mathit{m}}_{\mathit{\sigma }}{\tilde{\mathit{k}}}_{\mathit{i}}{\tilde{\mathit{\sigma }}}_{\mathit{i},\mathit{ max}}}{{\mathit{\sigma }}_{\mathit{i},\mathit{ max}}},(11)$

where ${\mathit{k}}_{\mathit{i}}$ are the true scale k coefficients, index i represents the various components, $m_{\sigma }$ is the normalization coefficient of the total uncertainty $\sigma \left(\lambda \right)$, ${\tilde{\mathit{k}}}_{\mathit{i}}$ are the normalized k coefficients determined by the optimization, and ${\tilde{\sigma }}_{i, max}$ and ${\sigma }_{i, max}$ are the spectral maxima of the normalized and true scale versions of the i^th $\sigma \left(\lambda \right)$ component, respectively.

These coefficients, once used to scale their corresponding radiance components, were interpreted as an indication of the major contributing components to the total observed R_rs variance ${\sigma }^2\left(\lambda \right)$.

Visible Infrared Imaging Radiometer Suite Satellite and AERONET-OC Data

VIIRS Data

Satellite imagery was downloaded for the period from January 2012 to October 2021 for the area of the Marine Optical BuoY (MOBY) in Hawaii and eight Aerosol Robotic Network for Ocean Color (AERONET-OC) sites (Figure 1): the University of South California (USC), Venise, Gloria, the Martha’s Vineyard Coastal Observatory (MVCO), COVE, WaveCIS, the Long Island Sound Coastal Observatory (LISCO), and the Helsinki Lighthouse (HLT).

FIGURE 1

FIGURE 1. Areas of study: global map showing the MOBY area and all the AERONET-OC stations.

VIIRS’s Satellite Level 2 imagery, version 2018.0, was downloaded from the NASA Ocean Color website https://oceancolor.gsfc.nasa.gov (Gordon and Wang 1994; Siegel et al., 2000; Bailey et al., 2010). Standard NASA Level 2 data files for the VIIRS include geophysical products of the atmosphere and ocean, such as aerosol optical thickness, remote sensing reflectance, R_rs(λ), in the visible wavelengths 410, 443, 486, 551, and 671 nm, and the level 2 quality flags. However, Sun zenith angle, sensor viewing angle, sensor azimuth angle, scattering angle, total radiance, and aerosol radiance are obtained from SeaDAS version 7.5.3 after processing VIIRS’s Satellite Level 1A imagery for those files that passed the matchup selection.

Pixels flagged by at least one of the following conditions were excluded: land, cloud, failure in the atmospheric correction, stray light (except for LISCO), bad navigation quality, high or moderate glint, negative Rayleigh-corrected radiance, negative water-leaving radiance, viewing angle larger than 60°, and solar zenith angle larger than 70°.

The VIIRS’s pixel resolution for the reflectance bands at nadir is 750 m. A file is selected if at least half of the pixels in the set plus one was flag-free. Pixels used for matchup comparison were averaged over 3 spatial resolutions: 2250, 3750, and 5250 m (3 × 3, 5 × 5, and 7 × 7 pixel boxes), centered at the AERONET site (Hlaing et al., 2013). Average R_rs(λ) and the standard deviation among pixels, geometry, and radiance were recorded.

In addition, aerosol radiances were also downloaded from SeaDAS for the comparison with aerosol radiances simulated based on parameters from AERONET-OC sites.

AERONET-OC Data

The ocean color component of the Aerosol Robotic Network (AERONET-OC) was implemented to support long-term ocean color investigations by collecting normalized water-leaving radiance and aerosol optical depth data using the SeaPRISM autonomous radiometer systems deployed on offshore fixed platforms (Zibordi et al., 2009; Zibordi et al., 2020). The SeaPRISM system is a CIMEL Electronique CE-318/CE-318T sunphotometer, used to retrieve atmospheric optical thickness and other atmospheric parameters, and modified to perform radiance measurements with a full-angle field of view of 1.2° to determine the total radiance from the sea surface, $L_t^A\left(\text{λ}\right)$, and the sky radiance, $L_{sky}^A\left(\text{λ}\right)$, as a function of solar zenith angle, sensor viewing angle, and relative azimuth with respect to the sun (Zibordi et al., 2009). From these, the normalized water-leaving radiance, $L_w^A\left(\text{λ}\right)$, and the remote sensing reflectance were determined. The wavelengths of the visible spectrum used in this analysis are 412, 443, 488, 551, and 667 nm from CE-318.

The aerosol optical depth, aerosol inversions, and ocean color data used in this analysis are version 3 level 1.5 data, which has been cloud-screened and quality-controlled to ensure the accuracy of the data. All matchups were observed within a ±2 h window between the satellite overpass and in situ observation (Zibordi et al., 2009; Zibordi et al., 2020). The spectrum was classified as coastal water if R_rs (412) ≤ 0.006, and as open ocean water if R_rs (412) > 0.006. More detailed information on AERONET-OC sites is listed in Table 1. At the Venise site, there were a large number of observations with both R_rs (412) ≤ 0.006 and R_rs (412) > 0.006, and they are presented separately. Due to the site location, Venise waters with R_rs (412) > 0.006 are far from clear sea waters, but the term remained to formally separate different types of the spectra. At a few other sites, there were a small number of spectra with R_rs (412) > 0.006; these were not considered because there were not enough data for reliable averaging and fitting procedures.

TABLE 1

TABLE 1. Location and parameters of the AERONET-OC sites.

The quality of R_rs(λ) data depends on the wind speed (Zibordi et al., 2009; Zibordi et al., 2020). Due to the accepted algorithm to minimize the impact of Sun glint on the measurements of the total above water radiance, which takes into account the lowest two out of 11 measurements, with increasing wind speeds above 5 m/s, water-leaving radiance becomes slightly lower. As can be seen in the following text, for several stations, that increased $\sigma \left(\lambda \right)$ for wind speed W>5 m/s and affected the bias in R_rs(λ), but there were exceptions from this trend, and this effect did not significantly change the whole picture of the uncertainties’ composition. Consequently, the results are presented mostly for the data averaged over all wind speed ranges; a few details are further discussed in the following section.

Marine Optical BuoY Data

Marine Optical BuoY (MOBY) radiometry data are used by the NASA-OBPG as part of ocean color validation and vicarious calibration activities (Clark et al., 1997). MOBY is an autonomous anchored buoy offshore of Lanai, Hawaii. On each day of deployment, it collects several measurements of upwelling radiance from sensors on its underwater arms (at approximately 1, 5, and 9 m depth) and downwelling irradiance from sensors on its underwater arms as well as at the surface (Voss et al., 2017).

From the MOBY “gold” directory, the MOBY data that matched the bands from the VIIRS were collected when the existing data were matched up with ±2 h of the satellite overpass.

The main atmospheric parameters at the studied sites determined from the AC processing and AERONET retrievals are provided in Table 2. Absorbing aerosols were noticeable at several sites with the average ${\omega }_0$ values even below 0.9 at LISCO, but the spectral dependence of ${\omega }_0$ was small.

TABLE 2

TABLE 2. Average atmospheric parameters at the sites of the study determined from satellite and AERONET retrievals.

Results

Mean R_rs(λ) spectra from the MOBY site and eight AERONET-OC stations are shown in Figure 2A. They represent water areas from very clear (MOBY and USC), moderate coastal (Venise, Gloria, MVCO, WaveCIS, and COVE), and very coastal (LISCO and HLT) waters with R_rs(λ) standard deviations shown in Figure 2B. Corresponding $\sigma \left(\lambda \right)$ spectra are shown in Figure 2C. They have different shapes with the highest values in the blue, which likely suggests the presence of both different and common spectral components in these total spectra.

FIGURE 2

FIGURE 2. (A) Mean R_rs (λ) spectra from the AERONET-OC, (B) standard deviations of R_rs (λ), and (C) corresponding uncertainties spectra $\sigma \left(\text{λ}\right)$ for all the areas of study.

First, the results of processing Eq. 5 without the ${\sigma }_{vc}\left(\lambda \right)$ term are reported. As explained previously, the fitting procedure included the sum of all the spectral components in Eq. 5 with the corresponding coefficients k. Typical spectra for all the components involved, normalized to their maxima, are shown for the open ocean and coastal water stations in Figure 3. In accordance with Eq. 4, ${\sigma }_R$, ${\sigma }_a$, and ${\sigma }_g$ are divided by the spectrum of the diffuse transmittance t for the propagation of light from the surface to TOA. The normalized spectra for the Rayleigh scattering and glint are the same for all the stations, and the spectra for the total $\sigma \left(\lambda \right)$, for aerosols, surface effects, and R_rs(λ) components are different. The Rayleigh scattering and surface effects spectra are both related to the sky spectra, but the former is divided by the spectrum of the diffuse transmittance t and the latter was simulated based on the composition of the Rayleigh and aerosol scattering, which makes these spectra distinct from each other.

FIGURE 3

FIGURE 3. Normalized spectra for the MOBY, USC, Venise, and LISCO sites (3 × 3 pixels, 2250 m resolution).

Results of the fitting are presented in Figure 4 for the case when ${\sigma }_a\left(\lambda \right)$ was considered proportional to the aerosol radiance spectra from Eq. 6b, and in Figure 5 for the case when ${\sigma }_a\left(\lambda \right)$ was considered proportional to the difference of aerosol radiances $\Delta L_a\left(\lambda \right)$ from SeaDAS and in situ (normalized by the downwelling irradiance). It was found that this normalized difference was about 5 times greater than the total $\sigma \left(\lambda \right)$, meaning that this difference does not fully represent uncertainties in aerosol radiances in the atmospheric correction process, so the coefficient ${\mathit{k}}_{\mathit{a}}$ was allowed to vary in a wide range. Note that the smaller number of measurements included in the averaging in Figure 5 in comparison with Figure 4 was due to the unavailability of some atmospheric parameters necessary to model aerosol radiances based on in situ data.

FIGURE 4

FIGURE 4. Results of fitting for all areas of the study. $\sigma \left(\text{λ}\right)$ was considered proportional to the aerosol radiance spectra from Eq. 6b.

FIGURE 5

FIGURE 5. Results of fitting for all areas of the study. $\sigma \left(\text{λ}\right)$ was considered proportional to the differences between the aerosol radiances from the SeaDAS and AERONET-OC sites (except for the MOBY site).

While there are some differences in the spectral components for these two cases, the main results remain the same. The main uncertainty for all the areas of AERONET-OC stations from clear waters to very coastal comes from the Rayleigh component, which is typically assumed to be well defined in the atmospheric correction process (Mobley et al., 2016) and to be pre-calculated in a very accurate manner. The presence of the strong Rayleigh component in the total $\sigma \left(\lambda \right)$ budget explains well the maximum of uncertainties in the blue part of the spectrum. For most cases, the maximum of ${\sigma }_R$ in the blue part is in a relative narrow range of 0.8–1.4 × 10^–3 sr⁻¹, leading to the conjecture that it is related to the specific physical effect.

The second major impact on the total $\sigma \left(\lambda \right)$ comes from the ${\sigma }_{water}$ term, which is due to the variability of water parameters spectrally proportional to R_rs(λ), and this brings a significant difference in the spectral shape of the total $\sigma \left(\lambda \right)$ depending on the site and the water type; surface effects ${\sigma }_{surface}$ are also noticeable at some stations. The contribution of the aerosol uncertainty component was found to be virtually null at almost all sites despite the fact that all data were used in the processing, including outliers such as R_rs(λ) spectra with negative values at 412 nm. It should be reminded that the fitting procedures were carried out with radiances corresponding to the mean values for each station in Table 2, where all ${\tau }_a\left(443\right)$ are below 0.2, which does not exclude higher contribution of aerosol radiance uncertainties in specific cases with higher ${\tau }_a\left(443\right)$ values. It should also be noted that in Figures 4, 5, standard deviations are presented, while the fitting process is carried out in variances, so small differences in the values of spectral components in Figures 4, 5 are amplified in the fitting procedure. From this point of view, the contribution of glint to the total $\sigma \left(\lambda \right)$ is mostly small in comparison with that of other components previously discussed. As previously mentioned, there was no filtering of any individual spectra of uncertainties, which contributed to $\sigma \left(\lambda \right)$ in Eq. 8 as long as all the necessary data were available; that includes outliers, which induced sharp features in some of the $\sigma \left(\lambda \right)$ spectra (WaveCIS, Gloria) and resulted in worse fitting results than for other stations.

Another important thing to be noted is that the optimization algorithm cannot easily distinguish between the components with a similar spectral behavior, particularly in the case of Rayleigh scattering and surface effects, and, to a lesser degree, aerosol scattering (cf. Figure 3). By the very definition of optimization, the procedure maximizes the efficiency in reconstructing the total uncertainty, which may include the singular utilization of one component at the expense of all others. This too may explain the general small contribution of the aerosol component, as well as sudden switches between the Rayleigh scattering and surface contributions like in the case of the WaveCIS site when the formulation of $\Delta L_a\left(\lambda \right)$ as the difference between SeaDAS and in situ values is used (Figure 5). Nevertheless, while possibly falling short of offering a true estimation of the real relative magnitudes of the individual contributions to $\sigma \left(\lambda \right)$, the consistent identification of the Rayleigh scattering as the component best capable of explaining $\sigma \left(\lambda \right)$ in the 400–500 nm range appears strongly suggestive of a much more critical role for this contribution than usually assumed in the literature.

The spectral composition of $\sigma \left(\lambda \right)$ at the MOBY site is different from the composition at AERONET-OC sites. The vicarious calibration of VIIRS and other NASA OC sensors is carried out at this site, meaning that all differences between the atmospherically corrected water-leaving radiance spectra and measured by MOBY are corrected by the calibration gains. Nevertheless, a strong Rayleigh component also exists in the uncertainty budget at this site.

As was shown before (Herrera- Estrella et al., 2021), there is almost no spatial variability of water measured by R_rs(λ) in the MOBY area. However, the second main spectral component in $\sigma \left(\lambda \right)$ at the MOBY site is the spectrum proportional to R_rs(λ), which is probably due to the temporal variability of the water parameters at this site. Inaccuracy of the radiances in the atmospheric correction at the MOBY site is most likely compensated by the ${\tau }_a\left(\lambda \right)$ values, which are usually higher than expected from field measurements.

Spectral biases of R_rs(λ) at all the sites are shown in Figure 6, with the global mean bias shown as a dashed line. This was calculated from all spectra excluding the ones for the MOBY site and from the LISCO site since the stray light flag was suspended for LISCO. This bias spectrum resembles the bias at the MOBY site, which will be smaller in the blue bands if the mean bias is subtracted. Bias at MOBY was already reduced in the VIIRS reprocessing 2018 in comparison with 2014 reprocessing (Franz et al., 2018), which improved the quality of NASA products.

FIGURE 6

FIGURE 6. Biases of R_rs(λ) in areas of study; dashed line is the mean bias.

As was mentioned before, the results of the fitting were slightly different for different wind speed brackets. These differences are demonstrated in Figure 7 where fitting results are presented for Venise and LISCO sites for three wind speed intervals as well as averaged over all cases.

FIGURE 7

FIGURE 7. Separation of the fitting results by the wind speed, Venise, and LISCO sites: Venise (top row), LISCO (middle row), $\sigma \left(\text{λ}\right)$ and bias for Venise and LISCO (bottom row). The last graph in the first and second rows represents data averaged over all wind speed intervals.

Surface effects were small at 3 < W ≤5 m/s and more noticeable at W > 5 m/s. This pattern was typical for most of the AERONET-OC sites. Some presence of the surface effects at low wind speeds W < 3 m/s is probably associated with the specifics of processing of the SeaPRISM data and/or accuracy of Cox–Munk surface slopes model (Cox and Munk, 1954) for the coastal sites in the satellite data processing. The uncertainty of the Rayleigh component is the highest in most cases, especially in the blue part of the spectra followed by water variability components and sometimes surface effects. $\sigma \left(\lambda \right)$ is higher at high wind speed W > 5 m/s for both stations, which can be partially due to the AERONET-OC processing algorithm also pushing biases’ spectra slightly higher for this wind speed range.

Fitting coefficients k from Eq. 5 for all stations are provided in Figure 8 as a function of the spatial resolution based on 3 × 3, 5 × 5, and 7 × 7 pixel processing by the model. There are no significant changes in the coefficients with an increase of the averaged area. The coefficients have different ranges of values for different parameters: they are smaller for $L_r$ and $L_g$, which are measured at the TOA level, and greater for surface effects and water variability, which at the surface level have radiances about 10 times smaller than $L_r$. The ${\mathit{k}}_{\mathit{a}}$ coefficients are shown for two different processing schemes, and the other coefficients for these two schemes remained in the same range. The ${\mathit{k}}_{\mathit{R}}$ coefficient related to the Rayleigh uncertainties is grouped in a very small range around 0.015, meaning that it is related to a specific effect common to all the stations, which further inspired processing with all terms in Eq. 5 including ${\sigma }_{vc}\left(\lambda \right)$. The ranges of the ${\mathit{k}}_{\mathit{a}}$ coefficient related to aerosols are slightly different for the two processing schemes because of different aerosol radiances represented in the fitting process and mostly represent a small contribution of the aerosol component as is seen in Figures 4, 5. ${\mathit{k}}_{\mathit{g}}$ varied in a broad range but as is seen in Figure 4, the contribution of glint effects was not significant at all the studied sites. Surface effects, characterized by the ${\mathit{k}}_{\mathit{S}}$ coefficient, were noticeable at many stations but with mostly small contribution: the highest value was found at the WaveCIS site in the second processing case (Figure 5, bottom row). For this station, the optimized solution surface effects probably replaced the Rayleigh component, which was dominant in the first processing case (Figure 4, bottom row). ${\mathit{k}}_{\mathit{Rrs}}$, which is related to water variability, was small in clear USC waters, about 0.1 at the MOBY site probably due to temporal variability and varied in a small 0.1–0.3 range at other sites, typically with a small change in the spatial resolution. The water variability component was prominent at all the sites and affected the shape of the uncertainties, as described previously.

FIGURE 8

FIGURE 8. k coefficients vs. resolution for all the stations.

Total uncertainties $\sigma \left(\lambda \right)$ and uncertainties related to the R_rs spatial variability [defined as ${\sigma }_{spat}\left(\lambda \right)$ in the study by Herrera-Estrella et al. (2021)] are well spectrally correlated in the open ocean (R² = 0.95–0.98), almost independently of the number of pixels averaged. In that study, ${\sigma }_{spat}\left(\lambda \right)$ was found to be 4–5 times lower than $\sigma \left(\lambda \right)$. These uncertainties are much less correlated in coastal waters with high water variability (typically R² = 0.70–0.80), and averaging over at least about 50 pixels is required to reach this R² level as the components are not directly related to the spatial variability to gradually propagate to ${\sigma }_{spat}\left(\lambda \right)$ statistics.

The results of data processing with ${\sigma }_{vc}\left(\lambda \right)$ included are further reported and discussed in the following. Normalized spectra for all the components are shown in Figure 9 similar to Figure 3 for the case without ${\sigma }_{vc}\left(\lambda \right)$. It can be seen that at almost at all sites, ${\sigma }_{vc}\left(\lambda \right)$ and the Rayleigh scattering spectra are very similar.

FIGURE 9

FIGURE 9. Normalized spectra for all sites (3 × 3 pixels, 2250 m resolution) with ${\sigma }_{vc}\left(\text{λ}\right)$ included.

The results of fitting with the inclusion of ${\sigma }_{vc}\left(\lambda \right)$ for the case where ${\sigma }_a\left(\lambda \right)$ is considered proportional to the difference $\Delta L_a\left(\lambda \right)$ of SeaDAS, and in situ aerosol radiances (normalized by the downwelling irradiance) are shown in Figure 10 for all the stations (similar to Figure 5). The fitting results for the first processing case (analogous to Figure 4) were less accurate and are not shown.

FIGURE 10

FIGURE 10. Results of fitting for all areas of the study with ${\sigma }_{vc}\left(\text{λ}\right)$ included: ${\sigma }_a\left(\text{λ}\right)$ was considered proportional to the differences between the aerosol radiances from the SeaDAS and AERONET-OC sites (except for the MOBY site).

Now the main uncertainty for all the AERONET-OC stations from clear water to very coastal water comes from the ${\sigma }_{vc}\left(\lambda \right)$ component, which results from the vicarious calibration of the sensor at the MOBY site (Franz et al., 2007). This component is slightly different at different sites because of differences in the average $L_t\left(\lambda \right)$ and $E_d\left(\lambda \right)$ values used in the transformation of ${\sigma }_{gains}\left(\lambda \right)$ to ${\sigma }_{vc}\left(\lambda \right)$. As aforementioned, the presence of a strong water variability term and weaker surface effects, aerosol, and glint components is observed.

CV spectra for all sites are shown in Figure 11 for both in situ and satellite data, with a typical increase of CV in the blue bands reaching about 2–2.5 at the very coastal sites where R_rs values are especially low. As expected, CVs are also high in the red part of the spectra for the clear water sites like USC because of the low R_rs values at those wavelengths.

FIGURE 11

FIGURE 11. CV spectra for all sites of the study: the in situ data (left) and the VIIRS data (right).

SeaPRISM measurements together with data processing introduce their own uncertainties, which need to be considered in the uncertainties’ budget. According to Gergely and Zibordi (2014), the CV for the Venise site at 412 nm is about 5% and at the Helsinki site is about 27% with the differences mostly due to the different R_rs values. Considering the R_rs spectra shown in Figure 2A, this corresponds to $\sigma \approx 2\times {10}^{-4}$ sr⁻¹, and thus makes a small contribution to the total $\sigma$. However, if some method is found to reduce the impact of the Rayleigh component on the total uncertainties, the contribution can become more significant.

Distributions of the individual spectra $\Delta R_{rs}=R_{rssat}^i-R_{rsin-situ}^i$ for the three stations (USC, Venise, and MVCO) are shown in the first column of Figure 12. Modified spectra after the removal of bias are shown in the central column of the same figure, demonstrating a strong symmetry of the spectra against the zero line. As a reminder, no outliers were removed at any of the stations as long as all the necessary data were fully available. In the third column of Figure 12, histograms of $\Delta R_{rs}\left(412\right)$—bias are shown, with distributions very close to normal in all the three cases. A similar effect was noticed in IOCCG (2019) for the comparison of data from MODIS-Aqua at the Venise site. It was also mentioned that uncertainties at different wavelengths are well correlated. All these features are consistent with the presence of the Rayleigh-type component as the main component in the total $\sigma \left(\lambda \right)$.

FIGURE 12

FIGURE 12. Analysis of the $\Delta R_{rs}$ distribution for MVCO, USC, and Venise sites. First column: all $\Delta R_{rs}$, second column: all $\Delta R_{rs}$—mean bias for the station, third column: histogram of all $\Delta R_{rs}$—mean bias for the station at 412 nm.

Finally, the $\Delta R_{rs}\left(412\right)$ distribution was checked in terms of correlations with $R_{rsin-situ}\left(412\right)$ for the in situ data. Almost no correlation (R² < 0.01) was found, indicating that $\Delta Rrs$ and total $\sigma \left(\lambda \right)$ are not related to the uncertainties of the in situ measurements. $\Delta R_{rs}\left(412\right)$ were instead correlated with $R_{rssat}$, where correlations were up to R² = 0.93 at the LISCO site.

Discussion and Conclusion

The model developed for the separation of remote sensing reflectance uncertainties into their spectral components is applied to the uncertainties’ spectra from the matchups of satellite and AERONET-OC data and MOBY measurements. It was shown that the main component in R_rs uncertainties at all the AERONET-OC sites is the Rayleigh-type component at the level of (0.8–1.4) ×10^–3 sr⁻¹ at 412 nm based on Figure 5, which with $E_d\left(412\right)\approx 100$ mWcm⁻²μm⁻¹sr⁻¹ corresponds to a standard deviation in radiance of about 0.12 mWcm⁻²μm⁻¹sr⁻¹. This uncertainty for $L_t\left(412\right)=$ 8–10 mWcm⁻²μm⁻¹sr⁻¹ is about ±1.2–1.5% of the total TOA radiance. The contribution of this component to total uncertainties is not constant, and the simple inclusion of ${\sigma }_{vc}\left(\lambda \right)$ does not account for the full range of uncertainties as is shown in Figure 10. This component is different at different stations.

A comparison of ${\sigma }_{gains}\left(\lambda \right)$ applied to the MOBY conditions from different satellite sensors: SeaWIFS, MODIS (Moore et al., 2015), Sentinel 3A OLCI (Lamquin et al., 2017), and VIIRS (NASA data) show significant similarity of these spectra, as demonstrated in Figure 13. It is unlikely that such similarity would exist for different sensors with different designs if ${\sigma }_{gains}\left(\lambda \right)$ were due to an imperfect sensing process. It is more likely that this variability is due to the variability in the atmosphere, with the main variability due to Rayleigh scattering. With ${\sigma }_{gains}\left(\lambda \right)$ applied not to $L_r\left(\lambda \right)$ but to $L_t\left(\lambda \right)$ at the MOBY site, these ${\sigma }_{gains}\left(\lambda \right)$ cannot be directly applied to other stations with different contributions of $L_r\left(\lambda \right)$ to $L_t\left(\lambda \right)$, primarily due to the different $L_w\left(\lambda \right)$ component at MOBY and other areas, as is visible in Figures 9, 10.

FIGURE 13

FIGURE 13. Comparison of ${\sigma }_{vc}\left(\text{λ}\right)$ for different satellite sensors.

The main parameter in $L_r\left(\lambda \right)$ that can force $L_r\left(\lambda \right)$ to fluctuate in this manner is the Rayleigh optical thickness ${\tau }_R\left(\lambda \right)$. In the current atmospheric correction processing, ${\tau }_R\left(\lambda \right)$ is calculated based on the study by Bodhaine et al. (1999). Several approaches in the calculation of ${\tau }_R\left(\lambda \right)$ were considered earlier (Teillet, 1990) with the differences of a few percent. Moreover, the natural variability of ${\tau }_R\left(\lambda \right)$ can be imagined due to the variability of the concentrations of the main gaseous components and their vertical distribution, as well as differences between the actual and estimated pressure and temperature values, that are each typically on the order of 0.2–0.3% of their actual values (Smith et al., 2001). This looks consistent with the close to normal distribution of ΔR_rs in Figure 12. The standard deviation of ${\tau }_R\left(\lambda \right)$, which may be considered as a measure of this variability, is about 1.5% of ${\tau }_R\left(\lambda \right)$, assuming $L_r\left(412\right)\approx 6.0$ mWcm⁻²μm⁻¹sr⁻¹ and transmittance coefficient $t\left(412\right)=0.8$. Rayleigh component uncertainties can also come from the interaction in the scattering process of molecules and aerosols. However, as previously mentioned, this impact should be small.

Variability of $L_r\left(\lambda \right)$ due to the changes in ${\tau }_R\left(\lambda \right)$ has a direct impact on the atmospheric correction. While the changes are too small in the NIR to affect the selection of aerosol models (Wang and Gordon, 2002), changes of ${\tau }_R\left(\lambda \right)$ in the blue part of the spectrum are stronger than the effect of absorbing aerosols. Since the Rayleigh component is pre-calculated, negative changes in the actual Rayleigh contribution will result in overcorrection, and this effect together with the low R_rs (412) values in coastal waters is probably the main factor that creates negative R_rs (412) in the atmospheric correction processing. In the study by Ransibrahmanakul and Stumpf (2006), negative values of R_rs (412) from SeaWiFS in Long Island Sound are shown to exist about 50% of the time, which is consistent with the probability of negative changes in the Rayleigh component but not with the frequency of very strong absorbing aerosols in the region. A power law-like artifact with exponent −6 was identified and corrected in R_rs (λ) in that work to eliminate an impact of negative R_rs (λ). Based on ${\omega }_0$ values in Table 2, there are no spectral changes of ${\omega }_0$ that can create a spectral component of R_rs error with power −6, but it was observed that uncertainties due to the Rayleigh component with power ∼ −4.5 (Rayleigh spectra divided by the diffuse transmittance) give similar results.

Thus, negative values of the Rayleigh component uncertainty, possibly together with the effect of absorbing aerosols at low actual R_rs (412) values, create negative R_rs (412) in the atmospheric correction processing. Similar effects with both positive and negative uncertainties most likely exist at higher R_rs (412) as well, inducing variability in R_rs (λ), but they are not explicitly visible as errors. They can just slightly affect the estimation of chlorophyll concentrations and water inherent optical properties.

The time series of $\Delta R_{rs}\left(412\right)$, which is in accordance with the abovementioned results can be considered as a proxy to the time series of ${\tau }_R\left(412\right)$ for the studied sites, are presented in Figure 14.

FIGURE 14

FIGURE 14. Time series of $\Delta R_{rs}\left(412\right)$ for the areas of the study.

It can be seen that $\Delta R_{rs}\left(412\right)$ uncertainties are not random. Rather, they represent gradual changes in the atmosphere if the measurements are close to each other in time, which is consistent with the assumption that they are mostly due to the changes in the total air column integrated mostly in the Rayleigh optical thickness, not in inaccuracies in the atmospheric correction process. Variability of the Rayleigh component also probably contributes to the recently reported seasonal bias between satellite observations and measurements at the MOBY site (Bisson et al., 2021).

It is possible that the methodology applied in this work underestimates the contribution of aerosols to the total uncertainties because of the difficulty to determine actual spectra of aerosol uncertainties and their spectral variability in various conditions. This contribution should be further analyzed using different approaches. Part of the uncertainties can come from inaccuracies in the computations of the Rayleigh and other components by vector radiative transfer codes. While current VRT codes calculate the Rayleigh component with uncertainties of a fraction of a percent (typically below 0.1%) (Kokhanovsky et al., 2010), related uncertainties should be carefully monitored.

One source of additional uncertainties in $\Delta R_{rs}\left(412\right)$ was recently discussed (Gilerson et al., 2018), and it is due to the combination of the Rayleigh component with surface effects, which is pre-calculated based on the Fresnel reflectance coefficient of the sky light reflected from the ocean surface. In the presence of aerosols in the sky, parameters of which are determined in the further steps of the atmospheric correction after the subtraction of the Rayleigh component, the reflectance coefficient changes. This effect is not accounted for in the Rayleigh component computations, thus creating circular relationships, that could lead to the abovementioned uncertainties. Simple estimations show that $L_{sky}$ and $L_r$ are of the same order. In the surface effects, $L_{sky}$ is multiplied by reflectance coefficient ρ, which is equal to 0.025–0.030 at viewing angles of 40⁰ and smaller, and is greater for larger viewing angles so that changes of ρ by 20–30% can account for about 1% of $L_r$. Similar uncertainties can appear because of deviations of the wave slopes distribution in coastal areas from the distribution of Cox and Munk (1954), which leads to changes in the reflectance coefficient. Surface effects were pronounced at several stations as can be seen in Figures 4, 5, 10. Spectral shapes of the Rayleigh component and of surface effects are close to each other (in the former, the Rayleigh spectra are divided by the transmittance, and in the latter, sky spectra include small effects of aerosols), and they can be complementary in the fitting procedure and thus not well distinguishable from each other.

As mentioned earlier, uncertainties can also include other errors due to detectors, polarization effects, stray light, etc., and such errors were not included in the model. Based on the similarity of the spectra of standard deviations of the vicarious gains from several satellite sensors shown in Figure 13, the contributions of these effects appear to be small for these sensors, which does not exclude the possibility that these terms will be larger for other sensors with different designs or radiometric performance.

Without the uncertainty due to changes of ${\tau }_R\left(\lambda \right)$, the total uncertainty could be reduced by several times, for example, at the USC site in the blue bands, and to a smaller degree but also significantly at the other coastal sites. To minimize this uncertainty, approaches based on the fitting of reflectance spectra (Steinmetz et al., 2011) should be explored more thoroughly.

Based on the analysis of biases of R_rs in Figure 6, the bias at the MOBY site is the highest and is positive. At the AERONET-OC sites, biases are positive at a few stations and negative at a few others. The reasons for such variability should be further studied. Biases in the 412–489 nm range can be partially improved by the slight change of gain values at these wavelengths.

In addition to the MOBY site, the USC site can also be considered for the vicarious calibration. It has minimal water variability and provides full information on aerosol parameters. This information together with the measured water-leaving radiance can be used for the calibration through the simulation of the total radiances (Bailey et al., 2010; Hlaing et al., 2014) in addition to the current approaches of the vicarious calibration (Franz et al., 2007; Werdell et al., 2007).

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

AG, EH-E and RF formulated the original concept of the model. EH-E processed satellite and AERONET-OC data. All the authors participated in the uncertainty analysis from different components and provided critical feedback to the final manuscript.

Funding

The funding of this study was provided by the NOAA JPSS Cal/Val and JPSS PGRR programs, CESSRST Grant, NA16SEC4810008, and the NASA OBB Grant, 80NSSC21K0562.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank NASA OBPG for the satellite imagery, and NASA AERONET Group for processing of the satellite data and support of the operation at the LISCO site. We are grateful to Kenneth Voss and the MOBY team for the MOBY data, PIs of AERONET-OC sites: Giuseppe Zibordi (Venise, Gloria, HLT), Brent Holben (COVE), Heidi Sosik and Hui Feng (MVCO), Robert Arnone, Alan Wiedemann, Bill Gibson, and Sherwin Ladner (WaveCIS). We also thank the Graduate Center of the City University of New York, The City College of New York, NOAA Center for Earth System Sciences and Remote Sensing Technologies, and NOAA Office of Education, Educational Partnership Program for fellowship support for EH-E. The statements contained within the research article are not the opinions of the funding agency or the U.S. government but reflect the authors’ opinions. We are grateful to two reviewers whose recommendations led to the improvement of the manuscript.

References

Ahmad, Z., Franz, B. A., McClain, C. R., Kwiatkowska, E. J., Werdell, J., Shettle, E. P., et al. (2010). New Aerosol Models for the Retrieval of Aerosol Optical Thickness and Normalized Water-Leaving Radiances from the SeaWiFS and MODIS Sensors over Coastal Regions and Open Oceans. Appl. Opt. 49, 5545. doi:10.1364/AO.49.005545

PubMed Abstract | CrossRef Full Text | Google Scholar

Antoine, D., d'Ortenzio, F., Hooker, S. B., Bécu, G., Gentili, B., Tailliez, D., et al. (2008). Assessment of Uncertainty in the Ocean Reflectance Determined by Three Satellite Ocean Color Sensors (MERIS, SeaWiFS and MODIS-A) at an Offshore Site in the Mediterranean Sea (BOUSSOLE Project). J. Geophys. Res. 113, C07013. doi:10.1029/2007JC004472

CrossRef Full Text | Google Scholar

Bailey, S. W., Franz, B. A., and Werdell, P. J. (2010). Estimation of Near-Infrared Water-Leaving Reflectance for Satellite Ocean Color Data Processing. Opt. Express 18, 7521–7527. doi:10.1364/OE.18.007521

PubMed Abstract | CrossRef Full Text | Google Scholar

Bisson, K. M., Boss, E., Werdell, P. J., Ibrahim, A., Frouin, R., and Behrenfeld, M. J. (2021). Seasonal Bias in Global Ocean Color Observations. Appl. Opt. 60 (23), 6978–6988. doi:10.1364/ao.426137

PubMed Abstract | CrossRef Full Text | Google Scholar

Bodhaine, B. A., Wood, N. B., Dutton, E. G., and Slusser, J. R. (1999). On Rayleigh Optical Depth Calculations. J. Atmos. Oceanic Technol. 16, 1854–1861. doi:10.1175/1520-0426(1999)016<1854:orodc>2.0.co;2

CrossRef Full Text | Google Scholar

Carrizo, C., Gilerson, A., Foster, R., Golovin, A., and El-Habashi, A. (2019). Characterization of Radiance from the Ocean Surface by Hyperspectral Imaging. Opt. Express 27 (2), 1750–1768. doi:10.1364/oe.27.001750

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, D. K., Gordon, H. R., Voss, K. J., Ge, Y., Broenkow, W., and Trees, C. (1997). Validation of Atmospheric Correction over the Oceans. J. Geophys. Res. 102, 17209–17217. doi:10.1029/96JD03345

CrossRef Full Text | Google Scholar

Coleman, T. F., and Li, Y. (1996). An Interior Trust Region Approach for Nonlinear Minimization Subject to Bounds. SIAM J. Optim. 6 (2), 418–445. doi:10.1137/0806023

CrossRef Full Text | Google Scholar

Coleman, T. F., and Li, Y. (1994). On the Convergence of interior-reflective Newton Methods for Nonlinear Minimization Subject to Bounds. Math. Programming 67 (2), 189–224. doi:10.1007/BF01582221

CrossRef Full Text | Google Scholar

Cox, C., and Munk, W. (1954). Measurement of the Roughness of the Sea Surface from Photographs of the Sun's Glitter. J. Opt. Soc. Am. 44 (11), 838–850. doi:10.1364/josa.44.000838

CrossRef Full Text | Google Scholar

El-Habashi, A., Ahmed, S., Ondrusek, M., and Lovko, V. (2019). Analyses of Satellite Ocean Color Retrievals Show Advantage of Neural Network Approaches and Algorithms that Avoid Deep Blue Bands. J. Appl. Rem. Sens. 13 (2), 1. doi:10.1117/1.JRS.13.024509

CrossRef Full Text | Google Scholar

Estrella, E. H., Gilerson, A., Foster, R., and Groetsch, P. (2021). Spectral Decomposition of Remote Sensing Reflectance Variance Due to the Spatial Variability from Ocean Color and High-Resolution Satellite Sensors. J. Appl. Rem. Sens. 15 (2), 024522. doi:10.1117/1.JRS.15.024522

CrossRef Full Text | Google Scholar

Fan, Y., Li, W., Chen, N., Ahn, J.-H., Park, Y.-J., Schroeder, T., et al. (2021). OC-SMART: A Machine Learning Based Data Analysis Platform for Satellite Ocean Color Sensors. Remote Sensing Environ. 253, 112236. doi:10.1016/j.rse.2020.112236

CrossRef Full Text | Google Scholar

IOCCG (2019). “Uncertainties in Ocean Colour Remote Sensing,” in Reports No. 18 of the International Ocean-Colour Coordinating Group. Editor F. Mélin (Dartmouth, NS: IOCCG). doi:10.25607/OBP-696)

CrossRef Full Text | Google Scholar

Franz, B. A., Bailey, S. W., Eplee, R. E., Lee, S., Patt, F. S., Proctor, C., et al. (2018) NASA Multi-Mission Ocean Color Reprocessing 2018.0.”in Proc. Of Ocean Optics XXIV. Dubrovnik, Croatia.

Google Scholar

Franz, B. A., Bailey, S. W., Werdell, P. J., and McClain, C. R. (2007). Sensor-independent Approach to the Vicarious Calibration of Satellite Ocean Color Radiometry. Appl. Opt. 46, 5068. doi:10.1364/AO.46.005068

PubMed Abstract | CrossRef Full Text | Google Scholar

Franz, B. A., and Karaköylü, E. M. (2016). “Estimating Uncertainty in the Retrieval of Water-Leaving Reflectance from Spaceborne Ocean Color Sensors: Effect of Instrument Noise.” in Proc. of Ocean Opties. Victoria, B.C.. Available at: https://www.researchgate.net/profile/Bryan-Franz-2/publication/309313055_Estimating_uncertainty_in_the_retrieval_of_water-leaving_reflectance_from_spaceborne_ocean_color_sensors_effect_of_instrument_noise/links/5809267408ae993dc0509fe1/Estimating-uncertainty-in-the-retrieval-of-water-leaving-reflectance-from-spaceborne-ocean-color-sensors-effect-of-instrument-noise.pdf

Google Scholar

Frouin, R. J., Franz, B. A., Ibrahim, A., Knobelspiesse, K., Ahmad, Z., Cairns, B., et al. (2019). Atmospheric Correction of Satellite Ocean-Color Imagery during the PACE Era. Front. Earth Sci. 7, 145. doi:10.3389/feart.2019.00145

PubMed Abstract | CrossRef Full Text | Google Scholar

Frouin, R., Schwindling, M., and Deschamps, P.-Y. (1996). Spectral Reflectance of Sea Foam in the Visible and Near-Infrared: In Situ Measurements and Remote Sensing Implications. J. Geophys. Res. 101 (C6), 14361–14371. doi:10.1029/96jc00629

CrossRef Full Text | Google Scholar

Gergely, M., and Zibordi, G. (2014). Assessment of AERONET-OC LWN Uncertainties. Metrologia 51, 40–47. doi:10.1088/0026-1394/51/1/40

CrossRef Full Text | Google Scholar

Gilerson, A., Carrizo, C., Foster, R., and Harmel, T. (2018). Variability of the Reflectance Coefficient of Skylight from the Ocean Surface and its Implications to Ocean Color. Opt. Express 26 (8), 9615–9633. doi:10.1364/oe.26.009615

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilerson, A., Malinowski, M., Herrera, E., Tomlinson, M. C., Stumpf, R. P., and Ondrusek, M. E. (2021). Estimation of Chlorophyll-A Concentration in Complex Coastal Waters from Satellite Imagery. Proc. SPIE 11752ocean Sensing Monit. XIII. doi:10.1117/12.2588004

CrossRef Full Text | Google Scholar

Gordon, H. R., Du, T., and Zhang, T. (1997). Remote Sensing of Ocean Color and Aerosol Properties: Resolving the Issue of Aerosol Absorption. Appl. Opt. 36, 8670. doi:10.1364/AO.36.008670

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, H. R., and Morel, A. Y. (1983). “In-water Algorithms,” in Remote Assessment of Ocean Color for Interpretation of Satellite Visible Imagery: A Review. (Berlin, Germany: Springer-Verlag). doi:10.1007/978-1-4684-6280-710.1007/978-1-4684-6280-7_3

CrossRef Full Text | Google Scholar

Gordon, H. R., and Wang, M. (1994). Retrieval of Water-Leaving Radiance and Aerosol Optical Thickness over the Oceans with SeaWiFS: a Preliminary Algorithm. Appl. Opt. 33, 443–452. doi:10.1364/ao.33.000443

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, H. R., and Wang, M. (1992). Surface-roughness Considerations for Atmospheric Correction of Ocean Color Sensors 1: The Rayleigh-Scattering Component. Appl. Opt. 31 (21), 4247–4260. doi:10.1364/ao.31.004247

PubMed Abstract | CrossRef Full Text | Google Scholar

Groetsch, P. M. M., Foster, R., and Gilerson, A. (2020). Exploring the Limits for Sky and Sun Glint Correction of Hyperspectral Above-Surface Reflectance Observations. Appl. Opt. 59 (9), 2942–2954. doi:10.1364/ao.385853

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrera-Estrella, E., Grotsch, P., Gilerson, A., Malinowski, M., and Ahmed, S. (2020). Blue Band Reflectance Uncertainties in Coastal Waters and Their Impact on Retrieval Algorithms. Proc. SPIE 11420 Ocean Sensing Monit. XII, 1142006. doi:10.1117/12.2559895

CrossRef Full Text | Google Scholar

Hlaing, S., Gilerson, A., Foster, R., Wang, M., Arnone, R., and Ahmed, S. (2014). Radiometric Calibration of Ocean Color Satellite Sensors Using AERONET-OC Data. Opt. Express 22, 23385–23401. doi:10.1364/OE.22.023385

PubMed Abstract | CrossRef Full Text | Google Scholar

Hlaing, S., Harmel, T., Gilerson, A., Foster, R., Weidemann, A., Arnone, R., et al. (2013). Evaluation of the VIIRS Ocean Color Monitoring Performance in Coastal Regions. Remote Sensing Environ. 139, 398–414. doi:10.1016/j.rse.2013.08.013

CrossRef Full Text | Google Scholar

Hu, C., and Carder, K. L. (2002). Atmospheric Correction for Airborne Sensors: Comment on a Scheme Used for CASI. Remote Sens. Environ. 79 (1), 134–137. doi:10.1016/S0034-4257(01)00232-2

CrossRef Full Text | Google Scholar

Hu, C., Feng, L., and Lee, Z. (2013). Uncertainties of SeaWiFS and MODIS Remote Sensing Reflectance: Implications from clear Water Measurements. Remote Sensing Environ. 133, 168–182. doi:10.1016/j.rse.2013.02.012

CrossRef Full Text | Google Scholar

Kokhanovsky, A. A., Budak, V. P., Cornet, C., Duan, M., Emde, C., Katsev, I. L., et al. (2010). Benchmark Results in Vector Atmospheric Radiative Transfer. J. Quantitative Spectrosc. Radiative Transfer 111 (12-13), 1931–1946. doi:10.1016/j.jqsrt.2010.03.005

CrossRef Full Text | Google Scholar

Lamquin, N., Bourg, L., Lerebourg, C., Martin-Lauzer, F. R., Kwiatkowska, E., and Dransfeld, S. (2017). System Vicarious Calibration of Sentinel-3 OLCI. Conference of Characterization and Radiometric Calibration for Remote Sensing. USA: CALCON. Utah.

Google Scholar

Li, J., Jamet, C., Zhu, J., Han, B., Li, T., Yang, A., et al. (2019). Error Budget in the Validation of Radiometric Products Derived from OLCI Around the China Sea from Open Ocean to Coastal Waters Compared with MODIS and VIIRS. Remote Sensing 11 (20), 2400. doi:10.3390/rs11202400

CrossRef Full Text | Google Scholar

Mobley, C. D. (1994). Light and Water: Radiative Transfer in Natural Waters. San Diego, CA: Academic Press.

Google Scholar

Mobley, C. D. (2022). The Oceanic Optics Book. Dartmouth, NS, Canada: IOCCG. doi:10.25607/OBP-1710

CrossRef Full Text | Google Scholar

Mobley, C. D., Werdell, J., Franz, B., Ahmad, Z., and Bailey, S. (2016). Atmospheric Correction for Satellite Ocean Color Radiometry - A Tutorial and Documentation of the Algorithms Used by the NASA Ocean Biology Processing Group. Greenbelt, MD, 73. NASA/TM-2016-217551. doi:10.13140/RG.2.2.23016.78081

CrossRef Full Text | Google Scholar

Moore, T. S., Campbell, J. W., and Feng, H. (2015). Characterizing the Uncertainties in Spectral Remote Sensing Reflectance for SeaWiFS and MODIS-Aqua Based on Global In Situ Matchup Data Sets. Remote Sensing Environ. 159, 14–27. doi:10.1016/j.rse.2014.11.025

CrossRef Full Text | Google Scholar

IOCCG (2010). “Atmospheric Correction for Remotely-Sensed Ocean-Colour Products,” in Reports No. 10 of the International Ocean-Colour Coordinating Group. Editor M. Wang (Dartmouth, NS: IOCCG).

Google Scholar

Oo, M., Vargas, M., Gilerson, A., Gross, B., Moshary, F., and Ahmed, S. (2008). Improving Atmospheric Correction for Highly Productive Coastal Waters Using the Short Wave Infrared Retrieval Algorithm with Water-Leaving Reflectance Constraints at 412 Nm. Appl. Opt. 47 (21), 3846–3859. doi:10.1364/AO.47.003846

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, L., Lee, Z., Hu, C., and Wang, M. (2017). Requirement of Minimal Signal-To-Noise Ratios of Ocean Color Sensors and Uncertainties of Ocean Color Products. J. Geophys. Res. Oceans 122 (3), 2595–2611. doi:10.1002/2016jc012558

CrossRef Full Text | Google Scholar

Ransibrahmanakul, V., and Stumpf, R. P. (2006). Correcting Ocean Colour Reflectance for Absorbing Aerosols. Int. J. Remote Sensing 27 (9), 1759–1774. doi:10.1080/01431160500380604

CrossRef Full Text | Google Scholar

Sayer, A. M., Govaerts, Y., Kolmonen, P., Lipponen, A., Luffarelli, M., Mielonen, T., et al. (2020). A Review and Framework for the Evaluation of Pixel-Level Uncertainty Estimates in Satellite Aerosol Remote Sensing. Atmos. Meas. Tech. 13, 373–404. doi:10.5194/amt-13-373-2020

CrossRef Full Text | Google Scholar

IOCCG (2021). “Observation of Harmful Algal Blooms with Ocean Colour Radiometry,” in Reports No. 20 of the International Ocean-Colour Coordinating Group. Editors S. Bernard, R. Kudela, L. Robertson Lain, and G. C. Pitcher (Dartmouth, NS: IOCCG). doi:10.25607/OBP-1042)

CrossRef Full Text | Google Scholar

Shi, W., and Wang, M. (2007). Detection of Turbid Waters and Absorbing Aerosols for the MODIS Ocean Color Data Processing. Remote Sensing Environ. 110 (2), 149161–161. doi:10.1016/j.rse.2007.02.013

CrossRef Full Text | Google Scholar

Siegel, D. A., Wang, M., Maritorena, S., and Robinson, W. (2000). Atmospheric Correction of Satellite Ocean Color Imagery: the Black Pixel assumption. Appl. Opt. 39 (21), 3582–3591. doi:10.1364/ao.39.003582

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, S. R., Legler, D. M., and Verzone, K. V. (2001). Quantifying Uncertainties in NCEP Reanalyses Using High-Quality Research Vessel Observations. J. Clim. 14 (20), 4062–4072. doi:10.1175/1520-0442(2001)014<4062:quinru>2.0.co;2

CrossRef Full Text | Google Scholar

Steinmetz, F., Deschamps, P.-Y., and Ramon, D. (2011). Atmospheric Correction in Presence of Sun Glint: Application to MERIS. Opt. Express 19 (10), 9783–9800. doi:10.1364/OE.19.009783

PubMed Abstract | CrossRef Full Text | Google Scholar

Teillet, P. M. (1990). Rayleigh Optical Depth Comparisons from Various Sources. Appl. Opt. 29 (13), 1897–1900. doi:10.1364/AO.29.001897

PubMed Abstract | CrossRef Full Text | Google Scholar

IOCCG (2008). “Why Ocean Colour? the Societal Benefits of Ocean-Colour Technology,” in Report No. 7 of the International Ocean-Colour Coordinating Group. Editors T. Platt, N. Hoepffner, V. Stuart, and C. Brown (Dartmouth, NS: IOCCG). doi:10.25607/OBP-97)

CrossRef Full Text | Google Scholar

Tynes, H. H., Kattawar, G. W., Zege, E. P., Katsev, I. L., Prikhach, A. S., and Chaikovskaya, L. I. (2001). Monte Carlo and Multicomponent Approximation Methods for Vector Radiative Transfer by Use of Effective Mueller Matrix Calculations. Appl. Opt. 40 (3), 400–412. doi:10.1364/ao.40.000400

PubMed Abstract | CrossRef Full Text | Google Scholar

Voss, K. J., Gordon, H. R., Flora, S., Johnson, B. C., Yarbrough, M., Feinholz, M., et al. (2017). A Method to Extrapolate the Diffuse Upwelling Radiance Attenuation Coefficient to the Surface as Applied to the Marine Optical Buoy (MOBY). J. Atmos. Ocean Technol. 34 (7), 1423–1432. doi:10.1175/JTECH-D-16-0235.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M. (1999). Atmospheric Correction of Ocean Color Sensors: Computing Atmospheric Diffuse Transmittance. Appl. Opt. 38, 451–455. doi:10.1364/AO.38.000451

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., and Bailey, S. W. (2001). Correction of Sun Glint Contamination on the SeaWiFS Ocean and Atmosphere Products. Appl. Opt. 40 (27), 4790–4798. doi:10.1364/ao.40.004790

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., and Gordon, H. R. (2002). Calibration of Ocean Color Scanners: How Much Error Is Acceptable in the Near Infrared? Rem. Sens 82 (2-3), 497–504. doi:10.1016/s0034-4257(02)00072-x

CrossRef Full Text | Google Scholar

Wei, J., Yu, X., Lee, Z., Wang, M., and Jiang, L. (2020). Improving Low-Quality Satellite Remote Sensing Reflectance at Blue Bands over Coastal and Inland Waters. Remote Sensing Environ. 250, 112029. doi:10.1016/j.rse.2020.112029

CrossRef Full Text | Google Scholar

Werdell, P. J., Bailey, S., Fargion, G., Pietras, C., Knobelspiesse, K., Feldman, G., et al. (2003). Unique Data Repository Facilitates Ocean Color Satellite Validation. Eos Trans. AGU 84 (38), 377–387. doi:10.1029/2003EO380001

CrossRef Full Text | Google Scholar

Werdell, P. J., Bailey, S. W., Franz, B. A., Morel, A., and McClain, C. R. (2007). On-orbit Vicarious Calibration of Ocean Color Sensors Using an Ocean Surface Reflectance Model. Appl. Opt. 46 (23), 5649–5666. doi:10.1364/AO.46.005649

PubMed Abstract | CrossRef Full Text | Google Scholar

Werdell, P. J., Behrenfeld, M. J., Bontempi, P. S., Boss, E., Cairns, B., Davis, G. T., et al. (2019). The Plankton, Aerosol, Cloud, Ocean Ecosystem Mission: Status, Science, Advances. Bull. Am. Meteorol. Soc. 100 (9), 1775–1794. doi:10.1175/BAMS-D-18-0056.1

CrossRef Full Text | Google Scholar

Xiong, X., Angal, A., Chang, T., Chiang, K., Lei, N., Li, Y., et al. (2020). MODIS and VIIRS Calibration and Characterization in Support of Producing Long-Term High-Quality Data Products. Remote Sensing 12 (19), 3167. doi:10.3390/rs12193167

CrossRef Full Text | Google Scholar

Zhang, M., Hu, C., and Barnes, B. B. (2019). Performance of POLYMER Atmospheric Correction of Ocean Color Imagery in the Presence of Absorbing Aerosols. IEEE Trans. Geosci. Remote Sensing 57 (9), 6666–6674. doi:10.1109/TGRS.2019.2907884

CrossRef Full Text | Google Scholar

Zibordi, G., Holben, B. N., Talone, M., D’Alimonte, D., Slutsker, I., Giles, D. M., et al. (2021). Advances in the Ocean Color Component of the Aerosol Robotic Network (AERONET-OC). J. Atmos. Ocean Technol. 38 (4), 725–746. doi:10.1175/JTECH-D-20-0085.1

CrossRef Full Text | Google Scholar

Zibordi, G., Mélin, F., Berthon, J.-F., Holben, B., Slutsker, I., Giles, D., et al. (2009). AERONET-OC: A Network for the Validation of Ocean Color Primary Products. J. Atmos. Ocean Technol. 26 (8), 1634–1651. doi:10.1175/2009jtecho654.1

CrossRef Full Text | Google Scholar