Sentinel-2/Landsat-8 product consistency and implications for monitoring aquatic systems

Sentinel-2 and Landsat data products when combined open opportunities for capturing the dynamics of near- shore coastal and inland waters at rates that have never been possible before. Recognizing the differences in their spectral and spatial sampling, to generate a seamless data record for global water quality monitoring, it is critical to quantify how well the derived data products agree under various atmospheric and aquatic conditions. This study provides an extensive quantitative assessment of how Landsat-8 and Sentinel-2A/B equivalent data pro- ducts compare and discusses implications on differences in downstream products generated via the SeaWiFS Data Analysis System (SeaDAS). These products include the top-of-atmosphere (TOA) reflectance ( ρ t ), the re- mote-sensing reflectance ( R rs ), as well as biogeochemical properties, such as the total suspended solids (TSS). The analyses are conducted a) for Landsat-8 and Sentinel-2A/B near-simultaneous nadir overpasses (n-SNO) and b) over several highly turbid/eutrophic inland/nearshore waters. Following the implementation of vicarious gains for Sentinel-2A, the n-SNO analyses indicated that Landsat-8 and Sentinel-2A agree within±1% in ρ t and±5% in R rs products across the visible and near-infrared (NIR) bands. Similar evaluations with preliminary vicarious gains for Sentinel-2B showed±2% in ρ t and±7% in R rs products. Considering Landsat-8-derived R rs products as a reference, we found<5% difference in Sentinel-2A and -2B R rs products. Analyses of combined TSS and R rs time-series products over several aquatic systems further corroborated these results and demon- strated the remarkable value of combined products. Occasional negative retrievals of R rs products over hyper-eutrophic and highly turbid waters suggest the need for improvements in the atmospheric correction procedure to empower science/application community to fully explore Landsat-Sentinel-2 products. With very similar absolute radiometric observations and products, the science community should consider developments of sui- table biogeochemical algorithms to maximize the utility of merged Landsat-Sentinel-2 products.


Introduction
Coastal waters are under pressure due to direct influence by human activities and climate change. Nearly 50% of the world's population lives in coastal counties, and all coastal cities depend on both a reliable source of high-quality drinking water to ensure human health as well as reliably clean receiving waters to ensure environmental health and to sustain high coastal ecosystem values (Vörösmarty et al., 2000). Further inland, freshwaters are critical sources of drinking waters, major hubs for recreational activities, and prime providers of ecosystem services to local communities. Human exploitations of land and extreme weather T 2017b; Pahlevan et al., 2017c) enabling limnologists, coastal oceanographers, aquatic ecologists, and water resource managers to enhance their monitoring efforts. The 2-3-day revisit time is readily achieved globally when Landsat-9 becomes operational in early 2021. Such frequent revisit times are essential to capture dynamical nearshore coastal and inland waters considering cloud coverage. With four visible channels, i.e., two in the blue and one in each of the green and red regions of the spectrum, and adequate radiometric performances, the OLI and MSI are capable of monitoring certain water quality indicators (Pahlevan et al., 2014;Pahlevan and Schott, 2013). In addition to the four visible bands, MSI has three additional near-infrared (NIR) bands that allow for using alternative algorithms for the retrievals of concentrations of chlorophyll-a (Chl) or other pigments in intense bloom conditions (Gower et al., 2005;Gower et al., 2008a;Moses et al., 2009). This certainly makes MSI more powerful in identifying areas impacted by algal blooms; however, this study primarily focuses on the OLI and MSI's similar (i.e., closely spaced) spectral bands, i.e., the visible bands, one NIR band (i.e., 865 nm), and the two shortwave-infrared (SWIR) channels.
Recent publications have demonstrated the utility of Landsat-8 and/ or Sentinel-2A data for water resource monitoring. Gernez et al. (2017) used Sentinel-2A to study the influence of variations in the concentration of total suspended solids (TSS) and Chl on oyster physiological response, and highlighted the utility of Sentinel-2A's near-infrared (NIR) bands to quantify TSS (Novoa et al., 2017). Manzo et al. (2015) performed a sensitivity analysis on the contributions of different optically active components in the spectral bands of Landsat-8, Sentinel-2, and Sentinel-3 and asserted difficulties in the retrievals of colored dissolved organic matter (CDOM) absorption using Landsat-8 and Sentinel-2 band configurations. Watanabe et al. (2017) attempted empirical methods for the retrievals of Chl in a hyper-eutrophic reservoir in Brazil and highlighted the utility of Sentinel-2A's near-infrared (NIR) bands (Gower et al., 2005). Lavrova et al. (2016) used Landsat-8 and Sentinel-2A imagery to study fine-scale hydrodynamic processes, i.e., internal waves, around river plumes in the eastern Black Sea and Mediterranean region. Similar studies have either attempted algorithm developments (Toming et al., 2016) or demonstrated the applicability of Sentinel-2A data for water quality (Liu et al., 2017) or bottom mapping (Dörnhöfer et al., 2016). Nonetheless, no attempt has been made to demonstrate the suitability of combined Sentinel-2 and Landsat-8 imagery and leverage their improved temporal coverage at 10-60 m spatial sampling (resolution).
To create a consistent, seamless multi-sensor data record for monitoring global aquatic systems, it is vital to ensure that OLI and MSI are consistent at both Level-1, i.e., calibrated top-of-atmosphere (TOA) reflectance (ρ t ), and Level-2 products, e.g., remote-sensing reflectance (R rs ) defined as the ratio of upwelling water-leaving radiance to the total downwelling irradiance just above water. Such exercises have extensively been practiced for ocean colour missions over open-ocean waters (Franz et al., 2005;Maritorena and Siegel, 2005;Morel et al., 2007;Müller et al., 2015) to gauge how differences in TOA measurements propagate to R rs products (Gordon, 1998) and other downstream physical/geophysical variables, including the inherent optical properties (IOPs). This manuscript evaluates and quantifies consistency in OLI, MSI-A, and MSI-B products through a) an extensive evaluation of nearsimultaneous nadir overpasses (n-SNOs) and b) analyses of time-series data records over several inland/coastal waters. While the former provides insights into the cross-mission consistency under ideal atmospheric conditions, i.e., low aerosol loading, the latter permits assessing relative performances under various atmospheric and aquatic conditions. This study regards OLI as the reference sensor to evaluate MSI-A and MSI-B's relative performances. This paper is, therefore, structured as follows. In Section 2, the material and methods are provided followed by the results presented in Section 4, which elaborates on the intercomparison analysis and the analyses of TSS and R rs over several bodies of waters. The discussions and conclusions are presented in Sections 4 and 5.

