Atmospheric Correction of OLCI Imagery over Extremely Turbid Waters Based on the Red , NIR and 1016 nm Bands and a New Baseline Residual Technique

A common approach to the pixel-by-pixel atmospheric correction of satellite water colour imagery is to calculate aerosol and water reflectance at two spectral bands, typically in the near infra-red (NIR, 700–1000 nm) or the short-wave-infra-red (SWIR, 1000–3000 nm), and then extrapolate aerosol reflectance to shorter wavelengths. For clear waters, this can be achieved simply for NIR bands, where the water reflectance can be assumed negligible i.e., the “black water” assumption. For moderately turbid waters, either the NIR water reflectance, which is non-negligible, must be modelled or longer wavelength SWIR bands, with negligible water reflectance, must be used. For extremely turbid waters, modelling of non-zero NIR water reflectance becomes uncertain because the spectral slopes of water and aerosol reflectance in the NIR become similar, making it difficult to distinguish between them. In such waters the use of SWIR bands is definitely preferred and the use of the MODIS bands at 1240 nm and 2130 nm is clearly established although, on many sensors such as the Ocean and Land Colour Instrument (OLCI), such SWIR bands are not included. Instead, a new, cheaper SWIR band at 1016 nm is available on OLCI with potential for much better atmospheric correction over extremely turbid waters. That potential is tested here. In this work, we demonstrate that for spectrally-close band triplets (such as OLCI bands at 779–865–1016 nm), the Rayleigh-corrected reflectance of the triplet’s “middle” band after baseline subtraction (or baseline residual, BLR) is essentially independent of the atmospheric conditions. We use the three BLRs defined by three consecutive band triplets of the group of bands 620–709–779–865–1016 nm to calculate water reflectance and hence aerosol reflectance at these wavelengths. Comparison with standard atmospheric correction algorithms shows similar performance in moderately turbid and clear waters and a considerable improvement in extremely turbid waters.


Introduction
The National Aeronautics and Space Administration (NASA) SeaDAS standard atmospheric correction algorithm applied to MODIS (Moderate Resolution Imaging Spectroradiometer) imagery uses two bands in the near infra-red (NIR) at 765 and 865 nm, assuming zero water-leaving radiance contribution ("black water" hypothesis) at these bands (Gordon and Wang 1994 [1]) or, for turbid waters, modelling the NIR water reflectance with an iterative process (Stumpf et al., 2003 [2]).The aerosol reflectance at the NIR bands is then used to estimate aerosol properties and hence extrapolate aerosol reflectance (and diffuse transmittance) to the visible (VIS) bands.For extremely turbid waters, (e.g., ρ w (865 nm) > 0.02 with suspended particulate matter (SPM) concentration above 100 g/m 3 ), this procedure may fail because of one or more problems: P1.The NIR water reflectance spectrum becomes flatter (Doxaran et al., 2002 [3]) and hence more similar to aerosol reflectance spectra (Ruddick and Vanhellemont 2015 [4]); P2.At high reflectance the simple analytic models of NIR water reflectance as polynomial expansions of the Gordon parameter (the ratio of b b a+b b ) are no longer accurate (Ruddick et al., 2006 [5]) and better modelling is needed (Lee et al., 2016 [6], Luo et al., 2018 [7]); P3.Top of atmosphere (TOA) reflectance may even exceed the photodetector saturation values of certain ocean colour sensors rendering these wavelengths unusable (Dogliotti et al., 2011 [8]).
In these cases, improvements to the standard "NIR-extrapolation" procedure are needed.One approach is to apply the black water assumption at bands in the short-wave-infra-red (SWIR) (e.g., 1240, 1640 and 2130 nm in MODIS), where water absorption is much higher giving negligible water-leaving radiance even in highly turbid waters (Wang and Shi 2005 [9]).Although the European Space Agency's (ESA) Ocean and Land Colour Instrument (OLCI) on Sentinel-3 lacks far-SWIR bands (e.g., 1240 nm, 2130 nm in MODIS), it incorporates a new SWIR band at 1016 nm (nominally 1020 nm, although here the value of 1016 nm is preferred, which corresponds to the mean wavelength of the spectral response function of the band).While the black water assumption does not hold at 1016 nm for extremely turbid water (Knaeps et al., 2012 [10]), the pure water absorption is much higher at this band than at the NIR, meaning that the flattening of the NIR water reflectance spectrum in extremely turbid waters (P1) does not extend to 1016 nm-with inclusion of this band, water-reflectance spectra have clearly different shape from aerosol reflectance spectra and decomposition of the TOA (or Rayleigh-corrected) reflectance into water and aerosol reflectance components can be achieved.NIR/SWIR extrapolative algorithms are the most conventional but clearly not the only possible approach to atmospheric correction.Considerable success has also been achieved by full spectrum coupled water-atmosphere optical models which use an optimization or iterative approach to obtain a best-fitting combination of water and atmosphere unknowns, such as in, e.g., Doerffer and Schiller 2007 [11], Chomko and Gordon 2001 [12].Such approaches have even achieved extraordinary success in removing strong but spectrally simple sunglint (e.g., POLYMER, Steinmetz et al., 2011 [13]) and tend to always provide a water reflectance spectrum that looks like water, avoiding for example negative reflectances, and proving very robust even for data with calibration bias.The operational potential of such approaches is fully recognized, but the extrapolative approach is preferred here because each step of the process has a clear physical basis, facilitating troubleshooting and improvement.For example, a poor water reflectance model or inaccurate calibration or inappropriate aerosol model will be clearly made evident with an extrapolative approach but may be less obvious with a full spectrum coupled water-atmosphere approach.In this work, we design a pixel-by-pixel aerosol correction algorithm for turbid waters by identifying, for the OLCI Red-NIR-SWIR (RNS) bands, "atmospheric invariant" quantities which are sensitive to water reflectance but insensitive to aerosols, e.g., Philpot 1991 [14].This approach is based on the fact that in the RNS region, aerosol reflectance is spectrally smoother than water reflectance (Ruddick and Vanhellemont 2015 [4]).To calculate directly the water reflectance, assuming only that aerosol reflectance is spectrally smooth, a "baseline residual" (BLR) algorithm is proposed here for Rayleigh-corrected reflectances at 3 RNS bands.A baseline is formed between the shortest and longest wavelength of the triplet and is subtracted from the middle wavelength to give the BLR.This approach is similar to that adopted, for example, in the fluorescence line height product (Letelier and Abott 1996 [15]) and gives a quantity, the BLR, that is invariant under spectrally-linear (or sufficiently smooth) signal addition.By calculating three different atmospherically "quasi-invariant" BLRs, defined by the Rayleigh-corrected (RC) reflectances at triplets (620-709-779 nm), (709-779-865 nm) and (779-865-1016 nm), and by applying over them an "equivalent transmittance" correction factor, it is then possible to calculate water reflectance in the RNS on a pixel-by-pixel basis even in the case of extremely turbid waters.Once the RNS water reflectance is calculated, the atmospheric correction can then revert to standard NIR/SWIR extrapolative algorithms by removing water reflectance from the Rayleigh-corrected reflectance to give aerosol reflectance in the RNS and then estimating an aerosol type and concentration in the NIR/SWIR and using tabulated models to extrapolate aerosol reflectance to all wavelengths as in Gordon and Wang 1994 [1].Alternatively, if a full spectrum atmospheric correction is not required, it is possible to use the RNS water reflectance directly to estimate products such as turbidity (Dogliotti et al., 2015 [16], Dogliotti et al., 2016 [17]) or SPM (Nechad et al., 2010 [18], Shen et al., 2010 [19], Moreira et al., 2013 [20]), or chlorophyll-a concentration in turbid waters (Gilerson et al., 2010 [21], Gitelson 1992 [22], Gons et al., 2005 [23]).

