Evaluating Atmospheric Correction Algorithms Applied to OLCI Sentinel-3 Data of Chesapeake Bay Waters

Abstract

Satellite remote sensing permits large-scale monitoring of coastal waters through synoptic measurements of water-leaving radiance that can be scaled to relevant water quality metrics and in turn help inform local and regional responses to a variety of stressors. As both the incident and water-leaving radiance are affected by interactions with the intervening atmosphere, the efficacy of atmospheric correction algorithms is essential to derive accurate water-leaving radiometry. Modern ocean color satellite sensors such as the Ocean and Land Colour Instrument (OLCI) onboard the Copernicus Sentinel-3A and -3B satellites are providing unprecedented operational data at the higher spatial, spectral, and temporal resolution that is necessary to resolve optically complex coastal water quality. Validating these satellite-based radiance measurements with vicarious in situ radiometry, especially in optically complex coastal waters, is a critical step in not only evaluating atmospheric correction algorithm performance but ultimately providing accurate water quality metrics for stakeholders. In this study, a regional in situ dataset from the Chesapeake Bay was used to evaluate the performance of four atmospheric correction algorithms applied to OLCI Level-1 data. Images of the Chesapeake Bay are processed through a neural-net based algorithm (C2RCC), a spectral optimization-based algorithm (POLYMER), an iterative two-band bio-optical-based algorithm (L2gen), and compared to the standard Level-2 OLCI data (BAC). Performance was evaluated through a matchup analysis to in situ remote sensing reflectance data. Statistical metrics demonstrated that C2RCC had the best performance, particularly in the longer wavelengths (>560 nm) and POLYMER contained the most clear day coverage (fewest flagged data). This study provides a framework with associated uncertainties and recommendations to utilize OLCI ocean color data to monitor the water quality and biogeochemical dynamics in Chesapeake Bay.

1. Introduction

At the top of the atmosphere, ocean observing satellites measure electromagnetic radiation reflected from the surface layers of water. Embedded in these measurements is water leaving radiance (L_w) which is regulated by the in-water interactions of incident irradiance (E_d) and the concentration and type of optically active constituents in the water [1]. Remote sensing reflectance (R_rs = L_w/E_d) can therefore be used to derive optical properties of the water that include phytoplankton pigment concentrations [2], dissolved organic matter [3] water clarity [4], and total suspended matter [5,6]. These data can significantly enhance water quality monitoring and provide a better understanding of global and local environmental processes in a changing climate.

Since the radiation source (i.e., the sun) and the satellite detector are both ‘remote’, the radiance incident on a given sensor is affected not only by the reflection of light from within and off the surface of water, but also by interactions with the intervening atmosphere. It is widely recognized that 90–99% of the total signal measured by a satellite sensor is due to confounding influences from the atmosphere and varying water surface facets [7]. Accurately removing these influences, termed atmospheric correction, enables the science of aquatic remote sensing.

Top of atmosphere radiance (L_t) can be partitioned into distinct contributions from the atmosphere and ocean surface as shown in Equation (1). L_r(λ) is radiance due to Rayleigh scattering by non-absorbing air molecules, L_a(λ) is radiance due to scattering by aerosols, L_ra(λ) is the multiple interaction term between scattering by air molecules and aerosols, L_wc(λ) and L_g(λ) are radiance due to whitecaps and specular reflection of sunlight off the water surface (e.g., glint), respectively, t(λ) and t₀(λ) are diffuse transmittances of the atmosphere from the sun to the surface and from the surface to the sensor, respectively, T(λ) is the direct transmittance from the surface to the sensor, and cosθ₀[L_w(λ)]_N is the normalized water leaving radiance where θ₀ is the solar zenith angle. Note that path geometry (i.e., solar and viewing zenith and azimuth angles, θ, Φ) is not shown for brevity, and λ represents wavelength.

$$\begin{gathered} {\mathrm{L}}_{\mathrm{t}}\left(\mathsf{\lambda }\right)={ \mathrm{L}}_{\mathrm{r}}\left(\mathsf{\lambda }\right)+{\mathrm{L}}_{\mathrm{a}}\left(\mathsf{\lambda }\right)+{ \mathrm{L}}_{\mathrm{ra}}\left(\mathsf{\lambda }\right)+{\mathrm{tL}}_{\mathrm{wc}}\left(\mathsf{\lambda }\right)+\mathrm{T}\left(\mathsf{\lambda }\right){\mathrm{L}}_{\mathrm{g}}\left(\mathsf{\lambda }\right) \\ +\mathrm{t}\left(\mathsf{\lambda }\right){\mathrm{t}}_0\left(\mathsf{\lambda }\right){\cos \mathsf{\theta }}_0{\left[{\mathrm{L}}_{\mathrm{w}}\left(\mathsf{\lambda }\right)\right]}_{\mathrm{N}} \end{gathered}$$

(1)

Atmospheric correction schemes aim to solve a set of deterministic radiative transfer models that resolve the individual terms in Equation (1) from L_t (λ) measurements [7,8]. Most of the partitioned radiances can be estimated with high confidence using standard radiative transfer theory and incorporating ancillary data (i.e., ozone, sea surface wind speed, atmospheric pressure measurements) [8]. In coastal waters, the main challenge is determining the aerosol contribution (L_a + L_ra) since aerosol concentrations, chemical composition, and morphology can be highly variable [7,8,9]. Heritage atmospheric correction schemes use aerosol models [9,10] to compute L_a + L_ra as a function of the aerosol optical thickness for a variety of sensor-sun geometries [7]. This process relies on the “black pixel assumption” which assumes L_w in the near-infrared (NIR) is negligible and entirely caused by atmospheric effects [11,12]. Spectral variation at two NIR wavebands is used to aid in the selection of an aerosol model and computed radiances are extrapolated through the visible bands [7,8,11,13,14]

To date, ocean color remote sensing has been successfully exploited in open ocean waters whose optical properties are primarily determined by phytoplankton and co-varying inherent optical properties. Atmospheric correction algorithms based on the heritage Gordon (1978) approach have been used to operationally process data from most ocean color satellite sensors; however, the application of such methods has been inherently more challenging in coastal waters. While the assumption of zero NIR L_w is typically valid in the open ocean where NIR absorption by water dwarfs particle scattering, it is often invalid in productive and turbid waters where concentrations of light-scattering particles are high [12,15,16]. Thus, the process of extrapolating aerosol path radiance into visible bands can result in an overestimation of aerosol scattering, negative reflectance at visible wavelengths, and ultimately erroneous satellite-derived water quality products [15,17,18].

