Performance of the Landsat 8 Provisional Aquatic Reflectance Product for Inland Waters

Recently, the United States Geological Service (USGS) released a new provisional product which estimates aquatic reflectance from Landsat 8 Operational Land Imager (OLI), called Landsat 8 Provisional Aquatic Reflectance (L8PAR). However, as indicated in the product guide, the use of this product for inland waters needs further verification and improvements. The goal of this study was to determine how the novel L8PAR product performs for different small turbid and eutrophic lakes in Northern Germany compared to in situ measurements of above water remote sensing reflectance (Rrs). For several recent scenes during our monitoring, the L8PAR product failed to produce full data for the lakes of our interest. For the best scene with in situ spectra, L8PAR was not able to retrieve any information for band 1 and not all information for bands 2, 3 and 4. The pixels with valid values for reflectance showed a weak relationship for band 2 (R2 of 0.24) and a medium relationship for bands 3 and 4 (R2 of 0.68 and 0.72, respectively). Compared to other atmospheric correction routines (ACOLITE, C2RCC, C2X, iCOR and L8SR), L8PAR was the only product which was not able to retrieve Rrs for all match up samples. This work provides an evaluation of the L8PAR product for inland waterbodies. Although more analysis and validation need to be conducted, our study suggests that the L8PAR product cannot be used for small inland lakes in its current state and has to be used with care for inland waters in general.


Introduction
The Operational Land Imager (OLI) onboard Landsat 8 is the most recent optical sensor from the Landsat satellite family [1]. Although it was intended primarily for monitoring land use/land cover, OLI images have been used for lake color studies for the estimation of water transparency, colored dissolved organic matter [2], identification of algal blooms [3], total suspended sediments (TSS) [4], chlorophyll-a [5], Secchi disk depth [6] and turbidy [7]. The major advantage of OLI over existing ocean color missions for inland waters is related to its spatial resolution of 30 m, which provides the capability to map optically active components (OACs) in small-sized aquatic systems (<5 km 2 ).
However, in aquatic environments the above-water remote sensing reflectance (R rs ) can be on the order of 1% at the peak of reflectance [8]. Because of the low intensity of the signal from aquatic Remote Sens. 2020, 12, 2410 2 of 12 environments, atmospheric correction routines for R rs retrieval require a very high accuracy [9]. For highly heterogeneous inland waters, this task is more complicated due to additional effects on the retrieval of R rs , such as water turbidity which can increase the signal in the near-infrared spectrum, whitecaps, sun-glint and adjacency effects from the surrounding land [10].
On 1 April 2020, the USGS/ESPA launched the Landsat 8 Provisional Aquatic Reflectance (L8PAR) Product [18]. This provisional scientific product for the retrieval of aquatic reflectance based on the estimation of R rs was validated using the coastal waters spectra from the Ocean Color Aerosol Robotic Network (OC-AERONET) [19,20]. Thus, the performance evaluation of the provisional L8PAR product over highly turbid and eutrophic inland waterbodies is still needed [18].
The goal of this study is to provide a radiometric and flags evaluation of the applicability of the L8PAR product for the retrieval of R rs for small inland waters (surface areas ranging from 1.05 to 2.58 km 2 ). The analysis compares the retrieved spectra from the L8PAR product to the in situ R rs measurements. Additionally, the L8PAR product was further evaluated by comparing it with the performance of other atmospheric correction routines using the same dataset. The evaluation of each atmospheric correction routine was computed via statistical error metrics in comparison to these in situ observations.