Materials and Methods
The following subsections describe the different datasets used to calibrate and validate the algorithm.

In Situ Radiometric Measurements
In the present study, in situ above-surface radiometric measurements were used: (i) as input to radiative transfer simulations to compute a transmittance factor to be applied to BLRs (see Sections 2 and 2.6) to compare between in situ and OLCI-derived BLRs (see Section 3.1).A total of 105 water reflectance measurements were selected from several field campaigns (2012-2015) performed in the Río de la Plata (RdP) (Argentina), a funnel-shaped estuary located at Eastern-Central South America (Figure 1), with mean suspended particulate matter concentrations (SPM) in the range 100-200 g/m 3 , reaching up to 400 g/m 3 (Framiñan and Brown 1996 [24]).The Analytic Spectral Devices (ASD FieldSpec FR) spectroradiometer works in the spectral range of 350-2500 nm (step 1 nm) and was used to measure upwelling and downwelling radiances, just above surface (L 0+ u and L 0+ d , resp.) at relative azimuth to sun of either 135°or 90°and reciprocal zenithal angles of θ = ±40°(angles in accordance with Mueller et al., 2003 [25]), and downwelling irradiance (E 0+ d ) inferred from measuring nadir radiance at a calibrated quasi-lambertian plaque and multiplying it by π.These are combined to compute the above-surface water reflectance, ρ 0+ w : where ρ M (W) is the surface radiance reflectance factor taken from Mobley 1999 [26], which depends on the wind speed, W. The measurement protocol is detailed in Knaeps et al., 2012 [10], and follows the generic "Abovewater method 2" of the NASA 2003 Ocean Optics Protocols (Mueller et al., 2003 [25]).

