Multiband Atmospheric Correction Algorithm for Ocean Color Retrievals

NASA’s current Atmospheric Correction (AC) algorithm for ocean color utilizes two bands and their ratio in the Near Infrared (NIR) to estimate aerosol reflectance and aerosol type. The algorithm then extrapolates the spectral dependence of aerosol reflectance to the visible wavelengths based on modeled spectral dependence of the identified aerosol type. Future advanced ocean color sensors, such as the Ocean Color Instrument (OCI) that will be carried on the Plankton, Aerosol, Cloud, and ocean Ecosystem (PACE) satellite, will be capable of measuring the hyperspectral radiance from 340 to 890 nm at 5-nm spectral resolution and at 7 discrete Short-wave Infrared (SWIR) channels: 940, 1038, 1250, 1378, 1615, 2130, and 2260 nm. To optimally employ this unprecedented instrument capability, we propose an improved AC algorithm that utilizes all atmospheric-window channels in the NIR to SWIR spectral range to reduce the uncertainty in the AC process. A theoretical uncertainty analysis of this, namely Multi-Band AC (MBAC), indicates that the algorithm can reduce the uncertainty in remote sensing reflectance (Rrs) retrievals of the ocean caused by sensor random noise. Furthermore, in optically complex waters, where the NIR signal is affected by contributions from highly-reflective turbid waters, the MBAC algorithm can be adaptively weighted to the strongly-absorbing SWIR channels to enable improved ocean color retrievals in coastal waters. We provide here a description of the algorithm and demonstrate the improved performance in ocean color retrievals, relative to the current NASA standard AC algorithm, through comparison with field measurements and assessment of propagated uncertainties in applying the MBAC algorithm to MODIS and simulated PACE OCI data.


