Atmospheric correction of Sentinel-3/OLCI data for mapping of suspended particulate matter and chlorophyll-a concentration in Belgian turbid coastal waters

The performance of different atmospheric correction algorithms for the Ocean and Land Colour Instrument (OLCI) on board of Sentinel-3 (S3) is evaluated for retrieval of water-leaving radiance reflectance, and derived parameters chlorophyll-a concentration and turbidity in turbid coastal waters in the Belgian Coastal Zone (BCZ). This is performed using in situ measurements from an autonomous pan-and-tilt hyperspectral radiometer system (PANTHYR). The PANTHYR provides validation data for any satellite band between 400 and 900 nm, with the deployment in the BCZ of particular interest due to the wide range of observed Near-InfraRed (NIR) reflectance. The Dark Spectrum Fitting (DSF) atmospheric correction algorithm is adapted for S3/OLCI processing in ACO-LITE, and its performance and that of 5 other processing algorithms (L2-WFR, POLYMER, C2RCC, SeaDAS, and SeaDAS-ALT) is compared to the in situ measured reflectances. Water turbidities across the matchups in the Belgian Coastal Zone are about 20 – 100 FNU, and the overall performance is best for ACOLITE and L2-WFR, with the former providing lowest relative (Mean Absolute Relative Difference, MARD 7 – 27%) and absolute errors (Mean Average Difference, MAD -0.002, Root Mean Squared Difference, RMSD 0.01 – 0.016) in the bands between 442 and 681 nm. L2-WFR provides the lowest errors at longer NIR wavelengths (754 – 885 nm). The algorithms that assume a water reflectance model, i.e. POLYMER and C2RCC, are at present not very suitable for processing imagery over the turbid Belgian coastal waters, with especially the latter introducing problems in the 665 and 709 nm bands, and hence the chlorophyll-a and turbidity retrievals. This may be caused by their internal model and/or training dataset not being well adapted to the waters encountered in the BCZ. The 1020 nm band is used most frequently by ACOLITE/DSF for the estimation of the atmospheric path reflectance (67% of matchups), indicating its usefulness for turbid water atmospheric correction. Turbidity retrieval using a single band algo-rithm showed good performance for L2-WFR and ACOLITE compared to PANTHYR for e.g. the 709 nm band (MARD 15 and 17%), where their reflectances were also very close to the in situ observations (MARD 11%). For the retrieval of chlorophyll-a, all methods except C2RCC gave similar performance, due to the RedEdge band-ratio algorithm being robust to typical spectrally flat atmospheric correction errors. C2RCC does not retain the spectral relationship in the Red and RedEdge bands, and hence its chlorophyll-a concentration retrieval is not at all reliable in Belgian coastal waters. L2-WFR and ACOLITE show similar performance compared to in situ radiometry, but due to the assumption of spatially consistent aerosols, ACOLITE provides less noisy products. With the superior performance of ACOLITE in the 490 – 681 nm wavelength range, and smoother output products, it can be recommended for processing of S3/OLCI data in turbid waters similar to those encountered in the BCZ. The ACOLITE processor for OLCI and the in situ matchup dataset used here are made available under an open source license.


Introduction
Ocean colour remote sensing has become a mature science, with operational products for the open ocean produced daily on a global scale for a number of sensors, at spatial resolutions of a few hundred metres to a few kilometres.The standard algorithms used for the atmospheric correction (AC) have a long heritage, developed in the 1980's and 1990's for satellite sensors such as the Coastal Zone Colour Scanner (CZCS), the Sea-viewing Wide Field-of-view Sensor (SeaWiFS), the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Medium Resolution Imaging Spectrometer (MERIS).These AC algorithms assume for the open ocean that the signal observed in the Near-InfraRed (NIR) part of the spectrum is entirely caused by atmospheric and airwater interface effects, and extrapolate the signal to the visible wavelength bands using a suite of aerosol models (Gordon and Wang, 1994;Antoine and Morel, 1999).Large parts of the ocean colour observations are contaminated by severe sun glint, and by a new approach to the problem, i.e. modeling the water signal and fitting the atmosphere and interface effects, POLYMER achieved excellent performance and dramatically increased data availability over the oceans, by effectively removing sun glint (Steinmetz et al., 2011).POLYMER has been adapted as the AC of choice for the Ocean Colour Climate Change Initiative, which aims to make long term time-series from multiple sensors consistent and inter-comparable, in particular for the derived chlorophyll-a concentration (Müller et al., 2015;Sathyendranath et al., 2019).
It is well known that in turbid or highly productive waters, the assumption of negligible water reflectance in the NIR breaks down, and a number of approaches have been suggested to deal with non-zero NIR water reflectance.Typically, ocean colour processing schemes adopt an assumption of spatially constant aerosol type and a simple water model (Ruddick et al., 2000;Hu et al., 2000;Goyens et al., 2013), or an iterative estimation of non-zero NIR (Stumpf et al., 2003;Bailey et al., 2010) for processing turbid waters.A Bright Pixel Atmospheric Correction (BPAC) has been integrated in the ground segment processor for MERIS and the Ocean and Land Colour Instrument (OLCI) processing, based on the work of Moore and Lavender (1999) and Lavender et al. (2005).The use of combined atmosphere and water retrievals using Neural Network (NN) approaches has been frequently suggested for more complex water or aerosol types (Chomko and Gordon, 1998;Schiller and Doerffer, 1999;Jamet et al., 2005;Brajard et al., 2006;Doerffer and Schiller, 2007;Schroeder et al., 2007;Hieronymi et al., 2017).Methods that have an explicit water model can however be challenged by target water types that are not in their training dataset or model, and careful assessment may be required for specific applications.For the processing of extremely turbid waters, the use of ShortWave InfraRed (SWIR) bands at 1.6 and 2.2 μm has been quite successful (Wang and Shi, 2005;Gao et al., 2007;Vanhellemont and Ruddick, 2015;Ibrahim et al., 2019), albeit at the cost of introducing additional noise (Werdell et al., 2010).Not all sensors have these SWIR bands, they are notably absent on MERIS and OLCI, and other methods may be needed for processing of extremely turbid waters.Gossn et al. (2019) proposed an alternative method to process OLCI data over extremely turbid waters, using the BaseLine Residual (BLR) from three pairs of spectrally close bands to separate aerosol and water contributions, and taking advantage of the new 1020 nm SWIR band.
Recently, the Dark Spectrum Fitting (DSF) atmospheric correction has been proposed for metre-and decametre-scale satellite sensors over turbid and inland waters (Vanhellemont and Ruddick, 2018;Vanhellemont, 2019aVanhellemont, , 2019bVanhellemont, , 2020)).The DSF does not assume which band has negligible water reflectance, and may even use non-water targets to estimate the atmospheric path reflectance, similar to the maximum aerosol optical depth estimated by Guanter et al. (2010).Depending on the aerosol model, the impact of residual target reflectance in the NIR can in fact have larger impacts than residual target reflectance in the Blue bands, due to the magnification of NIR errors towards the visible bands by the extrapolative fitting methods, and the general shape of atmospheric path reflectance.ACOLITE/DSF for decametre scale sensors, specifically those on Landsat 8 and Sentinel-2, has found many applications in coastal and inland waters, such as the estimation of turbidity (Braga et al., 2020;Balasubramanian et al., 2020;Ciancia et al., 2020), bathymetry, (Bué et al., 2020;Caballero and Stumpf, 2020) and algal bloom monitoring (Molkov et al., 2019;Saberioon et al., 2020).
In the present paper, the adaptation of the DSF for S3/OLCI processing of turbid coastal and inland waters is introduced, and the performance of a number of S3/OLCI products is evaluated using an autonomous pan-and-tilt hyperspectral radiometer (PANTHYR, Vansteenwegen et al. (2019)) deployed in Belgian coastal waters.The processing of S3/OLCI is of particular interest to remote sensing of optically complex coastal waters due to its very complete spectral coverage at moderate spatial resolution (300 m), with 21 bands of which about 16 will provide useful surface reflectances between 400 and 1020 nm.Full resolution products from the standard water processing, including the Bright Pixel correction (L2-WFR), as well as top-of-atmosphere L1 data processed with freely available processors, POLYMER (Steinmetz et al., 2011), C2RCC (Brockmann et al., 2016), SeaDAS/l2gen (Gordon and Wang, 1994;Bailey et al., 2010), and ACOLITE/DSF (present paper), are evaluated using in situ matchups.The overall performance of the different algorithms is discussed, and the error budgets and their implications are explored for typical turbid water applications such as the retrieval of turbidity and chlorophyll-a concentration.

