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

The accuracy of remote-sensing reflectance (Rrs) estimated from ocean color imagery through the atmospheric correction step is essential in conducting quantitative estimates of the inherent optical properties and biogeochemical parameters of seawater. Therefore, finding the main source of error is the first step toward improving the accuracy of Rrs. However, the classic validation exercises provide only the total error of the retrieved Rrs. They do not reveal the error sources. Moreover, how to effectively improve this satellite algorithm remains unknown. To better understand and improve various aspects of the satellite atmospheric correction algorithm, the error budget in the validation is required. Here, to find the primary error source from the OLCI Rrs, we evaluated the OLCI Rrs product with in-situ data around the China Sea from open ocean to coastal waters and compared them with the MODIS-AQUA and VIIRS products. The results show that the performances of OLCI are comparable to those MODIS-AQUA. The average percentage difference (APD) in Rrs is lowest at 490 nm (18%), and highest at 754 nm (79%). A more detailed analysis reveals that open ocean and coastal waters show opposite results: compared to coastal waters the satellite Rrs in open seas are higher than the in-situ measured values. An error budget for the three satellite-derived Rrs products is presented, showing that the primary error source in the China Sea was the aerosol estimation and the error on the Rayleigh-corrected radiance for OLCI, as well as for MODIS and VIIRS. This work suggests that to improve the accuracy of Sentinel-3A in the coastal waters of China, the accuracy of aerosol estimation in atmospheric correction must be improved.