Radiative Transfer Simulations
Simulations were performed to compare the behavior of BLRs computed from above-water reflectances, BLR(ρ w ), and BLRs computed from RC reflectances, BLR(ρ RC ), and also to establish an equivalent transmittance factor to be applied to BLRs computed from RC reflectances from OLCI scenes (see Section 2.6).In situ radiometric measurements were used as an input to the CNES-SOS v5.0 (Centre National d'Études Spatiales-Successive Orders of Scattering) radiative transfer (RT) code (Lenoble et al., 2007 [27]) to test how water reflectance baseline residuals (defined in Equations ( 5)-( 6)) relate to Rayleigh-corrected baseline residuals taking account of atmospheric transmittance.CNES-SOS was developed to solve the vector RT equation for marine and terrestrial environments using a plane-parallel assumption.In our simulations, the surface reflectance is composed of: (i) a component that accounts for air-water interface effects, modelled as a function of surface wind speed (Cox and Munk 1954 [28]), and (ii) a component, which is assumed to be Lambertian, that accounts for radiation that interacts with the in-water constituents (also called above-water reflectance, ρ w , Equation (1), or water-leaving radiance reflectance).
Combining all the in situ spectra with a set of multiple simulated atmospheric conditions, a total of 70875 Rayleigh-corrected (RC) reflectance spectra (Table 1) were computed with CNES-SOS RT code, corresponding to: 3 solar and viewing zenith angles at 0°, 30°and 60°(interpolated from the Gaussian quadrature angles [27]), 5 sun-relative azimuth angles in the range 0-180°(step 45°) (where 0°represents viewing into sunglint), 3 aerosol granulometries and refractive indices of type "Continental" (WMO-C), "Maritime" (WMO-M) and "Urban" (WMO-U) (from the World Meteorological Organization aerosol optical models [29]), 5 aerosol optical depths at 500 nm in the range 0-0.4 (step 0.1) and 105 in situ water reflectances (Section 2.1).The relative refractive index of the water-air interface was fixed at 1.334 throughout the spectrum and the wind speed fixed at 3 m/s.Rayleigh scattering from air molecules was set according to the values of Rayleigh optical thickness and depolarization factor reported by Bodhaine et al., 1999 [30], see also Mobley et al., 2016 [31].The molecular and aerosol vertical profiles were set to follow exponential decay functions of e-folding scales of 8 km and 2 km, respectively.The code was set to not exceed 20 orders of scattering and to include the effects of polarization.To obtain the Rayleigh-corrected reflectances of the set, the code was run twice to compute the Top-of-Atmosphere reflectance (ρ SOS TOA ): (i) using the input values specified by Table 1, and (ii) assuming a black marine environment (ρ w = 0), and without aerosol content (ρ a = 0).The second set was subtracted from the first to compute the RC signal, as it would be calculated in the atmospheric correction step of a satellite data processor:

Ocean and Land Colour Instrument (OLCI) Imagery
OLCI L1B and L2 imagery were downloaded from the Copernicus Online Data Access (CODA) system on 1 February 2018 [32].The L1B and L2 imagery used in this work correspond to the processing baseline v2.23 [33].The detailed list of scenes from Río de la Plata (RdP), Bahía Blanca (BBl, Argentina) and Belgian Coast (BE) used to calibrate and test performance of the algorithm is shown in Table 2. Before entering the BLR calculations, two preliminary corrections were applied to the L1B images.Firstly, a "PPE correction" is made to remove occasional pixel artifacts caused by prompt particle events (PPEs) following the algorithm proposed by Gossn 2018 [34].PPEs are caused by magnetospheric particles hitting OLCI's charge-coupled devices (CCDs) and generating across-track stripes of anomalous radiance values.PPEs are highly likely to occur over the South Atlantic Anomaly where Río de la Plata is located (D'Amico et al., 2015 [35]).Secondly, a "Rayleigh correction" is applied, using SeaDAS v7.5 (Bailey et al., 2010 [36], Stumpf et al., 2003 [2]), to remove the path radiance expected from a pure Rayleigh atmosphere with black sea.This gives the Rayleigh-corrected (RC) reflectance at all OLCI bands of interest, the ρ S product in SeaDAS.The NASA/SeaDAS and ESA standard atmospheric corrections were also applied to OLCI imagery to compare performances with the present scheme: (i) l2gen script from NASA's SeaDAS (Bayley et al., 2010 [36], Stumpf et al., 2003 [2]) was run with aerosol option "−2", i.e., assuming an iterative procedure to find the optimum aerosol model using bands 865 nm and 1016 nm (termed SeaDAS-2 (865, 1016) hereafter), and (ii) ESA's standard atmospheric correction (as performed for L2 imagery, processing baseline v2.23), that combines the baseline atmospheric correction (BAC, Antoine and Morel 1999 [37]), which is essentially a NIR-based black water approach, together with the bright pixel atmospheric correction (BPAC, Moore et al., 1999 [38]), which calculates first the NIR water reflectance from an iterative approach.