Modifications to the heritage atmospheric correction scheme have been developed to improve atmospheric correction in coastal waters. The utilization of longer wavelengths (i.e., shortwave infrared, SWIR) when present on sensors can be used to estimate and remove atmospheric signals since the absorption of water is greater and water-leaving radiance is smaller in that region [19,20,21,22,23,24]. Other modifications to the heritage correction scheme include applying a spatially constant aerosol type, utilizing water reflectance in the red and NIR [25,26], using a bio-optical model to conduct an iterative estimation of non-zero NIR [27,28], and exploiting differences in the spectral shape of aerosols and water reflectance [29]. Alternative deterministic and statistical approaches aimed at inverting L_t in a single step have also been proposed, such as spectral optimization and neural network techniques [30,31,32,33,34]. Some of these approaches are developed, trained, and validated with regional datasets that may perform erroneously in other coastal water bodies.

In this study, the performance of four freely available atmospheric correction processors applied to data collected from 2016 to 2020 by the Ocean and Land Colour Instrument (OLCI) onboard the Sentinel-3 satellites are evaluated against in situ radiometric measurements collected in the optically complex waters of the Chesapeake Bay. OLCI currently contains the most spectral bands and the highest spatial resolution (300 m) of the current suite of missions designed for ocean color. The constellation of OLCI also provides nearly daily coverage of Chesapeake Bay, making it an ideal sensor to study complex, small-scale dynamic processes typical of a coastal water body.

Several studies have assessed the efficacy of various atmospheric correction schemes applied to OLCI data collected in regional coastal waters by evaluating satellite retrievals against in situ R_rs measurements [23,35,36,37,38,39,40]. Results from these studies demonstrate that the selection of a best performing atmospheric correction algorithm varies due to the variability of optical properties found in coastal waters globally. For example, Mograne et al. (2019) and Vanhellemont and Ruddick (2021) applied the same atmospheric correction algorithms in French and Belgian coastal waters, respectively, and obtained contrasting performance results likely due to differences in turbidity concentrations. The present study contributes to this body of literature by assessing atmospheric correction algorithms to the coastal and estuarine waters of Chesapeake Bay.

2. Materials and Methods

2.1. Study Area

The Chesapeake Bay is a highly productive estuary located in the Mid-Atlantic region of the United States. The Bay is a shallow, partially mixed, temperate estuary spanning ~300 km in length and receives inputs of freshwater from an extensive 165,000 km² watershed [41]. The Bay follows a salinity gradient, with the Upper and Lower Bay largely oligohaline and polyhaline, respectively. The winter-spring freshet from the Susquehanna River typically regulates the spring phytoplankton bloom [42], while other non-algal constituents correlate with river discharge [43]. Water properties can also be influenced by estuarine circulation, frontal features, sediment resuspension, and tidal cycles. Consequently, optical properties of the Chesapeake Bay waters are strongly influenced by varying concentrations of phytoplankton, total suspended matter (TSM), and chromophoric dissolved organic matter (CDOM) [44]. R_rs at 865 nm (NIR wavelength) demonstrates the wide optical variability and can be used to infer a TSM gradient in tributaries (Figure 1).

2.2. In Situ Radiometry

Hyperspectral in situ radiometry was collected at a total of 172 stations in the Upper and Middle Bay and Choptank River, a tributary located on the eastern shore. Data were collected in different seasons and conditions between 2016 and 2020 using two methods (Figure 1). One method used a set of intercalibrated hyperspectral radiometers (TriOS RAMSES, Rastete DE) deployed on a mini-catamaran float platform [45] with a skylight-blocked approach [46,47,48]. The float held a downwelling irradiance sensor (E_d) and an inverted radiance sensor equipped with a plastic cone that extended fully beneath the air-water interface while the sensor fore-optics remained in the air, effectively blocking the reflected skylight contribution (L_w). The TriOS radiometers collect 256 wavelength bands at 3.3 nm intervals within the calibrated 320–950 nm range. This method collected L_w and E_d at 151 locations; in situ R_rs was calculated as L_w/E_d. The median R_rs was calculated from replicate measurements and any spectrum deviating more than ±10% was excluded in order to eliminate outliers due to variable measurement and illumination conditions.

The other method involved collecting underwater radiometry based on vertical profiling and propagating measurements of upwelling radiance across the water-air interface. Hyperspectral measurements were collected using a Satlantic HyperPro II (Satlantic is now Sea-Bird) free-falling profiler equipped with hyperspectral radiometers (Satlantic/Sea-Bird HyperOCRs) at 21 locations located in the Upper Bay [6,49]. The profiler was equipped with a downward directed radiance sensor to measure upwelling radiance (L_u) and an inverted upward irradiance sensor (E_d). An upward directed irradiance sensor was mounted on a gimbal at the top of the vessel to measure reference water surface irradiance (E_s) and avoided any structure shadow. The HyperOCR radiometers collect 256 wavelength bands at 3.3 nm intervals within the calibrated 350–800 nm range. The mean value of replicate casts was taken at each station. L_w was calculated just above the surface using Equation (2) where ⍴(λ,θ) is the Fresnel reflectance index of seawater and was set to 0.021 and n²_w is the Fresnel refractive index of seawater set to 1.345.

L_w(0⁺, λ) = L_u(0, λ) ∗ [(1 − ⍴(λ,θ))/n²_w(λ)]

(2)

All in situ hyperspectral measurements were interpolated at 1 nm intervals and spectrally weighted to match either the Sentinel-3A and 3B OLCI spectral response functions (OLCI Spectral Response Function Data. Accessed online 13 March 2022, https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-3-olci/olci-instrument/spectral-response-function-data), depending on satellite matchup (see Section 2.5). Due to differing approaches currently considered for the bidirectional reflectance distribution reflectance (BRDF) in coastal waters, in situ R_rs data were not BRDF corrected (see Section 2.4). Solar zenith angles were on average 33.7° (±14.3°) during in situ data collection.

2.3. Satellite Data Collection

The OLCI sensor onboard the European Commission Copernicus programmes Sentinel-3A (launched in February 2016) and Sentinel-3B (launched in April 2018) was built based upon the heritage of the Medium Resolution Imaging Spectrophotometer (MERIS). OLCI is a multispectral push-broom imaging spectrometer with five camera modules that collect top-of-atmosphere (TOA) radiance in 21 wavelengths in the 400 to 1020 nm range at a 300 m spatial resolution [50]. Both satellites sun-synchronously orbit at a mean altitude of 815 km, providing global coverage approximately every two days [51]. OLCI provides significant improvements to MERIS including an increase in the number of spectral bands, improved signal-to-noise, sun-glint mitigation through camera tilt, and higher spatial, temporal, and radiometric resolution [50].

OLCI imagery overlapping the Chesapeake Bay region (36.8, −77.35, 39.65, −75.6) from April 2016 to April 2021 were obtained from various data portals (Copernicus Online Data Access REProcessed, EUMETSAT Data Center Long-Term Archive, and Copernicus Online Data Access). A total of 2152 scenes of TOA, full resolution (EFR), non-time critical Level-1 data and 2119 scenes of full resolution (WFR), non-time critical Level-2 data were downloaded and stored locally. Full Level-1 and Level-2 OLCI scenes were clipped to the Chesapeake Bay region using the Sentinel Application Platform’s (SNAP) Graph Processing Tools (GPT, ‘Subset’).