INTRODUCTION
Ocean color retrieval algorithms require an atmospheric correction (AC) process to separate the radiometric contribution of the atmosphere from the ocean, given the radiance measured at the top of atmosphere (TOA). The National Aeronautics and Space Administration (NASA) standard approach to the AC is a two-step process: the atmosphere and surface contributions are first removed using minimal assumptions about the water optical properties (Gordon and Wang, 1994;Franz et al., 2007;Ahmad et al., 2010;Bailey et al., 2010;Mobley et al., 2016), and the resulting spectral water-leaving radiances are then used to infer water column optical properties and constituent concentrations (e.g., Gordon et al., 1988;O'Reilly et al., 1998;Morel and Maritorena, 2001;Hu et al., 2012;Werdell et al., 2013). AC involves finding a solution to a set of deterministic models, enabling the removal of atmospheric path and surface effects from the TOA signal. The primary challenge is determining the contribution of aerosols (a primary source of uncertainty in the AC) to the atmospheric path radiance, which is highly variable and, thus, must be inferred from observation. This approach takes advantage of the strong absorption of water in the Near Infrared to Shortwave Infrared (NIR-SWIR) (longward of 700 nm) to separate the atmospheric and oceanic signals (Gordon and Wang, 1994;Wang et al., 2009;Ahmad et al., 2010). This AC process aims for retrievals that meet ocean color requirements (i.e., 5% uncertainty on water-leaving radiance at 443 nm) over the open ocean and in the absence of strongly absorbing aerosols (i.e., dust, organic, and black carbon) (Gordon and Wang, 1994;Mobley et al., 2016). In coastal and optically complex waters, the water-leaving radiance retrievals often fail in meeting that requirement due to non-negligible water-leaving radiance in the NIR bands used for the AC, with water-leaving radiance often underestimated or even negative in value. Iterative AC techniques have been developed to model the ocean reflectance in the NIR, which mitigate the problem effectively in moderately turbid and productive waters (Siegel et al., 2000;Bailey et al., 2010). Although this iterative approach has been used in the NASA standard processing , it has been shown that such methods can fail in highly turbid waters Goyens et al., 2013). Other methods, based on the SWIR bands, have been shown effective in separating the aerosol signal from the TOA signal over highly turbid waters, due to the increased water absorption at these longer wavelengths Shi, 2007, 2012;Jiang and Wang, 2014;Vanhellemont and Ruddick, 2015;Pahlevan et al., 2017). The SWIR-based methods has limitations, however, for two reasons: (1) the signal-to-noise ratio (SNR) on SWIR bands is low for current ocean color sensors [i.e., MODIS (Moderate Resolution Imaging Spectroradiometer) and VIIRS (Visible Infrared Imaging Radiometer Suite)] and (2) the approach generally requires discrete switching from the NIR AC bands to SWIR bands over coastal waters, thus producing spatial inconsistency or artifacts in the retrievals. For example, Shi and Wang (2007) utilized 1,240 and 2,130 nm for MODIS to perform the AC using the 2-band Gordon and Wang (1994) AC approach. The algorithm switches from NIR to SWIR based on a turbidity index threshold Wang and Shi, 2007;Wang et al., 2009). Werdell et al. (2010) have shown the effectiveness of the 2-band SWIR AC approach, but with increased uncertainties in the ocean color retrievals due to low SWIR SNR of MODIS-Aqua. Still, other studies have shown failure in the NIR-SWIR switching method due to the NIR sensor saturation at moderate to high turbidity (i.e., >35 g m −3 ) where the NIR bands on MODIS-Aqua saturate at TOA radiances >2.8-3.45 mW cm −2 µm −1 sr −1 for 748 nm and >1.9-2.45 mW cm −2 µm −1 sr −1 for 869 nm and improper cloud masking (Aurin et al., 2013). A recent study by Liu et al. (2019) indicated that the turbidity index, responsible for the switching, varies with the aerosol size distribution, optical depth, and observing geometry, thus rendering the effectiveness of the NIR-SWIR switching questionable.
One way for reducing the uncertainty in ocean color retrievals using the SWIR bands is to use more spectral information to reduce sensor random noise impact produced by low SNR. Gao et al. (2000Gao et al. ( , 2007 used a spectral matching technique to fit the measured TOA signal with an aerosol lookup table (LUT) for bands >860 nm and demonstrated improved retrieval of ocean reflectance in turbid and shallow waters from MODIS-Aqua and Airborne Visible/Infrared Imaging Spectrometer (AVIRIS). They also showed a sensitivity study indicating that SWIR bands are sensitive to variations in aerosol size distributions including some fine-mode dominant types, although the sensitivity is weaker than that of the visible-NIR spectral regime. Other methods that utilize more spectral information for the AC, such as machine learning techniques (i.e., neural network) and coupled atmosphere-ocean simultaneous retrievals, have shown promising results (Gordon et al., 1997;Gordon, 1998, 2001;Chomko et al., 2003;Stamnes et al., 2003;Jamet et al., 2005;Brajard et al., 2006;Schroeder et al., 2007;Spurr et al., 2007;Kuchinke et al., 2009;Steinmetz et al., 2011;He et al., 2012;Frouin and Gross-Colzy, 2016;Fan et al., 2017;Gossn et al., 2019). Some of these algorithms were more successful over the open ocean, while others were specifically tuned for coastal cases, and they often require the use of a biooptical ocean model with inherent assumptions about the marine optical properties.
Cloud detection and masking in turbid waters can also be challenging. Extremely turbid waters are masked as clouds when a simple NIR threshold technique is used, as is the case for the standard NASA processing. Wang and Shi (2006) proposed an improved cloud masking in coastal regions using SWIR channels. Other cloud masking methods that utilize more spectral bands have shown improved discrimination between clouds and water in highly turbid conditions (Nordkvist et al., 2009).
NASA is planning to launch the Plankton, Aerosol, Cloud, and ocean Ecosystem (PACE) satellite in 2022. The PACE satellite will host a state-of-the-art hyperspectral ocean color sensor (the Ocean Color Instrument, OCI). OCI will measure the hyperspectral radiance at the TOA from 340 to 890 nm at 5-nm spectral resolution and at seven discrete channels in the SWIR: 940, 1,038, 1,250, 1,378, 1,615, 2,130, and 2,260 nm, of which two are water vapor bands (940 and 1,378 nm) and five are window bands (i.e., no major gas absorption features). In fact, OCI will host the first instrument with a set of SWIR bands primarily designed for the ocean AC. The SWIR bands will have sufficient radiometric performance (i.e., high SNR with low systematic bias) to enable their use for AC over open oceans as well as coastal water conditions (Ibrahim and McKinna, 2018). The analysis throughout the manuscript assumed the following SNRs for OCI (420 at 1,038 nm, 280 at 1,250 nm, 220 at 1,615 nm, and 76 at 2,260 nm, excluding 2,130 nm as cloud band), for the typical ocean scene, which are based on current best estimates since the design is currently subject to slight changes. The AC algorithm proposed here attempts to exploit the unprecedented capabilities of OCI to improve the quality of ocean color retrievals from the PACE mission.

MULTIBAND AC ALGORITHM
The heritage aerosol correction algorithm utilizes two bands either in the NIR or SWIR spectral range by calculating a band ratio of the Rayleigh (+ surface) and gas corrected reflectance, namely, epsilon, as follows: where λ s and λ l correspond to the short and long wavelength either in the NIR or SWIR range (e.g., for MODIS, λ s = 748 nm and λ l = 869 nm) and ρ a (λ) is the aerosol reflectance defined as: F 0 is the extraterrestrial solar irradiance and µ 0 is the cosine of the solar zenith angle. Since the aerosol radiance term is a strong function of the viewing/illumination geometry, the LUT is generated as a function of the solar and sensor zenith angle and relative azimuth angle. NASA's current approach to estimate the aerosol reflectance term is by finding the best match of the observed single scattering epsilon value to the model single scattering epsilon value. Since the LUT is calculated for a set of deterministic aerosol models based on Ahmad et al. (2010), the algorithm interpolates the single scattering aerosol reflectance from the LUT to match the observed reflectance. Accounting for multiple scattering interactions between the aerosol and molecules, is crucial when extrapolating the reflectance to shorter wavelengths. Given the aerosol type estimated from the single scattering epsilon in the longer wavelengths, the multiple scattered aerosol radiance is then estimated through a polynomial relationship derived through vector radiative transfer simulations (Gordon and Wang, 1994), and this multiscattering estimate of aerosol path reflectance is then used in the aerosol correction. The process of deriving single scattering epsilon and then converting it into multiple scattering radiance is unnecessary since currently computationally feasible multiple scattering tables can be easily generated and irreversible since there is no 1-1 relationship between multiple and single scattering when different aerosol types are mixed. Ahmad and Franz (2018) proposed an alternative approach, where the aerosol model and optical properties determination is done in multiple scattering space. It is based on the Ahmad et al. (2010) aerosol models, for which the TOA reflectance due to aerosols can be described by three parameters: relative humidity (RH), Angstrom exponent (α), or the analogous fine-mode fraction (f ), and optical thickness of the aerosols (τ ) in the model atmosphere. Briefly, for any relative humidity suite, the dependence of multiple-scattering epsilon ε (λ) = ρ(λ) ρ(869) on ρ(869) would look like the x-y plot shown in Figure 1. Here, ρ refers to TOA aerosol reflectance and the subscript numbers refer to wavelength bands.
In Figure 1, the aerosol optical depth (τ ) varies along the "horizontal" axis, and the Angstrom exponent (α) and FIGURE 1 | A graph of multiple-scattering epsilon (ε 748 = ρ 748 /ρ 869 ) vs. ρ 869 for solar zenith, θ 0 = 36 • , view zenith, θ = 30 • , relative azimuth angle, ϕ = 120 • , and relative humidity (RH) = 80% from MODIS LUT . extinction coefficients (ξ ) vary along the "vertical" axis. This parametrization of the LUT allows one to readily calculate the aerosol reflectance at any grid point. One advantage of this approach over the current NASA standard algorithm is that it provides a mathematical framework for the aerosol reflectance determination that facilitates the application of standard error propagation techniques.
The multiband AC (MBAC) algorithm is based on this multiscattering epsilon approach by Ahmad and Franz (2018), but it utilizes all the spectral window channels (i.e., not affected by strongly absorbing gases such as water vapor, oxygen, and carbon dioxide) available at which the water absorption is high to separate the aerosol signal (i.e., black-pixel). MBAC determines the best aerosol model by fitting the magnitude and the spectrum of the aerosol reflectance at all window NIR to SWIR bands (e.g., excluding 940-and 1,378-nm channels for OCI). The aerosol reflectance at each wavelength is stored in a LUT as a polynomial function of the geometry (not shown in Equation 3 for simplicity), aerosol optical depth, relative humidity, and fine-mode fraction: where a, b, and c are the fitting coefficients stored in the aerosol LUT. For any spectral band in the NIR-SWIR range at which the water-leaving radiance can be assumed negligible, the atmospheric path radiance can be derived from the observed TOA radiance. Assuming a proper correction for absorbing gases (discussed in the section Correction for Absorbing Gases) and Rayleigh scattering has been done, aerosol reflectance can be determined and the aerosol optical depth can then be estimated by solving the quadratic equation, as For consistency with the system vicarious calibration procedure, which is performed relative to one wavelength in the NIR (e.g., 869 nm for MODIS), the optical depth is estimated at that wavelength.
For an observed aerosol reflectance, the optical depth can be calculated for each aerosol model set (RH, f ) as shown in Figure 2.
The algorithm then calculates the aerosol optical depth and reflectance spectrally for all aerosol models as: where ξ λ, RH, f is the aerosol extinction coefficient for each aerosol model. Then, the spectral aerosol reflectance for each aerosol model set (RH, f ) is recalculated from the polynomial equation.
The MBAC algorithm fits the aerosol reflectance in the NIR and SWIR for all models calculated in the previous step. However, since multiple solutions can exist, especially at low optical depth over the open ocean, we constrain the model selection to consider only those models with the two relative humidities bracketing the observed RH (as obtained from ancillary meteorological data). The maximum likelihood method is used to find the closest solution by minimizing the cost function: where ρ obs (λ) is the observed aerosol reflectance and ρ a λ, RH, f is the aerosol reflectance from the LUT for each FIGURE 2 | The aerosol reflectance as a function of the aerosol optical depth at 869 nm for aerosol models with relative humidity RH = 75% and 80% and for fine-mode fraction f = 30% and 50%. relative humidity and fine-mode fraction. Note that the geometry is omitted from Equation 6 for simplicity; however, it depends on the solar and viewing zenith and relative azimuth angles. DOF is the degree of freedom of the observed signal. Assuming each band in the NIR and SWIR is an independent observations, DOF is the number of bands from the shortest band in the NIR, λ s , to the longest band in the SWIR, λ l . σ 2 (λ) is the squared measurement uncertainty (variance assuming gaussian random noise) of the sensor, as determined from the sensor-specific noise model for a given radiance level.
SW(λ) is a spectral weight (SW) that is set to change the relative weighting in the cost function for the spectral bands in different environmental conditions. For example, setting SW (λ) to 1 in all bands will minimize the cost function to fit the aerosol reflectance to all NIR/SWIR bands. The method will then be skewed to bands with better radiometric performance (i.e., the NIR). And when the weight is set 0 for NIR bands, the cost function will be skewed toward longer wavelengths in the SWIR, where the ocean signal is still negligible. In this implementation, we vary the SW, adaptively, as a function of the number of iterations in the iterative scheme of Bailey et al. (2010), which is currently used in the standard NASA AC algorithm to estimate any non-negligible water-leaving reflectance contribution in the NIR bands. The Bailey algorithm works by modeling the ocean reflectance in the NIR based on a bio-optical model for productive waters. The algorithm is initiated by a first guess of the NIR ocean reflectance [i.e., R rs (NIR) = 0], where the AC is performed, and then the R rs (NIR) is estimated based on a bio-optical model that is propagated to the TOA to perform another AC and this process is repeated until the changes in the retrieved ocean reflectance is <2%. Similarly, the MBAC algorithm performs the AC at every iteration while bio-optical model parameters and the SW are tuned until the convergence criterion of 2% change in ocean reflectance is met. Thus, more iterations indicate more difficulty in fitting the bio-optical model possibly due to a higher optical complexity in turbid waters. The SW is then calculated as , where i is the number of iterations in the current step and i max is the maximum number of iterations set in the Bailey algorithm (operationally i max = 11). Note that the actual number of iterations is reduced by 2 since the algorithm always performs at least two iterations in the AC: the first one with Rayleigh correction only (to estimate glint), and the second one for the aerosol correction.β (λ) is an empirically estimated value that determines the slope of the exponential function for a given sensor by processing various scenes to inspect for image artifacts. We adopted an exponential function for the SW to ensure a quick transition toward the SWIR bands once non-negligible NIR ocean signal is detected since the SWIR bands will have an effect only when the weight on the NIR bands is largely reduced. This is because measurement uncertainty on the NIR is nearly an order of magnitude lower than the SWIR bands (at least for MODIS and VIIRS sensors). Higher value of the slopeβ (λ) > 10 leads to a quick switching to the SWIR bands in coastal waters, including for low/moderate turbidity, while smaller values <5 give more dependence on the NIR AC in the same water conditions. For MODISA,β (λ) = 7 for NIR bands and 0 for SWIR bands, as this provides for a smooth transition from open ocean AC, where we want an equal SW on all bands to a smaller weight on the NIR in more turbid waters. In the open ocean, NIR bands will dominate the cost function; thus, the AC will be similar to the standard approach, while in turbid waters, the AC will depend more on the SWIR bands. The slopeβ (λ) will primarily impact the number of iterations required in the convergence; thus, it is not critically important to have an exact optimal value. What really matters is the bio-optical model used in the R rs convergence in the first few iterations where the AC is affected by the modeled NIR reflectance. After three or four iterations, the bio-optical model becomes irrelevant since the cost function will be mostly skewed toward the SWIR bands. With every iteration, a new cost function will be calculated with new aerosol models, and the NIR-SWIR ocean radiance will be estimated until the convergence criteria in Bailey et al. (2010) is met. If the bio-optical model does not converge in sediment-dominant waters, the SW in those cases will skew all of the cost function toward the SWIR bands, thus improving the AC.
At each iteration of the AC, the following steps are performed. Assuming a smooth unimodal shape to the cost function χ 2 RH, f , the minimum of the function at which the best solution exists can be calculated as: For the sake of simplicity in the operational code implementation, the 2-d partial derivate can be converted into 1-d minimization problem by finding the minimum of the cost function for a given RH set (i.e., RH 1 ) from National Centers for Environmental Prediction (NCEP). However, ideally, the cost function minimization could be done for optical depth, finemode fraction, and relative humidity [i.e., ∇χ 2 τ , RH, f = 0]. The cost function will be minimized for a set of fine-mode fractions of aerosols. The smallest two χ 2 values in the set will be: The observed aerosol reflectance typically falls in-between the two closest aerosol reflectance, ρ a λ, RH 1 , f min 1 and ρ a λ, RH 1 , f min 2 from the model set, and are chosen based on the χ 2 min 1 f and χ 2 min 2 f , similar to the current operational algorithm. Since the fine-mode fractions in the table are discrete values with a relatively coarse step, an interpolation or a smoothing scheme of the aerosol reflectance is necessary. One way of doing this is to weight all of the aerosol models within the set for all fine-mode fractions with the χ 2 .
where n is the number of models to be mixed such that n starts from minimum χ 2 and ends with the maximum χ 2 . However, there is no strong justification to mix all models since a mixture of all the models is not what happens in the atmosphere. A more physical description is to find the most immediate models around the observations, of which χ 2 min 1 and χ 2 min 2 can be used as a weight to the aerosol reflectance, where n is 2.
This process is repeated for the second set of RH in the case that the ancillary RH falls between two values in the aerosol LUT to estimate ρ a λ, RH 2 , f mix . Finally, a linear interpolation of the two estimated aerosol reflectance is calculated as follows: The final ρ a λ, RH obs , f mix should fit through the observed aerosol reflectance in the NIR and SWIR channels, while the extrapolated reflectance to the visible channels based on the model mixtures are used for the aerosol correction.
To demonstrate how the fitting of the spectral information improves the determination of the aerosol spectral dependence, we show in Figure 3 an example of two retrievals of the aerosol reflectance for a 50% fine-mode fraction case and coarsemode dominant aerosol case, 10% fine-mode fraction, at 77.5% RH, using the 2-band and the 6-band MBAC algorithms. The analysis is done by propagating an assumed ocean reflectance for Chlor-a = 0.03 mg m −3 , from The International Ocean-Color Coordinating Group (IOCCG) report 10 for MODIS bands, to the TOA and adding the Rayleigh (+ glint) and the aerosol reflectance at one geometry (solar zenith θ 0 = 25 • , view zenith, θ = 26 • , relative azimuth angle, ϕ = 90 • ) (IOCCG, 2010). We assume an optical depth of 0.2 at 869 nm for the two aerosol cases. The summation of all these signals at TOA is the total radiance (reflectance), of which we add sensor random + systematic noise (as described in the section Error Assessment). This would simulate the expected measurement uncertainty for these specific observations. The Rayleigh (+ glint) reflectance is then removed from the total signal, and the resultant signal is the noised observed TOA reflectance (blue curve), which is the summation of the water and aerosol reflectance at TOA needed for the aerosol correction algorithms. We also show the true (input) aerosol reflectance that we aim to retrieve, shown as the black curve. The red and green curves show the retrieved aerosol reflectance from the observed (noised) reflectance, for the 2-band standard AC and the 6-band MBAC algorithm, respectively.
Since the 2-band algorithm aims to match the spectral dependence of the aerosol as a two-band ratio in the NIR, extrapolating that information to the visible (and SWIR) will depend on the radiometric quality of these two NIR bands and the linear mixing of the aerosol types in the LUT necessary for the visible correction. As can be seen from Figure 3, the spectral dependence of the retrieved aerosol reflectance using the 2-band algorithm (red curve) is similar to the observed reflectance (blue curve) within the NIR range, which is expected. However, that does not necessarily retrieve the correct reflectance in the visible and SWIR, of which the retrieval departs from the input reference reflectance (black curve) in both fine-mode fraction cases. On the FIGURE 3 | Example to demonstrate the aerosol reflectance retrievals for MODISA using the 2-band (red curve) and 6-band (green curve) algorithms. The input model reflectance (black curve) is before adding noise; meanwhile, the observed reflectance (blue curve) is the noised reflectance after removing the Rayleigh + glint reflectance (i.e., the observed reflectance is the aerosol + water signal). The left panel is for aerosol with 50% fine-mode fraction and the right panel is for 10% fine-mode fraction, and in both cases for 0.2 optical depth at 869 nm. Solar zenith θ 0 = 25 • , view zenith, θ = 26 • , relative azimuth angle, ϕ = 90 • .
other hand, the 6-band algorithm fits the spectral information for the whole NIR-SWIR range and thus is less affected by the sensor noise as more spectral information improves the chance to fit the right model. This is also demonstrated as the 6-band algorithm retrievals (green curve) overlap with the input reference model reflectance (black curve) throughout the whole spectrum. This demonstration, however, is qualitative; thus, a more thorough analysis is shown in the section Error Assessment, including the method to add noise to the reflectance.

CORRECTION FOR ABSORBING GASES
Ocean color bands are designed to avoid detection in spectral regions contaminated by strongly absorbing gases such as water vapor, oxygen, methane, and carbon dioxide. However, due to the imperfect spectral response of the detectors (i.e., broad bandwidth and out-of-band response), the measured TOA radiance is modulated by strongly absorbing gases. AC requires a correction for strongly absorbing gases as well as broadly absorbing gases, such as ozone and nitrogen dioxide that mostly absorb in the visible. The current NASA algorithm uses the method of Gordon (1995) that was implemented for SeaWiFS and then extended to all ocean color sensors. The method corrects for strongly absorbing gases, primarily for water vapor, by correcting the epsilon ratio using an empirical polynomial function for a given column gas (water vapor) concentration and air mass. To avoid this step in our proposed AC algorithm, we generate a LUT of transmittances of all strongly absorbing gases that include oxygen, water vapor, methane, and carbon dioxide using HITRAN 2012 LBL absorption database (Rothman et al., 2013;. The water vapor LUT is a function of column water vapor following the ATREM method of Gao et al. (2000). The other gases are corrected based on ancillary or climatology column concentrations, while the transmittance is adjusted for air mass per pixel. MODIS and VIIRS sensors have strong out-of-band effects from absorbing gases in the SWIR range; thus, the LUT correction method allows for a proper correction to utilize all window SWIR bands (i.e., excluding water vapor bands) in the MBAC algorithm.

ERROR ASSESSMENT Algorithm
The aerosol LUT in the MBAC algorithm is pre-computed at discretized optical properties of the aerosols; thus, the interpolation step between different fine-mode fractions and relative humidity increases the uncertainty in the retrieval of the remote sensing reflectance. To assess the algorithm error, we introduce a set of aerosol optical properties that are not part of the discretized LUT. We do so by simulating the aerosol reflectance set for an RH of 77.5%, which falls between the 75% and 80% RH in the operational LUT. The aerosol properties and reflectance are then calculated using the model set of Ahmad et al. (2010) (i.e., 10 fine-mode fractions), but with an RH of 77.5%. To define the uncertainty in the AC, we estimate the root mean squared error (RMSE) in ρ w for the normalized water-leaving radiance. We do so by assuming a reference ocean reflectance from IOCCG report 10 for case I clear waters with chlorophyll concentration, Chlor-a, 0.03 mg m −3 (IOCCG, 2010). That reference reflectance is then propagated to the TOA assuming different aerosol types and optical thicknesses for different geometries to perform the uncertainty analysis. The aerosol models are based on Ahmad et al. (2010) operational LUTs. Figure 4 shows the polar plot of ρ w at 443 nm for all view geometries in the LUT, where 0 • relative azimuth is specular reflection plane, and for three solar angles, 10, 20, 30, and 60 • .
The aerosol models encompass all aerosol types, as defined in Ahmad et al. (2010), from fine-mode to coarse-mode dominant cases with equal weights, and the aerosol reflectances were calculated for optical depth ranging from 0.05 to 0.35 at 869 nm. The uncertainty level in the retrievals is generally smaller than the scale midpoint at 0.0015 and increases with increasing solar angle from 10 to 60 • . There are some hotspots in the geometry plot that indicated a higher uncertainty in the retrievals. These regions could be due to the decreased sensitivity of the aerosol models, where the algorithm finds it difficult to determine the type due to the ambiguity in the phase function of aerosols at certain scattering angles, or geometry midpoint interpolation errors, but further analysis is needed. The uncertainty in the aerosol reflectance at 443 nm, FIGURE 4 | (A-D) are the algorithm uncertainty in water reflectance, ρ w , at 443 nm as a function of the viewing geometry for solar zenith angles of 10, 20, 30, and 60 • , respectively. The simulations are done for all aerosol models with relative humidity RH = 77.5% at optical depth at 869 nm from 0.05 to 0.35, and the retrieval is done using models with RH = 75% and 80%.
FIGURE 5 | The normalized density histogram of the percentage difference for the aerosol reflectance and optical depth retrievals at 443 and 869 nm, respectively, and the Angstrom coefficient for all cases of different aerosol types, optical depths, and geometries. optical depth retrievals at 869 nm, and the Angstrom coefficient expressed as a percent difference are shown in Figure 5 for all simulated cases of 1,037,400 (i.e., 12 solar angle × 35 sensor angle × 19 azimuth angle × 10 models × 13 optical depth).
The mean bias in the aerosol reflectance at 443 nm is −0.04% and the standard deviation is 3.5%, while optical depth retrieval is 0.05% with a standard deviation of 1.35%, and the Angstrom coefficient retrieval shows a mean bias of −1.9% with a standard deviation of 15%. The small bias and standard deviation indicate that the algorithm can well-retrieve the optical depth; however, there is a more significant uncertainty in the Angstrom coefficient retrievals with tendency to retrieve coarser aerosol types than the input data. This could be an inherent error due to the definition of the Angstrom coefficient that assumes a power law relationship of the aerosol spectral dependence from 443 to 869 nm. Since we are including all cases of aerosol types with equal representation, the uncertainty is overestimated. In the cases at the extreme ends of the aerosol's LUT, there could be large errors in the retrievals because of the necessity to extrapolate to models that are not in the table. However, these extreme ends of the LUT are intentionally added to reduce that necessity for the extrapolation since these aerosol conditions are rarely observed based on AERONET data (Ahmad et al., 2010). Finally, note that the uncertainty increases at the highest two values of finemode fractions (i.e., f = 0.95 and f = 0.80). Since most of the aerosol absorption is in the fine-mode aerosols, these two models are more absorbing than others; thus, the difference in the single scattering albedo is >0.05 in the blue wavelengths, which can induce more than 3% errors in the aerosol reflectance at TOA by linear mixing (Abdou et al., 1997;Ahmad et al., 2010). Additionally, NIR and SWIR bands used in the cost function are less sensitive to fine-mode aerosols; thus, fine-mode dominant cases will be erroneously observed as coarser aerosols, which is the case for the operational algorithm. The only method to reduce that uncertainty is by inferring the aerosol properties at shorter wavelengths while minimizing the impact of the marine reflectance. For example, observations in the UV could improve the MBAC algorithm when non-absorbing fine-mode aerosols are present over coastal (absorbing) waters where the UV signal is negligible. OCI will measure in the UV, which will be a subject of future work.

Random Noise
The uncertainty in the TOA radiance measurements due to sensor random noise will impact the retrieval uncertainty of the ocean reflectance through influence on the AC process. The uncertainty propagation of the sensor noise is simulated here using the Monte Carlo method, where the sensor noise at a given SNR value is generated for 300 iterations for all aerosol model sets mentioned above. We performed the sensitivity analysis using the NIR 2-band AC (standard approach) and a 6-band MBAC (i.e.,for MODIS 748,859,869,1,240,1,640,and 2,130 nm). The theoretical analysis for MODIS-Aqua is done for the six bands; however, we drop 1,640 nm in our real retrievals because of defective detectors at that band. The noise model here is based on the MODIS model with the noise-equivalent radiance of: where C 0 (λ) and C 1 (λ) are linear fit coefficients of the noise model from Xiong et al. (2010), and S(λ) is the spectraldependent spatial weight to bring all bands to a common spatial resolution since MODIS land bands have finer spatial resolution than ocean bands. We also performed an analysis at PACE-OCI radiometric performance levels, where the MODIS noise model was scaled to OCI noise levels such that the noise equivalent radiance, NE L, response to the dynamic range is consistent for both instruments. The random noise for each iteration is generated as follows: and the Gaussian random noise is generated as: The result of the analysis shown in Figures 6A,B indicates that the geometrically averaged uncertainty in reflectance, ρ w , is closer to 0.0016 levels for the 6-band MBAC algorithm, while it is slightly higher near 0.0021 levels for the 2-band AC for MODIS-Aqua noise levels. The improvement in the uncertainty between the two algorithms is shown in Figure 6C as the ratio of ρ w (6 − band) to the ρ w (2 − band) AC. The average reduction in uncertainty is nearly 24% using the 6-band MBAC for the MODIS-Aqua sensor; however, when the same analysis for the 6-band MBAC is done using the expected OCI performance, the average reduction in the uncertainty approaches 29% as shown in Figure 6D. The relatively modest improvement for MODIS-Aqua using the MBAC algorithm is attributed to the reduced performance of the SWIR radiometric detection. The high NEL(λ) for MODIS SWIR channels reduces the weight in the cost function in determining the best aerosol models, relative to the higher radiometric quality of the NIR bands. The MBAC algorithm, however, shows more improvement for the OCI type of radiometric performance due to the increased SNR of the SWIR detectors relative to its NIR channels, although the visible performance is kept the same as MODIS for relative comparison. It is apparent that the merit of the MBAC algorithm will prevail when the estimated aerosol reflectance at the NIR bands is compromised. In fact, that is the case for nearly 10% of the global ocean, where the non-negligible reflectance from the turbid/coastal ocean leads to systematic underestimation and often negative reflectance retrievals due to overcorrection for aerosols. More details will be discussed in the section Application to Ocean Color Instrument.

Systematic Errors
To have a complete understanding of the potential sources of uncertainty for the MBAC algorithm, we also performed an analysis to understand the impact of systematic uncertainty on the retrievals. Systematic errors in measurements are very difficult to characterize post-launch due to the lack of an accurate absolute calibration apparatus on-orbit. Typically, a solar diffuser is used as a calibration reference; however, due to the immediate degradation of the diffuser, the absolute calibration with a high degree of certainty cannot generally be achieved. Thus, the ocean color community depends on a system vicarious calibration, which aims to remove any systematic bias due to the sensor or algorithm errors. For all ocean color sensors processed by NASA GSFC, we use the vicarious calibration method by Franz et al. (2007) using the MOBY in situ measurements as ground truth. The vicarious calibration procedure reduces the systematic uncertainty to 0.5% for all bands relative to the 869nm band for MODIS sensors. However, the absolute calibration at 869 nm and SWIR bands is assumed to be 2%. Typically, systematic errors will depend on the environment and geometry at which the instrument is observing the Earth. Systematic bias can be influenced by various observing conditions and sensor characteristics, such as the operating temperature, instrument polarization sensitivity, stray light, and optical or electronic FIGURE 6 | (A,B) are the algorithm + random uncertainty in water reflectance, ρ w , at 443 nm as a function of the viewing geometry for solar zenith angles of 30 • , using six bands MBAC and two NIR bands AC, respectively. (C,D) are the ratio between the reflectance uncertainty using six bands AC with MODIS-Aqua noise levels, and PACE OCI noise levels. The simulations are done for all aerosol models with relative humidity RH = 77.5% at optical depth at 869 nm from 0.05 to 0.35, and the retrieval is done using models with RH = 75 and 80%.
cross-talk. The behavior of the errors on a per-pixel basis cannot be easily predicted; thus, it remains unknown; however, the bias can be reduced through prelaunch characterization and derived correction algorithms (e.g., polarization correction of the instrument, correction for out-of-band response, estimation, and correction for straylight). The spectral correlation of systematic bias can play an important role in determining the uncertainty of the retrievals. For example, the spectral correlation of the two-band ratio algorithm will suffer when the systematic bias is random and independent at the two bands. Since such a spectral covariance in bias cannot be directly measured on-orbit, we assume a conservative scenario in our analysis by assuming a 10% correlation inter-bands, although the spectral correlation between bands could increase/decrease depending on the spectral distance. We ran the MC analysis assuming a random bias at the levels mentioned above for MODIS-Aqua observations. The covariance matrix, Cov (λ n , λ m ), of the system errors is equal to 0.1 for off-diagonal elements and 1 for diagonal ones. The correlated random bias is then calculated as follows: The bias is generated for the number of iterations in the MC analysis for each band. U(−1, 1) is uniform distribution of random numbers from −1 to 1, while L (λ n , λ m ) is Lower Cholesky factorization of the covariance matrix such that: L * (λ n , λ m ) is the conjugate value of the matrix and D(λ n , λ m ) is the diagonal matrix. The Cholesky matrix allows the generation of random errors with inter-bands correlation specified in the covariance matrix. The bias levels are then re-adjusted to match the level of the instrument. In Figures 7A,B, the uncertainty, ρ w , at 443 nm is due to algorithm + random + systematic errors for the 6-band MBAC and NIR 2-band AC, respectively. The average uncertainty for the 2-band AC is 0.0036, while for the 6-band MBAC algorithm, the average uncertainty is 0.0024, with an improvement of 33% over the 2-band standard algorithm as shown in Figure 7C. For the OCI configuration, at higher SNRs, the improvement is only slight at 37% relative to the MODIS-Aqua configuration in Figure 7D. This can be attributed to the dominant impact of systematic uncertainty over random noise; however, in practice, OCI is expected to have lower systematic errors at 1.6% levels, FIGURE 7 | (A,B) are the algorithm + random + systematic uncertainty in water reflectance, ρ w , at 443 nm as a function of the viewing geometry for solar zenith angles of 30 • , using six bands MBAC and NIR two bands AC, respectively. (C,D) are the ratio between the reflectance uncertainty using 6-bands AC with MODIS-Aqua noise levels and PACE OCI noise levels. The simulations are done for all aerosol models with relative humidity RH = 77.5% at optical depth at 869 nm from 0.05 to 0.35, and the retrieval is done using models with RH = 75 and 80%.
which is not simulated here. Additionally, the high SNR on OCI's SWIR will have a more significant impact when retrieving in coastal waters, when the NIR signal is contaminated by nonnegligible water-leaving radiance.

APPLICATION TO OCEAN COLOR INSTRUMENT
Unlike MODIS or VIIRS, the SWIR detection assembly on OCI will include channels specifically tuned for good radiometric performance over the relatively dark ocean. This significant advancement will provide an opportunity to improve ocean color observations in complex water environments (i.e., coastal ocean) due to the increased sensitivity to the aerosols and less to the water signal. With the MBAC algorithm, the combined utilization of all SWIR channels (excluding water vapor bands) and NIR channels will offer an improved AC over coastal waters and open ocean by utilizing more spectral bands to reduce the effect of sensor noise and with adaptive weighting toward the bands that are less influenced by non-negligible waterleaving radiance. Figure 8 shows the spectral uncertainty, ρ w , for proxy OCI wavelengths from the UV to red wavelengths for open ocean and coastal turbid waters (with high total suspended material, TSM). The open ocean reflectance was simulated with Hydrolight assuming Chlor-a = 0.03 mg/m 3 and TSM = 0 g m −3 while turbid waters assume TSM = 10 g m −3 that has been propagated to the TOA for each atmospheric composition and geometry [more details are in Ibrahim and McKinna (2018)]. The retrieval was done with the 2-band NIR-only algorithm and the SWIRonly MBAC algorithm. The uncertainty is calculated for the aerosol set as in the section Error Assessment except for one specific geometry in the operational LUTs (i.e., solar zenith, θ 0 = 25 • , view zenith, θ = 26 • , relative azimuth angle, ϕ = 90 • ). For an OCI-like instrument, the SWIR MBAC algorithm shows a slight reduction in the uncertainty as compared to the 2band NIR AC over open ocean conditions. In contrast, over turbid waters, the non-negligible water-leaving radiance at the NIR bands shows a large uncertainty (systematic bias) in the reflectance retrievals from the 2-band NIR AC, which is substantially reduced using the SWIR MBAC algorithm. As previously discussed, the SW, SW(λ), in the cost function is based on the NIR iteration scheme of Bailey et al. (2010) to reduce the impact of contaminated NIR bands in the AC over productive waters. The algorithm sequentially iterates over a bio-optical model until a convergence criterion is met. The SW is directly dependent on the convergence process such that the weight on FIGURE 8 | The spectral uncertainty, ρ w , for a proxy OCI wavelengths from the UV to red wavelengths for open ocean and coastal turbid waters. The uncertainty is calculated at one specific geometry (i.e., solar zenith θ 0 = 25 • , view zenith, θ = 26 • , relative azimuth angle, ϕ = 90 • ). The analysis assumes a 2% systematic bias at NIR and SWIR bands.
the contaminated NIR bands is quickly (exponentially) reduced with increasing difficulty in the bio-optical model fit. The SW is proportional to the number of iterations; thus, it translates linearly. At the maximum of the iterations range, the SW is set to zero on the NIR bands; thus, the algorithm transitions into a full SWIR MBAC. This adaptive weight allows for a smooth transition in the processing between clear waters and productive waters; thus, it does not require switching between two or more algorithms. This allows for producing more consistent products, while maximizing the information content from all observations at all bands.

APPLICATION TO MODIS RETRIEVALS
The MBAC algorithm has been implemented into NASA's L2gen processing code. L2gen (Level-2 generator) within the SeaDAS software package is the multi-sensor Level-1 to Level-2 processing code developed and maintained by NASA's Ocean Biology Processing Group (OBPG) that is capable of retrieving ocean color products from TOA radiances for a host of sensors. L2gen supports multiple AC methods and variations that can be applied to a variety of ocean color sensors (Gordon and Wang, 1994;Ruddick et al., 2000;Wang et al., 2009;Ahmad et al., 2010;Bailey et al., 2010;. Figure 9 shows the Chlor-a retrieval from MODISA observations for both the operational algorithm and the MBAC algorithm (748, 859, 869, 1,240, and 2,130 nm, excluding 1,640 nm for Aqua because of bad detectors), the average percent difference (APD), and the SW at NIR over the East Coast of the United States. Qualitatively, there is general agreement between the two algorithms for the Chlor-a product and through the APD histogram. This agreement confirms the feasibility of implementation in an operational environment. Differences are evident, however, in coastal waters such as the turbid Chesapeake Bay (39.5 • -37 • N, 75.5 • -76 • W) or the shallow Pamlico Sound (35 • N, 76 • W) as shown in the APD figure. There are some apparent artifacts in the upper part of the scene from the MBAC retrievals that could be due to strong NIR/SWIR calibration artifacts at the edge of the scan; however, a further assessment of SWIR calibration for MODIS is necessary for reliable retrievals. The SW map shows lower values (i.e., closer to 0) near the coast where the Bailey NIR iterative algorithm was triggered, thus relying more on the SWIR bands with each increased iteration. The SW is, however, higher over more clear water over the open Atlantic waters; thus, the algorithm relies on both the NIR and SWIR bands. We believe that the MBAC algorithm should be validated in the more common moderate turbid waters and highly productive waters. That scene includes regions with high turbidity such as the upper bay regions and inland rivers as well as shallow water regions such as Pamlico Sound, highly productive waters such as middle Chesapeake Bay region, highly absorbing waters such as Long Island Sound, and the less productive to more clear waters in the open Atlantic waters. Testing the MBAC algorithm over a large dynamic range of water conditions and more complex atmospheric conditions (i.e., including the large range of Angstrom coefficients from coarse to fine aerosols) provides more insight.
We also show in Figure 10 the aerosol's Angstrom coefficient retrieval over these coastal regions for both the operational and the MBAC algorithm.
The operational algorithm shows an increased Angstrom over turbid or shallow regions near the coast that can be attributed to the failure in the NIR iterative algorithm to mitigate the non-negligible water-leaving radiance in the NIR bands. This artifact in the aerosol properties typically leads to the overcorrection of aerosol radiance and the potential retrieval of negative radiances. However, using the adaptive SW in the MBAC algorithm, it quickly damps that NIR contamination by decreasing the weight of the NIR bands in the cost function, and thus, the Angstrom coefficient is reduced to a range that agrees better with the surrounding pixels. Notice that the southwesternmost part of the scene shows an increase in the Angstrom coefficient with the MBAC algorithm, which could be a more realistic retrieval or an artifact due to a bias in the SWIR calibration. However, this could indicate that the MBAC algorithm is sensitive to fine-mode as much as coarse-mode aerosols. Better quantitative evidence is further discussed in the next section.

VALIDATION
Using NASA's SeaBASS dataset of in situ R rs measurements, we performed match-up to MODIS retrievals at eight bands in the visible spectral range near 412, 443, 488, 531, 547, 555, 667, and 678 nm, using both the 2-band operation algorithm and the MBAC algorithm for MODIS-Aqua. Many of the SeaBASS data points are coastal, however less frequent than AERONET-OC data. Table 1 below shows the spectral uncertainty, R rs (λ), R 2 , and the mean bias in the retrievals compared to in situ observations for MODISA for both the operational and the MBAC algorithm. The retrieval uncertainty, R rs (λ), for the 6band MBAC algorithm is reduced as compared to the operational algorithm, especially in the green wavelengths at 531, 547, and 555 nm, where the uncertainty is reduced by 52, 55, and 29%, respectively. There is a slight increase in uncertainty at 678 nm, FIGURE 10 | MODISA retrieval of the aerosols' Angstrom coefficient of a scene over the east coast of the US using the operational 2-band and the MBAC algorithm. which could be due to the lack of high-quality in situ observations at longer wavelengths. The coefficient of determination, however, is slightly lower for 412 and 443 for the MBAC algorithm relative to the operational one. There is a more significant correlation for the green channels, 531, 547, and 555, using the MBAC algorithm. Overall, the uncertainty is reduced using the MBAC algorithm due to the decreased mean bias in retrievals.
Because of the extrapolation of the aerosol information from the NIR/SWIR to the visible, at shorter wavelengths, the retrievals uncertainty will increase specifically in coastal waters where the ocean blue reflectance is small. Thus, retrievals in the blue bands are slightly noisier due to the impact of low-SNR SWIR bands due to the spectral extrapolation of the selected aerosol models. A further investigation into a proper estimation of systematic uncertainty and its spectral correlation is necessary to provide uncertainties that are similar to in situ validation. It is important to note that the validation analysis has been done with an initial SWIR vicarious calibration following Franz et al. (2007). The reduction of systematic bias is essential to provide reliable MBAC retrievals, which will require a thorough assessment of the SWIR calibration in a future study.

CONCLUSION
The MBAC algorithm shows several potential merits over NASA's current operational Gordon and Wang (1994) 2-band ratio aerosols correction algorithm. The utilization of multibands from the NIR to SWIR in the AC reduces the uncertainty in the retrieval of ocean color reflectance due to sensor random and systematic noise. The adaptive SW reduces the weight on the NIR bands over turbid or highly productive waters; thus, the contamination from bright ocean signal in the NIR is reduced. This algorithm enables the production of spatially consistent ocean color data in turbid and open ocean conditions. The limitations of the algorithm are mainly tied to the radiometric performance in the SWIR bands. A significantly decreased performance in the SWIR for some sensors leads to noisy retrievals in the coastal regions. For the OCI on PACE, the SWIR bands are expected to have a much higher SNR compared to MODISA within the radiance range observed over oceans; thus, a significant improvement using the MBAC algorithm is expected, especially in coastal regions. The sensitivity analysis shows an overall reduction in the reflectance uncertainty by 33% over open ocean conditions and over 300% over turbid water conditions in the blue spectral region. In addition, the SeaBASS in situ match-ups show an improvement in the retrievals using the MBAC algorithm relative to the 2-band algorithm, especially in the green spectral range where the reduction in the uncertainty is nearly 50% at 531 nm. An improvement in the system vicarious calibration of the SWIR bands for current ocean color sensors shall increase the value of the MBAC algorithm in producing more consistent products.

AUTHOR CONTRIBUTIONS
AI performed the analysis in the manuscript for the sections Error Assessment, Application to Ocean Color Instrument, Application to MODIS Retrievals, and Validation. BF contributed to editing and guiding the manuscript discussion. ZA provided Figure 1 and the aerosol LUTs. SB provided the SWIR vicarious calibration gains necessary for the MBAC algorithm.
FUNDING Support for this work was provided by a NASA Ocean Biology and Biogeochemistry PACE Science Team and Aqua, Terra, and Sumi-NPP Science Team award.