Algorithm Theoretical Basis
Consider the following expression for the total reflectance at TOA (cf.Gordon and Wang 1994, Equation (1) [1]): where L TOA is the total radiance at TOA, d the Sun-Earth distance in Astronomical Units, E TOA is the solar extraterrestrial irradiance at TOA and θ s is the solar zenith angle.A first step common to all atmospheric correction schemes is to estimate the Rayleigh-corrected reflectance (ρ RC ), which is obtained by subtracting the path radiance expected from a pure Rayleigh atmosphere with black sea (ρ R ) produced by air molecules from the total signal at TOA, Mobley et al., 2016 [31]: The aerosol correction scheme presented here is based on a pixel-by-pixel approach called here baseline residual (BLR), which is the difference between the signal at a "middle" (M) band (i.e., of intermediate wavelength in a triplet of bands) and the value of the baseline formed by the signals at the left (L) and right (R) bands at this middle wavelength, as illustrated in Figure 2. In this work, we use BLRs computed using Rayleigh-corrected (RC) reflectances, so we define BLR(ρ RC ) as follows: where λ L , λ M and λ R stands for the left, middle and right wavelength of the triplet, ρ RC (λ M ) is the RC reflectance at λ M and BL(ρ RC )(λ M |λ L , λ R ) is the baseline term, given by: Other precursor algorithms in the ocean colour field that use BLRs, include the fluorescence line height (FLH [15]), defined for MODIS as the BLR of the normalized water-leaving radiances at bands 665 nm, 677 nm and 746 nm, i.e., FLH = BLR(nL w )(665, 677, 746); the floating algal index (FAI, Hu 2009 [39]), defined also for MODIS as the BLR of Rayleigh-corrected reflectances at 645 nm, 859 nm and 1240 nm, i.e., FAI = BLR(ρ RC )(645, 859, 1240), the maximum chlorophyll index for MERIS (MCI, Gower and King 2008 [40]), which can be expressed as MCI = BLR(L TOA )(681, 709, 754), or the synthetic chlorophyll index (SCI, Shen et al., 2010 [19]), which is computed using remote-sensing reflectances as SCI = −BLR(R rs )(620, 665, 681) − BLR(R rs )(560, 620, 681) (these bands correspond to the cases of MERIS and OLCI).In some cases these approaches (FLH, SCI) are implemented after an aerosol correction, while in other cases (FAI, MCI) no aerosol correction is applied.Nevertheless, all such approaches are similar in the sense that BLR indexes are essentially unaffected by undesirable signal such as aerosol and/or moderate sun glint (and hence also to typical errors in the process of removal of this signal), since these components are generally spectrally smoother than the in-water signal, especially in the RNS region.Mathematically, if we consider the atmospheric signal "ρ atm " inside the spectral range determined by the band triplet (λ L , λ M , λ R ) to be as "smooth" as possible, i.e., of linear wavelength dependence, The atmospheric correction algorithm presented here was designed taking into consideration that the dependence of BLRs on atmospheric components is best minimized by computing BLRs from RC reflectances at spectrally-close bands, and for longer wavelengths (λ > 600 nm), to reduce the impact of uncertainties in Rayleigh correction, including coupled aerosol-Rayleigh scattering (see Section 2.5 for further reasons).This can be illustrated by considering the following decomposition of RC reflectance, (cf.Gordon and Wang 1994, Equation (2) [1]): where each term accounts for photons that arrive at the sensor after: having interacted with aerosols or with both aerosols and molecules (expressed in a summarized form as ρ a ); being specularly reflected by the water surface (sunglint) (T(λ)ρ g , where T(λ) is the direct atmospheric transmittance), and having interacted with the in-water components (t(λ)ρ w (λ), where t(λ) is the diffuse atmospheric transmittance).Based on typical aerosol types and concentrations such as the modes reported by the World Meteorological Organization (WMO) [29], or by Shettle and Fenn 1979 [41], aerosol reflectances can, in most cases, be modelled as exponential functions of wavelength if considering a sufficiently short bandwidth (see Gordon and Wang 1994 [1], Figure 1): where c is related to the aerosol type and ρ a (λ L ) is the amplitude at λ L .Typically (but not always) c > 0, describing a monotonically decreasing wavelength dependence.In particular, in a sufficiently short spectral range, e.g., 250 nm in the RNS, the aerosol reflectance can be approximated as a linear function of wavelength, at least compared to the stronger spectral variability of turbid water reflectance: If Equation ( 9) holds for the spectral range involved in the calculation of BLR (i.e., from λ L to λ R ), the aerosol term does not contribute to BLR(ρ RC )(λ L , λ M , λ R ), i.e.,: BLR(ρ a )(λ L , λ M , λ R ) ≈ 0 Finally, the sunglint term, at least in a moderate regime, is essentially a white term, especially in short spectral ranges in the RNS region, where we can consider ∂T(λ) ∂λ ≈ 0, i.e., negligible contribution to BLR(ρ RC ).Thus, considering any triplet of spectrally close bands in the RNS implies a near-linear spectral dependence of the atmospheric-interface terms in the RC reflectance decomposition.Thus the BLR(ρ RC ) depends mainly on the in-water term:

Band Choice
In this study, we have chosen to use the RC reflectances computed from the OLCI bands at 620, 709, 779, 865 and 1016 nm.These 5 bands span a R 5 space, from which we obtained three successive linearly-independent BLRs: BLR(ρ RC )(620, 709, 779), BLR(ρ RC )(709, 779, 865) and BLR(ρ RC )(779, 865, 1016).These bands have been chosen in order to i) maximize the impact of the (turbid) water signal on BLR(ρ RC ), and simultaneously ii) minimize the impact of the atmospheric components.Other OLCI bands inside this spectral region were not considered, e.g., those in the range 760-770 nm and the band centered at 940 nm, as they are strongly affected by absorption of atmospheric components such as O 2 (oxygen) and H 2 Ov (water vapour).OLCI bands in the range 660-690 nm are also avoided because the water reflectances there can be affected by chlorophyll absorption and fluorescence.

