Field data were derived from autonomous above-water radiometric measurements from SeaWiFS¹ Photometer Revision for Incident Surface Measurements (SeaPRISM or simply PRS hereafter) systems operated on off-shore structures part of AERONET-OC (Zibordi et al., 2021). The sites included in this study (see map in Figure 1) are the Gustav Dalen Lighthouse Tower (GDLT, 58.594N, 17.467E), the Helsinki Lighthouse Tower (HLT, 59.949N,24.926E) and the Irbe Lighthouse Tower (IRLT, 57.751N, 21.723E) in the Baltic Sea, the Gloria (GLR, 44.600N, 29.360E) and Galata (GLT, 43.45N, 28.193E) platforms on the western shelf of the Black Sea, the Acqua Alta Oceanographic Tower (AAOT, 45.314N, 12.508E) in the northern Adriatic Sea and the Casablanca Platform (CSP, 40.717N, 1.358E) in the western Mediterranean Sea. In August 2019, the Gloria site was substituted by a system operating on the Section-7 platform (S7, 44.546N, 29.447E); considering their proximity, the data associated with the two sites were aggregated as one time series (GLRS7). The sites in the northern Adriatic and Black Sea are representative of coastal areas with moderately turbid conditions (see average R_RS spectra in Figure 2) while the Baltic sites are decidedly more influenced by higher levels of Chromophoric Dissolved Organic Matter (CDOM) and show lower R_RS, particularly in the blue (Zibordi et al., 2009, 2021). The CSP location, associated with lower turbidity and conditions often more typical of Case-1 waters (i.e., with optical properties defined by phytoplankton and derivatives), show R_RS generally decreasing from the blue to low values in the red.

This study relied on AERONET-OC (version 3) data of multi-spectral normalized water-leaving radiance corrected for bidirectional effects using an IOP-based approach where the IOPs were computed using regional empirical algorithms (Zibordi et al., 2004, 2009). These data products were then converted in remote-sensing reflectance R_RS by division by extra-terrestrial solar irradiance (Thuillier et al., 2003). Only level-2 R_RS data having undergone the complete AERONET-OC quality control were used in this work. Relative uncertainties obtained for R_RS at a site such as AAOT are ~5% at all bands except the red (~7%) (Zibordi et al., 2004) and tend to increase when R_RS is in a lower range of values such as in Baltic waters (Gergely and Zibordi, 2014). Uncertainty values used in this work for PRS data are those expressed in radiometric units as defined by Gergely and Zibordi (2014).

This study relied on the VIIRS data collected by the mission Suomi National Polar-orbiting Partnership (SNPP) launched in October 2011, and Joint Polar Satellite System (JPSS) -1 launched in November 2017 (Cao et al., 2013; Goldberg et al., 2013), which allows a current overlap of 3.5 years. Satellite Level-1A data for both VIIRS missions were acquired from the Ocean Biology Distributed Active Archive Center (OB.DAAC) of the National Aeronautics and Space Administration (NASA) and processed by the SeaWiFS Data Analysis System (SeaDAS, Fu et al., 1998) standard atmospheric correction (Gordon and Wang, 1994; Franz et al., 2007; Ahmad et al., 2010) according the latest NASA reprocessing option R2018. Besides being processed with the same code, both missions are also subject to the same strategy for system vicarious calibration (Franz et al., 2007), conditions that favor agreement between the products from the two missions. Nominal center-wavelengths for VIIRS on SNPP are given at 410, 443, 486, 551, and 671 nm, while they are 411, 445, 489, 556, and 667 nm for JPSS1.Thereafter, VIIRS-SNPP and VIIRS-JPSS1 are merely referred to as SNPP and JPSS1, respectively.

Macro-pixels of 3x3 pixels centered on each AERONET-OC site were then extracted from the Level-2 imagery. For each macro-pixel, field data registered less than 2 h from the satellite overpass, if any, were first selected. If field measurement records were available before and after the satellite overpass, then a weighted average was computed for the time of overpass using the two closest measurements found before and after the overpass; otherwise the closest measurement was adopted as the PRS field value. The maximum time difference allowed for the analysis (2 h) is in the continuity of past studies at the AERONET-OC sites (e.g., Zibordi et al., 2009). Considering their coastal, possibly more dynamic, character, it is lower than those adopted in studies more focused on the open ocean (e.g., 3 h, Bailey and Werdell, 2006).