Materials and methods
This section is comprised of a) a brief description of the atmospheric correction along with vicarious gain estimations, b) the n-SNO analysis for OLI-MSI intercomparisons, and c) the time-series analysis approach. Here, we refer to both Sentinel-2 payloads as MSI but will be specific, i.e., MSI-A and MSI-B, when necessary.

Atmospheric correction
The atmospheric correction implemented in the SeaWiFS Data Analysis System (SeaDAS) has been fully described in various publications (Gordon and Wang, 1994;Mobley et al., 2016;Pahlevan et al., 2017b); however, here, we briefly describe a theoretical background for completeness. The goal is to reduce OLI and MSI scene pairs used in the intercomparison analysis to R rs products, referred to as Level-2 products. Briefly, over bodies of water, the total signal recorded at TOA (ρ t (λ)) is modeled as a summation of different radiometric components as below: where t is the diffuse transmission, ρ r is the Rayleigh reflectance in the absence of aerosol, ρ a is the aerosol radiance, and ρ ar is the radiance arising from Rayleigh-aerosol multiple scattering, and ρ w is the waterleaving reflectance just above water, which can readily be converted to R rs (Mobley et al., 2016). Note that the contributions by whitecaps and sun-glint are discarded. Amongst all the unknown components, estimating the total aerosol contribution is the most challenging task and small errors may lead to high uncertainties in R rs retrievals over coastal and inland waters (IOCCG, 2010). To do so, Rayleigh-corrected radiance in the NIR (865 nm) and SWIR bands are supplied to infer dominant aerosol type and optical thickness. Owing to OLI and MSI's relatively low signal-to-noise ratios in these bands (Del Castillo et al., 2012;Pahlevan et al., 2017b), spatial aggregations are considered in the aerosol correction  to minimize noise in products. Since our processing system currently produces MSI products at 20-m resolution, 6 × 6-and 9 × 9-element averaging windows were applied for OLI and MSI images, respectively. Throughout this study, the aerosol corrections for OLI and MSI processing were conducted using the 864-and~1611-nm bands (see Figs. 4 and 6). It should be noted that MSI-A data collected within the 2015-2017 period are processed using original Relative Spectral Responses (RSRs) (Fig. 1) whereas the data acquired after 2018 are processed using updated RSRs (Clerc et al., 2018). The latest SeaDAS update (i.e., SeaDAS v7.5) applies the updated RSRs.
The processing system allows for vicariously calibrating MSI data Pahlevan et al., 2017b). This is performed using high-quality in situ data collected autonomously at the Marine Optical Buoy (MOBY) (Clark et al., 2003) and the BOUSSOLE (Antoine et al., 2008). Briefly, the atmospheric correction over these sites is conducted to obtain water-leaving radiance and the corresponding atmospheric parameters (Eq. (1)). The estimated water-leaving radiance is replaced with recorded field measurements, which are propagated back to TOA with the estimated atmospheric parameters. The in situ data collected under satellite overpasses within the 2015-2017 period were applied for the vicarious calibration. These gains are slightly different than those provided in Pahlevan et al. (2017b). The differences can be attributed to the availability of more matchups yielding more accurate estimations. The median of all the gains derived for the available images over the BOUSSOLE (n = 13) and MOBY (n = 4) sites were independently computed for MSI-A. The average of the two sites are reported in Table 1. For MSI-B, only four cloud-free images over the MOBY site were identified; hence, the derived gains are treated as preliminary. Due to major updates in MSI-A's 492-nm band implemented on January 18th 2018 (Fig. 1), a unit gain (similar to that for MSI-B) is recommended for the current data.