Transmittance Factor Treatment
The diffuse transmittance factor in Equation (10), accounts for the fact that photons may interact with atmospheric constituents on the path from sun to water and/or from water to sensor.Since it is a second order correction factor, we will assume an equivalent spectrally white diffuse transmittance inside each band triplet, whose expression might depend on geometric conditions and aerosol properties, i.e., as an approximation of Equation (10): where θ s and θ v account for the solar and sensor zenith angles, resp.; ∆φ is the sun-sensor azimuth angle difference, and AER stands for the dependence on aerosol type and concentration.Generally, we expect t BLR to be less than (but close to) unity.
To show how the atmosphere affects BLRs computed from the selected bands, Figure 3 shows the relation between BLR(ρ w ) computed from in situ measurements, (applying OLCI spectral response functions [42] and Equation ( 5)) and the corresponding BLR(ρ RC ), calculated from simulations made with the CNES-SOS radiative transfer code (see Section 2.2, Table 1).In Figure 3b-d, for each BLR(ρ w ) computed from the associated input water reflectance, a vertical line (inter-percentile range between percentiles 5 and 95, IPR(5,95)) is shown together with the extreme values (represented by " * "), corresponding to the set of computed RC reflectances for the combination of all possible atmospheric conditions described in Table 1.As a general rule, BLR(ρ RC ) tend to present lower absolute values than BLR(ρ w ), as suggested by the slopes of the linear regressions (i.e., m < 1).In terms of spectral shape, lower absolute BLRs are associated with spectrally smoother signals, which is consistent with an equivalent transmittance factor smaller than 1 (see Equations ( 10) and ( 11)).All scenarios associated to the min/max values, represented by " * ", correspond to cases where the sensor is positioned at reciprocal angles from the sun (i.e., direct sun glint).Apart from these exceptional cases, the low offsets computed for the linear regressions (b), are consistent with the assumptions taken to reach Equation (10), i.e., BLR(ρ a + Tρ g ) ≈ 0. It was observed, from the whole set of simulations, that the dependence of t BLR on the geometric conditions could be reduced to a single variable, the air mass factor, µ = 1 cos(θ s ) + 1 cos(θ v ) , since, except for the cases of direct sun glint, the relative azimuth was observed to have a very small effect on t BLR .This fact is well illustrated in Figure 4, where t BLR and a bias are plotted vs. the air mass factor µ considering the following expression: In Figure 4, each point represents the result of t BLR (upper insets) and bias (lower insets) of a linear regression of the form of Equation ( 12) done over subsets defined by a specific set of values for θ s , θ v and ∆φ, i.e., for specific geometric conditions.It can be observed that: (i) biases are usually smaller than 0.001, and will be neglected in the overall correction for transmittance, and (ii) the dependence of t BLR on µ can be considered linear in the range µ ∈ [2; 4].Similar results are observed for increasing aerosol optical thicknesses.11) and ( 12)) and bias vs. air mass factor, µ = 1 cos(θ s ) + 1 cos(θ v ) , for each of the BLRs used in this work.Turquoise (violet) dots represent subsets of simulations corresponding to different sun-view geometries with (without) direct sunglint.Dashed lines represent linear regressions.
On the basis of these results, the BLRs computed over RC OLCI scenes used to estimate the water signal are divided by the corresponding equivalent transmittance factors for each pixel (Figure 4), to estimate the BLRs from OLCI that will be associated to the water signal in OLCI scenes, i.e.,:

Relation between Baseline Residuals and Water Reflectances
The previous sections described how BLR(ρ w )(620, 709, 779), BLR(ρ w )(709, 779, 865) and BLR(ρ w )(779, 865, 1016) (Equations ( 5), ( 10) and ( 11)) represent a convenient set of quantities to estimate water reflectance at any of the considered bands.In this study, we focused on establishing a relation between the BLRs and water reflectances of at least two bands: 865 nm and 1016 nm; since this is the minimum required to eventually extrapolate the atmospheric signal to shorter wavelengths.We decided to use a subset of OLCI imagery to achieve a plausible calibration dataset, in order to avoid systematic errors induced by calibrating the algorithm with data coming from other sources (such as in situ measurements or radiative transfer simulations).To build this dataset, a total of 13 sub-regions from 9 cloud-, sunglint-and haze-free scenes were selected from Río de la Plata (see Table 2, magenta boxes at Figure 5).Over these subregions a non-operational atmospheric correction was applied to infer water reflectance at bands 865 nm and 1016 nm, based on estimating the atmospheric component from manually selected fixed clear water windows, close to the windows of interest, of 15 × 15 pixels, referred to herein as "Fix-AC" (see Table 2, white boxes at Figure 5), and assuming horizontally homogeneous aerosols.These "clear windows" were chosen to be as close as possible to the subregion of interest and totally free from water signal in the RNS (determined by visual inspection of the ρ RC rasters).This simple aerosol removal assumes a spatially uniform atmospheric signal which is subtracted from the whole subregion of interest as follows: where a simple expression is used for the two-way transmittance factor, t(λ, µ), (as applied by NASA OBPG in the absence of aerosol information and molecular absorption, Mobley et al., 2016 [31]) which depends on the air mass factor, µ, (see Section 2.6) and the Rayleigh optical thickness, τ Ray (λ) (Bodhaine et al., 1999 [30]): Once this process is performed, the calibration dataset formed by the mentioned subregions is used to fit a calibration surface in BLR 3-dimensional sub-space, which is spanned by BLR(620, 709, 779), BLR(709, 779, 865) and BLR(779, 865, 1016) (see Figures 5 and 6).The calibration surface is generated by binning on a 2-dimensional mesh-grid of "X" (BLR(620, 709, 779)) in the range (−0.0100; 0.0350) (step 0.0005) and "Y" (BLR(709, 779, 865)) in the range (−0.0300; 0.0150) (step 0.0005).The "Z" (BLR(779, 865, 1016)) values at each point of the grid are taken as the median of the Z values taken by the calibration dataset in the corresponding (X, Y) pair.This can be done in this way since the BLR surface does not present ambiguity in Z, i.e., the BLR surface can be considered as a function of (X, Y): Z = f (X, Y).The calibration point is discarded if less than 10 data fall inside the associated (X, Y) range.The BLRs(ρ w ) triplet computed at each input pixel is associated to the closest (in the Euclidean sense) BLRs(ρ w ) triplet that corresponds to the calibration curve, and the corresponding water reflectance at 865 nm and 1016 nm is assigned to the given pixel.

