Atmospheric correction algorithm over coastal and inland waters based on the red and NIR bands: application to Landsat-8/OLI and VNREDSat-1/NAOMI observations.

Water pixel extraction and correction of the atmospheric signal represent prerequisite steps prior to applying algorithms dedicated to the assessment of water quality of natural surface water bodies. The recent multiplication of medium spatial resolution sensors (10-60 m) provides the required constellation to monitoring bio-optical and biogeochemical parameters of surface waters at the relevant spatial-temporal scales. Here we present a new approach to identify water pixels and to extract the atmospheric contribution to the top of atmosphere signal measured by the NAOMI sensor on board the first Vietnamese satellite, VNREDSat-1. After verifying the TOA calibration of NAOMI through a vicarious calibration exercise, we adapt a recent water pixel extraction algorithm (WiPE) to NAOMI, and develop a new atmospheric correction algorithm (referred to as red-NIR) based on the use of the red and NIR bands (the only bands available for that purpose on NAOMI) and spectral relationships. The evaluation of red-NIR with a match-up data set gathering remote sensing reflectance, Rrs, measurements performed at the AERONET-OC stations in moderately turbid waters indicates excellent performance in the blue and green part of the spectrum (similar to the performances reached by the SeaDAS NIR-SWIR algorithms) and lower accuracy in the red. Intercomparison of simultaneous images collected by NAOMI and OLI over a more turbid water body shows an excellent agreement between the NAOMI-Rrs estimated by the present processing, and the OLI-Rrs estimated from the ACOLITE algorithm. This approach will allow sensors that do not have SWIR bands, such as SPOT-6 and -7, to be processed, making their data exploitation available for long-term temporal analyses.