2.4. Atmospheric Correction Algorithms

Four atmospheric correction algorithms were evaluated against in situ measurements: The Baseline Atmospheric Correction (BAC) v7.0, the Case 2 Regional Coast Color (C2RCC) v7.0, the Polynomial-based algorithm applied to MERIS (POLYMER) v4.12, and NASA’s standard Level-2 generator atmospheric correction algorithm (L2gen) v7.5.3 (Figure 2). Downloaded Level-2 data were pre-processed with BAC, and Level-1 data were locally processed with the other atmospheric correction algorithms. Additional information on processing parameters for each algorithm can be found in Figure 2 and Tables S1–S3. Due to differing approaches currently considered for the BRDF in coastal waters and following methods in similar optical complex R_rs matchup studies [23,40], BRDF correction was excluded from all atmospheric correction schemes.

The BAC algorithm is based on the basic principles for the MERIS atmospheric correction scheme which was designed for open ocean waters and uses a multiple scattering algorithm to assess the aerosol contribution [7,14]. BAC integrates a Bright Pixel Atmospheric Correction (BPAC) to account for when NIR water leaving radiance is not negligible, such as in high scattering waters [29,51,52]. A detected NIR signal partitions the TOA reflectance into components due to aerosols and TSM using a coupled hydrological and atmospheric model [52]. Validation of this approach was based on MERIS match-ups with in situ data provided by the MERIS Matchup In situ Database (MERMAID) which includes compiled global in situ radiometry data [52,53].

The C2RCC algorithm, originally developed by Doerffer and Schiller (2007) for MERIS, is based on an artificial neural network (ANN) inversion model to derive water leaving radiance and in-water optical properties. The ANN is trained on a large database of simulated angle-dependent water leaving reflectances and TOA radiances generated by the radiative transfer model HydroLight [54]. TOA radiances at fifteen wavelengths are corrected for absorbing gasses and the smile effect and used together with ancillary data (i.e., wind speed, salinity, temperature, viewing geometry) to perform an inversion of water-leaving radiance and atmospheric parameters [31,55]. C2RCC v7.0 was applied locally to Level-1 OLCI data using SNAP GPT. C2RCC v7.0 is updated with new default neural nets (i.e., NNv2) that fully supersedes the older neural net (Sentinel-3 Product Notice. Accessed 13 March 2022, https://www.eumetsat.int/media/48139). The most updated system vicarious calibration (SVC) gains were applied (Table S2) and the derived retrievals were not corrected for BRDF effects.

The POLYMER algorithm developed by HYGEOS (Lille, France) was also initially developed for MERIS and consists of a full-spectrum coupled spectral matching algorithm [33]. The algorithm relies on a bio-optical reflectance model of Park and Ruddick (2005) and can be modified with known chlorophyll a concentrations and particle backscattering coefficients to represent various oceanic and coastal waters [56]. A polynomial expression is used to represent the coupled atmosphere and water-leaving reflectance and an iterative process is used to optimize parameters to obtain the best spectral fit of water leaving reflectance [33,56,57]. This algorithm determines aerosol contributions from a linear combination of reflectance terms, rather than from a specific aerosol model, and can characterize complex atmospheric and surface effects such as residual sun glint. POLYMER was trained on in situ data collected from the North Sea [58] and has been validated using global AERONET-OC sites [57]. In this study, POLYMER was locally processed using Python v3.7. All POLYMER defaults were kept, including the hard-coded SVC gains (Table S2); however, BRDF correction was disabled (Table S1).

L2gen is the standard atmospheric correction algorithm that the NASA Ocean Biology Processing Group (OBPG) employs which builds on several decades of atmospheric correction research [9,11,27,28]. The current version of L2gen consists of using precomputed radiative transfer simulations, models, and ancillary information in order to correct TOA radiance for atmospheric terms such as gaseous transmittance, white caps, Rayleigh scattering, and sun glint [8]. The aerosol contribution is estimated using two NIR wavelengths where water leaving radiance can be accurately extrapolated to the visible wavelengths [28]. For waters with non-negligible water-leaving reflectance in the NIR, an optical model is used to estimate NIR and an iterative bio-optical modeling approach is used to extrapolate aerosol contribution to the visible wavelengths [28]. In this study, L2gen was locally processed using Python v3.7 with the L2gen aerosol mode option selected as multi-scattering with 2-band model selection and iterative NIR correction using OLCI wavelengths 865 nm and 1012 nm (Table S1). SVC gains were applied according to OBPG recommendations and BRDF effects were not corrected for (Tables S1 and S2).

All atmospherically corrected products were flagged using recommended quality flags for each atmospheric correction processor (Figure 2, Table S3). The Identification of Pixels (IdePix) pixel identification tool was run during the processing of C2RCC and associated flags were applied to all data for comprehensive image masking (Figure 2, Table S3). IdePix is a SNAP operator that calculates pixel-by-pixel cloud probability; however, other pixel types such as cloud shadow, coastline, or land can also be masked. IdePix is recognized to reliably mask cloudy pixels or pixels of mixed surface types and was applied to all atmospheric correction processors to enhance flagging results, similar to other atmospheric correction evaluations [58,59]. The WQSF_lsb data contain quality flags that are relevant to both the BAC and C2RCC processing; therefore, relevant WQSF_lsb quality flags were applied to each, in line with EUMETSAT match-up protocol guidance. In addition to IdePix, all the recommended quality flags were used for POLYMER and L2gen.

2.5. Match-Up Analysis

Satellite-derived R_rs measurements were co-located in time and space in comparison with in situ measurements following EUMETSAT recommendations [60]. Satellite scenes corresponding to the same date and within 3 hours of an in situ measurement were obtained. The nearest satellite pixel located ≤250 m from an in situ measurement location was identified and R_rs was extracted from a 3 × 3 full resolution pixel window (900 m²) centered on the nearest pixel location. R_rs values from these matchups were included if 50% + 1 pixels in the 3 × 3 pixel window (5 out of 9) were valid (not flagged). Pixel outliers were identified as mean satellite R_rs measurements greater or less than 1.5 times the standard deviation of the set of valid pixels inside the pixel window and were removed from the analysis. To ensure homogeneity, matchups were also discarded if the coefficient of variation at 560 nm was greater than 20%. Median and standard deviation R_rs values were extracted and compared to in situ values using various matchup statistics following EUMETSAT recommendations [60]. Statistical metrics included linear regression statistics (slope, intercept, coefficient of determination (r²)) as well as root mean square error (RMSE), mean deviation (MD) and mean percentage deviation (MPD) to investigate prediction errors, bias and dispersion, respectively.

$$\mathrm{MD}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{Rrs}}_{\mathrm{Sat}}-{\mathrm{Rrs}}_{\mathrm{in} \mathrm{situ}})}{\mathrm{n}}$$