Estimation of Aerosol Reflectance at Bands 865 nm and 1016 nm
The aerosol reflectance, ρ a , at 865 and 1016 nm is computed to test the spatial correlation with the estimated marine signal -which is expected to be low-using the following simplified expression: where ρ w,j is the water reflectance value assigned to the nearest Euclidean neighbour from the calibrations surface at the BLR space.The transmittance factor is calculated in the same way as in Equation ( 15).This same expression was applied to intercompare with the aerosol reflectances yielded by other preexistent schemes (BAC/BPAC and SeaDAS-2 (865, 1016)).A last correction is performed on ρ a (865) to restrain the derived a (865, 1016) = ρ a (865) ρ a (1016) inside the range of (0.85; 1.25).These bounds were determined as the extreme values taken over a set of 82 different selected windows of size 15px × 15px from OLCI-A scenes of clear water regions close to Río de la Plata, Bahía Blanca, North Sea, Yellow Sea, Amazonas and North Australia.Also these bounds are consistent with what was obtained over the CNES-SOS simulations.This correction is performed by imposing the following condition on a pixel-by-pixel basis: and then, consequently correcting the retrieved ρ w (865) values.This constraint is further discussed in Section 4.

Atmospheric Correction Scheme: Summary
The following list (shown in Figure 7) summarizes the atmospheric correction algorithm developed in this study: 1.
For each pixel, the computed BLRs are matched to the BLRs from the calibration surface that minimize the Euclidean distance in BLR space.The corresponding water reflectance at 865 nm and 1016 nm is assigned to the pixel (see Section 2.7).6.
The atmospheric residual at these bands is obtained by subtracting the assigned water signal to the RC reflectance (see Section 2.8, Equation ( 16)). 7.
A final constraint is applied to ρ a (865) to limit a (865, 1016) inside the reasonable range of (0.85; 1.25) (see Section 2.8, Equation ( 17)).This approach may be extended to the whole range of bands of interest by subtracting the water signal retrieved and applying the clear pixel assumption to extrapolate the aerosol signal to shorter wavelengths (Gordon and Wang 1994 [1], Stumpf et al., 2003 [2]).

Baseline Residuals in Clear and Turbid Waters
The fact that BLRs of the considered triplets approach near-zero values in clear waters is expressed through Equation (10), assuming that very clear waters appear as black in the RNS.This property is clearly seen from OLCI RC reflectances, and is evident from Figure 8 (as well as Figures 2 and 5), where BLR distributions are intercompared from OLCI scenes from (i) windows selected to perform the BLR-AC calibration (which contain mainly turbid waters, calibration windows in Table 2), and (ii) windows selected to perform the Fix-AC scheme (which contain very clear waters, Fix-AC windows in Table 2).It is clearly seen from Figure 8 that BLR mean values and interquartile ranges (IQR) over the Fix-AC windows are close to 0 (always smaller than 10 −3 ) and the Fix-AC window IQR are at least 9 times smaller than the computed IQRs for the calibration windows.This is in accordance with what has been predicted and shows that for this approach the black water condition is translated to the condition BLR(ρ w ) = 0 (also evident by observing the clear water regions at Figure 5).This means that, given a pixel where simultaneously BLR(ρ w ) = 0 for the three considered triplets, the algorithm yields ρ w (865) = ρ w (1016) = 0, i.e., the algorithm returns naturally to the standard AC approach for clear waters.2) and clear water windows (blue, Fix-AC windows in Table 2), showing that BLRs are close to 0 in clear waters.This is also evident from Figure 9, where three different sources of the BLR vs. ρ w relation are intercompared to show how the baseline residuals at each of the three considered triplets of wavelengths vary according to water reflectances at 865 nm and 1016 nm.Blue dots show data from the calibration windows (Table 2); in magenta, BLRs computed from in situ data collected from Río de la Plata (Section 2.1), and a family of colour-mapped curves show the BLR vs. ρ w relation using a quasi-single scattering approximation (qSSA) reflectance model, described in Appendix A, where all specific IOPs, except specific particle absorption at a fixed wavelength (443 nm), were fixed in ranges consistent with previous publications (cf.Luo et al., 2018, Table 1 [7]) and where SPM varied logarithmically between 0.001 g/m 3 and 10,000 g/m 3 .Although there is an evident high resemblance between the three sources, it must be noticed that the qSSA curves depart from (remotely and in situ) measured water reflectances mainly because the underlying analytical reflectance model is less reliable at high reflectance (see Appendix A).The variations over the BLR vs. ρ w relation observed between the different qSSA curves indicate that BLRs might be very sensitive to particle absorption; which means that they might be used as indicative of different specific optical properties of particles in highly turbid waters.Nonetheless, this hypothesis must be ascertained with simultaneous reflectance and particle absorption field data, to test and validate optical closure, in the same way as expressed in Luo et al., 2018 [7].
In the three triplets considered, modelled BLR(ρ w ) behave similarly with increasing reflectance.At ρ w = 0, they all approach values towards 0, recovering the black water condition BLR(ρ w = 0) = 0 seen in Figure 8.As reflectances increase, BLRs take negative values, which can be understood if we consider ρ w ∝ b * b,p a w for very low turbidities (low sediment content) (cf.Ruddick et al., 2006, Equation (14) [5]).Given the simple qSSA reflectance model described in Appendix A, with specific IOPs taken from reported values in Babin et al., 2003 (a and b) [43,44], this magnitude is convex for the considered OLCI bands, i.e., negative BLRs (Figure 10a).On the contrary, for very high turbidities, water reflectances tend to achieve a constant value (saturation) that depends on the specific IOPs of the particulate content of the water: Dogliotti et al., 2015, Equation (A2) [16]).In this case, this magnitude is concave for the considered OLCI bands in the qSSA relfectance model, i.e., positive BLRs (Figure 10b).