Introduction
The tight societal and economical connections of Vietnam with the sea, mainly due to its 3260 km long coastline, makes this country very sensitive to natural and anthropogenic disasters. Vietnam is one of the countries that will be severely impacted by climate change through the modification of precipitation patterns and sea level rise [1,2]. The abrupt changes occurring in the Mekong Delta area [3][4][5], the recurrent flooding or pollution events observed in different provinces [6,7], the potential link between water quality and aquaculture activities [8], are, among others, phenomena stressing the urgent need of a better regional monitoring of coastal and connected inland waters. In that context, satellite remote sensing represents a valuable technique for a better assessment of the status and evolution of Vietnamese coastal and inland waters. While standard ocean color observations provide relevant information about the bio-optical status of surface coastal waters at a daily coverage (in absence of clouds), their spatial resolution (between 300 m and 1 km) may not fully be adapted for the survey of specific areas such as sea harbors, aquaculture areas, industrial and urban coastal areas, and river outlets.
The recent radiometric improvement of medium spatial resolution (from 10 to 60 m) sensors now allows the monitoring of water surface quality at the required spatial resolution [9][10][11][12][13][14][15][16]. However, the relatively low revisit time (between 10 to 16 days under clear days) of these sensors is not adapted to monitor the complex physical and biological interactions occurring in coastal waters. Combining these sensors into a constellation allows the temporal resolution to be greatly increased, and thus the required time-scales to be better covered. For instance, the Sentinel-2 constellation provides a revisit time of 5 days at the equator in cloud-free conditions (instead of 10 days when the two Sentinel-2 are considered separately), and 2.9 days if combined with Landsat-8 [17]. To reduce the impact of cloud coverage on the surface water remote sensing monitoring, it is therefore essential to increase the number of observations. For that purpose, the first multi-spectral optical satellite in Vietnam, VNREDSat-1 has been launched on May 7, 2013. VNREDSat-1 [18] carries the NAOMI (New AstroSat Optical Modular Instrument) sensor [19], which is similar to SPOT-6 and -7, and which includes four multi-spectral bands and one panchromatic band at a spatial resolution of 10 and 2.5 m, respectively (Table 1). In order to use NAOMI images for aquatic applications over coastal and inland waters, water pixels have to be identified, while the top of atmosphere signal recorded by the sensor above these water surfaces has to be corrected from atmospheric effects. This is the objective of the present study. The extraction of water pixels on NAOMI images is performed in the present study by adapting the recent water pixels extraction algorithm (WiPE) developed for Landsat-8 Operational Land Imager (OLI) and Sentinel-2 Multispectral Instrument (MSI) [20]. Numerous approaches have already been developed to perform atmospheric correction over water pixels delivered by medium spatial resolution sensors observations [9,[21][22][23][24][25][26][27]. Depending on the number of NIR bands, and the availability of a SWIR band, such as on HSV (SPOT) or ETM+ (Landsat), the atmospheric correction procedures over water surfaces are based on different approaches. First, some directly apply a radiative transfer code such as 6S, combined with external information on the composition of the atmosphere, such as the aerosol optical thickness, the atmospheric profile, and the aerosol model. Others use the combination of radiative transfer calculations (6S) and the black target approach [28], as used by [29]. For the OLI sensor on Landsat-8 or MSI on Sentinel-2 (A and B) the different existing approaches essentially differ according to the use or not of several near-infrared (NIR) or short-wave infrared (SWIR) bands, and by their inherent assumptions. Different validation exercises have recently been performed to evaluate their respective performances over coastal and inland waters [25,[30][31][32][33]. However, despites these numerous inter-comparison exercises, a consensus on the best atmospheric correction method over coastal waters has not been reached yet. The purpose of the present paper is to develop a self-consistent (i.e. without any external input on aerosols) atmospheric correction algorithm (referred to as Red-NIR) for NAOMI, which has only one NIR band and no SWIR information. To reduce uncertainties related to the use of various in situ data sets, the performance of the present algorithm is evaluated using the ocean color component of the Aerosol Robotic Network (AERONET-OC) [34]. Because of the very limited number (N = 6) of coincident in situ measurements and VNREDSat-1 overpass, the Red-NIR algorithm is applied to OLI which has four spectral bands in common with NAOMI (Table 1), and for which the number of available match-up data points (N = 67) is more suitable for a validation exercise. A recent study emphasized that, among two current atmospheric correction algorithms developed for OLI (SeaDASS and ACOLITE), and tested over the AERONET-OC data set (i.e. over moderately turbid waters), SeaDAS performs better [33]. For that purpose, the performance of Red-NIR will be compared to that obtained from atmospheric correction algorithms implemented in SeaDAS and using different band combinations in the NIR and SWIR. The objective is here to show how the Red-NIR algorithm behaves in comparison with similar approaches using SWIR bands, and not to compare to all available atmospheric correction algorithms available for OLI. While the AERONET-OC match-up data set only allows evaluating Red-NIR over moderately turbid waters, the unique scene sampled almost simultaneously by VNREDSat-1 and Landsat-8 over turbid area will be used to assess its performance for higher turbidity levels. For that purpose, the ACOLITE [27] algorithm, which shows relatively good performances over highly turbid waters, will be used.
The radiometric in situ data set used in the algorithm development, and the match-up data set used for its validation are first presented. Due to the relatively low contribution of the water leaving radiance to the top of atmosphere (TOA) signal, it is essential to verify the absolute calibration of the sensor for ocean color applications, which has still not been performed for NAOMI over water surfaces. For that purpose, in situ remote sensing reflectance measurements performed at different stations of the AERONET network were propagated through the atmosphere using a radiative transfer code [35], and compared to the TOA reflectance measured by NAOMI. Then, the adaptation and validation of the WiPE algorithm to NAOMI is presented, and the new atmospheric correction method based on the use of the red and NIR bands available on NAOMI is described. Validation as described above is then provided.

Match-up data set
The match-up data set gathers 67 Landsat-8 scenes collected over 10 Aeronet-OC [34,[36][37] stations, and 6 VNREDSat-1 scenes over Venice and Lucinda Aeronet-OC stations (Fig. 1). The Landsat-8 scenes were downloaded from the Glovis portal (https://glovis.usgs.gov) of U.S. Geological Survey (USGS), while the VNREDSat-1 scenes were received from the center for control and exploitation of small satellite at the Space Technology Institute (STI) of the Vietnam Academy of Science and Technology (VAST). In situ AERONET-OC [37] radiometric data were collected at three spectral bands centered at 443, 551, and 667 nm. The match-up data selection is based on the same criteria as those applied in [25] for coastal waters and high spatial resolution observations. After application of these match-up criteria 67 and 6 match-up data points remained for Landsat-8 ( Table 2) and VNREDSat-1 (Table 3), respectively. This data set gathers relatively clear to eutrophic waters, with 26% of the in situ spectra presenting a maximum at 443 nm, and 74% at 531 nm. This data set is then characteristic of clear to moderately turbid waters. The wind speed values recorded at the AERONET stations are generally lower than 5 m s −1 .

Radiometric in situ data set used for trans-spectral relationships
The atmospheric algorithm developed for VNREDSat-1 in the present study requires the use of spectral relationships between R rs (865) and R rs (655), and R rs (655) and R rs (561) (see section 6.2). For that purpose, a large in situ data-set of remote sensing reflectance has been gathered from radiomeric measurements performed in the southern North Sea and English Channel [38,39], the Celtic Sea, the Ligurian Sea, the Adriatic Sea and in the Atlantic Ocean along the coasts of Portugal and French Guiana [40,41]. All measurements were performed using the TriOS-RAMSES hyper-spectral sensors. The measurement protocols and associated processing are fully described in [39] and [42]. This data set covers various optically contrasted coastal waters. While the clearest waters, characterized by a maximum of R rs (λ) in the blue, represent about 39% of the data set, green (with a maximum of R rs around 561 nm), turbid (maximum of R rs in the red) and very turbid waters (with a maximum of R rs at 700 nm) represent about 55.8, 2.4, and 2.8% of the total reflectance spectra, respectively. The suspended particulate matter (estimated based on [38] at 865 nm), SPM, concentration varies between 1.5 and 726.2 mg.l −1 , with a mean (median) and standard deviation value of 36.8 (8.2) ± 87.0 mg.l −1 .

Statistical indicators of model performance
Graphical comparisons of observations with model predictions are complemented with quantitative statistical metrics in order to evaluate model performance over the match-up data set. Statistical indicators that are typically utilized in the assessment of atmospheric correction model accuracy [25] are also calculated. These indicators include the root-mean-square difference (Eq. (1)), RMSD, the median percentage difference (Eq. (2)), MPD, and the mean bias (Eq. (3)), MB: where k is the match-up index, and N is the numbers of selected match-up.

Verification of the TOA calibration of NAOMI/V1
The top of atmosphere (TOA) radiance of NAOMI, L TOA (λ), is calculated through Eq. (4).
where DN(λ), bias(λ), gain(λ) are the digital number, the bias, and the gain spectral values provided for each NAOMI image. In a similar way, the top of atmosphere (TOA) radiance of OLI, L TOA (λ) (in W m −2 sr −1 µm −1 ), which will be used to test the present atmospheric correction algorithm due to the large match-up data set for OLI, is defined as follows: where M L (λ) is the radiance multiplicative scaling factor (i.e. gain) for each band, Q cal (λ) the Level-1 pixel value in digital number, and A L (λ) the radiance additive scaling factor (i.e. offset) for each band.
The TOA reflectance, ρ TOA (λ) for OLI and NAOMI are then calculated from L TOA (λ) using Eq. (6): where d, F 0 , and θ 0 are the Earth-Sun distance in Astronomical Units, the band average solar irradiance [43], and the sun zenith angle, respectively. A comparison between measured and simulated NAOMI TOA reflectance is performed to check the radiometric calibration of the sensor. Simulations are carried out using the Ocean Successive Orders with Atmosphere -Advanced (OSOAA V1.4) radiative transfer code [35]. The OSOAA code simulates the propagation of the radiation through the Ocean-Atmosphere including coupling terms, and the interactions with the wind-ruffled sea surface. The code outputs (I out ) are radiances normalized to an extraterrestrial solar irradiance set to π: Equation (6) and (7) lead to Eq. (8), that converts I out into TOA reflectance to compare directly with NAOMI measurements: Simulations are carried out for the exact sun-sensor geometry of the considered scenes presented in Table 3 and for λ from 420 to 920 nm, with a 1 nm step. For each channel, the spectral TOA reflectance is convoluted with the spectral response function of NAOMI (Fig. 2) and is finally integrated to generate the simulated bandpass average. The comparison analysis is realized over the Venice (Italy) and Lucinda (Australia) Aeronet-OC sites. The AERONET-OC ground-based measurements (the atmospheric pressure P 0 and the wind speed W at sea level, the aerosol optical thickness AOT(λ), and the normalized water-leaving radiance L WN (λ)) are used as code inputs. Coastal aerosols are defined according to [44]. Note that the normalized water-leaving radiance and aerosol optical thickness (from 0.02 and 0.25 in the NIR depending on the image) from AERONET-OC are interpolated over the wavelength to provide inputs for the code with a 1 nm resolution. The wavy sea surface is modeled in OSOAA using the Cox-Munk wave slope probability density depending on the surface wind speed [45]. Reflection by individual whitecaps is neglected, which is adequate for the purpose of our study as the wind speed is low (W < 2 m.s −1 ). The pressure at sea level (P 0 ) is used to derive the molecular optical thickness needed to simulate the Rayleigh scattering. Concerning the water body, the in-water part of the code has not been considered as the normalized water-leaving radiance, measured at the AERONET-OC site and converted into reflectance (Eq. (5)), has directly been injected into OSOAA (considering a lambertian surface): Figure 3 displays the NAOMI measured and simulated TOA reflectance. Considering the 6 match-up data points, the relative difference values calculated between the measured and simulated TOA reflectances, are 10% in the blue, 3.3% in the green, and 15% in the red. This number are similar to those reported for a similar exercise performed for the GOCI sensor for which the difference is between 1.05 and 8.14% [46]. The relative high value in the red is explained by the low absolute value of the signal, and to the signal over noise ratio of NAOMI.
No specific relationship has been found between the aerosol optical thickness and these TOA relative differences. Supplementary match-up data points are however required to confirm these results and to derive statistically robust vicarious gains.  (Table 3) and using the measured normalized water-leaving radiance.

The water pixel extraction procedure for VNREDSat-1
Detection and extraction of surface water pixels is a prerequisite step prior to applying atmospheric correction algorithms dedicated to the estimation of R rs (λ) from L TOA (λ). For that purpose, the water pixel extraction algorithm (WiPE) recently developed for the Landsat-8 OLI and Sentinel-2 MSI sensors [20] has been adapted to VNREDSat-1 for the present study. This algorithm which requires the Rayleigh-corrected reflectance spectrum, ρ rc (λ), as input paramet, is based on the application of two successive steps. First, different spectral criteria are applied to mainly disregard clouds, vegetation, barren land, and construction pixels. Due to the complex radiative interactions occurring between thin clouds or cloud shadow and the overlaid water pixels, the identification of these two latter objects cannot properly be performed using the spectral shape analysis developed in the first step. For that purpose, the ρ rc (λ) spectra are transferred into the Hue, Saturation Value (HSV) color space to improve the distinction between water pixels and thin cloud and shadow pixels over water areas. As the complete description of the WiPE algorithm and its performances are provided in [20], we only provide here the main modifications for the adaptation of WiPE to VNREDSat-1.
In the present study, 17 VNREDSat-1 images, collected over various bio-optical environments, similar to the ones described in [20], have been gathered for the development data set. This spectral database accounts for seven different objects: water (43,843 pixels), cloud (779 pixels), thin cloud (783 pixels), vegetation (3,122 pixels), construction (360 pixels), barren land (390 pixels), and shadow (1,503 pixels). Due to the lower number of spectral bands of NAOMI compared to OLI, the two different steps of the WiPE algorithm for NAOMI are based on a restricted number of spectral criteria (Fig. 4). The spectral analysis based on the development data set shows that if ρ rc (1) is greater than -0.12 ρ rc (4)/ρ rc (3) + 0.228 and if ρ rc (4)/ρ rc (3) is greater than 1.14 then cloud, some thin cloud and land pixels are removed (Fig. 5). In addition to removing the totality of the vegetation, cloud, barren land, and construction pixels, the application of this criteria also removed many thin clouds (91% of total thin clouds pixels) and some shadow pixels (7% of total shadow pixels). The remaining clouds shadow and thin cloud pixels may be partly removed by the HSV analysis where different criteria are applied depending on the wavelengths at which V max is reached (Fig. 4). At the end of the whole processing 2% of thin clouds and 82% of cloud shadows remain on the development data-set. The large percentage of remaining cloud shadow pixels is due to the unavailability of SWIR bands, in contrast to the WiPE version developed for OLI [20]. The WiPE algorithm adapted to NAOMI is validated using 5 images not included in the development data set and collected over contrasted areas (Table 4). For each image, a reference water pixel map (referred to as the reference WPM hereafter) was derived as in [20] using the QGIS software. The performance of the WiPE algorithm developed for NAOMI is then quantified using the Mean Absolute Percentage Difference (MAPD in %) calculated between the reference number of water pixels (WPM) and the number of water pixels generated by WiPE. The MAPD varies between 0.07% to 5.6% over the five validation scenes, the highest MAPD value (5.6%) being found for the scene collected over the Tonle Sap area (Fig. 6(c), 6(d)) which presents the highest cloud coverage (25%) and cloud shadow. The algorithm allows us to remarkably  (3) band ratio for the seven different objects of the development data set (each color corresponds to a different object as indicated). The black lines correspond to the limit adopted to remove vegetation, clouds, construction, barren land, and thin clouds (step 1.1). The number of pixels before (panel a) and after (panel b) the application of the criterion adopted in steps 1.1 are provided in each panel.
extract water pixels over very turbid waters, without confusion with land or clouds pixels, and also in complex environments such as aquaculture areas (Fig. 7(c), 7(d)), but seems however to underestimate the number of cloud shadow pixels. This latter pattern, which mainly explains why WiPE slightly over estimates the number of water pixels (Table 4), is due to the restricted number of pertinent bands to perform a proper correction based on the HSV analysis [20].

Basics of atmospheric correction
The TOA reflectance can be written as shown in Eq. (10) [47][48][49]: with ρ R (λ), the reflectance due to multiple scattering of a purely air molecules atmosphere; ρ a (λ), the reflectance of multiple scattering by aerosols in a pure aerosol atmosphere; and ρ Ra (λ), the interaction term in a real atmosphere containing both molecules and aerosols. T(λ) and t(λ) are the direct and diffuse transmittance of the atmosphere, respectively. ρ g (λ), ρ wc (λ), and ρ w (λ) are the reflectances due to sun glint, whitecaps, and in-water optically significant components, respectively. In the present study, the reflectance of whitecaps is estimated using the wind speed and the geometry [47,50,51]. The reflectance due to Rayleigh scattering in a purely molecular atmosphere is removed using the air pressure and the geometry [47]. Only glint-free images, identified from the geometry of observation and wind speed, are considered in this study. Correction of glitter images can however be performed following the different available methods using NIR and summarized in [52]. Finally, the atmospheric correction algorithm developed for NAOMI is based on the Rayleigh-corrected reflectance ρ rc (λ) defined as follows: The goal of the atmospheric correction is to determine ρ a and ρ Ra so the water-leaving reflectance, ρ w , can be estimated as: t(λ) is computed knowing the diffuse transmittance for the sun-sea (t 0 ) and sea-sensor (t v ) directions. τ r and τ oz are the band average Rayleigh and Ozone optical thickness and are estimated for a standard atmosphere and using the air pressure. θ 0 , θ v are sun and sensor zenith angles, respectively. At first approximation, the aerosol transmittance is set to 1.0 in the present study. Over open ocean waters, the black pixel assumption in the NIR (i.e., the ocean being totally absorbent) is used to estimate the aerosol models and optical properties that are then used to estimate ρ a and ρ Ra [48]. Over optically-complex waters, this assumption is not valid anymore and it is necessary to deal with a non-zero ρ w [21,[53][54][55][56][57][58]. This is the purpose of the following section, where an atmospheric correction algorithm over coastal waters is developed taking into account the limited number of NAOMI spectral bands.

Description of the atmospheric correction scheme
The atmospheric correction scheme presented here is decomposed into several steps (Fig. 8). The first step is to use the WiPE algorithm to identify the water pixels (section 5). Then the clearest water pixel is extracted over the region of interest. If blue pixels, defined as pixels for which ρ rc (482)>ρ rc (561)>ρ rc (655), exist, then the pixel with the combined maximum value of the ρ rc (482)/ρ rc (655) Rayleigh corrected reflectance ratio, and minimum ρ rc (865) value is selected. If no pixel fulfills the first conditions, the one with the minimum ρ rc (865) value is considered as the clearest pixel. Once the clearest pixel is extracted, we consider that ρ a (865) is equal to ρ rc (865) for this pixel, as ρ Ra (865) can be ignored in first approximation. Knowing ρ a (865) at the clearest pixel, led to the determination of ρ a (λ L ) at shorter wavelengths using the epsilon parameter [48]: ρ a (λ L ) = (ε 655,865 ) n ρ a (865) (14) with n = 865−λ L 865−655 and ε 655,865 = ρ a (655) ρ a (865) (15) Combining Eqs. (14) and (15), a first guess of the water leaving reflectance in the green band (561 nm) can then be estimated as follows: As mentioned previously, we considered at the first step that ρ a (865) and ρ a (655) were equal to ρ rc (865) and ρ rc (665) at the clearest pixel, respectively. At the second step, these values are re-calculated as we considered that the water leaving signal was negligible in the first step of the algorithm. The new values of ρ a (655) and ρ a (865) at the clearest pixel can be estimated as follows: ρ a (655) = ρ rc (655) − t(655).ρ w (655) (17) ρ a (865) = ρ rc (865) − t(865).ρ w (865) (18) At this point, we only have a first guess of ρ w (561), that will be used to estimate ρ w (665) using an empirical relationship developed from the in situ data set presented in section 2.2 ( Fig. 9(a)): Similarly, ρ w (865) is estimated from ρ w (655) using the following empirical relationship ( Fig. 9(b)): The new values of ρ a (655) and ρ a (865), considered equal to ρ rc (655) and ρ rc (865), respectively, at the first step of the algorithm, can therefore be calculated knowing ρ w (655) and ρ w (865) as follows: ρ a (865) = ρ rc (865) − t(865)(25.1ρ w 3 (655) − 1.09ρ w 2 (655) + 0.107ρ w (655) − 0.0000237) (22) These two new values of ρ a in the red and NIR are then used to update the atmospheric correction parameters ε i,j for a given scene, implying that the aerosols optical thickness and types are even over the whole scene, which is similar to the approach developed by [54]. The ρ a reflectances at the shorter bands λ L are then calculated by using Eq. (14), as explained previously (Eq. (16) has been written to compute ρ w (561), however the same equation is used for ρ w (482) where 561 should be replaced by 482). Finally, the water-leaving reflectances at the NAOMI spectral bands are estimated using Eq. (16), and the whole process is repeated one time when ρ a (655) and ρ a (865) are re-calculated from Eq. (17) and Eq. (18) using the new values of ρ w (655) and ρ w (865), respectively. We verified, using the match-up data set, that this iterative procedure allows us to improve the R rs (λ) retrieval.

Match-up exercise for OLI based on the NAOMI configuration
The performance of the Red-NIR algorithm is evaluated from the 67 match-up data points ( Fig. 9 and Table 5). On average, considering the whole spectrum ( Fig. 9(a)), the Red-NIR algorithm is able to assess R rs (λ) with a RMSD of 0.0017 sr −1 , a MPD of 13.3%, and a MB of 5.3 10 −4 sr −1 . However, this average good pattern (r 2 =0.75) masks some important spectral differences. Among the three available visible bands, the green band presents the best R rs retrieval accuracy (RMSD = 0.0015 sr −1 ; MPD=-3.21%; MB=-3.7 10 −4 sr −1 ), whereas R rs in the red is estimated with the lowest accuracy (RMSD = 0.0012 sr −1 ; MPD = 54%; MB = 5.5 10 −4 sr −1 ). This latter pattern is likely linked to uncertainties of Eq. (20) for low level of turbidity, and to the low level of signal at this band (lower than 0.003 sr −1 ) making the inversion more challenging. As a matter of fact, almost all data points have a R rs (655) value lower than 0.003 sr −1 , which corresponds to suspended particulate matter concentration lower than about 3.5 g. m −3 [59].
The performance of the Red-NIR algorithm is now compared to the ones achieved using the three different NIR and SWIR band combinations available in SeaDAS, which performances have been already tested for OLI in [25] (Table 5, Fig. 11). While the Red-NIR algorithm ( Fig. 11(a)) provides better results when all bands are pulled together (Table 5), it is also clearly visible that the use of NIR-SWIR bands (Figs. 11(b), 11(c)) allows R rs to be estimated with a much better accuracy in the red part of the spectrum. For instance, the RMSD value is 54% by the Red-NIR model, and -19% and -22% for the NIR-SWIR1 and NIR-SWIR2 SeaDAS combinations, respectively. For the blue and green bands, the Red-NIR algorithm and the NIR-SWIR1 combination provides about the same results, while the NIR-SWIR2 combination is slightly lower in terms of performance. Note however, that the scatter around the 1:1 line is always lower for the NIR-SWIR combinations (Figs. 11(b), 11(c)) than with the Red-NIR model ( Fig. 11(a)). At last, the combination of two SWIR bands (Fig. 11(d)) always provides lower performances compared to the NIR-SWIR combinations, in agreement with the results of [25], but also when compared to the Red-NIR model.
The match-up data set allowed the performance of the Red-NIR model to be evaluated and compared with the SeaDAS algorithms only over slightly turbid waters. The performance of the WiPE + Red-NIR process for turbid to very turbid waters is illustrated by a scene sampled by OLI over the Hangzhou Bay in China on February 8, 2015 (Fig. 12). This scene clearly stresses the interest of WiPE (Fig. 12(b)) for the identification of water pixels over very turbid waters characterized by a maximum of R rs (λ) in the NIR (not shown). This is consistent with recent  results showing that WiPE is able to assess water pixels over very turbid water surfaces [20], in contrast to other approaches such as the approach used in SeaDAS (Fig. 12(c)). A restricted inland area, characterized by a relatively high turbidity level with maximum R rs (λ) values of about 0.02 sr −1 in the green or red channels (depending on the considered pixel), allows the Red-NIR R rs (λ) and NIR1-SWIR R rs (λ) estimations to be compared (Figs. 12(d), 12(e), 12(f)). For this level of turbidity, an excellent agreement between the two processings is found in the red ( Fig. 12(g)) and green (Fig. 12(h)) bands, and systematic differences are observed in the blue, with higher R rs values estimated by NIR1-SWIR than by Red-NIR (Fig. 12(i)).

Match-up exercise for VNREDSat-1
The present algorithm is now validated based on the restricted match-up VNREDSat-1 data set presented in section 2.1. For that purpose, the in situ remote sensing reflectance values measured at the Venice (4 scenes) and Lucinda (2 scenes) are compared to the NAOMI-R rs (λ) at the center wavelength of NAOMI, using the spectral response function provided in Fig. 2. These latter bands are 488 nm, 565 nm, and 655 nm. Due to the restricted number of match-up data, the statistical results have to be interpreted with caution. However, the similar patterns (about the same RMSD and MPD values) as the one previously described for OLI can be observed, with the best retrieval achieved at 565 nm (Table 6 and Fig. 13(a)). Accounting for the new TOA calibration factors, estimated from the exercise performed in section 4, slightly improves the R rs retrieval, especially in the blue part of the spectrum (Fig. 13(b) and Table 6), This later result should however be interpreted with caution as more match-up data points are needed to confirm these TOA calibration factors.  the Red-NIR slightly under-estimate R rs (λ) in the visible compared to ACOLITE-DSF. If this pattern is confirmed from other images, it may be explained by the NAOMI TOA calibration factors, the adopted spectral relationships (Eq. (19) and Eq. (20)) which do not properly account for such high level of turbidity, and the aerosol transmittance which is fixed to 1.0 in the model. These three aspects should be addressed in the next version of the algorithm. The relatively narrow range of variability of the R rs (λ) spectral values in the Fig. 14 density plots reflects the homogeneous distribution of SPM, mainly due to the calm conditions encountered during the acquisition day of the image (low wind speed, and low river outflow characteristics of the dry season). For this reason, the two algorithms provide very similar R rs (λ) spatial distribution at the three considered wavelengths (not shown). The scatter points around the 1:1 lines are likely explained by the different spatial resolutions of the sensor and the presence of white-caps partly identified by WiPE (Fig. 14(c)).

Concluding remarks
Due to the unavailability of SWIR bands on NAOMI, we developed a new atmospheric correction algorithm (Red-NIR) dedicated to the estimation of R rs (λ) over coastal and inland waters. This iterative algorithm is based on the assumption of spatial homogeneity in the aerosol signal over the whole scene, and on the use of spectral relationships, similarly to the atmospheric correction algorithms of POLDER2 [60] and GOCI [61] between R rs (λ) in the NIR and visible parts of the spectrum. The different results, based on match-up analysis (Figs. 10 and 12), emphasize that the WiPE-Red-NIR algorithm applied to NAOMI is well adapted for the monitoring of R rs from clear to moderately turbid waters. Despite the absence of SWIR band on NAOMI, preliminary results based on the respective application of NIR-SWIR and ACOLITE/DSF [27] on simultaneous images collected by NAOMI and OLI over a more turbid water area, showed that these two algorithms performed similarly. Match-up exercise performed over very turbid waters should however be conducted in a near future to confirm this result, and to verify the performance of Red-NIR over such water areas. Specifically, we should test how the algorithm behaves for waters presenting very different spectral behavior as the one described by Eq. (19) and (20). This study also stresses the need of additional match-up data points for a robust vicarious calibration of NAOMI.