(3)

$$\mathrm{MPD}=\frac{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}100({\mathrm{Rrs}}_{\mathrm{Sat}}-{\mathrm{Rrs}}_{\mathrm{in} \mathrm{situ}})}{{\mathrm{Rrs}}_{\mathrm{in} \mathrm{situ}}}}{\mathrm{n}}$$

(4)

Performance metrics between in situ and satellite-derived R_rs for each atmospheric correction algorithm were also analyzed using Taylor Diagrams [61]. These polar coordinate diagrams provide a compact summary of how three complementary model performance statistics vary simultaneously [62]. The correlation coefficient, standard deviation, and RMSE are related to one another through the Law of Cosines and can be plotted together on one 2D graph [61]. The magnitude of the variability is indicated by computed standard deviation and illustrated as the radial distance from the origin of the plot with points farthest from the origin containing higher variability. Correlation coefficients are shown on the arc of the coordinate plot with points closest to the x-axis containing the highest correlation. RMSE is indicated by the concentric dashed lines emanating from the origin where points farthest from the observed value have high RMSE.

2.6. OLCI Composites

To investigate occurrence and spatial variability of quality flags, composites of clear (unflagged) pixels were created for each atmospheric correction. The number of unflagged pixels in all calendar year 2020 data were summed and mapped to illustrate the spatial variability of ‘clear day’ spatial coverage across Chesapeake Bay. Composites were produced by merging 2020 OLCI scenes processed by each atmospheric correction algorithm (n = 336, 307, 353, 355 for BAC, C2RCC, POLYMER, L2gen, respectively). Data were re-projected to the World Geodetic System 1984 coordinate system (SNAP GPT ‘Reproject’) and the quality flags used in the matchup analysis were applied (Figure 2, Table S3). The ‘xarray’ Python module was used to create climatology datasets. Seasonal R_rs climatologies for individual wavelengths were produced by averaging pixel values across each season (Dec–Feb, Mar–May, Jun–Aug, Sep–Nov).

3. Results

3.1. In Situ and Satellite-Derived Radiometry

In general, hyperspectral in situ R_rs spectra were representative of a eutrophic system (Figure 3) [63,64,65]. Due to absorption of CDOM and chlorophyll a in lower wavelengths, R_rs is low in the blue region (400–500 nm) and higher in the green region (500–600 nm). The red region (600–700 nm) decreases and contains a second minimum around 675 nm corresponding to the red chlorophyll a absorption maximum [64]. The rise in R_rs between 690 and 715 nm is the result of high backscattering and a minimum in absorption by all optically active components besides pure water [64]. This secondary peak can also be attributed to chlorophyll fluorescence emission [65]. Some R_rs spectra are high in the NIR (700–800 nm) due to highly scattering particulate matter [65].

3.2. Matchup Analysis

A time difference of ≤3 h with flagging applied resulted in 27 matchups for BAC, 41 for C2RCC, 48 for POLYMER, and 46 for L2gen. Comparisons of the matchup results to in situ R_rs are presented as mean spectral comparisons, scatterplots, and Taylor Diagrams in Figure 4, Figure 5 and Figure 6 respectively. Averaged across the matchup dataset, C2RCC and L2gen data are most similar to the in situ R_rs spectra across all OLCI wavebands. Conversely, BAC and POLYMER data overestimate in situ R_rs, particularly in the green and red spectral region (Figure 4). All atmospherically corrected R_rs tend to agree towards the longer wavelengths (681–779 nm) with lower standard errors (Figure 4).

Figure 5 and

Table 1. Statistical matchup results for satellite-derived R_rs processed by each atmospheric correction algorithm against in situ R_rs.

BAC (n = 27)							C2RCC (n = 41)
λ (nm)	Slope	Int	r2	RMSE (sr−1)	MD (sr−1)	MPD	Slope	Int	r2	RMSE (sr−1)	MD (sr−1)	MPD
400	−7.97	0.009	0.07	0.009	0.002	1470	0.03	0.002	0	0.001	0.001	303
412	−5.73	0.008	0.05	0.008	0.001	1130	0.02	0.002	0	0.001	0.001	236
443	−2.8	0.009	0.03	0.008	0.002	548	0.01	0.002	0	0.001	0.001	119
490	0.47	0.009	0	0.010	0.006	397	0	0.003	0	0.002	0.001	64
510	1.85	0.008	0.04	0.013	0.008	365	0.11	0.003	0.01	0.002	0.001	47
560	4.01	0.001	0.25	0.018	0.014	333	0.3	0.004	0.08	0.002	0	24
620	5.89	−0.006	0.71	0.015	0.009	311	0.28	0.003	0.16	0.001	0	29
665	5.43	−0.003	0.82	0.012	0.006	230	0.33	0.002	0.32	0.001	0	34
681	4.93	−0.003	0.78	0.011	0.006	298	0.32	0.002	0.3	0.001	0	30
709	5.98	−0.002	0.72	0.013	0.005	396	0.48	0.002	0.47	0.001	0.001	51
754	8.22	0.001	0.5	0.007	0.002	1043	0.94	0.001	0.72	0.001	0	249
779	4.52	0.001	0.48	0.005	0.002	774	0.95	0.001	0.7	0.001	0.001	378
Mean	2.07	0.003	0.37	0.011	0.010	608	0.31	0.002	0.23	0.001	0.001	130
POLYMER (n = 48)							L2gen (n = 46)
λ (nm)	Slope	Int	r2	RMSE (sr−1)	MD (sr−1)	MPD	Slope	Int	r2	RMSE (sr−1)	MD (sr−1)	MPD
400	−36.28	0.054	0.04	0.060	0.01	9809	−4.26	0.008	0.11	0.006	0.005	1514
412	−30.19	0.041	0.05	0.046	0.008	4483	−3.14	0.007	0.08	0.005	0.004	790
443	−8.32	0.024	0.05	0.024	0.006	1403	−0.34	0.006	0.01	0.005	0.004	394
490	1.13	0.008	0.04	0.010	0.008	395	0.5	0.004	0.04	0.004	0.003	154
510	2.08	0.005	0.43	0.009	0.008	277	0.73	0.003	0.13	0.004	0.002	103
560	2.54	0.003	0.6	0.013	0.012	233	0.98	0.002	0.36	0.004	0.003	60
620	2.8	−0.001	0.53	0.009	0.006	192	1.25	0	0.55	0.003	0.002	55
665	2.85	−0.001	0.49	0.007	0.004	205	1.21	0.001	0.55	0.002	0.001	59
681	2.71	−0.001	0.47	0.007	0.005	194	1.08	0.001	0.51	0.002	0.001	51
709	2.98	−0.002	0.49	0.006	0.003	196	1.33	0	0.61	0.002	0.001	667
754	4.12	0.001	0.29	0.004	0.002	1122	1.84	0.001	0.35	0.002	0.001	719
779	3.65	0.001	0.35	0.004	0.001	1059	1.54	0.001	0.39	0.002	0.001	757
Mean	−4.16	0.011	0.32	0.017	0.006	1631	0.23	0.003	0.31	0.003	0.002	444