Atmospheric Correction Performance
Given the lack of match-ups between in situ measurements and OLCI data for the studied regions (or any other extremely turbid waters), the Fix-AC approach (considered to be of high performance for the chosen scenes but not appropriate for a global automated processor) was used as a reference to validate the performance of the automated BLR-AC scheme (Figure 11).Pixels inside the marked coloured boxes, corresponding to the most turbid parts of the validation images listed in Table 2 were selected.The three regions were selected as representative of extremely turbid (Figure 11b) and moderately turbid (Figure 11a,c) waters.Water reflectances derived from both approaches show very similar patterns and very small root mean square difference (RMSD) values in both cases.In the case of the 620 nm band, ρ BLR−AC w values were obtained by using a simple linear aerosol extrapolation from the correction bands at 865 and 1016 nm.Another way of evaluating the performance of an atmospheric correction is through the analysis of the spatial (de)correlation between the retrieved water and aerosol reflectances.Figures 12-14 show water and aerosol reflectances at 1016 nm retrieved by BLR-AC from OLCI scenes over Río de la Plata (2017-01-21T13:24:42Z, 2016-06-08T13:09:51Z, 2017-12-12T12:59:01Z, Table 2), together with standard atmospheric corrections: ESA's BAC/BPAC and NASA's SeaDAS-2 (865, 1016) (see Section 2.3).The most challenging region for atmospheric correction over RdP is the turbidity front adjacent to Buenos Aires Province Coast (Argentina) at around 35°5 S; 57°0 W (marked as 0.8°× 0.8°black boxes), where SPM, turbidity and water reflectances at the RNS are maximum [8,24].It is evident from the higher slopes and R 2 values retrieved by linear regressions in Figure 12b,c how standard atmospheric corrections yield higher (unphysical) correlation between the water and aerosol signals.This is also evident by visually comparing the water and aerosol rasters in the BLR-AC case, compared to BAC/BPAC and SeaDAS-2 (865, 1016).In general terms, BAC/BPAC and SeaDAS-2 (865, 1016) underestimate water signal or simply do not converge to a numeric solution (i.e., NaNs, seen as magenta in Figures 12e,f,h,i, 13e,f,h,i and 14e,f,h,i).The image shown in Figure 12 is partially contaminated with moderate sunglint (toward the East edge), which is clearly removed as non-water signal for the BLR-AC and BAC/BPAC schemes.The image shown in Figure 14 is partially contaminated with thin clouds.BLR-AC is clearly capable of separating thin cloud from turbid water signal.