Satellite Image and Products, L8PAR and L8SR
A Landsat 8/OLI image from 26 July 2019 (Product ID: LC08_L1TP_193023_20190726_20190801_01_T1) was compared to the respective field data. The Landsat 8/OLI Level 1 product (LC08DN, Figure 1A) was downloaded from the USGS Earth Explorer website (earthexplorer.usgs.gov) and was used for the computation of the atmospheric correction routines (ACOLITE, C2RCC, C2X and iCOR). L8SR and L8PAR products were ordered via the USGS EROS ESPA on-demand interface (https://espa.cr.usgs.gov/). Version 1.0 of the L8PAR was used to retrieve R rs values from a set of inland waters. This product computes the aquatic reflectance which was divided by π for the computation of R rs . The product is based on the R rs Algorithm Theoretical Basis Document (ATBD) developed by the NASA's Ocean Biology Processing Group (OBPG) [19,20]. More information about this product can be found in the L8PAR Algorithm Description Document [21] and in Pahlevan et al. [19].
The L8SR product is generated using the Land Surface Reflectance Code (LaSRC) algorithm (Version 1.4.1), which uses band 1 (centered at 443 nm) to perform aerosol inversion tests, auxiliary climate data from the Moderate-Resolution Imaging Spectroradiometer (MODIS), and a radiative transfer model [17]. This product was developed for the retrieval of reflectance for land applications. However, we included it here for a comparison of performances among land and aquatic atmospheric correction routines to showcase the differences between both methodologies: one that considered the reflected skylight and one that does not. Therefore, the surface reflectance values from the L8SR product were divided by π (in this case the atmospheric correction routine considered that the upwelling radiance is the same as water leaving radiance and there is no influence of the reflected skylight).

Product Evaluation
For the evaluation of atmosphericly corrected products, statistical metrics were computed to assess the accuracy of the retrieved Rrs signals. Table 1 summarizes the list of error metrics used here.

Atmospheric Correction Routines
R rs retrieval from L8PAR and L8SR products was compared to the results of different atmospheric correction routines. ACOLITE is an atmospheric correction routine that uses the dark spectrum fitting approach by default, however, it can be configured to use the exponential extrapolation approach [11,12]. In this study, ACOLITE version Python 20190326 was used with its default configuration. C2RCC and C2X [13][14][15] rely on a database of radiative transfer simulations of water-leaving radiance and related Top-Of-Atmosphere (TOA) radiances via neural networks. The difference between them is that C2RCC is designed for low-to-medium concentration of OACs and typical ranges of inherent optical properties (IOPs), while C2X was designed for high concentration of OACs and extreme ranges of IOPs [13,15]. The Landsat 8/OLI image was processed according to the default configuration via the external tool available (version 1.2) on the Sentinel Application Platform (SNAP, version 7.0). iCOR uses MODTRAN 5 Look Up Tables to perform the atmospheric correction and needs information about the solar and viewing angles (Sun Zenith Angle, View Zenith Angle and Relative Azimuth Angle) and a digital elevation model [16]. The level 1 image was processed according to default processing parameters via the external tool available for SNAP (v. 7.0). For each product, R rs was retrieved from a single pixel, matching-up with in situ sampling.

Field Data
Field data was collected from two boats by two teams within ±3 h from the Landsat 8/OLI image acquisition in six stratified lakes (surface areas ranging from 1.05 to 2.58 km 2 and mean depth ranging from 6 to 10.8 m) located about 100 km north of Berlin, Germany (red box on Figure 1A). The sampled lakes ( Figure 1B) have been previously classified as eutrophic lakes with index values between 2.8 and 3.3 (data kindly provided by the Ministry of Agriculture and the Environment Mecklenburg-Vorpommern, Germany) for the trophic index of the German working group on water issues (LAWA) [22]. The values of the LAWA Index (trophic index) for lakes are based on chlorophyll-a measurements, algal community composition and abundance of certain algal indicator species [22]. We have chosen a set of eutrophic lakes, since the L8PAR product requires explicitely further verification in eutrophic inland waterbodies [15]. During this field campaign chlorophyll-a concentrations ranged from 4 to 46 mg/m 3 , with higher values found in Lake Zotzen (ZOT, Figure 1C) and Lake Priepert (PRI, Figure 1F), and lower values in Lake Ellbogen (ELL, Figure 1G) and Lake Ziern (ZIE, Figure 1H).
In situ radiometric data were collected to calculate the proximal hyperspectral R rs directly above the water surface using two spectroradiometers: the Analytical Spectral Devices FieldSpec Pro (Analytik Ltd., Cambridge, UK) and the JETI Specbos 1211 UV (JETI (Jena Technische Instrumente GmbH), Jena, Germany). Both spectroradiometers provide data resampled to 1 nm wavelength step. We used two different spectroradiometers to increase the number of match-up points which were collected from two boats by two teams in different locations where there was no influence by bottom reflectance (mean depth of these lakes ranging from 6 to 10.8 m, while maximum depth ranges from 12 m in ZIE to 32 m in PRI) and away from the shore line (at least 6 pixels-180 m away from the shoreline). ASD FieldSpec Pro was used for data in ZOT (concomitant to image acquisition), VIL (Lake Vilz, +1 h of image acquisition) and LAB (Lake Labus, +3 h of image acquisition) while JETI was used for data acquisition in ZIE (−0.5 h of image acquisition), ELL (+1.5 h of image acquisition) and PRI (+3 h of image acquisition). Data acquisition and R rs calculation were based on the viewing angles and equations proposed by Mobley [23]. For each sampling point, ten measurements were averaged for the computation of the R rs spectrum from 380 to 780 nm. A diffuse reflectance plate with 20% reflectivity (Zenith lite, SG3145, SphereOptics, GmbH, Herrsching, Germany) was used to obtain the irradiance reference. The final R rs spectra were normalized by subtracting R rs~7 50 nm following the protocol described in Mobley [23]. A total of 32 hyperspectral R rs spectra were taken from these six lakes ( Figure 1C).

Product Evaluation
For the evaluation of atmosphericly corrected products, statistical metrics were computed to assess the accuracy of the retrieved R rs signals. Table 1 summarizes the list of error metrics used here. Table 1. List of error metrics.

Metric
Abbreviation Formula 1,2 1 X ret -retrieved value of R rs from each atmosphericly corrected product; 2 X obs -observed value of in situ R rs .

L8PAR Performance
The L8PAR product was used to retrieve R rs for the corresponding pixels where in situ samples were taken (see Section 2.3). The L8PAR product only provides R rs information for the first 4 bands of Landsat 8/OLI image which are in the visible range. For band 1 (centered at 443 nm), the L8PAR product was not able to retrieve any R rs value (Figure 2A), since the product calculated only negative values for aquatic reflectance which is then classified as invalid. For spectral band 2 (centered at 483 nm), the product could only retrieve 9 of the 32 data points and no relationship to in situ R rs was observed ( Figure 2B). From spectral band 3 (centered at 561 nm), the product could retrieve data for 23 sampling points and for spectral band 4 (centered at 655 nm), the product could retrieve R rs values from 16 locations. Both, spectral bands 3 and 4 showed a linear relationship to in situ R rs , however R rs values were almost always underestimated ( Figure 2C,D). The closest match-up between L8PAR and in situ R rs was observed in band 4 ( Figure 2D). Overall, the L8PAR product was not able to accurately retrieve information from these six lakes.
The L8PAR product was used to retrieve Rrs for the corresponding pixels where in situ samples were taken (see Section 2.3). The L8PAR product only provides Rrs information for the first 4 bands of Landsat 8/OLI image which are in the visible range. For band 1 (centered at 443 nm), the L8PAR product was not able to retrieve any Rrs value (Figure 2A), since the product calculated only negative values for aquatic reflectance which is then classified as invalid. For spectral band 2 (centered at 483 nm), the product could only retrieve 9 of the 32 data points and no relationship to in situ Rrs was observed ( Figure 2B). From spectral band 3 (centered at 561 nm), the product could retrieve data for 23 sampling points and for spectral band 4 (centered at 655 nm), the product could retrieve Rrs values from 16 locations. Both, spectral bands 3 and 4 showed a linear relationship to in situ Rrs, however Rrs values were almost always underestimated ( Figure 2C,D). The closest match-up between L8PAR and in situ Rrs was observed in band 4 ( Figure 2D). Overall, the L8PAR product was not able to accurately retrieve information from these six lakes.  Figure 3 shows the reflectance retrieved from different L8 products (L8PAR and L8SR) and different atmospheric correction routines (ACOLITE, C2RCC, C2X, iCOR) in comparison to the in situ hyperspectral Rrs spectra for the six different lakes and 32 measurement points. In this Figure, each lake was plotted with a different color ZOT (blue), VIL (red) LAB (green), ZIE (pink), ELL (yellow and gray) and PRI (dark cyan). For the L8PAR product divided by π, an underestimation of Rrs was noticed for all wavelengths ( Figure 3A). The L8SR divided by π overestimated the retrieval of Rrs for bands 1, 2 and 5 (centered at 865 nm). The high deviation around the near-infrared spectral band (centered at 865 nm) indicated the interference of adjacency effects as well as the effect of the reflected skylight which were not accounted in this standard land atmospheric correction routine ( Figure 3B). Therefore, considering that adjacency effects are most visible in the near-infrared [24], the L8SR product should not be used for aquatic studies, especially because it also does not consider the effect of the reflected skylight. ACOLITE got a good retrieval of Rrs in the visible range, however, for the near-infrared band it also overestimated the intensity of the Rrs ( Figure 3C). C2RCC and C2X  Figure 3 shows the reflectance retrieved from different L8 products (L8PAR and L8SR) and different atmospheric correction routines (ACOLITE, C2RCC, C2X, iCOR) in comparison to the in situ hyperspectral R rs spectra for the six different lakes and 32 measurement points. In this Figure, each lake was plotted with a different color ZOT (blue), VIL (red) LAB (green), ZIE (pink), ELL (yellow and gray) and PRI (dark cyan). For the L8PAR product divided by π, an underestimation of R rs was noticed for all wavelengths ( Figure 3A). The L8SR divided by π overestimated the retrieval of R rs for bands 1, 2 and 5 (centered at 865 nm). The high deviation around the near-infrared spectral band (centered at 865 nm) indicated the interference of adjacency effects as well as the effect of the reflected skylight which were not accounted in this standard land atmospheric correction routine ( Figure 3B). Therefore, considering that adjacency effects are most visible in the near-infrared [24], the L8SR product should not be used for aquatic studies, especially because it also does not consider the effect of the reflected skylight. ACOLITE got a good retrieval of R rs in the visible range, however, for the near-infrared band it also overestimated the intensity of the R rs ( Figure 3C). C2RCC and C2X had similar performances and both routines underestimated the spectral band 3 intensity ( Figure 3D and 3E). Nevertheless, they did not overestimate the spectral band 5. iCOR overestimated all spectral bands ( Figure 3F). Overall, ACOLITE was the atmospheric correction routine that got better results for the visible range, while C2RCC and C2X got the best result for the near-infrared and similar results for spectral bands 1, 2 and 4. Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 12

Further Assessment of L8PAR for Small Inland Waters
Retrieval of Rrs from L8PAR for OLI's band 1 and band 2 were poor for the image from 26 of July, 2019 ( Figure 2). The RGB image from this date ( Figure 4A) showed no cloud cover over the scene, and the Level 2 flag of the product ( Figure 4B) showed that pixels over inland waters were classified as 6161536 (in yellow) which includes the following flags: ATMWARN, MODGLINT, MAXAERITER, RRSWARN, COCCOLITH and SEADAS_CLOUD. These flags are set when the value is invalid because: (i) the sun glint coefficient is higher than the minimum; (ii) the aerosol corrections exceed the maximum threshold, (iii) Rrs is out of range, (iv) of presence of coccolithophore blooms, and/or (v) of cloud cover. However, some misflaggings were observed such as: (1) this scene does not contain clouds and the product is identifying clouds, (2) the Rrs is not out of range but flagged as out of range by the L8PAR product. Especially Rrs was flagged as out of range by L8PAR despite the fact that Rrs was able to be retrieved by ACOLITE, C2RCC and C2X (which also assess for these issues, including

Error Analysis for Each Spectral Band
To quantify the performances of each atmosphericly corrected product, statistical errors were evaluated per spectral band per single pixel in the sampling location. Band 1 (centered at 443 nm) was best estimated using ACOLITE, C2RCC and C2X with bias, MAD and RMSD of ±0.002 sr −1 , ACOLITE got the lowest NRMSD while C2RCC and C2X got the lowest MAPD (Table S1-in Supplement 1). C2RCC and C2X also scored the best performance for band 2 (centered at 483 nm) with bias and MAD of ±0.001 sr −1 , a RMSD of ±0.002 sr −1 , NRMSD of 31.08% and 32.02% and MAPD of 29.4% and 29.6% (Table S2 in Supplement 1). ACOLITE was the best to retrieve the R rs at 561 nm (band 3) with a bias and MAD of ±0.001 sr −1 , a RMSD of ±0.002 sr −1 , NRMSD of 8.33% and a MAPD of 7.87% (Table S3 in Supplement 1). Band 4, in the red spectrum, was best estimated when ACOLITE, C2RCC and C2X was used with bias, MAD and RMSD of ±0.002 sr −1 , C2RCC and C2X got a NRMSD of 17.99% and 17.82% while the lowest MAPD was from L8PAR with 25.07% (Table S4 in Supplement 1). C2RCC and C2X got the best performance for near-infrared band (band 5) with bias, MAD and RMSD of <0.001 sr −1 and NRMSD of 3.52% and 3.59% while they got a MAPD of 27.22% and 27.77% respectively (Table S5 in Supplement 1). Additionally, we used a 3 × 3 pixel window to compute the same estimators and noticed that there were few significant changes in atmospheric correction routine performances. For band 1 (centered at 443 nm) C2RCC improved the estimation from a NRMSD of 48.40 to 38.12% (Table S6-in Supplement 1). Also for band 2 (centered at 483 nm) some significant improvement of C2RCC with a NRMSD changing from 31.08% to 26.61% was shown (Table S7-in Supplement 1). For band 3 (centered at 561 nm) a deterioration in the performance of L8PAR was observed which changed from a NRMSD of 8.32% to 542.03% while a significant improvement was observed for L8SR with a NRMSD changing from 19.92% to 11.17% (Table S8-in Supplement 1). The same pattern was observed for band 4 (centered at 655 nm) which showed a deterioration of the performance of L8PAR (NRMSD changing from 29.94% to 552.43%) and a small improvement of L8SR (NRMSD changing from 33.85% to 28.92%). For band 5 in the near-infrared the performance was similar with a small deterioration of C2RCC which got a NRMSD of 5.02% (in the single pixel analysis it was 3.52%). Therefore, similar to the result for single pixels also in the case of using a 3 × 3 window, the use of ACOLITE, C2RCC and C2X were more consistently showing good performances for different bands. However, the use os a 3 × 3 window decreased the performance of the L8PAR product, which coud be related to the large number of invalid pixels in this product.

Further Assessment of L8PAR for Small Inland Waters
Retrieval of R rs from L8PAR for OLI's band 1 and band 2 were poor for the image from 26 of July, 2019 ( Figure 2). The RGB image from this date ( Figure 4A) showed no cloud cover over the scene, and the Level 2 flag of the product ( Figure 4B) showed that pixels over inland waters were classified as 6161536 (in yellow) which includes the following flags: ATMWARN, MODGLINT, MAXAERITER, RRSWARN, COCCOLITH and SEADAS_CLOUD. These flags are set when the value is invalid because: (i) the sun glint coefficient is higher than the minimum; (ii) the aerosol corrections exceed the maximum threshold, (iii) R rs is out of range, (iv) of presence of coccolithophore blooms, and/or (v) of cloud cover. However, some misflaggings were observed such as: (1) this scene does not contain clouds and the product is identifying clouds, (2) the R rs is not out of range but flagged as out of range by the L8PAR product. Especially R rs was flagged as out of range by L8PAR despite the fact that R rs was able to be retrieved by ACOLITE, C2RCC and C2X (which also assess for these issues, including the sun glint as recently demonstrated by Caballero et al. [25,26]) and the fact that the retrieved R rs was well within the range of 0-0.3142 specified in the L8PAR product guide. The L8PAR aquatic reflectance for band 1 ( Figure 4C) and for band 2 ( Figure 4D) showed that most of the inland waters do not have positive aquatic reflectance values (in white, except for few pixels of mostly large aquatic systems such as the Szczecin Lagoon-green box, and Lake Müritz-blue box on Figure 4). These negative values (in black) seem to be related to non-valid pixels due to the applied flags when using L8PAR. Band 3 ( Figure 4E) is the band which was able to retrieve more pixels (positive values in white) although for the small inland water the retrieval was limited. Band 4 had less positive pixels than band 3 ( Figure 4F).
showed some misflags, for example, the identification of sea ice in the image. The pixel quality assessment for all images were in the range of valid values . Therefore, although the pixel quality assessment for all seven L8PAR products over our study site showed valid pixels, flags and aquatic reflectance values were not accurate. Considering that most of the pixels over inland waters for all images (Table S6 in Supplement 2) got a ATMWARN flag, the L8PAR product should not be used for inland water studies in its current state.  In addition to the scene used in the presented study, we qualitatively analyzed six additional L8PAR scenes with the same path and row from different periods of the year (Table S6 in Supplement 2) to extend this assessment. L8PAR products from October 14 and 30, 2019, and May 9, 2020, showed that the Level 2 flag classified all inland water pixels in the scene as 4,915,776.75 (Figures S1, S2 and S6 in Supplement 2). This class involves the following flags: ATMWARN, MAXAERITER, NAVWARN, CLOUD and COASTZ. These flags are set when the value is invalid because: (i) the aerosol corrections exceed the maximum threshold, (ii) of cloud cover and iii) of bathymetry (if it has less than 30 m depth). Scenes from March and April, 2020 ( Figures S3, S4 and S5 in Supplement 2) showed that the level 2 flag classified inland water pixels as 23,857,216. This class involves the following flags: SEAICE, ATMWARN, CHLWARN, MAXAERITER, RRSWARN, TURBIDW and COASTZ. These flags are set when the value is invalid because: (i) of ice cover, (ii) the aerosol corrections exceed the maximum threshold, (iii) chlorophyll-a values are out of the range, (iv) R rs is out of the range, (v) water is turbid, and/or (vi) of bathymetry [18]. These additional scenes also showed some misflags, for example, the identification of sea ice in the image. The pixel quality assessment for all images were in the range of valid values . Therefore, although the pixel quality assessment for all seven L8PAR products over our study site showed valid pixels, flags and aquatic reflectance values were not accurate. Considering that most of the pixels over inland waters for all images (Table S6 in Supplement 2) got a ATMWARN flag, the L8PAR product should not be used for inland water studies in its current state.