demonstrate that the accuracy of atmospherically corrected R_rs is wavelength-dependent as the correlation coefficient decreases with wavelength across all AC algorithms (Figure 5, Table 1). Regression slopes are low or negative and correlations of determination are weak in the shorter wavelengths (400–443 nm, slope = −36.28–0.01, r² = 0.00–0.08) and increase in the longer wavelengths (665–779 nm, slope = 0.33–8.22, r² = 0.28–0.78). C2RCC and L2gen are the only processors that do not result in negative regression slopes. BAC and C2RCC contained the highest correlations of determination (BAC at 665 nm, r² = 0.82, C2RCC at 754 nm, r² = 0.72). Mean RMSE and mean deviation was highest in POLYMER and BAC across all wavelengths (POLYMR mean RMSE = 0.017, POLYMER mean MD = 0.006, BAC mean RMSE = 0.011, BAC mean MD = 0.005). C2RCC had the lowest mean RMSE and mean deviation (C2RCC mean RMSE = 0.001, mean MD = 0.001). Mean percentage deviation was overall higher for POLYMER and BAC across all wavelengths (POLYMER mean MPD = 1630, BAC mean MPD = 608) and lowest for C2RCC (mean MPD = 131).

The location of points in the Taylor Diagrams in Figure 6 demonstrate how C2RCC had the lowest standard deviation at all wavelengths (points are close to the origin), contained the highest correlation across all wavelengths (points are close to x-axis), and had the smallest unbiased RMSE (points are close to in situ value). Taylor plots portray L2gen performing second best with POLYMER and BAC performing the poorest across all bands.

3.3. Pixel Flagging

The spatial variability of pixel flagging for each atmospheric correction processor demonstrates that more pixels are flagged (reduction in clear day pixels) in the Upper Bay and along the coastline (Figure 7). Across all 2020 OLCI scenes of only Chesapeake Bay water pixels, IdePix flagged 46% of pixels as clouds and 59% as associated cloud pixels (i.e., cloud margin, cloud shadow, etc.), 39% as bright pixels, and 4% as land although land was previously masked out. The BAC quality flags (WQSF_lsb) flagged 57% of water pixels as cloud-associated flags, 8% as high glint, and 3% as invalid or ‘AC fail’. The BAC science flags (WQSF_msb) flagged 2–11% of water pixels as negative water-leaving radiance for wavebands 412–665 nm, with more flagging at the lower wavelengths. POLYMER quality flags flagged 49% of water pixels as clouds, 4% as invalid and 1% as thick aerosol. L2gen flags flagged 61% as clouds, 9% as atmospheric correction warning, 6% as high glint, and 14% as land although the land was previously masked out. The default land mask for L2gen is approximately 500 m which leads to masking out water pixels in narrow tributaries and could be rectified by passing a null landmask or defining a custom mask (Sean Bailey, personal communication).

Clear day coverage is lowest for the L2gen processing workflow and highest for the POLYMER processing workflow (Figure 7). BAC flagging produced a mean of 52 clear days with higher coverage in the Middle and Lower Bay and less in the Upper Bay (Figure 7). C2RCC flagging produced an average of 58 clear days and is distributed across the entire bay. POLYMER flagging produced an average of 85 clear days which, like BAC, has higher coverage in the Middle and Lower Bay. L2gen flagging produced an average of 47 clear days which is more prominent in the Lower Bay (Figure 7). More pixels are flagged and removed in the Upper Bay and tributaries likely due to increased challenges in atmospheric correction and more optically complex waters. Pixels adjacent to the Bay shoreline are also increasingly flagged, likely due to the adjacency effect where nearby bright pixels are scattered into the field of view of neighboring water pixels.

3.4. R_rs Composites

R_rs composites of the 665 nm OLCI waveband demonstrate differences in R_rs across time and space for each atmospheric correction algorithm (Figure 8). In general, BAC contains the highest R_rs followed by POLYMER and C2RCC while L2gen produced lower overall R_rs (Figure 8). In all atmospheric correction algorithms, R_rs at 665 nm tends to be high in the Upper Bay and tributaries, especially during winter and spring owing to elevated suspended matter concentrations.

4. Discussion

This study is the first attempt to evaluate atmospheric correction performance applied to Sentinel-3A and 3B OLCI scenes of Chesapeake Bay waters. Overall, C2RCC performed best when compared to in situ data collected in the middle and upper Bay. Results have broad implications for improving satellite water quality retrievals in Chesapeake Bay, however, several aspects and caveats of this study merit additional discussion.

4.1. Performance of Atmospheric Correction Processors in Chesapeake Bay

In general, the accuracy of R_rs retrievals from the atmospheric correction processors was wavelength-dependent, with performance increasing at longer wavelengths (Table 1, Figure 4 and Figure 5). The spectral shapes of R_rs averaged across all valid matchup locations were similar to in situ spectra and typical reflectance spectra collected in coastal waters [64]. However, the magnitude of mean R_rs spectra varied across atmospheric correction algorithms (Figure 4). The range of R_rs in the green waveband (560 nm) demonstrates high variability in scattering by inorganic and organic particles (e.g., suspended sediment and phytoplankton). All algorithms produced a distinct secondary peak at 681 nm associated with the chlorophyll absorption maximum and potentially fluorescence, except for C2RCC (Figure 4). Vanhellemont and Ruddick (2021) hypothesize that the C2RCC neural net is limited by its training dataset, and in high reflectance waters, resorts to using spectra containing an absorption dip at 681 nm. Some in situ locations may have been in areas of high reflectance, prompting the neural net to use spectra that, when averaged, resulted in flattening at 681 nm. Therefore, water quality products derived using this wavelength could potentially be erroneous. The BAC and L2gen algorithms that use a two-band iterative approach likely had higher errors in the blue wavelengths associated with extrapolating the aerosol contribution from a NIR band to the visible [7,66]. Similar performance was reported in previous studies using L2gen on Chesapeake Bay scenes collected by NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS) [22,67]. These authors concluded that low signal-to-noise on NIR and SWIR bands propagates into the blue bands, leading to an overestimation of aerosols and negative or low R_rs retrievals [22,67]. The POLYMER algorithm, which uses a spectral optimization approach, does not use a specific aerosol model like BAC and L2gen, but instead, uses the whole spectral range from the blue to NIR to decouple the atmospheric and surface components of the TOA signal from the water leaving reflectance [57]. This spectral matching approach is based on chlorophyll a concentration and backscattering and does not contain a parameter directly related to non-varying CDOM absorption. Therefore, POLYMER likely overestimated R_rs in the blue bands due to the exclusion of highly absorbing CDOM found in Chesapeake Bay waters.