Intercomparisons at near-simultaneous nadir overpasses (n-SNO)
Landsat-8 and Sentinel-2A/B with equatorial crossing times at 10:00 a.m. ± 15 min and 10:30 AM and a revisit cycle of 16 and 5 days are placed in orbits that allow occasional near-simultaneous collections with < 30 min time difference over multiple locations on the Earth surface. Due to similarities in atmospheric and aquatic conditions, examining products at n-SNO events minimizes uncertainties in the intercomparison analysis. In addition, we choose to utilize scene pairs with low aerosol loading; therefore, differences in products, to a large extent, can be attributed to differences in sensors' absolute radiometric responses. Moreover, to reduce uncertainties in spectral band adjustments (Pahlevan et al., 2017d) inherent to sensor-to-senor intercomparisons, the analyses were limited to moderately eutrophic/turbid coastal waters only.
The TOA intercomparison is carried out for the visible, the 865-nm ( Fig. 1), and the SWIR bands. The n-SNO events were selected automatically by conducting an exhaustive search on the metadata and  (Clerc et al., 2018). The RSRs for the blue band is now centered at 492 nm, which was at 497 nm within the 2015-2017 timeframe. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
The updated and preliminary vicarious gains for MSI-A (n = 17) and MSI-B (n = 4), respectively, derived from in situ observations within the 2015-2017 timeframe at the BOUSSOLE and/or MOBY sites. The values correspond to the median and mean values for each center wavelength (CW). The gains in thẽ 865-nm and SWIR bands were kept unity. ⁎ This gain is applicable to the data collected in 2015-2017 (with the band centered at 497 nm). A unit gain for MSI-A's 492-nm band is recommended for 2018 and onward. limiting the difference in overpass times to < 30 min. This threshold resulted in n-SNOs distributed globally across different latitudes (see Appendix A. for Supplementary data). To perform the analyses, Level-1T OLI and Level-1C (L1C) MSI data were obtained from the United States Geological Survey (USGS) web portal (https://earthexplorer. usgs.gov) and ESA's science hub (https://scihub.copernicus.eu). Given the calibration coefficients in the metadata, the images were converted to unitless ρ t , which were applied in the TOA intercomparisons. In this study, TOA radiance products (L t ) are not analyzed and the vicariously calibrated OLI ρ t data products are treated as the reference (Pahlevan et al., 2017c) to assess MSI ρ t . The viewing zenith and azimuth angles of MSI observations were re-gridded to 10-, 20-, and 60-m resolution to match science pixel grids (Pahlevan et al., 2017b). To minimize the effects of the bidirectional reflectance distribution factor, intercomparisons were performed for pixels with view zenith angles (VZA) within the ± 5°range. Several tests and thresholds were applied to ρ t (1611) to avert intercomparisons under turbid aerosol conditions, in sun-glint regions, and to ensure similar aerosol conditions within the overpass time window (Pahlevan et al., 2017c). Based on previous research (Pahlevan et al., 2014) and after running a few experiments, we set L t (1600) < 0.24 wm −2 sr −1 μ as the threshold, which was converted to ρ t given the metadata information. Furthermore, spatial heterogeneity tests were carried out for 6 × 6-and 9 × 9-element (180 × 180 m 2 ) of OLI and MSI ρ t (560) image products to exclude highly eutrophic or turbid waters. In this study, sensors' point spread functions were assumed squares (Pahlevan et al., 2016). To avoid the large computational burden, for each OLI-MSI-A and OLI-MSI-B image pairs, only top 60 matchup areas (ranked based on low to high ρ t (1611)) within each scene pair were incorporated. An example of the candidate matchup sites is highlighted in Fig. 2. For each matchup, the differences in RSRs (Fig. 1) were accounted for using the method described in Pahlevan et al. (2017c). Briefly, we simulate hyperspectral TOA reflectance (ρ t s ) using the MODerate resolution atmospheric TRANsmission (MODTRAN) code (Berk et al., 2006) for each matchup site. The sun-sensor geometry, per-site ancillary data, including O 3 amount, wind speed, water vapor, NO 2 , and surface pressure were supplied to MODTRAN. In addition, hyperspectral water-leaving radiances required for TOA simulations were determined from OLI products through spectral fitting (Mobley et al., 2005). This was performed at each matchup site assuming in-water and atmospheric optical properties are known. However, recognizing that the aerosol optical properties derived from OLI data processing may carry large uncertainties, the derived aerosol optical thickness at 865nm band, i.e., AOT(865), was scaled up to ± 80% to model a variety of TOA reflectance spectra. Amongst the 80 different simulated TOA cases (i.e., 2 bounding aerosol models × 40 AOT(865)), the curve that best fits the OLI ρ t was identified. The ρ t s spectrum corresponding to the fitted curve was then convolved with the spectral responses of the sensors to compute the differences in RSRs. Adjustments were then applied to MSI-A and MSI-B ρ t observations on a band-by-band basis. The maximum differences were found in the red channels (i.e., 655 versus 664 nm). Out of all the n-SNO events for OLI-MSI-A and OLI-MSI-B, 103 and 24 scene pairs, respectively, had at least 10 pixels passing all the thresholding criteria. For a given scene pair, the median ratios of OLI and MSI ρ t for all the matchups represent the interconsistency metric. The median of all the median ratios of the scene pairs allude to the overall cross-mission consistency. The median AOT(865) for the valid OLI-MSI-A and OLI-MSI-B intercomparisons (estimated as byproducts in SeaDAS) were 0.039 and 0.042, respectively, suggesting, on average, a very low aerosol loading. To evaluate the effectiveness of vicarious gains listed in Table 1, the TOA intercomparisons are provided before and after applying the gains. The intercomparisons of R rs are also performed (independent of the TOA analyses) for the scene pairs and the associated matchups selected above. The MSI ρ t are supplied to SeaDAS two times; once with unity gains applied and once by incorporating gains tabulated in Table 1. To account for the differences in RSRs, we applied a machine learning technique (Pahlevan et al., 2017d), which has been demonstrated to outperform other existing spectral band adjustment methods. This approach uses a trained deep neural network (DNN) to predict MSI-derived R rs from those derived from OLI.

Time-series analysis
The n-SNO analysis reveals product consistency under near-ideal atmospheric/aquatic conditions enabling a rigorous assessment of how instrument calibration differences manifest in R rs products. Yet, it does not address how products compare under challenging atmospheric conditions and in aquatic systems with diverse optical regimes. Thus, we perform a time-series assessment of relative differences in OLI-and MSI-derived products, i.e., R rs and TSS, from mid 2015 when Sentinel-2A became operational. Here, no spectral band adjustment is implemented for R rs intercomparisons. For the temporal analyses of TSS products, we use a semi-empirical approach (Nechad et al., 2010), which has been adopted by the optical oceanography community since the mid-90's and has been found to work well, even with satellite sensors that have not been designed for water retrievals (Neukermans et al., 2009). Note that, here, the relative variations in products are of interest and their accuracy/precision is beyond the scope of this research.
Nine different aquatic systems ( Fig. 3) with various environmental conditions were selected and processed to R rs , from which TSS products were derived. Some of these systems (e.g., Apopka Lake in Florida) have occasionally experienced harmful algal blooms. Other lakes, reservoirs, or bays include Upper-Klamath reservoir (Oregon, USA), Lake Tana (Ethiopia), Lake Mohave (Nevada, USA), San-Francisco Bay (California, USA), Ratzeburger See (Germany), Lake Kyoga (Uganda), Lake Bangweulu (Zambia), and Phewa Lake (Nepal). San Francisco Bay and Lake Tana are mostly sediment dominated and some others like Ratzeburger See are considered mesotrophic waters. Upper Klamath Lake, a large and shallow eutrophic freshwater lake east of the Cascade Range, has been experiencing frequent blue-green algae blooms starting late spring into the summer season (Snyder, 2018). Lake Mohave is a reservoir on the Colorado River in Nevada and Arizona and provides recreational opportunities. In spring and summer of 2014-2015, the regional/local stakeholders issued health advisory for blue-green algae in the lake (Blunt et al., 2018). Lake Bangweulu is a shallow turbid freshwater system which is under pressure due to overexploitation and the changing climate. The lake water level changes drastically from wet to dry season (Kolding and van Zwieten, 2012). In addition, these sites span a wide range of atmospheric conditions ranging from various aerosol types (e.g., rural and urban) to different humidity levels.

Performance metrics
The differences in products are expressed in terms of the root-meansquare differences (RMSD): where λ i indicates band i. For the n-SNO analyses, N is the total number of scene pairs (Tables 2 & 3); whereas for the time-series analyses N is the total number of OLI-MSI same-day overpasses. The median difference (MD) is also calculated as below: The median is preferred over mean due to presence of noise in the analyses. These metrics along with slope and intercept allow us to   gauge interconsistency in OLI and MSI products.

Results
The results are presented in three subsections, including the intercomparisons of ρ t and R rs at n-SNO events (Sections 3.1 and 3.2) and the time-series analyses of R rs and TSS products over lake and coastal waters (Section 3.3). These assessments will collectively provide insights into how consistent Landsat-8 and Sentinel-2A/B products are when processed through SeaDAS.

Top-of-atmosphere reflectance (ρ t )
The intercomparisons of MSI-A and OLI ρ t products before (left panel) and after applying vicarious calibration gains to MSI-A products are shown in Fig. 4. The data points represent the median ratios of MSI-A to OLI (n = 103) and error bars indicate one-standard deviation. It is inferred that OLI and MSI-A products "as-is" (i.e., no gains applied to MSI-A ρ t ) agree to within 3% in the blue bands and to 1% in the green and red channels.
The NIR bands are found very consistent with < 0.1% difference. Following the vicarious calibration, the differences in ρ t in the 443-and 497-nm bands have significantly been reduced, increasing consistency in merged products. However, there is also 0.8% bias remaining in the red channel, which is likely due to uncertainties in the RSR adjustments in the red channels. A closer look at Fig. 1 indicates that MSI-A and OLI red channels do not fully overlap in this part of the spectrum where Chl absorption spectrum may significantly vary (Pahlevan et al., 2017d).
Thus, small errors in modeling this component results in inaccurate spectral band adjustments performed for ρ t products. Similar analyses for the OLI and MSI-B are illustrated in Fig. 5. The plot indicates that the two sensors are consistent within ± 2.5% (left panel) with the NIR bands agree the most, i.e.,~0.2%. The application of preliminary vicarious gains reduces the differences in the 442-and 559-nm channels; whereas, the differences remain nearly unchanged in the green and red channels.
By comparing Figs. 4 & 5 and the vicarious gains (Table 1), it is inferred that MSI-A and MSI-B differ by~1% in the visible bands except for the~492-nm band. Furthermore, the TOA observations differ on the order of~2% in the 704-and~740-nm while they agree well in both the 783-and the~865-nm channels (Tables 1 & 2). Further statistics for the intercomparisons are provided in Table 2 where average ρ t and L t are also included (i.e., t and L t ). The largest differences are found in the 492-nm band for both intercomparisons. In Section 3.2, we will further evaluate these differences for R rs products to better understand the uncertainties. Included in Table 2 are the differences in the SWIR bands for which no vicarious calibrations were attempted. It is found that the OLI and MSI's SWIR bands agree within ± 3%. Such small differences in the~1611-nm channels used in the aerosol removal are not expected to yield notable differences in R rs products (Pahlevan et al., 2017a). This interpretation also applies to relative differences in the MSI-A and MSI-B products.

Remote-sensing reflectance
The median ratios of vicariously calibrated MSI-and OLI-derived R rs Fig. 4. The intercomparison of OLI and MSI-A ρ t products before and after applying vicarious corrections (Table 1) for the visible and the 865-nm bands. The data points represent median values of MSI-A to OLI ratios computed from 103 scene pairs at n-SNO events. The error bars represent one-standard deviation. Note that the NIR band has not been vicariously adjusted. ⁎ The analysis associated with the 497-nm band applies to the 2017-2015 timeframe. Following the change in RSRs, this channel is centered at 492 nm. Fig. 5. The intercomparison of OLI and MSI-B ρ t products before (left) and after applying vicarious corrections using preliminary gains (Table 1) for the visible and the 865-nm bands. The data points represent median values of MSI-B to OLI ratios computed from 24 scene pairs collected at n-SNO events. The error bars represent one-standard deviation. Note that the NIR band has not been vicariously adjusted. products are shown in Fig. 6. Overall, the average differences in products are estimated to range from 1 to 7% indicating relatively consistent products. In particular, R rs (443) products agree to within < 1.5%. This is partially contradictory to the results shown for MSI-B's ρ t (442) (Fig. 5) where MSI-B was found~1.5% larger in magnitude than that of OLI. In principle, this difference is anticipated to translate to 15% difference in R rs (443) over open oceanic waters (Gordon, 1997). This unexpected deviation from theoretical realizations can be attributed either to uncertainties arising from differences in the spectral band adjustments (Pahlevan and Schott, 2012) or to the lack of rigor in the statistical analyses (n = 24).
Similar evaluations apply to the 492-nm bands where > 1.5% and < 4% differences were found in OLI-MSI-B ρ t and R rs , respectively. The differences in the red channels, however, comply with theoretical expectations, i.e., OLI-derived R rs products are larger in magnitude. The one-standard-deviation (i.e.,~2%) shown in forms of error bars in Fig. 6 indicate a robust analysis in Level-2 products. Furthermore, MSI-A and MSI-B products are found to differ < 2% in the visible bands except in the green channel where~6% difference exists. It should also be noted that our recently developed method (Pahlevan et al., 2017d) accounting for differences in spectral bands showed improvements when compared to the spectral matching method implemented in (Pahlevan et al., 2017c). More details on the relative differences and mean R rs (R rs ) are listed in Table 3. The RMSDs and MDs for both intercomparisons are well below acceptable uncertainties for ocean colour products, i.e., 5e −4 , (Del Castillo et al., 2012) except in the 492nm band for which uncertainties due to spectral band adjustments may be higher. Moreover, across-track nonuniformity may contribute to increasing the variability and median differences. The deviations from the unity in the slopes observed for the red channels are, in general, common due to narrow dynamic ranges in R rs in this region (Zibordi et al., 2009). Further matchups are required to provide more robust analyses.

Time-series analyses
In this section, we demonstrate the utility of combined OLI-MSI products by assessing time-series of TSS and vicariously corrected R rs products over several lakes, reservoirs, and bays under various atmospheric and climatic conditions. Figs. 7, 8, and 9 illustrates the timeseries plots exhibiting valid retrievals beginning from July 2015. These sites show various seasonal variability and dynamic ranges in R rs (665) and TSS. The star symbols mark a same-day overpass, i.e., < one-hour time difference, at these sites. The temporal data are extracted from arbitrarily select areas indicated with a star on the TSS maps provided on the left panel. The extraction was conducted by taking the median over 6 × 6-and 9 × 9-element windows (equivalent of 180 × 180 m 2 ) from OLI and MSI data, respectively. This is to minimize random noise and/or artifacts. Pixels adjacent to land were avoided in data extraction as frequent negative retrievals were found (Santer and Schmechtig, 2000). The analyses over various aquatic systems characterized with different TSS and R rs (665) magnitudes highlight the value of combined Landsat-8 and Sentinel-2 data products. It is clearly seen that MSI data fill the gap in OLI data products available every 16 (or 8) days. For instance, there is a major gap in OLI products during the wet season over the central basin of Lake Bangweulu (Fig. 7). This has been partially filled by MSI-A data products. Furthermore, Fig. 7 indicates that OLI produces slightly larger TSS products than those generated from MSI-A. This is, however, the case in 4 of the 7 same-say overpasses; the other three samples agree. The triplet of the sensors clearly capture the seasonal variability in Lake Apopka with peaks in TSS in the month of July (Fig. 8). Relatively frequent OLI, MSI-A, and MSI-B products towards the end of 2017 and in 2018 further demonstrate Landsat-Sentinel-2 capability in monitoring global inland and nearshore coastal waters. In San Francisco Bay (Fig. 8), a Gaussian-like feature is also captured in the winter-spring of 2017-2018. Likewise, a sharp peak is observed in Lake Kyoga in July 2017 (Fig. 7). These two features would have not been fully captured with each individual sensor. Another gap in time-series data occurs in the winter of 2015-2016 in Upper Klamath Lake, where OLI-MSI-A data products are very scarce due to cloud coverage. This gap has partially been filled in with MSI-B products in 2017. The differences in OLI-MSI-A and OLI-MSI-B are further highlighted with OLI giving out larger signals. High suspended solid loads are also captured in time-series plots of Upper Klamath Lake and Lake Mohave in January 2018 and July 2016, respectively. The quantitative analyses of interconsistency in OLI-and MSI-A-derived TSS and all R rs products are provided in Table 4 (OLI-MSI-B matchups were inadequate for statistical analyses). The RMSDs, MDs (Eqs. (2) and (3)) and mean R rs (R rs ) associated with the intercomparisons are listed. Only 19 valid matchups (2015-2017) associated with same-day overpasses were identified for this analysis.
Similar to the results presented in Section 3.2 (Fig. 6), these analyses verify that, on average, the MSI-A-derived R rs in the 497-and 665nm bands are slightly lower in magnitude than those of OLI. The difference in the red channel, however, is relatively small (i.e., < 5%) compared to that in the 497-nm, which is~9% and is much larger than what was reported in Fig. 6. While the MD and RMSD in the 443-nm channel are negligible (i.e., < 2%), the differences in the green channel is relatively high in magnitude. For these mostly turbid waters with high peaks in the green ( = R sr (559) 0.0143 [1/ ] rs ), these differences may  Table 1). The data points and the error bars represent median statistics and one standard deviations for all the available matchups, i.e., n-SNO events. The MSI products agree with those of OLI within ± 7% in all the visible bands. The differences in MSI-A and MSI-B R rs products are on the order of 2% except in the green channel where 5% difference in observed. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) not be significant (i.e., < 6%). The MSI-A's lower TSS (i.e., < 0.3 g m −3 ) is attributed to the differences in R rs (665). Further research must be dedicated to assessing other TSS algorithms or to developing new algorithms to reduce the RMSD (i.e., 0.41 g m −3 ) over a larger number of sites. The TSS algorithm fails when R rs (665) retrievals are invalid. Occasional negative R rs retrievals (from both OLI and MSI) in hyper-trophic and turbid waters further emphasize the need for improved aerosol removal in areas where NIR water-leaving radiance is relatively high (Bailey et al., 2010). For instance, the retrievals over Lake Tana were mostly invalid (not shown here). Invalid retrievals in Lake Kyoga, January 23 rd 2016 Lake Bangweulu, July 23 rd 2016 Fig. 7. Examples of TSS distributions together with time-series plots of an arbitrary location (marked with a filled star) in Lake Banglweulu and Lake Kyoga. The products shown are R rs (665) and TSS for OLI (blue circles), MSI-A (red triangle), and MSI-B (green triangle). The same-day OLI-MSI-A overpasses (denoted with solid and outlined stars) indicate that the OLI-derived products are larger in magnitude than those of MSI-A. The retrievals for same-day overpasses exhibit larger-thannormal magnitudes in both R rs (665) and TSS. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Apopka Lake, April 7 th San Francisco Bay, June 14 th 2017 Fig. 8. The time-series plots shown for an arbitrary location (marked with a filled star) in Apopka Lake and San Francisco Bay. Examples of TSS distribution are shown to the left. The relatively frequent data derived from the triplet of the sensors capture low-frequency seasonal variability for the two water bodies. The timeseries is more cluttered for San Francisco Bay due to strong influence of tidal forcing. (665) were also observed for Phewa Lake and Ratzeburger See. The Ratzeburger See, in particular, represents relatively clear waters where TSS < 0.5 g m −3 . Phewa Lake showed similar results with occasional high TSS loads. In such bodies of waters, R rs retrievals are often influenced by image artifacts or residual calibration biases.