Atmospheric Correction Routines Comparison
The study of the performance for R rs retrieval of two different Landsat 8 products (L8PAR and L8SR) and four atmospheric correction routines (ACOLITE; C2RCC, C2X and iCOR) showed that ACOLITE, C2RCC and C2X scored the best performance for the Landsat 8/OLI bands. L8PAR failed to retrieve any data for band 1 (443 nm) and was not able to produce data in the other three bands for several pixels. Overall, L8PAR strongly underestimated R rs for almost all calculated data points and a poor relationship with measured in situ R rs data was observed for band 2 and 3, while performance for band 4 was best and close to ACOLITE, C2RCC and C2X.
Although this study only used one single cloud-free scene, the performances of the atmospheric correction routines were consistent to previous studies of atmospheric correction for inland waters at different atmospheric conditions. Pereira-Sandoval et al. [27] evaluated different atmospheric correction routines for Spanish reservoirs from 2015 to 2018 using 21 cloud-free images from Sentinel 2/Multispectral Instrument (MSI). Similarly to our results, they also observed a poor performance of iCOR for different water types. In their study, C2RCC and C2X also showed an overall good performance, suggesting that the look-up-tables methods of simulated water leaving radiance is providing good results over semi-analytical methods. In contrast to our study, C2RCC and C2X showed slightly different performances in different bands, which we attribute to the more turbid and more productive waters investigated in their study. The main difference to our study was the performance of ACOLITE which got a poor overall performance for the Spanish reservoirs, suggesting that ACOLITE is not performing well for extremely turbid waters. Warren et al. [28] evaluated the application of different atmospheric correction routines for coastal and inland waters using Sentinel 2/MSI imagery. For inland waters, they observed that C2RCC got the lowest errors for most of Sentinel 2/MSI spectral bands and that ACOLITE got the lowest error for 705 nm band. Additionally, they also observed that iCOR had a poor performance for all spectral bands. These findings from 163 sampling points from different European inland waters also agree with our results shown in Figure 3 and Tables S3-S7.
For the Landsat 8 products (L8PAR and L8SR), which were developed for open ocean, coastal waters and land applications, results were expected to be less accurate than the products specifically tailored for inland waters. The major source of errors contributing to the uncertainties for the atmospheric correction routines for open ocean is typically the use of two or more near-infrared bands where the signal is assumed to be zero [29]. This assumption is valid for clear waters where pure water absorption is very high in this spectral region and there is no contribution from suspended particles. However, for turbid inland waters, where the concentration of suspended particles is high, the signal in the near-infrared region is not negligible [10]. Since turbid waters were not considered in the look-up-tables for the development of the product, the spectral information content from these waters may not be captured in the modeled spectra [19].
The strong deviation between the L8PAR product and in situ radiometric data at the sampling points obtained in our study suggests that further improvements are needed to make L8PAR applicable for inland turbid waters. Since the algorithm used to generate the L8PAR product was based on datasets from coastal waters which did not include extremely eutrophic, turbid or CDOM-rich waters [18,19], it was expected that the product would not perform well for turbid inland waters (especially when turbid waters is a flag of the product as well as depths lower than 30 m). However, the fact that the L8SR product was performing better than or equally to the L8PAR was a surprise since L8SR does not consider the effects of reflected skylight (Figure 2 and Tables S4-S6).
Overall, L8PAR did not perform well for retrieval of/failed to retrieve water reflectance for the lakes investigated in this study and was not able to provide valuable data for the lakes investigated in this study. Specifically, L8PAR could not successfully retrieve information for all sampling points in any of the four spectral bands that are available in the product. The reason for the underperformance of the product is probably caused by sensitivity to aerosol type estimation and aerosol extrapolation to these spectral bands. Additionally, commonly present in inland waters were not considered in the L8PAR product. Further studies are necessary and should cover larger data sets for different trophic conditions and different sizes of inland waters, as well as different levels and types of adjacency effects (forest × agriculture) and different weather (or atmospheric) conditions. While our results are not positive for the use of L8PAR for continental waters, a recent indirect evaluation of the L8PAR was conducted by Nazeer et al. [30] for coastal waters in the region of Hong Kong, China. In that study authors used L8PAR to calculate chlorophyll-a concentrations using an empirically calibrated bio-optical algorithm. Differently than the results presented in our study, L8PAR was able to retrieve data from the bands used in their algorithm for the coastal waters of Hong Kong, China. The better performance observed for coastal waters of Hong Kong could be related to the fact that L8PAR was calibrated using coastal waters dataset [18,19] and the use of flags that are appropriated for coastal and oceanic waters such as the bathmetry flag which flag pixels in areas where the depth is lower than 30 m. Nevertheless, in the study of Hong Kong coastal waters, L8PAR was not directly evaluated since no remote sensing reflectance measument was used, only chlorophyll-a concentration.
The development of analysis ready data (ARD) such as the L8PAR product are highly needed since its spatial resolution of 30 m makes it potentially applicable to the monitoring of small-sized inland waters which are of great abundance worldwide [31]. Therefore, the improvement of the L8PAR product will be a great contribution for studying inland water quality and to have a ready-to-use product for the development of bio-optical algorithms. Additionally, in combination with Sentinel 2 MSI, satellite imagery could be used to provide spatial and temporal information to support water quality monitoring programs [16]. To guide further development, this evaluation of L8PAR for inland waters aims to highlight that: (1) the L8PAR product in its current state should be useed with caution, especially for small and turbid inland waters; and (2) further evaluation of the product for inland waters under different atmospheric conditions is needed.