Study area
The Belgian Coastal Zone (BCZ) is a relatively shallow (<50 m) and well-mixed shelf sea, connected to the North Sea to the north, and the English Channel to the west (Ruddick and Lacroix, 2006).It is characterised by a relatively high suspended sediment concentration, with a gradient from several hundreds of g ⋅ m − 3 nearshore to <1 g ⋅ m − 3 in the offshore waters, inversely related to the bathymetry.Turbidity in Formazine Nephalometric Units (FNU), shows about a 1:1 relationship to the suspended sediment concentration, and both measurements are of interest to users in the BCZ (Nechad et al., 2009(Nechad et al., , 2010;;Neukermans et al., 2012).Strong tidal currents occur (exceeding 1 ms − 1 ), with the tidal ellipse oriented along the coast from south-west to north-east, with currents flowing north-east during flood tide.The tidal resuspension of sediments is the main cause of the high turbidity in the nearshore area (Fettweis and Van den Eynde, 2003;Fettweis et al., 2007).In near-bed measurements, the sediment concentration reaches the highest concentration at the end of ebb and the beginning of flood, as a result of resuspension during current maxima (Baeye et al., 2011).The same pattern is observed during neap and spring tides, although at lower concentration during neap tide (Baeye et al., 2011).Annually recurring spring and summer phytoplankton blooms, generally consisting of diatoms are observed, with important spatial variability throughout the BCZ (Rousseau et al., 2006;Muylaert et al., 2006;Lacroix et al., 2007).In recent years, the blooms have been occurring earlier, likely in response to sea surface temperature increases and changes in nutrient outputs (Desmit et al., 2020), and the two bloom peaks may have merged to one long growing season (Raitsos et al., 2014).Blooms of Phaeocystis globosa generally occur in between the diatom blooms, with a 4-13 week duration, and in most years, the biomass of the Phaeocystis blooms greatly exceeds that of the diatom blooms (Breton et al., 2006;Gypens et al., 2007).The use of multispectral reflectance information is not sufficient to distinguish which species is causing a bloom (Astoreca et al., 2009), but this can perhaps be done qualitatively (Lubac et al., 2008), or through evaluation of other characteristics, such as the presence of foam lines on higher resolution satellite imagery.Absolute chlorophyll-a concentration can be retrieved from satellite remote sensing data however, and is already important for Belgium's eutrophication reporting requirements for the European Commission's Marine Strategy Framework Directive (MSFD, 2008/56/EC and addenda).This reporting requires the estimate of the 90th percentile of chlorophyll-a in a European nation's waters, something that is generally only achievable through the use of satellite remote sensing data (Van der Zande et al., 2019;Gohin et al., 2019).In the BCZ, in situ measurements are performed every few months for about ten reference stations, and hence remotely sensed information provides crucial spatial and temporal components for Belgium's MSFD reporting.

In situ data
In situ data was collected using a prototype PANTHYR ( Vansteenwegen et al., 2019)

deployed on the Blue Innovation Platform
Research Tower 1 (RT1), just in front of the port of Oostende (51.2464 • N, 2.9193 • E).Fig. 1 shows the location of RT1 in the turbid coastal waters around Oostende, and the structure of the tower.PAN-THYR has two TriOS RAMSES radiometers mounted on a pan-and-tilt head, one for up-and downwelling spectral radiances, and one with a cosine collector to measure spectral irradiance.The PANTHYR measures autonomously every 20 min at programmed relative azimuth angles to the sun.In the present study, measurements were made at a 270 • azimuth angle.A measurement cycle consists of sequential scans of 3× spectral irradiance (E d ), 3× downwelling radiance (L d ), 11× upwelling radiance (L u ), 3× L d and 3× E d .Such a measurement cycle for a single azimuth angle takes about 1 min.Measurements are calibrated, dark current corrected, and resampled to a common wavelength grid, from 350 to 900 nm at 2.5 nm steps.Individual scans are subjected to a quality control as in Ruddick et al. (2006) and then averaged to provide a single equivalent measurement, if at least 9/11 E d and 5/6 of the L d and L u scans pass the quality control.The water-leaving radiance reflectance, ρ w is then computed using: where L w is the water-leaving radiance: where ρ f is the effective Fresnel coefficient linearly interpolated from (Mobley, 1999) for the sun zenith angle at the time of measurement, and the wind speed retrieved from the National Centers for Environmental Prediction (NCEP) 1 degree global model.No NIR correction is   performed and no water reflectance model is assumed in the in situ data processing.ρ w data were finally convolved to the relative spectral response (RSR) functions of the OLCI instruments on Sentinel-3 A and B (Sentinel-3 CalVal Team, 2016).In the period from 2019-12-11 to 2020-07-15, 1058 PANTHYR measurements made at 270 • relative azimuth passing quality control were available (Figs. 2 and 3).Waters at RT1 are rather turbid, with the 5-95 percentiles of ρ w 0.053-0.137at 560 nm, 0.024-0.109at 665 nm, and 0.021-0.090at 709 nm, with median values respectively 0.094, 0.060, and 0.052.The spectral percentiles are plotted in Fig. 2 for the full data archive and the obtained matchups (see further).
The surface sediment concentration in the near-shore Belgian coastal zone is mainly driven by tidal resuspension and horizontal advection.The tide measured in the port of Oostende is plotted in the background of Fig. 3 to show the neap-spring cycle.Some effects of the neap-spring tidal cycle can be seen in the in situ measured reflectance, as can impacts of strong wind events.Lower reflectances are observed just before, and higher reflectances just after neap tide, related to the resuspension of material that settled during the neap tide (Fettweis and Van den Eynde, 2003).Some reflectance peaks are associated with strong wind events, for example the mid February, end of March, and mid May storms (Fig. 4) all cause local peaks in the water reflectance.
A significant advantage of the PANTHYR/RT1 dataset, compared to AERONET-OC datasets (Zibordi et al., 2009) normally used for validation, is that the hyperspectral instrument allows validation of all OLCI VNIR bands in the range 400-900 nm, for OLCI specifically including the 400 nm and the various NIR bands that are not available from AERONET-OC and hence not yet validated.The RT1 location is in sufficiently turbid waters for a measurable water reflectance in the range 700-900 nm (Figs. 2 and 3).

Satellite data
The Ocean and Land Colour Instrument (OLCI) is a multispectral radiometer carried on board Sentinel-3A (launched in 2016) and B (launched in 2018) with 21 bands in the 400-1020 nm spectral range with high signal-to-noise ratio at an approximately 300 m spatial resolution.The two satellites in polar orbit can provide a revisit time of less than two days at the equator, with higher overpass frequencies at higher latitudes.S3/OLCI imagery for the RT1 site was retrieved as top-ofatmosphere (L1) full resolution (FR, 300 m) data from the Copernicus Open Access Hub (https://scihub.copernicus.eu),as processed with Instrument Processing Facility (IPF) IPF-OL-1-EO version 06.08, and baseline water products (L2-WFR) were retrieved from the Copernicus Online Data Access (CODA) hosted by EUMETSAT (coda.eumetsat.int).The following processing algorithms were evaluated:

Matchups
Matchups were identified as satellite overpasses with in situ data measured within 30 min, either (1) with measurements bounding the overpass, or if no bounding measurements available, (2) the closest measurement.In situ data from (1) were linearly interpolated from the bounding measurement times to the overpass time, data from (2) were used as is.Due to the high temporal variability of suspended sediments at the site, the reflectance can change rapidly, and the use of interpolated data is preferred.Due to large spatial variability around the site at short spatial scales, a single pixel containing the station coordinates was extracted from the satellite data.For completeness the ACOLITE matchups are also repeated for the 3 × 3 box mean in Supplementary Data 1, Table S1.Automated quality control was performed using the These metrics emphasise that both the satellite and in situ data are estimations of the true value.Noise in the satellite reflectance was computed using a Laplace operator using approximated second derivatives: where ρ is the reflectance image, and i and j are the pixel coordinates in image space.The operator was not computed for pixels where any of the required pixels were masked.The Laplace operator gives the sum of second derivatives in both directions, and is not sensitive to smooth gradients in the data, e.g. a near-to-offshore turbidity gradient.It does pick up pixel-to-pixel variability, e.g.sharp fronts and noise, and is commonly used in image processing and computer vision for noise (Tai and Yang, 2008) and edge detection (Van Vliet et al., 1989;Wang, 2007).It can be used to detect noise in remotely sensed images after e.g.subtracting an edge map (Corner et al., 2003).Since the noise levels are quite different across the processors, a consistent edge map could not be generated.For the present application of estimating noise levels across many images, the Laplace operator is used directly for a relatively clear location, away from the edges resulting from turbidity fronts.

Parameter retrieval algorithms
Water turbidity (T) is computed using the algorithm of Nechad et al. (2009): with ρ w at 709 nm, and calibration parameters taken at 710 nm from Nechad et al. (2016): A T = 498.52(FNU), and C T = 0.1892.The concentration of turbid water chlorophyll-a (C) is computed using the commonly used RedEdge algorithm of Gons et al. (2005): where R M is the 709/665 ρ w ratio, ϕ a * the chlorophyll-a specific absorption at 665 nm, here taken as 0.016 m 2 ⋅ mg − 1 , and b b the total backscatter (m − 1 ) as estimated from ρ w at 779 nm:

Matchups
In the deployment period between 2019-12-11 and 2020-07-15, 392 images were available from the Data Hub, about equally distributed among  and .According to the criteria listed above, 46 common matchups (S3A: 18, S3B: 28) were identified for L2-WFR and ACOLITE/DSF.Of these 46 matchups, 10 had negative reflectance at 400 nm from either L2-WFR (5) or ACOLITE (6), and 7 with negative reflectance at 412 nm (L2-WFR: 2, ACOLITE: 5).These spectra are not removed as it could bias the matchups to those with conditions where one processor may outperform others.Including the flagging from POLYMER and C2RCC reduced these to 39 matchups (S3A: 17, S3B: 22), and by including SeaDAS/l2gen the common subset was further reduced to 27 matchups (S3A: 9, S3B: 18).The common dataset contains only 13 matchups where none of the processors retrieve negatives anywhere in the spectrum (S3A: 6, S3B: 7).The subset of matchups is rather representative of the observed reflectance range at RT1, with a slight bias to the more turbid waters (Fig. 2).Selected matchups between PANTHYR and the six processors are shown in Fig. 5, a full set of matchups and the data are provided in Supplementary Data 1 (Fig. S1, matchup plots) and 2 (CSV database with convolved PANTHYR data).The panels in Fig. 5 show from top to bottom, a general decrease in reflectance due to a decrease in suspended sediments from winter to summer.High reflectances are still observed in March (2020-03-03 and 2020-03-20), with a Green peak reflectance at about 0.15.The third and fourth panels show spectra near the height of the spring bloom (2020-04-22) and in a second, more patchy bloom, observed at the end of June (2020-06-20), with Green peak reflectance respectively about 0.08 and 0.04.Both show a strong Red band chlorophyll-a absorption feature, and a RedEdge reflectance peak almost equaling the Green.For the matchups, the ranges of turbidity are about 20-100 FNU, and chlorophyll-a concentration of 0-60 mg ⋅ m − 3 , consistent across the different subsets of matchups.As expected, the band automatically selected for τ a determination across the 46 ACOLITE/DSF matchups, tended to the longer wavelengths, with a visible band (i.e. the 443 nm band) used in only 4 cases (9%).For the remaining 42 matchups bands >750 nm were used: 768 nm (6 cases, 13%), 885 nm (5 cases, 11%) and 1020 nm (31 cases, 67%).

Reflectance retrieval
Scatter plots for the full matchup dataset (46 scenes) for ACOLITE/ DSF are shown in Figs. 6 and 7.An increase of the RMSD is found towards the Blue and UltraViolet (UV) wavelengths, i.e. from about 0.005 at 865 nm to about 0.02 at 400 nm.These errors are consistent with the increase of atmospheric path reflectance from the NIR to the Blue.The MARD is <15% for all bands between 490 and 709 nm, reaching <8% for bands at 560 and 620 nm with the largest observed reflectance range.MARD are highest (>50%) for the shortest (400 and 412 nm) and longest wavelengths (865 and 885 nm), as a result of low data ranges, and, in the Blue, the difficult atmospheric correction.The remaining bands (442, 754, and 779 nm) show <30% MARD.Except for the bands at 709, 754, and 779 nm, the MAD are in general negative indicating a slight overestimation of the atmospheric path reflectance, and hence an underestimation of the water reflectance.
A comparison of all six processors is shown in Fig. 8 for a common subset of 27 matchups for three of the generally best performing bands, 560, 665, and 709 nm, which are typically used for retrieval of turbidity and turbid water chlorophyll-a.For these bands, the best performances are found for L2-WFR (10.4-13.6%MARD) and ACOLITE (7.4-11.2%MARD), with ACOLITE giving better results at 560 and 665 nm, and both having nearly equal performance at 709 nm.The ACOLITE results for this subset are close to the full 46 matchup results (Figs. 6-7), with a slightly better performance for the 27 matchup subset, likely as more challenging scenes are filtered out here.Both configurations of SeaDAS give decent results, with MARD ≤25%, and SeaDAS-ALT outperforming the default SeaDAS configuration for these turbid waters, likely because of the inclusion of the 1020 nm band.C2RCC gives a MARD between 25 and 30%.Note that the slope of the C2RCC regression line is well below 1 at 665 nm, indicating an underestimation of Red reflectance, which is especially pronounced at high reflectances.This effect can also be observed in the second panel of Fig. 5 (2020− 03− 21), where no Red band dip -characteristic of chlorophyll-a absorption -is observed in the in situ data, or in any of the processors except C2RCC.This is an indication that the C2RCC NN is limited by its training dataset, and that for these turbid waters, in order to minimise fitting errors, the NN has to resort to using the highest reflectance waters in its training dataset, which happen to have this Red absorption dip.This effect may explain the poor performance of the "case 2" water products in the MERIS 3rd reprocessing in turbid waters, where "winter blooms" were observed in the turbid coastal waters of the southern North Sea (Vanhellemont, 2012).Finally, POLYMER shows the highest MARD (40-57%), and gives significantly lower reflectances compared to in situ measurements, and indicates that the water reflectance model in POLYMER is unable to represent the higher reflectances of the turbid waters in the BCZ.Q. Vanhellemont and K. Ruddick The error metrics are represented spectrally for the common matchup datasets of 39 scenes (L2-WFR, POLYMER, C2RCC, ACOLITE, solid lines) and 27 scenes (including SeaDAS and SeaDAS-ALT, dashed lines) in Fig. 9. Two general shapes of the spectral errors are observed, one from the processors that constrain only weakly the water reflectance (L2-WFR, SeaDAS, and ACOLITE), and another from the processors that impose strong constraints on the water reflectance through a full VIS-NIR water reflectance model (POLYMER, C2RCC).For the first class, the RMSD increases from the NIR to the Blue, consistent with higher atmospheric signals to be removed in the shorter wavelength bands.The MARD shows a minimum in the visible bands, typically between 490 and 709 nm, i.e. in the bands with the highest observed reflectances, with sharp increases of relative errors to the Blue (high atmospheric signal) and the NIR (low water-leaving signal).For the second class, the RMSD peaks in the visible bands, with largest (absolute) MAD also found in the visible bands.The MARD is a bit flatter for this class, with only slight increases towards the NIR and Blue/UV.All three metrics have pronounced spectral features for the second class, where the water model in the processor presumably deviates from the observed reflectances, e.g. a lower Green reflectance peak for POLYMER, and differences in the Red/ RedEdge bands for C2RCC.In terms of MARD, ACOLITE gives the best performance across the bands between 442 and 681 nm, while the L2-WFR has lower relative errors between 709 and 885 nm.The RMSD shows a similar pattern, with ACOLITE giving more scatter at 400 nm compared to L2-WFR.All processors give negative MAD across the spectrum, apart from a handful of bands, indicating a tendency of overcorrecting for atmospheric effects, or a too low top-of-atmosphere calibration.ACOLITE gives lowest (absolute) MAD of around 0.002, with SeaDAS-ALT and C2RCC having a lower bias at 400 nm, and L2-WFR and SeaDAS-ALT having very comparable biases in the NIR.Note that system vicarious calibration gains are applied only for L2-WFR/S3A processing (EUMETSAT, 2019), and no gains are applied for L2-WFR/ S3B, SeaDAS and ACOLITE processing.POLYMER applies by default the gains obtained from S3A processing, and for C2RCC a common set of vicarious calibration gains are applied for S3A and B. For these processors constraining the water reflectance, the relative band-to-band calibration is of more importance than the absolute calibration.A summary of the two best performing methods, L2-WFR and ACOLITE is given in Table 1 (400-674 nm) and 2 (681-885 nm), for the full 46 matchup database that both methods provided results for.Summary statistics for the 13 scenes in the common dataset where no processor gave any negatives are provided in Supplementary Data 1, Tables S2 and  S3.

Derived parameters
Fig. 10 shows time-series of ρ w 665 and 709 nm, two bands that are frequently used to derive water turbidity or total suspended matter (Nechad et al., 2009(Nechad et al., , 2010)), and chlorophyll-a (Gons et al., 2005;Gilerson et al., 2010;Moses et al., 2012).Due to the high observation density of OLCI from combining S3A/B data, effects of tidal and wind driven resuspension are preserved in the satellite data.Generally the reflectance is higher at 665 nm than at 709 nm, except in the period from April to mid-May, where there is a sharp decrease in 665 nm reflectance, while the 709 nm reflectance remains at similar levels to before -this effect can also be seen more clearly in Fig. 3.This pattern is caused by the spring phytoplankton bloom, with the associated increase of chlorophyll-a absorption.Fig. 11 shows time-series of chlorophyll-a concentration, as derived by the RedEdge algorithm of Gons et al. (2005), showing a phytoplankton bloom occurring from about April to mid-May, and a second, short but intense bloom at the end of June, as is systematically observed in these waters, e.g.(Breton et al., 2006).The bloom causing the peak near the end of the time-series (2020-06-22) is visible as a dark patch on the last panel of Fig. 5, and the full extent of the bloom is plotted in the Graphical Abstract to this paper.
The two processors constrained by a water model (POLYMER and C2RCC) show a significant underestimation of the Red and RedEdge reflectance (Figs. 8,9,and,10).Despite the underestimation of reflectance, with biases at 665 and 709 nm of about − 0.03 and − 0.02, POLYMER retrieved the chlorophyll-a absorption feature, and was able to track the bloom development.It gives slightly higher chlorophyll-a values during the peak of the spring bloom, caused by the different magnitude of the biases at 665 and 709 nm "enhancing" the Red absorption feature.C2RCC provided significant chlorophyll-a values before the spring bloom, and during the blooms, the RedEdge peak as returned by C2RCC is reduced, and hence it gives a less pronounced difference in retrieved chlorophyll-a between non-bloom and bloom periods.For these turbid waters, even in the absence of a chlorophyll-a absorption feature, C2RCC provides a water reflectance spectrum including such a feature.This is presumably a result of the nature of the training dataset and the error minimisation procedure inherent to the NN.In the presence of a strong absorption feature, C2RCC will give a reduced and flattened RedEdge to the chlorophyll-a absorption dip.These effects can be seen in the example matchups, respectively in the 2nd and 3rd panel of Fig. 5, and are also clearly visible in the time-series of Fig. 11, where C2RCC consistently overestimates the chlorophyll-a, except during the April-May and end of June blooms, where it retrieves similar values to before the bloom.
These results are summarised in the scatter plots shown in Fig. 12 and Tables 3 and 4 for turbidity (Nechad et al., 2009) and chlorophyll-a (Gons et al., 2005) estimated from the reflectance data.For turbidity, L2-WFR and ACOLITE give the best performance, with around 15% MARD, and biases of 1-2 FNU and RMA slopes of 0.96-1.01.This is no surprise, as these processors gave the performance for the retrieved reflectance in the band itself (Fig. 8).The other processors show a significant underestimation of turbidity with biases of − 6 to − 20 FNU and RMA slopes of 0.39 (POLYMER) and 0.54-0.71(others).For chlorophyll-a, the matchups were further reduced from 27 to 17, as POLYMER retrieved a negative ρ w 779 nm, making the b b impossible to estimate (Eq.( 9)), and SeaDAS-ALT produced 9 matchups with negative ρ w 779 nm.C2RCC retrieves a RMA slope of 0.32, with an offset of 12.5 mg ⋅ m − 3 and a bias of − 4 mg ⋅ m − 3 , showing significant overestimation at low and underestimation at high chlorophyll-a concentrations.L2-WFR gives rather high slope (1.45) and bias (11 mg ⋅ m − 3 ), and the other processors have slopes between 1.1 and 1.3 and biases of 5-11 mg ⋅ m − 3 , indicating a general overestimation of chlorophyll-a compared to the values estimated from the in situ measurements.SeaDAS and SeaDAS-ALT show the lowest MAD and MARD, likely due to their better performance at lower chlorophyll-a concentration, as only 5 of the common matchups have chlorophyll-a > 35 mg ⋅ m − 3 .L2-WFR and ACOLITE are however able to provide double the number of observations with similar biases to SeaDAS.The different reflectance biases of POLYMER at 665 and 709 nm observed before, result in a rather high bias (11 mg ⋅ m − 3 ) and offset (4 mg ⋅ m − 3 ) in the RedEdge chlorophyll-a retrieval.

Spatial noise
Noise can be introduced or amplified in per-pixel processing, due to extrapolation of noise from longer wavelength bands, and from model fitting differences in adjacent pixels.Since ACOLITE uses a spatially fixed τ a , or a smoothed interpolated τ a to retrieve per-pixel atmospheric parameters, its results are smoother, and closer to the band signal-tonoise specifications, than other processors.In this section the noise level in the products from the two best performing processors in the matchup analysis is examined, i.e.ACOLITE and L2-WFR.An example of the L2-WFR and ACOLITE outputs are given in Fig. 13, showing speckled noise in the L2-WFR product.A comparison of noise levels was performed using a Laplace operator (Eq.( 6)) and extracted from the pixel containing 51.357 • N, 2.835 • E, a location with relatively clear waters, avoiding turbidity fronts.Scenes not covering this location, or having ρ t 1020 nm > 0.03 anywhere in the 3 × 3 pixel box were excluded, leaving 36 of the 46 matchup scenes.Fig. 14 shows an estimate of noiseequivalent reflectance in 16 OLCI bands at top-of-atmosphere (black), after L2-WFR (vermillion) and ACOLITE (grey) processing.ACOLITE results are only slightly noisier than the top-of-atmosphere input reflectances, where L2-WFR has considerably higher noise levels, especially for the bands <560 nm, and increasing with atmospheric path reflectance.ACOLITE/DSF does not smooth or amplify the noise present in the TOA data, which can also be observed in the bottom panel of Fig. 13.L2-WFR has slightly lower noise than the input reflectances at 885 and 1020 nm, possibly as a result of assumptions made on these bands during atmospheric correction.

Perspectives
In the present study, the baseline Sentinel-3 OLCI product (L2-WFR) has been confirmed to work well in the turbid Belgian coastal waters, as compared to 46 matchups with hyperspectral in situ measurements across winter and summer seasons (Tables 1 and 2).ACOLITE gave lower errors in the visible wavelengths compared to L2-WFR, and provides products with significantly less noise.POLYMER and C2RCC did not show good performance, the former due to the water model giving too low reflectances for these turbid waters, and the latter due to poor performance of the Red and RedEdge bands, i.e. by introducing artificial chlorophyll-a absorption features, or flattening the features that are present.Even though the OLCI processing in SeaDAS has not been formally evaluated by NASA/OBPG, the SeaDAS performance is reasonably good, albeit outperformed by L2-WFR and ACOLITE in these waters.The performance of a single band turbidity algorithm was of course closely related to that band's specific results for each processor, while a RedEdge band-ratio chlorophyll-a algorithm proved to be more robust to various (typical) atmospheric correction errors (Fig. 12 and Tables 3 and 4).Only the C2RCC chlorophyll-a retrieval was problematic due to (1) unrealistic Red band absorption features, and (2) flattening of the RedEdge peak being introduced by the NN fitting.Results presented here show that ACOLITE/DSF gives best performance for visible band reflectances (7-14% MARD from 490 to 681 nm), and is tied with L2-WFR for providing turbidity based on 709 nm reflectance (15-17% MARD), and does provide robust chlorophyll-a concentration (MARD 37-57%) estimates.
Results obtained by Renosh et al. (2020), show that the baseline

Table 1
Comparison of in situ and satellite ρ w for the full set of matchups from 400 to 674 nm.m and b are the RMA regression slope and offset.

Conclusions
• The Dark Spectrum Fitting (DSF) atmospheric correction algorithm in ACOLITE has been adapted and automated for processing of Sentinel-3/OLCI data, and it is here demonstrated and evaluated for the turbid waters in the Belgian Coastal Zone.   3 and 4.
Q. Vanhellemont and K. Ruddick Q. Vanhellemont and K. Ruddick • A single automated hyperspectral PANTHYR system, deployed at Research Tower 1 near Oostende, provided 46 matchups with the OLCI on Sentinel-3A ( 18) and Sentinel-3B ( 28), in about 7 months time, and data for all bands from 400 to 900 nm.This is an unprecedented dataset for evaluating OLCI performance in turbid coastal waters, with median reflectances observed over the 46 matchups of 0.094 at 560 nm, 0.060 at 665 nm, and 0.052 at 709 nm.• Six processors have been evaluated for mapping the turbid Belgian coastal waters using Sentinel-3/OLCI data, as compared to a number of autonomous hyperspectral in situ measurements.The best results are obtained for the processors that only weakly constrain the water reflectance, in particular the baseline L2-WFR (with a water model only in the NIR) and ACOLITE (only assuming non-zero water reflectance in the darkest pixels), with the latter providing the best performance in the visible bands (<15% MARD with a 0.01 RMSD and a near-zero MAD), with the lowest noise levels.• The two processors based on a full visible to NIR water reflectance model (POLYMER and C2RCC) showed a significant underestimation of the visible reflectance.C2RCC gave poor performance in particular for the 665 and 709 nm reflectance, either by introducing artificial Red band chlorophyll-a absorption, or by reducing the presence of a real RedEdge peak.These effects are both problematic for typical eutrophication monitoring applications in turbid waters.At present these processors cannot be recommended for turbid water monitoring using OLCI, and will likely show better performance with better adapted specific water models and training data.• Using the turbid water Red/RedEdge chlorophyll-a algorithm of Gons et al. (2005), all processors except C2RCC could reliably track in time the development of two blooms, a spring bloom in April-May and a summer bloom at the end of June.This algorithm was found to be robust to the absolute performance of the atmospheric correction, as long as the relative signal in both Red and RedEdge bands is maintained.• For mapping of turbidity or suspended particulate matter concentration in moderately turbid waters, the use of 560-709 nm bands is recommended (Nechad et al., 2009(Nechad et al., , 2010;;Novoa et al., 2017) and hence ACOLITE and L2-WFR can be recommended for this application, both giving comparable MARD at 709 nm (11%) and ACOLITE performing best at 560 nm (7%) and 665 nm (10%), and both processors giving similar scatter (RMSD of around 0.01).Slightly higher errors are found for retrieved turbidity due to the non-linearity of the turbidity algorithm, i.e. 15 and 17% for 709 nm derived turbidity.The lower spatial noise of ACOLITE is also seen in the turbidity product.• Due to the high combined observation density of Sentinel-3A/B, the turbidity changes in the Belgian coastal waters during neap-spring tide cycles and as a result of strong wind events, could be observed using satellite remote sensing.Sentinel-3A/B is the only satellite    mission that at present can provide spatially and temporally resolved information about coastal phytoplankton blooms in the BCZ.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix B. Processing with fixed or tiled τ a
Most ocean colour atmospheric correction algorithms are typically applied on a per-pixel basis, correcting a single observed ρ t spectrum at a time.
In contrast, the DSF works on a spatial subset of data in order to retrieve representative atmospheric parameters, which is especially useful for correcting difficult targets, e.g.extremely turbid or inland waters, and to reduce the amplification of noise during the atmospheric correction.The exploitation by DSF of the fact that atmospheric variability generally has a longer length scale than a single pixel has already proven very useful for higher resolution sensors with low signal:noise (Vanhellemont and Ruddick, 2018;Vanhellemont, 2019aVanhellemont, , 2019bVanhellemont, , 2020)), but is also useful for the better-specified OLCI sensor.The DSF can be applied using a spatially uniform or tiled estimation of aerosol optical thickness (τ a at 550 nm), a choice that depends on the size of the region of interest (ROI) and application.In this Appendix, the fixed (i.e.spatially uniform) and tiled options are described, and the performance analysed for various settings, in order to make an optimal choice for processing Sentinel-3/OLCI imagery for the turbid waters in the BCZ.In the spatially uniform τ a processing mode, a single model and τ a is retrieved for the ROI, and the atmospheric correction parameters are then retrieved from the LUT using the per-pixel geometry.Fig. A1 shows the MARD with 46 matchups with in situ PANTHYR measurements for different ROI sizes using fixed processing (F-prefix).Different subsets (6 × 6, 12 × 12, 24 × 24, 36 × 36, and 48 × 48 km) centred on RT1 provided quite consistent results, with a minimum MARD for the 12 and 24 km ROIs, with a larger MARD for the 6, 36 and 48 km ROIs.The larger MARD are caused by the selection of less representative ρ dark in the DSF, i.e. brighter pixels for the 6 km ROI, and darker pixels for the 36 and 48 km ROIs.The presence of darker pixels in the larger ROI is not necessarily bad, but they may provide a less representative (i.e.lower) estimate of the path reflectance.The 12 and 24 km ROI seem to provide the optimal presence of dark targets, while the likelihood of even darker targets increases in a larger search radius, causing a potential underestimation of the path reflectance in the larger ROIs.Compared to the 12 and 24 km ROIs, the 6 km ROI is too small to provide similar targets with negligible surface signal, and hence the path reflectance is overestimated.Note that the 6, 36, and 48 km ROIs still provide very good results, their performance is just not as good compared to the 12 and 24 km ROIs.
For tiled processing, the scene is divided into square tiles of given dimensions, and the DSF applied individually to each tile, retrieving an aerosol model and τ a per tile.For tiles where no result is found, the results from the closest neighbouring tile are used.To avoid tile edges in the final product, the retrieved τ a per model, and the selected model per tile, are linearly interpolated from the tile centres to each pixel, and smoothed using a Gaussian filter with a 24 pixel (7200 m) kernel standard deviation.This provides a per-pixel estimate of τ a per model, as well as a per-pixel model weighting factor.The atmospheric parameters are then retrieved using the per-pixel geometry, the retrieved τ a for both aerosol models, and are finally linearly weighted using the per-pixel model weighting factor.MARD for different tile sizes (12, 24 and 36 km) are also plotted in Fig. A1 (results with the Vprefix).An example of the tiled processing for a subscene over the southern North Sea is provided in Fig. A2.In comparison with 46 RT1 matchups, the 36 km tiled processing gives the best results, in terms of MARD (Fig. A1), and (absolute) MAD (not shown separately).The RMSD (i.e.scatter in the matchup comparison) is slightly higher for the tiled processing compared to the fixed processing at wavelengths >490 nm (not shown separately), as a result of the model mixing and spatial interpolation of atmospheric parameters.Small differences are found between the 24 km and 36 tiling grids, while for the 12 km tiling grid higher MAD and RMSD were retrieved, with about the same MARD.For processing of complete full resolution top-of-atmosphere (L1-FR) scenes or large areas of interest, the tiled processing with a 36 km tiling grid is recommended.For local studies, e.g. at estuarine scale, a fixed processing may provide better (less noisy) results.In the main paper a ROI of 36 × 36 km was processed, and hence the fixed τ a option was used, even though here lower errors were retrieved for the 36 km tiled processing.Both processing options are fully automated in ACOLITE and the choice may depend on a specific application.

Appendix C. Supplementary data
Supplementary data to this article can be found online at https://doi.org/10.1016/j.rse.2021.112284.

Fig. 1 .
Fig. 1.Left: the location of the Blue Innovation Platform RT1 just in front of Oostende plotted on a 12 × 12 km subset of the Sentinel-2A/MSI image taken 2020-03-31, and right: photos of the platform and the PANTHYR installation provided by Dieter Vansteenwegen (VLIZ).

Fig. 3 .
Fig. 3. Timeseries of PANTHYR data convolved to the OLCI RSR for four bands.Black/grey triangles denote the position of the 46/27 matchup subsets.Tide height above the TAW reference level as provided by Meetnet Vlaamse Banken (at station'Ostend Harbor') is plotted in grey.

•
L2-WFR: Baseline product from EUMETSAT/CODA as processed with IPF-OL-2 version 06.13 (EUMETSAT, 2019).Standard masking was used, i.e. excluding INVALID, LAND, CLOUD and CLOUD_AMBIG-UOUS pixels.• POLYMER: POLYMER (v.4.13) with current default settings, i.e. using the Park and Ruddick (2005) water model.The bitmask dataset was used to keep only pixels passing the recommended quality check, i.e. pixels with no mask (bitmask = 0) or pixels where flags 10 (CASE2) and 11 (INCONSISTENCY) were set.• C2RCC-ALT: The C2RCC ALTERNATIVE NN as provided in SNAP 7.0 with default settings.The c2rcc_flags dataset was used to reject pixels where flags 1, 3, or 4 (Rtosa_OOR, Cloud_risk, Iop_OOR) were set.• SeaDAS: SeaDAS/l2gen (9.5.0-V2019.3,git clone dated 2020-07-28) was used with default settings for OLCI, i.e. using the 2 band 779/ 865 multi-scattering algorithm with Relative Humidity based model selection (Ahmad et al., 2010) and iterative NIR correction (Bailey et al., 2010).Standard flagging was used, i.e. an external land mask and a threshold on the NIR reflectance for cloud masking (ρ t 865 nm > 0.027).• SeaDAS-ALT: An alternative processing for SeaDAS/l2gen, using the 2 band 865/1020 multi-scattering algorithm.Standard flagging was used -see above.Disabling the chlorophyll-a based iterative NIR and BRDF corrections (aer_opt = − 1, and brdf_opt = 3) did not provide outputs at 865 nm, so these options were kept as default (aer_opt = − 2, and brdf_opt = 7).• ACOLITE: A new OLCI version of the DSF was implemented in ACOLITE (details in Appendix A) and imagery was processed to water reflectances using a fixed aerosol optical thickness τ a retrieved from a 36 × 36 km region of interest (ROI) centred on RT1.Automated full scene processing is supported, and an analysis of other processing settings regarding the spatial variability of τ a is given in Appendix B.

Fig. 4 .
Fig. 4. Timeseries of average wind speed at 10 m (red) and wind direction (blue) from Meetnet Vlaamse Banken station'Ostend -weather station'.10 min data is plotted in the background, with solid foreground lines the 48 h moving average.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5.A few selected matchups of Sentinel-3/OLCI and the PANTHYR deployed at RT1.The RGB composites use surface reflectances at 665, 560, 490 nm scaled from 0 to 0.15 in the 8 bit R, G, B channels.The open red circle on the RGB composites is centred on RT1.The title for each plot shows the aerosol optical thickness at 550 nm (τ a ) and whether the Continental (C) or Maritime (M) aerosol model was used for ACOLITE/DSF processing.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.Time-series plots of ρ w at 665 nm (top) and 709 nm (bottom) for the different processors.In situ PANTHYR data is plotted in red.The grey background shows the tide as measured by Meetnet Vlaamse Banken at station'Ostend Harbor'.The second panel in each plot shows the difference between the satellite and the in situ measured ρ w .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12.Comparison of turbidity and chlorophyll-a concentration computed from satellite and in situ reflectances for 27 (turbidity) and 17 (chlorophyll-a) matchups.Statistics are shown in Tables3 and 4.

Fig. 13 .
Fig. 13.Top: ρ w 709 nm from the baseline L2-WFR product (left) and as processed by ACOLITE (right) for a subset of the S3A image of 2020-03-22.Middle: Noise calculated using a Laplacian (absolute value) for L2-WFR and ACOLITE.Bottom: Noise for the TOA and ACOLITE results.Pixels are masked using a 0.03 threshold on ρ t 1020 nm.

Fig. 14 .
Fig. 14.Spectral noise-equivalent reflectance as computed using a Laplacian on a relatively clear location in the BCZ.Solid lines show the median noise from 36 scenes, the dashed lines the interquartile range.

Fig. A1 .
Fig. A1.MARD between ACOLITE/DSF processed OLCI and in situ measurements for different ROI sizes for the fixed (F-) and tiled (V-) DSF processing for 46 matchups.The black solid line shows the 36 × 36km ROI with fixed processing as used in the main paper.

Fig. A2 .
Fig. A2.An example of the 36 × 36 km tiled processing of an OLCI scene (S3B, 2020-04-20) over the southern North Sea, showing (left) the lowest retrieved τ a at 550 nm, (middle) the interpolated and smoothed ρ path at 560 nm, and (right) the RGB surface reflectance composite.

Table 2
Same as Table 1 for bands 681-885 nm.

Table 4
Same as Table3but for chlorophyll-a concentration, with 36 common matchups (with ρ w 779 nm > 0) from the full dataset also given for L2-WFR and ACOLITE.

Table 3
Comparison of in situ and satellite derived turbidity computed from the 709 nm band for the common subset of 27 matchups.For L2-WFR and ACOLITE the full matchup dataset performance is also provided.