Discussions and future directions
Throughout this study, we have shown that OLI and MSI products, including ρ t , R rs , and TSS, are sufficiently consistent. Our approach to direct inter-comparisons of ρ t products has proven to provide robust results. For example, our findings for ρ t (443) products agree well with the corresponding results provided in . The residual differences in ρ t after vicarious calibration (in particular for OLI and MSI-B) can be improved when more Landsat-8 and Sentinel-2 data become available over the MOBY and BOUSSOLE sites. This will directly enhance consistency in R rs and downstream products. Further improvements may involve improving the spectral band adjustments in the blue and red portion of the spectrum in both ρ t and R rs domains. Some of the uncertainties in direct intercomparisons also arise from potential MSI cross-track non-uniformities (Pahlevan et al., 2017c). Here, the analyses were based upon image products processed via SeaDAS. The heritage physics-based, atmospheric correction algorithm built into SeaDAS with well-known strengths and sources of uncertainties (IOCCG, 2010; Mobley et al., 2016; Pahlevan et al., 2017a) enables straightforward predications of how image artifacts are translated to R rs products. Further, the relatively small differences in OLI and MSI's ρ t (NIR) and ρ t (SWIR) ( Table 2) are expected to make minimal contributions to the differences in R rs . Nonetheless, it would be valuable to evaluate the relative consistency in the OLI-MSI products processed through other available processing systems. This is, in particular, required for CDOM-rich or hyper-eutrophic waters.
The combined TSS and R rs (665) time-series products at various sites further emphasize the consistency in Landsat-8 and Sentinel-2 missions, however, research has to continue to enable retrievals in extremely eutrophic/turbid conditions and/or CDOM-rich waters where SeaDAS retrievals often fail. From the statistics tabulated in Table 4, it is found that a median difference on the order of 3.8E-4 [1/sr] in R rs (665) results in 0.27 g m −3 (Nechad et al., 2010). This difference in magnitude is equivalent of~10% difference in R rs (665). As pointed out, improving the vicarious gains in this channel can diminish these discrepancies, which may be regarded negligible for applications in more turbid waters (with R rs (665) > 0.0096 [1/sr]). Alternatively, one may use the 865-nm channel for TSS retrievals where OLI, MSI-A, and MSI-B are most consistent, i.e., < 0.3% in ρ t , (Fig. 5). Although vicarious gains provided for MSI-A and MSI-B's 704-and 740-nm bands (Table 1) are anticipated to enhance interconsistency in corresponding products, further research is required to extensively assess how they compare for various applications.
Future studies should focus on improvements in the retrievals of Chl or IOPs. In general, for a consistent product record for a Landsat-Sentinel-2 constellation, it is vital to use consistent algorithms. Given MSI's additional NIR bands, which permit more alternatives for Chl retrievals in hyper-eutrophic systems (Gower et al., 2008b), there must be a compromise in product consistency. Perhaps, MSI-derived Chl record can be produced independent of OLI over hyper-eutrophic waters while another set of Chl product record can be produced from both OLI and MSI using consistent algorithms.

