Determining the primary sources of uncertainty in the retrieval of marine remote sensing reflectance from satellite ocean color sensors II. Sentinel 3 OLCI sensors

Uncertainties in remote sensing reflectance R r s for the Ocean Color sensors strongly affect the quality of the retrieval of concentrations of chlorophyll-a and water properties. By comparison of data from SNPP VIIRS and several AERONET-OC stations and MOBY, it was recently shown that the main uncertainties come from the Rayleigh-type spectral component (Gilerson et al., 2022), which was associated with small variability in the Rayleigh optical thickness in the atmosphere and/or its calculation. In addition, water variability spectra proportional to R r s were found to play a significant role in coastal waters, while other components including radiances from aerosols and glint were small. This work expands on the previous study, following a similar procedure and applying the same model for the characterization of uncertainties to the Sentinel-3A and B OLCI sensors. It is shown that the primary sources of uncertainties are the same as for VIIRS, i.e., dominated by the Rayleigh-type component, with the total uncertainties for OLCI sensors typically higher in coastal areas than for VIIRS.


Introduction
Ocean Color (OC) satellite sensors provide radiometric information in multiple bands, allowing to determine spectra of absorption and backscattering of water components and through them, concentrations of chlorophyll-a, colored dissolved organic matter (CDOM) and non-algal particles (Mobley, 2022). These quantities are monitored over global and regional scales and are used in multiple applications including fisheries, ecological models (IOCCG, 2021) and for the detection of algal blooms. Thus, OC is indicative of ocean health and biochemistry and for that reason is listed as an essential climate variable (ECV) (IOCCG, 2008). Atmospheric correction is the critical part in the process of derivation of water parameters, since the water leaving radiance is about 10% of the total radiance measured by the sensor at the top of the atmosphere (TOA) in the blue/green spectral range over typical clear waters, with the other 90% coming from scattering and absorption processes in the atmosphere, including from molecules and aerosol particles, and from reflection from the water surface (Gordon and Morel, 1983;IOCCG, 2010;IOCCG, 2019).
Uncertainties in the water leaving radiance and corresponding remote sensing reflectance R rs after atmospheric correction are assumed to come mostly from the non-ideal modeling of aerosol scattering spectra and reflection of Sun and sky light from the windroughened water surface (Gordon and Wang, 1992;Gordon and Wang, 1994;Frouin et al., 1996;Wang and Bailey, 2001;Ahmad et al., 2010;Frouin et al., 2019) from the presence of absorbing aerosols Ransibrahmanakul and Stumpf, 2006;Shi and Wang, 2007) and are mostly pronounced in the blue bands with small R rs . It is assumed that these effects lead to the negative values of R rs in the blue bands (Frouin et al., 2019) that are often observed in coastal areas. More complex processing in the atmospheric correction Oo et al., 2008), removal of the uncertainty as a power law-like "artefact" with exponent −6 (Ransibrahmanakul and Stumpf, 2006), Neural Network approaches (Hieronymi et al., 2017;Fan et al., 2021), utilization of AC algorithms based on the fitting of Rayleigh spectra (Steinmetz et al., 2011;Zhang et al., 2019), and avoidance of the blue bands in algorithms for retrieval of water parameters (El-Habashi et al., 2019;Gilerson et al., 2021) represent only partial solutions to the problem. One of the goals of future satellite missions like NASA Plankton, Aerosol, Cloud, ocean Ecosystem, PACE (Werdell et al., 2019), and EUMETSAT 3MI (Fougnie et al., 2018) with polarimeters on board is to better characterize aerosol parameters, especially in the blue part of the spectrum. However, with three Visible Infrared Imaging Radiometer Suite (VIIRS) sensors operated by NOAA and two European Commission's Copernicus Programme Ocean and Land Colour Instruments (OLCI) currently in space and with more launches planned, these sensors will provide most of OC information in the next decade and beyond. Therefore, it is important for uncertainties in the products from these sensors to be well characterized.
The estimation of uncertainties can be carried out by comparison of the parameters determined from satellite imagery after atmospheric correction with in-situ values, which include data from the AERONET-OC Network (Zibordi et al., 2009;Zibordi et al., 2021), from buoys like Marine Optical BuoY (MOBY) (Clark et al., 1997) and from ships (Moore et al., 2015). Specifically, in Moore et al., 2015 uncertainties in R rs were estimated from such sources for 7 optical water types (OWT), and it was found that R rs uncertainties are generally highest in the blue part of the spectrum in both clear and coastal waters.
In Zibordi et al., 2022 a comprehensive assessment of OLCI products and related uncertainties was carried out over European water areas through the comparison of satellite and AERONET-OC data and the additional comparison with similar VIIRS products, showing that OLCI typically displays higher uncertainties than VIIRS in coastal waters. A special case of OLCI-B and OLCI-A called "Tandem Phase" was analyzed, with Sentinel-3B and Sentinel-3A flying 30 s apart on the same orbit, demonstrating much smaller uncertainties than with satellites on the same orbit but with a different phase, which is the default operational configuration of the satellites.
In Herrera-Estrella et al., 2021, a first version of the model used in this study was developed to evaluate the spectral composition of R rs uncertainties and was applied to characterize uncertainties due to the R rs spatial distribution in images from the VIIRS sensor on the SNPP platform and the Landsat-8 Operational Land Imager (OLI) at different spatial resolutions. Most of those uncertainties were attributed to the surface effects and water variability conditions.
In Gilerson et al., 2022, the model was expanded and elaborated to determine spectral components of R rs uncertainties, calculated from the comparison of SNPP-VIIRS data and data from MOBY and 8 AERONET-OC stations in the US and Europe. It was shown that the main R rs uncertainties come from Rayleigh-type spectral components (Rayleigh scattering and surface effects) strongly peaked towards the blue, as well as from water variability with spectra proportional to R rs spectra and mostly pronounced in coastal waters. It was hypothesized that Rayleigh uncertainties are associated with a 1%-1.5% inaccuracy in the estimation of the Rayleigh optical thickness (ROT) due to the discrepancy between modeled and measured surface pressure and with the possible variability of the vertical distribution of the gaseous component in the atmosphere. It was also shown that Rayleigh uncertainties at all stations are well characterized by the spectra of the standard deviations of gains measured in the process of the vicarious calibration of the sensor at the MOBY site, meaning that uncertainties are similar at MOBY and other sites and are related to the atmospheric processes and/or their processing in the atmospheric correction model.
In this work, the same model is now applied to estimate spectral components of the uncertainties in R rs retrieval by comparing Sentinel 3A and 3B OLCI satellite data and in situ data from the MOBY site and a similar group of AERONET-OC stations in US and European waters (Table 1), with results compared to corresponding SNPP-VIIRS-AERONET-OC matchups.
In Section 2 the model (Gilerson et al., 2022) is briefly discussed, in Section 3 satellite and in situ data are described and results are presented in Section 4. Conclusions are provided in Section 5.