All atmospheric correction algorithms had better performance in the longer wavelengths (560–865 nm) with C2RCC containing the lowest RMSE, MD, and MPD (Table 1, Figure 4 and Figure 5). POLYMER retrieved a few unrealistic negative values in bands 665–865 nm, likely caused by the spectral matching approach. The applicability of the assumptions made in the in-water model used in POLYMER may not be suitable for eutrophic waters, where different IOP models produce contrasting outputs with varying phytoplankton concentrations (Lain et al., 2014). R_rs at 665–709 nm are wavelengths that are frequently used to derive water turbidity or total suspended matter [5], and chlorophyll a in coastal waters [68,69]; therefore, water quality products derived from these bands are likely to be more accurate than band-ratio algorithms incorporating blue wavelengths.

In situ measurements in this study cover a dynamic range for the upper and middle Bay as well as locations in the Choptank River (Figure 3). These measurements, conducted throughout the year, are likely to be representative of waters spanning the entire Chesapeake Bay; however, the authors acknowledge the limitation in the number of samples and the need for further radiometric validation with a larger and more extensive dataset. Future studies should consider incorporating in situ radiometric data in the Lower Bay and measurements collected from the recently installed radiometer included in the Ocean Color Component of the NASA Aerosol Robotic Network (AERONET-OC) located at 39.13N, −76.32W [70].

4.2. Comparison to Similar Studies

The results of this study can be compared to other OLCI matchup analyses in coastal waters found in the literature (

Table 2. Summary of existing OLCI atmospheric correction evaluation literature. Table information includes the atmospheric correction algorithms used, study site, in situ R_rs matchup sample size, and performance metric(s) for four wavelengths, if available. BAC: EUMETSAT’s Baseline Atmospheric correction algorithm, C2RCC: EUMETSAT’s Case 2 Regional Coast Color, C2R-CCAltNets: alternative neural nets of C2RCC, POLYMER: POLYnomial based algorithm applied to MERIS, L2gen: NASA’s standard atmospheric correction algorithm, ACOLITE/DSF: ACOLITE atmospheric correction with dark spectrum fitting (DSF). MD: Mean deviation, MPD: Mean Percentage Deviation.

		443			560			665			754
	n	r2	MD	MPD	r2	MD	MPD	r2	MD	MPD	r2	MD	MPD
This study
BAC	27	0.03	0.002	548	0.25	0.014	333	0.82	0.006	230	0.50	0.002	1043
C2RCC	41	0	0.001	119	0.08	0	24	0.32	0	34	0.72	0	249
POLYMER	48	0.05	0.006	1403	0.60	0.012	233	0.49	0.004	205	0.29	0.002	1122
L2gen	46	0.01	0.004	394	0.36	0.003	60	0.55	0.001	59	0.35	0.001	719
Mograne et al., 2019
BAC	18	0.38 *	--	−40 *	0.95 *	--	−10 *	0.95 *	--	−40 *	--	--	--
C2RCC	16	0.63 *	--	−5 *	0.93 *	--	0 *	0.95 *	--	−20 *	--	--	--
C2R-CCAltNets	19	0.84 *	--	15 *	0.94 *	--	10 *	0.90 *	--	5 *	--	--	--
POLYMER	35	0.54 *	--	−18 *	0.86 *	--	−12 *	0.93 *	--	−18 *	--	--	--
L2gen	17	0.50 *	--	22 *	0.98 *	--	−4 *	0.98 *	--	−12 *	--	--	--
Alikas et al., 2020
BAC	22	0.028	--	−192	0.59	--	−2	--	--	12	--	--	58
C2RCC	22	0.29	--	64	0.6	--	16	--	--	39	--	--	122
C2RCC ALTNN	22	0.003	--	50	0.37	--	7	--	--	41	--	--	72
POLYMER	22	0.82	--	6	0.96	--	14	--	--	20	--	--	55
Vanhellemont & Ruddick, 2021
BAC	46	0.61	−0.013	--	0.86	−0.008	--	0.90	0.000	--	0.73	−0.001	--
POLYMER	27	--	−0.025 *	--	--	−0.033	--	--	−0.025	--	--	−0.015 *	--
C2RCC	27	--	−0.005 *	--	--	−0.018	--	--	−0.017	--	--	−0.008 *	--
L2gen	27	--	−0.005 *	--	--	−0.012	--	--	−0.007	--	--	0.000 *	--
ACOLITE/DSF	46	0.37	−0.005	--	0.86	−0.004	--	0.90	−0.003	--	0.67	0.001	--
Giannini et al., 2021
C2RCC	24	0.74	--	−33.6	0.83	--	−28.8	0.76	--	−36.3	--	--	--
C2RCC altNN	24	0.36	--	−36.9	0.72	--	−18.2	0.66	--	−12	--	--	--
POLYMER	24	0.73	--	−41.5	0.83	--	−23.6	0.71	--	−28.4	--	--	--

* estimated values from spectra plots of derived statistics.

) [23,34,35,40,71]. However, statistical reporting across these studies is not consistent, which makes it challenging to compare atmospheric correction performance across different study sites. EUMETSAT provides guidelines for a common matchup approach to achieve a consistent validation baseline to compare across different studies [60], and were used in this study. Adopting these matchup statistics will provide comparable values and allow for better investigation of the performance of atmospheric correction processors across various regional coastal water bodies.

Mograne et al. (2019) conducted an R_rs matchup analysis focused on the coastal waters of the Eastern English Channel and French Guiana (n = 18, 19, 35, 17 for BAC, C2RCC, POLYMER and L2gen, respectively). For low wavelengths (400–510 nm), BAC was always negatively biased and contained the highest RMSE while C2RCC was always positively biased and contained the lowest RMSE. For longer wavelengths (560–754 nm), C2RCC also contained the lowest RMSE while POLYMER contained high RMSE. The authors recommended the use of POLYMER or C2RCC to process OLCI images over their region of interest. An R_rs matchup analysis focusing on the coastal waters of British Columbia (n = 24) demonstrated that C2RCC performed best since it allowed for large variability in water-leaving radiance for regions impacted by the presence of CDOM [40]. POLYMER showed slightly inferior radiometric performance; however, it was the most accurate when comparing Level-2 operational products (e.g., chlorophyll a and TSM) products to in situ measurements. The authors hypothesize this was due to the preservation of the spectral shape in R_rs compared to C2RCC [40]. Hieronymi et al. (2017) developed the OLCI Neural Network Swarm (ONNS) which consists of several specialized neural networks to derive ocean color remote sensing products in a wide range of optical properties. This study applied ONNS to C2RCC-corrected OLCI data and compared derived water quality concentrations to in situ measurements (n = 43). The results of this study show relatively good agreement to in situ data, particularly for inorganic suspended matter and CDOM (r = 0.86 and 0.73, respectively).