Summary and conclusion
This study demonstrates the consistency in Landsat-8 and Sentinel-2A/B aquatic products provided via NASA's processing software suite, Upper Klamath Lake, May 7 th 2018 Lake Mohave, June 22 nd 2015 Fig. 9. The time-series plot shown for an arbitrary location (marked with a filled star) for Upper Klamath Lake and Lake Mohave. Examples of TSS distribution are shown to the left. A pronounced peak in Lake Kyoga is captured by OLI and MSI-A in the June-July 2017 time-frame. Major gaps are also found in data products earlier period (i.e., 2016) when MSI-A data had not gone into full global acquisition mode.

Table 4
The intercomparison metrics for OLI and MSI-A products over the sites in which valid TSS [gm 3 ] retrievals were available. (July 2015-January 2018). The results are consistent with the findings in Fig. 6 i.e., SeaDAS. We showed that the OLI and MSI data when vicariously calibrated agree to within ± 7% in R rs products, which are central to generating higher level aquatic products. This demonstration performed at near-simultaneous nadir overpasses under near-ideal atmospheric conditions was further verified over select highly turbid/eutrophic waters. The time-series analyses showed that the total suspended solids (TSS) (Nechad et al., 2010) derived from OLI and MSI were, on average, in good agreements, i.e., ∆TSS = 0.27 g m −3 , with MSI producing lower TSS products. While there is more data needed to evaluate MSI-B, the initial analyses of OLI-MSI-A and -MSI-B suggests that MSI-A-and MSI-B-derived R rs products are well in agreements in the visible bands (i.e., < 6%). Further investigations are required to affirm these differences. With two-to-three-day revisit time enabled by combined Landsat-8/9 and Sentinel-2A/B data, the aquatic science and end-user community can benefit from high-quality and consistent products. While the science community continues to develop algorithms to improve products, it is now possible to exploit this virtual constellation for operational puposes.