Discussion
The new atmospheric correction approach presented here has the advantage that it is designed on the basis of clear physical principles and thus highlights (e.g., Figure 6) how factors such as natural variability of specific inherent optical properties of water, TOA satellite data calibration and the choice of spectral bands may impact algorithm performance.BLRs computed from in situ measurements and OLCI data behave in a similar way and in accordance with a quasi-single scattering approximation (qSSA) reflectance model, although the qSSA model departs from measured reflectances at extreme turbidities, because of the underlying qSSA approximations (Figure 9).
While the new algorithm already gives a good performance for retrieval of Red/NIR/SWIR water reflectance, particularly in extremely turbid waters, and can thus be used directly for estimation of turbidity, SPM and chlorophyll-a concentration (in turbid waters using Red/NIR algorithms), further improvements may be necessary for retrieval of water reflectance at shorter wavelengths in turbid waters, if the latter is actually needed for subsequent water constituent retrieval algorithms.In extremely turbid waters the Red/NIR Rayleigh-corrected reflectance is dominated by water reflectance and retrieved aerosol reflectances and, in particular, the retrieved aerosol type (through a (865, 1016)) can be highly impacted by small errors in the water reflectance estimates-in fact, in this situation it is precisely the water reflectance that can be better estimated in relative terms than the aerosol reflectance.If more reliable a (865, 1016) are required in extremely turbid waters for extrapolation to shorter wavelengths then the current approach could be improved in the future by integration of extra information.For example, extra constraints could be applied to a (865, 1016) by the use of ultraviolet/violet bands (He et al., 2012 [45]) or by spatial smoothing of a (865, 1016) assuming that the horizontal length scale for variation of aerosol type is larger than the pixel size (Ruddick et al., 2000 [46]).In the present study no such horizontal smoothing has been applied, enabling us to better understand the capabilities and limitations of an automated pixel-by-pixel approach and hence identify, for example, the intrinsic value of the various spectral bands in separating aerosol and water reflectance in the RNS and the potential sources of uncertainty in this decomposition.In a future version of the algorithm, the constraint on a (865, 1016) applied at the last step (Section 2.9, step 7) of the processing might be replaced by either or both of these spectral/spatial constraints, and might address BLR optimization with a more refined method.
While the new algorithm relies on a BLR calibration surface (Section 2.9, step 5) constructed with Río de la Plata OLCI imagery, the AC scheme has been also tested in many other sites with sediment-dominated waters, like the Amazonas Plume, Bohai Sea, Belgian Coast, Bahía Blanca, (see Figure 11), using this calibration surface, so we currently consider that it will not be necessary to apply different calibrations to different regions.

Conclusions
A new atmospheric correction algorithm for highly turbid waters is presented for S3/OLCI.The key assumption is that the spectral convexity of Rayleigh-corrected reflectance in the RNS is determined by water reflectance and is essentially unaffected by aerosols.This spectral convexity is quantified here through baseline residuals (BLRs, Equation ( 5)) of the three consecutive triplets of the OLCI bands centered at 620, 709, 779, 865 and 1016 nm.An "equivalent transmittance" applicable to BLRs that depends on the air mass factor was derived from in situ measurements and radiative transfer simulations.The algorithm that relates BLRs to water reflectances at 865 and 1016 nm was calibrated by using OLCI data from selected subregions of images from Río de la Plata (Argentina), which were atmospherically-corrected by a simple fix window atmospheric correction (Fix-AC) that assumes spatial homogeneity of the atmosphere (given very clear sky conditions and small horizontal distances) and computes the atmospheric reflectance from fixed windows corresponding to clear water regions.The performance of the BLR-AC algorithm compares favourably with existing atmospheric correction algorithms (the standard ESA/OLCI processor and the SeaDAS processor implemented with the OLCI 865 nm and 1016 nm bands), particularly in extremely turbid waters, where existing algorithms either fail or show unphysical correlation between aerosol and water reflectance images (Figures 12-14).The performance improvement is achieved due to use of multiple red/NIR/SWIR bands, including the new OLCI 1016 nm band, which facilitates separation of water and aerosol reflectances even for extremely turbid waters.Although the BLR approach is designed here for OLCI imagery, it might be easily expandable to other sensors that have spectrally-close triplets of bands in the Red/NIR/SWIR, such as the Argentine-Brazilian joint mission SABIA-Mar, or hyperspectral sensors such as HICO or CHRIS/PROBA.

Figure 1 .
Figure 1.Location of the set of above-water reflectance measurements (red dots) taken from several field campaigns in the Río de la Plata (Argentina) in the period 2012-2015.

Figure 5 .
Figure 5. Example of a scene which forms part of the dataset used to calibrate the BLR(ρ w ) vs. ρ w relation.The magenta box surrounds the selected subset of the image added to the calibration dataset, while the aerosol signal was subtracted using the Fix Clear Window Atmospheric Correction Scheme (Fix-AC, Equation (14)) from clear water pixels (15 × 15 pixels window), indicated by the white boxes.

Figure 6 .
Figure 6.BLR(ρ w ) 3D sub-space, formed by the three linearly independent BLRs defined by the three consecutive triplets of the five OLCI bands 620, 709, 779, 865 and 1016 nm.Small dots: OLCI-A calibration dataset.Big colour-mapped dots: Calibration surface obtained from the OLCI-A calibration dataset, whose colour indicates water reflectance at 865 nm.Magenta dots are in situ data.The origin is indicated with an "X" and corresponds to "clear waters".

Figure 7 .
Figure 7. BLR atmospheric correction scheme, as listed in the summary.

Figure 10 .
Figure 10.Theoretical water reflectances at BLR bands yielded by the quasi-single scattering approximation (qSSA) at: (a) very low and (b) very high particulate concentration.

Figure 11 .
Figure 11.Intercomparison between BLR-AC and Fix-AC schemes.RGBs of scenes at Bahía Blanca (ARG), Río de la Plata and Belgian Coast from which the datasets are acquired (turquoise, violet and ocher yellow boxes at insets (a-c), resp.).ρ BLR−AC w

Table 1 .
Atmospheric, surface, geometric and spectral parameters used as input to the set of CNES-SOS simulations.

Table 2 .
List of Ocean and Land Colour Instrument (OLCI) scenes from Río de la Plata (L1B, processing baseline v2.23) used to calibrate and validate the algorithm, described by date of acquisition.