Macro-pixels were excluded from the analysis if one pixel was affected by the standard level-2 flags of the SeaDAS processor applied when building level-3 products², to the exception of the flag signaling coccolithophores, the presence of which can be recurrent in the Black Sea but that does not appear to degrade the performance of the atmospheric correction (Cazzaniga et al., 2021). Standard flags excluded all conditions leading to a failure or non-application of the atmospheric correction (such as cloudy conditions, high glint, influence of stray light) and large zenith angles for the illumination and observation directions. Furthermore, match-ups were ignored if the coefficient of variation (CV, ratio of standard deviation and average over the macro-pixel) calculated on R_RS at selected wavelengths (486/489 and 551/556 nm) is larger than 0.2 to avoid spatially heterogeneous conditions. When comparing data from both VIIRS missions, the protocol selecting valid macro-pixels was similar, keeping inter-mission match-ups when the satellite overpass did not exceed 2 h (in practice all inter-mission match-ups were separated by less than 1 h).

Besides slightly different acquisition times, the comparison between satellite and field data is beset by a mismatch in spatial scales (pixel-size vs. point-based), which leads to a representation error (Oke and Sakov, 2008). As said above, heterogeneous conditions at the scale of the 3x3 macro-pixels were excluded with a test on CV. To further reduce the representation error, the satellite value compared with field data was computed at the position of the measurement site with a bilinear interpolation of the values associated with the four closest pixels surrounding its location, and the associated local variability was quantified by the standard deviation found for the four pixel values. Using Landsat-8 high-resolution data, Pahlevan et al. (2016) recommended using the closest pixel, while recognizing that other factors (striping or noise) could then adversely affect the comparison, so that the adopted approach appears as a good compromise. In practice, validation statistics are barely affected when using a 3x3-pixel average or the 2x2-pixel interpolation.

As far as a satellite-to-satellite comparison is concerned, Pahlevan et al. (2016) recommended a 7x7-km window size with size actually varying according to the selected area, but their decision criterion was fairly strict. Moreover the window size might be diversely translated in number of pixels for VIIRS (pixel size at nadir of ~0.75-km) with respect to a sensor such as the Moderate Resolution Spectroradiometer (MODIS, ~1-km). Considering also that the studied sites are in coastal waters, selecting a 3x3-pixel window rather than a larger one seemed safer to avoid potential adjacency effects (Bulgarelli and Zibordi, 2018).

Comparison between field and VIIRS data or between the two VIIRS missions was conducted at corresponding center-wavelengths but these may somewhat differ. Differences in center-wavelengths were corrected by a band-shifting scheme based on regional bio-optical relationships following Zibordi et al. (2009). In practice field data were expressed at the closest VIIRS center-wavelengths. The inter-mission comparison was made with three options, expressing R_RS from SNPP at the nominal center-wavelengths of JPSS1, expressing R_RS from JPSS1 at the nominal center-wavelengths of SNPP, and without performing any band shifting. It is anticipated that the choice of option barely affected comparison statistics.

Considering a set of N match-ups with (_{xi)i = 1, N} and (_{yi)i = 1, N}, the PRS field data and satellite products, respectively, the following comparison metrics were introduced (Mélin and Franz, 2014; IOCCG, 2019):

with the overline indicating an average value. Δ, the root-mean-square (RMS) difference between x and y, can be decomposed by a development of squares into the average difference δ and a centered RMS difference Δ_c that is free from systematic effects. Δ and δ were also computed for the differences between SNPP and JPSS1 and are noted Δ_VRS and δ_VRS (JPSS1 minus SNPP).

Relative differences were expressed with the following definitions:

for the median absolute relative difference |ψ|_m and the median relative difference ψ_m, where the choice of the “median” operator was made to avoid the impact of outliers when the denominator was nearing 0.

Similar equations were used to compute differences between the two satellite products (SNPP and JPSS1). When applied to satellite data, relative differences were written in their unbiased (also called symmetric) form (taking the average of x and y as a reference, Armstrong, 1985; Mélin and Franz, 2014):

Assuming a general approximate linear relationship between the distributions of the field data x and their satellite equivalent y (as supported by validation results, e.g., Bailey and Werdell, 2006; Zibordi et al., 2009; Mélin et al., 2011; Moore et al., 2015), a linear error model was adopted:

The field data x were written as the sum of a target reference state t and a zero-mean random error term ξ (random being here understood as uncorrelated with other quantities). The relationship between the reference t and the true value is undefined with the proviso that non-systematic effects are captured by ξ; in practice t only served as a link between x and y. Similarly, the satellite data y were written as a function of t with additive and multiplicative biases, α and β, respectively, and a zero-mean random error ϵ. It can be noted that the bias δ introduced by Equation (2) is equal to (β-1)t¯+α.

In this framework, the standard deviation of (_{ξi)i = 1, N}, σ_ξ, was assumed equal to the standard uncertainty defined for the field data (and therefore known, see section 2.1), and the objective was to define the terms characterizing the uncertainty budget of the satellite data y, particularly the term related to random effects quantified by the standard deviation of (ϵ_i)_{i = 1, N}, σ_ϵ.

Writing the variance and covariance terms σx2, σy2 and σ_xy from Equations (8) and (9) leads to:

taking advantage of the fact that covariance terms with ξ or ϵ are equal to 0. A similar framework was applied by Mélin et al. (2016) to the case of two satellite data sets, where the system of equations was solved by assuming known the ratio between σ_ξ and σ_ϵ. Here instead, the system can be solved to calculate σ_ϵ knowing σ_ξ, which leads to:

From this equation, σ_ϵ is a decreasing function of σ_ξ, having as maximum σy2-σxy2σx2, i.e., σy2(1-r2) where r is the Pearson correlation coefficient between x and y. Equation (13) is only valid with a data set having a significant range of variability (σ_x>σ_ξ).

For completeness, the value of the slope β is also given (model II regression, Legendre and Legendre, 1998):

noting γ=σ_ϵ/σ_ξ. If γ is very large (in situ data considered error-free and σ_ξ=0), β is the slope of an ordinary least-square regression σ_xy/σx2. Differently if γ=1 (field and satellite data characterized by the same level of non-systematic effects in their uncertainty), β is the slope of a major axis regression (Legendre and Legendre, 1998).

The centered RMS difference Δ_c between field and satellite data can also be linked to the terms β, σ_ϵ, and σ_ξ (by considering that Δc2 is the variance of (y_i − x_i)_{i = 1, N} and therefore a combination of the three terms defined by Equations 10 to 12):

It can be noticed that Δc2 is close to σξ2+σϵ2, i.e., the sum, in variance space, of the uncertainty terms associated with random effects, only if β is close to 1. In the same way, putting aside systematic differences (biases) and assuming field data as error-free (σ_ξ=0), then the difference between field and satellite data Δ_c is equal to ((β-1)2σx2+σϵ2)1/2, and can be equated to the uncertainties associated with satellite data (as represented by σ_ϵ) only if β is close to 1.

The comparison between VIIRS and PRS field data is shown for each site in Figure 3 for selected wavelengths while comparison statistics are given for the entire spectrum in Figures 4A,B for |ψ|_m and ψ_m, and in Figures 5, 6 for Δ, Δ_c, and δ. More detailed statistics are given in Supplementary Material. The number of match-ups varies extensively from 684 at AAOT to 58 at IRLT for SNPP, and from 245 at AAOT to 32 at HLT for JPSS1. Besides the fact that SNPP has been in orbit for a longer period than JPSS1, the number of match-ups obviously depends on the environmental conditions encountered at each site as well as on the periods of operation of the PRS systems. For instance there are no match-ups for the Baltic sites in the winter period and the site of HLT was not active in 2018 and 2020, explaining the low number of match-ups for JPSS1.

Validation results appear well consistent for both VIIRS missions (as well as with previous missions; e.g., Zibordi et al., 2009, 2015; Mélin et al., 2012), with match-up points fairly well-distributed along the 1:1 line when all sites are considered (Figure 3). The maximum values of R_RS are associated with conditions of coccolithophore blooms in the Black Sea (Cazzaniga et al., 2021). The median absolute relative difference |ψ|_m show the expected U shape for this type of metrics (e.g., Mélin and Franz, 2014) with larger values for the blue or red bands when the R_RS values are low (Figures 4A,B). This is obviously the case for Baltic waters that are highly absorbing with low R_RS in the blue (see Figure 2) and |ψ|_m above 60% at 410/411 nm for both SNPP and JPSS1. To the contrary, the CSP site with fairly clear waters is characterized by low R_RS in the red (671/667 nm) and the highest |ψ|_m in that spectral region (44% and 35% for SNPP and JPSS1, respectively). For the two Black Sea sites (GLRS7 and GLT), |ψ|_m is approximately 30 and 20% for the first two bands while it is lower at AAOT and CSP (~21% and 12–14%). For SNPP at 551 nm, |ψ|_m varies in the interval 7–10% for all sites (except 12% at CSP), and for JPSS1 at 556 nm in the interval 7–11% (except 4% at HLT). At 486/489 nm, |ψ|_m is also lower than 10% for AAOT, CSP and the Black Sea sites. In the red, |ψ|_m is in the interval 9–15% for the Baltic sites (except 20% at GDLT for SNPP). Considering the similar properties and validation results and for ease of presentation, the match-ups of the three Baltic Sea sites are grouped to compute uncertainty statistics; the same procedure is carried out with the two Black Sea sites. This also reinforces the statistical significance of the results by increasing the sample size.