Results of this study also contrast with similar matchup studies focused in coastal waters [23,71]. An OLCI R_rs matchup in the Belgian coastal zone (n = 46) demonstrated poor performance in POLYMER and C2RCC while BAC performed relatively well [23]. Vanhellemont and Ruddick (2021) hypothesize that POLYMER and C2RCC are constrained by the bio-optical models and are unable to represent higher reflectance in turbid water, resulting in an underestimation of R_rs at red wavelengths. These results are not consistent with the present study; however, in situ locations were also not located in highly turbid waters. Future work should test the accuracy of C2RCC and POLYMER in more turbid conditions in Chesapeake Bay. Additionally, a matchup analysis focused on Estonian inland waters and Baltic Sea coastal waters (n = 22) demonstrated that POLYMER was the most suitable and provided the most amount of unflagged data and C2RCC was not suitable for highly absorbing waters [71]

The contrasting radiometric performance across various regional coastal water bodies (Table 2) demonstrates the challenge of developing a robust atmospheric correction algorithm for global coastal waters. Extensive validation with regional in situ data is required to ensure water-leaving retrievals are reliable for input into inversion models to derive accurate water quality components used for research and management. Comparison of performance in different regions can also yield vital information to refine assumptions made in different algorithms, with regard to the inherent optical properties of different regimes.

It is also important to note that there are other atmospheric correction algorithms that were not evaluated in this study that could potentially perform better than the processors applied in this study. For example, Vanhellemont and Ruddick (2021) apply ACOLITE, an atmospheric correction processor developed by the Royal Belgian Institute of Natural Sciences for OLCI imagery of coastal Belgian waters. The adaptation of the ‘Dark Spectrum Fitting’ approach was shown to provide the best performance in the visible bands in comparison to other atmospheric correction processors including C2RCC and POLYMER. Xue et al. (2019) applied the 6SV atmospheric correction (the vector version of the Second Simulation of the Satellite Signal in the Solar Spectrum correction scheme) to OLCI imagery of Chinese inland waters which demonstrated greater performance than C2RCC and POLYMER [72]. Schroeder et al. (2022) developed a new coastal atmospheric correction algorithm for OLCI that consists of an ensemble of ANNs and obtained lower mean absolute percentage error and band-averaged bias when compared to C2RCC [73]. This algorithm also provides pixel-based estimation of model inversion uncertainty and sensor noise propagation which is critical to determine the quality with which ocean color products can be estimated from [73].

4.3. Matchup Uncertainties

Although radiometric in situ measurements are typically used to assess satellite-derived measurements, they can contain variability and uncertainties due to instrument calibration, instrument or ship shadowing, environmental factors, and methods to address surface reflection factors [74,75,76,77,78]. Quantifying uncertainties from in situ measurements is an important step to obtain a reliable fiducial reference measurement and has been performed in previous studies [77,79,80]. An intercomparison of in situ radiometers and computing an uncertainty budget were not the main objectives for this study and therefore were not included. However, quality control measures were taken to eliminate outliers in in situ spectra. The agreement between R_rs derived from the two in situ measurements techniques (Figure 3) provides assurance that measurements were collected with methods that minimize uncertainties. Nonetheless, additional measures to reduce and compute uncertainty would potentially provide higher precision and accuracy of in situ R_rs.

Radiometric data can also be influenced by viewing and illumination geometries (i.e., solar zenith angle, viewing zenith angle, relative azimuth angle between the Sun and the sensor). This geometry dependence, known as BRDF, aims to characterize how water-leaving radiance varies depending on viewing and illumination angles [81]. A BRDF correction converts measured water-leaving radiance at any viewing geometry to a conceptual “exact normalized” water-leaving radiance measurement at a nadir direction using values obtained through a look-up table derived by Morel et al. (2002). This correction is based on atmospheric and bio-optical models developed for open ocean waters and requires a priori knowledge of chlorophyll a concentration. This BRDF correction has been validated against open ocean measurements to show good agreement [82]; however, it can fail in regions with different atmospheric conditions and confounding inherent optical properties, such as coastal waters. Methods to improve the BRDF correction in coastal waters is an active area of research [83,84]; however, due to the current uncertainty, this study did not apply BRDF corrections to in situ or satellite-derived R_rs. A well performing BRDF correction for coastal waters could potentially improve matchup uncertainties.