Conclusions
This study provides an evaluation of the novel L8PAR product which derives aquatic reflectance. The L8PAR product showed a poor relationship to the in situ R rs from six turbid/eutrophic German lakes for a single scene with best performance. The main problem of L8PAR was the estimation of negative values for aquatic reflectance in inland waters. Compared to other atmosphericly corrected products which were able to estimate R rs values for all 32 sampling locations, the L8PAR in its current state (version 1.0) is not suitable for realtively small, turbid inland waters. Based on the comparison with other atmospheric correction routines, the ACOLITE routine performed well for spectral bands in the visible range, while C2RCC and C2X performed well when considering the near infrared band. Therefore, for a retrieval of spectral bands in the visible and near-infrared, we woud for now recommend the use of C2X and C2RCC.
This study validates the L8PAR product in inland waters. It is important to emphasize that the results of the algorithms comparison presented here originate from a single scene. However, our study further showed that for multiple other scenes, the L8PAR product could not retrive aquatic reflectance from our study sites (Figure 4 and Figures S1-S6) and from large waterbodies also presented in the same image (Szczecin Lagoon and Lake Müritz). This happened due to misflagging of the L8PAR product which are using flags that are not appropriated for inland waters. We would like to highlight that the use of L8PAR should be carefully evaluated before its use, especially for eutrophic, turbid inland waters. Future studies need to be conducted for different types and sizes of inland waterbodies, as well as different atmospheric conditions. Although Landsat 8 is limited by the low revisit time (16 days), the harmonization of OLI imagery with Sentinel 2/MSI bears great potential for inland water remote sensing, where an improved L8PAR product that works reliably also for small and turbid lakes, like the ones used in this study, is highly desired.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/15/2410/s1, Table S1: Statistical error metrics for the estimation of Rrs for band 1 (443 nm). Shaded areas indicate best performances, Table S2: Statistical error metrics for the estimation of Rrs for band 2 (483 nm). Shaded areas indicate best performances, Table S3: Statistical error metrics for the estimation of Rrs for band 3 (561 nm). Shaded areas indicate best performances, Table S4: Statistical error metrics for the estimation of R rs for band 4 (655 nm). Shaded areas indicate best performances, Table S5: Statistical error metrics for the estimation of R rs for band 5 (865 nm). Shaded areas indicate best performances, Table S6: Statistical error metrics for the estimation of R rs for band 1 (443 nm) using a 3 × 3 window. Shaded areas indicate best performances, Table S7: Statistical error metrics for the estimation of R rs for band 2 (483 nm) using a 3 × 3 window. Shaded areas indicate best performances, Table S8: Statistical error metrics for the estimation of R rs for band 3 (561 nm) using a 3 × 3 window. Shaded areas indicate best performances, Table S9: Statistical error metrics for the estimation of R rs for band 4 (655 nm) using a 3 × 3 window. Shaded areas indicate best performances, Table S10: Statistical error metrics for the estimation of R rs for band 5 (865 nm) using a 3 × 3 window. Shaded areas indicate best performances. Figure