Introduction
The first global chlorophyll-a concentration image from the NASA Coastal Zone Color Scanner (CZCS, 1978(CZCS, -1986) launched efforts to begin conducting optical remote sensing of the world's oceans.Many ocean color sensors have been launched since CZCS, including the NASA Sea-viewing Wide Field-of-view Sensor (SeaWiFS, 1997(SeaWiFS, -2010)), the NASA Moderate Resolution Imaging Spectroradiometer (MODIS-T, 1999-present on board the Terra platform, and MODIS-A, 2002-present on board the Aqua platform), the ESA Medium Resolution Imaging Spectrometer (MERIS, 2001(MERIS, -2012)), NOAA Visible Infrared Imaging Radiometer Suite (VIIRS, 2011-present on board the Suomi NPP) and the most recent The latest generation of ocean color remote sensor is the Ocean and Land Color Imager (OLCI), which is on board the Sentinel-3A (S3A) satellite, which was successfully launched by the European Space Agency (ESA) on February 16th, 2016.This satellite features many new capabilities, such as (1) global coverage at 300 m (Full Resolution, FR) or 1200 m (Reduced Resolution, RR) resolution; (2) 21 spectral bands from 400 nm to 1020 nm; (3) sun glint minimization by tilting the sensor from the nadir.OLCI will quickly become the main remote sensor for studying the open ocean and coastal waters [17,18] because four successive versions of the sensor are planned for launch over the next 10-15 years (OLCI-B, on board the Sentinel-3B has been in space since 25 April 2018).
Since OLCI was launched, only a few works have been performed to validate its products [19][20][21][22].Shen et al. (2017) [19] developed a dual band ratio algorithm to calculate the downwelling diffuse attenuation coefficient at 490 nm (Kd(490)) for the waters of Lake Taihu with OLCI data and showed that the new OLCI product has a smoother spatial distribution and finer textural characteristics than does the MODIS product and it contained notably higher-quality data.Zibordi et al. (2018) [20] summarized a regional assessment of radiometric data products from OLCI with in-situ data from the ocean color component of the Aerosol Robotic Network (AERONET-OC) and the bio-optical mapping of marine properties (BiOMaP) program and revealed that there R rs was systematically underestimated while the aerosol optical thickness was overestimated, explainable by biases in calibration coefficients or poor performance of bright pixel correction.Mograne et al. (2019) [21] validated the OLCI water-leaving reflectance products over two contrasted French coastal waters obtained by five different atmospheric correction algorithms (AC), and discovered that the polymer and C2R-CCAltNets algorithms obtain high performances.Gossn et al. (2019) [22] developed a new atmospheric correction algorithm (BLR-AC) for turbid waters based on the red, near-infrared (NIR) and 1016 nm bands of OLCI.They presented a comparison with the NASA/SeaDAS and ESA standard atmospheric correction algorithms, showing that the BLR-AC is better than the NASA/SeaDAS and ESA AC, particularly over extremely turbid waters.However, these published OLCI validation studies described only the overall error of OLCI in the region of interest; they did not analyze the error sources.Therefore, there is a need to validate OLCI products in other regions of interest (in our case, the China Sea) and, more importantly, to further analyze the sources of their overall error.
Here, we present a validation of the OLCI-R rs product and compare it with the MODIS and VIIRS-derived R rs (λ) products around the China Sea based on in situ measurements and match-up analysis.Then, we further assess the uncertainties, the error budget, and the factors that influence the accuracy of the R rs product.

Description of Atmospheric Correction
The goal of atmospheric correction is to estimate and then remove the atmospheric path radiance contribution.The water-leaving radiance is at most 10-20% of the total top-of-atmosphere (TOA) radiance in the visible bands (VIS) over open ocean waters, and it can reach 50% in the red bands over turbid waters.Therefore, atmospheric correction is a critical step for remotely sensed data [11].For the ocean-atmosphere system, after pre-correcting for gas absorption, whitecaps, and sun glitter on the sea surface, the radiance (L t (λ)) measured by the remote sensor can be decomposed as follows [11,23,24] where λ is the wavelength.L r (λ) is the radiance due to Rayleigh scattering, L a (λ) + L ra (λ) is the contribution of the aerosols to scattering and the scattering between aerosols and air molecules, t 0 (λ) is the diffuse transmittance of the atmosphere from the surface to the sensor, θ 0 is the viewing direction, F 0 is the solar irradiance at the mean Earth-Sun distance, and R rs (λ) is the remote-sensing reflectance.
Over open ocean waters, the L A (λ) estimation is based on the hypothesis of a black ocean in the near-infrared (NIR) bands [23].Over turbid waters, this hypothesis is no longer valid; thus, the contribution of the ocean to the NIR must be estimated to accurately estimate R rs in the visible (VIS) bands [11,12,33].
In this study, we focused on the standard OLCI radiometric product obtained with the standard OLCI atmospheric correction algorithm and the MODIS-AQUA and VIIRS radiometric products obtained with the NASA standard atmospheric correction algorithm [24].The latter algorithm has been extensively validated worldwide [13,14,34,35].The hypotheses for the standard OLCI atmospheric correction [36,37] are similar to those of NASA MODIS/VIIRS atmospheric correction [23,24,38,39].
The standard OLCI atmospheric correction algorithm first retrieves the suspended particulate matter (SPM) by assuming the water reflectance (ρ w (λ)) and the single scattering aerosol reflectance (ρ as (λ)) in the NIR bands based on the black ocean hypothesis.Then, the initial SPM is used to get an initial estimate of ρ w (λ) and ρ as (λ) in the NIR bands via an empirical relationship [37].Finally, the ρ as (λ) in the NIR bands is used to get ρ w (λ) in the VIS bands by accounting for the multiple scattering of air molecules and aerosols [36] using Equation (3).
To compare the OLCI product with the MODIS-AQUA and VIIRS R rs , we first present a brief description of the standard NASA atmospheric correction algorithm.This algorithm is described in Bailey et al. (2010) [24].First, the black-pixel assumption is used for both NIR bands to obtain the first initial R rs (λ) estimation [23].Second, Equation ( 3) is used to obtain the initial estimate of R rs (λ) in the NIR bands.Third, the initial R rs (λ) is used to obtain an initial estimate of the chlorophyll concentration by using a bio-optical model.Fourth, this chlorophyll concentration is used to obtain the absorption (a(λ)) and backscatter coefficients (b b (λ)) at the NIR bands via an empirical relationship [24,38].Fifth, a(λ) and b b (λ) in the NIR bands are then used to obtain R rs (λ) in the NIR bands and these quantities are used to remove the non-zero R rs (NIR) contribution to L rc (λ) from the NIR bands.Finally, this process is repeated until R rs (NIR) convergence to obtain the R rs (λ) in all bands.
To summarize, the general flow of these algorithms is as follows: Step 1. Estimate R rs (λ) (or equivalently ρ w (λ)) at the NIR reference wavelengths using the iterative model [24,38,40] so that the non-zero water-leaving radiance can be removed from the TOA signal, leaving only the aerosol reflectance as the contribution to L A (λ).
Step 2. The aerosol reflectance at the NIR bands is used to estimate the aerosol properties and extrapolate aerosol reflectance to the VIS bands.Then the aerosol reflectance can be removed from the TOA signal, leaving the R rs (λ) (or equivalently ρ w (λ)) in the VIS bands.
In the iterative model, the difference between the real R rs (NIR) and the estimated R rs (NIR) is called the error of the iterative model, and it stems from the error that occurs when estimating L A (NIR) and is later transferred to L A (λ) in the VIS bands (L A (VIS)), leading to inaccurate R rs (λ) in the VIS bands (R rs (VIS)).
The aerosol lookup tables (LUTs) used in the atmospheric correction are obtained from simulations of the radiative transfer and take a given number of aerosol models into account [23,33,36].These estimated aerosol LUTs are ideal, but they may differ from the values observed over the ocean, especially over coastal waters.The difference between the real L A (VIS) and the estimated L A (VIS), which we term the error of the aerosols LUTs, also has an impact on the final R rs (VIS) estimation.
Thus, Step 1 of the algorithms depends on the accuracy of the iterative model used to estimate the aerosol properties in the NIR bands, while Step 2 depends on the accuracy of the aerosol LUTs, which are used to extrapolate the aerosol reflectance to the VIS bands.Note that the iterative model and the aerosol LUTs algorithms differ between ESA and NASA, although their general principles are similar.

Error Budget
To improve the product accuracy for ocean color satellites, it is necessary to not only determine the total error but also to find the most algorithm component that contributes the most to the total error.The error budget decomposes the total error into sub-error categories and then compares the sub-errors to obtain the maximal error contribution.Therefore, calculating the error budget is an essential part of the validation.
If the measurement result is determined by a functional relation from other quantities, rather than directly measured, that function or model should represent not only the physical/bio-optical laws but also the measurement process.In particular, it should include all the quantities that may have a significant impact on the uncertainty of the measurement result [41].Here, we introduce an error budget model that represents both the physical/bio-optical laws and the measurement processes and apply it to estimate the error budget of the R rs -derived product [42].
If the result (C) includes two independent parts, such as A and B. The error of C is If the function is given as q = x + . . .+ z − (u + . . .+ w) then the error budget is δq ≈ δx + . . .
where q is the result; x, z, u, and w are function quantities; δq is the error of the result; and δx, δz, δu, and δw are the error of quantities [42].
If the function is given as with B, being a constant.Then the error budget is where q is the result, x is the quantity, δq is the error of the result, and δx is the error of the quantity [42].

Total Error of the Satellite Product
With the in-situ data and match-up procedures, the total difference in R rs (λ) between the in-situ data and the satellite product can be obtained in the VIS bands, for example, the absolute percentage difference (APD, Equation (25)) and the bias (Equation ( 24)).The in-situ data and the satellite product are independent, according to Equation (4); therefore, the difference between in-situ data and the satellite product can be obtained by where δE insitu−satellite is the difference between the in-situ data and a satellite-derived product, δE insitu is the error of the in-situ data (also named the uncertainty on the measurement of the in-situ data), δE satellite is the error of the satellite-derived product.Then, we can calculate the error of the satellite-derived product in the VIS and NIR bands.
Because the error of in-situ data is independent of the satellite measurement, the total error of the satellite-derived product can be computed as δR rs−satellite = δR rs−insitu−satellite − δR rs−insitu (10) where δR rs−insitu−satellite is the difference between the in situ and satellite-derived R rs , δR rs−insitu is the error of the in situ R rs , and δR rs−satellite is the error of satellite-derived R rs .

Decomposition of the Total Error
For NASA and OLCI atmospheric correction algorithms, R rs (λ) can be derived using the equation For a given band, when the observation geometry and aerosol type are known, then t 0 (λ)cosθ 0 F 0 (λ) is also known, and can be considered a constant [43].Then, the error of R rs (λ) can be derived using Equations ( 5)-( 8), leading to Thus, the error of R rs (λ) includes two parts: the error of L rc (λ), which is called the error of the Rayleigh-corrected radiance and the error of L A (λ), which includes the error of the iterative model and the aerosol LUTs.

Error of the Iterative Model
For the NASA and OLCI AC algorithms, R rs (NIR) is estimated using an iterative model, and the error of R rs (λ) is the difference between the R rs (NIR) from the iterative model and the R rs (NIR) from the in-situ data.Thus, the error of R rs (NIR) is where δR rs (NIR) is the error of R rs (λ) at NIR bands, R rs−ite (NIR), the R rs (λ) from the iterative model in the NIR bands, and R rs−true (NIR), the true/in situ R rs (λ) in the NIR bands.
Because the error budget relies on the calculation/iterative process, by determination the δR rs (NIR), the error budget of the error of the iterative model is calculated as follows: (a) The error on R rs (NIR) is passed to L A (NIR) using Equation (11), and the error of L A (NIR) due to the iterative model (δL (b) Considering the laws of aerosol radiance in atmospheric correction algorithms [23,33,44], L A (λ) can be derived as where α is the Ångström exponent or Ångström coefficient, which is derived from aerosol LUTs using the aerosol optical thickness τ.
For a given band and a known aerosol type, , and the error of L A (VIS) from the iterative model (δL A−ite (NIR)) can be derived using Equations ( 4)-( 7) (10), and the error of R rs (VIS) from the iterative model (δR rs−ite (VIS)) is where δR rs−ite (λ) is the error of the iterative model in the VIS bands.Because t o (λ) is not included in the L2 standard products, it can be obtained from [45] Remote Sens. 2019, 11, 2400 where τ r (λ), τ A (λ), and τ O are the Rayleigh, aerosol, and ozone optical thicknesses, respectively; ω A (λ) is the aerosol single scattering albedo; and η R (λ) and η A (λ) are the Rayleigh and aerosol forward scattering probabilities, respectively.The τ r (λ) value was computed using Bodhaine et al. (1999) [27].
The value of τ O was taken as a constant equal to 0.008.The value η R (λ) was also taken as a constant equal to 0.5, while η A (λ)ω A (λ) was allowed to spectrally vary from 0.89 at 412 nm to 0.86 at 670 nm.The value θ 0 was extracted from the OLCI L2 radiometric product, and from the orbitons for VIIRS and MODIS.The value of τ A (λ) was extracted from the L2 standard products.

Error of the Aerosol LUTs and the Rayleigh-Corrected Radiance
Given the total error of the satellite-derived product (Equation ( 10)) and the error of the iterative model (Equation ( 17)), the error of the aerosol LUTs and the error of the Rayleigh-corrected radiance can be calculated using Equation ( 4).The error of the aerosol LUTs and the Rayleigh-corrected radiance is where δR rs−LUTs (λ) is the error of the aerosol LUTs, and δL rc (λ) is the error of Rayleigh-corrected radiance.The error of L rc (λ) includes the errors of both L t (λ) and L r (λ).

In-Situ Data
The in-situ data were collected in the coastal and offshore waters of the China Sea (see Figure 1 and Table 2).The stations used during the four campaigns were conducted around the Bo Sea (BoS), Yellow Sea (YS), East China Sea (ECS), the Pearl River Estuary (PRE), and South China Sea (SCS); these stations cover most of the China Sea.The R rs (λ) were obtained using two above-water optical instruments, SAS and ASD, following the NASA optical protocols [46].The SAS and ASD were calibrated before and after each campaign, and the R rs (λ) (unit = sr −1 ) were estimated from the in situ radiometric parameters using the R06 approach [47,48] R rs (λ) = L t (λ) − ρL sky (λ) where L t (λ) is the total radiance, L sky (λ) is the sky irradiance, E s (λ) is the total irradiance, λ is the wavelength of incident light, and ρ is the Fresnel reflectance of the air-water interface (ρ is a function of wind speed and cloud cover).Finally, a residual sun glint or white offset correction ([a*R rs (780) − R rs (720)]/a−) was implemented with a = 2.35 [48].The spectral backscattering coefficient (b b , m −1 ) using Hydroscatt-6 (HS6) [49] and the spectral absorption coefficient (a, m −1 ) using AC9 [49] were collected during only two campaigns (YS, ECS&SCS).

Quality Control of the In Situ Data
Quality control of the in-situ data was essential because the in-situ data were considered as the sea-truth and because we sampled the R rs values using two different radiometers.

Consistency of Multiple Measurements
Each in situ R rs was measured three times (three replicates) at each sampling station.The samples were excluded if the coefficient of variability (CV = standard deviation (R rs )/mean (R rs )) of the three R rs spectra at 490 nm at any given station was higher than 5% (Figure 2).Based on this threshold, 18 stations were removed.

Removal of the Surface-Reflected Radiance
The surface-reflected radiance (mainly due to sky glint because direct sun glint is usually weak with measurement geometry) is recognized as a major sources of uncertainty [50,51]; thus, a station  The surface-reflected radiance (mainly due to sky glint because direct sun glint is usually weak with measurement geometry) is recognized as a major sources of uncertainty [50,51]; thus, a station was discarded when the ρL sky /L t exceeded 50%, [50,52].Based on this threshold, eight stations were removed.

Comparison with an IOP Model
R rs can also be estimated from the IOPs of seawater, i.e., the backscattering (obtained here using the HS6 instrument) and absorption coefficients (using the AC9 instrument).The R rs can be modeled using the IOPs, and the modeled R rs values can be compared to the ASD/SAS in situ R rs measurements [53,54] using the equation The absolute difference between the in situ R rs and the modeled R rs is calculated as  The surface-reflected radiance (mainly due to sky glint because direct sun glint is usually weak with measurement geometry) is recognized as a major sources of uncertainty [50,51]; thus, a station The coefficients of the model presented in Equation ( 22) are from the GSM model [55,56].If di f f was larger than 30%, the corresponding in situ measurement was discarded.The PRE and BoS campaigns were not included in this closure exercise because no IOPs measurements were collected for these areas due to the lack of AC9 and HS6 data.The R rs values from 17 stations were removed based on this method.The remaining R rs values are depicted in Figure 3.Some stations were kept when only one band exceeded the threshold, which is the case for stations where the di f f was above 30% at 650 and 676 nm.The mean absolute percentage difference (APD, Equation ( 25)) and the relative percentage difference (RPD, Equation ( 26)) between the modeled and in situ R rs values are shown in Figure 3. Overall, there is good consistency between the in situ and modelled R rs in the VIS; the APD varies from 11.82% (412 nm) to 25.35% (676 nm), and the values increase with increasing wavelengths [57,58].
The total number of the discarded stations is 43 (some values were discarded because there were several R rs for a given station).After undergoing this quality control process, the data from 167 stations remained and were used for the match-up exercise.

Match-Up Procedures
The match-up analyses and quality control were performed for the R rs values from VIIRS, MODIS, and OLCI radiometric products.First, strict rules for time interval and spatial distance to the match-up data were applied.We selected the mean R rs values over a 0.034×0.034degree (VIIRS), 0.05 × 0.05 degree (MODIS and OLCI RR) rectangle, from the L2 product within a time window of ±3 h.The spatial distance the same as in-situ data, corresponding to about a 5 × 5-pixel box near the nadir, but this pixel box size was adjusted at the scan edge to address spatial heterogeneity as suggested by Barnes et al. (2019) [59].Second, views and solar zenith angles above 60 • and 70 • , respectively, were removed for all match-up data.Moreover, the L2 flags of the match-up data were applied.The L2 flags for OLCI include INVALID, LAND, CLOUD, CLOUD_AMBIGUOUS, CLOUD_MARGIN, SNOW_ICE, SUSPECT, HISOLZEN, SATURATED, HIGHGLINT, WHITECAPS, AC_FAIL, OC4ME_FAIL, ANNOT_TAU06, RWNEG_O2, RWNEG_O3, RWNEG_O4, RWNEG_O5, RWNEG_O6, RWNEG_O7, RWNEG_O8, and the 490 nm reflectance of is below 0.02.The L2 flags for MODIS and VIIRS include ATMFAIL, LAND, HIGLINT, HILT, HISATZEN, STRAYLIGHT, CLDICE, COCCOLITH, HISOLZEN, LOWLW, CHLFAIL, NAVWARN, ABSAER, MAXAERITER, ATMWARN, NAVFAIL [60].Third, because the accuracy of R rs (488/490) is usually the highest [59], a coefficient of variability CV (CV = STD/mean) was calculated for all boxes at 490 (488) nm in the space distance was calculated and a match-up was considered valid when the CV was smaller than 0.15 at 490 nm [60].Finally, the percentage of valid pixels in each box was checked, and when this percentage was no less than 50%, the mean values of the valid pixels in the box was calculated and compared to the in-situ data.
The absolute difference between the in situ  and the modeled  is calculated as The coefficients of the model presented in Equation ( 22) are from the GSM model [55,56].If  was larger than 30%, the corresponding in situ measurement was discarded.The PRE and BoS campaigns were not included in this closure exercise because no IOPs measurements were collected for these areas due to the lack of AC9 and HS6 data.The  values from 17 stations were removed based on this method.The remaining  values are depicted in Figure 3.Some stations were kept when only one band exceeded the threshold, which is the case for stations where the  was above 30% at 650 and 676 nm.The mean absolute percentage difference (APD, Equation ( 25)) and the relative percentage difference (RPD, Equation ( 26)) between the modeled and in situ  values are shown in Figure 3. Overall, there is good consistency between the in situ and modelled  in the VIS; the APD varies from 11.82% (412 nm) to 25.35% (676 nm), and the values increase with increasing wavelengths [57,58].
The total number of the discarded stations is 43 (some values were discarded because there were several  for a given station).After undergoing this quality control process, the data from 167 stations remained and were used for the match-up exercise.

Statistical Method
Several statistical parameters were used to evaluate the match-up results.The definitions of these parameters are given below: The bias: The absolute percentage difference (APD): The relative percentage difference (RPD): The root mean square error (RMSE): The unbiased RMS difference (u ∆ ) [61]: with y i , the ith satellite-retrieved value, x i , the ith in-situ value, and n, the number of match-up data points.

Variability of the In-Situ R rs Data
The total dataset of R rs for the four campaigns are shown in Figure 4.The R rs spectra include both open and coastal waters.SCS is mainly composed of open ocean waters [62], while YS and BoS are SPM-dominated waters [63].We classified the R rs values following the method developed in Wei et al. ( 2016) [64].The R rs of the in-situ dataset corresponds to the classes 2-14 and 16 from Wei et al. ( 2016) [64].Figure 4 shows a color code that depends on the classes defined by Wei et al. ( 2016) [64].The blue plots are typical of open ocean waters (classes 2-4), whose R rs values depend only on phytoplankton.The green spectra correspond to more optically complex waters, including those with CDOM and phytoplankton in the water (classes 5-9).The black spectra correspond to SPM-dominated waters (classes 10-14 and 16).For open ocean waters (blue spectra), we can observe the well-known shape of the spectra with high values of R rs at 410 and 443 nm and a decreasing magnitude from the blue bands to the NIR.For the intermediate waters (green spectra), we can observe that the highest magnitude is observed in the green bands and that their magnitudes are low in the NIR bands.Due to the high CDOM absorbance in the purple/blue bands, the R rs at 410 and 443 nm are lower than those in the green bands.In turbid waters (black spectra), the highest magnitude is observed in the yellow bands and high magnitudes also exist in the NIR bands (compared to the open ocean and intermediate waters) due to the strong SPM backscatter.In conclusion, our R rs dataset covers a large set of water types in the China Sea.
The spatial location and the R rs spectra of the match-ups are shown in Figure 5.The number of match-up stations for OLCI, MODIS, and VIIRS are 13, 13, and 15, respectively.Both the intermediate waters (green) and the turbid waters (black) are coastal waters.Figure 5 shows the match-up stations of the three sensors covering the coastal and open ocean waters.The match-up stations of MODIS and VIIRS cover SCS, YS, and BoS.The match-up stations of OLCI cover SCS, ECS, and YS.The R rs from the match-up stations of OLCI and MODIS included all three types of waters.The R rs from the match-up stations of VIIRS included only two water types (open ocean and turbid waters).The match-up stations of three sensors cover most of the China Sea and most of the water types.
observe that the highest magnitude is observed in the green bands and that their magnitudes are low in the NIR bands.Due to the high CDOM absorbance in the purple/blue bands, the  at 410 and 443 nm are lower than those in the green bands.In turbid waters (black spectra), the highest magnitude is observed in the yellow bands and high magnitudes also exist in the NIR bands (compared to the open ocean and intermediate waters) due to the strong SPM backscatter.In conclusion, our  dataset covers a large set of water types in the China Sea.

Validation Results
Scatterplots showing the relationship between the in situ and satellite-derived R rs are shown in Figure 6 (OLCI) and Figure 7 (VIIRS and MODIS) and the statistical parameters are listed in Table 3.For OLCI, the points of the scatterplots are mainly distributed around the 1:1 line.The scatterplots show that the satellite-derived values are in good agreement with the in situ R rs values, except at 400 and 410 nm, which corresponds to the campaign in PRE (coastal waters; the red circle in Figure 6).The bias values show that the satellite-derived R rs values are under-estimated compared to the in-situ R rs , except at 400, 443, 490, and 510 nm.Then, in the visible spectra, the R rs are mainly over-estimated and under-estimated in the NIR (negative bias values).The RMSE and u ∆ are low for all bands, with values below 0.003 sr −1 , which means that the satellite-derived values are in good agreement with the in situ R rs .The highest RMSE value is obtained at 400 nm (0.003 sr −1 ).The statistical parameters showed that the APDs at 400, 410, 443, 490, 510, 560, 620, 665, 674, 681, 708, and 754 nm are 43%, 30%, 23%, 18%, 20%, 21%, 68%, 57%, 77%, 66%, 56%, and 79% for OLCI, respectively.The R rs values at 490 nm show the lowest APD making it the most accurate band.
waters (green) and the turbid waters (black) are coastal waters.Figure 5 shows the match-up stations of the three sensors covering the coastal and open ocean waters.The match-up stations of MODIS and VIIRS cover SCS, YS, and BoS.The match-up stations of OLCI cover SCS, ECS, and YS.The  from the match-up stations of OLCI and MODIS included all three types of waters.The  from the match-up stations of VIIRS included only two water types (open ocean and turbid waters).The match-up stations of three sensors cover most of the China Sea and most of the water types.

Validation Results
Scatterplots showing the relationship between the in situ and satellite-derived  are shown in Figure 6 (OLCI) and Figure 7 (VIIRS and MODIS) and the statistical parameters are listed in Table 3.For OLCI, the points of the scatterplots are mainly distributed around the 1:1 line.The scatterplots show that the satellite-derived values are in good agreement with the in situ  values, except at 400 and 410 nm, which corresponds to the campaign in PRE (coastal waters; the red circle in Figure 6).The bias values show that the satellite-derived  values are under-estimated compared to the in-situ  , except at 400, 443, 490, and 510 nm.Then, in the visible spectra, the  are mainly overestimated and under-estimated in the NIR (negative bias values).The RMSE and  ∆ are low for all The MODIS/VIIRS derived-R rs scatterplots are shown in Figure 7, and the statistical parameters are listed in Table 3.For OLCI, the points of the scatterplots are mainly distributed around the 1:1 line and show that the satellite-derived values are in agreement with the in situ R rs for MODIS/ VIIRS, other than a few some outliers (red circle), where the in situ R rs were collected in BOS (coastal waters).The Bias values showed that the satellite-derived R rs values were under-estimated compared to the in situ R rs in all bands (except at 678 nm for MODIS, perhaps due to chlorophyll fluorescence).The RMSE and u ∆ are low in all bands, and have values below 0.002 sr −1 , which means that the satellite-derived values are in good agreement with the in situ R rs .The highest RMSE values were obtained at 412 nm (MODIS) and at 443 nm (VIIRS).The statistical parameters show that the APDs at 412, 443, 488, 531, 555, 645, 667, and 678 nm are 39%, 30%, 26%, 19%, 19%, 48%, 32%, and 40% for MODIS, respectively, while the APDs at 410, 443, 486, 510, and 671 nm were 46%, 42%, 30%, 27%, and 36% for VIIRS, respectively.The R rs (531) and R rs (551) are the most accurate for MODIS and VIIRS, respectively.Although PRE and BOS are located in different regions of the China Sea, their water types (which belong to coastal waters) are similar, and the aerosols are affected by terrestrial sources, which might be the main source of the outliers for these three sensors.
MODIS, and VIIRS in the bands between 490 nm and 754 nm.The UV band (400 nm) is the least accurate for OLCI.The APD value of OLCI is lower than that of MODIS, and VIIRS in the 412, 443, 490, and 510 nm bands.The APD of OLCI is higher than for MODIS and VIIRS in the NIR bands (from 600 nm to 678 nm).Although their match-up stations are different, the performance of OLCI is similar to that of MODIS in the visible bands.MODIS, and VIIRS in the bands between 490 nm and 754 nm.The UV band (400 nm) is the least accurate for OLCI.The APD value of OLCI is lower than that of MODIS, and VIIRS in the 412, 443, 490, and 510 nm bands.The APD of OLCI is higher than for MODIS and VIIRS in the NIR bands (from 600 nm to 678 nm).Although their match-up stations are different, the performance of OLCI is similar to that of MODIS in the visible bands.The spectral variation of the Bias, RMSE, APD, and u ∆ for OLCI, MODIS and VIIRS are shown in Figure 8.The trend of RMSE, APD, and u ∆ are the same for OLCI, MODIS and VIIRS.The bias of OLCI is positive in the 400, 443, 490, and 510 nm bands, which differs from MODIS and VIIRS.The absolute value of the bias for OLCI is lower than that for MODIS and VIIRS in most of the bands (from 410 nm to 620 nm).R rs (400) is the least accurate band for OLCI.Among the three sensors, the RMSE and u ∆ are consistent in term of both trend and value.The u ∆ of OLCI is lower than that of MODIS, and VIIRS in the bands between 490 nm and 754 nm.The UV band (400 nm) is the least accurate for OLCI.
The APD value of OLCI is lower than that of MODIS, and VIIRS in the 412, 443, 490, and 510 nm bands.The APD of OLCI is higher than for MODIS and VIIRS in the NIR bands (from 600 nm to 678 nm).Although their match-up stations are different, the performance of OLCI is similar to that of MODIS in the visible bands.

Difference between Open Ocean and Coastal Waters
Although the mean bias or APD values differ between the satellites (see Figure 8), we noticed that the satellite-derived R rs has a small difference in the open ocean and a large difference in the coastal water (see Figures 6 and 7).The reason may be the surface-reflected radiance which is one of the main sources of uncertainties in satellite R rs retrievals [65,66].These effects are small in the open ocean with small aerosol loads and much more pronounced in coastal waters with higher aerosol loads, because they strongly depend on the aerosol optical thickness [67].The phenomenon that the APD of three sensors decreases with increasing band is also the same influencing factor, because the sky radiance reflected from the water surface, which carries uncertainties to the satellite-derived R rs , is highest in the blue part of the spectrum [59].
Although the mean bias or APD values differ between the satellites (see Figure 8), we noticed that the satellite-derived R rs were higher than the in situ values in the visible bands over open ocean waters, which is a different situation than coastal waters (Figures 6 and 7).The RPD values from OLCI, MODIS, and VIIRS are shown in Figure 9, showing that the three satellites have similar results and, except for a few points, the RPD values of the R rs are positive over open ocean but negative in coastal waters.

Difference between Open Ocean and Coastal Waters
Although the mean bias or APD values differ between the satellites (see Figure 8), we noticed that the satellite-derived  has a small difference in the open ocean and a large difference in the coastal water (see Figures 6 and 7).The reason may be the surface-reflected radiance which is one of the main sources of uncertainties in satellite  retrievals [65,66].These effects are small in the open ocean with small aerosol loads and much more pronounced in coastal waters with higher aerosol loads, because they strongly depend on the aerosol optical thickness [67].The phenomenon that the APD of three sensors decreases with increasing band is also the same influencing factor, because the sky radiance reflected from the water surface, which carries uncertainties to the satellite-derived  , is highest in the blue part of the spectrum [59].
Although the mean bias or APD values differ between the satellites (see Figure 8), we noticed that the satellite-derived  were higher than the in situ values in the visible bands over open ocean waters, which is a different situation than coastal waters (Figures 6 and 7).The RPD values from OLCI, MODIS, and VIIRS are shown in Figure 9, showing that the three satellites have similar results and, except for a few points, the RPD values of the  are positive over open ocean but negative in coastal waters.

Discussion
Few works exist that the error of  () from in-situ data [34,50,68].The error of  () from BOUSSOLE was reported to be 6% in the VIS bands [68], and the error of in-situ data from AERONET-OC was approximately 5% from 412 nm to 551 nm, and 8% at 667 nm [34].The difference between directly measured  () and  () obtain using the optimization approach (bio-optical model) is approximately 11% in all bands [50].The data in Figure 3 show that the difference in  () between the GSM model and the in-situ data is approximately 14% in the VIS bands and 26% in the NIR bands.We assume that the error of  () from the GSM model is the same as the error of the in-situ data.Then, using Equation ( 3), the error of in-situ  () is 7% in the VIS bands and 13% in the NIR bands.
As discussed in Section 2.2, we estimated an error budget.We used 754, 678, and 671 nm as the

Discussion
Few works exist that the error of R rs (λ) from in-situ data [34,50,68].The error of R rs (λ) from BOUSSOLE was reported to be 6% in the VIS bands [68], and the error of in-situ data from AERONET-OC was approximately 5% from 412 nm to 551 nm, and 8% at 667 nm [34].The difference between directly measured R rs (λ) and R rs (λ) obtain using the optimization approach (bio-optical model) is approximately 11% in all bands [50].The data in Figure 3 show that the difference in R rs (λ) between the GSM model and the in-situ data is approximately 14% in the VIS bands and 26% in the NIR bands.We assume that the error of R rs (λ) from the GSM model is the same as the error of the in-situ data.Then, using Equation (3), the error of in-situ R rs (λ) is 7% in the VIS bands and 13% in the NIR bands.
As discussed in Section 2.2, we estimated an error budget.We used 754, 678, and 671 nm as the NIR bands for OLCI, MODIS, and VIIRS, respectively, due to the lack of R rs (NIR) in standard MODIS and VIIRS products.We first calculated the error of R rs and the error of the iterative model; then, we calculated the percentage of the error of the iterative model in the error of R rs .Finally, we calculated the average result according to the location of stations (open ocean or coastal waters).We investigated the factors that impacted the errors on the estimation of R rs in the open ocean and coastal waters using the equation presented in Section 2.2.

Influencing Factors in the Open Ocean
The mean error of R rs for the three sensors in the open ocean is shown in Figure 10a, revealing that the three sensors have similar R rs errors, all of which are relatively high in the blue band (approximately 0.0008 sr −1 ) and then decrease as the band increases (falling to approximately 0.0002 sr −1 for the red band).In terms of numerical values, VIIRS shows the highest error, followed by OLCI and MODIS, which show the lowest errors at 443 and 486 nm, respectively.The 400 nm band of OLCI has the largest error.

Influencing Factors in the Coastal Waters
Similar operations were conducted for coastal waters.The mean error of  for the three sensors in the coastal waters is shown in Figure 11a.The figure shows that the  error trends of the three satellites are similar to those in the open ocean, but the blue bands show different behaviors.In the bands shorter than 500 nm, the  error of OLCI is the highest, and the  error of MODIS is the smallest.For bands above 510 nm, the  error of OLCI and MODIS are both less than the  error of VIIRS.We can also observe that the  error of OLCI increases at the 410 and 443 nm bands.The  errors of the three satellites over coastal waters are higher than those over open ocean, especially for OLCI in bands shorter than 450 nm.
The mean percentage of the error of the iterative model in the error of  in the coastal waters is shown in Figure 11b, where the iterative algorithm error of MODIS has the highest proportion in the total error of  except at 531 nm, and it is less than 50% from 440 to 555 nm, which are the main VIS bands, while the ratio of the iterative algorithm error of OLCI and VIIRS is less than 50% in the total error of  in all bands.Although the impact of the iterative algorithm on the  error is higher in MODIS, the aerosol estimation and the error of the Rayleigh-corrected radiance in the coastal waters remain the major error sources for all three sensors over coastal waters.The mean percentage of the iterative model error in the error of R rs over the open ocean is shown in Figure 10b, showing that the percentages of the iterative model error in the error of R rs from the three sensors have similar trends, all have high percentages in the red bands and lower values in the blue band, except for 410 nm (OLCI) and 412 nm (MODIS).The ratio of the iteration error of OLCI is less than that MODIS and VIIRS in the bands before 500 nm.For the open ocean region of the China Sea, the ratio of the iteration error of the three satellites from 400 to 600 nm in the total error is less than 50%, which indicates that the main influencing factors of the three satellites are not the iterative algorithm, but the aerosols LUTs (δL A (λ)) [65] and the Rayleigh-corrected radiance (δL rc (λ)).This conclusion was also reached by Zibordi et al. (2018) [20].

Influencing Factors in the Coastal Waters
Similar operations were conducted for coastal waters.The mean error of R rs for the three sensors in the coastal waters is shown in Figure 11a.The figure shows that the R rs error trends of the three satellites are similar to those in the open ocean, but the blue bands show different behaviors.In the bands shorter than 500 nm, the R rs error of OLCI is the highest, and the R rs error of MODIS is the smallest.For bands above 510 nm, the R rs error of OLCI and MODIS are both less than the R rs error of VIIRS.We can also observe that the R rs error of OLCI increases at the 410 and 443 nm bands.The R rs errors of the three satellites over coastal waters are higher than those over open ocean, especially for OLCI in bands shorter than 450 nm.

Influencing Factors in the Coastal Waters
Similar operations were conducted for coastal waters.The mean error of  for the three sensors in the coastal waters is shown in Figure 11a.The figure shows that the  error trends of the three satellites are similar to those in the open ocean, but the blue bands show different behaviors.In the bands shorter than 500 nm, the  error of OLCI is the highest, and the  error of MODIS is the smallest.For bands above 510 nm, the  error of OLCI and MODIS are both less than the  error of VIIRS.We can also observe that the  error of OLCI increases at the 410 and 443 nm bands.The  errors of the three satellites over coastal waters are higher than those over open ocean, especially for OLCI in bands shorter than 450 nm.
The mean percentage of the error of the iterative model in the error of  in the coastal waters is shown in Figure 11b, where the iterative algorithm error of MODIS has the highest proportion in the total error of  except at 531 nm, and it is less than 50% from 440 to 555 nm, which are the main VIS bands, while the ratio of the iterative algorithm error of OLCI and VIIRS is less than 50% in the total error of  in all bands.Although the impact of the iterative algorithm on the  error is higher in MODIS, the aerosol estimation and the error of the Rayleigh-corrected radiance in the coastal waters remain the major error sources for all three sensors over coastal waters.The mean percentage of the error of the iterative model in the error of R rs in the coastal waters is shown in Figure 11b, where the iterative algorithm error of MODIS has the highest proportion in the total error of R rs except at 531 nm, and it is less than 50% from 440 to 555 nm, which are the main VIS bands, while the ratio of the iterative algorithm error of OLCI and VIIRS is less than 50% in the total error of R rs in all bands.Although the impact of the iterative algorithm on the R rs error is higher in MODIS, the aerosol estimation and the error of the Rayleigh-corrected radiance in the coastal waters remain the major error sources for all three sensors over coastal waters.

Summary
The OLCI, MODIS, and VIIRS R rs products were evaluated against in situ measurements through a match-up analysis in the China Sea.This analysis showed that the satellite-derived-R rs are in good agreement with the in-situ data and that the R rs from the three sensors are overestimated over the open ocean and underestimated over the coastal waters around China.The highest uncertainty was observed in the NIR bands for the three missions.The bands with the best performances for the three missions were all between 488 and 560 nm, which agree with the results of previous studies.The performance of OLCI is good (the APD of R rs at 490 nm is the most accurate with values of 18%) and close to or better than MODIS (the APD of R rs at 531 nm is the most accurate with values of 19%) around the China Sea.
The calculated error budget showed that the iterative model used in the atmospheric correction algorithm of OLCI has a smaller impact on the error of R rs compared to the NASA iterative model used to process the MODIS and VIIRS data over the open ocean around China.The main influencing factors on the error of R rs for OLCI, MODIS, and VIIRS are the aerosol LUTs and the Rayleigh-corrected radiance (the error of system vicarious calibration, the error of Rayleigh scattering, whitecaps, gas and

Figure 1 .
Figure 1.Stations of the China Sea used in the study.

Figure 1 .
Figure 1.Stations of the China Sea used in the study.

Figure 1 .
Figure 1.Stations of the China Sea used in the study.

Figure 3 .
Figure 3. Difference between the in-situ  and the  estimated from the GSM model using the in-situ IOPs.

Figure 3 .
Figure 3. Difference between the in-situ R rs and the R rs estimated from the GSM model using the in-situ IOPs.

Figure 4 .
Figure 4. In-situ  spectra at OLCI bands: (a) in the BoS, (b) in the YS, (c) in the PRE, (d) in the ECS&SCS.The blue spectra correspond to the open ocean, the black spectra to the turbid waters, the green corresponds to intermediate waters between open ocean and turbid waters.

Figure 4 .
Figure 4. In-situ R rs spectra at OLCI bands: (a) in the BoS, (b) in the YS, (c) in the PRE, (d) in the ECS&SCS.The blue spectra correspond to the open ocean, the black spectra to the turbid waters, the green corresponds to intermediate waters between open ocean and turbid waters.

Figure 5 .
Figure 5. In-situ  spectra and the map of match-up stations: (a) OLCI, (b)VIIRS, (c)MODIS, (d) the geographical map with the location of the match-up for each sensor.The color code is same as Figure 4.

Figure 5 .
Figure 5. In-situ R rs spectra and the map of match-up stations: (a) OLCI, (b)VIIRS, (c) MODIS, (d) the geographical map with the location of the match-up for each sensor.The color code is same as Figure 4.

Figure 6 .
Figure 6.Scatterplot of OLCI-derived versus in situ remote-sensing reflectance in the China Sea, total R 2 = 0.83.The blue scatterplots correspond to the open ocean, the black scatterplots to the coastal waters.The red circles correspond to the outliers for the PRE region (the APD of the outlier points is larger than the mean APD +/− 3*sigma of the other points.).

Figure 7 .
Figure 7. Scatterplot of the satellite-derived versus in situ remote-sensing reflectance in the SCS (left: MODIS-AQUA, total R 2 = 0.76; right: VIIRS, total R 2 = 0.80).The color code and the outliers are similar as in Figure 6.The red circles correspond to the outliers for the BoS region.

Figure 6 .
Figure 6.Scatterplot of OLCI-derived versus in situ remote-sensing reflectance in the China Sea, total R 2 = 0.83.The blue scatterplots correspond to the open ocean, the black scatterplots to the coastal waters.The red circles correspond to the outliers for the PRE region (the APD of the outlier points is larger than the mean APD +/− 3*sigma of the other points).

Figure 6 .
Figure 6.Scatterplot of OLCI-derived versus in situ remote-sensing reflectance in the China Sea, total R 2 = 0.83.The blue scatterplots correspond to the open ocean, the black scatterplots to the coastal waters.The red circles correspond to the outliers for the PRE region (the APD of the outlier points is larger than the mean APD +/− 3*sigma of the other points.).

Figure 7 .
Figure 7. Scatterplot of the satellite-derived versus in situ remote-sensing reflectance in the SCS (left: MODIS-AQUA, total R 2 = 0.76; right: VIIRS, total R 2 = 0.80).The color code and the outliers are similar as in Figure 6.The red circles correspond to the outliers for the BoS region.

Figure 7 .
Figure 7. Scatterplot of the satellite-derived versus in situ remote-sensing reflectance in the SCS (left: MODIS-AQUA, total R 2 = 0.76; right: VIIRS, total R 2 = 0.80).The color code and the outliers are similar as in Figure 6.The red circles correspond to the outliers for the BoS region.

Figure 9 .
Figure 9. RPD of R rs at each station from OLCI (a), MODIS (b), VIIRS (c).The vertical red dashed line shows the limit between the open and coastal waters.

21 Figure 10 .
Figure 10.Mean error of  from OLCI, MODIS, and VIIRS in open ocean (a) and the mean percentage of the error of the iterative model in the error of  from OLCI, MODIS, and VIIRS in the open ocean (SCS) (b).

Figure 10 .
Figure 10.Mean error of R rs from OLCI, MODIS, and VIIRS in open ocean (a) and the mean percentage of the error of the iterative model in the error of R rs from OLCI, MODIS, and VIIRS in the open ocean (SCS) (b).

Figure 11 .
Figure 11.Mean error of R rs from OLCI, MODIS, and VIIRS in coastal water (a) and the mean percentage of the error of the iterative model in the error of R rs from OLCI, MODIS, and VIIRS in coastal water (b).

Table 1 .
Symbols and definitions.

Table 2 .
Location, period, number, and source of R rs samples used in the study.

Table 3 .
Statistical parameters of OLCI, MODIS, and VIIRS around the China Sea.