Ocean color sensors are pre-calibrated before launch; however, residual calibration errors develop throughout time, requiring a calibration process to ensure satellite retrievals meet water-leaving radiance requirements (i.e., L_w within 5% uncertainties) [8,85]. This process includes measuring the difference between satellite-retrieved TOA radiance to radiance estimated by back propagating high quality in situ water leaving radiance to the TOA [8,85]. Therefore, calibration gain factors are specific to a particular sensor and atmospheric correction algorithm. Gain factors, often termed system vicarious calibration (SVC) gains, are recomputed by conducting a series of matchup comparisons of satellite and in situ data and are applied in routine atmospheric correction processing chains. Recommended SVC gains and sources for each atmospheric correction algorithm can be found in Table S2. Publicly available Level 2 OLCI data (BAC) is processed with the most current SVC gains [60]. In this study, BAC data were downloaded at different times throughout 2019 to 2021 and contained varying SVC gains which could potentially impact matchup results. It is recommended to download data once to ensure SVC gains are consistent throughout. C2RCC was processed with SVC gains listed in the 2021 EUMETSAT Level 2 report for baseline collection [60]. POLYMER was processed with the default SVC gains hardcoded in POLYMER v4.12 which were updated in 2019 (POLYMER changelog, 2021-12-17. Accessed online: 13 March 2022, https://www.hygeos.com/static/polymer/CHANGELOG.TXT). OLCI-A gains were applied by default to OLCI-B, which could potentially impact the accuracy of OLCI-B retrievals since gains should be sensor specific. L2gen was processed with SVC gains derived following procedures of Franz et al. (2007) and Werdell et al. (2007) with the use of a model based on a climatology of chlorophyll at the Marine Optical Buoy (MOBY) site (NASA Ocean Color Forum Post. Accessed online: 13 March 2022, https://forum.earthdata.nasa.gov/viewtopic.php?f=7&t=1273&sid=5e1209c2e4f00bec7683116cc017d6fb). Operational L2gen SVC gains are derived using open ocean calibration targets and typical maritime atmospheric conditions which may not be representative of a coastal water body [85]. Therefore, regionally specific vicarious calibration may be necessary in coastal and inland waters and would potentially improve the matchup analysis in this study.

4.4. Pixel Flagging

Many solutions to atmospheric correction challenges have been operationally implemented in processing schemes; however, there are some issues where satisfactory solutions have not been developed, causing atmospheric correction to fail, or flag pixels as suspect quality such as absorbing aerosols, whitecaps, and the adjacency effect [66]. POLYMER produced more clear day spatial coverage than the other atmospheric correction algorithms, likely due to POLYMER being capable of retrieving data under sun-glint conditions. L2gen resulted in the highest amount of flagged pixels, likely due to the restrictive filtering of stray light and sun glint, which has been shown to be the primary quality flags affecting L2gen data quantity [86] and the 500 m default land mask.

More pixels were flagged in the Upper Bay for all atmospheric correction algorithms likely due to high reflectance from scattering suspended mineral particles associated with terrestrial runoff or potential errors associated with nonuniform aerosols. The Chesapeake Bay watershed contains aerosol influences from nearby industrial cities and surrounding agriculture; however, studies have shown that aerosol characteristics are consistent and related to relative humidity [9]. Future work should focus on correcting for the land adjacency effect as demonstrated in [87] and relaxing filtering criteria to observe differences in data quality and quantity.

4.5. Management Implications

Satellite remote sensing provides water quality observations at spatial and temporal scales that exceed current field monitoring techniques and can be used to enhance the understanding of coastal processes and monitoring of marine resources. Water quality concentrations such as chlorophyll a provide indicators of nutrient enrichment and can be used to predict eutrophic conditions [88]. Anomalous high levels of chlorophyll a concentration can also be used to identify harmful algal blooms which are only expected to increase in abundance and severity with climate change [89,90]. Measuring satellite-derived TSM concentrations can provide insight into nutrient loads, salinity, and light levels [6].

Satellite remote sensing has been used to study water quality trends in Chesapeake Bay with various ocean color sensors [6,22,43,63,67,90,91,92,93]. Several of these studies demonstrate erroneous water-leaving radiance due to the narrowing shape of the estuary, high optical complexity, and increased urban aerosols that confound the atmospheric correction process [22,67]. These uncertainties have prompted the use of utilizing R_rhos or R_rc instead of R_rs to derive water quality parameters [90,91,92]. R_rhos is derived by a partial atmospheric correction that corrects TOA radiance for gaseous absorption (mainly ozone) and Rayleigh scattering effects, resulting in Rayleigh-corrected reflectance [91]. However, since this method does not account for aerosol contributions, it is likely that R_rhos has aerosol-associated uncertainties. Other studies have used custom workflows for MODIS Aqua data based on a threshold value at R_rs at 645 nm that switches to SWIR radiances, unavailable with OLCI, in highly turbid waters [93,94].

OLCI currently contains the most spectral bands of any ocean color satellite sensor and provides nearly daily coverage at 300 m spatial resolution; characteristics which make it an important sensor for near real time retrieval of ocean color products. A reliable and well validated atmospheric correction algorithm for Chesapeake Bay waters is necessary to utilize the benefits OLCI contains for reliable water quality monitoring. The C2RCC atmospheric correction is the best performing algorithm for Chesapeake Bay waters in this study and agencies can consider incorporating R_rs and water quality products produced by C2RCC into operational data portals (i.e., NOAA CoastWatch, NASA OceanColor Web) to help end-users to select the most appropriate and accurate data for their needs.

4.6. Future Approaches

Atmospheric correction is done on a pixel by pixel basis, which can be computationally intensive due to the large number of pixels contained in multiple satellite images captured by numerous satellite sensors per day. Therefore, algorithms need to be very computationally fast while maintaining accuracy to produce data in near real-time. Algorithms will likely advance and change in accordance with computational power, new ocean color sensors, and innovative ideas, much like they have over the past several decades. There have been several new proposed atmospheric correction algorithms that have the potential to enhance the upcoming Plankton, Aerosol, Cloud, and ocean Ecosystem (PACE) mission, expected to launch in 2024. PACE will host the first hyperspectral ocean color sensor (the Ocean Color Instrument, OCI) which will measure unprecedented TOA hyperspectral radiance from 340 to 890 nm at 5-nm spectral resolution [95]. It will also contain seven discrete channels in the SWIR with high radiometric performance to enable atmospheric correction over open ocean and coastal waters [24]. Ibrahim et al. (2019) describes a multi-band atmospheric correction (MBAC) which utilizes six bands from the NIR and SWIR region as compared to the heritage two-band ratio aerosol correction algorithm. This multiband approach reduces uncertainty caused by sensor noise using adaptive spectral weights to decrease the contamination of a highly scattering signal [24]. Preliminary analyses show significant improvement in radiometric retrievals when compared to NASA’s two-step atmospheric correction procedure [24].

Quantifying polarizing properties of reflected sunlight also has the potential to improve atmospheric correction. Solar radiation is unpolarized when it enters the Earth’s atmosphere; as it interacts with aerosols and molecules and refracted at the atmosphere-ocean interface, it can become partially polarized. Studies have shown that polarization of incident light on the sea surface can lead to significant variations in water-leaving radiances [96]. Accounting for polarization can enhance the contribution of water-leaving reflectance, improve retrievals impacted by sun glint, and better characterize aerosol properties [66]. It has also been demonstrated that measuring polarization at multiple viewing angles can significantly improve aerosol radiation properties which can assist in determining appropriate aerosol types or models in atmospheric correction schemes [66]. Multi-angle polarimetry is an emerging field of active research and will be included in the upcoming PACE mission to improve atmospheric correction [95].

5. Conclusions

Coastal and estuarine ecosystems are dynamic environments where tidal forces, sediment resuspension, and riverine and terrigenous inputs continuously work to alter the relative concentrations of optically active constituents. The OLCI sensor can resolve optical data in coastal water bodies at higher spatial (300 m), spectral (21 bands), and temporal (2-day) resolutions than any other operational ocean color satellite sensor. Since water-leaving radiance is at most 10% of the total radiance OLCI collects, accurate atmospheric correction is essential. The challenge of atmospheric correction is exacerbated in coastal waters due to the optically complex water-leaving reflectance signal and the presence of strongly absorbing aerosols with differing composition and vertical distributions. This study conducted a radiometric evaluation of four different atmospheric correction algorithms compared to a regional in situ dataset from Chesapeake Bay waters. A matchup analysis demonstrated that C2RCC had the best performance, particularly in the longer wavelengths. The results from this study both agreed and contrasted with similar studies in regional coastal water bodies, depicting the challenge of developing a robust atmospheric correction algorithm that works in all global coastal waters. Future research will include comparing biophysical products derived from C2RCC processed data to an extensive dataset of in situ water quality concentrations collected in the Chesapeake Bay to evaluate performance and analyze water quality variability. This study provides a framework with associated uncertainties and recommendations to utilize OLCI ocean color data to monitor the biogeochemical dynamics in the Chesapeake Bay.