Main relationships in the estimation of uncertainties
A brief description of the model from Gilerson et al., 2022, is given here for the convenience of the reader. 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 (λ) L w (λ)/E d (λ), where L w (λ) is the water leaving radiance, E d (λ) is the downwelling irradiance and λ is the wavelength. The main relationship for the top of the atmosphere (TOA) total radiance, L t (λ), can be represented as (Gordon and Wang, 1994) where L r (λ) is due to the Rayleigh scattering and surface effects, L a (λ) is due to the aerosol scattering and Rayleigh-aerosol Frontiers in Remote Sensing frontiersin.org 02 interaction, L g (λ) is due to the Sun glint from the water surface, L w (λ) is due to the water leaving radiance and t(λ) is the diffuse transmittance of light from the water surface to the TOA. L t (λ) measured at the satellite sensor at 3 × 3 or more pixels has uncertainties due to all of these components in addition to vicarious calibration and sensor noise. It is important to emphasize that in the recording process L t (λ) already contains uncertainties from all its components in Eq. 1a because of the spatial variability of these components in the pixels of study. These include uncertainties from the water spatial variability, which as was shown by Gilerson et al., 2022, are not negligible in coastal areas. In the process of retrieval of the water leaving radiance, L r (λ), L a (λ), L g (λ) and t(λ) are modeled and radiances are subtracted from the measured L t (λ), which introduces another set of uncertainties between actual and modeled radiances and transmittance coefficients: In the model, L r (λ) was further divided into radiance from the Rayleigh scattering L R (λ) in the atmosphere and reflectance from the ocean surface: where L surf (λ) L sky (λ)*ρ, L sky (λ) is the sky radiance and ρ is the reflectance coefficient from the water surface. While usually in the satellite atmospheric correction procedure averaged surface effects are considered as a part of L r (λ) (Cox and Munk, 1954;Gordon and Wang, 1992), in this model they are considered separately because of differences in the spectra of L R (λ) and L surf (λ). Additionally, 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. Uncertainties from all components included in Eq. 1a, 1b, 1c in the recording of the signal and in the retrieval process need to be taken into account. After normalizing radiances by the downwelling irradiance, E d (λ), the uncertainty in remote sensing reflectance σ in sr −1 can be determined from: Variances for the quantities at TOA σ 2 t , σ 2 R , σ 2 a , σ 2 g are divided by t 2 in accordance with Eq. 1a, 1b, 1c. It was shown (Herrera-Estrella, et al., 2021) that estimated σ noise for the VIIRS sensor (Qi, et al., 2017;Xiong, et al., 2020) is significantly smaller than the total uncertainties σ(λ) (Moore et al., 2015). OLCI σ noise (based on S3 OLCI Cyclic Performance Report, 2021) is of the same order (0.4*10 -4 sr −1 at 412 nm and slightly lower for other bands) and will thus not be considered further. As before for VIIRS, OLCI σ 2 R , σ 2 a , σ 2 g and σ 2 surf contain uncertainties due to both natural variability inside the set of pixels and uncertainties due to modeling inaccuracies, while σ 2 t is at least partially due to the vicarious calibration. It can also include other systematic errors due to detectors, polarization effects, stray light, etc., but these errors are not included in the model.
As one of the main assumptions of the model, all standard deviation components in Eq. 2 except σ t were, as a first approximation, considered proportional to the corresponding mean values of the normalized radiances with proportionality coefficient k: where σ vc (λ) = σ gains (λ)L t (λ)/E d (λ), and σ gains (λ) is the standard deviation of the system vicarious calibration gains (unitless) for OLCI from EUMETSAT processing. Additionally, In Eq. 4a F 0 (λ) is the extraterrestrial irradiance, θ is the sensor zenith angle, Θ is the scattering angle, the angle between solar and viewing directions; in Eq. 4b ω 0 (λ) is the single scattering albedo, P a is the scattering function for aerosols; in (5c) T 0 (λ), T(λ) 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 −1 (Wang and Bailey, 2001); in Eq. 4d θ 0 and t 0 (λ) are the Sun zenith angle and the corresponding diffuse transmittance; in Eq. 4e 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 θ 0 42°at a viewing zenith angle of 40°, with τ a (443) and Angstrom coefficient average values for each specific area based on the satellite processing values given below in Table 2. Aerosol radiance in Eq.4b was determined assuming ω 0 = 1 for all wavelengths, with aerosol optical thickness values and the Angstrom coefficient derived from satellite imagery, and with the phase function (PF) equal to 0.3 corresponding to the scattering angle around 120°. The water component σ water was expressed directly proportional to the remote sensing reflectance R rs .
For each available matchup between satellite and AERONET-OC measurements, all radiance spectra in Eqs 4a, 4b, 4c, 4d 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. 3, and on mean R rs (λ) to determine the contribution of each component to the total σ 2 (λ). Because of the shorter period of operation of OLCI on S3A and especially on S3B than of VIIRS, the number of matchups was significantly smaller than for VIIRS and Frontiers in Remote Sensing frontiersin.org 03 differences in wind speeds were not considered. In addition, outliers in uncertainties spectra or unusual spectra not typical for respective water areas from the AERONET-OC (in-situ) data were removed from matchups since they had a significant impact on results given the small number of total data points.
σ(λ) was calculated from the comparison of satellite data with the corresponding AERONET-OC station data in all bands as with σ(λ) = RMSD. Biases were also calculated as Each of radiances in Eq. 4, has individual spectral variable(s) not repeated in other equations, so all spectra of uncertainties, which were considered proportional to these radiances, should be independent of each other. An optimization procedure was carried out in the same manner as in Gilerson et al., 2022 in 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 σ(λ) components and their individual contribution to the total observed R rs variance σ 2 (λ) as described in Eq. 3. Spectral components were normalized before being used as inputs in the procedure, and results were averaged over 100 runs with random initial conditions to test the robustness of the solution against possible local minima. The normalization did not affect spectral shapes of components and was used to avoid small values of uncertainties, which can be not treated correctly in the fitting process. Once a solution was reached, all normalized k coefficients were then scaled back to their true scale values and were interpreted as an indication of the major contributing components to the total observed R rs variance σ 2 (λ).
Equation 3 can be alternatively presented as where normalized radiance spectra L R /E d , etc. are divided by their values at 412 nm (at 560 nm for R rs ) and corresponding k coefficients are multipled by these values. Thus k coefficients are replaced by the standard deviation for each term at the reference wavelength, 412 and 560 respectively. In this case standard deviations at the reference wavelengths have clear physical interpretation and they can be conveniently compared for different terms.
3 Satellite and AERONET-OC data 3.1 OLCI data OLCI S3A and S3B Collection 3 level 2 Full Resolution imagery (300 m spatial resolution by pixel, EUMETSAT, 2021) was downloaded from the EUMETSAT website for the period of January 2017 to July 2022 for S3A and from May 2018 to July 2022 for S3B for the areas of MOBY in Lanai, Hawaii and seven Aerosol Robotic Network for Ocean Color (AERONET-OC) sites: Casablanca, University of South California (USC), Venise, Gloria, WaveCIS, the Long Island Sound Coastal Observatory (LISCO) and the Helsinki Lighthouse (HLT) with locations shown in Figure 1.
Collection 3 was obtained using the extraction data base (EDB) workflow scripts (https://github.com/juanchossn/ ThoMaS), and Level 2 SeaBASS -Ocean Colour Database (OCDB) format images were retrieved with a window size of 25 × 25 pixels centered on the previously named sites. Each level 2 file includes geophysical products of the atmosphere and ocean, such as aerosol optical thickness and Angstrom exponent at 865 nm, water-leaving reflectance at 412.5, 443.5, 490, 560, and 665 nm among others, sensor zenith angle, solar zenith angle, and quality flags. R rs (λ) is calculated dividing the water-leaving reflectance spectra by π.
OLCI Collection 3 level 2 operational water reflectance products are not corrected for the bidirectional reflectance distribution function (BRDF). The BRDF correction is not applied for historical reasons as the main users have typically been interested in coastal and inland waters where the standard open-ocean BRDF approach does not apply. Nevertheless, the BRDF correction was applied externally in this study to OLCI L2 products as it is necessary for the current analyses (Morel et al., 2002).
Pixels flagged by at least one of the following conditions were excluded: invalid flag, land, cloud (including ambiguous and marginal), coastline, high solar zenith angle larger than 70°, saturated flag, moderate or high glint, whitecaps, and atmospheric correction fail. It should be noted that this set of flags is slightly different from the set recommended for OLCI by EUMETSAT . 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 5 spatial resolutions: 900, 1,500, 2,100, 3,900, and 5,100 m (3 × 3, 5 × 5, 7 × 7, 13 × 13, and 17 × 17 pixel boxes), centered at the AERONET-OC site (Hlaing et al., 2013). Average R rs (λ) and the standard deviation between pixels and geometry were recorded.

VIIRS data
VIIRS's Satellite Level 2-Version 2018.0 imagery was downloaded, for the same areas, 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 VIIRS (with a pixel resolution of 750 m at nadir) include geophysical products of the atmosphere and ocean, such as aerosol optical thickness at 862 nm, remote sensing reflectance, R rs (λ), in the visible wavelengths 410, 443, 486, 551, and 671 nm, and the level 2 quality flags.
Pixels flagged by at least one of the following conditions were excluded: land, cloud, failure in 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°. A file is selected if at least half of the pixels in the set plus one was flag-free. Pixels used for matchup comparison Frontiers in Remote Sensing frontiersin.org were averaged over 3 spatial resolutions: 2,250, 3,750, and 5,250 m (3 × 3, 5 × 5, and 7 × 7 pixel boxes), centered at the AERONET-OC site (Hlaing et al., 2013;Gilerson et al., 2022). Average R rs (λ) and the standard deviation between pixels, geometry, and radiance were recorded.

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 FIGURE 2 (A) Mean R rs (λ) spectra from AERONET-OC, (B) standard deviations of R rs (λ) for matchups with S3A. (C, D, and E) corresponding uncertainties spectra σ(λ) for all areas of study, S3A, 3B and VIIRS respectively, (F, G, and H) biases for S3A, S3B and VIIRS.

Frontiers in Remote Sensing
frontiersin.org 06 radiance and aerosol optical depth data using the SeaPRISM autonomous radiometer systems CE-318/CE-318T (CIMEL Electronique, France) deployed on offshore fixed platforms (Zibordi et al., 2009;Zibordi et al., 2021). SeaPRISM systems are 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 A t (λ), and the sky radiance, L A sky (λ), as a function of solar zenith angle, sensor viewing angle, and azimuth angle relative to the Sun (Zibordi et al., 2009).
The normalized water leaving radiances, L A w , which were measured and then calculated in accordance with AERONET-OC protocols with the BRDF correction based on the open ocean approach (Zibordi et al., 2009;Zibordi et al., 2021), were downloaded for the sites mentioned above from the AERONET-OC website and were converted into remote sensing reflectance at the wavelengths 412, 443, 488, 551(or 560 when available), and 667 nm available on the CE-318 sunphotometers.
Several AERONET-OC stations had 551 nm band for the whole period of OLCI operation, some of them had new sensor heads with 560 nm band for some part of this period. Considering water type variability between sites and at sites themselves, introduction of some algorithm of band shift most likely would have introduced more uncertainties, so available data at AERONET-OC stations at 551 (or 560) nm were used without changes. SNPP VIIRS also has 551 nm band. Tests were conducted by changing the value of R rs (551) by 5% to simulate a band-shift effect, and the results of the spectral decomposition were not noticeably affected.
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 satellite overpass and in situ observation (Zibordi et al., 2009;Zibordi et al., 2021). More detailed information on AERONET-OC sites is listed in Table 1.

The Marine Optical BuoY (MOBY) data
MOBY (Marine Optical BuoY) is an autonomous buoy anchored offshore of Lanai, Hawaii. Its radiometry data is used by the NASA-OBPG and EUMETSAT as part of their ocean color validation and vicarious calibration activities (Clark et al., 1997). 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, data from deployments 266 to 272 (2019)(2020)(2021) was collected for OLCI S3A and S3B, and 249 to 270 (2012-2021) for VIIRS-NPP. Normalized water-leaving radiance calculated from the top sensor and bottom sensor corrected for BRDF is used here because it has the larger available dataset. "Good" and "questionable" data were used to match up with satellite data. Matchups with more than 5 percent differences were excluded from the analysis. MOBY data that matched the bands from satellite sensors was collected with ±2 h of the satellite overpass.
Main atmospheric parameters at the studied sites determined from AC processing and AERONET retrievals are provided in Table 2. Absorbing aerosols were noticeable at several sites with

Frontiers in Remote Sensing
frontiersin.org 07 the average ω 0 values usually above 0.9 and even below 0.9 at LISCO, but spectral dependence was small.

Results: Spectral components of uncertainties
Mean R rs (λ) spectra from the MOBY site and AERONET-OC stations based on in-situ measurements are shown in Figure 2A; waters were clear or moderately clear at MOBY, Casablanca and USC with R rs (443)> R rs (551), moderately eutrophic at Venise, Gloria and WaveCIS with R rs (443)< R rs (551) and more eutrophic at LISCO and HLT with also low R rs (412) R rs (λ). standard deviations are shown in Figure 2B. Corresponding σ(λ) spectra are shown in Figures 2C-E Figure 3 and these spectra were used in the fitting procedure. In accordance with Eq. 3 σ R , σ a , σ 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 stations, while spectra for the total σ(λ), for aerosols and R rs (λ) components are different.
Spectra for surface effects are slightly different for different stations because of the small difference in the aerosol parameters used for their calculations. 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 Rayleigh and aerosol scattering, which makes these spectra distinct from each other.
Results of fitting are presented in Figure 4 for S3A and S3B. Because of the similarity of the total uncertainties' spectra for S3A and S3B, it is not surprising that spectral components are also similar. For most of the stations fitted spectra (dash black lines) match well solid lines of total uncertainties spectra confirming that most likely all relevant components were taken into account. For all stations the main component of σ(λ) is Rayleigh scattering. Glint and surface effects are small, and aerosol contributions are noticeable for LISCO only. With some differences, water variability terms (green solid lines) are quite pronounced, especially at coastal stations, which is usually not considered in the uncertainties budget.
Results are very similar for SNPP VIIRS with relatively stable σ 412 R (0.9 − 1.4)*10 −3 sr −1 as seen in Figure 5. It should be emphasized that in Figures 4, 5 standard deviations are presented, while the fitting process is carried out on variances, so even small differences in the values of spectral components in Figures 4, 5 are amplified in the fitting procedure.
Since the validity of the assumption that the uncertainties spectra are proportional to the spectra of the components themselves can be less accurate for aerosols than for other components, it is possible that the optimization procedure underestimates the contribution of the aerosol component to the total uncertainties. Generally, however, the aerosol uncertainty levels are consistent with their estimations or slightly below them (Wang, 2007;Gao et al., 2022), where they are at the level of 10 -4 -4*10 -4 sr −1 . LISCO station has a pronounced aerosol component for both OLCI S3A and S3B, which is not visible in VIIRS. It can be partially due to the different atmospheric correction processes in OLCI and VIIRS, as well as different number of matchups for these sensors. Some inaccuracies in spectral composition can also be present after optimization because of the similarity of some of the spectra (cf. Rayleigh and surface effects spectra).
For both OLCI and VIIRS sensors the water variability term (solid green line) is very pronounced at the MOBY site. This can be partially due to the similarity of R rs and Rayleigh spectra but is probably also related to the temporal variability of R rs at the site since, as was shown in Herrera-Estrella et al., 2021, spatial variability of R rs is very small in such waters. Also, Brown et al. (2007) estimated uncertainties for upwelling radiance L u (412) for "good" days at 2.4% and 4.7% in general. Assuming similar uncertainties for the downwelling irradiance and R rs (412) = 0.012 sr −1 , related σ 412 can be about 6*10 -4 sr −1 , which at least partially justifies large uncertainty proportional to the R rs at the MOBY site, well visible in Figures 4, 5 as a green solid line. While results have many similarities between VIIRS and OLCI sensors, some differences are present due to higher uncertainties at OLCI, differences in atmospheric correction algorithms, imperfectness of the fitting approach for the separation of components, especially with similar spectra, and the different number of matchups. Thus, for example, at the LISCO site the shape of the surface and aerosol spectra seemed to switch. For S3A/ B, aerosol spectra are close to exponential with σ a (412) ≈ 10 -3 , and the surface spectra are almost flat and near zero. For VIIRS, the exact opposite happens: the surface spectrum is close to exponential with a magnitude of 10 -3 and the aerosol spectrum is flat and near zero. Interestingly, at the Venise site there is a big surface contribution for S3B, but not for S3A. Glint effects are more visible in VIIRS than in OLCI S3A/B. Standard deviations σ 412 and σ 551 retrieved in the fitting procedure for both OLCI sensors and VIIRS for all areas are shown in Figure 6 as a function of spatial resolution, starting from 3 × 3 pixels (900 × 900 m) to 17 × 17 pixels (5100 × 5100 m) for OLCI and from 3 × 3 pixels (2250 × 2250 m) to 7 × 7 pixels (5250 × 5250 m) for VIIRS. Fully in accordance with Figures 4, 5, with exception of Venise for VIIRS, the standard deviations for aerosol, glint and surface effects components are much smaller than for the Rayleigh scattering component, while water variability through σ 551 is pronounced for most of coastal sites. For aerosols, glint, surface effects and water variability σ 412 , σ 551 are in the same range for both OLCI and VIIRS; however, they are different for the Rayleigh component, as will be discussed in further detail below. Some dependence of these standard deviations on the spatial resolution as well as some differences between OLCI S3A and S3B are not systematic and are most likely due to the different number and timing of matchups, which are not properly averaged over the relatively small number of observations. Rayleigh component σ 412 R for OLCI and VIIRS are compared in Figure 7. It can be clearly seen that these standard deviations differ in a small range for VIIRS and in a much larger range for OLCI sensors. At the MOBY site OLCI σ 412 R values are similar to VIIRS σ 412 R values but are instead higher at other stations and consistently so for both OLCI 3A and 3B. It can be also noticed that these OLCI σ 412 R values, with some exceptions, tend to vary in accordance with the latitude of the stations, as can be seen roughly in Figure 1; Table 1, which leads to their possible dependence on the solar zenith angle, mean values for which are provided in Table 3 together with sensor zenith angle. In Zibordi et al., 2022 it was also shown that dependence of uncertainties exists for the sensor zenith angle as well; however, this is not apparent in VIIRS data. Both these effects merit further study. It should be also noticed that comparison of processing of OLCI data using EUMETSAT and NOAA MSL12 algorithms (Mikelsons et al., 2022) showed smaller R rs uncertainties for MSL12, which can assume that larger OLCI uncertainties than for VIIRS shown above are associated with the processing algorithm and not with design features of the instruments.
In Gilerson et al., 2022 it was shown that σ 412 R = 1*10 -3 sr −1 corresponds approximately to the uncertainty in the Rayleigh optical thickness (ROT) of about 1.5%. The uncertainty in ROT was estimated as the ratio of determined Rayleigh uncertainties spectra to the average Rayleigh radiance component at TOA normalized by the downwelling irradiance. This number remains the same for OLCI in case of the MOBY site. Higher uncertainties for coastal sites are most likely associated with other reasons related to the processing algorithms (for example, higher uncertainties in modeling of the surface pressure, which is directly related to the estimation of ROT).
Following the hypothesis in Gilerson et al., 2022 that Rayleigh components of uncertainties are related to the standard deviation of gains σ gains (λ) during the vicarious calibration at the MOBY site, the fitting procedure was repeated with vicarious calibration terms included as it was given in Eq. 3 and σ gains (λ) below it, with results shown in Figure 8 for both OLCI and VIIRS sensors. As discussed in Gilerson et al., 2022, the standard deviation of gains σ gains (λ) is mostly related to the variability of atmospheric parameters and thus should not depend on the sensor. OLCI S3B σ gains (λ) are very close to the ones from VIIRS, while OLCI S3A σ gains (λ) are about 15% higher. Based on the similarity of the total uncertainties' spectra and their spectral components for OLCI S3A and S3B as shown in Figures 2, 4, this difference in σ gains (λ) (see values for S3A and S3B in Section 2) does not have an immediate explanation. Results of fitting for all areas of study for OLCI S3A, S3B and VIIRS with VC term included.
Frontiers in Remote Sensing frontiersin.org 12 While for VIIRS σ vc /t replaced most of Rayleigh components, for the OLCI S3B sensor an additional Rayleigh component remained, which correlated with higher total uncertainties for OLCI than for VIIRS. The reason for additional OLCI uncertainties is not obvious, but it is probably related to some features in the processing Analysis of the S3A ΔR rs distribution for MOBY, Venise and LISCO sites. First column: all ΔR rs , second column: all ΔR rs -mean bias for the station, third column: histogram of all ΔR rsmean bias for the station at 412 nm.
Frontiers in Remote Sensing frontiersin.org 13 algorithms for these sensors. For OLCI S3A, as for VIIRS, σ vc /t replaced most of Rayleigh components shown in Figure 4, with the obvious exception of the Casablanca station where, for reasons yet unclear, the fitted curve was higher than the spectrum of the total uncertainties. Other components remain similar to how they were in Figure 4 for both OLCI sensors.
Melin, 2022 compared satellite and AERONET-OC data for the European sites and for several sensors launched before OLCI and found that for specific stations the spectra of uncertainties from different sensors are close to each other, thus supporting the hypothesis in Gilerson et., 2022 that such uncertainties are mostly due to the variability of atmospheric parameters (or their calculation).
Zibordi et al., 2022 evaluated uncertainties for the European sites over the period of the so-called "Tandem Phase," with Sentinel-3B and Sentinel-3A flying just 30 s apart on the same orbit, and showed spectral median percent differences lower than ±1% in the 412-560 nm interval, much lower than when the sensors were on the same orbit but with the nominal phase difference. This is also very consistent with the hypothesis of a dependence of the uncertainties on the atmospheric parameters (or their calculation), since these should be almost the same for both satellites in the Tandem Phase configuration.
As discussed in Gilerson et al., 2022, uncertainties related to AERONET-OC measurements (Gergely and Zibordi, 2014), are about 4 times smaller than the total uncertainties at all stations and thus play a small role in the total uncertainties budget. With OLCI uncertainties typically higher than for VIIRS, this is also true for the evaluation of OLCI sensors data. However, if it will be possible to make the Rayleigh component contribution to the total uncertainties smaller, remaining uncertainties can be similar to the ones from AERONET-OC and the issue will require more attention.
In Figure 9 the spectra of ΔR rs R rs i sat − R rs i in−situ are shown for MOBY and two AERONET-OC sites in the first column, with bias subtracted in the second column and the histogram of the second column data at 412 nm in the third column. Satellite data are from OLCI S3A. After correction for biases the ΔR rs distribution is clearly almost symmetrical and close to normal for the open ocean MOBY and coastal Venise and LISCO sites, confirming that negative uncertainties values are for the most part not associated with absorbing aerosols.
In Figure 10 σ vc spectra for OLCI S3A/B and VIIRS calculated as σ vc (λ) = σ gains (λ)L t (λ)/E d (λ) at the MOBY site are compared with Rayleigh components determined from the fitting procedure for these sensors also for the MOBY site and shown in Figures 4, 5. For OLCI S3B and VIIRS σ vc spectra are very similar to each other and to Rayleigh components. The differences are larger for S3A, which are most likely due to the inaccuracies of measurements, different number of matchups, etc., considering that all other spectral characterisitics related to R rs and uncertanties are very similar. This again confirms that Rayleigh uncertainties spectra determined in this work are close to σ vc at the MOBY site and represent a significant part of Rayleigh uncertainties at other sites with remaining differences can be probably mitigated with advanced processing algorithms. To remind, as was shown in Gilerson et al., 2022, for VIIRS σ vc were very close to Rayleigh components for all coastal sites as well.

Conclusion
Uncertainties in remote sensing reflectance for S3A and S3B OLCI sensors were evaluated by comparisons with MOBY and AERONET-OC data for several US and European sites. We applied a previously developed model for the decomposition of uncertainties spectra for OC satellite sensors to OLCI sensors and compared the results with the uncertainties from SNPP VIIRS. The uncertainties for OLCI and VIIRS are found to be spectrally similar, but at the coastal sites OLCI uncertainties are about 50% higher.
It is shown that as previously for VIIRS, the main component in R rs uncertainties spectra from OLCI is the Rayleigh component. For coastal water areas, the water variability term proportional to R rs spectra is also pronounced. All other components like aerosol, surface and glint effects are much smaller. Since the spectra of uncertainties from aerosols are not known, it is possible that the contribution of these uncertainties is underestimated. The effect of absorbing aerosols can be more visible in biases, but their spectral dependence is small at the sites considered here.
As in Gilerson et al., 2022 it is assumed that uncertainties in the Rayleigh component are mostly associated with the variability of the atmospheric parameters (ROT) or with their estimation in the current atmospheric correction algorithms (inaccuracies in modeling of the surface pressure), and/or with variability in the vertical distribution of atmospheric gaseous components. This uncertainty corresponds to about 1.5% of ROT at the MOBY site and higher at other stations. The Rayleigh spectral component of the uncertainties was equal to the standard deviation of gains in the vicarious calibration process after conversion to R rs uncertainties and is close to similar terms for other OC sensors. The reasons for some small differences remain unclear because of the differences in the standard deviations of the spectra of gains for OLCI S3A and S3B. Rayleigh-type uncertainties appear to be the main reason for the negative R rs values in coastal waters, where R rs themselves are FIGURE 10 Standard deviations during system vicarious calibration for OLCI and VIIRS σ vc compared with Rayleigh component spectra for these sensors at the MOBY site.
Frontiers in Remote Sensing frontiersin.org typically small. Positive uncertainties in coastal areas and all uncertainties in the open ocean create some inaccuracies in R rs retrievals, but they are not explicitly visible as errors at high R rs . The reason for the higher uncertainties of OLCI, primarily with Rayleigh spectral shape, in comparison with VIIRS should be further studied. They may be partially related to some specifics in the atmospheric correction processing. In a preliminary manner, we associate them here with the latitude of the site and thus with the solar zenith angle, where in previous studies (Zibordi et al., 2022) dependence on the viewing angle was also noticed. Such dependencies do not exist for the SNPP VIIRS sensor.
To mitigate Rayleigh-type uncertainties atmospheric correction schemes should probably include additional constraints (Steinmetz et al., 2011), which can be even different for open ocean and costal water areas.

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