When expressed in radiometric units, statistics tend to be closer among the sites (Figures 5, 6), with the RMS difference Δ decreasing regularly with wavelength from an interval among all sites of 0.76–1.7 10⁻³ sr⁻¹ at 410/411 nm to 0.14–0.41 10⁻³ sr⁻¹ at 671/667 nm. The AAOT and CSP sites are characterized by the highest Δ in the blue, while the Baltic sites have the lowest values. When looking at the bias δ, a clear agreement between sites and missions can be noted: δ tends to be negative at 410/411 nm (with a value as large as −0.001 sr⁻¹ for JPSS1 at CSP) and positive at 443/445 nm, illustrating a certain spectral inconsistency for these two bands. Above 486 nm, δ is in the interval ±0.3 10⁻³ sr⁻¹ (and most often ±0.2 10⁻³ sr⁻¹). The centered RMS difference Δ_c is mostly very close to Δ except for the cases where δ is largely negative, at CSP in the blue for both VIIRS missions, and at AAOT for JPSS1.

As far as the multiplicative bias is concerned, β is mostly in the interval 0.8–1.1, except lower values at 410 nm for SNPP for the Baltic sites (0.35) and in the red at CSP (~0.6). This indicates that the first term in the expression of Δ_c [proportional to (β−1)², Equation 15] is usually fairly small.

Using the period of mission overlap, Figure 7 shows the direct comparison between SNPP and JPSS1 R_RS that confirms the fair agreement between the two missions anticipated from the validation statistics. As mentioned in section 2.2, the comparison between the two missions was made with three cases: without band-shifting (keeping R_RS unchanged) and after application of a band-shifting expressing SNPP R_RS at the JPSS1 bands, or JPSS1 R_RS at the SNPP bands. Comparison statistics change very little with these three cases. For instance, Figures 5, 6 show the RMS difference Δ_VRS as an envelop defined by the maximum and minimum Δ_VRS obtained for the three cases at each wavelength: it is barely distinguishable from a single thick line. From now onward, only the case without band-shifting is discussed. It is however stressed that the application of a band-shifting scheme might be required for other pairs of sensors and/or other water types.

Figures 4C,D provides the median (unbiased) relative differences |ψ_u|_m and ψ_{u, m} for the various sites. Again, a U shape is seen for |ψ_u|_m with high values in the blue or red bands when R_RS is low: at 410/411 nm, |ψ_u|_m is higher than 30% for the Baltic sites (up to 48% at HLT), and 23% at 671/667 nm for CSP. To the contrary, the lowest |ψ_u|_m in the red is found at IRLT (7%), GLR (8%), and HLT (9%), and at 410/411 nm at CSP (13%). The agreement is encouraging at 486/489 and 551/556 nm, with |ψ_u|_m always lower than 10% (and often lower than 7%). The relative bias (JPSS1 values minus SNPP values) tends to increase with wavelength, from usually between −4 and −10% at 410/411 and 443/445 nm to positive values in the red. When expressed in radiometric units, the mean difference δ_VRS is also usually negative in the blue (i.e., JPSS1 R_RS lower) and lower (in modulus) than 0.2 10⁻³ sr⁻¹, except at 443/445 nm where δ_VRS can be as low as −0.5 10⁻³ sr⁻¹ (Figures 5, 6).

