Atmospheric correction of Landsat-8/OLI and Sentinel-2/MSI data using iCOR algorithm: validation for coastal and inland waters

ABSTRACT Image correction for atmospheric effects (iCOR) is an atmospheric correction tool that can process satellite data collected over coastal, inland or transitional waters and land. The tool is adaptable with minimal effort to hyper- or multi-spectral radiometric sensors. By using a single atmospheric correction implementation for land and water, discontinuities in reflectance within one scene are reduced. iCOR derives aerosol optical thickness from the image and allows for adjacency correction, which is SIMilarity Environmental Correction (SIMEC) over water. This paper illustrates the performance of iCOR for Landsat-8 OLI and Sentinel-2 MSI data acquired over water. An intercomparison of water leaving reflectance between iCOR and Aerosol Robotic Network – Ocean Color provided a quantitative assessment of performance and produced coefficient of determination (R2) higher than 0.88 in all wavebands except the 865 nm band. For inland waters, the SIMEC adjacency correction improved results in the red-edge and near-infrared region in relation to optical in situ measurements collected during field campaigns.


Introduction
Remote sensing has been proven useful in many terrestrial and aquatic applications such as crop monitoring, mapping of invasive species or obtaining water quality information. However, deriving reliable information from remotely sensed data is challenging as the signal detected by the sensor is subject to atmospheric influences. The scattering and absorption of molecules and aerosols present in the atmosphere modify the pure reflectance originating from the target object. Many atmospheric correction algorithms have been developed to remove these unwanted effects, however multiple correction procedures are sensor specific and either developed for land (e.g. Conel, Green, Vane, Bruegge, & Alley, 1987;Hadjimitsis et al., 2010;Kruse, 1988) or water (Gao, Montes, Davis, & Goetz, 2009;Gordon & Wang, 1994;Ioccg, 2009;Mobley, Werdell, Franz, Ahmad, & Bailey, 2016). Common methods for ocean colour atmospheric correction do not look at elevation or adjacency effects and make assumptions on reflectance in the near-infrared (NIR) (e.g. Gordon & Wang, 1994) or short-wave infrared (SWIR) (e.g. Vanhellemont & Ruddick, 2015) wavelengths. These assumptions are not valid for land targets. Land atmospheric corrections often consider a Lambertian surface, while the air-water interface has a specular reflection (Gao et al., 2009). We present an alternative atmospheric correction tool image correction for atmospheric effects (iCOR) designed to work over inland, coastal or transitional waters and land.
iCOR depends on auxillary data for the atmospheric correction. Auxiliary data can originate from external sources or be derived from the image itself and include: digital elevation model (DEM), solar and viewing angles and the atmospheric composition. The satellite overpass time, sensor and sun position provide information on solar and viewing angles. The atmospheric composition is more difficult to retrieve and consists of molecular and aerosol scattering and absorption by gases, such as water vapor, ozone and oxygen. The aerosol contribution can be described by a combination of an aerosol model (e.g. Urban, Rural, Maritime or Desert Model (Carr, 2005)) and the aerosol optical thickness (AOT). The aerosol model describes the optical properties of the aerosols, i.d. single scattering albedo, phase function, extinction and scattering coefficients (Carr, 2005). The AOT indicates how much direct sunlight is prevented from reaching the ground by these aerosol particles (van Donkelaar et al., 2010). An important, but often neglected contributor to the TOA radiance is the adjacency effect (Richter, Bachmann, Dorigo, & Müller, 2006a;Vermote, Tanré, Deuzé, Herman, & Morcrette, 1997) which refers to light originating from neighbouring pixels and scattered into the target-sensor path. If adjacency effects are not properly considered, the overall apparent surface contrast decreases as bright pixels will be darkened and dark pixels brightened (Lyapustin & Kaufman, 2001). When atmospheric components are known or estimated, they can be imported into radiative transfer models (e.g. Moderate-Resolution Atmospheric Radiance and Transmittance Modelversion 5 (MODTRAN5) (Berk et al., 2006) or Second Simulation of a Satellite Signal in the Solar Spectrum (6S) (Vermote, Tanré, Deuzé, Herman, & Morcrette, 1994)) to retrieve the surface reflectance from the TOA radiance. iCOR uses MODTRAN5 Look-Up-Tables (LUTs) to remove atmospheric effects and the implemented AOT retrieval and adjacency correction over water are supported by image-based information. The next two paragraphs will discuss both topics more in detail.
Image-based AOT retrieval algorithms typically make assumptions about dark objects in the image. These dark objects are more sensitive to atmospheric effects as the atmospheric contribution to the overall signal is higher when compared to bright targets. Over land one of the approaches is the Dark Dense Vegetation method (DDV) (Kaufman & Sendra, 1988). This method begins with the assumption that vegetation is sufficiently dark and that the ratio between bottom-of-atmosphere reflectance at different wavelengths is constant (Ouaidrari & Vermote, 1999). This method is included in the SEN2COR atmospheric correction (Müller-Wilm, 2016). Over water, a black-pixel assumption (Gordon & Wang, 1994) is often used. This assumption expects that water absorption in the NIR wavelength region is so high that reflectance detected in this spectral region is only caused by atmospheric effects. For highly turbid waters, where this assumption is violated, one can switch to an adapted SWIR variant (Vanhellemont & Ruddick, 2015;Wang, Shi, & Tang, 2011). This black pixel approach for NIR or SWIR is implemented in ACOLITE . In iCOR, an adapted version of the Self Contained Atmospheric Parameters Estimation from MERIS data (SCAPE-M) method, developed by Guanter (2007), is implemented. The method exploits the spectral variability within a subset, making it less dependent on the presence of DDV and it does not pose any requirements on the water composition nor its spectral characteristics. SCAPE-M has been proven successful over land on MERIS data (L. Guanter, Del Carmen González-Sanpedro, & Moreno, 2007), Chris/Proba hyperspectral data (Guanter, Alonso, & Moreno, 2005b) and for the correction of inland waters (Guanter et al., 2010). A known shortcoming of SCAPE-M is the inability to retrieve information on the aerosol type over land (Ramon & Santer, 2005;Santer, Zagolski, & Gilson, 2007).
Adjacency effects can be corrected by including spectral information on the environment of the target in the calculations. However, the extent of the environment that influences the target observation is often unknown. Alternatively, a fixed window range can be considered which is valid for a horizontally homogeneous underlying surface and produces good results for land targets (Minomura, Kuze, & Takeuchi, 2001;Richter, Schläpfer, & Müller, 2006b). In coastal or inland waters, this approach is less favourable as the sliding average will vary from typical land to typical water value (Kiselev, Bulgarelli, & Heege, 2015). Alternative methods developed to tackle the influence of adjacent light include improve contrast between ocean & land (ICOL) (Santer et al., 2007), C-Wombat-C (Brando & Dekker, 2003) or methods based on the point spread function (Kiselev et al., 2015). The adjacency correction options available in iCOR are a fixed window range for land and the SIMilarity Environmental Correction (SIMEC) method developed by Sterckx, Knaeps, Kratzer, and Ruddick (2014) for water. SIMEC estimates the background contribution iteratively by checking against the NIR similarity spectrum (Ruddick, De Cauwer, Park, & Moore, 2006). This paper describes the methodology of iCOR, previously known as OPERA , and shows preliminary results for the Operational Land Imager (OLI) on board of Landsat-8 (L8) and the MultiSpectral Imager (MSI) onboard of Sentinel-2 (S2). Although iCOR works over land and water, the primary focus will be the validation of coastal and inland waters. As indicated by Palmer, Kutser, and Hunter (2015) the inland water community often has to choose between the coarse spatial resolution of ocean colour sensors or the sub-optimal spectral characteristics and radiometric sensitivity of land missions. L8 and S2 are both land missions, but different studies (e.g. Pahlevan & Schott, 2013; show their applicability for inland and coastal water quality monitoring thanks to their high spatial resolution and high signal-to-noise ratio (SNR). For validation data from the Aerosol Robotic Network -Ocean Color (AERONET-OC) of the Belgian coast and optical in situ measurements collected from field campaigns in multiple inland lakes (and a lagoon) across Europe are used. Specific attention will be placed on adjacency effects in inland lakes.

Icor workflow
The overall workflow of iCOR is presented in Figure 1. This schematic can be summarised into the following steps, (i) land and water pixels are distinguished, (ii) land pixels are used to derive the AOT based on an adapted version of the method developed by Guanter (2007), (iii) an adjacency correction is performed using SIMEC (Sterckx et al., 2014) over water and fixed background ranges over land targets (Berk et al., 2006), and (iv) the radiative transfer equation is solved. iCOR uses MODTRAN5 (Berk et al., 2006) LUT to perform the atmospheric correction and uses information about the solar and viewing angles (Sun zenith angle (SZA), view zenith angle (VZA) and relative azimuth angle (RAA)) and a DEM. The next paragraphs focus on different components of the iCOR, first a small introduction on the radiative transfer function is provided, next the implemented AOT retrieval technique and the adjacency correction algorithms are described.

Radiative transfer and atmospheric correction
The radiance received by the sensor (ðL rs target Þ consists of atmospheric path radiance L atmÀpath (at sensor level), background path radiance L backgrÀpath and ground reflected radiance L target (Kaufman, 1984): With • L atmÀpath photon scattered into the sensor's instantaneous field-of-view, without having ground contact. • L backgrÀpath the reflected radiation from the neighbourhood of the target pixel and diffusely scattered into the field-of-view. This term is responsible for the adjacency effect. • L target the ground reflected radiance. Over water, this component consists of the transmitted water leaving radiance and the transmitted sky glint radiance (scattered light reflecting from the surface). Sun glint is ignored.
In the atmospheric correction process, the signal is corrected for the first two terms and the surface reflectance is then retrieved from L target . The three radiance components can be written as: where F 0 is the extra-terrestrial solar irradiance, R atm θ v ; θ s ; ϕ v À ϕ s À Á is a coefficient describing the reflection of light by the atmosphere, θ s and ϕ s are respectively the sun zenith and azimuth angle, θ v and ϕ v are respectively the viewing zenith and azimuth angle, E d is downwelling irradiance just above the surface, t dif and t dir are respectively the diffuse and direct ground-tosensor transmittance, τ is the optical depth in the atmosphere, ρ s is the surface reflectance and ρ backgr is the background reflectance, defined in a similar way as ρ s . The downwelling irradiance above the surface E d can be written as: Þ) and s Ã the spherical albedo of the atmosphere. Sterckx, Knaeps, and Ruddick (2011) described how to derive, starting from these equations, surface reflectance ρ s , and derived the following equation: Over water surfaces, an extra correction for reflected sky glint (d 1 Þ results in the water leaving reflectance ρ w : where • L sky L sky the downwelling sky radiance The parameters (c 1...5 ; d 1 ) depend on a large number of geometric and atmospheric parameters (solar and viewing zenith angle, ground altitude, aerosol density, water content and aerosol type) which vary over the scenes. A MODTRAN5 interrogation technique, as described in detail by De Haan and Kokke (1996), can be used to derive values for these parameters. To minimize computing time, a MODTRAN5 LUT approach is used as described in Biesemans et al. (2007). This approach assumes that the pixel is not affected by direct sun glint and that wind speed is low so that foam or white cap radiance can be ignored. Light is considered to be unpolarised and the sky radiance to be distributed in a uniform way.

Aerosol optical thickness
One of the initially unknown input parameters of the radiative transfer equation is the AOT. Various methods exist for estimating AOT, either using external sources (e.g. AERONET stations) or image-based. In iCOR we implemented an adapted version of the land based AOT retrieval technique described by Guanter et al. (2005b). The flow chart is included in Figure 1. This AOT retrieval algorithm makes use of the spectral variability of the land pixels within an image.
In a preliminary step cloud and water pixels are masked using a simple single band threshold for water and multiple thresholding levels defined for clouds (based on Guanter, Alonso, Moreno, and Member (2005a).
• The average top-of-atmosphere (TOA) reflectance for all visible and near-infrared (VNIR) bands is calculated and compared with a predefined "average" threshold. If this threshold is exceeded and the reflectance in the blue band is higher than a minimum threshold, the pixel is considered as cloud. • For sensors with a cirrus band, an additional threshold is set on this band for the detection of cirrus clouds.
The cloud mask is dilated with an extra surrounding border to ensure that pixels affected by clouds but undetected are not considered in the AOT retrieval phase. The raw TOA image is subdivided into tiles of about 15 × 15 km which are small enough to assume atmospheric homogeneity and large enough to include high spectral variation (Guanter et al., 2010). For each tile, the lowest radiance value within the tile is selected for each spectral band and the corresponding path radiance of this approximated dark target spectrum is retrieved using the pre-calculated MODTRAN5 LUT. The AOT value leading to the path radiance closest to the dark spectrum becomes the tile upper AOT boundary, preventing path radiance to be higher than the dark spectrum in any of the spectral bands. In the next step, the initial AOT estimation is refined through a multiparameter endmember inversion technique. Five pixels with high spectral contrast (selected based on the NDVI values from TOA reflectance) are represented by a linear combination of two predefined default vegetation spectra and a soil spectrum (endmembers) to estimate the surface reflectance: With ρ s , ρ veg and ρ soil the surface reflectances of respectively the reference pixels and predefined vegetation and soil spectra. The C v;s parameters are left free in the inversion. This results in an 11-D inversion, with two parameters for every five pixels and AOT as degrees of freedom. The multi-parameter endmember inversion is performed though the minimisation of the Merit function : L SIM is simulated TOA radiance, stored in and retrieved from MODTRAN5 LUT, L SENS is the measured TOA radiance, λ i is the center wavelength of the ii-th band and ω pix is the weighting factor, which is 2.0 for pure vegetation, 1.5 for mixed and 1.0 for pure soil pixels to enhance the sensitivity in vegetation targets to aerosol loading. The function is weighted by λ À2 i to drive the inversion towards the smaller wavelengths, where the aerosol effect is much larger than the reflectance of most natural surfaces. Minimisation of the merit function is done by the Powell's minimisation method (Press, Flannery, Teukolosky, & Vetterling, 1986).
The last step in the AOT retrieval scheme is the interpolation of missing pixels and the smoothing of the resulting mosaic. The missing data cells, deselected due to cloudiness or the land/water mask, are interpolated from neighbouring cells. To scale from cell image to per-pixel image, the cubic convolution method is used.
There are two important restrictions of this approach. Firstly, surface reflectance should be representable by a linear combination of two pure green vegetation and a bare soil endmember. Open oceans, deserts or snow landscapes do not meet this requirement. The second limitation is the difficulty in deriving a reliable estimate of aerosol model over land (Ramon and Santer (2005), Grey, North, and Los (2006)). iCOR considers a fixed rural aerosol model as default, similar as proposed by Guanter et al. (2005b) for the CHRIS/PROBA mission.

Adjacency correction
For land targets, adjacency correction assumes that environmental influences can be approximated by averaging (weighted average) the spectra of the neighbouring pixels over a fixed distance, which can be chosen arbitrarily. This weighted average spectrum can then be inserted into the radiative transfer equation as background radiance (L rs backgr in Equation 6). For water targets, the SIMEC adjacency correction (Sterckx et al., 2014), originally developed for hyperspectral airborne imagery (Sterckx et al., 2011), is implemented. This starts from the NIR similarity (NIR sim ) assumption (Ruddick et al., 2006) which states that the shape of the water spectrum in the NIR region is invariant. After normalizing the water leaving reflectance at 780 nm, the value should fall within a predefined range: With ρ t w and ρ r w the retrieved water leaving reflectance for respectively a "test" spectral band and the "reference" spectral band situated near 780 nm. When this requirement is not fulfilled, pixels are assumed to be influenced by adjacency effects. The background radiance (L rs backgr in equation 6) for range N is calculated as a weighted average of the pixel radiance values surrounding the target pixel. In an iterative manner, the optimal range (N) for defining the environmental influences is determined as shown in Figure 1. Sterckx et al. (2011) defined optimal waveband settings for the "test" band as; (i) minimum influence of gaseous absorption, (ii) not lower than 690 nm because of increasing uncertainty, (iii) preferably located in the red-edge region of the spectrum where contrast between water and land is more pronounced. The availability of spectral bands is a restrictive factor in multispectral sensors. The band setting of Landsat-8 is not ideal because spectral bands at 655 nm and 865 nm are selected respectively as "test" and "reference". Sentinel-2 allows a better spectral band selection with 705 nm and 783 nm selected respectively. Sterckx et al. (2014) defined restrictions for the use of SIMEC which tends to fail in (i) high turbid waters where the NIR reflectance is flattened (Doron, Bélanger, Doxaran, & Babin, 2011;Goyens, Jamet, & Ruddick, 2013), (ii) in waters with macrophyte growth or specific algae blooms and in (iii) areas where bottom effects are significant in the NIR (optically shallow waters).

Coastal waters
The performance of iCOR in coastal water was validated by comparing the atmospherically corrected data with AERONET-OC data. AERONET is developed by the National Aeronautics and Space Administration (NASA) to sustain atmospheric studies at various scales with worldwide autonomous sun-photometer measurements (Zibordi et al., 2009). This network has been extended in AERONET-OC to support marine applications, by providing additional capability of measuring the radiance emerging from the sea (i.e. water-leaving radiance) with modified sun-photometers installed on offshore platforms. AERONET-OC is instrumental in satellite ocean colour validation activities through standardized measurements (i) performed at different sites with a single measuring system and protocol, (ii) calibrated with an identical reference source and method, and (iii) processed with the same code (Zibordi et al., 2009).
Two AERONET-OC stations situated in the North Sea near the Belgian coast were selected for validation; Zeebrugge-MOW1, located in coastal nearshore (3.65 km from land) waters (51.36°N, 3.12°E), and Thornton_C-Power, situated further offshore (26.25 km from land) in clearer waters (51.53°N, 2.96°E), see Figure 2. The AERONET-OC stations have a SeaWiFS Photometer Revision for Incident Surface Measurements (SeaPrism) installed, which autonomously performs multiple sky-and sea-radiance observations at programmable viewing and azimuth angles at eight (nine in the 2006 instrument release) center wavelengths in the 412-1020 nm spectral range (Zibordi et al., 2006).
Level 2.0 Quality Assured, or in absence Level 1.5 Real Time Cloud Screened, Normalized Water Leaving Radiance (Lwn) from AERONET-OC was downloaded from the AERONET website (http://aero net.gsfc.nasa.gov/). The maximum allowed time difference between AERONET-OC measurements and the image acquisition time was set at 1 h, since water can be highly dynamic in space and time . If more than one AERONET-OC measurement fulfilled this requirement, the data were interpolated to the corresponding satellite overpass time. The validated Landsat-8 and Sentinel-2 scenes are listed in Table 1. The Lwn data [mW/(cm 2 sr um)] from AERONET-OC was converted to water leaving reflectance (ρ w ) through following formula: With F0 λ ð Þ the exo-atmospheric solar irradiance [mW/(cm 2 sr um] at wavelength λ ð Þ from Thuillier et al. (2003). The differences in the spectral response curves between Landsat-8, Sentinel-2 and AERONET-OC were disregarded in this analysis.
From the satellite imagery, the median ρ w value of a 300×300 m window surrounding the AERONET-OC station was used to exclude noise coming from the station or its shadow. For Landsat-8, this resulted in an 11 × 11 window for the five bands in the VNIR. For Sentinel-2 the window was 5 × 5 pixels for B1 (443nm at 60m), 15 × 15 pixels for B5 (865 nm at 20m), and 30 × 30 pixels for B2, B3 and B4 (resp. 490, 560 and 665nm at 10m). Since the AERONET-OC stations are located at a certain distance from the coast, the effect of the surrounding land on the target pixels is assumed to be negligible and no SIMEC corrected was performed.

Inland and transitional waters
In situ data for validation of Landsat-8 and Sentinel-2 imagery were gathered during multiple field campaigns (2014-2016) in the following five European lakes: Balaton, Marken, Mantua, Garda and Geneva. Moreover field data from Curonian Lagoon was also included in the validation exercise. The location of the lakes is shown in Figure 3. An overview of the test sites is included in this paper, a full detailed description of the field campaigns can be found in Reusen et al. (2014) and Reusen et al. (2016); • Lake Balaton is Europe's largest shallow lake at 592 km 2 , located in Hungary. In spite of its large surface area, it is very shallow with a mean depth of approximately 3.2 m. The lake has historically been subdivided into four basins (west to east): Keszthely, Szigliget, Szemes and Siófok. The main inflow to the lake is the River Zala in the western Keszthely basin and the only outflow is via a regulated canal in the eastern Siófok basin. The trophic state of Lake Balaton varies between mesotrophic and eutrophic regimes. • Lake Marken (Markermeer) is a 700 km 2 large, artificial, shallow, eutrophic lake in The Netherlands. It is part of a former brackish, 4000 km 2 inland sea (Zuiderzee) that was dammed and turned into Lake IJssel in 1932. After the construction of the Houtribdijk dam between the cities Lelystad and Enkhuizen, Lake Marken has been separated by sluices from the major inputs of River IJssel.
• The lakes of Mantua are composed by three small and shallow eutrophic fluvial lakes located in northern Italy. The Mantua Lakes have features typical of a shallow lentic environment due to low water velocity and limited depth, and of a wetland hosting dense macrophyte meadows, and therefore can be defined as a fluvial lake system. • Lake Garda, located in the subalpine ecoregion of Italy has an area of 370 km 2 , a water volume of 50 km 3 and a maximum depth of 346 m. From an ecological point of view, the basin can be categorised as a warm, monomictic and oligomictic basin, classified as oligo-mesotropic. • Lake Geneva is one of the largest lakes in Western Europe with a surface area of approximately 580 km 2 . The lake is extremely  deep, with a mean depth of 152 m and a maximum depth of 310 m. The main inflow into the lake is the Rhone River in the east. The lake is divided into three basins: The Upper Lake (Haut Lac) in the east receives water from the surrounding Alps via the Rhone River and as such often exhibits relatively high concentrations of mineral particles; the Large Lake (Grand Lac) is the largest and deepest basin, while the Small Lake (Petit Lac) in the south-west is shallower. Historically the lake has suffered from eutrophication, but reductions in the nutrient load have returned the lake towards an oligo-mesotrophic statealthough phytoplankton blooms can occur on the lake during summer. • Curonian Lagoon, which is the largest in Europe, is a shallow water body (total area 1584 km 2 , mean depth 3.8 m) located along the south eastern coast of the Baltic Sea and geographically positioned between Lithuania and the Russian Federation. The waters are considered eutrophic or hyper-eutrophic.
Measurements of downwelling irradiance ðE d Þ, skylight radiance (L sky ) and total upwelling radiance from the water (L u ) for computation of ρ w were made using various optical radiometric devices; ASD FieldSpec FR, Satlantic HyperSAS, Trios RAMSES, Spectral Evolution or WISP-3 (Hommerson et al., 2012) and used as input in following formula (Mobley, 1999): with r as the air-water interface reflection coefficient. The HyperSAS and RAMSES data were processed and quality controlled using the "fingerprint" method described in Simis and Olsson (2013). In situ ρ w data, collected during Landsat-8 and Sentinel-2 overpasses were compared with iCOR derived ρ w values through scatterplots, boxplots and correlation coefficients. The list of Landsat-8 images is shown in Table 2. In total 56 matchups with Landsat-8 and 27 matchup points of Sentinel-2 were used for validation, from which the hyperspectral in situ reflectance data were resampled to fit the Landsat-8 and Sentinel-2 spectral response curves. The mean time difference between in situ measurements and image acquisition was less then 2 h.
Since inland waters can be affected by adjacency effects (e.g. Odermatt, Kiselev, Heege, Kneubühler, & Itten, 2008) two iCOR runs were performed with and without SIMEC adjacency correction, which will be referred to as the SIMEC and BASIC run respectively. In contrast to the coastal water validation, where the median value in a window of 300×300 m was selected, the inland water validation considered the average ρ w value between the target pixel and its immediate surrounding pixels. This is to reduce the chance of including land pixels in the analysis. For Landsat-8 the mean value of a 3 × 3 pixel window was selected. For Sentinel-2, a 5 × 5 pixel window was selected for the 10 m bands, a 3 × 3 pixel window for the 20 m bands and only one pixel for the 60 m band. Figure 4 shows the intercomparison between iCOR and AERONET-OC data for different bands in the VNIR regions. The four visible (VIS) bands show reasonable results with a slope varying from 0.99 to 1.17, an R 2 between 0.88 and 0.98, the root mean square error (RMSE) between 8.5e-3 and 1.4e-2 and the mean absolute percentage error (MAPE) between 13% and 65%. Band 560 nm produced the best results. Both L8 and S2 follow the same trend in these bands. The performance of the NIR band is less accurate with a slope of 2.25, and R 2 of 0.54 and RMSE of 8.9e-3. iCOR tends to overcorrect the ρ w value in this band, which is Table 2. List of inland lakes with the selected Landsat-8 and Sentinel-2 images and the corresponding date and number (N°) of matchups during which in-situ measurements were taken. Scene Info contains information on the Path and Row of the Landsat-8 and on Orbit and Tile of the Sentinel-2 image acquisition. particularly evident at CPower station where AERONET-OC data is always lower than 4.5e-4 while iCOR ranges to 1.2e-2. The different spectral band characteristics of the sensors (neglected in this study) might explain the poorer performance at 865 nm where water absorption is high.

Inland and transitional waters
Scatterplots of iCOR ρ w values and in situ collected data are shown in Figure 5 for L8 data and Figure 6 for S2 data which presents the results for the BASIC and SIMEC runs. For Landsat-8, good agreement is demonstrated from the intercomparison exercise, especially in the VIS region where: the slope varies between 1.07 and 1.17, the R 2 value is higher than 0.87 and MAPE values are lower than 77%. The results of the 865 nm band are poorer with a slope of 2.21, R 2 of 0.6 and MAPE of 836%. Correcting for adjacency effects with SIMEC has a slightly positive effect in the VIS regions (R 2 and MAPE values remained or improved) and caused a significant amelioration in the 865 nm band: the slope dropped to 1.54, R 2 increased to 0.83 and MAPE dropped to 237%. This means that, despite the non-ideal bandsettting for Landsat-8, running SIMEC in inland waters can improve the retrieved reflectance signal.
For Sentinel-2 promising results are also retrieved in the VIS region where the slope varies between 0.8 and 0.88, R 2 is higher than 0.7 and the MAPE values are lower than 100% and with exclusion of the aerosol band (443 nm) which was lower   than 40%. In the NIR region, the performance of the BASIC iCOR run is less satisfying with slope values high (1.31-1.7), the R 2 varying between 0.5 and 0.63 and MAPE values exceeding 100% (even 386% at 865 nm). Applying SIMEC lowered the slope in all bands. In the VIS bands this resulted is a slight deterioration of the performance with a drop in R 2 of maximum 0.1 (560 nm), RMSE increase of 0.002 and MAPE increase of up to 8%. However, from 705 nm onwards SIMEC improved the results significantly. The largest effect was noted in band 783 nm where R 2 increased from 0.63 to 0.75, the RMSE dropped from 9.6e-3 to 3.3e-3 and the MAPE value decreased from 248% to 47%. Figure 7 shows the subtraction of in situ collected data from iCOR retrieved ρ w data in boxplots. For Landsat-8, the bias is slightly lower than zero for all bands, both for the BASIC and SIMEC run. For Sentinel-2 this average approaches zero for the VIS bands in the BASIC run and are negative for the NIR bands. SIMEC seems to increase the average value, resulting in positive values in the VIS and close to zero values in the NIR. An overcorrection by SIMEC in the VIS might be related to a slight overcorrection of the AOT or caused by the affected spectral characterisation of S2A e.g. bands B01 and B02 have been identified as being affected by a measurement defect (ESA, 2017). To better understand the effect of SIMEC on the overall spectral shape of a water pixel, Pearson Correlation coefficients between in situ and iCOR ρ w values are plotted in boxplot format in Figure 8. Here it is clear that SIMEC has a positive effect on the spectral shape for both Landsat-8 and Sentinel-2 as the median Pearson correction value after SIMEC is closer to 1 and also the outliers have a higher correlation value.
The spatial extent of adjacency effects is illustrated in Figure 9 for a Landsat-8 image of Lake Mantua (23/09/2014) and in Figure 10 for a Sentinel-2 image of Lake Marken (15/09/2016). In these figures the iCOR retrieved ρ w values at 865 nm are displayed with a colour legend. From each image three match up points are selected and their spectra (in situ, BASIC and SIMEC) are plotted in Figures 11 and  12 for Mantua Lakes and Lake Marken respectively. Figure 9 shows that the reflectance in NIR without SIMEC correction is higher and the sediment plume (slightly east of Point 3) becomes more blurred. Figure 11 illustrates that the spectral shape of ρ w is reasonably in line with the in situ results after performing SIMEC. For Point 2 the reflectance is slightly overcorrected in the VIS. In Lake Marken, adjacency effects are especially pronounced at the border of the lake near West side.what has the highest presence of vegetation. After performing SIMEC, border effects are reduced. The spectra shown in Figure 12 indicate that SIMEC had no impact in the VIS region, but improved the results in the NIR.

Discussion
iCOR is developed to be applicable to land and inland, transitional or coastal water scenes. The advantages and limitations of iCOR are linked to the different steps in the process.
A MODTRAN5 LUT approach solves the radiative transfer equation and takes terrain elevation and sky glint into account. However, the current version of iCOR does not correct for sun glint effects. This is an important limitation, because the viewing angles of the OLI and MSI sensor make images vulnerable to Figure 9. Water leaving reflectance at 865 nm for the Landsat-8 image of Lake Mantua (23/09/2014) without SIMEC adjacency correction (a) and with SIMEC adjacency correction (b). Three locations of insitu sampling are highlighted from which the water leaving reflectance spectra are plotted in Figure 11. © U.S. Geological Survey Landsat data, processed by VITO. Figure 10. Water leaving reflectance at 842 nm for the Sentinel-2 image of Lake Marken (15/09/2016) without SIMEC adjacency correction (a) and with SIMEC adjacency correction (b). Three locations of in-situ sampling are highlighted from which the water leaving reflectance spectra are plotted in Figure 12. © Copernicus Sentinel data processed by VITO.
sun glint contamination (Harmel, Chami, Tormos, Reynaud, & Danis, 2018). As pointed out by Kay, Hedley, and Lavender (2009) methods using statistical sea surface models, commonly used for processing of ocean colour images (e.g. Cox & Munk, 1954), are not suitable for higher resolution imaging sensors and an image-based solution is required. Removal or correction of sun glint effects will be considered in a next phase of iCOR.
With the implementation of a land-based AOT retrieval, no assumptions are made on the spectral shape of the water signal. This is a significant advantage for inland, coastal and transitional waters, since the constituents of the water are often not known and there is no dependency on SWIR bands in highly turbid waters (Knaeps, Dogliotti, Raymaekers, Ruddick, & Sterckx, 2012) which have lower SNR, especially for terrestrial missions such as L8 and S2. The main disadvantage of the current methodology is that land pixels should be present and should be representable by a linear combination of green vegetation and bare soil. When this is not the case, the user can set a default AOT value. The validation results presented in this paper show that the proximity of land allows using a continental rural aerosol model which is set as default since it is difficult to retrieve the aerosol type over land pixels (Grey et al., 2006;Ramon & Santer, 2005). In a next phase of iCOR, we will include and validate water based AOT retrieval in combination with the current land based implementation. It is expected that this will reduce errors caused by extrapolation of AOT over large water bodies.
iCOR allows the user to perform an adjacency correction over land and water. For water the SIMEC method (Sterckx et al., 2014) is implemented, which has been previously validated for MERIS data. In this paper SIMEC was validated on L8 and S2 data, and although the bandsetting of L8 was not ideal, adjacency effects were reduced after running SIMEC. But, as pointed out by Sterckx et al. (2014) some caution is needed when running SIMEC in high turbid waters, in waters with macrophyte growth or specific algae blooms or in areas where bottom effects are significant in the NIR (optically shallow waters).
In general, the assessment of iCOR for inland and coastal waters shows reasonable results for L8 and S2.
The retrieved values are in line with ACOLITE results for the coastal AERONET-OC stations (Van der Zande, Vanhellemont, De Keukelaere, Knaeps, & Ruddick, 2016). iCOR results over land were tested and validated in the Atmospheric Correction Intercomparison Exercise (ACIX) organised by ESA-NASA (Doxani et al., 2018).

Conclusion
The iCOR atmospheric correction workflow is MODTRAN5 based, written in C++, and includes (i) land/water masking, (ii) land-based AOT retrieval, (iii) adjacency correction for both land (fixed ranges) and water (SIMEC) and (iv) final atmospheric correction. In this work, the iCOR performance has been tested using L8 and S2 imagery for inland, transitional and coastal waters. The validation exercise was two-fold; (i) a comparison was performed with AERONET-OC data over coastal water and (ii) with optical in situ data collected during dedicated field campaigns in multiple inland lakes and a lagoon around Europe. In the latter case, additional attention was paid to adjacency effects and the role of SIMEC.
Imagery acquired over the Belgian North Sea had been compared with the AERONET-OC data of Zeebrugge_MOW1, located in coastal nearshore waters, and Thornton_C-Power, further offshore in clearer waters. In these coastal waters iCOR performed well, despite the use of a fixed rural aerosol model in the atmospheric correction instead of a marine aerosol model. This was particularly true in the VIS region with an R 2 value higher than 0.88. The NIR band (865 nm) showed less accurate performance with high MAPE values and lower R 2 values. Similar trends were observed in inland lakes. The SIMEC adjacency correction showed a negligible or slightly negative effect on the ρ w values in the VIS range (R 2 drop lower than 0.03 and RMSE increase of max 0.002), but a significant improvement in the NIR region. Here the R 2 increased with 0.23 for the 865 nm Landsat-8 band and with 0.12 for the 783 nm Sentinel-2 band. As a consequence, the overall spectral shape of the water pixels improved after running SIMEC. Also for the inland waters, band 865 nm showed the poorest results.
The iCOR atmospheric correction for L8 and S2 is freely available for users as a plugin in the SNAP Sentinel toolbox. iCOR is user-friendly and requires minimal manual interaction as the tool attempts to retrieve all required information from the image itself. Neverthess, the user is free to select various options, such as cloud or water detecting settings, or enabling the correction of adjacency effects over land and/or water.