The RMS difference Δ_VRS between SNPP and JPSS1 R_RS, illustrated in Figures 5, 6, shows a regular decrease with wavelength, fairly parallel to Δ. At AAOT, CSP and the Black Sea sites, Δ_VRS is 0.9–1 10⁻³ sr⁻¹, and is slightly lower at the Baltic site (0.76 10⁻³ sr⁻¹). In the red, Δ_VRS is increasing from CSP (0.06 10⁻³ sr⁻¹), the Baltic sites (0.17 10⁻³ sr⁻¹), the Black Sea sites (0.22 10⁻³ sr⁻¹) to AAOT (0.29 10⁻³ sr⁻¹). In general, Δ_VRS is clearly lower than Δ (RMS difference with respect to field data) except at the Baltic sites when Δ is computed for JPSS1 (Figure 6D).

Following the framework introduced in section 2.4, the uncertainty term associated with random effects was computed with Equation (13) and is here noted σ_VRS for the VIIRS products (or σ_SNPP and σ_JPSS1 when specific to a sensor) and σ_PRS for the PRS field data (equivalent to σ_ϵ and σ_ξ in section 2.4).

As seen in Figures 5, 6 and Supplementary Material, σ_VRS decreases with wavelengths for both VIIRS missions and all regions. At 410/411 nm, it is found in the interval 1.2–1.4 10⁻³ sr⁻¹ at 410/411 nm for the northern Adriatic and western Mediterranean areas (AAOT and CSP), slightly lower for the Black Sea sites (~1.1 10⁻³ sr⁻¹) and lowest for the Baltic sites (0.8–0.9 10⁻³ sr⁻¹). In the green, AAOT and the Black Sea sites have the highest values of σ_VRS (0.5–0.8 10⁻³ sr⁻¹), while it is lower for CSP and the Baltic sites (~0.3 10⁻³ sr⁻¹). In the red, it is highest for AAOT (0.2–0.3 10⁻³ sr⁻¹), lower for Baltic and Black Sea sites (0.15–0.25 10⁻³ sr⁻¹) and lowest (<0.1 10⁻³ sr⁻¹) at CSP, a site where R_RS is low. An important result is that σ_SNPP and σ_JPSS1 are close to each other. Considering the four regions and the five bands, the ratio σ_JPSS1/σ_SNPP is mostly in the interval 0.8–1.05 (except in the green and red bands at AAOT, where it is ~0.7), with σ_JPSS1 lower in almost all cases. Besides possible differences between the two missions, differences in σ_VRS can be explained by the differing number of match-ups found for each mission. On the other hand, the general agreement (compare Figures 5, 6) is favored by the similar characteristics of the sensors and by the use of common algorithms to process their data.

The spectra of σ_VRS appear fairly close to Δ_c, which, considering Equation (15), is explained by the relatively low values of σ_PRS and values of β not far from 1. σ_PRS is always lower than 0.3 10⁻³ sr⁻¹, a value reached in the green at AAOT and in the blue for CSP (in association with relatively high values of R_RS). Considering that σ_VRS is clearly higher than σ_PRS and that the uncertainties associated with PRS data are at least 5% when expressed in relative terms (Gergely and Zibordi, 2014), the VIIRS R_RS products at the considered sites do not comply with an objective of a 5% uncertainty advocated by the Global Climate Observing System (GCOS, 2011). It is however fair to recall that the AERONET-OC sites are located in coastal regions whereas this 5% objective applies to open ocean waters, and that even in these conditions, it might not be currently fulfilled for all bands (Hu et al., 2013; Mélin et al., 2016).

There are no readily available independent estimates of uncertainty for R_RS that could serve as a point of comparison with the present results. Hu et al. (2013) provide uncertainties for SeaWiFS and MODIS (Moderate Resolution Imaging Spectroradiometer) R_RS data by looking at the deviations of R_RS with respect to reference values defined in relation to a reference chlorophyll-a algorithm (Hu et al., 2012). This approach is only applicable in oligotrophic waters, conditions which are usually not observed at the AERONET-OC sites. However, considering that the CSP site is often associated with Case-1 waters, a comparison can be performed for completeness. For the highest chlorophyll-a concentration tabulated by Hu et al. (2013) (0.2 mg m⁻³ in the North Atlantic), uncertainty estimates in the blue (~412 nm) are much lower than σ_VRS at CSP: 0.68 10⁻³ sr⁻¹ and 0.79 10⁻³ sr⁻¹ for MODIS and SeaWiFS, respectively, vs. ~1.3 10⁻³ sr⁻¹ for σ_VRS; on the other hand they are of the same order at the green band, ~0.3 10⁻³ sr⁻¹.